ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

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

Comparing UserCode/cbrown/Development/Plotting/Modules/SUSYScan.C (file contents):
Revision 1.3 by buchmann, Thu Mar 15 20:34:15 2012 UTC vs.
Revision 1.17 by buchmann, Thu May 17 19:55:38 2012 UTC

#	Line 4 | Line 4
4		#include <fstream>
5		#include <pthread.h>
6		#include <assert.h>
7	+	#include <cctype>
8
9		#include <TCut.h>
10		#include <TROOT.h>
#	Line 19 | Line 20
20		#include <TProfile.h>
21		#include <TKey.h>
22
23	+	#include "ShapeDroplet.C"
24	+	#include "ShapeLimit.C"
25	+
26		//#include "TTbar_stuff.C"
27		using namespace std;
28
29		using namespace PlottingSetup;
30
31	+	void scan_SUSY_parameter_space(string mcjzb,string datajzb,vector<float> jzb_cut,bool requireZ, float peakerror, float peakerrordata, float njobs, float jobnumber, bool systonly, bool effonly,bool shapeanalysis);
32	+
33	+	bool isValidShapeResult(ShapeDroplet droplet) {
34	+	if(droplet.isValid==false) return false;
35	+	if(droplet.observed>50000) {
36	+	// error code of ShapeLimit/CreateModel.sh, invoked if there were too many errors during computation
37	+	write_warning(__FUNCTION__,"Shape limit was ridiculously high ("+any2string(droplet.observed)+"); this is a typical error signature of the combine tool. Shape result will be discarded, and c&c will be carried out instead.");
38	+	return false;
39	+	}
40	+	if(droplet.observed==-12345) {
41	+	write_warning(__FUNCTION__,"Shape limit was ("+any2string(droplet.observed)+"); this is the error signature of CreateModel.sh / ShapeDroplet. Shape result will be discarded, and c&c will be carried out instead.");
42	+	return false;
43	+	}
44	+	if(droplet.expected>-1&&droplet.expected==0&&droplet.observed>5000) {
45	+	write_warning(__FUNCTION__,"Observed shape limit was ("+any2string(droplet.observed)+") while expected limit was "+any2string(droplet.expected)+"; this is a typical error signature of the combine tool. Shape result will be discarded, and c&c will be carried out instead.");
46	+	return false;
47	+	}
48	+	if(droplet.expected>-1&&droplet.expected/droplet.observed>50) {
49	+	write_warning(__FUNCTION__,"Observed shape limit was ("+any2string(droplet.observed)+") while expected limit was "+any2string(droplet.expected)+"; such a high ratio between expected & observed limit is a typical error signature of the combine tool. Shape result will be discarded, and c&c will be carried out instead.");
50	+	return false;
51	+	}
52	+	if(droplet.expected>-1&&droplet.observed/droplet.expected>50) {
53	+	write_warning(__FUNCTION__,"Observed shape limit was ("+any2string(droplet.observed)+") while expected limit was "+any2string(droplet.expected)+"; such a high ratio between observed & expected limit is a typical error signature of the combine tool. Shape result will be discarded, and c&c will be carried out instead.");
54	+	return false;
55	+	}
56	+	return true;
57	+	}
58	+
59		void susy_scan_axis_labeling(TH2F *histo) {
60		histo->GetXaxis()->SetTitle("m_{#Chi_{2}^{0}}-m_{LSP}");
61		histo->GetXaxis()->CenterTitle();
#	Line 61 | Line 93 | void set_SUSY_style() {
93
94		gStyle->SetPalette(100, MyPalette);
95
64	–	float rightmargin=gStyle->GetPadRightMargin();
96		gStyle->SetPadRightMargin(0.15);
97
98		}
99
100		bool delete_any_cached_scans() {
70	–	//This should only be called via CRAB
71	–	cout << "   Deleting all cached files." << endl;
72	–	gSystem->Exec("ls -ltrh | grep root");
101		char hostname[1023];
102		gethostname(hostname,1023);
103	<	if((Contains(hostname,"t3ui")||Contains(hostname,"t3wn"))) return false;
104	<	else {
103	>
104	>
105	>	/*
106	>	//This should only be called via CRAB
107	>	dout << "   Deleting all cached files." << endl;
108	>	gSystem->Exec("ls -ltrh | grep root");
109	>	if(!((Contains(hostname,"t3ui")||Contains(hostname,"t3wn"))))
110	>	{
111		vector<string> all_files;
112		char currentpath[1024];
113		TString directory=TString(getcwd(currentpath,1024));
#	Line 86 | Line 120 | bool delete_any_cached_scans() {
120		gSystem->Exec(((string)"rm "+all_files[ifile]).c_str());
121		}
122		}
123	<	return true;
123	>	return true;
124		}
125	+	*/
126	+	if(!((Contains(hostname,"t3ui")||Contains(hostname,"t3wn")))) {
127	+	dout << "Going to purge files in local_storage that have been loaded!" << endl;
128	+	for(int i=0;i<(int)SUSYScanSpace::loaded_files.size();i++) {
129	+	gSystem->Exec(((string)"rm "+SUSYScanSpace::loaded_files[i]).c_str());
130	+	dout << "      Purging : Deleted file " << SUSYScanSpace::loaded_files[i] << endl;
131	+	}
132	+	}
133	+	return true;
134		}
135
136		bool initialized_t2=false;
137
138		int srmcpretries=0;
139
140	<	void load_scan_sample(int a, int b, int &scanfileindex,int scantype,bool isretry=false) {
141	<	/*
142	<	// There is no need to define your sample here. That is done in Setup.C where you define the loading directory!
100	<
101	<	dout << "Going to load file with Xzone=" << a << " and  Yzone=" << b << endl;
102	<	stringstream filetoload;
103	<	char hostname[1023];
104	<	gethostname(hostname,1023);
105	<
106	<	if((Contains(hostname,"t3ui")||Contains(hostname,"t3wn"))) { // local case
140	>	void load_local_scan_sample(int a, int b, int &scanfileindex,int scantype,int mglu, int mlsp, bool isretry=false) {
141	>	stringstream filetoload;
142	>	string samplename;
143		filetoload << "/shome/buchmann/ntuples/"<<PlottingSetup::ScanSampleDirectory;
144	<	if(scantype==mSUGRA) filetoload << "/mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
145	<	if(scantype==SMS) filetoload << "/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
146	<	if(scantype==GMSB) filetoload << "/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
144	>	if(scantype==mSUGRA) {
145	>	//filetoload << "/mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
146	>	filetoload.str("");
147	>	filetoload << "/shome/lbaeni/jzb/" << PlottingSetup::ScanSampleDirectory  << "/mSUGRA_M0_" << mglu << "__M12_" << mlsp << ".root";
148	>	samplename="mSUGRA";
149	>	//      filetoload << "/mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
150	>	}
151	>	if(scantype==SMS) {
152	>	filetoload << "/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
153	>	samplename="SMS";
154	>	}
155	>	if(scantype==GMSB) {
156	>	filetoload << "/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
157	>	samplename="GMSB";
158	>	}
159	>
160	>	if(scantype==SMS && Contains(PlottingSetup::ScanSampleDirectory,"T1lh")) {
161	>	filetoload.str("");
162	>	filetoload << "/shome/lbaeni/jzb/" << PlottingSetup::ScanSampleDirectory  << "/SMS_MassGlu_" << mglu << "__MassLSP_" << mlsp << ".root";
163	>	samplename="SMS";
164	>	}
165	>
166		if(scansample.collection.size()<1||!Contains(((scansample.collection)[(scansample.collection).size()-1]).filename,"_"+any2string(a)+"_"+any2string(b))) {
167		dout << "The last sample is NOT the same one as the current one, possibly popping off last one and adding the new one." << endl;
168	<	if((scansample.collection).size()>1) {
168	>	if((scansample.collection).size()>=1) {
169		scansample.RemoveLastSample();
170		}
171		scanfileindex=(scansample.collection).size();
172		//New: Loading file when necessary, not before (avoiding high memory usage and startup times)
173	<	if(scanfileindex!=0) scansample.AddSample(filetoload.str(),"scansample",1,1,false,true,scanfileindex,kRed);
173	>
174	>	scansample.AddSample(filetoload.str(),samplename,1,1,false,true,scanfileindex,kRed);
175	>	dout << "Just added the following file: " << filetoload.str() << endl;
176		} else {
177		dout << "Last sample is the same as the current one. Recycling it." << endl;
178		scanfileindex=(scansample.collection).size()-1;
179		}
180		dout << " Going to use the following file: " << filetoload.str() << endl;
181	<	} // end of t3 case
182	<	else {
181	>	}
182	>
183	>	void load_remote_scan_sample(int a, int b, int &scanfileindex,int scantype,int mglu, int mlsp, bool isretry=false) {
184		write_info(__FUNCTION__,"Hello, CRAB! Might need to load files ...");
185	+	string samplename;
186		PlottingSetup::limitpatience=PlottingSetup::limitpatienceCRAB;
187		stringstream copyfile;
188	+	gSystem->Exec("mkdir -p local_storage");
189		copyfile << "lcg-cp -b -T srmv2 srm://storage01.lcg.cscs.ch:8443/srm/managerv2?SFN=/pnfs/lcg.cscs.ch/cms/trivcat/store/user/buchmann/bussola/" << PlottingSetup::ScanSampleDirectory;
190	<	if(scantype==mSUGRA) copyfile << "/mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root " << "mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root ";
191	<	if(scantype==SMS) copyfile << "/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root " << "SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root ";
192	<	if(scantype==GMSB) copyfile << "/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root " << "GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root ";
190	>	if(scantype==mSUGRA) {
191	>	copyfile<< "/mSUGRA_M0_" << mglu << "__M12_" << mlsp << ".root local_storage/mSUGRA_M0_" << mglu << "__M12_" << mlsp << ".root ";
192	>	samplename="mSUGRA";
193	>	}
194	>	if(scantype==SMS) {
195	>	copyfile << "/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root " << "local_storage/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root ";
196	>	samplename="SMS";
197	>	}
198	>	if(scantype==GMSB) {
199	>	copyfile << "/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root " << "local_storage/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root ";
200	>	samplename="GMSB";
201	>	}
202		stringstream newfilename;
203	<	if(scantype==mSUGRA) newfilename << "mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
204	<	if(scantype==SMS) newfilename << "SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
205	<	if(scantype==GMSB) newfilename << "GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
203	>	//    if(scantype==mSUGRA) newfilename << "local_storage/mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
204	>	if(scantype==mSUGRA) newfilename << "local_storage/mSUGRA_M0_" << mglu << "__M12_" << mlsp << ".root";
205	>	if(scantype==SMS) newfilename << "local_storage/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
206	>	if(scantype==GMSB) newfilename << "local_storage/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
207	>
208	>	SUSYScanSpace::SUSYscantype=scantype;
209
210	<	if(!Contains(((scansample.collection)[(scansample.collection).size()-1]).filename,"_"+any2string(a)+"_"+any2string(b))||(a==0&&b==0)&&!initialized_t2) {
210	>	if((scansample.collection).size()==0||!Contains(((scansample.collection)[(scansample.collection).size()-1]).filename,"_"+any2string(a)+"_"+any2string(b))||(a==0&&b==0)&&!initialized_t2) {
211		dout << "The last sample is NOT the same one as the current one, possibly popping off last one and downloading as well as adding the new one." << endl;
212	<	if((scansample.collection).size()>1) {
212	>	while((scansample.collection).size()>0) {
213		scansample.RemoveLastSample();
214		}
215	+	dout << "New scanfileindex stems from scansample size: " << (scansample.collection).size() << endl;
216		scanfileindex=(scansample.collection).size();
217		//New: Loading file when necessary, not before (avoiding high memory usage and startup times)
218	<	if(scanfileindex!=0||isretry) {
219	<	delete_any_cached_scans();
220	<	dout << "Going to download the scan file with the following copy command: " << copyfile.str() << endl;
221	<	gSystem->Exec(copyfile.str().c_str());
218	>	if(1) {
219	>	if(!doesROOTFileExist(newfilename.str())) {
220	>	dout << "Going to download the scan file with the following copy command: " << copyfile.str() << endl;
221	>	int retcode = gSystem->Exec(copyfile.str().c_str()); // download it if it hasn't been downloaded before
222	>	if(retcode==0) SUSYScanSpace::loaded_files.push_back(newfilename.str()); // adding it to the list of loaded files (will be deleted once we finish running
223	>	}
224		if(doesROOTFileExist(newfilename.str())) {
225		initialized_t2=true;
226	<	scansample.AddSample(newfilename.str(),"scansample",1,1,false,true,scanfileindex,kRed);
226	>	dout << "Going to add the following sample: " << newfilename.str() << endl;
227	>	gSystem->Exec(("ls -ltrh "+newfilename.str()).c_str());
228	>	scansample.AddSample(newfilename.str(),samplename,1,1,false,true,scanfileindex,kRed);
229		} else {
230		srmcpretries++;
231		if(srmcpretries<5) {
232		dout << "The file could not be loaded correctly - retrying!" << endl;
233		sleep(5);
234	<	load_scan_sample(a,b,scanfileindex,scantype,true);
234	>	load_remote_scan_sample(a,b,scanfileindex,scantype,mglu,mlsp,true);
235		} else {
236		dout << "Have tried 5 times to load this sample. Giving up now and failing the program execution" << endl;
237	<	assert(0);
237	>	assert(srmcpretries<5);
238		dout << "The command " << copyfile.str() << " failed to copy the file to the local storage (saving it as " << newfilename.str() << ")" << endl;
239		}
240		}
#	Line 166 | Line 243 | void load_scan_sample(int a, int b, int
243		dout << "Last sample is the same as the current one. Recycling it." << endl;
244		scanfileindex=(scansample.collection).size()-1;
245		}
246	<
247	<	}
248	<	*/
246	>	}
247	>
248	>	void load_scan_sample(int a, int b, int &scanfileindex,int scantype,int mglu, int mlsp, bool isretry=false) {
249	>
250	>	// There is no need to define your sample here. That is done in Setup.C where you define the loading directory!
251
252	<	// WATCH OUT THIS LINE NEEDS TO BE REMOVED
253	<	scanfileindex=0;
252	>	SUSYScanSpace::SUSYscantype=scantype;
253	>	dout << "Going to load file with Xzone=" << a << " and  Yzone=" << b << " for scantype (mSugra: " << (bool)(scantype==mSUGRA) << " , SMS: " << (bool)(scantype==SMS) << " , GMSB: " << (bool)(scantype==GMSB) << ")" << endl;
254	>	char hostname[1023];
255	>	gethostname(hostname,1023);
256	>	if((Contains(hostname,"t3ui")||Contains(hostname,"t3wn"))) load_local_scan_sample(a,b,scanfileindex,scantype,mglu,mlsp,isretry);
257	>	else load_remote_scan_sample(a,b,scanfileindex,scantype,mglu,mlsp,isretry);
258		}
259
260		/*float get_xs(float &altxs, float mglu, float mlsp, string massgluname, string massLSPname, map <  pair<float, float>, map<string, float>  >  &xsec, string mcjzb, bool requireZ)  {
#	Line 205 | Line 288 | void establish_SUSY_limits(string mcjzb,
288		while ((key = (TKey*)nextkey()))
289		{
290		TObject *obj = key->ReadObj();
291	<	if(Contains((string)(obj->GetName()),"mSUGRA")) scantype=mSUGRA;
291	>	if(Contains((string)(obj->GetName()),"SUGRA")) scantype=mSUGRA;
292		if(Contains((string)(obj->GetName()),"GMSB")) scantype=GMSB;
293		}
294
#	Line 215 | Line 298 | void establish_SUSY_limits(string mcjzb,
298		massgluname="M0"; // this is the "x axis" in the limit plot (like the gluino in the SMS case)
299		massLSPname="M12"; // this is the "y axis" in the limit plot (like the LSP in the SMS case)
300		mglustart=m0start;
301	<	xsec=getXsec(PlottingSetup::cbafbasedir+"/Plotting/Modules/external/scale_xsection_nlo1.0_m0_m12_10_0_1v1.txt");
301	>	xsec=getXsec(PlottingSetup::cbafbasedir+"/"+mSUGRAxsFile);
302	>
303		mgluend=m0end;
304		mglustep=m0step;
305		mLSPstart=m12start;
#	Line 223 | Line 307 | void establish_SUSY_limits(string mcjzb,
307		mLSPstep=m12step;
308		prefix="mSUGRA_";
309		dout << "mSUGRA scan has been set up." << endl;
226	–	//    xsec=getXsec("/scratch/buchmann/C/scale_xsection_nlo1.0_m0_m12_10_0_1v1.txt");
310		}
311		if(scantype==SMS) {
312		dout << "SMS scan has been set up." << endl;
#	Line 260 | Line 343 | void establish_SUSY_limits(string mcjzb,
343		TH2F *exp2mlimitmap  = new TH2F((prefix+"exp2mlimitmap"+any2string(jzbSel)).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,(int)((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);//expected limit - 2 sigma
344
345		TH2F *exclmap  = new TH2F((prefix+"exclusionmap"+any2string(jzb_cut[ibin])).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,(int)((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
263	–	TH2F *xsmap  = new TH2F((prefix+"crosssectionmap"+any2string(jzb_cut[ibin])).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,(int)((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
264	–	TH2F *totxsmap  = new TH2F((prefix+"absolutecrosssectionmap"+any2string(jzb_cut[ibin])).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,(int)((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
346		TH2F *flipmap  = new TH2F((prefix+"limitflipmap"+any2string(jzb_cut[ibin])).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,(int)((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
347	+	TH2F *asymptoticmap  = new TH2F((prefix+"asymptoticmap"+any2string(jzb_cut[ibin])).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,(int)((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
348
349
350		if(!fullerr || !mceff || !NEvents) {
#	Line 294 | Line 376 | void establish_SUSY_limits(string mcjzb,
376		continue;
377		}
378		vector<float> sigmas;
379	<	if(scantype!=mSUGRA) sigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped);
379	>	if(scantype!=mSUGRA) {
380	>	sigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped,true);
381	>	write_warning(__FUNCTION__,"Temporarily doing asymptotic limits!!!!");
382	>	}
383		else {
384		//this is a bit trickier - what we want to do is compute the limit fast, and if it is in [0.5 xs, 2xs], compute it again but "correctly"!
385		vector<float> asigmas;
386		asigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped,true); // asymptotic limit first
387		float strength=asigmas[0]/XSmap->GetBinContent(GlobalBin);
388	<	if(strength>0.5&&strength<2) sigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped); // asymptotic limit first
389	<	else sigmas=asigmas;
388	>	/*      if(strength>0.5&&strength<2) {
389	>	sigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped); // asymptotic limit first
390	>	asymptoticmap->SetBinContent(GlobalBin,0);
391	>	} else {*/
392	>	sigmas=asigmas;
393	>	asymptoticmap->SetBinContent(GlobalBin,1);
394	>	/*      }*/
395		exclmap->SetBinContent(GlobalBin,strength);
306	–	xsmap->SetBinContent(GlobalBin,XSmap->GetBinContent(GlobalBin));
396		}
397
309	–
310	–
311	–
398		if(sigmas[0]>0) limitmap->SetBinContent(GlobalBin,sigmas[0]); //anything else is an error code
399		if(sigmas.size()>1) {
400		explimitmap->SetBinContent(GlobalBin,sigmas[1]);
#	Line 323 | Line 409 | void establish_SUSY_limits(string mcjzb,
409		}
410
411		prepare_scan_axis(limitmap,scantype);
412	+	gSystem->Exec("mkdir -p output");
413		TFile *outputfile=new TFile(("output/DistributedLimitsFromSystematics_job"+string(any2string(jobnumber))+"_of_"+string(any2string(njobs))+".root").c_str(),"UPDATE");//needs to be "UPDATE" as we can get to this point for different JZB cuts and don't want to erase the previous data :-)
414		outputfile->cd();
415		limitmap->Write();
416	+	XSmap->Write();
417		flipmap->Write();
418	+	asymptoticmap->Write();
419		if(doexpected) {
420		explimitmap->Write();
421		exp1plimitmap->Write();
#	Line 334 | Line 423 | void establish_SUSY_limits(string mcjzb,
423		exp2plimitmap->Write();
424		exp2mlimitmap->Write();
425		}
426	<	if(scantype==mSUGRA) {
338	<	exclmap->Write();
339	<	xsmap->Write();
340	<	totxsmap->Write();
341	<	}
426	>	if(scantype==mSUGRA) exclmap->Write();
427		outputfile->Close();
428		delete limcanvas;
429		}
#	Line 350 | Line 435 | void establish_SUSY_limits(string mcjzb,
435		write_error(__FUNCTION__,"You provided an invalid systematics file. Please change it ... ");
436		return;
437		}
438	<	for(int ibin=0;ibin<jzb_cut.size();ibin++) {
438	>	for(int ibin=0;ibin<(int)jzb_cut.size();ibin++) {
439		establish_SUSY_limits(mcjzb,datajzb,jzb_cut,requireZ, peakerror, fsyst, ibin,njobs, jobnumber);
440		}
441		}
442
443
444	+	string SimpleFormat(string input) {
445	+	string output;
446	+	for (unsigned int i=0; i<input.length();i++) {
447	+	if (isalnum(input[i])) //If current character is alpha-numeric
448	+	output+= input[i]; //Add the character to strConverted
449	+
450	+	}
451	+	return output;
452	+	}
453
454
455	<	void scan_SUSY_parameter_space(string mcjzb,string datajzb,vector<float> jzb_cut,bool requireZ, float peakerror, int ibin,float njobs=-1, float jobnumber=-1, bool systematicsonly=false,bool efficiencyonly=false) {
455	>	void scan_SUSY_parameter_space_for_one_cut(string mcjzb,string datajzb,vector<float> jzb_cut,bool requireZ, float peakerror, float peakerrordata, int ibin,float njobs=-1, float jobnumber=-1, bool systematicsonly=false,bool efficiencyonly=false, bool shapeanalysis=false) {
456	>
457	>	if(shapeanalysis&&ibin>0) return; // pointless since we're doing a shape analysis :-)
458	>
459		bool runninglocally=true;
460		if(njobs>-1&&jobnumber>-1) {
461		runninglocally=false;
#	Line 370 | Line 467 | void scan_SUSY_parameter_space(string mc
467		}
468
469		jzbSel=jzb_cut[ibin];
470	+	if(shapeanalysis) jzbSel=0;
471		geqleq="geq";
472		automatized=true;
473		mcjzbexpression=mcjzb;
#	Line 380 | Line 478 | void scan_SUSY_parameter_space(string mc
478		string prefix="SMS_";
479		// up to here, everything is set up for SMS; now we need to switch stuff around if we're dealing with an mSUGRA scan!
480		int scantype=SMS;
481	<	if(Contains((scansample.collection)[0].samplename,"mSUGRA")) scantype=mSUGRA;
482	<	if(Contains((scansample.collection)[0].samplename,"GMSB")) scantype=GMSB;
481	>	if(Contains((PlottingSetup::ScanSampleDirectory),"SUGRA")) scantype=mSUGRA;
482	>	if(Contains((PlottingSetup::ScanSampleDirectory),"GMSB")) scantype=GMSB;
483
484		if(scantype==mSUGRA) {
485		massgluname="M0"; // this is the "x axis" in the limit plot (like the gluino in the SMS case)
#	Line 415 | Line 513 | void scan_SUSY_parameter_space(string mc
513
514		gStyle->SetPalette(100, MyPalette);
515
516	+	TH2F *rawlimitmap  = new TH2F((prefix+"rawlimitmap"+any2string(jzbSel)).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);//observed limit
517	+	TH2F *rawelimitmap  = new TH2F((prefix+"rawelimitmap"+any2string(jzbSel)).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);//observed limit
518	+
519		TH2F *limitmap  = new TH2F((prefix+"limitmap"+any2string(jzbSel)).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);//observed limit
520		TH2F *explimitmap  = new TH2F((prefix+"explimitmap"+any2string(jzbSel)).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);//expected limit
521		TH2F *exp1plimitmap  = new TH2F((prefix+"exp1plimitmap"+any2string(jzbSel)).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);//expected limit + 1 sigma
#	Line 441 | Line 542 | mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLS
542		TH2F *systotmap      = new TH2F((prefix+"systotmap"+any2string(jzbSel)).c_str(),"",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
543		TH2F *sysstatmap     = new TH2F((prefix+"sysstatmap"+any2string(jzbSel)).c_str(),"",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
544		TH2F *XSmap     = new TH2F((prefix+"XS"+any2string(jzbSel)).c_str(),"",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
545	+	TH2F *absXSmap     = new TH2F((prefix+"absXS"+any2string(jzbSel)).c_str(),"",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
546	+	TH2F *FilterEff     = new TH2F((prefix+"FilterEfficiency"+any2string(jzbSel)).c_str(),"",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
547	+
548	+	TH2F *ScanIdentifier = new TH2F(("ScanIdentifier"+any2string(jzbSel)+SimpleFormat(PlottingSetup::ScanSampleDirectory)).c_str(),("ScanIdentifier"+any2string(jzbSel)+SimpleFormat(PlottingSetup::ScanSampleDirectory)).c_str(),1,0,1,1,0,1);
549	+
550	+
551	+	TH2F *asymptoticmap  = new TH2F((prefix+"asymptoticmap"+any2string(jzbSel)).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,(int)((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
552	+	// only used for asymptoticmap
553
554		TH2F *imposedxmap;
555		TH2F *realxmap;
#	Line 452 | Line 561 | mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLS
561
562		TH2F *flipmap        = new TH2F((prefix+"flipmap"+any2string(jzbSel)).c_str(),"",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
563
455	–	TH2F *nZmap    = new TH2F((prefix+"nZmap"+any2string(jzb_cut[ibin])).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,(int)((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
456	–
564		if(ibin==0) {
565		centermap    = new TH2F((prefix+"centermap").c_str(),"",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
566		widthmap       = new TH2F((prefix+"widthmap").c_str(),"",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
#	Line 462 | Line 569 | mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLS
569		realxmap       = new TH2F((prefix+"realxmap").c_str(),"",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
570		}
571
572	+
573	+	if(shapeanalysis) jzbSel=jzb_cut[ibin];
574		vector<int> susy_Zdecay_origin;
575		vector<TH2F*> h_susy_Zdecay_origin;
576		// we only generate these histograms IF we're dealing with the lowest cut! (we don't need them n times, they're the same for all of them)
#	Line 483 | Line 592 | mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLS
592
593		write_warning(__FUNCTION__,"CURRENTLY SWITCHING AUTOMATIZED MODE OFF!");automatized=false;
594
486	–	float rightmargin=gStyle->GetPadRightMargin();
595		gStyle->SetPadRightMargin(0.15);
596
597		TCanvas *limcanvas = new TCanvas("limcanvas","Limit canvas");
#	Line 498 | Line 606 | write_warning(__FUNCTION__,"CURRENTLY SW
606
607		map <  pair<float, float>, map<string, float>  >  xsec;
608		if(scantype==mSUGRA) {
609	<	xsec=getXsec(PlottingSetup::mSUGRAxsFile);
609	>	xsec=getXsec(PlottingSetup::cbafbasedir+"/"+mSUGRAxsFile);
610		}
611
612		int ipoint=-1;
#	Line 508 | Line 616 | write_warning(__FUNCTION__,"CURRENTLY SW
616		ipoint++;
617		int GlobalBin=efficiencymap->FindBin(mglu,mlsp);
618		if(!runninglocally&&!do_this_point(ipoint,(int)Npoints,(int)jobnumber,(int)njobs)) continue;
619	+
620		int flipped=0;
621		float result=-987,resulterr=-987;
622		float m0trees = PlottingSetup::mgluend;
#	Line 515 | Line 624 | write_warning(__FUNCTION__,"CURRENTLY SW
624		int a = int((PlottingSetup::ScanXzones*mlsp)/(m12trees+1));
625		int b = int((PlottingSetup::ScanYzones*mglu)/(m0trees+1));
626		int scanfileindex=PlottingSetup::ScanYzones*a+b;
627	<	load_scan_sample(a,b,scanfileindex,scantype);
627	>	load_scan_sample(a,b,scanfileindex,scantype,mglu,mlsp);
628
629		clock_t start,finish;
630		start = clock(); // starting the clock to measure how long the computation takes!
#	Line 529 | Line 638 | write_warning(__FUNCTION__,"CURRENTLY SW
638		continue;
639		} else {
640		dout << "OK! This point ("<<ipoint<<") with configuration ("<<massgluname<<"="<<mglu<<" , "<<massLSPname<<"="<<mlsp << ") contains " << nevents << " and will therefore \033[1;32m not be skipped \033[0m. " << endl;
641	<	if(PlottingSetup::RestrictToMassPeak==true) cout << "Analysis type: on-peak" << endl;
642	<	if(PlottingSetup::RestrictToMassPeak==false) cout << "Analysis type: offpeak" << endl;
641	>	if(PlottingSetup::RestrictToMassPeak==true) dout << "Analysis type: on-peak" << endl;
642	>	if(PlottingSetup::RestrictToMassPeak==false) dout << "Analysis type: offpeak" << endl;
643	>	if(shapeanalysis) dout << "Shape analysis? Yes." << endl;
644	>	else dout << "Shape analysis? No." << endl;
645		}
646		float hitrecord=0;
647		int hitrecordholder=0;
648		float hitsecondplace=0;
649		int hitsecondplaceholder=0;
650	<	for(int imo=0;imo<susy_Zdecay_origin.size()&&ibin==0;imo++) {
650	>	for(int imo=0;imo<(int)susy_Zdecay_origin.size()&&ibin==0;imo++) {
651		if(ibin>0) continue;
652		int MotherParticlePDGid = susy_Zdecay_origin[imo];
653		(scansample.collection)[scanfileindex].events->Draw("eventNum",("("+addcut.str()+")&&(abs(SourceOfZ-"+any2string(MotherParticlePDGid)+")<0.5)").c_str(),"goff");
#	Line 578 | Line 689 | write_warning(__FUNCTION__,"CURRENTLY SW
689		}
690
691
692	<	if(nevents!=0&&(efficiencyonly||systematicsonly)) {
692	>	if(scantype==mSUGRA) {
693	>	float absxs=0;
694	>	for(int i=0;i<12;i++) absxs+=GetXSecForPointAndChannel(mglu,mlsp,xsec,i);
695	>	TFile *FilterEffFile = new TFile(PlottingSetup::FilterEfficiencyFile.c_str());
696	>	TH2F *FilterEfficiency = (TH2F*)FilterEffFile->Get("FilterEfficiency");
697	>	float filtereff = FilterEfficiency->GetBinContent(FilterEfficiency->FindBin(mglu,mlsp));
698	>	FilterEff->SetBinContent(GlobalBin,filtereff);
699	>	FilterEffFile->Close();
700	>	absXSmap->SetBinContent(GlobalBin,absxs);
701	>	absxs*=filtereff;
702	>	XSmap->SetBinContent(GlobalBin,absxs);
703	>	}
704	>
705	>	if(nevents!=0) {
706		Value effwosigcont = MCefficiency((scansample.collection)[scanfileindex].events,result,resulterr,flipped,mcjzb,requireZ,(int)nevents,scantype,xsec,addcut.str(),-1);
707		if(result<0&&allowflipping) {
708		write_info(__FUNCTION__,"Found negative efficiency (signal distribution populating JZB<0 more than JZB>0) - flipping analysis!");
#	Line 615 | Line 739 | write_warning(__FUNCTION__,"CURRENTLY SW
739		Neventsmap->SetBinContent(GlobalBin,nevents);
740		ipointmap->SetBinContent(GlobalBin,ipoint);
741		if(efficiencyonly) continue;
742	<
743	<	do_systematics_for_one_file((scansample.collection)[scanfileindex].events,(int)nevents,"SUSY SCAN", systematics,flipped, xsec, mcjzb,datajzb,peakerror,requireZ, addcut.str(),true,scantype);
744	<	float JZBcutat = systematics[0][0];
745	<	float mceff    = systematics[0][1];
746	<	float mcefferr = systematics[0][2];//MC stat error
747	<	float toterr   = systematics[0][4];
748	<	float sys_jes  = systematics[0][5]; // Jet Energy Scale
749	<	float sys_jsu  = systematics[0][6]; // JZB scale uncertainty
750	<	float sys_res  = systematics[0][7]; // resolution
751	<	float mcwoscef = systematics[0][8]; // efficiency without signal contamination
752	<	float mcwoscefr= systematics[0][9]; // error on efficiency without signal contamination
753	<	float sys_pdf   = 0;
754	<	if(systematics[0].size()>10) sys_pdf = systematics[0][10]; // PDF
755	<
756	<	if(mceff!=mceff||toterr!=toterr||mceff<0 && (!systematicsonly&&!efficiencyonly)) {
757	<	dout << "Limits can't be calculated for this configuration (mglu="<<mglu<<" , mlsp="<<mlsp << ") as either the efficiency or its error are not positive numbers! (mceff="<<mceff<<" and toterr="<<toterr<<")"<< endl;
758	<	continue;
759	<	} else {
760	<	if(!systematicsonly&&!efficiencyonly) {
761	<	dout << "Calculating limit now for "<<massgluname<<"="<<mglu<<" , "<<massLSPname<<"="<<mlsp <<endl;
762	<	vector<float> sigmas;
763	<	string plotfilename=(string)(TString((scansample.collection)[scanfileindex].samplename)+TString(massgluname)+TString(any2string(mglu))+TString("__")+TString(massLSPname)+TString(any2string(mlsp))+TString(".png"));
764	<	sigmas=compute_one_upper_limit(mceff,toterr,ibin,mcjzb,plotfilename,true,flipped);
765	<	//      do_limit_wrapper(mceff,toterr,ibin,mcjzb,sigmas,plotfilename);
766	<	dout << "back in " << __FUNCTION__ << endl;
767	<	if(sigmas[0]>-0.5) { // negative sigmas are the error signature of do_limit_wrapper, so we want to exclude them.
768	<	limitmap->SetBinContent(GlobalBin,sigmas[0]);
769	<	if(sigmas.size()>1) {
770	<	explimitmap->SetBinContent(GlobalBin,sigmas[1]);
771	<	exp1plimitmap->SetBinContent(GlobalBin,sigmas[2]);
772	<	exp1mlimitmap->SetBinContent(GlobalBin,sigmas[3]);
773	<	exp2plimitmap->SetBinContent(GlobalBin,sigmas[4]);
774	<	exp2mlimitmap->SetBinContent(GlobalBin,sigmas[5]);
742	>
743	>	if(!shapeanalysis) {
744	>	do_systematics_for_one_file((scansample.collection)[scanfileindex].events,(int)nevents,"SUSY SCAN", systematics,flipped, xsec, mcjzb,datajzb,peakerror,requireZ, addcut.str(),true,scantype);
745	>	//      float JZBcutat = systematics[0][0];
746	>	float mceff    = systematics[0][1];
747	>	float mcefferr = systematics[0][2];//MC stat error
748	>	float toterr   = systematics[0][4];
749	>	float sys_jes  = systematics[0][5]; // Jet Energy Scale
750	>	float sys_jsu  = systematics[0][6]; // JZB scale uncertainty
751	>	float sys_res  = systematics[0][7]; // resolution
752	>	//      float mcwoscef = systematics[0][8]; // efficiency without signal contamination
753	>	//      float mcwoscefr= systematics[0][9]; // error on efficiency without signal contamination
754	>	float sys_pdf   = 0;
755	>	if(systematics[0].size()>10) sys_pdf = systematics[0][10]; // PDF
756	>	dout << "Going to store: total error is: " << toterr << " and the relative tot error, i.e. total error / mc efficiency is " << toterr/mceff << endl;
757	>	if(mceff!=mceff||toterr!=toterr||mceff<0 && (!systematicsonly&&!efficiencyonly)) {
758	>	dout << "Limits can't be calculated for this configuration (" << massgluname <<"="<<mglu<<" , " << massLSPname << "="<<mlsp << ") as either the efficiency or its error are not positive numbers! (mceff="<<mceff<<"and toterr="<<toterr<<")"<< endl;
759	>	continue;
760	>	} else {
761	>	// Systematics and efficiencies make sense (i.e. non-zero, and all numbers)
762	>	if(!systematicsonly&&!efficiencyonly) {
763	>	dout << "Calculating limit now for "<<massgluname<<"="<<mglu<<" , "<<massLSPname<<"="<<mlsp <<endl;
764	>	vector<float> sigmas;
765	>	string plotfilename=(string)(TString((scansample.collection)[scanfileindex].samplename)+TString(massgluname)+TString(any2string(mglu))+TString("__")+TString(massLSPname)+TString(any2string(mlsp))+TString(".png"));
766	>	sigmas=compute_one_upper_limit(mceff,toterr,ibin,mcjzb,plotfilename,true,flipped);
767	>	//          do_limit_wrapper(mceff,toterr,ibin,mcjzb,sigmas,plotfilename);
768	>	if(sigmas[0]>-0.5) { // negative sigmas are the error signature of do_limit_wrapper, so we want to exclude them.
769	>	limitmap->SetBinContent(GlobalBin,sigmas[0]);
770	>	if(sigmas.size()>1) {
771	>	explimitmap->SetBinContent(GlobalBin,sigmas[1]);
772	>	exp1plimitmap->SetBinContent(GlobalBin,sigmas[2]);
773	>	exp1mlimitmap->SetBinContent(GlobalBin,sigmas[3]);
774	>	exp2plimitmap->SetBinContent(GlobalBin,sigmas[4]);
775	>	exp2mlimitmap->SetBinContent(GlobalBin,sigmas[5]);
776	>	}
777	>	sysjesmap->SetBinContent(GlobalBin,sys_jes);
778	>	sysjsumap->SetBinContent(GlobalBin,sys_jsu);
779	>	sysresmap->SetBinContent(GlobalBin,sys_res);
780	>	syspdfmap->SetBinContent(GlobalBin,sys_pdf);
781	>	systotmap->SetBinContent(GlobalBin,toterr/mceff);//total relative (!) error
782	>	sysstatmap->SetBinContent(GlobalBin,mcefferr);//total relative (!) error
783	>	dout << "A limit has been added at " << sigmas[0] << " for m_{glu}="<<mglu << " and m_{lsp}="<<mlsp<<endl;
784	>	}//end of if sigma is positive
785	>	}//end of not systematics only condition
786	>	if(systematicsonly) {
787	>	sysjesmap->SetBinContent(GlobalBin,sys_jes);
788	>	sysjsumap->SetBinContent(GlobalBin,sys_jsu);
789	>	sysresmap->SetBinContent(GlobalBin,sys_res);
790	>	syspdfmap->SetBinContent(GlobalBin,sys_pdf);
791	>	systotmap->SetBinContent(GlobalBin,toterr/mceff);//total relative (!) error
792	>	sysstatmap->SetBinContent(GlobalBin,mcefferr);//total relative (!) error
793	>	}
794	>	}//efficiency is valid
795	>	}
796	>	if(shapeanalysis) {
797	>	//shape analysis
798	>	dout << "This is the shape analysis. Have fun with it!" << endl;
799	>	float quickmceff,quickmcefferr;
800	>	MCefficiency((scansample.collection)[scanfileindex].events,quickmceff,quickmcefferr,flipped,mcjzb,requireZ,(int)nevents,scantype,xsec,addcut.str(),-1);
801	>	if(quickmceff==0) {
802	>	write_warning(__FUNCTION__,"The MC efficiency is zero - need to skip this point.");
803	>	continue;
804	>	}
805	>	stringstream PointName;
806	>	PointName << massgluname << "_" << mglu << "__" << massLSPname << "_" << mlsp;
807	>	SUSYScanSpace::SavedMLSP=mlsp;
808	>	SUSYScanSpace::SavedMGlu=mglu;
809	>	SUSYScanSpace::SavedMLSPname=massLSPname;
810	>	SUSYScanSpace::SavedMGluname=massgluname;
811	>	float referenceXS=0;
812	>	if(scantype==mSUGRA) referenceXS=XSmap->GetBinContent(GlobalBin);
813	>	else referenceXS=getSMSxs(mlsp,mglu);
814	>	float low_fullCLs=referenceXS*0.5;
815	>	float high_fullCLs=referenceXS*2.5;
816	>	if(scantype!=mSUGRA) {
817	>	low_fullCLs=-1;
818	>	high_fullCLs=10e10;
819	>	}
820	>	ShapeDroplet droplet = LimitsFromShapes(low_fullCLs,high_fullCLs, (scansample.collection)[scanfileindex].events,addcut.str().c_str(),PointName.str(),mcjzb,datajzb,jzb_cut,peakerror,peakerror);
821	>	if(!isValidShapeResult(droplet)) {
822	>	write_error(__FUNCTION__,"Something went horribly wrong with the shape analysis - ignoring this point!!!");
823	>	//        scan_SUSY_parameter_space(mcjzb,datajzb,jzb_cut,requireZ, peakerror, peakerrordata, (float)Npoints, (float)ipoint, systematicsonly, efficiencyonly,false); // don't uncomment this line; it would lead to multiple limitmap histos ...
824	>	continue;
825		}
826
827	<	sysjesmap->SetBinContent(GlobalBin,sys_jes);
828	<	sysjsumap->SetBinContent(GlobalBin,sys_jsu);
829	<	sysresmap->SetBinContent(GlobalBin,sys_res);
830	<	syspdfmap->SetBinContent(GlobalBin,sys_pdf);
831	<	systotmap->SetBinContent(GlobalBin,toterr/mceff);//total relative (!) error
832	<	sysstatmap->SetBinContent(GlobalBin,mcefferr);//total relative (!) error
833	<	dout << "A limit has been added at " << sigmas[0] << " for m_{glu}="<<mglu << " and m_{lsp}="<<mlsp<<endl;
834	<	} //end of if sigma is positive
835	<
836	<	//end of not systematics only condition
837	<	}
838	<	if(systematicsonly) {
839	<	sysjesmap->SetBinContent(GlobalBin,sys_jes);
840	<	sysjsumap->SetBinContent(GlobalBin,sys_jsu);
841	<	sysresmap->SetBinContent(GlobalBin,sys_res);
842	<	syspdfmap->SetBinContent(GlobalBin,sys_pdf);
843	<	systotmap->SetBinContent(GlobalBin,toterr/mceff);//total relative (!) error
844	<	sysstatmap->SetBinContent(GlobalBin,mcefferr);//total relative (!) error
845	<	}
672	<	}//efficiency is valid
827	>	float excludedXSobs = droplet.observed/droplet.SignalIntegral;
828	>	absXSmap->SetBinContent(GlobalBin,referenceXS);
829	>	write_info(__FUNCTION__,"Established limit at "+any2string(excludedXSobs)+" , ref XS is "+any2string(referenceXS));
830	>
831	>	rawlimitmap->SetBinContent(GlobalBin,droplet.observed);
832	>	rawelimitmap->SetBinContent(GlobalBin,droplet.expected);
833	>	limitmap->SetBinContent(GlobalBin,droplet.observed/droplet.SignalIntegral);
834	>	explimitmap->SetBinContent(GlobalBin,droplet.expected/droplet.SignalIntegral);
835	>	exp1plimitmap->SetBinContent(GlobalBin,droplet.expectedPlus1Sigma/droplet.SignalIntegral);
836	>	exp1mlimitmap->SetBinContent(GlobalBin,droplet.expectedMinus1Sigma/droplet.SignalIntegral);
837	>	exp2plimitmap->SetBinContent(GlobalBin,droplet.expectedPlus2Sigma/droplet.SignalIntegral);
838	>	exp2mlimitmap->SetBinContent(GlobalBin,droplet.expectedMinus2Sigma/droplet.SignalIntegral);
839	>	dout <<"Saved the following shape limit: " << droplet.observed/droplet.SignalIntegral << " (observed) , " << droplet.expected/droplet.SignalIntegral << " (expected) " << endl;
840	>	sysjesmap->SetBinContent(GlobalBin,droplet.JES);
841	>	sysjsumap->SetBinContent(GlobalBin,droplet.JSU);
842	>	syspdfmap->SetBinContent(GlobalBin,droplet.PDF);
843	>	systotmap->SetBinContent(GlobalBin,droplet.toterr);
844	>	sysstatmap->SetBinContent(GlobalBin,droplet.staterr);
845	>	}
846		finish = clock();
847		timemap->SetBinContent(GlobalBin,((float(finish)-float(start))/CLOCKS_PER_SEC));
848		}
849		}
677	–
850		prepare_scan_axis(limitmap,scantype);
851		prepare_scan_axis(sysjesmap,scantype);
852		prepare_scan_axis(sysjsumap,scantype);
#	Line 685 | Line 857 | write_warning(__FUNCTION__,"CURRENTLY SW
857		prepare_scan_axis(timemap,scantype);
858
859		if(!systematicsonly&&!efficiencyonly) {
860	<	limcanvas->cd();
860	>	TObject* old=gDirectory->GetList()->FindObject("limcanvas");
861	>	if(old) gDirectory->GetList()->Remove(old);
862	>	TCanvas *limcanvas = new TCanvas("limcanvas","Limit canvas");
863		limitmap->Draw("COLZ");
864		string titletobewritten="Limits in LSP-Glu plane";
865		if(scantype==mSUGRA) titletobewritten="Limits in m_{1/2}-m_{0} plane";
#	Line 699 | Line 873 | write_warning(__FUNCTION__,"CURRENTLY SW
873		if(systematicsonly) outputfile=new TFile(("output/DistributedSystematics_job"+string(any2string(jobnumber))+"_of_"+string(any2string(njobs))+".root").c_str(),"UPDATE");//needs to be "UPDATE" as we can get to this point for different JZB cuts and don't want to erase the previous data :-)
874		else outputfile=new TFile(("output/DistributedLimits_job"+string(any2string(jobnumber))+"_of_"+string(any2string(njobs))+".root").c_str(),"UPDATE");//needs to be "UPDATE" as we can get to this point for
875		limitmap->Write();
876	+	rawlimitmap->Write();
877	+	rawelimitmap->Write();
878		if(doexpected) {
879		explimitmap->Write();
880		exp1plimitmap->Write();
#	Line 724 | Line 900 | write_warning(__FUNCTION__,"CURRENTLY SW
900		systotmap->Write();
901		sysstatmap->Write();
902		XSmap->Write();
903	+	absXSmap->Write();
904	+	FilterEff->Write();
905	+	if(shapeanalysis) asymptoticmap->Write();
906		if(ibin==0) imposedxmap->Write();
907		if(ibin==0) realxmap->Write();
908		Neventsmap->Write();
909	+	ScanIdentifier->Write();
910		ipointmap->Write();
911		timemap->Write();
912	<	for(int i=0;i<susy_Zdecay_origin.size();i++) {
912	>	for(int i=0;i<(int)susy_Zdecay_origin.size();i++) {
913		h_susy_Zdecay_origin[i]->Write();
914		}
915		outputfile->Close();
#	Line 782 | Line 962 | write_warning(__FUNCTION__,"CURRENTLY SW
962		noscefficiencymap->Write();
963		efficiencymapContam->Write();
964		Neventsmap->Write();
965	+	ScanIdentifier->Write();
966		ipointmap->Write();
967		syspdfmap->Write();
968		systotmap->Write();
969		sysstatmap->Write();
970	+	absXSmap->Write();
971		XSmap->Write();
972	+	FilterEff->Write();
973	+	if(shapeanalysis) asymptoticmap->Write();
974		if(ibin==0) imposedxmap->Write();
975		if(ibin==0) realxmap->Write();
976		timemap->Write();
977	<	for(int i=0;i<susy_Zdecay_origin.size();i++) {
977	>	for(int i=0;i<(int)susy_Zdecay_origin.size();i++) {
978		h_susy_Zdecay_origin[i]->Write();
979		}
980		outputfile->Close();
#	Line 857 | Line 1041 | write_warning(__FUNCTION__,"CURRENTLY SW
1041		TFile *outputfile=new TFile(("output/DistributedSystematics_job"+string(any2string(jobnumber))+"_of_"+string(any2string(njobs))+".root").c_str(),"UPDATE");
1042		ipointmap->Write();
1043		Neventsmap->Write();
1044	+	ScanIdentifier->Write();
1045		noscefficiencymap->Write();
1046		efficiencymapContam->Write();
1047		efficiencymap->Write();
#	Line 872 | Line 1057 | write_warning(__FUNCTION__,"CURRENTLY SW
1057		outputfile->Close();
1058		}
1059		}//end of efficiencies only
1060	<
1061	<	for(int i=0;i<susy_Zdecay_origin.size();i++) {
1060	>
1061	>	for(int i=0;i<(int)susy_Zdecay_origin.size();i++) {
1062		delete h_susy_Zdecay_origin[i];
1063		}
1064
#	Line 895 | Line 1080 | write_warning(__FUNCTION__,"CURRENTLY SW
1080		delete efficiencymapJets;
1081		delete efficiencymapSignal;
1082		delete Neventsmap;
1083	+	delete ScanIdentifier;
1084		delete ipointmap;
1085		delete syspdfmap;
1086		delete systotmap;
1087		delete sysstatmap;
1088		delete XSmap;
1089	<	delete limcanvas;
1089	>	delete absXSmap;
1090	>	TObject* old=gDirectory->GetList()->FindObject("limcanvas");
1091	>	if(old) gDirectory->GetList()->Remove(old);
1092		}
1093
1094	<	void scan_SUSY_parameter_space(string mcjzb,string datajzb,vector<float> jzb_cut,bool requireZ, float peakerror, float njobs=-1, float jobnumber=-1, bool systonly=false, bool effonly=false) {
1094	>	void scan_SUSY_parameter_space(string mcjzb,string datajzb,vector<float> jzb_cut,bool requireZ, float peakerror, float peakerrordata, float njobs=-1, float jobnumber=-1, bool systonly=false, bool effonly=false,bool shapeanalysis=false) {
1095		dout << "Starting the SUSY scan now with all " << jzb_cut.size() << " bin(s)" << endl;
1096	<	for(int ibin=0;ibin<jzb_cut.size();ibin++) {
1097	<	scan_SUSY_parameter_space(mcjzb,datajzb,jzb_cut,requireZ, peakerror, ibin, njobs, jobnumber,systonly,effonly);
1096	>	for(int ibin=0;ibin<(int)jzb_cut.size();ibin++) {
1097	>	scan_SUSY_parameter_space_for_one_cut(mcjzb,datajzb,jzb_cut,requireZ, peakerror, peakerrordata, ibin, njobs, jobnumber,systonly,effonly,shapeanalysis);
1098		}
1099	+	delete_any_cached_scans();// tidy up ...
1100		}
1101

Diff Legend

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