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.5 by buchmann, Sun Mar 18 17:53:29 2012 UTC

#	Line 94 | Line 94 | bool initialized_t2=false;
94
95		int srmcpretries=0;
96
97	<	void load_scan_sample(int a, int b, int &scanfileindex,int scantype,bool isretry=false) {
97	>	void load_scan_sample(int a, int b, int &scanfileindex,int scantype,int mglu, int mlsp, bool isretry=false) {
98
99		// 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;
101	>	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;
102		stringstream filetoload;
103		char hostname[1023];
104		gethostname(hostname,1023);
105	<
105	>	string samplename;
106		if((Contains(hostname,"t3ui")||Contains(hostname,"t3wn"))) { // local case
107		filetoload << "/shome/buchmann/ntuples/"<<PlottingSetup::ScanSampleDirectory;
108	<	if(scantype==mSUGRA) filetoload << "/mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
109	<	if(scantype==SMS) filetoload << "/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
110	<	if(scantype==GMSB) filetoload << "/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
111	<	if(!Contains(((scansample.collection)[(scansample.collection).size()-1]).filename,"_"+any2string(a)+"_"+any2string(b))) {
108	>	if(scantype==mSUGRA) {
109	>	//filetoload << "/mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
110	>	filetoload.str("");
111	>	filetoload << "/shome/lbaeni/jzb/" << PlottingSetup::ScanSampleDirectory  << "/mSUGRA_M0_" << mglu << "__M12_" << mlsp << ".root";
112	>	samplename="mSUGRA";
113	>	//      filetoload << "/mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
114	>	}
115	>	if(scantype==SMS) {
116	>	filetoload << "/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
117	>	samplename="SMS";
118	>	}
119	>	if(scantype==GMSB) {
120	>	filetoload << "/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
121	>	samplename="GMSB";
122	>	}
123	>	if(scansample.collection.size()<1||!Contains(((scansample.collection)[(scansample.collection).size()-1]).filename,"_"+any2string(a)+"_"+any2string(b))) {
124		dout << "The last sample is NOT the same one as the current one, possibly popping off last one and adding the new one." << endl;
125	<	if((scansample.collection).size()>1) {
125	>	if((scansample.collection).size()>=1) {
126		scansample.RemoveLastSample();
127		}
128		scanfileindex=(scansample.collection).size();
129		//New: Loading file when necessary, not before (avoiding high memory usage and startup times)
130	<	if(scanfileindex!=0) scansample.AddSample(filetoload.str(),"scansample",1,1,false,true,scanfileindex,kRed);
130	>
131	>	scansample.AddSample(filetoload.str(),samplename,1,1,false,true,scanfileindex,kRed);
132	>	dout << "Just added the following file: " << filetoload.str() << endl;
133		} else {
134		dout << "Last sample is the same as the current one. Recycling it." << endl;
135		scanfileindex=(scansample.collection).size()-1;
#	Line 123 | Line 137 | void load_scan_sample(int a, int b, int
137		dout << " Going to use the following file: " << filetoload.str() << endl;
138		} // end of t3 case
139		else {
140	+	//CRAB case
141		write_info(__FUNCTION__,"Hello, CRAB! Might need to load files ...");
142		PlottingSetup::limitpatience=PlottingSetup::limitpatienceCRAB;
143		stringstream copyfile;
144	+	gSystem->Exec("mkdir -p local_storage");
145		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;
146	<	if(scantype==mSUGRA) copyfile << "/mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root " << "mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root ";
147	<	if(scantype==SMS) copyfile << "/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root " << "SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root ";
148	<	if(scantype==GMSB) copyfile << "/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root " << "GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root ";
146	>	if(scantype==mSUGRA) {
147	>	copyfile<< "/mSUGRA_M0_" << mglu << "__M12_" << mlsp << ".root local_storage/mSUGRA_M0_" << mglu << "__M12_" << mlsp << ".root ";
148	>	samplename="mSUGRA";
149	>	}
150	>	if(scantype==SMS) {
151	>	copyfile << "/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root " << "local_storage/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root ";
152	>	samplename="SMS";
153	>	}
154	>	if(scantype==GMSB) {
155	>	copyfile << "/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root " << "local_storage/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root ";
156	>	samplename="GMSB";
157	>	}
158		stringstream newfilename;
159	<	if(scantype==mSUGRA) newfilename << "mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
160	<	if(scantype==SMS) newfilename << "SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
161	<	if(scantype==GMSB) newfilename << "GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
159	>	//    if(scantype==mSUGRA) newfilename << "local_storage/mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
160	>	if(scantype==mSUGRA) newfilename << "local_storage/mSUGRA_M0_" << mglu << "__M12_" << mlsp << ".root";
161	>	if(scantype==SMS) newfilename << "local_storage/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
162	>	if(scantype==GMSB) newfilename << "local_storage/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
163
164		if(!Contains(((scansample.collection)[(scansample.collection).size()-1]).filename,"_"+any2string(a)+"_"+any2string(b))||(a==0&&b==0)&&!initialized_t2) {
165		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;
166	<	if((scansample.collection).size()>1) {
166	>	while((scansample.collection).size()>0) {
167		scansample.RemoveLastSample();
168		}
169	+	cout << "New scanfileindex stems from scansample size: " << (scansample.collection).size() << endl;
170		scanfileindex=(scansample.collection).size();
171		//New: Loading file when necessary, not before (avoiding high memory usage and startup times)
172	<	if(scanfileindex!=0||isretry) {
173	<	delete_any_cached_scans();
174	<	dout << "Going to download the scan file with the following copy command: " << copyfile.str() << endl;
175	<	gSystem->Exec(copyfile.str().c_str());
172	>	if(1) {
173	>	//        delete_any_cached_scans();
174	>	if(!doesROOTFileExist(newfilename.str())) {
175	>	dout << "Going to download the scan file with the following copy command: " << copyfile.str() << endl;
176	>	gSystem->Exec(copyfile.str().c_str()); // download it if it hasn't been downloaded before
177	>	}
178		if(doesROOTFileExist(newfilename.str())) {
179		initialized_t2=true;
180	<	scansample.AddSample(newfilename.str(),"scansample",1,1,false,true,scanfileindex,kRed);
180	>	cout << "Going to add the following sample: " << newfilename.str() << endl;
181	>	gSystem->Exec(("ls -ltrh "+newfilename.str()).c_str());
182	>	scansample.AddSample(newfilename.str(),samplename,1,1,false,true,scanfileindex,kRed);
183		} else {
184		srmcpretries++;
185		if(srmcpretries<5) {
186		dout << "The file could not be loaded correctly - retrying!" << endl;
187		sleep(5);
188	<	load_scan_sample(a,b,scanfileindex,scantype,true);
188	>	load_scan_sample(a,b,scanfileindex,scantype,mglu,mlsp,true);
189		} else {
190		dout << "Have tried 5 times to load this sample. Giving up now and failing the program execution" << endl;
191	<	assert(0);
191	>	assert(srmcpretries<5);
192		dout << "The command " << copyfile.str() << " failed to copy the file to the local storage (saving it as " << newfilename.str() << ")" << endl;
193		}
194		}
#	Line 168 | Line 199 | void load_scan_sample(int a, int b, int
199		}
200
201		}
202	+
203	+
204	+	// WATCH OUT THIS LINE NEEDS TO BE REMOVED
205	+	scanfileindex=0;
206		}
207
208	<	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)  {
208	>	/*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)  {
209		altxs=0;
210	<	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	<	}
210	>	for(int iproc=1;iproc<=10;iproc++) altxs+=GetXSecForPointAndChannel(mglu,mlsp,xsec,iproc);
211		return altxs;
212	<	}
212	>	}*/
213
214		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) {
215		bool runninglocally=true;
#	Line 206 | Line 236 | void establish_SUSY_limits(string mcjzb,
236		while ((key = (TKey*)nextkey()))
237		{
238		TObject *obj = key->ReadObj();
239	<	if(Contains((string)(obj->GetName()),"mSUGRA")) scantype=mSUGRA;
239	>	if(Contains((string)(obj->GetName()),"SUGRA")) scantype=mSUGRA;
240		if(Contains((string)(obj->GetName()),"GMSB")) scantype=GMSB;
241		}
242
#	Line 251 | Line 281 | void establish_SUSY_limits(string mcjzb,
281		TH2F *fullerr   = (TH2F*)fsyst->Get((prefix+"systotmap"+any2string(jzb_cut[ibin])).c_str());
282		TH2F *NEvents   = (TH2F*)fsyst->Get((prefix+"Neventsmap"+any2string(jzb_cut[ibin])).c_str());
283		TH2F *lflip     = (TH2F*)fsyst->Get((prefix+"flipmap"+any2string(jzb_cut[ibin])).c_str());
284	+	TH2F *XSmap     = (TH2F*)fsyst->Get((prefix+"XS"+any2string(jzb_cut[ibin])).c_str());
285
286		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
287		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 291 | void establish_SUSY_limits(string mcjzb,
291		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
292
293		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);
294		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);
295	+	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);
296
297
298		if(!fullerr || !mceff || !NEvents) {
#	Line 294 | Line 324 | void establish_SUSY_limits(string mcjzb,
324		continue;
325		}
326		vector<float> sigmas;
327	<	sigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped);
328	<
327	>	if(scantype!=mSUGRA) sigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped);
328	>	else {
329	>	//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"!
330	>	vector<float> asigmas;
331	>	asigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped,true); // asymptotic limit first
332	>	float strength=asigmas[0]/XSmap->GetBinContent(GlobalBin);
333	>	if(strength>0.5&&strength<2) {
334	>	sigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped); // asymptotic limit first
335	>	asymptoticmap->SetBinContent(GlobalBin,0);
336	>	}
337	>	else {
338	>	sigmas=asigmas;
339	>	asymptoticmap->SetBinContent(GlobalBin,1);
340	>	}
341	>	exclmap->SetBinContent(GlobalBin,strength);
342	>	}
343
344		if(sigmas[0]>0) limitmap->SetBinContent(GlobalBin,sigmas[0]); //anything else is an error code
345		if(sigmas.size()>1) {
#	Line 307 | Line 351 | void establish_SUSY_limits(string mcjzb,
351		}
352
353		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	–	}
354		}
355		}
356
#	Line 325 | Line 359 | void establish_SUSY_limits(string mcjzb,
359		outputfile->cd();
360		limitmap->Write();
361		flipmap->Write();
362	+	asymptoticmap->Write();
363		if(doexpected) {
364		explimitmap->Write();
365		exp1plimitmap->Write();
#	Line 332 | Line 367 | void establish_SUSY_limits(string mcjzb,
367		exp2plimitmap->Write();
368		exp2mlimitmap->Write();
369		}
370	<	if(scantype==mSUGRA) {
336	<	exclmap->Write();
337	<	xsmap->Write();
338	<	totxsmap->Write();
339	<	}
370	>	if(scantype==mSUGRA) exclmap->Write();
371		outputfile->Close();
372		delete limcanvas;
373		}
#	Line 378 | Line 409 | void scan_SUSY_parameter_space(string mc
409		string prefix="SMS_";
410		// up to here, everything is set up for SMS; now we need to switch stuff around if we're dealing with an mSUGRA scan!
411		int scantype=SMS;
412	<	if(Contains((scansample.collection)[0].samplename,"mSUGRA")) scantype=mSUGRA;
412	>	if(Contains((scansample.collection)[0].samplename,"SUGRA")) scantype=mSUGRA;
413		if(Contains((scansample.collection)[0].samplename,"GMSB")) scantype=GMSB;
414
415		if(scantype==mSUGRA) {
#	Line 430 | Line 461 | void scan_SUSY_parameter_space(string mc
461		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);
462		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);
463		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);
464	<	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);
464	>	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-
465	>	mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
466		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);
467		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);
468		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);
469		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);
470		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);
471		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);
472	+	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);
473
474		TH2F *imposedxmap;
475		TH2F *realxmap;
476
477	+	TH2F *centermap;
478	+	TH2F *widthmap;
479	+
480		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);
481
482		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 484 | void scan_SUSY_parameter_space(string mc
484		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);
485
486		if(ibin==0) {
487	+	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);
488	+	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);
489	+
490		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);
491		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);
492		}
#	Line 485 | Line 524 | write_warning(__FUNCTION__,"CURRENTLY SW
524		for(int mglu=(int)mglustart;mglu<=mgluend;mglu+=(int)mglustep) {
525		for (int mlsp=(int)mLSPstart;mlsp<=mLSPend&&(scantype==mSUGRA||mlsp<=mglu);mlsp+=(int)mLSPstep) Npoints++;
526		}
527	<
527	>
528	>	map <  pair<float, float>, map<string, float>  >  xsec;
529	>	if(scantype==mSUGRA) {
530	>	xsec=getXsec(PlottingSetup::mSUGRAxsFile);
531	>	}
532	>
533		int ipoint=-1;
534		for(int mglu=(int)mglustart;mglu<=mgluend;mglu+=(int)mglustep) {
535		for (int mlsp=(int)mLSPstart;mlsp<=mLSPend&&(scantype==mSUGRA||mlsp<=mglu);mlsp+=(int)mLSPstep)
#	Line 500 | Line 544 | write_warning(__FUNCTION__,"CURRENTLY SW
544		int a = int((PlottingSetup::ScanXzones*mlsp)/(m12trees+1));
545		int b = int((PlottingSetup::ScanYzones*mglu)/(m0trees+1));
546		int scanfileindex=PlottingSetup::ScanYzones*a+b;
547	<	load_scan_sample(a,b,scanfileindex,scantype);
547	>	load_scan_sample(a,b,scanfileindex,scantype,mglu,mlsp);
548
549		clock_t start,finish;
550		start = clock(); // starting the clock to measure how long the computation takes!
#	Line 510 | Line 554 | write_warning(__FUNCTION__,"CURRENTLY SW
554		float nevents = (scansample.collection)[scanfileindex].events->GetSelectedRows();
555		vector<vector<float> > systematics;
556		if(nevents<10) {
557	<	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;
557	>	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;
558		continue;
559		} else {
560		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 596 | write_warning(__FUNCTION__,"CURRENTLY SW
596		(scansample.collection)[scanfileindex].events->Draw("imposedx>>realxhisto",addcut.str().c_str(),"goff");
597		imposedxmap->SetBinContent(GlobalBin,realxhisto->GetMean());
598		imposedxmap->SetBinError(GlobalBin,realxhisto->GetRMS());
599	+	stringstream specialexpression;
600	+	specialexpression << mcjzb << ">>jzbhisto";
601	+	TH1F *jzbhisto = new TH1F("jzbhisto","jzbhisto",100,-500,500);
602	+	(scansample.collection)[scanfileindex].events->Draw(specialexpression.str().c_str(),addcut.str().c_str(),"goff");
603	+	centermap->SetBinContent(GlobalBin,jzbhisto->GetMean());
604	+	widthmap->SetBinContent(GlobalBin,jzbhisto->GetRMS());
605		delete realxhisto;
606	+	delete jzbhisto;
607		}
608
609
610	+	if(scantype==mSUGRA) {
611	+	float absxs=0;
612	+	for(int i=0;i<12;i++) absxs+=GetXSecForPointAndChannel(mglu,mlsp,xsec,i);
613	+	XSmap->SetBinContent(GlobalBin,absxs);
614	+	}
615	+
616		if(nevents!=0&&(efficiencyonly||systematicsonly)) {
617	<	Value effwosigcont = MCefficiency((scansample.collection)[scanfileindex].events,result,resulterr,flipped,mcjzb,requireZ,(int)nevents,addcut.str(),-1);
617	>	Value effwosigcont = MCefficiency((scansample.collection)[scanfileindex].events,result,resulterr,flipped,mcjzb,requireZ,(int)nevents,scantype,xsec,addcut.str(),-1);
618		if(result<0&&allowflipping) {
619		write_info(__FUNCTION__,"Found negative efficiency (signal distribution populating JZB<0 more than JZB>0) - flipping analysis!");
620		flipped=1;
621	<	MCefficiency((scansample.collection)[scanfileindex].events,result,resulterr,flipped,mcjzb,requireZ,(int)nevents,addcut.str(),-1);
621	>	MCefficiency((scansample.collection)[scanfileindex].events,result,resulterr,flipped,mcjzb,requireZ,(int)nevents,scantype,xsec,addcut.str(),-1);
622		}
623		efficiencymap->SetBinContent(GlobalBin,result);
624		efficiencymap->SetBinError(GlobalBin,resulterr);
#	Line 594 | Line 651 | write_warning(__FUNCTION__,"CURRENTLY SW
651		ipointmap->SetBinContent(GlobalBin,ipoint);
652		if(efficiencyonly) continue;
653
654	<	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
654	>	do_systematics_for_one_file((scansample.collection)[scanfileindex].events,(int)nevents,"SUSY SCAN", systematics,flipped, xsec, mcjzb,datajzb,peakerror,requireZ, addcut.str(),true,scantype);
655		float JZBcutat = systematics[0][0];
656		float mceff    = systematics[0][1];
657		float mcefferr = systematics[0][2];//MC stat error
#	Line 606 | Line 663 | write_warning(__FUNCTION__,"CURRENTLY SW
663		float mcwoscefr= systematics[0][9]; // error on efficiency without signal contamination
664		float sys_pdf   = 0;
665		if(systematics[0].size()>10) sys_pdf = systematics[0][10]; // PDF
666	<
666	>	cout << "Going to store: total error is: " << toterr << " and the relative tot error, i.e. total error / mc efficiency is " << toterr/mceff << endl;
667		if(mceff!=mceff||toterr!=toterr||mceff<0 && (!systematicsonly&&!efficiencyonly)) {
668		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;
669		continue;
#	Line 688 | Line 745 | write_warning(__FUNCTION__,"CURRENTLY SW
745		sysjsumap->Write();
746		sysresmap->Write();
747		efficiencymap->Write();
748	+	if(ibin==0) centermap->Write();
749	+	if(ibin==0) widthmap->Write();
750		flipmap->Write();
751		efficiencymapmet->Write();
752		efficiencymapAcc->Write();
#	Line 699 | Line 758 | write_warning(__FUNCTION__,"CURRENTLY SW
758		syspdfmap->Write();
759		systotmap->Write();
760		sysstatmap->Write();
761	+	XSmap->Write();
762		if(ibin==0) imposedxmap->Write();
763		if(ibin==0) realxmap->Write();
764		Neventsmap->Write();
765		ipointmap->Write();
766		timemap->Write();
707	–	outputfile->Close();
767		for(int i=0;i<susy_Zdecay_origin.size();i++) {
768		h_susy_Zdecay_origin[i]->Write();
769		}
770	+	outputfile->Close();
771		}
772		}
773		if(systematicsonly) { // systematics only :
#	Line 746 | Line 806 | write_warning(__FUNCTION__,"CURRENTLY SW
806		sysjsumap->Write();
807		sysresmap->Write();
808		flipmap->Write();
809	+	if(ibin==0) centermap->Write();
810	+	if(ibin==0) widthmap->Write();
811		efficiencymap->Write();
812		efficiencymapmet->Write();
813		efficiencymapAcc->Write();
#	Line 759 | Line 821 | write_warning(__FUNCTION__,"CURRENTLY SW
821		syspdfmap->Write();
822		systotmap->Write();
823		sysstatmap->Write();
824	+	XSmap->Write();
825		if(ibin==0) imposedxmap->Write();
826		if(ibin==0) realxmap->Write();
827		timemap->Write();
765	–	outputfile->Close();
828		for(int i=0;i<susy_Zdecay_origin.size();i++) {
829		h_susy_Zdecay_origin[i]->Write();
830		}
831	+	outputfile->Close();
832		}
833		}//end of systematics only
834		if(efficiencyonly) {
#	Line 832 | Line 895 | write_warning(__FUNCTION__,"CURRENTLY SW
895		noscefficiencymap->Write();
896		efficiencymapContam->Write();
897		efficiencymap->Write();
898	+	if(ibin==0) centermap->Write();
899	+	if(ibin==0) widthmap->Write();
900		flipmap->Write();
901		efficiencymapmet->Write();
902		efficiencymapAcc->Write();
#	Line 869 | Line 934 | write_warning(__FUNCTION__,"CURRENTLY SW
934		delete syspdfmap;
935		delete systotmap;
936		delete sysstatmap;
937	+	delete XSmap;
938		delete limcanvas;
939		}
940

Diff Legend

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