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.5 by buchmann, Wed Jul 20 14:34:31 2011 UTC vs.
Revision 1.31 by buchmann, Thu Nov 24 08:19:17 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>
17 + #include <fstream>
18  
19   #include <TCut.h>
20   #include <TROOT.h>
# Line 14 | 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 167 | 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 >    char hostname[1023];
245 >    gethostname(hostname,1023);
246 >    if((!((Contains(hostname,"t3ui")||Contains(hostname,"t3wn"))))&&retval==256) {
247 >        //running via CRAB and encountered the same problem too often: place a problem file to mark this problem!
248 >        stringstream markproblem;
249 >        markproblem << "touch " << PlottingSetup::basedirectory << "/exchange/problemswhilesettinglimits.txt";
250 >        gSystem->Exec(markproblem.str().c_str());      
251 >    }
252 >    LimitDroplet limres;
253 >    limres.readDroplet(repname.str());
254 >    dout << limres << endl;
255 >    remove(repname.str().c_str());
256 >    sigma95=limres.observed;
257 >
258 >    
259 >    ///------------------------------------------ < /NEW > ----------------------------------------------------------
260    vector<float> sigmas;
261    sigmas.push_back(sigma95);
262 <  sigmas.push_back(sigma95A);
262 >  if(doexpected) {
263 >    sigmas.push_back(limres.expected);
264 >    sigmas.push_back(limres.upper68);
265 >    sigmas.push_back(limres.lower68);
266 >    sigmas.push_back(limres.upper95);
267 >    sigmas.push_back(limres.lower95);
268 >  }
269 >  
270    return sigmas;
271 +  
272 +
273 +  }//end of mc efficiency is ok
274   }
275  
276 < void compute_upper_limits_from_counting_experiment(vector<vector<float> > uncertainties,vector<float> jzbcuts, string mcjzb, bool doobserved) {
276 > void compute_upper_limits_from_counting_experiment(vector<vector<float> > uncertainties,vector<float> jzbcuts, string mcjzb, bool doexpected, int flipped) {
277    dout << "Doing counting experiment ... " << endl;
278    vector<vector<string> > limits;
279    vector<vector<float> > vlimits;
# Line 200 | Line 291 | void compute_upper_limits_from_counting_
291        float staterr=uncertainties[isample*jzbcuts.size()+ibin][2];
292        float systerr=uncertainties[isample*jzbcuts.size()+ibin][3];
293        float toterr =uncertainties[isample*jzbcuts.size()+ibin][4];
294 <      float observed,null,result;
204 <      fill_result_histos(observed, null,null,null,null,null,null,null,mcjzb,JZBcutat,(int)5,result,(signalsamples.FindSample(signalsamples.collection[isample].filename)),signalsamples);
205 <      observed-=result;//this is the actual excess we see!
206 <      float expected=observed/luminosity;
294 >      float observed,observederr,null,result;
295        
296 +      string plotfilename=(string)(TString(signalsamples.collection[isample].samplename)+TString("___JZB_geq_")+TString(any2string(JZBcutat))+TString(".png"));
297        dout << "Sample: " << signalsamples.collection[isample].samplename << ", JZB>"<<JZBcutat<< " : " << mceff << " +/- " << staterr << " (stat) +/- " << systerr << " (syst) --> toterr = " << toterr << endl;
298 <      vector<float> sigmas = compute_one_upper_limit(mceff,toterr,ibin,mcjzb,doobserved);
299 <      
300 <      if(doobserved) {
301 <        rows.push_back(any2string(sigmas[0])+";"+any2string(sigmas[1])+";"+"("+any2string(expected)+")");
298 >      vector<float> sigmas = compute_one_upper_limit(mceff,toterr,ibin,mcjzb,plotfilename,doexpected,flipped);
299 >
300 >      tripple LibraryUpperLimits;
301 >      LibraryUpperLimits.name=signalsamples.collection[isample].samplename;
302 >      LibraryUpperLimits.first=mceff*signalsamples.collection[isample].xs * PlottingSetup::luminosity;
303 >      LibraryUpperLimits.second=staterr*signalsamples.collection[isample].xs * PlottingSetup::luminosity;
304 >      int resultindex=PlottingSetup::allresults.Find(jzbcuts[ibin]);
305 >      (allresults.predictions[resultindex]).SignalYield.push_back(LibraryUpperLimits);
306 >
307 >      if(doexpected) {
308 >        rows.push_back(any2string(sigmas[0])+";"+any2string(sigmas[1])+";"+"("+any2string(signalsamples.collection[isample].xs)+")");
309          vrows.push_back(sigmas[0]);
310          vrows.push_back(sigmas[1]);
311 <        vrows.push_back(expected);
311 >        vrows.push_back(signalsamples.collection[isample].xs);
312        }
313        else {
314 <        rows.push_back(any2string(sigmas[0])+"("+any2string(expected)+")");
314 >        rows.push_back(any2string(sigmas[0]));
315          vrows.push_back(sigmas[0]);
316 <        vrows.push_back(expected);
316 >        vrows.push_back(signalsamples.collection[isample].xs);
317        }
318      }//end of bin loop
319      limits.push_back(rows);
320      vlimits.push_back(vrows);
321    }//end of sample loop
322 <  dout << endl << endl << "PAS table 3: " << endl << endl;
322 >
323 >  dout << endl << endl << "_______________________________________________________________________________________" << endl;
324 >  dout << "Going to store upper limit on event yield in result library: " << endl;
325 >  for(int ibin=0;ibin<jzbcuts.size();ibin++) {
326 >      int resultindex=PlottingSetup::allresults.Find(jzbcuts[ibin]);
327 >      vector<float> Normsigmas = compute_one_upper_limit(1.0,0.0, resultindex, mcjzb, "UPPERLIMIT", false, 0);
328 >      (allresults.predictions[resultindex]).UpperLimit=Normsigmas[0]*PlottingSetup::luminosity;
329 >  }
330 >  dout << endl << "_______________________________________________________________________________________" << endl;
331 >
332 >  dout << endl << endl << endl << "_________________________________________________________________________________________________" << endl << endl;
333 >  dout << endl << endl << "PAS table 3:   (notation: limit [95%CL])" << endl << endl;
334    dout << "\t";
335    for (int irow=0;irow<jzbcuts.size();irow++) {
336      dout << jzbcuts[irow] << "\t";
# Line 231 | Line 338 | void compute_upper_limits_from_counting_
338    dout << endl;
339    for(int irow=0;irow<limits.size();irow++) {
340      for(int ientry=0;ientry<limits[irow].size();ientry++) {
341 <      dout << limits[irow][ientry] << "\t";
341 >      if (limits[irow][ientry]>0) dout << limits[irow][ientry] << "\t";
342 >      else dout << " (N/A) \t";
343      }
344      dout << endl;
345    }
346    
347 <  if(!doobserved) {
348 <    dout << endl << endl << "LIMITS: " << endl;
349 <    dout << "\t";
347 >  if(!doexpected) {
348 >    dout << endl << endl << "LIMITS: (Tex)" << endl;
349 >    tout << "\\begin{table}[hbtp]" << endl;
350 >    tout << "\\renewcommand{\\arraystretch}{1.3}" << endl;
351 >    tout << "\\begin{center}" << endl;
352 >    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;
353 >    tout << "" << endl;
354 >    tout << "\\begin{tabular}{ | l | ";
355 >    for (int irow=0;irow<jzbcuts.size();irow++) tout << " l |";
356 >    tout << "} " << endl << " \\hline " << endl << "& \t ";
357      for (int irow=0;irow<jzbcuts.size();irow++) {
358 <      dout << jzbcuts[irow] << "\t";
358 >      tout << "JZB $>$ " << jzbcuts[irow] << " GeV & \t ";
359      }
360 <    dout << endl;
360 >    tout << " \\\\ \\hline " << endl;
361      for(int irow=0;irow<limits.size();irow++) {
362 <      dout << limits[irow][0] << "\t";
362 >      tout << limits[irow][0] << " \t";
363        for(int ientry=0;ientry<jzbcuts.size();ientry++) {
364 <        dout << Round(vlimits[irow][2*ientry] / vlimits[irow][2*ientry+1],3)<< "\t";
364 >        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";
365 >        else tout << " & ( N / A ) \t";
366 > //      dout << Round(vlimits[irow][2*ientry],3) << " / " << Round(vlimits[irow][2*ientry+1],3)<< "\t";
367        }
368 <      dout << endl;
368 >      tout << " \\\\ \\hline " << endl;
369      }
370 +      tout << "\\end{tabular}" << endl;
371 +      tout << "      \\end{tabular}"<< endl;
372 +      tout << "\\end{center}"<< endl;
373 +      tout << "\\end{table} "<< endl;
374 +
375    }//do observed
376    
377    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;
378 <  dout << "Scenario \t Efficiency [%] \t Upper limits [pb] \t Prediction [pb]" << endl;
378 >  dout << "Scenario \t Efficiency [%] \t Upper limits [pb] \t \\sigma [pb]" << endl;
379    for(int icut=0;icut<jzbcuts.size();icut++) {
380 <    dout << "Region with JZB>" << jzbcuts[icut] << endl;
380 >    dout << "Region with JZB>" << jzbcuts[icut] << (ConsiderSignalContaminationForLimits?"  (accounting for signal contamination)":"  (not accounting for signal contamination)") << endl;
381      for(int isample=0;isample<signalsamples.collection.size();isample++) {
382 <      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;
382 >      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;
383      }
384      dout << endl;
385    }
386 +  allresults.Print();
387 +
388    
265 –  write_warning("compute_upper_limits_from_counting_experiment","Still need to update the script");
389   }
390  
268 – void susy_scan_axis_labeling(TH2F *histo) {
269 –  histo->GetXaxis()->SetTitle("#Chi_{2}^{0}-LSP");
270 –  histo->GetXaxis()->CenterTitle();
271 –  histo->GetYaxis()->SetTitle("m_{#tilde{q}}");
272 –  histo->GetYaxis()->CenterTitle();
273 – }
391  
392 < void scan_susy_space(string mcjzb, string datajzb) {
393 <  TCanvas *c3 = new TCanvas("c3","c3");
277 <  vector<float> binning;
278 <  binning=allsamples.get_optimal_binsize(mcjzb,cutmass&&cutOSSF&&cutnJets,20,50,800);
279 <  float arrbinning[binning.size()];
280 <  for(int i=0;i<binning.size();i++) arrbinning[i]=binning[i];
281 <  TH1F *puredata   = allsamples.Draw("puredata",  datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
282 <  puredata->SetMarkerSize(DataMarkerSize);
283 <  TH1F *allbgs   = allsamples.Draw("allbgs",  "-"+datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
284 <  TH1F *allbgsb   = allsamples.Draw("allbgsb",  "-"+datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data,luminosity);
285 <  TH1F *allbgsc   = allsamples.Draw("allbgsc",  datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data,luminosity);
286 <  allbgs->Add(allbgsb,-1);
287 <  allbgs->Add(allbgsc);
288 <  int ndata=puredata->Integral();
289 <  ofstream myfile;
290 <  myfile.open ("susyscan_log.txt");
291 <  TFile *susyscanfile = new TFile("/scratch/fronga/SMS/T5z_SqSqToQZQZ_38xFall10.root");
292 <  TTree *suevents = (TTree*)susyscanfile->Get("events");
293 <  TH2F *exclusionmap = new TH2F("exclusionmap","",20,0,500,20,0,1000);
294 <  TH2F *exclusionmap1s = new TH2F("exclusionmap1s","",20,0,500,20,0,1000);
295 <  TH2F *exclusionmap2s = new TH2F("exclusionmap2s","",20,0,500,20,0,1000);
296 <  TH2F *exclusionmap3s = new TH2F("exclusionmap3s","",20,0,500,20,0,1000);
297 <  
298 <  susy_scan_axis_labeling(exclusionmap);
299 <  susy_scan_axis_labeling(exclusionmap1s);
300 <  susy_scan_axis_labeling(exclusionmap2s);
301 <  susy_scan_axis_labeling(exclusionmap3s);
302 <  
303 <  Int_t MyPalette[100];
304 <  Double_t r[]    = {0., 0.0, 1.0, 1.0, 1.0};
305 <  Double_t g[]    = {0., 0.0, 0.0, 1.0, 1.0};
306 <  Double_t b[]    = {0., 1.0, 0.0, 0.0, 1.0};
307 <  Double_t stop[] = {0., .25, .50, .75, 1.0};
308 <  Int_t FI = TColor::CreateGradientColorTable(5, stop, r, g, b, 100);
309 <  for (int i=0;i<100;i++) MyPalette[i] = FI+i;
310 <  
311 <  gStyle->SetPalette(100, MyPalette);
312 <  
313 <  for(int m23=50;m23<500;m23+=25) {
314 <    for (int m0=(2*(m23-50)+150);m0<=1000;m0+=50)
315 <    {
316 <      c3->cd();
317 <      stringstream drawcondition;
318 <      drawcondition << "pfJetGoodNum>=3&&(TMath::Abs(masses[0]-"<<m0<<")<10&&TMath::Abs(masses[2]-masses[3]-"<<m23<<")<10)&&mll>5&&id1==id2";
319 <      TH1F *puresignal = new TH1F("puresignal","puresignal",binning.size()-1,arrbinning);
320 <      TH1F *puresignall= new TH1F("puresignall","puresignal",binning.size()-1,arrbinning);
321 <      stringstream drawvar,drawvar2;
322 <      drawvar<<mcjzb<<">>puresignal";
323 <      drawvar2<<"-"<<mcjzb<<">>puresignall";
324 <      suevents->Draw(drawvar.str().c_str(),drawcondition.str().c_str());
325 <      suevents->Draw(drawvar2.str().c_str(),drawcondition.str().c_str());
326 <      if(puresignal->Integral()<60) {
327 <        delete puresignal;
328 <        continue;
329 <      }
330 <      puresignal->Add(puresignall,-1);//we need to correct for the signal contamination - we effectively only see (JZB>0)-(JZB<0) !!
331 <      puresignal->Scale(ndata/(20*puresignal->Integral()));//normalizing it to 5% of the data
332 <      stringstream saveas;
333 <      saveas<<"Model_Scan/m0_"<<m0<<"__m23_"<<m23;
334 <      dout << "PLEASE KEEP IN MIND THAT SIGNAL CONTAMINATION IS NOT REALLY TAKEN CARE OF YET DUE TO LOW STATISTICS! SHOULD BE SOMETHING LIKE THIS : "<< endl;
335 < //        TH1F *signalpredlo   = allsamples.Draw("signalpredlo",  "-"+mcjzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc,  luminosity,allsamples.FindSample("LM4"));
336 < //        TH1F *signalpredro   = allsamples.Draw("signalpredro",      mcjzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc,  luminosity,allsamples.FindSample("LM4"));
337 < //        TH1F *puredata       = allsamples.Draw("puredata",          datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
338 < //        signalpred->Add(signalpredlo,-1);
339 < //        signalpred->Add(signalpredro);
340 < //        puresignal->Add(signalpred,-1);//subtracting signal contamination
341 < //---------------------
342 < //      dout << "(m0,m23)=("<<m0<<","<<m23<<") contains " << puresignal->Integral() << endl;
343 < //    TH1F *puresignal = allsamples.Draw("puresignal",mcjzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,mc,  luminosity,allsamples.FindSample("LM4"));
344 <      vector<float> results=establish_upper_limits(puredata,allbgs,puresignal,saveas.str(),myfile);  
345 <      if(results.size()==0) {
346 <        delete puresignal;
347 <        continue;
348 <      }
349 <      exclusionmap->Fill(m23,m0,results[0]);
350 <      exclusionmap1s->Fill(m23,m0,results[1]);
351 <      exclusionmap2s->Fill(m23,m0,results[2]);
352 <      exclusionmap3s->Fill(m23,m0,results[3]);
353 <      delete puresignal;
354 <      dout << "(m0,m23)=("<<m0<<","<<m23<<") : 3 sigma at " << results[3] << endl;
355 <    }
356 <  }//end of model scan for loop
357 <  
358 <  dout << "Exclusion Map contains" << exclusionmap->Integral() << " (integral) and entries: " << exclusionmap->GetEntries() << endl;
359 <  c3->cd();
360 <  exclusionmap->Draw("CONTZ");
361 <  CompleteSave(c3,"Model_Scan/CONT/Model_Scan_Mean_values");
362 <  exclusionmap->Draw("COLZ");
363 <  CompleteSave(c3,"Model_Scan/COL/Model_Scan_Mean_values");
364 <  
365 <  exclusionmap1s->Draw("CONTZ");
366 <  CompleteSave(c3,"Model_Scan/CONT/Model_Scan_1sigma_values");
367 <  exclusionmap1s->Draw("COLZ");
368 <  CompleteSave(c3,"Model_Scan/COL/Model_Scan_1sigma_values");
369 <  
370 <  exclusionmap2s->Draw("CONTZ");
371 <  CompleteSave(c3,"Model_Scan/CONT/Model_Scan_2sigma_values");
372 <  exclusionmap2s->Draw("COLZ");
373 <  CompleteSave(c3,"Model_Scan/COL/Model_Scan_2sigma_values");
374 <  
375 <  exclusionmap3s->Draw("CONTZ");
376 <  CompleteSave(c3,"Model_Scan/CONT/Model_Scan_3sigma_values");
377 <  exclusionmap3s->Draw("COLZ");
378 <  CompleteSave(c3,"Model_Scan/COL/Model_Scan_3sigma_values");
379 <  
380 <  TFile *exclusion_limits = new TFile("exclusion_limits.root","RECREATE");
381 <  exclusionmap->Write();
382 <  exclusionmap1s->Write();
383 <  exclusionmap2s->Write();
384 <  exclusionmap3s->Write();
385 <  exclusion_limits->Close();
386 <  susyscanfile->Close();
387 <  
388 <  myfile.close();
389 < }
392 >
393 > /********************************************************************** new : Limits using SHAPES ***********************************
394  
395  
396 +   SSSSSSSSSSSSSSS hhhhhhh                                                                                      
397 + SS:::::::::::::::Sh:::::h                                                                                      
398 + S:::::SSSSSS::::::Sh:::::h                                                                                      
399 + S:::::S     SSSSSSSh:::::h                                                                                      
400 + S:::::S             h::::h hhhhh         aaaaaaaaaaaaa  ppppp   ppppppppp       eeeeeeeeeeee        ssssssssss  
401 + S:::::S             h::::hh:::::hhh      a::::::::::::a p::::ppp:::::::::p    ee::::::::::::ee    ss::::::::::s  
402 + S::::SSSS          h::::::::::::::hh    aaaaaaaaa:::::ap:::::::::::::::::p  e::::::eeeee:::::eess:::::::::::::s
403 +  SS::::::SSSSS     h:::::::hhh::::::h            a::::app::::::ppppp::::::pe::::::e     e:::::es::::::ssss:::::s
404 +    SSS::::::::SS   h::::::h   h::::::h    aaaaaaa:::::a p:::::p     p:::::pe:::::::eeeee::::::e s:::::s  ssssss
405 +       SSSSSS::::S  h:::::h     h:::::h  aa::::::::::::a p:::::p     p:::::pe:::::::::::::::::e    s::::::s      
406 +            S:::::S h:::::h     h:::::h a::::aaaa::::::a p:::::p     p:::::pe::::::eeeeeeeeeee        s::::::s  
407 +            S:::::S h:::::h     h:::::ha::::a    a:::::a p:::::p    p::::::pe:::::::e           ssssss   s:::::s
408 + SSSSSSS     S:::::S h:::::h     h:::::ha::::a    a:::::a p:::::ppppp:::::::pe::::::::e          s:::::ssss::::::s
409 + S::::::SSSSSS:::::S h:::::h     h:::::ha:::::aaaa::::::a p::::::::::::::::p  e::::::::eeeeeeee  s::::::::::::::s
410 + S:::::::::::::::SS  h:::::h     h:::::h a::::::::::aa:::ap::::::::::::::pp    ee:::::::::::::e   s:::::::::::ss  
411 + SSSSSSSSSSSSSSS    hhhhhhh     hhhhhhh  aaaaaaaaaa  aaaap::::::pppppppp        eeeeeeeeeeeeee    sssssssssss    
412 +                                                         p:::::p                                                
413 +                                                         p:::::p                                                
414 +                                                        p:::::::p                                                
415 +                                                        p:::::::p                                                
416 +                                                        p:::::::p                                                
417 +                                                        ppppppppp                                                
418 +                                                                                                                
419 +
420 + *********************************************************************** new : Limits using SHAPES ***********************************/
421  
393 –
394 – //********************************************************************** new : Limits using SHAPES ***********************************
422  
423   void limit_shapes_for_systematic_effect(TFile *limfile, string identifier, string mcjzb, string datajzb, int JES,vector<float> binning, TCanvas *limcan) {
424    dout << "Creatig shape templates ... ";
425    if(identifier!="") dout << "for systematic called "<<identifier;
426    dout << endl;
427    int dataormc=mcwithsignal;//this is only for tests - for real life you want dataormc=data !!!
428 <  if(dataormc!=data) write_warning("limit_shapes_for_systematic_effect","WATCH OUT! Not using data for limits!!!! this is ok for tests, but not ok for anything official!");
428 >  if(dataormc!=data) write_warning(__FUNCTION__,"WATCH OUT! Not using data for limits!!!! this is ok for tests, but not ok for anything official!");
429    
430    TCut limitnJetcut;
431    if(JES==noJES) limitnJetcut=cutnJets;
# Line 411 | Line 438 | void limit_shapes_for_systematic_effect(
438    TH1F *ZOSSFN = allsamples.Draw("ZOSSFN","-"+datajzb,binning, "JZB4limits", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,dataormc,luminosity);
439    TH1F *ZOSOFN = allsamples.Draw("ZOSOFN","-"+datajzb,binning, "JZB4limits", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,dataormc,luminosity);
440    
441 <  TH1F *SBOSSFP = allsamples.Draw("SBOSSFP",datajzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
442 <  TH1F *SBOSOFP = allsamples.Draw("SBOSOFP",datajzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
443 <  TH1F *SBOSSFN = allsamples.Draw("SBOSSFN","-"+datajzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
444 <  TH1F *SBOSOFN = allsamples.Draw("SBOSOFN","-"+datajzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
441 >  TH1F *SBOSSFP;
442 >  TH1F *SBOSOFP;
443 >  TH1F *SBOSSFN;
444 >  TH1F *SBOSOFN;
445    
446    TH1F *LZOSSFP = allsamples.Draw("LZOSSFP",mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("LM4"));
447    TH1F *LZOSOFP = allsamples.Draw("LZOSOFP",mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("LM4"));
448    TH1F *LZOSSFN = allsamples.Draw("LZOSSFN","-"+mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("LM4"));
449    TH1F *LZOSOFN = allsamples.Draw("LZOSOFN","-"+mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("LM4"));
450    
451 <  TH1F *LSBOSSFP = allsamples.Draw("LSBOSSFP",mcjzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
452 <  TH1F *LSBOSOFP = allsamples.Draw("LSBOSOFP",mcjzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
453 <  TH1F *LSBOSSFN = allsamples.Draw("LSBOSSFN","-"+mcjzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
454 <  TH1F *LSBOSOFN = allsamples.Draw("LSBOSOFN","-"+mcjzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
451 >  TH1F *LSBOSSFP;
452 >  TH1F *LSBOSOFP;
453 >  TH1F *LSBOSSFN;
454 >  TH1F *LSBOSOFN;
455 >  
456 >  flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
457 >  if(PlottingSetup::RestrictToMassPeak) {
458 >      SBOSSFP = allsamples.Draw("SBOSSFP",datajzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
459 >      SBOSOFP = allsamples.Draw("SBOSOFP",datajzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
460 >      SBOSSFN = allsamples.Draw("SBOSSFN","-"+datajzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
461 >      SBOSOFN = allsamples.Draw("SBOSOFN","-"+datajzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
462 >      
463 >      LSBOSSFP = allsamples.Draw("LSBOSSFP",mcjzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
464 >      LSBOSOFP = allsamples.Draw("LSBOSOFP",mcjzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
465 >      LSBOSSFN = allsamples.Draw("LSBOSSFN","-"+mcjzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
466 >      LSBOSOFN = allsamples.Draw("LSBOSOFN","-"+mcjzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
467 >  }
468    
469    string obsname="data_obs";
470    string predname="background";
# Line 435 | Line 475 | void limit_shapes_for_systematic_effect(
475      signalname="signal_"+identifier;
476    }
477    
478 <  TH1F *obs = (TH1F*)ZOSSFP->Clone();
478 >  TH1F *obs = (TH1F*)ZOSSFP->Clone("observation");
479    obs->SetName(obsname.c_str());
480    obs->Write();
481 <  TH1F *pred = (TH1F*)ZOSSFN->Clone();
482 <  pred->Add(ZOSOFP,1.0/3);
483 <  pred->Add(ZOSOFN,-1.0/3);
484 <  pred->Add(SBOSSFP,1.0/3);
485 <  pred->Add(SBOSSFN,-1.0/3);
486 <  pred->Add(SBOSOFP,1.0/3);
487 <  pred->Add(SBOSOFN,-1.0/3);
481 >  TH1F *pred = (TH1F*)ZOSSFN->Clone("prediction");
482 >    flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
483 >  if(PlottingSetup::RestrictToMassPeak) {
484 >    pred->Add(ZOSOFP,1.0/3);
485 >    pred->Add(ZOSOFN,-1.0/3);
486 >    pred->Add(SBOSSFP,1.0/3);
487 >    pred->Add(SBOSSFN,-1.0/3);
488 >    pred->Add(SBOSOFP,1.0/3);
489 >    pred->Add(SBOSOFN,-1.0/3);
490 >  } else {
491 >    pred->Add(ZOSOFP,1.0);
492 >    pred->Add(ZOSOFN,-1.0);
493 >  }
494 >    
495    pred->SetName(predname.c_str());
496    pred->Write();
497    
# Line 455 | Line 502 | void limit_shapes_for_systematic_effect(
502    TH1F *Lpred = new TH1F("Lpred","Lpred",binning.size()-1,&binning[0]);
503    Lobs->Add(LZOSSFP);
504    Lpred->Add(LZOSSFN);
505 <  Lpred->Add(LZOSOFP,1.0/3);
506 <  Lpred->Add(LZOSOFN,-1.0/3);
507 <  Lpred->Add(LSBOSSFP,1.0/3);
508 <  Lpred->Add(LSBOSSFN,-1.0/3);
509 <  Lpred->Add(LSBOSOFP,1.0/3);
510 <  Lpred->Add(LSBOSOFN,-1.0/3);
505 >    flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
506 >  if(PlottingSetup::RestrictToMassPeak) {
507 >    Lpred->Add(LZOSOFP,1.0/3);
508 >    Lpred->Add(LZOSOFN,-1.0/3);
509 >    Lpred->Add(LSBOSSFP,1.0/3);
510 >    Lpred->Add(LSBOSSFN,-1.0/3);
511 >    Lpred->Add(LSBOSOFP,1.0/3);
512 >    Lpred->Add(LSBOSOFN,-1.0/3);
513 >  } else {
514 >    Lpred->Add(LZOSOFP,1.0);
515 >    Lpred->Add(LZOSOFN,-1.0);
516 >  }
517 >
518    TH1F *signal = (TH1F*)Lobs->Clone();
519    signal->Add(Lpred,-1);
520    signal->SetName(signalname.c_str());
# Line 474 | Line 528 | void limit_shapes_for_systematic_effect(
528    delete ZOSSFN;
529    delete ZOSOFN;
530    
531 <  delete SBOSSFP;
532 <  delete SBOSOFP;
533 <  delete SBOSSFN;
534 <  delete SBOSOFN;
531 >  if(PlottingSetup::RestrictToMassPeak) {
532 >    delete SBOSSFP;
533 >    delete SBOSOFP;
534 >    delete SBOSSFN;
535 >    delete SBOSOFN;
536 >  }
537    
538    delete LZOSSFP;
539    delete LZOSOFP;
540    delete LZOSSFN;
541    delete LZOSOFN;
542    
543 <  delete LSBOSSFP;
544 <  delete LSBOSOFP;
545 <  delete LSBOSSFN;
546 <  delete LSBOSOFN;
543 >  if(PlottingSetup::RestrictToMassPeak) {
544 >    delete LSBOSSFP;
545 >    delete LSBOSOFP;
546 >    delete LSBOSSFN;
547 >    delete LSBOSOFN;
548 >  }
549  
550   }
551  
# Line 498 | Line 556 | void prepare_datacard(TFile *f) {
556  
557   ofstream datacard;
558   ensure_directory_exists(get_directory()+"/limits");
559 < datacard.open ((get_directory()+"/limits/susylm4datacard.txt").c_str());
559 > datacard.open ((get_directory()+"/limits/susydatacard.txt").c_str());
560   datacard << "Writing this to a file.\n";
561   datacard << "imax 1\n";
562   datacard << "jmax 1\n";
# Line 534 | Line 592 | void prepare_limits(string mcjzb, string
592    limit_shapes_for_systematic_effect(limfile,"JESDown",mcjzb,datajzb,JESdown,jzbbins,limcan);
593    
594    prepare_datacard(limfile);
537 –  
538 –  write_info("prepare_limits","limitfile.root and datacard.txt have been generated. You can now use them to calculate limits!");
595    limfile->Close();
596 +  write_info("prepare_limits","limitfile.root and datacard.txt have been generated. You can now use them to calculate limits!");
597    
598 < }
598 > }

Diff Legend

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