| 183 |
|
|
| 184 |
|
} |
| 185 |
|
|
| 186 |
< |
vector<float> compute_one_upper_limit(float mceff,float mcefferr, int ibin, string mcjzb, string plotfilename, bool doexpected) { |
| 186 |
> |
vector<float> compute_one_upper_limit(float mceff,float mcefferr, int ibin, string mcjzb, string plotfilename, bool doexpected, int flipped) { |
| 187 |
|
float sigma95=-9.9,sigma95A=-9.9; |
| 188 |
|
/* |
| 189 |
|
USAGE OF ROOSTATS_CL95 |
| 226 |
|
- JZB cut [9] |
| 227 |
|
- plot name [10]*/ |
| 228 |
|
|
| 229 |
< |
dout << "Calling limit capsule instead of calling : CL95(" << luminosity << "," << lumiuncert*luminosity << "," << mceff << "," << mcefferr << "," << Npred[ibin] << "," << Nprederr[ibin] << "," << Nobs[ibin] << "," << false << "," << nuisancemodel<< ") " << endl; |
| 229 |
> |
if(flipped==0) dout << "Calling limit capsule instead of calling : CL95(" << luminosity << "," << lumiuncert*luminosity << "," << mceff << "," << mcefferr << "," << Npred[ibin] << "," << Nprederr[ibin] << "," << Nobs[ibin] << "," << false << "," << nuisancemodel<< ") " << endl; |
| 230 |
> |
if(flipped>0) dout << "Calling limit capsule instead of calling : CL95(" << luminosity << "," << lumiuncert*luminosity << "," << mceff << "," << mcefferr << "," << flippedNpred[ibin] << "," << flippedNprederr[ibin] << "," << flippedNobs[ibin] << "," << false << "," << nuisancemodel<< ") " << endl; |
| 231 |
|
|
| 232 |
|
stringstream command; |
| 233 |
< |
command << PlottingSetup::cbafbasedir << "/DistributedModelCalculations/Limits/TimedLimitCapsule.exec " << repname.str() << " " << luminosity << " " << luminosity*lumiuncert << " " << mceff << " " << mcefferr << " " << Npred[ibin] << " " << Nprederr[ibin] << " " << Nobs[ibin] << " " << -1 << " " << PlottingSetup::basedirectory << "/" << plotfilename << " " << doexpected; |
| 233 |
> |
if(flipped==0) command << PlottingSetup::cbafbasedir << "/DistributedModelCalculations/Limits/TimedLimitCapsule.exec " << repname.str() << " " << luminosity << " " << luminosity*lumiuncert << " " << mceff << " " << mcefferr << " " << Npred[ibin] << " " << Nprederr[ibin] << " " << Nobs[ibin] << " " << -1 << " " << PlottingSetup::basedirectory << "/" << plotfilename << " " << doexpected; |
| 234 |
> |
if(flipped>0) command << PlottingSetup::cbafbasedir << "/DistributedModelCalculations/Limits/TimedLimitCapsule.exec " << repname.str() << " " << luminosity << " " << luminosity*lumiuncert << " " << mceff << " " << mcefferr << " " << flippedNpred[ibin] << " " << flippedNprederr[ibin] << " " << flippedNobs[ibin] << " " << -1 << " " << PlottingSetup::basedirectory << "/" << plotfilename << " " << doexpected; |
| 235 |
|
dout << command.str() << endl; |
| 236 |
|
|
| 237 |
|
int retval = 256; |
| 238 |
|
int attempts=0; |
| 239 |
|
while(!(retval==0||attempts>=3)) {//try up to 3 times |
| 240 |
|
attempts++; |
| 241 |
< |
dout << "Starting limit calculation (LimitCapsule) now : Attempt " << attempts << endl; |
| 241 |
> |
dout << "Starting limit calculation (TimedLimitCapsule) now : Attempt " << attempts << endl; |
| 242 |
|
retval=gSystem->Exec(command.str().c_str()); |
| 243 |
|
} |
| 244 |
< |
|
| 244 |
> |
char hostname[1023]; |
| 245 |
> |
gethostname(hostname,1023); |
| 246 |
> |
if((!((Contains(hostname,"t3ui")||Contains(hostname,"t3wn"))))&&retval==256) { |
| 247 |
> |
//running via CRAB and encountered the same problem too often: place a problem file to mark this problem! |
| 248 |
> |
stringstream markproblem; |
| 249 |
> |
markproblem << "touch " << PlottingSetup::basedirectory << "/exchange/problemswhilesettinglimits.txt"; |
| 250 |
> |
gSystem->Exec(markproblem.str().c_str()); |
| 251 |
> |
} |
| 252 |
|
LimitDroplet limres; |
| 253 |
|
limres.readDroplet(repname.str()); |
| 254 |
|
dout << limres << endl; |
| 273 |
|
}//end of mc efficiency is ok |
| 274 |
|
} |
| 275 |
|
|
| 276 |
< |
void compute_upper_limits_from_counting_experiment(vector<vector<float> > uncertainties,vector<float> jzbcuts, string mcjzb, bool doexpected) { |
| 276 |
> |
void compute_upper_limits_from_counting_experiment(vector<vector<float> > uncertainties,vector<float> jzbcuts, string mcjzb, bool doexpected, int flipped) { |
| 277 |
|
dout << "Doing counting experiment ... " << endl; |
| 278 |
|
vector<vector<string> > limits; |
| 279 |
|
vector<vector<float> > vlimits; |
| 293 |
|
float toterr =uncertainties[isample*jzbcuts.size()+ibin][4]; |
| 294 |
|
float observed,observederr,null,result; |
| 295 |
|
|
| 287 |
– |
// fill_result_histos(observed,observederr, null,null,null,null,null,null,null,mcjzb,JZBcutat,14000,(int)5,result,(signalsamples.FindSample(signalsamples.collection[isample].filename)),signalsamples); |
| 288 |
– |
// observed-=result;//this is the actual excess we see! |
| 289 |
– |
// float expected=observed/luminosity; |
| 296 |
|
string plotfilename=(string)(TString(signalsamples.collection[isample].samplename)+TString("___JZB_geq_")+TString(any2string(JZBcutat))+TString(".png")); |
| 297 |
|
dout << "Sample: " << signalsamples.collection[isample].samplename << ", JZB>"<<JZBcutat<< " : " << mceff << " +/- " << staterr << " (stat) +/- " << systerr << " (syst) --> toterr = " << toterr << endl; |
| 298 |
< |
vector<float> sigmas = compute_one_upper_limit(mceff,toterr,ibin,mcjzb,plotfilename,doexpected); |
| 299 |
< |
|
| 298 |
> |
vector<float> sigmas = compute_one_upper_limit(mceff,toterr,ibin,mcjzb,plotfilename,doexpected,flipped); |
| 299 |
> |
|
| 300 |
> |
tripple LibraryUpperLimits; |
| 301 |
> |
LibraryUpperLimits.name=signalsamples.collection[isample].samplename; |
| 302 |
> |
LibraryUpperLimits.first=mceff*signalsamples.collection[isample].xs * PlottingSetup::luminosity; |
| 303 |
> |
LibraryUpperLimits.second=staterr*signalsamples.collection[isample].xs * PlottingSetup::luminosity; |
| 304 |
> |
int resultindex=PlottingSetup::allresults.Find(jzbcuts[ibin]); |
| 305 |
> |
(allresults.predictions[resultindex]).SignalYield.push_back(LibraryUpperLimits); |
| 306 |
> |
|
| 307 |
|
if(doexpected) { |
| 295 |
– |
// rows.push_back(any2string(sigmas[0])+";"+any2string(sigmas[1])+";"+"("+any2string(expected)+")"); |
| 308 |
|
rows.push_back(any2string(sigmas[0])+";"+any2string(sigmas[1])+";"+"("+any2string(signalsamples.collection[isample].xs)+")"); |
| 309 |
|
vrows.push_back(sigmas[0]); |
| 310 |
|
vrows.push_back(sigmas[1]); |
| 299 |
– |
// vrows.push_back(expected); |
| 311 |
|
vrows.push_back(signalsamples.collection[isample].xs); |
| 312 |
|
} |
| 313 |
|
else { |
| 303 |
– |
// rows.push_back(any2string(sigmas[0])+"("+any2string(expected)+")"); |
| 314 |
|
rows.push_back(any2string(sigmas[0])); |
| 315 |
|
vrows.push_back(sigmas[0]); |
| 316 |
|
vrows.push_back(signalsamples.collection[isample].xs); |
| 307 |
– |
// vrows.push_back(expected); |
| 317 |
|
} |
| 318 |
|
}//end of bin loop |
| 319 |
|
limits.push_back(rows); |
| 320 |
|
vlimits.push_back(vrows); |
| 321 |
|
}//end of sample loop |
| 322 |
+ |
|
| 323 |
+ |
dout << endl << endl << "_______________________________________________________________________________________" << endl; |
| 324 |
+ |
dout << "Going to store upper limit on event yield in result library: " << endl; |
| 325 |
+ |
for(int ibin=0;ibin<jzbcuts.size();ibin++) { |
| 326 |
+ |
int resultindex=PlottingSetup::allresults.Find(jzbcuts[ibin]); |
| 327 |
+ |
vector<float> Normsigmas = compute_one_upper_limit(1.0,0.0, resultindex, mcjzb, "UPPERLIMIT", false, 0); |
| 328 |
+ |
(allresults.predictions[resultindex]).UpperLimit=Normsigmas[0]*PlottingSetup::luminosity; |
| 329 |
+ |
} |
| 330 |
+ |
dout << endl << "_______________________________________________________________________________________" << endl; |
| 331 |
+ |
|
| 332 |
|
dout << endl << endl << endl << "_________________________________________________________________________________________________" << endl << endl; |
| 333 |
|
dout << endl << endl << "PAS table 3: (notation: limit [95%CL])" << endl << endl; |
| 334 |
|
dout << "\t"; |
| 383 |
|
} |
| 384 |
|
dout << endl; |
| 385 |
|
} |
| 386 |
+ |
allresults.Print(); |
| 387 |
+ |
|
| 388 |
+ |
|
| 389 |
|
} |
| 390 |
|
|
| 391 |
|
|
