170 |
|
|
171 |
|
} |
172 |
|
|
173 |
< |
vector<float> compute_one_upper_limit(float mceff,float mcefferr, int ibin, string mcjzb, string plotfilename, bool doobserved) { |
173 |
> |
vector<float> compute_one_upper_limit(float mceff,float mcefferr, int ibin, string mcjzb, string plotfilename, bool doexpected) { |
174 |
|
float sigma95=-9.9,sigma95A=-9.9; |
175 |
|
/* |
176 |
|
USAGE OF ROOSTATS_CL95 |
216 |
|
dout << "Calling limit capsule instead of calling : CL95(" << luminosity << "," << lumiuncert*luminosity << "," << mceff << "," << mcefferr << "," << Npred[ibin] << "," << Nprederr[ibin] << "," << Nobs[ibin] << "," << false << "," << nuisancemodel<< ") " << endl; |
217 |
|
|
218 |
|
stringstream command; |
219 |
< |
command << PlottingSetup::cbafbasedir << "/DistributedModelCalculations/Limits/LimitCapsule.exec " << repname.str() << " " << luminosity << " " << luminosity*lumiuncert << " " << mceff << " " << mcefferr << " " << Npred[ibin] << " " << Nprederr[ibin] << " " << Nobs[ibin] << " " << -1 << " " << plotfilename; |
219 |
> |
command << PlottingSetup::cbafbasedir << "/DistributedModelCalculations/Limits/NewLimitCapsule.exec " << repname.str() << " " << luminosity << " " << luminosity*lumiuncert << " " << mceff << " " << mcefferr << " " << Npred[ibin] << " " << Nprederr[ibin] << " " << Nobs[ibin] << " " << -1 << " " << PlottingSetup::basedirectory << "/" << plotfilename << " " << doexpected; |
220 |
|
dout << command.str() << endl; |
221 |
|
|
222 |
|
int retval = 256; |
231 |
|
limres.readDroplet(repname.str()); |
232 |
|
dout << limres << endl; |
233 |
|
remove(repname.str().c_str()); |
234 |
< |
sigma95=limres.expected; |
234 |
> |
sigma95=limres.observed; |
235 |
|
|
236 |
|
|
237 |
|
///------------------------------------------ < /NEW > ---------------------------------------------------------- |
238 |
– |
//if(doobserved) nlimittoysused=nlimittoys; |
239 |
– |
// nlimittoysused=nlimittoys; |
240 |
– |
// dout << "Now calling : CL95(" << luminosity << "," << lumiuncert*luminosity << "," << mceff << "," << mcefferr << "," << Npred[ibin] << "," << Nprederr[ibin] << "," << Nobs[ibin] << "," << false << "," << nuisancemodel<< ") " << endl; |
241 |
– |
// sigma95 = CL95(luminosity, lumiuncert*luminosity, mceff, mcefferr, Npred[ibin], Nprederr[ibin], Nobs[ibin], false, nuisancemodel); |
242 |
– |
// sigma95 = CL95(luminosity, lumiuncert*luminosity, mceff, mcefferr, Npred[ibin], Nprederr[ibin], Nobs[ibin], false, nuisancemodel); |
243 |
– |
|
244 |
– |
/* dout << "Now calling : roostats_cl95(" << luminosity << "," << lumiuncert*luminosity << ","<<mceff <<","<<mcefferr<<","<<Npred[ibin]<<","<<Nprederr[ibin] << ",n=" << nlimittoysused << ",gauss=" << false << ",nuisanceModel="<<nuisancemodel<<",method="<<limitmethod<<",plotfilename="<<plotfilename<<",seed=0) " << endl; |
245 |
– |
/* dout << "Now calling : roostats_limit(" << luminosity << "," << lumiuncert*luminosity << ","<<mceff <<","<<mcefferr<<","<<Npred[ibin]<<","<<Nprederr[ibin] << ",n=" << nlimittoysused << ",gauss=" << false << ", nuisanceModel="<<nuisancemodel<<",method="<<limitmethod<<",plotfilename="<<plotfilename<<",seed=1) " << endl; |
246 |
– |
LimitResult limit = roostats_limit(luminosity,lumiuncert*luminosity,mceff,mcefferr,Npred[ibin],Nprederr[ibin],nlimittoysused,false,nuisancemodel,limitmethod,plotfilename,0); |
247 |
– |
dout << "Now interpreting and saving results ... " << endl; |
248 |
– |
vector<float> sigmas; |
249 |
– |
sigmas.push_back(limit.GetExpectedLimit());//expected |
250 |
– |
sigmas.push_back(limit.GetObservedLimit());//observed |
251 |
– |
//up to here for backward compatibility |
252 |
– |
sigmas.push_back(limit.GetOneSigmaHighRange());//expected, up |
253 |
– |
sigmas.push_back(limit.GetTwoSigmaHighRange());//expected, 2 up |
254 |
– |
sigmas.push_back(limit.GetOneSigmaLowRange());//expected, down |
255 |
– |
sigmas.push_back(limit.GetTwoSigmaLowRange());//expected, 2 down |
256 |
– |
*/ |
257 |
– |
// float limit = roostats_cl95(luminosity,lumiuncert*luminosity,mceff,mcefferr,Npred[ibin],Nprederr[ibin],nlimittoysused,false,nuisancemodel,limitmethod,plotfilename,0); |
258 |
– |
if(doobserved) { |
259 |
– |
write_warning(__FUNCTION__,"OBSERVED LIMITS HAVE BEEN SWITCHED OFF TEMPORARILY"); |
260 |
– |
// dout << "Now calling : CLA(" << luminosity << "," << lumiuncert*luminosity << "," << mceff << "," << mcefferr << "," << Npred[ibin] << "," << Nprederr[ibin] << "," << nuisancemodel<< ") " << endl; |
261 |
– |
// sigma95A = CLA(luminosity, lumiuncert*luminosity, mceff, mcefferr, Npred[ibin], Nprederr[ibin], nuisancemodel); |
262 |
– |
} |
263 |
– |
// vector<float> sigmas; |
264 |
– |
// sigmas.push_back(limit); |
238 |
|
vector<float> sigmas; |
239 |
|
sigmas.push_back(sigma95); |
240 |
< |
sigmas.push_back(sigma95A); |
240 |
> |
if(doexpected) { |
241 |
> |
sigmas.push_back(limres.expected); |
242 |
> |
sigmas.push_back(limres.upper68); |
243 |
> |
sigmas.push_back(limres.lower68); |
244 |
> |
sigmas.push_back(limres.upper95); |
245 |
> |
sigmas.push_back(limres.lower95); |
246 |
> |
} |
247 |
> |
|
248 |
|
return sigmas; |
249 |
|
|
250 |
|
|
251 |
< |
} |
272 |
< |
write_warning(__FUNCTION__,"STILL MISSING SIGMAS, LIMITS, EVERYTHING ..."); |
251 |
> |
}//end of mc efficiency is ok |
252 |
|
} |
253 |
|
|
254 |
< |
void compute_upper_limits_from_counting_experiment(vector<vector<float> > uncertainties,vector<float> jzbcuts, string mcjzb, bool doobserved) { |
254 |
> |
void compute_upper_limits_from_counting_experiment(vector<vector<float> > uncertainties,vector<float> jzbcuts, string mcjzb, bool doexpected) { |
255 |
|
dout << "Doing counting experiment ... " << endl; |
256 |
|
vector<vector<string> > limits; |
257 |
|
vector<vector<float> > vlimits; |
276 |
|
// float expected=observed/luminosity; |
277 |
|
string plotfilename=(string)(TString(signalsamples.collection[isample].samplename)+TString("___JZB_geq_")+TString(any2string(JZBcutat))+TString(".png")); |
278 |
|
dout << "Sample: " << signalsamples.collection[isample].samplename << ", JZB>"<<JZBcutat<< " : " << mceff << " +/- " << staterr << " (stat) +/- " << systerr << " (syst) --> toterr = " << toterr << endl; |
279 |
< |
vector<float> sigmas = compute_one_upper_limit(mceff,toterr,ibin,mcjzb,plotfilename,doobserved); |
279 |
> |
vector<float> sigmas = compute_one_upper_limit(mceff,toterr,ibin,mcjzb,plotfilename,doexpected); |
280 |
|
|
281 |
< |
if(doobserved) { |
281 |
> |
if(doexpected) { |
282 |
|
// rows.push_back(any2string(sigmas[0])+";"+any2string(sigmas[1])+";"+"("+any2string(expected)+")"); |
283 |
|
rows.push_back(any2string(sigmas[0])+";"+any2string(sigmas[1])+";"+"("+any2string(signalsamples.collection[isample].xs)+")"); |
284 |
|
vrows.push_back(sigmas[0]); |
312 |
|
dout << endl; |
313 |
|
} |
314 |
|
|
315 |
< |
if(!doobserved) { |
315 |
> |
if(!doexpected) { |
316 |
|
dout << endl << endl << "LIMITS: (Tex)" << endl; |
317 |
|
tout << "\\begin{table}[hbtp]" << endl; |
318 |
|
tout << "\\renewcommand{\\arraystretch}{1.3}" << endl; |
351 |
|
} |
352 |
|
dout << endl; |
353 |
|
} |
375 |
– |
|
376 |
– |
write_warning(__FUNCTION__,"Still need to update the script"); |
354 |
|
} |
355 |
|
|
356 |
|
|