ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/SUSYScan.C
Revision: 1.6
Committed: Wed Mar 21 22:00:49 2012 UTC (13 years, 1 month ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.5: +26 -5 lines
Log Message: 
Made scans more economic: files are on SE, point we want to consider gets downloaded, analyzed and the local copy is deleted; that way we significantly reduce the local storage required for the scans

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