ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/ShapeLimit.C

Comparing UserCode/cbrown/Development/Plotting/Modules/ShapeLimit.C (file contents):
Revision 1.10 by buchmann, Mon May 7 16:44:27 2012 UTC vs.
Revision 1.13 by buchmann, Thu May 17 19:52:55 2012 UTC

#	Line 18 | Line 18
18		#include <TSystem.h>
19		#include <TRandom3.h>
20
21	–	//#include "TTbar_stuff.C"
21		using namespace std;
22
23		using namespace PlottingSetup;
#	Line 44 | Line 43 | void EliminateNegativeEntries(TH1F *hist
43		}
44		}
45
46	<	void prepare_limit_datacard(string RunDirectory, TFile *f, float uJSU, float uPDF) {
47	<	TH1F *dataob = (TH1F*)f->Get("data_obs");
48	<	TH1F *signal = (TH1F*)f->Get("signal");
49	<	TH1F *Tbackground = (TH1F*)f->Get("TTbarBackground");
50	<	TH1F *Zbackground  = (TH1F*)f->Get("ZJetsBackground");
51	<
52	<	assert(dataob);
53	<	assert(signal);
54	<	assert(Tbackground);
55	<	assert(Zbackground);
56	<
46	>	vector<float> get_expected_points(float observed,int npoints, string RunDirectory, bool doPlot=false) {
47	>	TF1 *gaus = new TF1("gaus","gaus(0)",0,10*observed);
48	>	gaus->SetParameters(1,observed,2*observed);
49	>	gaus->SetParameter(0,1/(gaus->Integral(0,10*observed)));
50	>	float currentpoint=0.01;
51	>	float stepsize=observed/100;
52	>	vector<float> points;
53	>	while(currentpoint<25*observed && points.size()<npoints) {
54	>
55	>	if(gaus->Integral(0,currentpoint)>((points.size()+1.0)/(npoints+2.0))) {
56	>	if(Verbosity>0) dout << "Added points for calculation at " << currentpoint << " (integral is " << gaus->Integral(0,currentpoint) << ")" << endl;
57	>	points.push_back(currentpoint);
58	>	if(points.size()==npoints) break;
59	>	}
60	>	currentpoint+=stepsize;
61	>	}
62	>	if(doPlot) {
63	>	gaus->SetName("Expected limit computation points");
64	>	gaus->SetTitle("Expected limit computation points");
65	>	TCanvas *can = new TCanvas("can","can");
66	>	gaus->Draw();
67	>	TLine *lines[npoints];
68	>	for(int i=0;i<npoints;i++) {
69	>	lines[i] = new TLine(points[i],0,points[i],gaus->GetMaximum());
70	>	lines[i]->SetLineColor(kBlue);
71	>	lines[i]->Draw();
72	>	}
73	>	can->SaveAs((RunDirectory+"/DistributionOfComputedPoints.png").c_str());
74	>	delete can;
75	>	for(int i=0;i<npoints;i++) delete lines[i];
76	>	}
77	>	delete gaus;
78	>
79	>	return points;
80	>	}
81	>
82	>	void prepare_MC_datacard(string RunDirectory, TFile *f, float uJSU, float uPDF) {
83	>	TH1F *dataob = (TH1F*)f->Get("data_obs");
84	>	TH1F *signal = (TH1F*)f->Get("signal");
85	>	assert(dataob);
86	>	assert(signal);
87	>
88	>	write_warning(__FUNCTION__,"The following two defs for th1's are fake");TH1F *Tbackground; TH1F *Zbackground;
89	>	write_warning(__FUNCTION__,"Not correctly implemented yet for MC");
90	>	ofstream datacard;
91	>	write_warning(__FUNCTION__,"Need to rethink systematics we want to use !");
92	>	datacard.open ((RunDirectory+"/susydatacard.txt").c_str());
93	>	dout << "Writing this to a file.\n";
94	>	datacard << "imax 1\n"; // number of channels
95	>	datacard << "jmax 2\n"; // number of backgrounds (Z+Jets prediction and TTbar prediction)
96	>	datacard << "kmax *\n"; // number of nuisance parameters (sources of systematic uncertainties)
97	>	datacard << "---------------\n";
98	>	datacard << "shapes * * limitfile.root $PROCESS $PROCESS_$SYSTEMATIC\n";
99	>	datacard << "---------------\n";
100	>	datacard << "bin 1\n";
101	>	datacard << "observation "<<dataob->Integral()<<"\n";
102	>	datacard << "------------------------------\n";
103	>	datacard << "bin             1          1          1\n";
104	>	datacard << "process         signal     TTbarBackground     ZJetsBackground\n";
105	>	datacard << "process         0          1          2\n";
106	>	datacard << "rate            "<<signal->Integral()<<"         "<<Tbackground->Integral()<<"         "<<Zbackground->Integral()<<"\n";
107	>	datacard << "--------------------------------\n";
108	>	datacard << "lumi     lnN    " << 1+ PlottingSetup::lumiuncert << "       -       -    luminosity uncertainty\n"; // only affects MC -> signal!
109	>	datacard << "Stat   shape    -          1          -    statistical uncertainty (ttbar)\n";
110	>	datacard << "Stat   shape    -          -          1    statistical uncertainty (zjets)\n";
111	>	datacard << "Sys    shape    -          1          -    systematic uncertainty on ttbar\n";
112	>	datacard << "Sys    shape    -          -          1    systematic uncertainty on zjets\n";
113	>	datacard << "Stat   shape    1          -          -    statistical uncertainty (signal)\n";
114	>	datacard << "JES    shape    1          -          -    uncertainty on Jet Energy Scale\n";
115	>	datacard << "JSU      lnN    " << 1+uJSU << "          -          -    JZB Scale Uncertainty\n";
116	>	if(uPDF>0) datacard << "PDF      lnN    " << 1+uPDF << "          -          -    uncertainty from PDFs\n";
117	>	datacard.close();
118	>	}
119	>
120	>	void prepare_datadriven_datacard(string RunDirectory, TFile *f, float uJSU, float uPDF) {
121	>	TH1F *dataob = (TH1F*)f->Get("data_obs");
122	>	TH1F *signal = (TH1F*)f->Get("signal");
123	>	TH1F *Tbackground = (TH1F*)f->Get("TTbarBackground");
124	>	TH1F *Zbackground  = (TH1F*)f->Get("ZJetsBackground");
125	>
126	>	assert(dataob);
127	>	assert(signal);
128	>	assert(Tbackground);
129	>	assert(Zbackground);
130
131	<	if(FullMCAnalysis) {
60	<	//dealing with MC analysis - need to use shapes.
61	<	write_warning(__FUNCTION__,"Not correctly implemented yet for MC");
62	<	ofstream datacard;
63	<	write_warning(__FUNCTION__,"Need to rethink systematics we want to use !");
64	<	datacard.open ((RunDirectory+"/susydatacard.txt").c_str());
65	<	datacard << "Writing this to a file.\n";
66	<	datacard << "imax 1\n"; // number of channels
67	<	datacard << "jmax 2\n"; // number of backgrounds (Z+Jets prediction and TTbar prediction)
68	<	datacard << "kmax *\n"; // number of nuisance parameters (sources of systematic uncertainties)
69	<	datacard << "---------------\n";
70	<	datacard << "shapes * * limitfile.root $PROCESS $PROCESS_$SYSTEMATIC\n";
71	<	datacard << "---------------\n";
72	<	datacard << "bin 1\n";
73	<	datacard << "observation "<<dataob->Integral()<<"\n";
74	<	datacard << "------------------------------\n";
75	<	datacard << "bin             1          1          1\n";
76	<	datacard << "process         signal     TTbarBackground     ZJetsBackground\n";
77	<	datacard << "process         0          1          2\n";
78	<	datacard << "rate            "<<signal->Integral()<<"         "<<Tbackground->Integral()<<"         "<<Zbackground->Integral()<<"\n";
79	<	datacard << "--------------------------------\n";
80	<	datacard << "lumi     lnN    " << 1+ PlottingSetup::lumiuncert << "       -       -    luminosity uncertainty\n"; // only affects MC -> signal!
81	<	datacard << "Stat   shape    -          1          -    statistical uncertainty (ttbar)\n";
82	<	datacard << "Stat   shape    -          -          1    statistical uncertainty (zjets)\n";
83	<	datacard << "Sys    shape    -          1          -    systematic uncertainty on ttbar\n";
84	<	datacard << "Sys    shape    -          -          1    systematic uncertainty on zjets\n";
85	<	datacard << "Stat   shape    1          -          -    statistical uncertainty (signal)\n";
86	<	datacard << "JES    shape    1          -          -    uncertainty on Jet Energy Scale\n";
87	<	datacard << "JSU      lnN    " << 1+uJSU << "          -          -    JZB Scale Uncertainty\n";
88	<	if(uPDF>0) datacard << "PDF      lnN    " << 1+uPDF << "          -          -    uncertainty from PDFs\n";
89	<	datacard.close();
90	<	} else {
91	<	//doing mutibin analysis
131	>
132		ofstream datacard;
133		write_warning(__FUNCTION__,"Need to rethink systematics we want to use !");
134		datacard.open ((RunDirectory+"/susydatacard.txt").c_str());
#	Line 212 | Line 252 | void prepare_limit_datacard(string RunDi
252		}
253
254		datacard.close();
255	+	}
256	+
257	+	void prepare_limit_datacard(string RunDirectory, TFile *f, float uJSU, float uPDF) {
258	+
259	+	if(FullMCAnalysis) {
260	+	//dealing with MC analysis - need to use shapes.
261	+	prepare_MC_datacard(RunDirectory,f,uJSU,uPDF);
262	+	} else {
263	+	//doing mutibin analysis
264	+	prepare_datadriven_datacard(RunDirectory,f,uJSU,uPDF);
265		}
266
267		}
#	Line 261 | Line 311 | TH1F* Multiply(TH1F *h1, TH1F *h2) {
311		return h;
312		}
313
314	+
315	+	string ReduceNumericHistoName(string origname) {
316	+	int pos = origname.find("__h_");
317	+	if(pos == -1) return origname;
318	+	return origname.substr(0,pos);
319	+	}
320	+
321	+
322	+
323	+	void generate_mc_shapes(bool dotree, TTree *signalevents, string mcjzb,string datajzb,vector<float> binning, TCut weight, TCut JetCut, string addcut, string identifier, map <  pair<float, float>, map<string, float>  > &xsec) {
324	+	//weight can be varied to take care of efficiency errors (weightEffUp, weightEffDown)
325	+	vector<float> mbinning;
326	+	mbinning.push_back(-8000);
327	+	for(int i=binning.size()-1;i>=0;i--) mbinning.push_back(-binning[i]);
328	+	mbinning.push_back(0);
329	+	for(int i=0;i<binning.size();i++) mbinning.push_back(binning[i]);
330	+	mbinning.push_back(8000);
331	+
332	+
333	+	TCanvas *shapecan = new TCanvas("shapecan","shapecan");
334	+	TH1F* Signal;
335	+
336	+	stringstream sInverseCutWeight;
337	+	sInverseCutWeight << "(1.0/" << (const char*)cutWeight<<")";
338	+	TCut InverseCutWeight(sInverseCutWeight.str().c_str());
339	+	if(weight==cutWeight) InverseCutWeight = TCut("1.0");
340	+	if(!dotree) {
341	+	//dealing with ALL MC samples (when we save the standard file)
342	+	THStack McPredicted = allsamples.DrawStack("McPred",mcjzb,mbinning, "JZB [GeV]", "events", ((cutmass&&cutOSSF&&JetCut)*(weight*InverseCutWeight)),mc,luminosity);
343	+	int nhists = McPredicted.GetHists()->GetEntries();
344	+	for (int i = 0; i < nhists; i++) {
345	+	TH1* hist = static_cast<TH1*>(McPredicted.GetHists()->At(i));
346	+	//          cout << hist->GetName() << " : " << hist->Integral() << endl;
347	+	//          cout << "     " << ReduceNumericHistoName(hist->GetName()) << endl;
348	+	if(identifier=="StatUp") {
349	+	float appliedweight = hist->Integral() / hist->GetEntries();
350	+	for(int ibin=1;ibin<=hist->GetNbinsX();ibin++) hist->SetBinContent(ibin,hist->GetBinContent(ibin)+appliedweight*TMath::Sqrt(hist->GetBinContent(ibin)/appliedweight));
351	+	}
352	+	if(identifier=="StatDown") {
353	+	float appliedweight = hist->Integral() / hist->GetEntries();
354	+	for(int ibin=1;ibin<=hist->GetNbinsX();ibin++) hist->SetBinContent(ibin,hist->GetBinContent(ibin)-appliedweight*TMath::Sqrt(hist->GetBinContent(ibin)/appliedweight));
355	+	}
356	+	hist->SetName(("mc_"+ReduceNumericHistoName(hist->GetName())+"_"+identifier).c_str());
357	+	hist->Write();
358	+	}
359	+	} else {
360	+	string histoname="mc_signal";
361	+	if(identifier!="") histoname=histoname+"_"+identifier;
362	+	Signal = QuickDraw(signalevents,histoname,mcjzb,binning,"JZB", "events",((cutmass&&cutOSSF&&JetCut&&basiccut))*(weight*InverseCutWeight),addcut,mc,luminosity,xsec);
363	+	if(identifier=="StatUp") {
364	+	float appliedweight = Signal->Integral() / Signal->GetEntries();
365	+	for(int ibin=1;ibin<=Signal->GetNbinsX();ibin++) Signal->SetBinContent(ibin,Signal->GetBinContent(ibin)+appliedweight*TMath::Sqrt(Signal->GetBinContent(ibin)/appliedweight));
366	+	}
367	+	if(identifier=="StatDown") {
368	+	float appliedweight = Signal->Integral() / Signal->GetEntries();
369	+	for(int ibin=1;ibin<=Signal->GetNbinsX();ibin++) Signal->SetBinContent(ibin,Signal->GetBinContent(ibin)-appliedweight*TMath::Sqrt(Signal->GetBinContent(ibin)/appliedweight));
370	+	}
371	+	Signal->Write();
372	+	}
373	+
374	+	/*
375	+	*
376	+	* Jet energy scale (easy, use different JetCut)
377	+	* JZB Scale Uncertainty (will be a number - signal only)
378	+	* Efficiency (will be weightEffUp, weightEffDown)
379	+	* PDFs (signal only) -> Will be a number yet again, computed in the traditional way
380	+	*/
381	+	delete shapecan;
382	+
383	+	}
384	+
385	+
386		void generate_shapes_for_systematic(bool signalonly, bool dataonly,TFile *limfile, TTree *signalevents, string identifier, string mcjzb, string datajzb, int JES,vector<float> binning, TCanvas *limcan, string addcut, map <  pair<float, float>, map<string, float>  > &xsec) {
387
388		binning.push_back(8000);
#	Line 824 | Line 946 | void generate_shapes_for_systematic(bool
946
947		}
948
949	<	ShapeDroplet LimitsFromShapes(bool asymptotic, float firstGuess, TTree *events, string addcut, string name, string mcjzb, string datajzb, vector<float> jzbbins, float jzbpeakerrormc, float jzbpeakerrordata) {
949	>	void DoFullCls(string RunDirectory,float firstGuess) {
950	>	vector<float> epoints = get_expected_points(firstGuess,25,RunDirectory,true);//get 20 points;
951	>	ofstream RealScript;
952	>	dout << "Going to write the full CLs script to " << RunDirectory << "/LimitScript.sh" << endl;
953	>	RealScript.open ((RunDirectory+"/LimitScript.sh").c_str());
954	>	RealScript << "#!/bin/bash \n";
955	>	RealScript << "\n";
956	>	RealScript << "RunDirectory=" << RunDirectory << "\n";
957	>	RealScript << "CMSSWDirectory=" << PlottingSetup::CMSSW_dir << "\n";
958	>	RealScript << "ORIGTMP=$TMP\n";
959	>	RealScript << "ORIGTMPDIR=$TMPDIR\n";
960	>	RealScript << "export TMP=" << RunDirectory << "\n";
961	>	RealScript << "export TMPDIR=" << RunDirectory << "\n";
962	>	RealScript << "\n";
963	>	RealScript << "echo \"Going to compute limit in $RunDirectory\"\n";
964	>	RealScript << "ORIGDIR=`pwd`\n";
965	>	RealScript << "\n";
966	>	RealScript << "pushd $CMSSWDirectory\n";
967	>	RealScript << "cd ~/final_production_2011/CMSSW_4_2_8/src/HiggsAnalysis/\n";
968	>	RealScript << "origscramarch=$SCRAM_ARCH\n";
969	>	RealScript << "origbase=$CMSSW_BASE\n";
970	>	RealScript << "export SCRAM_ARCH=slc5_amd64_gcc434\n";
971	>	RealScript << "eval `scram ru -sh`\n";
972	>	RealScript << "\n";
973	>	RealScript << "mkdir -p $RunDirectory\n";
974	>	RealScript << "cd $RunDirectory\n";
975	>	RealScript << "echo `pwd`\n";
976	>	RealScript << "mkdir -p FullCLs\n";
977	>	RealScript << "cd FullCLs\n";
978	>	RealScript << "\n";
979	>	RealScript << "combineEXE=\"${CMSSW_BASE}/bin/${SCRAM_ARCH}/combine\"\n";
980	>	RealScript << "\n";
981	>	RealScript << "dobreak=0\n";
982	>	RealScript << "npoints=0\n";
983	>	RealScript << "time for i in {";
984	>	RealScript << epoints[0]/4 << ","; // adding a VERY VERY low point just to be totally safe
985	>	RealScript << epoints[0]/2 << ","; // adding a VERY low point just to be safe
986	>	for(int i=0;i<epoints.size();i++) RealScript << epoints[i] << ",";
987	>	RealScript << 2*epoints[epoints.size()-1] << ","; // adding a VERY high point just to be safe
988	>	RealScript << 4*epoints[epoints.size()-1]; // adding a VERY VERY high point just to be totally safe
989	>	RealScript << "}; do \n";
990	>	RealScript << "let \"npoints = npoints +1\"\n";
991	>	RealScript << "if [[ $dobreak -gt 0 ]]; then \n";
992	>	RealScript << "    continue \n";
993	>	RealScript << "fi \n";
994	>	RealScript << "echo -e \" Going to compute CLs for $i \\n\"\n";
995	>	RealScript << "echo $(date +%d%b%Y,%H:%M:%S)\n";
996	>	RealScript << "time " << PlottingSetup::cbafbasedir << "/DistributedModelCalculations/Utilities/TimeOut \"$combineEXE ../susydatacard.txt -M HybridNew --freq -T 500 --clsAcc 0 -v 0 -n TestPoint --saveHybridResult --saveToys -s -1 -i 8 --singlePoint $i >/dev/null 2>/dev/null\"\n";
997	>	//  RealScript << "$combineEXE ../susydatacard.txt -M HybridNew --freq -T 500 --clsAcc 0 -v 0 -n TestPoint --saveHybridResult --saveToys -s -1 -i 8 --singlePoint $i >/dev/null 2>/dev/null \n";
998	>	RealScript << "errorsize=$?\n";
999	>	RealScript << "rm " << RunDirectory << "/rstat* 2>/dev/null \n";
1000	>	RealScript << "if [[ $errorsize -gt 0 ]]; then \n";
1001	>	RealScript << "    echo \"Apparently something went wrong (had to be aborted)\"\n";
1002	>	RealScript << "    if [[ $npoints -gt 10 ]]; then \n";
1003	>	RealScript << "        dobreak=1 \n";
1004	>	RealScript << "    fi \n";
1005	>	RealScript << "fi \n";
1006	>	RealScript << "done\n";
1007	>
1008	>	RealScript << "\n";
1009	>	RealScript << "hadd -f FullCLs.root higgsCombine*.root\n";
1010	>	RealScript << "mkdir originalFiles\n";
1011	>	RealScript << "mv higgsCombine* originalFiles\n";
1012	>	RealScript << "$combineEXE ../susydatacard.txt -M HybridNew --freq --grid=FullCLs.root\n";
1013	>	RealScript << "$combineEXE ../susydatacard.txt -M HybridNew --freq --grid=FullCLs.root --expectedFromGrid 0.5\n";
1014	>	RealScript << "$combineEXE ../susydatacard.txt -M HybridNew --freq --grid=FullCLs.root --expectedFromGrid 0.84\n";
1015	>	RealScript << "$combineEXE ../susydatacard.txt -M HybridNew --freq --grid=FullCLs.root --expectedFromGrid 0.16\n";
1016	>	RealScript << "$combineEXE ../susydatacard.txt -M HybridNew --freq --grid=FullCLs.root --expectedFromGrid 0.975\n";
1017	>	RealScript << "$combineEXE ../susydatacard.txt -M HybridNew --freq --grid=FullCLs.root --expectedFromGrid 0.025\n";
1018	>	RealScript << "\n";
1019	>	RealScript << "hadd FinalResult.root higgsCombineTest.HybridNew*\n";
1020	>	RealScript << "\n";
1021	>	RealScript << "g++ " << PlottingSetup::cbafbasedir << "/DistributedModelCalculations/ShapeLimits/ReadAndSave.C -o ReadAndSave.exec `root-config --glibs --cflags` \n";
1022	>	RealScript << "./ReadAndSave.exec FinalResult.root " << RunDirectory << "/FullCLsShapeDropletResult.txt \n";
1023	>	RealScript << "\n";
1024	>	RealScript << "rm " << RunDirectory << "/rstat* \n";
1025	>	RealScript << "\n";
1026	>	RealScript << "#####\n";
1027	>	RealScript << "echo \"Finalizing ...\"\n";
1028	>	RealScript << "cd $ORIGDIR\n";
1029	>	RealScript << "echo $ORIGDIR\n";
1030	>	RealScript << "ls -ltrh ../SandBox/ | grep combine\n";
1031	>	RealScript << "export SCRAM_ARCH=$origscramarch\n";
1032	>	RealScript << "export TMP=$ORIGTMP\n";
1033	>	RealScript << "export TMPDIR=$ORIGTMPDIR\n";
1034	>	RealScript << "exit 0\n";
1035	>	RealScript << "\n";
1036	>	RealScript.close();
1037	>	gSystem->Exec(("bash "+RunDirectory+"/LimitScript.sh").c_str());
1038	>	}
1039	>
1040	>
1041	>	ShapeDroplet LimitsFromShapes(float low_fullCLs, float high_CLs, TTree *events, string addcut, string name, string mcjzb, string datajzb, vector<float> jzbbins, float jzbpeakerrormc, float jzbpeakerrordata) {
1042		map <  pair<float, float>, map<string, float>  >  xsec;
1043		if(SUSYScanSpace::SUSYscantype==mSUGRA) xsec=getXsec(PlottingSetup::cbafbasedir+"/"+mSUGRAxsFile);
1044
#	Line 929 | Line 1143 | ShapeDroplet LimitsFromShapes(bool asymp
1143		limfile->Close();
1144		dout << "Info: Shape root file and datacard have been generated in " << RunDirectory << endl;
1145		stringstream command;
1146	<	if(asymptotic) {
1147	<	if(firstGuess>0) command << "bash " << PlottingSetup::cbafbasedir<< "/DistributedModelCalculations/ShapeLimits/CreateModel.sh " << RunDirectory << " susydatacard.txt" << " 1 0 0 " << firstGuess; // ASYMPTOTIC LIMITS WITH FIRST GUESS
1148	<	else command << "bash " << PlottingSetup::cbafbasedir<< "/DistributedModelCalculations/ShapeLimits/CreateModel.sh " << RunDirectory << " susydatacard.txt" << " 0 0 0 0"; // ASYMPTOTIC LIMITS
1149	<	}
1150	<	else command << "bash " << PlottingSetup::cbafbasedir<< "/DistributedModelCalculations/ShapeLimits/CreateModel.sh " << RunDirectory << " susydatacard.txt" << " 2 0 " << int(0.333 * firstGuess) << " " << int(firstGuess); // ASYMPTOTIC LIMITS
1146	>	string DropletLocation;
1147	>	int CreatedModelFileExitCode;
1148	>	//----------------------------------------------------------------
1149	>	//do an asymptotic limit first
1150	>	command << "bash " << PlottingSetup::cbafbasedir<< "/DistributedModelCalculations/ShapeLimits/CreateModel.sh " << RunDirectory << " susydatacard.txt" << " 0 0 0 0"; // ASYMPTOTIC LIMITS
1151		dout <<"Going to run : " << command.str() << endl;
1152	<	int CreatedModelFileExitCode = gSystem->Exec(command.str().c_str());
1152	>	CreatedModelFileExitCode = gSystem->Exec(command.str().c_str());
1153		dout << "exit code of limit algorithm (CreateModel.sh) : " << CreatedModelFileExitCode << endl;
1154		if(!(CreatedModelFileExitCode==0)) {
1155	<	write_warning(__FUNCTION__,"Something bad happened. It looks like a shape analysis is not the way to go. ");
1156	<	ShapeDroplet alpha;
1157	<	alpha.observed=-12345; // this is the flag to say watch out something went wrong with the signal ...
1158	<	alpha.SignalIntegral=1;
1159	<	return alpha;
1155	>	write_warning(__FUNCTION__,"Something bad happened. It looks like a shape analysis is not the way to go. ");
1156	>	ShapeDroplet beta;
1157	>	beta.observed=-12345; // this is the flag to say watch out something went wrong with the signal ...
1158	>	beta.SignalIntegral=1;
1159	>	return beta;
1160	>	}
1161	>	DropletLocation=RunDirectory+"/ShapeDropletResult.txt";
1162	>	ShapeDroplet beta;
1163	>	beta.readDroplet(DropletLocation);
1164	>	float obsLimit = beta.observed/SignalIntegral;
1165	>	dout << "Checking if the obtained limit (" << beta.observed << " / " << SignalIntegral << " = " << beta.observed/SignalIntegral << " is in [" << low_fullCLs << " , " << high_CLs << "]" << endl;
1166	>	if(obsLimit<high_CLs && obsLimit>low_fullCLs) {
1167	>	//if(PlottingSetup::scantype!=mSUGRA||(obsLimit<high_CLs && obsLimit>low_fullCLs)) {
1168	>	dout << " It is! Full CLs ahead!" << endl;
1169	>	DoFullCls(RunDirectory,beta.observed);
1170	>	DropletLocation=RunDirectory+"/FullCLsShapeDropletResult.txt";
1171	>	} else {
1172	>	dout << " It is not ... happy with asymptotic limits." << endl;
1173		}
1174	+	//----------------------------------------------------------------
1175		ShapeDroplet alpha;
1176	<	alpha.readDroplet(RunDirectory+"/ShapeDropletResult.txt");
1176	>	alpha.readDroplet(DropletLocation);
1177		alpha.PDF=PDFuncert;
1178		alpha.JSU=JZBscale;
1179		alpha.SignalIntegral=SignalIntegral;
1180		dout << "Done reading limit information from droplet" << endl;
1181		dout << alpha << endl;
1182	<
1182	>
1183		dout << "Everything is saved in " << RunDirectory << endl;
1184		dout << "Cleaning up ... " << std::flush;
1185	<	/*  dout << "   1) Make sure models directory exists ... " << std::flush;
1186	<	gSystem->Exec("mkdir -p models/");
959	<	dout << " ok!" << endl;
960	<	dout << "   2) Deleting any previous model files with the same name ... " << std::flush;
961	<	stringstream delcommand;
962	<	delcommand << "rm models/model_" << name << ".root";
963	<	gSystem->Exec(delcommand.str().c_str());
964	<	stringstream delcommand2;
965	<	delcommand2 << "rm models/model_" << name << ".txt";
966	<	gSystem->Exec(delcommand2.str().c_str());
967	<	dout << " ok!" << endl;
968	<	dout << "   3) Transfering the new model files ... " << std::flush;
969	<	stringstream copycommand;
970	<	copycommand << "cp " << RunDirectory << "/model.root models/model_" << name << ".root";
971	<	gSystem->Exec(copycommand.str().c_str());
972	<	stringstream copycommand2;
973	<	copycommand2 << "cp " << RunDirectory << "/susydatacard.txt models/model_" << name << ".txt";
974	<	gSystem->Exec(copycommand2.str().c_str());
975	<	stringstream copycommand3;
976	<	copycommand3 << "cp " << RunDirectory << "/limitfile.root models/model_" << name << "_histo.root";
977	<	gSystem->Exec(copycommand3.str().c_str());
978	<	dout << " ok!" << endl;
979	<	dout << "   4) Removing original working directory (" << RunDirectory << ") ... " << std::flush;*/
980	<	gSystem->Exec(("rm -r "+RunDirectory).c_str());
1185	>	write_warning(__FUNCTION__,"NOT CLEANING UP");
1186	>	gSystem->Exec(("rm -rf "+RunDirectory).c_str());
1187		dout << " ok!" << endl;
1188		delete limcan;
1189		return alpha;
1190		}
1191
1192	<	void PrepareDataShapes(string datajzb, vector<float> jzbbins, float jzbpeakerrordata) {
1193	<	//  dout << "Generating data shape templates ... " << endl;
1192	>	void PrepareDataShapes(string mcjzb, string datajzb, vector<float> jzbbins, float jzbpeakerrordata) {
1193	>
1194		map <  pair<float, float>, map<string, float>  >  xsec;
1195		TFile *datafile = new TFile("../DistributedModelCalculations/StoredShapes.root","RECREATE");
1196		TCanvas *limcan = new TCanvas("limcan","Canvas for calculating limits");
1197		TTree *faketree;
1198		bool dataonly=true;
1199		bool signalonly=false;
1200	<	string mcjzb="JZBforMCinPrepareDataShapes";//this string is not used.
995	<	//  float jzbpeakerrormc=0;
1200	>
1201		generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"",mcjzb,datajzb,noJES,jzbbins,limcan,"",xsec);
1202	<	// don't need these effects for obs & pred, only for signal!
1203	<	//  generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"peakUp",newjzbexpression(mcjzb,jzbpeakerrormc),newjzbexpression(datajzb,jzbpeakerrordata),noJES,jzbbins,limcan,"",xsec);
1204	<	//  generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"peakDown",newjzbexpression(mcjzb,-jzbpeakerrormc),newjzbexpression(datajzb,-jzbpeakerrordata),noJES,jzbbins,limcan,"",xsec);
1205	<	//  generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"JESUp",mcjzb,datajzb,JESup,jzbbins,limcan,"",xsec);
1206	<	//  generate_shapes_for_systematic(signalonly,dataonly,datafile,faketree,"JESDown",mcjzb,datajzb,JESdown,jzbbins,limcan,"",xsec);
1202	>
1203	>	generate_mc_shapes(false,faketree,mcjzb,datajzb,jzbbins,cutWeight,cutnJets,"","",xsec);
1204	>	generate_mc_shapes(false,faketree,mcjzb,datajzb,jzbbins,cutWeight,cutnJetsJESup,"","JESUp",xsec);
1205	>	generate_mc_shapes(false,faketree,mcjzb,datajzb,jzbbins,cutWeight,cutnJetsJESdown,"","JESDown",xsec);
1206	>	generate_mc_shapes(false,faketree,mcjzb,datajzb,jzbbins,TCut("weightEffUp"),cutnJets,"","EffUp",xsec);
1207	>	generate_mc_shapes(false,faketree,mcjzb,datajzb,jzbbins,TCut("weightEffDown"),cutnJets,"","EffDown",xsec);
1208	>	generate_mc_shapes(false,faketree,mcjzb,datajzb,jzbbins,cutWeight,cutnJets,"","StatUp",xsec);
1209	>	generate_mc_shapes(false,faketree,mcjzb,datajzb,jzbbins,cutWeight,cutnJets,"","StatDown",xsec);
1210	>	/*
1211		datafile->Close();
1212	+	*/
1213		}
1214

Diff Legend

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