| 24 |
|
|
| 25 |
|
using namespace PlottingSetup; |
| 26 |
|
|
| 27 |
+ |
namespace SUSYScanSpace { |
| 28 |
+ |
vector<string> loaded_files; |
| 29 |
+ |
int SUSYscantype; |
| 30 |
+ |
} |
| 31 |
+ |
|
| 32 |
|
void susy_scan_axis_labeling(TH2F *histo) { |
| 33 |
|
histo->GetXaxis()->SetTitle("m_{#Chi_{2}^{0}}-m_{LSP}"); |
| 34 |
|
histo->GetXaxis()->CenterTitle(); |
| 72 |
|
} |
| 73 |
|
|
| 74 |
|
bool delete_any_cached_scans() { |
| 75 |
+ |
|
| 76 |
+ |
|
| 77 |
+ |
/* |
| 78 |
|
//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 |
< |
if((Contains(hostname,"t3ui")||Contains(hostname,"t3wn"))) return false; |
| 84 |
< |
else { |
| 83 |
> |
if(!((Contains(hostname,"t3ui")||Contains(hostname,"t3wn")))) |
| 84 |
> |
{ |
| 85 |
|
vector<string> all_files; |
| 86 |
|
char currentpath[1024]; |
| 87 |
|
TString directory=TString(getcwd(currentpath,1024)); |
| 94 |
|
gSystem->Exec(((string)"rm "+all_files[ifile]).c_str()); |
| 95 |
|
} |
| 96 |
|
} |
| 97 |
< |
return true; |
| 97 |
> |
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 |
+ |
return true; |
| 108 |
|
} |
| 109 |
|
|
| 110 |
|
bool initialized_t2=false; |
| 111 |
|
|
| 112 |
|
int srmcpretries=0; |
| 113 |
|
|
| 114 |
< |
void load_scan_sample(int a, int b, int &scanfileindex,int scantype,bool isretry=false) { |
| 115 |
< |
/* |
| 114 |
> |
void load_scan_sample(int a, int b, int &scanfileindex,int scantype,int mglu, int mlsp, bool isretry=false) { |
| 115 |
> |
|
| 116 |
|
// There is no need to define your sample here. That is done in Setup.C where you define the loading directory! |
| 117 |
|
|
| 118 |
< |
dout << "Going to load file with Xzone=" << a << " and Yzone=" << b << endl; |
| 118 |
> |
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 |
|
stringstream filetoload; |
| 120 |
|
char hostname[1023]; |
| 121 |
|
gethostname(hostname,1023); |
| 122 |
< |
|
| 123 |
< |
if((Contains(hostname,"t3ui")||Contains(hostname,"t3wn"))) { // local case |
| 122 |
> |
string samplename; |
| 123 |
> |
if((Contains(hostname,"t3ui")||Contains(hostname,"t3wn"))&&scantype!=mSUGRA) { // local case |
| 124 |
|
filetoload << "/shome/buchmann/ntuples/"<<PlottingSetup::ScanSampleDirectory; |
| 125 |
< |
if(scantype==mSUGRA) filetoload << "/mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root"; |
| 126 |
< |
if(scantype==SMS) filetoload << "/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root"; |
| 127 |
< |
if(scantype==GMSB) filetoload << "/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root"; |
| 125 |
> |
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 |
|
if(scansample.collection.size()<1||!Contains(((scansample.collection)[(scansample.collection).size()-1]).filename,"_"+any2string(a)+"_"+any2string(b))) { |
| 141 |
|
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 |
< |
if((scansample.collection).size()>1) { |
| 142 |
> |
if((scansample.collection).size()>=1) { |
| 143 |
|
scansample.RemoveLastSample(); |
| 144 |
|
} |
| 145 |
|
scanfileindex=(scansample.collection).size(); |
| 146 |
|
//New: Loading file when necessary, not before (avoiding high memory usage and startup times) |
| 147 |
< |
if(scanfileindex!=0) scansample.AddSample(filetoload.str(),"scansample",1,1,false,true,scanfileindex,kRed); |
| 147 |
> |
|
| 148 |
> |
scansample.AddSample(filetoload.str(),samplename,1,1,false,true,scanfileindex,kRed); |
| 149 |
> |
dout << "Just added the following file: " << filetoload.str() << endl; |
| 150 |
|
} else { |
| 151 |
|
dout << "Last sample is the same as the current one. Recycling it." << endl; |
| 152 |
|
scanfileindex=(scansample.collection).size()-1; |
| 154 |
|
dout << " Going to use the following file: " << filetoload.str() << endl; |
| 155 |
|
} // end of t3 case |
| 156 |
|
else { |
| 157 |
+ |
//CRAB case |
| 158 |
|
write_info(__FUNCTION__,"Hello, CRAB! Might need to load files ..."); |
| 159 |
|
PlottingSetup::limitpatience=PlottingSetup::limitpatienceCRAB; |
| 160 |
|
stringstream copyfile; |
| 161 |
+ |
gSystem->Exec("mkdir -p local_storage"); |
| 162 |
|
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 |
< |
if(scantype==mSUGRA) copyfile << "/mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root " << "mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root "; |
| 164 |
< |
if(scantype==SMS) copyfile << "/SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root " << "SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root "; |
| 165 |
< |
if(scantype==GMSB) copyfile << "/GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root " << "GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root "; |
| 163 |
> |
if(scantype==mSUGRA) { |
| 164 |
> |
copyfile<< "/mSUGRA_M0_" << mglu << "__M12_" << mlsp << ".root local_storage/mSUGRA_M0_" << mglu << "__M12_" << mlsp << ".root "; |
| 165 |
> |
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 |
|
stringstream newfilename; |
| 176 |
< |
if(scantype==mSUGRA) newfilename << "mSUGRA_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root"; |
| 177 |
< |
if(scantype==SMS) newfilename << "SMS_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root"; |
| 178 |
< |
if(scantype==GMSB) newfilename << "GMSB_clean_splitup_" << any2string(a) << "_" << any2string(b) << ".root"; |
| 176 |
> |
// 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 |
> |
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 |
> |
|
| 181 |
> |
SUSYScanSpace::SUSYscantype=scantype; |
| 182 |
|
|
| 183 |
|
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 |
< |
if((scansample.collection).size()>1) { |
| 185 |
> |
while((scansample.collection).size()>0) { |
| 186 |
|
scansample.RemoveLastSample(); |
| 187 |
|
} |
| 188 |
+ |
cout << "New scanfileindex stems from scansample size: " << (scansample.collection).size() << endl; |
| 189 |
|
scanfileindex=(scansample.collection).size(); |
| 190 |
|
//New: Loading file when necessary, not before (avoiding high memory usage and startup times) |
| 191 |
< |
if(scanfileindex!=0||isretry) { |
| 192 |
< |
delete_any_cached_scans(); |
| 193 |
< |
dout << "Going to download the scan file with the following copy command: " << copyfile.str() << endl; |
| 194 |
< |
gSystem->Exec(copyfile.str().c_str()); |
| 191 |
> |
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 |
> |
SUSYScanSpace::loaded_files.push_back(newfilename.str()); // adding it to the list of loaded files (will be deleted once we finish running |
| 196 |
> |
} |
| 197 |
|
if(doesROOTFileExist(newfilename.str())) { |
| 198 |
|
initialized_t2=true; |
| 199 |
< |
scansample.AddSample(newfilename.str(),"scansample",1,1,false,true,scanfileindex,kRed); |
| 199 |
> |
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 |
|
} else { |
| 203 |
|
srmcpretries++; |
| 204 |
|
if(srmcpretries<5) { |
| 205 |
|
dout << "The file could not be loaded correctly - retrying!" << endl; |
| 206 |
|
sleep(5); |
| 207 |
< |
load_scan_sample(a,b,scanfileindex,scantype,true); |
| 207 |
> |
load_scan_sample(a,b,scanfileindex,scantype,mglu,mlsp,true); |
| 208 |
|
} else { |
| 209 |
|
dout << "Have tried 5 times to load this sample. Giving up now and failing the program execution" << endl; |
| 210 |
< |
assert(0); |
| 210 |
> |
assert(srmcpretries<5); |
| 211 |
|
dout << "The command " << copyfile.str() << " failed to copy the file to the local storage (saving it as " << newfilename.str() << ")" << endl; |
| 212 |
|
} |
| 213 |
|
} |
| 218 |
|
} |
| 219 |
|
|
| 220 |
|
} |
| 221 |
< |
*/ |
| 221 |
> |
|
| 222 |
|
|
| 223 |
|
// WATCH OUT THIS LINE NEEDS TO BE REMOVED |
| 224 |
|
scanfileindex=0; |
| 255 |
|
while ((key = (TKey*)nextkey())) |
| 256 |
|
{ |
| 257 |
|
TObject *obj = key->ReadObj(); |
| 258 |
< |
if(Contains((string)(obj->GetName()),"mSUGRA")) scantype=mSUGRA; |
| 258 |
> |
if(Contains((string)(obj->GetName()),"SUGRA")) scantype=mSUGRA; |
| 259 |
|
if(Contains((string)(obj->GetName()),"GMSB")) scantype=GMSB; |
| 260 |
|
} |
| 261 |
|
|
| 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); |
| 263 |
– |
TH2F *xsmap = new TH2F((prefix+"crosssectionmap"+any2string(jzb_cut[ibin])).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,(int)((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep); |
| 264 |
– |
TH2F *totxsmap = new TH2F((prefix+"absolutecrosssectionmap"+any2string(jzb_cut[ibin])).c_str(), "",(int)((mgluend-mglustart)/mglustep+1),mglustart-0.5*mglustep,mgluend+0.5*mglustep,(int)((mLSPend-mLSPstart)/mLSPstep+1),mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep); |
| 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 |
+ |
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 |
|
|
| 316 |
|
|
| 317 |
|
if(!fullerr || !mceff || !NEvents) { |
| 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 |
< |
if(strength>0.5&&strength<2) sigmas=compute_one_upper_limit(currmceff,currtoterr,ibin,mcjzb,plotfilename,true, flipped); // asymptotic limit first |
| 353 |
< |
else sigmas=asigmas; |
| 352 |
> |
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 |
|
exclmap->SetBinContent(GlobalBin,strength); |
| 306 |
– |
xsmap->SetBinContent(GlobalBin,XSmap->GetBinContent(GlobalBin)); |
| 361 |
|
} |
| 362 |
|
|
| 309 |
– |
|
| 310 |
– |
|
| 311 |
– |
|
| 363 |
|
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]); |
| 374 |
|
} |
| 375 |
|
|
| 376 |
|
prepare_scan_axis(limitmap,scantype); |
| 377 |
+ |
gSystem->Exec("mkdir -p output"); |
| 378 |
|
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 |
+ |
asymptoticmap->Write(); |
| 383 |
|
if(doexpected) { |
| 384 |
|
explimitmap->Write(); |
| 385 |
|
exp1plimitmap->Write(); |
| 387 |
|
exp2plimitmap->Write(); |
| 388 |
|
exp2mlimitmap->Write(); |
| 389 |
|
} |
| 390 |
< |
if(scantype==mSUGRA) { |
| 338 |
< |
exclmap->Write(); |
| 339 |
< |
xsmap->Write(); |
| 340 |
< |
totxsmap->Write(); |
| 341 |
< |
} |
| 390 |
> |
if(scantype==mSUGRA) exclmap->Write(); |
| 391 |
|
outputfile->Close(); |
| 392 |
|
delete limcanvas; |
| 393 |
|
} |
| 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 |
< |
if(Contains((scansample.collection)[0].samplename,"mSUGRA")) scantype=mSUGRA; |
| 432 |
> |
if(Contains((scansample.collection)[0].samplename,"SUGRA")) scantype=mSUGRA; |
| 433 |
|
if(Contains((scansample.collection)[0].samplename,"GMSB")) scantype=GMSB; |
| 434 |
|
|
| 435 |
|
if(scantype==mSUGRA) { |
| 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 |
< |
load_scan_sample(a,b,scanfileindex,scantype); |
| 567 |
> |
load_scan_sample(a,b,scanfileindex,scantype,mglu,mlsp); |
| 568 |
|
|
| 569 |
|
clock_t start,finish; |
| 570 |
|
start = clock(); // starting the clock to measure how long the computation takes! |
| 627 |
|
} |
| 628 |
|
|
| 629 |
|
|
| 630 |
+ |
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 |
|
if(nevents!=0&&(efficiencyonly||systematicsonly)) { |
| 637 |
|
Value effwosigcont = MCefficiency((scansample.collection)[scanfileindex].events,result,resulterr,flipped,mcjzb,requireZ,(int)nevents,scantype,xsec,addcut.str(),-1); |
| 638 |
|
if(result<0&&allowflipping) { |
| 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 |
< |
|
| 686 |
> |
cout << "Going to store: total error is: " << toterr << " and the relative tot error, i.e. total error / mc efficiency is " << toterr/mceff << endl; |
| 687 |
|
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; |
| 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 |
+ |
delete_any_cached_scans();// tidy up ... |
| 967 |
|
} |
| 968 |
|
|
