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.2 by buchmann, Fri Mar 2 08:40:26 2012 UTC vs.
Revision 1.7 by buchmann, Fri Mar 23 12:04:29 2012 UTC

#	Line 24 | Line 24 | using namespace std;
24
25		using namespace PlottingSetup;
26
27	+	namespace SUSYScanSpace {
28	+	vector<string> loaded_files;
29	+	int SUSYscantype;
30	+	}
31	+
32		void susy_scan_axis_labeling(TH2F *histo) {
33		histo->GetXaxis()->SetTitle("m_{#Chi_{2}^{0}}-m_{LSP}");
34		histo->GetXaxis()->CenterTitle();
#	Line 67 | Line 72 | void set_SUSY_style() {
72		}
73
74		bool delete_any_cached_scans() {
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);
83	<	if((Contains(hostname,"t3ui")||Contains(hostname,"t3wn"))) return false;
84	<	else {
83	>	if(!((Contains(hostname,"t3ui")||Contains(hostname,"t3wn"))))
84	>	{
85		vector<string> all_files;
86		char currentpath[1024];
87		TString directory=TString(getcwd(currentpath,1024));
#	Line 86 | Line 94 | bool delete_any_cached_scans() {
94		gSystem->Exec(((string)"rm "+all_files[ifile]).c_str());
95		}
96		}
97	<	return true;
97	>	return true;
98	>	}
99	>	*/
100	>	if(SUSYScanSpace::SUSYscantype==mSUGRA) {
101	>	cout << "Going to purge files in local_storage that have been loaded!" << endl;
102	>	for(int i=0;i<SUSYScanSpace::loaded_files.size();i++) {
103	>	gSystem->Exec(((string)"rm "+SUSYScanSpace::loaded_files[i]).c_str());
104	>	dout << "      Purging : Deleted file " << SUSYScanSpace::loaded_files[i] << endl;
105	>	}
106		}
107	+	return true;
108		}
109
110		bool initialized_t2=false;
111
112		int srmcpretries=0;
113
114	<	void load_scan_sample(int a, int b, int &scanfileindex,int scantype,bool isretry=false) {
114	>	void load_scan_sample(int a, int b, int &scanfileindex,int scantype,int mglu, int mlsp, bool isretry=false) {
115
116		// There is no need to define your sample here. That is done in Setup.C where you define the loading directory!
117
118	<	dout << "Going to load file with Xzone=" << a << " and  Yzone=" << b << endl;
118	>	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;
119		stringstream filetoload;
120		char hostname[1023];
121		gethostname(hostname,1023);
122	<
123	<	if((Contains(hostname,"t3ui")||Contains(hostname,"t3wn"))) { // local case
122	>	string samplename;
123	>	if((Contains(hostname,"t3ui")||Contains(hostname,"t3wn"))&&scantype!=mSUGRA) { // local case
124		filetoload << "/shome/buchmann/ntuples/"<<PlottingSetup::ScanSampleDirectory;
125	<	if(scantype==mSUGRA) filetoload << "/mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
126	<	if(scantype==SMS) filetoload << "/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
127	<	if(scantype==GMSB) filetoload << "/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
128	<	if(!Contains(((scansample.collection)[(scansample.collection).size()-1]).filename,"_"+any2string(a)+"_"+any2string(b))) {
125	>	if(scantype==mSUGRA) {
126	>	//filetoload << "/mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
127	>	filetoload.str("");
128	>	filetoload << "/shome/lbaeni/jzb/" << PlottingSetup::ScanSampleDirectory  << "/mSUGRA_M0_" << mglu << "__M12_" << mlsp << ".root";
129	>	samplename="mSUGRA";
130	>	//      filetoload << "/mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
131	>	}
132	>	if(scantype==SMS) {
133	>	filetoload << "/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
134	>	samplename="SMS";
135	>	}
136	>	if(scantype==GMSB) {
137	>	filetoload << "/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
138	>	samplename="GMSB";
139	>	}
140	>	if(scansample.collection.size()<1||!Contains(((scansample.collection)[(scansample.collection).size()-1]).filename,"_"+any2string(a)+"_"+any2string(b))) {
141		dout << "The last sample is NOT the same one as the current one, possibly popping off last one and adding the new one." << endl;
142	<	if((scansample.collection).size()>1) {
142	>	if((scansample.collection).size()>=1) {
143		scansample.RemoveLastSample();
144		}
145		scanfileindex=(scansample.collection).size();
146		//New: Loading file when necessary, not before (avoiding high memory usage and startup times)
147	<	if(scanfileindex!=0) scansample.AddSample(filetoload.str(),"scansample",1,1,false,true,scanfileindex,kRed);
147	>
148	>	scansample.AddSample(filetoload.str(),samplename,1,1,false,true,scanfileindex,kRed);
149	>	dout << "Just added the following file: " << filetoload.str() << endl;
150		} else {
151		dout << "Last sample is the same as the current one. Recycling it." << endl;
152		scanfileindex=(scansample.collection).size()-1;
#	Line 123 | Line 154 | void load_scan_sample(int a, int b, int
154		dout << " Going to use the following file: " << filetoload.str() << endl;
155		} // end of t3 case
156		else {
157	+	//CRAB case
158		write_info(__FUNCTION__,"Hello, CRAB! Might need to load files ...");
159		PlottingSetup::limitpatience=PlottingSetup::limitpatienceCRAB;
160		stringstream copyfile;
161	+	gSystem->Exec("mkdir -p local_storage");
162		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;
163	<	if(scantype==mSUGRA) copyfile << "/mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root " << "mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root ";
164	<	if(scantype==SMS) copyfile << "/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root " << "SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root ";
165	<	if(scantype==GMSB) copyfile << "/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root " << "GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root ";
163	>	if(scantype==mSUGRA) {
164	>	copyfile<< "/mSUGRA_M0_" << mglu << "__M12_" << mlsp << ".root local_storage/mSUGRA_M0_" << mglu << "__M12_" << mlsp << ".root ";
165	>	samplename="mSUGRA";
166	>	}
167	>	if(scantype==SMS) {
168	>	copyfile << "/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root " << "local_storage/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root ";
169	>	samplename="SMS";
170	>	}
171	>	if(scantype==GMSB) {
172	>	copyfile << "/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root " << "local_storage/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root ";
173	>	samplename="GMSB";
174	>	}
175		stringstream newfilename;
176	<	if(scantype==mSUGRA) newfilename << "mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
177	<	if(scantype==SMS) newfilename << "SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
178	<	if(scantype==GMSB) newfilename << "GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
176	>	//    if(scantype==mSUGRA) newfilename << "local_storage/mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
177	>	if(scantype==mSUGRA) newfilename << "local_storage/mSUGRA_M0_" << mglu << "__M12_" << mlsp << ".root";
178	>	if(scantype==SMS) newfilename << "local_storage/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
179	>	if(scantype==GMSB) newfilename << "local_storage/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
180
181	<	if(!Contains(((scansample.collection)[(scansample.collection).size()-1]).filename,"_"+any2string(a)+"_"+any2string(b))||(a==0&&b==0)&&!initialized_t2) {
181	>	SUSYScanSpace::SUSYscantype=scantype;
182	>
183	>	if((scansample.collection).size()==0||!Contains(((scansample.collection)[(scansample.collection).size()-1]).filename,"_"+any2string(a)+"_"+any2string(b))||(a==0&&b==0)&&!initialized_t2) {
184		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;
185	<	if((scansample.collection).size()>1) {
185	>	while((scansample.collection).size()>0) {
186		scansample.RemoveLastSample();
187		}
188	+	cout << "New scanfileindex stems from scansample size: " << (scansample.collection).size() << endl;
189		scanfileindex=(scansample.collection).size();
190		//New: Loading file when necessary, not before (avoiding high memory usage and startup times)
191	<	if(scanfileindex!=0||isretry) {
192	<	delete_any_cached_scans();
193	<	dout << "Going to download the scan file with the following copy command: " << copyfile.str() << endl;
194	<	gSystem->Exec(copyfile.str().c_str());
191	>	if(1) {
192	>	if(!doesROOTFileExist(newfilename.str())) {
193	>	dout << "Going to download the scan file with the following copy command: " << copyfile.str() << endl;
194	>	int retcode = gSystem->Exec(copyfile.str().c_str()); // download it if it hasn't been downloaded before
195	>	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
196	>	}
197		if(doesROOTFileExist(newfilename.str())) {
198		initialized_t2=true;
199	<	scansample.AddSample(newfilename.str(),"scansample",1,1,false,true,scanfileindex,kRed);
199	>	cout << "Going to add the following sample: " << newfilename.str() << endl;
200	>	gSystem->Exec(("ls -ltrh "+newfilename.str()).c_str());
201	>	scansample.AddSample(newfilename.str(),samplename,1,1,false,true,scanfileindex,kRed);
202		} else {
203		srmcpretries++;
204		if(srmcpretries<5) {
205		dout << "The file could not be loaded correctly - retrying!" << endl;
206		sleep(5);
207	<	load_scan_sample(a,b,scanfileindex,scantype,true);
207	>	load_scan_sample(a,b,scanfileindex,scantype,mglu,mlsp,true);
208		} else {
209		dout << "Have tried 5 times to load this sample. Giving up now and failing the program execution" << endl;
210	<	assert(0);
210	>	assert(srmcpretries<5);
211		dout << "The command " << copyfile.str() << " failed to copy the file to the local storage (saving it as " << newfilename.str() << ")" << endl;
212		}
213		}
#	Line 168 | Line 218 | void load_scan_sample(int a, int b, int
218		}
219
220		}
221	+
222	+
223	+	// WATCH OUT THIS LINE NEEDS TO BE REMOVED
224	+	scanfileindex=0;
225		}
226
227	<	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)  {
227	>	/*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)  {
228		altxs=0;
229	<	bool HaveCorrectlyImplementedkFactorsFormSUGRA=false; // still need to work on k factors for mSUGRA!
176	<	assert(HaveCorrectlyImplementedkFactorsFormSUGRA);
177	<	for(int iproc=1;iproc<=10;iproc++) {
178	<	float process_xs = GetXSecForPointAndChannel(mglu,mlsp,xsec,iproc);
179	<	altxs+=process_xs;
180	<	}
229	>	for(int iproc=1;iproc<=10;iproc++) altxs+=GetXSecForPointAndChannel(mglu,mlsp,xsec,iproc);
230		return altxs;
231	<	}
231	>	}*/
232
233		void establish_SUSY_limits(string mcjzb,string datajzb,vector<float> jzb_cut,bool requireZ, float peakerror, TFile *fsyst, int ibin,float njobs=-1, float jobnumber=-1) {
234		bool runninglocally=true;
#	Line 206 | Line 255 | void establish_SUSY_limits(string mcjzb,
255		while ((key = (TKey*)nextkey()))
256		{
257		TObject *obj = key->ReadObj();
258	<	if(Contains((string)(obj->GetName()),"mSUGRA")) scantype=mSUGRA;
258	>	if(Contains((string)(obj->GetName()),"SUGRA")) scantype=mSUGRA;
259		if(Contains((string)(obj->GetName()),"GMSB")) scantype=GMSB;
260		}
261
#	Line 251 | Line 300 | void establish_SUSY_limits(string mcjzb,
300		TH2F *fullerr   = (TH2F*)fsyst->Get((prefix+"systotmap"+any2string(jzb_cut[ibin])).c_str());
301		TH2F *NEvents   = (TH2F*)fsyst->Get((prefix+"Neventsmap"+any2string(jzb_cut[ibin])).c_str());
302		TH2F *lflip     = (TH2F*)fsyst->Get((prefix+"flipmap"+any2string(jzb_cut[ibin])).c_str());
303	+	TH2F *XSmap     = (TH2F*)fsyst->Get((prefix+"XS"+any2string(jzb_cut[ibin])).c_str());
304
305		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
306		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
#	Line 260 | Line 310 | void establish_SUSY_limits(string mcjzb,
310		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
311
312		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);
313		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);
314	+	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);
315
316
317		if(!fullerr || !mceff || !NEvents) {
#	Line 294 | Line 343 | void establish_SUSY_limits(string mcjzb,
343		continue;
344		}
345		vector<float> sigmas;
346	<	sigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped);
347	<
346	>	if(scantype!=mSUGRA) sigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped);
347	>	else {
348	>	//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"!
349	>	vector<float> asigmas;
350	>	asigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped,true); // asymptotic limit first
351	>	float strength=asigmas[0]/XSmap->GetBinContent(GlobalBin);
352	>	if(strength>0.5&&strength<2) {
353	>	sigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped); // asymptotic limit first
354	>	asymptoticmap->SetBinContent(GlobalBin,0);
355	>	}
356	>	else {
357	>	sigmas=asigmas;
358	>	asymptoticmap->SetBinContent(GlobalBin,1);
359	>	}
360	>	exclmap->SetBinContent(GlobalBin,strength);
361	>	}
362
363		if(sigmas[0]>0) limitmap->SetBinContent(GlobalBin,sigmas[0]); //anything else is an error code
364		if(sigmas.size()>1) {
#	Line 307 | Line 370 | void establish_SUSY_limits(string mcjzb,
370		}
371
372		dout << "An upper limit has been added for this point ( " << massgluname << " = " << mglu << " and " << massLSPname << " = " << mlsp << " ) at " << sigmas[0] << endl;
310	–	if(scantype==mSUGRA) { // for SMS this is a bit easier at the moment - we have a reference XS file which we use when plotting
311	–	dout << "Computing exclusion status" << endl;
312	–	float rel_limit=0;
313	–	float totxs;
314	–	float xs=get_xs(totxs,mglu,mlsp,massgluname,massLSPname,xsec,mcjzb,requireZ);
315	–	if(xs>0) rel_limit=sigmas[0]/xs;
316	–	exclmap->SetBinContent(GlobalBin,rel_limit);
317	–	xsmap->SetBinContent(GlobalBin,xs);
318	–	totxsmap->SetBinContent(GlobalBin,totxs);
319	–	}
373		}
374		}
375
376		prepare_scan_axis(limitmap,scantype);
377	+	gSystem->Exec("mkdir -p output");
378		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 :-)
379		outputfile->cd();
380		limitmap->Write();
381		flipmap->Write();
382	+	asymptoticmap->Write();
383		if(doexpected) {
384		explimitmap->Write();
385		exp1plimitmap->Write();
#	Line 332 | Line 387 | void establish_SUSY_limits(string mcjzb,
387		exp2plimitmap->Write();
388		exp2mlimitmap->Write();
389		}
390	<	if(scantype==mSUGRA) {
336	<	exclmap->Write();
337	<	xsmap->Write();
338	<	totxsmap->Write();
339	<	}
390	>	if(scantype==mSUGRA) exclmap->Write();
391		outputfile->Close();
392		delete limcanvas;
393		}
#	Line 378 | Line 429 | void scan_SUSY_parameter_space(string mc
429		string prefix="SMS_";
430		// up to here, everything is set up for SMS; now we need to switch stuff around if we're dealing with an mSUGRA scan!
431		int scantype=SMS;
432	<	if(Contains((scansample.collection)[0].samplename,"mSUGRA")) scantype=mSUGRA;
432	>	if(Contains((scansample.collection)[0].samplename,"SUGRA")) scantype=mSUGRA;
433		if(Contains((scansample.collection)[0].samplename,"GMSB")) scantype=GMSB;
434
435		if(scantype==mSUGRA) {
#	Line 430 | Line 481 | void scan_SUSY_parameter_space(string mc
481		TH2F *efficiencymapID   =  new TH2F((prefix+"efficiencymapID"+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);
482		TH2F *efficiencymapJets =  new TH2F((prefix+"efficiencymapJets"+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);
483		TH2F *efficiencymapSignal = new TH2F((prefix+"efficiencymapSignal"+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);
484	<	TH2F *efficiencymapContam = new TH2F((prefix+"efficiencymapContam"+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);
484	>	TH2F *efficiencymapContam = new TH2F((prefix+"efficiencymapContam"+any2string(jzbSel)).c_str(),"",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-
485	>	mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
486		TH2F *noscefficiencymap = new TH2F((prefix+"noscefficiencymap"+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);
487		TH2F *Neventsmap     = new TH2F((prefix+"Neventsmap"+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);
488		TH2F *ipointmap      = new TH2F((prefix+"ipointmap"+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);
489		TH2F *syspdfmap      = new TH2F((prefix+"syspdfmap"+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);
490		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);
491		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);
492	+	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);
493
494		TH2F *imposedxmap;
495		TH2F *realxmap;
496
497	+	TH2F *centermap;
498	+	TH2F *widthmap;
499	+
500		TH2F *timemap        = new TH2F((prefix+"timemap"+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
502		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);
#	Line 448 | Line 504 | void scan_SUSY_parameter_space(string mc
504		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);
505
506		if(ibin==0) {
507	+	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);
508	+	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);
509	+
510		imposedxmap    = new TH2F((prefix+"imposedxmap").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);
511		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);
512		}
#	Line 485 | Line 544 | write_warning(__FUNCTION__,"CURRENTLY SW
544		for(int mglu=(int)mglustart;mglu<=mgluend;mglu+=(int)mglustep) {
545		for (int mlsp=(int)mLSPstart;mlsp<=mLSPend&&(scantype==mSUGRA||mlsp<=mglu);mlsp+=(int)mLSPstep) Npoints++;
546		}
547	<
547	>
548	>	map <  pair<float, float>, map<string, float>  >  xsec;
549	>	if(scantype==mSUGRA) {
550	>	xsec=getXsec(PlottingSetup::mSUGRAxsFile);
551	>	}
552	>
553		int ipoint=-1;
554		for(int mglu=(int)mglustart;mglu<=mgluend;mglu+=(int)mglustep) {
555		for (int mlsp=(int)mLSPstart;mlsp<=mLSPend&&(scantype==mSUGRA||mlsp<=mglu);mlsp+=(int)mLSPstep)
#	Line 500 | Line 564 | write_warning(__FUNCTION__,"CURRENTLY SW
564		int a = int((PlottingSetup::ScanXzones*mlsp)/(m12trees+1));
565		int b = int((PlottingSetup::ScanYzones*mglu)/(m0trees+1));
566		int scanfileindex=PlottingSetup::ScanYzones*a+b;
567	<	load_scan_sample(a,b,scanfileindex,scantype);
567	>	load_scan_sample(a,b,scanfileindex,scantype,mglu,mlsp);
568
569		clock_t start,finish;
570		start = clock(); // starting the clock to measure how long the computation takes!
#	Line 510 | Line 574 | write_warning(__FUNCTION__,"CURRENTLY SW
574		float nevents = (scansample.collection)[scanfileindex].events->GetSelectedRows();
575		vector<vector<float> > systematics;
576		if(nevents<10) {
577	<	dout << "This point ("<<ipoint<<") with configuration ("<<massgluname<<"="<<mglu<<" , "<<massLSPname<<"="<<mlsp << ") does not contain enough events and will be \033[1;31m skipped \033[0m."<< endl;
577	>	dout << "This point ("<<ipoint<<" / " << Npoints << ") with configuration ("<<massgluname<<"="<<mglu<<" , "<<massLSPname<<"="<<mlsp << ") does not contain enough events and will be \033[1;31m skipped \033[0m."<< endl;
578		continue;
579		} else {
580		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;
#	Line 552 | Line 616 | write_warning(__FUNCTION__,"CURRENTLY SW
616		(scansample.collection)[scanfileindex].events->Draw("imposedx>>realxhisto",addcut.str().c_str(),"goff");
617		imposedxmap->SetBinContent(GlobalBin,realxhisto->GetMean());
618		imposedxmap->SetBinError(GlobalBin,realxhisto->GetRMS());
619	+	stringstream specialexpression;
620	+	specialexpression << mcjzb << ">>jzbhisto";
621	+	TH1F *jzbhisto = new TH1F("jzbhisto","jzbhisto",100,-500,500);
622	+	(scansample.collection)[scanfileindex].events->Draw(specialexpression.str().c_str(),addcut.str().c_str(),"goff");
623	+	centermap->SetBinContent(GlobalBin,jzbhisto->GetMean());
624	+	widthmap->SetBinContent(GlobalBin,jzbhisto->GetRMS());
625		delete realxhisto;
626	+	delete jzbhisto;
627		}
628
629
630	+	if(scantype==mSUGRA) {
631	+	float absxs=0;
632	+	for(int i=0;i<12;i++) absxs+=GetXSecForPointAndChannel(mglu,mlsp,xsec,i);
633	+	XSmap->SetBinContent(GlobalBin,absxs);
634	+	}
635	+
636		if(nevents!=0&&(efficiencyonly||systematicsonly)) {
637	<	Value effwosigcont = MCefficiency((scansample.collection)[scanfileindex].events,result,resulterr,flipped,mcjzb,requireZ,(int)nevents,addcut.str(),-1);
637	>	Value effwosigcont = MCefficiency((scansample.collection)[scanfileindex].events,result,resulterr,flipped,mcjzb,requireZ,(int)nevents,scantype,xsec,addcut.str(),-1);
638		if(result<0&&allowflipping) {
639		write_info(__FUNCTION__,"Found negative efficiency (signal distribution populating JZB<0 more than JZB>0) - flipping analysis!");
640		flipped=1;
641	<	MCefficiency((scansample.collection)[scanfileindex].events,result,resulterr,flipped,mcjzb,requireZ,(int)nevents,addcut.str(),-1);
641	>	MCefficiency((scansample.collection)[scanfileindex].events,result,resulterr,flipped,mcjzb,requireZ,(int)nevents,scantype,xsec,addcut.str(),-1);
642		}
643		efficiencymap->SetBinContent(GlobalBin,result);
644		efficiencymap->SetBinError(GlobalBin,resulterr);
#	Line 594 | Line 671 | write_warning(__FUNCTION__,"CURRENTLY SW
671		ipointmap->SetBinContent(GlobalBin,ipoint);
672		if(efficiencyonly) continue;
673
674	<	do_systematics_for_one_file((scansample.collection)[scanfileindex].events,(int)nevents,"SUSY SCAN", systematics,flipped, mcjzb,datajzb,peakerror,requireZ, addcut.str(),true);//mSUGRA is now always true here because we always want to compute PDF systematics
674	>	do_systematics_for_one_file((scansample.collection)[scanfileindex].events,(int)nevents,"SUSY SCAN", systematics,flipped, xsec, mcjzb,datajzb,peakerror,requireZ, addcut.str(),true,scantype);
675		float JZBcutat = systematics[0][0];
676		float mceff    = systematics[0][1];
677		float mcefferr = systematics[0][2];//MC stat error
#	Line 606 | Line 683 | write_warning(__FUNCTION__,"CURRENTLY SW
683		float mcwoscefr= systematics[0][9]; // error on efficiency without signal contamination
684		float sys_pdf   = 0;
685		if(systematics[0].size()>10) sys_pdf = systematics[0][10]; // PDF
686	<
686	>	cout << "Going to store: total error is: " << toterr << " and the relative tot error, i.e. total error / mc efficiency is " << toterr/mceff << endl;
687		if(mceff!=mceff||toterr!=toterr||mceff<0 && (!systematicsonly&&!efficiencyonly)) {
688		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;
689		continue;
#	Line 688 | Line 765 | write_warning(__FUNCTION__,"CURRENTLY SW
765		sysjsumap->Write();
766		sysresmap->Write();
767		efficiencymap->Write();
768	+	if(ibin==0) centermap->Write();
769	+	if(ibin==0) widthmap->Write();
770		flipmap->Write();
771		efficiencymapmet->Write();
772		efficiencymapAcc->Write();
#	Line 699 | Line 778 | write_warning(__FUNCTION__,"CURRENTLY SW
778		syspdfmap->Write();
779		systotmap->Write();
780		sysstatmap->Write();
781	+	XSmap->Write();
782		if(ibin==0) imposedxmap->Write();
783		if(ibin==0) realxmap->Write();
784		Neventsmap->Write();
785		ipointmap->Write();
786		timemap->Write();
707	–	outputfile->Close();
787		for(int i=0;i<susy_Zdecay_origin.size();i++) {
788		h_susy_Zdecay_origin[i]->Write();
789		}
790	+	outputfile->Close();
791		}
792		}
793		if(systematicsonly) { // systematics only :
#	Line 746 | Line 826 | write_warning(__FUNCTION__,"CURRENTLY SW
826		sysjsumap->Write();
827		sysresmap->Write();
828		flipmap->Write();
829	+	if(ibin==0) centermap->Write();
830	+	if(ibin==0) widthmap->Write();
831		efficiencymap->Write();
832		efficiencymapmet->Write();
833		efficiencymapAcc->Write();
#	Line 759 | Line 841 | write_warning(__FUNCTION__,"CURRENTLY SW
841		syspdfmap->Write();
842		systotmap->Write();
843		sysstatmap->Write();
844	+	XSmap->Write();
845		if(ibin==0) imposedxmap->Write();
846		if(ibin==0) realxmap->Write();
847		timemap->Write();
765	–	outputfile->Close();
848		for(int i=0;i<susy_Zdecay_origin.size();i++) {
849		h_susy_Zdecay_origin[i]->Write();
850		}
851	+	outputfile->Close();
852		}
853		}//end of systematics only
854		if(efficiencyonly) {
#	Line 832 | Line 915 | write_warning(__FUNCTION__,"CURRENTLY SW
915		noscefficiencymap->Write();
916		efficiencymapContam->Write();
917		efficiencymap->Write();
918	+	if(ibin==0) centermap->Write();
919	+	if(ibin==0) widthmap->Write();
920		flipmap->Write();
921		efficiencymapmet->Write();
922		efficiencymapAcc->Write();
#	Line 869 | Line 954 | write_warning(__FUNCTION__,"CURRENTLY SW
954		delete syspdfmap;
955		delete systotmap;
956		delete sysstatmap;
957	+	delete XSmap;
958		delete limcanvas;
959		}
960
#	Line 877 | Line 963 | void scan_SUSY_parameter_space(string mc
963		for(int ibin=0;ibin<jzb_cut.size();ibin++) {
964		scan_SUSY_parameter_space(mcjzb,datajzb,jzb_cut,requireZ, peakerror, ibin, njobs, jobnumber,systonly,effonly);
965		}
966	+	delete_any_cached_scans();// tidy up ...
967		}
968

Diff Legend

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