ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/ShapeLimit.C

Comparing UserCode/cbrown/Development/Plotting/Modules/ShapeLimit.C (file contents):
Revision 1.4 by buchmann, Wed Apr 25 12:59:53 2012 UTC vs.
Revision 1.10 by buchmann, Mon May 7 16:44:27 2012 UTC

#	Line 38 | Line 38 | const char* strip_flip_away(string flipp
38		return flipped_name.c_str();
39		}
40
41	<	void prepare_limit_datacard(string RunDirectory, TFile *f, float uJSU, float uPDF, int flipped) {
42	<	TH1F *dataob;
43	<	if(flipped) dataob = (TH1F*)f->Get("flipped_data_obs");
44	<	else dataob = (TH1F*)f->Get("data_obs");
45	<	TH1F *signal;
46	<
47	<	if(flipped) signal = (TH1F*)f->Get("flipped_signal");
48	<	else signal = (TH1F*)f->Get("signal");
49	<
50	<	TH1F *background1;
51	<	if(flipped) background1 = (TH1F*)f->Get("flipped_TTbarBackground");
52	<	else background1 = (TH1F*)f->Get("TTbarBackground");
53	<
54	<	TH1F *background2;
55	<	if(flipped) background2 = (TH1F*)f->Get("flipped_ZJetsBackground");
56	<	else background2 = (TH1F*)f->Get("ZJetsBackground");
41	>	void EliminateNegativeEntries(TH1F *histo) {
42	>	for(int i=1;i<=histo->GetNbinsX();i++) {
43	>	if(histo->GetBinContent(i)<0) histo->SetBinContent(i,0);
44	>	}
45	>	}
46	>
47	>	void prepare_limit_datacard(string RunDirectory, TFile *f, float uJSU, float uPDF) {
48	>	TH1F *dataob = (TH1F*)f->Get("data_obs");
49	>	TH1F *signal = (TH1F*)f->Get("signal");
50	>	TH1F *Tbackground = (TH1F*)f->Get("TTbarBackground");
51	>	TH1F *Zbackground  = (TH1F*)f->Get("ZJetsBackground");
52
53		assert(dataob);
54		assert(signal);
55	<	assert(background1);
56	<	assert(background2);
55	>	assert(Tbackground);
56	>	assert(Zbackground);
57
63	–	ofstream datacard;
64	–	write_warning(__FUNCTION__,"Need to rethink systematics we want to use !");
65	–	datacard.open ((RunDirectory+"/susydatacard.txt").c_str());
66	–	datacard << "Writing this to a file.\n";
67	–	datacard << "imax 1\n"; // number of channels
68	–	datacard << "jmax 2\n"; // number of backgrounds (Z+Jets prediction and TTbar prediction)
69	–	datacard << "kmax *\n"; // number of nuisance parameters (sources of systematic uncertainties)
70	–	datacard << "---------------\n";
71	–	datacard << "shapes * * limitfile.root $PROCESS $PROCESS_$SYSTEMATIC\n";
72	–	datacard << "---------------\n";
73	–	datacard << "bin 1\n";
74	–	datacard << "observation "<<dataob->Integral()<<"\n";
75	–	datacard << "------------------------------\n";
76	–	datacard << "bin             1          1          1\n";
77	–	datacard << "process         signal     TTbarBackground     ZJetsBackground\n";
78	–	datacard << "process         0          1          2\n";
79	–	datacard << "rate            "<<signal->Integral()<<"         "<<background1->Integral()<<"         "<<background2->Integral()<<"\n";
80	–	datacard << "--------------------------------\n";
81	–	datacard << "lumi     lnN    " << 1+ PlottingSetup::lumiuncert << "       -       -    luminosity uncertainty\n"; // only affects MC -> signal!
82	–	datacard << "Stat   shape    -          1          -    statistical uncertainty (ttbar)\n";
83	–	datacard << "Stat   shape    -          -          1    statistical uncertainty (zjets)\n";
84	–	datacard << "Sys    shape    -          1          -    systematic uncertainty on ttbar\n";
85	–	datacard << "Sys    shape    -          -          1    systematic uncertainty on zjets\n";
86	–	datacard << "Stat   shape    1          -          -    statistical uncertainty (signal)\n";
87	–	datacard << "JES    shape    1          -          -    uncertainty on Jet Energy Scale\n";
88	–	datacard << "JSU      lnN    " << 1+uJSU << "          -          -    JZB Scale Uncertainty\n";
89	–	if(uPDF>0) datacard << "PDF      lnN    " << 1+uPDF << "          -          -    uncertainty from PDFs\n";
58
59	<	// datacard << "peak   shape    1          1    uncertainty on signal resolution.
60	<	datacard.close();
59	>	if(FullMCAnalysis) {
60	>	//dealing with MC analysis - need to use shapes.
61	>	write_warning(__FUNCTION__,"Not correctly implemented yet for MC");
62	>	ofstream datacard;
63	>	write_warning(__FUNCTION__,"Need to rethink systematics we want to use !");
64	>	datacard.open ((RunDirectory+"/susydatacard.txt").c_str());
65	>	datacard << "Writing this to a file.\n";
66	>	datacard << "imax 1\n"; // number of channels
67	>	datacard << "jmax 2\n"; // number of backgrounds (Z+Jets prediction and TTbar prediction)
68	>	datacard << "kmax *\n"; // number of nuisance parameters (sources of systematic uncertainties)
69	>	datacard << "---------------\n";
70	>	datacard << "shapes * * limitfile.root $PROCESS $PROCESS_$SYSTEMATIC\n";
71	>	datacard << "---------------\n";
72	>	datacard << "bin 1\n";
73	>	datacard << "observation "<<dataob->Integral()<<"\n";
74	>	datacard << "------------------------------\n";
75	>	datacard << "bin             1          1          1\n";
76	>	datacard << "process         signal     TTbarBackground     ZJetsBackground\n";
77	>	datacard << "process         0          1          2\n";
78	>	datacard << "rate            "<<signal->Integral()<<"         "<<Tbackground->Integral()<<"         "<<Zbackground->Integral()<<"\n";
79	>	datacard << "--------------------------------\n";
80	>	datacard << "lumi     lnN    " << 1+ PlottingSetup::lumiuncert << "       -       -    luminosity uncertainty\n"; // only affects MC -> signal!
81	>	datacard << "Stat   shape    -          1          -    statistical uncertainty (ttbar)\n";
82	>	datacard << "Stat   shape    -          -          1    statistical uncertainty (zjets)\n";
83	>	datacard << "Sys    shape    -          1          -    systematic uncertainty on ttbar\n";
84	>	datacard << "Sys    shape    -          -          1    systematic uncertainty on zjets\n";
85	>	datacard << "Stat   shape    1          -          -    statistical uncertainty (signal)\n";
86	>	datacard << "JES    shape    1          -          -    uncertainty on Jet Energy Scale\n";
87	>	datacard << "JSU      lnN    " << 1+uJSU << "          -          -    JZB Scale Uncertainty\n";
88	>	if(uPDF>0) datacard << "PDF      lnN    " << 1+uPDF << "          -          -    uncertainty from PDFs\n";
89	>	datacard.close();
90	>	} else {
91	>	//doing mutibin analysis
92	>	ofstream datacard;
93	>	write_warning(__FUNCTION__,"Need to rethink systematics we want to use !");
94	>	datacard.open ((RunDirectory+"/susydatacard.txt").c_str());
95	>	datacard << "Writing this to a file.\n";
96	>	datacard << "imax " << dataob->GetNbinsX() << "\n"; // number of channels
97	>	datacard << "jmax 2\n"; // number of backgrounds (Z+Jets prediction and TTbar prediction)
98	>	datacard << "kmax *\n"; // number of nuisance parameters (sources of systematic uncertainties)
99	>	datacard << "---------------\n";
100	>	datacard << "bin\t";
101	>	for(int i=1;i<=dataob->GetNbinsX();i++) datacard << " " << dataob->GetBinLowEdge(i) << "to" << dataob->GetBinLowEdge(i)+dataob->GetBinWidth(i) << " \t";
102	>	datacard << "\n";
103	>	datacard << "observation\t";
104	>	for(int i=1;i<=dataob->GetNbinsX();i++) datacard << " " << dataob->GetBinContent(i) << " \t";
105	>	datacard<<"\n";
106	>	datacard << "------------------------------\n";
107	>	datacard << "bin\t";
108	>	for(int i=1;i<=dataob->GetNbinsX();i++) {
109	>	datacard << " " << dataob->GetBinLowEdge(i) << "to" << dataob->GetBinLowEdge(i)+dataob->GetBinWidth(i) << "\t" <<
110	>	" " << dataob->GetBinLowEdge(i) << "to" << dataob->GetBinLowEdge(i)+dataob->GetBinWidth(i) << "\t" <<
111	>	" " << dataob->GetBinLowEdge(i) << "to" << dataob->GetBinLowEdge(i)+dataob->GetBinWidth(i) << "\t";
112	>	}
113	>	datacard << "\n";
114	>	datacard << "process\t";
115	>	for(int i=1;i<=dataob->GetNbinsX();i++) datacard << "\t signal \t TTbarBackground \t ZJetsBackground";
116	>	datacard << "\n";
117	>	datacard << "process\t";
118	>	for(int i=1;i<=dataob->GetNbinsX();i++) datacard << "\t 0 \t 1 \t 2";
119	>	datacard << "\n";
120	>
121	>
122	>	datacard << "rate\t";
123	>	for(int i=1;i<=dataob->GetNbinsX();i++) datacard << "\t " << signal->GetBinContent(i) << " \t " << Tbackground->GetBinContent(i) << " \t " << Zbackground->GetBinContent(i) << "\t";
124	>	datacard<<"\n";
125	>
126	>	datacard << "--------------------------------\n";
127	>
128	>
129	>	datacard << "lumi     lnN    \t";
130	>	for(int i=1;i<=dataob->GetNbinsX();i++) datacard << " " << 1+ PlottingSetup::lumiuncert << "\t - \t -";
131	>	datacard << "       luminosity uncertainty\n"; // only affects MC -> signal!
132	>
133	>	// Statistical uncertainty
134	>	datacard << "Stat     lnN    \t";
135	>	for(int i=1;i<=dataob->GetNbinsX();i++) {
136	>	//Signal
137	>	float central = signal->GetBinContent(i);
138	>	float up      = ((TH1F*)f->Get("signal_StatDown"))->GetBinContent(i);
139	>	float down    = ((TH1F*)f->Get("signal_StatUp"))->GetBinContent(i);
140	>	float suncert=0;
141	>	if(central>0) {
142	>	if(abs(up-central)>abs(down-central)) suncert=abs(up-central)/central;
143	>	else suncert=abs(central-down)/central;
144	>	}
145	>
146	>	//TTbar
147	>	central = Tbackground->GetBinContent(i);
148	>	up      = ((TH1F*)f->Get("TTbarBackground_StatUp"))->GetBinContent(i);
149	>	down    = ((TH1F*)f->Get("TTbarBackground_StatDown"))->GetBinContent(i);
150	>	float tuncert=0;
151	>	if(central>0) {
152	>	if(abs(up-central)>abs(down-central)) tuncert=abs(up-central)/central;
153	>	else tuncert=abs(central-down)/central;
154	>	}
155	>	//ZJets
156	>	central = Zbackground->GetBinContent(i);
157	>	up      = ((TH1F*)f->Get("ZJetsBackground_StatUp"))->GetBinContent(i);
158	>	down    = ((TH1F*)f->Get("ZJetsBackground_StatDown"))->GetBinContent(i);
159	>	float zuncert=0;
160	>	if(central>0) {
161	>	if(abs(up-central)>abs(down-central)) zuncert=abs(up-central)/central;
162	>	else zuncert=abs(central-down)/central;
163	>	}
164	>	datacard << " " << 1+suncert << " \t " << 1+tuncert << "\t " << 1+zuncert << " \t ";
165	>	}
166	>
167	>	datacard << "       statistical uncertainty\n";
168	>
169	>	// Statistical uncertainty
170	>	datacard << "Sys     lnN    \t";
171	>	for(int i=1;i<=dataob->GetNbinsX();i++) {
172	>	float central = Tbackground->GetBinContent(i);
173	>	float up      = ((TH1F*)f->Get("TTbarBackground_SysUp"))->GetBinContent(i);
174	>	float down    = ((TH1F*)f->Get("TTbarBackground_SysDown"))->GetBinContent(i);
175	>	float tuncert=0;
176	>	if(central>0) {
177	>	if(abs(up-central)>abs(down-central)) tuncert=abs(up-central)/central;
178	>	else tuncert=abs(central-down)/central;
179	>	}
180	>	central = Zbackground->GetBinContent(i);
181	>	up      = ((TH1F*)f->Get("ZJetsBackground_SysUp"))->GetBinContent(i);
182	>	down    = ((TH1F*)f->Get("ZJetsBackground_SysDown"))->GetBinContent(i);
183	>	float zuncert=0;
184	>	if(central>0) {
185	>	if(abs(up-central)>abs(down-central)) zuncert=abs(up-central)/central;
186	>	else zuncert=abs(central-down)/central;
187	>	}
188	>	datacard << " - \t " << 1+tuncert << "\t " << 1+zuncert << " \t ";
189	>	}
190	>
191	>	datacard << "       systematic uncertainty\n";
192	>
193	>
194	>	float JESup = ((TH1F*)f->Get("signal_JESUp"))->Integral();
195	>	float JESdn = ((TH1F*)f->Get("signal_JESDown"))->Integral();
196	>	float central = signal->Integral();
197	>	float uJES=0;
198	>	if(abs(JESup-central)>abs(JESdn-central)) uJES=abs(JESup-central)/central;
199	>	else uJES=abs(central-JESdn)/central;
200	>	datacard << "JES    lnN";
201	>	for(int i=1;i<=dataob->GetNbinsX();i++) datacard << " " << 1+uJES << "\t - \t  - \t";
202	>	datacard << "uncertainty on Jet Energy Scale\n";
203	>
204	>	datacard << "JSU      lnN    ";
205	>	for(int i=1;i<=dataob->GetNbinsX();i++) datacard << " " << 1+uJSU << "\t - \t  - \t";
206	>	datacard << "JZB Scale Uncertainty\n";
207	>
208	>	if(uPDF>0) {
209	>	datacard << "PDF      lnN    ";
210	>	for(int i=1;i<=dataob->GetNbinsX();i++) datacard << " " << 1+uPDF << "\t - \t - \t";
211	>	datacard << "uncertainty from PDFs\n";
212	>	}
213	>
214	>	datacard.close();
215	>	}
216	>
217		}
218
219		float QuickDrawNevents=0;
#	Line 122 | Line 246 | TH1F* QuickDraw(TTree *events, string hn
246		return histo;
247		}
248
125	–	/*void do_stat_up(TH1F *h, float sign=1.0) {
126	–	for(int i=1;i<=h->GetNbinsX();i++) {
127	–	h->SetBinContent(i,h->GetBinContent(i)+sign*h->GetBinError(i));
128	–	}
129	–	h->Write();
130	–	}
131	–
132	–	void do_stat_dn(TH1F *h) {
133	–	do_stat_up(h,-1.0);
134	–	}*/
249
250		void SQRT(TH1F *h) {
251		for (int i=1;i<=h->GetNbinsX();i++) {
#	Line 176 | Line 290 | void generate_shapes_for_systematic(bool
290		}
291
292		if(signalonly) {
293	<	cout << "Processing a signal with mcjzb: " << mcjzb << " (identifier: " << identifier << ")" << endl;
293	>	dout << "Processing a signal with mcjzb: " << mcjzb << " (identifier: '" << identifier << "')" << endl;
294		TH1F *ZOSSFP = QuickDraw(signalevents,"ZOSSFP",mcjzb,binning,     "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
181	–	TH1F *ZOSOFP = QuickDraw(signalevents,"ZOSOFP",mcjzb,binning,     "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
295		TH1F *ZOSSFN = QuickDraw(signalevents,"ZOSSFN","-"+mcjzb,binning, "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
296	<	TH1F *ZOSOFN = QuickDraw(signalevents,"ZOSOFN","-"+mcjzb,binning, "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
296	>	TH1F *ZOSOFP;
297	>	TH1F *ZOSOFN;
298	>
299	>	if(!PlottingSetup::FullMCAnalysis) {
300	>	ZOSOFP = QuickDraw(signalevents,"ZOSOFP",mcjzb,binning,     "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
301	>	ZOSOFN = QuickDraw(signalevents,"ZOSOFN","-"+mcjzb,binning, "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
302	>	}
303
304		TH1F *SBOSSFP;
305		TH1F *SBOSOFP;
306		TH1F *SBOSSFN;
307		TH1F *SBOSOFN;
308
309	<	if(PlottingSetup::RestrictToMassPeak) {
309	>	if(PlottingSetup::RestrictToMassPeak&&!PlottingSetup::FullMCAnalysis) {
310		SBOSSFP = QuickDraw(signalevents,"SBOSSFP",mcjzb,binning, "JZB", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
311		SBOSOFP = QuickDraw(signalevents,"SBOSOFP",mcjzb,binning, "JZB", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
312		SBOSSFN = QuickDraw(signalevents,"SBOSSFN","-"+mcjzb,binning, "JZB", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
#	Line 206 | Line 325 | void generate_shapes_for_systematic(bool
325		TH1F *flippedLpred = new TH1F("flippedLpred","flippedLpred",binning.size()-1,&binning[0]);
326
327		Lobs->Add(ZOSSFP);
328	<	Lpred->Add(ZOSSFN);
328	>	if(!PlottingSetup::FullMCAnalysis) Lpred->Add(ZOSSFN);
329
330	<	cout << "SITUATION FOR SIGNAL: " << endl;
331	<	cout << " OSSF     JZB> 0 : " << ZOSSFP->Integral() << "     JZB < 0 :" << ZOSSFN->Integral() << endl;
332	<	cout << " OSOF     JZB> 0 : " << ZOSOFP->Integral() << "     JZB < 0 :" << ZOSOFN->Integral() << endl;
333	<	if(PlottingSetup::RestrictToMassPeak) {
334	<	cout << " OSSF SB  JZB> 0 : " << SBOSSFP->Integral() << "     JZB < 0 :" << SBOSSFN->Integral() << endl;
335	<	cout << " OSOF SB  JZB> 0 : " << SBOSOFP->Integral() << "     JZB < 0 :" << SBOSOFN->Integral() << endl;
330	>	dout << "SITUATION FOR SIGNAL: " << endl;
331	>
332	>
333	>	if(PlottingSetup::FullMCAnalysis) {
334	>	dout << " OSSF     JZB> 0 : " << ZOSSFP->Integral() << endl;
335	>	} else {
336	>	dout << " OSSF     JZB> 0 : " << ZOSSFP->Integral() << "     JZB < 0 :" << ZOSSFN->Integral() << endl;
337	>	dout << " OSOF     JZB> 0 : " << ZOSOFP->Integral() << "     JZB < 0 :" << ZOSOFN->Integral() << endl;
338	>	if(PlottingSetup::RestrictToMassPeak&&!PlottingSetup::FullMCAnalysis) {
339	>	dout << " OSSF SB  JZB> 0 : " << SBOSSFP->Integral() << "     JZB < 0 :" << SBOSSFN->Integral() << endl;
340	>	dout << " OSOF SB  JZB> 0 : " << SBOSOFP->Integral() << "     JZB < 0 :" << SBOSOFN->Integral() << endl;
341	>	}
342		}
343
344
345		flippedLobs->Add(ZOSSFN);
346	<	flippedLpred->Add(ZOSSFP);
346	>	if(!PlottingSetup::FullMCAnalysis) flippedLpred->Add(ZOSSFP);
347
348	<	if(PlottingSetup::RestrictToMassPeak) {
349	<	Lpred->Add(ZOSOFP,1.0/3);
350	<	Lpred->Add(ZOSOFN,-1.0/3);
351	<	Lpred->Add(SBOSSFP,1.0/3);
352	<	Lpred->Add(SBOSSFN,-1.0/3);
353	<	Lpred->Add(SBOSOFP,1.0/3);
354	<	Lpred->Add(SBOSOFN,-1.0/3);
355	<
356	<	//flipped prediction
357	<	flippedLpred->Add(ZOSOFP,-1.0/3);
358	<	flippedLpred->Add(ZOSOFN,1.0/3);
359	<	flippedLpred->Add(SBOSSFP,-1.0/3);
360	<	flippedLpred->Add(SBOSSFN,1.0/3);
361	<	flippedLpred->Add(SBOSOFP,-1.0/3);
362	<	flippedLpred->Add(SBOSOFN,1.0/3);
363	<
364	<	} else {
365	<	Lpred->Add(ZOSOFP,1.0);
366	<	Lpred->Add(ZOSOFN,-1.0);
367	<
368	<	//flipped prediction
369	<	flippedLpred->Add(ZOSOFP,-1.0);
370	<	flippedLpred->Add(ZOSOFN,1.0);
348	>	if(!PlottingSetup::FullMCAnalysis) {
349	>	if(PlottingSetup::RestrictToMassPeak) {
350	>	Lpred->Add(ZOSOFP,1.0/3);
351	>	Lpred->Add(ZOSOFN,-1.0/3);
352	>	Lpred->Add(SBOSSFP,1.0/3);
353	>	Lpred->Add(SBOSSFN,-1.0/3);
354	>	Lpred->Add(SBOSOFP,1.0/3);
355	>	Lpred->Add(SBOSOFN,-1.0/3);
356	>
357	>	//flipped prediction
358	>	flippedLpred->Add(ZOSOFP,-1.0/3);
359	>	flippedLpred->Add(ZOSOFN,1.0/3);
360	>	flippedLpred->Add(SBOSSFP,-1.0/3);
361	>	flippedLpred->Add(SBOSSFN,1.0/3);
362	>	flippedLpred->Add(SBOSOFP,-1.0/3);
363	>	flippedLpred->Add(SBOSOFN,1.0/3);
364	>	} else {
365	>	Lpred->Add(ZOSOFP,1.0);
366	>	Lpred->Add(ZOSOFN,-1.0);
367	>
368	>	//flipped prediction
369	>	flippedLpred->Add(ZOSOFP,-1.0);
370	>	flippedLpred->Add(ZOSOFN,1.0);
371	>	}
372		}
373
374		TH1F *signal = (TH1F*)Lobs->Clone("signal");
375	<	signal->Add(Lpred,-1);
375	>	if(!PlottingSetup::FullMCAnalysis) signal->Add(Lpred,-1);
376		signal->SetName(signalname.c_str());
377		signal->SetTitle(signalname.c_str());
378		signal->Write();
379
380		TH1F *flippedsignal = (TH1F*)flippedLobs->Clone();
381	<	flippedsignal->Add(flippedLpred,-1);
381	>	if(!PlottingSetup::FullMCAnalysis) flippedsignal->Add(flippedLpred,-1);
382		flippedsignal->SetName(("flipped_"+signalname).c_str());
383		flippedsignal->Write();
384
#	Line 266 | Line 392 | void generate_shapes_for_systematic(bool
392		signalStatDn->SetTitle(((string)signal->GetTitle()+"_StatDown").c_str());
393
394		for(int i=1;i<=signalStatDn->GetNbinsX();i++) {
395	<	float staterr = TMath::Sqrt(Lpred->GetBinContent(i) + Lobs->GetBinContent(i));
396	<	cout << "Stat err in bin " << i << " : " << staterr << endl;
397	<	cout << "    prediction: " << Lpred->GetBinContent(i) << " , observation: " << Lobs->GetBinContent(i) << " --> signal: " << signal->GetBinContent(i) << endl;
398	<	cout << "    we obtain : " << signal->GetBinContent(i)-staterr << " , " << signal->GetBinContent(i)+staterr << endl;
399	<	signalStatDn->SetBinContent(i,signal->GetBinContent(i)-staterr);
395	>	float staterr;
396	>	if(!PlottingSetup::FullMCAnalysis) staterr = TMath::Sqrt(Lpred->GetBinContent(i) + Lobs->GetBinContent(i));
397	>	else staterr = TMath::Sqrt(Lobs->GetBinContent(i));
398	>
399	>	if(!PlottingSetup::FullMCAnalysis) {
400	>	dout << "Stat err in bin " << i << " : " << staterr << endl;
401	>	dout << "    prediction: " << Lpred->GetBinContent(i) << " , observation: " << Lobs->GetBinContent(i) << " --> signal: " << signal->GetBinContent(i) << endl;
402	>	dout << "    we obtain : " << signal->GetBinContent(i)-staterr << " , " << signal->GetBinContent(i)+staterr << endl;
403	>	}
404	>	if(signal->GetBinContent(i)-staterr>0) signalStatDn->SetBinContent(i,signal->GetBinContent(i)-staterr);
405	>	else signalStatDn->SetBinContent(i,0);
406		signalStatUp->SetBinContent(i,signal->GetBinContent(i)+staterr);
407		signal->SetBinError(i,staterr);
408		}
#	Line 290 | Line 422 | void generate_shapes_for_systematic(bool
422		flippedsignalStatDn->SetTitle(((string)flippedsignal->GetTitle()+"_StatDown").c_str());
423
424		for(int i=1;i<=flippedsignalStatDn->GetNbinsX();i++) {
425	<	float staterr = TMath::Sqrt(flippedLpred->GetBinContent(i) + flippedLobs->GetBinContent(i));
426	<	flippedsignalStatDn->SetBinContent(i,flippedsignal->GetBinContent(i)-staterr);
425	>	float staterr;
426	>	if(!PlottingSetup::FullMCAnalysis) staterr = TMath::Sqrt(flippedLpred->GetBinContent(i) + flippedLobs->GetBinContent(i));
427	>	else staterr = TMath::Sqrt(flippedLobs->GetBinContent(i));
428	>	if(flippedsignal->GetBinContent(i)-staterr>0) flippedsignalStatDn->SetBinContent(i,flippedsignal->GetBinContent(i)-staterr);
429	>	else flippedsignalStatDn->SetBinContent(i,0);
430		flippedsignalStatUp->SetBinContent(i,flippedsignal->GetBinContent(i)+staterr);
431		flippedsignal->SetBinError(i,staterr);
432		}
#	Line 323 | Line 458 | void generate_shapes_for_systematic(bool
458
459
460		if(dataonly) {
461	<	cout << "Processing data with datajzb: " << datajzb << endl;
461	>	dout << "Processing data with datajzb: " << datajzb << endl;
462		TH1F *ZOSSFP = allsamples.Draw("ZOSSFP",datajzb,binning, "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,data,luminosity);
463		TH1F *ZOSOFP = allsamples.Draw("ZOSOFP",datajzb,binning, "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,data,luminosity);
464		TH1F *ZOSSFN = allsamples.Draw("ZOSSFN","-"+datajzb,binning, "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,data,luminosity);
#	Line 457 | Line 592 | void generate_shapes_for_systematic(bool
592		SQRT(predstaterr);
593		TH1F *bgStatUp = (TH1F*)pred->Clone("background_StatUp");
594		bgStatUp->Add(predstaterr);
595	+	EliminateNegativeEntries(bgStatUp);
596		bgStatUp->Write();
597		TH1F *bgStatDn = (TH1F*)pred->Clone("background_StatDown");
598		bgStatDn->Add(predstaterr,-1);
599	+	EliminateNegativeEntries(bgStatDn);
600		bgStatDn->Write();
601		//      delete bgStatDn;
602		//      delete bgStatUp;
#	Line 476 | Line 613 | void generate_shapes_for_systematic(bool
613		SQRT(flippedpredstaterr);
614		TH1F *fbgStatUp = (TH1F*)flippedpred->Clone("flipped_background_StatUp");
615		fbgStatUp->Add(predstaterr);
616	+	EliminateNegativeEntries(fbgStatUp);
617		fbgStatUp->Write();
618		TH1F *fbgStatDn = (TH1F*)flippedpred->Clone("flipped_background_StatDown");
619		fbgStatDn->Add(predstaterr,-1);
620	+	EliminateNegativeEntries(fbgStatDn);
621		fbgStatDn->Write();
622		//      delete fbgStatDn;
623		//      delete fbgStatUp;
#	Line 491 | Line 630 | void generate_shapes_for_systematic(bool
630		SQRT(Tpredstaterr);
631		TH1F *TpredStatUp = (TH1F*)Tpred->Clone("TTbarBackground_StatUp");
632		TpredStatUp->Add(Tpredstaterr);
633	+	EliminateNegativeEntries(TpredStatUp);
634		TpredStatUp->Write();
635		TH1F *TpredStatDn = (TH1F*)Tpred->Clone("TTbarBackground_StatDown");
636		TpredStatDn->Add(Tpredstaterr,-1);
637	+	EliminateNegativeEntries(TpredStatDn);
638		TpredStatDn->Write();
639		//      delete TpredStatDn;
640		//      delete TpredStatUp;
#	Line 506 | Line 647 | void generate_shapes_for_systematic(bool
647		SQRT(fTpredstaterr);
648		TH1F *fTpredStatUp = (TH1F*)flippedTpred->Clone("flipped_TTbarBackground_StatUp");
649		fTpredStatUp->Add(fTpredstaterr);
650	+	EliminateNegativeEntries(fTpredStatUp);
651		fTpredStatUp->Write();
652		TH1F *fTpredStatDn = (TH1F*)flippedTpred->Clone("flipped_TTbarBackground_StatDown");
653		fTpredStatDn->Add(fTpredstaterr,-1);
654	+	EliminateNegativeEntries(fTpredStatDn);
655		fTpredStatDn->Write();
656		//      delete fTpredStatDn;
657		//      delete fTpredStatUp;
#	Line 522 | Line 665 | void generate_shapes_for_systematic(bool
665		SQRT(Zpredstaterr);
666		TH1F *ZpredStatUp = (TH1F*)Zpred->Clone("ZJetsBackground_StatUp");
667		ZpredStatUp->Add(Zpredstaterr);
668	+	EliminateNegativeEntries(ZpredStatUp);
669		ZpredStatUp->Write();
670		TH1F *ZpredStatDn = (TH1F*)Zpred->Clone("ZJetsBackground_StatDown");
671		ZpredStatDn->Add(Zpredstaterr,-1);
672	+	EliminateNegativeEntries(ZpredStatDn);
673		ZpredStatDn->Write();
674		//      delete ZpredStatDn;
675		//      delete ZpredStatUp;
#	Line 537 | Line 682 | void generate_shapes_for_systematic(bool
682		SQRT(fTpredstaterr);
683		TH1F *fZpredStatUp = (TH1F*)flippedZpred->Clone("flipped_ZJetsBackground_StatUp");
684		fZpredStatUp->Add(fZpredstaterr);
685	+	EliminateNegativeEntries(fZpredStatUp);
686		fZpredStatUp->Write();
687		TH1F *fZpredStatDn = (TH1F*)flippedZpred->Clone("flipped_ZJetsBackground_StatDown");
688		fZpredStatDn->Add(fZpredstaterr,-1);
689	+	EliminateNegativeEntries(fZpredStatDn);
690		fZpredStatDn->Write();
691		//      delete fZpredStatDn;
692		//      delete fZpredStatUp;
#	Line 559 | Line 706 | void generate_shapes_for_systematic(bool
706		SQRT(predsyserr);
707		TH1F *bgSysUp = (TH1F*)pred->Clone("background_SysUp");
708		bgSysUp->Add(predsyserr);
709	+	EliminateNegativeEntries(bgSysUp);
710		bgSysUp->Write();
711		TH1F *bgSysDn = (TH1F*)pred->Clone("background_SysDown");
712		bgSysDn->Add(predsyserr,-1);
713	+	EliminateNegativeEntries(bgSysDn);
714		bgSysDn->Write();
715		delete predsyserr;
716
#	Line 577 | Line 726 | void generate_shapes_for_systematic(bool
726		SQRT(fpredsyserr);
727		TH1F *fbgSysUp = (TH1F*)flippedpred->Clone("flipped_background_SysUp");
728		fbgSysUp->Add(fpredsyserr);
729	+	EliminateNegativeEntries(fbgSysUp);
730		fbgSysUp->Write();
731		TH1F *fbgSysDn = (TH1F*)flippedpred->Clone("flipped_background_SysDown");
732		fbgSysDn->Add(fpredsyserr,-1);
733	+	EliminateNegativeEntries(fbgSysDn);
734		fbgSysDn->Write();
735		delete fpredsyserr;
736
#	Line 592 | Line 743 | void generate_shapes_for_systematic(bool
743		SQRT(Tpredsyserr);
744		TH1F *TpredSysUp = (TH1F*)Tpred->Clone("TTbarBackground_SysUp");
745		TpredSysUp->Add(Tpredsyserr);
746	+	EliminateNegativeEntries(TpredSysUp);
747		TpredSysUp->Write();
748		TH1F *TpredSysDn = (TH1F*)Tpred->Clone("TTbarBackground_SysDown");
749		TpredSysDn->Add(Tpredsyserr,-1);
750	+	EliminateNegativeEntries(TpredSysDn);
751		TpredSysDn->Write();
752		delete Tpredsyserr;
753
#	Line 606 | Line 759 | void generate_shapes_for_systematic(bool
759		SQRT(fTpredsyserr);
760		TH1F *fTpredSysUp = (TH1F*)flippedTpred->Clone("flipped_TTbarBackground_SysUp");
761		fTpredSysUp->Add(fTpredsyserr);
762	+	EliminateNegativeEntries(fTpredSysUp);
763		fTpredSysUp->Write();
764		TH1F *fTpredSysDn = (TH1F*)flippedTpred->Clone("flipped_TTbarBackground_SysDown");
765		fTpredSysDn->Add(fTpredsyserr,-1);
766	+	EliminateNegativeEntries(fTpredSysDn);
767		fTpredSysDn->Write();
768		delete fTpredsyserr;
769
#	Line 623 | Line 778 | void generate_shapes_for_systematic(bool
778		SQRT(Zpredsyserr);
779		TH1F *ZpredSysUp = (TH1F*)Zpred->Clone("ZJetsBackground_SysUp");
780		ZpredSysUp->Add(Zpredsyserr);
781	+	EliminateNegativeEntries(ZpredSysUp);
782		ZpredSysUp->Write();
783		TH1F *ZpredSysDn = (TH1F*)Zpred->Clone("ZJetsBackground_SysDown");
784		ZpredSysDn->Add(Zpredsyserr,-1);
785	+	EliminateNegativeEntries(ZpredSysDn);
786		ZpredSysDn->Write();
787		delete Zpredsyserr;
788
#	Line 638 | Line 795 | void generate_shapes_for_systematic(bool
795		SQRT(fZpredsyserr);
796		TH1F *fZpredSysUp = (TH1F*)flippedZpred->Clone("flipped_ZJetsBackground_SysUp");
797		fZpredSysUp->Add(fZpredsyserr);
798	+	EliminateNegativeEntries(fZpredSysUp);
799		fZpredSysUp->Write();
800		TH1F *fZpredSysDn = (TH1F*)flippedZpred->Clone("flipped_ZJetsBackground_SysDown");
801		fZpredSysDn->Add(fZpredsyserr,-1);
802	+	EliminateNegativeEntries(fZpredSysDn);
803		fZpredSysDn->Write();
804		delete fZpredsyserr;
805		}
806
807		/*if(identifier=="") {
808		for(int i=0;i<binning.size()-1;i++) {
809	<	cout << "[ " << binning[i] << " , " << binning[i+1] << "] : O " << obs->GetBinContent(i+1) << " P " << pred->GetBinContent(i+1) << " (Z: " << Zpred->GetBinContent(i+1) << " , T: " << Tpred->GetBinContent(i+1) << ")" << endl;
809	>	dout << "[ " << binning[i] << " , " << binning[i+1] << "] : O " << obs->GetBinContent(i+1) << " P " << pred->GetBinContent(i+1) << " (Z: " << Zpred->GetBinContent(i+1) << " , T: " << Tpred->GetBinContent(i+1) << ")" << endl;
810		}
811		}*/
812		delete ZOSSFP;
#	Line 679 | Line 838 | ShapeDroplet LimitsFromShapes(bool asymp
838
839		ensure_directory_exists(RunDirectory);
840
841	<	TFile *datafile = new TFile("../StoredShapes.root","READ");
841	>	TFile *datafile = new TFile((PlottingSetup::cbafbasedir+"/DistributedModelCalculations/StoredShapes.root").c_str(),"READ");
842		if(datafile->IsZombie()) {
843		write_error(__FUNCTION__,"Fatal error: The stored shapes are not available!");
844		assert(!datafile->IsZombie());
845		}
846	<	cout << "Run Directory: " << RunDirectory << endl;
846	>	dout << "Run Directory: " << RunDirectory << endl;
847		TFile *limfile = new TFile((RunDirectory+"/PRElimitfile.root").c_str(),"RECREATE");
848
849		TIter nextkey(datafile->GetListOfKeys());
#	Line 702 | Line 861 | ShapeDroplet LimitsFromShapes(bool asymp
861		bool dataonly=false;
862
863		generate_shapes_for_systematic(signalonly,dataonly,limfile,events,"",mcjzb,datajzb,noJES,jzbbins,limcan,addcut,xsec);
705	–	//  generate_shapes_for_systematic(signalonly,dataonly,limfile,events,"peakUp",newjzbexpression(mcjzb,jzbpeakerrormc),newjzbexpression(datajzb,jzbpeakerrordata),noJES,jzbbins,limcan,addcut,xsec);
706	–	//  generate_shapes_for_systematic(signalonly,dataonly,limfile,events,"peakDown",newjzbexpression(mcjzb,-jzbpeakerrormc),newjzbexpression(datajzb,-jzbpeakerrordata),noJES,jzbbins,limcan,addcut,xsec);
864		generate_shapes_for_systematic(signalonly,dataonly,limfile,events,"JESUp",mcjzb,datajzb,JESup,jzbbins,limcan,addcut,xsec);
865		generate_shapes_for_systematic(signalonly,dataonly,limfile,events,"JESDown",mcjzb,datajzb,JESdown,jzbbins,limcan,addcut,xsec);
866
#	Line 717 | Line 874 | ShapeDroplet LimitsFromShapes(bool asymp
874		bool docomplicatedmSUGRAxsreweighting=false; //if you modify this value please also adapt it in Systematics.C not only here in ShapeLimit.C
875
876		MCefficiency(events,mceff,mcefferr,flipped,mcjzb,requireZ,Neventsinfile,SUSYScanSpace::SUSYscantype,xsec,addcut,-1);
877	<	if(mceff<0) flipped=1;
877	>	if(mceff<0) {
878	>	flipped=1;
879	>	write_info(__FUNCTION__,"Doing flipping!");
880	>	}
881		doJZBscale(events,SUSYScanSpace::SUSYscantype==mSUGRA&&docomplicatedmSUGRAxsreweighting,xsec,scaledown,scaleup,scalesyst,JZBScaleUncert,informalname,flipped,requireZ,addcut);
882		float PDFuncert=0;
883		int NPdfs = get_npdfs(events);
#	Line 763 | Line 923 | ShapeDroplet LimitsFromShapes(bool asymp
923		obj->Write();
924		}
925
926	<	prepare_limit_datacard(RunDirectory,final_limfile,JZBscale,PDFuncert,flipped);
926	>	prepare_limit_datacard(RunDirectory,final_limfile,JZBscale,PDFuncert);
927
928		final_limfile->Close();
929		limfile->Close();
930		dout << "Info: Shape root file and datacard have been generated in " << RunDirectory << endl;
931		stringstream command;
932		if(asymptotic) {
933	<	if(firstGuess>0) command << "bash CreateModel.sh " << RunDirectory << " susydatacard.txt" << " 1 0 0 " << firstGuess; // ASYMPTOTIC LIMITS WITH FIRST GUESS
934	<	else command << "bash CreateModel.sh " << RunDirectory << " susydatacard.txt" << " 0 0 0 0"; // ASYMPTOTIC LIMITS
933	>	if(firstGuess>0) command << "bash " << PlottingSetup::cbafbasedir<< "/DistributedModelCalculations/ShapeLimits/CreateModel.sh " << RunDirectory << " susydatacard.txt" << " 1 0 0 " << firstGuess; // ASYMPTOTIC LIMITS WITH FIRST GUESS
934	>	else command << "bash " << PlottingSetup::cbafbasedir<< "/DistributedModelCalculations/ShapeLimits/CreateModel.sh " << RunDirectory << " susydatacard.txt" << " 0 0 0 0"; // ASYMPTOTIC LIMITS
935		}
936	<	else command << "bash CreateModel.sh " << RunDirectory << " susydatacard.txt" << " 2 0 " << int(0.5 * firstGuess) << " " << int(2*firstGuess); // ASYMPTOTIC LIMITS
936	>	else command << "bash " << PlottingSetup::cbafbasedir<< "/DistributedModelCalculations/ShapeLimits/CreateModel.sh " << RunDirectory << " susydatacard.txt" << " 2 0 " << int(0.333 * firstGuess) << " " << int(firstGuess); // ASYMPTOTIC LIMITS
937		dout <<"Going to run : " << command.str() << endl;
938		int CreatedModelFileExitCode = gSystem->Exec(command.str().c_str());
939	<	cout << "exit code of limit algorithm (CreateModel.sh) : " << CreatedModelFileExitCode << endl;
940	<	assert(CreatedModelFileExitCode==0);
939	>	dout << "exit code of limit algorithm (CreateModel.sh) : " << CreatedModelFileExitCode << endl;
940	>	if(!(CreatedModelFileExitCode==0)) {
941	>	write_warning(__FUNCTION__,"Something bad happened. It looks like a shape analysis is not the way to go. ");
942	>	ShapeDroplet alpha;
943	>	alpha.observed=-12345; // this is the flag to say watch out something went wrong with the signal ...
944	>	alpha.SignalIntegral=1;
945	>	return alpha;
946	>	}
947		ShapeDroplet alpha;
948		alpha.readDroplet(RunDirectory+"/ShapeDropletResult.txt");
949		alpha.PDF=PDFuncert;
#	Line 787 | Line 953 | ShapeDroplet LimitsFromShapes(bool asymp
953		dout << alpha << endl;
954
955		dout << "Everything is saved in " << RunDirectory << endl;
956	<	dout << "Will transfer model and datacard over for possible post-processing" << endl;
957	<	dout << "   1) Make sure models directory exists ... " << std::flush;
956	>	dout << "Cleaning up ... " << std::flush;
957	>	/*  dout << "   1) Make sure models directory exists ... " << std::flush;
958		gSystem->Exec("mkdir -p models/");
959		dout << " ok!" << endl;
960		dout << "   2) Deleting any previous model files with the same name ... " << std::flush;
#	Line 810 | Line 976 | ShapeDroplet LimitsFromShapes(bool asymp
976		copycommand3 << "cp " << RunDirectory << "/limitfile.root models/model_" << name << "_histo.root";
977		gSystem->Exec(copycommand3.str().c_str());
978		dout << " ok!" << endl;
979	<	dout << "   4) Removing original working directory (" << RunDirectory << ") ... " << std::flush;
980	<	write_warning(__FUNCTION__,"Watch out : need to uncomment the line below to remove the original working directory again");
815	<	//  gSystem->Exec(("rm -r "+RunDirectory).c_str());
979	>	dout << "   4) Removing original working directory (" << RunDirectory << ") ... " << std::flush;*/
980	>	gSystem->Exec(("rm -r "+RunDirectory).c_str());
981		dout << " ok!" << endl;
982		delete limcan;
983		return alpha;
#	Line 827 | Line 992 | void PrepareDataShapes(string datajzb, v
992		bool dataonly=true;
993		bool signalonly=false;
994		string mcjzb="JZBforMCinPrepareDataShapes";//this string is not used.
995	<	float jzbpeakerrormc=0;
995	>	//  float jzbpeakerrormc=0;
996		generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"",mcjzb,datajzb,noJES,jzbbins,limcan,"",xsec);
997		// don't need these effects for obs & pred, only for signal!
998		//  generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"peakUp",newjzbexpression(mcjzb,jzbpeakerrormc),newjzbexpression(datajzb,jzbpeakerrordata),noJES,jzbbins,limcan,"",xsec);

Diff Legend

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