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.7 by buchmann, Fri Mar 23 12:04:29 2012 UTC vs.
Revision 1.13 by buchmann, Fri Apr 27 08:46:00 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
28		using namespace PlottingSetup;
29
30	<	namespace SUSYScanSpace {
31	<	vector<string> loaded_files;
32	<	int SUSYscantype;
30	>	void scan_SUSY_parameter_space(string mcjzb,string datajzb,vector<float> jzb_cut,bool requireZ, float peakerror, float peakerrordata, float njobs, float jobnumber, bool systonly, bool effonly,bool shapeanalysis);
31	>
32	>	bool isValidShapeResult(ShapeDroplet droplet) {
33	>	if(droplet.observed>50000) {
34	>	// error code of ShapeLimit/CreateModel.sh, invoked if there were too many errors during computation
35	>	write_warning(__FUNCTION__,"Shape limit was ridiculously high ("+any2string(droplet.observed)+"); this is a typical error signature of the combine tool. Shape result will be discarded, and c&c will be carried out instead.");
36	>	return false;
37	>	}
38	>	if(droplet.observed==-12345) {
39	>	write_warning(__FUNCTION__,"Shape limit was ("+any2string(droplet.observed)+"); this is the error signature of CreateModel.sh / ShapeDroplet. Shape result will be discarded, and c&c will be carried out instead.");
40	>	return false;
41	>	}
42	>	if(droplet.expected>-1&&droplet.expected==0&&droplet.observed>5000) {
43	>	write_warning(__FUNCTION__,"Observed shape limit was ("+any2string(droplet.observed)+") while expected limit was "+any2string(droplet.expected)+"; this is a typical error signature of the combine tool. Shape result will be discarded, and c&c will be carried out instead.");
44	>	return false;
45	>	}
46	>	if(droplet.expected>-1&&droplet.expected/droplet.observed>50) {
47	>	write_warning(__FUNCTION__,"Observed shape limit was ("+any2string(droplet.observed)+") while expected limit was "+any2string(droplet.expected)+"; such a high ratio between expected & observed limit is a typical error signature of the combine tool. Shape result will be discarded, and c&c will be carried out instead.");
48	>	return false;
49	>	}
50	>	if(droplet.expected>-1&&droplet.observed/droplet.expected>50) {
51	>	write_warning(__FUNCTION__,"Observed shape limit was ("+any2string(droplet.observed)+") while expected limit was "+any2string(droplet.expected)+"; such a high ratio between observed & expected limit is a typical error signature of the combine tool. Shape result will be discarded, and c&c will be carried out instead.");
52	>	return false;
53	>	}
54	>	return true;
55		}
56
57		void susy_scan_axis_labeling(TH2F *histo) {
#	Line 72 | Line 97 | void set_SUSY_style() {
97		}
98
99		bool delete_any_cached_scans() {
100	+	char hostname[1023];
101	+	gethostname(hostname,1023);
102
103
104		/*
105		//This should only be called via CRAB
106		cout << "   Deleting all cached files." << endl;
107		gSystem->Exec("ls -ltrh | grep root");
81	–	char hostname[1023];
82	–	gethostname(hostname,1023);
108		if(!((Contains(hostname,"t3ui")||Contains(hostname,"t3wn"))))
109		{
110		vector<string> all_files;
#	Line 97 | Line 122 | bool delete_any_cached_scans() {
122		return true;
123		}
124		*/
125	<	if(SUSYScanSpace::SUSYscantype==mSUGRA) {
125	>	if(!((Contains(hostname,"t3ui")||Contains(hostname,"t3wn")))) {
126		cout << "Going to purge files in local_storage that have been loaded!" << endl;
127		for(int i=0;i<SUSYScanSpace::loaded_files.size();i++) {
128		gSystem->Exec(((string)"rm "+SUSYScanSpace::loaded_files[i]).c_str());
#	Line 111 | Line 136 | bool initialized_t2=false;
136
137		int srmcpretries=0;
138
139	<	void load_scan_sample(int a, int b, int &scanfileindex,int scantype,int mglu, int mlsp, bool isretry=false) {
140	<
141	<	// 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 << " 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	<	string samplename;
123	<	if((Contains(hostname,"t3ui")||Contains(hostname,"t3wn"))&&scantype!=mSUGRA) { // local case
139	>	void load_local_scan_sample(int a, int b, int &scanfileindex,int scantype,int mglu, int mlsp, bool isretry=false) {
140	>	stringstream filetoload;
141	>	string samplename;
142		filetoload << "/shome/buchmann/ntuples/"<<PlottingSetup::ScanSampleDirectory;
143		if(scantype==mSUGRA) {
144		//filetoload << "/mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
#	Line 137 | Line 155 | void load_scan_sample(int a, int b, int
155		filetoload << "/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
156		samplename="GMSB";
157		}
158	+
159	+	if(scantype==SMS && Contains(PlottingSetup::ScanSampleDirectory,"T1lh")) {
160	+	filetoload.str("");
161	+	filetoload << "/shome/lbaeni/jzb/" << PlottingSetup::ScanSampleDirectory  << "/SMS_MassGlu_" << mglu << "__MassLSP_" << mlsp << ".root";
162	+	samplename="SMS";
163	+	}
164	+
165		if(scansample.collection.size()<1||!Contains(((scansample.collection)[(scansample.collection).size()-1]).filename,"_"+any2string(a)+"_"+any2string(b))) {
166		dout << "The last sample is NOT the same one as the current one, possibly popping off last one and adding the new one." << endl;
167		if((scansample.collection).size()>=1) {
#	Line 152 | Line 177 | void load_scan_sample(int a, int b, int
177		scanfileindex=(scansample.collection).size()-1;
178		}
179		dout << " Going to use the following file: " << filetoload.str() << endl;
180	<	} // end of t3 case
181	<	else {
182	<	//CRAB case
180	>	}
181	>
182	>	void load_remote_scan_sample(int a, int b, int &scanfileindex,int scantype,int mglu, int mlsp, bool isretry=false) {
183		write_info(__FUNCTION__,"Hello, CRAB! Might need to load files ...");
184	+	string samplename;
185		PlottingSetup::limitpatience=PlottingSetup::limitpatienceCRAB;
186		stringstream copyfile;
187		gSystem->Exec("mkdir -p local_storage");
#	Line 204 | Line 230 | void load_scan_sample(int a, int b, int
230		if(srmcpretries<5) {
231		dout << "The file could not be loaded correctly - retrying!" << endl;
232		sleep(5);
233	<	load_scan_sample(a,b,scanfileindex,scantype,mglu,mlsp,true);
233	>	load_remote_scan_sample(a,b,scanfileindex,scantype,mglu,mlsp,true);
234		} else {
235		dout << "Have tried 5 times to load this sample. Giving up now and failing the program execution" << endl;
236		assert(srmcpretries<5);
#	Line 216 | Line 242 | void load_scan_sample(int a, int b, int
242		dout << "Last sample is the same as the current one. Recycling it." << endl;
243		scanfileindex=(scansample.collection).size()-1;
244		}
245	<
246	<	}
247	<
245	>	}
246	>
247	>	void load_scan_sample(int a, int b, int &scanfileindex,int scantype,int mglu, int mlsp, bool isretry=false) {
248
249	<	// WATCH OUT THIS LINE NEEDS TO BE REMOVED
250	<	scanfileindex=0;
249	>	// There is no need to define your sample here. That is done in Setup.C where you define the loading directory!
250	>
251	>	SUSYScanSpace::SUSYscantype=scantype;
252	>	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;
253	>	char hostname[1023];
254	>	gethostname(hostname,1023);
255	>	if((Contains(hostname,"t3ui")||Contains(hostname,"t3wn"))) load_local_scan_sample(a,b,scanfileindex,scantype,mglu,mlsp,isretry);
256	>	else load_remote_scan_sample(a,b,scanfileindex,scantype,mglu,mlsp,isretry);
257		}
258
259		/*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 265 | Line 297 | void establish_SUSY_limits(string mcjzb,
297		massgluname="M0"; // this is the "x axis" in the limit plot (like the gluino in the SMS case)
298		massLSPname="M12"; // this is the "y axis" in the limit plot (like the LSP in the SMS case)
299		mglustart=m0start;
300	<	xsec=getXsec(PlottingSetup::cbafbasedir+"/Plotting/Modules/external/scale_xsection_nlo1.0_m0_m12_10_0_1v1.txt");
300	>	xsec=getXsec(PlottingSetup::cbafbasedir+"/"+mSUGRAxsFile);
301	>
302		mgluend=m0end;
303		mglustep=m0step;
304		mLSPstart=m12start;
#	Line 273 | Line 306 | void establish_SUSY_limits(string mcjzb,
306		mLSPstep=m12step;
307		prefix="mSUGRA_";
308		dout << "mSUGRA scan has been set up." << endl;
276	–	//    xsec=getXsec("/scratch/buchmann/C/scale_xsection_nlo1.0_m0_m12_10_0_1v1.txt");
309		}
310		if(scantype==SMS) {
311		dout << "SMS scan has been set up." << endl;
#	Line 343 | Line 375 | void establish_SUSY_limits(string mcjzb,
375		continue;
376		}
377		vector<float> sigmas;
378	<	if(scantype!=mSUGRA) sigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped);
378	>	if(scantype!=mSUGRA) {
379	>	sigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped,true);
380	>	write_warning(__FUNCTION__,"Temporarily doing asymptotic limits!!!!");
381	>	}
382		else {
383		//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"!
384		vector<float> asigmas;
385		asigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped,true); // asymptotic limit first
386		float strength=asigmas[0]/XSmap->GetBinContent(GlobalBin);
387	<	if(strength>0.5&&strength<2) {
387	>	/*      if(strength>0.5&&strength<2) {
388		sigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped); // asymptotic limit first
389		asymptoticmap->SetBinContent(GlobalBin,0);
390	<	}
356	<	else {
390	>	} else {*/
391		sigmas=asigmas;
392		asymptoticmap->SetBinContent(GlobalBin,1);
393	<	}
393	>	/*      }*/
394		exclmap->SetBinContent(GlobalBin,strength);
395		}
396
#	Line 378 | Line 412 | void establish_SUSY_limits(string mcjzb,
412		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 :-)
413		outputfile->cd();
414		limitmap->Write();
415	+	XSmap->Write();
416		flipmap->Write();
417		asymptoticmap->Write();
418		if(doexpected) {
#	Line 407 | Line 442 | void establish_SUSY_limits(string mcjzb,
442
443
444
445	<	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) {
445	>	void scan_SUSY_parameter_space_for_one_cut(string mcjzb,string datajzb,vector<float> jzb_cut,bool requireZ, float peakerror, float peakerrordata, int ibin,float njobs=-1, float jobnumber=-1, bool systematicsonly=false,bool efficiencyonly=false, bool shapeanalysis=false) {
446	>
447	>	if(shapeanalysis&&ibin>0) return; // pointless since we're doing a shape analysis :-)
448	>
449		bool runninglocally=true;
450		if(njobs>-1&&jobnumber>-1) {
451		runninglocally=false;
#	Line 419 | Line 457 | void scan_SUSY_parameter_space(string mc
457		}
458
459		jzbSel=jzb_cut[ibin];
460	+	if(shapeanalysis) jzbSel=0;
461		geqleq="geq";
462		automatized=true;
463		mcjzbexpression=mcjzb;
#	Line 429 | Line 468 | void scan_SUSY_parameter_space(string mc
468		string prefix="SMS_";
469		// up to here, everything is set up for SMS; now we need to switch stuff around if we're dealing with an mSUGRA scan!
470		int scantype=SMS;
471	<	if(Contains((scansample.collection)[0].samplename,"SUGRA")) scantype=mSUGRA;
472	<	if(Contains((scansample.collection)[0].samplename,"GMSB")) scantype=GMSB;
471	>	if(Contains((PlottingSetup::ScanSampleDirectory),"SUGRA")) scantype=mSUGRA;
472	>	if(Contains((PlottingSetup::ScanSampleDirectory),"GMSB")) scantype=GMSB;
473
474		if(scantype==mSUGRA) {
475		massgluname="M0"; // this is the "x axis" in the limit plot (like the gluino in the SMS case)
#	Line 464 | Line 503 | void scan_SUSY_parameter_space(string mc
503
504		gStyle->SetPalette(100, MyPalette);
505
506	+	TH2F *rawlimitmap  = new TH2F((prefix+"rawlimitmap"+any2string(jzbSel)).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);//observed limit
507	+	TH2F *rawelimitmap  = new TH2F((prefix+"rawelimitmap"+any2string(jzbSel)).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);//observed limit
508	+
509		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
510		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
511		TH2F *exp1plimitmap  = new TH2F((prefix+"exp1plimitmap"+any2string(jzbSel)).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);//expected limit + 1 sigma
#	Line 490 | Line 532 | mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLS
532		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);
533		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);
534		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);
535	+	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);
536	+	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);
537	+
538	+	TH2F *asymptoticmap  = new TH2F((prefix+"asymptoticmap"+any2string(jzbSel)).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,(int)((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
539	+	// only used for asymptoticmap
540
541		TH2F *imposedxmap;
542		TH2F *realxmap;
#	Line 501 | Line 548 | mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLS
548
549		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);
550
551	<	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);
551	>	TH2F *nZmap    = new TH2F((prefix+"nZmap"+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);
552
553		if(ibin==0) {
554		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);
#	Line 511 | Line 558 | mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLS
558		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);
559		}
560
561	+
562	+	if(shapeanalysis) jzbSel=jzb_cut[ibin];
563		vector<int> susy_Zdecay_origin;
564		vector<TH2F*> h_susy_Zdecay_origin;
565		// we only generate these histograms IF we're dealing with the lowest cut! (we don't need them n times, they're the same for all of them)
#	Line 547 | Line 596 | write_warning(__FUNCTION__,"CURRENTLY SW
596
597		map <  pair<float, float>, map<string, float>  >  xsec;
598		if(scantype==mSUGRA) {
599	<	xsec=getXsec(PlottingSetup::mSUGRAxsFile);
599	>	xsec=getXsec(PlottingSetup::cbafbasedir+"/"+mSUGRAxsFile);
600		}
601
602		int ipoint=-1;
#	Line 557 | Line 606 | write_warning(__FUNCTION__,"CURRENTLY SW
606		ipoint++;
607		int GlobalBin=efficiencymap->FindBin(mglu,mlsp);
608		if(!runninglocally&&!do_this_point(ipoint,(int)Npoints,(int)jobnumber,(int)njobs)) continue;
609	+
610		int flipped=0;
611		float result=-987,resulterr=-987;
612		float m0trees = PlottingSetup::mgluend;
#	Line 578 | Line 628 | write_warning(__FUNCTION__,"CURRENTLY SW
628		continue;
629		} else {
630		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;
631	<	if(PlottingSetup::RestrictToMassPeak==true) cout << "Analysis type: on-peak" << endl;
632	<	if(PlottingSetup::RestrictToMassPeak==false) cout << "Analysis type: offpeak" << endl;
631	>	if(PlottingSetup::RestrictToMassPeak==true) dout << "Analysis type: on-peak" << endl;
632	>	if(PlottingSetup::RestrictToMassPeak==false) dout << "Analysis type: offpeak" << endl;
633	>	if(shapeanalysis) dout << "Shape analysis? Yes." << endl;
634	>	else dout << "Shape analysis? No." << endl;
635		}
636		float hitrecord=0;
637		int hitrecordholder=0;
#	Line 630 | Line 682 | write_warning(__FUNCTION__,"CURRENTLY SW
682		if(scantype==mSUGRA) {
683		float absxs=0;
684		for(int i=0;i<12;i++) absxs+=GetXSecForPointAndChannel(mglu,mlsp,xsec,i);
685	+	TFile *FilterEffFile = new TFile(PlottingSetup::FilterEfficiencyFile.c_str());
686	+	TH2F *FilterEfficiency = (TH2F*)FilterEffFile->Get("FilterEfficiency");
687	+	float filtereff = FilterEfficiency->GetBinContent(FilterEfficiency->FindBin(mglu,mlsp));
688	+	FilterEff->SetBinContent(GlobalBin,filtereff);
689	+	FilterEffFile->Close();
690	+	absXSmap->SetBinContent(GlobalBin,absxs);
691	+	absxs*=filtereff;
692		XSmap->SetBinContent(GlobalBin,absxs);
693		}
694
#	Line 670 | Line 729 | write_warning(__FUNCTION__,"CURRENTLY SW
729		Neventsmap->SetBinContent(GlobalBin,nevents);
730		ipointmap->SetBinContent(GlobalBin,ipoint);
731		if(efficiencyonly) continue;
732	<
733	<	do_systematics_for_one_file((scansample.collection)[scanfileindex].events,(int)nevents,"SUSY SCAN", systematics,flipped, xsec, mcjzb,datajzb,peakerror,requireZ, addcut.str(),true,scantype);
734	<	float JZBcutat = systematics[0][0];
735	<	float mceff    = systematics[0][1];
736	<	float mcefferr = systematics[0][2];//MC stat error
737	<	float toterr   = systematics[0][4];
738	<	float sys_jes  = systematics[0][5]; // Jet Energy Scale
739	<	float sys_jsu  = systematics[0][6]; // JZB scale uncertainty
740	<	float sys_res  = systematics[0][7]; // resolution
741	<	float mcwoscef = systematics[0][8]; // efficiency without signal contamination
742	<	float mcwoscefr= systematics[0][9]; // error on efficiency without signal contamination
743	<	float sys_pdf   = 0;
744	<	if(systematics[0].size()>10) sys_pdf = systematics[0][10]; // PDF
745	<	cout << "Going to store: total error is: " << toterr << " and the relative tot error, i.e. total error / mc efficiency is " << toterr/mceff << endl;
746	<	if(mceff!=mceff||toterr!=toterr||mceff<0 && (!systematicsonly&&!efficiencyonly)) {
747	<	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;
748	<	continue;
749	<	} else {
750	<	if(!systematicsonly&&!efficiencyonly) {
751	<	dout << "Calculating limit now for "<<massgluname<<"="<<mglu<<" , "<<massLSPname<<"="<<mlsp <<endl;
752	<	vector<float> sigmas;
753	<	string plotfilename=(string)(TString((scansample.collection)[scanfileindex].samplename)+TString(massgluname)+TString(any2string(mglu))+TString("__")+TString(massLSPname)+TString(any2string(mlsp))+TString(".png"));
754	<	sigmas=compute_one_upper_limit(mceff,toterr,ibin,mcjzb,plotfilename,true,flipped);
755	<	//      do_limit_wrapper(mceff,toterr,ibin,mcjzb,sigmas,plotfilename);
756	<	dout << "back in " << __FUNCTION__ << endl;
757	<	if(sigmas[0]>-0.5) { // negative sigmas are the error signature of do_limit_wrapper, so we want to exclude them.
758	<	limitmap->SetBinContent(GlobalBin,sigmas[0]);
759	<	if(sigmas.size()>1) {
760	<	explimitmap->SetBinContent(GlobalBin,sigmas[1]);
761	<	exp1plimitmap->SetBinContent(GlobalBin,sigmas[2]);
762	<	exp1mlimitmap->SetBinContent(GlobalBin,sigmas[3]);
763	<	exp2plimitmap->SetBinContent(GlobalBin,sigmas[4]);
764	<	exp2mlimitmap->SetBinContent(GlobalBin,sigmas[5]);
732	>
733	>	if(!shapeanalysis) {
734	>	do_systematics_for_one_file((scansample.collection)[scanfileindex].events,(int)nevents,"SUSY SCAN", systematics,flipped, xsec, mcjzb,datajzb,peakerror,requireZ, addcut.str(),true,scantype);
735	>	float JZBcutat = systematics[0][0];
736	>	float mceff    = systematics[0][1];
737	>	float mcefferr = systematics[0][2];//MC stat error
738	>	float toterr   = systematics[0][4];
739	>	float sys_jes  = systematics[0][5]; // Jet Energy Scale
740	>	float sys_jsu  = systematics[0][6]; // JZB scale uncertainty
741	>	float sys_res  = systematics[0][7]; // resolution
742	>	float mcwoscef = systematics[0][8]; // efficiency without signal contamination
743	>	float mcwoscefr= systematics[0][9]; // error on efficiency without signal contamination
744	>	float sys_pdf   = 0;
745	>	if(systematics[0].size()>10) sys_pdf = systematics[0][10]; // PDF
746	>	cout << "Going to store: total error is: " << toterr << " and the relative tot error, i.e. total error / mc efficiency is " << toterr/mceff << endl;
747	>	if(mceff!=mceff||toterr!=toterr||mceff<0 && (!systematicsonly&&!efficiencyonly)) {
748	>	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;
749	>	continue;
750	>	} else {
751	>	// Systematics and efficiencies make sense (i.e. non-zero, and all numbers)
752	>	if(!systematicsonly&&!efficiencyonly) {
753	>	dout << "Calculating limit now for "<<massgluname<<"="<<mglu<<" , "<<massLSPname<<"="<<mlsp <<endl;
754	>	vector<float> sigmas;
755	>	string plotfilename=(string)(TString((scansample.collection)[scanfileindex].samplename)+TString(massgluname)+TString(any2string(mglu))+TString("__")+TString(massLSPname)+TString(any2string(mlsp))+TString(".png"));
756	>	sigmas=compute_one_upper_limit(mceff,toterr,ibin,mcjzb,plotfilename,true,flipped);
757	>	//          do_limit_wrapper(mceff,toterr,ibin,mcjzb,sigmas,plotfilename);
758	>	if(sigmas[0]>-0.5) { // negative sigmas are the error signature of do_limit_wrapper, so we want to exclude them.
759	>	limitmap->SetBinContent(GlobalBin,sigmas[0]);
760	>	if(sigmas.size()>1) {
761	>	explimitmap->SetBinContent(GlobalBin,sigmas[1]);
762	>	exp1plimitmap->SetBinContent(GlobalBin,sigmas[2]);
763	>	exp1mlimitmap->SetBinContent(GlobalBin,sigmas[3]);
764	>	exp2plimitmap->SetBinContent(GlobalBin,sigmas[4]);
765	>	exp2mlimitmap->SetBinContent(GlobalBin,sigmas[5]);
766	>	}
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	>	dout << "A limit has been added at " << sigmas[0] << " for m_{glu}="<<mglu << " and m_{lsp}="<<mlsp<<endl;
774	>	}//end of if sigma is positive
775	>	}//end of not systematics only condition
776	>	if(systematicsonly) {
777	>	sysjesmap->SetBinContent(GlobalBin,sys_jes);
778	>	sysjsumap->SetBinContent(GlobalBin,sys_jsu);
779	>	sysresmap->SetBinContent(GlobalBin,sys_res);
780	>	syspdfmap->SetBinContent(GlobalBin,sys_pdf);
781	>	systotmap->SetBinContent(GlobalBin,toterr/mceff);//total relative (!) error
782	>	sysstatmap->SetBinContent(GlobalBin,mcefferr);//total relative (!) error
783	>	}
784	>	}//efficiency is valid
785	>	}
786	>	if(shapeanalysis) {
787	>	//shape analysis
788	>	cout << "This is the shape analysis. Have fun with it!" << endl;
789	>	float quickmceff,quickmcefferr;
790	>	MCefficiency((scansample.collection)[scanfileindex].events,quickmceff,quickmcefferr,flipped,mcjzb,requireZ,(int)nevents,scantype,xsec,addcut.str(),-1);
791	>	if(quickmceff==0) {
792	>	write_warning(__FUNCTION__,"The MC efficiency is zero - need to skip this point.");
793	>	continue;
794	>	}
795	>	stringstream PointName;
796	>	PointName << massgluname << "_" << mglu << "__" << massLSPname << "_" << mlsp;
797	>	SUSYScanSpace::SavedMLSP=mlsp;
798	>	SUSYScanSpace::SavedMGlu=mglu;
799	>	SUSYScanSpace::SavedMLSPname=massLSPname;
800	>	SUSYScanSpace::SavedMGluname=massgluname;
801	>	bool asymptoticonly=true; // do asymptotic computation first for everything, then (way below) if necessary switch on complete computation
802	>	ShapeDroplet droplet = LimitsFromShapes(asymptoticonly, 0, (scansample.collection)[scanfileindex].events,addcut.str().c_str(),PointName.str(),mcjzb,datajzb,jzb_cut,peakerror,peakerror);
803	>	if(!isValidShapeResult(droplet)) {
804	>	write_error(__FUNCTION__,"Something went horribly wrong with the shape analysis - skipping this point!!!");
805	>	//        scan_SUSY_parameter_space(mcjzb,datajzb,jzb_cut,requireZ, peakerror, peakerrordata, (float)Npoints, (float)ipoint, systematicsonly, efficiencyonly,false); // don't uncomment this line; it would lead to multiple limitmap histos ...
806	>	continue;
807	>	}
808	>
809	>	droplet = LimitsFromShapes(asymptoticonly, 25*droplet.expected, (scansample.collection)[scanfileindex].events,addcut.str().c_str(),PointName.str(),mcjzb,datajzb,jzb_cut,peakerror,peakerror);
810	>	if(!isValidShapeResult(droplet)) {
811	>	write_error(__FUNCTION__,"Something went horribly wrong with the shape analysis - skipping this point !!!");
812	>	//        scan_SUSY_parameter_space(mcjzb,datajzb,jzb_cut,requireZ, peakerror, peakerrordata, (float)Npoints, (float)ipoint, systematicsonly, efficiencyonly,false); // don't uncomment this line; it would lead to multiple limitmap histos ...
813	>	continue;
814		}
815
816	<	sysjesmap->SetBinContent(GlobalBin,sys_jes);
817	<	sysjsumap->SetBinContent(GlobalBin,sys_jsu);
818	<	sysresmap->SetBinContent(GlobalBin,sys_res);
819	<	syspdfmap->SetBinContent(GlobalBin,sys_pdf);
820	<	systotmap->SetBinContent(GlobalBin,toterr/mceff);//total relative (!) error
821	<	sysstatmap->SetBinContent(GlobalBin,mcefferr);//total relative (!) error
822	<	dout << "A limit has been added at " << sigmas[0] << " for m_{glu}="<<mglu << " and m_{lsp}="<<mlsp<<endl;
823	<	} //end of if sigma is positive
824	<
825	<	//end of not systematics only condition
826	<	}
827	<	if(systematicsonly) {
828	<	sysjesmap->SetBinContent(GlobalBin,sys_jes);
829	<	sysjsumap->SetBinContent(GlobalBin,sys_jsu);
830	<	sysresmap->SetBinContent(GlobalBin,sys_res);
831	<	syspdfmap->SetBinContent(GlobalBin,sys_pdf);
832	<	systotmap->SetBinContent(GlobalBin,toterr/mceff);//total relative (!) error
833	<	sysstatmap->SetBinContent(GlobalBin,mcefferr);//total relative (!) error
834	<	}
835	<	}//efficiency is valid
816	>	float referenceXS=0;
817	>	if(scantype==mSUGRA) referenceXS=XSmap->GetBinContent(GlobalBin);
818	>	else referenceXS=getSMSxs(mlsp,mglu);
819	>	float excludedXSobs = droplet.observed/droplet.SignalIntegral;
820	>	absXSmap->SetBinContent(GlobalBin,referenceXS);
821	>	write_info(__FUNCTION__,"Established limit at "+any2string(excludedXSobs)+" , ref XS is "+any2string(referenceXS));
822	>	asymptoticmap->SetBinContent(GlobalBin,1);
823	>	if(excludedXSobs>0.5*referenceXS && excludedXSobs<2.5*referenceXS||scantype!=mSUGRA) {
824	>	// for mSUGRA we only do an in-depth scan of a band around the exclusion, for SMSes we want to know the full CLs UL everywhere
825	>	write_warning(__FUNCTION__,"NOT DOING FULL CLS RIGHT NOW");
826	>	droplet = LimitsFromShapes(false, droplet.observed, (scansample.collection)[scanfileindex].events,addcut.str().c_str(),PointName.str(),mcjzb,datajzb,jzb_cut,peakerror,peakerror);
827	>	if(!isValidShapeResult(droplet)) {
828	>	write_error(__FUNCTION__,"Something went horribly wrong with the shape analysis - skipping this point !!!");
829	>	//          scan_SUSY_parameter_space(mcjzb,datajzb,jzb_cut,requireZ, peakerror, peakerrordata, (float)Npoints, (float)ipoint, systematicsonly, efficiencyonly,false); // don't uncomment this line; it would lead to multiple limitmap histos ...
830	>	continue;
831	>	}
832	>	asymptoticmap->SetBinContent(GlobalBin,0);
833	>	}
834	>	rawlimitmap->SetBinContent(GlobalBin,droplet.observed);
835	>	rawelimitmap->SetBinContent(GlobalBin,droplet.expected);
836	>	limitmap->SetBinContent(GlobalBin,droplet.observed/droplet.SignalIntegral);
837	>	explimitmap->SetBinContent(GlobalBin,droplet.expected/droplet.SignalIntegral);
838	>	exp1plimitmap->SetBinContent(GlobalBin,droplet.expectedPlus1Sigma/droplet.SignalIntegral);
839	>	exp1mlimitmap->SetBinContent(GlobalBin,droplet.expectedMinus1Sigma/droplet.SignalIntegral);
840	>	exp2plimitmap->SetBinContent(GlobalBin,droplet.expectedPlus2Sigma/droplet.SignalIntegral);
841	>	exp2mlimitmap->SetBinContent(GlobalBin,droplet.expectedMinus2Sigma/droplet.SignalIntegral);
842	>
843	>	sysjesmap->SetBinContent(GlobalBin,droplet.JES);
844	>	sysjsumap->SetBinContent(GlobalBin,droplet.JSU);
845	>	syspdfmap->SetBinContent(GlobalBin,droplet.PDF);
846	>	systotmap->SetBinContent(GlobalBin,droplet.toterr);
847	>	sysstatmap->SetBinContent(GlobalBin,droplet.staterr);
848	>	}
849		finish = clock();
850		timemap->SetBinContent(GlobalBin,((float(finish)-float(start))/CLOCKS_PER_SEC));
851		}
852		}
732	–
853		prepare_scan_axis(limitmap,scantype);
854		prepare_scan_axis(sysjesmap,scantype);
855		prepare_scan_axis(sysjsumap,scantype);
#	Line 740 | Line 860 | write_warning(__FUNCTION__,"CURRENTLY SW
860		prepare_scan_axis(timemap,scantype);
861
862		if(!systematicsonly&&!efficiencyonly) {
863	<	limcanvas->cd();
863	>	TObject* old=gDirectory->GetList()->FindObject("limcanvas");
864	>	if(old) gDirectory->GetList()->Remove(old);
865	>	TCanvas *limcanvas = new TCanvas("limcanvas","Limit canvas");
866		limitmap->Draw("COLZ");
867		string titletobewritten="Limits in LSP-Glu plane";
868		if(scantype==mSUGRA) titletobewritten="Limits in m_{1/2}-m_{0} plane";
#	Line 754 | Line 876 | write_warning(__FUNCTION__,"CURRENTLY SW
876		if(systematicsonly) outputfile=new TFile(("output/DistributedSystematics_job"+string(any2string(jobnumber))+"_of_"+string(any2string(njobs))+".root").c_str(),"UPDATE");//needs to be "UPDATE" as we can get to this point for different JZB cuts and don't want to erase the previous data :-)
877		else outputfile=new TFile(("output/DistributedLimits_job"+string(any2string(jobnumber))+"_of_"+string(any2string(njobs))+".root").c_str(),"UPDATE");//needs to be "UPDATE" as we can get to this point for
878		limitmap->Write();
879	+	rawlimitmap->Write();
880	+	rawelimitmap->Write();
881		if(doexpected) {
882		explimitmap->Write();
883		exp1plimitmap->Write();
#	Line 779 | Line 903 | write_warning(__FUNCTION__,"CURRENTLY SW
903		systotmap->Write();
904		sysstatmap->Write();
905		XSmap->Write();
906	+	absXSmap->Write();
907	+	FilterEff->Write();
908	+	if(shapeanalysis) asymptoticmap->Write();
909		if(ibin==0) imposedxmap->Write();
910		if(ibin==0) realxmap->Write();
911		Neventsmap->Write();
#	Line 841 | Line 968 | write_warning(__FUNCTION__,"CURRENTLY SW
968		syspdfmap->Write();
969		systotmap->Write();
970		sysstatmap->Write();
971	+	absXSmap->Write();
972		XSmap->Write();
973	+	FilterEff->Write();
974	+	if(shapeanalysis) asymptoticmap->Write();
975		if(ibin==0) imposedxmap->Write();
976		if(ibin==0) realxmap->Write();
977		timemap->Write();
#	Line 927 | Line 1057 | write_warning(__FUNCTION__,"CURRENTLY SW
1057		outputfile->Close();
1058		}
1059		}//end of efficiencies only
1060	<
1060	>
1061		for(int i=0;i<susy_Zdecay_origin.size();i++) {
1062		delete h_susy_Zdecay_origin[i];
1063		}
#	Line 955 | Line 1085 | write_warning(__FUNCTION__,"CURRENTLY SW
1085		delete systotmap;
1086		delete sysstatmap;
1087		delete XSmap;
1088	<	delete limcanvas;
1088	>	delete absXSmap;
1089	>	TObject* old=gDirectory->GetList()->FindObject("limcanvas");
1090	>	if(old) gDirectory->GetList()->Remove(old);
1091		}
1092
1093	<	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) {
1093	>	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) {
1094		dout << "Starting the SUSY scan now with all " << jzb_cut.size() << " bin(s)" << endl;
1095		for(int ibin=0;ibin<jzb_cut.size();ibin++) {
1096	<	scan_SUSY_parameter_space(mcjzb,datajzb,jzb_cut,requireZ, peakerror, ibin, njobs, jobnumber,systonly,effonly);
1096	>	scan_SUSY_parameter_space_for_one_cut(mcjzb,datajzb,jzb_cut,requireZ, peakerror, peakerrordata, ibin, njobs, jobnumber,systonly,effonly,shapeanalysis);
1097	>	if(shapeanalysis) scan_SUSY_parameter_space_for_one_cut(mcjzb,datajzb,jzb_cut,requireZ, peakerror, peakerrordata, ibin, njobs, jobnumber,systonly,effonly,!shapeanalysis); // if we run the shape analysis we also run cut & count - we pick the best result.
1098		}
1099		delete_any_cached_scans();// tidy up ...
1100		}

Diff Legend

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