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.21 by buchmann, Mon Jul 9 12:32:45 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 "TTbar_stuff.C"
23	>	#include "ShapeDroplet.C"
24	>	#include "ShapeLimit.C"
25	>	#include "EdgeLimit.C"
26	>
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	<	// WATCH OUT THIS LINE NEEDS TO BE REMOVED
251	<	scanfileindex=0;
250	>	// There is no need to define your sample here. That is done in Setup.C where you define the loading directory!
251	>
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 238 | Line 321 | void establish_SUSY_limits(string mcjzb,
321		}
322
323		set_SUSY_style();
324	+	std::cout << __LINE__ << std::endl;
325
326		TCanvas *limcanvas = new TCanvas("limcanvas","Limit canvas");
327
#	Line 260 | Line 344 | void establish_SUSY_limits(string mcjzb,
344		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
345
346		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);
347		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);
348	+	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);
349
350
351		if(!fullerr || !mceff || !NEvents) {
#	Line 294 | Line 377 | void establish_SUSY_limits(string mcjzb,
377		continue;
378		}
379		vector<float> sigmas;
380	<	if(scantype!=mSUGRA) sigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped);
380	>	if(scantype!=mSUGRA) {
381	>	sigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped,true);
382	>	write_warning(__FUNCTION__,"Temporarily doing asymptotic limits!!!!");
383	>	}
384		else {
385		//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"!
386		vector<float> asigmas;
387		asigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped,true); // asymptotic limit first
388		float strength=asigmas[0]/XSmap->GetBinContent(GlobalBin);
389	<	if(strength>0.5&&strength<2) sigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped); // asymptotic limit first
390	<	else sigmas=asigmas;
389	>	/*      if(strength>0.5&&strength<2) {
390	>	sigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped); // asymptotic limit first
391	>	asymptoticmap->SetBinContent(GlobalBin,0);
392	>	} else {*/
393	>	sigmas=asigmas;
394	>	asymptoticmap->SetBinContent(GlobalBin,1);
395	>	/*      }*/
396		exclmap->SetBinContent(GlobalBin,strength);
306	–	xsmap->SetBinContent(GlobalBin,XSmap->GetBinContent(GlobalBin));
397		}
398
309	–
310	–
311	–
399		if(sigmas[0]>0) limitmap->SetBinContent(GlobalBin,sigmas[0]); //anything else is an error code
400		if(sigmas.size()>1) {
401		explimitmap->SetBinContent(GlobalBin,sigmas[1]);
#	Line 323 | Line 410 | void establish_SUSY_limits(string mcjzb,
410		}
411
412		prepare_scan_axis(limitmap,scantype);
413	+	gSystem->Exec("mkdir -p output");
414		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 :-)
415		outputfile->cd();
416		limitmap->Write();
417	+	XSmap->Write();
418		flipmap->Write();
419	+	asymptoticmap->Write();
420		if(doexpected) {
421		explimitmap->Write();
422		exp1plimitmap->Write();
#	Line 334 | Line 424 | void establish_SUSY_limits(string mcjzb,
424		exp2plimitmap->Write();
425		exp2mlimitmap->Write();
426		}
427	<	if(scantype==mSUGRA) {
338	<	exclmap->Write();
339	<	xsmap->Write();
340	<	totxsmap->Write();
341	<	}
427	>	if(scantype==mSUGRA) exclmap->Write();
428		outputfile->Close();
429		delete limcanvas;
430		}
#	Line 350 | Line 436 | void establish_SUSY_limits(string mcjzb,
436		write_error(__FUNCTION__,"You provided an invalid systematics file. Please change it ... ");
437		return;
438		}
439	<	for(int ibin=0;ibin<jzb_cut.size();ibin++) {
439	>	for(int ibin=0;ibin<(int)jzb_cut.size();ibin++) {
440		establish_SUSY_limits(mcjzb,datajzb,jzb_cut,requireZ, peakerror, fsyst, ibin,njobs, jobnumber);
441		}
442		}
443
444
445	+	string SimpleFormat(string input) {
446	+	string output;
447	+	for (unsigned int i=0; i<input.length();i++) {
448	+	if (isalnum(input[i])) //If current character is alpha-numeric
449	+	output+= input[i]; //Add the character to strConverted
450	+
451	+	}
452	+	return output;
453	+	}
454	+
455
456	+	vector<float> 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,string specialfilename="") {
457	+	//note: njobs==-2 means that we're running a benchmark point
458	+	dout << "Starting the SUSY scan now with " << jzb_cut[ibin] << " (bin #" << ibin << " )" << endl;
459	+	dout << "   mcjzb: " << mcjzb << endl;
460	+	dout << "   datajzb: " << datajzb << endl;
461	+	dout << "   jzb_cuts: ";
462	+	for(int ib=0;ib<(int)jzb_cut.size()-1;ib++) dout << jzb_cut[ib] << ", ";
463	+	dout << jzb_cut[jzb_cut.size()-1] << endl;
464	+	dout << "   requireZ: " << requireZ << endl;
465	+	dout << "   Peak error: " << peakerror << endl;
466	+	dout << "   Data Peak e: " << peakerrordata << endl;
467	+	dout << "   NJobs : " << njobs << endl;
468	+	dout << "   Job num: " << jobnumber << endl;
469	+	dout << "   Syst only? " << systematicsonly << endl;
470	+	dout << "   Eff only? " << efficiencyonly << endl;
471	+	dout << "   Shape analysis?"  << shapeanalysis << endl;
472	+	dout << "   Special file name? \"" << specialfilename<< "\"" << endl;
473	+
474	+	vector<float> empty;
475	+
476	+	if(shapeanalysis&&ibin>0) return empty; // pointless since we're doing a shape analysis :-)
477
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) {
478		bool runninglocally=true;
479		if(njobs>-1&&jobnumber>-1) {
480		runninglocally=false;
#	Line 366 | Line 482 | void scan_SUSY_parameter_space(string mc
482		dout << "Running on the GRID (this is job " << jobnumber << "/" << njobs << ") for a jzb cut at " << jzb_cut[ibin] << endl;
483		} else {
484		dout << "Running locally " << endl;
485	<	dout << "This will take a really, really long time - if  you want to see the results THIS month try running the LimitWorkerScript on the grid (DistributedModelCalculation/Limits/)" << endl;
485	>	if(njobs>-2) dout << "If you are running more than one point and you want to see the results THIS month try running scripts from DistributedModelCalculation" << endl;
486		}
487
488		jzbSel=jzb_cut[ibin];
489	+	if(shapeanalysis) jzbSel=0;
490		geqleq="geq";
491		automatized=true;
492		mcjzbexpression=mcjzb;
#	Line 380 | Line 497 | void scan_SUSY_parameter_space(string mc
497		string prefix="SMS_";
498		// up to here, everything is set up for SMS; now we need to switch stuff around if we're dealing with an mSUGRA scan!
499		int scantype=SMS;
500	<	if(Contains((scansample.collection)[0].samplename,"mSUGRA")) scantype=mSUGRA;
501	<	if(Contains((scansample.collection)[0].samplename,"GMSB")) scantype=GMSB;
500	>	write_info(__FUNCTION__,PlottingSetup::ScanSampleDirectory);
501	>	if(Contains((PlottingSetup::ScanSampleDirectory),"SUGRA")) scantype=mSUGRA;
502	>	if(Contains((PlottingSetup::ScanSampleDirectory),"GMSB")) scantype=GMSB;
503
504		if(scantype==mSUGRA) {
505		massgluname="M0"; // this is the "x axis" in the limit plot (like the gluino in the SMS case)
#	Line 415 | Line 533 | void scan_SUSY_parameter_space(string mc
533
534		gStyle->SetPalette(100, MyPalette);
535
536	+	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
537	+	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
538	+
539		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
540		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
541		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 562 | mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLS
562		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);
563		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);
564		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);
565	+	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);
566	+	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);
567	+
568	+	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);
569	+
570	+
571	+	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);
572	+	// only used for asymptoticmap
573
574		TH2F *imposedxmap;
575		TH2F *realxmap;
#	Line 452 | Line 581 | mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLS
581
582		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);
583
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	–
584		if(ibin==0) {
585		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);
586		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 589 | mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLS
589		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);
590		}
591
592	+
593	+	if(shapeanalysis) jzbSel=jzb_cut[ibin];
594		vector<int> susy_Zdecay_origin;
595		vector<TH2F*> h_susy_Zdecay_origin;
596		// 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 612 | mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLS
612
613		write_warning(__FUNCTION__,"CURRENTLY SWITCHING AUTOMATIZED MODE OFF!");automatized=false;
614
486	–	float rightmargin=gStyle->GetPadRightMargin();
615		gStyle->SetPadRightMargin(0.15);
616
617		TCanvas *limcanvas = new TCanvas("limcanvas","Limit canvas");
#	Line 498 | Line 626 | write_warning(__FUNCTION__,"CURRENTLY SW
626
627		map <  pair<float, float>, map<string, float>  >  xsec;
628		if(scantype==mSUGRA) {
629	<	xsec=getXsec(PlottingSetup::mSUGRAxsFile);
629	>	xsec=getXsec(PlottingSetup::cbafbasedir+"/"+mSUGRAxsFile);
630		}
631
632		int ipoint=-1;
#	Line 508 | Line 636 | write_warning(__FUNCTION__,"CURRENTLY SW
636		ipoint++;
637		int GlobalBin=efficiencymap->FindBin(mglu,mlsp);
638		if(!runninglocally&&!do_this_point(ipoint,(int)Npoints,(int)jobnumber,(int)njobs)) continue;
639	+
640		int flipped=0;
641		float result=-987,resulterr=-987;
642		float m0trees = PlottingSetup::mgluend;
#	Line 515 | Line 644 | write_warning(__FUNCTION__,"CURRENTLY SW
644		int a = int((PlottingSetup::ScanXzones*mlsp)/(m12trees+1));
645		int b = int((PlottingSetup::ScanYzones*mglu)/(m0trees+1));
646		int scanfileindex=PlottingSetup::ScanYzones*a+b;
647	<	load_scan_sample(a,b,scanfileindex,scantype);
647	>	if(njobs>-2) load_scan_sample(a,b,scanfileindex,scantype,mglu,mlsp);
648	>	else {
649	>	if((scansample.collection).size()>=1) scansample.RemoveLastSample();
650	>	scanfileindex=(scansample.collection).size();
651	>	scansample.AddSample(specialfilename,"LMpoint",1,1,false,true,scanfileindex,kRed);
652	>	}
653
654		clock_t start,finish;
655		start = clock(); // starting the clock to measure how long the computation takes!
656		stringstream addcut;
657		addcut << "(TMath::Abs("<<massgluname<<"-"<<mglu<<")<5)&&(TMath::Abs("<<massLSPname<<"-"<<mlsp<<")<5)";
658	+	if(njobs==-2) {
659	+	dout << "Processing a signal sample, not a scan : setting addcut to \"eventNum>0||eventNum<1\" (empty cut)" << endl;
660	+	addcut.str("(eventNum>0||eventNum<1)");
661	+	}
662	+
663		(scansample.collection)[scanfileindex].events->Draw("eventNum",addcut.str().c_str(),"goff");
664		float nevents = (scansample.collection)[scanfileindex].events->GetSelectedRows();
665		vector<vector<float> > systematics;
#	Line 529 | Line 668 | write_warning(__FUNCTION__,"CURRENTLY SW
668		continue;
669		} else {
670		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;
671	<	if(PlottingSetup::RestrictToMassPeak==true) cout << "Analysis type: on-peak" << endl;
672	<	if(PlottingSetup::RestrictToMassPeak==false) cout << "Analysis type: offpeak" << endl;
671	>	if(PlottingSetup::RestrictToMassPeak==true) dout << "Analysis type: on-peak" << endl;
672	>	if(PlottingSetup::RestrictToMassPeak==false) dout << "Analysis type: offpeak" << endl;
673	>	if(shapeanalysis) dout << "Shape analysis? Yes." << endl;
674	>	else dout << "Shape analysis? No." << endl;
675		}
676		float hitrecord=0;
677		int hitrecordholder=0;
678		float hitsecondplace=0;
679		int hitsecondplaceholder=0;
680	<	for(int imo=0;imo<susy_Zdecay_origin.size()&&ibin==0;imo++) {
680	>	for(int imo=0;imo<(int)susy_Zdecay_origin.size()&&ibin==0&&njobs>-2;imo++) {
681		if(ibin>0) continue;
682		int MotherParticlePDGid = susy_Zdecay_origin[imo];
683		(scansample.collection)[scanfileindex].events->Draw("eventNum",("("+addcut.str()+")&&(abs(SourceOfZ-"+any2string(MotherParticlePDGid)+")<0.5)").c_str(),"goff");
#	Line 559 | Line 700 | write_warning(__FUNCTION__,"CURRENTLY SW
700		SecondaryZSource->SetBinContent(GlobalBin,hitsecondplaceholder);
701		}
702
703	<	if(ibin==0) {
703	>	if(ibin==0&&njobs>-2) {
704		TH1F *realxhisto = new TH1F("realxhisto","realxhisto",1000,0,1);
705		(scansample.collection)[scanfileindex].events->Draw("realx>>realxhisto",addcut.str().c_str(),"goff");
706		realxmap->SetBinContent(GlobalBin,realxhisto->GetMean());
#	Line 578 | Line 719 | write_warning(__FUNCTION__,"CURRENTLY SW
719		}
720
721
722	<	if(nevents!=0&&(efficiencyonly||systematicsonly)) {
722	>	if(scantype==mSUGRA) {
723	>	float absxs=0;
724	>	for(int i=0;i<12;i++) absxs+=GetXSecForPointAndChannel(mglu,mlsp,xsec,i);
725	>	TFile *FilterEffFile = new TFile(PlottingSetup::FilterEfficiencyFile.c_str());
726	>	TH2F *FilterEfficiency = (TH2F*)FilterEffFile->Get("FilterEfficiency");
727	>	float filtereff = FilterEfficiency->GetBinContent(FilterEfficiency->FindBin(mglu,mlsp));
728	>	FilterEff->SetBinContent(GlobalBin,filtereff);
729	>	FilterEffFile->Close();
730	>	absXSmap->SetBinContent(GlobalBin,absxs);
731	>	absxs*=filtereff;
732	>	XSmap->SetBinContent(GlobalBin,absxs);
733	>	}
734	>
735	>	if(nevents!=0) {
736		Value effwosigcont = MCefficiency((scansample.collection)[scanfileindex].events,result,resulterr,flipped,mcjzb,requireZ,(int)nevents,scantype,xsec,addcut.str(),-1);
737		if(result<0&&allowflipping) {
738		write_info(__FUNCTION__,"Found negative efficiency (signal distribution populating JZB<0 more than JZB>0) - flipping analysis!");
#	Line 615 | Line 769 | write_warning(__FUNCTION__,"CURRENTLY SW
769		Neventsmap->SetBinContent(GlobalBin,nevents);
770		ipointmap->SetBinContent(GlobalBin,ipoint);
771		if(efficiencyonly) continue;
772	+
773	+	if(!shapeanalysis) {
774	+	do_systematics_for_one_file((scansample.collection)[scanfileindex].events,(int)nevents,"SUSY SCAN", systematics,flipped, xsec, mcjzb,datajzb,peakerror,requireZ, addcut.str(),true,scantype);
775	+	//      float JZBcutat = systematics[0][0];
776	+	float mceff    = systematics[0][1];
777	+	float mcefferr = systematics[0][2];//MC stat error
778	+	float toterr   = systematics[0][4];
779	+	float sys_jes  = systematics[0][5]; // Jet Energy Scale
780	+	float sys_jsu  = systematics[0][6]; // JZB scale uncertainty
781	+	float sys_res  = systematics[0][7]; // resolution
782	+	//      float mcwoscef = systematics[0][8]; // efficiency without signal contamination
783	+	//      float mcwoscefr= systematics[0][9]; // error on efficiency without signal contamination
784	+	float sys_pdf   = 0;
785	+	if(systematics[0].size()>10) sys_pdf = systematics[0][10]; // PDF
786	+	dout << "Going to store: total error is: " << toterr << " and the relative tot error, i.e. total error / mc efficiency is " << toterr/mceff << endl;
787	+	if(mceff!=mceff||toterr!=toterr||mceff<0 && (!systematicsonly&&!efficiencyonly)) {
788	+	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;
789	+	continue;
790	+	} else {
791	+	// Systematics and efficiencies make sense (i.e. non-zero, and all numbers)
792	+	if(!systematicsonly&&!efficiencyonly) {
793	+	dout << "Calculating limit now for "<<massgluname<<"="<<mglu<<" , "<<massLSPname<<"="<<mlsp <<endl;
794	+	vector<float> sigmas;
795	+	string plotfilename=(string)(TString((scansample.collection)[scanfileindex].samplename)+TString(massgluname)+TString(any2string(mglu))+TString("__")+TString(massLSPname)+TString(any2string(mlsp))+TString(".png"));
796	+	sigmas=compute_one_upper_limit(mceff,toterr,ibin,mcjzb,plotfilename,true,flipped);
797	+	//          do_limit_wrapper(mceff,toterr,ibin,mcjzb,sigmas,plotfilename);
798	+	if(sigmas[0]>-0.5) { // negative sigmas are the error signature of do_limit_wrapper, so we want to exclude them.
799	+	limitmap->SetBinContent(GlobalBin,sigmas[0]);
800	+	if(sigmas.size()>1) {
801	+	explimitmap->SetBinContent(GlobalBin,sigmas[1]);
802	+	exp1plimitmap->SetBinContent(GlobalBin,sigmas[2]);
803	+	exp1mlimitmap->SetBinContent(GlobalBin,sigmas[3]);
804	+	exp2plimitmap->SetBinContent(GlobalBin,sigmas[4]);
805	+	exp2mlimitmap->SetBinContent(GlobalBin,sigmas[5]);
806	+	}
807	+	sysjesmap->SetBinContent(GlobalBin,sys_jes);
808	+	sysjsumap->SetBinContent(GlobalBin,sys_jsu);
809	+	sysresmap->SetBinContent(GlobalBin,sys_res);
810	+	syspdfmap->SetBinContent(GlobalBin,sys_pdf);
811	+	systotmap->SetBinContent(GlobalBin,toterr/mceff);//total relative (!) error
812	+	sysstatmap->SetBinContent(GlobalBin,mcefferr);//total relative (!) error
813	+	dout << "A limit has been added at " << sigmas[0] << " for m_{glu}="<<mglu << " and m_{lsp}="<<mlsp<<endl;
814	+	}//end of if sigma is positive
815	+	}//end of not systematics only condition
816	+	if(systematicsonly) {
817	+	sysjesmap->SetBinContent(GlobalBin,sys_jes);
818	+	sysjsumap->SetBinContent(GlobalBin,sys_jsu);
819	+	sysresmap->SetBinContent(GlobalBin,sys_res);
820	+	syspdfmap->SetBinContent(GlobalBin,sys_pdf);
821	+	systotmap->SetBinContent(GlobalBin,toterr/mceff);//total relative (!) error
822	+	sysstatmap->SetBinContent(GlobalBin,mcefferr);//total relative (!) error
823	+	}
824	+	}//efficiency is valid
825	+	}
826	+	if(shapeanalysis) {
827	+	//shape analysis
828	+	dout << "This is the shape analysis. Have fun with it!" << endl;
829	+	float quickmceff,quickmcefferr;
830	+	MCefficiency((scansample.collection)[scanfileindex].events,quickmceff,quickmcefferr,flipped,mcjzb,requireZ,(int)nevents,scantype,xsec,addcut.str(),-1);
831	+	if(quickmceff==0) {
832	+	write_warning(__FUNCTION__,"The MC efficiency is zero - need to skip this point.");
833	+	continue;
834	+	}
835	+	stringstream PointName;
836	+	PointName << massgluname << "_" << mglu << "__" << massLSPname << "_" << mlsp;
837	+	SUSYScanSpace::SavedMLSP=mlsp;
838	+	SUSYScanSpace::SavedMGlu=mglu;
839	+	SUSYScanSpace::SavedMLSPname=massLSPname;
840	+	SUSYScanSpace::SavedMGluname=massgluname;
841	+	float referenceXS=0;
842	+	if(scantype==mSUGRA) referenceXS=XSmap->GetBinContent(GlobalBin);
843	+	else referenceXS=getSMSxs(mlsp,mglu);
844	+	float low_fullCLs=referenceXS*0.5;
845	+	float high_fullCLs=referenceXS*2.5;
846	+	if(scantype!=mSUGRA) {
847	+	low_fullCLs=-1;
848	+	high_fullCLs=10e10;
849	+	}
850	+	ShapeDroplet droplet = LimitsFromShapes(low_fullCLs,high_fullCLs, (scansample.collection)[scanfileindex].events,addcut.str().c_str(),PointName.str(),mcjzb,datajzb,jzb_cut,peakerror,peakerror);
851	+	if(!isValidShapeResult(droplet)) {
852	+	write_error(__FUNCTION__,"Something went horribly wrong with the shape analysis - ignoring this point!!!");
853	+	//        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 ...
854	+	if(njobs==-2) {
855	+	delete limcanvas;
856	+	return empty;
857	+	}
858	+	continue;
859
619	–	do_systematics_for_one_file((scansample.collection)[scanfileindex].events,(int)nevents,"SUSY SCAN", systematics,flipped, xsec, mcjzb,datajzb,peakerror,requireZ, addcut.str(),true,scantype);
620	–	float JZBcutat = systematics[0][0];
621	–	float mceff    = systematics[0][1];
622	–	float mcefferr = systematics[0][2];//MC stat error
623	–	float toterr   = systematics[0][4];
624	–	float sys_jes  = systematics[0][5]; // Jet Energy Scale
625	–	float sys_jsu  = systematics[0][6]; // JZB scale uncertainty
626	–	float sys_res  = systematics[0][7]; // resolution
627	–	float mcwoscef = systematics[0][8]; // efficiency without signal contamination
628	–	float mcwoscefr= systematics[0][9]; // error on efficiency without signal contamination
629	–	float sys_pdf   = 0;
630	–	if(systematics[0].size()>10) sys_pdf = systematics[0][10]; // PDF
631	–
632	–	if(mceff!=mceff||toterr!=toterr||mceff<0 && (!systematicsonly&&!efficiencyonly)) {
633	–	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;
634	–	continue;
635	–	} else {
636	–	if(!systematicsonly&&!efficiencyonly) {
637	–	dout << "Calculating limit now for "<<massgluname<<"="<<mglu<<" , "<<massLSPname<<"="<<mlsp <<endl;
638	–	vector<float> sigmas;
639	–	string plotfilename=(string)(TString((scansample.collection)[scanfileindex].samplename)+TString(massgluname)+TString(any2string(mglu))+TString("__")+TString(massLSPname)+TString(any2string(mlsp))+TString(".png"));
640	–	sigmas=compute_one_upper_limit(mceff,toterr,ibin,mcjzb,plotfilename,true,flipped);
641	–	//      do_limit_wrapper(mceff,toterr,ibin,mcjzb,sigmas,plotfilename);
642	–	dout << "back in " << __FUNCTION__ << endl;
643	–	if(sigmas[0]>-0.5) { // negative sigmas are the error signature of do_limit_wrapper, so we want to exclude them.
644	–	limitmap->SetBinContent(GlobalBin,sigmas[0]);
645	–	if(sigmas.size()>1) {
646	–	explimitmap->SetBinContent(GlobalBin,sigmas[1]);
647	–	exp1plimitmap->SetBinContent(GlobalBin,sigmas[2]);
648	–	exp1mlimitmap->SetBinContent(GlobalBin,sigmas[3]);
649	–	exp2plimitmap->SetBinContent(GlobalBin,sigmas[4]);
650	–	exp2mlimitmap->SetBinContent(GlobalBin,sigmas[5]);
860		}
861
862	<	sysjesmap->SetBinContent(GlobalBin,sys_jes);
863	<	sysjsumap->SetBinContent(GlobalBin,sys_jsu);
864	<	sysresmap->SetBinContent(GlobalBin,sys_res);
865	<	syspdfmap->SetBinContent(GlobalBin,sys_pdf);
866	<	systotmap->SetBinContent(GlobalBin,toterr/mceff);//total relative (!) error
867	<	sysstatmap->SetBinContent(GlobalBin,mcefferr);//total relative (!) error
868	<	dout << "A limit has been added at " << sigmas[0] << " for m_{glu}="<<mglu << " and m_{lsp}="<<mlsp<<endl;
869	<	} //end of if sigma is positive
870	<
871	<	//end of not systematics only condition
872	<	}
873	<	if(systematicsonly) {
874	<	sysjesmap->SetBinContent(GlobalBin,sys_jes);
875	<	sysjsumap->SetBinContent(GlobalBin,sys_jsu);
876	<	sysresmap->SetBinContent(GlobalBin,sys_res);
877	<	syspdfmap->SetBinContent(GlobalBin,sys_pdf);
878	<	systotmap->SetBinContent(GlobalBin,toterr/mceff);//total relative (!) error
879	<	sysstatmap->SetBinContent(GlobalBin,mcefferr);//total relative (!) error
880	<	}
881	<	}//efficiency is valid
862	>	float excludedXSobs = droplet.observed/droplet.SignalIntegral;
863	>	absXSmap->SetBinContent(GlobalBin,referenceXS);
864	>	write_info(__FUNCTION__,"Established limit at "+any2string(excludedXSobs)+" , ref XS is "+any2string(referenceXS));
865	>	if(njobs==-2) {
866	>	//this is the non-scan case
867	>	vector<float> sigmas;
868	>	sigmas.push_back(droplet.observed/droplet.SignalIntegral);
869	>	sigmas.push_back(droplet.expected/droplet.SignalIntegral);
870	>	sigmas.push_back(droplet.expectedPlus1Sigma/droplet.SignalIntegral);
871	>	sigmas.push_back(droplet.expectedMinus1Sigma/droplet.SignalIntegral);
872	>	sigmas.push_back(droplet.expectedPlus2Sigma/droplet.SignalIntegral);
873	>	sigmas.push_back(droplet.expectedMinus2Sigma/droplet.SignalIntegral);
874	>	delete limcanvas;
875	>	return sigmas;
876	>	}
877	>
878	>	rawlimitmap->SetBinContent(GlobalBin,droplet.observed);
879	>	rawelimitmap->SetBinContent(GlobalBin,droplet.expected);
880	>	limitmap->SetBinContent(GlobalBin,droplet.observed/droplet.SignalIntegral);
881	>	explimitmap->SetBinContent(GlobalBin,droplet.expected/droplet.SignalIntegral);
882	>	exp1plimitmap->SetBinContent(GlobalBin,droplet.expectedPlus1Sigma/droplet.SignalIntegral);
883	>	exp1mlimitmap->SetBinContent(GlobalBin,droplet.expectedMinus1Sigma/droplet.SignalIntegral);
884	>	exp2plimitmap->SetBinContent(GlobalBin,droplet.expectedPlus2Sigma/droplet.SignalIntegral);
885	>	exp2mlimitmap->SetBinContent(GlobalBin,droplet.expectedMinus2Sigma/droplet.SignalIntegral);
886	>	dout <<"Saved the following shape limit: " << droplet.observed/droplet.SignalIntegral << " (observed) , " << droplet.expected/droplet.SignalIntegral << " (expected) " << endl;
887	>	sysjesmap->SetBinContent(GlobalBin,droplet.JES);
888	>	sysjsumap->SetBinContent(GlobalBin,droplet.JSU);
889	>	syspdfmap->SetBinContent(GlobalBin,droplet.PDF);
890	>	systotmap->SetBinContent(GlobalBin,droplet.toterr);
891	>	sysstatmap->SetBinContent(GlobalBin,droplet.staterr);
892	>	}
893		finish = clock();
894		timemap->SetBinContent(GlobalBin,((float(finish)-float(start))/CLOCKS_PER_SEC));
895		}
896		}
677	–
897		prepare_scan_axis(limitmap,scantype);
898		prepare_scan_axis(sysjesmap,scantype);
899		prepare_scan_axis(sysjsumap,scantype);
#	Line 685 | Line 904 | write_warning(__FUNCTION__,"CURRENTLY SW
904		prepare_scan_axis(timemap,scantype);
905
906		if(!systematicsonly&&!efficiencyonly) {
907	<	limcanvas->cd();
907	>	TObject* old=gDirectory->GetList()->FindObject("limcanvas");
908	>	if(old) gDirectory->GetList()->Remove(old);
909	>	TCanvas *limcanvas = new TCanvas("limcanvas","Limit canvas");
910		limitmap->Draw("COLZ");
911		string titletobewritten="Limits in LSP-Glu plane";
912		if(scantype==mSUGRA) titletobewritten="Limits in m_{1/2}-m_{0} plane";
#	Line 699 | Line 920 | write_warning(__FUNCTION__,"CURRENTLY SW
920		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 :-)
921		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
922		limitmap->Write();
923	+	rawlimitmap->Write();
924	+	rawelimitmap->Write();
925		if(doexpected) {
926		explimitmap->Write();
927		exp1plimitmap->Write();
#	Line 724 | Line 947 | write_warning(__FUNCTION__,"CURRENTLY SW
947		systotmap->Write();
948		sysstatmap->Write();
949		XSmap->Write();
950	+	absXSmap->Write();
951	+	FilterEff->Write();
952	+	if(shapeanalysis) asymptoticmap->Write();
953		if(ibin==0) imposedxmap->Write();
954		if(ibin==0) realxmap->Write();
955		Neventsmap->Write();
956	+	ScanIdentifier->Write();
957		ipointmap->Write();
958		timemap->Write();
959	<	for(int i=0;i<susy_Zdecay_origin.size();i++) {
959	>	for(int i=0;i<(int)susy_Zdecay_origin.size();i++) {
960		h_susy_Zdecay_origin[i]->Write();
961		}
962		outputfile->Close();
#	Line 782 | Line 1009 | write_warning(__FUNCTION__,"CURRENTLY SW
1009		noscefficiencymap->Write();
1010		efficiencymapContam->Write();
1011		Neventsmap->Write();
1012	+	ScanIdentifier->Write();
1013		ipointmap->Write();
1014		syspdfmap->Write();
1015		systotmap->Write();
1016		sysstatmap->Write();
1017	+	absXSmap->Write();
1018		XSmap->Write();
1019	+	FilterEff->Write();
1020	+	if(shapeanalysis) asymptoticmap->Write();
1021		if(ibin==0) imposedxmap->Write();
1022		if(ibin==0) realxmap->Write();
1023		timemap->Write();
1024	<	for(int i=0;i<susy_Zdecay_origin.size();i++) {
1024	>	for(int i=0;i<(int)susy_Zdecay_origin.size();i++) {
1025		h_susy_Zdecay_origin[i]->Write();
1026		}
1027		outputfile->Close();
#	Line 857 | Line 1088 | write_warning(__FUNCTION__,"CURRENTLY SW
1088		TFile *outputfile=new TFile(("output/DistributedSystematics_job"+string(any2string(jobnumber))+"_of_"+string(any2string(njobs))+".root").c_str(),"UPDATE");
1089		ipointmap->Write();
1090		Neventsmap->Write();
1091	+	ScanIdentifier->Write();
1092		noscefficiencymap->Write();
1093		efficiencymapContam->Write();
1094		efficiencymap->Write();
#	Line 872 | Line 1104 | write_warning(__FUNCTION__,"CURRENTLY SW
1104		outputfile->Close();
1105		}
1106		}//end of efficiencies only
1107	<
1108	<	for(int i=0;i<susy_Zdecay_origin.size();i++) {
1107	>
1108	>	for(int i=0;i<(int)susy_Zdecay_origin.size();i++) {
1109		delete h_susy_Zdecay_origin[i];
1110		}
1111
#	Line 895 | Line 1127 | write_warning(__FUNCTION__,"CURRENTLY SW
1127		delete efficiencymapJets;
1128		delete efficiencymapSignal;
1129		delete Neventsmap;
1130	+	delete ScanIdentifier;
1131		delete ipointmap;
1132		delete syspdfmap;
1133		delete systotmap;
1134		delete sysstatmap;
1135		delete XSmap;
1136	<	delete limcanvas;
1136	>	delete absXSmap;
1137	>	TObject* old=gDirectory->GetList()->FindObject("limcanvas");
1138	>	if(old) gDirectory->GetList()->Remove(old);
1139	>	return empty;
1140		}
1141
1142	<	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) {
1142	>	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) {
1143		dout << "Starting the SUSY scan now with all " << jzb_cut.size() << " bin(s)" << endl;
1144	<	for(int ibin=0;ibin<jzb_cut.size();ibin++) {
1145	<	scan_SUSY_parameter_space(mcjzb,datajzb,jzb_cut,requireZ, peakerror, ibin, njobs, jobnumber,systonly,effonly);
1144	>	dout << "   mcjzb: " << mcjzb << endl;
1145	>	dout << "   datajzb: " << datajzb << endl;
1146	>	dout << "   jzb_cut: ";
1147	>	for(int ib=0;ib<(int)jzb_cut.size()-1;ib++) dout << jzb_cut[ib] << ", ";
1148	>	dout << jzb_cut[jzb_cut.size()-1] << endl;
1149	>	dout << "   requireZ: " << requireZ << endl;
1150	>	dout << "   Peak error: " << peakerror << endl;
1151	>	dout << "   Data Peak e: " << peakerrordata << endl;
1152	>	dout << "   NJobs : " << njobs << endl;
1153	>	dout << "   Job num: " << jobnumber << endl;
1154	>	dout << "   Syst only?" << systonly << endl;
1155	>	dout << "   Eff only? " << effonly << endl;
1156	>	dout << "   Shape analysis?" << shapeanalysis << endl;
1157	>
1158	>	for(int ibin=0;ibin<(int)jzb_cut.size();ibin++) {
1159	>	scan_SUSY_parameter_space_for_one_cut(mcjzb,datajzb,jzb_cut,requireZ, peakerror, peakerrordata, ibin, njobs, jobnumber,systonly,effonly,shapeanalysis);
1160		}
1161	+	delete_any_cached_scans();// tidy up ...
1162		}
1163
1164	+
1165	+
1166	+	void compute_upper_limits_from_shapes(string mcjzb,string datajzb,vector<float> jzb_cut,bool requireZ, float peakerror, float peakerrordata) {
1167	+	dout << "Starting the SUSY scan now with all " << jzb_cut.size() << " bin(s)" << endl;
1168	+	vector<vector<float> > limits;
1169	+	vector<string> samplename;
1170	+	for(int isample=0;isample<(int)signalsamples.collection.size();isample++) {
1171	+	limits.push_back(scan_SUSY_parameter_space_for_one_cut(mcjzb,datajzb,jzb_cut,requireZ, peakerror, peakerrordata, 0, -2, -2,false,false,true,signalsamples.collection[isample].filename));// the "-2" stands for "load from customized path"
1172	+	}
1173	+	for(int isample=0;isample<(int)signalsamples.collection.size();isample++) {
1174	+	dout << signalsamples.collection[isample].samplename << " :" << endl;
1175	+	dout << "       observed :   " << limits[isample][0] << endl;
1176	+	dout << "       expected :   " << limits[isample][1] << endl;
1177	+	dout << "       exp. band:   [ " << limits[isample][5] << " , " << limits[isample][3] << " , \033[1;34m " << limits[isample][1] << "\033[0m , " << limits[isample][2] << " , " << limits[isample][4] << "] " << endl;
1178	+	}
1179	+	}

Diff Legend

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