ForestErrorAnal Lab7.C

void BGDiff_Plot(Int_t Min, Int_t Max, double p)
{
  // k= number of successful trials                                                    
  //p = probability of success                                                         
  //N= number of trials                                                                
  // Example: the number of heads = k after flipping a coin n times, p=0.5             
  TH1F *Binomial=new TH1F("Binomial","Binomial",Max-Min,Min,Max);

  Binomial->GetXaxis()->SetTitle("x");
  Binomial->GetXaxis()->SetTitleSize(0.06);
  Binomial->GetXaxis()->CenterTitle();
  Binomial->GetYaxis()->SetTitle("P(x)");
  Binomial->GetYaxis()->SetTitleOffset(1.2);
  Binomial->GetYaxis()->SetTitleSize(0.06);
  Binomial->SetLineColor(4);
  Binomial->SetLineWidth(3);
  Binomial->SetLineStyle(2);

  TH1F *Gauss=new TH1F("Gauss","Gauss",Max-Min,Min,Max);
  Gauss->SetLineColor(2);
  Gauss->SetLineWidth(3);
  Gauss->SetLineStyle(9);
  double x,Mean, sigma,Binomial_x;

  TH1F *BGDiff=new TH1F("BGDiff","BGDiff",Max-Min,Min,Max);
  TH1F *BGDiffSqrd=new TH1F("BGDiffSqrd","BGDiffSqrd",Max-Min,Min,Max);

  gStyle->SetOptStat(1111111);

  for(Int_t i=Min;i<Max;i++)
    {
      Binomial_x=ROOT::Math::binomial_pdf(i,p,Max-Min);
      Binomial->Fill(i,Binomial_x);
      //printf("Step \t %d \n",i);                                                     
      Mean=p*(Max-Min);
      sigma=sqrt((Max-Min)/4);
      x=ROOT::Math::gaussian_pdf(i,sigma,Mean);
      Gauss->Fill(i,x);

      BGDiff->Fill(i,Binomial_x-x);
      BGDiffSqrd->Fill(i,(Binomial_x-x)*(Binomial_x-x));

 }
  //   BGDiff->Draw();                                                                 
  BGDiffSqrd->Draw();
  //    Binomial->Draw();                                                              
  // Gauss->Draw("same");                                                              


}

void PGDiff_Plot(Int_t Min, Int_t Max, double p)
{
  // k= number of successful trials                                                    
  //p = probability of success                                                         
  //N= number of trials                                                                
  // Example: the number of heads = k after flipping a coin n times, p=0.5             

  TH1F *Poisson=new TH1F("Poisson","Poisson",Max-Min,Min,Max);
  Poisson->GetXaxis()->SetTitle("x");
  Poisson->GetXaxis()->SetTitleSize(0.06);
  Poisson->GetXaxis()->CenterTitle();
  Poisson->GetYaxis()->SetTitle("P(x)");
  Poisson->GetYaxis()->SetTitleOffset(1.2);
  Poisson->GetYaxis()->SetTitleSize(0.06);
  Poisson->SetLineColor(3);
  Poisson->SetLineWidth(3);
  Poisson->SetLineStyle(4);

  TH1F *Gauss=new TH1F("Gauss","Gauss",Max-Min,Min,Max);
  Gauss->SetLineColor(2);
  Gauss->SetLineWidth(3);
  Gauss->SetLineStyle(9);
  double x,Mean, sigma,Poisson_x;

  TH1F *PGDiff=new TH1F("PGDiff","PGDiff",Max-Min,Min,Max);
  TH1F *PGDiffSqrd=new TH1F("PGDiffSqrd","PGDiffSqrd",Max-Min,Min,Max);

  gStyle->SetOptStat(1111111);

  for(Int_t i=Min;i<Max;i++)
    {
      //printf("Step \t %d \n",i);                                                     
            Mean=p*(Max-Min);
            Poisson_x=ROOT::Math::poisson_pdf(i,Mean);
            Poisson->Fill(i,x);
            sigma=sqrt(Mean);

            x=ROOT::Math::gaussian_pdf(i,sigma,Mean);
            Gauss->Fill(i,x);
            PGDiff->Fill(i,Poisson_x-x);
            PGDiffSqrd->Fill(i,(Poisson_x-x)*(Poisson_x-x));

 }
  //   Poisson->Draw();                                                                
  //   Gauss->Draw("same");                                                            
  PGDiffSqrd->Draw();


}