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 |
|
|
|