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

Comparing UserCode/cbrown/Development/Plotting/Modules/LimitCalculation.C (file contents):
Revision 1.2 by buchmann, Thu Mar 15 20:35:39 2012 UTC vs.
Revision 1.7 by buchmann, Wed Aug 15 13:52:12 2012 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 – ****/
1   #include <iostream>
2   #include <vector>
3   #include <sys/stat.h>
# Line 55 | Line 42 | void rediscover_the_top(string mcjzb, st
42    TH1F *puresignal     = allsamples.Draw("puresignal",    datajzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data,  luminosity);
43   //  TH1F *puresignal     = allsamples.Draw("puresignal",        mcjzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,mc,  luminosity,allsamples.FindSample("TTJets"));
44    TH1F *observed       = allsamples.Draw("observed",          datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
45 <  /*
45 >  
46    ofstream myfile;
47    TH1F *ratio = (TH1F*)observed->Clone();
48    ratio->Divide(dataprediction);
# Line 148 | Line 135 | ratio_binning.push_back(80);
135    TH1F *datapredictiontbo = allsamples.Draw("datapredictiontbo",    "-"+datajzb, ratio_binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data,  luminosity);
136    datapredictiontb->Add(datapredictiontbo,-1);
137    TH1F *analytical_background_predictionb = allsamples.Draw("analytical_background_predictionb",datajzb, ratio_binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets&&"mll<2",data,  luminosity);
138 <  for(int i=0;i<=ratio_binning.size();i++) {
138 >  for(int i=0;i<=(int)ratio_binning.size();i++) {
139      analytical_background_predictionb->SetBinContent(i+1,functions[1]->Eval(analytical_background_predictionb->GetBinCenter(i)));
140      analytical_background_predictionb->SetBinError(i+1,TMath::Sqrt(functions[1]->Eval(analytical_background_predictionb->GetBinCenter(i))));
141    }
# Line 156 | Line 143 | ratio_binning.push_back(80);
143    TH1F *ratio = (TH1F*) observedtb->Clone();
144    ratio->Divide(datapredictiontb);
145    
146 <  for (int i=0;i<=ratio_binning.size();i++) {
146 >  for (int i=0;i<=(int)ratio_binning.size();i++) {
147      dout << observedtb->GetBinLowEdge(i+1) << ";"<<observedtb->GetBinContent(i+1) << ";" << datapredictiontb->GetBinContent(i+1) << " --> " << ratio->GetBinContent(i+1) << "+/-" << ratio->GetBinError(i+1) << endl;
148    }
149    
# Line 184 | Line 171 | ratio_binning.push_back(80);
171   }
172  
173   vector<float> compute_one_upper_limit(float mceff,float mcefferr, int ibin, string mcjzb, string plotfilename, bool doexpected, int flipped, bool doasymptotic=false) {
174 <  float sigma95=-9.9,sigma95A=-9.9;
174 >  float sigma95=-9.9;
175   /*
176   USAGE OF ROOSTATS_CL95
177   " Double_t             limit = roostats_cl95(ilum, slum, eff, seff, bck, sbck, n, gauss = false, nuisanceModel, method, plotFileName, seed); \n"
# Line 207 | Line 194 | USAGE OF ROOSTATS_CL95
194      sigmas.push_back(-1);
195      return sigmas;
196    } else {
210 –    int nlimittoysused=1;
211 –    
197      ///------------------------------------------ < NEW > ----------------------------------------------------------
198      
199      int secondssince1970=time(NULL);
# Line 228 | Line 213 | USAGE OF ROOSTATS_CL95
213  
214    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;
215    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;
216 +  
217 +  stringstream command;
218 +  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 << " " << doasymptotic;
219 +  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<< " " << doasymptotic;
220 +  
221 +  dout << command.str() << endl;
222      
223 <    stringstream command;
224 <    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 << " " << doasymptotic;
225 <    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<< " " << doasymptotic;
226 <    dout << command.str() << endl;
227 <    
237 <    int retval = 256;
238 <    int attempts=0;
239 <    while(!(retval==0||attempts>=3)) {//try up to 3 times
223 >  if(doasymptotic) write_warning(__FUNCTION__, "DOING ASYMPTOTIC LIMIT!");
224 >  
225 >  int retval = 256;
226 >  int attempts=0;
227 >  while(!(retval==0||attempts>=3)) {//try up to 3 times
228          attempts++;
229          dout << "Starting limit calculation (TimedLimitCapsule) now : Attempt " << attempts << endl;
230          retval=gSystem->Exec(command.str().c_str());
# Line 255 | Line 243 | USAGE OF ROOSTATS_CL95
243      remove(repname.str().c_str());
244      sigma95=limres.observed;
245  
258 –    
246      ///------------------------------------------ < /NEW > ----------------------------------------------------------
247    vector<float> sigmas;
248    sigmas.push_back(sigma95);
# Line 279 | Line 266 | vector<float> compute_upper_limits_from_
266    vector<vector<float> > vlimits;
267    
268    
269 <  for(int isample=0;isample<signalsamples.collection.size();isample++) {
269 >  for(int isample=0;isample<(int)signalsamples.collection.size();isample++) {
270      vector<string> rows;
271      vector<float> vrows;
272      dout << "Considering sample " << signalsamples.collection[isample].samplename << endl;
273      rows.push_back(signalsamples.collection[isample].samplename);
274 <    for(int ibin=0;ibin<jzbcuts.size();ibin++) {
274 >    for(int ibin=0;ibin<(int)jzbcuts.size();ibin++) {
275        dout << "_________________________________________________________________________________" << endl;
276        float JZBcutat=uncertainties[isample*jzbcuts.size()+ibin][0];
277        float mceff=uncertainties[isample*jzbcuts.size()+ibin][1];
278        float staterr=uncertainties[isample*jzbcuts.size()+ibin][2];
279        float systerr=uncertainties[isample*jzbcuts.size()+ibin][3];
280        float toterr =uncertainties[isample*jzbcuts.size()+ibin][4];
281 <      float observed,observederr,null,result;
281 > //      float observed,observederr,null,result;
282        
283        string plotfilename=(string)(TString(signalsamples.collection[isample].samplename)+TString("___JZB_geq_")+TString(any2string(JZBcutat))+TString(".png"));
284        dout << "Sample: " << signalsamples.collection[isample].samplename << ", JZB>"<<JZBcutat<< " : " << mceff << " +/- " << staterr << " (stat) +/- " << systerr << " (syst) --> toterr = " << toterr << endl;
# Line 322 | Line 309 | vector<float> compute_upper_limits_from_
309  
310    dout << endl << endl << "_______________________________________________________________________________________" << endl;
311    dout << "Going to store upper limit on event yield in result library: " << endl;
312 <  for(int ibin=0;ibin<jzbcuts.size();ibin++) {
312 >  for(int ibin=0;ibin<(int)jzbcuts.size();ibin++) {
313        int resultindex=PlottingSetup::allresults.Find(jzbcuts[ibin]);
314        vector<float> Normsigmas = compute_one_upper_limit(1.0,0.0, resultindex, mcjzb, "UPPERLIMIT", false, 0);
315        (allresults.predictions[resultindex]).UpperLimit=Normsigmas[0]*PlottingSetup::luminosity;
# Line 332 | Line 319 | vector<float> compute_upper_limits_from_
319    dout << endl << endl << endl << "_________________________________________________________________________________________________" << endl << endl;
320    dout << endl << endl << "PAS table 3:   (notation: limit [95%CL])" << endl << endl;
321    dout << "\t";
322 <  for (int irow=0;irow<jzbcuts.size();irow++) {
322 >  for (int irow=0;irow<(int)jzbcuts.size();irow++) {
323      dout << jzbcuts[irow] << "\t";
324    }
325    dout << endl;
326 <  for(int irow=0;irow<limits.size();irow++) {
327 <    for(int ientry=0;ientry<limits[irow].size();ientry++) {
326 >  for(int irow=0;irow<(int)limits.size();irow++) {
327 >    for(int ientry=0;ientry<(int)limits[irow].size();ientry++) {
328        if (limits[irow][ientry]>0) dout << limits[irow][ientry] << "\t";
329        else dout << " (N/A) \t";
330      }
# Line 352 | Line 339 | vector<float> compute_upper_limits_from_
339      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;
340      tout << "" << endl;
341      tout << "\\begin{tabular}{ | l | ";
342 <    for (int irow=0;irow<jzbcuts.size();irow++) tout << " l |";
342 >    for (int irow=0;irow<(int)jzbcuts.size();irow++) tout << " l |";
343      tout << "} " << endl << " \\hline " << endl << "& \t ";
344 <    for (int irow=0;irow<jzbcuts.size();irow++) {
344 >    for (int irow=0;irow<(int)jzbcuts.size();irow++) {
345        tout << "JZB $>$ " << jzbcuts[irow] << " GeV & \t ";
346      }
347      tout << " \\\\ \\hline " << endl;
348 <    for(int irow=0;irow<limits.size();irow++) {
348 >    for(int irow=0;irow<(int)limits.size();irow++) {
349        tout << limits[irow][0] << " \t";
350 <      for(int ientry=0;ientry<jzbcuts.size();ientry++) {
350 >      for(int ientry=0;ientry<(int)jzbcuts.size();ientry++) {
351          if(vlimits[irow][2*ientry]>0) tout << " & " << Round(vlimits[irow][2*ientry],2) << " \t (" << Round(vlimits[irow][2*ientry] / vlimits[irow][2*ientry+1],3)<< "x \\sigma ) \t";
352          else tout << " & ( N / A ) \t";
353   //      dout << Round(vlimits[irow][2*ientry],3) << " / " << Round(vlimits[irow][2*ientry+1],3)<< "\t";
# Line 376 | Line 363 | vector<float> compute_upper_limits_from_
363    
364    dout << endl << endl << "Final selection efficiencies with total statistical and systematic errors, and corresponding observed and expected upper limits (UL) on ($\\sigma\\times$  BR $\\times$ acceptance) for the LM4 and LM8 scenarios, in the different regions. The last column contains the predicted ($\\sigma \\times $BR$\\times$ acceptance) at NLO obtained from Monte Carlo simulation." << endl;
365    dout << "Scenario \t Efficiency [%] \t Upper limits [pb] \t \\sigma [pb]" << endl;
366 <  for(int icut=0;icut<jzbcuts.size();icut++) {
366 >  for(int icut=0;icut<(int)jzbcuts.size();icut++) {
367      dout << "Region with JZB>" << jzbcuts[icut] << (ConsiderSignalContaminationForLimits?"  (accounting for signal contamination)":"  (not accounting for signal contamination)") << endl;
368 <    for(int isample=0;isample<signalsamples.collection.size();isample++) {
368 >    for(int isample=0;isample<(int)signalsamples.collection.size();isample++) {
369        dout << limits[isample][0] << "\t" << Round(100*uncertainties[isample*jzbcuts.size()+icut][1],3) << "+/-" << Round(100*uncertainties[isample*jzbcuts.size()+icut][2],3) << " (stat) +/- " << Round(100*uncertainties[isample*jzbcuts.size()+icut][3],3) << " (syst) \t" << Round((vlimits[isample][2*icut]),3) << "\t" << Round(vlimits[isample][2*icut+1],3) << endl;
370      }
371      dout << endl;
# Line 388 | Line 375 | vector<float> compute_upper_limits_from_
375    //---------------------------------------------
376    
377    vector<float> lowestULs;
378 <  for(int isample=0;isample<signalsamples.collection.size();isample++) {
378 >  for(int isample=0;isample<(int)signalsamples.collection.size();isample++) {
379      float lowestUL=-1;
380 <    for(int icut=0;icut<jzbcuts.size();icut++) {
380 >    for(int icut=0;icut<(int)jzbcuts.size();icut++) {
381        float currUL=Round((vlimits[isample][2*icut]),3);
382        if(currUL>0) {
383          if(lowestUL<0) lowestUL=currUL;
# Line 470 | Line 457 | void limit_shapes_for_systematic_effect(
457    TH1F *LSBOSOFN;
458    
459    flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
460 <  if(PlottingSetup::RestrictToMassPeak) {
460 >  if(PlottingSetup::RestrictToMassPeak&&PlottingSetup::UseSidebandsForcJZB) {
461        SBOSSFP = allsamples.Draw("SBOSSFP",datajzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
462        SBOSOFP = allsamples.Draw("SBOSOFP",datajzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
463        SBOSSFN = allsamples.Draw("SBOSSFN","-"+datajzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
# Line 496 | Line 483 | void limit_shapes_for_systematic_effect(
483    obs->Write();
484    TH1F *pred = (TH1F*)ZOSSFN->Clone("prediction");
485      flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
486 <  if(PlottingSetup::RestrictToMassPeak) {
486 >  if(PlottingSetup::RestrictToMassPeak&&PlottingSetup::UseSidebandsForcJZB) {
487      pred->Add(ZOSOFP,1.0/3);
488      pred->Add(ZOSOFN,-1.0/3);
489      pred->Add(SBOSSFP,1.0/3);
# Line 519 | Line 506 | void limit_shapes_for_systematic_effect(
506    Lobs->Add(LZOSSFP);
507    Lpred->Add(LZOSSFN);
508      flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
509 <  if(PlottingSetup::RestrictToMassPeak) {
509 >  if(PlottingSetup::RestrictToMassPeak&&PlottingSetup::UseSidebandsForcJZB) {
510      Lpred->Add(LZOSOFP,1.0/3);
511      Lpred->Add(LZOSOFN,-1.0/3);
512      Lpred->Add(LSBOSSFP,1.0/3);
# Line 544 | Line 531 | void limit_shapes_for_systematic_effect(
531    delete ZOSSFN;
532    delete ZOSOFN;
533    
534 <  if(PlottingSetup::RestrictToMassPeak) {
534 >  if(PlottingSetup::RestrictToMassPeak&&PlottingSetup::UseSidebandsForcJZB) {
535      delete SBOSSFP;
536      delete SBOSOFP;
537      delete SBOSSFN;
# Line 556 | Line 543 | void limit_shapes_for_systematic_effect(
543    delete LZOSSFN;
544    delete LZOSOFN;
545    
546 <  if(PlottingSetup::RestrictToMassPeak) {
546 >  if(PlottingSetup::RestrictToMassPeak&&PlottingSetup::UseSidebandsForcJZB) {
547      delete LSBOSSFP;
548      delete LSBOSOFP;
549      delete LSBOSSFN;

Diff Legend

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