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.17 by buchmann, Thu Sep 1 16:23:46 2011 UTC vs.
Revision 1.27 by buchmann, Mon Oct 24 15:26:14 2011 UTC

# Line 1 | Line 1
1 + /****
2 +
3 + Off peak status (RestrictToMassPeak) :
4 +
5 + x  Necessary adaptations identified
6 + x  Started working on necessary adaptations
7 + x  Necessary adaptations implemented
8 + x  Necessary adaptations tested
9 +
10 + DONE!
11 +
12 +
13 + ****/
14   #include <iostream>
15   #include <vector>
16   #include <sys/stat.h>
# Line 15 | Line 28
28   #include <TF1.h>
29   #include <TSQLResult.h>
30   #include <TProfile.h>
31 + #include <TSystem.h>
32 + #include "LimitDroplet.C"
33  
34   //#include "TTbar_stuff.C"
35   using namespace std;
# Line 168 | Line 183 | ratio_binning.push_back(80);
183    
184   }
185  
186 < vector<float> compute_one_upper_limit(float mceff,float mcefferr, int ibin, string mcjzb, string plotfilename, bool doobserved) {
186 > vector<float> compute_one_upper_limit(float mceff,float mcefferr, int ibin, string mcjzb, string plotfilename, bool doexpected) {
187    float sigma95=-9.9,sigma95A=-9.9;
173 –  int nuisancemodel=0;
188   /*
189   USAGE OF ROOSTATS_CL95
190   " Double_t             limit = roostats_cl95(ilum, slum, eff, seff, bck, sbck, n, gauss = false, nuisanceModel, method, plotFileName, seed); \n"
# Line 194 | Line 208 | USAGE OF ROOSTATS_CL95
208      return sigmas;
209    } else {
210      int nlimittoysused=1;
211 <    if(doobserved) nlimittoysused=nlimittoys;
212 < /*  dout << "Now calling : CL95(" << luminosity << "," <<  lumiuncert*luminosity << "," << mceff << "," << mcefferr << "," << Npred[ibin] << "," << Nprederr[ibin] << "," << Nobs[ibin] << "," << false << "," << nuisancemodel<< ") " << endl;
213 <  sigma95 = CL95(luminosity, lumiuncert*luminosity, mceff, mcefferr, Npred[ibin], Nprederr[ibin], Nobs[ibin], false, nuisancemodel);
214 <  */
215 <  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;
216 <  LimitResult limit = roostats_limit(luminosity,lumiuncert*luminosity,mceff,mcefferr,Npred[ibin],Nprederr[ibin],nlimittoysused,false,nuisancemodel,limitmethod,plotfilename,1);
217 <  
218 <  vector<float> sigmas;
219 <  sigmas.push_back(limit.GetExpectedLimit());//expected
220 <  sigmas.push_back(limit.GetObservedLimit());//observed
221 <  //up to here for backward compatibility
222 <  sigmas.push_back(limit.GetOneSigmaHighRange());//expected, up
223 <  sigmas.push_back(limit.GetTwoSigmaHighRange());//expected, 2 up
224 <  sigmas.push_back(limit.GetOneSigmaLowRange());//expected, down
225 <  sigmas.push_back(limit.GetTwoSigmaLowRange());//expected, 2 down
226 < //  if(doobserved) {
227 < //    dout << "Now calling : CLA(" << luminosity << "," <<  lumiuncert*luminosity << "," << mceff << "," << mcefferr << "," << Npred[ibin] << "," << Nprederr[ibin] << "," << nuisancemodel<< ") " << endl;
228 < //    sigma95A = CLA(luminosity, lumiuncert*luminosity, mceff, mcefferr, Npred[ibin], Nprederr[ibin], nuisancemodel);
229 < //  }
211 >    
212 >    ///------------------------------------------ < NEW > ----------------------------------------------------------
213 >    
214 >    int secondssince1970=time(NULL);
215 >    stringstream repname;
216 >    repname << PlottingSetup::cbafbasedir << "/exchange/report_" << secondssince1970 << "_"<<plotfilename<< "__"<< ".txt";
217 >    
218 >      /* - report filename [1]
219 >         - luminosity [2]
220 >         - lumi uncert [3]
221 >         - MC efficiency [4]
222 >         - MC efficiency error [5]
223 >         - Npred [6]
224 >         - Nprederr [7]
225 >         - Nobs [8]
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;
230 >    
231 >    stringstream command;
232 >    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 >    dout << command.str() << endl;
234 >    
235 >    int retval = 256;
236 >    int attempts=0;
237 >    while(!(retval==0||attempts>=3)) {//try up to 3 times
238 >        attempts++;
239 >        dout << "Starting limit calculation (LimitCapsule) now : Attempt " << attempts << endl;
240 >        retval=gSystem->Exec(command.str().c_str());
241 >    }
242 >    
243 >    LimitDroplet limres;
244 >    limres.readDroplet(repname.str());
245 >    dout << limres << endl;
246 >    remove(repname.str().c_str());
247 >    sigma95=limres.observed;
248  
249 < /*  vector<float> sigmas;
249 >    
250 >    ///------------------------------------------ < /NEW > ----------------------------------------------------------
251 >  vector<float> sigmas;
252    sigmas.push_back(sigma95);
253 <  sigmas.push_back(sigma95A);*/
253 >  if(doexpected) {
254 >    sigmas.push_back(limres.expected);
255 >    sigmas.push_back(limres.upper68);
256 >    sigmas.push_back(limres.lower68);
257 >    sigmas.push_back(limres.upper95);
258 >    sigmas.push_back(limres.lower95);
259 >  }
260 >  
261    return sigmas;
262    
263  
264 <  }
224 <  write_warning(__FUNCTION__,"STILL MISSING SIGMAS, LIMITS, EVERYTHING ...");
264 >  }//end of mc efficiency is ok
265   }
266  
267 < void compute_upper_limits_from_counting_experiment(vector<vector<float> > uncertainties,vector<float> jzbcuts, string mcjzb, bool doobserved) {
267 > void compute_upper_limits_from_counting_experiment(vector<vector<float> > uncertainties,vector<float> jzbcuts, string mcjzb, bool doexpected) {
268    dout << "Doing counting experiment ... " << endl;
269    vector<vector<string> > limits;
270    vector<vector<float> > vlimits;
# Line 249 | Line 289 | void compute_upper_limits_from_counting_
289   //      float expected=observed/luminosity;
290        string plotfilename=(string)(TString(signalsamples.collection[isample].samplename)+TString("___JZB_geq_")+TString(any2string(JZBcutat))+TString(".png"));
291        dout << "Sample: " << signalsamples.collection[isample].samplename << ", JZB>"<<JZBcutat<< " : " << mceff << " +/- " << staterr << " (stat) +/- " << systerr << " (syst) --> toterr = " << toterr << endl;
292 <      vector<float> sigmas = compute_one_upper_limit(mceff,toterr,ibin,mcjzb,plotfilename,doobserved);
292 >      vector<float> sigmas = compute_one_upper_limit(mceff,toterr,ibin,mcjzb,plotfilename,doexpected);
293        
294 <      if(doobserved) {
294 >      if(doexpected) {
295   //      rows.push_back(any2string(sigmas[0])+";"+any2string(sigmas[1])+";"+"("+any2string(expected)+")");
296          rows.push_back(any2string(sigmas[0])+";"+any2string(sigmas[1])+";"+"("+any2string(signalsamples.collection[isample].xs)+")");
297          vrows.push_back(sigmas[0]);
# Line 285 | Line 325 | void compute_upper_limits_from_counting_
325      dout << endl;
326    }
327    
328 <  if(!doobserved) {
328 >  if(!doexpected) {
329      dout << endl << endl << "LIMITS: (Tex)" << endl;
330      tout << "\\begin{table}[hbtp]" << endl;
331 <    tout << "\\renewcommand{\arraystretch}{1.3}" << endl;
331 >    tout << "\\renewcommand{\\arraystretch}{1.3}" << endl;
332      tout << "\\begin{center}" << endl;
333      tout << "\\caption{Observed upper limits on the cross section of different LM benchmark points " << (ConsiderSignalContaminationForLimits?"  (accounting for signal contamination)":"  (not accounting for signal contamination)") << "}\\label{tab:lmresults}" << endl;
334      tout << "" << endl;
# Line 324 | Line 364 | void compute_upper_limits_from_counting_
364      }
365      dout << endl;
366    }
327 –  
328 –  write_warning(__FUNCTION__,"Still need to update the script");
367   }
368  
369  
# Line 378 | Line 416 | void limit_shapes_for_systematic_effect(
416    TH1F *ZOSSFN = allsamples.Draw("ZOSSFN","-"+datajzb,binning, "JZB4limits", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,dataormc,luminosity);
417    TH1F *ZOSOFN = allsamples.Draw("ZOSOFN","-"+datajzb,binning, "JZB4limits", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,dataormc,luminosity);
418    
419 <  TH1F *SBOSSFP = allsamples.Draw("SBOSSFP",datajzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
420 <  TH1F *SBOSOFP = allsamples.Draw("SBOSOFP",datajzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
421 <  TH1F *SBOSSFN = allsamples.Draw("SBOSSFN","-"+datajzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
422 <  TH1F *SBOSOFN = allsamples.Draw("SBOSOFN","-"+datajzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
419 >  TH1F *SBOSSFP;
420 >  TH1F *SBOSOFP;
421 >  TH1F *SBOSSFN;
422 >  TH1F *SBOSOFN;
423    
424    TH1F *LZOSSFP = allsamples.Draw("LZOSSFP",mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("LM4"));
425    TH1F *LZOSOFP = allsamples.Draw("LZOSOFP",mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("LM4"));
426    TH1F *LZOSSFN = allsamples.Draw("LZOSSFN","-"+mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("LM4"));
427    TH1F *LZOSOFN = allsamples.Draw("LZOSOFN","-"+mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("LM4"));
428    
429 <  TH1F *LSBOSSFP = allsamples.Draw("LSBOSSFP",mcjzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
430 <  TH1F *LSBOSOFP = allsamples.Draw("LSBOSOFP",mcjzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
431 <  TH1F *LSBOSSFN = allsamples.Draw("LSBOSSFN","-"+mcjzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
432 <  TH1F *LSBOSOFN = allsamples.Draw("LSBOSOFN","-"+mcjzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
429 >  TH1F *LSBOSSFP;
430 >  TH1F *LSBOSOFP;
431 >  TH1F *LSBOSSFN;
432 >  TH1F *LSBOSOFN;
433 >  
434 >  flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
435 >  if(PlottingSetup::RestrictToMassPeak) {
436 >      SBOSSFP = allsamples.Draw("SBOSSFP",datajzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
437 >      SBOSOFP = allsamples.Draw("SBOSOFP",datajzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
438 >      SBOSSFN = allsamples.Draw("SBOSSFN","-"+datajzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
439 >      SBOSOFN = allsamples.Draw("SBOSOFN","-"+datajzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
440 >      
441 >      LSBOSSFP = allsamples.Draw("LSBOSSFP",mcjzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
442 >      LSBOSOFP = allsamples.Draw("LSBOSOFP",mcjzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
443 >      LSBOSSFN = allsamples.Draw("LSBOSSFN","-"+mcjzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
444 >      LSBOSOFN = allsamples.Draw("LSBOSOFN","-"+mcjzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
445 >  }
446    
447    string obsname="data_obs";
448    string predname="background";
# Line 402 | Line 453 | void limit_shapes_for_systematic_effect(
453      signalname="signal_"+identifier;
454    }
455    
456 <  TH1F *obs = (TH1F*)ZOSSFP->Clone();
456 >  TH1F *obs = (TH1F*)ZOSSFP->Clone("observation");
457    obs->SetName(obsname.c_str());
458    obs->Write();
459 <  TH1F *pred = (TH1F*)ZOSSFN->Clone();
460 <  pred->Add(ZOSOFP,1.0/3);
461 <  pred->Add(ZOSOFN,-1.0/3);
462 <  pred->Add(SBOSSFP,1.0/3);
463 <  pred->Add(SBOSSFN,-1.0/3);
464 <  pred->Add(SBOSOFP,1.0/3);
465 <  pred->Add(SBOSOFN,-1.0/3);
459 >  TH1F *pred = (TH1F*)ZOSSFN->Clone("prediction");
460 >    flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
461 >  if(PlottingSetup::RestrictToMassPeak) {
462 >    pred->Add(ZOSOFP,1.0/3);
463 >    pred->Add(ZOSOFN,-1.0/3);
464 >    pred->Add(SBOSSFP,1.0/3);
465 >    pred->Add(SBOSSFN,-1.0/3);
466 >    pred->Add(SBOSOFP,1.0/3);
467 >    pred->Add(SBOSOFN,-1.0/3);
468 >  } else {
469 >    pred->Add(ZOSOFP,1.0);
470 >    pred->Add(ZOSOFN,-1.0);
471 >  }
472 >    
473    pred->SetName(predname.c_str());
474    pred->Write();
475    
# Line 422 | Line 480 | void limit_shapes_for_systematic_effect(
480    TH1F *Lpred = new TH1F("Lpred","Lpred",binning.size()-1,&binning[0]);
481    Lobs->Add(LZOSSFP);
482    Lpred->Add(LZOSSFN);
483 <  Lpred->Add(LZOSOFP,1.0/3);
484 <  Lpred->Add(LZOSOFN,-1.0/3);
485 <  Lpred->Add(LSBOSSFP,1.0/3);
486 <  Lpred->Add(LSBOSSFN,-1.0/3);
487 <  Lpred->Add(LSBOSOFP,1.0/3);
488 <  Lpred->Add(LSBOSOFN,-1.0/3);
483 >    flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
484 >  if(PlottingSetup::RestrictToMassPeak) {
485 >    Lpred->Add(LZOSOFP,1.0/3);
486 >    Lpred->Add(LZOSOFN,-1.0/3);
487 >    Lpred->Add(LSBOSSFP,1.0/3);
488 >    Lpred->Add(LSBOSSFN,-1.0/3);
489 >    Lpred->Add(LSBOSOFP,1.0/3);
490 >    Lpred->Add(LSBOSOFN,-1.0/3);
491 >  } else {
492 >    Lpred->Add(LZOSOFP,1.0);
493 >    Lpred->Add(LZOSOFN,-1.0);
494 >  }
495 >
496    TH1F *signal = (TH1F*)Lobs->Clone();
497    signal->Add(Lpred,-1);
498    signal->SetName(signalname.c_str());
# Line 441 | Line 506 | void limit_shapes_for_systematic_effect(
506    delete ZOSSFN;
507    delete ZOSOFN;
508    
509 <  delete SBOSSFP;
510 <  delete SBOSOFP;
511 <  delete SBOSSFN;
512 <  delete SBOSOFN;
509 >  if(PlottingSetup::RestrictToMassPeak) {
510 >    delete SBOSSFP;
511 >    delete SBOSOFP;
512 >    delete SBOSSFN;
513 >    delete SBOSOFN;
514 >  }
515    
516    delete LZOSSFP;
517    delete LZOSOFP;
518    delete LZOSSFN;
519    delete LZOSOFN;
520    
521 <  delete LSBOSSFP;
522 <  delete LSBOSOFP;
523 <  delete LSBOSSFN;
524 <  delete LSBOSOFN;
521 >  if(PlottingSetup::RestrictToMassPeak) {
522 >    delete LSBOSSFP;
523 >    delete LSBOSOFP;
524 >    delete LSBOSSFN;
525 >    delete LSBOSOFN;
526 >  }
527  
528   }
529  

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines