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.11 by buchmann, Mon May 14 15:14:03 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");
57	<
58	<	assert(dataob);
59	<	assert(signal);
60	<	assert(background1);
61	<	assert(background2);
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_MC_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	>	assert(dataob);
51	>	assert(signal);
52	>
53	>	write_warning(__FUNCTION__,"The following two defs for th1's are fake");TH1F *Tbackground; TH1F *Zbackground;
54	>	write_warning(__FUNCTION__,"Not correctly implemented yet for MC");
55	>	ofstream datacard;
56	>	write_warning(__FUNCTION__,"Need to rethink systematics we want to use !");
57	>	datacard.open ((RunDirectory+"/susydatacard.txt").c_str());
58	>	dout << "Writing this to a file.\n";
59	>	datacard << "imax 1\n"; // number of channels
60	>	datacard << "jmax 2\n"; // number of backgrounds (Z+Jets prediction and TTbar prediction)
61	>	datacard << "kmax *\n"; // number of nuisance parameters (sources of systematic uncertainties)
62	>	datacard << "---------------\n";
63	>	datacard << "shapes * * limitfile.root $PROCESS $PROCESS_$SYSTEMATIC\n";
64	>	datacard << "---------------\n";
65	>	datacard << "bin 1\n";
66	>	datacard << "observation "<<dataob->Integral()<<"\n";
67	>	datacard << "------------------------------\n";
68	>	datacard << "bin             1          1          1\n";
69	>	datacard << "process         signal     TTbarBackground     ZJetsBackground\n";
70	>	datacard << "process         0          1          2\n";
71	>	datacard << "rate            "<<signal->Integral()<<"         "<<Tbackground->Integral()<<"         "<<Zbackground->Integral()<<"\n";
72	>	datacard << "--------------------------------\n";
73	>	datacard << "lumi     lnN    " << 1+ PlottingSetup::lumiuncert << "       -       -    luminosity uncertainty\n"; // only affects MC -> signal!
74	>	datacard << "Stat   shape    -          1          -    statistical uncertainty (ttbar)\n";
75	>	datacard << "Stat   shape    -          -          1    statistical uncertainty (zjets)\n";
76	>	datacard << "Sys    shape    -          1          -    systematic uncertainty on ttbar\n";
77	>	datacard << "Sys    shape    -          -          1    systematic uncertainty on zjets\n";
78	>	datacard << "Stat   shape    1          -          -    statistical uncertainty (signal)\n";
79	>	datacard << "JES    shape    1          -          -    uncertainty on Jet Energy Scale\n";
80	>	datacard << "JSU      lnN    " << 1+uJSU << "          -          -    JZB Scale Uncertainty\n";
81	>	if(uPDF>0) datacard << "PDF      lnN    " << 1+uPDF << "          -          -    uncertainty from PDFs\n";
82	>	datacard.close();
83	>	}
84	>
85	>	void prepare_datadriven_datacard(string RunDirectory, TFile *f, float uJSU, float uPDF) {
86	>	TH1F *dataob = (TH1F*)f->Get("data_obs");
87	>	TH1F *signal = (TH1F*)f->Get("signal");
88	>	TH1F *Tbackground = (TH1F*)f->Get("TTbarBackground");
89	>	TH1F *Zbackground  = (TH1F*)f->Get("ZJetsBackground");
90	>
91	>	assert(dataob);
92	>	assert(signal);
93	>	assert(Tbackground);
94	>	assert(Zbackground);
95
96	<	ofstream datacard;
97	<	write_warning(__FUNCTION__,"Need to rethink systematics we want to use !");
98	<	datacard.open ((RunDirectory+"/susydatacard.txt").c_str());
99	<	datacard << "Writing this to a file.\n";
100	<	datacard << "imax 1\n"; // number of channels
101	<	datacard << "jmax 2\n"; // number of backgrounds (Z+Jets prediction and TTbar prediction)
102	<	datacard << "kmax *\n"; // number of nuisance parameters (sources of systematic uncertainties)
103	<	datacard << "---------------\n";
104	<	datacard << "shapes * * limitfile.root $PROCESS $PROCESS_$SYSTEMATIC\n";
105	<	datacard << "---------------\n";
106	<	datacard << "bin 1\n";
107	<	datacard << "observation "<<dataob->Integral()<<"\n";
108	<	datacard << "------------------------------\n";
109	<	datacard << "bin             1          1          1\n";
110	<	datacard << "process         signal     TTbarBackground     ZJetsBackground\n";
111	<	datacard << "process         0          1          2\n";
112	<	datacard << "rate            "<<signal->Integral()<<"         "<<background1->Integral()<<"         "<<background2->Integral()<<"\n";
113	<	datacard << "--------------------------------\n";
114	<	datacard << "lumi     lnN    " << 1+ PlottingSetup::lumiuncert << "       -       -    luminosity uncertainty\n"; // only affects MC -> signal!
115	<	datacard << "Stat   shape    -          1          -    statistical uncertainty (ttbar)\n";
116	<	datacard << "Stat   shape    -          -          1    statistical uncertainty (zjets)\n";
117	<	datacard << "Sys    shape    -          1          -    systematic uncertainty on ttbar\n";
118	<	datacard << "Sys    shape    -          -          1    systematic uncertainty on zjets\n";
119	<	datacard << "Stat   shape    1          -          -    statistical uncertainty (signal)\n";
120	<	datacard << "JES    shape    1          -          -    uncertainty on Jet Energy Scale\n";
121	<	datacard << "JSU      lnN    " << 1+uJSU << "          -          -    JZB Scale Uncertainty\n";
122	<	if(uPDF>0) datacard << "PDF      lnN    " << 1+uPDF << "          -          -    uncertainty from PDFs\n";
96	>
97	>	ofstream datacard;
98	>	write_warning(__FUNCTION__,"Need to rethink systematics we want to use !");
99	>	datacard.open ((RunDirectory+"/susydatacard.txt").c_str());
100	>	datacard << "Writing this to a file.\n";
101	>	datacard << "imax " << dataob->GetNbinsX() << "\n"; // number of channels
102	>	datacard << "jmax 2\n"; // number of backgrounds (Z+Jets prediction and TTbar prediction)
103	>	datacard << "kmax *\n"; // number of nuisance parameters (sources of systematic uncertainties)
104	>	datacard << "---------------\n";
105	>	datacard << "bin\t";
106	>	for(int i=1;i<=dataob->GetNbinsX();i++) datacard << " " << dataob->GetBinLowEdge(i) << "to" << dataob->GetBinLowEdge(i)+dataob->GetBinWidth(i) << " \t";
107	>	datacard << "\n";
108	>	datacard << "observation\t";
109	>	for(int i=1;i<=dataob->GetNbinsX();i++) datacard << " " << dataob->GetBinContent(i) << " \t";
110	>	datacard<<"\n";
111	>	datacard << "------------------------------\n";
112	>	datacard << "bin\t";
113	>	for(int i=1;i<=dataob->GetNbinsX();i++) {
114	>	datacard << " " << dataob->GetBinLowEdge(i) << "to" << dataob->GetBinLowEdge(i)+dataob->GetBinWidth(i) << "\t" <<
115	>	" " << dataob->GetBinLowEdge(i) << "to" << dataob->GetBinLowEdge(i)+dataob->GetBinWidth(i) << "\t" <<
116	>	" " << dataob->GetBinLowEdge(i) << "to" << dataob->GetBinLowEdge(i)+dataob->GetBinWidth(i) << "\t";
117	>	}
118	>	datacard << "\n";
119	>	datacard << "process\t";
120	>	for(int i=1;i<=dataob->GetNbinsX();i++) datacard << "\t signal \t TTbarBackground \t ZJetsBackground";
121	>	datacard << "\n";
122	>	datacard << "process\t";
123	>	for(int i=1;i<=dataob->GetNbinsX();i++) datacard << "\t 0 \t 1 \t 2";
124	>	datacard << "\n";
125	>
126	>
127	>	datacard << "rate\t";
128	>	for(int i=1;i<=dataob->GetNbinsX();i++) datacard << "\t " << signal->GetBinContent(i) << " \t " << Tbackground->GetBinContent(i) << " \t " << Zbackground->GetBinContent(i) << "\t";
129	>	datacard<<"\n";
130	>
131	>	datacard << "--------------------------------\n";
132	>
133	>
134	>	datacard << "lumi     lnN    \t";
135	>	for(int i=1;i<=dataob->GetNbinsX();i++) datacard << " " << 1+ PlottingSetup::lumiuncert << "\t - \t -";
136	>	datacard << "       luminosity uncertainty\n"; // only affects MC -> signal!
137	>
138	>	// Statistical uncertainty
139	>	datacard << "Stat     lnN    \t";
140	>	for(int i=1;i<=dataob->GetNbinsX();i++) {
141	>	//Signal
142	>	float central = signal->GetBinContent(i);
143	>	float up      = ((TH1F*)f->Get("signal_StatDown"))->GetBinContent(i);
144	>	float down    = ((TH1F*)f->Get("signal_StatUp"))->GetBinContent(i);
145	>	float suncert=0;
146	>	if(central>0) {
147	>	if(abs(up-central)>abs(down-central)) suncert=abs(up-central)/central;
148	>	else suncert=abs(central-down)/central;
149	>	}
150	>
151	>	//TTbar
152	>	central = Tbackground->GetBinContent(i);
153	>	up      = ((TH1F*)f->Get("TTbarBackground_StatUp"))->GetBinContent(i);
154	>	down    = ((TH1F*)f->Get("TTbarBackground_StatDown"))->GetBinContent(i);
155	>	float tuncert=0;
156	>	if(central>0) {
157	>	if(abs(up-central)>abs(down-central)) tuncert=abs(up-central)/central;
158	>	else tuncert=abs(central-down)/central;
159	>	}
160	>	//ZJets
161	>	central = Zbackground->GetBinContent(i);
162	>	up      = ((TH1F*)f->Get("ZJetsBackground_StatUp"))->GetBinContent(i);
163	>	down    = ((TH1F*)f->Get("ZJetsBackground_StatDown"))->GetBinContent(i);
164	>	float zuncert=0;
165	>	if(central>0) {
166	>	if(abs(up-central)>abs(down-central)) zuncert=abs(up-central)/central;
167	>	else zuncert=abs(central-down)/central;
168	>	}
169	>	datacard << " " << 1+suncert << " \t " << 1+tuncert << "\t " << 1+zuncert << " \t ";
170	>	}
171	>
172	>	datacard << "       statistical uncertainty\n";
173	>
174	>	// Statistical uncertainty
175	>	datacard << "Sys     lnN    \t";
176	>	for(int i=1;i<=dataob->GetNbinsX();i++) {
177	>	float central = Tbackground->GetBinContent(i);
178	>	float up      = ((TH1F*)f->Get("TTbarBackground_SysUp"))->GetBinContent(i);
179	>	float down    = ((TH1F*)f->Get("TTbarBackground_SysDown"))->GetBinContent(i);
180	>	float tuncert=0;
181	>	if(central>0) {
182	>	if(abs(up-central)>abs(down-central)) tuncert=abs(up-central)/central;
183	>	else tuncert=abs(central-down)/central;
184	>	}
185	>	central = Zbackground->GetBinContent(i);
186	>	up      = ((TH1F*)f->Get("ZJetsBackground_SysUp"))->GetBinContent(i);
187	>	down    = ((TH1F*)f->Get("ZJetsBackground_SysDown"))->GetBinContent(i);
188	>	float zuncert=0;
189	>	if(central>0) {
190	>	if(abs(up-central)>abs(down-central)) zuncert=abs(up-central)/central;
191	>	else zuncert=abs(central-down)/central;
192	>	}
193	>	datacard << " - \t " << 1+tuncert << "\t " << 1+zuncert << " \t ";
194	>	}
195	>
196	>	datacard << "       systematic uncertainty\n";
197	>
198	>
199	>	float JESup = ((TH1F*)f->Get("signal_JESUp"))->Integral();
200	>	float JESdn = ((TH1F*)f->Get("signal_JESDown"))->Integral();
201	>	float central = signal->Integral();
202	>	float uJES=0;
203	>	if(abs(JESup-central)>abs(JESdn-central)) uJES=abs(JESup-central)/central;
204	>	else uJES=abs(central-JESdn)/central;
205	>	datacard << "JES    lnN";
206	>	for(int i=1;i<=dataob->GetNbinsX();i++) datacard << " " << 1+uJES << "\t - \t  - \t";
207	>	datacard << "uncertainty on Jet Energy Scale\n";
208	>
209	>	datacard << "JSU      lnN    ";
210	>	for(int i=1;i<=dataob->GetNbinsX();i++) datacard << " " << 1+uJSU << "\t - \t  - \t";
211	>	datacard << "JZB Scale Uncertainty\n";
212	>
213	>	if(uPDF>0) {
214	>	datacard << "PDF      lnN    ";
215	>	for(int i=1;i<=dataob->GetNbinsX();i++) datacard << " " << 1+uPDF << "\t - \t - \t";
216	>	datacard << "uncertainty from PDFs\n";
217	>	}
218	>
219	>	datacard.close();
220	>	}
221	>
222	>	void prepare_limit_datacard(string RunDirectory, TFile *f, float uJSU, float uPDF) {
223
224	<	// datacard << "peak   shape    1          1    uncertainty on signal resolution.
225	<	datacard.close();
224	>	if(FullMCAnalysis) {
225	>	//dealing with MC analysis - need to use shapes.
226	>	prepare_MC_datacard(RunDirectory,f,uJSU,uPDF);
227	>	} else {
228	>	//doing mutibin analysis
229	>	prepare_datadriven_datacard(RunDirectory,f,uJSU,uPDF);
230	>	}
231	>
232		}
233
234		float QuickDrawNevents=0;
#	Line 122 | Line 261 | TH1F* QuickDraw(TTree *events, string hn
261		return histo;
262		}
263
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	–	}*/
264
265		void SQRT(TH1F *h) {
266		for (int i=1;i<=h->GetNbinsX();i++) {
#	Line 147 | Line 276 | TH1F* Multiply(TH1F *h1, TH1F *h2) {
276		return h;
277		}
278
279	+
280	+	string ReduceNumericHistoName(string origname) {
281	+	int pos = origname.find("__h_");
282	+	if(pos == -1) return origname;
283	+	return origname.substr(0,pos);
284	+	}
285	+
286	+
287	+
288	+	void generate_mc_shapes(bool dotree, TTree *signalevents, string mcjzb,string datajzb,vector<float> binning, TCut weight, TCut JetCut, string addcut, string identifier, map <  pair<float, float>, map<string, float>  > &xsec) {
289	+	//weight can be varied to take care of efficiency errors (weightEffUp, weightEffDown)
290	+	vector<float> mbinning;
291	+	mbinning.push_back(-8000);
292	+	for(int i=binning.size()-1;i>=0;i--) mbinning.push_back(-binning[i]);
293	+	mbinning.push_back(0);
294	+	for(int i=0;i<binning.size();i++) mbinning.push_back(binning[i]);
295	+	mbinning.push_back(8000);
296	+
297	+
298	+	TCanvas *shapecan = new TCanvas("shapecan","shapecan");
299	+	TH1F* Signal;
300	+
301	+	stringstream sInverseCutWeight;
302	+	sInverseCutWeight << "(1.0/" << (const char*)cutWeight<<")";
303	+	TCut InverseCutWeight(sInverseCutWeight.str().c_str());
304	+	if(weight==cutWeight) InverseCutWeight = TCut("1.0");
305	+	if(!dotree) {
306	+	//dealing with ALL MC samples (when we save the standard file)
307	+	THStack McPredicted = allsamples.DrawStack("McPred",mcjzb,mbinning, "JZB [GeV]", "events", ((cutmass&&cutOSSF&&JetCut)*(weight*InverseCutWeight)),mc,luminosity);
308	+	int nhists = McPredicted.GetHists()->GetEntries();
309	+	for (int i = 0; i < nhists; i++) {
310	+	TH1* hist = static_cast<TH1*>(McPredicted.GetHists()->At(i));
311	+	//          cout << hist->GetName() << " : " << hist->Integral() << endl;
312	+	//          cout << "     " << ReduceNumericHistoName(hist->GetName()) << endl;
313	+	if(identifier=="StatUp") {
314	+	float appliedweight = hist->Integral() / hist->GetEntries();
315	+	for(int ibin=1;ibin<=hist->GetNbinsX();ibin++) hist->SetBinContent(ibin,hist->GetBinContent(ibin)+appliedweight*TMath::Sqrt(hist->GetBinContent(ibin)/appliedweight));
316	+	}
317	+	if(identifier=="StatDown") {
318	+	float appliedweight = hist->Integral() / hist->GetEntries();
319	+	for(int ibin=1;ibin<=hist->GetNbinsX();ibin++) hist->SetBinContent(ibin,hist->GetBinContent(ibin)-appliedweight*TMath::Sqrt(hist->GetBinContent(ibin)/appliedweight));
320	+	}
321	+	hist->SetName(("mc_"+ReduceNumericHistoName(hist->GetName())+"_"+identifier).c_str());
322	+	hist->Write();
323	+	}
324	+	} else {
325	+	string histoname="mc_signal";
326	+	if(identifier!="") histoname=histoname+"_"+identifier;
327	+	Signal = QuickDraw(signalevents,histoname,mcjzb,binning,"JZB", "events",((cutmass&&cutOSSF&&JetCut&&basiccut))*(weight*InverseCutWeight),addcut,mc,luminosity,xsec);
328	+	if(identifier=="StatUp") {
329	+	float appliedweight = Signal->Integral() / Signal->GetEntries();
330	+	for(int ibin=1;ibin<=Signal->GetNbinsX();ibin++) Signal->SetBinContent(ibin,Signal->GetBinContent(ibin)+appliedweight*TMath::Sqrt(Signal->GetBinContent(ibin)/appliedweight));
331	+	}
332	+	if(identifier=="StatDown") {
333	+	float appliedweight = Signal->Integral() / Signal->GetEntries();
334	+	for(int ibin=1;ibin<=Signal->GetNbinsX();ibin++) Signal->SetBinContent(ibin,Signal->GetBinContent(ibin)-appliedweight*TMath::Sqrt(Signal->GetBinContent(ibin)/appliedweight));
335	+	}
336	+	Signal->Write();
337	+	}
338	+
339	+	/*
340	+	*
341	+	* Jet energy scale (easy, use different JetCut)
342	+	* JZB Scale Uncertainty (will be a number - signal only)
343	+	* Efficiency (will be weightEffUp, weightEffDown)
344	+	* PDFs (signal only) -> Will be a number yet again, computed in the traditional way
345	+	*/
346	+	delete shapecan;
347	+
348	+	}
349	+
350	+
351		void generate_shapes_for_systematic(bool signalonly, bool dataonly,TFile *limfile, TTree *signalevents, string identifier, string mcjzb, string datajzb, int JES,vector<float> binning, TCanvas *limcan, string addcut, map <  pair<float, float>, map<string, float>  > &xsec) {
352
353		binning.push_back(8000);
#	Line 176 | Line 377 | void generate_shapes_for_systematic(bool
377		}
378
379		if(signalonly) {
380	<	cout << "Processing a signal with mcjzb: " << mcjzb << " (identifier: " << identifier << ")" << endl;
380	>	dout << "Processing a signal with mcjzb: " << mcjzb << " (identifier: '" << identifier << "')" << endl;
381		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);
382		TH1F *ZOSSFN = QuickDraw(signalevents,"ZOSSFN","-"+mcjzb,binning, "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
383	<	TH1F *ZOSOFN = QuickDraw(signalevents,"ZOSOFN","-"+mcjzb,binning, "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
383	>	TH1F *ZOSOFP;
384	>	TH1F *ZOSOFN;
385	>
386	>	if(!PlottingSetup::FullMCAnalysis) {
387	>	ZOSOFP = QuickDraw(signalevents,"ZOSOFP",mcjzb,binning,     "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
388	>	ZOSOFN = QuickDraw(signalevents,"ZOSOFN","-"+mcjzb,binning, "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
389	>	}
390
391		TH1F *SBOSSFP;
392		TH1F *SBOSOFP;
393		TH1F *SBOSSFN;
394		TH1F *SBOSOFN;
395
396	<	if(PlottingSetup::RestrictToMassPeak) {
396	>	if(PlottingSetup::RestrictToMassPeak&&!PlottingSetup::FullMCAnalysis) {
397		SBOSSFP = QuickDraw(signalevents,"SBOSSFP",mcjzb,binning, "JZB", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
398		SBOSOFP = QuickDraw(signalevents,"SBOSOFP",mcjzb,binning, "JZB", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
399		SBOSSFN = QuickDraw(signalevents,"SBOSSFN","-"+mcjzb,binning, "JZB", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
#	Line 206 | Line 412 | void generate_shapes_for_systematic(bool
412		TH1F *flippedLpred = new TH1F("flippedLpred","flippedLpred",binning.size()-1,&binning[0]);
413
414		Lobs->Add(ZOSSFP);
415	<	Lpred->Add(ZOSSFN);
415	>	if(!PlottingSetup::FullMCAnalysis) Lpred->Add(ZOSSFN);
416
417	<	cout << "SITUATION FOR SIGNAL: " << endl;
418	<	cout << " OSSF     JZB> 0 : " << ZOSSFP->Integral() << "     JZB < 0 :" << ZOSSFN->Integral() << endl;
419	<	cout << " OSOF     JZB> 0 : " << ZOSOFP->Integral() << "     JZB < 0 :" << ZOSOFN->Integral() << endl;
420	<	if(PlottingSetup::RestrictToMassPeak) {
421	<	cout << " OSSF SB  JZB> 0 : " << SBOSSFP->Integral() << "     JZB < 0 :" << SBOSSFN->Integral() << endl;
422	<	cout << " OSOF SB  JZB> 0 : " << SBOSOFP->Integral() << "     JZB < 0 :" << SBOSOFN->Integral() << endl;
417	>	dout << "SITUATION FOR SIGNAL: " << endl;
418	>
419	>
420	>	if(PlottingSetup::FullMCAnalysis) {
421	>	dout << " OSSF     JZB> 0 : " << ZOSSFP->Integral() << endl;
422	>	} else {
423	>	dout << " OSSF     JZB> 0 : " << ZOSSFP->Integral() << "     JZB < 0 :" << ZOSSFN->Integral() << endl;
424	>	dout << " OSOF     JZB> 0 : " << ZOSOFP->Integral() << "     JZB < 0 :" << ZOSOFN->Integral() << endl;
425	>	if(PlottingSetup::RestrictToMassPeak&&!PlottingSetup::FullMCAnalysis) {
426	>	dout << " OSSF SB  JZB> 0 : " << SBOSSFP->Integral() << "     JZB < 0 :" << SBOSSFN->Integral() << endl;
427	>	dout << " OSOF SB  JZB> 0 : " << SBOSOFP->Integral() << "     JZB < 0 :" << SBOSOFN->Integral() << endl;
428	>	}
429		}
430
431
432		flippedLobs->Add(ZOSSFN);
433	<	flippedLpred->Add(ZOSSFP);
433	>	if(!PlottingSetup::FullMCAnalysis) flippedLpred->Add(ZOSSFP);
434
435	<	if(PlottingSetup::RestrictToMassPeak) {
436	<	Lpred->Add(ZOSOFP,1.0/3);
437	<	Lpred->Add(ZOSOFN,-1.0/3);
438	<	Lpred->Add(SBOSSFP,1.0/3);
439	<	Lpred->Add(SBOSSFN,-1.0/3);
440	<	Lpred->Add(SBOSOFP,1.0/3);
441	<	Lpred->Add(SBOSOFN,-1.0/3);
442	<
443	<	//flipped prediction
444	<	flippedLpred->Add(ZOSOFP,-1.0/3);
445	<	flippedLpred->Add(ZOSOFN,1.0/3);
446	<	flippedLpred->Add(SBOSSFP,-1.0/3);
447	<	flippedLpred->Add(SBOSSFN,1.0/3);
448	<	flippedLpred->Add(SBOSOFP,-1.0/3);
449	<	flippedLpred->Add(SBOSOFN,1.0/3);
450	<
451	<	} else {
452	<	Lpred->Add(ZOSOFP,1.0);
453	<	Lpred->Add(ZOSOFN,-1.0);
454	<
455	<	//flipped prediction
456	<	flippedLpred->Add(ZOSOFP,-1.0);
457	<	flippedLpred->Add(ZOSOFN,1.0);
435	>	if(!PlottingSetup::FullMCAnalysis) {
436	>	if(PlottingSetup::RestrictToMassPeak) {
437	>	Lpred->Add(ZOSOFP,1.0/3);
438	>	Lpred->Add(ZOSOFN,-1.0/3);
439	>	Lpred->Add(SBOSSFP,1.0/3);
440	>	Lpred->Add(SBOSSFN,-1.0/3);
441	>	Lpred->Add(SBOSOFP,1.0/3);
442	>	Lpred->Add(SBOSOFN,-1.0/3);
443	>
444	>	//flipped prediction
445	>	flippedLpred->Add(ZOSOFP,-1.0/3);
446	>	flippedLpred->Add(ZOSOFN,1.0/3);
447	>	flippedLpred->Add(SBOSSFP,-1.0/3);
448	>	flippedLpred->Add(SBOSSFN,1.0/3);
449	>	flippedLpred->Add(SBOSOFP,-1.0/3);
450	>	flippedLpred->Add(SBOSOFN,1.0/3);
451	>	} else {
452	>	Lpred->Add(ZOSOFP,1.0);
453	>	Lpred->Add(ZOSOFN,-1.0);
454	>
455	>	//flipped prediction
456	>	flippedLpred->Add(ZOSOFP,-1.0);
457	>	flippedLpred->Add(ZOSOFN,1.0);
458	>	}
459		}
460
461		TH1F *signal = (TH1F*)Lobs->Clone("signal");
462	<	signal->Add(Lpred,-1);
462	>	if(!PlottingSetup::FullMCAnalysis) signal->Add(Lpred,-1);
463		signal->SetName(signalname.c_str());
464		signal->SetTitle(signalname.c_str());
465		signal->Write();
466
467		TH1F *flippedsignal = (TH1F*)flippedLobs->Clone();
468	<	flippedsignal->Add(flippedLpred,-1);
468	>	if(!PlottingSetup::FullMCAnalysis) flippedsignal->Add(flippedLpred,-1);
469		flippedsignal->SetName(("flipped_"+signalname).c_str());
470		flippedsignal->Write();
471
#	Line 266 | Line 479 | void generate_shapes_for_systematic(bool
479		signalStatDn->SetTitle(((string)signal->GetTitle()+"_StatDown").c_str());
480
481		for(int i=1;i<=signalStatDn->GetNbinsX();i++) {
482	<	float staterr = TMath::Sqrt(Lpred->GetBinContent(i) + Lobs->GetBinContent(i));
483	<	cout << "Stat err in bin " << i << " : " << staterr << endl;
484	<	cout << "    prediction: " << Lpred->GetBinContent(i) << " , observation: " << Lobs->GetBinContent(i) << " --> signal: " << signal->GetBinContent(i) << endl;
485	<	cout << "    we obtain : " << signal->GetBinContent(i)-staterr << " , " << signal->GetBinContent(i)+staterr << endl;
486	<	signalStatDn->SetBinContent(i,signal->GetBinContent(i)-staterr);
482	>	float staterr;
483	>	if(!PlottingSetup::FullMCAnalysis) staterr = TMath::Sqrt(Lpred->GetBinContent(i) + Lobs->GetBinContent(i));
484	>	else staterr = TMath::Sqrt(Lobs->GetBinContent(i));
485	>
486	>	if(!PlottingSetup::FullMCAnalysis) {
487	>	dout << "Stat err in bin " << i << " : " << staterr << endl;
488	>	dout << "    prediction: " << Lpred->GetBinContent(i) << " , observation: " << Lobs->GetBinContent(i) << " --> signal: " << signal->GetBinContent(i) << endl;
489	>	dout << "    we obtain : " << signal->GetBinContent(i)-staterr << " , " << signal->GetBinContent(i)+staterr << endl;
490	>	}
491	>	if(signal->GetBinContent(i)-staterr>0) signalStatDn->SetBinContent(i,signal->GetBinContent(i)-staterr);
492	>	else signalStatDn->SetBinContent(i,0);
493		signalStatUp->SetBinContent(i,signal->GetBinContent(i)+staterr);
494		signal->SetBinError(i,staterr);
495		}
#	Line 290 | Line 509 | void generate_shapes_for_systematic(bool
509		flippedsignalStatDn->SetTitle(((string)flippedsignal->GetTitle()+"_StatDown").c_str());
510
511		for(int i=1;i<=flippedsignalStatDn->GetNbinsX();i++) {
512	<	float staterr = TMath::Sqrt(flippedLpred->GetBinContent(i) + flippedLobs->GetBinContent(i));
513	<	flippedsignalStatDn->SetBinContent(i,flippedsignal->GetBinContent(i)-staterr);
512	>	float staterr;
513	>	if(!PlottingSetup::FullMCAnalysis) staterr = TMath::Sqrt(flippedLpred->GetBinContent(i) + flippedLobs->GetBinContent(i));
514	>	else staterr = TMath::Sqrt(flippedLobs->GetBinContent(i));
515	>	if(flippedsignal->GetBinContent(i)-staterr>0) flippedsignalStatDn->SetBinContent(i,flippedsignal->GetBinContent(i)-staterr);
516	>	else flippedsignalStatDn->SetBinContent(i,0);
517		flippedsignalStatUp->SetBinContent(i,flippedsignal->GetBinContent(i)+staterr);
518		flippedsignal->SetBinError(i,staterr);
519		}
#	Line 323 | Line 545 | void generate_shapes_for_systematic(bool
545
546
547		if(dataonly) {
548	<	cout << "Processing data with datajzb: " << datajzb << endl;
548	>	dout << "Processing data with datajzb: " << datajzb << endl;
549		TH1F *ZOSSFP = allsamples.Draw("ZOSSFP",datajzb,binning, "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,data,luminosity);
550		TH1F *ZOSOFP = allsamples.Draw("ZOSOFP",datajzb,binning, "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,data,luminosity);
551		TH1F *ZOSSFN = allsamples.Draw("ZOSSFN","-"+datajzb,binning, "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,data,luminosity);
#	Line 457 | Line 679 | void generate_shapes_for_systematic(bool
679		SQRT(predstaterr);
680		TH1F *bgStatUp = (TH1F*)pred->Clone("background_StatUp");
681		bgStatUp->Add(predstaterr);
682	+	EliminateNegativeEntries(bgStatUp);
683		bgStatUp->Write();
684		TH1F *bgStatDn = (TH1F*)pred->Clone("background_StatDown");
685		bgStatDn->Add(predstaterr,-1);
686	+	EliminateNegativeEntries(bgStatDn);
687		bgStatDn->Write();
688		//      delete bgStatDn;
689		//      delete bgStatUp;
#	Line 476 | Line 700 | void generate_shapes_for_systematic(bool
700		SQRT(flippedpredstaterr);
701		TH1F *fbgStatUp = (TH1F*)flippedpred->Clone("flipped_background_StatUp");
702		fbgStatUp->Add(predstaterr);
703	+	EliminateNegativeEntries(fbgStatUp);
704		fbgStatUp->Write();
705		TH1F *fbgStatDn = (TH1F*)flippedpred->Clone("flipped_background_StatDown");
706		fbgStatDn->Add(predstaterr,-1);
707	+	EliminateNegativeEntries(fbgStatDn);
708		fbgStatDn->Write();
709		//      delete fbgStatDn;
710		//      delete fbgStatUp;
#	Line 491 | Line 717 | void generate_shapes_for_systematic(bool
717		SQRT(Tpredstaterr);
718		TH1F *TpredStatUp = (TH1F*)Tpred->Clone("TTbarBackground_StatUp");
719		TpredStatUp->Add(Tpredstaterr);
720	+	EliminateNegativeEntries(TpredStatUp);
721		TpredStatUp->Write();
722		TH1F *TpredStatDn = (TH1F*)Tpred->Clone("TTbarBackground_StatDown");
723		TpredStatDn->Add(Tpredstaterr,-1);
724	+	EliminateNegativeEntries(TpredStatDn);
725		TpredStatDn->Write();
726		//      delete TpredStatDn;
727		//      delete TpredStatUp;
#	Line 506 | Line 734 | void generate_shapes_for_systematic(bool
734		SQRT(fTpredstaterr);
735		TH1F *fTpredStatUp = (TH1F*)flippedTpred->Clone("flipped_TTbarBackground_StatUp");
736		fTpredStatUp->Add(fTpredstaterr);
737	+	EliminateNegativeEntries(fTpredStatUp);
738		fTpredStatUp->Write();
739		TH1F *fTpredStatDn = (TH1F*)flippedTpred->Clone("flipped_TTbarBackground_StatDown");
740		fTpredStatDn->Add(fTpredstaterr,-1);
741	+	EliminateNegativeEntries(fTpredStatDn);
742		fTpredStatDn->Write();
743		//      delete fTpredStatDn;
744		//      delete fTpredStatUp;
#	Line 522 | Line 752 | void generate_shapes_for_systematic(bool
752		SQRT(Zpredstaterr);
753		TH1F *ZpredStatUp = (TH1F*)Zpred->Clone("ZJetsBackground_StatUp");
754		ZpredStatUp->Add(Zpredstaterr);
755	+	EliminateNegativeEntries(ZpredStatUp);
756		ZpredStatUp->Write();
757		TH1F *ZpredStatDn = (TH1F*)Zpred->Clone("ZJetsBackground_StatDown");
758		ZpredStatDn->Add(Zpredstaterr,-1);
759	+	EliminateNegativeEntries(ZpredStatDn);
760		ZpredStatDn->Write();
761		//      delete ZpredStatDn;
762		//      delete ZpredStatUp;
#	Line 537 | Line 769 | void generate_shapes_for_systematic(bool
769		SQRT(fTpredstaterr);
770		TH1F *fZpredStatUp = (TH1F*)flippedZpred->Clone("flipped_ZJetsBackground_StatUp");
771		fZpredStatUp->Add(fZpredstaterr);
772	+	EliminateNegativeEntries(fZpredStatUp);
773		fZpredStatUp->Write();
774		TH1F *fZpredStatDn = (TH1F*)flippedZpred->Clone("flipped_ZJetsBackground_StatDown");
775		fZpredStatDn->Add(fZpredstaterr,-1);
776	+	EliminateNegativeEntries(fZpredStatDn);
777		fZpredStatDn->Write();
778		//      delete fZpredStatDn;
779		//      delete fZpredStatUp;
#	Line 559 | Line 793 | void generate_shapes_for_systematic(bool
793		SQRT(predsyserr);
794		TH1F *bgSysUp = (TH1F*)pred->Clone("background_SysUp");
795		bgSysUp->Add(predsyserr);
796	+	EliminateNegativeEntries(bgSysUp);
797		bgSysUp->Write();
798		TH1F *bgSysDn = (TH1F*)pred->Clone("background_SysDown");
799		bgSysDn->Add(predsyserr,-1);
800	+	EliminateNegativeEntries(bgSysDn);
801		bgSysDn->Write();
802		delete predsyserr;
803
#	Line 577 | Line 813 | void generate_shapes_for_systematic(bool
813		SQRT(fpredsyserr);
814		TH1F *fbgSysUp = (TH1F*)flippedpred->Clone("flipped_background_SysUp");
815		fbgSysUp->Add(fpredsyserr);
816	+	EliminateNegativeEntries(fbgSysUp);
817		fbgSysUp->Write();
818		TH1F *fbgSysDn = (TH1F*)flippedpred->Clone("flipped_background_SysDown");
819		fbgSysDn->Add(fpredsyserr,-1);
820	+	EliminateNegativeEntries(fbgSysDn);
821		fbgSysDn->Write();
822		delete fpredsyserr;
823
#	Line 592 | Line 830 | void generate_shapes_for_systematic(bool
830		SQRT(Tpredsyserr);
831		TH1F *TpredSysUp = (TH1F*)Tpred->Clone("TTbarBackground_SysUp");
832		TpredSysUp->Add(Tpredsyserr);
833	+	EliminateNegativeEntries(TpredSysUp);
834		TpredSysUp->Write();
835		TH1F *TpredSysDn = (TH1F*)Tpred->Clone("TTbarBackground_SysDown");
836		TpredSysDn->Add(Tpredsyserr,-1);
837	+	EliminateNegativeEntries(TpredSysDn);
838		TpredSysDn->Write();
839		delete Tpredsyserr;
840
#	Line 606 | Line 846 | void generate_shapes_for_systematic(bool
846		SQRT(fTpredsyserr);
847		TH1F *fTpredSysUp = (TH1F*)flippedTpred->Clone("flipped_TTbarBackground_SysUp");
848		fTpredSysUp->Add(fTpredsyserr);
849	+	EliminateNegativeEntries(fTpredSysUp);
850		fTpredSysUp->Write();
851		TH1F *fTpredSysDn = (TH1F*)flippedTpred->Clone("flipped_TTbarBackground_SysDown");
852		fTpredSysDn->Add(fTpredsyserr,-1);
853	+	EliminateNegativeEntries(fTpredSysDn);
854		fTpredSysDn->Write();
855		delete fTpredsyserr;
856
#	Line 623 | Line 865 | void generate_shapes_for_systematic(bool
865		SQRT(Zpredsyserr);
866		TH1F *ZpredSysUp = (TH1F*)Zpred->Clone("ZJetsBackground_SysUp");
867		ZpredSysUp->Add(Zpredsyserr);
868	+	EliminateNegativeEntries(ZpredSysUp);
869		ZpredSysUp->Write();
870		TH1F *ZpredSysDn = (TH1F*)Zpred->Clone("ZJetsBackground_SysDown");
871		ZpredSysDn->Add(Zpredsyserr,-1);
872	+	EliminateNegativeEntries(ZpredSysDn);
873		ZpredSysDn->Write();
874		delete Zpredsyserr;
875
#	Line 638 | Line 882 | void generate_shapes_for_systematic(bool
882		SQRT(fZpredsyserr);
883		TH1F *fZpredSysUp = (TH1F*)flippedZpred->Clone("flipped_ZJetsBackground_SysUp");
884		fZpredSysUp->Add(fZpredsyserr);
885	+	EliminateNegativeEntries(fZpredSysUp);
886		fZpredSysUp->Write();
887		TH1F *fZpredSysDn = (TH1F*)flippedZpred->Clone("flipped_ZJetsBackground_SysDown");
888		fZpredSysDn->Add(fZpredsyserr,-1);
889	+	EliminateNegativeEntries(fZpredSysDn);
890		fZpredSysDn->Write();
891		delete fZpredsyserr;
892		}
893
894		/*if(identifier=="") {
895		for(int i=0;i<binning.size()-1;i++) {
896	<	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;
896	>	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;
897		}
898		}*/
899		delete ZOSSFP;
#	Line 679 | Line 925 | ShapeDroplet LimitsFromShapes(bool asymp
925
926		ensure_directory_exists(RunDirectory);
927
928	<	TFile *datafile = new TFile("../StoredShapes.root","READ");
928	>	//  write_warning(__FUNCTION__,"Modified the shape output file! PLEASE RESTORE!!!");  TFile *datafile = new TFile((PlottingSetup::cbafbasedir+"/DistributedModelCalculations/TEST___StoredShapes.root").c_str(),"READ");
929	>	TFile *datafile = new TFile((PlottingSetup::cbafbasedir+"/DistributedModelCalculations/StoredShapes.root").c_str(),"READ");
930		if(datafile->IsZombie()) {
931		write_error(__FUNCTION__,"Fatal error: The stored shapes are not available!");
932		assert(!datafile->IsZombie());
933		}
934	<	cout << "Run Directory: " << RunDirectory << endl;
934	>	dout << "Run Directory: " << RunDirectory << endl;
935		TFile *limfile = new TFile((RunDirectory+"/PRElimitfile.root").c_str(),"RECREATE");
936
937		TIter nextkey(datafile->GetListOfKeys());
#	Line 702 | Line 949 | ShapeDroplet LimitsFromShapes(bool asymp
949		bool dataonly=false;
950
951		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);
952		generate_shapes_for_systematic(signalonly,dataonly,limfile,events,"JESUp",mcjzb,datajzb,JESup,jzbbins,limcan,addcut,xsec);
953		generate_shapes_for_systematic(signalonly,dataonly,limfile,events,"JESDown",mcjzb,datajzb,JESdown,jzbbins,limcan,addcut,xsec);
954
#	Line 717 | Line 962 | ShapeDroplet LimitsFromShapes(bool asymp
962		bool docomplicatedmSUGRAxsreweighting=false; //if you modify this value please also adapt it in Systematics.C not only here in ShapeLimit.C
963
964		MCefficiency(events,mceff,mcefferr,flipped,mcjzb,requireZ,Neventsinfile,SUSYScanSpace::SUSYscantype,xsec,addcut,-1);
965	<	if(mceff<0) flipped=1;
965	>	if(mceff<0) {
966	>	flipped=1;
967	>	write_info(__FUNCTION__,"Doing flipping!");
968	>	}
969		doJZBscale(events,SUSYScanSpace::SUSYscantype==mSUGRA&&docomplicatedmSUGRAxsreweighting,xsec,scaledown,scaleup,scalesyst,JZBScaleUncert,informalname,flipped,requireZ,addcut);
970		float PDFuncert=0;
971		int NPdfs = get_npdfs(events);
#	Line 763 | Line 1011 | ShapeDroplet LimitsFromShapes(bool asymp
1011		obj->Write();
1012		}
1013
1014	<	prepare_limit_datacard(RunDirectory,final_limfile,JZBscale,PDFuncert,flipped);
1014	>	prepare_limit_datacard(RunDirectory,final_limfile,JZBscale,PDFuncert);
1015
1016		final_limfile->Close();
1017		limfile->Close();
1018		dout << "Info: Shape root file and datacard have been generated in " << RunDirectory << endl;
1019		stringstream command;
1020		if(asymptotic) {
1021	<	if(firstGuess>0) command << "bash CreateModel.sh " << RunDirectory << " susydatacard.txt" << " 1 0 0 " << firstGuess; // ASYMPTOTIC LIMITS WITH FIRST GUESS
1022	<	else command << "bash CreateModel.sh " << RunDirectory << " susydatacard.txt" << " 0 0 0 0"; // ASYMPTOTIC LIMITS
1021	>	if(firstGuess>0) command << "bash " << PlottingSetup::cbafbasedir<< "/DistributedModelCalculations/ShapeLimits/CreateModel.sh " << RunDirectory << " susydatacard.txt" << " 1 0 0 " << firstGuess; // ASYMPTOTIC LIMITS WITH FIRST GUESS
1022	>	else command << "bash " << PlottingSetup::cbafbasedir<< "/DistributedModelCalculations/ShapeLimits/CreateModel.sh " << RunDirectory << " susydatacard.txt" << " 0 0 0 0"; // ASYMPTOTIC LIMITS
1023		}
1024	<	else command << "bash CreateModel.sh " << RunDirectory << " susydatacard.txt" << " 2 0 " << int(0.5 * firstGuess) << " " << int(2*firstGuess); // ASYMPTOTIC LIMITS
1024	>	else command << "bash " << PlottingSetup::cbafbasedir<< "/DistributedModelCalculations/ShapeLimits/CreateModel.sh " << RunDirectory << " susydatacard.txt" << " 2 0 " << int(0.333 * firstGuess) << " " << int(firstGuess); // ASYMPTOTIC LIMITS
1025		dout <<"Going to run : " << command.str() << endl;
1026		int CreatedModelFileExitCode = gSystem->Exec(command.str().c_str());
1027	<	cout << "exit code of limit algorithm (CreateModel.sh) : " << CreatedModelFileExitCode << endl;
1028	<	assert(CreatedModelFileExitCode==0);
1027	>	dout << "exit code of limit algorithm (CreateModel.sh) : " << CreatedModelFileExitCode << endl;
1028	>	if(!(CreatedModelFileExitCode==0)) {
1029	>	write_warning(__FUNCTION__,"Something bad happened. It looks like a shape analysis is not the way to go. ");
1030	>	ShapeDroplet alpha;
1031	>	alpha.observed=-12345; // this is the flag to say watch out something went wrong with the signal ...
1032	>	alpha.SignalIntegral=1;
1033	>	return alpha;
1034	>	}
1035		ShapeDroplet alpha;
1036		alpha.readDroplet(RunDirectory+"/ShapeDropletResult.txt");
1037		alpha.PDF=PDFuncert;
#	Line 785 | Line 1039 | ShapeDroplet LimitsFromShapes(bool asymp
1039		alpha.SignalIntegral=SignalIntegral;
1040		dout << "Done reading limit information from droplet" << endl;
1041		dout << alpha << endl;
1042	<
1042	>	if(asymptotic) {
1043	>	dout << "Doing asymptotic limits, therefore adding an additional round!" << endl;
1044	>	command.str("");
1045	>	command << "bash " << PlottingSetup::cbafbasedir<< "/DistributedModelCalculations/ShapeLimits/CreateModel.sh " << RunDirectory << " susydatacard.txt" << " 1 0 0 " << alpha.expected; // ASYMPTOTIC LIMITS WITH FIRST GUESS
1046	>	CreatedModelFileExitCode = gSystem->Exec(command.str().c_str());
1047	>	dout << "exit code of limit algorithm (CreateModel.sh) : " << CreatedModelFileExitCode << endl;
1048	>	if(!(CreatedModelFileExitCode==0)) {
1049	>	write_warning(__FUNCTION__,"Something bad happened. It looks like a shape analysis is not the way to go. ");
1050	>	ShapeDroplet alpha;
1051	>	alpha.observed=-12345; // this is the flag to say watch out something went wrong with the signal ...
1052	>	alpha.SignalIntegral=1;
1053	>	return alpha;
1054	>	}
1055	>	alpha.readDroplet(RunDirectory+"/ShapeDropletResult.txt");
1056	>	alpha.PDF=PDFuncert;
1057	>	alpha.JSU=JZBscale;
1058	>	alpha.SignalIntegral=SignalIntegral;
1059	>	dout << "Done reading updated limit information from droplet" << endl;
1060	>	dout << alpha << endl;
1061	>	}
1062	>
1063		dout << "Everything is saved in " << RunDirectory << endl;
1064	<	dout << "Will transfer model and datacard over for possible post-processing" << endl;
1065	<	dout << "   1) Make sure models directory exists ... " << std::flush;
1064	>	dout << "Cleaning up ... " << std::flush;
1065	>	/*  dout << "   1) Make sure models directory exists ... " << std::flush;
1066		gSystem->Exec("mkdir -p models/");
1067		dout << " ok!" << endl;
1068		dout << "   2) Deleting any previous model files with the same name ... " << std::flush;
#	Line 810 | Line 1084 | ShapeDroplet LimitsFromShapes(bool asymp
1084		copycommand3 << "cp " << RunDirectory << "/limitfile.root models/model_" << name << "_histo.root";
1085		gSystem->Exec(copycommand3.str().c_str());
1086		dout << " ok!" << endl;
1087	<	dout << "   4) Removing original working directory (" << RunDirectory << ") ... " << std::flush;
814	<	write_warning(__FUNCTION__,"Watch out : need to uncomment the line below to remove the original working directory again");
1087	>	dout << "   4) Removing original working directory (" << RunDirectory << ") ... " << std::flush;*/
1088		//  gSystem->Exec(("rm -r "+RunDirectory).c_str());
1089		dout << " ok!" << endl;
1090	+	write_warning(__FUNCTION__,"Not cleaning up right now!");
1091		delete limcan;
1092		return alpha;
1093		}
1094
1095	<	void PrepareDataShapes(string datajzb, vector<float> jzbbins, float jzbpeakerrordata) {
1096	<	//  dout << "Generating data shape templates ... " << endl;
1095	>	void PrepareDataShapes(string mcjzb, string datajzb, vector<float> jzbbins, float jzbpeakerrordata) {
1096	>
1097		map <  pair<float, float>, map<string, float>  >  xsec;
1098		TFile *datafile = new TFile("../DistributedModelCalculations/StoredShapes.root","RECREATE");
1099		TCanvas *limcan = new TCanvas("limcan","Canvas for calculating limits");
1100		TTree *faketree;
1101		bool dataonly=true;
1102		bool signalonly=false;
1103	<	string mcjzb="JZBforMCinPrepareDataShapes";//this string is not used.
830	<	float jzbpeakerrormc=0;
1103	>
1104		generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"",mcjzb,datajzb,noJES,jzbbins,limcan,"",xsec);
1105	<	// don't need these effects for obs & pred, only for signal!
1106	<	//  generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"peakUp",newjzbexpression(mcjzb,jzbpeakerrormc),newjzbexpression(datajzb,jzbpeakerrordata),noJES,jzbbins,limcan,"",xsec);
1107	<	//  generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"peakDown",newjzbexpression(mcjzb,-jzbpeakerrormc),newjzbexpression(datajzb,-jzbpeakerrordata),noJES,jzbbins,limcan,"",xsec);
1108	<	//  generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"JESUp",mcjzb,datajzb,JESup,jzbbins,limcan,"",xsec);
1109	<	//  generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"JESDown",mcjzb,datajzb,JESdown,jzbbins,limcan,"",xsec);
1105	>
1106	>	generate_mc_shapes(false,faketree,mcjzb,datajzb,jzbbins,cutWeight,cutnJets,"","",xsec);
1107	>	generate_mc_shapes(false,faketree,mcjzb,datajzb,jzbbins,cutWeight,cutnJetsJESup,"","JESUp",xsec);
1108	>	generate_mc_shapes(false,faketree,mcjzb,datajzb,jzbbins,cutWeight,cutnJetsJESdown,"","JESDown",xsec);
1109	>	generate_mc_shapes(false,faketree,mcjzb,datajzb,jzbbins,TCut("weightEffUp"),cutnJets,"","EffUp",xsec);
1110	>	generate_mc_shapes(false,faketree,mcjzb,datajzb,jzbbins,TCut("weightEffDown"),cutnJets,"","EffDown",xsec);
1111	>	generate_mc_shapes(false,faketree,mcjzb,datajzb,jzbbins,cutWeight,cutnJets,"","StatUp",xsec);
1112	>	generate_mc_shapes(false,faketree,mcjzb,datajzb,jzbbins,cutWeight,cutnJets,"","StatDown",xsec);
1113	>	/*
1114		datafile->Close();
1115	+	*/
1116		}
1117

Diff Legend

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