ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/AnalysisFramework/Plotting/Modules/SUSYScan.C
Revision: 1.59
Committed: Wed Nov 16 13:41:30 2011 UTC (13 years, 5 months ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.58: +1 -0 lines
Log Message: 
Two new features: the time limit for the limit capsule algorithm is raised for the crab case (1h); if the algorithm runs into this limit three times (i.e. fails completely) a file is placed in 'exchange' which is then detected by the run script, as the algorithm counts as failed at this point, and the CRAB jobs will be considered failed

File Contents

# User Rev Content
1 buchmann 1.1 #include <iostream>
2     #include <vector>
3     #include <sys/stat.h>
4     #include <fstream>
5 buchmann 1.7 #include <pthread.h>
6 buchmann 1.51 #include <assert.h>
7 buchmann 1.1
8     #include <TCut.h>
9     #include <TROOT.h>
10     #include <TCanvas.h>
11     #include <TMath.h>
12     #include <TColor.h>
13     #include <TPaveText.h>
14     #include <TRandom.h>
15     #include <TH1.h>
16     #include <TH2.h>
17     #include <TF1.h>
18     #include <TSQLResult.h>
19     #include <TProfile.h>
20 buchmann 1.31 #include <TKey.h>
21 buchmann 1.1
22     //#include "TTbar_stuff.C"
23     using namespace std;
24    
25     using namespace PlottingSetup;
26    
27     void susy_scan_axis_labeling(TH2F *histo) {
28 buchmann 1.37 histo->GetXaxis()->SetTitle("m_{#Chi_{2}^{0}}-m_{LSP}");
29 buchmann 1.1 histo->GetXaxis()->CenterTitle();
30     histo->GetYaxis()->SetTitle("m_{#tilde{q}}");
31     histo->GetYaxis()->CenterTitle();
32     }
33    
34 buchmann 1.22 void prepare_scan_axis(TH2 *variablemap,bool ismSUGRA) {
35 buchmann 1.9 variablemap->GetXaxis()->SetTitle("m_{glu}");
36     variablemap->GetYaxis()->SetTitle("m_{LSP}");
37 buchmann 1.35
38 buchmann 1.22 if(ismSUGRA) {
39     variablemap->GetXaxis()->SetTitle("m_{0}");
40     variablemap->GetYaxis()->SetTitle("m_{1/2}");
41     }
42 buchmann 1.35
43 buchmann 1.9 variablemap->GetXaxis()->CenterTitle();
44     variablemap->GetYaxis()->CenterTitle();
45     }
46    
47 buchmann 1.31 void set_SUSY_style() {
48     Int_t MyPalette[100];
49     Double_t r[] = {0., 0.0, 1.0, 1.0, 1.0};
50     Double_t g[] = {0., 0.0, 0.0, 1.0, 1.0};
51     Double_t b[] = {0., 1.0, 0.0, 0.0, 1.0};
52     Double_t stop[] = {0., .25, .50, .75, 1.0};
53     Int_t FI = TColor::CreateGradientColorTable(5, stop, r, g, b, 110);
54     for (int i=0;i<100;i++) MyPalette[i] = FI+i;
55 buchmann 1.35
56 buchmann 1.31 gStyle->SetPalette(100, MyPalette);
57 buchmann 1.35
58 buchmann 1.31 float rightmargin=gStyle->GetPadRightMargin();
59     gStyle->SetPadRightMargin(0.15);
60    
61     }
62    
63 buchmann 1.53 bool delete_any_cached_scans() {
64 buchmann 1.50 //This should only be called via CRAB
65 buchmann 1.57 char hostname[1023];
66 buchmann 1.53 string completehostname=GetCompleteHostname();
67 buchmann 1.57 gethostname(hostname,1023);
68     if((Contains(hostname,"t3ui")||Contains(hostname,"t3wn"))) return false;
69 buchmann 1.53 else {
70     vector<string> all_files;
71     char currentpath[1024];
72     TString directory=TString(getcwd(currentpath,1024));
73     directory=directory+"/";
74     process_directory(directory,all_files);
75     for(int ifile=0;ifile<all_files.size();ifile++)
76 buchmann 1.50 {
77     if(Contains(all_files[ifile],"_clean_splitup_")) {
78     dout << "About to delete " << all_files[ifile] << endl;
79     gSystem->Exec(((string)"rm "+all_files[ifile]).c_str());
80     }
81     }
82 buchmann 1.53 return true;
83     }
84 buchmann 1.50 }
85    
86     bool initialized_t2=false;
87    
88 buchmann 1.51 int srmcpretries=0;
89    
90 buchmann 1.50 void load_scan_sample(int a, int b, int &scanfileindex,bool ismSUGRA) {
91 buchmann 1.55
92     // There is no need to define your sample here. That is done in Setup.C where you define the loading directory!
93    
94 buchmann 1.50 dout << "Going to load file with Xzone=" << a << " and Yzone=" << b << endl;
95 buchmann 1.48 stringstream filetoload;
96 buchmann 1.57 char hostname[1023];
97     string completehostname=GetCompleteHostname();
98     gethostname(hostname,1023);
99 buchmann 1.50
100 buchmann 1.57 if((Contains(hostname,"t3ui")||Contains(hostname,"t3wn"))) { // local case
101 buchmann 1.55 filetoload << "/shome/buchmann/ntuples/"<<PlottingSetup::ScanSampleDirectory;
102     if(ismSUGRA) filetoload << "/mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
103     else filetoload << "/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
104 buchmann 1.50 if(!Contains(((scansample.collection)[(scansample.collection).size()-1]).filename,"_"+any2string(a)+"_"+any2string(b))) {
105     dout << "The last sample is NOT the same one as the current one, possibly popping off last one and adding the new one." << endl;
106 buchmann 1.48 if((scansample.collection).size()>1) {
107     scansample.RemoveLastSample();
108     }
109     scanfileindex=(scansample.collection).size();
110     //New: Loading file when necessary, not before (avoiding high memory usage and startup times)
111     if(scanfileindex!=0) scansample.AddSample(filetoload.str(),"scansample",1,1,false,true,scanfileindex,kRed);
112 buchmann 1.50 } else {
113     dout << "Last sample is the same as the current one. Recycling it." << endl;
114     scanfileindex=(scansample.collection).size()-1;
115     }
116 buchmann 1.58 dout << " Going to use the following file: " << filetoload.str() << endl;
117 buchmann 1.50 } // end of t3 case
118     else {
119     write_info(__FUNCTION__,"Hello, CRAB! Might need to load files ...");
120 buchmann 1.59 PlottingSetup::limitpatience=PlottingSetup::limitpatienceCRAB;
121 buchmann 1.50 stringstream copyfile;
122 buchmann 1.57 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;
123 buchmann 1.55 if(ismSUGRA) copyfile << "/mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root " << "mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root ";
124     else copyfile << "/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root " << "SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root ";
125 buchmann 1.50 stringstream newfilename;
126     if(ismSUGRA) newfilename << "mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
127     else newfilename << "SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
128    
129     if(!Contains(((scansample.collection)[(scansample.collection).size()-1]).filename,"_"+any2string(a)+"_"+any2string(b))||(a==0&&b==0)&&!initialized_t2) {
130     initialized_t2=true;
131     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;
132     if((scansample.collection).size()>1) {
133     scansample.RemoveLastSample();
134     }
135     scanfileindex=(scansample.collection).size();
136     //New: Loading file when necessary, not before (avoiding high memory usage and startup times)
137     if(scanfileindex!=0) {
138     delete_any_cached_scans();
139     dout << "Going to download the scan file with the following copy command: " << copyfile.str() << endl;
140     gSystem->Exec(copyfile.str().c_str());
141 buchmann 1.51 if(doesROOTFileExist(newfilename.str())) {
142     scansample.AddSample(newfilename.str(),"scansample",1,1,false,true,scanfileindex,kRed);
143     } else {
144     srmcpretries++;
145     if(srmcpretries<5) {
146     dout << "The file could not be loaded correctly - retrying!" << endl;
147 buchmann 1.52 sleep(5);
148 buchmann 1.51 load_scan_sample(a,b,scanfileindex,ismSUGRA);
149     } else {
150     dout << "Have tried 5 times to load this sample. Giving up now and failing the program execution" << endl;
151     assert(0);
152     dout << "The command " << copyfile.str() << " failed to copy the file to the local storage (saving it as " << newfilename.str() << ")" << endl;
153     }
154     }
155 buchmann 1.50 }
156     } else {
157     dout << "Last sample is the same as the current one. Recycling it." << endl;
158 buchmann 1.48 scanfileindex=(scansample.collection).size()-1;
159 buchmann 1.50 }
160    
161 buchmann 1.48 }
162 buchmann 1.50 }
163 buchmann 1.48
164 buchmann 1.50 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) {
165     altxs=0;
166 buchmann 1.46 for(int iproc=1;iproc<=10;iproc++) {
167 buchmann 1.35 float process_xs = GetXSecForPointAndChannel(mglu,mlsp,xsec,iproc);
168 buchmann 1.49 altxs+=process_xs;
169 buchmann 1.46 }
170 buchmann 1.50 return altxs;
171 buchmann 1.35 }
172    
173 buchmann 1.31 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) {
174     bool runninglocally=true;
175     if(njobs>-1&&jobnumber>-1) {
176     runninglocally=false;
177     dout << "Running on the GRID (this is job " << jobnumber << "/" << njobs << ") for a jzb cut at " << jzb_cut[ibin] << endl;
178     } else {
179     dout << "Running locally " << endl;
180 buchmann 1.37 dout << "This will take a really, really long time - if you want to see the results within hours instead of weeks try running the worker script on the grid (DistributedModelCalculation/Limits/)" << endl;
181 buchmann 1.31 }
182 buchmann 1.35
183 buchmann 1.31 string massgluname="MassGlu";
184     string massLSPname="MassLSP";
185 buchmann 1.39 jzbSel=jzb_cut[ibin];
186     geqleq="geq";
187     automatized=true;
188     mcjzbexpression=mcjzb;
189    
190 buchmann 1.31 string prefix="SMS_";
191     // up to here, everything is set up for SMS; now we need to switch stuff around if we're dealing with an mSUGRA scan!
192     bool ismSUGRA=false;
193     TIter nextkey(fsyst->GetListOfKeys());
194     TKey *key;
195     while ((key = (TKey*)nextkey()))
196     {
197     TObject *obj = key->ReadObj();
198     if(Contains((string)(obj->GetName()),"mSUGRA")) ismSUGRA=true;
199     }
200 buchmann 1.35
201     map < pair<float, float>, map<string, float> > xsec;
202 buchmann 1.31
203     if(ismSUGRA) {
204     massgluname="M0"; // this is the "x axis" in the limit plot (like the gluino in the SMS case)
205     massLSPname="M12"; // this is the "y axis" in the limit plot (like the LSP in the SMS case)
206     mglustart=m0start;
207 buchmann 1.48 xsec=getXsec(PlottingSetup::cbafbasedir+"/Plotting/Modules/external/scale_xsection_nlo1.0_m0_m12_10_0_1v1.txt");
208 buchmann 1.31 mgluend=m0end;
209     mglustep=m0step;
210     mLSPstart=m12start;
211     mLSPend=m12end;
212     mLSPstep=m12step;
213     prefix="mSUGRA_";
214 buchmann 1.50 dout << "mSUGRA scan has been set up." << endl;
215 buchmann 1.31 } else {
216 buchmann 1.50 dout << "SMS scan has been set up." << endl;
217 buchmann 1.35 write_warning(__FUNCTION__,"Don't have the correct XS file yet");
218     xsec=getXsec("/scratch/buchmann/C/scale_xsection_nlo1.0_m0_m12_10_0_1v1.txt");
219 buchmann 1.31 }
220 buchmann 1.35
221 buchmann 1.31 set_SUSY_style();
222 buchmann 1.35
223 buchmann 1.31 TCanvas *limcanvas = new TCanvas("limcanvas","Limit canvas");
224 buchmann 1.35
225 buchmann 1.31 int Npoints=0;
226 buchmann 1.56 for(int mglu=(int)mglustart;mglu<=mgluend;mglu+=(int)mglustep) {
227     for (int mlsp=(int)mLSPstart;mlsp<=mLSPend&&(ismSUGRA||mlsp<=mglu);mlsp+=(int)mLSPstep) Npoints++;
228 buchmann 1.31 }
229 buchmann 1.37 TH2F *mceff = (TH2F*)fsyst->Get((prefix+"efficiencymap"+any2string(jzb_cut[ibin])).c_str());
230     //write_warning(__FUNCTION__,"Currently the efficiencymap name was switched to Pablo's convention. This NEEDS to be switched BACK!");TH2F *mceff = (TH2F*)fsyst->Get(("efficiency_jzbdiff"+any2string(jzb_cut[ibin])).c_str());
231 buchmann 1.31 TH2F *fullerr = (TH2F*)fsyst->Get((prefix+"systotmap"+any2string(jzb_cut[ibin])).c_str());
232 buchmann 1.33 TH2F *NEvents = (TH2F*)fsyst->Get((prefix+"Neventsmap"+any2string(jzb_cut[ibin])).c_str());
233 buchmann 1.52 TH2F *lflip = (TH2F*)fsyst->Get((prefix+"flipmap"+any2string(jzb_cut[ibin])).c_str());
234 buchmann 1.35
235 buchmann 1.56 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
236     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
237     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
238     TH2F *exp2plimitmap = new TH2F((prefix+"exp2plimitmap"+any2string(jzbSel)).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,(int)((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);//expected limit + 2 sigma
239     TH2F *exp1mlimitmap = new TH2F((prefix+"exp1mlimitmap"+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
240     TH2F *exp2mlimitmap = new TH2F((prefix+"exp2mlimitmap"+any2string(jzbSel)).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,(int)((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);//expected limit - 2 sigma
241 buchmann 1.37
242 buchmann 1.56 TH2F *exclmap = new TH2F((prefix+"exclusionmap"+any2string(jzb_cut[ibin])).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,(int)((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
243     TH2F *xsmap = new TH2F((prefix+"crosssectionmap"+any2string(jzb_cut[ibin])).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,(int)((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
244     TH2F *totxsmap = new TH2F((prefix+"absolutecrosssectionmap"+any2string(jzb_cut[ibin])).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,(int)((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
245     TH2F *flipmap = new TH2F((prefix+"limitflipmap"+any2string(jzb_cut[ibin])).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,(int)((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
246 buchmann 1.52
247    
248 buchmann 1.35
249 buchmann 1.32 if(!fullerr || !mceff || !NEvents) {
250 buchmann 1.31 write_error(__FUNCTION__,"The supplied systematics file did not contain the correct histograms - please check the file");
251 buchmann 1.34 dout << "mc eff address: " << mceff << " , error address: " << fullerr << " , NEvents address: " << NEvents << endl;
252 buchmann 1.31 delete limcanvas;
253     return;
254     }
255 buchmann 1.37
256     bool doexpected=true;
257 buchmann 1.35
258 buchmann 1.31 int ipoint=-1;
259 buchmann 1.56 for(int mglu=(int)mglustart;mglu<=mgluend;mglu+=(int)mglustep) {
260     for (int mlsp=(int)mLSPstart;mlsp<=mLSPend&&(ismSUGRA||mlsp<=mglu);mlsp+=(int)mLSPstep)
261 buchmann 1.31 {
262     ipoint++;
263 buchmann 1.56 if(!runninglocally&&!do_this_point(ipoint,(int)Npoints,(int)jobnumber,(int)njobs)) continue;
264 buchmann 1.31 float currmceff=mceff->GetBinContent(mceff->FindBin(mglu,mlsp));
265 buchmann 1.33 float currtoterr=(fullerr->GetBinContent(fullerr->FindBin(mglu,mlsp)))*currmceff;
266 buchmann 1.32 float nevents=NEvents->GetBinContent(NEvents->FindBin(mglu,mlsp));
267 buchmann 1.56 int flipped = (int)(lflip->GetBinContent(lflip->FindBin(mglu,mlsp)));
268 buchmann 1.52 flipmap->Fill(mglu,mlsp,flipped);
269 buchmann 1.34 dout << "Looking at point " << ipoint << " / " << Npoints << " with masses " << massgluname << " = " << mglu << " and " << massLSPname << " = " << mlsp << endl;
270 buchmann 1.35 dout << "Found : MCeff=" << currmceff << " and total error=" << currtoterr << " and Nevents=" << nevents << endl;
271 buchmann 1.31 string plotfilename=(string)(TString(massgluname)+TString(any2string(mglu))+TString("__")+TString(massLSPname)+TString(any2string(mlsp))+TString(".png"));
272 buchmann 1.34 if(currmceff<=0||currtoterr<=0||nevents==0) {
273     dout << " Nothing to work with, skipping this point." << endl;
274     continue;
275     }
276 buchmann 1.31 vector<float> sigmas;
277 buchmann 1.52 sigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped);
278 buchmann 1.48
279    
280 buchmann 1.34 if(sigmas[0]>0) limitmap->Fill(mglu,mlsp,sigmas[0]); //anything else is an error code
281 buchmann 1.37 if(sigmas.size()>1) {
282     explimitmap->Fill(mglu,mlsp,sigmas[1]);
283     exp1plimitmap->Fill(mglu,mlsp,sigmas[2]);
284     exp1mlimitmap->Fill(mglu,mlsp,sigmas[3]);
285     exp2plimitmap->Fill(mglu,mlsp,sigmas[4]);
286     exp2mlimitmap->Fill(mglu,mlsp,sigmas[5]);
287     }
288    
289 buchmann 1.35 dout << "An upper limit has been added for this point ( " << massgluname << " = " << mglu << " and " << massLSPname << " = " << mlsp << " ) at " << sigmas[0] << endl;
290 buchmann 1.39 if(ismSUGRA) { // for SMS this is a bit easier at the moment - we have a reference XS file which we use when plotting
291 buchmann 1.35 dout << "Computing exclusion status" << endl;
292     float rel_limit=0;
293 buchmann 1.49 float totxs;
294     float xs=get_xs(totxs,mglu,mlsp,massgluname,massLSPname,xsec,mcjzb,requireZ);
295 buchmann 1.35 if(xs>0) rel_limit=sigmas[0]/xs;
296     exclmap->Fill(mglu,mlsp,rel_limit);
297     xsmap->Fill(mglu,mlsp,xs);
298 buchmann 1.49 totxsmap->Fill(mglu,mlsp,totxs);
299 buchmann 1.32 }
300 buchmann 1.31 }
301     }
302 buchmann 1.35
303 buchmann 1.31 prepare_scan_axis(limitmap,ismSUGRA);
304 buchmann 1.32 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 :-)
305     outputfile->cd();
306 buchmann 1.31 limitmap->Write();
307 buchmann 1.52 flipmap->Write();
308 buchmann 1.37 if(doexpected) {
309     explimitmap->Write();
310     exp1plimitmap->Write();
311     exp1mlimitmap->Write();
312     exp2plimitmap->Write();
313     exp2mlimitmap->Write();
314     }
315 buchmann 1.32 if(ismSUGRA) {
316     exclmap->Write();
317     xsmap->Write();
318 buchmann 1.49 totxsmap->Write();
319 buchmann 1.32 }
320     outputfile->Close();
321 buchmann 1.31 delete limcanvas;
322     }
323    
324     void establish_SUSY_limits(string mcjzb,string datajzb,vector<float> jzb_cut,bool requireZ, float peakerror, string fsystfilename, float njobs=-1, float jobnumber=-1) {
325     dout << "Starting the SUSY scan from systematics file now with all " << jzb_cut.size() << " bin(s)" << " and using " << fsystfilename << endl;
326     TFile *fsyst = new TFile(fsystfilename.c_str());
327     if(!fsyst) {
328     write_error(__FUNCTION__,"You provided an invalid systematics file. Please change it ... ");
329     return;
330     }
331     for(int ibin=0;ibin<jzb_cut.size();ibin++) {
332     establish_SUSY_limits(mcjzb,datajzb,jzb_cut,requireZ, peakerror, fsyst, ibin,njobs, jobnumber);
333     }
334     }
335    
336    
337    
338    
339 buchmann 1.13 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) {
340 buchmann 1.6 bool runninglocally=true;
341     if(njobs>-1&&jobnumber>-1) {
342     runninglocally=false;
343 buchmann 1.13 dout << "Running on the GRID (this is job " << jobnumber << "/" << njobs << ") for a jzb cut at " << jzb_cut[ibin] << endl;
344 buchmann 1.6 } else {
345 buchmann 1.7 dout << "Running locally " << endl;
346 buchmann 1.40 dout << "This will take a really, really long time - if you want to see the results THIS month try running the LimitWorkerScript on the grid (DistributedModelCalculation/Limits/)" << endl;
347 buchmann 1.6 }
348 buchmann 1.35
349 buchmann 1.4 jzbSel=jzb_cut[ibin];
350     geqleq="geq";
351     automatized=true;
352     mcjzbexpression=mcjzb;
353 buchmann 1.35
354 buchmann 1.4 string massgluname="MassGlu";
355     string massLSPname="MassLSP";
356 buchmann 1.35
357 buchmann 1.23 string prefix="SMS_";
358 buchmann 1.22 // up to here, everything is set up for SMS; now we need to switch stuff around if we're dealing with an mSUGRA scan!
359     bool ismSUGRA=false;
360 buchmann 1.24 if(Contains((scansample.collection)[0].samplename,"mSUGRA")) ismSUGRA=true;
361 buchmann 1.22 if(ismSUGRA) {
362     massgluname="M0"; // this is the "x axis" in the limit plot (like the gluino in the SMS case)
363     massLSPname="M12"; // this is the "y axis" in the limit plot (like the LSP in the SMS case)
364     mglustart=m0start;
365     mgluend=m0end;
366     mglustep=m0step;
367     mLSPstart=m12start;
368     mLSPend=m12end;
369     mLSPstep=m12step;
370 buchmann 1.23 prefix="mSUGRA_";
371 buchmann 1.50 dout << "mSUGRA scan has been set up." << endl;
372 buchmann 1.24 } else {
373 buchmann 1.50 dout << "SMS scan has been set up." << endl;
374 buchmann 1.24 }
375 buchmann 1.35
376 buchmann 1.4 Int_t MyPalette[100];
377     Double_t r[] = {0., 0.0, 1.0, 1.0, 1.0};
378     Double_t g[] = {0., 0.0, 0.0, 1.0, 1.0};
379     Double_t b[] = {0., 1.0, 0.0, 0.0, 1.0};
380     Double_t stop[] = {0., .25, .50, .75, 1.0};
381     Int_t FI = TColor::CreateGradientColorTable(5, stop, r, g, b, 110);
382     for (int i=0;i<100;i++) MyPalette[i] = FI+i;
383 buchmann 1.35
384 buchmann 1.4 gStyle->SetPalette(100, MyPalette);
385 buchmann 1.35
386 buchmann 1.56 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
387     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
388     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
389     TH2F *exp2plimitmap = new TH2F((prefix+"exp2plimitmap"+any2string(jzbSel)).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);//expected limit + 2 sigma
390     TH2F *exp1mlimitmap = new TH2F((prefix+"exp1mlimitmap"+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
391     TH2F *exp2mlimitmap = new TH2F((prefix+"exp2mlimitmap"+any2string(jzbSel)).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,int((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);//expected limit - 2 sigma
392 buchmann 1.37
393 buchmann 1.56 TH2F *exclmap = new TH2F((prefix+"exclusionmap"+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);
394     TH2F *sysjesmap = new TH2F((prefix+"sysjes"+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);
395     TH2F *sysjsumap = new TH2F((prefix+"sysjsu"+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);
396     TH2F *sysresmap = new TH2F((prefix+"sysresmap"+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);
397     TH2F *efficiencymap = new TH2F((prefix+"efficiencymap"+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);
398     TH2F *efficiencymapAcc = new TH2F((prefix+"efficiencymapAcc"+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);
399     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);
400     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);
401     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);
402     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);
403     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);
404     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);
405     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);
406     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);
407     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);
408 buchmann 1.27
409 buchmann 1.56 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);
410 buchmann 1.36
411 buchmann 1.56 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);
412 buchmann 1.52
413 buchmann 1.29 write_warning(__FUNCTION__,"CURRENTLY SWITCHING AUTOMATIZED MODE OFF!");automatized=false;
414 buchmann 1.35
415 buchmann 1.4 float rightmargin=gStyle->GetPadRightMargin();
416     gStyle->SetPadRightMargin(0.15);
417 buchmann 1.35
418 buchmann 1.4 TCanvas *limcanvas = new TCanvas("limcanvas","Limit canvas");
419 buchmann 1.35
420 buchmann 1.37
421     bool doexpected=true;
422    
423 buchmann 1.6 int Npoints=0;
424 buchmann 1.56 for(int mglu=(int)mglustart;mglu<=mgluend;mglu+=(int)mglustep) {
425     for (int mlsp=(int)mLSPstart;mlsp<=mLSPend&&(ismSUGRA||mlsp<=mglu);mlsp+=(int)mLSPstep) Npoints++;
426 buchmann 1.6 }
427 buchmann 1.35
428 buchmann 1.6 int ipoint=-1;
429 buchmann 1.56 for(int mglu=(int)mglustart;mglu<=mgluend;mglu+=(int)mglustep) {
430     for (int mlsp=(int)mLSPstart;mlsp<=mLSPend&&(ismSUGRA||mlsp<=mglu);mlsp+=(int)mLSPstep)
431 buchmann 1.4 {
432 buchmann 1.6 ipoint++;
433 buchmann 1.56 if(!runninglocally&&!do_this_point(ipoint,(int)Npoints,(int)jobnumber,(int)njobs)) continue;
434 buchmann 1.52 int flipped=0;
435 buchmann 1.19 float result=-987,resulterr=-987;
436 buchmann 1.50 float m0trees = PlottingSetup::mgluend;
437     float m12trees = PlottingSetup::mLSPend;
438 buchmann 1.42 int a = int((PlottingSetup::ScanXzones*mlsp)/(m12trees+1));
439     int b = int((PlottingSetup::ScanYzones*mglu)/(m0trees+1));
440 buchmann 1.50 int scanfileindex=PlottingSetup::ScanYzones*a+b;
441     load_scan_sample(a,b,scanfileindex,ismSUGRA);
442 buchmann 1.42
443 buchmann 1.36 clock_t start,finish;
444     start = clock(); // starting the clock to measure how long the computation takes!
445 buchmann 1.4 stringstream addcut;
446     addcut << "(TMath::Abs("<<massgluname<<"-"<<mglu<<")<5)&&(TMath::Abs("<<massLSPname<<"-"<<mlsp<<")<5)";
447 buchmann 1.40 (scansample.collection)[scanfileindex].events->Draw("eventNum",addcut.str().c_str(),"goff");
448     float nevents = (scansample.collection)[scanfileindex].events->GetSelectedRows();
449 buchmann 1.4 vector<vector<float> > systematics;
450 buchmann 1.30 if(nevents<10) {
451 buchmann 1.40 dout << "This point ("<<ipoint<<") with configuration ("<<massgluname<<"="<<mglu<<" , "<<massLSPname<<"="<<mlsp << ") does not contain enough events and will be skipped."<< endl;
452     continue;
453 buchmann 1.6 } else {
454 buchmann 1.48 dout << "OK! This point ("<<ipoint<<") with configuration ("<<massgluname<<"="<<mglu<<" , "<<massLSPname<<"="<<mlsp << ") contains " << nevents << " and will therefore not be skipped. " << endl;
455     write_analysis_type(PlottingSetup::RestrictToMassPeak);
456 buchmann 1.22 }
457 buchmann 1.50 if(nevents!=0&&(efficiencyonly||systematicsonly)) {
458 buchmann 1.56 Value effwosigcont = MCefficiency((scansample.collection)[scanfileindex].events,result,resulterr,flipped,mcjzb,requireZ,(int)nevents,addcut.str(),-1);
459 buchmann 1.52 if(result<0) {
460     write_info(__FUNCTION__,"Found negative efficiency (signal distribution populating JZB<0 more than JZB>0) - flipping analysis!");
461     flipped=1;
462 buchmann 1.56 MCefficiency((scansample.collection)[scanfileindex].events,result,resulterr,flipped,mcjzb,requireZ,(int)nevents,addcut.str(),-1);
463 buchmann 1.52 }
464 buchmann 1.36 efficiencymap->Fill(mglu,mlsp,result);
465 buchmann 1.47 efficiencymap->SetBinError(mglu,mlsp,resulterr);
466 buchmann 1.52 flipmap->Fill(mglu,mlsp,flipped);
467 buchmann 1.37 noscefficiencymap->Fill(mglu,mlsp,effwosigcont.getValue());
468     noscefficiencymap->SetBinError(mglu,mlsp,effwosigcont.getError());
469 buchmann 1.47 //Partial efficiencies
470     float resultAcc, resultID, resultJets, resultSignal;
471     float resulterrAcc, resulterrID, resulterrJets, resulterrSignal;
472 buchmann 1.56 MCPartialefficiency((scansample.collection)[scanfileindex].events,resultAcc,resulterrAcc,flipped,mcjzb,requireZ,(int)nevents,addcut.str(),-1, 1);
473     MCPartialefficiency((scansample.collection)[scanfileindex].events,resultID,resulterrID,flipped,mcjzb,requireZ,(int)nevents,addcut.str(),-1, 2);
474     MCPartialefficiency((scansample.collection)[scanfileindex].events,resultJets,resulterrJets,flipped,mcjzb,requireZ,(int)nevents,addcut.str(),-1, 3);
475     MCPartialefficiency((scansample.collection)[scanfileindex].events,resultSignal,resulterrSignal,flipped,mcjzb,requireZ,(int)nevents,addcut.str(),-1, 4);
476 buchmann 1.47 efficiencymapAcc->Fill(mglu,mlsp,resultAcc);
477     efficiencymapAcc->SetBinError(mglu,mlsp,resulterrAcc);
478     efficiencymapID->Fill(mglu,mlsp,resultID);
479     efficiencymapAcc->SetBinError(mglu,mlsp,resulterrID);
480     efficiencymapJets->Fill(mglu,mlsp,resultJets);
481     efficiencymapAcc->SetBinError(mglu,mlsp,resulterrJets);
482     efficiencymapSignal->Fill(mglu, mlsp, resultSignal);
483     efficiencymapAcc->SetBinError(mglu,mlsp,resulterrSignal);
484 buchmann 1.36 finish = clock();
485 buchmann 1.13 }
486 buchmann 1.17 Neventsmap->Fill(mglu,mlsp,nevents);
487 buchmann 1.24 ipointmap->Fill(mglu,mlsp,ipoint);
488 buchmann 1.17 if(efficiencyonly) continue;
489    
490 buchmann 1.56 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
491 buchmann 1.9 float JZBcutat = systematics[0][0];
492     float mceff = systematics[0][1];
493 buchmann 1.22 float mcefferr = systematics[0][2];//MC stat error
494 buchmann 1.9 float toterr = systematics[0][4];
495     float sys_jes = systematics[0][5]; // Jet Energy Scale
496     float sys_jsu = systematics[0][6]; // JZB scale uncertainty
497     float sys_res = systematics[0][7]; // resolution
498 buchmann 1.37 float mcwoscef = systematics[0][8]; // efficiency without signal contamination
499     float mcwoscefr= systematics[0][9]; // error on efficiency without signal contamination
500 buchmann 1.27 float sys_pdf = 0;
501 buchmann 1.37 if(systematics[0].size()>10) sys_pdf = systematics[0][10]; // PDF
502 buchmann 1.52 // efficiencymap->Fill(mglu,mlsp,mceff);
503     // efficiencymap->SetBinError(efficiencymap->FindBin(mglu,mlsp),mcefferr);
504 buchmann 1.37 noscefficiencymap->Fill(mglu,mlsp,mcwoscef);
505     noscefficiencymap->SetBinError(efficiencymap->FindBin(mglu,mlsp),mcwoscefr);
506 buchmann 1.35
507 buchmann 1.47 if(mceff!=mceff||toterr!=toterr||mceff<0 && (!systematicsonly&&!efficiencyonly)) {
508 buchmann 1.35 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;
509     continue;
510 buchmann 1.5 } else {
511 buchmann 1.35 if(!systematicsonly&&!efficiencyonly) {
512     dout << "Calculating limit now for "<<massgluname<<"="<<mglu<<" , "<<massLSPname<<"="<<mlsp <<endl;
513     vector<float> sigmas;
514 buchmann 1.40 string plotfilename=(string)(TString((scansample.collection)[scanfileindex].samplename)+TString(massgluname)+TString(any2string(mglu))+TString("__")+TString(massLSPname)+TString(any2string(mlsp))+TString(".png"));
515 buchmann 1.52 sigmas=compute_one_upper_limit(mceff,toterr,ibin,mcjzb,plotfilename,true,flipped);
516 buchmann 1.46 // do_limit_wrapper(mceff,toterr,ibin,mcjzb,sigmas,plotfilename);
517 buchmann 1.50 dout << "back in " << __FUNCTION__ << endl;
518 buchmann 1.35 if(sigmas[0]>-0.5) { // negative sigmas are the error signature of do_limit_wrapper, so we want to exclude them.
519     limitmap->Fill(mglu,mlsp,sigmas[0]);
520 buchmann 1.37 if(sigmas.size()>1) {
521     explimitmap->Fill(mglu,mlsp,sigmas[1]);
522     exp1plimitmap->Fill(mglu,mlsp,sigmas[2]);
523     exp1mlimitmap->Fill(mglu,mlsp,sigmas[3]);
524     exp2plimitmap->Fill(mglu,mlsp,sigmas[4]);
525     exp2mlimitmap->Fill(mglu,mlsp,sigmas[5]);
526 buchmann 1.38 }
527 buchmann 1.37
528 buchmann 1.35 sysjesmap->Fill(mglu,mlsp,sys_jes);
529     sysjsumap->Fill(mglu,mlsp,sys_jsu);
530     sysresmap->Fill(mglu,mlsp,sys_res);
531     syspdfmap->Fill(mglu,mlsp,sys_pdf);
532     systotmap->Fill(mglu,mlsp,toterr/mceff);//total relative (!) error
533     sysstatmap->Fill(mglu,mlsp,mcefferr);//total relative (!) error
534     dout << "A limit has been added at " << sigmas[0] << " for m_{glu}="<<mglu << " and m_{lsp}="<<mlsp<<endl;
535     } //end of if sigma is positive
536 buchmann 1.36
537 buchmann 1.35 //end of not systematics only condition
538     }
539     if(systematicsonly) {
540     sysjesmap->Fill(mglu,mlsp,sys_jes);
541     sysjsumap->Fill(mglu,mlsp,sys_jsu);
542     sysresmap->Fill(mglu,mlsp,sys_res);
543     syspdfmap->Fill(mglu,mlsp,sys_pdf);
544     systotmap->Fill(mglu,mlsp,toterr/mceff);//total relative (!) error
545     sysstatmap->Fill(mglu,mlsp,mcefferr);//total relative (!) error
546     }
547 buchmann 1.15 }//efficiency is valid
548 buchmann 1.52 finish = clock();
549     timemap->Fill(mglu,mlsp,((float(finish)-float(start))/CLOCKS_PER_SEC));
550 buchmann 1.4 }
551     }
552 buchmann 1.52
553 buchmann 1.22 prepare_scan_axis(limitmap,ismSUGRA);
554     prepare_scan_axis(sysjesmap,ismSUGRA);
555     prepare_scan_axis(sysjsumap,ismSUGRA);
556     prepare_scan_axis(sysresmap,ismSUGRA);
557 buchmann 1.27 prepare_scan_axis(syspdfmap,ismSUGRA);
558     prepare_scan_axis(systotmap,ismSUGRA);
559     prepare_scan_axis(sysstatmap,ismSUGRA);
560 buchmann 1.36 prepare_scan_axis(timemap,ismSUGRA);
561 buchmann 1.35
562 buchmann 1.13 if(!systematicsonly&&!efficiencyonly) {
563 buchmann 1.9 limcanvas->cd();
564     limitmap->Draw("COLZ");
565 buchmann 1.22 string titletobewritten="Limits in LSP-Glu plane";
566     if(ismSUGRA) titletobewritten="Limits in m_{1/2}-m_{0} plane";
567     TText *title = write_title(titletobewritten);
568 buchmann 1.9 title->Draw();
569     if(runninglocally) {
570     CompleteSave(limcanvas,"SUSYScan/Limits_JZB_geq"+any2string(jzb_cut[ibin]));
571     } else {
572 buchmann 1.12 TFile *outputfile;
573     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 :-)
574 buchmann 1.35 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
575 buchmann 1.9 limitmap->Write();
576 buchmann 1.37 if(doexpected) {
577     explimitmap->Write();
578     exp1plimitmap->Write();
579     exp1mlimitmap->Write();
580     exp2plimitmap->Write();
581     exp2mlimitmap->Write();
582     }
583 buchmann 1.18 sysjesmap->Write();
584     sysjsumap->Write();
585 buchmann 1.28 sysresmap->Write();
586 buchmann 1.18 efficiencymap->Write();
587 buchmann 1.52 flipmap->Write();
588 buchmann 1.47 efficiencymapAcc->Write();
589     efficiencymapID->Write();
590     efficiencymapJets->Write();
591     efficiencymapSignal->Write();
592 buchmann 1.37 noscefficiencymap->Write();
593 buchmann 1.27 syspdfmap->Write();
594     systotmap->Write();
595     sysstatmap->Write();
596 buchmann 1.18 Neventsmap->Write();
597 buchmann 1.24 ipointmap->Write();
598 buchmann 1.36 timemap->Write();
599 buchmann 1.9 outputfile->Close();
600     }
601 buchmann 1.35 }
602 buchmann 1.13 if(systematicsonly) { // systematics only :
603 buchmann 1.9 limcanvas->cd();
604     sysjesmap->Draw("COLZ");
605 buchmann 1.22 string titletobewritten="Jet Energy scale in LSP-Glu plane";
606     if(ismSUGRA) titletobewritten="Limits in m_{1/2}-m_{0} plane";
607     TText *title = write_title(titletobewritten);
608 buchmann 1.9 title->Draw();
609     if(runninglocally) {
610     CompleteSave(limcanvas,"SUSYScan/JES_geq"+any2string(jzb_cut[ibin]));
611     }
612     sysjsumap->Draw("COLZ");
613 buchmann 1.22 titletobewritten="JZB Scale Uncertainty in LSP-Glu plane";
614     if(ismSUGRA) titletobewritten="JZB Scale Uncertainty in m_{1/2}-m_{0} plane";
615     TText *title2 = write_title(titletobewritten);
616 buchmann 1.9 title2->Draw();
617     if(runninglocally) {
618     CompleteSave(limcanvas,"SUSYScan/JSU_geq"+any2string(jzb_cut[ibin]));
619     }
620     sysresmap->Draw("COLZ");
621 buchmann 1.22 titletobewritten="Resolution in LSP-Glu plane";
622     if(ismSUGRA) titletobewritten="Resolution in m_{1/2}-m_{0} plane";
623     TText *title3 = write_title(titletobewritten);
624 buchmann 1.9 title3->Draw();
625     if(runninglocally) {
626     CompleteSave(limcanvas,"SUSYScan/Resolution_geq"+any2string(jzb_cut[ibin]));
627     }
628 buchmann 1.35
629 buchmann 1.9 if(!runninglocally) {
630     TFile *outputfile=new TFile(("output/DistributedSystematics_job"+string(any2string(jobnumber))+"_of_"+string(any2string(njobs))+".root").c_str(),"UPDATE");
631     sysjesmap->Write();
632     sysjsumap->Write();
633     sysresmap->Write();
634 buchmann 1.52 flipmap->Write();
635 buchmann 1.18 efficiencymap->Write();
636 buchmann 1.47 efficiencymapAcc->Write();
637     efficiencymapID->Write();
638     efficiencymapJets->Write();
639     efficiencymapSignal->Write();
640 buchmann 1.37 noscefficiencymap->Write();
641 buchmann 1.18 Neventsmap->Write();
642 buchmann 1.24 ipointmap->Write();
643 buchmann 1.27 syspdfmap->Write();
644     systotmap->Write();
645     sysstatmap->Write();
646 buchmann 1.36 timemap->Write();
647 buchmann 1.9 outputfile->Close();
648     }
649     }//end of systematics only
650 buchmann 1.13 if(efficiencyonly) {
651     limcanvas->cd();
652 buchmann 1.22 string titletobewritten="Efficiencies in LSP-Glu plane";
653     if(ismSUGRA) titletobewritten="Efficiencies in m_{1/2}-m_{0} plane";
654     TText *title = write_title(titletobewritten);
655 buchmann 1.47 efficiencymap->Draw("COLZ");
656 buchmann 1.13 title->Draw();
657     if(runninglocally) {
658     CompleteSave(limcanvas,"SUSYScan/Efficiency_geq"+any2string(jzb_cut[ibin]));
659     }
660     limcanvas->cd();
661 buchmann 1.47 efficiencymapAcc->Draw("COLZ");
662     title->Draw();
663     if(runninglocally) {
664     CompleteSave(limcanvas,"SUSYScan/EfficiencyAcc_geq"+any2string(jzb_cut[ibin]));
665     }
666     limcanvas->cd();
667     efficiencymapID->Draw("COLZ");
668     title->Draw();
669     if(runninglocally) {
670     CompleteSave(limcanvas,"SUSYScan/EfficiencyID_geq"+any2string(jzb_cut[ibin]));
671     }
672     limcanvas->cd();
673     efficiencymapJets->Draw("COLZ");
674     title->Draw();
675     if(runninglocally) {
676     CompleteSave(limcanvas,"SUSYScan/EfficiencyJets_geq"+any2string(jzb_cut[ibin]));
677     }
678     limcanvas->cd();
679     efficiencymapSignal->Draw("COLZ");
680     title->Draw();
681     if(runninglocally) {
682     CompleteSave(limcanvas,"SUSYScan/EfficiencySignal_geq"+any2string(jzb_cut[ibin]));
683     }
684     limcanvas->cd();
685     noscefficiencymap->Draw("COLZ");
686     title->Draw();
687     if(runninglocally) {
688     CompleteSave(limcanvas,"SUSYScan/NoEfficiency_geq"+any2string(jzb_cut[ibin]));
689     }
690     limcanvas->cd();
691 buchmann 1.13 sysjesmap->Draw("COLZ");
692 buchmann 1.22 titletobewritten="N(events) in LSP-Glu plane";
693     if(ismSUGRA) titletobewritten="N(events) in m_{1/2}-m_{0} plane";
694     TText *title2 = write_title(titletobewritten);
695 buchmann 1.13 title2->Draw();
696     if(runninglocally) {
697     CompleteSave(limcanvas,"SUSYScan/Nevents_geq"+any2string(jzb_cut[ibin]));
698     }
699     if(!runninglocally) {
700     TFile *outputfile=new TFile(("output/DistributedSystematics_job"+string(any2string(jobnumber))+"_of_"+string(any2string(njobs))+".root").c_str(),"UPDATE");
701 buchmann 1.24 ipointmap->Write();
702 buchmann 1.13 Neventsmap->Write();
703 buchmann 1.37 noscefficiencymap->Write();
704 buchmann 1.13 efficiencymap->Write();
705 buchmann 1.52 flipmap->Write();
706 buchmann 1.47 efficiencymapAcc->Write();
707     efficiencymapID->Write();
708     efficiencymapJets->Write();
709     efficiencymapSignal->Write();
710 buchmann 1.36 timemap->Write();
711 buchmann 1.13 outputfile->Close();
712     }
713     }//end of efficiencies only
714 buchmann 1.35
715 buchmann 1.31 delete limitmap;
716     delete exclmap;
717     delete sysjesmap;
718     delete sysjsumap;
719     delete sysresmap;
720 buchmann 1.37 delete noscefficiencymap;
721 buchmann 1.31 delete efficiencymap;
722 buchmann 1.47 delete efficiencymapAcc;
723     delete efficiencymapID;
724     delete efficiencymapJets;
725     delete efficiencymapSignal;
726 buchmann 1.31 delete Neventsmap;
727     delete ipointmap;
728     delete syspdfmap;
729     delete systotmap;
730     delete sysstatmap;
731 buchmann 1.13 delete limcanvas;
732 buchmann 1.1 }
733 buchmann 1.4
734 buchmann 1.13 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) {
735 buchmann 1.10 dout << "Starting the SUSY scan now with all " << jzb_cut.size() << " bin(s)" << endl;
736 buchmann 1.4 for(int ibin=0;ibin<jzb_cut.size();ibin++) {
737 buchmann 1.13 scan_SUSY_parameter_space(mcjzb,datajzb,jzb_cut,requireZ, peakerror, ibin, njobs, jobnumber,systonly,effonly);
738 buchmann 1.4 }
739 buchmann 1.6 }
740 buchmann 1.35