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.11 by buchmann, Wed Apr 18 09:22:22 2012 UTC

#	Line 19 | Line 19
19		#include <TProfile.h>
20		#include <TKey.h>
21
22	+	#include "ShapeDroplet.C"
23	+	#include "ShapeLimit.C"
24	+
25		//#include "TTbar_stuff.C"
26		using namespace std;
27
#	Line 67 | Line 70 | void set_SUSY_style() {
70		}
71
72		bool delete_any_cached_scans() {
73	+	char hostname[1023];
74	+	gethostname(hostname,1023);
75	+
76	+
77	+	/*
78		//This should only be called via CRAB
79		cout << "   Deleting all cached files." << endl;
80		gSystem->Exec("ls -ltrh | grep root");
81	<	char hostname[1023];
82	<	gethostname(hostname,1023);
75	<	if((Contains(hostname,"t3ui")||Contains(hostname,"t3wn"))) return false;
76	<	else {
81	>	if(!((Contains(hostname,"t3ui")||Contains(hostname,"t3wn"))))
82	>	{
83		vector<string> all_files;
84		char currentpath[1024];
85		TString directory=TString(getcwd(currentpath,1024));
#	Line 86 | Line 92 | bool delete_any_cached_scans() {
92		gSystem->Exec(((string)"rm "+all_files[ifile]).c_str());
93		}
94		}
95	<	return true;
95	>	return true;
96	>	}
97	>	*/
98	>	if(!((Contains(hostname,"t3ui")||Contains(hostname,"t3wn")))) {
99	>	cout << "Going to purge files in local_storage that have been loaded!" << endl;
100	>	for(int i=0;i<SUSYScanSpace::loaded_files.size();i++) {
101	>	gSystem->Exec(((string)"rm "+SUSYScanSpace::loaded_files[i]).c_str());
102	>	dout << "      Purging : Deleted file " << SUSYScanSpace::loaded_files[i] << endl;
103	>	}
104		}
105	+	return true;
106		}
107
108		bool initialized_t2=false;
109
110		int srmcpretries=0;
111
112	<	void load_scan_sample(int a, int b, int &scanfileindex,int scantype,bool isretry=false) {
113	<	/*
114	<	// 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
112	>	void load_local_scan_sample(int a, int b, int &scanfileindex,int scantype,int mglu, int mlsp, bool isretry=false) {
113	>	stringstream filetoload;
114	>	string samplename;
115		filetoload << "/shome/buchmann/ntuples/"<<PlottingSetup::ScanSampleDirectory;
116	<	if(scantype==mSUGRA) filetoload << "/mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
117	<	if(scantype==SMS) filetoload << "/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
118	<	if(scantype==GMSB) filetoload << "/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
116	>	if(scantype==mSUGRA) {
117	>	//filetoload << "/mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
118	>	filetoload.str("");
119	>	filetoload << "/shome/lbaeni/jzb/" << PlottingSetup::ScanSampleDirectory  << "/mSUGRA_M0_" << mglu << "__M12_" << mlsp << ".root";
120	>	samplename="mSUGRA";
121	>	//      filetoload << "/mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
122	>	}
123	>	if(scantype==SMS) {
124	>	filetoload << "/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
125	>	samplename="SMS";
126	>	}
127	>	if(scantype==GMSB) {
128	>	filetoload << "/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
129	>	samplename="GMSB";
130	>	}
131	>
132	>	if(scantype==SMS && Contains(PlottingSetup::ScanSampleDirectory,"T1lh")) {
133	>	filetoload.str("");
134	>	filetoload << "/shome/lbaeni/jzb/" << PlottingSetup::ScanSampleDirectory  << "/SMS_MassGlu_" << mglu << "__MassLSP_" << mlsp << ".root";
135	>	samplename="SMS";
136	>	}
137	>
138		if(scansample.collection.size()<1||!Contains(((scansample.collection)[(scansample.collection).size()-1]).filename,"_"+any2string(a)+"_"+any2string(b))) {
139		dout << "The last sample is NOT the same one as the current one, possibly popping off last one and adding the new one." << endl;
140	<	if((scansample.collection).size()>1) {
140	>	if((scansample.collection).size()>=1) {
141		scansample.RemoveLastSample();
142		}
143		scanfileindex=(scansample.collection).size();
144		//New: Loading file when necessary, not before (avoiding high memory usage and startup times)
145	<	if(scanfileindex!=0) scansample.AddSample(filetoload.str(),"scansample",1,1,false,true,scanfileindex,kRed);
145	>
146	>	scansample.AddSample(filetoload.str(),samplename,1,1,false,true,scanfileindex,kRed);
147	>	dout << "Just added the following file: " << filetoload.str() << endl;
148		} else {
149		dout << "Last sample is the same as the current one. Recycling it." << endl;
150		scanfileindex=(scansample.collection).size()-1;
151		}
152		dout << " Going to use the following file: " << filetoload.str() << endl;
153	<	} // end of t3 case
154	<	else {
153	>	}
154	>
155	>	void load_remote_scan_sample(int a, int b, int &scanfileindex,int scantype,int mglu, int mlsp, bool isretry=false) {
156		write_info(__FUNCTION__,"Hello, CRAB! Might need to load files ...");
157	+	string samplename;
158		PlottingSetup::limitpatience=PlottingSetup::limitpatienceCRAB;
159		stringstream copyfile;
160	+	gSystem->Exec("mkdir -p local_storage");
161		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;
162	<	if(scantype==mSUGRA) copyfile << "/mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root " << "mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root ";
163	<	if(scantype==SMS) copyfile << "/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root " << "SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root ";
164	<	if(scantype==GMSB) copyfile << "/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root " << "GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root ";
162	>	if(scantype==mSUGRA) {
163	>	copyfile<< "/mSUGRA_M0_" << mglu << "__M12_" << mlsp << ".root local_storage/mSUGRA_M0_" << mglu << "__M12_" << mlsp << ".root ";
164	>	samplename="mSUGRA";
165	>	}
166	>	if(scantype==SMS) {
167	>	copyfile << "/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root " << "local_storage/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root ";
168	>	samplename="SMS";
169	>	}
170	>	if(scantype==GMSB) {
171	>	copyfile << "/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root " << "local_storage/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root ";
172	>	samplename="GMSB";
173	>	}
174		stringstream newfilename;
175	<	if(scantype==mSUGRA) newfilename << "mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
176	<	if(scantype==SMS) newfilename << "SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
177	<	if(scantype==GMSB) newfilename << "GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
175	>	//    if(scantype==mSUGRA) newfilename << "local_storage/mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
176	>	if(scantype==mSUGRA) newfilename << "local_storage/mSUGRA_M0_" << mglu << "__M12_" << mlsp << ".root";
177	>	if(scantype==SMS) newfilename << "local_storage/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
178	>	if(scantype==GMSB) newfilename << "local_storage/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
179
180	<	if(!Contains(((scansample.collection)[(scansample.collection).size()-1]).filename,"_"+any2string(a)+"_"+any2string(b))||(a==0&&b==0)&&!initialized_t2) {
180	>	SUSYScanSpace::SUSYscantype=scantype;
181	>
182	>	if((scansample.collection).size()==0||!Contains(((scansample.collection)[(scansample.collection).size()-1]).filename,"_"+any2string(a)+"_"+any2string(b))||(a==0&&b==0)&&!initialized_t2) {
183		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;
184	<	if((scansample.collection).size()>1) {
184	>	while((scansample.collection).size()>0) {
185		scansample.RemoveLastSample();
186		}
187	+	cout << "New scanfileindex stems from scansample size: " << (scansample.collection).size() << endl;
188		scanfileindex=(scansample.collection).size();
189		//New: Loading file when necessary, not before (avoiding high memory usage and startup times)
190	<	if(scanfileindex!=0||isretry) {
191	<	delete_any_cached_scans();
192	<	dout << "Going to download the scan file with the following copy command: " << copyfile.str() << endl;
193	<	gSystem->Exec(copyfile.str().c_str());
190	>	if(1) {
191	>	if(!doesROOTFileExist(newfilename.str())) {
192	>	dout << "Going to download the scan file with the following copy command: " << copyfile.str() << endl;
193	>	int retcode = gSystem->Exec(copyfile.str().c_str()); // download it if it hasn't been downloaded before
194	>	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
195	>	}
196		if(doesROOTFileExist(newfilename.str())) {
197		initialized_t2=true;
198	<	scansample.AddSample(newfilename.str(),"scansample",1,1,false,true,scanfileindex,kRed);
198	>	cout << "Going to add the following sample: " << newfilename.str() << endl;
199	>	gSystem->Exec(("ls -ltrh "+newfilename.str()).c_str());
200	>	scansample.AddSample(newfilename.str(),samplename,1,1,false,true,scanfileindex,kRed);
201		} else {
202		srmcpretries++;
203		if(srmcpretries<5) {
204		dout << "The file could not be loaded correctly - retrying!" << endl;
205		sleep(5);
206	<	load_scan_sample(a,b,scanfileindex,scantype,true);
206	>	load_remote_scan_sample(a,b,scanfileindex,scantype,mglu,mlsp,true);
207		} else {
208		dout << "Have tried 5 times to load this sample. Giving up now and failing the program execution" << endl;
209	<	assert(0);
209	>	assert(srmcpretries<5);
210		dout << "The command " << copyfile.str() << " failed to copy the file to the local storage (saving it as " << newfilename.str() << ")" << endl;
211		}
212		}
#	Line 166 | Line 215 | void load_scan_sample(int a, int b, int
215		dout << "Last sample is the same as the current one. Recycling it." << endl;
216		scanfileindex=(scansample.collection).size()-1;
217		}
218	<
219	<	}
220	<	*/
218	>	}
219	>
220	>	void load_scan_sample(int a, int b, int &scanfileindex,int scantype,int mglu, int mlsp, bool isretry=false) {
221
222	<	// WATCH OUT THIS LINE NEEDS TO BE REMOVED
223	<	scanfileindex=0;
222	>	// There is no need to define your sample here. That is done in Setup.C where you define the loading directory!
223	>
224	>	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;
225	>	char hostname[1023];
226	>	gethostname(hostname,1023);
227	>	if((Contains(hostname,"t3ui")||Contains(hostname,"t3wn"))) load_local_scan_sample(a,b,scanfileindex,scantype,mglu,mlsp,isretry);
228	>	else load_remote_scan_sample(a,b,scanfileindex,scantype,mglu,mlsp,isretry);
229		}
230
231		/*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 259 | void establish_SUSY_limits(string mcjzb,
259		while ((key = (TKey*)nextkey()))
260		{
261		TObject *obj = key->ReadObj();
262	<	if(Contains((string)(obj->GetName()),"mSUGRA")) scantype=mSUGRA;
262	>	if(Contains((string)(obj->GetName()),"SUGRA")) scantype=mSUGRA;
263		if(Contains((string)(obj->GetName()),"GMSB")) scantype=GMSB;
264		}
265
#	Line 215 | Line 269 | void establish_SUSY_limits(string mcjzb,
269		massgluname="M0"; // this is the "x axis" in the limit plot (like the gluino in the SMS case)
270		massLSPname="M12"; // this is the "y axis" in the limit plot (like the LSP in the SMS case)
271		mglustart=m0start;
272	<	xsec=getXsec(PlottingSetup::cbafbasedir+"/Plotting/Modules/external/scale_xsection_nlo1.0_m0_m12_10_0_1v1.txt");
272	>	xsec=getXsec(PlottingSetup::cbafbasedir+"/"+mSUGRAxsFile);
273	>
274		mgluend=m0end;
275		mglustep=m0step;
276		mLSPstart=m12start;
#	Line 223 | Line 278 | void establish_SUSY_limits(string mcjzb,
278		mLSPstep=m12step;
279		prefix="mSUGRA_";
280		dout << "mSUGRA scan has been set up." << endl;
226	–	//    xsec=getXsec("/scratch/buchmann/C/scale_xsection_nlo1.0_m0_m12_10_0_1v1.txt");
281		}
282		if(scantype==SMS) {
283		dout << "SMS scan has been set up." << endl;
#	Line 260 | Line 314 | void establish_SUSY_limits(string mcjzb,
314		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
315
316		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);
317		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);
318	+	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);
319
320
321		if(!fullerr || !mceff || !NEvents) {
#	Line 294 | Line 347 | void establish_SUSY_limits(string mcjzb,
347		continue;
348		}
349		vector<float> sigmas;
350	<	if(scantype!=mSUGRA) sigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped);
350	>	if(scantype!=mSUGRA) {
351	>	sigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped,true);
352	>	write_warning(__FUNCTION__,"Temporarily doing asymptotic limits!!!!");
353	>	}
354		else {
355		//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"!
356		vector<float> asigmas;
357		asigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped,true); // asymptotic limit first
358		float strength=asigmas[0]/XSmap->GetBinContent(GlobalBin);
359	<	if(strength>0.5&&strength<2) sigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped); // asymptotic limit first
360	<	else sigmas=asigmas;
359	>	/*      if(strength>0.5&&strength<2) {
360	>	sigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped); // asymptotic limit first
361	>	asymptoticmap->SetBinContent(GlobalBin,0);
362	>	} else {*/
363	>	sigmas=asigmas;
364	>	asymptoticmap->SetBinContent(GlobalBin,1);
365	>	/*      }*/
366		exclmap->SetBinContent(GlobalBin,strength);
306	–	xsmap->SetBinContent(GlobalBin,XSmap->GetBinContent(GlobalBin));
367		}
368
309	–
310	–
311	–
369		if(sigmas[0]>0) limitmap->SetBinContent(GlobalBin,sigmas[0]); //anything else is an error code
370		if(sigmas.size()>1) {
371		explimitmap->SetBinContent(GlobalBin,sigmas[1]);
#	Line 323 | Line 380 | void establish_SUSY_limits(string mcjzb,
380		}
381
382		prepare_scan_axis(limitmap,scantype);
383	+	gSystem->Exec("mkdir -p output");
384		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 :-)
385		outputfile->cd();
386		limitmap->Write();
387	+	XSmap->Write();
388		flipmap->Write();
389	+	asymptoticmap->Write();
390		if(doexpected) {
391		explimitmap->Write();
392		exp1plimitmap->Write();
#	Line 334 | Line 394 | void establish_SUSY_limits(string mcjzb,
394		exp2plimitmap->Write();
395		exp2mlimitmap->Write();
396		}
397	<	if(scantype==mSUGRA) {
338	<	exclmap->Write();
339	<	xsmap->Write();
340	<	totxsmap->Write();
341	<	}
397	>	if(scantype==mSUGRA) exclmap->Write();
398		outputfile->Close();
399		delete limcanvas;
400		}
#	Line 358 | Line 414 | void establish_SUSY_limits(string mcjzb,
414
415
416
417	<	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) {
417	>	void scan_SUSY_parameter_space(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) {
418	>
419	>	if(shapeanalysis&&ibin>0) return; // pointless since we're doing a shape analysis :-)
420	>
421		bool runninglocally=true;
422		if(njobs>-1&&jobnumber>-1) {
423		runninglocally=false;
#	Line 380 | Line 439 | void scan_SUSY_parameter_space(string mc
439		string prefix="SMS_";
440		// up to here, everything is set up for SMS; now we need to switch stuff around if we're dealing with an mSUGRA scan!
441		int scantype=SMS;
442	<	if(Contains((scansample.collection)[0].samplename,"mSUGRA")) scantype=mSUGRA;
443	<	if(Contains((scansample.collection)[0].samplename,"GMSB")) scantype=GMSB;
442	>	if(Contains((PlottingSetup::ScanSampleDirectory),"SUGRA")) scantype=mSUGRA;
443	>	if(Contains((PlottingSetup::ScanSampleDirectory),"GMSB")) scantype=GMSB;
444
445		if(scantype==mSUGRA) {
446		massgluname="M0"; // this is the "x axis" in the limit plot (like the gluino in the SMS case)
#	Line 441 | Line 500 | mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLS
500		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);
501		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);
502		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);
503	+	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);
504	+	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);
505
506		TH2F *imposedxmap;
507		TH2F *realxmap;
#	Line 515 | Line 576 | write_warning(__FUNCTION__,"CURRENTLY SW
576		int a = int((PlottingSetup::ScanXzones*mlsp)/(m12trees+1));
577		int b = int((PlottingSetup::ScanYzones*mglu)/(m0trees+1));
578		int scanfileindex=PlottingSetup::ScanYzones*a+b;
579	<	load_scan_sample(a,b,scanfileindex,scantype);
579	>	load_scan_sample(a,b,scanfileindex,scantype,mglu,mlsp);
580
581		clock_t start,finish;
582		start = clock(); // starting the clock to measure how long the computation takes!
#	Line 578 | Line 639 | write_warning(__FUNCTION__,"CURRENTLY SW
639		}
640
641
642	+	if(scantype==mSUGRA) {
643	+	float absxs=0;
644	+	for(int i=0;i<12;i++) absxs+=GetXSecForPointAndChannel(mglu,mlsp,xsec,i);
645	+	TFile *FilterEffFile = new TFile(PlottingSetup::FilterEfficiencyFile.c_str());
646	+	TH2F *FilterEfficiency = (TH2F*)FilterEffFile->Get("FilterEfficiency");
647	+	float filtereff = FilterEfficiency->GetBinContent(FilterEfficiency->FindBin(mglu,mlsp));
648	+	FilterEff->SetBinContent(GlobalBin,filtereff);
649	+	FilterEffFile->Close();
650	+	absXSmap->SetBinContent(GlobalBin,absxs);
651	+	absxs*=filtereff;
652	+	XSmap->SetBinContent(GlobalBin,absxs);
653	+	}
654	+
655		if(nevents!=0&&(efficiencyonly||systematicsonly)) {
656		Value effwosigcont = MCefficiency((scansample.collection)[scanfileindex].events,result,resulterr,flipped,mcjzb,requireZ,(int)nevents,scantype,xsec,addcut.str(),-1);
657		if(result<0&&allowflipping) {
#	Line 615 | Line 689 | write_warning(__FUNCTION__,"CURRENTLY SW
689		Neventsmap->SetBinContent(GlobalBin,nevents);
690		ipointmap->SetBinContent(GlobalBin,ipoint);
691		if(efficiencyonly) continue;
692	<
693	<	do_systematics_for_one_file((scansample.collection)[scanfileindex].events,(int)nevents,"SUSY SCAN", systematics,flipped, xsec, mcjzb,datajzb,peakerror,requireZ, addcut.str(),true,scantype);
694	<	float JZBcutat = systematics[0][0];
695	<	float mceff    = systematics[0][1];
696	<	float mcefferr = systematics[0][2];//MC stat error
697	<	float toterr   = systematics[0][4];
698	<	float sys_jes  = systematics[0][5]; // Jet Energy Scale
699	<	float sys_jsu  = systematics[0][6]; // JZB scale uncertainty
700	<	float sys_res  = systematics[0][7]; // resolution
701	<	float mcwoscef = systematics[0][8]; // efficiency without signal contamination
702	<	float mcwoscefr= systematics[0][9]; // error on efficiency without signal contamination
703	<	float sys_pdf   = 0;
704	<	if(systematics[0].size()>10) sys_pdf = systematics[0][10]; // PDF
705	<
706	<	if(mceff!=mceff||toterr!=toterr||mceff<0 && (!systematicsonly&&!efficiencyonly)) {
707	<	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;
708	<	continue;
692	>
693	>	if(!shapeanalysis) {
694	>	do_systematics_for_one_file((scansample.collection)[scanfileindex].events,(int)nevents,"SUSY SCAN", systematics,flipped, xsec, mcjzb,datajzb,peakerror,requireZ, addcut.str(),true,scantype);
695	>	float JZBcutat = systematics[0][0];
696	>	float mceff    = systematics[0][1];
697	>	float mcefferr = systematics[0][2];//MC stat error
698	>	float toterr   = systematics[0][4];
699	>	float sys_jes  = systematics[0][5]; // Jet Energy Scale
700	>	float sys_jsu  = systematics[0][6]; // JZB scale uncertainty
701	>	float sys_res  = systematics[0][7]; // resolution
702	>	float mcwoscef = systematics[0][8]; // efficiency without signal contamination
703	>	float mcwoscefr= systematics[0][9]; // error on efficiency without signal contamination
704	>	float sys_pdf   = 0;
705	>	if(systematics[0].size()>10) sys_pdf = systematics[0][10]; // PDF
706	>	cout << "Going to store: total error is: " << toterr << " and the relative tot error, i.e. total error / mc efficiency is " << toterr/mceff << endl;
707	>	if(mceff!=mceff||toterr!=toterr||mceff<0 && (!systematicsonly&&!efficiencyonly)) {
708	>	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;
709	>	continue;
710	>	} else {
711	>	// Systematics and efficiencies make sense (i.e. non-zero, and all numbers)
712	>	if(!systematicsonly&&!efficiencyonly) {
713	>	dout << "Calculating limit now for "<<massgluname<<"="<<mglu<<" , "<<massLSPname<<"="<<mlsp <<endl;
714	>	vector<float> sigmas;
715	>	string plotfilename=(string)(TString((scansample.collection)[scanfileindex].samplename)+TString(massgluname)+TString(any2string(mglu))+TString("__")+TString(massLSPname)+TString(any2string(mlsp))+TString(".png"));
716	>	sigmas=compute_one_upper_limit(mceff,toterr,ibin,mcjzb,plotfilename,true,flipped);
717	>	//          do_limit_wrapper(mceff,toterr,ibin,mcjzb,sigmas,plotfilename);
718	>	if(sigmas[0]>-0.5) { // negative sigmas are the error signature of do_limit_wrapper, so we want to exclude them.
719	>	limitmap->SetBinContent(GlobalBin,sigmas[0]);
720	>	if(sigmas.size()>1) {
721	>	explimitmap->SetBinContent(GlobalBin,sigmas[1]);
722	>	exp1plimitmap->SetBinContent(GlobalBin,sigmas[2]);
723	>	exp1mlimitmap->SetBinContent(GlobalBin,sigmas[3]);
724	>	exp2plimitmap->SetBinContent(GlobalBin,sigmas[4]);
725	>	exp2mlimitmap->SetBinContent(GlobalBin,sigmas[5]);
726	>	}
727	>	sysjesmap->SetBinContent(GlobalBin,sys_jes);
728	>	sysjsumap->SetBinContent(GlobalBin,sys_jsu);
729	>	sysresmap->SetBinContent(GlobalBin,sys_res);
730	>	syspdfmap->SetBinContent(GlobalBin,sys_pdf);
731	>	systotmap->SetBinContent(GlobalBin,toterr/mceff);//total relative (!) error
732	>	sysstatmap->SetBinContent(GlobalBin,mcefferr);//total relative (!) error
733	>	dout << "A limit has been added at " << sigmas[0] << " for m_{glu}="<<mglu << " and m_{lsp}="<<mlsp<<endl;
734	>	}//end of if sigma is positive
735	>	}//end of not systematics only condition
736	>	if(systematicsonly) {
737	>	sysjesmap->SetBinContent(GlobalBin,sys_jes);
738	>	sysjsumap->SetBinContent(GlobalBin,sys_jsu);
739	>	sysresmap->SetBinContent(GlobalBin,sys_res);
740	>	syspdfmap->SetBinContent(GlobalBin,sys_pdf);
741	>	systotmap->SetBinContent(GlobalBin,toterr/mceff);//total relative (!) error
742	>	sysstatmap->SetBinContent(GlobalBin,mcefferr);//total relative (!) error
743	>	}
744	>	}//efficiency is valid
745		} else {
746	<	if(!systematicsonly&&!efficiencyonly) {
747	<	dout << "Calculating limit now for "<<massgluname<<"="<<mglu<<" , "<<massLSPname<<"="<<mlsp <<endl;
748	<	vector<float> sigmas;
749	<	string plotfilename=(string)(TString((scansample.collection)[scanfileindex].samplename)+TString(massgluname)+TString(any2string(mglu))+TString("__")+TString(massLSPname)+TString(any2string(mlsp))+TString(".png"));
750	<	sigmas=compute_one_upper_limit(mceff,toterr,ibin,mcjzb,plotfilename,true,flipped);
751	<	//      do_limit_wrapper(mceff,toterr,ibin,mcjzb,sigmas,plotfilename);
752	<	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]);
746	>	//shape analysis
747	>	cout << "This is the shape analysis. Have fun with it!" << endl;
748	>	float quickmceff,quickmcefferr;
749	>	MCefficiency((scansample.collection)[scanfileindex].events,quickmceff,quickmcefferr,flipped,mcjzb,requireZ,(int)nevents,scantype,xsec,addcut.str(),-1);
750	>	if(quickmceff==0) {
751	>	write_warning(__FUNCTION__,"The MC efficiency is zero - need to skip this point.");
752	>	continue;
753		}
754	<
755	<	sysjesmap->SetBinContent(GlobalBin,sys_jes);
756	<	sysjsumap->SetBinContent(GlobalBin,sys_jsu);
757	<	sysresmap->SetBinContent(GlobalBin,sys_res);
758	<	syspdfmap->SetBinContent(GlobalBin,sys_pdf);
759	<	systotmap->SetBinContent(GlobalBin,toterr/mceff);//total relative (!) error
760	<	sysstatmap->SetBinContent(GlobalBin,mcefferr);//total relative (!) error
761	<	dout << "A limit has been added at " << sigmas[0] << " for m_{glu}="<<mglu << " and m_{lsp}="<<mlsp<<endl;
762	<	} //end of if sigma is positive
763	<
764	<	//end of not systematics only condition
765	<	}
766	<	if(systematicsonly) {
767	<	sysjesmap->SetBinContent(GlobalBin,sys_jes);
768	<	sysjsumap->SetBinContent(GlobalBin,sys_jsu);
769	<	sysresmap->SetBinContent(GlobalBin,sys_res);
770	<	syspdfmap->SetBinContent(GlobalBin,sys_pdf);
771	<	systotmap->SetBinContent(GlobalBin,toterr/mceff);//total relative (!) error
772	<	sysstatmap->SetBinContent(GlobalBin,mcefferr);//total relative (!) error
773	<	}
672	<	}//efficiency is valid
754	>	stringstream PointName;
755	>	PointName << massgluname << "_" << mglu << "__" << massLSPname << "_" << mlsp;
756	>	SUSYScanSpace::SavedMLSP=mlsp;
757	>	SUSYScanSpace::SavedMGlu=mglu;
758	>	SUSYScanSpace::SavedMLSPname=massLSPname;
759	>	SUSYScanSpace::SavedMGluname=massgluname;
760	>	ShapeDroplet droplet = LimitsFromShapes((scansample.collection)[scanfileindex].events,addcut.str().c_str(),PointName.str(),mcjzb,datajzb,jzb_cut,peakerror,peakerror);
761	>	limitmap->SetBinContent(GlobalBin,droplet.observed);
762	>	explimitmap->SetBinContent(GlobalBin,droplet.expected);
763	>	exp1plimitmap->SetBinContent(GlobalBin,droplet.expectedPlus1Sigma);
764	>	exp1mlimitmap->SetBinContent(GlobalBin,droplet.expectedMinus1Sigma);
765	>	exp2plimitmap->SetBinContent(GlobalBin,droplet.expectedPlus2Sigma);
766	>	exp2mlimitmap->SetBinContent(GlobalBin,droplet.expectedMinus2Sigma);
767	>
768	>	sysjesmap->SetBinContent(GlobalBin,droplet.JES);
769	>	sysjsumap->SetBinContent(GlobalBin,droplet.JSU);
770	>	syspdfmap->SetBinContent(GlobalBin,droplet.PDF);
771	>	systotmap->SetBinContent(GlobalBin,droplet.toterr);
772	>	sysstatmap->SetBinContent(GlobalBin,droplet.staterr);
773	>	}
774		finish = clock();
775		timemap->SetBinContent(GlobalBin,((float(finish)-float(start))/CLOCKS_PER_SEC));
776		}
#	Line 724 | Line 825 | write_warning(__FUNCTION__,"CURRENTLY SW
825		systotmap->Write();
826		sysstatmap->Write();
827		XSmap->Write();
828	+	absXSmap->Write();
829	+	FilterEff->Write();
830		if(ibin==0) imposedxmap->Write();
831		if(ibin==0) realxmap->Write();
832		Neventsmap->Write();
#	Line 786 | Line 889 | write_warning(__FUNCTION__,"CURRENTLY SW
889		syspdfmap->Write();
890		systotmap->Write();
891		sysstatmap->Write();
892	+	absXSmap->Write();
893		XSmap->Write();
894	+	FilterEff->Write();
895		if(ibin==0) imposedxmap->Write();
896		if(ibin==0) realxmap->Write();
897		timemap->Write();
#	Line 900 | Line 1005 | write_warning(__FUNCTION__,"CURRENTLY SW
1005		delete systotmap;
1006		delete sysstatmap;
1007		delete XSmap;
1008	+	delete absXSmap;
1009		delete limcanvas;
1010		}
1011
1012	<	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) {
1012	>	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) {
1013		dout << "Starting the SUSY scan now with all " << jzb_cut.size() << " bin(s)" << endl;
1014		for(int ibin=0;ibin<jzb_cut.size();ibin++) {
1015	<	scan_SUSY_parameter_space(mcjzb,datajzb,jzb_cut,requireZ, peakerror, ibin, njobs, jobnumber,systonly,effonly);
1015	>	scan_SUSY_parameter_space(mcjzb,datajzb,jzb_cut,requireZ, peakerror, peakerrordata, ibin, njobs, jobnumber,systonly,effonly,shapeanalysis);
1016		}
1017	+	delete_any_cached_scans();// tidy up ...
1018		}
1019

Diff Legend

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