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.3 by buchmann, Thu Mar 15 20:34:15 2012 UTC

#	Line 95 | Line 95 | bool initialized_t2=false;
95		int srmcpretries=0;
96
97		void load_scan_sample(int a, int b, int &scanfileindex,int scantype,bool isretry=false) {
98	<
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;
#	Line 108 | Line 108 | void load_scan_sample(int a, int b, int
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))) {
111	>	if(scansample.collection.size()<1||!Contains(((scansample.collection)[(scansample.collection).size()-1]).filename,"_"+any2string(a)+"_"+any2string(b))) {
112		dout << "The last sample is NOT the same one as the current one, possibly popping off last one and adding the new one." << endl;
113		if((scansample.collection).size()>1) {
114		scansample.RemoveLastSample();
#	Line 168 | Line 168 | void load_scan_sample(int a, int b, int
168		}
169
170		}
171	+	*/
172	+
173	+	// WATCH OUT THIS LINE NEEDS TO BE REMOVED
174	+	scanfileindex=0;
175		}
176
177	<	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)  {
177	>	/*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)  {
178		altxs=0;
179	<	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	<	}
179	>	for(int iproc=1;iproc<=10;iproc++) altxs+=GetXSecForPointAndChannel(mglu,mlsp,xsec,iproc);
180		return altxs;
181	<	}
181	>	}*/
182
183		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) {
184		bool runninglocally=true;
#	Line 251 | Line 250 | void establish_SUSY_limits(string mcjzb,
250		TH2F *fullerr   = (TH2F*)fsyst->Get((prefix+"systotmap"+any2string(jzb_cut[ibin])).c_str());
251		TH2F *NEvents   = (TH2F*)fsyst->Get((prefix+"Neventsmap"+any2string(jzb_cut[ibin])).c_str());
252		TH2F *lflip     = (TH2F*)fsyst->Get((prefix+"flipmap"+any2string(jzb_cut[ibin])).c_str());
253	+	TH2F *XSmap     = (TH2F*)fsyst->Get((prefix+"XS"+any2string(jzb_cut[ibin])).c_str());
254
255		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
256		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 294 | void establish_SUSY_limits(string mcjzb,
294		continue;
295		}
296		vector<float> sigmas;
297	<	sigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped);
297	>	if(scantype!=mSUGRA) sigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped);
298	>	else {
299	>	//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"!
300	>	vector<float> asigmas;
301	>	asigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped,true); // asymptotic limit first
302	>	float strength=asigmas[0]/XSmap->GetBinContent(GlobalBin);
303	>	if(strength>0.5&&strength<2) sigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped); // asymptotic limit first
304	>	else sigmas=asigmas;
305	>	exclmap->SetBinContent(GlobalBin,strength);
306	>	xsmap->SetBinContent(GlobalBin,XSmap->GetBinContent(GlobalBin));
307	>	}
308	>
309	>
310
311
312		if(sigmas[0]>0) limitmap->SetBinContent(GlobalBin,sigmas[0]); //anything else is an error code
#	Line 307 | Line 319 | void establish_SUSY_limits(string mcjzb,
319		}
320
321		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	–	}
322		}
323		}
324
#	Line 430 | Line 432 | void scan_SUSY_parameter_space(string mc
432		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);
433		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);
434		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);
435	<	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);
435	>	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-
436	>	mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
437		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);
438		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);
439		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);
440		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);
441		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);
442		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);
443	+	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);
444
445		TH2F *imposedxmap;
446		TH2F *realxmap;
447
448	+	TH2F *centermap;
449	+	TH2F *widthmap;
450	+
451		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);
452
453		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 455 | void scan_SUSY_parameter_space(string mc
455		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);
456
457		if(ibin==0) {
458	+	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);
459	+	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);
460	+
461		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);
462		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);
463		}
#	Line 485 | Line 495 | write_warning(__FUNCTION__,"CURRENTLY SW
495		for(int mglu=(int)mglustart;mglu<=mgluend;mglu+=(int)mglustep) {
496		for (int mlsp=(int)mLSPstart;mlsp<=mLSPend&&(scantype==mSUGRA||mlsp<=mglu);mlsp+=(int)mLSPstep) Npoints++;
497		}
498	<
498	>
499	>	map <  pair<float, float>, map<string, float>  >  xsec;
500	>	if(scantype==mSUGRA) {
501	>	xsec=getXsec(PlottingSetup::mSUGRAxsFile);
502	>	}
503	>
504		int ipoint=-1;
505		for(int mglu=(int)mglustart;mglu<=mgluend;mglu+=(int)mglustep) {
506		for (int mlsp=(int)mLSPstart;mlsp<=mLSPend&&(scantype==mSUGRA||mlsp<=mglu);mlsp+=(int)mLSPstep)
#	Line 510 | Line 525 | write_warning(__FUNCTION__,"CURRENTLY SW
525		float nevents = (scansample.collection)[scanfileindex].events->GetSelectedRows();
526		vector<vector<float> > systematics;
527		if(nevents<10) {
528	<	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;
528	>	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;
529		continue;
530		} else {
531		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 567 | write_warning(__FUNCTION__,"CURRENTLY SW
567		(scansample.collection)[scanfileindex].events->Draw("imposedx>>realxhisto",addcut.str().c_str(),"goff");
568		imposedxmap->SetBinContent(GlobalBin,realxhisto->GetMean());
569		imposedxmap->SetBinError(GlobalBin,realxhisto->GetRMS());
570	+	stringstream specialexpression;
571	+	specialexpression << mcjzb << ">>jzbhisto";
572	+	TH1F *jzbhisto = new TH1F("jzbhisto","jzbhisto",100,-500,500);
573	+	(scansample.collection)[scanfileindex].events->Draw(specialexpression.str().c_str(),addcut.str().c_str(),"goff");
574	+	centermap->SetBinContent(GlobalBin,jzbhisto->GetMean());
575	+	widthmap->SetBinContent(GlobalBin,jzbhisto->GetRMS());
576		delete realxhisto;
577	+	delete jzbhisto;
578		}
579
580
581		if(nevents!=0&&(efficiencyonly||systematicsonly)) {
582	<	Value effwosigcont = MCefficiency((scansample.collection)[scanfileindex].events,result,resulterr,flipped,mcjzb,requireZ,(int)nevents,addcut.str(),-1);
582	>	Value effwosigcont = MCefficiency((scansample.collection)[scanfileindex].events,result,resulterr,flipped,mcjzb,requireZ,(int)nevents,scantype,xsec,addcut.str(),-1);
583		if(result<0&&allowflipping) {
584		write_info(__FUNCTION__,"Found negative efficiency (signal distribution populating JZB<0 more than JZB>0) - flipping analysis!");
585		flipped=1;
586	<	MCefficiency((scansample.collection)[scanfileindex].events,result,resulterr,flipped,mcjzb,requireZ,(int)nevents,addcut.str(),-1);
586	>	MCefficiency((scansample.collection)[scanfileindex].events,result,resulterr,flipped,mcjzb,requireZ,(int)nevents,scantype,xsec,addcut.str(),-1);
587		}
588		efficiencymap->SetBinContent(GlobalBin,result);
589		efficiencymap->SetBinError(GlobalBin,resulterr);
#	Line 594 | Line 616 | write_warning(__FUNCTION__,"CURRENTLY SW
616		ipointmap->SetBinContent(GlobalBin,ipoint);
617		if(efficiencyonly) continue;
618
619	<	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
619	>	do_systematics_for_one_file((scansample.collection)[scanfileindex].events,(int)nevents,"SUSY SCAN", systematics,flipped, xsec, mcjzb,datajzb,peakerror,requireZ, addcut.str(),true,scantype);
620		float JZBcutat = systematics[0][0];
621		float mceff    = systematics[0][1];
622		float mcefferr = systematics[0][2];//MC stat error
#	Line 688 | Line 710 | write_warning(__FUNCTION__,"CURRENTLY SW
710		sysjsumap->Write();
711		sysresmap->Write();
712		efficiencymap->Write();
713	+	if(ibin==0) centermap->Write();
714	+	if(ibin==0) widthmap->Write();
715		flipmap->Write();
716		efficiencymapmet->Write();
717		efficiencymapAcc->Write();
#	Line 699 | Line 723 | write_warning(__FUNCTION__,"CURRENTLY SW
723		syspdfmap->Write();
724		systotmap->Write();
725		sysstatmap->Write();
726	+	XSmap->Write();
727		if(ibin==0) imposedxmap->Write();
728		if(ibin==0) realxmap->Write();
729		Neventsmap->Write();
730		ipointmap->Write();
731		timemap->Write();
707	–	outputfile->Close();
732		for(int i=0;i<susy_Zdecay_origin.size();i++) {
733		h_susy_Zdecay_origin[i]->Write();
734		}
735	+	outputfile->Close();
736		}
737		}
738		if(systematicsonly) { // systematics only :
#	Line 746 | Line 771 | write_warning(__FUNCTION__,"CURRENTLY SW
771		sysjsumap->Write();
772		sysresmap->Write();
773		flipmap->Write();
774	+	if(ibin==0) centermap->Write();
775	+	if(ibin==0) widthmap->Write();
776		efficiencymap->Write();
777		efficiencymapmet->Write();
778		efficiencymapAcc->Write();
#	Line 759 | Line 786 | write_warning(__FUNCTION__,"CURRENTLY SW
786		syspdfmap->Write();
787		systotmap->Write();
788		sysstatmap->Write();
789	+	XSmap->Write();
790		if(ibin==0) imposedxmap->Write();
791		if(ibin==0) realxmap->Write();
792		timemap->Write();
765	–	outputfile->Close();
793		for(int i=0;i<susy_Zdecay_origin.size();i++) {
794		h_susy_Zdecay_origin[i]->Write();
795		}
796	+	outputfile->Close();
797		}
798		}//end of systematics only
799		if(efficiencyonly) {
#	Line 832 | Line 860 | write_warning(__FUNCTION__,"CURRENTLY SW
860		noscefficiencymap->Write();
861		efficiencymapContam->Write();
862		efficiencymap->Write();
863	+	if(ibin==0) centermap->Write();
864	+	if(ibin==0) widthmap->Write();
865		flipmap->Write();
866		efficiencymapmet->Write();
867		efficiencymapAcc->Write();
#	Line 869 | Line 899 | write_warning(__FUNCTION__,"CURRENTLY SW
899		delete syspdfmap;
900		delete systotmap;
901		delete sysstatmap;
902	+	delete XSmap;
903		delete limcanvas;
904		}
905

Diff Legend

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