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.16 by buchmann, Fri May 4 11:48:39 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.observed>50000) {
35	+	// error code of ShapeLimit/CreateModel.sh, invoked if there were too many errors during computation
36	+	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.");
37	+	return false;
38	+	}
39	+	if(droplet.observed==-12345) {
40	+	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.");
41	+	return false;
42	+	}
43	+	if(droplet.expected>-1&&droplet.expected==0&&droplet.observed>5000) {
44	+	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.");
45	+	return false;
46	+	}
47	+	if(droplet.expected>-1&&droplet.expected/droplet.observed>50) {
48	+	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.");
49	+	return false;
50	+	}
51	+	if(droplet.expected>-1&&droplet.observed/droplet.expected>50) {
52	+	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.");
53	+	return false;
54	+	}
55	+	return true;
56	+	}
57	+
58		void susy_scan_axis_labeling(TH2F *histo) {
59		histo->GetXaxis()->SetTitle("m_{#Chi_{2}^{0}}-m_{LSP}");
60		histo->GetXaxis()->CenterTitle();
#	Line 61 | Line 92 | void set_SUSY_style() {
92
93		gStyle->SetPalette(100, MyPalette);
94
64	–	float rightmargin=gStyle->GetPadRightMargin();
95		gStyle->SetPadRightMargin(0.15);
96
97		}
98
99		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");
100		char hostname[1023];
101		gethostname(hostname,1023);
102	<	if((Contains(hostname,"t3ui")||Contains(hostname,"t3wn"))) return false;
103	<	else {
102	>
103	>
104	>	/*
105	>	//This should only be called via CRAB
106	>	dout << "   Deleting all cached files." << endl;
107	>	gSystem->Exec("ls -ltrh | grep root");
108	>	if(!((Contains(hostname,"t3ui")||Contains(hostname,"t3wn"))))
109	>	{
110		vector<string> all_files;
111		char currentpath[1024];
112		TString directory=TString(getcwd(currentpath,1024));
#	Line 86 | Line 119 | bool delete_any_cached_scans() {
119		gSystem->Exec(((string)"rm "+all_files[ifile]).c_str());
120		}
121		}
122	<	return true;
122	>	return true;
123	>	}
124	>	*/
125	>	if(!((Contains(hostname,"t3ui")||Contains(hostname,"t3wn")))) {
126	>	dout << "Going to purge files in local_storage that have been loaded!" << endl;
127	>	for(int i=0;i<(int)SUSYScanSpace::loaded_files.size();i++) {
128	>	gSystem->Exec(((string)"rm "+SUSYScanSpace::loaded_files[i]).c_str());
129	>	dout << "      Purging : Deleted file " << SUSYScanSpace::loaded_files[i] << endl;
130	>	}
131		}
132	+	return true;
133		}
134
135		bool initialized_t2=false;
136
137		int srmcpretries=0;
138
139	<	void load_scan_sample(int a, int b, int &scanfileindex,int scantype,bool isretry=false) {
140	<	/*
141	<	// 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
139	>	void load_local_scan_sample(int a, int b, int &scanfileindex,int scantype,int mglu, int mlsp, bool isretry=false) {
140	>	stringstream filetoload;
141	>	string samplename;
142		filetoload << "/shome/buchmann/ntuples/"<<PlottingSetup::ScanSampleDirectory;
143	<	if(scantype==mSUGRA) filetoload << "/mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
144	<	if(scantype==SMS) filetoload << "/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
145	<	if(scantype==GMSB) filetoload << "/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
143	>	if(scantype==mSUGRA) {
144	>	//filetoload << "/mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
145	>	filetoload.str("");
146	>	filetoload << "/shome/lbaeni/jzb/" << PlottingSetup::ScanSampleDirectory  << "/mSUGRA_M0_" << mglu << "__M12_" << mlsp << ".root";
147	>	samplename="mSUGRA";
148	>	//      filetoload << "/mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
149	>	}
150	>	if(scantype==SMS) {
151	>	filetoload << "/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
152	>	samplename="SMS";
153	>	}
154	>	if(scantype==GMSB) {
155	>	filetoload << "/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
156	>	samplename="GMSB";
157	>	}
158	>
159	>	if(scantype==SMS && Contains(PlottingSetup::ScanSampleDirectory,"T1lh")) {
160	>	filetoload.str("");
161	>	filetoload << "/shome/lbaeni/jzb/" << PlottingSetup::ScanSampleDirectory  << "/SMS_MassGlu_" << mglu << "__MassLSP_" << mlsp << ".root";
162	>	samplename="SMS";
163	>	}
164	>
165		if(scansample.collection.size()<1||!Contains(((scansample.collection)[(scansample.collection).size()-1]).filename,"_"+any2string(a)+"_"+any2string(b))) {
166		dout << "The last sample is NOT the same one as the current one, possibly popping off last one and adding the new one." << endl;
167	<	if((scansample.collection).size()>1) {
167	>	if((scansample.collection).size()>=1) {
168		scansample.RemoveLastSample();
169		}
170		scanfileindex=(scansample.collection).size();
171		//New: Loading file when necessary, not before (avoiding high memory usage and startup times)
172	<	if(scanfileindex!=0) scansample.AddSample(filetoload.str(),"scansample",1,1,false,true,scanfileindex,kRed);
172	>
173	>	scansample.AddSample(filetoload.str(),samplename,1,1,false,true,scanfileindex,kRed);
174	>	dout << "Just added the following file: " << filetoload.str() << endl;
175		} else {
176		dout << "Last sample is the same as the current one. Recycling it." << endl;
177		scanfileindex=(scansample.collection).size()-1;
178		}
179		dout << " Going to use the following file: " << filetoload.str() << endl;
180	<	} // end of t3 case
181	<	else {
180	>	}
181	>
182	>	void load_remote_scan_sample(int a, int b, int &scanfileindex,int scantype,int mglu, int mlsp, bool isretry=false) {
183		write_info(__FUNCTION__,"Hello, CRAB! Might need to load files ...");
184	+	string samplename;
185		PlottingSetup::limitpatience=PlottingSetup::limitpatienceCRAB;
186		stringstream copyfile;
187	+	gSystem->Exec("mkdir -p local_storage");
188		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;
189	<	if(scantype==mSUGRA) copyfile << "/mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root " << "mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root ";
190	<	if(scantype==SMS) copyfile << "/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root " << "SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root ";
191	<	if(scantype==GMSB) copyfile << "/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root " << "GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root ";
189	>	if(scantype==mSUGRA) {
190	>	copyfile<< "/mSUGRA_M0_" << mglu << "__M12_" << mlsp << ".root local_storage/mSUGRA_M0_" << mglu << "__M12_" << mlsp << ".root ";
191	>	samplename="mSUGRA";
192	>	}
193	>	if(scantype==SMS) {
194	>	copyfile << "/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root " << "local_storage/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root ";
195	>	samplename="SMS";
196	>	}
197	>	if(scantype==GMSB) {
198	>	copyfile << "/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root " << "local_storage/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root ";
199	>	samplename="GMSB";
200	>	}
201		stringstream newfilename;
202	<	if(scantype==mSUGRA) newfilename << "mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
203	<	if(scantype==SMS) newfilename << "SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
204	<	if(scantype==GMSB) newfilename << "GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
202	>	//    if(scantype==mSUGRA) newfilename << "local_storage/mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
203	>	if(scantype==mSUGRA) newfilename << "local_storage/mSUGRA_M0_" << mglu << "__M12_" << mlsp << ".root";
204	>	if(scantype==SMS) newfilename << "local_storage/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
205	>	if(scantype==GMSB) newfilename << "local_storage/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
206
207	<	if(!Contains(((scansample.collection)[(scansample.collection).size()-1]).filename,"_"+any2string(a)+"_"+any2string(b))||(a==0&&b==0)&&!initialized_t2) {
207	>	SUSYScanSpace::SUSYscantype=scantype;
208	>
209	>	if((scansample.collection).size()==0||!Contains(((scansample.collection)[(scansample.collection).size()-1]).filename,"_"+any2string(a)+"_"+any2string(b))||(a==0&&b==0)&&!initialized_t2) {
210		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;
211	<	if((scansample.collection).size()>1) {
211	>	while((scansample.collection).size()>0) {
212		scansample.RemoveLastSample();
213		}
214	+	dout << "New scanfileindex stems from scansample size: " << (scansample.collection).size() << endl;
215		scanfileindex=(scansample.collection).size();
216		//New: Loading file when necessary, not before (avoiding high memory usage and startup times)
217	<	if(scanfileindex!=0||isretry) {
218	<	delete_any_cached_scans();
219	<	dout << "Going to download the scan file with the following copy command: " << copyfile.str() << endl;
220	<	gSystem->Exec(copyfile.str().c_str());
217	>	if(1) {
218	>	if(!doesROOTFileExist(newfilename.str())) {
219	>	dout << "Going to download the scan file with the following copy command: " << copyfile.str() << endl;
220	>	int retcode = gSystem->Exec(copyfile.str().c_str()); // download it if it hasn't been downloaded before
221	>	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
222	>	}
223		if(doesROOTFileExist(newfilename.str())) {
224		initialized_t2=true;
225	<	scansample.AddSample(newfilename.str(),"scansample",1,1,false,true,scanfileindex,kRed);
225	>	dout << "Going to add the following sample: " << newfilename.str() << endl;
226	>	gSystem->Exec(("ls -ltrh "+newfilename.str()).c_str());
227	>	scansample.AddSample(newfilename.str(),samplename,1,1,false,true,scanfileindex,kRed);
228		} else {
229		srmcpretries++;
230		if(srmcpretries<5) {
231		dout << "The file could not be loaded correctly - retrying!" << endl;
232		sleep(5);
233	<	load_scan_sample(a,b,scanfileindex,scantype,true);
233	>	load_remote_scan_sample(a,b,scanfileindex,scantype,mglu,mlsp,true);
234		} else {
235		dout << "Have tried 5 times to load this sample. Giving up now and failing the program execution" << endl;
236	<	assert(0);
236	>	assert(srmcpretries<5);
237		dout << "The command " << copyfile.str() << " failed to copy the file to the local storage (saving it as " << newfilename.str() << ")" << endl;
238		}
239		}
#	Line 166 | Line 242 | void load_scan_sample(int a, int b, int
242		dout << "Last sample is the same as the current one. Recycling it." << endl;
243		scanfileindex=(scansample.collection).size()-1;
244		}
245	<
246	<	}
247	<	*/
245	>	}
246	>
247	>	void load_scan_sample(int a, int b, int &scanfileindex,int scantype,int mglu, int mlsp, bool isretry=false) {
248	>
249	>	// There is no need to define your sample here. That is done in Setup.C where you define the loading directory!
250
251	<	// WATCH OUT THIS LINE NEEDS TO BE REMOVED
252	<	scanfileindex=0;
251	>	SUSYScanSpace::SUSYscantype=scantype;
252	>	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;
253	>	char hostname[1023];
254	>	gethostname(hostname,1023);
255	>	if((Contains(hostname,"t3ui")||Contains(hostname,"t3wn"))) load_local_scan_sample(a,b,scanfileindex,scantype,mglu,mlsp,isretry);
256	>	else load_remote_scan_sample(a,b,scanfileindex,scantype,mglu,mlsp,isretry);
257		}
258
259		/*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 287 | void establish_SUSY_limits(string mcjzb,
287		while ((key = (TKey*)nextkey()))
288		{
289		TObject *obj = key->ReadObj();
290	<	if(Contains((string)(obj->GetName()),"mSUGRA")) scantype=mSUGRA;
290	>	if(Contains((string)(obj->GetName()),"SUGRA")) scantype=mSUGRA;
291		if(Contains((string)(obj->GetName()),"GMSB")) scantype=GMSB;
292		}
293
#	Line 215 | Line 297 | void establish_SUSY_limits(string mcjzb,
297		massgluname="M0"; // this is the "x axis" in the limit plot (like the gluino in the SMS case)
298		massLSPname="M12"; // this is the "y axis" in the limit plot (like the LSP in the SMS case)
299		mglustart=m0start;
300	<	xsec=getXsec(PlottingSetup::cbafbasedir+"/Plotting/Modules/external/scale_xsection_nlo1.0_m0_m12_10_0_1v1.txt");
300	>	xsec=getXsec(PlottingSetup::cbafbasedir+"/"+mSUGRAxsFile);
301	>
302		mgluend=m0end;
303		mglustep=m0step;
304		mLSPstart=m12start;
#	Line 223 | Line 306 | void establish_SUSY_limits(string mcjzb,
306		mLSPstep=m12step;
307		prefix="mSUGRA_";
308		dout << "mSUGRA scan has been set up." << endl;
226	–	//    xsec=getXsec("/scratch/buchmann/C/scale_xsection_nlo1.0_m0_m12_10_0_1v1.txt");
309		}
310		if(scantype==SMS) {
311		dout << "SMS scan has been set up." << endl;
#	Line 260 | Line 342 | void establish_SUSY_limits(string mcjzb,
342		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
343
344		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);
345		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);
346	+	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);
347
348
349		if(!fullerr || !mceff || !NEvents) {
#	Line 294 | Line 375 | void establish_SUSY_limits(string mcjzb,
375		continue;
376		}
377		vector<float> sigmas;
378	<	if(scantype!=mSUGRA) sigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped);
378	>	if(scantype!=mSUGRA) {
379	>	sigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped,true);
380	>	write_warning(__FUNCTION__,"Temporarily doing asymptotic limits!!!!");
381	>	}
382		else {
383		//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"!
384		vector<float> asigmas;
385		asigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped,true); // asymptotic limit first
386		float strength=asigmas[0]/XSmap->GetBinContent(GlobalBin);
387	<	if(strength>0.5&&strength<2) sigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped); // asymptotic limit first
388	<	else sigmas=asigmas;
387	>	/*      if(strength>0.5&&strength<2) {
388	>	sigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped); // asymptotic limit first
389	>	asymptoticmap->SetBinContent(GlobalBin,0);
390	>	} else {*/
391	>	sigmas=asigmas;
392	>	asymptoticmap->SetBinContent(GlobalBin,1);
393	>	/*      }*/
394		exclmap->SetBinContent(GlobalBin,strength);
306	–	xsmap->SetBinContent(GlobalBin,XSmap->GetBinContent(GlobalBin));
395		}
396
309	–
310	–
311	–
397		if(sigmas[0]>0) limitmap->SetBinContent(GlobalBin,sigmas[0]); //anything else is an error code
398		if(sigmas.size()>1) {
399		explimitmap->SetBinContent(GlobalBin,sigmas[1]);
#	Line 323 | Line 408 | void establish_SUSY_limits(string mcjzb,
408		}
409
410		prepare_scan_axis(limitmap,scantype);
411	+	gSystem->Exec("mkdir -p output");
412		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 :-)
413		outputfile->cd();
414		limitmap->Write();
415	+	XSmap->Write();
416		flipmap->Write();
417	+	asymptoticmap->Write();
418		if(doexpected) {
419		explimitmap->Write();
420		exp1plimitmap->Write();
#	Line 334 | Line 422 | void establish_SUSY_limits(string mcjzb,
422		exp2plimitmap->Write();
423		exp2mlimitmap->Write();
424		}
425	<	if(scantype==mSUGRA) {
338	<	exclmap->Write();
339	<	xsmap->Write();
340	<	totxsmap->Write();
341	<	}
425	>	if(scantype==mSUGRA) exclmap->Write();
426		outputfile->Close();
427		delete limcanvas;
428		}
#	Line 350 | Line 434 | void establish_SUSY_limits(string mcjzb,
434		write_error(__FUNCTION__,"You provided an invalid systematics file. Please change it ... ");
435		return;
436		}
437	<	for(int ibin=0;ibin<jzb_cut.size();ibin++) {
437	>	for(int ibin=0;ibin<(int)jzb_cut.size();ibin++) {
438		establish_SUSY_limits(mcjzb,datajzb,jzb_cut,requireZ, peakerror, fsyst, ibin,njobs, jobnumber);
439		}
440		}
441
442
443	+	string SimpleFormat(string input) {
444	+	string output;
445	+	for (unsigned int i=0; i<input.length();i++) {
446	+	if (isalnum(input[i])) //If current character is alpha-numeric
447	+	output+= input[i]; //Add the character to strConverted
448	+
449	+	}
450	+	return output;
451	+	}
452	+
453
454	+	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) {
455	+
456	+	if(shapeanalysis&&ibin>0) return; // pointless since we're doing a shape analysis :-)
457
361	–	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) {
458		bool runninglocally=true;
459		if(njobs>-1&&jobnumber>-1) {
460		runninglocally=false;
#	Line 370 | Line 466 | void scan_SUSY_parameter_space(string mc
466		}
467
468		jzbSel=jzb_cut[ibin];
469	+	if(shapeanalysis) jzbSel=0;
470		geqleq="geq";
471		automatized=true;
472		mcjzbexpression=mcjzb;
#	Line 380 | Line 477 | void scan_SUSY_parameter_space(string mc
477		string prefix="SMS_";
478		// up to here, everything is set up for SMS; now we need to switch stuff around if we're dealing with an mSUGRA scan!
479		int scantype=SMS;
480	<	if(Contains((scansample.collection)[0].samplename,"mSUGRA")) scantype=mSUGRA;
481	<	if(Contains((scansample.collection)[0].samplename,"GMSB")) scantype=GMSB;
480	>	if(Contains((PlottingSetup::ScanSampleDirectory),"SUGRA")) scantype=mSUGRA;
481	>	if(Contains((PlottingSetup::ScanSampleDirectory),"GMSB")) scantype=GMSB;
482
483		if(scantype==mSUGRA) {
484		massgluname="M0"; // this is the "x axis" in the limit plot (like the gluino in the SMS case)
#	Line 415 | Line 512 | void scan_SUSY_parameter_space(string mc
512
513		gStyle->SetPalette(100, MyPalette);
514
515	+	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
516	+	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
517	+
518		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
519		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
520		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 541 | mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLS
541		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);
542		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);
543		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);
544	+	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);
545	+	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);
546	+
547	+	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);
548	+
549	+
550	+	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);
551	+	// only used for asymptoticmap
552
553		TH2F *imposedxmap;
554		TH2F *realxmap;
#	Line 452 | Line 560 | mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLS
560
561		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);
562
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	–
563		if(ibin==0) {
564		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);
565		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 568 | mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLS
568		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);
569		}
570
571	+
572	+	if(shapeanalysis) jzbSel=jzb_cut[ibin];
573		vector<int> susy_Zdecay_origin;
574		vector<TH2F*> h_susy_Zdecay_origin;
575		// 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 591 | mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLS
591
592		write_warning(__FUNCTION__,"CURRENTLY SWITCHING AUTOMATIZED MODE OFF!");automatized=false;
593
486	–	float rightmargin=gStyle->GetPadRightMargin();
594		gStyle->SetPadRightMargin(0.15);
595
596		TCanvas *limcanvas = new TCanvas("limcanvas","Limit canvas");
#	Line 498 | Line 605 | write_warning(__FUNCTION__,"CURRENTLY SW
605
606		map <  pair<float, float>, map<string, float>  >  xsec;
607		if(scantype==mSUGRA) {
608	<	xsec=getXsec(PlottingSetup::mSUGRAxsFile);
608	>	xsec=getXsec(PlottingSetup::cbafbasedir+"/"+mSUGRAxsFile);
609		}
610
611		int ipoint=-1;
#	Line 508 | Line 615 | write_warning(__FUNCTION__,"CURRENTLY SW
615		ipoint++;
616		int GlobalBin=efficiencymap->FindBin(mglu,mlsp);
617		if(!runninglocally&&!do_this_point(ipoint,(int)Npoints,(int)jobnumber,(int)njobs)) continue;
618	+
619		int flipped=0;
620		float result=-987,resulterr=-987;
621		float m0trees = PlottingSetup::mgluend;
#	Line 515 | Line 623 | write_warning(__FUNCTION__,"CURRENTLY SW
623		int a = int((PlottingSetup::ScanXzones*mlsp)/(m12trees+1));
624		int b = int((PlottingSetup::ScanYzones*mglu)/(m0trees+1));
625		int scanfileindex=PlottingSetup::ScanYzones*a+b;
626	<	load_scan_sample(a,b,scanfileindex,scantype);
626	>	load_scan_sample(a,b,scanfileindex,scantype,mglu,mlsp);
627
628		clock_t start,finish;
629		start = clock(); // starting the clock to measure how long the computation takes!
#	Line 529 | Line 637 | write_warning(__FUNCTION__,"CURRENTLY SW
637		continue;
638		} else {
639		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;
640	<	if(PlottingSetup::RestrictToMassPeak==true) cout << "Analysis type: on-peak" << endl;
641	<	if(PlottingSetup::RestrictToMassPeak==false) cout << "Analysis type: offpeak" << endl;
640	>	if(PlottingSetup::RestrictToMassPeak==true) dout << "Analysis type: on-peak" << endl;
641	>	if(PlottingSetup::RestrictToMassPeak==false) dout << "Analysis type: offpeak" << endl;
642	>	if(shapeanalysis) dout << "Shape analysis? Yes." << endl;
643	>	else dout << "Shape analysis? No." << endl;
644		}
645		float hitrecord=0;
646		int hitrecordholder=0;
647		float hitsecondplace=0;
648		int hitsecondplaceholder=0;
649	<	for(int imo=0;imo<susy_Zdecay_origin.size()&&ibin==0;imo++) {
649	>	for(int imo=0;imo<(int)susy_Zdecay_origin.size()&&ibin==0;imo++) {
650		if(ibin>0) continue;
651		int MotherParticlePDGid = susy_Zdecay_origin[imo];
652		(scansample.collection)[scanfileindex].events->Draw("eventNum",("("+addcut.str()+")&&(abs(SourceOfZ-"+any2string(MotherParticlePDGid)+")<0.5)").c_str(),"goff");
#	Line 578 | Line 688 | write_warning(__FUNCTION__,"CURRENTLY SW
688		}
689
690
691	<	if(nevents!=0&&(efficiencyonly||systematicsonly)) {
691	>	if(scantype==mSUGRA) {
692	>	float absxs=0;
693	>	for(int i=0;i<12;i++) absxs+=GetXSecForPointAndChannel(mglu,mlsp,xsec,i);
694	>	TFile *FilterEffFile = new TFile(PlottingSetup::FilterEfficiencyFile.c_str());
695	>	TH2F *FilterEfficiency = (TH2F*)FilterEffFile->Get("FilterEfficiency");
696	>	float filtereff = FilterEfficiency->GetBinContent(FilterEfficiency->FindBin(mglu,mlsp));
697	>	FilterEff->SetBinContent(GlobalBin,filtereff);
698	>	FilterEffFile->Close();
699	>	absXSmap->SetBinContent(GlobalBin,absxs);
700	>	absxs*=filtereff;
701	>	XSmap->SetBinContent(GlobalBin,absxs);
702	>	}
703	>
704	>	if(nevents!=0) {
705		Value effwosigcont = MCefficiency((scansample.collection)[scanfileindex].events,result,resulterr,flipped,mcjzb,requireZ,(int)nevents,scantype,xsec,addcut.str(),-1);
706		if(result<0&&allowflipping) {
707		write_info(__FUNCTION__,"Found negative efficiency (signal distribution populating JZB<0 more than JZB>0) - flipping analysis!");
#	Line 615 | Line 738 | write_warning(__FUNCTION__,"CURRENTLY SW
738		Neventsmap->SetBinContent(GlobalBin,nevents);
739		ipointmap->SetBinContent(GlobalBin,ipoint);
740		if(efficiencyonly) continue;
741	<
742	<	do_systematics_for_one_file((scansample.collection)[scanfileindex].events,(int)nevents,"SUSY SCAN", systematics,flipped, xsec, mcjzb,datajzb,peakerror,requireZ, addcut.str(),true,scantype);
743	<	float JZBcutat = systematics[0][0];
744	<	float mceff    = systematics[0][1];
745	<	float mcefferr = systematics[0][2];//MC stat error
746	<	float toterr   = systematics[0][4];
747	<	float sys_jes  = systematics[0][5]; // Jet Energy Scale
748	<	float sys_jsu  = systematics[0][6]; // JZB scale uncertainty
749	<	float sys_res  = systematics[0][7]; // resolution
750	<	float mcwoscef = systematics[0][8]; // efficiency without signal contamination
751	<	float mcwoscefr= systematics[0][9]; // error on efficiency without signal contamination
752	<	float sys_pdf   = 0;
753	<	if(systematics[0].size()>10) sys_pdf = systematics[0][10]; // PDF
754	<
755	<	if(mceff!=mceff||toterr!=toterr||mceff<0 && (!systematicsonly&&!efficiencyonly)) {
756	<	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;
757	<	continue;
758	<	} else {
759	<	if(!systematicsonly&&!efficiencyonly) {
760	<	dout << "Calculating limit now for "<<massgluname<<"="<<mglu<<" , "<<massLSPname<<"="<<mlsp <<endl;
761	<	vector<float> sigmas;
762	<	string plotfilename=(string)(TString((scansample.collection)[scanfileindex].samplename)+TString(massgluname)+TString(any2string(mglu))+TString("__")+TString(massLSPname)+TString(any2string(mlsp))+TString(".png"));
763	<	sigmas=compute_one_upper_limit(mceff,toterr,ibin,mcjzb,plotfilename,true,flipped);
764	<	//      do_limit_wrapper(mceff,toterr,ibin,mcjzb,sigmas,plotfilename);
765	<	dout << "back in " << __FUNCTION__ << endl;
766	<	if(sigmas[0]>-0.5) { // negative sigmas are the error signature of do_limit_wrapper, so we want to exclude them.
767	<	limitmap->SetBinContent(GlobalBin,sigmas[0]);
768	<	if(sigmas.size()>1) {
769	<	explimitmap->SetBinContent(GlobalBin,sigmas[1]);
770	<	exp1plimitmap->SetBinContent(GlobalBin,sigmas[2]);
771	<	exp1mlimitmap->SetBinContent(GlobalBin,sigmas[3]);
772	<	exp2plimitmap->SetBinContent(GlobalBin,sigmas[4]);
773	<	exp2mlimitmap->SetBinContent(GlobalBin,sigmas[5]);
741	>
742	>	if(!shapeanalysis) {
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	>	dout << "Going to store: total error is: " << toterr << " and the relative tot error, i.e. total error / mc efficiency is " << toterr/mceff << endl;
756	>	if(mceff!=mceff||toterr!=toterr||mceff<0 && (!systematicsonly&&!efficiencyonly)) {
757	>	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;
758	>	continue;
759	>	} else {
760	>	// Systematics and efficiencies make sense (i.e. non-zero, and all numbers)
761	>	if(!systematicsonly&&!efficiencyonly) {
762	>	dout << "Calculating limit now for "<<massgluname<<"="<<mglu<<" , "<<massLSPname<<"="<<mlsp <<endl;
763	>	vector<float> sigmas;
764	>	string plotfilename=(string)(TString((scansample.collection)[scanfileindex].samplename)+TString(massgluname)+TString(any2string(mglu))+TString("__")+TString(massLSPname)+TString(any2string(mlsp))+TString(".png"));
765	>	sigmas=compute_one_upper_limit(mceff,toterr,ibin,mcjzb,plotfilename,true,flipped);
766	>	//          do_limit_wrapper(mceff,toterr,ibin,mcjzb,sigmas,plotfilename);
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]);
775	>	}
776	>	sysjesmap->SetBinContent(GlobalBin,sys_jes);
777	>	sysjsumap->SetBinContent(GlobalBin,sys_jsu);
778	>	sysresmap->SetBinContent(GlobalBin,sys_res);
779	>	syspdfmap->SetBinContent(GlobalBin,sys_pdf);
780	>	systotmap->SetBinContent(GlobalBin,toterr/mceff);//total relative (!) error
781	>	sysstatmap->SetBinContent(GlobalBin,mcefferr);//total relative (!) error
782	>	dout << "A limit has been added at " << sigmas[0] << " for m_{glu}="<<mglu << " and m_{lsp}="<<mlsp<<endl;
783	>	}//end of if sigma is positive
784	>	}//end of not systematics only condition
785	>	if(systematicsonly) {
786	>	sysjesmap->SetBinContent(GlobalBin,sys_jes);
787	>	sysjsumap->SetBinContent(GlobalBin,sys_jsu);
788	>	sysresmap->SetBinContent(GlobalBin,sys_res);
789	>	syspdfmap->SetBinContent(GlobalBin,sys_pdf);
790	>	systotmap->SetBinContent(GlobalBin,toterr/mceff);//total relative (!) error
791	>	sysstatmap->SetBinContent(GlobalBin,mcefferr);//total relative (!) error
792	>	}
793	>	}//efficiency is valid
794	>	}
795	>	if(shapeanalysis) {
796	>	//shape analysis
797	>	dout << "This is the shape analysis. Have fun with it!" << endl;
798	>	float quickmceff,quickmcefferr;
799	>	MCefficiency((scansample.collection)[scanfileindex].events,quickmceff,quickmcefferr,flipped,mcjzb,requireZ,(int)nevents,scantype,xsec,addcut.str(),-1);
800	>	if(quickmceff==0) {
801	>	write_warning(__FUNCTION__,"The MC efficiency is zero - need to skip this point.");
802	>	continue;
803	>	}
804	>	stringstream PointName;
805	>	PointName << massgluname << "_" << mglu << "__" << massLSPname << "_" << mlsp;
806	>	SUSYScanSpace::SavedMLSP=mlsp;
807	>	SUSYScanSpace::SavedMGlu=mglu;
808	>	SUSYScanSpace::SavedMLSPname=massLSPname;
809	>	SUSYScanSpace::SavedMGluname=massgluname;
810	>	bool asymptoticonly=true; // do asymptotic computation first for everything, then (way below) if necessary switch on complete computation
811	>	ShapeDroplet droplet = LimitsFromShapes(asymptoticonly, 0, (scansample.collection)[scanfileindex].events,addcut.str().c_str(),PointName.str(),mcjzb,datajzb,jzb_cut,peakerror,peakerror);
812	>	if(!isValidShapeResult(droplet)) {
813	>	write_error(__FUNCTION__,"Something went horribly wrong with the shape analysis - ignoring this point!!!");
814	>	//        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 ...
815	>	continue;
816	>	}
817	>
818	>	droplet = LimitsFromShapes(asymptoticonly, 25*droplet.expected, (scansample.collection)[scanfileindex].events,addcut.str().c_str(),PointName.str(),mcjzb,datajzb,jzb_cut,peakerror,peakerror);
819	>	if(!isValidShapeResult(droplet)) {
820	>	write_error(__FUNCTION__,"Something went horribly wrong with the shape analysis - ignoring this point !!!");
821	>	//        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 ...
822	>	continue;
823		}
824
825	<	sysjesmap->SetBinContent(GlobalBin,sys_jes);
826	<	sysjsumap->SetBinContent(GlobalBin,sys_jsu);
827	<	sysresmap->SetBinContent(GlobalBin,sys_res);
828	<	syspdfmap->SetBinContent(GlobalBin,sys_pdf);
829	<	systotmap->SetBinContent(GlobalBin,toterr/mceff);//total relative (!) error
830	<	sysstatmap->SetBinContent(GlobalBin,mcefferr);//total relative (!) error
831	<	dout << "A limit has been added at " << sigmas[0] << " for m_{glu}="<<mglu << " and m_{lsp}="<<mlsp<<endl;
832	<	} //end of if sigma is positive
833	<
834	<	//end of not systematics only condition
835	<	}
836	<	if(systematicsonly) {
837	<	sysjesmap->SetBinContent(GlobalBin,sys_jes);
838	<	sysjsumap->SetBinContent(GlobalBin,sys_jsu);
839	<	sysresmap->SetBinContent(GlobalBin,sys_res);
840	<	syspdfmap->SetBinContent(GlobalBin,sys_pdf);
841	<	systotmap->SetBinContent(GlobalBin,toterr/mceff);//total relative (!) error
842	<	sysstatmap->SetBinContent(GlobalBin,mcefferr);//total relative (!) error
843	<	}
844	<	}//efficiency is valid
825	>	float referenceXS=0;
826	>	if(scantype==mSUGRA) referenceXS=XSmap->GetBinContent(GlobalBin);
827	>	else referenceXS=getSMSxs(mlsp,mglu);
828	>	float excludedXSobs = droplet.observed/droplet.SignalIntegral;
829	>	absXSmap->SetBinContent(GlobalBin,referenceXS);
830	>	write_info(__FUNCTION__,"Established limit at "+any2string(excludedXSobs)+" , ref XS is "+any2string(referenceXS));
831	>	asymptoticmap->SetBinContent(GlobalBin,1);
832	>	if(excludedXSobs>0.5*referenceXS && excludedXSobs<2.5*referenceXS||scantype!=mSUGRA) {
833	>	// for mSUGRA we only do an in-depth scan of a band around the exclusion, for SMSes we want to know the full CLs UL everywhere
834	>	droplet = LimitsFromShapes(false, droplet.expected, (scansample.collection)[scanfileindex].events,addcut.str().c_str(),PointName.str(),mcjzb,datajzb,jzb_cut,peakerror,peakerror);
835	>	if(!isValidShapeResult(droplet)) {
836	>	write_error(__FUNCTION__,"Something went horribly wrong with the shape analysis - ignoring this point !!!");
837	>	//          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 ...
838	>	continue;
839	>	}
840	>	asymptoticmap->SetBinContent(GlobalBin,0);
841	>	}
842	>	rawlimitmap->SetBinContent(GlobalBin,droplet.observed);
843	>	rawelimitmap->SetBinContent(GlobalBin,droplet.expected);
844	>	limitmap->SetBinContent(GlobalBin,droplet.observed/droplet.SignalIntegral);
845	>	explimitmap->SetBinContent(GlobalBin,droplet.expected/droplet.SignalIntegral);
846	>	exp1plimitmap->SetBinContent(GlobalBin,droplet.expectedPlus1Sigma/droplet.SignalIntegral);
847	>	exp1mlimitmap->SetBinContent(GlobalBin,droplet.expectedMinus1Sigma/droplet.SignalIntegral);
848	>	exp2plimitmap->SetBinContent(GlobalBin,droplet.expectedPlus2Sigma/droplet.SignalIntegral);
849	>	exp2mlimitmap->SetBinContent(GlobalBin,droplet.expectedMinus2Sigma/droplet.SignalIntegral);
850	>	dout <<"Saved the following shape limit: " << droplet.observed/droplet.SignalIntegral << " (observed) , " << droplet.expected/droplet.SignalIntegral << " (expected) " << endl;
851	>	sysjesmap->SetBinContent(GlobalBin,droplet.JES);
852	>	sysjsumap->SetBinContent(GlobalBin,droplet.JSU);
853	>	syspdfmap->SetBinContent(GlobalBin,droplet.PDF);
854	>	systotmap->SetBinContent(GlobalBin,droplet.toterr);
855	>	sysstatmap->SetBinContent(GlobalBin,droplet.staterr);
856	>	}
857		finish = clock();
858		timemap->SetBinContent(GlobalBin,((float(finish)-float(start))/CLOCKS_PER_SEC));
859		}
860		}
677	–
861		prepare_scan_axis(limitmap,scantype);
862		prepare_scan_axis(sysjesmap,scantype);
863		prepare_scan_axis(sysjsumap,scantype);
#	Line 685 | Line 868 | write_warning(__FUNCTION__,"CURRENTLY SW
868		prepare_scan_axis(timemap,scantype);
869
870		if(!systematicsonly&&!efficiencyonly) {
871	<	limcanvas->cd();
871	>	TObject* old=gDirectory->GetList()->FindObject("limcanvas");
872	>	if(old) gDirectory->GetList()->Remove(old);
873	>	TCanvas *limcanvas = new TCanvas("limcanvas","Limit canvas");
874		limitmap->Draw("COLZ");
875		string titletobewritten="Limits in LSP-Glu plane";
876		if(scantype==mSUGRA) titletobewritten="Limits in m_{1/2}-m_{0} plane";
#	Line 699 | Line 884 | write_warning(__FUNCTION__,"CURRENTLY SW
884		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 :-)
885		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
886		limitmap->Write();
887	+	rawlimitmap->Write();
888	+	rawelimitmap->Write();
889		if(doexpected) {
890		explimitmap->Write();
891		exp1plimitmap->Write();
#	Line 724 | Line 911 | write_warning(__FUNCTION__,"CURRENTLY SW
911		systotmap->Write();
912		sysstatmap->Write();
913		XSmap->Write();
914	+	absXSmap->Write();
915	+	FilterEff->Write();
916	+	if(shapeanalysis) asymptoticmap->Write();
917		if(ibin==0) imposedxmap->Write();
918		if(ibin==0) realxmap->Write();
919		Neventsmap->Write();
920	+	ScanIdentifier->Write();
921		ipointmap->Write();
922		timemap->Write();
923	<	for(int i=0;i<susy_Zdecay_origin.size();i++) {
923	>	for(int i=0;i<(int)susy_Zdecay_origin.size();i++) {
924		h_susy_Zdecay_origin[i]->Write();
925		}
926		outputfile->Close();
#	Line 782 | Line 973 | write_warning(__FUNCTION__,"CURRENTLY SW
973		noscefficiencymap->Write();
974		efficiencymapContam->Write();
975		Neventsmap->Write();
976	+	ScanIdentifier->Write();
977		ipointmap->Write();
978		syspdfmap->Write();
979		systotmap->Write();
980		sysstatmap->Write();
981	+	absXSmap->Write();
982		XSmap->Write();
983	+	FilterEff->Write();
984	+	if(shapeanalysis) asymptoticmap->Write();
985		if(ibin==0) imposedxmap->Write();
986		if(ibin==0) realxmap->Write();
987		timemap->Write();
988	<	for(int i=0;i<susy_Zdecay_origin.size();i++) {
988	>	for(int i=0;i<(int)susy_Zdecay_origin.size();i++) {
989		h_susy_Zdecay_origin[i]->Write();
990		}
991		outputfile->Close();
#	Line 857 | Line 1052 | write_warning(__FUNCTION__,"CURRENTLY SW
1052		TFile *outputfile=new TFile(("output/DistributedSystematics_job"+string(any2string(jobnumber))+"_of_"+string(any2string(njobs))+".root").c_str(),"UPDATE");
1053		ipointmap->Write();
1054		Neventsmap->Write();
1055	+	ScanIdentifier->Write();
1056		noscefficiencymap->Write();
1057		efficiencymapContam->Write();
1058		efficiencymap->Write();
#	Line 872 | Line 1068 | write_warning(__FUNCTION__,"CURRENTLY SW
1068		outputfile->Close();
1069		}
1070		}//end of efficiencies only
1071	<
1072	<	for(int i=0;i<susy_Zdecay_origin.size();i++) {
1071	>
1072	>	for(int i=0;i<(int)susy_Zdecay_origin.size();i++) {
1073		delete h_susy_Zdecay_origin[i];
1074		}
1075
#	Line 895 | Line 1091 | write_warning(__FUNCTION__,"CURRENTLY SW
1091		delete efficiencymapJets;
1092		delete efficiencymapSignal;
1093		delete Neventsmap;
1094	+	delete ScanIdentifier;
1095		delete ipointmap;
1096		delete syspdfmap;
1097		delete systotmap;
1098		delete sysstatmap;
1099		delete XSmap;
1100	<	delete limcanvas;
1100	>	delete absXSmap;
1101	>	TObject* old=gDirectory->GetList()->FindObject("limcanvas");
1102	>	if(old) gDirectory->GetList()->Remove(old);
1103		}
1104
1105	<	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) {
1105	>	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) {
1106		dout << "Starting the SUSY scan now with all " << jzb_cut.size() << " bin(s)" << endl;
1107	<	for(int ibin=0;ibin<jzb_cut.size();ibin++) {
1108	<	scan_SUSY_parameter_space(mcjzb,datajzb,jzb_cut,requireZ, peakerror, ibin, njobs, jobnumber,systonly,effonly);
1107	>	for(int ibin=0;ibin<(int)jzb_cut.size();ibin++) {
1108	>	scan_SUSY_parameter_space_for_one_cut(mcjzb,datajzb,jzb_cut,requireZ, peakerror, peakerrordata, ibin, njobs, jobnumber,systonly,effonly,shapeanalysis);
1109	>	if(shapeanalysis) scan_SUSY_parameter_space_for_one_cut(mcjzb,datajzb,jzb_cut,requireZ, peakerror, peakerrordata, ibin, njobs, jobnumber,systonly,effonly,!shapeanalysis); // if we run the shape analysis we also run cut & count - we pick the best result.
1110		}
1111	+	delete_any_cached_scans();// tidy up ...
1112		}
1113

Diff Legend

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