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.4 by buchmann, Sun Mar 18 11:53:56 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 126 | Line 140 | void load_scan_sample(int a, int b, int
140		write_info(__FUNCTION__,"Hello, CRAB! Might need to load files ...");
141		PlottingSetup::limitpatience=PlottingSetup::limitpatienceCRAB;
142		stringstream copyfile;
143	+	gSystem->Exec("mkdir -p local_storage");
144		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;
145	<	if(scantype==mSUGRA) copyfile << "/mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root " << "mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root ";
146	<	if(scantype==SMS) copyfile << "/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root " << "SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root ";
147	<	if(scantype==GMSB) copyfile << "/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root " << "GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root ";
145	>	if(scantype==mSUGRA) {
146	>	copyfile << "/mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root " << "local_storage/mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root ";
147	>	samplename="mSUGRA";
148	>	}
149	>	if(scantype==SMS) {
150	>	copyfile << "/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root " << "local_storage/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root ";
151	>	samplename="SMS";
152	>	}
153	>	if(scantype==GMSB) {
154	>	copyfile << "/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root " << "local_storage/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root ";
155	>	samplename="GMSB";
156	>	}
157		stringstream newfilename;
158	<	if(scantype==mSUGRA) newfilename << "mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
159	<	if(scantype==SMS) newfilename << "SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
160	<	if(scantype==GMSB) newfilename << "GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
158	>	if(scantype==mSUGRA) newfilename << "local_storage/mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
159	>	if(scantype==SMS) newfilename << "local_storage/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
160	>	if(scantype==GMSB) newfilename << "local_storage/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
161
162		if(!Contains(((scansample.collection)[(scansample.collection).size()-1]).filename,"_"+any2string(a)+"_"+any2string(b))||(a==0&&b==0)&&!initialized_t2) {
163		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;
#	Line 143 | Line 167 | void load_scan_sample(int a, int b, int
167		scanfileindex=(scansample.collection).size();
168		//New: Loading file when necessary, not before (avoiding high memory usage and startup times)
169		if(scanfileindex!=0||isretry) {
170	<	delete_any_cached_scans();
170	>	//        delete_any_cached_scans();
171		dout << "Going to download the scan file with the following copy command: " << copyfile.str() << endl;
172	<	gSystem->Exec(copyfile.str().c_str());
172	>	if(!doesROOTFileExist(newfilename.str())) gSystem->Exec(copyfile.str().c_str()); // download it if it hasn't been downloaded before
173		if(doesROOTFileExist(newfilename.str())) {
174		initialized_t2=true;
175		scansample.AddSample(newfilename.str(),"scansample",1,1,false,true,scanfileindex,kRed);
#	Line 154 | Line 178 | void load_scan_sample(int a, int b, int
178		if(srmcpretries<5) {
179		dout << "The file could not be loaded correctly - retrying!" << endl;
180		sleep(5);
181	<	load_scan_sample(a,b,scanfileindex,scantype,true);
181	>	load_scan_sample(a,b,scanfileindex,scantype,mglu,mlsp,true);
182		} else {
183		dout << "Have tried 5 times to load this sample. Giving up now and failing the program execution" << endl;
184	<	assert(0);
184	>	assert(srmcpretries<5);
185		dout << "The command " << copyfile.str() << " failed to copy the file to the local storage (saving it as " << newfilename.str() << ")" << endl;
186		}
187		}
#	Line 168 | Line 192 | void load_scan_sample(int a, int b, int
192		}
193
194		}
195	+
196	+
197	+	// WATCH OUT THIS LINE NEEDS TO BE REMOVED
198	+	scanfileindex=0;
199		}
200
201	<	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)  {
201	>	/*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)  {
202		altxs=0;
203	<	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	<	}
203	>	for(int iproc=1;iproc<=10;iproc++) altxs+=GetXSecForPointAndChannel(mglu,mlsp,xsec,iproc);
204		return altxs;
205	<	}
205	>	}*/
206
207		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) {
208		bool runninglocally=true;
#	Line 206 | Line 229 | void establish_SUSY_limits(string mcjzb,
229		while ((key = (TKey*)nextkey()))
230		{
231		TObject *obj = key->ReadObj();
232	<	if(Contains((string)(obj->GetName()),"mSUGRA")) scantype=mSUGRA;
232	>	if(Contains((string)(obj->GetName()),"SUGRA")) scantype=mSUGRA;
233		if(Contains((string)(obj->GetName()),"GMSB")) scantype=GMSB;
234		}
235
#	Line 251 | Line 274 | void establish_SUSY_limits(string mcjzb,
274		TH2F *fullerr   = (TH2F*)fsyst->Get((prefix+"systotmap"+any2string(jzb_cut[ibin])).c_str());
275		TH2F *NEvents   = (TH2F*)fsyst->Get((prefix+"Neventsmap"+any2string(jzb_cut[ibin])).c_str());
276		TH2F *lflip     = (TH2F*)fsyst->Get((prefix+"flipmap"+any2string(jzb_cut[ibin])).c_str());
277	+	TH2F *XSmap     = (TH2F*)fsyst->Get((prefix+"XS"+any2string(jzb_cut[ibin])).c_str());
278
279		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
280		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 294 | Line 318 | void establish_SUSY_limits(string mcjzb,
318		continue;
319		}
320		vector<float> sigmas;
321	<	sigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped);
321	>	if(scantype!=mSUGRA) sigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped);
322	>	else {
323	>	//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"!
324	>	vector<float> asigmas;
325	>	asigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped,true); // asymptotic limit first
326	>	float strength=asigmas[0]/XSmap->GetBinContent(GlobalBin);
327	>	if(strength>0.5&&strength<2) sigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped); // asymptotic limit first
328	>	else sigmas=asigmas;
329	>	exclmap->SetBinContent(GlobalBin,strength);
330	>	xsmap->SetBinContent(GlobalBin,XSmap->GetBinContent(GlobalBin));
331	>	}
332	>
333	>
334
335
336		if(sigmas[0]>0) limitmap->SetBinContent(GlobalBin,sigmas[0]); //anything else is an error code
#	Line 307 | Line 343 | void establish_SUSY_limits(string mcjzb,
343		}
344
345		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	–	}
346		}
347		}
348
#	Line 378 | Line 404 | void scan_SUSY_parameter_space(string mc
404		string prefix="SMS_";
405		// up to here, everything is set up for SMS; now we need to switch stuff around if we're dealing with an mSUGRA scan!
406		int scantype=SMS;
407	<	if(Contains((scansample.collection)[0].samplename,"mSUGRA")) scantype=mSUGRA;
407	>	if(Contains((scansample.collection)[0].samplename,"SUGRA")) scantype=mSUGRA;
408		if(Contains((scansample.collection)[0].samplename,"GMSB")) scantype=GMSB;
409
410		if(scantype==mSUGRA) {
#	Line 430 | Line 456 | void scan_SUSY_parameter_space(string mc
456		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);
457		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);
458		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);
459	<	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);
459	>	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-
460	>	mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
461		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);
462		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);
463		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);
464		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);
465		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);
466		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);
467	+	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);
468
469		TH2F *imposedxmap;
470		TH2F *realxmap;
471
472	+	TH2F *centermap;
473	+	TH2F *widthmap;
474	+
475		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);
476
477		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 479 | void scan_SUSY_parameter_space(string mc
479		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);
480
481		if(ibin==0) {
482	+	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);
483	+	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);
484	+
485		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);
486		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);
487		}
#	Line 485 | Line 519 | write_warning(__FUNCTION__,"CURRENTLY SW
519		for(int mglu=(int)mglustart;mglu<=mgluend;mglu+=(int)mglustep) {
520		for (int mlsp=(int)mLSPstart;mlsp<=mLSPend&&(scantype==mSUGRA||mlsp<=mglu);mlsp+=(int)mLSPstep) Npoints++;
521		}
522	<
522	>
523	>	map <  pair<float, float>, map<string, float>  >  xsec;
524	>	if(scantype==mSUGRA) {
525	>	xsec=getXsec(PlottingSetup::mSUGRAxsFile);
526	>	}
527	>
528		int ipoint=-1;
529		for(int mglu=(int)mglustart;mglu<=mgluend;mglu+=(int)mglustep) {
530		for (int mlsp=(int)mLSPstart;mlsp<=mLSPend&&(scantype==mSUGRA||mlsp<=mglu);mlsp+=(int)mLSPstep)
#	Line 500 | Line 539 | write_warning(__FUNCTION__,"CURRENTLY SW
539		int a = int((PlottingSetup::ScanXzones*mlsp)/(m12trees+1));
540		int b = int((PlottingSetup::ScanYzones*mglu)/(m0trees+1));
541		int scanfileindex=PlottingSetup::ScanYzones*a+b;
542	<	load_scan_sample(a,b,scanfileindex,scantype);
542	>	load_scan_sample(a,b,scanfileindex,scantype,mglu,mlsp);
543
544		clock_t start,finish;
545		start = clock(); // starting the clock to measure how long the computation takes!
#	Line 510 | Line 549 | write_warning(__FUNCTION__,"CURRENTLY SW
549		float nevents = (scansample.collection)[scanfileindex].events->GetSelectedRows();
550		vector<vector<float> > systematics;
551		if(nevents<10) {
552	<	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;
552	>	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;
553		continue;
554		} else {
555		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 591 | write_warning(__FUNCTION__,"CURRENTLY SW
591		(scansample.collection)[scanfileindex].events->Draw("imposedx>>realxhisto",addcut.str().c_str(),"goff");
592		imposedxmap->SetBinContent(GlobalBin,realxhisto->GetMean());
593		imposedxmap->SetBinError(GlobalBin,realxhisto->GetRMS());
594	+	stringstream specialexpression;
595	+	specialexpression << mcjzb << ">>jzbhisto";
596	+	TH1F *jzbhisto = new TH1F("jzbhisto","jzbhisto",100,-500,500);
597	+	(scansample.collection)[scanfileindex].events->Draw(specialexpression.str().c_str(),addcut.str().c_str(),"goff");
598	+	centermap->SetBinContent(GlobalBin,jzbhisto->GetMean());
599	+	widthmap->SetBinContent(GlobalBin,jzbhisto->GetRMS());
600		delete realxhisto;
601	+	delete jzbhisto;
602		}
603
604
605		if(nevents!=0&&(efficiencyonly||systematicsonly)) {
606	<	Value effwosigcont = MCefficiency((scansample.collection)[scanfileindex].events,result,resulterr,flipped,mcjzb,requireZ,(int)nevents,addcut.str(),-1);
606	>	Value effwosigcont = MCefficiency((scansample.collection)[scanfileindex].events,result,resulterr,flipped,mcjzb,requireZ,(int)nevents,scantype,xsec,addcut.str(),-1);
607		if(result<0&&allowflipping) {
608		write_info(__FUNCTION__,"Found negative efficiency (signal distribution populating JZB<0 more than JZB>0) - flipping analysis!");
609		flipped=1;
610	<	MCefficiency((scansample.collection)[scanfileindex].events,result,resulterr,flipped,mcjzb,requireZ,(int)nevents,addcut.str(),-1);
610	>	MCefficiency((scansample.collection)[scanfileindex].events,result,resulterr,flipped,mcjzb,requireZ,(int)nevents,scantype,xsec,addcut.str(),-1);
611		}
612		efficiencymap->SetBinContent(GlobalBin,result);
613		efficiencymap->SetBinError(GlobalBin,resulterr);
#	Line 594 | Line 640 | write_warning(__FUNCTION__,"CURRENTLY SW
640		ipointmap->SetBinContent(GlobalBin,ipoint);
641		if(efficiencyonly) continue;
642
643	<	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
643	>	do_systematics_for_one_file((scansample.collection)[scanfileindex].events,(int)nevents,"SUSY SCAN", systematics,flipped, xsec, mcjzb,datajzb,peakerror,requireZ, addcut.str(),true,scantype);
644		float JZBcutat = systematics[0][0];
645		float mceff    = systematics[0][1];
646		float mcefferr = systematics[0][2];//MC stat error
#	Line 606 | Line 652 | write_warning(__FUNCTION__,"CURRENTLY SW
652		float mcwoscefr= systematics[0][9]; // error on efficiency without signal contamination
653		float sys_pdf   = 0;
654		if(systematics[0].size()>10) sys_pdf = systematics[0][10]; // PDF
655	<
655	>	cout << "Going to store: total error is: " << toterr << " and the relative tot error, i.e. total error / mc efficiency is " << toterr/mceff << endl;
656		if(mceff!=mceff||toterr!=toterr||mceff<0 && (!systematicsonly&&!efficiencyonly)) {
657		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;
658		continue;
#	Line 688 | Line 734 | write_warning(__FUNCTION__,"CURRENTLY SW
734		sysjsumap->Write();
735		sysresmap->Write();
736		efficiencymap->Write();
737	+	if(ibin==0) centermap->Write();
738	+	if(ibin==0) widthmap->Write();
739		flipmap->Write();
740		efficiencymapmet->Write();
741		efficiencymapAcc->Write();
#	Line 699 | Line 747 | write_warning(__FUNCTION__,"CURRENTLY SW
747		syspdfmap->Write();
748		systotmap->Write();
749		sysstatmap->Write();
750	+	XSmap->Write();
751		if(ibin==0) imposedxmap->Write();
752		if(ibin==0) realxmap->Write();
753		Neventsmap->Write();
754		ipointmap->Write();
755		timemap->Write();
707	–	outputfile->Close();
756		for(int i=0;i<susy_Zdecay_origin.size();i++) {
757		h_susy_Zdecay_origin[i]->Write();
758		}
759	+	outputfile->Close();
760		}
761		}
762		if(systematicsonly) { // systematics only :
#	Line 746 | Line 795 | write_warning(__FUNCTION__,"CURRENTLY SW
795		sysjsumap->Write();
796		sysresmap->Write();
797		flipmap->Write();
798	+	if(ibin==0) centermap->Write();
799	+	if(ibin==0) widthmap->Write();
800		efficiencymap->Write();
801		efficiencymapmet->Write();
802		efficiencymapAcc->Write();
#	Line 759 | Line 810 | write_warning(__FUNCTION__,"CURRENTLY SW
810		syspdfmap->Write();
811		systotmap->Write();
812		sysstatmap->Write();
813	+	XSmap->Write();
814		if(ibin==0) imposedxmap->Write();
815		if(ibin==0) realxmap->Write();
816		timemap->Write();
765	–	outputfile->Close();
817		for(int i=0;i<susy_Zdecay_origin.size();i++) {
818		h_susy_Zdecay_origin[i]->Write();
819		}
820	+	outputfile->Close();
821		}
822		}//end of systematics only
823		if(efficiencyonly) {
#	Line 832 | Line 884 | write_warning(__FUNCTION__,"CURRENTLY SW
884		noscefficiencymap->Write();
885		efficiencymapContam->Write();
886		efficiencymap->Write();
887	+	if(ibin==0) centermap->Write();
888	+	if(ibin==0) widthmap->Write();
889		flipmap->Write();
890		efficiencymapmet->Write();
891		efficiencymapAcc->Write();
#	Line 869 | Line 923 | write_warning(__FUNCTION__,"CURRENTLY SW
923		delete syspdfmap;
924		delete systotmap;
925		delete sysstatmap;
926	+	delete XSmap;
927		delete limcanvas;
928		}
929

Diff Legend

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