ROOT_528-00b_version: tutorials/fit/TestBinomial.C Source File

00001 // Perform a fit to a set of data with binomial errors 
00002 // like those derived from the division of two histograms. 
00003 // Three different fits are performed and compared:  
00004 //
00005 //   -  simple least square fit to the divided histogram obtained 
00006 //      from TH1::Divide with option b
00007 //   -  least square fit to the TGraphAsymmErrors obtained from 
00008 //      TGraphAsymmErrors::BayesDivide
00009 //   -  likelihood fit performed on the dividing histograms using 
00010 //      binomial statistics with the TBinomialEfficiency class
00011 // 
00012 // The first two methods are biased while the last one  is statistical correct. 
00013 // Running the script passing an integer value n larger than 1, n fits are 
00014 // performed and the bias are also shown.   
00015 // To run the script : 
00016 // 
00017 //  to show the bias performing 100 fits for 1000 events per "experiment"
00018 //  root[0]: .x TestBinomial.C+
00019 // 
00020 //  to show the bias performing 100 fits for 1000 events per "experiment"
00021 //           .x TestBinomial.C+(100, 1000)
00022 // 
00023 //
00024 #include "TBinomialEfficiencyFitter.h"
00025 #include "TVirtualFitter.h"
00026 #include "TH1.h"
00027 #include "TRandom3.h"
00028 #include "TF1.h"
00029 #include "TStyle.h"
00030 #include "TCanvas.h"
00031 #include "TLegend.h"
00032 #include "TPaveStats.h"
00033 #include <cassert>
00034 #include <iostream>
00035 
00036 void TestBinomial(int nloop = 100, int nevts = 100, bool plot=false)
00037 {
00038    gStyle->SetMarkerStyle(20);
00039    gStyle->SetLineWidth(2.0);
00040    gStyle->SetOptStat(11);
00041 
00042    TObjArray hbiasNorm;
00043    hbiasNorm.Add(new TH1D("h0Norm", "Bias Histogram fit",100,-5,5));
00044    hbiasNorm.Add(new TH1D("h1Norm","Bias Binomial fit",100,-5,5));
00045    TObjArray hbiasThreshold;
00046    hbiasThreshold.Add(new TH1D("h0Threshold", "Bias Histogram fit",100,-5,5));
00047    hbiasThreshold.Add(new TH1D("h1Threshold","Bias Binomial fit",100,-5,5));
00048    TObjArray hbiasWidth;
00049    hbiasWidth.Add(new TH1D("h0Width", "Bias Histogram fit",100,-5,5));
00050    hbiasWidth.Add(new TH1D("h1Width","Bias Binomial fit",100,-5,5));
00051    TH1D* hChisquared = new TH1D("hChisquared", 
00052       "#chi^{2} probability (Baker-Cousins)", 200, 0.0, 1.0);
00053 
00054    //TVirtualFitter::SetDefaultFitter("Minuit2");
00055 
00056    // Note: in order to be able to use TH1::FillRandom() to generate
00057    //       pseudo-experiments, we use a trick: generate "selected"
00058    //       and "non-selected" samples independently. These are
00059    //       statistically independent and therefore can be safely
00060    //       added to yield the "before selection" sample.
00061 
00062 
00063    // Define (arbitrarily?) a distribution of input events.
00064    // Here: assume a x^(-2) distribution. Boundaries: [10, 100].
00065 
00066    Double_t xmin =10, xmax = 100;
00067    TH1D* hM2D = new TH1D("hM2D", "x^(-2) denominator distribution", 
00068       45, xmin, xmax);
00069    TH1D* hM2N = new TH1D("hM2N", "x^(-2) numerator distribution",   
00070       45, xmin, xmax);
00071    TH1D* hM2E = new TH1D("hM2E", "x^(-2) efficiency",               
00072       45, xmin, xmax);
00073 
00074    TF1*  fM2D = new TF1("fM2D", "(1-[0]/(1+exp(([1]-x)/[2])))/(x*x)", 
00075       xmin, xmax);
00076    TF1*  fM2N = new TF1("fM2N", "[0]/(1+exp(([1]-x)/[2]))/(x*x)",     
00077       xmin, xmax);
00078    TF1*  fM2Fit = new TF1("fM2Fit", "[0]/(1+exp(([1]-x)/[2]))",       
00079       xmin, xmax);
00080    TF1*  fM2Fit2 = 0;
00081 
00082    TRandom3 rb(0);
00083 
00084    // First try: use a single set of parameters.
00085    // For each try, we need to find the overall normalization
00086 
00087    Double_t norm = 0.80;
00088    Double_t threshold = 25.0;
00089    Double_t width = 5.0;
00090 
00091    fM2D->SetParameter(0, norm);
00092    fM2D->SetParameter(1, threshold);
00093    fM2D->SetParameter(2, width);
00094    fM2N->SetParameter(0, norm);
00095    fM2N->SetParameter(1, threshold);
00096    fM2N->SetParameter(2, width);
00097    Double_t integralN = fM2N->Integral(xmin, xmax);
00098    Double_t integralD = fM2D->Integral(xmin, xmax);
00099    Double_t fracN = integralN/(integralN+integralD);
00100    Int_t nevtsN = rb.Binomial(nevts, fracN);
00101    Int_t nevtsD = nevts - nevtsN;
00102 
00103    gStyle->SetOptFit(1111);
00104 
00105    // generate many times to see the bias
00106    for (int iloop = 0; iloop < nloop; ++iloop) {
00107 
00108      // generate pseudo-experiments
00109      hM2D->Reset();
00110      hM2N->Reset();
00111      hM2D->FillRandom(fM2D->GetName(), nevtsD);
00112      hM2N->FillRandom(fM2N->GetName(), nevtsN);
00113      hM2D->Add(hM2N);
00114 
00115      // construct the "efficiency" histogram
00116      hM2N->Sumw2();
00117      hM2E->Divide(hM2N, hM2D, 1, 1, "b");
00118 
00119      // Fit twice, using the same fit function.
00120      // In the first (standard) fit, initialize to (arbitrary) values.
00121      // In the second fit, use the results from the first fit (this
00122      // makes it easier for the fit -- but the purpose here is not to
00123      // show how easy or difficult it is to obtain results, but whether
00124      // the CORRECT results are obtained or not!).
00125 
00126      fM2Fit->SetParameter(0, 0.5);
00127      fM2Fit->SetParameter(1, 15.0);
00128      fM2Fit->SetParameter(2, 2.0);
00129      fM2Fit->SetParError(0, 0.1);
00130      fM2Fit->SetParError(1, 1.0);
00131      fM2Fit->SetParError(2, 0.2);
00132      TCanvas* cEvt;
00133      if (plot) {
00134        cEvt = new TCanvas(Form("cEnv%d",iloop),
00135                           Form("plots for experiment %d", iloop),
00136                           1000, 600);
00137        cEvt->Divide(1,2);
00138        cEvt->cd(1);
00139        hM2D->DrawCopy("HIST");
00140        hM2N->SetLineColor(kRed);
00141        hM2N->DrawCopy("HIST SAME");
00142        cEvt->cd(2);
00143      }
00144      for (int fit = 0; fit < 2; ++fit) {
00145        Int_t status = 0;
00146        switch (fit) {
00147        case 0:
00148          hM2E->Fit(fM2Fit, "RN");
00149          if (plot) {
00150            hM2E->DrawCopy("E");
00151            fM2Fit->DrawCopy("SAME");
00152          }
00153          break;
00154        case 1:
00155          if (fM2Fit2) delete fM2Fit2;
00156          fM2Fit2 = dynamic_cast<TF1*>(fM2Fit->Clone("fM2Fit2"));
00157          if (fM2Fit2->GetParameter(0) >= 1.0)
00158            fM2Fit2->SetParameter(0, 0.95);
00159          fM2Fit2->SetParLimits(0, 0.0, 1.0);
00160          TBinomialEfficiencyFitter bef(hM2N, hM2D);
00161          status = bef.Fit(fM2Fit2,"RI");
00162          if (status!=0) {
00163            std::cerr << "Error performing binomial efficiency fit, result = "
00164                      << status << std::endl;
00165            continue;
00166          }
00167          if (plot) {
00168            fM2Fit2->SetLineColor(kRed);
00169            fM2Fit2->DrawCopy("SAME");
00170          }
00171        }
00172 
00173        if (status != 0) break;
00174 
00175        Double_t fnorm = fM2Fit->GetParameter(0);
00176        Double_t enorm = fM2Fit->GetParError(0);
00177        Double_t fthreshold = fM2Fit->GetParameter(1);
00178        Double_t ethreshold = fM2Fit->GetParError(1);
00179        Double_t fwidth = fM2Fit->GetParameter(2);
00180        Double_t ewidth = fM2Fit->GetParError(2);
00181        if (fit == 1) {
00182          fnorm = fM2Fit2->GetParameter(0);
00183          enorm = fM2Fit2->GetParError(0);
00184          fthreshold = fM2Fit2->GetParameter(1);
00185          ethreshold = fM2Fit2->GetParError(1);
00186          fwidth = fM2Fit2->GetParameter(2);
00187          ewidth = fM2Fit2->GetParError(2);
00188          hChisquared->Fill(fM2Fit2->GetProb());
00189        }
00190 
00191        TH1D* h = dynamic_cast<TH1D*>(hbiasNorm[fit]);
00192        h->Fill((fnorm-norm)/enorm);
00193        h = dynamic_cast<TH1D*>(hbiasThreshold[fit]);
00194        h->Fill((fthreshold-threshold)/ethreshold);
00195        h = dynamic_cast<TH1D*>(hbiasWidth[fit]);
00196        h->Fill((fwidth-width)/ewidth);
00197      }
00198    }
00199    
00200 
00201    TCanvas* c1 = new TCanvas("c1",
00202       "Efficiency fit biases",10,10,1000,800); 
00203    c1->Divide(2,2);
00204 
00205    TH1D *h0, *h1;
00206    c1->cd(1);
00207    h0 = dynamic_cast<TH1D*>(hbiasNorm[0]);
00208    h0->Draw("HIST");
00209    h1 = dynamic_cast<TH1D*>(hbiasNorm[1]);
00210    h1->SetLineColor(kRed);
00211    h1->Draw("HIST SAMES");
00212    TLegend* l1 = new TLegend(0.1, 0.75, 0.5, 0.9, 
00213       "plateau parameter", "ndc");
00214    l1->AddEntry(h0, Form("histogram: mean = %4.2f RMS = \
00215       %4.2f", h0->GetMean(), h0->GetRMS()), "l");
00216    l1->AddEntry(h1, Form("binomial : mean = %4.2f RMS = \
00217       %4.2f", h1->GetMean(), h1->GetRMS()), "l");
00218    l1->Draw();
00219 
00220    c1->cd(2);
00221    h0 = dynamic_cast<TH1D*>(hbiasThreshold[0]);
00222    h0->Draw("HIST");
00223    h1 = dynamic_cast<TH1D*>(hbiasThreshold[1]);
00224    h1->SetLineColor(kRed);
00225    h1->Draw("HIST SAMES");
00226    TLegend* l2 = new TLegend(0.1, 0.75, 0.5, 0.9, 
00227       "threshold parameter", "ndc");
00228    l2->AddEntry(h0, Form("histogram: mean = %4.2f RMS = \
00229       %4.2f", h0->GetMean(), h0->GetRMS()), "l");
00230    l2->AddEntry(h1, Form("binomial : mean = %4.2f RMS = \
00231       %4.2f", h1->GetMean(), h1->GetRMS()), "l");
00232    l2->Draw();
00233 
00234    c1->cd(3);
00235    h0 = dynamic_cast<TH1D*>(hbiasWidth[0]);
00236    h0->Draw("HIST");
00237    h1 = dynamic_cast<TH1D*>(hbiasWidth[1]);
00238    h1->SetLineColor(kRed);
00239    h1->Draw("HIST SAMES");
00240    TLegend* l3 = new TLegend(0.1, 0.75, 0.5, 0.9, "width parameter", "ndc");
00241    l3->AddEntry(h0, Form("histogram: mean = %4.2f RMS = \
00242       %4.2f", h0->GetMean(), h0->GetRMS()), "l");
00243    l3->AddEntry(h1, Form("binomial : mean = %4.2f RMS = \
00244       %4.2f", h1->GetMean(), h1->GetRMS()), "l");
00245    l3->Draw();
00246 
00247    c1->cd(4);
00248    hChisquared->Draw("HIST");
00249 }