ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/AnalysisFramework/Plotting/Modules/Systematics.C
(Generate patch)

Comparing UserCode/cbrown/AnalysisFramework/Plotting/Modules/Systematics.C (file contents):
Revision 1.29 by buchmann, Wed Sep 7 06:42:21 2011 UTC vs.
Revision 1.42 by buchmann, Mon Oct 24 15:05:37 2011 UTC

# Line 1 | Line 1
1   #include <iostream>
2   #include <vector>
3   #include <sys/stat.h>
4 + #include <algorithm>
5 + #include <cmath>
6  
7   #include <TMath.h>
8   #include <TColor.h>
# Line 100 | Line 102 | float allcontributionsplot(TTree* events
102          hname=GetNumericHistoName();
103          TH1F* hosofn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
104          events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kMassCut&&JZBNegCut&&cutOSOF,"goff");
103
104        hname=GetNumericHistoName();
105        TH1F* sbhossfp = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
106        events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kSidebandCut&&JZBPosCut&&cutOSSF,"goff");
107        hname=GetNumericHistoName();
108        TH1F* sbhossfn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
109        events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kSidebandCut&&JZBNegCut&&cutOSSF,"goff");
110
111        hname=GetNumericHistoName();
112        TH1F* sbhosofp = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
113        events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kSidebandCut&&JZBPosCut&&cutOSOF,"goff");
114        hname=GetNumericHistoName();
115        TH1F* sbhosofn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
116        events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kSidebandCut&&JZBNegCut&&cutOSOF,"goff");
105          
106 <        float obs = hossfp->Integral();
107 <        float pred= hossfn->Integral() + (1.0/3)*( hosofp->Integral() - hosofn->Integral() + sbhossfp->Integral() - sbhossfn->Integral() + sbhosofp->Integral() - sbhosofn->Integral());
106 >        float obs=0;
107 >        float pred=0;
108 >        flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
109 >        if(PlottingSetup::RestrictToMassPeak) {
110 >          hname=GetNumericHistoName();
111 >          TH1F* sbhossfp = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
112 >          events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kSidebandCut&&JZBPosCut&&cutOSSF,"goff");
113 >          hname=GetNumericHistoName();
114 >          TH1F* sbhossfn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
115 >          events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kSidebandCut&&JZBNegCut&&cutOSSF,"goff");
116 >          
117 >          hname=GetNumericHistoName();
118 >          TH1F* sbhosofp = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
119 >          events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kSidebandCut&&JZBPosCut&&cutOSOF,"goff");
120 >          hname=GetNumericHistoName();
121 >          TH1F* sbhosofn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
122 >          events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kSidebandCut&&JZBNegCut&&cutOSOF,"goff");
123 >          
124 >          obs = hossfp->Integral();
125 >          pred= hossfn->Integral() + (1.0/3)*( hosofp->Integral() - hosofn->Integral() + sbhossfp->Integral() - sbhossfn->Integral() + sbhosofp->Integral() - sbhosofn->Integral());
126 >          delete sbhossfp,sbhossfn,sbhosofp,sbhosofn;
127 >        } else {
128 >          obs = hossfp->Integral();
129 >          pred= hossfn->Integral() + (hosofp->Integral() - hosofn->Integral());
130 >        }
131          
132          delete hossfp,hossfn,hosofp,hosofn;
122        delete sbhossfp,sbhossfn,sbhosofp,sbhosofn;
133          return obs-pred;
134   }
135  
# Line 188 | Line 198 | void master_formula(std::vector<float> e
198  
199   //________________________________________________________________________________________
200   // Get normalization factor for the PDFs
201 < float get_norm_pdf_factor(TTree *events, int k) {
201 > float get_norm_pdf_factor(TTree *events, int k, string addcut) {
202  
203    TH1F *haux = new TH1F("haux", "", 10000, 0, 5);
204    char nameVar[20];
205    sprintf(nameVar, "pdfW[%d]", k);
206 <  events->Project("haux", nameVar);
206 >  events->Project("haux", nameVar, addcut.c_str());
207    float thisW = haux->Integral();
208    events->Project("haux", "pdfW[0]");
209    float normW = haux->Integral();
# Line 215 | Line 225 | float pileup(TTree *events, bool require
225          jzbMax = myJzbMax;
226          
227          // Acceptance cuts
228 <        TCut kbase("abs(genMll-91.2)<20&&genNjets>2&&genZPt>0&&abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
228 >        flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
229 >        TCut kbase(PlottingSetup::genMassCut&&"genNjets>2&&genZPt>0"&&cutmass&&cutOSSF);
230          if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
231          
232 <        if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
232 >        if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
233          TH1F* hLM4 = plotEff(events,kbase,informalname);
234          hLM4->SetMinimum(0.);
235          
# Line 241 | Line 252 | float pileup(TTree *events, bool require
252   //____________________________________________________________________________________
253   // Effect of peak shifting
254   void PeakError(TTree *events,float &result, string mcjzb, float peakerr,string addcut="") {
255 +    //Note: the cut used here is something like (JZBEXPRESSION+(peakerr)>50) without all the other cuts, to increase statistics (particularly for scans)
256          TString peakup("("+TString(mcjzb)+"+"+TString(any2string(TMath::Abs(peakerr)))+")"+geq_or_leq()+TString(any2string(jzbSel)));
257          TString peakdown("("+TString(mcjzb)+"-"+TString(any2string(TMath::Abs(peakerr)))+")"+geq_or_leq()+TString(any2string(jzbSel)));
258          TString peakcentral("("+TString(mcjzb)+")"+geq_or_leq()+TString(any2string(jzbSel)));
259          TString npeakup("("+TString(mcjzb)+"+"+TString(any2string(TMath::Abs(peakerr)))+")"+ngeq_or_leq()+"-"+TString(any2string(jzbSel)));
260          TString npeakdown("("+TString(mcjzb)+"-"+TString(any2string(TMath::Abs(peakerr)))+")"+ngeq_or_leq()+"-"+TString(any2string(jzbSel)));
261          TString npeakcentral("("+TString(mcjzb)+")"+ngeq_or_leq()+"-"+TString(any2string(jzbSel)));
250        
262          nBins = 1;
263          string informalname="PeakErrorCalculation";
264          float resup,resdown,rescent;
# Line 270 | Line 281 | void PeakError(TTree *events,float &resu
281            else if(i==1) resdown=res;
282            else if(i==2) resup=res;
283          }
284 <        if(TMath::Abs(rescent-resup)>TMath::Abs(rescent-resdown)) result=(TMath::Abs(rescent-resup)/rescent);
285 <        else result=(TMath::Abs(rescent-resdown)/rescent);
284 >        if(TMath::Abs(rescent-resup)>TMath::Abs(rescent-resdown)) result=(TMath::Abs(rescent-resup)/(float)rescent);
285 >        else result=(TMath::Abs(rescent-resdown)/(float)rescent);
286   }
287  
288   //____________________________________________________________________________________
289   // Total selection efficiency (MC)
290 < void MCefficiency(TTree *events,float &result, float &resulterr,string mcjzb,bool requireZ,int Neventsinfile, string addcut="", int k = 0) {
290 > //returns the efficiency WITHOUT signal contamination, and the result and resulterr contain the result and the corresponding error
291 > Value MCefficiency(TTree *events,float &result, float &resulterr,string mcjzb,bool requireZ,int Neventsinfile, string addcut="", int k = 0) {
292 >        write_warning(__FUNCTION__,"Setting automatized to off!"); automatized=false;
293 >        if(!events) {
294 >          write_error(__FUNCTION__,"Tree passed for efficiency calculation is invalid!");
295 >          result=0;
296 >          resulterr=0;
297 >          return Value(0,0);
298 >        }
299          
300          char jzbSelStr[256]; sprintf(jzbSelStr,"%f",jzbSel);
301          // All acceptance cuts at gen. level
302          //TCut kbase("abs(genMll-91.2)<20&&genNjets>2&&genZPt>0&&genJZB"+geq_or_leq()+TString(jzbSelStr)+"&&genId1==-genId2");
303          TCut kbase("");
304 <        if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
304 >        
305 >        flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
306 >        if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
307          if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
308          // Corresponding reco. cuts
309 <        TCut ksel("pfJetGoodNum>2&&abs(mll-91.2)<20&&id1==id2&&"+TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
310 <        TCut ksel2("pfJetGoodNum>2&&abs(mll-91.2)<20&&id1==id2&&"+TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr));
309 >        
310 >        flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
311 >        TCut ksel;//("pfJetGoodNum>2"&&cutmass&&"id1==id2&&"+TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
312 >        TCut ksel2;//("pfJetGoodNum>2"&&cutmass&&"id1==id2&&"+TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr));
313 >        flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
314 >        if(PlottingSetup::RestrictToMassPeak||!ConsiderSignalContaminationForLimits) {
315 >          ksel=TCut("pfJetGoodNum>2"&&cutmass&&"id1==id2&&"+TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
316 >          ksel2=TCut("pfJetGoodNum>2"&&cutmass&&"id1==id2&&"+TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr));
317 >        } else {
318 >          //for off peak analysis we don't use the OSSF condition here yet so we can recycle these two cuts for the em condition!
319 >          ksel=TCut("pfJetGoodNum>2"&&cutmass&&(TString(mcjzb)+geq_or_leq()+TString(jzbSelStr)));
320 >          ksel2=TCut("pfJetGoodNum>2"&&cutmass&&(TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr)));
321 >        }
322 >            
323          TCut posSide = kbase&&ksel;
324          TCut negSide = kbase&&ksel2;
325          string sposSide(posSide);
326          string snegSide(negSide);
327          char var[20];
328          sprintf(var, "pdfW[%d]", k);
329 +        if(k==-1) sprintf(var,"1.0");//case in which we don't want to evaluate PDFs
330          string svar(var);
331 <        string newPosSide = "(" + sposSide + ")*" + svar;
332 <        string newNegSide = "(" + snegSide + ")*" + svar;
331 >        string newPosSide = "((id1==id2)&&(" + sposSide + "))*" + svar;
332 >        string newNegSide = "((id1==id2)&&(" + snegSide + "))*" + svar;
333 >        string emnewPosSide = "((id1!=id2)&&(" + sposSide + "))*" + svar; // only used for off peak analysis
334 >        string emnewNegSide = "((id1!=id2)&&(" + snegSide + "))*" + svar; // only used for off peak analysis
335  
336          TH1F *effh= new TH1F("effh","effh",1,-14000,14000);
337          if(k>=0)events->Draw((mcjzbexpression+">>effh").c_str(), newPosSide.c_str(),"goff");
338 <        else events->Draw((mcjzbexpression+">>effh").c_str(), sposSide.c_str(),"goff");
338 >        else events->Draw((mcjzbexpression+">>effh").c_str(), (sposSide+"&&(id1==id2)").c_str(),"goff");//the OSSF condition is added for the offpeak analysis, in onpeak case it's there already but doesn't change anything.
339          Float_t sel = effh->Integral();
340          Float_t nsel=0;
341 +        
342 +        ///----------------------------------------------- THIS PART REQUIRES STUDYING! -------------------------
343 +        
344          if(ConsiderSignalContaminationForLimits) {
345 <          if(k>=0)events->Draw((mcjzbexpression+">>effh").c_str(), newNegSide.c_str(),"goff");
346 <          else events->Draw((mcjzbexpression+">>effh").c_str(), snegSide.c_str(),"goff");
347 <          nsel = effh->Integral();
345 >          flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
346 >          if(PlottingSetup::RestrictToMassPeak) {
347 >            events->Draw((mcjzbexpression+">>effh").c_str(), newNegSide.c_str(),"goff");
348 >            nsel += effh->Integral();
349 >          } else {
350 >            events->Draw((mcjzbexpression+">>effh").c_str(), newNegSide.c_str(),"goff");
351 >            nsel += effh->Integral();
352 >            events->Draw((mcjzbexpression+">>effh").c_str(), emnewPosSide.c_str(),"goff");
353 >            nsel += effh->Integral();
354 >            events->Draw((mcjzbexpression+">>effh").c_str(), emnewNegSide.c_str(),"goff");
355 >            nsel -= effh->Integral();
356 >          }
357          }
358 +
359          //Corrections due to normalization in the PDF. This has to be applied as well to the number of events in a file if the definition changes at some point.
360          float normFactor = 1;
361 <        if(k>=0) get_norm_pdf_factor(events, k);
361 >        if(k>=0) get_norm_pdf_factor(events, k, addcut);
362          sel = sel/normFactor;
363          nsel = nsel/normFactor;
364  
# Line 317 | Line 366 | void MCefficiency(TTree *events,float &r
366   //      Float_t tot = events->GetSelectedRows();
367          Float_t tot = Neventsinfile;
368          
369 +        Value result_wo_signalcont;
370 +
371          if(ConsiderSignalContaminationForLimits) {
372            result=(sel-nsel)/tot;
373            resulterr=(1.0/tot)*TMath::Sqrt(sel+nsel+(sel-nsel)*(sel-nsel)/tot);
374 +          result_wo_signalcont=Value(sel/tot,TMath::Sqrt(sel/tot*(1+sel/tot)/tot));
375          } else {//no signal contamination considered:
376            result=(sel)/tot;
377            resulterr=TMath::Sqrt(sel/tot*(1+sel/tot)/tot);
378 +          result_wo_signalcont=Value(result,resulterr);
379          }
380 <        if(!automatized) dout << "  MC efficiency: " << result << "+-" << resulterr << "  ( JZB>" << jzbSel << " : " << sel << " , JZB<-" << jzbSel << " : " << nsel << " and nevents=" << tot << ") with normFact=" << normFactor << std::endl;
380 >        if(!automatized && k>0 ) dout << "PDF assessment: ";
381 >        if(!automatized) dout << "  MC efficiency: " << result << "+-" << resulterr << "  ( JZB>" << jzbSel << " : " << sel << " , signal contamination : " << nsel << " and nevents=" << tot << ") with normFact=" << normFactor << std::endl;
382          delete effh;
383 +        return result_wo_signalcont;
384   }
385  
386  
387 +
388   //____________________________________________________________________________________
389 < // Total selection efficiency (MC)
389 > // Selection efficiency for one process (MC)
390   vector<float> processMCefficiency(TTree *events,string mcjzb,bool requireZ,int Neventsinfile, string addcut) {
391    vector<float> process_efficiencies;
392    for(int iprocess=0;iprocess<=10;iprocess++) {
# Line 345 | Line 401 | vector<float> processMCefficiency(TTree
401          
402  
403   void JZBefficiency(TTree *events, string informalname, float &jzbeff, float &jzbefferr, bool requireZ, string addcut="") {
404 <        TCut kbase("abs(genMll-91.2)<20&&genNjets>2&&genZPt>0&&abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
404 >        TCut kbase(genMassCut&&"genNjets>2&&genZPt>0"&&cutmass&&cutOSSF);
405          if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
406 <        if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
406 >        if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
407          TH1F* hLM4 = plotEff(events,kbase,informalname);
408          Int_t bin = hLM4->FindBin(jzbSel); // To get the error
409          jzbeff=Interpolate(jzbSel,hLM4);
# Line 359 | Line 415 | void JZBefficiency(TTree *events, string
415   //________________________________________________________________________
416   // Effect of energy scale on efficiency
417   void JZBjetScale(TTree *events, float &jesdown, float &jesup, string informalname,bool requireZ,string addcut="",float syst=0.1, Float_t jzbSelection=-1, TString plotName = "" ) {
418 <        TCut kbase("abs(genMll-91.2)<20&&genZPt>0");
418 >        TCut kbase(genMassCut&&"genZPt>0");
419          if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
420 <        if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
420 >        flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
421 >        if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
422  
423 <        TCut ksel("abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
423 >        TCut ksel(cutmass&&cutOSSF);
424          TCut nJets("pfJetGoodNum>2");
425          stringstream down,up;
426          down << "pfJetGoodNum"<<30*(1-syst)<<">=3";
# Line 397 | Line 454 | void JZBjetScale(TTree *events, float &j
454   // Effect of energy scale on JZB efficiency
455   void doJZBscale(TTree *events, float &down, float &up, float &syst, float systematic, string informalname, bool requireZ, string addcut) {
456          
457 <        TCut kbase("abs(genMll-91.2)<20&&genZPt>0&&genNjets>2");
457 >        TCut kbase(genMassCut&&"genZPt>0&&genNjets>2");
458          if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
459 <        if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
460 <        TCut ksel("abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
459 >        flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
460 >        if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
461 >        TCut ksel(cutmass&&cutOSSF);
462          
463          nBins =    50;
464          jzbMin =   0.5*jzbSel;
# Line 412 | Line 470 | void doJZBscale(TTree *events, float &do
470          Float_t eff  = Interpolate(jzbSel,hist);
471          Float_t effp = Interpolate(jzbSel*(1.+systematic),hist);
472          Float_t effm = Interpolate(jzbSel*(1.-systematic),hist);
473 <        if(!automatized) dout << "  efficiency at JZB==" << jzbSel*(1.+systematic)  << "(-"<<syst*100<<"%) : " << effp << " (" << ((effp-eff)/eff)*100. << "%)"  << std::endl;
473 >        if(!automatized) dout << "  efficiency at JZB==" << jzbSel*(1.+systematic)  << "(-"<<systematic*100<<"%) : " << effp << " (" << ((effp-eff)/eff)*100. << "%)"  << std::endl;
474          if(!automatized) dout << "  efficiency at JZB==" << jzbSel  << ": " << eff << std::endl;
475 <        if(!automatized) dout << "  efficiency at JZB==" << jzbSel*(1.-systematic)  << "(-"<<syst*100<<"%) : " << effm << " (" << ((effm-eff)/eff)*100. << "%)"  << std::endl;
475 >        if(!automatized) dout << "  efficiency at JZB==" << jzbSel*(1.-systematic)  << "(-"<<systematic*100<<"%) : " << effm << " (" << ((effm-eff)/eff)*100. << "%)"  << std::endl;
476          up=((effp-eff)/eff);
477          down=((effm-eff)/eff);
478   }
# Line 424 | Line 482 | void doJZBscale(TTree *events, float &do
482   void JZBresponse(TTree *events, bool requireZ, float &resp, float &resperr, string addcut="",bool isMET = kFALSE, Float_t myJzbMax = 200., Int_t nPeriods = 9 ) {
483          
484          jzbMin = 20;
485 <        TCut kbase("abs(genMll-91.2)<20&&genZPt>0&&genNjets>2");
485 >        flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
486 >        TCut kbase(genMassCut&&"genZPt>0&&genNjets>2");
487          if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
488 <        if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
489 <        TCut ksel("abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
488 >        flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
489 >        if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
490 >        flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
491 >        TCut ksel(cutmass&&cutOSSF);
492          
493          TProfile* hJzbResp = new TProfile("hJzbResp","JZB response  ; JZB true (GeV/c); JZB reco. / JZB true", nPeriods, jzbMin, myJzbMax, "" );
494          
# Line 438 | Line 499 | void JZBresponse(TTree *events, bool req
499          hJzbResp->SetMinimum(0.2);
500          hJzbResp->Fit("pol0","Q");
501          TF1 *fittedfunction = hJzbResp->GetFunction("pol0");
502 <        resp=fittedfunction->GetParameter(0);
503 <        resperr=fittedfunction->GetParError(0);
504 <        if(!automatized) dout << "  Response: " << resp << " +/- " << resperr << endl;
502 >        if(!fittedfunction) {
503 >                // in case there are not enough points passing our selection
504 >                cout << "OOPS response function invalid, assuming 100% error !!!!" << endl;
505 >                resp=1;
506 >                resperr=1;
507 >        } else {
508 >                resp=fittedfunction->GetParameter(0);
509 >                resperr=fittedfunction->GetParError(0);
510 >                if(!automatized) dout << "  Response: " << resp << " +/- " << resperr << endl;
511 >        }
512          delete hJzbResp;
513   }
514  
# Line 478 | Line 546 | void do_systematics_for_one_file(TTree *
546    float triggereff=5.0/100;// in range [0,1]
547    dout << "Trigger efficiency not implemented in this script  yet, still using external one" << endl;
548    float leptonseleff=2.0/100;// in range [0,1]
549 +  leptonseleff=TMath::Sqrt(leptonseleff*leptonseleff+leptonseleff*leptonseleff); // because the 2% is per lepton
550    dout << "Lepton selection efficiency not implemented in this script  yet, still using external one" << endl;
551    
552    int NPdfs=0;
# Line 485 | Line 554 | void do_systematics_for_one_file(TTree *
554    
555    float mceff,mcefferr,jzbeff,jzbefferr;
556    if(!automatized) dout << "MC efficiencies:" << endl;
557 <  MCefficiency(events,mceff,mcefferr,mcjzb,requireZ,Neventsinfile,addcut,-1);
558 <  JZBefficiency(events,informalname,jzbeff,jzbefferr,requireZ,addcut);
557 >  Value mceff_nosigcont = MCefficiency(events,mceff,mcefferr,mcjzb,requireZ,Neventsinfile,addcut,-1);
558 >  if(!automatized) cout << "   Without signal contamination, we find an efficiency of " << mceff_nosigcont << endl;
559 >
560 >  if(PlottingSetup::computeJZBefficiency) JZBefficiency(events,informalname,jzbeff,jzbefferr,requireZ,addcut);
561    if(!automatized) dout << "JZB efficiency: " << jzbeff << "+/-" << jzbefferr << endl;
562    
563    if(!automatized) dout << "Error from Peak position:" << endl;
# Line 503 | Line 574 | void do_systematics_for_one_file(TTree *
574    
575    if(!automatized) dout << "JZB response: " << std::endl;
576    float resp,resperr;
577 <  JZBresponse(events,requireZ,resp,resperr,addcut);
577 >  if(PlottingSetup::computeJZBresponse) {
578 >        if(!automatized) dout << "JZB response: " << std::endl;
579 >        JZBresponse(events,requireZ,resp,resperr,addcut);
580 >  }
581  
582    if(!automatized) dout << "Pileup: " << std::endl;
583 <  float resolution=pileup(events,requireZ,informalname,addcut);
583 >  float resolution;
584 >  resolution=pileup(events,requireZ,informalname,addcut);
585  
586    float PDFuncert=0;
587 +  if(!automatized) dout << "Assessing PDF uncertainty: " << std::endl;
588    if(ismSUGRA) PDFuncert = get_pdf_uncertainty(events, mcjzb, requireZ, Neventsinfile, NPdfs, addcut);
589  
590    dout << "_______________________________________________" << endl;
# Line 523 | Line 599 | void do_systematics_for_one_file(TTree *
599    dout << "Resolution : " << resolution << endl; // in range [0,1]
600    dout << "From peak : " << sysfrompeak << endl; // in range [0,1]
601    if(ismSUGRA) dout << "PDF uncertainty  : " << PDFuncert << endl; // in range [0,1]
602 <  dout << "JZB efficiency: " << jzbeff << "+/-" << jzbefferr << " (not yet included below) " << endl; // in range [0,1]
603 <  dout << "JZB response  : " << resp << " +/-" << resperr << " (not yet included below) " << endl; // in range [0,1]
602 >  if(PlottingSetup::computeJZBefficiency) dout << "JZB efficiency: " << jzbeff << "+/-" << jzbefferr << " (not yet included below) " << endl; // in range [0,1]
603 >  if(PlottingSetup::computeJZBresponse)dout << "JZB response  : " << resp << " +/-" << resperr << " (not yet included below) " << endl; // in range [0,1]
604    
605    float toterr=0;
606    toterr+=(triggereff)*(triggereff);
# Line 540 | Line 616 | void do_systematics_for_one_file(TTree *
616    
617    dout << "FINAL RESULT : " << 100*mceff << " +/- "<< 100*mcefferr << " (stat) +/- " << 100*systerr << " (syst)   %" << endl;
618    dout << "     we thus use the sqrt of the sum of the squares of the stat & syst err, which is : " << 100*toterr << endl;
619 +  dout << "_______________________________________________" << endl;
620    
621    //Do not modify the lines below or mess with the order; this order is expected by all limit calculating functions!
622    vector<float> res;
# Line 551 | Line 628 | void do_systematics_for_one_file(TTree *
628    if(fabs(jesup)>fabs(jesdown)) res.push_back(fabs(jesup)); else res.push_back(fabs(jesdown));
629    if(fabs(scaleup)>fabs(scaledown)) res.push_back(fabs(scaleup)); else res.push_back(fabs(scaledown));
630    res.push_back(fabs(resolution));
631 +  res.push_back(mceff_nosigcont.getValue());
632 +  res.push_back(mceff_nosigcont.getError());
633    if(ismSUGRA) res.push_back(PDFuncert);
634    results.push_back(res);
635   }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines