Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:33:22

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