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.8 by buchmann, Thu Jul 21 08:26:12 2011 UTC vs.
Revision 1.18 by buchmann, Wed Aug 17 17:07:18 2011 UTC

# Line 158 | Line 158 | TH1F* plotEff(TTree* events, TCut kbase,
158  
159   //________________________________________________________________________________________
160   // Pile-up efficiency
161 < float pileup(TTree *events, bool requireZ, string informalname, Float_t myJzbMax = 140. ) {
161 > float pileup(TTree *events, bool requireZ, string informalname, string addcut="",Float_t myJzbMax = 140. ) {
162          nBins = 16;
163          jzbMax = myJzbMax;
164          
165          // Acceptance cuts
166          TCut kbase("abs(genMll-91.2)<20&&genNjets>2&&genZPt>0&&abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
167 +        if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
168 +        
169          if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
170          TH1F* hLM4 = plotEff(events,kbase,informalname);
171          hLM4->SetMinimum(0.);
# Line 180 | Line 182 | float pileup(TTree *events, bool require
182          if(!automatized) dout << "  PU: " << funcUp->Eval(jzbSel) << " " <<  func->Eval(jzbSel)
183          << "(" << (funcUp->Eval(jzbSel)-func->Eval(jzbSel))/func->Eval(jzbSel)*100. << "%)" << std::endl;
184          
185 <        return (funcUp->Eval(jzbSel)-func->Eval(jzbSel))/func->Eval(jzbSel)*100.;
185 >        return (funcUp->Eval(jzbSel)-func->Eval(jzbSel))/func->Eval(jzbSel);
186          
187   }
188  
189   //____________________________________________________________________________________
190   // Effect of peak shifting
191 < void PeakError(TTree *events,float &result, string mcjzb, float peakerr) {
191 > void PeakError(TTree *events,float &result, string mcjzb, float peakerr,string addcut="") {
192          TString peakup("("+TString(mcjzb)+"+"+TString(any2string(TMath::Abs(peakerr)))+")"+geq_or_leq()+TString(any2string(jzbSel)));
193          TString peakdown("("+TString(mcjzb)+"-"+TString(any2string(TMath::Abs(peakerr)))+")"+geq_or_leq()+TString(any2string(jzbSel)));
194          TString peakcentral("("+TString(mcjzb)+")"+geq_or_leq()+TString(any2string(jzbSel)));
# Line 209 | Line 211 | void PeakError(TTree *events,float &resu
211              poscut=peakup;
212              negcut=npeakup;
213            }
214 <          float res=allcontributionsplot(events,cutnJets,cutmass,sidebandcut,poscut.c_str(),negcut.c_str());
214 >          float res;
215 >          if(addcut=="") res=allcontributionsplot(events,cutnJets,cutmass,sidebandcut,poscut.c_str(),negcut.c_str());
216 >          else res=allcontributionsplot(events,cutnJets&&addcut.c_str(),cutmass,sidebandcut,poscut.c_str(),negcut.c_str());
217            if(i==0) rescent=res;
218            else if(i==1) resdown=res;
219            else if(i==2) resup=res;
220          }
221 <        if(TMath::Abs(rescent-resup)>TMath::Abs(rescent-resdown)) result=(TMath::Abs(rescent-resup)/rescent)*100;
222 <        else result=(TMath::Abs(rescent-resdown)/rescent)*100;
221 >        if(TMath::Abs(rescent-resup)>TMath::Abs(rescent-resdown)) result=(TMath::Abs(rescent-resup)/rescent);
222 >        else result=(TMath::Abs(rescent-resdown)/rescent);
223   }
224  
225   //____________________________________________________________________________________
226   // Total selection efficiency (MC)
227 < void MCefficiency(TTree *events,float &result, float &resulterr,string mcjzb,bool requireZ) {
227 > void MCefficiency(TTree *events,float &result, float &resulterr,string mcjzb,bool requireZ,int Neventsinfile, string addcut="") {
228          
229          char jzbSelStr[256]; sprintf(jzbSelStr,"%f",jzbSel);
230          // All acceptance cuts at gen. level
231 <        TCut kbase("abs(genMll-91.2)<20&&genNjets>2&&genZPt>0&&genJZB"+geq_or_leq()+TString(jzbSelStr)+"&&genId1==-genId2");
231 >        //TCut kbase("abs(genMll-91.2)<20&&genNjets>2&&genZPt>0&&genJZB"+geq_or_leq()+TString(jzbSelStr)+"&&genId1==-genId2");
232 >        TCut kbase("");
233          if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
234 +        if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
235          // Corresponding reco. cuts
236 <        TCut ksel("abs(mll-91.2)<20&&id1==id2&&"+TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
237 <        
236 >        TCut ksel("pfJetGoodNum>2&&abs(mll-91.2)<20&&id1==id2&&"+TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
237 >        TCut ksel2("pfJetGoodNum>2&&abs(mll-91.2)<20&&id1==id2&&"+TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr));
238          events->Draw(mcjzbexpression.c_str(),kbase&&ksel,"goff");
239          Float_t sel = events->GetSelectedRows();
240 <        events->Draw(mcjzbexpression.c_str(),kbase,"goff");
241 <        Float_t tot = events->GetSelectedRows();
240 >        events->Draw(mcjzbexpression.c_str(),kbase&&ksel2,"goff");
241 >        Float_t nsel = events->GetSelectedRows();
242 > //      events->Draw(mcjzbexpression.c_str(),kbase,"goff");
243 > //      Float_t tot = events->GetSelectedRows();
244 >        Float_t tot = Neventsinfile;
245          
246 <        result=sel/tot;
246 >        result=(sel-nsel)/tot;
247          resulterr=TMath::Sqrt(sel/tot*(1-sel/tot)/tot);
248 <        if(!automatized) dout << "  MC efficiency: " << result << "+-" << resulterr << std::endl;
248 >        dout << "  MC efficiency: " << result << "+-" << resulterr << "  ( JZB>" << jzbSel << " : " << sel << " , JZB<-" << jzbSel << " : " << nsel << " and nevents=" << tot << ")" << std::endl;
249   }
250  
251 < void JZBefficiency(TTree *events, string informalname, float &jzbeff, float &jzbefferr, bool requireZ) {
251 > void JZBefficiency(TTree *events, string informalname, float &jzbeff, float &jzbefferr, bool requireZ, string addcut="") {
252          TCut kbase("abs(genMll-91.2)<20&&genNjets>2&&genZPt>0&&abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
253 +        cout << "Getting started with JZB efficiency" << endl;
254 +        if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
255          if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
256          TH1F* hLM4 = plotEff(events,kbase,informalname);
257          Int_t bin = hLM4->FindBin(jzbSel); // To get the error
# Line 252 | Line 263 | void JZBefficiency(TTree *events, string
263  
264   //________________________________________________________________________
265   // Effect of energy scale on efficiency
266 < void JZBjetScale(TTree *events, float &jesdown, float &jesup, string informalname,bool requireZ,float syst=0.1, Float_t jzbSelection=-1, TString plotName = "" ) {
266 > void JZBjetScale(TTree *events, float &jesdown, float &jesup, string informalname,bool requireZ,string addcut="",float syst=0.1, Float_t jzbSelection=-1, TString plotName = "" ) {
267          TCut kbase("abs(genMll-91.2)<20&&genZPt>0");
268 +        if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
269          if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
270  
271          TCut ksel("abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
# Line 282 | Line 294 | void JZBjetScale(TTree *events, float &j
294          if(!automatized) dout << "    JESup: " << effp << " (" << (effp-eff)/eff*100. << "%)" << std::endl;
295          if(!automatized) dout << "    central:  " << eff << std::endl;
296          if(!automatized) dout << "    JESdown: " << effm << " (" << (effm-eff)/eff*100. << "%)" << std::endl;
297 <        jesup=(effp-eff)/eff*100.;
298 <        jesdown=(effm-eff)/eff*100.;
297 >        jesup=(effp-eff)/eff;
298 >        jesdown=(effm-eff)/eff;
299   }
300  
301   //________________________________________________________________________
302   // Effect of energy scale on JZB efficiency
303 < void doJZBscale(TTree *events, float &down, float &up, float &syst, float systematic, string informalname, bool requireZ) {
303 > void doJZBscale(TTree *events, float &down, float &up, float &syst, float systematic, string informalname, bool requireZ, string addcut) {
304          
305          TCut kbase("abs(genMll-91.2)<20&&genZPt>0&&genNjets>2");
306 +        if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
307          if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
308          TCut ksel("abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
309          
# Line 307 | Line 320 | void doJZBscale(TTree *events, float &do
320          if(!automatized) dout << "  efficiency at JZB==" << jzbSel*(1.+systematic)  << "(-"<<syst*100<<"%) : " << effp << " (" << ((effp-eff)/eff)*100. << "%)"  << std::endl;
321          if(!automatized) dout << "  efficiency at JZB==" << jzbSel  << ": " << eff << std::endl;
322          if(!automatized) dout << "  efficiency at JZB==" << jzbSel*(1.-systematic)  << "(-"<<syst*100<<"%) : " << effm << " (" << ((effm-eff)/eff)*100. << "%)"  << std::endl;
323 <        up=((effp-eff)/eff)*100;
324 <        down=((effm-eff)/eff)*100;
323 >        up=((effp-eff)/eff);
324 >        down=((effm-eff)/eff);
325   }
326  
327   //________________________________________________________________________
328   // JZB response (true/reco. vs. true)
329 < void JZBresponse(TTree *events, bool requireZ, float &resp, float &resperr, bool isMET = kFALSE, Float_t myJzbMax = 200., Int_t nPeriods = 9 ) {
329 > void JZBresponse(TTree *events, bool requireZ, float &resp, float &resperr, string addcut="",bool isMET = kFALSE, Float_t myJzbMax = 200., Int_t nPeriods = 9 ) {
330          
331          jzbMin = 20;
332          TCut kbase("abs(genMll-91.2)<20&&genZPt>0&&genNjets>2");
333 +        if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
334          if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
335          TCut ksel("abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
336          
# Line 336 | Line 350 | void JZBresponse(TTree *events, bool req
350   }
351  
352  
353 < void do_systematics_for_one_file(TTree *events,string informalname, vector<vector<float> > &results,string mcjzb,string datajzb,float peakerror,bool requireZ=false) {
353 > void do_systematics_for_one_file(TTree *events,int Neventsinfile,string informalname, vector<vector<float> > &results,string mcjzb,string datajzb,float peakerror,bool requireZ=false, string addcut="") {
354    
355    float JetEnergyScaleUncert=0.1;
356    float JZBScaleUncert=0.1;
357    mcjzbexpression=mcjzb;
358 <  
345 <  float triggereff=4;//percent!
358 >  float triggereff=4.0/100;// in range [0,1]
359    dout << "Trigger efficiency not implemented in this script  yet, still using external one" << endl;
360 <  float leptonseleff=2;//percent!
360 >  float leptonseleff=2.0/100;// in range [0,1]
361    dout << "Lepton selection efficiency not implemented in this script  yet, still using external one" << endl;
362    
363    float mceff,mcefferr,jzbeff,jzbefferr;
364    if(!automatized) dout << "MC efficiencies:" << endl;
365 <  MCefficiency(events,mceff,mcefferr,mcjzb,requireZ);
366 <  JZBefficiency(events,informalname,jzbeff,jzbefferr,requireZ);
365 >  MCefficiency(events,mceff,mcefferr,mcjzb,requireZ,Neventsinfile,addcut);
366 >  JZBefficiency(events,informalname,jzbeff,jzbefferr,requireZ,addcut);
367    if(!automatized) dout << "JZB efficiency: " << jzbeff << "+/-" << jzbefferr << endl;
368    
369    if(!automatized) dout << "Error from Peak position:" << endl;
370    float sysfrompeak=0;
371 <  PeakError(events,sysfrompeak,mcjzb,peakerror);
371 >  PeakError(events,sysfrompeak,mcjzb,peakerror,addcut);
372      
373    if(!automatized) dout << "Jet energy scale: " << std::endl;
374    float jesup,jesdown;
375 <  JZBjetScale(events,jesdown,jesup,informalname,JetEnergyScaleUncert,requireZ);
375 >  JZBjetScale(events,jesdown,jesup,informalname,requireZ,addcut,JetEnergyScaleUncert);
376    
377    if(!automatized) dout << "JZB scale: " << std::endl;
378    float scaleup,scaledown,scalesyst;
379 <  doJZBscale(events,scaledown,scaleup,scalesyst,JZBScaleUncert,informalname,requireZ);
379 >  doJZBscale(events,scaledown,scaleup,scalesyst,JZBScaleUncert,informalname,requireZ,addcut);
380    
381    if(!automatized) dout << "JZB response: " << std::endl;
382    float resp,resperr;
383 <  JZBresponse(events,requireZ,resp,resperr);
383 >  JZBresponse(events,requireZ,resp,resperr,addcut);
384  
385    if(!automatized) dout << "Pileup: " << std::endl;
386 <  float resolution=pileup(events,requireZ,informalname);
386 >  float resolution=pileup(events,requireZ,informalname,addcut);
387    
388    dout << "_______________________________________________" << endl;
389 <  dout << "                 SUMMARY FOR " << informalname << " with JZB>" << jzbSel << "  (all in %) " << endl;
390 <  dout << "MC efficiency: " << 100*mceff << "+/-" << 100*mcefferr << endl;
391 <  dout << "Trigger efficiency: " << triggereff << endl;
392 <  dout << "Lepton Sel Eff: " << leptonseleff << endl;
393 <  dout << "Jet energy scale: " << jesup << " " << jesdown << endl;
394 <  dout << "JZB Scale Uncert: " << scaledown << " " << scaleup << endl;
395 <  dout << "Resolution : " << resolution << endl;
396 <  dout << "From peak : " << sysfrompeak << endl;
397 <  dout << "JZB efficiency: " << jzbeff << "+/-" << jzbefferr << " (not yet included below) " << endl;
398 <  dout << "JZB response  : " << resp << " +/-" << resperr << " (not yet included below) " << endl;
389 >  dout << "                 SUMMARY FOR " << informalname << " with JZB>" << jzbSel << "  (all in %) ";
390 >  if(addcut!="") dout << "With additional cut: " << addcut;
391 >  dout << endl;
392 >  dout << "MC efficiency: " << mceff << "+/-" << mcefferr << endl; // in range [0,1]
393 >  dout << "Trigger efficiency: " << triggereff << endl; // in range [0,1]
394 >  dout << "Lepton Sel Eff: " << leptonseleff << endl; // in range [0,1]
395 >  dout << "Jet energy scale: " << jesup << " " << jesdown << endl; // in range [0,1]
396 >  dout << "JZB Scale Uncert: " << scaledown << " " << scaleup << endl; // in range [0,1]
397 >  dout << "Resolution : " << resolution << endl; // in range [0,1]
398 >  dout << "From peak : " << sysfrompeak << endl; // in range [0,1]
399 >  dout << "JZB efficiency: " << jzbeff << "+/-" << jzbefferr << " (not yet included below) " << endl; // in range [0,1]
400 >  dout << "JZB response  : " << resp << " +/-" << resperr << " (not yet included below) " << endl; // in range [0,1]
401    
402    float toterr=0;
403 <  toterr+=(triggereff/100)*(triggereff/100);
404 <  toterr+=(leptonseleff/100)*(leptonseleff/100);
405 <  if(fabs(jesup)>fabs(jesdown)) toterr+=(jesup/100)*(jesup/100); else toterr+=(jesdown/100)*(jesdown/100);
406 <  if(fabs(scaleup)>fabs(scaledown)) toterr+=(scaleup/100)*(scaleup/100); else toterr+=(scaledown/100)*(scaledown/100);
407 <  toterr+=(resolution/100)*(resolution/100);
408 <  toterr+=(sysfrompeak/100)*(sysfrompeak/100);
409 <  toterr=TMath::Sqrt(toterr);
403 >  toterr+=(triggereff)*(triggereff);
404 >  toterr+=(leptonseleff)*(leptonseleff);
405 >  if(fabs(jesup)>fabs(jesdown)) toterr+=(jesup*jesup); else toterr+=(jesdown*jesdown);
406 >  if(fabs(scaleup)>fabs(scaledown)) toterr+=(scaleup*scaleup); else toterr+=(scaledown*scaledown);
407 >  toterr+=(resolution*resolution);
408 >  toterr+=(sysfrompeak*sysfrompeak);
409 >  dout << "TOTAL SYSTEMATICS: " << TMath::Sqrt(toterr) << " --> " << TMath::Sqrt(toterr)*mceff << endl;
410 >  toterr=TMath::Sqrt(toterr)*mceff;
411    dout << "FINAL RESULT : " << 100*mceff << " +/- "<< 100*mcefferr << " (stat) +/- " << 100*toterr << " (syst)   %" << endl;
412 <  dout << "     we thus use the sqrt of the sum of the squares which is : " << 100*TMath::Sqrt(mcefferr*mcefferr+(toterr*toterr)) << endl;
412 >  dout << "     we thus use the sqrt of the sum of the squares which is : " << 100*toterr << endl;
413    vector<float> res;
414    res.push_back(jzbSel);
415    res.push_back(mceff);
416    res.push_back(mcefferr);
417    res.push_back(toterr);
418    res.push_back(TMath::Sqrt((mcefferr)*(mcefferr)+(toterr*toterr)));
419 <  
419 >  if(fabs(jesup)>fabs(jesdown)) res.push_back(fabs(jesup)); else res.push_back(fabs(jesdown));
420 >  if(fabs(scaleup)>fabs(scaledown)) res.push_back(fabs(scaleup)); else res.push_back(fabs(scaledown));
421 >  res.push_back(fabs(resolution));
422    results.push_back(res);
423   }
424  
# Line 412 | Line 430 | vector<vector<float> > compute_systemati
430        for(int ibin=0;ibin<bins.size();ibin++) {
431          jzbSel=bins[ibin];
432          geqleq="geq";
433 <        do_systematics_for_one_file((signalsamples.collection)[isignal].events,(signalsamples.collection)[isignal].samplename,systematics,mcjzb,datajzb,mcpeakerror,requireZ);
433 >        do_systematics_for_one_file((signalsamples.collection)[isignal].events,(signalsamples.collection)[isignal].Nentries,(signalsamples.collection)[isignal].samplename,systematics,mcjzb,datajzb,mcpeakerror,requireZ);
434        }//end of bin loop
435    }//end of signal loop
436    return systematics;

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines