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.6 by buchmann, Fri Jul 22 10:15:34 2011 UTC vs.
Revision 1.28 by buchmann, Mon Nov 7 15:24:59 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, bool doobserved=false) {
187 <  float sigma95=0.0,sigma95A=0.0;
188 <  int nuisancemodel=1;
189 <  dout << "Now calling : CL95(" << luminosity << "," <<  lumiuncert*luminosity << "," << mceff << "," << mcefferr << "," << Npred[ibin] << "," << Nprederr[ibin] << "," << Nobs[ibin] << "," << false << "," << nuisancemodel<< ") " << endl;
190 <  sigma95 = CL95(luminosity, lumiuncert*luminosity, mceff, mcefferr, Npred[ibin], Nprederr[ibin], Nobs[ibin], false, nuisancemodel);
191 <  if(doobserved) {
192 <    dout << "Now calling : CLA(" << luminosity << "," <<  lumiuncert*luminosity << "," << mceff << "," << mcefferr << "," << Npred[ibin] << "," << Nprederr[ibin] << "," << nuisancemodel<< ") " << endl;
193 <    sigma95A = CLA(luminosity, lumiuncert*luminosity, mceff, mcefferr, Npred[ibin], Nprederr[ibin], nuisancemodel);
194 <  }
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
190 > " Double_t             limit = roostats_cl95(ilum, slum, eff, seff, bck, sbck, n, gauss = false, nuisanceModel, method, plotFileName, seed); \n"
191 > " LimitResult expected_limit = roostats_clm(ilum, slum, eff, seff, bck, sbck, ntoys, nuisanceModel, method, seed); \n"
192 > " Double_t     average_limit = roostats_cla(ilum, slum, eff, seff, bck, sbck, nuisanceModel, method, seed); \n"
193 > "                                                                     \n"
194 > "
195 > " Double_t obs_limit = limit.GetObservedLimit();                      \n"
196 > " Double_t exp_limit = limit.GetExpectedLimit();                      \n"
197 > " Double_t exp_up    = limit.GetOneSigmaHighRange();                  \n"
198 > " Double_t exp_down  = limit.GetOneSigmaLowRange();                   \n"
199 > " Double_t exp_2up   = limit.GetTwoSigmaHighRange();                  \n"
200 > " Double_t exp_2down = limit.GetTwoSigmaLowRange();                   \n"
201 > */
202 >  if(mceff<=0) {
203 >    write_warning(__FUNCTION__,"Cannot compute upper limit in this configuration as the efficiency is negative:");
204 >    dout << "mc efficiency=" << mceff << " +/- " << mcefferr;
205 >    vector<float> sigmas;
206 >    sigmas.push_back(-1);
207 >    sigmas.push_back(-1);
208 >    return sigmas;
209 >  } else {
210 >    int nlimittoysused=1;
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 >  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 >    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 (TimedLimitCapsule) now : Attempt " << attempts << endl;
242 >        retval=gSystem->Exec(command.str().c_str());
243 >    }
244 >    
245 >    LimitDroplet limres;
246 >    limres.readDroplet(repname.str());
247 >    dout << limres << endl;
248 >    remove(repname.str().c_str());
249 >    sigma95=limres.observed;
250 >
251 >    
252 >    ///------------------------------------------ < /NEW > ----------------------------------------------------------
253    vector<float> sigmas;
254    sigmas.push_back(sigma95);
255 <  sigmas.push_back(sigma95A);
255 >  if(doexpected) {
256 >    sigmas.push_back(limres.expected);
257 >    sigmas.push_back(limres.upper68);
258 >    sigmas.push_back(limres.lower68);
259 >    sigmas.push_back(limres.upper95);
260 >    sigmas.push_back(limres.lower95);
261 >  }
262 >  
263    return sigmas;
264 +  
265 +
266 +  }//end of mc efficiency is ok
267   }
268  
269 < void compute_upper_limits_from_counting_experiment(vector<vector<float> > uncertainties,vector<float> jzbcuts, string mcjzb, bool doobserved) {
269 > void compute_upper_limits_from_counting_experiment(vector<vector<float> > uncertainties,vector<float> jzbcuts, string mcjzb, bool doexpected, int flipped) {
270    dout << "Doing counting experiment ... " << endl;
271    vector<vector<string> > limits;
272    vector<vector<float> > vlimits;
# Line 201 | Line 284 | void compute_upper_limits_from_counting_
284        float staterr=uncertainties[isample*jzbcuts.size()+ibin][2];
285        float systerr=uncertainties[isample*jzbcuts.size()+ibin][3];
286        float toterr =uncertainties[isample*jzbcuts.size()+ibin][4];
287 <      float observed,null,result;
205 <      fill_result_histos(observed, null,null,null,null,null,null,null,mcjzb,JZBcutat,(int)5,result,(signalsamples.FindSample(signalsamples.collection[isample].filename)),signalsamples);
206 <      observed-=result;//this is the actual excess we see!
207 <      float expected=observed/luminosity;
287 >      float observed,observederr,null,result;
288        
289 + //      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);
290 + //      observed-=result;//this is the actual excess we see!
291 + //      float expected=observed/luminosity;
292 +      string plotfilename=(string)(TString(signalsamples.collection[isample].samplename)+TString("___JZB_geq_")+TString(any2string(JZBcutat))+TString(".png"));
293        dout << "Sample: " << signalsamples.collection[isample].samplename << ", JZB>"<<JZBcutat<< " : " << mceff << " +/- " << staterr << " (stat) +/- " << systerr << " (syst) --> toterr = " << toterr << endl;
294 <      vector<float> sigmas = compute_one_upper_limit(mceff,toterr,ibin,mcjzb,doobserved);
294 >      vector<float> sigmas = compute_one_upper_limit(mceff,toterr,ibin,mcjzb,plotfilename,doexpected,flipped);
295        
296 <      if(doobserved) {
297 <        rows.push_back(any2string(sigmas[0])+";"+any2string(sigmas[1])+";"+"("+any2string(expected)+")");
296 >      if(doexpected) {
297 > //      rows.push_back(any2string(sigmas[0])+";"+any2string(sigmas[1])+";"+"("+any2string(expected)+")");
298 >        rows.push_back(any2string(sigmas[0])+";"+any2string(sigmas[1])+";"+"("+any2string(signalsamples.collection[isample].xs)+")");
299          vrows.push_back(sigmas[0]);
300          vrows.push_back(sigmas[1]);
301 <        vrows.push_back(expected);
301 > //      vrows.push_back(expected);
302 >        vrows.push_back(signalsamples.collection[isample].xs);
303        }
304        else {
305 <        rows.push_back(any2string(sigmas[0])+"("+any2string(expected)+")");
305 > //      rows.push_back(any2string(sigmas[0])+"("+any2string(expected)+")");
306 >        rows.push_back(any2string(sigmas[0]));
307          vrows.push_back(sigmas[0]);
308 <        vrows.push_back(expected);
308 >        vrows.push_back(signalsamples.collection[isample].xs);
309 > //      vrows.push_back(expected);
310        }
311      }//end of bin loop
312      limits.push_back(rows);
313      vlimits.push_back(vrows);
314    }//end of sample loop
315 <  dout << endl << endl << "PAS table 3: " << endl << endl;
315 >  dout << endl << endl << endl << "_________________________________________________________________________________________________" << endl << endl;
316 >  dout << endl << endl << "PAS table 3:   (notation: limit [95%CL])" << endl << endl;
317    dout << "\t";
318    for (int irow=0;irow<jzbcuts.size();irow++) {
319      dout << jzbcuts[irow] << "\t";
# Line 232 | Line 321 | void compute_upper_limits_from_counting_
321    dout << endl;
322    for(int irow=0;irow<limits.size();irow++) {
323      for(int ientry=0;ientry<limits[irow].size();ientry++) {
324 <      dout << limits[irow][ientry] << "\t";
324 >      if (limits[irow][ientry]>0) dout << limits[irow][ientry] << "\t";
325 >      else dout << " (N/A) \t";
326      }
327      dout << endl;
328    }
329    
330 <  if(!doobserved) {
331 <    dout << endl << endl << "LIMITS: " << endl;
332 <    dout << "\t";
330 >  if(!doexpected) {
331 >    dout << endl << endl << "LIMITS: (Tex)" << endl;
332 >    tout << "\\begin{table}[hbtp]" << endl;
333 >    tout << "\\renewcommand{\\arraystretch}{1.3}" << endl;
334 >    tout << "\\begin{center}" << endl;
335 >    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;
336 >    tout << "" << endl;
337 >    tout << "\\begin{tabular}{ | l | ";
338 >    for (int irow=0;irow<jzbcuts.size();irow++) tout << " l |";
339 >    tout << "} " << endl << " \\hline " << endl << "& \t ";
340      for (int irow=0;irow<jzbcuts.size();irow++) {
341 <      dout << jzbcuts[irow] << "\t";
341 >      tout << "JZB $>$ " << jzbcuts[irow] << " GeV & \t ";
342      }
343 <    dout << endl;
343 >    tout << " \\\\ \\hline " << endl;
344      for(int irow=0;irow<limits.size();irow++) {
345 <      dout << limits[irow][0] << "\t";
345 >      tout << limits[irow][0] << " \t";
346        for(int ientry=0;ientry<jzbcuts.size();ientry++) {
347 <        dout << Round(vlimits[irow][2*ientry] / vlimits[irow][2*ientry+1],3)<< "\t";
347 >        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";
348 >        else tout << " & ( N / A ) \t";
349 > //      dout << Round(vlimits[irow][2*ientry],3) << " / " << Round(vlimits[irow][2*ientry+1],3)<< "\t";
350        }
351 <      dout << endl;
351 >      tout << " \\\\ \\hline " << endl;
352      }
353 +      tout << "\\end{tabular}" << endl;
354 +      tout << "      \\end{tabular}"<< endl;
355 +      tout << "\\end{center}"<< endl;
356 +      tout << "\\end{table} "<< endl;
357 +
358    }//do observed
359    
360    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;
361 <  dout << "Scenario \t Efficiency [%] \t Upper limits [pb] \t Prediction [pb]" << endl;
361 >  dout << "Scenario \t Efficiency [%] \t Upper limits [pb] \t \\sigma [pb]" << endl;
362    for(int icut=0;icut<jzbcuts.size();icut++) {
363 <    dout << "Region with JZB>" << jzbcuts[icut] << endl;
363 >    dout << "Region with JZB>" << jzbcuts[icut] << (ConsiderSignalContaminationForLimits?"  (accounting for signal contamination)":"  (not accounting for signal contamination)") << endl;
364      for(int isample=0;isample<signalsamples.collection.size();isample++) {
365 <      dout << limits[icut][0] << "\t" << Round(100*uncertainties[isample*jzbcuts.size()+icut][1],1) << "+/-" << Round(100*uncertainties[isample*jzbcuts.size()+icut][2],1) << " (stat) +/- " << Round(100*uncertainties[isample*jzbcuts.size()+icut][3],1) << " (syst) \t" << Round((vlimits[isample][2*icut]),3) << "\t" << Round(vlimits[isample][2*icut+1],3) << endl;
365 >      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;
366      }
367      dout << endl;
368    }
265  
266  write_warning("compute_upper_limits_from_counting_experiment","Still need to update the script");
369   }
370  
371  
372  
373 < //********************************************************************** new : Limits using SHAPES ***********************************
373 > /********************************************************************** new : Limits using SHAPES ***********************************
374 >
375 >
376 >   SSSSSSSSSSSSSSS hhhhhhh                                                                                      
377 > SS:::::::::::::::Sh:::::h                                                                                      
378 > S:::::SSSSSS::::::Sh:::::h                                                                                      
379 > S:::::S     SSSSSSSh:::::h                                                                                      
380 > S:::::S             h::::h hhhhh         aaaaaaaaaaaaa  ppppp   ppppppppp       eeeeeeeeeeee        ssssssssss  
381 > S:::::S             h::::hh:::::hhh      a::::::::::::a p::::ppp:::::::::p    ee::::::::::::ee    ss::::::::::s  
382 > S::::SSSS          h::::::::::::::hh    aaaaaaaaa:::::ap:::::::::::::::::p  e::::::eeeee:::::eess:::::::::::::s
383 >  SS::::::SSSSS     h:::::::hhh::::::h            a::::app::::::ppppp::::::pe::::::e     e:::::es::::::ssss:::::s
384 >    SSS::::::::SS   h::::::h   h::::::h    aaaaaaa:::::a p:::::p     p:::::pe:::::::eeeee::::::e s:::::s  ssssss
385 >       SSSSSS::::S  h:::::h     h:::::h  aa::::::::::::a p:::::p     p:::::pe:::::::::::::::::e    s::::::s      
386 >            S:::::S h:::::h     h:::::h a::::aaaa::::::a p:::::p     p:::::pe::::::eeeeeeeeeee        s::::::s  
387 >            S:::::S h:::::h     h:::::ha::::a    a:::::a p:::::p    p::::::pe:::::::e           ssssss   s:::::s
388 > SSSSSSS     S:::::S h:::::h     h:::::ha::::a    a:::::a p:::::ppppp:::::::pe::::::::e          s:::::ssss::::::s
389 > S::::::SSSSSS:::::S h:::::h     h:::::ha:::::aaaa::::::a p::::::::::::::::p  e::::::::eeeeeeee  s::::::::::::::s
390 > S:::::::::::::::SS  h:::::h     h:::::h a::::::::::aa:::ap::::::::::::::pp    ee:::::::::::::e   s:::::::::::ss  
391 > SSSSSSSSSSSSSSS    hhhhhhh     hhhhhhh  aaaaaaaaaa  aaaap::::::pppppppp        eeeeeeeeeeeeee    sssssssssss    
392 >                                                         p:::::p                                                
393 >                                                         p:::::p                                                
394 >                                                        p:::::::p                                                
395 >                                                        p:::::::p                                                
396 >                                                        p:::::::p                                                
397 >                                                        ppppppppp                                                
398 >                                                                                                                
399 >
400 > *********************************************************************** new : Limits using SHAPES ***********************************/
401 >
402  
403   void limit_shapes_for_systematic_effect(TFile *limfile, string identifier, string mcjzb, string datajzb, int JES,vector<float> binning, TCanvas *limcan) {
404    dout << "Creatig shape templates ... ";
# Line 288 | Line 418 | void limit_shapes_for_systematic_effect(
418    TH1F *ZOSSFN = allsamples.Draw("ZOSSFN","-"+datajzb,binning, "JZB4limits", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,dataormc,luminosity);
419    TH1F *ZOSOFN = allsamples.Draw("ZOSOFN","-"+datajzb,binning, "JZB4limits", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,dataormc,luminosity);
420    
421 <  TH1F *SBOSSFP = allsamples.Draw("SBOSSFP",datajzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
422 <  TH1F *SBOSOFP = allsamples.Draw("SBOSOFP",datajzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
423 <  TH1F *SBOSSFN = allsamples.Draw("SBOSSFN","-"+datajzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
424 <  TH1F *SBOSOFN = allsamples.Draw("SBOSOFN","-"+datajzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
421 >  TH1F *SBOSSFP;
422 >  TH1F *SBOSOFP;
423 >  TH1F *SBOSSFN;
424 >  TH1F *SBOSOFN;
425    
426    TH1F *LZOSSFP = allsamples.Draw("LZOSSFP",mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("LM4"));
427    TH1F *LZOSOFP = allsamples.Draw("LZOSOFP",mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("LM4"));
428    TH1F *LZOSSFN = allsamples.Draw("LZOSSFN","-"+mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("LM4"));
429    TH1F *LZOSOFN = allsamples.Draw("LZOSOFN","-"+mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("LM4"));
430    
431 <  TH1F *LSBOSSFP = allsamples.Draw("LSBOSSFP",mcjzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
432 <  TH1F *LSBOSOFP = allsamples.Draw("LSBOSOFP",mcjzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
433 <  TH1F *LSBOSSFN = allsamples.Draw("LSBOSSFN","-"+mcjzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
434 <  TH1F *LSBOSOFN = allsamples.Draw("LSBOSOFN","-"+mcjzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
431 >  TH1F *LSBOSSFP;
432 >  TH1F *LSBOSOFP;
433 >  TH1F *LSBOSSFN;
434 >  TH1F *LSBOSOFN;
435 >  
436 >  flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
437 >  if(PlottingSetup::RestrictToMassPeak) {
438 >      SBOSSFP = allsamples.Draw("SBOSSFP",datajzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
439 >      SBOSOFP = allsamples.Draw("SBOSOFP",datajzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
440 >      SBOSSFN = allsamples.Draw("SBOSSFN","-"+datajzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
441 >      SBOSOFN = allsamples.Draw("SBOSOFN","-"+datajzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
442 >      
443 >      LSBOSSFP = allsamples.Draw("LSBOSSFP",mcjzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
444 >      LSBOSOFP = allsamples.Draw("LSBOSOFP",mcjzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
445 >      LSBOSSFN = allsamples.Draw("LSBOSSFN","-"+mcjzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
446 >      LSBOSOFN = allsamples.Draw("LSBOSOFN","-"+mcjzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
447 >  }
448    
449    string obsname="data_obs";
450    string predname="background";
# Line 312 | Line 455 | void limit_shapes_for_systematic_effect(
455      signalname="signal_"+identifier;
456    }
457    
458 <  TH1F *obs = (TH1F*)ZOSSFP->Clone();
458 >  TH1F *obs = (TH1F*)ZOSSFP->Clone("observation");
459    obs->SetName(obsname.c_str());
460    obs->Write();
461 <  TH1F *pred = (TH1F*)ZOSSFN->Clone();
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);
461 >  TH1F *pred = (TH1F*)ZOSSFN->Clone("prediction");
462 >    flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
463 >  if(PlottingSetup::RestrictToMassPeak) {
464 >    pred->Add(ZOSOFP,1.0/3);
465 >    pred->Add(ZOSOFN,-1.0/3);
466 >    pred->Add(SBOSSFP,1.0/3);
467 >    pred->Add(SBOSSFN,-1.0/3);
468 >    pred->Add(SBOSOFP,1.0/3);
469 >    pred->Add(SBOSOFN,-1.0/3);
470 >  } else {
471 >    pred->Add(ZOSOFP,1.0);
472 >    pred->Add(ZOSOFN,-1.0);
473 >  }
474 >    
475    pred->SetName(predname.c_str());
476    pred->Write();
477    
# Line 332 | Line 482 | void limit_shapes_for_systematic_effect(
482    TH1F *Lpred = new TH1F("Lpred","Lpred",binning.size()-1,&binning[0]);
483    Lobs->Add(LZOSSFP);
484    Lpred->Add(LZOSSFN);
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);
485 >    flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
486 >  if(PlottingSetup::RestrictToMassPeak) {
487 >    Lpred->Add(LZOSOFP,1.0/3);
488 >    Lpred->Add(LZOSOFN,-1.0/3);
489 >    Lpred->Add(LSBOSSFP,1.0/3);
490 >    Lpred->Add(LSBOSSFN,-1.0/3);
491 >    Lpred->Add(LSBOSOFP,1.0/3);
492 >    Lpred->Add(LSBOSOFN,-1.0/3);
493 >  } else {
494 >    Lpred->Add(LZOSOFP,1.0);
495 >    Lpred->Add(LZOSOFN,-1.0);
496 >  }
497 >
498    TH1F *signal = (TH1F*)Lobs->Clone();
499    signal->Add(Lpred,-1);
500    signal->SetName(signalname.c_str());
# Line 351 | Line 508 | void limit_shapes_for_systematic_effect(
508    delete ZOSSFN;
509    delete ZOSOFN;
510    
511 <  delete SBOSSFP;
512 <  delete SBOSOFP;
513 <  delete SBOSSFN;
514 <  delete SBOSOFN;
511 >  if(PlottingSetup::RestrictToMassPeak) {
512 >    delete SBOSSFP;
513 >    delete SBOSOFP;
514 >    delete SBOSSFN;
515 >    delete SBOSOFN;
516 >  }
517    
518    delete LZOSSFP;
519    delete LZOSOFP;
520    delete LZOSSFN;
521    delete LZOSOFN;
522    
523 <  delete LSBOSSFP;
524 <  delete LSBOSOFP;
525 <  delete LSBOSSFN;
526 <  delete LSBOSOFN;
523 >  if(PlottingSetup::RestrictToMassPeak) {
524 >    delete LSBOSSFP;
525 >    delete LSBOSOFP;
526 >    delete LSBOSSFN;
527 >    delete LSBOSOFN;
528 >  }
529  
530   }
531  
# Line 414 | Line 575 | void prepare_limits(string mcjzb, string
575    limfile->Close();
576    write_info("prepare_limits","limitfile.root and datacard.txt have been generated. You can now use them to calculate limits!");
577    
578 < }
578 > }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines