ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/cbrown/AnalysisFramework/Plotting/Modules/Systematics.C

Comparing UserCode/cbrown/AnalysisFramework/Plotting/Modules/Systematics.C (file contents):
Revision 1.41 by buchmann, Wed Oct 19 10:36:31 2011 UTC vs.
Revision 1.46 by buchmann, Wed Nov 9 10:55:25 2011 UTC

#	Line 85 | Line 85 | Double_t Interpolate(Double_t x, TH1 *hi
85
86		//____________________________________________________________________________________
87		// Plotting with all contributions, i.e. sidebands, peak, osof,ossf ... (for a systematic)
88	<	float allcontributionsplot(TTree* events, TCut kBaseCut, TCut kMassCut, TCut kSidebandCut, TCut JZBPosCut, TCut JZBNegCut) {
88	>	float allcontributionsplot(TTree* events, TCut kBaseCut, TCut kMassCut, TCut kSidebandCut, TCut JZBPosCut, TCut JZBNegCut, int flipped) {
89		iplot++;
90		int count=iplot;
91	+	string locmcjzbexpression=mcjzbexpression;
92		// Define new histogram
93		string hname=GetNumericHistoName();
94		TH1F* hossfp = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
95	<	events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kMassCut&&JZBPosCut&&cutOSSF,"goff");
95	>	events->Draw("("+TString(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kMassCut&&JZBPosCut&&cutOSSF,"goff");
96		hname=GetNumericHistoName();
97		TH1F* hossfn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
98	<	events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kMassCut&&JZBNegCut&&cutOSSF,"goff");
98	>	events->Draw("("+TString(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kMassCut&&JZBNegCut&&cutOSSF,"goff");
99
100		hname=GetNumericHistoName();
101		TH1F* hosofp = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
102	<	events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kMassCut&&JZBPosCut&&cutOSOF,"goff");
102	>	events->Draw("("+TString(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kMassCut&&JZBPosCut&&cutOSOF,"goff");
103		hname=GetNumericHistoName();
104		TH1F* hosofn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
105	<	events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kMassCut&&JZBNegCut&&cutOSOF,"goff");
105	<
106	<	hname=GetNumericHistoName();
107	<	TH1F* sbhossfp = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
108	<	events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kSidebandCut&&JZBPosCut&&cutOSSF,"goff");
109	<	hname=GetNumericHistoName();
110	<	TH1F* sbhossfn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
111	<	events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kSidebandCut&&JZBNegCut&&cutOSSF,"goff");
112	<
113	<	hname=GetNumericHistoName();
114	<	TH1F* sbhosofp = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
115	<	events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kSidebandCut&&JZBPosCut&&cutOSOF,"goff");
116	<	hname=GetNumericHistoName();
117	<	TH1F* sbhosofn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
118	<	events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kSidebandCut&&JZBNegCut&&cutOSOF,"goff");
105	>	events->Draw("("+TString(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kMassCut&&JZBNegCut&&cutOSOF,"goff");
106
107	<	float obs = hossfp->Integral();
108	<	float pred= hossfn->Integral() + (1.0/3)*( hosofp->Integral() - hosofn->Integral() + sbhossfp->Integral() - sbhossfn->Integral() + sbhosofp->Integral() - sbhosofn->Integral());
107	>	float obs=0;
108	>	float pred=0;
109	>	flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
110	>	if(PlottingSetup::RestrictToMassPeak) {
111	>	hname=GetNumericHistoName();
112	>	TH1F* sbhossfp = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
113	>	events->Draw("("+TString(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kSidebandCut&&JZBPosCut&&cutOSSF,"goff");
114	>	hname=GetNumericHistoName();
115	>	TH1F* sbhossfn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
116	>	events->Draw("("+TString(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kSidebandCut&&JZBNegCut&&cutOSSF,"goff");
117	>
118	>	hname=GetNumericHistoName();
119	>	TH1F* sbhosofp = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
120	>	events->Draw("("+TString(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kSidebandCut&&JZBPosCut&&cutOSOF,"goff");
121	>	hname=GetNumericHistoName();
122	>	TH1F* sbhosofn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
123	>	events->Draw("("+TString(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kSidebandCut&&JZBNegCut&&cutOSOF,"goff");
124	>
125	>	obs = hossfp->Integral();
126	>	pred= hossfn->Integral() + (1.0/3)*( hosofp->Integral() - hosofn->Integral() + sbhossfp->Integral() - sbhossfn->Integral() + sbhosofp->Integral() - sbhosofn->Integral());
127	>
128	>	if(flipped>0) {
129	>	obs = hossfn->Integral();
130	>	pred= hossfp->Integral() - (1.0/3)*( hosofp->Integral() - hosofn->Integral() + sbhossfp->Integral() - sbhossfn->Integral() + sbhosofp->Integral() - sbhosofn->Integral());
131	>	}
132	>	delete sbhossfp,sbhossfn,sbhosofp,sbhosofn;
133	>	} else {
134	>	obs = hossfp->Integral();
135	>	pred= hossfn->Integral() + (hosofp->Integral() - hosofn->Integral());
136	>	if(flipped>0) {
137	>	obs = hossfn->Integral();
138	>	pred= hossfp->Integral() - (hosofp->Integral() - hosofn->Integral());;
139	>	}
140	>	}
141
142		delete hossfp,hossfn,hosofp,hosofn;
124	–	delete sbhossfp,sbhossfn,sbhosofp,sbhosofn;
143		return obs-pred;
144		}
145
146
147		//____________________________________________________________________________________
148		// Efficiency plot
149	<	TH1F* plotEff(TTree* events, TCut kbase, TString informalname) {
149	>	TH1F* plotEff(TTree* events, TCut kbase, TString informalname, int flipped) {
150		iplot++;
151		int count=iplot;
152		// Define new histogram
#	Line 137 | Line 155 | TH1F* plotEff(TTree* events, TCut kbase,
155		nBins,jzbMin,jzbMax);
156		Float_t step = (jzbMax-jzbMin)/static_cast<Float_t>(nBins);
157
158	<	events->Draw(mcjzbexpression.c_str(),"genJZB>-400"&&kbase,"goff");
158	>	if(flipped==0) events->Draw(mcjzbexpression.c_str(),"genJZB>-400"&&kbase,"goff");
159	>	else events->Draw(("-"+mcjzbexpression).c_str(),"genJZB>-14000"&&kbase,"goff");
160		Float_t maxEff = events->GetSelectedRows();
161		if(verbose>0) dout << hname << " (" << informalname <<") " << maxEff <<  std::endl;
162
#	Line 190 | Line 209 | void master_formula(std::vector<float> e
209
210		//________________________________________________________________________________________
211		// Get normalization factor for the PDFs
212	<	float get_norm_pdf_factor(TTree *events, int k) {
212	>	float get_norm_pdf_factor(TTree *events, int k, string addcut) {
213
214		TH1F *haux = new TH1F("haux", "", 10000, 0, 5);
215		char nameVar[20];
216		sprintf(nameVar, "pdfW[%d]", k);
217	<	events->Project("haux", nameVar);
217	>	events->Project("haux", nameVar, addcut.c_str());
218		float thisW = haux->Integral();
219		events->Project("haux", "pdfW[0]");
220		float normW = haux->Integral();
#	Line 212 | Line 231 | float get_norm_pdf_factor(TTree *events,
231
232		//________________________________________________________________________________________
233		// Pile-up efficiency
234	<	float pileup(TTree *events, bool requireZ, string informalname, string addcut="",Float_t myJzbMax = 140. ) {
234	>	float pileup(TTree *events, bool requireZ, string informalname, int flipped, string addcut="",Float_t myJzbMax = 140. ) {
235		nBins = 16;
236		jzbMax = myJzbMax;
237
238		// Acceptance cuts
239	<	TCut kbase("abs(genMll-91.2)<20&&genNjets>2&&genZPt>0&&abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
239	>	flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
240	>	TCut kbase(PlottingSetup::genMassCut&&"genNjets>2&&genZPt>0"&&cutmass&&cutOSSF);
241		if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
242
243	<	if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
244	<	TH1F* hLM4 = plotEff(events,kbase,informalname);
243	>	if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
244	>	TH1F* hLM4 = plotEff(events,kbase,informalname,flipped);
245		hLM4->SetMinimum(0.);
246
247		// Nominal function
#	Line 242 | Line 262 | float pileup(TTree *events, bool require
262
263		//____________________________________________________________________________________
264		// Effect of peak shifting
265	<	void PeakError(TTree *events,float &result, string mcjzb, float peakerr,string addcut="") {
265	>	void PeakError(TTree *events,float &result, string mcjzb, float peakerr,int flipped,string addcut="") {
266	>	//Note: the cut used here is something like (JZBEXPRESSION+(peakerr)>50) without all the other cuts, to increase statistics (particularly for scans)
267		TString peakup("("+TString(mcjzb)+"+"+TString(any2string(TMath::Abs(peakerr)))+")"+geq_or_leq()+TString(any2string(jzbSel)));
268		TString peakdown("("+TString(mcjzb)+"-"+TString(any2string(TMath::Abs(peakerr)))+")"+geq_or_leq()+TString(any2string(jzbSel)));
269		TString peakcentral("("+TString(mcjzb)+")"+geq_or_leq()+TString(any2string(jzbSel)));
#	Line 265 | Line 286 | void PeakError(TTree *events,float &resu
286		negcut=npeakup;
287		}
288		float res;
289	<	if(addcut=="") res=allcontributionsplot(events,cutnJets,cutmass,sidebandcut,poscut.c_str(),negcut.c_str());
290	<	else res=allcontributionsplot(events,cutnJets&&addcut.c_str(),cutmass,sidebandcut,poscut.c_str(),negcut.c_str());
289	>	if(addcut=="") res=allcontributionsplot(events,cutnJets,cutmass,sidebandcut,poscut.c_str(),negcut.c_str(),flipped);
290	>	else res=allcontributionsplot(events,cutnJets&&addcut.c_str(),cutmass,sidebandcut,poscut.c_str(),negcut.c_str(),flipped);
291		if(i==0) rescent=res;
292		else if(i==1) resdown=res;
293		else if(i==2) resup=res;
#	Line 275 | Line 296 | void PeakError(TTree *events,float &resu
296		else result=(TMath::Abs(rescent-resdown)/(float)rescent);
297		}
298
299	+
300	+	void  MCPartialefficiency(TTree *events,float &result, float &resulterr,int flipped, string mcjzb,bool requireZ,int Neventsinfile, string addcut="", int k = 0, int type = 0) {
301	+	if(!events) {
302	+	write_error(__FUNCTION__,"Tree passed for efficiency calculation is invalid!");
303	+	result=0;resulterr=0;
304	+	return;
305	+	}
306	+
307	+	char jzbSelStr[256]; sprintf(jzbSelStr,"%f",jzbSel);
308	+	// All acceptance cuts at gen. level
309	+	//TCut kbase("abs(genMll-91.2)<20&&genNjets>2&&genZPt>0&&genJZB"+geq_or_leq()+TString(jzbSelStr)+"&&genId1==-genId2");
310	+	TCut kbase("");
311	+	if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
312	+	if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
313	+	// Corresponding reco. cuts
314	+
315	+	TCut acceptance("genPt2 != 0");
316	+	TCut massId(cutmass&&cutOSSF);
317	+	TCut njets(cutnJets);
318	+	TCut jzbp;
319	+	TCut jzbn;
320	+	if(flipped==0) {
321	+	jzbp=TCut((TString(mcjzb)+geq_or_leq()+TString(jzbSelStr)));
322	+	jzbn=TCut((TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr)));
323	+	} else {
324	+	jzbp=TCut(TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr));
325	+	jzbn=TCut(TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
326	+	}
327	+	float ntotal = events->Draw("pt1", kbase, "goff");
328	+	TCut theCut;
329	+	switch(type) {
330	+	case 1:
331	+	theCut = kbase+acceptance;
332	+	break;
333	+	case 2:
334	+	theCut = kbase+massId;
335	+	break;
336	+	case 3:
337	+	theCut = kbase+massId+njets;
338	+	break;
339	+	case 4:
340	+	theCut = kbase+massId+njets+jzbn;
341	+	break;
342	+	default:
343	+	theCut = kbase+massId+njets+jzbn;
344	+	break;
345	+	}
346	+
347	+	string stheCut(theCut);
348	+	char var[20];
349	+	sprintf(var, "pdfW[%d]", k);
350	+
351	+	string svar(var);
352	+	string newtheCut;
353	+	if(k>0) newtheCut = "(" + stheCut + ")*" + svar;
354	+	else newtheCut = "(" + stheCut + ")"; // for k==0 or even k==-1 we don't need to evaluate PDFs
355	+
356	+	TH1F *effh= new TH1F("effh","effh",1,-14000,14000);
357	+	if(k>=0) events->Draw((mcjzbexpression+">>effh").c_str(), newtheCut.c_str(),"goff");
358	+	else events->Draw((mcjzbexpression+">>effh").c_str(), theCut,"goff");
359	+	Float_t sel = effh->Integral();
360	+	Float_t nsel=0;
361	+	//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.
362	+	float normFactor = 1;
363	+	if(k>=0) get_norm_pdf_factor(events, k, addcut);
364	+	sel = sel/normFactor;
365	+
366	+	result=(sel)/ntotal;
367	+	resulterr=TMath::Sqrt(sel/ntotal*(1+sel/ntotal)/ntotal);
368	+
369	+	}
370	+
371		//____________________________________________________________________________________
372		// Total selection efficiency (MC)
373		//returns the efficiency WITHOUT signal contamination, and the result and resulterr contain the result and the corresponding error
374	<	Value MCefficiency(TTree *events,float &result, float &resulterr,string mcjzb,bool requireZ,int Neventsinfile, string addcut="", int k = 0) {
374	>	Value MCefficiency(TTree *events,float &result, float &resulterr, int flipped, string mcjzb,bool requireZ,int Neventsinfile, string addcut="", int k = 0) {
375		if(!events) {
376	<	write_error(__FUNCTION__,"Tree passed for efficiency calculation is invalid!");
377	<	result=0;resulterr=0;
378	<	return Value(0,0);
376	>	write_error(__FUNCTION__,"Tree passed for efficiency calculation is invalid!");
377	>	result=0;
378	>	resulterr=0;
379	>	return Value(0,0);
380		}
381	+
382		char jzbSelStr[256]; sprintf(jzbSelStr,"%f",jzbSel);
383		// All acceptance cuts at gen. level
384		//TCut kbase("abs(genMll-91.2)<20&&genNjets>2&&genZPt>0&&genJZB"+geq_or_leq()+TString(jzbSelStr)+"&&genId1==-genId2");
385		TCut kbase("");
386	<	if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
386	>
387	>	flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
388	>	if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
389		if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
390		// Corresponding reco. cuts
391	<	TCut ksel("pfJetGoodNum>2&&abs(mll-91.2)<20&&id1==id2&&"+TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
392	<	TCut ksel2("pfJetGoodNum>2&&abs(mll-91.2)<20&&id1==id2&&"+TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr));
391	>
392	>	flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
393	>	TCut ksel;//("pfJetGoodNum>2"&&cutmass&&"id1==id2&&"+TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
394	>	TCut ksel2;//("pfJetGoodNum>2"&&cutmass&&"id1==id2&&"+TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr));
395	>	flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
396	>	if(PlottingSetup::RestrictToMassPeak||!ConsiderSignalContaminationForLimits) {
397	>	ksel=TCut("pfJetGoodNum>2"&&cutmass&&"id1==id2&&"+TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
398	>	ksel2=TCut("pfJetGoodNum>2"&&cutmass&&"id1==id2&&"+TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr));
399	>	if(flipped>0) {
400	>	ksel=TCut("pfJetGoodNum>2"&&cutmass&&"id1==id2&&"+TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr));
401	>	ksel2=TCut("pfJetGoodNum>2"&&cutmass&&"id1==id2&&"+TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
402	>	}
403	>	} else {
404	>	//for off peak analysis we don't use the OSSF condition here yet so we can recycle these two cuts for the em condition!
405	>	ksel=TCut("pfJetGoodNum>2"&&cutmass&&(TString(mcjzb)+geq_or_leq()+TString(jzbSelStr)));
406	>	ksel2=TCut("pfJetGoodNum>2"&&cutmass&&(TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr)));
407	>	if(flipped>0) {
408	>	ksel=TCut("pfJetGoodNum>2"&&cutmass&&(TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr)));
409	>	ksel2=TCut("pfJetGoodNum>2"&&cutmass&&(TString(mcjzb)+geq_or_leq()+TString(jzbSelStr)));
410	>	}
411	>	}
412	>
413		TCut posSide = kbase&&ksel;
414		TCut negSide = kbase&&ksel2;
415		string sposSide(posSide);
416		string snegSide(negSide);
417		char var[20];
418		sprintf(var, "pdfW[%d]", k);
419	+	if(k==-1) sprintf(var,"1.0");//case in which we don't want to evaluate PDFs
420		string svar(var);
421	<	string newPosSide = "(" + sposSide + ")*" + svar;
422	<	string newNegSide = "(" + snegSide + ")*" + svar;
421	>	string newPosSide = "((id1==id2)&&(" + sposSide + "))*" + svar;
422	>	string newNegSide = "((id1==id2)&&(" + snegSide + "))*" + svar;
423	>	string emnewPosSide = "((id1!=id2)&&(" + sposSide + "))*" + svar; // only used for off peak analysis
424	>	string emnewNegSide = "((id1!=id2)&&(" + snegSide + "))*" + svar; // only used for off peak analysis
425
426		TH1F *effh= new TH1F("effh","effh",1,-14000,14000);
427		if(k>=0)events->Draw((mcjzbexpression+">>effh").c_str(), newPosSide.c_str(),"goff");
428	<	else events->Draw((mcjzbexpression+">>effh").c_str(), sposSide.c_str(),"goff");
428	>	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.
429		Float_t sel = effh->Integral();
430		Float_t nsel=0;
431	+
432	+	///----------------------------------------------- THIS PART REQUIRES STUDYING! -------------------------
433	+
434		if(ConsiderSignalContaminationForLimits) {
435	<	if(k>=0)events->Draw((mcjzbexpression+">>effh").c_str(), newNegSide.c_str(),"goff");
436	<	else events->Draw((mcjzbexpression+">>effh").c_str(), snegSide.c_str(),"goff");
437	<	nsel = effh->Integral();
435	>	flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
436	>	if(PlottingSetup::RestrictToMassPeak) {
437	>	events->Draw((mcjzbexpression+">>effh").c_str(), newNegSide.c_str(),"goff");
438	>	nsel += effh->Integral();
439	>	} else {
440	>	events->Draw((mcjzbexpression+">>effh").c_str(), newNegSide.c_str(),"goff");
441	>	nsel += effh->Integral();
442	>	events->Draw((mcjzbexpression+">>effh").c_str(), emnewPosSide.c_str(),"goff");
443	>	nsel += effh->Integral();
444	>	events->Draw((mcjzbexpression+">>effh").c_str(), emnewNegSide.c_str(),"goff");
445	>	nsel -= effh->Integral();
446	>	}
447		}
448	+
449		//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.
450		float normFactor = 1;
451	<	if(k>=0) get_norm_pdf_factor(events, k);
451	>	if(k>=0) get_norm_pdf_factor(events, k, addcut);
452		sel = sel/normFactor;
453		nsel = nsel/normFactor;
454
#	Line 334 | Line 467 | Value MCefficiency(TTree *events,float &
467		resulterr=TMath::Sqrt(sel/tot*(1+sel/tot)/tot);
468		result_wo_signalcont=Value(result,resulterr);
469		}
470	<	if(!automatized && k>0 ) dout << "PDF assessment (" << k << ") : ";
471	<	if(!automatized) dout << "  MC efficiency: " << result << "+-" << resulterr << "  ( JZB>" << jzbSel << " : " << sel << " , JZB<-" << jzbSel << " : " << nsel << " and nevents=" << tot << ") with normFact=" << normFactor << std::endl;
470	>	if(!automatized && k>0 ) dout << "PDF assessment [" << k << "] : ";
471	>	if(!automatized) dout << "  MC efficiency: " << result << "+-" << resulterr << "  ( JZB>" << jzbSel << " : " << sel << " , signal contamination : " << nsel << " and nevents=" << tot << ") with normFact=" << normFactor << std::endl;
472		delete effh;
473		return result_wo_signalcont;
474		}
475
476
477	+
478		//____________________________________________________________________________________
479		// Selection efficiency for one process (MC)
480	<	vector<float> processMCefficiency(TTree *events,string mcjzb,bool requireZ,int Neventsinfile, string addcut) {
480	>	// not in use anymore.
481	>	/*
482	>	vector<float> processMCefficiency(TTree *events,int flipped, string mcjzb,bool requireZ,int Neventsinfile, string addcut) {
483		vector<float> process_efficiencies;
484		for(int iprocess=0;iprocess<=10;iprocess++) {
485		float this_process_efficiency,efferr;
486		stringstream addcutplus;
487		addcutplus<<addcut<<"&&(process=="<<iprocess<<")";
488	<	MCefficiency(events,this_process_efficiency, efferr,mcjzb,requireZ,Neventsinfile, addcutplus.str(),-1);
488	>	MCefficiency(events,this_process_efficiency, efferr,flipped,mcjzb,requireZ,Neventsinfile, addcutplus.str(),-1);
489		process_efficiencies.push_back(this_process_efficiency);
490		}
491		return process_efficiencies;
492		}
493	<
493	>	*/
494
495	<	void JZBefficiency(TTree *events, string informalname, float &jzbeff, float &jzbefferr, bool requireZ, string addcut="") {
496	<	TCut kbase("abs(genMll-91.2)<20&&genNjets>2&&genZPt>0&&abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
495	>	void JZBefficiency(TTree *events, string informalname, float &jzbeff, float &jzbefferr, int flipped, bool requireZ, string addcut="") {
496	>	TCut kbase(genMassCut&&"genNjets>2&&genZPt>0"&&cutmass&&cutOSSF);
497		if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
498	<	if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
499	<	TH1F* hLM4 = plotEff(events,kbase,informalname);
498	>	if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
499	>	TH1F* hLM4 = plotEff(events,kbase,informalname,flipped);
500		Int_t bin = hLM4->FindBin(jzbSel); // To get the error
501		jzbeff=Interpolate(jzbSel,hLM4);
502		jzbefferr=hLM4->GetBinError(bin);
#	Line 370 | Line 506 | void JZBefficiency(TTree *events, string
506
507		//________________________________________________________________________
508		// Effect of energy scale on efficiency
509	<	void JZBjetScale(TTree *events, float &jesdown, float &jesup, string informalname,bool requireZ,string addcut="",float syst=0.1, Float_t jzbSelection=-1, TString plotName = "" ) {
510	<	TCut kbase("abs(genMll-91.2)<20&&genZPt>0");
509	>	void JZBjetScale(TTree *events, float &jesdown, float &jesup, string informalname, int flipped, bool requireZ,string addcut="",float syst=0.1, Float_t jzbSelection=-1, TString plotName = "" ) {
510	>	TCut kbase(genMassCut&&"genZPt>0");
511		if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
512	<	if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
512	>	flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
513	>	if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
514
515	<	TCut ksel("abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
515	>	TCut ksel(cutmass&&cutOSSF);
516		TCut nJets("pfJetGoodNum>2");
517		stringstream down,up;
518		down << "pfJetGoodNum"<<30*(1-syst)<<">=3";
#	Line 387 | Line 524 | void JZBjetScale(TTree *events, float &j
524		if ( !(plotName.Length()>1) ) plotName = informalname;
525
526		nBins = 1; jzbMin = jzbSel*0.95; jzbMax = jzbSel*1.05;
527	<	TH1F* hist = plotEff(events,(kbase&&ksel&&nJets),informalname);
527	>	TH1F* hist = plotEff(events,(kbase&&ksel&&nJets),informalname,flipped);
528
529	<	TH1F* histp = plotEff(events,(kbase&&ksel&&nJetsP),informalname);
529	>	TH1F* histp = plotEff(events,(kbase&&ksel&&nJetsP),informalname,flipped);
530
531	<	TH1F* histm = plotEff(events,(kbase&&ksel&&nJetsM),informalname);
531	>	TH1F* histm = plotEff(events,(kbase&&ksel&&nJetsM),informalname,flipped);
532
533		// Dump some information
534		Float_t eff  = Interpolate(jzbSel,hist);
#	Line 407 | Line 544 | void JZBjetScale(TTree *events, float &j
544
545		//________________________________________________________________________
546		// Effect of energy scale on JZB efficiency
547	<	void doJZBscale(TTree *events, float &down, float &up, float &syst, float systematic, string informalname, bool requireZ, string addcut) {
547	>	void doJZBscale(TTree *events, float &down, float &up, float &syst, float systematic, string informalname, int flipped, bool requireZ, string addcut) {
548
549	<	TCut kbase("abs(genMll-91.2)<20&&genZPt>0&&genNjets>2");
549	>	TCut kbase(genMassCut&&"genZPt>0&&genNjets>2");
550		if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
551	<	if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
552	<	TCut ksel("abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
551	>	flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
552	>	if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
553	>	TCut ksel(cutmass&&cutOSSF);
554
555		nBins =    50;
556		jzbMin =   0.5*jzbSel;
557		jzbMax =   2.0*jzbSel;
558
559	<	TH1F* hist = plotEff(events,kbase&&ksel,informalname);
559	>	TH1F* hist = plotEff(events,kbase&&ksel,informalname,flipped);
560
561		// Dump some information
562		Float_t eff  = Interpolate(jzbSel,hist);
#	Line 433 | Line 571 | void doJZBscale(TTree *events, float &do
571
572		//________________________________________________________________________
573		// JZB response (true/reco. vs. true)
574	<	void JZBresponse(TTree *events, bool requireZ, float &resp, float &resperr, string addcut="",bool isMET = kFALSE, Float_t myJzbMax = 200., Int_t nPeriods = 9 ) {
574	>	void JZBresponse(TTree *events, bool requireZ, float &resp, float &resperr, int flipped, string addcut="", bool isMET = kFALSE, Float_t myJzbMax = 200., Int_t nPeriods = 9 ) {
575
576		jzbMin = 20;
577	<	TCut kbase("abs(genMll-91.2)<20&&genZPt>0&&genNjets>2");
577	>	flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
578	>	TCut kbase(genMassCut&&"genZPt>0&&genNjets>2");
579		if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
580	<	if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
581	<	TCut ksel("abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
580	>	flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
581	>	if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
582	>	flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
583	>	TCut ksel(cutmass&&cutOSSF);
584
585		TProfile* hJzbResp = new TProfile("hJzbResp","JZB response  ; JZB true (GeV/c); JZB reco. / JZB true", nPeriods, jzbMin, myJzbMax, "" );
586
587	<	if (!isMET) events->Project("hJzbResp","("+TString(mcjzbexpression)+")/genJZB:genJZB",kbase&&ksel);
587	>	string locmcjzbexpression=mcjzbexpression;
588	>	if(flipped>0) locmcjzbexpression="-"+locmcjzbexpression;
589	>	string possibleminus="";
590	>	if(flipped>0) possibleminus="-";
591	>	if (!isMET) events->Project("hJzbResp","("+TString(locmcjzbexpression)+")/("+possibleminus+"genJZB):("+possibleminus+"genJZB)",kbase&&ksel);
592		else events->Project("hJzbResp","met[4]/genMET:genMET",kbase&&ksel);
593
594		hJzbResp->SetMaximum(1.2);
#	Line 466 | Line 611 | void JZBresponse(TTree *events, bool req
611
612		//________________________________________________________________________________________
613		// PDF uncertainty
614	<	float get_pdf_uncertainty(TTree *events, string mcjzb, bool requireZ, int Neventsinfile, int NPdfs, string addcut="") {
614	>	float get_pdf_uncertainty(TTree *events, int flipped, string mcjzb, bool requireZ, int Neventsinfile, int NPdfs, string addcut="") {
615		std::vector<float> efficiency;
616		for(int k = 1; k < NPdfs; k++) {
617		float result, resulterr;
618	<	MCefficiency(events, result, resulterr, mcjzb, requireZ, Neventsinfile, addcut, k);
618	>	Value flipval;
619	>	MCefficiency(events, result, resulterr, flipped, mcjzb, requireZ, Neventsinfile, addcut, k);
620		efficiency.push_back(result);
621		}
622		float errHi, errLow,err;
#	Line 490 | Line 636 | int get_npdfs(TTree *events) {
636		}
637
638
639	<	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="", bool ismSUGRA=false) {
639	>	void do_systematics_for_one_file(TTree *events,int Neventsinfile,string informalname, vector<vector<float> > &results,int flipped, string mcjzb,string datajzb,float peakerror,bool requireZ=false, string addcut="", bool ismSUGRA=false) {
640		float JetEnergyScaleUncert=0.1;
641		float JZBScaleUncert=0.1;
642		mcjzbexpression=mcjzb;
#	Line 505 | Line 651 | void do_systematics_for_one_file(TTree *
651
652		float mceff,mcefferr,jzbeff,jzbefferr;
653		if(!automatized) dout << "MC efficiencies:" << endl;
654	<	Value mceff_nosigcont = MCefficiency(events,mceff,mcefferr,mcjzb,requireZ,Neventsinfile,addcut,-1);
654	>	Value flipefficiency;
655	>	Value mceff_nosigcont = MCefficiency(events,mceff,mcefferr,flipped,mcjzb,requireZ,Neventsinfile,addcut,-1);
656		if(!automatized) cout << "   Without signal contamination, we find an efficiency of " << mceff_nosigcont << endl;
657
658	<	if(PlottingSetup::computeJZBefficiency) JZBefficiency(events,informalname,jzbeff,jzbefferr,requireZ,addcut);
658	>	if(PlottingSetup::computeJZBefficiency) JZBefficiency(events,informalname,jzbeff,jzbefferr,flipped,requireZ,addcut);
659		if(!automatized) dout << "JZB efficiency: " << jzbeff << "+/-" << jzbefferr << endl;
660
514	–
515	–	if(!(mceff>0)) write_warning(__FUNCTION__,"Efficiency is zero - the systematics will not be computed!");
516	–
517	–
661		if(!automatized) dout << "Error from Peak position:" << endl;
662		float sysfrompeak=0;
663	<	if(mceff>0) PeakError(events,sysfrompeak,mcjzb,peakerror,addcut);
521	<	else dout << "Not computed." << endl;
663	>	PeakError(events,sysfrompeak,mcjzb,peakerror,flipped,addcut);
664
665		if(!automatized) dout << "Jet energy scale: " << std::endl;
666	<	float jesup=0,jesdown=0;
667	<	if(mceff>0) JZBjetScale(events,jesdown,jesup,informalname,requireZ,addcut,JetEnergyScaleUncert);
526	<	else dout << "Not computed." << endl;
666	>	float jesup,jesdown;
667	>	JZBjetScale(events,jesdown,jesup,informalname,flipped,requireZ,addcut,JetEnergyScaleUncert);
668
669		if(!automatized) dout << "JZB scale: " << std::endl;
670	<	float scaleup=0,scaledown=0,scalesyst=0;
671	<	if(mceff>0) doJZBscale(events,scaledown,scaleup,scalesyst,JZBScaleUncert,informalname,requireZ,addcut);
531	<	else dout << "Not computed." << endl;
670	>	float scaleup,scaledown,scalesyst;
671	>	doJZBscale(events,scaledown,scaleup,scalesyst,JZBScaleUncert,informalname,flipped,requireZ,addcut);
672
673		if(!automatized) dout << "JZB response: " << std::endl;
674	<	float resp=0,resperr=0;
675	<	if(PlottingSetup::computeJZBresponse && mceff>0) {
674	>	float resp,resperr;
675	>	if(PlottingSetup::computeJZBresponse) {
676		if(!automatized) dout << "JZB response: " << std::endl;
677	<	JZBresponse(events,requireZ,resp,resperr,addcut);
677	>	JZBresponse(events,requireZ,resp,resperr,flipped,addcut);
678		}
679
680		if(!automatized) dout << "Pileup: " << std::endl;
681	<	float resolution=0;
682	<	if(mceff>0) resolution=pileup(events,requireZ,informalname,addcut);
543	<	else dout << "Not computed." << endl;
681	>	float resolution;
682	>	resolution=pileup(events,requireZ,informalname,flipped,addcut);
683
684		float PDFuncert=0;
685	<	if(!automatized&&mceff>0) dout << "Assessing PDF uncertainty: " << std::endl;
686	<	if(ismSUGRA&&mceff>0) PDFuncert = get_pdf_uncertainty(events, mcjzb, requireZ, Neventsinfile, NPdfs, addcut);
685	>	if(!automatized) dout << "Assessing PDF uncertainty: " << std::endl;
686	>	if(ismSUGRA) PDFuncert = get_pdf_uncertainty(events, flipped, mcjzb, requireZ, Neventsinfile, NPdfs, addcut);
687
688		dout << "_______________________________________________" << endl;
689		dout << "                 SUMMARY FOR " << informalname << " with JZB>" << jzbSel << "  (all in %) ";
#	Line 593 | Line 732 | void do_systematics_for_one_file(TTree *
732		results.push_back(res);
733		}
734
735	<	vector<vector<float> > compute_systematics(string mcjzb, float mcpeakerror, string datajzb, samplecollection &signalsamples, vector<float> bins, bool requireZ=false) {
735	>	vector<vector<float> > compute_systematics(string mcjzb, float mcpeakerror, int flipped, string datajzb, samplecollection &signalsamples, vector<float> bins, bool requireZ=false) {
736		automatized=true;
737		vector< vector<float> > systematics;
738		for (int isignal=0; isignal<signalsamples.collection.size();isignal++) {
#	Line 601 | Line 740 | vector<vector<float> > compute_systemati
740		for(int ibin=0;ibin<bins.size();ibin++) {
741		jzbSel=bins[ibin];
742		geqleq="geq";
743	<	do_systematics_for_one_file((signalsamples.collection)[isignal].events,(signalsamples.collection)[isignal].Nentries,(signalsamples.collection)[isignal].samplename,systematics,mcjzb,datajzb,mcpeakerror,requireZ);
743	>	do_systematics_for_one_file((signalsamples.collection)[isignal].events,(signalsamples.collection)[isignal].Nentries,(signalsamples.collection)[isignal].samplename,systematics,flipped,mcjzb,datajzb,mcpeakerror,requireZ);
744		}//end of bin loop
745		}//end of signal loop
746		return systematics;

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines