ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/SUSYScan.C
Revision: 1.11
Committed: Wed Apr 18 09:22:22 2012 UTC (13 years ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.10: +6 -0 lines
Log Message:
Added a check for shapes to see if the efficiency is non-zero

File Contents

# User Rev Content
1 buchmann 1.1 #include <iostream>
2     #include <vector>
3     #include <sys/stat.h>
4     #include <fstream>
5     #include <pthread.h>
6     #include <assert.h>
7    
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     #include <TKey.h>
21    
22 buchmann 1.9 #include "ShapeDroplet.C"
23     #include "ShapeLimit.C"
24    
25 buchmann 1.1 //#include "TTbar_stuff.C"
26     using namespace std;
27    
28     using namespace PlottingSetup;
29    
30     void susy_scan_axis_labeling(TH2F *histo) {
31     histo->GetXaxis()->SetTitle("m_{#Chi_{2}^{0}}-m_{LSP}");
32     histo->GetXaxis()->CenterTitle();
33     histo->GetYaxis()->SetTitle("m_{#tilde{q}}");
34     histo->GetYaxis()->CenterTitle();
35     }
36    
37     void prepare_scan_axis(TH2 *variablemap,int scantype) {
38     //assuming regular SMS
39     variablemap->GetXaxis()->SetTitle("m_{glu}");
40     variablemap->GetYaxis()->SetTitle("m_{LSP}");
41    
42     if(scantype==mSUGRA) {
43     variablemap->GetXaxis()->SetTitle("m_{0}");
44     variablemap->GetYaxis()->SetTitle("m_{1/2}");
45     }
46    
47     if(scantype==GMSB) {
48     variablemap->GetXaxis()->SetTitle("m_{glu}");
49     variablemap->GetYaxis()->SetTitle("m_{NLSP}");
50     }
51    
52     variablemap->GetXaxis()->CenterTitle();
53     variablemap->GetYaxis()->CenterTitle();
54     }
55    
56     void set_SUSY_style() {
57     Int_t MyPalette[100];
58     Double_t r[] = {0., 0.0, 1.0, 1.0, 1.0};
59     Double_t g[] = {0., 0.0, 0.0, 1.0, 1.0};
60     Double_t b[] = {0., 1.0, 0.0, 0.0, 1.0};
61     Double_t stop[] = {0., .25, .50, .75, 1.0};
62     Int_t FI = TColor::CreateGradientColorTable(5, stop, r, g, b, 110);
63     for (int i=0;i<100;i++) MyPalette[i] = FI+i;
64    
65     gStyle->SetPalette(100, MyPalette);
66    
67     float rightmargin=gStyle->GetPadRightMargin();
68     gStyle->SetPadRightMargin(0.15);
69    
70     }
71    
72     bool delete_any_cached_scans() {
73 buchmann 1.8 char hostname[1023];
74     gethostname(hostname,1023);
75 buchmann 1.6
76    
77     /*
78 buchmann 1.1 //This should only be called via CRAB
79     cout << " Deleting all cached files." << endl;
80     gSystem->Exec("ls -ltrh | grep root");
81 buchmann 1.6 if(!((Contains(hostname,"t3ui")||Contains(hostname,"t3wn"))))
82     {
83 buchmann 1.1 vector<string> all_files;
84     char currentpath[1024];
85     TString directory=TString(getcwd(currentpath,1024));
86     directory=directory+"/";
87     process_directory(directory,all_files);
88     for(int ifile=0;ifile<all_files.size();ifile++)
89     {
90     if(Contains(all_files[ifile],"_clean_splitup_")) {
91     dout << " About to delete " << all_files[ifile] << endl;
92     gSystem->Exec(((string)"rm "+all_files[ifile]).c_str());
93     }
94     }
95 buchmann 1.6 return true;
96     }
97     */
98 buchmann 1.8 if(!((Contains(hostname,"t3ui")||Contains(hostname,"t3wn")))) {
99 buchmann 1.6 cout << "Going to purge files in local_storage that have been loaded!" << endl;
100     for(int i=0;i<SUSYScanSpace::loaded_files.size();i++) {
101     gSystem->Exec(((string)"rm "+SUSYScanSpace::loaded_files[i]).c_str());
102     dout << " Purging : Deleted file " << SUSYScanSpace::loaded_files[i] << endl;
103     }
104     }
105 buchmann 1.1 return true;
106     }
107    
108     bool initialized_t2=false;
109    
110     int srmcpretries=0;
111    
112 buchmann 1.8 void load_local_scan_sample(int a, int b, int &scanfileindex,int scantype,int mglu, int mlsp, bool isretry=false) {
113     stringstream filetoload;
114     string samplename;
115 buchmann 1.1 filetoload << "/shome/buchmann/ntuples/"<<PlottingSetup::ScanSampleDirectory;
116 buchmann 1.4 if(scantype==mSUGRA) {
117     //filetoload << "/mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
118     filetoload.str("");
119     filetoload << "/shome/lbaeni/jzb/" << PlottingSetup::ScanSampleDirectory << "/mSUGRA_M0_" << mglu << "__M12_" << mlsp << ".root";
120     samplename="mSUGRA";
121     // filetoload << "/mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
122     }
123     if(scantype==SMS) {
124     filetoload << "/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
125     samplename="SMS";
126     }
127     if(scantype==GMSB) {
128     filetoload << "/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
129     samplename="GMSB";
130     }
131 buchmann 1.9
132     if(scantype==SMS && Contains(PlottingSetup::ScanSampleDirectory,"T1lh")) {
133     filetoload.str("");
134     filetoload << "/shome/lbaeni/jzb/" << PlottingSetup::ScanSampleDirectory << "/SMS_MassGlu_" << mglu << "__MassLSP_" << mlsp << ".root";
135     samplename="SMS";
136     }
137    
138 buchmann 1.3 if(scansample.collection.size()<1||!Contains(((scansample.collection)[(scansample.collection).size()-1]).filename,"_"+any2string(a)+"_"+any2string(b))) {
139 buchmann 1.1 dout << "The last sample is NOT the same one as the current one, possibly popping off last one and adding the new one." << endl;
140 buchmann 1.4 if((scansample.collection).size()>=1) {
141 buchmann 1.1 scansample.RemoveLastSample();
142     }
143     scanfileindex=(scansample.collection).size();
144     //New: Loading file when necessary, not before (avoiding high memory usage and startup times)
145 buchmann 1.4
146     scansample.AddSample(filetoload.str(),samplename,1,1,false,true,scanfileindex,kRed);
147     dout << "Just added the following file: " << filetoload.str() << endl;
148 buchmann 1.1 } else {
149     dout << "Last sample is the same as the current one. Recycling it." << endl;
150     scanfileindex=(scansample.collection).size()-1;
151     }
152     dout << " Going to use the following file: " << filetoload.str() << endl;
153 buchmann 1.8 }
154    
155     void load_remote_scan_sample(int a, int b, int &scanfileindex,int scantype,int mglu, int mlsp, bool isretry=false) {
156 buchmann 1.1 write_info(__FUNCTION__,"Hello, CRAB! Might need to load files ...");
157 buchmann 1.8 string samplename;
158 buchmann 1.1 PlottingSetup::limitpatience=PlottingSetup::limitpatienceCRAB;
159     stringstream copyfile;
160 buchmann 1.4 gSystem->Exec("mkdir -p local_storage");
161 buchmann 1.1 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;
162 buchmann 1.4 if(scantype==mSUGRA) {
163 buchmann 1.5 copyfile<< "/mSUGRA_M0_" << mglu << "__M12_" << mlsp << ".root local_storage/mSUGRA_M0_" << mglu << "__M12_" << mlsp << ".root ";
164 buchmann 1.4 samplename="mSUGRA";
165     }
166     if(scantype==SMS) {
167     copyfile << "/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root " << "local_storage/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root ";
168     samplename="SMS";
169     }
170     if(scantype==GMSB) {
171     copyfile << "/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root " << "local_storage/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root ";
172     samplename="GMSB";
173     }
174 buchmann 1.1 stringstream newfilename;
175 buchmann 1.5 // if(scantype==mSUGRA) newfilename << "local_storage/mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
176     if(scantype==mSUGRA) newfilename << "local_storage/mSUGRA_M0_" << mglu << "__M12_" << mlsp << ".root";
177 buchmann 1.4 if(scantype==SMS) newfilename << "local_storage/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
178     if(scantype==GMSB) newfilename << "local_storage/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root";
179 buchmann 1.1
180 buchmann 1.6 SUSYScanSpace::SUSYscantype=scantype;
181    
182 buchmann 1.7 if((scansample.collection).size()==0||!Contains(((scansample.collection)[(scansample.collection).size()-1]).filename,"_"+any2string(a)+"_"+any2string(b))||(a==0&&b==0)&&!initialized_t2) {
183 buchmann 1.1 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;
184 buchmann 1.5 while((scansample.collection).size()>0) {
185 buchmann 1.1 scansample.RemoveLastSample();
186     }
187 buchmann 1.5 cout << "New scanfileindex stems from scansample size: " << (scansample.collection).size() << endl;
188 buchmann 1.1 scanfileindex=(scansample.collection).size();
189     //New: Loading file when necessary, not before (avoiding high memory usage and startup times)
190 buchmann 1.5 if(1) {
191     if(!doesROOTFileExist(newfilename.str())) {
192     dout << "Going to download the scan file with the following copy command: " << copyfile.str() << endl;
193 buchmann 1.7 int retcode = gSystem->Exec(copyfile.str().c_str()); // download it if it hasn't been downloaded before
194     if(retcode==0) SUSYScanSpace::loaded_files.push_back(newfilename.str()); // adding it to the list of loaded files (will be deleted once we finish running
195 buchmann 1.5 }
196 buchmann 1.1 if(doesROOTFileExist(newfilename.str())) {
197     initialized_t2=true;
198 buchmann 1.5 cout << "Going to add the following sample: " << newfilename.str() << endl;
199     gSystem->Exec(("ls -ltrh "+newfilename.str()).c_str());
200     scansample.AddSample(newfilename.str(),samplename,1,1,false,true,scanfileindex,kRed);
201 buchmann 1.1 } else {
202     srmcpretries++;
203     if(srmcpretries<5) {
204     dout << "The file could not be loaded correctly - retrying!" << endl;
205     sleep(5);
206 buchmann 1.8 load_remote_scan_sample(a,b,scanfileindex,scantype,mglu,mlsp,true);
207 buchmann 1.1 } else {
208     dout << "Have tried 5 times to load this sample. Giving up now and failing the program execution" << endl;
209 buchmann 1.4 assert(srmcpretries<5);
210 buchmann 1.1 dout << "The command " << copyfile.str() << " failed to copy the file to the local storage (saving it as " << newfilename.str() << ")" << endl;
211     }
212     }
213     }
214     } else {
215     dout << "Last sample is the same as the current one. Recycling it." << endl;
216     scanfileindex=(scansample.collection).size()-1;
217     }
218 buchmann 1.8 }
219    
220     void load_scan_sample(int a, int b, int &scanfileindex,int scantype,int mglu, int mlsp, bool isretry=false) {
221    
222     // There is no need to define your sample here. That is done in Setup.C where you define the loading directory!
223 buchmann 1.3
224 buchmann 1.8 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;
225     char hostname[1023];
226     gethostname(hostname,1023);
227     if((Contains(hostname,"t3ui")||Contains(hostname,"t3wn"))) load_local_scan_sample(a,b,scanfileindex,scantype,mglu,mlsp,isretry);
228     else load_remote_scan_sample(a,b,scanfileindex,scantype,mglu,mlsp,isretry);
229 buchmann 1.1 }
230    
231 buchmann 1.3 /*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) {
232 buchmann 1.1 altxs=0;
233 buchmann 1.3 for(int iproc=1;iproc<=10;iproc++) altxs+=GetXSecForPointAndChannel(mglu,mlsp,xsec,iproc);
234 buchmann 1.1 return altxs;
235 buchmann 1.3 }*/
236 buchmann 1.1
237     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) {
238     bool runninglocally=true;
239     if(njobs>-1&&jobnumber>-1) {
240     runninglocally=false;
241     dout << "Running on the GRID (this is job " << jobnumber << "/" << njobs << ") for a jzb cut at " << jzb_cut[ibin] << endl;
242     } else {
243     dout << "Running locally " << endl;
244     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;
245     }
246    
247     string massgluname="MassGlu";
248     string massLSPname="MassLSP";
249     jzbSel=jzb_cut[ibin];
250     geqleq="geq";
251     automatized=true;
252     mcjzbexpression=mcjzb;
253    
254     string prefix="SMS_";
255     // up to here, everything is set up for SMS; now we need to switch stuff around if we're dealing with an mSUGRA scan!
256     int scantype=SMS;
257     TIter nextkey(fsyst->GetListOfKeys());
258     TKey *key;
259     while ((key = (TKey*)nextkey()))
260     {
261     TObject *obj = key->ReadObj();
262 buchmann 1.4 if(Contains((string)(obj->GetName()),"SUGRA")) scantype=mSUGRA;
263 buchmann 1.1 if(Contains((string)(obj->GetName()),"GMSB")) scantype=GMSB;
264     }
265    
266     map < pair<float, float>, map<string, float> > xsec;
267    
268     if(scantype==mSUGRA) {
269     massgluname="M0"; // this is the "x axis" in the limit plot (like the gluino in the SMS case)
270     massLSPname="M12"; // this is the "y axis" in the limit plot (like the LSP in the SMS case)
271     mglustart=m0start;
272 buchmann 1.10 xsec=getXsec(PlottingSetup::cbafbasedir+"/"+mSUGRAxsFile);
273    
274 buchmann 1.1 mgluend=m0end;
275     mglustep=m0step;
276     mLSPstart=m12start;
277     mLSPend=m12end;
278     mLSPstep=m12step;
279     prefix="mSUGRA_";
280     dout << "mSUGRA scan has been set up." << endl;
281     }
282     if(scantype==SMS) {
283     dout << "SMS scan has been set up." << endl;
284     write_warning(__FUNCTION__,"Don't have the correct XS file yet");
285     }
286     if(scantype==GMSB) {
287     string massgluname="MassGMSBGlu";
288     string massLSPname="MassGMSBChi";
289     prefix="GMSB_";
290     dout << "GMSB scan has been set up." << endl;
291     write_warning(__FUNCTION__,"Don't have the correct XS file yet");
292     }
293    
294     set_SUSY_style();
295    
296     TCanvas *limcanvas = new TCanvas("limcanvas","Limit canvas");
297    
298     int Npoints=0;
299     for(int mglu=(int)mglustart;mglu<=mgluend;mglu+=(int)mglustep) {
300     for (int mlsp=(int)mLSPstart;mlsp<=mLSPend&&(scantype==mSUGRA||mlsp<=mglu);mlsp+=(int)mLSPstep) Npoints++;
301     }
302     TH2F *mceff = (TH2F*)fsyst->Get((prefix+"efficiencymap"+any2string(jzb_cut[ibin])).c_str());
303     //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());
304     TH2F *fullerr = (TH2F*)fsyst->Get((prefix+"systotmap"+any2string(jzb_cut[ibin])).c_str());
305     TH2F *NEvents = (TH2F*)fsyst->Get((prefix+"Neventsmap"+any2string(jzb_cut[ibin])).c_str());
306     TH2F *lflip = (TH2F*)fsyst->Get((prefix+"flipmap"+any2string(jzb_cut[ibin])).c_str());
307 buchmann 1.3 TH2F *XSmap = (TH2F*)fsyst->Get((prefix+"XS"+any2string(jzb_cut[ibin])).c_str());
308 buchmann 1.1
309     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
310     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
311     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
312     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
313     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
314     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
315    
316     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);
317     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);
318 buchmann 1.5 TH2F *asymptoticmap = new TH2F((prefix+"asymptoticmap"+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);
319 buchmann 1.1
320    
321     if(!fullerr || !mceff || !NEvents) {
322     write_error(__FUNCTION__,"The supplied systematics file did not contain the correct histograms - please check the file");
323     dout << "mc eff address: " << mceff << " , error address: " << fullerr << " , NEvents address: " << NEvents << endl;
324     delete limcanvas;
325     return;
326     }
327    
328     bool doexpected=true;
329    
330     int ipoint=-1;
331     for(int mglu=(int)mglustart;mglu<=mgluend;mglu+=(int)mglustep) {
332     for (int mlsp=(int)mLSPstart;mlsp<=mLSPend&&(scantype==mSUGRA||mlsp<=mglu);mlsp+=(int)mLSPstep)
333     {
334     ipoint++;
335     if(!runninglocally&&!do_this_point(ipoint,(int)Npoints,(int)jobnumber,(int)njobs)) continue;
336     int GlobalBin=flipmap->FindBin(mglu,mlsp);
337     float currmceff=mceff->GetBinContent(mceff->FindBin(mglu,mlsp));
338     float currtoterr=(fullerr->GetBinContent(fullerr->FindBin(mglu,mlsp)))*currmceff;
339     float nevents=NEvents->GetBinContent(NEvents->FindBin(mglu,mlsp));
340     int flipped = (int)(lflip->GetBinContent(lflip->FindBin(mglu,mlsp)));
341     flipmap->SetBinContent(GlobalBin,flipped);
342     dout << "Looking at point " << ipoint << " / " << Npoints << " with masses " << massgluname << " = " << mglu << " and " << massLSPname << " = " << mlsp << endl;
343     dout << "Found : MCeff=" << currmceff << " and total error=" << currtoterr << " and Nevents=" << nevents << endl;
344     string plotfilename=(string)(TString(massgluname)+TString(any2string(mglu))+TString("__")+TString(massLSPname)+TString(any2string(mlsp))+TString(".png"));
345     if(currmceff<=0||currtoterr<=0||nevents==0) {
346     dout << " Nothing to work with, skipping this point." << endl;
347     continue;
348     }
349     vector<float> sigmas;
350 buchmann 1.9 if(scantype!=mSUGRA) {
351     sigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped,true);
352     write_warning(__FUNCTION__,"Temporarily doing asymptotic limits!!!!");
353     }
354 buchmann 1.3 else {
355     //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"!
356     vector<float> asigmas;
357     asigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped,true); // asymptotic limit first
358     float strength=asigmas[0]/XSmap->GetBinContent(GlobalBin);
359 buchmann 1.9 /* if(strength>0.5&&strength<2) {
360 buchmann 1.5 sigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped); // asymptotic limit first
361     asymptoticmap->SetBinContent(GlobalBin,0);
362 buchmann 1.9 } else {*/
363 buchmann 1.5 sigmas=asigmas;
364     asymptoticmap->SetBinContent(GlobalBin,1);
365 buchmann 1.9 /* }*/
366 buchmann 1.3 exclmap->SetBinContent(GlobalBin,strength);
367     }
368    
369 buchmann 1.1 if(sigmas[0]>0) limitmap->SetBinContent(GlobalBin,sigmas[0]); //anything else is an error code
370     if(sigmas.size()>1) {
371     explimitmap->SetBinContent(GlobalBin,sigmas[1]);
372     exp1plimitmap->SetBinContent(GlobalBin,sigmas[2]);
373     exp1mlimitmap->SetBinContent(GlobalBin,sigmas[3]);
374     exp2plimitmap->SetBinContent(GlobalBin,sigmas[4]);
375     exp2mlimitmap->SetBinContent(GlobalBin,sigmas[5]);
376     }
377    
378     dout << "An upper limit has been added for this point ( " << massgluname << " = " << mglu << " and " << massLSPname << " = " << mlsp << " ) at " << sigmas[0] << endl;
379     }
380     }
381    
382     prepare_scan_axis(limitmap,scantype);
383 buchmann 1.6 gSystem->Exec("mkdir -p output");
384 buchmann 1.1 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 :-)
385     outputfile->cd();
386     limitmap->Write();
387 buchmann 1.8 XSmap->Write();
388 buchmann 1.1 flipmap->Write();
389 buchmann 1.5 asymptoticmap->Write();
390 buchmann 1.1 if(doexpected) {
391     explimitmap->Write();
392     exp1plimitmap->Write();
393     exp1mlimitmap->Write();
394     exp2plimitmap->Write();
395     exp2mlimitmap->Write();
396     }
397 buchmann 1.5 if(scantype==mSUGRA) exclmap->Write();
398 buchmann 1.1 outputfile->Close();
399     delete limcanvas;
400     }
401    
402     void establish_SUSY_limits(string mcjzb,string datajzb,vector<float> jzb_cut,bool requireZ, float peakerror, string fsystfilename, float njobs=-1, float jobnumber=-1) {
403     dout << "Starting the SUSY scan from systematics file now with all " << jzb_cut.size() << " bin(s)" << " and using " << fsystfilename << endl;
404     TFile *fsyst = new TFile(fsystfilename.c_str());
405     if(!fsyst) {
406     write_error(__FUNCTION__,"You provided an invalid systematics file. Please change it ... ");
407     return;
408     }
409     for(int ibin=0;ibin<jzb_cut.size();ibin++) {
410     establish_SUSY_limits(mcjzb,datajzb,jzb_cut,requireZ, peakerror, fsyst, ibin,njobs, jobnumber);
411     }
412     }
413    
414    
415    
416    
417 buchmann 1.9 void scan_SUSY_parameter_space(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) {
418    
419     if(shapeanalysis&&ibin>0) return; // pointless since we're doing a shape analysis :-)
420    
421 buchmann 1.1 bool runninglocally=true;
422     if(njobs>-1&&jobnumber>-1) {
423     runninglocally=false;
424 buchmann 1.2 dout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
425 buchmann 1.1 dout << "Running on the GRID (this is job " << jobnumber << "/" << njobs << ") for a jzb cut at " << jzb_cut[ibin] << endl;
426     } else {
427     dout << "Running locally " << endl;
428     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;
429     }
430    
431     jzbSel=jzb_cut[ibin];
432     geqleq="geq";
433     automatized=true;
434     mcjzbexpression=mcjzb;
435    
436     string massgluname="MassGlu";
437     string massLSPname="MassLSP";
438    
439     string prefix="SMS_";
440     // up to here, everything is set up for SMS; now we need to switch stuff around if we're dealing with an mSUGRA scan!
441     int scantype=SMS;
442 buchmann 1.9 if(Contains((PlottingSetup::ScanSampleDirectory),"SUGRA")) scantype=mSUGRA;
443     if(Contains((PlottingSetup::ScanSampleDirectory),"GMSB")) scantype=GMSB;
444 buchmann 1.1
445     if(scantype==mSUGRA) {
446     massgluname="M0"; // this is the "x axis" in the limit plot (like the gluino in the SMS case)
447     massLSPname="M12"; // this is the "y axis" in the limit plot (like the LSP in the SMS case)
448     mglustart=m0start;
449     mgluend=m0end;
450     mglustep=m0step;
451     mLSPstart=m12start;
452     mLSPend=m12end;
453     mLSPstep=m12step;
454     prefix="mSUGRA_";
455     dout << "mSUGRA scan has been set up." << endl;
456     }
457     if(scantype==SMS) {
458     dout << "SMS scan has been set up." << endl;
459     }
460     if(scantype==GMSB) {
461     dout << "GMSB scan has been set up." << endl;
462     prefix="GMSB_";
463     massgluname="MassGMSBGlu"; // this is the "x axis" in the limit plot (like the gluino in the SMS case)
464     massLSPname="MassGMSBChi"; // this is the "y axis" in the limit plot (like the LSP in the SMS case)
465     }
466    
467     Int_t MyPalette[100];
468     Double_t r[] = {0., 0.0, 1.0, 1.0, 1.0};
469     Double_t g[] = {0., 0.0, 0.0, 1.0, 1.0};
470     Double_t b[] = {0., 1.0, 0.0, 0.0, 1.0};
471     Double_t stop[] = {0., .25, .50, .75, 1.0};
472     Int_t FI = TColor::CreateGradientColorTable(5, stop, r, g, b, 110);
473     for (int i=0;i<100;i++) MyPalette[i] = FI+i;
474    
475     gStyle->SetPalette(100, MyPalette);
476    
477     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
478     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
479     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
480     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
481     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
482     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
483    
484 buchmann 1.2 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);
485     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);
486     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);
487     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);
488     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);
489     TH2F *efficiencymapmet = new TH2F((prefix+"efficiencymapmet"+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);
490     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);
491     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);
492 buchmann 1.1 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);
493     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);
494 buchmann 1.3 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-
495     mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
496 buchmann 1.1 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);
497 buchmann 1.2 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);
498     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);
499     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);
500     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);
501     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);
502 buchmann 1.3 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);
503 buchmann 1.8 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);
504 buchmann 1.9 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);
505 buchmann 1.1
506 buchmann 1.2 TH2F *imposedxmap;
507     TH2F *realxmap;
508 buchmann 1.1
509 buchmann 1.3 TH2F *centermap;
510     TH2F *widthmap;
511    
512 buchmann 1.2 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);
513    
514     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);
515 buchmann 1.1
516 buchmann 1.2 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);
517 buchmann 1.1
518 buchmann 1.2 if(ibin==0) {
519 buchmann 1.3 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);
520     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);
521    
522 buchmann 1.2 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);
523     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);
524     }
525    
526     vector<int> susy_Zdecay_origin;
527     vector<TH2F*> h_susy_Zdecay_origin;
528     // 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)
529     for(int i=1000022;i<1000039;i++) {
530     if(ibin>0) continue;
531     susy_Zdecay_origin.push_back(i);
532     TH2F *susy_dec = new TH2F((prefix+"SUSY_Decay_From_"+any2string(i)+"_to_Z").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     h_susy_Zdecay_origin.push_back(susy_dec);
534     }
535     for(int i=2000006;i<2000016;i++) {
536     if(ibin>0) continue;
537     susy_Zdecay_origin.push_back(i);
538     TH2F *susy_dec = new TH2F((prefix+"SUSY_Decay_From_"+any2string(i)+"_to_Z").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     h_susy_Zdecay_origin.push_back(susy_dec);
540     }
541    
542     TH2F *PrimaryZSource = new TH2F((prefix+"PrimaryZSource").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);
543     TH2F *SecondaryZSource = new TH2F((prefix+"SecondaryZSource").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);
544    
545 buchmann 1.1 write_warning(__FUNCTION__,"CURRENTLY SWITCHING AUTOMATIZED MODE OFF!");automatized=false;
546    
547     float rightmargin=gStyle->GetPadRightMargin();
548     gStyle->SetPadRightMargin(0.15);
549    
550     TCanvas *limcanvas = new TCanvas("limcanvas","Limit canvas");
551    
552    
553     bool doexpected=true;
554    
555     int Npoints=0;
556     for(int mglu=(int)mglustart;mglu<=mgluend;mglu+=(int)mglustep) {
557     for (int mlsp=(int)mLSPstart;mlsp<=mLSPend&&(scantype==mSUGRA||mlsp<=mglu);mlsp+=(int)mLSPstep) Npoints++;
558     }
559 buchmann 1.3
560     map < pair<float, float>, map<string, float> > xsec;
561     if(scantype==mSUGRA) {
562     xsec=getXsec(PlottingSetup::mSUGRAxsFile);
563     }
564    
565 buchmann 1.1 int ipoint=-1;
566     for(int mglu=(int)mglustart;mglu<=mgluend;mglu+=(int)mglustep) {
567     for (int mlsp=(int)mLSPstart;mlsp<=mLSPend&&(scantype==mSUGRA||mlsp<=mglu);mlsp+=(int)mLSPstep)
568     {
569     ipoint++;
570     int GlobalBin=efficiencymap->FindBin(mglu,mlsp);
571     if(!runninglocally&&!do_this_point(ipoint,(int)Npoints,(int)jobnumber,(int)njobs)) continue;
572     int flipped=0;
573     float result=-987,resulterr=-987;
574     float m0trees = PlottingSetup::mgluend;
575     float m12trees = PlottingSetup::mLSPend;
576     int a = int((PlottingSetup::ScanXzones*mlsp)/(m12trees+1));
577     int b = int((PlottingSetup::ScanYzones*mglu)/(m0trees+1));
578     int scanfileindex=PlottingSetup::ScanYzones*a+b;
579 buchmann 1.4 load_scan_sample(a,b,scanfileindex,scantype,mglu,mlsp);
580 buchmann 1.1
581     clock_t start,finish;
582     start = clock(); // starting the clock to measure how long the computation takes!
583     stringstream addcut;
584     addcut << "(TMath::Abs("<<massgluname<<"-"<<mglu<<")<5)&&(TMath::Abs("<<massLSPname<<"-"<<mlsp<<")<5)";
585     (scansample.collection)[scanfileindex].events->Draw("eventNum",addcut.str().c_str(),"goff");
586     float nevents = (scansample.collection)[scanfileindex].events->GetSelectedRows();
587     vector<vector<float> > systematics;
588     if(nevents<10) {
589 buchmann 1.3 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;
590 buchmann 1.1 continue;
591     } else {
592 buchmann 1.2 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;
593     if(PlottingSetup::RestrictToMassPeak==true) cout << "Analysis type: on-peak" << endl;
594     if(PlottingSetup::RestrictToMassPeak==false) cout << "Analysis type: offpeak" << endl;
595     }
596     float hitrecord=0;
597     int hitrecordholder=0;
598     float hitsecondplace=0;
599     int hitsecondplaceholder=0;
600     for(int imo=0;imo<susy_Zdecay_origin.size()&&ibin==0;imo++) {
601     if(ibin>0) continue;
602     int MotherParticlePDGid = susy_Zdecay_origin[imo];
603     (scansample.collection)[scanfileindex].events->Draw("eventNum",("("+addcut.str()+")&&(abs(SourceOfZ-"+any2string(MotherParticlePDGid)+")<0.5)").c_str(),"goff");
604     float hits = (scansample.collection)[scanfileindex].events->GetSelectedRows();
605     hits/=nevents;
606     h_susy_Zdecay_origin[imo]->SetBinContent(GlobalBin,hits);
607     if(hits>hitrecord) {
608     hitsecondplace=hitrecord;
609     hitsecondplaceholder=hitrecordholder;
610     hitrecord=hits;
611     hitrecordholder=imo;
612     }
613     if(hits<hitrecord&&hits>hitsecondplace) {
614     hitsecondplace=hits;
615     hitsecondplaceholder=imo;
616     }
617     }
618     if(ibin==0) {
619     PrimaryZSource->SetBinContent(GlobalBin,hitrecordholder);
620     SecondaryZSource->SetBinContent(GlobalBin,hitsecondplaceholder);
621     }
622    
623     if(ibin==0) {
624     TH1F *realxhisto = new TH1F("realxhisto","realxhisto",1000,0,1);
625     (scansample.collection)[scanfileindex].events->Draw("realx>>realxhisto",addcut.str().c_str(),"goff");
626     realxmap->SetBinContent(GlobalBin,realxhisto->GetMean());
627     realxmap->SetBinError(GlobalBin,realxhisto->GetRMS());
628     (scansample.collection)[scanfileindex].events->Draw("imposedx>>realxhisto",addcut.str().c_str(),"goff");
629     imposedxmap->SetBinContent(GlobalBin,realxhisto->GetMean());
630     imposedxmap->SetBinError(GlobalBin,realxhisto->GetRMS());
631 buchmann 1.3 stringstream specialexpression;
632     specialexpression << mcjzb << ">>jzbhisto";
633     TH1F *jzbhisto = new TH1F("jzbhisto","jzbhisto",100,-500,500);
634     (scansample.collection)[scanfileindex].events->Draw(specialexpression.str().c_str(),addcut.str().c_str(),"goff");
635     centermap->SetBinContent(GlobalBin,jzbhisto->GetMean());
636     widthmap->SetBinContent(GlobalBin,jzbhisto->GetRMS());
637 buchmann 1.2 delete realxhisto;
638 buchmann 1.3 delete jzbhisto;
639 buchmann 1.1 }
640    
641    
642 buchmann 1.5 if(scantype==mSUGRA) {
643     float absxs=0;
644     for(int i=0;i<12;i++) absxs+=GetXSecForPointAndChannel(mglu,mlsp,xsec,i);
645 buchmann 1.8 TFile *FilterEffFile = new TFile(PlottingSetup::FilterEfficiencyFile.c_str());
646     TH2F *FilterEfficiency = (TH2F*)FilterEffFile->Get("FilterEfficiency");
647     float filtereff = FilterEfficiency->GetBinContent(FilterEfficiency->FindBin(mglu,mlsp));
648     FilterEff->SetBinContent(GlobalBin,filtereff);
649     FilterEffFile->Close();
650     absXSmap->SetBinContent(GlobalBin,absxs);
651     absxs*=filtereff;
652 buchmann 1.5 XSmap->SetBinContent(GlobalBin,absxs);
653     }
654    
655 buchmann 1.1 if(nevents!=0&&(efficiencyonly||systematicsonly)) {
656 buchmann 1.3 Value effwosigcont = MCefficiency((scansample.collection)[scanfileindex].events,result,resulterr,flipped,mcjzb,requireZ,(int)nevents,scantype,xsec,addcut.str(),-1);
657 buchmann 1.2 if(result<0&&allowflipping) {
658 buchmann 1.1 write_info(__FUNCTION__,"Found negative efficiency (signal distribution populating JZB<0 more than JZB>0) - flipping analysis!");
659     flipped=1;
660 buchmann 1.3 MCefficiency((scansample.collection)[scanfileindex].events,result,resulterr,flipped,mcjzb,requireZ,(int)nevents,scantype,xsec,addcut.str(),-1);
661 buchmann 1.1 }
662     efficiencymap->SetBinContent(GlobalBin,result);
663     efficiencymap->SetBinError(GlobalBin,resulterr);
664     flipmap->SetBinContent(GlobalBin,flipped);
665     noscefficiencymap->SetBinContent(GlobalBin,effwosigcont.getValue());
666     noscefficiencymap->SetBinError(GlobalBin,effwosigcont.getError());
667 buchmann 1.2 efficiencymapContam->SetBinContent(GlobalBin,effwosigcont.getValue()-result);
668     efficiencymapContam->SetBinError(GlobalBin,TMath::Sqrt(effwosigcont.getError()+resulterr));
669 buchmann 1.1 //Partial efficiencies
670     float resultAcc, resultID, resultJets, resultSignal, resultmet;
671     float resulterrAcc, resulterrID, resulterrJets, resulterrSignal, resulterrmet;
672     MCPartialefficiency((scansample.collection)[scanfileindex].events,resultAcc,resulterrAcc,flipped,mcjzb,requireZ,(int)nevents,addcut.str(),-1, 1);
673     MCPartialefficiency((scansample.collection)[scanfileindex].events,resultID,resulterrID,flipped,mcjzb,requireZ,(int)nevents,addcut.str(),-1, 2);
674     MCPartialefficiency((scansample.collection)[scanfileindex].events,resultJets,resulterrJets,flipped,mcjzb,requireZ,(int)nevents,addcut.str(),-1, 3);
675     MCPartialefficiency((scansample.collection)[scanfileindex].events,resultSignal,resulterrSignal,flipped,mcjzb,requireZ,(int)nevents,addcut.str(),-1, 4);
676     MCPartialefficiency((scansample.collection)[scanfileindex].events,resultmet,resulterrmet,flipped,mcjzb,requireZ,(int)nevents,addcut.str(),-1, 5);
677     efficiencymapAcc->SetBinContent(GlobalBin,resultAcc);
678     efficiencymapAcc->SetBinError(GlobalBin,resulterrAcc);
679     efficiencymapmet->SetBinContent(GlobalBin,resultmet);
680     efficiencymapmet->SetBinError(GlobalBin,resulterrmet);
681     efficiencymapID->SetBinContent(GlobalBin,resultID);
682     efficiencymapID->SetBinError(GlobalBin,resulterrID);
683     efficiencymapJets->SetBinContent(GlobalBin,resultJets);
684     efficiencymapJets->SetBinError(GlobalBin,resulterrJets);
685     efficiencymapSignal->SetBinContent(GlobalBin,resultSignal);
686     efficiencymapSignal->SetBinError(GlobalBin,resulterrSignal);
687     finish = clock();
688     }
689     Neventsmap->SetBinContent(GlobalBin,nevents);
690     ipointmap->SetBinContent(GlobalBin,ipoint);
691     if(efficiencyonly) continue;
692 buchmann 1.9
693     if(!shapeanalysis) {
694     do_systematics_for_one_file((scansample.collection)[scanfileindex].events,(int)nevents,"SUSY SCAN", systematics,flipped, xsec, mcjzb,datajzb,peakerror,requireZ, addcut.str(),true,scantype);
695     float JZBcutat = systematics[0][0];
696     float mceff = systematics[0][1];
697     float mcefferr = systematics[0][2];//MC stat error
698     float toterr = systematics[0][4];
699     float sys_jes = systematics[0][5]; // Jet Energy Scale
700     float sys_jsu = systematics[0][6]; // JZB scale uncertainty
701     float sys_res = systematics[0][7]; // resolution
702     float mcwoscef = systematics[0][8]; // efficiency without signal contamination
703     float mcwoscefr= systematics[0][9]; // error on efficiency without signal contamination
704     float sys_pdf = 0;
705     if(systematics[0].size()>10) sys_pdf = systematics[0][10]; // PDF
706     cout << "Going to store: total error is: " << toterr << " and the relative tot error, i.e. total error / mc efficiency is " << toterr/mceff << endl;
707     if(mceff!=mceff||toterr!=toterr||mceff<0 && (!systematicsonly&&!efficiencyonly)) {
708     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;
709     continue;
710     } else {
711     // Systematics and efficiencies make sense (i.e. non-zero, and all numbers)
712     if(!systematicsonly&&!efficiencyonly) {
713     dout << "Calculating limit now for "<<massgluname<<"="<<mglu<<" , "<<massLSPname<<"="<<mlsp <<endl;
714     vector<float> sigmas;
715     string plotfilename=(string)(TString((scansample.collection)[scanfileindex].samplename)+TString(massgluname)+TString(any2string(mglu))+TString("__")+TString(massLSPname)+TString(any2string(mlsp))+TString(".png"));
716     sigmas=compute_one_upper_limit(mceff,toterr,ibin,mcjzb,plotfilename,true,flipped);
717     // do_limit_wrapper(mceff,toterr,ibin,mcjzb,sigmas,plotfilename);
718     if(sigmas[0]>-0.5) { // negative sigmas are the error signature of do_limit_wrapper, so we want to exclude them.
719     limitmap->SetBinContent(GlobalBin,sigmas[0]);
720     if(sigmas.size()>1) {
721     explimitmap->SetBinContent(GlobalBin,sigmas[1]);
722     exp1plimitmap->SetBinContent(GlobalBin,sigmas[2]);
723     exp1mlimitmap->SetBinContent(GlobalBin,sigmas[3]);
724     exp2plimitmap->SetBinContent(GlobalBin,sigmas[4]);
725     exp2mlimitmap->SetBinContent(GlobalBin,sigmas[5]);
726     }
727     sysjesmap->SetBinContent(GlobalBin,sys_jes);
728     sysjsumap->SetBinContent(GlobalBin,sys_jsu);
729     sysresmap->SetBinContent(GlobalBin,sys_res);
730     syspdfmap->SetBinContent(GlobalBin,sys_pdf);
731     systotmap->SetBinContent(GlobalBin,toterr/mceff);//total relative (!) error
732     sysstatmap->SetBinContent(GlobalBin,mcefferr);//total relative (!) error
733     dout << "A limit has been added at " << sigmas[0] << " for m_{glu}="<<mglu << " and m_{lsp}="<<mlsp<<endl;
734     }//end of if sigma is positive
735     }//end of not systematics only condition
736     if(systematicsonly) {
737     sysjesmap->SetBinContent(GlobalBin,sys_jes);
738     sysjsumap->SetBinContent(GlobalBin,sys_jsu);
739     sysresmap->SetBinContent(GlobalBin,sys_res);
740     syspdfmap->SetBinContent(GlobalBin,sys_pdf);
741     systotmap->SetBinContent(GlobalBin,toterr/mceff);//total relative (!) error
742     sysstatmap->SetBinContent(GlobalBin,mcefferr);//total relative (!) error
743     }
744     }//efficiency is valid
745 buchmann 1.1 } else {
746 buchmann 1.9 //shape analysis
747     cout << "This is the shape analysis. Have fun with it!" << endl;
748 buchmann 1.11 float quickmceff,quickmcefferr;
749     MCefficiency((scansample.collection)[scanfileindex].events,quickmceff,quickmcefferr,flipped,mcjzb,requireZ,(int)nevents,scantype,xsec,addcut.str(),-1);
750     if(quickmceff==0) {
751     write_warning(__FUNCTION__,"The MC efficiency is zero - need to skip this point.");
752     continue;
753     }
754 buchmann 1.9 stringstream PointName;
755     PointName << massgluname << "_" << mglu << "__" << massLSPname << "_" << mlsp;
756 buchmann 1.10 SUSYScanSpace::SavedMLSP=mlsp;
757     SUSYScanSpace::SavedMGlu=mglu;
758     SUSYScanSpace::SavedMLSPname=massLSPname;
759     SUSYScanSpace::SavedMGluname=massgluname;
760     ShapeDroplet droplet = LimitsFromShapes((scansample.collection)[scanfileindex].events,addcut.str().c_str(),PointName.str(),mcjzb,datajzb,jzb_cut,peakerror,peakerror);
761 buchmann 1.9 limitmap->SetBinContent(GlobalBin,droplet.observed);
762     explimitmap->SetBinContent(GlobalBin,droplet.expected);
763     exp1plimitmap->SetBinContent(GlobalBin,droplet.expectedPlus1Sigma);
764     exp1mlimitmap->SetBinContent(GlobalBin,droplet.expectedMinus1Sigma);
765     exp2plimitmap->SetBinContent(GlobalBin,droplet.expectedPlus2Sigma);
766     exp2mlimitmap->SetBinContent(GlobalBin,droplet.expectedMinus2Sigma);
767    
768     sysjesmap->SetBinContent(GlobalBin,droplet.JES);
769     sysjsumap->SetBinContent(GlobalBin,droplet.JSU);
770     syspdfmap->SetBinContent(GlobalBin,droplet.PDF);
771     systotmap->SetBinContent(GlobalBin,droplet.toterr);
772     sysstatmap->SetBinContent(GlobalBin,droplet.staterr);
773     }
774 buchmann 1.1 finish = clock();
775     timemap->SetBinContent(GlobalBin,((float(finish)-float(start))/CLOCKS_PER_SEC));
776     }
777     }
778    
779     prepare_scan_axis(limitmap,scantype);
780     prepare_scan_axis(sysjesmap,scantype);
781     prepare_scan_axis(sysjsumap,scantype);
782     prepare_scan_axis(sysresmap,scantype);
783     prepare_scan_axis(syspdfmap,scantype);
784     prepare_scan_axis(systotmap,scantype);
785     prepare_scan_axis(sysstatmap,scantype);
786     prepare_scan_axis(timemap,scantype);
787    
788     if(!systematicsonly&&!efficiencyonly) {
789     limcanvas->cd();
790     limitmap->Draw("COLZ");
791     string titletobewritten="Limits in LSP-Glu plane";
792     if(scantype==mSUGRA) titletobewritten="Limits in m_{1/2}-m_{0} plane";
793     if(scantype==GMSB) titletobewritten="Limits in Chi-Glu plane";
794     TText *title = write_title(titletobewritten);
795     title->Draw();
796     if(runninglocally) {
797     CompleteSave(limcanvas,"SUSYScan/Limits_JZB_geq"+any2string(jzb_cut[ibin]));
798     } else {
799     TFile *outputfile;
800     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 :-)
801     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
802     limitmap->Write();
803     if(doexpected) {
804     explimitmap->Write();
805     exp1plimitmap->Write();
806     exp1mlimitmap->Write();
807     exp2plimitmap->Write();
808     exp2mlimitmap->Write();
809     }
810     sysjesmap->Write();
811     sysjsumap->Write();
812     sysresmap->Write();
813     efficiencymap->Write();
814 buchmann 1.3 if(ibin==0) centermap->Write();
815     if(ibin==0) widthmap->Write();
816 buchmann 1.1 flipmap->Write();
817     efficiencymapmet->Write();
818     efficiencymapAcc->Write();
819     efficiencymapID->Write();
820     efficiencymapJets->Write();
821     efficiencymapSignal->Write();
822     noscefficiencymap->Write();
823 buchmann 1.2 efficiencymapContam->Write();
824 buchmann 1.1 syspdfmap->Write();
825     systotmap->Write();
826     sysstatmap->Write();
827 buchmann 1.3 XSmap->Write();
828 buchmann 1.8 absXSmap->Write();
829     FilterEff->Write();
830 buchmann 1.2 if(ibin==0) imposedxmap->Write();
831     if(ibin==0) realxmap->Write();
832 buchmann 1.1 Neventsmap->Write();
833     ipointmap->Write();
834     timemap->Write();
835 buchmann 1.2 for(int i=0;i<susy_Zdecay_origin.size();i++) {
836     h_susy_Zdecay_origin[i]->Write();
837     }
838 buchmann 1.3 outputfile->Close();
839 buchmann 1.1 }
840     }
841     if(systematicsonly) { // systematics only :
842     limcanvas->cd();
843     sysjesmap->Draw("COLZ");
844     string titletobewritten="Jet Energy scale in LSP-Glu plane";
845     if(scantype==mSUGRA) titletobewritten="Limits in m_{1/2}-m_{0} plane";
846     if(scantype==GMSB) titletobewritten="Limits in Chi-Glu plane";
847     TText *title = write_title(titletobewritten);
848     title->Draw();
849     if(runninglocally) {
850     CompleteSave(limcanvas,"SUSYScan/JES_geq"+any2string(jzb_cut[ibin]));
851     }
852     sysjsumap->Draw("COLZ");
853     titletobewritten="JZB Scale Uncertainty in LSP-Glu plane";
854     if(scantype==mSUGRA) titletobewritten="JZB Scale Uncertainty in m_{1/2}-m_{0} plane";
855     if(scantype==GMSB) titletobewritten="JZB Scale Uncertainty in Chi-Glu plane";
856     TText *title2 = write_title(titletobewritten);
857     title2->Draw();
858     if(runninglocally) {
859     CompleteSave(limcanvas,"SUSYScan/JSU_geq"+any2string(jzb_cut[ibin]));
860     }
861     sysresmap->Draw("COLZ");
862     titletobewritten="Resolution in LSP-Glu plane";
863     if(scantype==mSUGRA) titletobewritten="Resolution in m_{1/2}-m_{0} plane";
864     if(scantype==GMSB) titletobewritten="Resolution in Chi-Glu plane";
865     TText *title3 = write_title(titletobewritten);
866     title3->Draw();
867     if(runninglocally) {
868     CompleteSave(limcanvas,"SUSYScan/Resolution_geq"+any2string(jzb_cut[ibin]));
869     }
870    
871     if(!runninglocally) {
872     TFile *outputfile=new TFile(("output/DistributedSystematics_job"+string(any2string(jobnumber))+"_of_"+string(any2string(njobs))+".root").c_str(),"UPDATE");
873     sysjesmap->Write();
874     sysjsumap->Write();
875     sysresmap->Write();
876     flipmap->Write();
877 buchmann 1.3 if(ibin==0) centermap->Write();
878     if(ibin==0) widthmap->Write();
879 buchmann 1.1 efficiencymap->Write();
880     efficiencymapmet->Write();
881     efficiencymapAcc->Write();
882     efficiencymapID->Write();
883     efficiencymapJets->Write();
884     efficiencymapSignal->Write();
885     noscefficiencymap->Write();
886 buchmann 1.2 efficiencymapContam->Write();
887 buchmann 1.1 Neventsmap->Write();
888     ipointmap->Write();
889     syspdfmap->Write();
890     systotmap->Write();
891     sysstatmap->Write();
892 buchmann 1.8 absXSmap->Write();
893 buchmann 1.9 XSmap->Write();
894 buchmann 1.8 FilterEff->Write();
895 buchmann 1.2 if(ibin==0) imposedxmap->Write();
896     if(ibin==0) realxmap->Write();
897 buchmann 1.1 timemap->Write();
898 buchmann 1.2 for(int i=0;i<susy_Zdecay_origin.size();i++) {
899     h_susy_Zdecay_origin[i]->Write();
900     }
901 buchmann 1.3 outputfile->Close();
902 buchmann 1.1 }
903     }//end of systematics only
904     if(efficiencyonly) {
905     limcanvas->cd();
906     string titletobewritten="Efficiencies in LSP-Glu plane";
907     if(scantype==mSUGRA) titletobewritten="Efficiencies in m_{1/2}-m_{0} plane";
908     if(scantype==GMSB) titletobewritten="Efficiencies in Chi-Glu plane";
909     TText *title = write_title(titletobewritten);
910     efficiencymap->Draw("COLZ");
911     title->Draw();
912     if(runninglocally) {
913     CompleteSave(limcanvas,"SUSYScan/Efficiency_geq"+any2string(jzb_cut[ibin]));
914     }
915     limcanvas->cd();
916     efficiencymapmet->Draw("COLZ");
917     title->Draw();
918     if(runninglocally) {
919     CompleteSave(limcanvas,"SUSYScan/EfficiencyMet_geq"+any2string(jzb_cut[ibin]));
920     }
921     limcanvas->cd();
922     efficiencymapAcc->Draw("COLZ");
923     title->Draw();
924     if(runninglocally) {
925     CompleteSave(limcanvas,"SUSYScan/EfficiencyAcc_geq"+any2string(jzb_cut[ibin]));
926     }
927     limcanvas->cd();
928     efficiencymapID->Draw("COLZ");
929     title->Draw();
930     if(runninglocally) {
931     CompleteSave(limcanvas,"SUSYScan/EfficiencyID_geq"+any2string(jzb_cut[ibin]));
932     }
933     limcanvas->cd();
934     efficiencymapJets->Draw("COLZ");
935     title->Draw();
936     if(runninglocally) {
937     CompleteSave(limcanvas,"SUSYScan/EfficiencyJets_geq"+any2string(jzb_cut[ibin]));
938     }
939     limcanvas->cd();
940     efficiencymapSignal->Draw("COLZ");
941     title->Draw();
942     if(runninglocally) {
943     CompleteSave(limcanvas,"SUSYScan/EfficiencySignal_geq"+any2string(jzb_cut[ibin]));
944     }
945     limcanvas->cd();
946     noscefficiencymap->Draw("COLZ");
947     title->Draw();
948     if(runninglocally) {
949     CompleteSave(limcanvas,"SUSYScan/NoEfficiency_geq"+any2string(jzb_cut[ibin]));
950     }
951     limcanvas->cd();
952     sysjesmap->Draw("COLZ");
953     titletobewritten="N(events) in LSP-Glu plane";
954     if(scantype==mSUGRA) titletobewritten="N(events) in m_{1/2}-m_{0} plane";
955     if(scantype==GMSB) titletobewritten="Resolution in Chi-Glu plane";
956     TText *title2 = write_title(titletobewritten);
957     title2->Draw();
958     if(runninglocally) {
959     CompleteSave(limcanvas,"SUSYScan/Nevents_geq"+any2string(jzb_cut[ibin]));
960     }
961     if(!runninglocally) {
962     TFile *outputfile=new TFile(("output/DistributedSystematics_job"+string(any2string(jobnumber))+"_of_"+string(any2string(njobs))+".root").c_str(),"UPDATE");
963     ipointmap->Write();
964     Neventsmap->Write();
965     noscefficiencymap->Write();
966 buchmann 1.2 efficiencymapContam->Write();
967 buchmann 1.1 efficiencymap->Write();
968 buchmann 1.3 if(ibin==0) centermap->Write();
969     if(ibin==0) widthmap->Write();
970 buchmann 1.1 flipmap->Write();
971     efficiencymapmet->Write();
972     efficiencymapAcc->Write();
973     efficiencymapID->Write();
974     efficiencymapJets->Write();
975     efficiencymapSignal->Write();
976     timemap->Write();
977     outputfile->Close();
978     }
979     }//end of efficiencies only
980    
981 buchmann 1.2 for(int i=0;i<susy_Zdecay_origin.size();i++) {
982     delete h_susy_Zdecay_origin[i];
983     }
984    
985     delete PrimaryZSource;
986     delete SecondaryZSource;
987    
988    
989 buchmann 1.1 delete limitmap;
990     delete exclmap;
991     delete sysjesmap;
992     delete sysjsumap;
993     delete sysresmap;
994     delete noscefficiencymap;
995 buchmann 1.2 delete efficiencymapContam;
996 buchmann 1.1 delete efficiencymap;
997     delete efficiencymapmet;
998     delete efficiencymapAcc;
999     delete efficiencymapID;
1000     delete efficiencymapJets;
1001     delete efficiencymapSignal;
1002     delete Neventsmap;
1003     delete ipointmap;
1004     delete syspdfmap;
1005     delete systotmap;
1006     delete sysstatmap;
1007 buchmann 1.3 delete XSmap;
1008 buchmann 1.8 delete absXSmap;
1009 buchmann 1.1 delete limcanvas;
1010     }
1011    
1012 buchmann 1.9 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) {
1013 buchmann 1.1 dout << "Starting the SUSY scan now with all " << jzb_cut.size() << " bin(s)" << endl;
1014     for(int ibin=0;ibin<jzb_cut.size();ibin++) {
1015 buchmann 1.9 scan_SUSY_parameter_space(mcjzb,datajzb,jzb_cut,requireZ, peakerror, peakerrordata, ibin, njobs, jobnumber,systonly,effonly,shapeanalysis);
1016 buchmann 1.1 }
1017 buchmann 1.6 delete_any_cached_scans();// tidy up ...
1018 buchmann 1.1 }
1019