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.1 by buchmann, Wed Apr 11 15:37:24 2012 UTC vs.
Revision 1.11 by buchmann, Mon May 14 15:14:03 2012 UTC

#	Line 32 | Line 32 | namespace SUSYScanSpace {
32		string SavedMGluname;
33		}
34
35	+	const char* strip_flip_away(string flipped_name) {
36	+	int pos = flipped_name.find("flipped_");
37	+	if(pos>=0&&pos<1000) flipped_name=flipped_name.substr(8,1000);
38	+	return flipped_name.c_str();
39	+	}
40	+
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	+
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	+	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;
235	+
236		TH1F* QuickDraw(TTree *events, string hname, string variable, vector<float> binning, string xlabel, string ylabel, TCut cut, string addcut, bool dataormc, float luminosity, map <  pair<float, float>, map<string, float>  > &xsec) {
237		TH1F *histo = new TH1F(hname.c_str(),hname.c_str(),binning.size()-1,&binning[0]);
238		histo->Sumw2();
#	Line 47 | Line 248 | TH1F* QuickDraw(TTree *events, string hn
248		stringstream mSUGRAweight;
249		float xs = XSForProcessViaAddCutWrapper(addcut,xsec,i);
250		totalxs+=xs;
251	+	mSUGRAweight << "(0.5*(2*(process==" << i << "))*" << xs << ")";
252		events->Draw(drawcommand2.str().c_str(),(essentialcut&&cut&&addcut.c_str())*cutWeight*mSUGRAweight.str().c_str());
253		}
254	<	histo->Scale(1.0/totalxs);
254	>	histo->Scale(1.0/totalxs); // to get back to 1 /pb normalization
255		}
256		events->Draw("eventNum",addcut.c_str(),"goff");
257		float nevents = events->GetSelectedRows();
258	+	QuickDrawNevents=nevents;
259		histo->Scale(luminosity/nevents); // this will give a histogram which is normalized to 1 pb so that any limit will be with respect to 1 pb
260
261		return histo;
262		}
263
61	–	/*void do_stat_up(TH1F *h, float sign=1.0) {
62	–	for(int i=1;i<=h->GetNbinsX();i++) {
63	–	h->SetBinContent(i,h->GetBinContent(i)+sign*h->GetBinError(i));
64	–	}
65	–	h->Write();
66	–	}
67	–
68	–	void do_stat_dn(TH1F *h) {
69	–	do_stat_up(h,-1.0);
70	–	}*/
264
265		void SQRT(TH1F *h) {
266		for (int i=1;i<=h->GetNbinsX();i++) {
#	Line 83 | 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);
354		limfile->cd();
355	<	dout << "Creatig shape templates ";
355	>	dout << "Creating shape templates ";
356		if(identifier!="") dout << "for "<<identifier;
357		dout << " ... " << endl;
358
359		TCut limitnJetcut;
360	+
361		if(JES==noJES) limitnJetcut=cutnJets;
362		else {
363		if(JES==JESdown) limitnJetcut=cutnJetsJESdown;
364		if(JES==JESup) limitnJetcut=cutnJetsJESup;
365		}
366
100	–	stringstream mSUGRAweight;
101	–	if(SUSYScanSpace::SUSYscantype==mSUGRA) {
102	–	//for(int i<
103	–	mSUGRAweight << "bla";
104	–	} else mSUGRAweight << "(1)";
105	–
106	–	write_warning(__FUNCTION__,"Not done yet for mSUGRA");
107	–
367		float zjetsestimateuncert=zjetsestimateuncertONPEAK;
368		float emuncert=emuncertONPEAK;
369		float emsidebanduncert=emsidebanduncertONPEAK;
#	Line 118 | Line 377 | void generate_shapes_for_systematic(bool
377		}
378
379		if(signalonly) {
380	<	cout << "Processing a signal with mcjzb: " << mcjzb << 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);
123	–	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 148 | 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	<	flippedLobs->Add(ZOSSFN);
154	<	flippedLpred->Add(ZOSSFP);
417	>	dout << "SITUATION FOR SIGNAL: " << endl;
418
419	<	if(PlottingSetup::RestrictToMassPeak) {
420	<	Lpred->Add(ZOSOFP,1.0/3);
421	<	Lpred->Add(ZOSOFN,-1.0/3);
159	<	Lpred->Add(SBOSSFP,1.0/3);
160	<	Lpred->Add(SBOSSFN,-1.0/3);
161	<	Lpred->Add(SBOSOFP,1.0/3);
162	<	Lpred->Add(SBOSOFN,-1.0/3);
163	<
164	<	//flipped prediction
165	<	flippedLpred->Add(ZOSOFP,-1.0/3);
166	<	flippedLpred->Add(ZOSOFN,1.0/3);
167	<	flippedLpred->Add(SBOSSFP,-1.0/3);
168	<	flippedLpred->Add(SBOSSFN,1.0/3);
169	<	flippedLpred->Add(SBOSOFP,-1.0/3);
170	<	flippedLpred->Add(SBOSOFN,1.0/3);
171	<
419	>
420	>	if(PlottingSetup::FullMCAnalysis) {
421	>	dout << " OSSF     JZB> 0 : " << ZOSSFP->Integral() << endl;
422		} else {
423	<	Lpred->Add(ZOSOFP,1.0);
424	<	Lpred->Add(ZOSOFN,-1.0);
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	<	//flipped prediction
432	<	flippedLpred->Add(ZOSOFP,-1.0);
433	<	flippedLpred->Add(ZOSOFN,1.0);
431	>
432	>	flippedLobs->Add(ZOSSFN);
433	>	if(!PlottingSetup::FullMCAnalysis) flippedLpred->Add(ZOSSFP);
434	>
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();
462	<	signal->Add(Lpred,-1);
461	>	TH1F *signal = (TH1F*)Lobs->Clone("signal");
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
472		if(identifier=="") {
192	–	write_warning(__FUNCTION__,"Don't know right now how to define stat err on signal");
473		TH1F *signalStatUp = (TH1F*)signal->Clone();
474		signalStatUp->SetName(((string)signal->GetName()+"_StatUp").c_str());
475		signalStatUp->SetTitle(((string)signal->GetTitle()+"_StatUp").c_str());
#	Line 199 | 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	<	signalStatDn->SetBinContent(i,signal->GetBinContent(i)-signal->GetBinError(i));
483	<	signalStatUp->SetBinContent(i,signal->GetBinContent(i)+signal->GetBinError(i));
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		}
496
497		signalStatDn->Write();
#	Line 218 | 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	<	flippedsignalStatDn->SetBinContent(i,flippedsignal->GetBinContent(i)-flippedsignal->GetBinError(i));
513	<	flippedsignalStatUp->SetBinContent(i,flippedsignal->GetBinContent(i)+flippedsignal->GetBinError(i));
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		}
520
521		flippedsignalStatDn->Write();
#	Line 249 | 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 383 | 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_StatDn");
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 402 | 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_StatDn");
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 417 | 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_StatDn");
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 432 | 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_StatDn");
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 448 | 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_StatDn");
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 463 | 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_StatDn");
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 485 | 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_SysDn");
798	>	TH1F *bgSysDn = (TH1F*)pred->Clone("background_SysDown");
799		bgSysDn->Add(predsyserr,-1);
800	+	EliminateNegativeEntries(bgSysDn);
801		bgSysDn->Write();
802		delete predsyserr;
803
#	Line 503 | 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_SysDn");
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 518 | 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_SysDn");
835	>	TH1F *TpredSysDn = (TH1F*)Tpred->Clone("TTbarBackground_SysDown");
836		TpredSysDn->Add(Tpredsyserr,-1);
837	+	EliminateNegativeEntries(TpredSysDn);
838		TpredSysDn->Write();
839		delete Tpredsyserr;
840
#	Line 532 | 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_SysDn");
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 549 | 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_SysDn");
870	>	TH1F *ZpredSysDn = (TH1F*)Zpred->Clone("ZJetsBackground_SysDown");
871		ZpredSysDn->Add(Zpredsyserr,-1);
872	+	EliminateNegativeEntries(ZpredSysDn);
873		ZpredSysDn->Write();
874		delete Zpredsyserr;
875
#	Line 564 | 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_SysDn");
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=="") {
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	<	}
898	>	}*/
899		delete ZOSSFP;
900		delete ZOSOFP;
901		delete ZOSSFN;
#	Line 591 | Line 911 | if(identifier=="") {
911
912		}
913
914	<	ShapeDroplet LimitsFromShapes(TTree *events, string addcut, string name, string mcjzb, string datajzb, vector<float> jzbbins, float jzbpeakerrormc, float jzbpeakerrordata) {
914	>	ShapeDroplet LimitsFromShapes(bool asymptotic, float firstGuess, TTree *events, string addcut, string name, string mcjzb, string datajzb, vector<float> jzbbins, float jzbpeakerrormc, float jzbpeakerrordata) {
915		map <  pair<float, float>, map<string, float>  >  xsec;
916		if(SUSYScanSpace::SUSYscantype==mSUGRA) xsec=getXsec(PlottingSetup::cbafbasedir+"/"+mSUGRAxsFile);
917
598	–	write_info(__FUNCTION__,"Implementing the shape function");
918		time_t timestamp;
919		tm *now;
920		timestamp = time(0);
#	Line 606 | Line 925 | ShapeDroplet LimitsFromShapes(TTree *eve
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;
935	<	TFile *limfile = new TFile((RunDirectory+"/limitfile.root").c_str(),"RECREATE");
934	>	dout << "Run Directory: " << RunDirectory << endl;
935	>	TFile *limfile = new TFile((RunDirectory+"/PRElimitfile.root").c_str(),"RECREATE");
936
937		TIter nextkey(datafile->GetListOfKeys());
938		TKey *key;
#	Line 629 | Line 949 | ShapeDroplet LimitsFromShapes(TTree *eve
949		bool dataonly=false;
950
951		generate_shapes_for_systematic(signalonly,dataonly,limfile,events,"",mcjzb,datajzb,noJES,jzbbins,limcan,addcut,xsec);
632	–	generate_shapes_for_systematic(signalonly,dataonly,limfile,events,"peakUp",newjzbexpression(mcjzb,jzbpeakerrormc),newjzbexpression(datajzb,jzbpeakerrordata),noJES,jzbbins,limcan,addcut,xsec);
633	–	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
637	–	// JSU & PDF: global factors
638	–	write_warning(__FUNCTION__,"NEED TO IMPLEMENT JSU & PDF ");
639	–
955		float JZBScaleUncert=0.05; //if you modify this value please also adapt it in Systematics.C not only here in ShapeLimit.C
956		float scaledown=0,scaleup=0,scalesyst=0;
957		float mceff=0,mcefferr=0;
#	Line 647 | Line 962 | ShapeDroplet LimitsFromShapes(TTree *eve
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;
970	>	float PDFuncert=0;
971		int NPdfs = get_npdfs(events);
972		if(SUSYScanSpace::SUSYscantype!=mSUGRA) PDFuncert = get_pdf_uncertainty(events, flipped, mcjzb, requireZ, Neventsinfile, SUSYScanSpace::SUSYscantype, xsec, NPdfs, addcut);
973
974	+	float JZBscale=scaledown;
975	+	if(scaleup>scaledown) JZBscale=scaleup;
976	+	dout << endl;
977	+	dout << endl;
978	+	dout << "Will use the following scalar (non-shape) systematics : " << endl;
979	+	dout << "       JZB Scale (JSU) : " << JZBscale << endl;
980	+	dout << "       PDF             : " << PDFuncert << endl;
981	+
982	+	float SignalIntegral;
983	+	if(flipped) {
984	+	TH1F *flisi = (TH1F*)(limfile->Get("flipped_signal"));
985	+	SignalIntegral= flisi ->Integral();
986	+	delete flisi;
987	+	}
988	+	else {
989	+	TH1F *flisi = (TH1F*)(limfile->Get("signal"));
990	+	SignalIntegral= flisi ->Integral();
991	+	delete flisi;
992	+	}
993	+
994	+	TFile *final_limfile = new TFile((RunDirectory+"/limitfile.root").c_str(),"RECREATE");
995	+
996	+	TIter nnextkey(limfile->GetListOfKeys());
997	+	TKey *nkey;
998	+
999	+	if(SignalIntegral==0) SignalIntegral=1; //case when the integral is zero - don't try to normalize, just leave it as it is
1000	+
1001	+	while ((nkey = (TKey*)nnextkey()))
1002	+	{
1003	+	TH1F *obj =(TH1F*) nkey->ReadObj();
1004	+	if(flipped&&!Contains(obj->GetName(),"flipped")) continue;
1005	+	if(!flipped&&Contains(obj->GetName(),"flipped")) continue;
1006	+	if(flipped) obj->SetName(strip_flip_away(obj->GetName()));
1007	+	final_limfile->cd();
1008	+	if(Contains(obj->GetName(),"signal")) {
1009	+	obj->Scale(1.0/SignalIntegral);
1010	+	}
1011	+	obj->Write();
1012	+	}
1013
1014	<	prepare_datacard(limfile);
658	<
659	<	limfile->Close();
660	<	write_info(__FUNCTION__,"limitfile.root and datacard.txt have been generated");
661	<
662	<	write_error(__FUNCTION__,"Implementation of limit calculation is still missing");
1014	>	prepare_limit_datacard(RunDirectory,final_limfile,JZBscale,PDFuncert);
1015
1016	<	write_warning(__FUNCTION__,"Once you're done implementing the limit computation, please uncomment the following line to enable clean up.");
1017	<	dout << "Everything is saved in " << RunDirectory << endl;
1018	<	//  gSystem->Exec(("rm -r "+RunDirectory).c_str());
1019	<	write_warning(__FUNCTION__,"Once you're done implementing the limit computation, please store the information below.");
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 " << 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 " << 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	>	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.expected=3;
1032	<
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;
1038	+	alpha.JSU=JZBscale;
1039	+	alpha.SignalIntegral=SignalIntegral;
1040	+	dout << "Done reading limit information from droplet" << endl;
1041	+	dout << alpha << endl;
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 << "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;
1069	+	stringstream delcommand;
1070	+	delcommand << "rm models/model_" << name << ".root";
1071	+	gSystem->Exec(delcommand.str().c_str());
1072	+	stringstream delcommand2;
1073	+	delcommand2 << "rm models/model_" << name << ".txt";
1074	+	gSystem->Exec(delcommand2.str().c_str());
1075	+	dout << " ok!" << endl;
1076	+	dout << "   3) Transfering the new model files ... " << std::flush;
1077	+	stringstream copycommand;
1078	+	copycommand << "cp " << RunDirectory << "/model.root models/model_" << name << ".root";
1079	+	gSystem->Exec(copycommand.str().c_str());
1080	+	stringstream copycommand2;
1081	+	copycommand2 << "cp " << RunDirectory << "/susydatacard.txt models/model_" << name << ".txt";
1082	+	gSystem->Exec(copycommand2.str().c_str());
1083	+	stringstream copycommand3;
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;*/
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;
682	–	string mcjzb="JZBforMCinPrepareDataShapes";//this string is not used.
683	–	float jzbpeakerrormc=0;
684	–	generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"",mcjzb,datajzb,noJES,jzbbins,limcan,"",xsec);
685	–	// don't need these effects for obs & pred, only for signal!
686	–	//  generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"peakUp",newjzbexpression(mcjzb,jzbpeakerrormc),newjzbexpression(datajzb,jzbpeakerrordata),noJES,jzbbins,limcan,"",xsec);
687	–	//  generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"peakDown",newjzbexpression(mcjzb,-jzbpeakerrormc),newjzbexpression(datajzb,-jzbpeakerrordata),noJES,jzbbins,limcan,"",xsec);
688	–	//  generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"JESUp",mcjzb,datajzb,JESup,jzbbins,limcan,"",xsec);
689	–	//  generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"JESDown",mcjzb,datajzb,JESdown,jzbbins,limcan,"",xsec);
1103
1104	+	generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"",mcjzb,datajzb,noJES,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