AliPhysics  a8afd6c (a8afd6c)
MakePurityCalculationAndComparisons.C
Go to the documentation of this file.
1 
12 #if !defined(__CINT__) || defined(__MAKECINT__)
13 
14 #include <TFile.h>
15 #include <TDirectoryFile.h>
16 #include <TList.h>
17 #include <TString.h>
18 #include <TROOT.h>
19 #include <TStyle.h>
20 #include <TH1D.h>
21 #include <TH2F.h>
22 #include <TCanvas.h>
23 #include <TPad.h>
24 #include <TLegend.h>
25 #include <TObject.h>
26 #include <TAxis.h>
27 #include <TGaxis.h>
28 #include <TLine.h>
29 #include <TF1.h>
30 #include <TMath.h>
31 #include <TGraphErrors.h>
32 #include "PlotUtils.C"
33 #endif
34 
35 Int_t color [] = {1,4,2,6,8,9,kYellow-6,kCyan,kOrange+2,kViolet,kOrange-2,4,2,6,8,9,kYellow-6,kCyan,kOrange+2,kViolet,kOrange-2};
36 Int_t lineStyle[] = {1,1,1,1,1,1, 1, 1, 1, 1, 1,2,2,2,2,2, 2, 2, 2, 2, 2,2,2,2,2,2,2,2,2,2,2,2};
37 Int_t marker [] = {24,21,20,25,24,24,24,24,24,24,24};
38 TString format = ".eps";
39 
40 // Main input data file location, set them in MakePurityCalculationAndComparison method
41 TString opt = "100MeV_New";
42 TString dataFileDir = "data3";
43 
44 //------------------------------------------------------------------------------
65 //------------------------------------------------------------------------------
66 void MakeDoubleRatios
67 (
68  TString titleProd, TString gjProd
69  , Int_t nCases, TString * mcCasesPath
70  , TString * mcCases, TString * mcLeg
71  , Bool_t addGJ = 1
72  , Float_t gjFactor = 1
73  , Bool_t usePi0 = 0
74  , Int_t rebin = 0
75  , Float_t m02LowCutMin = 0.10
76  , Float_t m02LowCutMax = 0.27//0.27
77  , Float_t m02HigCutMin = 0.4//0.40
78  , Float_t m02HigCutMax = 3.00
79  , Int_t debug = kFALSE
80 )
81 {
82  // Do some checks and modify stuff
83  //
84 
85  if ( addGJ && titleProd.Contains("LHC") )
86  {
87  if ( debug )
88  {
89  printf("++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n");
90  printf("In case of data, not possible to add GJ simulation, do nothing, deactivate this option\n");
91  printf("++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n");
92  }
93  return;
94  }
95 
96  // Create histograms arrays
97  //
98  const Int_t nMCCases = nCases;
99 
100  TFile * fileJJ[nMCCases];
101  TFile * fileGJ[nMCCases];
102  Int_t firstMC = 0;
103 
104  TH2D* h2M02Iso [nMCCases];
105  TH2D* h2M02Not [nMCCases];
106  TH2D* h2M02IsoGJ[nMCCases];
107 
108  TH1D* hM02LowIs [nMCCases];
109  TH1D* hM02LowNo [nMCCases];
110  TH1D* hM02HigIs [nMCCases];
111  TH1D* hM02HigNo [nMCCases];
112 
113  TH1D* hM02LowIsRat[nMCCases];
114  TH1D* hM02LowNoRat[nMCCases];
115  TH1D* hM02HigIsRat[nMCCases];
116  TH1D* hM02HigNoRat[nMCCases];
117 
118  TH1D* hDoubleRatioCases[nMCCases];
119  TH1D* hDoubleRatio[nMCCases];
120 
121  // Init the arrays, just in case
122  for(Int_t icase = 0; icase < nMCCases; icase++)
123  {
124  fileJJ[icase] = 0;
125  fileGJ[icase] = 0;
126 
127  hM02LowIs [icase]=0;
128  hM02LowNo [icase]=0;
129  hM02HigIs [icase]=0;
130  hM02HigNo [icase]=0;
131 
132  hM02LowIsRat[icase]=0;
133  hM02LowNoRat[icase]=0;
134  hM02HigIsRat[icase]=0;
135  hM02HigNoRat[icase]=0;
136 
137  hDoubleRatio[icase]=0;
138  hDoubleRatioCases[icase]=0;
139  }
140 
141  // Define energy binning
142  //
143  // Rebin (>0)
144  const Int_t nBins = 14;
145  Double_t binE[]={8,9,10,11,12,13,14,16,18,20,25,30,40,50,60};
146  // Histogram limits
147  Float_t emin = 8;
148  Float_t emax = 60;
149  if(titleProd =="Low")
150  {
151  emin = 8;
152  emax = 30;
153  }
154  if(titleProd =="High")
155  {
156  emin = 18;
157  emax = 60;
158  }
159 
160  // String used later to name the output files that depends on the parameters
161  // provided
162  TString cuts = "";
163  if ( usePi0 )
164  cuts = Form("Low_%2.2f-%2.2f" ,m02LowCutMin, m02LowCutMax);
165  else
166  cuts = Form("Low_%2.2f-%2.2f_High_%2.2f-%2.2f",m02LowCutMin, m02LowCutMax, m02HigCutMin, m02HigCutMax);
167 
168  TString prodTitleOutput = titleProd;
169  if ( addGJ ) prodTitleOutput+=Form("_GJx%2.2f",gjFactor);
170  if ( usePi0 ) prodTitleOutput+="_MergedPi0";
171 
172  //---------------------------------------
173  // Open the files and get the histograms
174  //---------------------------------------
175 
176  if ( debug )
177  printf("Get files and histograms \n");
178  // Simu
179  for(Int_t icase = firstMC; icase < nMCCases; icase++)
180  {
181  if ( debug )
182  printf("\t Case %d - %s\n",icase, mcCases[icase].Data());
183  if ( titleProd.Contains("LHC") )
184  {
185  fileJJ[icase] = TFile::Open(Form("%s/%s/%s.root",opt.Data(),dataFileDir.Data(),titleProd.Data()));
186  }
187  else
188  fileJJ[icase] = TFile::Open(Form("%s/%s/%s/ScaledMerged_%s.root",opt.Data(),mcCasesPath[icase].Data(),titleProd.Data(), mcCases[icase].Data()));
189 
190  if ( addGJ )
191  fileGJ[icase] = TFile::Open(Form("%s/%s/%s/ScaledMerged_%s.root",opt.Data(),mcCasesPath[icase].Data(),gjProd.Data(), mcCases[icase].Data()));
192 
193  if ( debug )
194  {
195  printf("Data/MC file %s %p\n",titleProd.Data(),fileJJ[icase]);
196  printf("GJ file %s %p\n",gjProd.Data(),fileGJ[icase]);
197  }
198 
199  if ( !fileJJ[icase] )
200  {
201  h2M02Iso[icase] = 0;
202  h2M02Not[icase] = 0;
203 
204  continue;
205  }
206 
207  //--------------------------------------------------------------------------
208  // Get M02 vs E histograms
209  //--------------------------------------------------------------------------
210 
211  h2M02Not[icase] = (TH2D*) fileJJ[icase]->Get("AnaIsolPhoton_hPtLambda0NoIso");
212  if ( addGJ )
213  h2M02Not[icase]->Add((TH2D*) fileGJ[icase]->Get("AnaIsolPhoton_hPtLambda0NoIso"),gjFactor);
214 
215  h2M02Iso[icase] = (TH2D*) fileJJ[icase]->Get("AnaIsolPhoton_hPtLambda0Iso");
216  // P3 does not need GJ in photon region, only in pi0 region
217  //if ( addGJ )
218  // h2M02Iso[icase]->Add((TH2D*) fileGJ[icase]->Get("AnaIsolPhoton_hPtLambda0Iso"),gjFactor);
219 
220  if ( addGJ )
221  h2M02IsoGJ[icase] = (TH2D*) fileGJ[icase]->Get("AnaIsolPhoton_hPtLambda0Iso");
222 
223 // h2M02Not[icase] = (TH2D*) fileJJ[icase]->Get("AnaIsolPhoton_TM1_hPtLambda0NoIso");
224 // if ( addGJ )
225 // h2M02Not[icase]->Add((TH2D*) fileGJ[icase]->Get("AnaIsolPhoton_TM1_hPtLambda0NoIso"),gjFactor);
226 // h2M02Iso[icase] = (TH2D*) fileJJ[icase]->Get("AnaIsolPhoton_TM1_hPtLambda0Iso");
227 // // P3 does not need GJ in photon region, only in pi0 region
228 // //if ( addGJ )
229 // // h2M02Iso[icase]->Add((TH2D*) fileGJ[icase]->Get("AnaIsolPhoton_TM1_hPtLambda0Iso"),gjFactor);
230 //
231 // if ( addGJ )
232 // h2M02IsoGJ[icase] = (TH2D*) fileGJ[icase]->Get("AnaIsolPhoton_TM1_hPtLambda0Iso");
233 
234  //--------------------------------------------------------------------------
235  // Project M02 vs E plot in signal region
236  //--------------------------------------------------------------------------
237 
238  Int_t binm02LowCutMin = h2M02Iso[icase]->GetYaxis()->FindBin(m02LowCutMin);
239  Int_t binm02LowCutMax = h2M02Iso[icase]->GetYaxis()->FindBin(m02LowCutMax)-1;
240 
241  if ( debug )
242  {
243  printf("Signal region: M02 [%2.2f,%2.2f] - Bins [%d,%d] \n",
244  m02LowCutMin,m02LowCutMax,binm02LowCutMin,binm02LowCutMax);
245  }
246 
247  hM02LowNo[icase] =
248  (TH1D*) h2M02Not[icase]->ProjectionX(Form("LowM02_NoIso_%s_%s",titleProd.Data(),cuts.Data()),
249  binm02LowCutMin,binm02LowCutMax);
250  // if(hM02LowNo[icase])
251  // printf("Entries %e Integral %e\n",hM02LowNo[icase]->GetEntries(),hM02LowNo[icase]->Integral());
252 
253  hM02LowIs[icase] =
254  (TH1D*) h2M02Iso[icase]->ProjectionX(Form("LowM02_Iso_%s_%s",titleProd.Data(),cuts.Data()),
255  binm02LowCutMin,binm02LowCutMax);
256 
257  // if(hM02LowIs[icase])
258  // printf("Entries %e Integral %e\n",hM02LowIs[icase]->GetEntries(),hM02LowIs[icase]->Integral());
259 
260  //--------------------------------------------------------------------------
261  // Project M02 vs E plot in bkg region,
262  // if usePi0=kTRUE, recover the corresponding histogram
263  //--------------------------------------------------------------------------
264 
265  if ( usePi0 )
266  {
267  hM02HigNo[icase] = (TH1D*) fileJJ[icase]->Get("AnaIsolPi0SS_hPtNoIso");
268  if(addGJ)hM02HigNo[icase]->Add((TH2D*) fileGJ[icase]->Get("AnaIsolPi0SS_hPtNoIso"),gjFactor);
269 
270  //if(hM02HigNo[icase])
271  // printf("Entries %e Integral %e\n",hM02HigNo[icase]->GetEntries(),hM02HigNo[icase]->Integral());
272 
273  hM02HigIs[icase] = (TH1D*) fileJJ[icase]->Get("AnaIsolPi0SS_hE");
274  if(addGJ)hM02HigIs[icase]->Add((TH2D*) fileGJ[icase]->Get("AnaIsolPi0SS_hE"),gjFactor);
275 
276  // if(hM02HigIs[icase])
277  // printf("Entries %e Integral %e\n",hM02HigIs[icase]->GetEntries(),hM02HigIs[icase]->Integral());
278 
279  //hM02HigNo[icase] = (TH1D*) fileJJ[icase]->Get("AnaIsolPi0SS_TM1_hPtNoIso");
280  //if(addGJ)hM02HigNo[icase]->Add((TH2D*) fileGJ[icase]->Get("AnaIsolPi0SS_TM1_hPtNoIso"),gjFactor);
281  //hM02HigIs[icase] = (TH1D*) fileJJ[icase]->Get("AnaIsolPi0SS_TM1_hE");
282  //if(addGJ)hM02HigIs[icase]->Add((TH2D*) fileGJ[icase]->Get("AnaIsolPi0SS_TM1_hE"),gjFactor);
283  }
284  else
285  {
286  // Project M02 vs E plot
287 
288  Int_t binm02HigCutMin = h2M02Iso[icase]->GetYaxis()->FindBin(m02HigCutMin);
289  Int_t binm02HigCutMax = h2M02Iso[icase]->GetYaxis()->FindBin(m02HigCutMax)-1;
290 
291  if ( debug )
292  {
293  printf("Bkg region: M02 [%2.2f,%2.2f] - Bins [%d,%d] \n",
294  m02HigCutMin,m02HigCutMax,binm02HigCutMin,binm02HigCutMax);
295  }
296 
297  hM02HigNo[icase] =
298  (TH1D*) h2M02Not[icase]->ProjectionX(Form("HigM02_NoIso_%s_%s",titleProd.Data(),cuts.Data()),
299  binm02HigCutMin,binm02HigCutMax);
300  //if(hM02HigNo[icase])
301  // printf("Entries %e Integral %e\n",hM02HigNo[icase]->GetEntries(),hM02HigNo[icase]->Integral());
302 
303  hM02HigIs[icase] =
304  (TH1D*) h2M02Iso[icase]->ProjectionX(Form("HigM02_Iso_%s_%s",titleProd.Data(),cuts.Data()),
305  binm02HigCutMin,binm02HigCutMax);
306 
307  // P3 needs GJ on iso bkg photon region not in signal
308  if ( addGJ )
309  {
310  hM02HigIs[icase] ->Add
311  ( (TH1D*) h2M02IsoGJ[icase]->ProjectionX(Form("HigM02_GJ_Iso_%s_%s",titleProd.Data(),cuts.Data()),
312  binm02HigCutMin,binm02HigCutMax)
313  ,gjFactor);
314  }
315  //if(hM02HigIs[icase])
316  // printf("Entries %e Integral %e\n",hM02HigIs[icase]->GetEntries(),hM02HigIs[icase]->Integral());
317  }
318 
319  // Set histogram ranges and titles
320  //
321  hM02HigIs[icase]->SetMarkerStyle(marker[1]);
322  hM02HigIs[icase]->SetMarkerColor(color[icase]);
323  hM02HigIs[icase]->SetLineColor (color[icase]);
324  hM02HigIs[icase]->SetAxisRange(emin,emax,"X");
325  hM02HigIs[icase]->SetYTitle("d #sigma / d #it{p}_{T} (mb)");
326 
327  hM02HigNo[icase]->SetMarkerStyle(marker[3]);
328  hM02HigNo[icase]->SetMarkerColor(color[icase]);
329  hM02HigNo[icase]->SetLineColor (color[icase]);
330  hM02HigNo[icase]->SetAxisRange(emin,emax,"X");
331  hM02HigNo[icase]->SetYTitle("d #sigma / d #it{p}_{T} (mb)");
332 
333  hM02LowIs[icase]->SetMarkerStyle(marker[0]);
334  hM02LowIs[icase]->SetMarkerColor(color[icase]);
335  hM02LowIs[icase]->SetLineColor (color[icase]);
336  hM02LowIs[icase]->SetAxisRange(emin,emax,"X");
337  hM02LowIs[icase]->SetYTitle("d #sigma / d #it{p}_{T} (mb)");
338 
339  hM02LowNo[icase]->SetMarkerStyle(marker[2]);
340  hM02LowNo[icase]->SetMarkerColor(color[icase]);
341  hM02LowNo[icase]->SetLineColor (color[icase]);
342  hM02LowNo[icase]->SetAxisRange(emin,emax,"X");
343  hM02LowNo[icase]->SetYTitle("d #sigma / d #it{p}_{T} (mb)");
344 
345  hM02LowNo[icase]->SetTitle(Form("C: %2.2f < #sigma_{long} < %2.2f - #Sigma #it{p}_{T} > 1 GeV/c",
346  m02LowCutMin,m02LowCutMax));
347  hM02LowNo[icase]->SetTitleOffset(1.6,"Y");
348 
349  hM02LowIs[icase]->SetTitle(Form("A: %2.2f < #sigma_{long} < %2.2f - #Sigma #it{p}_{T} < 1 GeV/c",
350  m02LowCutMin,m02LowCutMax));
351  hM02LowIs[icase]->SetTitleOffset(1.6,"Y");
352 
353  if(usePi0)
354  {
355  hM02HigIs[icase]->SetTitle("B: Merged #pi^{0} - #Sigma #it{p}_{T} < 1 GeV/c");
356  hM02HigIs[icase]->SetTitleOffset(1.6,"Y");
357 
358  hM02HigNo[icase]->SetTitle("D: Merged #pi^{0} - #Sigma #it{p}_{T} > 1 GeV/c");
359  hM02HigNo[icase]->SetTitleOffset(1.6,"Y");
360  }
361  else
362  {
363  hM02HigNo[icase]->SetTitle(Form("D: %2.2f < #sigma_{long} < %2.2f - #Sigma #it{p}_{T} > 1 GeV/c",
364  m02HigCutMin,m02HigCutMax));
365  hM02HigNo[icase]->SetTitleOffset(1.6,"Y");
366 
367  hM02HigIs[icase]->SetTitle(Form("B: %2.2f < #sigma_{long} < %2.2f - #Sigma #it{p}_{T} < 1 GeV/c",
368  m02HigCutMin,m02HigCutMax));
369  hM02HigIs[icase]->SetTitleOffset(1.6,"Y");
370  }
371 
372  if ( titleProd.Contains("LHC") )
373  {
374  hM02LowIs[icase]->Sumw2();
375  hM02HigIs[icase]->Sumw2();
376  hM02LowNo[icase]->Sumw2();
377  hM02HigNo[icase]->Sumw2();
378  }
379 
380  // Rebin
381  if(rebin > 0)
382  {
383  hM02LowIs[icase]->Rebin(rebin);
384  hM02HigIs[icase]->Rebin(rebin);
385  hM02LowNo[icase]->Rebin(rebin);
386  hM02HigNo[icase]->Rebin(rebin);
387  }
388  else
389  {
390  //printf("Rebin with n bins %d\n",nBins);
391  hM02LowIs[icase] = (TH1D*) hM02LowIs[icase]->Rebin(nBins,Form("%s_Rebin",hM02LowIs[icase]->GetName()),binE);
392  hM02HigIs[icase] = (TH1D*) hM02HigIs[icase]->Rebin(nBins,Form("%s_Rebin",hM02HigIs[icase]->GetName()),binE);
393  hM02LowNo[icase] = (TH1D*) hM02LowNo[icase]->Rebin(nBins,Form("%s_Rebin",hM02LowNo[icase]->GetName()),binE);
394  hM02HigNo[icase] = (TH1D*) hM02HigNo[icase]->Rebin(nBins,Form("%s_Rebin",hM02HigNo[icase]->GetName()),binE);
395 
396  if ( debug )
397  {
398  printf("Rebinned ABCD names: %s, %s,\n %s, %s\n",
399  hM02LowIs[icase]->GetName(),hM02HigIs[icase]->GetName(),
400  hM02LowNo[icase]->GetName(),hM02HigNo[icase]->GetName());
401  }
402 
403  // From PlotUtils.C
404  ScaleBinBySize(hM02LowIs[icase]);
405  ScaleBinBySize(hM02LowNo[icase]);
406  ScaleBinBySize(hM02HigIs[icase]);
407  ScaleBinBySize(hM02HigNo[icase]);
408  }
409  }
410 
411  //
412  // Plot
413  //
414 
415  if ( debug )
416  printf("PLOT Double Ratios \n");
417 
418  TString fileName ;
419 
420  gStyle->SetPadRightMargin(0.01);
421  gStyle->SetPadLeftMargin(0.12);
422  gStyle->SetPadTopMargin(0.08);
423 
424  gStyle->SetTitleFontSize(0.06);
425  // gStyle->SetStatFontSize(0.06);
426  gStyle->SetOptStat(0);
427  gStyle->SetOptTitle(1);
428 
430  TCanvas * cABCD = new TCanvas(Form("cABCD_%s_%s" ,cuts.Data(),prodTitleOutput.Data()),
431  Form("ABCD, %s, %s",cuts.Data(),prodTitleOutput.Data()),1000,1000);
432 
433  cABCD->Divide(2,2);
434 
435  TLegend *lABCD = new TLegend(0.3,0.6,0.89,0.89);
436  lABCD->SetFillColor(0);
437  lABCD->SetFillStyle(0);
438  lABCD->SetLineColor(0);
439  lABCD->SetBorderSize(0);
440  lABCD->SetTextSize(0.05);
441 
443  cABCD->cd(1);
444 
445  gPad->SetLogy();
446 
447  for(Int_t icase = 0; icase < nMCCases; icase++)
448  {
449  //printf("icase %d %p\n",icase,hM02LowIs[icase]);
450  if(!hM02LowIs[icase]) continue;
451 
452  if ( icase==0 ) hM02LowIs[icase]->Draw("H");
453  else hM02LowIs[icase]->Draw("H same");
454 
455 // if(hM02LowIs[icase]->GetMaximum() > hM02LowIs[0]->GetMaximum())
456 // hM02LowIs[0]->SetMaximum(hM02LowIs[icase]->GetMaximum()*1.2);
457 
458  lABCD->AddEntry(hM02LowIs[icase],Form("%s",mcLeg[icase].Data()),"L");
459  }
460 
461  if ( nCases > 1 )
462  lABCD->Draw();
463 
465 
466  cABCD->cd(2);
467 
468  gPad->SetLogy();
469 
470  for(Int_t icase = 0; icase < nMCCases; icase++)
471  {
472  if(!hM02HigIs[icase]) continue;
473 
474  if ( icase==0 ) hM02HigIs[icase]->Draw("H");
475  else hM02HigIs[icase]->Draw("H same");
476 
477  }
478 
479  if ( nCases > 1 )
480  lABCD->Draw();
481 
483  cABCD->cd(3);
484 
485  gPad->SetLogy();
486 
487  for(Int_t icase = 0; icase < nMCCases; icase++)
488  {
489  if(!hM02LowNo[icase]) continue;
490 
491  if ( icase==0 ) hM02LowNo[icase]->Draw("H");
492  else hM02LowNo[icase]->Draw("H same");
493 
494  // if(hM02LowNo[icase]->GetMaximum() > hM02LowNo[0]->GetMaximum())
495  // hM02LowNo[0]->SetMaximum(hM02LowNo[icase]->GetMaximum()*1.2);
496  }
497 
498  if ( nCases > 1 )
499  lABCD->Draw();
500 
502 
503  cABCD->cd(4);
504 
505  gPad->SetLogy();
506 
507  for(Int_t icase = 0; icase < nMCCases; icase++)
508  {
509  if(!hM02HigNo[icase]) continue;
510 
511  if ( icase==0 ) hM02HigNo[icase]->Draw("H");
512  else hM02HigNo[icase]->Draw("H same");
513 
514  // if(hM02HigNo[icase]->GetMaximum() > hM02HigNo[0]->GetMaximum())
515  // hM02HigNo[0]->SetMaximum(hM02HigNo[icase]->GetMaximum()*1.2);
516  }
517 
518  if ( nCases > 1 )
519  lABCD->Draw();
520 
521  fileName = Form("%s/figures/ABCD_%s_pT_%s",opt.Data(),cuts.Data(),prodTitleOutput.Data());
522  fileName+=".eps";
523  cABCD->Print(fileName);
524 
526  if ( nMCCases > 1 )
527  {
528  TCanvas * cABCDRatio = new TCanvas(Form("cABCDRatio_%s_%s" ,cuts.Data(),prodTitleOutput.Data()),
529  Form("ABCD cases Ratio, %s, %s",cuts.Data(),prodTitleOutput.Data()),
530  1000,1000);
531 
532  cABCDRatio->Divide(2,2);
533 
534  TLegend *lABCDRat = new TLegend(0.3,0.6,0.89,0.89);
535  lABCDRat->SetFillColor(0);
536  lABCDRat->SetFillStyle(0);
537  lABCDRat->SetLineColor(0);
538  lABCDRat->SetBorderSize(0);
539  lABCDRat->SetTextSize(0.05);
540  //lABCDRat->SetHeader(Form("Den.: %s",mcLeg[nMCCases-1].Data()));
541  lABCDRat->SetHeader(Form("Den.: %s",mcLeg[0].Data()));
542 
544  cABCDRatio->cd(1);
545  gPad->SetGridy();
546 
547  //for(Int_t icase = 0; icase < nMCCases-1; icase++)
548  for(Int_t icase = 1; icase < nMCCases; icase++)
549  {
550  if(!hM02LowIs[icase]) continue;
551 
552  TString cloneName = Form("M02LowIsCasesRat%s_%s_%s",mcCases[icase].Data(),mcCasesPath[icase].Data(),cuts.Data());
553  // Careful, in case of data, take the data name and not the MC case name
554  if ( titleProd.Contains("LHC") ) cloneName = Form("M02LowIsCasesRat%s_%s_%s",titleProd.Data(),"",cuts.Data());
555 
556  //printf(">>>>>>>>>>>>>> Clone name %s\n",cloneName.Data());
557 
558  hM02LowIsRat[icase] = (TH1D*) hM02LowIs[icase]->Clone(cloneName);
559  //hM02LowIsRat[icase]->Divide(hM02LowIs[nMCCases-1]);
560  hM02LowIsRat[icase]->Divide(hM02LowIs[0]);
561 
562  hM02LowIsRat[icase]->SetYTitle("Ratio");
563  hM02LowIsRat[icase]->SetMaximum(2.5);
564  hM02LowIsRat[icase]->SetMinimum(0.5);
565 
566  if(icase == 0) hM02LowIsRat[icase]->Draw("H");
567  else hM02LowIsRat[icase]->Draw("H same");
568 
569  lABCDRat->AddEntry(hM02LowIsRat[icase],Form("Num.:%s",mcLeg[icase].Data()),"L");
570  }
571 
572  lABCDRat->Draw();
573 
575  cABCDRatio->cd(2);
576  gPad->SetGridy();
577 
578  for(Int_t icase = 1; icase < nMCCases; icase++)
579  //for(Int_t icase = 0; icase < nMCCases-1; icase++)
580  {
581  if(!hM02HigIs[icase]) continue;
582 
583  TString cloneName = Form("M02HighIsCasesRat%s_%s_%s",mcCases[icase].Data(),mcCasesPath[icase].Data(),cuts.Data());
584  // Careful, in case of data, take the data name and not the MC case name
585  if ( titleProd.Contains("LHC") ) cloneName = Form("M02HighIsCasesRat%s_%s_%s",titleProd.Data(),"",cuts.Data());
586 
587  hM02HigIsRat[icase] = (TH1D*) hM02HigIs[icase]->Clone(cloneName);
588  //hM02HigIsRat[icase]->Divide(hM02HigIs[nMCCases-1]);
589  hM02HigIsRat[icase]->Divide(hM02HigIs[0]);
590 
591  hM02HigIsRat[icase]->SetYTitle("Ratio");
592  hM02HigIsRat[icase]->SetMaximum(2.5);
593  hM02HigIsRat[icase]->SetMinimum(0.5);
594 
595  if(icase == 0) hM02HigIsRat[icase]->Draw("H");
596  else hM02HigIsRat[icase]->Draw("H same");
597  }
598 
599  lABCDRat->Draw();
600 
602  cABCDRatio->cd(3);
603  gPad->SetGridy();
604 
605  for(Int_t icase = 1; icase < nMCCases; icase++)
606  //for(Int_t icase = 0; icase < nMCCases-1; icase++)
607  {
608  if(!hM02LowNo[icase]) continue;
609 
610  TString cloneName = Form("M02LowNoCasesRat%s_%s_%s",mcCases[icase].Data(),mcCasesPath[icase].Data(),cuts.Data());
611  // Careful, in case of data, take the data name and not the MC case name
612  if ( titleProd.Contains("LHC") ) cloneName = Form("M02LowNoCasesRat%s_%s_%s",titleProd.Data(),"",cuts.Data());
613 
614  hM02LowNoRat[icase] = (TH1D*) hM02LowNo[icase]->Clone(cloneName);
615  hM02LowNoRat[icase]->Divide(hM02LowNo[0]);
616  //hM02LowNoRat[icase]->Divide(hM02LowNo[nMCCases-1]);
617 
618  hM02LowNoRat[icase]->SetYTitle("Ratio");
619  hM02LowNoRat[icase]->SetMaximum(2.5);
620  hM02LowNoRat[icase]->SetMinimum(0.5);
621 
622  if(icase == 0) hM02LowNoRat[icase]->Draw("H");
623  else hM02LowNoRat[icase]->Draw("H same");
624  }
625 
626  lABCDRat->Draw();
627 
629  cABCDRatio->cd(4);
630  gPad->SetGridy();
631 
632  for(Int_t icase = 1; icase < nMCCases; icase++)
633  //for(Int_t icase = 0; icase < nMCCases-1; icase++)
634  {
635  if(!hM02HigNo[icase]) continue;
636 
637  TString cloneName = Form("M02HighNoCasesRat%s_%s_%s",mcCases[icase].Data(),mcCasesPath[icase].Data(),cuts.Data());
638  // Careful, in case of data, take the data name and not the MC case name
639  if ( titleProd.Contains("LHC") ) cloneName = Form("M02HighNoCasesRat%s_%s_%s",titleProd.Data(),"",cuts.Data());
640 
641  hM02HigNoRat[icase] = (TH1D*) hM02HigNo[icase]->Clone(cloneName);
642  //hM02HigNoRat[icase]->Divide(hM02HigNo[nMCCases-1]);
643  hM02HigNoRat[icase]->Divide(hM02HigNo[0]);
644 
645  hM02HigNoRat[icase]->SetYTitle("Ratio");
646  hM02HigNoRat[icase]->SetMaximum(2.5);
647  hM02HigNoRat[icase]->SetMinimum(0.5);
648 
649  if(icase == 0) hM02HigNoRat[icase]->Draw("H");
650  else hM02HigNoRat[icase]->Draw("H same");
651  }
652 
653  lABCDRat->Draw();
654 
655 
656  fileName = Form("%s/figures/ABCD_%s_pT_CasesRatio_%s",opt.Data(),cuts.Data(),prodTitleOutput.Data());
657  fileName+=format;
658  cABCDRatio->Print(fileName);
659  }
660 
662  TCanvas * cABCD2Ratio = new TCanvas(Form("cABCD2Ratio_%s_%s" ,cuts.Data(),prodTitleOutput.Data()),
663  Form("ABCD Double Ratio, %s, %s",cuts.Data(),prodTitleOutput.Data()),
664  1000,1000);
665 
666  TLegend *lABCD2Rat = new TLegend(0.15,0.6,0.5,0.89);
667  lABCD2Rat->SetFillColor(0);
668  lABCD2Rat->SetFillStyle(0);
669  lABCD2Rat->SetLineColor(0);
670  lABCD2Rat->SetBorderSize(0);
671  lABCD2Rat->SetTextSize(0.05);
672 
673  gPad->SetGridy();
674 
675  for(Int_t icase = 0; icase < nMCCases; icase++)
676  {
677  if(!hM02LowIs[icase]) continue;
678 
679  TString cloneName = Form("DoubleRat%s_%s_%s",mcCases[icase].Data(),mcCasesPath[icase].Data(),cuts.Data());
680  // Careful, in case of data, take the data name and not the MC case name
681  if ( titleProd.Contains("LHC") ) cloneName = Form("DoubleRat%s_%s",titleProd.Data(),cuts.Data());
682 
683  //printf(">>>>> CloneName %s for %s\n",cloneName.Data(),titleProd.Data());
684 
685  hDoubleRatio[icase] = (TH1D*) hM02LowIs[icase]->Clone(cloneName);
686  hDoubleRatio[icase]->Divide (hM02LowNo[icase]);
687  hDoubleRatio[icase]->Divide (hM02HigIs[icase]);
688  hDoubleRatio[icase]->Multiply(hM02HigNo[icase]);
689 
690  hDoubleRatio[icase]->SetYTitle("(A*D)/(C*B)");
691  hDoubleRatio[icase]->SetTitle(Form("#Sigma #it{p}_{T} > 1 GeV/c - AC: %2.2f < #sigma_{long} < %2.2f - BD: %2.2f < #sigma_{long} < %2.2f",
692  m02LowCutMin,m02LowCutMax,m02HigCutMin,m02HigCutMax));
693 
694  if(usePi0) hDoubleRatio[icase]->SetTitle(Form("#Sigma #it{p}_{T} > 1 GeV/c - AC: %2.2f < #sigma_{long} < %2.2f - BD:Merged #pi^{0} band",
695  m02LowCutMin,m02LowCutMax));
696 
697  hDoubleRatio[icase]->SetMaximum(1.4);
698  hDoubleRatio[icase]->SetMinimum(0.7);
699 
700  if(titleProd.Contains("LHC"))
701  {
702  hDoubleRatio[icase]->SetMaximum(8);
703  hDoubleRatio[icase]->SetMinimum(0.8);
704  }
705 
706  if(icase == 0) hDoubleRatio[icase]->Draw("H");
707  else hDoubleRatio[icase]->Draw("H same");
708 
709  lABCD2Rat->AddEntry(hDoubleRatio[icase],Form("%s",mcLeg[icase].Data()),"L");
710  }
711 
712  lABCD2Rat->Draw();
713 
714 
715  fileName = Form("%s/figures/ABCD_%s_DoubleRatio_%s",opt.Data(),cuts.Data(),prodTitleOutput.Data());
716  fileName+=format;
717  cABCD2Ratio->Print(fileName);
718 
719  if(nMCCases > 1)
720  {
722  TCanvas * cABCD2RatioCases = new TCanvas(Form("cABCD2RatioCases_%s_%s" ,cuts.Data(),prodTitleOutput.Data()),
723  Form("ABCD Double Ratio Cases, %s, %s",cuts.Data(),prodTitleOutput.Data()),
724  1000,1000);
725 
726  TLegend *lABCD2RatCa = new TLegend(0.15,0.6,0.5,0.89);
727  lABCD2RatCa->SetFillColor(0);
728  lABCD2RatCa->SetFillStyle(0);
729  lABCD2RatCa->SetLineColor(0);
730  lABCD2RatCa->SetBorderSize(0);
731  lABCD2RatCa->SetTextSize(0.05);
732  //lABCD2RatCa->SetHeader(Form("Den.: %s",mcLeg[nMCCases-1].Data()));
733  lABCD2RatCa->SetHeader(Form("Den.: %s",mcLeg[0].Data()));
734 
735  gPad->SetGridy();
736 
737  //for(Int_t icase = 0; icase < nMCCases-1; icase++)
738  for(Int_t icase = 1; icase < nMCCases; icase++)
739  {
740  if(!hM02LowIs[icase]) continue;
741 
742  TString cloneName = Form("DoubleRatCases%s_%s_%s",mcCases[icase].Data(),mcCasesPath[icase].Data(),cuts.Data());
743  // Careful, in case of data, take the data name and not the MC case name
744  if ( titleProd.Contains("LHC") ) cloneName = Form("DoubleRatCases%s_%s_%s",titleProd.Data(),"",cuts.Data());
745 
746  hDoubleRatioCases[icase] = (TH1D*) hDoubleRatio[icase]->Clone(cloneName);
747  //hDoubleRatioCases[icase]->Divide(hDoubleRatio[nMCCases-1]);
748  hDoubleRatioCases[icase]->Divide(hDoubleRatio[0]);
749 
750  hDoubleRatioCases[icase]->SetTitle(Form("#Sigma #it{p}_{T} > 1 GeV/c - AC: %2.2f < #sigma_{long} < %2.2f - BD: %2.2f < #sigma_{long} < %2.2f",
751  m02LowCutMin,m02LowCutMax,m02HigCutMin,m02HigCutMax));
752 
753  if(usePi0) hDoubleRatioCases[icase]->SetTitle(Form("#Sigma #it{p}_{T} > 1 GeV/c - AC: %2.2f < #sigma_{long} < %2.2f - BD: Merged #pi^{0} band",
754  m02LowCutMin,m02LowCutMax));
755  hDoubleRatioCases[icase]->SetYTitle("Ratio");
756  hDoubleRatioCases[icase]->SetMaximum(1.4);
757  hDoubleRatioCases[icase]->SetMinimum(0.7);
758 
759  if(icase == 0) hDoubleRatioCases[icase]->Draw("");
760  else hDoubleRatioCases[icase]->Draw("same");
761 
762  lABCD2RatCa->AddEntry(hDoubleRatioCases[icase],Form("Num.: %s",mcLeg[icase].Data()),"L");
763  }
764 
765  lABCD2RatCa->Draw();
766 
767 
768  fileName = Form("%s/figures/ABCD_%s_DoubleRatio_CaseRatio_%s",opt.Data(),cuts.Data(),prodTitleOutput.Data());
769  fileName+=format;
770  cABCD2RatioCases->Print(fileName);
771  }
772 
773  //-----------------------------------------------
774  // Save the ratios in external file
775  //-----------------------------------------------
776  fileName = Form("%s/figures/ABCD_%s_DoubleRatios_%s.root",opt.Data(),cuts.Data(),prodTitleOutput.Data());
777 
778  TFile * fout = new TFile(fileName,"recreate");
779 
780  for(Int_t icase = 0; icase < nMCCases; icase++)
781  {
782  //printf("icase %d\n",icase);
783  if(hM02LowIs[icase])hM02LowIs [icase]->Write();
784  if(hM02LowNo[icase])hM02LowNo [icase]->Write();
785  if(hM02HigIs[icase])hM02HigIs [icase]->Write();
786  if(hM02HigNo[icase])hM02HigNo [icase]->Write();
787 
788  if(hM02LowIsRat[icase])hM02LowIsRat[icase]->Write();
789  if(hM02LowNoRat[icase])hM02LowNoRat[icase]->Write();
790  if(hM02HigIsRat[icase])hM02HigIsRat[icase]->Write();
791  if(hM02HigNoRat[icase])hM02HigNoRat[icase]->Write();
792 
793  if(hDoubleRatio [icase])hDoubleRatio [icase]->Write();
794  if(hDoubleRatioCases[icase])hDoubleRatioCases[icase]->Write();
795  }
796 
797  fout->Close();
798 }
799 
800 //------------------------------------------------------------------------------
819 //------------------------------------------------------------------------------
820 void MakePurity
821 (
822  TString dataName, TString simuName
823  , Int_t nCases, TString * mcCasesPath, TString * mcCases, TString * mcLeg
824  , TString gj = ""//"_GJx0.50","_GJx1.00","_GJx2.00"
825  , Bool_t usePi0 = kTRUE
826  , Float_t m02LowCutMin = 0.10
827  , Float_t m02LowCutMax = 0.27
828  , Float_t m02HigCutMin = 0.40
829  , Float_t m02HigCutMax = 3.00
830  , Bool_t debug = kFALSE
831  )
832 {
833  // Init histogram arrays and other things
834  //
835  const Int_t nMCCases = nCases;
836 
837  TH1D* hDoubleRatioData = 0;
838  TH1D* hDoubleRatioSimu [nMCCases];
839  TH1D* hPurityCases [nMCCases];
840  TH1D* hPurityCasesRatio[nMCCases];
841 
842  for(Int_t icase = 0; icase < nMCCases; icase++)
843  {
844  hDoubleRatioSimu [icase]=0;
845  hPurityCases [icase]=0;
846  hPurityCasesRatio[icase]=0;
847  }
848 
849  TString cuts = "";
850  TString merged = "";
851  if ( usePi0 )
852  {
853  merged = "_MergedPi0";
854  cuts = Form("Low_%2.2f-%2.2f" ,m02LowCutMin,m02LowCutMax);
855  }
856  else
857  cuts = Form("Low_%2.2f-%2.2f_High_%2.2f-%2.2f",m02LowCutMin,m02LowCutMax,m02HigCutMin,m02HigCutMax);
858 
859  // Histogram limits
860  Float_t emin = 8;
861  Float_t emax = 29;
862  if ( simuName.Contains("High") )
863  {
864  emin = 16;
865  emax = 45;
866  }
867 
868  // ----------
869  // Open files
870  // ----------
871  TFile* fileData = TFile::Open(Form("%s/figures/ABCD_%s_DoubleRatios_%s%s.root" ,opt.Data(),cuts.Data(),dataName.Data(),merged.Data()));
872  TFile* fileSimu = TFile::Open(Form("%s/figures/ABCD_%s_DoubleRatios_%s%s%s.root",opt.Data(),cuts.Data(),simuName.Data(),gj.Data(),merged.Data()));
873 
874  if ( debug )
875  {
876  printf("Get files (%p, %p) and histograms \n",fileData,fileSimu);
877  printf("\t %s/figures/ABCD_%s_DoubleRatios_%s%s.root\n" ,opt.Data(),cuts.Data(),dataName.Data(),merged.Data());
878  printf("\t %s/figures/ABCD_%s_DoubleRatios_%s%s%s.root\n",opt.Data(),cuts.Data(),simuName.Data(),gj.Data(),merged.Data());
879  }
880 
881  if ( !fileData || !fileSimu ) return;
882 
883  // Data
884  hDoubleRatioData = (TH1D*) fileData->Get(Form("DoubleRat%s_%s",dataName.Data(),cuts.Data()));
885 
886  if ( debug )
887  {
888  printf("\t data histo %p\n",hDoubleRatioData);
889  printf("\t \t DoubleRat%s_\n",dataName.Data());
890  }
891 
892  if ( !hDoubleRatioData ) return;
893 
894  // Simu
895  for(Int_t icase = 0; icase < nMCCases; icase++)
896  {
897  hDoubleRatioSimu[icase] = (TH1D*) fileSimu->Get(Form("DoubleRat%s_%s_%s",mcCases[icase].Data(),mcCasesPath[icase].Data(),cuts.Data()));
898  if ( debug )
899  {
900  printf("\t simu case %d histo %p\n",icase, hDoubleRatioSimu[icase]);
901  printf("\t \t DoubleRat%s_%s\n",mcCases[icase].Data(),mcCasesPath[icase].Data());
902  }
903 
904  if ( !hDoubleRatioSimu[icase] ) return;
905 
906  hPurityCases[icase] = (TH1D*) hDoubleRatioSimu[icase]->Clone(Form("Purity_Data_%s_MC_%s_%s_%s%s",
907  dataName.Data(),mcCases[icase].Data(),mcCasesPath[icase].Data(),cuts.Data(),gj.Data()));
908  hPurityCases[icase]->Divide(hPurityCases[icase],hDoubleRatioData,1,1,"");
909 
910  for(Int_t ibin = 1; ibin < hPurityCases[icase]->GetNbinsX(); ibin++)
911  {
912  // Double_t content = 1.-hDoubleRatioSimu[icase]->GetBinContent(ibin)/hDoubleRatioData->GetBinContent(ibin);
913  // Double_t error = content* TMath::Sqrt( hDoubleRatioSimu[icase]->GetBinError (ibin)*hDoubleRatioSimu[icase]->GetBinError (ibin)/
914  // (hDoubleRatioSimu[icase]->GetBinContent(ibin)*hDoubleRatioSimu[icase]->GetBinContent(ibin)) +
915  // hDoubleRatioData ->GetBinError (ibin)*hDoubleRatioData ->GetBinError (ibin)/
916  // (hDoubleRatioData ->GetBinContent(ibin)*hDoubleRatioData ->GetBinContent(ibin)));
917 
918  Double_t content2 = 1.-hPurityCases[icase]->GetBinContent(ibin);
919  Double_t error2 = hPurityCases[icase]->GetBinError(ibin);
920  // printf("ibin %d purity = 1-%2.2f(+-%2.2f)*%2.2f(+-%2.2f) = %2.2f +- %2.2f; p2 = %2.2f +- %2.2f \n",
921  // ibin,hDoubleRatioSimu[icase]->GetBinContent(ibin),hDoubleRatioSimu[icase]->GetBinError(ibin),
922  // hDoubleRatioData ->GetBinContent(ibin),hDoubleRatioData ->GetBinError(ibin),
923  // content,error,content2,error2);
924 
925  hPurityCases[icase]->SetBinContent(ibin,content2);
926  hPurityCases[icase]->SetBinError (ibin,error2 );
927  }
928 
929  hPurityCases[icase]->SetMarkerStyle(marker[1]);
930  hPurityCases[icase]->SetMarkerColor(color[icase]);
931  hPurityCases[icase]->SetLineColor (color[icase]);
932  hPurityCases[icase]->SetAxisRange(emin,emax,"X");
933  hPurityCases[icase]->SetYTitle("Purity");
934  } // simu case loop
935 
936  if ( debug )
937  printf("** Make purity plots ** \n");
938 
939  // Ratio to last case
940  for(Int_t icase = 1; icase < nMCCases; icase++)
941  //for(Int_t icase = 0; icase < nMCCases-1; icase++)
942  {
943  hPurityCasesRatio[icase] = (TH1D*) hPurityCases[icase]->Clone(Form("Purity_%s_%s_%s%s_Ratio",
944  mcCases[icase].Data(),mcCasesPath[icase].Data(),cuts.Data(),gj.Data()));
945  hPurityCasesRatio[icase]->Divide(hPurityCases[0] );
946  //hPurityCasesRatio[icase]->Divide(hPurityCases[nMCCases-1] );
947  //hPurityCases[icase]->SetYTitle(Form("Purity X / Purity %s",mcLeg[nMCCases-1].Data()));
948  hPurityCases[icase]->SetYTitle(Form("Purity X / Purity %s",mcLeg[0].Data()));
949  }
950 
951  gStyle->SetPadRightMargin(0.01);
952  gStyle->SetPadLeftMargin(0.12);
953  gStyle->SetPadTopMargin(0.08);
954 
955  gStyle->SetTitleFontSize(0.06);
956  // gStyle->SetStatFontSize(0.06);
957  gStyle->SetOptStat(0);
958  gStyle->SetOptTitle(1);
959 
961  TCanvas * cPurity = new TCanvas(Form("cPurity_%s_Data_%s_MC_%s%s%s",cuts.Data(),dataName.Data(),simuName.Data(),gj.Data(),merged.Data()),
962  Form("Purity: %s_Data_%s_MC_%s%s%s",cuts.Data(),dataName.Data(),simuName.Data(),gj.Data(),merged.Data())
963  ,1000,1000);
964 
965  TLegend *lPu = new TLegend(0.3,0.15,0.98,0.4);
966  lPu->SetFillColor(0);
967  lPu->SetFillStyle(0);
968  lPu->SetLineColor(0);
969  lPu->SetBorderSize(0);
970  lPu->SetTextSize(0.05);
971 
972  gPad->SetGridy();
973 
974  for(Int_t icase = 0; icase < nMCCases; icase++)
975  {
976  hPurityCases[icase]->SetTitle(Form("#Sigma #it{p}_{T} > 1 GeV/c - AC: %2.2f < #sigma_{long} < %2.2f - BD: %2.2f < #sigma_{long} < %2.2f ",
977  m02LowCutMin,m02LowCutMax,m02HigCutMin,m02HigCutMax));
978 
979  if(usePi0) hPurityCases[icase]->SetTitle(Form("#Sigma #it{p}_{T} > 1 GeV/c - AC: %2.2f < #sigma_{long} < %2.2f - BD: Merged #pi^{0} band",
980  m02LowCutMin,m02LowCutMax));
981 
982  hPurityCases[icase]->SetMaximum(0.85);
983  hPurityCases[icase]->SetMinimum(0.0);
984 
985  if(icase == 0) hPurityCases[icase]->Draw("H");
986  else hPurityCases[icase]->Draw("H same");
987 
988  lPu->AddEntry(hPurityCases[icase],Form("%s",mcLeg[icase].Data()),"L");
989  }
990 
991  lPu->Draw();
992 
993 
994  TString fileName = Form("%s/figures/Purity_%s_Data_%s_MC_%s%s%s",opt.Data(),cuts.Data(),dataName.Data(),simuName.Data(),gj.Data(),merged.Data());
995 
996  fileName+=format;
997  cPurity->Print(fileName);
998 
999  if(nMCCases > 1)
1000  {
1002  TCanvas * cPurityRat = new TCanvas(Form("cPurityRat_%s_Data_%s_MC_%s%s%s",cuts.Data(),dataName.Data(),simuName.Data(),gj.Data(),merged.Data()),
1003  Form("Purity Ratio Cases: %s_Data_%s_MC_%s%s%s",cuts.Data(),dataName.Data(),simuName.Data(),gj.Data(),merged.Data()),
1004  1000,1000);
1005 
1006  TLegend *lPuCa = new TLegend(0.3,0.7,0.98,0.89);
1007  lPuCa->SetFillColor(0);
1008  lPuCa->SetFillStyle(0);
1009  lPuCa->SetLineColor(0);
1010  lPuCa->SetBorderSize(0);
1011  lPuCa->SetTextSize(0.04);
1012  //lPuCa->SetHeader(Form("Den.: %s",mcLeg[nMCCases-1].Data()));
1013  lPuCa->SetHeader(Form("Den.: %s",mcLeg[0].Data()));
1014 
1015  gPad->SetGridy();
1016 
1017  for(Int_t icase = 1; icase < nMCCases; icase++)
1018  //for(Int_t icase = 0; icase < nMCCases-1; icase++)
1019  {
1020  hPurityCasesRatio[icase]->SetTitle(Form("#Sigma #it{p}_{T} > 1 GeV/c - AC: %2.2f < #sigma_{long} < %2.2f - BD: %2.2f < #sigma_{long} < %2.2f",
1021  m02LowCutMin,m02LowCutMax,m02HigCutMin,m02HigCutMax));
1022  if(usePi0) hPurityCasesRatio[icase]->SetTitle(Form("#Sigma #it{p}_{T} > 1 GeV/c - AC: %2.2f < #sigma_{long} < %2.2f - BD: Merged #pi^{0} band",
1023  m02LowCutMin,m02LowCutMax));
1024 
1025  hPurityCasesRatio[icase]->SetYTitle("Ratio");
1026 
1027  hPurityCasesRatio[icase]->SetMaximum(1.5);
1028  hPurityCasesRatio[icase]->SetMinimum(0.7);
1029 
1030  if(icase == 0) hPurityCasesRatio[icase]->Draw("");
1031  else hPurityCasesRatio[icase]->Draw("same");
1032 
1033  lPuCa->AddEntry(hPurityCasesRatio[icase],Form("Num.: %s",mcLeg[icase].Data()),"L");
1034  }
1035 
1036  lPuCa->Draw();
1037 
1038 
1039  fileName = Form("%s/figures/Purity_%s_CaseRatio_Data_%s_MC_%s%s%s",opt.Data(),cuts.Data(),dataName.Data(),simuName.Data(),gj.Data(),merged.Data());
1040  fileName+=format;
1041  cPurityRat->Print(fileName);
1042  }
1043 
1044 }
1045 
1046 //------------------------------------------------------------------------------
1053 //------------------------------------------------------------------------------
1054 void MakePurityCalculationAndComparisons
1055 (
1056  Bool_t doDoRat = kTRUE,
1057  Bool_t doPurity = kTRUE,
1058  Bool_t debug = kFALSE
1059 )
1060 {
1061  // Declare the input histograms, output of the analysis
1062  // First the data analysis, after all the MC analysis from different jet-jet productions
1063  // Do not merge the Gamma-jet to the jet-jet production, only internally.
1064  //
1065  Int_t nData = 3;
1066  TString titleDataMC[] = {"LHC11cd_EMC7","pp_7TeV_Dec_Low","pp_7TeV_Dec_High"};
1067  TString gjProd = "pp_7TeV_GJ";
1068 
1069  // In case gamma-jet is added internally, define the extra scaling factors for comparison
1070  // in a later stage.
1071  Int_t nfactors = 3;
1072  Float_t gjFactor [] = { 1, 2, 0.5 }; // Used MakeDoubleRatio
1073  TString gjFactorString [] = {"","_GJx1.00","_GJx2.00","_GJx0.50"}; // Used in MakePurity
1074 
1075  // Declare here the MC analysis with different analysis cuts/emulation/settings
1076  // and their corresponding title for the legends and output files
1077  //
1078  Int_t nCases = 2;
1079 
1080  TString mcCasesPath[] =
1081  {
1082  "mimic0_Scaled2_v3"
1083  , "mimic10c_EcellCut_Scaled2_v3"
1084  };
1085 
1086  TString mcCases[] =
1087  {
1088  "Default"
1089  , "TCardEmulation"
1090  };
1091 
1092  TString mcLeg[] =
1093  {
1094  "MC default"
1095  , "MC mimic"
1096  };
1097 
1098  // Declare here the shoer shape cuts for the ABCD regions
1099  //
1100  Int_t nShShCutsMaxLow = 3;
1101  Int_t nShShCutsMinHigh = 3;
1102  Float_t m02LowCutMin = 0.10;
1103  Float_t m02LowCutMax [] = {0.27, 0.30, 0.40};
1104  Float_t m02HigCutMin [] = {0.30, 0.40, 0.50};
1105  Float_t m02HigCutMax = 2.00;
1106 
1107  if ( debug )
1108  printf("N data %d, cases %d nShShCuts MaxLow %d MinHigh %d\n",
1109  nData, nCases, nShShCutsMaxLow, nShShCutsMinHigh);
1110 
1111  Int_t rebin = 0; // Use defined E binning
1112  Int_t npi0 = 2;
1113  Int_t ngj = 2;
1114 
1115 // nfactors = 1;
1116 // nShShCutsMaxLow = 1;
1117 // nShShCutsMinHigh = 1;
1118  //nCases = 1;
1119 
1120  // Calculate and plot double ratios
1121  //
1122  if ( doDoRat )
1123  {
1124  if ( debug )
1125  printf("=== Make double ratios === \n");
1126 
1127  for(Int_t igj = 0; igj < ngj; igj++)
1128  {
1129  for(Int_t ifactor = 0; ifactor < nfactors; ifactor++)
1130  {
1131  if ( igj == 0 && ifactor > 0 ) continue;
1132 
1133  for(Int_t idata = 0; idata < nData; idata++)
1134  {
1135  if( idata == 0 && igj > 0 ) continue;
1136 
1137  for(Int_t ish = 0; ish < nShShCutsMaxLow; ish++)
1138  {
1139  for(Int_t jsh = 0; jsh < nShShCutsMinHigh; jsh++)
1140  {
1141  if ( m02LowCutMax[ish] > m02HigCutMin[jsh] ) continue;
1142 
1143  for(Int_t ipi0 = 0; ipi0 < npi0; ipi0++)
1144  {
1145  if ( ipi0 && (ish > 1 || jsh > 0) )
1146  continue;
1147 
1148  if ( debug )
1149  printf("\t idata %d %s, gj %s, usepi0 %d, addGJ %d, GJ factor %d-%2.1f, ish %d %2.2f, jsh %d %2.2f\n",
1150  idata, titleDataMC[idata].Data(), gjProd.Data(),
1151  ipi0,igj,ifactor, gjFactor[ifactor],
1152  ish, m02LowCutMax[ish],
1153  jsh, m02HigCutMin[jsh]);
1154 
1155  Int_t nCases2 = nCases;
1156 
1157  // For data, there is only one case in this example
1158  if ( titleDataMC[idata].Contains("LHC") )
1159  nCases2 = 1;
1160 
1161  MakeDoubleRatios
1162  (titleDataMC[idata], gjProd,
1163  nCases2, mcCasesPath, mcCases, mcLeg,
1164  igj, gjFactor[ifactor], ipi0, rebin,
1165  m02LowCutMin , m02LowCutMax[ish],
1166  m02HigCutMin[jsh], m02HigCutMax,
1167  debug );
1168  } // use Pi0
1169  } // jsh
1170  } // ish
1171  } // idata
1172  } // GJ scale factor
1173  } // add GJ
1174  } // make double ratios
1175 
1176 
1177  // Calculate and plot the purity
1178  //
1179  //nfactors = 0;
1180  if ( doPurity )
1181  {
1182  if ( debug )
1183  printf("=== Make purity === \n");
1184 
1185  // Data and JJ low
1186  for(Int_t ifactor = 0; ifactor < nfactors+1; ifactor++)
1187  {
1188  for(Int_t ish = 0; ish < nShShCutsMaxLow; ish++)
1189  {
1190  for(Int_t jsh = 0; jsh < nShShCutsMinHigh; jsh++)
1191  {
1192  if ( m02LowCutMax[ish] > m02HigCutMin[jsh] ) continue;
1193 
1194  for(Int_t ipi0 = 0; ipi0 < npi0; ipi0++)
1195  {
1196  if ( ipi0 && (ish > 1 || jsh > 0) )
1197  continue;
1198 
1199  if ( debug )
1200  printf("\t usepi0 %d, GJ factor %d-%2.1f, ish %d %2.2f, jsh %d %2.2f\n",
1201  ipi0,ifactor, gjFactor[ifactor],
1202  ish, m02LowCutMax[ish],
1203  jsh, m02HigCutMin[jsh]);
1204 
1205  // JJ low
1206  MakePurity(titleDataMC[0],titleDataMC[1],
1207  nCases, mcCasesPath, mcCases, mcLeg,
1208  gjFactorString[ifactor], ipi0,
1209  m02LowCutMin , m02LowCutMax[ish],
1210  m02HigCutMin[jsh], m02HigCutMax, debug);
1211 
1212  // JJ High
1213  MakePurity(titleDataMC[0],titleDataMC[2],
1214  nCases, mcCasesPath, mcCases, mcLeg,
1215  gjFactorString[ifactor], ipi0,
1216  m02LowCutMin , m02LowCutMax[ish],
1217  m02HigCutMin[jsh], m02HigCutMax, debug);
1218  } // ipi0
1219  } // jsh
1220  }// ish
1221  } // GJ factor
1222  } // Make purity
1223 }
1224 
1225 
Double_t
double Double_t
Definition: External.C:58
marker
Int_t marker[]
Definition: MakePurityCalculationAndComparisons.C:37
ScaleBinBySize
static void ScaleBinBySize(TH1D *h)
Scale histogram bins by its size.
Definition: PlotUtils.C:285
MakePurityCalculationAndComparisons
void MakePurityCalculationAndComparisons(Bool_t doDoRat=kTRUE, Bool_t doPurity=kTRUE, Bool_t debug=kFALSE)
Definition: MakePurityCalculationAndComparisons.C:1055
fileName
TString fileName
Definition: ProjectCombinHFAndFit.C:25
Rebin
TH1 * Rebin(TH1 *h, Int_t rebin, Bool_t cutEdges) const
Definition: ExtractGSEs.C:167
opt
TString opt
Definition: MakePurityCalculationAndComparisons.C:41
dataFileDir
TString dataFileDir
Definition: MakePurityCalculationAndComparisons.C:42
PlotUtils.C
Plotting utilities.
format
TString format
Definition: MakePurityCalculationAndComparisons.C:38
MakePurity
void MakePurity(TString dataName, TString simuName, Int_t nCases, TString *mcCasesPath, TString *mcCases, TString *mcLeg, TString gj="", Bool_t usePi0=kTRUE, Float_t m02LowCutMin=0.10, Float_t m02LowCutMax=0.27, Float_t m02HigCutMin=0.40, Float_t m02HigCutMax=3.00, Bool_t debug=kFALSE)
Definition: MakePurityCalculationAndComparisons.C:821
Int_t
int Int_t
Definition: External.C:63
Float_t
float Float_t
Definition: External.C:68
TH2D
Definition: External.C:228
TH1D
Definition: External.C:212
Data
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
Definition: charmCutsOptimization.C:371
MakeDoubleRatios
void MakeDoubleRatios(TString titleProd, TString gjProd, Int_t nCases, TString *mcCasesPath, TString *mcCases, TString *mcLeg, Bool_t addGJ=1, Float_t gjFactor=1, Bool_t usePi0=0, Int_t rebin=0, Float_t m02LowCutMin=0.10, Float_t m02LowCutMax=0.27, Float_t m02HigCutMin=0.4, Float_t m02HigCutMax=3.00, Int_t debug=kFALSE)
Definition: MakePurityCalculationAndComparisons.C:67
rebin
Int_t rebin
Definition: charmCutsOptimization.C:32
color
Int_t color[]
Definition: MakePurityCalculationAndComparisons.C:35
binE
Double_t binE[]
Definition: CalculateShowerShapeFractionPerNCellCut.C:53
Bool_t
bool Bool_t
Definition: External.C:53
fout
TFile * fout
input train file
Definition: DrawAnaCaloTrackQA.C:90
lineStyle
Int_t lineStyle[]
Definition: MakePurityCalculationAndComparisons.C:36