Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:16:45 

0001 #include <iostream>
0002 #include <iomanip>
0003 
0004 
0005 
0006 void CutNCountElectronTagAndProbe(){
0007 
0008 /////// Are you running over data or MC ???
0009     TCut preCut("mcTrue");
0010     //TCut preCut("");
0011 
0012    TCut cleanSC("(( (abs(probe_eta)<1.5)"
0013    "&& ((probe_trkSumPtHollowConeDR03 + probe_ecalRecHitSumEtConeDR03 + probe_hcalTowerSumEtConeDR03)/probe_et < 0.07)"
0014    " && (probe_sigmaIetaIeta<0.01)"
0015    " && (probe_hadronicOverEm<0.15)"
0016    ")"
0017    " || ((abs(probe_eta)>1.5)"
0018    "&& ( (probe_trkSumPtHollowConeDR03 + probe_ecalRecHitSumEtConeDR03 + probe_hcalTowerSumEtConeDR03)/probe_et < 0.06 )"
0019    " && (probe_sigmaIetaIeta<0.03)"
0020    " && (probe_hadronicOverEm<0.15) "
0021    "))");
0022 
0023 
0024    TCut cleanGsf("(probe_gsfEle_HoverE<0.15) && (probe_gsfEle_HoverE<0.15)");
0025 
0026    TCut ID95("probe_passingId");
0027    TCut NotID95(cleanGsf && "!probe_passingId");
0028    TCut NotPPass("!probe_passing");
0029    TCut PassAll("probe_passingALL");
0030    TCut NotPassAll("!probe_passingALL");
0031    TCut ID80("probe_passingId80");
0032    TCut NotID80(cleanGsf && "!probe_passingId80");
0033    TCut EMINUS("probe_gsfEle_charge<0");
0034    TCut EMINUSSC("tag_gsfEle_charge>0");
0035    TCut EPLUS("probe_gsfEle_charge>0");
0036    TCut EPLUSSC("tag_gsfEle_charge<0");
0037    TCut BARREL("abs(probe_sc_eta)<1.4442");
0038    TCut BARRELSC("abs(probe_eta)<1.4442");
0039    TCut ENDCAPS("abs(probe_sc_eta)>1.566");
0040    TCut ENDCAPSSC("abs(probe_eta)>1.566");
0041 
0042 
0043 // //////////////////////////////////////////////////////////
0044    cout << "probe type" << "         efficiency " << "       Npass" <<  
0045       "       Nfail" << endl;
0046 // //////////////////////////////////////////////////////////
0047 
0048 
0049 
0050 // //////////////////////////////////////////////////////////
0051 //   //  Super cluster --> gsfElectron efficiency
0052 // //////////////////////////////////////////////////////////
0053 
0054    TCut GsfPass = preCut && cleanGsf;
0055    TCut GsfFail = cleanSC && preCut && NotPPass;
0056    ComputeEfficiency("Gsf", GsfPass, GsfFail);
0057    ComputeEfficiency("Gsf", GsfPass, GsfFail, "Barrel", BARREL, BARRELSC);
0058    ComputeEfficiency("Gsf", GsfPass, GsfFail, "Endcap", ENDCAPS, ENDCAPSSC);
0059    ComputeEfficiency("Gsf", GsfPass, GsfFail, "", "", "", "_eminus", EMINUS, EMINUSSC);
0060    ComputeEfficiency("Gsf", GsfPass, GsfFail, "", "", "", "_eplus", EPLUS, EPLUSSC);
0061    ComputeEfficiency("Gsf", GsfPass, GsfFail, "Barrel", BARREL, BARRELSC, "_eminus", EMINUS, EMINUSSC);
0062    ComputeEfficiency("Gsf", GsfPass, GsfFail, "Barrel",  BARREL, BARRELSC, "_eplus", EPLUS, EPLUSSC);
0063    ComputeEfficiency("Gsf", GsfPass, GsfFail, "Endcap", ENDCAPS, ENDCAPSSC, "_eminus", EMINUS, EMINUSSC);
0064    ComputeEfficiency("Gsf", GsfPass, GsfFail, "Endcap", ENDCAPS, ENDCAPSSC, "_eplus", EPLUS, EPLUSSC);
0065    cout << "########################################" << endl;
0066 
0067 
0068 
0069 // //////////////////////////////////////////////////////////
0070 //   //  gsfElectron --> WP-95 selection efficiency
0071 // //////////////////////////////////////////////////////////
0072 
0073 
0074    TCut Id95Pass = preCut && ID95;
0075    TCut Id95Fail = preCut && NotID95;
0076    ComputeEfficiency("Id95", Id95Pass, Id95Fail);
0077    ComputeEfficiency("Id95", Id95Pass, Id95Fail, "Barrel", BARREL, BARREL);
0078    ComputeEfficiency("Id95", Id95Pass, Id95Fail, "Endcap", ENDCAPS, ENDCAPS);
0079    ComputeEfficiency("Id95", Id95Pass, Id95Fail, "", "", "", "_eminus", EMINUS, EMINUS);
0080    ComputeEfficiency("Id95", Id95Pass, Id95Fail, "", "", "", "_eplus", EPLUS, EPLUS);
0081    ComputeEfficiency("Id95", Id95Pass, Id95Fail, "Barrel", BARREL, BARREL, "_eminus", EMINUS, EMINUS);
0082    ComputeEfficiency("Id95", Id95Pass, Id95Fail, "Barrel",  BARREL, BARREL, "_eplus", EPLUS, EPLUS);
0083    ComputeEfficiency("Id95", Id95Pass, Id95Fail, "Endcap", ENDCAPS, ENDCAPS, "_eminus", EMINUS, EMINUS);
0084    ComputeEfficiency("Id95", Id95Pass, Id95Fail, "Endcap", ENDCAPS, ENDCAPS, "_eplus", EPLUS, EPLUS);
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    cout << "########################################" << endl;
0106 
0107 
0108 
0109 // //////////////////////////////////////////////////////////
0110 //   //   WP-95 --> HLT triggering efficiency
0111 // //////////////////////////////////////////////////////////
0112 
0113 
0114    TCut HLT95Pass = preCut && ID95 && PassAll;
0115    TCut HLT95Fail = preCut && ID80 && NotPassAll;
0116    ComputeEfficiency("HLT95", HLT95Pass, HLT95Fail);
0117    ComputeEfficiency("HLT95", HLT95Pass, HLT95Fail, "Barrel", BARREL, BARREL);
0118    ComputeEfficiency("HLT95", HLT95Pass, HLT95Fail, "Endcap", ENDCAPS, ENDCAPS);
0119    ComputeEfficiency("HLT95", HLT95Pass, HLT95Fail, "", "", "", "_eminus", EMINUS, EMINUS);
0120    ComputeEfficiency("HLT95", HLT95Pass, HLT95Fail, "", "", "", "_eplus", EPLUS, EPLUS);
0121    ComputeEfficiency("HLT95", HLT95Pass, HLT95Fail, "Barrel", BARREL, BARREL, "_eminus", EMINUS, EMINUS);
0122    ComputeEfficiency("HLT95", HLT95Pass, HLT95Fail, "Barrel",  BARREL, BARREL, "_eplus", EPLUS, EPLUS);
0123    ComputeEfficiency("HLT95", HLT95Pass, HLT95Fail, "Endcap", ENDCAPS, ENDCAPS, "_eminus", EMINUS, EMINUS);
0124    ComputeEfficiency("HLT95", HLT95Pass, HLT95Fail, "Endcap", ENDCAPS, ENDCAPS, "_eplus", EPLUS, EPLUS);
0125    cout << "########################################" << endl;
0126 
0127 
0128 
0129 
0130 // //////////////////////////////////////////////////////////
0131 //   //   WP-80 --> HLT triggering efficiency
0132 // //////////////////////////////////////////////////////////
0133 
0134  
0135    TCut HLT80Pass = preCut && ID80 && PassAll;
0136    TCut HLT80Fail = preCut && ID80 && NotPassAll;
0137    ComputeEfficiency("HLT80", HLT80Pass, HLT80Fail);
0138    ComputeEfficiency("HLT80", HLT80Pass, HLT80Fail, "Barrel", BARREL, BARREL);
0139    ComputeEfficiency("HLT80", HLT80Pass, HLT80Fail, "Endcap", ENDCAPS, ENDCAPS);
0140    ComputeEfficiency("HLT80", HLT80Pass, HLT80Fail, "", "", "", "_eminus", EMINUS, EMINUS);
0141    ComputeEfficiency("HLT80", HLT80Pass, HLT80Fail, "", "", "", "_eplus", EPLUS, EPLUS);
0142    ComputeEfficiency("HLT80", HLT80Pass, HLT80Fail, "Barrel", BARREL, BARREL, "_eminus", EMINUS, EMINUS);
0143    ComputeEfficiency("HLT80", HLT80Pass, HLT80Fail, "Barrel",  BARREL, BARREL, "_eplus", EPLUS, EPLUS);
0144    ComputeEfficiency("HLT80", HLT80Pass, HLT80Fail, "Endcap", ENDCAPS, ENDCAPS, "_eminus", EMINUS, EMINUS);
0145    ComputeEfficiency("HLT80", HLT80Pass, HLT80Fail, "Endcap", ENDCAPS, ENDCAPS, "_eplus", EPLUS, EPLUS);
0146    cout << "########################################" << endl;
0147 }
0148 
0149 
0150 
0151 
0152 
0153 
0154 
0155 
0156 
0157 
0158 
0159 void ComputeEfficiency(char* primaryLabel, TCut primaryCutPass,
0160                        TCut primaryCutFail, char* secondaryLabel="", 
0161                        TCut secondaryCutPass="", TCut secondaryCutFail="",
0162                        char* tertiaryLabel="", 
0163                        TCut tertiaryCutPass = "", TCut tertiaryCutFail = "") 
0164 {
0165    // TFile* f = new TFile("newhistos.root");
0166    TFile* f = new TFile("Zee_new80.root");
0167 
0168   TTree* scTree = (TTree*) f->Get("PhotonToGsf/fitter_tree");
0169   TTree* gsfTree = (TTree*) f->Get("GsfToIso/fitter_tree");
0170   char namePass[50];
0171   char nameFail[50];
0172 
0173   sprintf(namePass,"Zmass%sPass%s%s",primaryLabel,secondaryLabel,tertiaryLabel);
0174   sprintf(nameFail,"Zmass%sFail%s%s",primaryLabel,secondaryLabel,tertiaryLabel);
0175   TH1F* histPass = createHistogram(namePass, 30);
0176   TH1F* histFail = createHistogram(nameFail, 12);
0177   gsfTree->Draw("mass>>"+TString(namePass), primaryCutPass && secondaryCutPass && 
0178   tertiaryCutPass,"goff");
0179   TString checkForSCTree(primaryLabel);
0180   TString plotVar = "mass>>"+TString(nameFail);
0181   if(checkForSCTree.Contains("Gsf"))
0182      scTree->Draw(plotVar,primaryCutFail && secondaryCutFail && tertiaryCutFail,"goff");
0183   else 
0184      gsfTree->Draw(plotVar,primaryCutFail && secondaryCutFail && tertiaryCutFail,"goff");
0185   ComputeEfficiency(*histPass, *histFail);
0186 
0187   delete histPass;
0188   delete histFail;
0189 }
0190 
0191 
0192 
0193 
0194 
0195 
0196 TH1F* createHistogram(char* name, int nbins=12) {
0197   TH1F* hist = new TH1F(name,name, nbins, 60, 120);
0198   hist->SetTitle("");
0199   char temp[100];
0200   sprintf(temp, "Events / %.1f GeV/c^{2}", 60./nbins);
0201   hist->GetXaxis()->SetTitle("m_{ee} (GeV/c^{2})");
0202   hist->GetYaxis()->SetTitle(temp);
0203   return hist;
0204 }
0205 
0206 
0207 
0208 
0209 
0210 
0211 void ComputeEfficiency( TH1& hist_pass, TH1& hist_fail)
0212 {
0213 
0214 
0215   TString effType = hist_pass.GetName();
0216   effType.ReplaceAll("Zmass","");
0217   effType.ReplaceAll("Pass","");
0218  
0219   double lowerLim = 0.0;
0220   double upperLim = 1.0;
0221   double numerator = hist_pass.Integral();
0222   double nfails    = hist_fail.Integral();
0223   double bkgFractionInNum = 0.0;
0224   double bkgFractionInFail = 0.0;
0225 
0226 //   if(effType.Contains("Gsf"))  { 
0227 //      bkgFractionInNum = 0.05;
0228 //      bkgFractionInFail = 0.02;
0229 //   }
0230 //   if(effType.Contains("Id"))  bkgFractionInFail= 0.05; 
0231 
0232   nfails  *= (1.0 - bkgFractionInFail); 
0233   numerator *= (1.0 - bkgFractionInNum); 
0234 
0235   double denominator = numerator + nfails;
0236 
0237   double eff = numerator/ denominator;
0238   //double err = sqrt(2*(1.0-eff)/totsig);
0239   ClopperPearsonLimits(numerator, denominator, lowerLim, upperLim);
0240   double errUp = upperLim -eff;
0241   double errLo = eff -  lowerLim;
0242 
0243 
0244   // ********** Make and save Canvas for the plots ********** //
0245   gROOT->ProcessLine(".L ~/tdrstyle.C");
0246   setTDRStyle();
0247   tdrStyle->SetErrorX(0.5);
0248   tdrStyle->SetPadLeftMargin(0.19);
0249   tdrStyle->SetPadRightMargin(0.10);
0250   tdrStyle->SetPadBottomMargin(0.15);
0251   tdrStyle->SetLegendBorderSize(0);
0252   tdrStyle->SetTitleYOffset(1.5);
0253   char temp[100];
0254 
0255   TString cname = TString("plot_") + effType;
0256   cname.ReplaceAll("Gsf","ScToGsf");
0257   cname.ReplaceAll("Id","GsfToWP");
0258   cname.Append("_Pass");
0259   TCanvas* c = new TCanvas(cname,cname,500,500);
0260   hist_pass.Draw("E1");
0261   TPaveText *plotlabel4 = new TPaveText(0.6,0.82,0.8,0.87,"NDC");
0262    plotlabel4->SetTextColor(kBlack);
0263    plotlabel4->SetFillColor(kWhite);
0264    plotlabel4->SetBorderSize(0);
0265    plotlabel4->SetTextAlign(12);
0266    plotlabel4->SetTextSize(0.03);
0267    sprintf(temp, "Signal = %.2f #pm %.2f", numerator, sqrt(numerator) );
0268    plotlabel4->AddText(temp);
0269   TPaveText *plotlabel5 = new TPaveText(0.6,0.77,0.8,0.82,"NDC");
0270    plotlabel5->SetTextColor(kBlack);
0271    plotlabel5->SetFillColor(kWhite);
0272    plotlabel5->SetBorderSize(0);
0273    plotlabel5->SetTextAlign(12);
0274    plotlabel5->SetTextSize(0.03);
0275    sprintf(temp, "Bkg fraction = %.2f", bkgFractionInNum);
0276    plotlabel5->AddText(temp);
0277   TPaveText *plotlabel6 = new TPaveText(0.6,0.87,0.8,0.92,"NDC");
0278    plotlabel6->SetTextColor(kBlack);
0279    plotlabel6->SetFillColor(kWhite);
0280    plotlabel6->SetBorderSize(0);
0281    plotlabel6->SetTextAlign(12);
0282    plotlabel6->SetTextSize(0.03);
0283    plotlabel6->AddText("Passing probes");
0284   TPaveText *plotlabel7 = new TPaveText(0.6,0.72,0.8,0.77,"NDC");
0285    plotlabel7->SetTextColor(kBlack);
0286    plotlabel7->SetFillColor(kWhite);
0287    plotlabel7->SetBorderSize(0);
0288    plotlabel7->SetTextAlign(12);
0289    plotlabel7->SetTextSize(0.03); 
0290    sprintf(temp, "Eff = %.3f + %.3f - %.3f", eff, errUp, errLo); 
0291    plotlabel7->AddText(temp);
0292   plotlabel4->Draw();
0293   plotlabel5->Draw();
0294   plotlabel6->Draw();
0295   plotlabel7->Draw();
0296   c->SaveAs( cname + TString(".eps"));
0297   c->SaveAs( cname + TString(".gif"));
0298   //c->SaveAs( cname + TString(".root"));
0299   delete c;
0300 
0301 
0302   cname.ReplaceAll("_Pass", "_Fail"); 
0303   TCanvas* c2 = new TCanvas(cname,cname,500,500);
0304   hist_fail.Draw("E1");
0305 
0306   TPaveText *plotlabel4 = new TPaveText(0.6,0.82,0.8,0.87,"NDC");
0307    plotlabel4->SetTextColor(kBlack);
0308    plotlabel4->SetFillColor(kWhite);
0309    plotlabel4->SetBorderSize(0);
0310    plotlabel4->SetTextAlign(12);
0311    plotlabel4->SetTextSize(0.03);
0312    sprintf(temp, "Signal = %.2f #pm %.2f", nfails, sqrt(nfails) );
0313    plotlabel4->AddText(temp);
0314   TPaveText *plotlabel5 = new TPaveText(0.6,0.77,0.8,0.82,"NDC");
0315    plotlabel5->SetTextColor(kBlack);
0316    plotlabel5->SetFillColor(kWhite);
0317    plotlabel5->SetBorderSize(0);
0318    plotlabel5->SetTextAlign(12);
0319    plotlabel5->SetTextSize(0.03);
0320    sprintf(temp, "Bkg fraction = %.2f", bkgFractionInFail);
0321    plotlabel5->AddText(temp);
0322   TPaveText *plotlabel6 = new TPaveText(0.6,0.87,0.8,0.92,"NDC");
0323    plotlabel6->SetTextColor(kBlack);
0324    plotlabel6->SetFillColor(kWhite);
0325    plotlabel6->SetBorderSize(0);
0326    plotlabel6->SetTextAlign(12);
0327    plotlabel6->SetTextSize(0.03);
0328    plotlabel6->AddText("Failing probes");
0329   TPaveText *plotlabel7 = new TPaveText(0.6,0.72,0.8,0.77,"NDC");
0330    plotlabel7->SetTextColor(kBlack);
0331    plotlabel7->SetFillColor(kWhite);
0332    plotlabel7->SetBorderSize(0);
0333    plotlabel7->SetTextAlign(12);
0334    plotlabel7->SetTextSize(0.03);
0335    sprintf(temp, "Eff = %.3f + %.3f - %.3f", eff, errUp, errLo);
0336    plotlabel7->AddText(temp);
0337    plotlabel4->Draw();
0338    plotlabel5->Draw();
0339    plotlabel6->Draw();
0340    plotlabel7->Draw();
0341 
0342      c2->SaveAs( cname + TString(".eps"));
0343     c2->SaveAs( cname + TString(".gif"));
0344   //c2->SaveAs( cname + TString(".root"));
0345   delete c2;
0346 
0347   // cout << "########################################" << endl;
0348   effType.ReplaceAll("Gsf","");
0349   effType.ReplaceAll("Id95_","");
0350   effType.ReplaceAll("Id95","");
0351   effType.ReplaceAll("Id80_","");
0352   effType.ReplaceAll("Id80","");
0353   effType.ReplaceAll("HLT95_","");
0354   effType.ReplaceAll("HLT95","");
0355   effType.ReplaceAll("HLT80_","");
0356   effType.ReplaceAll("HLT80","");
0357 
0358  char* effTypeToPrint = (char*) effType;
0359 
0360   cout << effTypeToPrint << "    " 
0361        << setiosflags(ios::fixed) << setprecision(4) << eff 
0362        << " + " << errUp << " - " <<  errLo 
0363        << setiosflags(ios::fixed) << setprecision(0)
0364        << "    " << numerator << "    "  << nfails << endl;
0365   //cout << "########################################" << endl;
0366 
0367 }
0368 
0369 
0370 double ErrorInProduct(double x, double errx, double y, 
0371                       double erry, double corr) {
0372    double xFrErr = errx/x;
0373    double yFrErr = erry/y;
0374    return sqrt(xFrErr**2 +yFrErr**2 + 2.0*corr*xFrErr*yFrErr)*x*y;
0375 }
0376 
0377 
0378 void ClopperPearsonLimits(double numerator, double denominator, 
0379 double &lowerLimit, double &upperLimit, const double CL_low=1.0, 
0380 const double CL_high=1.0) 
0381 {  
0382 //Confidence intervals are in the units of \sigma.
0383 
0384    double ratio = numerator/denominator;
0385    
0386 // first get the lower limit
0387    if(numerator==0)   lowerLimit = 0.0; 
0388    else { 
0389       double v=ratio/2; 
0390       double vsL=0; 
0391       double vsH=ratio; 
0392       double p=CL_low/100;
0393       while((vsH-vsL)>1e-5) { 
0394          if(BinP(denominator,v,numerator,denominator)>p) 
0395          { vsH=v; v=(vsL+v)/2; } 
0396          else { vsL=v; v=(v+vsH)/2; } 
0397       }
0398       lowerLimit = v; 
0399    }
0400    
0401 // now get the upper limit
0402    if(numerator==denominator) upperLimit = 1.0;
0403    else { 
0404       double v=(1+ratio)/2; 
0405       double vsL=ratio; 
0406       double vsH=1; 
0407       double p=CL_high/100;
0408       while((vsH-vsL)>1e-5) { 
0409          if(BinP(denominator,v,0,numerator)<p) { vsH=v; v=(vsL+v)/2; } 
0410          else { vsL=v; v=(v+vsH)/2; } 
0411       }
0412       upperLimit = v;
0413    }
0414 }
0415 
0416 
0417 
0418 
0419 double BinP(int N, double p, int x1, int x2) {
0420    double q=p/(1-p); 
0421    int k=0; 
0422    double v = 1; 
0423    double s=0; 
0424    double tot=0.0;
0425     while(k<=N) {
0426        tot=tot+v;
0427        if(k>=x1 & k<=x2) { s=s+v; }
0428        if(tot>1e30){s=s/1e30; tot=tot/1e30; v=v/1e30;}
0429        k=k+1; 
0430        v=v*q*(N+1-k)/k;
0431     }
0432     return s/tot;
0433 }
0434 
0435 
0436 
0437
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!