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.2 by buchmann, Thu Apr 12 14:58:25 2012 UTC vs.
Revision 1.12 by buchmann, Tue May 15 19:37:47 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);
62	<
63	<	ofstream datacard;
64	<
65	<	datacard.open ((RunDirectory+"/susydatacard.txt").c_str());
66	<	datacard << "Writing this to a file.\n";
67	<	datacard << "imax 1\n"; // number of channels
68	<	datacard << "jmax 2\n"; // number of backgrounds (Z+Jets prediction and TTbar prediction)
69	<	datacard << "kmax *\n"; // number of nuisance parameters (sources of systematic uncertainties)
70	<	datacard << "---------------\n";
71	<	datacard << "shapes * * limitfile.root $PROCESS $PROCESS_$SYSTEMATIC\n";
72	<	datacard << "---------------\n";
73	<	datacard << "bin 1\n";
74	<	datacard << "observation "<<dataob->Integral()<<"\n";
75	<	datacard << "------------------------------\n";
76	<	datacard << "bin             1          1          1\n";
77	<	datacard << "process         signal     TTbarBackground     ZJetsBackground\n";
78	<	datacard << "process         0          1          2\n";
79	<	datacard << "rate            "<<signal->Integral()<<"         "<<background1->Integral()<<"         "<<background2->Integral()<<"\n";
80	<	datacard << "--------------------------------\n";
81	<	datacard << "lumi     lnN    " << 1+ PlottingSetup::lumiuncert << "       1.0       1.0    luminosity uncertainty\n"; // only affects MC -> signal!
82	<	// datacard << "bgnorm   lnN    1.00       1.4  uncertainty on our prediction (40%)\n";
83	<	datacard << "Stat   shape    1          1          1    statistical uncertainty\n";
84	<	datacard << "Sys    shape    -          1          1    systematic uncertainty\n";
85	<	datacard << "JES    shape    1          -          -    uncertainty on Jet Energy Scale\n";
86	<	datacard << "JSU      lnN    " << 1+uJSU << "          -          -    JZB Scale Uncertainty\n";
87	<	if(uPDF>0) datacard << "PDF      lnN    " << 1+uPDF << "          -          -    uncertainty from PDFs\n";
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	>	vector<float> get_expected_points(float observed,int npoints, string RunDirectory, bool doPlot=false) {
48	>	TF1 *gaus = new TF1("gaus","gaus(0)",0,10*observed);
49	>	gaus->SetParameters(1,observed,2*observed);
50	>	gaus->SetParameter(0,1/(gaus->Integral(0,10*observed)));
51	>	float currentpoint=0.01;
52	>	float stepsize=observed/100;
53	>	vector<float> points;
54	>	while(currentpoint<25*observed && points.size()<npoints) {
55	>
56	>	if(gaus->Integral(0,currentpoint)>((points.size()+1.0)/(npoints+2.0))) {
57	>	if(Verbosity>0) dout << "Added points for calculation at " << currentpoint << " (integral is " << gaus->Integral(0,currentpoint) << ")" << endl;
58	>	points.push_back(currentpoint);
59	>	if(points.size()==npoints) break;
60	>	}
61	>	currentpoint+=stepsize;
62	>	}
63	>	if(doPlot) {
64	>	gaus->SetName("Expected limit computation points");
65	>	gaus->SetTitle("Expected limit computation points");
66	>	TCanvas *can = new TCanvas("can","can");
67	>	gaus->Draw();
68	>	TLine *lines[npoints];
69	>	for(int i=0;i<npoints;i++) {
70	>	lines[i] = new TLine(points[i],0,points[i],gaus->GetMaximum());
71	>	lines[i]->SetLineColor(kBlue);
72	>	lines[i]->Draw();
73	>	}
74	>	can->SaveAs((RunDirectory+"/DistributionOfComputedPoints.png").c_str());
75	>	delete can;
76	>	for(int i=0;i<npoints;i++) delete lines[i];
77	>	}
78	>	delete gaus;
79	>
80	>	return points;
81	>	}
82	>
83	>	void prepare_MC_datacard(string RunDirectory, TFile *f, float uJSU, float uPDF) {
84	>	TH1F *dataob = (TH1F*)f->Get("data_obs");
85	>	TH1F *signal = (TH1F*)f->Get("signal");
86	>	assert(dataob);
87	>	assert(signal);
88	>
89	>	write_warning(__FUNCTION__,"The following two defs for th1's are fake");TH1F *Tbackground; TH1F *Zbackground;
90	>	write_warning(__FUNCTION__,"Not correctly implemented yet for MC");
91	>	ofstream datacard;
92	>	write_warning(__FUNCTION__,"Need to rethink systematics we want to use !");
93	>	datacard.open ((RunDirectory+"/susydatacard.txt").c_str());
94	>	dout << "Writing this to a file.\n";
95	>	datacard << "imax 1\n"; // number of channels
96	>	datacard << "jmax 2\n"; // number of backgrounds (Z+Jets prediction and TTbar prediction)
97	>	datacard << "kmax *\n"; // number of nuisance parameters (sources of systematic uncertainties)
98	>	datacard << "---------------\n";
99	>	datacard << "shapes * * limitfile.root $PROCESS $PROCESS_$SYSTEMATIC\n";
100	>	datacard << "---------------\n";
101	>	datacard << "bin 1\n";
102	>	datacard << "observation "<<dataob->Integral()<<"\n";
103	>	datacard << "------------------------------\n";
104	>	datacard << "bin             1          1          1\n";
105	>	datacard << "process         signal     TTbarBackground     ZJetsBackground\n";
106	>	datacard << "process         0          1          2\n";
107	>	datacard << "rate            "<<signal->Integral()<<"         "<<Tbackground->Integral()<<"         "<<Zbackground->Integral()<<"\n";
108	>	datacard << "--------------------------------\n";
109	>	datacard << "lumi     lnN    " << 1+ PlottingSetup::lumiuncert << "       -       -    luminosity uncertainty\n"; // only affects MC -> signal!
110	>	datacard << "Stat   shape    -          1          -    statistical uncertainty (ttbar)\n";
111	>	datacard << "Stat   shape    -          -          1    statistical uncertainty (zjets)\n";
112	>	datacard << "Sys    shape    -          1          -    systematic uncertainty on ttbar\n";
113	>	datacard << "Sys    shape    -          -          1    systematic uncertainty on zjets\n";
114	>	datacard << "Stat   shape    1          -          -    statistical uncertainty (signal)\n";
115	>	datacard << "JES    shape    1          -          -    uncertainty on Jet Energy Scale\n";
116	>	datacard << "JSU      lnN    " << 1+uJSU << "          -          -    JZB Scale Uncertainty\n";
117	>	if(uPDF>0) datacard << "PDF      lnN    " << 1+uPDF << "          -          -    uncertainty from PDFs\n";
118	>	datacard.close();
119	>	}
120	>
121	>	void prepare_datadriven_datacard(string RunDirectory, TFile *f, float uJSU, float uPDF) {
122	>	TH1F *dataob = (TH1F*)f->Get("data_obs");
123	>	TH1F *signal = (TH1F*)f->Get("signal");
124	>	TH1F *Tbackground = (TH1F*)f->Get("TTbarBackground");
125	>	TH1F *Zbackground  = (TH1F*)f->Get("ZJetsBackground");
126	>
127	>	assert(dataob);
128	>	assert(signal);
129	>	assert(Tbackground);
130	>	assert(Zbackground);
131
132	+
133	+	ofstream datacard;
134	+	write_warning(__FUNCTION__,"Need to rethink systematics we want to use !");
135	+	datacard.open ((RunDirectory+"/susydatacard.txt").c_str());
136	+	datacard << "Writing this to a file.\n";
137	+	datacard << "imax " << dataob->GetNbinsX() << "\n"; // number of channels
138	+	datacard << "jmax 2\n"; // number of backgrounds (Z+Jets prediction and TTbar prediction)
139	+	datacard << "kmax *\n"; // number of nuisance parameters (sources of systematic uncertainties)
140	+	datacard << "---------------\n";
141	+	datacard << "bin\t";
142	+	for(int i=1;i<=dataob->GetNbinsX();i++) datacard << " " << dataob->GetBinLowEdge(i) << "to" << dataob->GetBinLowEdge(i)+dataob->GetBinWidth(i) << " \t";
143	+	datacard << "\n";
144	+	datacard << "observation\t";
145	+	for(int i=1;i<=dataob->GetNbinsX();i++) datacard << " " << dataob->GetBinContent(i) << " \t";
146	+	datacard<<"\n";
147	+	datacard << "------------------------------\n";
148	+	datacard << "bin\t";
149	+	for(int i=1;i<=dataob->GetNbinsX();i++) {
150	+	datacard << " " << dataob->GetBinLowEdge(i) << "to" << dataob->GetBinLowEdge(i)+dataob->GetBinWidth(i) << "\t" <<
151	+	" " << dataob->GetBinLowEdge(i) << "to" << dataob->GetBinLowEdge(i)+dataob->GetBinWidth(i) << "\t" <<
152	+	" " << dataob->GetBinLowEdge(i) << "to" << dataob->GetBinLowEdge(i)+dataob->GetBinWidth(i) << "\t";
153	+	}
154	+	datacard << "\n";
155	+	datacard << "process\t";
156	+	for(int i=1;i<=dataob->GetNbinsX();i++) datacard << "\t signal \t TTbarBackground \t ZJetsBackground";
157	+	datacard << "\n";
158	+	datacard << "process\t";
159	+	for(int i=1;i<=dataob->GetNbinsX();i++) datacard << "\t 0 \t 1 \t 2";
160	+	datacard << "\n";
161	+
162	+
163	+	datacard << "rate\t";
164	+	for(int i=1;i<=dataob->GetNbinsX();i++) datacard << "\t " << signal->GetBinContent(i) << " \t " << Tbackground->GetBinContent(i) << " \t " << Zbackground->GetBinContent(i) << "\t";
165	+	datacard<<"\n";
166	+
167	+	datacard << "--------------------------------\n";
168	+
169	+
170	+	datacard << "lumi     lnN    \t";
171	+	for(int i=1;i<=dataob->GetNbinsX();i++) datacard << " " << 1+ PlottingSetup::lumiuncert << "\t - \t -";
172	+	datacard << "       luminosity uncertainty\n"; // only affects MC -> signal!
173	+
174	+	// Statistical uncertainty
175	+	datacard << "Stat     lnN    \t";
176	+	for(int i=1;i<=dataob->GetNbinsX();i++) {
177	+	//Signal
178	+	float central = signal->GetBinContent(i);
179	+	float up      = ((TH1F*)f->Get("signal_StatDown"))->GetBinContent(i);
180	+	float down    = ((TH1F*)f->Get("signal_StatUp"))->GetBinContent(i);
181	+	float suncert=0;
182	+	if(central>0) {
183	+	if(abs(up-central)>abs(down-central)) suncert=abs(up-central)/central;
184	+	else suncert=abs(central-down)/central;
185	+	}
186	+
187	+	//TTbar
188	+	central = Tbackground->GetBinContent(i);
189	+	up      = ((TH1F*)f->Get("TTbarBackground_StatUp"))->GetBinContent(i);
190	+	down    = ((TH1F*)f->Get("TTbarBackground_StatDown"))->GetBinContent(i);
191	+	float tuncert=0;
192	+	if(central>0) {
193	+	if(abs(up-central)>abs(down-central)) tuncert=abs(up-central)/central;
194	+	else tuncert=abs(central-down)/central;
195	+	}
196	+	//ZJets
197	+	central = Zbackground->GetBinContent(i);
198	+	up      = ((TH1F*)f->Get("ZJetsBackground_StatUp"))->GetBinContent(i);
199	+	down    = ((TH1F*)f->Get("ZJetsBackground_StatDown"))->GetBinContent(i);
200	+	float zuncert=0;
201	+	if(central>0) {
202	+	if(abs(up-central)>abs(down-central)) zuncert=abs(up-central)/central;
203	+	else zuncert=abs(central-down)/central;
204	+	}
205	+	datacard << " " << 1+suncert << " \t " << 1+tuncert << "\t " << 1+zuncert << " \t ";
206	+	}
207	+
208	+	datacard << "       statistical uncertainty\n";
209	+
210	+	// Statistical uncertainty
211	+	datacard << "Sys     lnN    \t";
212	+	for(int i=1;i<=dataob->GetNbinsX();i++) {
213	+	float central = Tbackground->GetBinContent(i);
214	+	float up      = ((TH1F*)f->Get("TTbarBackground_SysUp"))->GetBinContent(i);
215	+	float down    = ((TH1F*)f->Get("TTbarBackground_SysDown"))->GetBinContent(i);
216	+	float tuncert=0;
217	+	if(central>0) {
218	+	if(abs(up-central)>abs(down-central)) tuncert=abs(up-central)/central;
219	+	else tuncert=abs(central-down)/central;
220	+	}
221	+	central = Zbackground->GetBinContent(i);
222	+	up      = ((TH1F*)f->Get("ZJetsBackground_SysUp"))->GetBinContent(i);
223	+	down    = ((TH1F*)f->Get("ZJetsBackground_SysDown"))->GetBinContent(i);
224	+	float zuncert=0;
225	+	if(central>0) {
226	+	if(abs(up-central)>abs(down-central)) zuncert=abs(up-central)/central;
227	+	else zuncert=abs(central-down)/central;
228	+	}
229	+	datacard << " - \t " << 1+tuncert << "\t " << 1+zuncert << " \t ";
230	+	}
231	+
232	+	datacard << "       systematic uncertainty\n";
233	+
234	+
235	+	float JESup = ((TH1F*)f->Get("signal_JESUp"))->Integral();
236	+	float JESdn = ((TH1F*)f->Get("signal_JESDown"))->Integral();
237	+	float central = signal->Integral();
238	+	float uJES=0;
239	+	if(abs(JESup-central)>abs(JESdn-central)) uJES=abs(JESup-central)/central;
240	+	else uJES=abs(central-JESdn)/central;
241	+	datacard << "JES    lnN";
242	+	for(int i=1;i<=dataob->GetNbinsX();i++) datacard << " " << 1+uJES << "\t - \t  - \t";
243	+	datacard << "uncertainty on Jet Energy Scale\n";
244	+
245	+	datacard << "JSU      lnN    ";
246	+	for(int i=1;i<=dataob->GetNbinsX();i++) datacard << " " << 1+uJSU << "\t - \t  - \t";
247	+	datacard << "JZB Scale Uncertainty\n";
248	+
249	+	if(uPDF>0) {
250	+	datacard << "PDF      lnN    ";
251	+	for(int i=1;i<=dataob->GetNbinsX();i++) datacard << " " << 1+uPDF << "\t - \t - \t";
252	+	datacard << "uncertainty from PDFs\n";
253	+	}
254	+
255	+	datacard.close();
256	+	}
257	+
258	+	void prepare_limit_datacard(string RunDirectory, TFile *f, float uJSU, float uPDF) {
259
260	<	// datacard << "peak   shape    1          1    uncertainty on signal resolution. Assume the histogram is a 2 sigma shift, \n";
261	<	// datacard << "#                                so divide the unit gaussian by 2 before doing the interpolation\n";
262	<	datacard.close();
260	>	if(FullMCAnalysis) {
261	>	//dealing with MC analysis - need to use shapes.
262	>	prepare_MC_datacard(RunDirectory,f,uJSU,uPDF);
263	>	} else {
264	>	//doing mutibin analysis
265	>	prepare_datadriven_datacard(RunDirectory,f,uJSU,uPDF);
266	>	}
267	>
268		}
269
270		float QuickDrawNevents=0;
#	Line 109 | Line 284 | TH1F* QuickDraw(TTree *events, string hn
284		stringstream mSUGRAweight;
285		float xs = XSForProcessViaAddCutWrapper(addcut,xsec,i);
286		totalxs+=xs;
287	+	mSUGRAweight << "(0.5*(2*(process==" << i << "))*" << xs << ")";
288		events->Draw(drawcommand2.str().c_str(),(essentialcut&&cut&&addcut.c_str())*cutWeight*mSUGRAweight.str().c_str());
289		}
290	<	histo->Scale(1.0/totalxs);
290	>	histo->Scale(1.0/totalxs); // to get back to 1 /pb normalization
291		}
292		events->Draw("eventNum",addcut.c_str(),"goff");
293		float nevents = events->GetSelectedRows();
#	Line 121 | Line 297 | TH1F* QuickDraw(TTree *events, string hn
297		return histo;
298		}
299
124	–	/*void do_stat_up(TH1F *h, float sign=1.0) {
125	–	for(int i=1;i<=h->GetNbinsX();i++) {
126	–	h->SetBinContent(i,h->GetBinContent(i)+sign*h->GetBinError(i));
127	–	}
128	–	h->Write();
129	–	}
130	–
131	–	void do_stat_dn(TH1F *h) {
132	–	do_stat_up(h,-1.0);
133	–	}*/
300
301		void SQRT(TH1F *h) {
302		for (int i=1;i<=h->GetNbinsX();i++) {
#	Line 146 | Line 312 | TH1F* Multiply(TH1F *h1, TH1F *h2) {
312		return h;
313		}
314
315	+
316	+	string ReduceNumericHistoName(string origname) {
317	+	int pos = origname.find("__h_");
318	+	if(pos == -1) return origname;
319	+	return origname.substr(0,pos);
320	+	}
321	+
322	+
323	+
324	+	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) {
325	+	//weight can be varied to take care of efficiency errors (weightEffUp, weightEffDown)
326	+	vector<float> mbinning;
327	+	mbinning.push_back(-8000);
328	+	for(int i=binning.size()-1;i>=0;i--) mbinning.push_back(-binning[i]);
329	+	mbinning.push_back(0);
330	+	for(int i=0;i<binning.size();i++) mbinning.push_back(binning[i]);
331	+	mbinning.push_back(8000);
332	+
333	+
334	+	TCanvas *shapecan = new TCanvas("shapecan","shapecan");
335	+	TH1F* Signal;
336	+
337	+	stringstream sInverseCutWeight;
338	+	sInverseCutWeight << "(1.0/" << (const char*)cutWeight<<")";
339	+	TCut InverseCutWeight(sInverseCutWeight.str().c_str());
340	+	if(weight==cutWeight) InverseCutWeight = TCut("1.0");
341	+	if(!dotree) {
342	+	//dealing with ALL MC samples (when we save the standard file)
343	+	THStack McPredicted = allsamples.DrawStack("McPred",mcjzb,mbinning, "JZB [GeV]", "events", ((cutmass&&cutOSSF&&JetCut)*(weight*InverseCutWeight)),mc,luminosity);
344	+	int nhists = McPredicted.GetHists()->GetEntries();
345	+	for (int i = 0; i < nhists; i++) {
346	+	TH1* hist = static_cast<TH1*>(McPredicted.GetHists()->At(i));
347	+	//          cout << hist->GetName() << " : " << hist->Integral() << endl;
348	+	//          cout << "     " << ReduceNumericHistoName(hist->GetName()) << endl;
349	+	if(identifier=="StatUp") {
350	+	float appliedweight = hist->Integral() / hist->GetEntries();
351	+	for(int ibin=1;ibin<=hist->GetNbinsX();ibin++) hist->SetBinContent(ibin,hist->GetBinContent(ibin)+appliedweight*TMath::Sqrt(hist->GetBinContent(ibin)/appliedweight));
352	+	}
353	+	if(identifier=="StatDown") {
354	+	float appliedweight = hist->Integral() / hist->GetEntries();
355	+	for(int ibin=1;ibin<=hist->GetNbinsX();ibin++) hist->SetBinContent(ibin,hist->GetBinContent(ibin)-appliedweight*TMath::Sqrt(hist->GetBinContent(ibin)/appliedweight));
356	+	}
357	+	hist->SetName(("mc_"+ReduceNumericHistoName(hist->GetName())+"_"+identifier).c_str());
358	+	hist->Write();
359	+	}
360	+	} else {
361	+	string histoname="mc_signal";
362	+	if(identifier!="") histoname=histoname+"_"+identifier;
363	+	Signal = QuickDraw(signalevents,histoname,mcjzb,binning,"JZB", "events",((cutmass&&cutOSSF&&JetCut&&basiccut))*(weight*InverseCutWeight),addcut,mc,luminosity,xsec);
364	+	if(identifier=="StatUp") {
365	+	float appliedweight = Signal->Integral() / Signal->GetEntries();
366	+	for(int ibin=1;ibin<=Signal->GetNbinsX();ibin++) Signal->SetBinContent(ibin,Signal->GetBinContent(ibin)+appliedweight*TMath::Sqrt(Signal->GetBinContent(ibin)/appliedweight));
367	+	}
368	+	if(identifier=="StatDown") {
369	+	float appliedweight = Signal->Integral() / Signal->GetEntries();
370	+	for(int ibin=1;ibin<=Signal->GetNbinsX();ibin++) Signal->SetBinContent(ibin,Signal->GetBinContent(ibin)-appliedweight*TMath::Sqrt(Signal->GetBinContent(ibin)/appliedweight));
371	+	}
372	+	Signal->Write();
373	+	}
374	+
375	+	/*
376	+	*
377	+	* Jet energy scale (easy, use different JetCut)
378	+	* JZB Scale Uncertainty (will be a number - signal only)
379	+	* Efficiency (will be weightEffUp, weightEffDown)
380	+	* PDFs (signal only) -> Will be a number yet again, computed in the traditional way
381	+	*/
382	+	delete shapecan;
383	+
384	+	}
385	+
386	+
387		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) {
388	+
389		binning.push_back(8000);
390		limfile->cd();
391	<	dout << "Creatig shape templates ";
391	>	dout << "Creating shape templates ";
392		if(identifier!="") dout << "for "<<identifier;
393		dout << " ... " << endl;
394
395		TCut limitnJetcut;
396	+
397		if(JES==noJES) limitnJetcut=cutnJets;
398		else {
399		if(JES==JESdown) limitnJetcut=cutnJetsJESdown;
400		if(JES==JESup) limitnJetcut=cutnJetsJESup;
401		}
402
163	–	stringstream mSUGRAweight;
164	–	if(SUSYScanSpace::SUSYscantype==mSUGRA) {
165	–	//for(int i<
166	–	mSUGRAweight << "bla";
167	–	} else mSUGRAweight << "(1)";
168	–
169	–	write_warning(__FUNCTION__,"Not done yet for mSUGRA");
170	–
403		float zjetsestimateuncert=zjetsestimateuncertONPEAK;
404		float emuncert=emuncertONPEAK;
405		float emsidebanduncert=emsidebanduncertONPEAK;
#	Line 181 | Line 413 | void generate_shapes_for_systematic(bool
413		}
414
415		if(signalonly) {
416	<	cout << "Processing a signal with mcjzb: " << mcjzb << " (identifier: " << identifier << ")" << endl;
416	>	dout << "Processing a signal with mcjzb: " << mcjzb << " (identifier: '" << identifier << "')" << endl;
417		TH1F *ZOSSFP = QuickDraw(signalevents,"ZOSSFP",mcjzb,binning,     "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
186	–	TH1F *ZOSOFP = QuickDraw(signalevents,"ZOSOFP",mcjzb,binning,     "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
418		TH1F *ZOSSFN = QuickDraw(signalevents,"ZOSSFN","-"+mcjzb,binning, "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
419	<	TH1F *ZOSOFN = QuickDraw(signalevents,"ZOSOFN","-"+mcjzb,binning, "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
419	>	TH1F *ZOSOFP;
420	>	TH1F *ZOSOFN;
421	>
422	>	if(!PlottingSetup::FullMCAnalysis) {
423	>	ZOSOFP = QuickDraw(signalevents,"ZOSOFP",mcjzb,binning,     "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
424	>	ZOSOFN = QuickDraw(signalevents,"ZOSOFN","-"+mcjzb,binning, "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,addcut,mc,luminosity,xsec);
425	>	}
426
427		TH1F *SBOSSFP;
428		TH1F *SBOSOFP;
429		TH1F *SBOSSFN;
430		TH1F *SBOSOFN;
431
432	<	if(PlottingSetup::RestrictToMassPeak) {
432	>	if(PlottingSetup::RestrictToMassPeak&&!PlottingSetup::FullMCAnalysis) {
433		SBOSSFP = QuickDraw(signalevents,"SBOSSFP",mcjzb,binning, "JZB", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
434		SBOSOFP = QuickDraw(signalevents,"SBOSOFP",mcjzb,binning, "JZB", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
435		SBOSSFN = QuickDraw(signalevents,"SBOSSFN","-"+mcjzb,binning, "JZB", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,addcut,mc,luminosity,xsec);
#	Line 211 | Line 448 | void generate_shapes_for_systematic(bool
448		TH1F *flippedLpred = new TH1F("flippedLpred","flippedLpred",binning.size()-1,&binning[0]);
449
450		Lobs->Add(ZOSSFP);
451	<	Lpred->Add(ZOSSFN);
451	>	if(!PlottingSetup::FullMCAnalysis) Lpred->Add(ZOSSFN);
452
453	<	cout << "SITUATION FOR SIGNAL: " << endl;
217	<	cout << " OSSF     JZB> 0 : " << ZOSSFP->Integral() << "     JZB < 0 :" << ZOSSFN->Integral() << endl;
218	<	cout << " OSOF     JZB> 0 : " << ZOSOFP->Integral() << "     JZB < 0 :" << ZOSOFN->Integral() << endl;
219	<	if(PlottingSetup::RestrictToMassPeak) {
220	<	cout << " OSSF SB  JZB> 0 : " << SBOSSFP->Integral() << "     JZB < 0 :" << SBOSSFN->Integral() << endl;
221	<	cout << " OSOF SB  JZB> 0 : " << SBOSOFP->Integral() << "     JZB < 0 :" << SBOSOFN->Integral() << endl;
222	<	}
453	>	dout << "SITUATION FOR SIGNAL: " << endl;
454
455
456	+	if(PlottingSetup::FullMCAnalysis) {
457	+	dout << " OSSF     JZB> 0 : " << ZOSSFP->Integral() << endl;
458	+	} else {
459	+	dout << " OSSF     JZB> 0 : " << ZOSSFP->Integral() << "     JZB < 0 :" << ZOSSFN->Integral() << endl;
460	+	dout << " OSOF     JZB> 0 : " << ZOSOFP->Integral() << "     JZB < 0 :" << ZOSOFN->Integral() << endl;
461	+	if(PlottingSetup::RestrictToMassPeak&&!PlottingSetup::FullMCAnalysis) {
462	+	dout << " OSSF SB  JZB> 0 : " << SBOSSFP->Integral() << "     JZB < 0 :" << SBOSSFN->Integral() << endl;
463	+	dout << " OSOF SB  JZB> 0 : " << SBOSOFP->Integral() << "     JZB < 0 :" << SBOSOFN->Integral() << endl;
464	+	}
465	+	}
466	+
467	+
468		flippedLobs->Add(ZOSSFN);
469	<	flippedLpred->Add(ZOSSFP);
469	>	if(!PlottingSetup::FullMCAnalysis) flippedLpred->Add(ZOSSFP);
470
471	<	if(PlottingSetup::RestrictToMassPeak) {
472	<	Lpred->Add(ZOSOFP,1.0/3);
473	<	Lpred->Add(ZOSOFN,-1.0/3);
474	<	Lpred->Add(SBOSSFP,1.0/3);
475	<	Lpred->Add(SBOSSFN,-1.0/3);
476	<	Lpred->Add(SBOSOFP,1.0/3);
477	<	Lpred->Add(SBOSOFN,-1.0/3);
478	<
479	<	//flipped prediction
480	<	flippedLpred->Add(ZOSOFP,-1.0/3);
481	<	flippedLpred->Add(ZOSOFN,1.0/3);
482	<	flippedLpred->Add(SBOSSFP,-1.0/3);
483	<	flippedLpred->Add(SBOSSFN,1.0/3);
484	<	flippedLpred->Add(SBOSOFP,-1.0/3);
485	<	flippedLpred->Add(SBOSOFN,1.0/3);
486	<
487	<	} else {
488	<	Lpred->Add(ZOSOFP,1.0);
489	<	Lpred->Add(ZOSOFN,-1.0);
490	<
491	<	//flipped prediction
492	<	flippedLpred->Add(ZOSOFP,-1.0);
493	<	flippedLpred->Add(ZOSOFN,1.0);
471	>	if(!PlottingSetup::FullMCAnalysis) {
472	>	if(PlottingSetup::RestrictToMassPeak) {
473	>	Lpred->Add(ZOSOFP,1.0/3);
474	>	Lpred->Add(ZOSOFN,-1.0/3);
475	>	Lpred->Add(SBOSSFP,1.0/3);
476	>	Lpred->Add(SBOSSFN,-1.0/3);
477	>	Lpred->Add(SBOSOFP,1.0/3);
478	>	Lpred->Add(SBOSOFN,-1.0/3);
479	>
480	>	//flipped prediction
481	>	flippedLpred->Add(ZOSOFP,-1.0/3);
482	>	flippedLpred->Add(ZOSOFN,1.0/3);
483	>	flippedLpred->Add(SBOSSFP,-1.0/3);
484	>	flippedLpred->Add(SBOSSFN,1.0/3);
485	>	flippedLpred->Add(SBOSOFP,-1.0/3);
486	>	flippedLpred->Add(SBOSOFN,1.0/3);
487	>	} else {
488	>	Lpred->Add(ZOSOFP,1.0);
489	>	Lpred->Add(ZOSOFN,-1.0);
490	>
491	>	//flipped prediction
492	>	flippedLpred->Add(ZOSOFP,-1.0);
493	>	flippedLpred->Add(ZOSOFN,1.0);
494	>	}
495		}
496
497		TH1F *signal = (TH1F*)Lobs->Clone("signal");
498	<	signal->Add(Lpred,-1);
498	>	if(!PlottingSetup::FullMCAnalysis) signal->Add(Lpred,-1);
499		signal->SetName(signalname.c_str());
500		signal->SetTitle(signalname.c_str());
501		signal->Write();
502
503		TH1F *flippedsignal = (TH1F*)flippedLobs->Clone();
504	<	flippedsignal->Add(flippedLpred,-1);
504	>	if(!PlottingSetup::FullMCAnalysis) flippedsignal->Add(flippedLpred,-1);
505		flippedsignal->SetName(("flipped_"+signalname).c_str());
506		flippedsignal->Write();
507
#	Line 271 | Line 515 | void generate_shapes_for_systematic(bool
515		signalStatDn->SetTitle(((string)signal->GetTitle()+"_StatDown").c_str());
516
517		for(int i=1;i<=signalStatDn->GetNbinsX();i++) {
518	<	float staterr = TMath::Sqrt(Lpred->GetBinContent(i) + Lobs->GetBinContent(i));
519	<	signalStatDn->SetBinContent(i,signal->GetBinContent(i)-staterr);
518	>	float staterr;
519	>	if(!PlottingSetup::FullMCAnalysis) staterr = TMath::Sqrt(Lpred->GetBinContent(i) + Lobs->GetBinContent(i));
520	>	else staterr = TMath::Sqrt(Lobs->GetBinContent(i));
521	>
522	>	if(!PlottingSetup::FullMCAnalysis) {
523	>	dout << "Stat err in bin " << i << " : " << staterr << endl;
524	>	dout << "    prediction: " << Lpred->GetBinContent(i) << " , observation: " << Lobs->GetBinContent(i) << " --> signal: " << signal->GetBinContent(i) << endl;
525	>	dout << "    we obtain : " << signal->GetBinContent(i)-staterr << " , " << signal->GetBinContent(i)+staterr << endl;
526	>	}
527	>	if(signal->GetBinContent(i)-staterr>0) signalStatDn->SetBinContent(i,signal->GetBinContent(i)-staterr);
528	>	else signalStatDn->SetBinContent(i,0);
529		signalStatUp->SetBinContent(i,signal->GetBinContent(i)+staterr);
530		signal->SetBinError(i,staterr);
531		}
#	Line 292 | Line 545 | void generate_shapes_for_systematic(bool
545		flippedsignalStatDn->SetTitle(((string)flippedsignal->GetTitle()+"_StatDown").c_str());
546
547		for(int i=1;i<=flippedsignalStatDn->GetNbinsX();i++) {
548	<	float staterr = TMath::Sqrt(flippedLpred->GetBinContent(i) + flippedLobs->GetBinContent(i));
549	<	flippedsignalStatDn->SetBinContent(i,flippedsignal->GetBinContent(i)-staterr);
548	>	float staterr;
549	>	if(!PlottingSetup::FullMCAnalysis) staterr = TMath::Sqrt(flippedLpred->GetBinContent(i) + flippedLobs->GetBinContent(i));
550	>	else staterr = TMath::Sqrt(flippedLobs->GetBinContent(i));
551	>	if(flippedsignal->GetBinContent(i)-staterr>0) flippedsignalStatDn->SetBinContent(i,flippedsignal->GetBinContent(i)-staterr);
552	>	else flippedsignalStatDn->SetBinContent(i,0);
553		flippedsignalStatUp->SetBinContent(i,flippedsignal->GetBinContent(i)+staterr);
554		flippedsignal->SetBinError(i,staterr);
555		}
#	Line 325 | Line 581 | void generate_shapes_for_systematic(bool
581
582
583		if(dataonly) {
584	<	cout << "Processing data with datajzb: " << datajzb << endl;
584	>	dout << "Processing data with datajzb: " << datajzb << endl;
585		TH1F *ZOSSFP = allsamples.Draw("ZOSSFP",datajzb,binning, "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,data,luminosity);
586		TH1F *ZOSOFP = allsamples.Draw("ZOSOFP",datajzb,binning, "JZB", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,data,luminosity);
587		TH1F *ZOSSFN = allsamples.Draw("ZOSSFN","-"+datajzb,binning, "JZB", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,data,luminosity);
#	Line 459 | Line 715 | void generate_shapes_for_systematic(bool
715		SQRT(predstaterr);
716		TH1F *bgStatUp = (TH1F*)pred->Clone("background_StatUp");
717		bgStatUp->Add(predstaterr);
718	+	EliminateNegativeEntries(bgStatUp);
719		bgStatUp->Write();
720		TH1F *bgStatDn = (TH1F*)pred->Clone("background_StatDown");
721		bgStatDn->Add(predstaterr,-1);
722	+	EliminateNegativeEntries(bgStatDn);
723		bgStatDn->Write();
724		//      delete bgStatDn;
725		//      delete bgStatUp;
#	Line 478 | Line 736 | void generate_shapes_for_systematic(bool
736		SQRT(flippedpredstaterr);
737		TH1F *fbgStatUp = (TH1F*)flippedpred->Clone("flipped_background_StatUp");
738		fbgStatUp->Add(predstaterr);
739	+	EliminateNegativeEntries(fbgStatUp);
740		fbgStatUp->Write();
741		TH1F *fbgStatDn = (TH1F*)flippedpred->Clone("flipped_background_StatDown");
742		fbgStatDn->Add(predstaterr,-1);
743	+	EliminateNegativeEntries(fbgStatDn);
744		fbgStatDn->Write();
745		//      delete fbgStatDn;
746		//      delete fbgStatUp;
#	Line 493 | Line 753 | void generate_shapes_for_systematic(bool
753		SQRT(Tpredstaterr);
754		TH1F *TpredStatUp = (TH1F*)Tpred->Clone("TTbarBackground_StatUp");
755		TpredStatUp->Add(Tpredstaterr);
756	+	EliminateNegativeEntries(TpredStatUp);
757		TpredStatUp->Write();
758		TH1F *TpredStatDn = (TH1F*)Tpred->Clone("TTbarBackground_StatDown");
759		TpredStatDn->Add(Tpredstaterr,-1);
760	+	EliminateNegativeEntries(TpredStatDn);
761		TpredStatDn->Write();
762		//      delete TpredStatDn;
763		//      delete TpredStatUp;
#	Line 508 | Line 770 | void generate_shapes_for_systematic(bool
770		SQRT(fTpredstaterr);
771		TH1F *fTpredStatUp = (TH1F*)flippedTpred->Clone("flipped_TTbarBackground_StatUp");
772		fTpredStatUp->Add(fTpredstaterr);
773	+	EliminateNegativeEntries(fTpredStatUp);
774		fTpredStatUp->Write();
775		TH1F *fTpredStatDn = (TH1F*)flippedTpred->Clone("flipped_TTbarBackground_StatDown");
776		fTpredStatDn->Add(fTpredstaterr,-1);
777	+	EliminateNegativeEntries(fTpredStatDn);
778		fTpredStatDn->Write();
779		//      delete fTpredStatDn;
780		//      delete fTpredStatUp;
#	Line 524 | Line 788 | void generate_shapes_for_systematic(bool
788		SQRT(Zpredstaterr);
789		TH1F *ZpredStatUp = (TH1F*)Zpred->Clone("ZJetsBackground_StatUp");
790		ZpredStatUp->Add(Zpredstaterr);
791	+	EliminateNegativeEntries(ZpredStatUp);
792		ZpredStatUp->Write();
793		TH1F *ZpredStatDn = (TH1F*)Zpred->Clone("ZJetsBackground_StatDown");
794		ZpredStatDn->Add(Zpredstaterr,-1);
795	+	EliminateNegativeEntries(ZpredStatDn);
796		ZpredStatDn->Write();
797		//      delete ZpredStatDn;
798		//      delete ZpredStatUp;
#	Line 539 | Line 805 | void generate_shapes_for_systematic(bool
805		SQRT(fTpredstaterr);
806		TH1F *fZpredStatUp = (TH1F*)flippedZpred->Clone("flipped_ZJetsBackground_StatUp");
807		fZpredStatUp->Add(fZpredstaterr);
808	+	EliminateNegativeEntries(fZpredStatUp);
809		fZpredStatUp->Write();
810		TH1F *fZpredStatDn = (TH1F*)flippedZpred->Clone("flipped_ZJetsBackground_StatDown");
811		fZpredStatDn->Add(fZpredstaterr,-1);
812	+	EliminateNegativeEntries(fZpredStatDn);
813		fZpredStatDn->Write();
814		//      delete fZpredStatDn;
815		//      delete fZpredStatUp;
#	Line 561 | Line 829 | void generate_shapes_for_systematic(bool
829		SQRT(predsyserr);
830		TH1F *bgSysUp = (TH1F*)pred->Clone("background_SysUp");
831		bgSysUp->Add(predsyserr);
832	+	EliminateNegativeEntries(bgSysUp);
833		bgSysUp->Write();
834		TH1F *bgSysDn = (TH1F*)pred->Clone("background_SysDown");
835		bgSysDn->Add(predsyserr,-1);
836	+	EliminateNegativeEntries(bgSysDn);
837		bgSysDn->Write();
838		delete predsyserr;
839
#	Line 579 | Line 849 | void generate_shapes_for_systematic(bool
849		SQRT(fpredsyserr);
850		TH1F *fbgSysUp = (TH1F*)flippedpred->Clone("flipped_background_SysUp");
851		fbgSysUp->Add(fpredsyserr);
852	+	EliminateNegativeEntries(fbgSysUp);
853		fbgSysUp->Write();
854		TH1F *fbgSysDn = (TH1F*)flippedpred->Clone("flipped_background_SysDown");
855		fbgSysDn->Add(fpredsyserr,-1);
856	+	EliminateNegativeEntries(fbgSysDn);
857		fbgSysDn->Write();
858		delete fpredsyserr;
859
#	Line 594 | Line 866 | void generate_shapes_for_systematic(bool
866		SQRT(Tpredsyserr);
867		TH1F *TpredSysUp = (TH1F*)Tpred->Clone("TTbarBackground_SysUp");
868		TpredSysUp->Add(Tpredsyserr);
869	+	EliminateNegativeEntries(TpredSysUp);
870		TpredSysUp->Write();
871		TH1F *TpredSysDn = (TH1F*)Tpred->Clone("TTbarBackground_SysDown");
872		TpredSysDn->Add(Tpredsyserr,-1);
873	+	EliminateNegativeEntries(TpredSysDn);
874		TpredSysDn->Write();
875		delete Tpredsyserr;
876
#	Line 608 | Line 882 | void generate_shapes_for_systematic(bool
882		SQRT(fTpredsyserr);
883		TH1F *fTpredSysUp = (TH1F*)flippedTpred->Clone("flipped_TTbarBackground_SysUp");
884		fTpredSysUp->Add(fTpredsyserr);
885	+	EliminateNegativeEntries(fTpredSysUp);
886		fTpredSysUp->Write();
887		TH1F *fTpredSysDn = (TH1F*)flippedTpred->Clone("flipped_TTbarBackground_SysDown");
888		fTpredSysDn->Add(fTpredsyserr,-1);
889	+	EliminateNegativeEntries(fTpredSysDn);
890		fTpredSysDn->Write();
891		delete fTpredsyserr;
892
#	Line 625 | Line 901 | void generate_shapes_for_systematic(bool
901		SQRT(Zpredsyserr);
902		TH1F *ZpredSysUp = (TH1F*)Zpred->Clone("ZJetsBackground_SysUp");
903		ZpredSysUp->Add(Zpredsyserr);
904	+	EliminateNegativeEntries(ZpredSysUp);
905		ZpredSysUp->Write();
906		TH1F *ZpredSysDn = (TH1F*)Zpred->Clone("ZJetsBackground_SysDown");
907		ZpredSysDn->Add(Zpredsyserr,-1);
908	+	EliminateNegativeEntries(ZpredSysDn);
909		ZpredSysDn->Write();
910		delete Zpredsyserr;
911
#	Line 640 | Line 918 | void generate_shapes_for_systematic(bool
918		SQRT(fZpredsyserr);
919		TH1F *fZpredSysUp = (TH1F*)flippedZpred->Clone("flipped_ZJetsBackground_SysUp");
920		fZpredSysUp->Add(fZpredsyserr);
921	+	EliminateNegativeEntries(fZpredSysUp);
922		fZpredSysUp->Write();
923		TH1F *fZpredSysDn = (TH1F*)flippedZpred->Clone("flipped_ZJetsBackground_SysDown");
924		fZpredSysDn->Add(fZpredsyserr,-1);
925	+	EliminateNegativeEntries(fZpredSysDn);
926		fZpredSysDn->Write();
927		delete fZpredsyserr;
928		}
929
930		/*if(identifier=="") {
931		for(int i=0;i<binning.size()-1;i++) {
932	<	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;
932	>	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;
933		}
934		}*/
935		delete ZOSSFP;
#	Line 667 | Line 947 | void generate_shapes_for_systematic(bool
947
948		}
949
950	<	ShapeDroplet LimitsFromShapes(TTree *events, string addcut, string name, string mcjzb, string datajzb, vector<float> jzbbins, float jzbpeakerrormc, float jzbpeakerrordata) {
950	>	void DoFullCls(string RunDirectory,float firstGuess) {
951	>	vector<float> epoints = get_expected_points(firstGuess,25,RunDirectory,true);//get 20 points;
952	>	ofstream RealScript;
953	>	dout << "Going to write the full CLs script to " << RunDirectory << "/LimitScript.sh" << endl;
954	>	RealScript.open ((RunDirectory+"/LimitScript.sh").c_str());
955	>	RealScript << "#!/bin/bash \n";
956	>	RealScript << "\n";
957	>	RealScript << "RunDirectory=" << RunDirectory << "\n";
958	>	RealScript << "CMSSWDirectory=" << PlottingSetup::CMSSW_dir << "\n";
959	>	RealScript << "\n";
960	>	RealScript << "echo \"Going to compute limit in $RunDirectory\"\n";
961	>	RealScript << "ORIGDIR=`pwd`\n";
962	>	RealScript << "\n";
963	>	RealScript << "pushd $CMSSWDirectory\n";
964	>	RealScript << "cd ~/final_production_2011/CMSSW_4_2_8/src/HiggsAnalysis/\n";
965	>	RealScript << "origscramarch=$SCRAM_ARCH\n";
966	>	RealScript << "origbase=$CMSSW_BASE\n";
967	>	RealScript << "export SCRAM_ARCH=slc5_amd64_gcc434\n";
968	>	RealScript << "eval `scram ru -sh`\n";
969	>	RealScript << "\n";
970	>	RealScript << "mkdir -p $RunDirectory\n";
971	>	RealScript << "cd $RunDirectory\n";
972	>	RealScript << "echo `pwd`\n";
973	>	RealScript << "mkdir -p FullCLs\n";
974	>	RealScript << "cd FullCLs\n";
975	>	RealScript << "\n";
976	>	RealScript << "combineEXE=\"${CMSSW_BASE}/bin/${SCRAM_ARCH}/combine\"\n";
977	>	RealScript << "\n";
978	>	RealScript << "time for i in {";
979	>	RealScript << epoints[0]/4 << ","; // adding a VERY VERY low point just to be totally safe
980	>	RealScript << epoints[0]/2 << ","; // adding a VERY low point just to be safe
981	>	for(int i=0;i<epoints.size();i++) RealScript << epoints[i] << ",";
982	>	RealScript << 2*epoints[epoints.size()-1] << ","; // adding a VERY high point just to be safe
983	>	RealScript << 4*epoints[epoints.size()-1] << ","; // adding a VERY VERY high point just to be totally safe
984	>	RealScript << "}; do time $combineEXE ../susydatacard.txt -M HybridNew --freq -T 500 --clsAcc 0 -v 0 -n TestPoint --saveHybridResult --saveToys -s -1 -i 8 --singlePoint $i; done\n"; // 500 toys, epoints.size+2 testpoints, 8 iterations
985	>	RealScript << "\n";
986	>	RealScript << "hadd FullCLs.root higgsCombine*.root\n";
987	>	RealScript << "mkdir originalFiles\n";
988	>	RealScript << "mv higgsCombine* originalFiles\n";
989	>	RealScript << "$combineEXE ../susydatacard.txt -M HybridNew --freq --grid=FullCLs.root\n";
990	>	RealScript << "$combineEXE ../susydatacard.txt -M HybridNew --freq --grid=FullCLs.root --expectedFromGrid 0.5\n";
991	>	RealScript << "$combineEXE ../susydatacard.txt -M HybridNew --freq --grid=FullCLs.root --expectedFromGrid 0.84\n";
992	>	RealScript << "$combineEXE ../susydatacard.txt -M HybridNew --freq --grid=FullCLs.root --expectedFromGrid 0.16\n";
993	>	RealScript << "$combineEXE ../susydatacard.txt -M HybridNew --freq --grid=FullCLs.root --expectedFromGrid 0.975\n";
994	>	RealScript << "$combineEXE ../susydatacard.txt -M HybridNew --freq --grid=FullCLs.root --expectedFromGrid 0.025\n";
995	>	RealScript << "\n";
996	>	RealScript << "hadd FinalResult.root higgsCombineTest.HybridNew*\n";
997	>	RealScript << "\n";
998	>	RealScript << "g++ $ORIGDIR/ReadAndSave.C -o ReadAndSave.exec `root-config --glibs --cflags` \n";
999	>	RealScript << "./ReadAndSave.exec FinalResult.root " << RunDirectory << "/FullCLsShapeDropletResult.txt \n";
1000	>	RealScript << "\n";
1001	>	RealScript << "#####\n";
1002	>	RealScript << "echo \"Finalizing ...\"\n";
1003	>	RealScript << "cd $ORIGDIR\n";
1004	>	RealScript << "echo $ORIGDIR\n";
1005	>	RealScript << "ls -ltrh ../SandBox/ | grep combine\n";
1006	>	RealScript << "export SCRAM_ARCH=$origscramarch\n";
1007	>	RealScript << "exit 0\n";
1008	>	RealScript << "\n";
1009	>	RealScript.close();
1010	>	gSystem->Exec(("bash "+RunDirectory+"/LimitScript.sh").c_str());
1011	>	}
1012	>
1013	>
1014	>	ShapeDroplet LimitsFromShapes(bool asymptotic, float firstGuess, TTree *events, string addcut, string name, string mcjzb, string datajzb, vector<float> jzbbins, float jzbpeakerrormc, float jzbpeakerrordata) {
1015		map <  pair<float, float>, map<string, float>  >  xsec;
1016		if(SUSYScanSpace::SUSYscantype==mSUGRA) xsec=getXsec(PlottingSetup::cbafbasedir+"/"+mSUGRAxsFile);
1017
674	–	write_info(__FUNCTION__,"Implementing the shape function");
1018		time_t timestamp;
1019		tm *now;
1020		timestamp = time(0);
#	Line 682 | Line 1025 | ShapeDroplet LimitsFromShapes(TTree *eve
1025
1026		ensure_directory_exists(RunDirectory);
1027
1028	<	TFile *datafile = new TFile("../StoredShapes.root","READ");
1028	>	TFile *datafile = new TFile((PlottingSetup::cbafbasedir+"/DistributedModelCalculations/StoredShapes.root").c_str(),"READ");
1029		if(datafile->IsZombie()) {
1030		write_error(__FUNCTION__,"Fatal error: The stored shapes are not available!");
1031		assert(!datafile->IsZombie());
1032		}
1033	<	cout << "Run Directory: " << RunDirectory << endl;
1033	>	dout << "Run Directory: " << RunDirectory << endl;
1034		TFile *limfile = new TFile((RunDirectory+"/PRElimitfile.root").c_str(),"RECREATE");
1035
1036		TIter nextkey(datafile->GetListOfKeys());
#	Line 705 | Line 1048 | ShapeDroplet LimitsFromShapes(TTree *eve
1048		bool dataonly=false;
1049
1050		generate_shapes_for_systematic(signalonly,dataonly,limfile,events,"",mcjzb,datajzb,noJES,jzbbins,limcan,addcut,xsec);
708	–	//  generate_shapes_for_systematic(signalonly,dataonly,limfile,events,"peakUp",newjzbexpression(mcjzb,jzbpeakerrormc),newjzbexpression(datajzb,jzbpeakerrordata),noJES,jzbbins,limcan,addcut,xsec);
709	–	//  generate_shapes_for_systematic(signalonly,dataonly,limfile,events,"peakDown",newjzbexpression(mcjzb,-jzbpeakerrormc),newjzbexpression(datajzb,-jzbpeakerrordata),noJES,jzbbins,limcan,addcut,xsec);
1051		generate_shapes_for_systematic(signalonly,dataonly,limfile,events,"JESUp",mcjzb,datajzb,JESup,jzbbins,limcan,addcut,xsec);
1052		generate_shapes_for_systematic(signalonly,dataonly,limfile,events,"JESDown",mcjzb,datajzb,JESdown,jzbbins,limcan,addcut,xsec);
1053
#	Line 720 | Line 1061 | ShapeDroplet LimitsFromShapes(TTree *eve
1061		bool docomplicatedmSUGRAxsreweighting=false; //if you modify this value please also adapt it in Systematics.C not only here in ShapeLimit.C
1062
1063		MCefficiency(events,mceff,mcefferr,flipped,mcjzb,requireZ,Neventsinfile,SUSYScanSpace::SUSYscantype,xsec,addcut,-1);
1064	<	if(mceff<0) flipped=1;
1064	>	if(mceff<0) {
1065	>	flipped=1;
1066	>	write_info(__FUNCTION__,"Doing flipping!");
1067	>	}
1068		doJZBscale(events,SUSYScanSpace::SUSYscantype==mSUGRA&&docomplicatedmSUGRAxsreweighting,xsec,scaledown,scaleup,scalesyst,JZBScaleUncert,informalname,flipped,requireZ,addcut);
1069	<	float PDFuncert;
1069	>	float PDFuncert=0;
1070		int NPdfs = get_npdfs(events);
1071	<	if(SUSYScanSpace::SUSYscantype!=mSUGRA) PDFuncert = get_pdf_uncertainty(events, flipped, mcjzb, requireZ, Neventsinfile, SUSYScanSpace::SUSYscantype, xsec, NPdfs, addcut);
1072	<
1071	>	write_warning(__FUNCTION__,"Deactivated PDFs temporarily for exp limits!!!");
1072	>	//  if(SUSYScanSpace::SUSYscantype!=mSUGRA) PDFuncert = get_pdf_uncertainty(events, flipped, mcjzb, requireZ, Neventsinfile, SUSYScanSpace::SUSYscantype, xsec, NPdfs, addcut);
1073
1074		float JZBscale=scaledown;
1075		if(scaleup>scaledown) JZBscale=scaleup;
#	Line 735 | Line 1079 | ShapeDroplet LimitsFromShapes(TTree *eve
1079		dout << "       JZB Scale (JSU) : " << JZBscale << endl;
1080		dout << "       PDF             : " << PDFuncert << endl;
1081
1082	<	prepare_limit_datacard(RunDirectory,limfile,JZBscale,PDFuncert,flipped);
1082	>	float SignalIntegral;
1083	>	if(flipped) {
1084	>	TH1F *flisi = (TH1F*)(limfile->Get("flipped_signal"));
1085	>	SignalIntegral= flisi ->Integral();
1086	>	delete flisi;
1087	>	}
1088	>	else {
1089	>	TH1F *flisi = (TH1F*)(limfile->Get("signal"));
1090	>	SignalIntegral= flisi ->Integral();
1091	>	delete flisi;
1092	>	}
1093	>
1094		TFile *final_limfile = new TFile((RunDirectory+"/limitfile.root").c_str(),"RECREATE");
1095
1096		TIter nnextkey(limfile->GetListOfKeys());
1097		TKey *nkey;
1098	+
1099	+	if(SignalIntegral==0) SignalIntegral=1; //case when the integral is zero - don't try to normalize, just leave it as it is
1100	+
1101		while ((nkey = (TKey*)nnextkey()))
1102		{
1103		TH1F *obj =(TH1F*) nkey->ReadObj();
#	Line 747 | Line 1105 | ShapeDroplet LimitsFromShapes(TTree *eve
1105		if(!flipped&&Contains(obj->GetName(),"flipped")) continue;
1106		if(flipped) obj->SetName(strip_flip_away(obj->GetName()));
1107		final_limfile->cd();
1108	+	if(Contains(obj->GetName(),"signal")) {
1109	+	obj->Scale(1.0/SignalIntegral);
1110	+	}
1111		obj->Write();
1112		}
1113
1114	+	prepare_limit_datacard(RunDirectory,final_limfile,JZBscale,PDFuncert);
1115	+
1116		final_limfile->Close();
1117		limfile->Close();
1118	<	write_info(__FUNCTION__,"Shape root file and datacard have been generated in "+RunDirectory);
1119	<
1120	<	int CreatedModelFileExitCode = gSystem->Exec(("bash CreateModel.sh "+RunDirectory+" susydatacard.txt").c_str());
1121	<	cout << "exit code of limit algorithm (CreateModel.sh) : " << CreatedModelFileExitCode << endl;
1122	<	assert(CreatedModelFileExitCode==0);
1118	>	dout << "Info: Shape root file and datacard have been generated in " << RunDirectory << endl;
1119	>	stringstream command;
1120	>	string DropletLocation;
1121	>	int CreatedModelFileExitCode;
1122	>	if(asymptotic) {
1123	>	if(firstGuess>0) {
1124	>	command << "bash " << PlottingSetup::cbafbasedir<< "/DistributedModelCalculations/ShapeLimits/CreateModel.sh " << RunDirectory << " susydatacard.txt" << " 1 0 0 " << firstGuess; // ASYMPTOTIC LIMITS WITH FIRST GUESS
1125	>	} else {
1126	>	command << "bash " << PlottingSetup::cbafbasedir<< "/DistributedModelCalculations/ShapeLimits/CreateModel.sh " << RunDirectory << " susydatacard.txt" << " 0 0 0 0"; // ASYMPTOTIC LIMITS
1127	>	}
1128	>	dout <<"Going to run : " << command.str() << endl;
1129	>	CreatedModelFileExitCode = gSystem->Exec(command.str().c_str());
1130	>	dout << "exit code of limit algorithm (CreateModel.sh) : " << CreatedModelFileExitCode << endl;
1131	>	if(!(CreatedModelFileExitCode==0)) {
1132	>	write_warning(__FUNCTION__,"Something bad happened. It looks like a shape analysis is not the way to go. ");
1133	>	ShapeDroplet alpha;
1134	>	alpha.observed=-12345; // this is the flag to say watch out something went wrong with the signal ...
1135	>	alpha.SignalIntegral=1;
1136	>	return alpha;
1137	>	}
1138	>	DropletLocation=RunDirectory+"/ShapeDropletResult.txt";
1139	>	}
1140	>	else {
1141	>	DoFullCls(RunDirectory,firstGuess);
1142	>	DropletLocation=RunDirectory+"/FullCLsShapeDropletResult.txt";
1143	>	}
1144		ShapeDroplet alpha;
1145	<	alpha.readDroplet(RunDirectory+"/ShapeDropletResult.txt");
1145	>	alpha.readDroplet(DropletLocation);
1146		alpha.PDF=PDFuncert;
1147		alpha.JSU=JZBscale;
1148	+	alpha.SignalIntegral=SignalIntegral;
1149		dout << "Done reading limit information from droplet" << endl;
1150		dout << alpha << endl;
1151	<
1152	<	write_warning(__FUNCTION__,"Once you're done implementing the limit computation, please uncomment the following line to enable clean up.");
1151	>	if(asymptotic) {
1152	>	dout << "Doing asymptotic limits, therefore adding an additional round!" << endl;
1153	>	command.str("");
1154	>	command << "bash " << PlottingSetup::cbafbasedir<< "/DistributedModelCalculations/ShapeLimits/CreateModel.sh " << RunDirectory << " susydatacard.txt" << " 1 0 0 " << alpha.expected; // ASYMPTOTIC LIMITS WITH FIRST GUESS
1155	>	CreatedModelFileExitCode = gSystem->Exec(command.str().c_str());
1156	>	dout << "exit code of limit algorithm (CreateModel.sh) : " << CreatedModelFileExitCode << endl;
1157	>	if(!(CreatedModelFileExitCode==0)) {
1158	>	write_warning(__FUNCTION__,"Something bad happened. It looks like a shape analysis is not the way to go. ");
1159	>	ShapeDroplet alpha;
1160	>	alpha.observed=-12345; // this is the flag to say watch out something went wrong with the signal ...
1161	>	alpha.SignalIntegral=1;
1162	>	return alpha;
1163	>	}
1164	>	alpha.readDroplet(RunDirectory+"/ShapeDropletResult.txt");
1165	>	alpha.PDF=PDFuncert;
1166	>	alpha.JSU=JZBscale;
1167	>	alpha.SignalIntegral=SignalIntegral;
1168	>	dout << "Done reading updated limit information from droplet" << endl;
1169	>	dout << alpha << endl;
1170	>	}
1171	>
1172		dout << "Everything is saved in " << RunDirectory << endl;
1173	<	//  gSystem->Exec(("rm -r "+RunDirectory).c_str());
1174	<
1175	<	return alpha;
1173	>	dout << "Cleaning up ... " << std::flush;
1174	>	/*  dout << "   1) Make sure models directory exists ... " << std::flush;
1175	>	gSystem->Exec("mkdir -p models/");
1176	>	dout << " ok!" << endl;
1177	>	dout << "   2) Deleting any previous model files with the same name ... " << std::flush;
1178	>	stringstream delcommand;
1179	>	delcommand << "rm models/model_" << name << ".root";
1180	>	gSystem->Exec(delcommand.str().c_str());
1181	>	stringstream delcommand2;
1182	>	delcommand2 << "rm models/model_" << name << ".txt";
1183	>	gSystem->Exec(delcommand2.str().c_str());
1184	>	dout << " ok!" << endl;
1185	>	dout << "   3) Transfering the new model files ... " << std::flush;
1186	>	stringstream copycommand;
1187	>	copycommand << "cp " << RunDirectory << "/model.root models/model_" << name << ".root";
1188	>	gSystem->Exec(copycommand.str().c_str());
1189	>	stringstream copycommand2;
1190	>	copycommand2 << "cp " << RunDirectory << "/susydatacard.txt models/model_" << name << ".txt";
1191	>	gSystem->Exec(copycommand2.str().c_str());
1192	>	stringstream copycommand3;
1193	>	copycommand3 << "cp " << RunDirectory << "/limitfile.root models/model_" << name << "_histo.root";
1194	>	gSystem->Exec(copycommand3.str().c_str());
1195	>	dout << " ok!" << endl;
1196	>	dout << "   4) Removing original working directory (" << RunDirectory << ") ... " << std::flush;*/
1197	>	gSystem->Exec(("rm -r "+RunDirectory).c_str());
1198	>	dout << " ok!" << endl;
1199	>	delete limcan;
1200	>	return alpha;
1201		}
1202
1203	<	void PrepareDataShapes(string datajzb, vector<float> jzbbins, float jzbpeakerrordata) {
1204	<	//  dout << "Generating data shape templates ... " << endl;
1203	>	void PrepareDataShapes(string mcjzb, string datajzb, vector<float> jzbbins, float jzbpeakerrordata) {
1204	>
1205		map <  pair<float, float>, map<string, float>  >  xsec;
1206		TFile *datafile = new TFile("../DistributedModelCalculations/StoredShapes.root","RECREATE");
1207		TCanvas *limcan = new TCanvas("limcan","Canvas for calculating limits");
1208		TTree *faketree;
1209		bool dataonly=true;
1210		bool signalonly=false;
1211	<	string mcjzb="JZBforMCinPrepareDataShapes";//this string is not used.
783	<	float jzbpeakerrormc=0;
1211	>
1212		generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"",mcjzb,datajzb,noJES,jzbbins,limcan,"",xsec);
1213	<	// don't need these effects for obs & pred, only for signal!
1214	<	//  generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"peakUp",newjzbexpression(mcjzb,jzbpeakerrormc),newjzbexpression(datajzb,jzbpeakerrordata),noJES,jzbbins,limcan,"",xsec);
1215	<	//  generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"peakDown",newjzbexpression(mcjzb,-jzbpeakerrormc),newjzbexpression(datajzb,-jzbpeakerrordata),noJES,jzbbins,limcan,"",xsec);
1216	<	//  generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"JESUp",mcjzb,datajzb,JESup,jzbbins,limcan,"",xsec);
1217	<	//  generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"JESDown",mcjzb,datajzb,JESdown,jzbbins,limcan,"",xsec);
1213	>
1214	>	generate_mc_shapes(false,faketree,mcjzb,datajzb,jzbbins,cutWeight,cutnJets,"","",xsec);
1215	>	generate_mc_shapes(false,faketree,mcjzb,datajzb,jzbbins,cutWeight,cutnJetsJESup,"","JESUp",xsec);
1216	>	generate_mc_shapes(false,faketree,mcjzb,datajzb,jzbbins,cutWeight,cutnJetsJESdown,"","JESDown",xsec);
1217	>	generate_mc_shapes(false,faketree,mcjzb,datajzb,jzbbins,TCut("weightEffUp"),cutnJets,"","EffUp",xsec);
1218	>	generate_mc_shapes(false,faketree,mcjzb,datajzb,jzbbins,TCut("weightEffDown"),cutnJets,"","EffDown",xsec);
1219	>	generate_mc_shapes(false,faketree,mcjzb,datajzb,jzbbins,cutWeight,cutnJets,"","StatUp",xsec);
1220	>	generate_mc_shapes(false,faketree,mcjzb,datajzb,jzbbins,cutWeight,cutnJets,"","StatDown",xsec);
1221	>	/*
1222		datafile->Close();
1223	+	*/
1224		}
1225

Diff Legend

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