test/utilities/SimplifiedElectronTagAndProbe.C

0001 #define LUMINOSITY 0.874 //(in pb^-1)
0002 #define NBINSPASS 24
0003 #define NBINSFAIL 6
0004 #include <iostream>
0005 #include <iomanip>
0006 #include <fstream>
0007
0008 ofstream effTextFile("efficiency.txt");
0009
0010
0011
0012 void SimplifiedElectronTagAndProbe(){
0013
0014 /////// Are you running over data or MC ???
0015    //TCut preCut("mcTrue");
0016     TCut preCut("");
0017     TCut cleanSC("probe_hadronicOverEm<0.15");
0018
0019    TCut cleanGsf("(probe_gsfEle_HoverE<0.15) && (probe_gsfEle_HoverE<0.15)");
0020    TCut ID95("probe_passingId");
0021    TCut NotID95(cleanGsf && "!probe_passingId");
0022    TCut NotPPass("!probe_passing");
0023    TCut PassAll("probe_passingALL");
0024    TCut NotPassAll("!probe_passingALL");
0025    TCut ID80("probe_passingId80");
0026    TCut NotID80(cleanGsf && "!probe_passingId80");
0027    TCut EMINUS("probe_gsfEle_charge<0");
0028    TCut EMINUSSC("tag_gsfEle_charge>0");
0029    TCut EPLUS("probe_gsfEle_charge>0");
0030    TCut EPLUSSC("tag_gsfEle_charge<0");
0031    TCut BARREL("abs(probe_sc_eta)<1.4442");
0032    TCut BARRELSC("abs(probe_eta)<1.4442");
0033    TCut ENDCAPS("abs(probe_sc_eta)>1.566");
0034    TCut ENDCAPSSC("abs(probe_eta)>1.566");
0035
0036
0037 // //////////////////////////////////////////////////////////
0038    effTextFile << "probe type" << "         efficiency " << "       Npass" <<
0039       "       Nfail" << endl;
0040 // //////////////////////////////////////////////////////////
0041
0042
0043
0044 // //////////////////////////////////////////////////////////
0045 //   //  Super cluster --> gsfElectron efficiency
0046 // //////////////////////////////////////////////////////////
0047
0048    TCut GsfPass = preCut && cleanGsf;
0049    TCut GsfFail = cleanSC && preCut && NotPPass;
0050    ComputeEfficiency("Gsf", GsfPass, GsfFail);
0051    ComputeEfficiency("Gsf", GsfPass, GsfFail, "Barrel", BARREL, BARRELSC);
0052    ComputeEfficiency("Gsf", GsfPass, GsfFail, "Endcap", ENDCAPS, ENDCAPSSC);
0053    ComputeEfficiency("Gsf", GsfPass, GsfFail, "", "", "", "_eminus", EMINUS, EMINUSSC);
0054    ComputeEfficiency("Gsf", GsfPass, GsfFail, "", "", "", "_eplus", EPLUS, EPLUSSC);
0055    ComputeEfficiency("Gsf", GsfPass, GsfFail, "Barrel", BARREL, BARRELSC, "_eminus", EMINUS, EMINUSSC);
0056    ComputeEfficiency("Gsf", GsfPass, GsfFail, "Barrel",  BARREL, BARRELSC, "_eplus", EPLUS, EPLUSSC);
0057    ComputeEfficiency("Gsf", GsfPass, GsfFail, "Endcap", ENDCAPS, ENDCAPSSC, "_eminus", EMINUS, EMINUSSC);
0058    ComputeEfficiency("Gsf", GsfPass, GsfFail, "Endcap", ENDCAPS, ENDCAPSSC, "_eplus", EPLUS, EPLUSSC);
0059    ComputeEfficiency("Gsf", GsfPass, GsfFail, "EndcapPlus", "probe_sc_eta>1.566", "probe_eta>1.566");
0060    ComputeEfficiency("Gsf", GsfPass, GsfFail, "EndcapMinus", "probe_sc_eta<-1.566", "probe_eta<-1.566");
0061
0062    cout << "########################################" << endl;
0063
0064
0065
0066 // //////////////////////////////////////////////////////////
0067 //   //  gsfElectron --> WP-95 selection efficiency
0068 // //////////////////////////////////////////////////////////
0069
0070
0071    TCut Id95Pass = preCut && ID95;
0072    TCut Id95Fail = preCut && NotID95;
0073    ComputeEfficiency("Id95", Id95Pass, Id95Fail);
0074    ComputeEfficiency("Id95", Id95Pass, Id95Fail, "Barrel", BARREL, BARREL);
0075    ComputeEfficiency("Id95", Id95Pass, Id95Fail, "Endcap", ENDCAPS, ENDCAPS);
0076    ComputeEfficiency("Id95", Id95Pass, Id95Fail, "", "", "", "_eminus", EMINUS, EMINUS);
0077    ComputeEfficiency("Id95", Id95Pass, Id95Fail, "", "", "", "_eplus", EPLUS, EPLUS);
0078    ComputeEfficiency("Id95", Id95Pass, Id95Fail, "Barrel", BARREL, BARREL, "_eminus", EMINUS, EMINUS);
0079    ComputeEfficiency("Id95", Id95Pass, Id95Fail, "Barrel",  BARREL, BARREL, "_eplus", EPLUS, EPLUS);
0080    ComputeEfficiency("Id95", Id95Pass, Id95Fail, "Endcap", ENDCAPS, ENDCAPS, "_eminus", EMINUS, EMINUS);
0081    ComputeEfficiency("Id95", Id95Pass, Id95Fail, "Endcap", ENDCAPS, ENDCAPS, "_eplus", EPLUS, EPLUS);
0082    ComputeEfficiency("Id95", Id95Pass, Id95Fail, "EndcapPlus", "probe_sc_eta>1.566", "probe_eta>1.566");
0083    ComputeEfficiency("Id95", Id95Pass, Id95Fail, "EndcapMinus", "probe_sc_eta<-1.566", "probe_eta<-1.566");
0084
0085    cout << "########################################" << endl;
0086
0087
0088
0089 // //////////////////////////////////////////////////////////
0090 //   //  gsfElectron --> WP-80 selection efficiency
0091 // //////////////////////////////////////////////////////////
0092
0093
0094    TCut Id80Pass = preCut && ID80;
0095    TCut Id80Fail = preCut && NotID80;
0096    ComputeEfficiency("Id80", Id80Pass, Id80Fail);
0097    ComputeEfficiency("Id80", Id80Pass, Id80Fail, "Barrel", BARREL, BARREL);
0098    ComputeEfficiency("Id80", Id80Pass, Id80Fail, "Endcap", ENDCAPS, ENDCAPS);
0099    ComputeEfficiency("Id80", Id80Pass, Id80Fail, "", "", "", "_eminus", EMINUS, EMINUS);
0100    ComputeEfficiency("Id80", Id80Pass, Id80Fail, "", "", "", "_eplus", EPLUS, EPLUS);
0101    ComputeEfficiency("Id80", Id80Pass, Id80Fail, "Barrel", BARREL, BARREL, "_eminus", EMINUS, EMINUS);
0102    ComputeEfficiency("Id80", Id80Pass, Id80Fail, "Barrel",  BARREL, BARREL, "_eplus", EPLUS, EPLUS);
0103    ComputeEfficiency("Id80", Id80Pass, Id80Fail, "Endcap", ENDCAPS, ENDCAPS, "_eminus", EMINUS, EMINUS);
0104    ComputeEfficiency("Id80", Id80Pass, Id80Fail, "Endcap", ENDCAPS, ENDCAPS, "_eplus", EPLUS, EPLUS);
0105    ComputeEfficiency("Id80", Id80Pass, Id80Fail, "EndcapPlus", "probe_sc_eta>1.566", "probe_eta>1.566");
0106    ComputeEfficiency("Id80", Id80Pass, Id80Fail, "EndcapMinus", "probe_sc_eta<-1.566", "probe_eta<-1.566");
0107
0108    cout << "########################################" << endl;
0109
0110
0111
0112 // //////////////////////////////////////////////////////////
0113 //   //   WP-95 --> HLT triggering efficiency
0114 // //////////////////////////////////////////////////////////
0115
0116
0117    TCut HLT95Pass = preCut && ID95 && PassAll;
0118    TCut HLT95Fail = preCut && ID80 && NotPassAll;
0119    ComputeEfficiency("HLT95", HLT95Pass, HLT95Fail);
0120    ComputeEfficiency("HLT95", HLT95Pass, HLT95Fail, "Barrel", BARREL, BARREL);
0121    ComputeEfficiency("HLT95", HLT95Pass, HLT95Fail, "Endcap", ENDCAPS, ENDCAPS);
0122    ComputeEfficiency("HLT95", HLT95Pass, HLT95Fail, "", "", "", "_eminus", EMINUS, EMINUS);
0123    ComputeEfficiency("HLT95", HLT95Pass, HLT95Fail, "", "", "", "_eplus", EPLUS, EPLUS);
0124    ComputeEfficiency("HLT95", HLT95Pass, HLT95Fail, "Barrel", BARREL, BARREL, "_eminus", EMINUS, EMINUS);
0125    ComputeEfficiency("HLT95", HLT95Pass, HLT95Fail, "Barrel",  BARREL, BARREL, "_eplus", EPLUS, EPLUS);
0126    ComputeEfficiency("HLT95", HLT95Pass, HLT95Fail, "Endcap", ENDCAPS, ENDCAPS, "_eminus", EMINUS, EMINUS);
0127    ComputeEfficiency("HLT95", HLT95Pass, HLT95Fail, "Endcap", ENDCAPS, ENDCAPS, "_eplus", EPLUS, EPLUS);
0128    ComputeEfficiency("HLT95", HLT95Pass, HLT95Fail, "EndcapPlus", "probe_sc_eta>1.566", "probe_eta>1.566");
0129    ComputeEfficiency("HLT95", HLT95Pass, HLT95Fail, "EndcapMinus", "probe_sc_eta<-1.566", "probe_eta<-1.566");
0130
0131    cout << "########################################" << endl;
0132
0133
0134
0135
0136 // //////////////////////////////////////////////////////////
0137 //   //   WP-80 --> HLT triggering efficiency
0138 // //////////////////////////////////////////////////////////
0139
0140
0141    TCut HLT80Pass = preCut && ID80 && PassAll;
0142    TCut HLT80Fail = preCut && ID80 && NotPassAll;
0143    ComputeEfficiency("HLT80", HLT80Pass, HLT80Fail);
0144    ComputeEfficiency("HLT80", HLT80Pass, HLT80Fail, "Barrel", BARREL, BARREL);
0145    ComputeEfficiency("HLT80", HLT80Pass, HLT80Fail, "Endcap", ENDCAPS, ENDCAPS);
0146    ComputeEfficiency("HLT80", HLT80Pass, HLT80Fail, "", "", "", "_eminus", EMINUS, EMINUS);
0147    ComputeEfficiency("HLT80", HLT80Pass, HLT80Fail, "", "", "", "_eplus", EPLUS, EPLUS);
0148    ComputeEfficiency("HLT80", HLT80Pass, HLT80Fail, "Barrel", BARREL, BARREL, "_eminus", EMINUS, EMINUS);
0149    ComputeEfficiency("HLT80", HLT80Pass, HLT80Fail, "Barrel",  BARREL, BARREL, "_eplus", EPLUS, EPLUS);
0150    ComputeEfficiency("HLT80", HLT80Pass, HLT80Fail, "Endcap", ENDCAPS, ENDCAPS, "_eminus", EMINUS, EMINUS);
0151    ComputeEfficiency("HLT80", HLT80Pass, HLT80Fail, "Endcap", ENDCAPS, ENDCAPS, "_eplus", EPLUS, EPLUS);
0152    ComputeEfficiency("HLT80", HLT80Pass, HLT80Fail, "EndcapPlus", "probe_sc_eta>1.566", "probe_eta>1.566");
0153    ComputeEfficiency("HLT80", HLT80Pass, HLT80Fail, "EndcapMinus", "probe_sc_eta<-1.566", "probe_eta<-1.566");
0154
0155    cout << "########################################" << endl;
0156
0157    effTextFile.close();
0158 }
0159
0160
0161
0162
0163
0164
0165
0166
0167
0168
0169
0170 void ComputeEfficiency(char* primaryLabel, TCut primaryCutPass,
0171                        TCut primaryCutFail, char* secondaryLabel="",
0172                        TCut secondaryCutPass="", TCut secondaryCutFail="",
0173                        char* tertiaryLabel="",
0174                        TCut tertiaryCutPass = "", TCut tertiaryCutFail = "")
0175 {
0176     TFile* f = new TFile("allTPtrees_836nb.root");
0177
0178   TTree* scTree = (TTree*) f->Get("PhotonToGsf/fitter_tree");
0179   TTree* gsfTree = (TTree*) f->Get("GsfToIso/fitter_tree");
0180   char namePass[50];
0181   char nameFail[50];
0182
0183   sprintf(namePass,"Zmass%sPass%s%s",primaryLabel,secondaryLabel,tertiaryLabel);
0184   sprintf(nameFail,"Zmass%sFail%s%s",primaryLabel,secondaryLabel,tertiaryLabel);
0185   TH1F* histPass = createHistogram(namePass, 30);
0186   TH1F* histFail = createHistogram(nameFail, 12);
0187   gsfTree->Draw("1.02*mass>>"+TString(namePass), primaryCutPass && secondaryCutPass &&
0188   tertiaryCutPass,"goff");
0189   TString checkForSCTree(primaryLabel);
0190   TString plotVar = "1.03*mass>>"+TString(nameFail);
0191   if(checkForSCTree.Contains("Gsf"))
0192      scTree->Draw(plotVar,primaryCutFail && secondaryCutFail && tertiaryCutFail,"goff");
0193   else
0194      gsfTree->Draw(plotVar,primaryCutFail && secondaryCutFail && tertiaryCutFail,"goff");
0195   ComputeEfficiency(*histPass, *histFail);
0196
0197   delete histPass;
0198   delete histFail;
0199 }
0200
0201
0202
0203
0204
0205
0206 TH1F* createHistogram(char* name, int nbins=12) {
0207   TH1F* hist = new TH1F(name,name, nbins, 60, 120);
0208   hist->SetTitle("");
0209   char temp[100];
0210   sprintf(temp, "Events / %.1f GeV/c^{2}", 60./nbins);
0211   hist->GetXaxis()->SetTitle("m_{ee} (GeV/c^{2})");
0212   hist->GetYaxis()->SetTitle(temp);
0213   return hist;
0214 }
0215
0216
0217
0218
0219
0220
0221 double ErrorInProduct(double x, double errx, double y,
0222                       double erry, double corr) {
0223    double xFrErr = errx/x;
0224    double yFrErr = erry/y;
0225    return sqrt(xFrErr**2 +yFrErr**2 + 2.0*corr*xFrErr*yFrErr)*x*y;
0226 }
0227
0228
0229 void ComputeEfficiency( TH1& hist_pass, TH1& hist_fail)
0230 {
0231
0232
0233   TString effType = hist_pass.GetName();
0234   effType.ReplaceAll("Zmass","");
0235   effType.ReplaceAll("Pass","");
0236
0237   // The fit variable - lepton invariant mass
0238   RooRealVar* rooMass_ = new RooRealVar("Mass","m_{ee}",60.0, 120.0, "GeV/c^{2}");
0239   RooRealVar Mass = *rooMass_;
0240
0241   // Make the category variable that defines the two fits,
0242   // namely whether the probe passes or fails the eff criteria.
0243   RooCategory sample("sample","") ;
0244   sample.defineType("Pass", 1) ;
0245   sample.defineType("Fail", 2) ;
0246
0247
0248   gROOT->cd();
0249
0250
0251   ///////// convert Histograms into RooDataHists
0252   RooDataHist* data_pass = new RooDataHist("data_pass","data_pass",
0253                       RooArgList(Mass), &hist_pass);
0254   RooDataHist* data_fail = new RooDataHist("data_fail","data_fail",
0255                       RooArgList(Mass), &hist_fail);
0256
0257   RooDataHist* data = new RooDataHist( "fitData","fitData",
0258   RooArgList(Mass),RooFit::Index(sample),
0259   RooFit::Import("Pass",*data_pass), RooFit::Import("Fail",*data_fail) );
0260
0261
0262
0263
0264  // Signal pdf
0265
0266   Zeelineshape_file =  new TFile("Zlineshapes.root", "READ");
0267   TH1* histbbpass = (TH1D*) Zeelineshape_file->Get("pass_BB");
0268   TH1* histebpass = (TH1D*) Zeelineshape_file->Get("pass_BE");
0269   TH1* histeepass = (TH1D*) Zeelineshape_file->Get("pass_EE");
0270   TH1D* th1 = (TH1D*) histbbpass->Clone("th1");
0271   th1->Add(histebpass);
0272   th1->Add(histeepass);
0273   RooDataHist* rdh = new RooDataHist("rdh","", Mass, th1);
0274   RooHistPdf* signalShapePdf = new RooHistPdf("signalShapePdf", "",
0275   RooArgSet(Mass), *rdh);
0276
0277
0278
0279 //   RooRealVar* Mean_   = new RooRealVar("Mean","Mean", 90.0, 86.2, 95.2);
0280 //   RooRealVar* Width_  = new RooRealVar("Width","Width", 2.5);
0281 //   RooRealVar* Resolution_  = new RooRealVar("Resolution","Resolution", 5.6, 0., 10.);
0282 //   RooRealVar* MeanBifG_   = new RooRealVar("MeanBifG","MeanBifG", 87.50, 86.2, 95.2);
0283 //   RooRealVar* WidthL_  = new RooRealVar("WidthL","WidthL", 3., 2., 6.);
0284 //   RooRealVar* WidthR_  = new RooRealVar("WidthR","WidthR", 6.);
0285 //   RooRealVar* BifurGaussFrac_ = new RooRealVar("BifurGaussFrac","",0.2, 0.0, 1.0);
0286
0287 //   // Voigtian
0288 //   RooAbsPdf* voigtianPdf= new RooVoigtian("voigtianPdf", "",
0289 //   Mass, *Mean_, *Width_, *Resolution_);
0290
0291 //   RooAbsPdf* signalShapePdf= voigtianPdf;
0292
0293
0294 //   // Bifurcated Gaussian
0295 //   RooAbsPdf* bifurGaussPdf_ = new RooBifurGauss("bifurGaussPdf", "",
0296 //   Mass, *Mean_, *WidthL_, *WidthL_);
0297
0298 //   // Signal PDF
0299 //   RooAddPdf* signalShapePdf= new RooAddPdf("signalShapePdf", "",
0300 //   *voigtianPdf, *bifurGaussPdf_,*BifurGaussFrac_);
0301
0302 //   // Signal PDF --> Failing sample
0303 //   RooRealVar* ResolutionFail_  = new RooRealVar("ResolutionFail","ResolutionFail", 5.6, 0., 10.);
0304 //   RooAbsPdf* signalShapePdfFail= new RooVoigtian("signalShapePdfFail", "",
0305 //   Mass, *Mean_, *Width_, *ResolutionFail_);
0306
0307
0308   // Background PDF
0309   RooRealVar* bkgShape = new RooRealVar("bkgShape","bkgShape",-0.2,-10.,0.);
0310   RooExponential* bkgShapePdf = new RooExponential("bkgShapePdf","bkgShapePdf",Mass, *bkgShape);
0311
0312
0313   // Now define some efficiency/yield variables
0314   RooRealVar* numSignal = new RooRealVar("numSignal","numSignal", 4000.0, 0.0, 100000.0);
0315   RooRealVar* eff = new RooRealVar("eff","eff", 0.9, 0.5, 1.0);
0316   RooFormulaVar* nSigPass = new RooFormulaVar("nSigPass", "eff*numSignal", RooArgList(*eff,*numSignal));
0317   RooFormulaVar* nSigFail = new RooFormulaVar("nSigFail", "(1.0-eff)*numSignal", RooArgList(*eff,*numSignal));
0318   RooRealVar* nBkgPass = new RooRealVar("nBkgPass","nBkgPass", 1000.0, 0.0, 10000000.0);
0319   RooRealVar* nBkgFail = new RooRealVar("nBkgFail","nBkgFail", 1000.0, 0.0, 10000000.0);
0320
0321
0322   RooArgList componentsPass(*signalShapePdf,*bkgShapePdf);
0323   RooArgList yieldsPass(*nSigPass, *nBkgPass);
0324   //RooArgList componentsFail(*signalShapePdfFail,*bkgShapePdf);
0325   RooArgList componentsFail(*signalShapePdf,*bkgShapePdf);
0326   RooArgList yieldsFail(*nSigFail, *nBkgFail);
0327
0328
0329    RooAddPdf pdfPass("pdfPass","extended sum pdf", componentsPass, yieldsPass);
0330    RooAddPdf pdfFail("pdfFail","extended sum pdf", componentsFail, yieldsFail);
0331
0332
0333
0334    // The total simultaneous fit ...
0335    RooSimultaneous totalPdf("totalPdf","totalPdf", sample);
0336    totalPdf.addPdf(pdfPass,"Pass");
0337    totalPdf.Print();
0338    totalPdf.addPdf(pdfFail,"Fail");
0339    totalPdf.Print();
0340
0341
0342   // ********* Do the Actual Fit ********** //
0343    RooFitResult *fitResult = totalPdf.fitTo(*data, RooFit::Save(true),
0344    RooFit::Extended(true), RooFit::PrintLevel(-1));
0345   fitResult->Print("v");
0346
0347   double numerator = nSigPass->getVal();
0348   double nfails    = nSigFail->getVal();
0349   double denominator = numerator + nfails;
0350
0351
0352
0353
0354   // ********** Make and save Canvas for the plots ********** //
0355   gROOT->ProcessLine(".L ~/tdrstyle.C");
0356   setTDRStyle();
0357   tdrStyle->SetErrorX(0.5);
0358   tdrStyle->SetPadLeftMargin(0.19);
0359   tdrStyle->SetPadRightMargin(0.10);
0360   tdrStyle->SetPadBottomMargin(0.15);
0361   tdrStyle->SetLegendBorderSize(0);
0362   tdrStyle->SetTitleYOffset(1.5);
0363   RooAbsData::ErrorType errorType = RooAbsData::Poisson;
0364
0365   TString cname = TString("fit_") + hist_pass.GetName();
0366   TCanvas* c = new TCanvas(cname,cname,500,500);
0367   RooPlot* frame1 = Mass.frame();
0368   frame1->SetMinimum(0);
0369   data_pass->plotOn(frame1,RooFit::DataError(errorType));
0370   pdfPass.plotOn(frame1,RooFit::ProjWData(*data_pass),
0371   RooFit::Components(*bkgShapePdf),RooFit::LineColor(kRed));
0372   pdfPass.plotOn(frame1,RooFit::ProjWData(*data_pass));
0373   frame1->Draw("e0");
0374
0375   TPaveText *plotlabel = new TPaveText(0.23,0.87,0.43,0.92,"NDC");
0376    plotlabel->SetTextColor(kBlack);
0377    plotlabel->SetFillColor(kWhite);
0378    plotlabel->SetBorderSize(0);
0379    plotlabel->SetTextAlign(12);
0380    plotlabel->SetTextSize(0.03);
0381    plotlabel->AddText("CMS Preliminary 2010");
0382   TPaveText *plotlabel2 = new TPaveText(0.23,0.82,0.43,0.87,"NDC");
0383    plotlabel2->SetTextColor(kBlack);
0384    plotlabel2->SetFillColor(kWhite);
0385    plotlabel2->SetBorderSize(0);
0386    plotlabel2->SetTextAlign(12);
0387    plotlabel2->SetTextSize(0.03);
0388    plotlabel2->AddText("#sqrt{s} = 7 TeV");
0389   TPaveText *plotlabel3 = new TPaveText(0.23,0.75,0.43,0.80,"NDC");
0390    plotlabel3->SetTextColor(kBlack);
0391    plotlabel3->SetFillColor(kWhite);
0392    plotlabel3->SetBorderSize(0);
0393    plotlabel3->SetTextAlign(12);
0394    plotlabel3->SetTextSize(0.03);
0395   char temp[100];
0396   sprintf(temp, "%.4f", LUMINOSITY);
0397   plotlabel3->AddText((string("#int#font[12]{L}dt = ") +
0398   temp + string(" pb^{ -1}")).c_str());
0399   TPaveText *plotlabel4 = new TPaveText(0.6,0.82,0.8,0.87,"NDC");
0400    plotlabel4->SetTextColor(kBlack);
0401    plotlabel4->SetFillColor(kWhite);
0402    plotlabel4->SetBorderSize(0);
0403    plotlabel4->SetTextAlign(12);
0404    plotlabel4->SetTextSize(0.03);
0405    double nsig = numSignal->getVal();
0406    double nErr = numSignal->getError();
0407    double e = eff->getVal();
0408    double eErr = eff->getError();
0409    double corr = fitResult->correlation(*eff, *numSignal);
0410    double err = ErrorInProduct(nsig, nErr, e, eErr, corr);
0411    sprintf(temp, "Signal = %.2f #pm %.2f", nSigPass->getVal(), err);
0412    plotlabel4->AddText(temp);
0413   TPaveText *plotlabel5 = new TPaveText(0.6,0.77,0.8,0.82,"NDC");
0414    plotlabel5->SetTextColor(kBlack);
0415    plotlabel5->SetFillColor(kWhite);
0416    plotlabel5->SetBorderSize(0);
0417    plotlabel5->SetTextAlign(12);
0418    plotlabel5->SetTextSize(0.03);
0419    sprintf(temp, "Bkg = %.2f #pm %.2f", nBkgPass->getVal(), nBkgPass->getError());
0420    plotlabel5->AddText(temp);
0421   TPaveText *plotlabel6 = new TPaveText(0.6,0.87,0.8,0.92,"NDC");
0422    plotlabel6->SetTextColor(kBlack);
0423    plotlabel6->SetFillColor(kWhite);
0424    plotlabel6->SetBorderSize(0);
0425    plotlabel6->SetTextAlign(12);
0426    plotlabel6->SetTextSize(0.03);
0427    plotlabel6->AddText("Passing probes");
0428   TPaveText *plotlabel7 = new TPaveText(0.6,0.72,0.8,0.77,"NDC");
0429    plotlabel7->SetTextColor(kBlack);
0430    plotlabel7->SetFillColor(kWhite);
0431    plotlabel7->SetBorderSize(0);
0432    plotlabel7->SetTextAlign(12);
0433    plotlabel7->SetTextSize(0.03);
0434    sprintf(temp, "Eff = %.3f #pm %.3f", eff->getVal(), eff->getErrorHi());
0435    plotlabel7->AddText(temp);
0436 //   TPaveText *plotlabel9 = new TPaveText(0.6,0.67,0.8,0.72,"NDC");
0437 //    plotlabel9->SetTextColor(kBlack);
0438 //    plotlabel9->SetFillColor(kWhite);
0439 //    plotlabel9->SetBorderSize(0);
0440 //    plotlabel9->SetTextAlign(12);
0441 //    plotlabel9->SetTextSize(0.03);
0442 //    sprintf(temp, "Resolution = %.2f #pm %.2f", Resolution_->getVal(), Resolution_->getError());
0443 //    plotlabel9->AddText(temp);
0444 //   plotlabel->Draw();
0445 //   plotlabel2->Draw();
0446 //   plotlabel3->Draw();
0447   plotlabel4->Draw();
0448   plotlabel5->Draw();
0449   plotlabel6->Draw();
0450   plotlabel7->Draw();
0451 //   plotlabel9->Draw();
0452
0453   c->SaveAs( cname + TString(".eps"));
0454   c->SaveAs( cname + TString(".gif"));
0455   c->SaveAs( cname + TString(".root"));
0456   delete c;
0457
0458
0459   cname = TString("fit_") + hist_fail.GetName();
0460   TCanvas* c2 = new TCanvas(cname,cname,500,500);
0461   RooPlot* frame2 = Mass.frame();
0462   frame2->SetMinimum(0);
0463   data_fail->plotOn(frame2,RooFit::DataError(errorType));
0464   pdfFail.plotOn(frame2,RooFit::ProjWData(*data_fail),
0465   RooFit::Components(*bkgShapePdf),RooFit::LineColor(kRed));
0466   pdfFail.plotOn(frame2,RooFit::ProjWData(*data_fail));
0467   frame2->Draw("e0");
0468
0469   TPaveText *plotlabel = new TPaveText(0.23,0.87,0.43,0.92,"NDC");
0470    plotlabel->SetTextColor(kBlack);
0471    plotlabel->SetFillColor(kWhite);
0472    plotlabel->SetBorderSize(0);
0473    plotlabel->SetTextAlign(12);
0474    plotlabel->SetTextSize(0.03);
0475    plotlabel->AddText("CMS Preliminary 2010");
0476   TPaveText *plotlabel2 = new TPaveText(0.23,0.82,0.43,0.87,"NDC");
0477    plotlabel2->SetTextColor(kBlack);
0478    plotlabel2->SetFillColor(kWhite);
0479    plotlabel2->SetBorderSize(0);
0480    plotlabel2->SetTextAlign(12);
0481    plotlabel2->SetTextSize(0.03);
0482    plotlabel2->AddText("#sqrt{s} = 7 TeV");
0483   TPaveText *plotlabel3 = new TPaveText(0.23,0.75,0.43,0.80,"NDC");
0484    plotlabel3->SetTextColor(kBlack);
0485    plotlabel3->SetFillColor(kWhite);
0486    plotlabel3->SetBorderSize(0);
0487    plotlabel3->SetTextAlign(12);
0488    plotlabel3->SetTextSize(0.03);
0489   char temp[100];
0490   sprintf(temp, "%.4f", LUMINOSITY);
0491   plotlabel3->AddText((string("#int#font[12]{L}dt = ") +
0492   temp + string(" pb^{ -1}")).c_str());
0493   TPaveText *plotlabel4 = new TPaveText(0.6,0.82,0.8,0.87,"NDC");
0494    plotlabel4->SetTextColor(kBlack);
0495    plotlabel4->SetFillColor(kWhite);
0496    plotlabel4->SetBorderSize(0);
0497    plotlabel4->SetTextAlign(12);
0498    plotlabel4->SetTextSize(0.03);
0499    double err = ErrorInProduct(nsig, nErr, 1.0-e, eErr, corr);
0500    sprintf(temp, "Signal = %.2f #pm %.2f", nSigFail->getVal(), err);
0501    plotlabel4->AddText(temp);
0502   TPaveText *plotlabel5 = new TPaveText(0.6,0.77,0.8,0.82,"NDC");
0503    plotlabel5->SetTextColor(kBlack);
0504    plotlabel5->SetFillColor(kWhite);
0505    plotlabel5->SetBorderSize(0);
0506    plotlabel5->SetTextAlign(12);
0507    plotlabel5->SetTextSize(0.03);
0508    sprintf(temp, "Bkg = %.2f #pm %.2f", nBkgFail->getVal(), nBkgFail->getError());
0509    plotlabel5->AddText(temp);
0510   TPaveText *plotlabel6 = new TPaveText(0.6,0.87,0.8,0.92,"NDC");
0511    plotlabel6->SetTextColor(kBlack);
0512    plotlabel6->SetFillColor(kWhite);
0513    plotlabel6->SetBorderSize(0);
0514    plotlabel6->SetTextAlign(12);
0515    plotlabel6->SetTextSize(0.03);
0516    plotlabel6->AddText("Failing probes");
0517   TPaveText *plotlabel7 = new TPaveText(0.6,0.72,0.8,0.77,"NDC");
0518    plotlabel7->SetTextColor(kBlack);
0519    plotlabel7->SetFillColor(kWhite);
0520    plotlabel7->SetBorderSize(0);
0521    plotlabel7->SetTextAlign(12);
0522    plotlabel7->SetTextSize(0.03);
0523    sprintf(temp, "Eff = %.3f #pm %.3f", eff->getVal(), eff->getErrorHi(), eff->getErrorLo());
0524    plotlabel7->AddText(temp);
0525
0526 //   plotlabel->Draw();
0527 //   plotlabel2->Draw();
0528 //   plotlabel3->Draw();
0529   plotlabel4->Draw();
0530   plotlabel5->Draw();
0531   plotlabel6->Draw();
0532   plotlabel7->Draw();
0533
0534   c2->SaveAs( cname + TString(".eps"));
0535   c2->SaveAs( cname + TString(".gif"));
0536   c2->SaveAs( cname + TString(".root"));
0537   delete c2;
0538
0539   // cout << "########################################" << endl;
0540   effType.ReplaceAll("Gsf","");
0541   effType.ReplaceAll("Id95_","");
0542   effType.ReplaceAll("Id95","");
0543   effType.ReplaceAll("Id80_","");
0544   effType.ReplaceAll("Id80","");
0545   effType.ReplaceAll("HLT95_","");
0546   effType.ReplaceAll("HLT95","");
0547   effType.ReplaceAll("HLT80_","");
0548   effType.ReplaceAll("HLT80","");
0549
0550  char* effTypeToPrint = (char*) effType;
0551
0552   effTextFile << effTypeToPrint << "    "
0553        << setiosflags(ios::fixed) << setprecision(4) << eff.getVal()
0554        << " + " << eff->getErrorHi() << " - " <<  eff->getErrorLo()
0555        << setiosflags(ios::fixed) << setprecision(2)
0556        << "    " << numerator << "    "  << nfails << endl;
0557   //cout << "########################################" << endl;
0558 }
0559
0560
0561 double ErrorInProduct(double x, double errx, double y,
0562                       double erry, double corr) {
0563    double xFrErr = errx/x;
0564    double yFrErr = erry/y;
0565    return sqrt(xFrErr**2 +yFrErr**2 + 2.0*corr*xFrErr*yFrErr)*x*y;
0566 }
0567
0568
0569
0570
0571
0572
0573 void ClopperPearsonLimits(double numerator, double denominator,
0574 double &lowerLimit, double &upperLimit, const double CL_low=1.0,
0575 const double CL_high=1.0)
0576 {
0577 //Confidence intervals are in the units of \sigma.
0578
0579    double ratio = numerator/denominator;
0580
0581 // first get the lower limit
0582    if(numerator==0)   lowerLimit = 0.0;
0583    else {
0584       double v=ratio/2;
0585       double vsL=0;
0586       double vsH=ratio;
0587       double p=CL_low/100;
0588       while((vsH-vsL)>1e-5) {
0589          if(BinP(denominator,v,numerator,denominator)>p)
0590          { vsH=v; v=(vsL+v)/2; }
0591          else { vsL=v; v=(v+vsH)/2; }
0592       }
0593       lowerLimit = v;
0594    }
0595
0596 // now get the upper limit
0597    if(numerator==denominator) upperLimit = 1.0;
0598    else {
0599       double v=(1+ratio)/2;
0600       double vsL=ratio;
0601       double vsH=1;
0602       double p=CL_high/100;
0603       while((vsH-vsL)>1e-5) {
0604          if(BinP(denominator,v,0,numerator)<p) { vsH=v; v=(vsL+v)/2; }
0605          else { vsL=v; v=(v+vsH)/2; }
0606       }
0607       upperLimit = v;
0608    }
0609 }
0610
0611
0612
0613
0614 double BinP(int N, double p, int x1, int x2) {
0615    double q=p/(1-p);
0616    int k=0;
0617    double v = 1;
0618    double s=0;
0619    double tot=0.0;
0620     while(k<=N) {
0621        tot=tot+v;
0622        if(k>=x1 & k<=x2) { s=s+v; }
0623        if(tot>1e30){s=s/1e30; tot=tot/1e30; v=v/1e30;}
0624        k=k+1;
0625        v=v*q*(N+1-k)/k;
0626     }
0627     return s/tot;
0628 }
0629
0630
0631