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.7 by buchmann, Wed Jul 20 13:58:05 2011 UTC vs.
Revision 1.17 by buchmann, Wed Aug 17 14:55:51 2011 UTC

# Line 158 | Line 158 | TH1F* plotEff(TTree* events, TCut kbase,
158  
159   //________________________________________________________________________________________
160   // Pile-up efficiency
161 < float pileup(TTree *events, 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.);
172          
# Line 186 | Line 188 | float pileup(TTree *events, string infor
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;
219        
223   }
224  
225   //____________________________________________________________________________________
226   // Total selection efficiency (MC)
227 < void MCefficiency(TTree *events,float &result, float &resulterr,string mcjzb) {
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 < float JZBefficiency(TTree *events, string informalname) {
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
258 +        jzbeff=Interpolate(jzbSel,hLM4);
259 +        jzbefferr=hLM4->GetBinError(bin);
260          if(!automatized) dout << "  Efficiency at JZB==" << jzbSel  << std::endl;
261 <        if(!automatized) dout << "    " << Interpolate(jzbSel,hLM4) << "+-" << hLM4->GetBinError(bin)  << std::endl;
248 <        return -1;
261 >        if(!automatized) dout << "    " << jzbeff << "+-" << jzbefferr  << std::endl;
262   }
263  
264   //________________________________________________________________________
265   // Effect of energy scale on efficiency
266 < void JZBjetScale(TTree *events, float &jesdown, float &jesup, string informalname="",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");
272          TCut nJets("pfJetGoodNum>2");
273          stringstream down,up;
# Line 284 | Line 300 | void JZBjetScale(TTree *events, float &j
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) {
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          
310          nBins =    50;
# Line 308 | Line 326 | void doJZBscale(TTree *events, float &do
326  
327   //________________________________________________________________________
328   // JZB response (true/reco. vs. true)
329 < void JZBresponse(TTree *events, 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          
337 <        TProfile* hJzbResp = new TProfile("hJzbResp","JZB response  ; JZB true (GeV/c); JZB reco. / JZB true",
318 <                                                                          nPeriods, jzbMin, myJzbMax, "" );
337 >        TProfile* hJzbResp = new TProfile("hJzbResp","JZB response  ; JZB true (GeV/c); JZB reco. / JZB true", nPeriods, jzbMin, myJzbMax, "" );
338          
339          if (!isMET) events->Project("hJzbResp","("+TString(mcjzbexpression)+")/genJZB:genJZB",kbase&&ksel);
340          else events->Project("hJzbResp","met[4]/genMET:genMET",kbase&&ksel);
# Line 324 | Line 343 | void JZBresponse(TTree *events, bool isM
343          hJzbResp->SetMinimum(0.2);
344          hJzbResp->Fit("pol0","Q");
345          TF1 *fittedfunction = hJzbResp->GetFunction("pol0");
346 <        if(!automatized) dout << "  Response: " << fittedfunction->GetParameter(0) << " +/- " << fittedfunction->GetParError(0) << endl;
346 >        resp=fittedfunction->GetParameter(0);
347 >        resperr=fittedfunction->GetParError(0);
348 >        if(!automatized) dout << "  Response: " << resp << " +/- " << resperr << endl;
349          delete hJzbResp;
350   }
351  
352  
353 < void do_systematics_for_one_file(TTree *events,string informalname, vector<vector<float> > &results,string mcjzb,string datajzb,float peakerror) {
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;
337  
358    float triggereff=4;//percent!
359    dout << "Trigger efficiency not implemented in this script  yet, still using external one" << endl;
360    float leptonseleff=2;//percent!
361    dout << "Lepton selection efficiency not implemented in this script  yet, still using external one" << endl;
362    
363 <  float mceff,mcefferr;
363 >  float mceff,mcefferr,jzbeff,jzbefferr;
364    if(!automatized) dout << "MC efficiencies:" << endl;
365 <  MCefficiency(events,mceff,mcefferr,mcjzb);
366 <  JZBefficiency(events,informalname);
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);
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);
379 >  doJZBscale(events,scaledown,scaleup,scalesyst,JZBScaleUncert,informalname,requireZ,addcut);
380    
381    if(!automatized) dout << "JZB response: " << std::endl;
382 <  JZBresponse(events);
382 >  float resp,resperr;
383 >  JZBresponse(events,requireZ,resp,resperr,addcut);
384  
385    if(!automatized) dout << "Pileup: " << std::endl;
386 <  float resolution=pileup(events,informalname);
386 >  float resolution=pileup(events,requireZ,informalname,addcut);
387    
388    dout << "_______________________________________________" << endl;
389 <  dout << "                 SUMMARY FOR " << informalname << " with JZB>" << jzbSel << "  (all in %) " << 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: " << 100*mceff << "+/-" << 100*mcefferr << endl;
393    dout << "Trigger efficiency: " << triggereff << endl;
394    dout << "Lepton Sel Eff: " << leptonseleff << endl;
# Line 372 | Line 396 | void do_systematics_for_one_file(TTree *
396    dout << "JZB Scale Uncert: " << scaledown << " " << scaleup << endl;
397    dout << "Resolution : " << resolution << endl;
398    dout << "From peak : " << sysfrompeak << endl;
399 <  
399 >  dout << "JZB efficiency: " << jzbeff << "+/-" << jzbefferr << " (not yet included below) " << endl;
400 >  dout << "JZB response  : " << resp << " +/-" << resperr << " (not yet included below) " << endl;
401    
402    float toterr=0;
403    toterr+=(triggereff/100)*(triggereff/100);
# Line 390 | Line 415 | void do_systematics_for_one_file(TTree *
415    res.push_back(mcefferr);
416    res.push_back(toterr);
417    res.push_back(TMath::Sqrt((mcefferr)*(mcefferr)+(toterr*toterr)));
418 <  
418 >  if(fabs(jesup)>fabs(jesdown)) res.push_back(fabs(jesup/100)); else res.push_back(fabs(jesdown)/100);
419 >  if(fabs(scaleup)>fabs(scaledown)) res.push_back(fabs(scaleup)/100); else res.push_back(fabs(scaledown)/100);
420 >  res.push_back(fabs(resolution)/100);
421    results.push_back(res);
422   }
423  
424 < vector<vector<float> > compute_systematics(string mcjzb, float mcpeakerror, string datajzb, samplecollection &signalsamples, vector<float> bins) {
424 > vector<vector<float> > compute_systematics(string mcjzb, float mcpeakerror, string datajzb, samplecollection &signalsamples, vector<float> bins, bool requireZ=false) {
425    automatized=true;
426    vector< vector<float> > systematics;
427    for (int isignal=0; isignal<signalsamples.collection.size();isignal++) {
# Line 402 | Line 429 | vector<vector<float> > compute_systemati
429        for(int ibin=0;ibin<bins.size();ibin++) {
430          jzbSel=bins[ibin];
431          geqleq="geq";
432 <        do_systematics_for_one_file((signalsamples.collection)[isignal].events,(signalsamples.collection)[isignal].samplename,systematics,mcjzb,datajzb,mcpeakerror);
432 >        do_systematics_for_one_file((signalsamples.collection)[isignal].events,(signalsamples.collection)[isignal].Nentries,(signalsamples.collection)[isignal].samplename,systematics,mcjzb,datajzb,mcpeakerror,requireZ);
433        }//end of bin loop
434    }//end of signal loop
435    return systematics;

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines