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.7 by buchmann, Fri Apr 27 14:07:39 2012 UTC vs.
Revision 1.11 by buchmann, Mon May 14 15:14:03 2012 UTC

#	Line 44 | Line 44 | void EliminateNegativeEntries(TH1F *hist
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(Tbackground);
56	<	assert(Zbackground);
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()<<"         "<<Tbackground->Integral()<<"         "<<Zbackground->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 117 | Line 261 | TH1F* QuickDraw(TTree *events, string hn
261		return histo;
262		}
263
120	–	/*void do_stat_up(TH1F *h, float sign=1.0) {
121	–	for(int i=1;i<=h->GetNbinsX();i++) {
122	–	h->SetBinContent(i,h->GetBinContent(i)+sign*h->GetBinError(i));
123	–	}
124	–	h->Write();
125	–	}
126	–
127	–	void do_stat_dn(TH1F *h) {
128	–	do_stat_up(h,-1.0);
129	–	}*/
264
265		void SQRT(TH1F *h) {
266		for (int i=1;i<=h->GetNbinsX();i++) {
#	Line 142 | 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 173 | Line 379 | void generate_shapes_for_systematic(bool
379		if(signalonly) {
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);
176	–	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 201 | 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		dout << "SITUATION FOR SIGNAL: " << endl;
418	<	dout << " OSSF     JZB> 0 : " << ZOSSFP->Integral() << "     JZB < 0 :" << ZOSSFN->Integral() << endl;
419	<	dout << " OSOF     JZB> 0 : " << ZOSOFP->Integral() << "     JZB < 0 :" << ZOSOFN->Integral() << endl;
420	<	if(PlottingSetup::RestrictToMassPeak) {
421	<	dout << " OSSF SB  JZB> 0 : " << SBOSSFP->Integral() << "     JZB < 0 :" << SBOSSFN->Integral() << endl;
422	<	dout << " OSOF SB  JZB> 0 : " << SBOSOFP->Integral() << "     JZB < 0 :" << SBOSOFN->Integral() << 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 261 | 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	<	dout << "Stat err in bin " << i << " : " << staterr << endl;
484	<	dout << "    prediction: " << Lpred->GetBinContent(i) << " , observation: " << Lobs->GetBinContent(i) << " --> signal: " << signal->GetBinContent(i) << endl;
485	<	dout << "    we obtain : " << signal->GetBinContent(i)-staterr << " , " << signal->GetBinContent(i)+staterr << endl;
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);
#	Line 286 | 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));
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);
#	Line 700 | Line 925 | ShapeDroplet LimitsFromShapes(bool asymp
925
926		ensure_directory_exists(RunDirectory);
927
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!");
#	Line 723 | 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);
726	–	//  generate_shapes_for_systematic(signalonly,dataonly,limfile,events,"peakUp",newjzbexpression(mcjzb,jzbpeakerrormc),newjzbexpression(datajzb,jzbpeakerrordata),noJES,jzbbins,limcan,addcut,xsec);
727	–	//  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 794 | Line 1018 | ShapeDroplet LimitsFromShapes(bool asymp
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		dout << "exit code of limit algorithm (CreateModel.sh) : " << CreatedModelFileExitCode << endl;
#	Line 815 | 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 << "Cleaning up ... " << std::flush;
1065		/*  dout << "   1) Make sure models directory exists ... " << std::flush;
#	Line 841 | Line 1085 | ShapeDroplet LimitsFromShapes(bool asymp
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());
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.
859	<	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