ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/AnalysisFramework/Plotting/Modules/LimitCalculation.C
(Generate patch)

Comparing UserCode/cbrown/AnalysisFramework/Plotting/Modules/LimitCalculation.C (file contents):
Revision 1.21 by buchmann, Wed Sep 14 12:02:01 2011 UTC vs.
Revision 1.25 by buchmann, Mon Sep 26 15:48:20 2011 UTC

# Line 170 | Line 170 | ratio_binning.push_back(80);
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
# Line 216 | Line 216 | 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;
# Line 231 | Line 231 | USAGE OF ROOSTATS_CL95
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;
# Line 297 | Line 276 | void compute_upper_limits_from_counting_
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]);
# Line 333 | Line 312 | void compute_upper_limits_from_counting_
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;
# Line 372 | Line 351 | void compute_upper_limits_from_counting_
351      }
352      dout << endl;
353    }
375  
376  write_warning(__FUNCTION__,"Still need to update the script");
354   }
355  
356  

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines