AliPhysics  eb0e5d9 (eb0e5d9)
CalculateShowerShapeFractionPerNCellCut.C
Go to the documentation of this file.
1 
16 #if !defined(__CINT__) || defined(__MAKECINT__)
17 
18 #include <TFile.h>
19 #include <TDirectoryFile.h>
20 #include <TList.h>
21 #include <TString.h>
22 #include <TH1D.h>
23 #include <TH2F.h>
24 #include <TH3F.h>
25 #include <TCanvas.h>
26 #include <TPad.h>
27 #include <TLegend.h>
28 #include <TObject.h>
29 #include <TAxis.h>
30 #include <TGaxis.h>
31 #include <TLine.h>
32 #include <TF1.h>
33 #include <TMath.h>
34 #include <TGraphErrors.h>
35 #include <TROOT.h>
36 #include <TStyle.h>
37 #include "PlotUtils.C"
38 
39 #endif
40 
41 //--------------------Global variables------------------------------------------
42 
43 // Cluster types
44 TString anaName [] = {"NoIso","Iso","All"};
45 TString anaTitle[] = {"NOT Isolated","Isolated","Inclusive"};
46 
47 // Cut cases, 0-no cut, 1-n cell > 4, 2-ncell <=4
48 const Int_t nHisto = 3;
49 TString histoLeg[] = {"#it{n}^{w}_{cell}>1", "#it{n}^{w}_{cell} > 4","#it{n}^{w}_{cell} #leq 4"};
50 
51 // Energy bins
52 const Int_t nEBins = 18;
53 Double_t binE [] = { 2, 3, 4, 5, 6, 7, 8,10,12,14,16,18,20,22, 25, 30, 35, 40,50};
54 Double_t binEErr[] = {0.5,0.5,0.5,0.5,0.5,0.5,0.5, 1, 1, 1, 1, 1, 1, 1,1.5, 2.5, 2.5, 2.5, 5};
55 
56 // Shower shape max cut array
57 const Int_t nShShMaxCut = 3;
58 Float_t maxShShCut[] = {0.27,0.3,0.4};
59 
60 //const Int_t nShShMaxCut = 1;
61 //Float_t maxShShCut[] = {0.3};
62 
63 Int_t shshRef = 1; // take 0.3 for some plots
64 
65 // Histogram settings
66 Int_t color [] = {1,4,2,6,8,9,kYellow-6,kCyan,kOrange+2,kViolet,kOrange-2,
67  4,2,6,8,9,kYellow-6,kCyan,kOrange+2,kViolet,kOrange-2};
68 Int_t lineStyle[] = {1,1,1,1,1,1, 1, 1, 1,
69  1, 1,2,2,2,2,2, 2, 2, 2, 2,
70  2,2,2,2,2,2,2,2,2,2,2,2};
71 Int_t marker [] = {20,20,20,21,24,24,24,24,24,24,24};
72 
73 // Plot file format
74 TString fileFormat = ".eps";
75 
76 //------------------------------------------------------------------------------
84 //------------------------------------------------------------------------------
85 void GetSMNumber(TString titleName,
86  Int_t & totalSM, Int_t & firstSM,
87  Int_t & ncol , Int_t & nrow)
88 {
89  totalSM = 10;
90  firstSM = 0;
91  ncol = 4;
92  nrow = 3;
93  GetCanvasColRowNumber(totalSM,ncol,nrow); // PlotUtils.C
94 
95  if(titleName.Contains("15") || titleName.Contains("16") ||
96  titleName.Contains("17") || titleName.Contains("18"))
97  {
98  totalSM = 12;
99  GetCanvasColRowNumber(totalSM,ncol,nrow); // PlotUtils.C
100  if(titleName.Contains("DCAL") || titleName.Contains("DMC") ||
101  titleName.Contains("DG") || titleName.Contains("DMG"))
102  {
103  totalSM = 20;
104  firstSM = 12;
105  GetCanvasColRowNumber(8,ncol,nrow); // PlotUtils.C
106  }
107  }
108 }
109 
110 //------------------------------------------------------------------------------
120 //------------------------------------------------------------------------------
121 void CalculateAndPlot
122 (TString titleName, TString filePath,
123  Int_t iana = 0,
124  TString tm = "_TMDep",
125  Bool_t plotRat = kFALSE,
126  Bool_t debug = kFALSE)
127 {
128  // ------------------------------------------------
129  // Declare histogram arrays and different settings
130  // ------------------------------------------------
131 
132  // Histogram settings
133  Int_t rebin = 2;
134 
135  Double_t xmin = 0.1;
136  Double_t xmax = 1.1;
137 
138  // Energy bin settings
139  // bins defined as global
140  Double_t lowE[nEBins];
141  Double_t higE[nEBins];
142  Double_t muE [nEBins];
143  Double_t errE[nEBins];
144  Int_t ncolE = 5;
145  Int_t nrowE = 4;
146  GetCanvasColRowNumber(nEBins,ncolE,nrowE); // PlotUtils.C
147 
148  // SM
149  //
150  Int_t totalSM = 10, firstSM = 0;
151  Int_t ncol = 4, nrow = 3;
152  GetSMNumber(titleName.Data(),totalSM,firstSM,ncol,nrow);
153  const Int_t nSM = totalSM;
154 
155  // histograms
156  TH1D* h [nHisto][nEBins][nSM];
157  TH1D* hA[nHisto][nEBins];
158 
159  // Shower shape max cut
160  Double_t frac[nEBins][nSM+1][nShShMaxCut];
161  Double_t efrac[nEBins][nSM+1][nShShMaxCut];
162 
163  // Merge good/bad SMs
164  const Int_t nGroup = 3;
165  TH1D* hGroup[nHisto][nEBins][nGroup];
166 
167  // ------------------------------------------------
168  // Open input file and fill histogram arrays
169  // ------------------------------------------------
170 
171  // Recover data file
172  TFile* file = TFile::Open(Form("%s.root",filePath.Data()));
173  if ( debug )
174  printf("Read: %s.root, %p\n",filePath.Data(),file);
175 
176  if ( !file ) return;
177 
178  //titleName.ReplaceAll("/","_");
179 
180  //
181  // Histograms from AliAnaEMCalClusterShape
182  // TH3 histogram, treat differently to AliAnaParticleIsolation
183  //
184  if ( iana == 2 )
185  {
186  TH3F* h3[nHisto];
187 
188  TString start = "Shape"; // Careful, it can change.
189  h3[0] = (TH3F*) file->Get(Form("%s%s_hSMM02NoCut_Neutral",start.Data(),tm.Data()));
190  h3[1] = (TH3F*) file->Get(Form("%s%s_hSMM02_Neutral" ,start.Data(),tm.Data()));
191  if ( debug )
192  printf("AliAnaClusterShapeStudies TH3F %p %p\n",h3[0],h3[1]);
193 
194  if(filePath.Contains("data"))
195  {
196  h3[0]->Sumw2();
197  h3[1]->Sumw2();
198  }
199 
200  if(nHisto > 2) h3[2] = (TH3F*) h3[0]->Clone("Histo2_Iso2");
201  h3[2]->Add(h3[1],-1);
202 
203  //
204  // Project
205  //
206  Double_t width = 0;
207  Int_t ebinMin =-1;
208  Int_t ebinMax =-1;
209  Int_t smbin =-1;
210  for(Int_t ihisto = 0; ihisto < nHisto; ihisto++)
211  {
212  for(Int_t iebin = 0; iebin < nEBins; iebin++)
213  {
214  ebinMin = h3[ihisto]->GetXaxis()->FindBin(binE[iebin ]);
215  ebinMax = h3[ihisto]->GetXaxis()->FindBin(binE[iebin+1])-1;
216  lowE[iebin] = h3[ihisto]->GetXaxis()->GetBinLowEdge(ebinMin);
217  width = h3[ihisto]->GetXaxis()->GetBinWidth (ebinMax);
218  higE[iebin] = h3[ihisto]->GetXaxis()->GetBinLowEdge(ebinMax)+width;
219  muE [iebin] = (lowE[iebin]+higE[iebin])/2.;
220  errE[iebin] = (higE[iebin]-lowE[iebin])/2.;
221 
222  for(Int_t ism = firstSM; ism < nSM; ism++)
223  {
224  smbin = h3[ihisto]->GetYaxis()->FindBin(ism);
225  h[ihisto][iebin][ism] =
226  (TH1D*) h3[ihisto]->ProjectionZ(Form("Histo%d_BinE%d_SM%d",ihisto,iebin,ism),ebinMin,ebinMax,smbin,smbin);
227 
228  if ( !h[ihisto][iebin][ism] ) continue;
229 
230 // printf("%s, histo %d, iebin %d, [%2.1f,%2.1f] GeV, muE %2.1f errE %2.1f,ism %d, smbin%d integral %e\n",
231 // h[ihisto][iebin][ism]->GetName(),ihisto,
232 // iebin,lowE[iebin],higE[iebin],muE[iebin],errE[iebin],
233 // ism,smbin,h[ihisto][iebin][ism]->Integral());
234 
235  h[ihisto][iebin][ism]->SetLineColor(color[ihisto]);
236 
237  h[ihisto][iebin][ism]->SetLineWidth(2);
238 
239  h[ihisto][iebin][ism]->SetLineStyle(lineStyle[ihisto]);
240 
241  h[ihisto][iebin][ism]->SetMarkerStyle(marker[ihisto]);
242 
243  h[ihisto][iebin][ism]->SetMarkerColor(color[ihisto]);
244 
245  h[ihisto][iebin][ism]->SetMarkerSize(0.5);
246 
247  if(rebin > 1)
248  h[ihisto][iebin][ism]->Rebin(rebin);
249 
250  h[ihisto][iebin][ism]->SetAxisRange(xmin,xmax,"X");
251  } // SM histo
252  } // E bin
253  } // cut histo
254  } // END Histograms from AliAnaEMCalClusterShape
255  //
256  // Histograms from AliAnaParticleIsolation
257  //
258  else
259  {
260  TH2F* h2[nHisto][nSM];
261  for(Int_t ism = firstSM; ism < nSM; ism++)
262  {
263  for(Int_t ihisto = 0; ihisto < nHisto; ihisto++)
264  {
265  if ( ihisto==0 ) h2[ihisto][ism] = (TH2F*) file->Get(Form("AnaIsolPhoton_hPtLambda0_SM%d_%s" ,ism,anaName[iana].Data()));
266  if ( ihisto==1 ) h2[ihisto][ism] = (TH2F*) file->Get(Form("AnaIsolPhoton_hPtLambda0_NCellCut_SM%d_%s",ism,anaName[iana].Data()));
267 
268  if ( ihisto==2 )
269  {
270  h2[ihisto][ism] = (TH2F*) h2[0][ism]->Clone(Form("Histo%d_Iso%d_SM%d",ihisto,iana,ism));
271  h2[ihisto][ism]->Add(h2[1][ism],-1);
272  }
273 
274  if ( ihisto < 2 )
275  {
276  if(filePath.Contains("data") ) h2[ihisto][ism]->Sumw2();
277  if(filePath.Contains("simu") && filePath.Contains("MB") ) h2[ihisto][ism]->Sumw2();
278  }
279 
280  } // cut histo
281  } // SM histo
282 
283  //
284  // Project
285  //
286  Double_t width = 0;
287  Int_t ebinMin =-1;
288  Int_t ebinMax =-1;
289  for(Int_t ism = firstSM; ism < nSM; ism++)
290  {
291  for(Int_t ihisto = 0; ihisto < nHisto; ihisto++)
292  {
293  for(Int_t iebin = 0; iebin < nEBins; iebin++)
294  {
295  ebinMin = h2[ihisto][ism]->GetXaxis()->FindBin(binE[iebin ]);
296  ebinMax = h2[ihisto][ism]->GetXaxis()->FindBin(binE[iebin+1])-1;
297  lowE[iebin] = h2[ihisto][ism]->GetXaxis()->GetBinLowEdge(ebinMin);
298  width = h2[ihisto][ism]->GetXaxis()->GetBinWidth (ebinMax);
299  higE[iebin] = h2[ihisto][ism]->GetXaxis()->GetBinLowEdge(ebinMax)+width;
300  muE [iebin] = (lowE[iebin]+higE[iebin])/2.;
301  errE[iebin] = (higE[iebin]-lowE[iebin])/2.;
302 
303  h[ihisto][iebin][ism] =
304  (TH1D*) h2[ihisto][ism]->ProjectionY(Form("Histo%d_BinE%d_SM%d",ihisto,iebin,ism),ebinMin,ebinMax);
305 
306  //printf("\t %p \n",h[ihisto][iebin][ism]);
307 
308  if ( !h[ihisto][iebin][ism] ) continue;
309 
310 // printf("%s, histo %d, iebin %d, [%2.1f,%2.1f] GeV, muE %2.1f errE %2.1f,ism %d, integral %e\n",
311 // h[ihisto][iebin][ism]->GetName(),ihisto,
312 // iebin,lowE[iebin],higE[iebin],muE[iebin],errE[iebin],
313 // ism,h[ihisto][iebin][ism]->Integral());
314 
315  h[ihisto][iebin][ism]->SetLineColor(color[ihisto]);
316 
317  h[ihisto][iebin][ism]->SetLineWidth(2);
318 
319  h[ihisto][iebin][ism]->SetLineStyle(lineStyle[ihisto]);
320 
321  h[ihisto][iebin][ism]->SetMarkerStyle(marker[ihisto]);
322 
323  h[ihisto][iebin][ism]->SetMarkerColor(color[ihisto]);
324 
325  h[ihisto][iebin][ism]->SetMarkerSize(0.5);
326 
327  if(rebin > 1)
328  h[ihisto][iebin][ism]->Rebin(rebin);
329 
330  h[ihisto][iebin][ism]->SetAxisRange(xmin,xmax,"X");
331  } // ie
332  } // ihisto
333  } // ism
334  } // END Histograms from AliAnaParticleIsolation
335 
336 
337  //
338  // Merged SM
339  //
340  for(Int_t ihisto = 0; ihisto < nHisto; ihisto++)
341  {
342  for(Int_t iebin = 0; iebin < nEBins; iebin++)
343  {
344  if ( h[ihisto][iebin][firstSM] )
345  {
346  hA [ihisto][iebin] = (TH1D*) h[ihisto][iebin][firstSM]->Clone(Form("SMAll_Cut%d_Ana%d" ,ihisto,iana));
347 
348  for(Int_t ism = firstSM+1; ism < nSM; ism++)
349  hA[ihisto][iebin]->Add( h[ihisto][iebin][ism] );
350  }
351 
352 
353  if( !titleName.Contains("DCAL") && !titleName.Contains("DMC") &&
354  !titleName.Contains("DG") && !titleName.Contains("DMG"))
355  {
356  hGroup[ihisto][iebin][0] = (TH1D*) h[ihisto][iebin][0]->Clone(Form("SM_Group0_Cut%d_Ana%d",ihisto, iana));
357  hGroup[ihisto][iebin][1] = (TH1D*) h[ihisto][iebin][1]->Clone(Form("SM_Group1_Cut%d_Ana%d",ihisto, iana));
358  hGroup[ihisto][iebin][2] = (TH1D*) h[ihisto][iebin][3]->Clone(Form("SM_Group2_Cut%d_Ana%d",ihisto, iana));
359 
360  hGroup[ihisto][iebin][0]->Add( h[ihisto][iebin][4] );
361  hGroup[ihisto][iebin][0]->Add( h[ihisto][iebin][5] );
362  hGroup[ihisto][iebin][0]->Add( h[ihisto][iebin][6] );
363  hGroup[ihisto][iebin][0]->Add( h[ihisto][iebin][8] );
364  hGroup[ihisto][iebin][0]->Add( h[ihisto][iebin][9] );
365 
366  hGroup[ihisto][iebin][1]->Add( h[ihisto][iebin][2] );
367 
368  hGroup[ihisto][iebin][2]->Add( h[ihisto][iebin][7] );
369 
370  hGroup[ihisto][iebin][0]->Scale(1./6.);
371  hGroup[ihisto][iebin][1]->Scale(1./2.);
372  hGroup[ihisto][iebin][2]->Scale(1./2.);
373  }
374 
375  } // ebins
376  } // prod bins
377 
378  // ---------------------------------------
379  // Plot
380  // ---------------------------------------
381  if ( debug )
382  printf("PLOT\n");
383 
384  TString fileName ;
385 
386  gStyle->SetPadRightMargin(0.01);
387  gStyle->SetPadLeftMargin(0.12);
388  gStyle->SetTitleFontSize(0.06);
389  // gStyle->SetStatFontSize(0.06);
390  gStyle->SetOptStat(0);
391 
392  //
393  // Plot comparison of 3 cut cases and calculate fractions per E bin
394  //
395  Int_t shshRef = 1; //Select a max Cut to be shown in the legend.
396 
397  // Plots per E bin, each par per SM
398  for(Int_t iebin = 0; iebin < nEBins; iebin++)
399  {
400  gStyle->SetOptTitle(1);
401  gStyle->SetPadTopMargin(0.1);
402 
403  TCanvas * c = new TCanvas(Form("c_ebin%d_iana%d_%s",
404  iebin,iana, titleName.Data()),
405  Form("iana %d, %2.1f < E < %2.1f, %s",
406  iana, lowE[iebin],higE[iebin], titleName.Data()),
407  ncol*2000,nrow*2000);
408  c->Divide(ncol,nrow);
409 
410  TLegend *l = new TLegend(0,0,1,1);
411  l->SetFillColor(0);
412  l->SetFillStyle(0);
413  l->SetLineColor(0);
414  l->SetBorderSize(0);
415  l->SetTextSize(0.07);
416  l->SetHeader(Form("%2.1f < #it{E} < %2.1f GeV", lowE[iebin],higE[iebin]));
417  l->AddEntry("",anaTitle[iana].Data(),"");
418 
419  for(Int_t ism = firstSM; ism < nSM; ism++)
420  {
421  c->cd(ism+1);
422 
423  //gPad->SetLogy();
424  if ( debug )
425  printf("iana %d iE %d ism %d: %p %p %p\n",
426  iana,iebin,ism,h[0][iebin][ism],h[1][iebin][ism],h[2][iebin][ism]);
427 
428  if(!h[0][iebin][ism]) continue;
429 
430  for(Int_t ihisto = 0; ihisto < nHisto; ihisto++)
431  {
432  //printf("\t prod %d %p\n",ihisto,h[ihisto][iebin][ism]);
433  if(!h[ihisto][iebin][ism]) continue;
434 
435  h[ihisto][iebin][ism]->SetTitle(Form("SM %d",ism));
436 
437  h[ihisto][iebin][ism]->SetTitleOffset(1.8,"Y");
438  //h[iebin][0]->SetLineColor(1);
439  //h[iebin][0]->SetAxisRange(0.1,2.5,"X");
440 
441  if ( ihisto == 0 ) h[ihisto][iebin][ism]->Draw("H");
442  else h[ihisto][iebin][ism]->Draw("H same");
443 
444  if(h[ihisto][iebin][ism]->GetMaximum() > h[0][iebin][ism]->GetMaximum())
445  h[0][iebin][ism]->SetMaximum(h[ihisto][iebin][ism]->GetMaximum()*1.2);
446 
447  if(ism==firstSM)
448  l->AddEntry(h[ihisto][iebin][ism],Form("%s",histoLeg[ihisto].Data()),"PL");
449  } // ihisto
450 
451  // *****
452  // Calculate the fraction of clusters with n cell > 4
453  // *****
454  Double_t integral0 = 0;
455  Double_t integral1 = 0;
456  Double_t integral0Err = 0;
457  Double_t integral1Err = 0;
458  Int_t binM02Min = h[0][iebin][ism]->FindBin(0.1);
459  for(Int_t imax = 0; imax < nShShMaxCut; imax++)
460  {
461  Int_t binM02Max = h[0][iebin][ism]->FindBin(maxShShCut[imax]);
462 
463  GetRangeIntegralAndError(h[0][iebin][ism], binM02Min, binM02Max, integral0, integral0Err );
464  GetRangeIntegralAndError(h[1][iebin][ism], binM02Min, binM02Max, integral1, integral1Err );
465 
466  Double_t ratio = 0;
467  if ( integral0 > 0 ) ratio = integral1 / integral0 ;
468  Double_t eratio = GetFractionError(integral1,integral0,integral1Err,integral0Err);
469 
470  frac [iebin][ism][imax] = ratio*100;
471  efrac[iebin][ism][imax] = eratio*100;
472 
473  if( debug )
474  printf("\t max %2.2f, frac %2.2f, err %2.2f\n",
475  maxShShCut[imax],frac [iebin][ism][imax],efrac[iebin][ism][imax]);
476  }
477 
478  TLegend *lsm = new TLegend(0.15,0.7,0.95,0.9);
479  lsm->SetFillColor(0);
480  lsm->SetFillStyle(0);
481  lsm->SetLineColor(0);
482  lsm->SetBorderSize(0);
483  lsm->SetTextSize(0.07);
484  lsm->SetHeader(Form(" For 0.10<#sigma_{long}^{2}<%2.2f",maxShShCut[shshRef]));
485  lsm->AddEntry("",Form("#it{R}_{#sigma_{long}}^{#it{n}^{w}_{cell}} = %2.1f #pm %2.1f %%",frac[iebin][ism][shshRef],efrac[iebin][ism][shshRef]),"") ;
486  lsm->Draw("same");
487  //printf("\t end\n");
488  } // param
489 
490  c->cd(ncol*nrow);
491  l->Draw();
492 
493  fileName = Form("figures/Comparison_M02_NCellCut_Ebin%d_%s_%s",
494  iebin,anaName[iana].Data(),titleName.Data());
495  fileName+=fileFormat;
496  c->Print(fileName);
497 
498  if ( !plotRat ) continue;
499 
500  TCanvas * cR = new TCanvas(Form("cR_ebin%d_iana%d_%s",
501  iebin,iana, titleName.Data()),
502  Form("Ratio iana %d, %2.1f < E < %2.1f, %s",
503  iana, lowE[iebin],higE[iebin], titleName.Data()),
504  ncol*2000,nrow*2000);
505  cR->Divide(ncol,nrow);
506 
507  for(Int_t ism = 0; ism < nSM; ism++)
508  {
509  cR->cd(ism+1);
510 
511  //gPad->SetLogy();
512 
513  //printf("iE %d ism %d\n",iebin,ism);
514  if(!h[1][iebin][ism]) continue;
515 
516  for(Int_t ihisto = 1; ihisto < nHisto; ihisto++)
517  {
518  //printf("\t prod %d %p\n",ihisto,h[ihisto][iebin][ism]);
519  if(!h[ihisto][iebin][ism]) continue;
520 
521  TH1F* hRat = (TH1F*) h[ihisto][iebin][ism]->Clone(Form("%s_Ratio",h[ihisto][iebin][ism]->GetName()));
522  hRat->Divide(h[0][iebin][ism]);
523 
524  if(ihisto==0) hRat->Draw("H");
525  else hRat->Draw("H same");
526 
527  hRat->SetYTitle("Ratio MC to Data");
528  hRat->SetMaximum(1.2);
529  hRat->SetMinimum(0);
530 
531  } // ihisto
532 
533  } // param
534 
535  cR->cd(ncol*nrow);
536  l->Draw();
537 
538  fileName = Form("figures/Comparison_Ratio_M02_NCellCut_Ebin%d_%s_%s",
539  iebin,anaName[iana].Data(),titleName.Data());
540  fileName+=fileFormat;
541 
542  cR->Print(fileName);
543 
544  } // e bin
545 
546  // Sum of all SM, each pad one energy
547  //
548  TCanvas * cA = new TCanvas(Form("cAllSM_iana%d_%s",iana, titleName.Data()),
549  Form("All SM, iana %d, %s",iana, titleName.Data()),
550  ncolE*2000,nrowE*2000);
551 
552  cA->Divide(ncolE,nrowE);
553 
554  TLegend *l2 = new TLegend(0,0.,1,1);
555  l2->SetFillColor(0);
556  l2->SetFillStyle(0);
557  l2->SetLineColor(0);
558  l2->SetBorderSize(0);
559  l2->SetTextSize(0.07);
560  l2->AddEntry("",Form("All SM, %s",anaTitle[iana].Data()),"");
561 
562  for(Int_t iebin = 0; iebin < nEBins; iebin++)
563  {
564  cA->cd(iebin+1);
565 
566  if ( debug )
567  printf("All SM, iana %d, ebin %d: %p, %p, %p\n",
568  iana,iebin,hA[0][iebin],hA[1][iebin],hA[2][iebin]);
569 
570  if(!hA[0][iebin]) continue;
571 
572  //gPad->SetGridy();
573  //gPad->SetGridx();
574  //gPad->SetLogy();
575 
576  gStyle->SetOptTitle(1);
577  gStyle->SetPadTopMargin(0.1);
578 
579  hA[0][iebin]->Draw("H");
580 
581  for(Int_t ihisto = 0; ihisto < nHisto; ihisto++)
582  {
583  if(!hA[ihisto][iebin]) continue;
584 
585  hA[ihisto][iebin]->SetTitleOffset(1.8,"Y");
586  //h[iebin][0]->SetAxisRange(0.1,2.5,"X");
587  //hA[ihisto][iebin]->SetMinimum(1e-2);
588  //printf("iebin %d, ism %d 0 %d %p\n",iebin,ism,0,h[0][iebin]);
589 
590  hA[ihisto][iebin]->SetTitle(Form("%2.1f < #it{E}^{clus} < %2.1f GeV",lowE[iebin],higE[iebin]));
591  if(higE[iebin] > 90) hA[ihisto][iebin]->SetTitle(Form("#it{E} > %2.1f GeV",lowE[iebin]));
592 
593  hA[ihisto][iebin]->Draw("H same");
594 
595  if(hA[ihisto][iebin]->GetMaximum() > hA[0][iebin]->GetMaximum())
596  hA[0][iebin]->SetMaximum(hA[ihisto][iebin]->GetMaximum()*1.2);
597 
598  if(iebin==0)
599  l2->AddEntry(hA[ihisto][iebin],Form("%s",histoLeg[ihisto].Data()),"PL");
600  }
601 
602  // *****
603  // Calculate the fraction of clusters with n cell > 4
604  // *****
605  Double_t integral0 = 0;
606  Double_t integral1 = 0;
607  Double_t integral0Err = 0;
608  Double_t integral1Err = 0;
609  Int_t binM02Min = hA[0][iebin]->FindBin(0.1);
610  for(Int_t imax = 0; imax < nShShMaxCut; imax++)
611  {
612  Int_t binM02Max = hA[0][iebin]->FindBin(maxShShCut[imax]);
613 
614  GetRangeIntegralAndError(hA[0][iebin], binM02Min, binM02Max, integral0, integral0Err );
615  GetRangeIntegralAndError(hA[1][iebin], binM02Min, binM02Max, integral1, integral1Err );
616 
617  Double_t ratio = 0;
618  if ( integral0 > 0 ) ratio = integral1 / integral0 ;
619  Double_t eratio = GetFractionError(integral1,integral0,integral1Err,integral0Err);
620 
621  // printf("A ratio %2.3f, eratio %2.3f;"
622  // " int0 %2.3e, int1 %2.3e;"
623  // " Eint0 %2.3e, Eint1 %2.3e;"
624  // " sqrt0 %2.3f sqrt1 %2.3f \n",
625  // ratio,eratio,
626  // integral0,integral1,
627  // integral0Err,integral1Err,
628  // TMath::Sqrt(integral0),TMath::Sqrt(integral1));
629 
630  frac [iebin][nSM][imax] = ratio*100;
631  efrac[iebin][nSM][imax] = eratio*100;
632 
633  if( debug )
634  printf("\t max %2.2f, frac %2.2f, err %2.2f\n",
635  maxShShCut[imax],frac [iebin][nSM][imax],efrac[iebin][nSM][imax]);
636  }
637 
638  TLegend *lE = new TLegend(0.15,0.7,0.95,0.9);
639  lE->SetFillColor(0);
640  lE->SetFillStyle(0);
641  lE->SetLineColor(0);
642  lE->SetBorderSize(0);
643  lE->SetTextSize(0.07);
644  lE->SetHeader(Form(" For 0.10<#sigma_{long}^{2}<%2.2f",maxShShCut[shshRef]));
645  lE->AddEntry("",Form("#it{R}_{#sigma_{long}}^{#it{n}^{w}_{cell}} = %2.1f #pm %2.1f %%",frac[iebin][nSM][shshRef],efrac[iebin][nSM][shshRef]),"") ;
646  lE->Draw("same");
647  }
648 
649  cA->cd(ncolE*nrowE);
650  l2->Draw("same");
651 
652  fileName = Form("figures/Comparison_M02_NCellCut_AllSM_%s_%s",
653  anaName[iana].Data(),titleName.Data());
654  fileName+=fileFormat;
655  cA->Print(fileName);
656 
657  if ( plotRat )
658  {
659  TCanvas * cAR = new TCanvas(Form("cR_AllSM_iana%d_%s",iana, titleName.Data()),
660  Form("All SM, iana %d, %s",iana, titleName.Data()),
661  ncolE*2000,nrowE*2000);
662 
663  cAR->Divide(ncolE,nrowE);
664 
665  for(Int_t iebin = 0; iebin < nEBins; iebin++)
666  {
667  cAR->cd(iebin+1);
668 
669  if(!hA[0][iebin]) continue;
670 
671  //if(quantity.Contains("ECell"))gPad->SetLogy();
672 
673  gStyle->SetOptTitle(1);
674  gStyle->SetPadTopMargin(0.1);
675 
676  hA[0][iebin]->Draw("H");
677 
678  for(Int_t ihisto = 1; ihisto < nHisto; ihisto++)
679  {
680  if(!hA[ihisto][iebin]) continue;
681 
682  TH1F* hRat = (TH1F*) hA[ihisto][iebin]->Clone(Form("%s_Ratio",hA[ihisto][iebin]->GetName()));
683  hRat->Divide(hA[0][iebin]);
684 
685  if(ihisto==1) hRat->Draw("H");
686  else hRat->Draw("H same");
687 
688  hRat->SetYTitle("Ratio MC to Data");
689  hRat->SetMaximum(1.2);
690  hRat->SetMinimum(0);
691  }
692  }
693 
694  cAR->cd(ncolE*nrowE);
695  l2->Draw("same");
696 
697  fileName = Form("figures/Comparison_Ratio_M02_NCellCut_AllSM_%s_%s",
698  anaName[iana].Data(),titleName.Data());
699 
700  fileName+=fileFormat;
701  cAR->Print(fileName);
702  }
703 
704  // Now plot each file per SM, per E bin in each pad
705  //
706  for(Int_t ism = firstSM; ism < nSM; ism++)
707  {
708  TCanvas * cSM = new TCanvas(Form("cSM%d_iana%d_%s",ism,iana,titleName.Data()),
709  Form("SM %d, iana %d, %s",ism,iana,titleName.Data()),
710  ncolE*2000,nrowE*2000);
711 
712  cSM->Divide(ncolE,nrowE);
713 
714  TLegend *l2 = new TLegend(0,0.,1,1);
715  l2->SetFillColor(0);
716  l2->SetFillStyle(0);
717  l2->SetLineColor(0);
718  l2->SetBorderSize(0);
719  l2->SetTextSize(0.07);
720  l2->AddEntry("",Form("SM %d, %s",ism,anaTitle[iana].Data()),"");
721 
722  for(Int_t iebin = 0; iebin < nEBins; iebin++)
723  {
724  cSM->cd(iebin+1);
725 
726  if(!h[0][iebin][ism]) continue;
727 
728  //gPad->SetLogy();
729 
730  gStyle->SetOptTitle(1);
731  gStyle->SetPadTopMargin(0.1);
732 
733  h[0][iebin][ism]->Draw("H");
734 
735  for(Int_t ihisto = 0; ihisto < nHisto; ihisto++)
736  {
737  if(!h[ihisto][iebin][ism]) continue;
738 
739  h[ihisto][iebin][ism]->SetTitle(Form("%2.1f < #it{E}^{clus} < %2.1f GeV",lowE[iebin],higE[iebin]));
740 
741  h[ihisto][iebin][ism]->Draw("H same");
742 
743  if(h[ihisto][iebin][ism]->GetMaximum() > h[0][iebin][ism]->GetMaximum())
744  h[0][iebin][ism]->SetMaximum(h[ihisto][iebin][ism]->GetMaximum()*1.2);
745 
746  if(iebin==0)
747  l2->AddEntry(h[ihisto][iebin][ism],Form("%s",histoLeg[ihisto].Data()),"PL");
748  }
749 
750  TLegend *lE = new TLegend(0.15,0.7,0.95,0.9);
751  lE->SetFillColor(0);
752  lE->SetFillStyle(0);
753  lE->SetLineColor(0);
754  lE->SetBorderSize(0);
755  lE->SetTextSize(0.07);
756  lE->SetHeader(Form(" For 0.10<#sigma_{long}^{2}<%2.2f",maxShShCut[shshRef]));
757  lE->AddEntry("",Form("#it{R}_{#sigma_{long}}^{#it{n}^{w}_{cell}} = %2.1f #pm %2.1f %%",frac[iebin][ism][shshRef],efrac[iebin][ism][shshRef]),"") ;
758  lE->Draw("same");
759  }
760 
761  cSM->cd(ncolE*nrowE);
762  l2->Draw("same");
763 
764  fileName = Form("figures/Comparison_M02_NCellCut_SM%d_%s_%s",
765  ism,anaName[iana].Data(),titleName.Data());
766  fileName+=fileFormat;
767  cSM->Print(fileName);
768  } // sm
769 
770  //=================================================================
771  // Put the calculated fractions in GRAPHS
772  //=================================================================
773  TGraphErrors* gFracPerEn[nSM+1] [nShShMaxCut];
774  TGraphErrors* gFracPerSM[nEBins][nShShMaxCut];
775 
776  Double_t fracEn[nEBins];
777  Double_t efracEn[nEBins];
778  Double_t fracSM[nSM+1];
779  Double_t efracSM[nSM+1];
780 
781  Double_t binSM [] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22};
782  Double_t binSMErr [] = {0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5};
783  for(Int_t imax = 0; imax < nShShMaxCut; imax++)
784  {
785  // Graphs frac vs SM per E bin
786  //
787  for(Int_t ie = 0; ie < nEBins; ie++)
788  {
789  for(Int_t ism = firstSM; ism < nSM+1; ism++)
790  {
791  fracSM[ism] = frac[ie][ism][imax];
792  efracSM[ism] = efrac[ie][ism][imax];
793  }
794 
795  gFracPerSM[ie][imax] = new TGraphErrors(nSM+1,binSM,fracSM,binSMErr,efracSM);
796  } // iE bin
797 
798  // Graphs frac vs E per SM
799  //
800  for(Int_t ism = firstSM; ism < nSM+1; ism++)
801  {
802  for(Int_t ie = 0; ie < nEBins; ie++)
803  {
804 
805  fracEn[ie] = frac[ie][ism][imax];
806  efracEn[ie] = efrac[ie][ism][imax];
807  }
808 
809  gFracPerEn[ism][imax] = new TGraphErrors(nEBins,muE,fracEn,errE,efracEn);
810 
811  // Set E ranges depending the data type
812  //
813  if( filePath.Contains("simu") )
814  {
815  if(filePath.Contains("High"))
816  {
817  gFracPerEn[ism][imax]->GetHistogram()->SetAxisRange(16,45,"X");
818  }
819  else if(filePath.Contains("Low"))
820  {
821  gFracPerEn[ism][imax]->GetHistogram()->SetAxisRange(7,20,"X");
822  }
823  else if(filePath.Contains("MB"))
824  {
825  gFracPerEn[ism][imax]->GetHistogram()->SetAxisRange(2,10,"X");
826  }
827  } // simu
828  else
829  {
830  if( filePath.Contains("EMC") || filePath.Contains("DMC") ||
831  filePath.Contains("EG") || filePath.Contains("DG") )
832  {
833  gFracPerEn[ism][imax]->GetHistogram()->SetAxisRange(6,50,"X");
834  }
835  else if ( filePath.Contains("MB") || filePath.Contains("INT") )
836  {
837  gFracPerEn[ism][imax]->GetHistogram()->SetAxisRange(2,10,"X");
838  if(filePath.Contains("LHC17"))
839  {
840  gFracPerEn[ism][imax]->GetHistogram()->SetAxisRange(2,14,"X");
841  }
842  }
843  } // simu
844  } // ism
845 
846  } // imax
847 
848  // --------------------------------------------------
849  // Plot fractions graphs, each pad one max Shower shape cut
850  // --------------------------------------------------
851 
852 // Float_t colorSM[] = {2,kBlue-3,kBlue+3,kRed-3,8,kYellow-2,kOrange-2,kRed+3,kCyan, kViolet,1};
853 // Float_t styleSM[] = {21,24,24,22,21,21,21,22,21,21,20};
854 
855  Int_t colorSM[]={1 , 1, 2, 2, 3, 3, 4, 4, 7, 7, 6, 6, 2, 3, 4, 7, 6, 2, 2, 3, 3, 4, 4, 6, 6};
856  Int_t styleSM[]={24,25,25,24,25,24,25,24,25,24,25,21,21,21,21,21,22,26,22,26,22,26,22,26};
857 
858  TCanvas * cGraphE = new TCanvas(Form("cGraphE_iana%d_%s",iana, titleName.Data()),
859  Form("Graph E bin, iana %d, %s",iana, titleName.Data()),
860  2*2000,2*2000);
861 
862  cGraphE->Divide(2,2);
863 
864  TLegend *lSM = new TLegend(0.1,0.1,1,1);
865  lSM->SetFillColor(0);
866  lSM->SetFillStyle(0);
867  lSM->SetLineColor(0);
868  lSM->SetBorderSize(0);
869  lSM->SetTextSize(0.07);
870 
871  for(Int_t imax = 0; imax < nShShMaxCut; imax++)
872  {
873  cGraphE->cd(imax+2);
874 
875  //gPad->SetGridy();
876  //gPad->SetGridx();
877 
878  gFracPerEn[nSM][imax]->SetTitle(Form("0.10<#sigma_{long}^{2}<%2.2f, %s",maxShShCut[imax],anaTitle[iana].Data()));
879  gFracPerEn[nSM][imax]->Draw("AP");
880 
881  for(Int_t ism = firstSM; ism < nSM+1; ism++)
882  {
883  gFracPerEn[ism][imax]->Draw("P");
884  gFracPerEn[ism][imax]->GetYaxis()->SetTitle("#it{R}_{#sigma_{long}}^{#it{n}^{w}_{cell}} (%)");
885  if ( iana < 2 ) gFracPerEn[ism][imax]->GetXaxis()->SetTitle("#it{E}_{T} (GeV)");
886  else gFracPerEn[ism][imax]->GetXaxis()->SetTitle("#it{E} (GeV)");
887  gFracPerEn[ism][imax]->SetMaximum(50);
888  gFracPerEn[ism][imax]->SetMinimum(0);
889 
890  if(filePath.Contains("Default"))
891  {
892  gFracPerEn[ism][imax]->SetMaximum(5);
893  gFracPerEn[ism][imax]->SetMinimum(0);
894  }
895 
896  gFracPerEn[ism][imax]->SetMarkerColor(colorSM[ism]);
897  gFracPerEn[ism][imax]->SetLineColor (colorSM[ism]);
898  gFracPerEn[ism][imax]->SetMarkerStyle(styleSM[ism]);
899  gFracPerEn[ism][imax]->SetMarkerSize (5);
900  if(ism!=nSM)
901  {
902  if ( imax == 0 ) lSM->AddEntry(gFracPerEn[ism][imax],Form("SM %d",ism),"PL");
903  }
904  else
905  {
906  gFracPerEn[ism][imax]->SetMarkerColor(1);
907  gFracPerEn[ism][imax]->SetLineColor (1);
908  gFracPerEn[ism][imax]->SetMarkerStyle(20);
909  if ( imax == 0 ) lSM->AddEntry(gFracPerEn[ism][imax],"All SM","PL");
910  }
911  }
912  }
913  cGraphE->cd(1);
914 
915  lSM->Draw();
916 
917  fileName = Form("figures/Comparison_Fraction_M02PhotonToAll_NCellCut_YaxisE_PerSM_%s_%s",
918  anaName[iana].Data(),titleName.Data());
919 
920  fileName+=fileFormat;
921  cGraphE->Print(fileName);
922 
923 
924  gStyle->SetOptTitle(1);
925 
926  // Plot each max cut individually
927  //
928  for(Int_t imax = 0; imax < nShShMaxCut; imax++)
929  {
930  TCanvas * cGraphEMax = new TCanvas(Form("cGraphE_Sh%d_iana%d_%s",imax,iana, titleName.Data()),
931  Form("Graph E bin, Sh < %2.2f, iana %d, %s", maxShShCut[imax], iana, titleName.Data()),
932  1*2000,1*2000);
933 
934  //gPad->SetGridy();
935  //gPad->SetGridx();
936 
937  //gFracPerEn[nSM][imax]->SetMaximum(25);
938  //gFracPerEn[nSM][imax]->SetMinimum(0);
939 
940  TLegend *lSM = new TLegend(0.85,0.2,0.98,0.85);
941  lSM->SetFillColor(1);
942  lSM->SetFillStyle(1);
943  lSM->SetLineColor(0);
944  lSM->SetBorderSize(0);
945  lSM->SetTextSize(0.03);
946 
947  TH1F * hAxis = new TH1F(Form("hAxis_Max%d",imax),
948  Form("0.1<#sigma_{long}^{2}<%2.2f, %s",maxShShCut[imax], anaTitle[iana].Data()),60,5,65);
949  hAxis->SetYTitle("#it{R}_{#sigma_{long}}^{#it{n}^{w}_{cell}} (%)");
950  if ( iana < 2 ) hAxis->SetXTitle("#it{E}_{T} (GeV)");
951  else hAxis->SetXTitle("#it{E} (GeV)");
952  //hAxis->SetTitle(Form("0.1<#sigma_{long}^{2}<0.3, %s",anaTitle[iana].Data()));
953  hAxis->SetMaximum(50);
954  hAxis->SetMinimum(0);
955  hAxis->Draw("");
956 
957  gFracPerEn[nSM][imax]->Draw("P");
958  lSM->AddEntry(gFracPerEn[nSM][imax],"All SM","PL");
959 
960  for(Int_t ism = firstSM; ism < nSM+1; ism++)
961  {
962  gFracPerEn[ism][imax]->Draw("P");
963  if(ism < nSM) lSM->AddEntry(gFracPerEn[ism][imax],Form("SM %d",ism),"PL");
964  }
965 
966  lSM->Draw();
967 
968  fileName = Form("figures/Comparison_Fraction_M02PhotonToAll_NCellCut_YaxisE_PerSM_ShMax%d_%s_%s",
969  imax, anaName[iana].Data(),titleName.Data());
970 
971  fileName+=fileFormat;
972  cGraphEMax->Print(fileName);
973  }
974 
975  // --------------------------------------------------
976  // Plot fractions graphs vs SM, each pad one max Shower shape cut
977  // --------------------------------------------------
978  TCanvas * cGraphSM = new TCanvas(Form("cGraphSM_iana%d_%s",iana, titleName.Data()),
979  Form("Graph E bin, iana %d, %s",iana, titleName.Data()),
980  2*2000,2*2000);
981 
982  cGraphSM->Divide(2,2);
983 
984  TLegend *lE = new TLegend(0.1,0.1,1,1);
985  lE->SetFillColor(0);
986  lE->SetFillStyle(0);
987  lE->SetLineColor(0);
988  lE->SetBorderSize(0);
989  lE->SetTextSize(0.07);
990 
991  for(Int_t imax = 0; imax < nShShMaxCut; imax++)
992  {
993  cGraphSM->cd(imax+2);
994 
995  gFracPerSM[0][imax]->SetTitle(Form("0.10<#sigma_{long}^{2}<%2.2f, %s",maxShShCut[imax],anaTitle[iana].Data()));
996  gFracPerSM[0][imax]->Draw("AP");
997 
998  for(Int_t ie = 0; ie < nEBins; ie++)
999  {
1000  gFracPerSM[ie][imax]->Draw("P");
1001  gFracPerSM[ie][imax]->GetYaxis()->SetTitle("Fraction %");
1002  gFracPerSM[ie][imax]->GetXaxis()->SetTitle("SM");
1003  gFracPerSM[ie][imax]->SetMaximum(50);
1004  gFracPerSM[ie][imax]->SetMinimum(0);
1005  if ( filePath.Contains("Default") )
1006  {
1007  gFracPerSM[ie][imax]->SetMaximum(5);
1008  gFracPerSM[ie][imax]->SetMinimum(0);
1009  }
1010  gFracPerSM[ie][imax]->SetMarkerColor(color[ie]);
1011  gFracPerSM[ie][imax]->SetLineColor (color[ie]);
1012  gFracPerSM[ie][imax]->SetMarkerStyle(24);
1013  gFracPerSM[ie][imax]->SetMarkerSize (5);
1014  if ( imax == 0 ) lE->AddEntry(gFracPerSM[ie][imax],Form("%2.f <E %2.1f GeV",lowE[ie],higE[ie]),"PL");
1015  }
1016  }
1017  cGraphSM->cd(1);
1018  lE->Draw();
1019 
1020  fileName = Form("figures/Comparison_Fraction_M02PhotonToAll_NCellCut_YaxisSM_PerE_%s_%s",
1021  anaName[iana].Data(),titleName.Data());
1022 
1023  fileName+=fileFormat;
1024  cGraphSM->Print(fileName);
1025 
1026  //****************************************************
1027  // Save fractions in histogram
1028  //****************************************************
1029  TFile *fout = new TFile(Form("figures/Fraction_M02PhotonToAll_NCellCut_%s_%s.root",
1030  anaName[iana].Data(),titleName.Data()),"recreate");
1031  for(Int_t imax = 0; imax < nShShMaxCut; imax++)
1032  {
1033  for(Int_t ism = firstSM; ism < nSM+1; ism++)
1034  {
1035  gFracPerEn[ism][imax]->SetName(Form("gFractionFrom%1.2f_SM%d_PerE",maxShShCut[imax], ism));
1036  gFracPerEn[ism][imax]->Write();
1037  }
1038 
1039  for(Int_t ie = 0; ie < nEBins; ie++)
1040  {
1041  gFracPerSM[ie][imax]->SetName(Form("gFractionFrom%1.2f_Ebin%d_PerSM",maxShShCut[imax],ie));
1042  gFracPerSM[ie][imax]->Write();
1043  }
1044  }
1045  fout->Print();
1046  fout->Close();
1047 
1048  file->Close();
1049  delete file;
1050 }
1051 
1052 
1053 //------------------------------------------------------------------------------
1063 //------------------------------------------------------------------------------
1064 void CompareProductions
1065 ( Int_t nProd,
1066  TString * fileName,
1067  TString * legName,
1068  Int_t iana = 0,
1069  Bool_t debug = kFALSE)
1070 {
1071  // SM
1072  //
1073  Int_t totalSM = 10, firstSM = 0;
1074  Int_t ncolSM = 4, nrowSM = 3;
1075  GetSMNumber(fileName[0].Data(),totalSM,firstSM,ncolSM,nrowSM);
1076  const Int_t nSM = totalSM+1; // Add last SM containing the sum of all
1077 
1078  if ( debug )
1079  printf("Compare %d productions for ana %d, nSM %d+1, first SM %d, ncol %d, nrow %d\n",
1080  nProd,iana,totalSM,firstSM,ncolSM,nrowSM);
1081 
1082  // Init arrays
1083  TFile * file[nProd];
1084  TGraphErrors * gFrac[nSM][nShShMaxCut][nProd];
1085  TH1F * hAxis[nShShMaxCut][nSM];
1086 
1087  // Get the files with graphs
1088  //
1089  for(Int_t iprod = 0; iprod < nProd; iprod++)
1090  {
1091  file[iprod] = new TFile(Form("figures/Fraction_M02PhotonToAll_NCellCut_%s_%s.root",
1092  anaName[iana].Data(),fileName[iprod].Data()),"read");
1093 
1094  if ( debug )
1095  {
1096  printf("figures/Fraction_M02PhotonToAll_NCellCut_%s_%s.root \n",
1097  anaName[iana].Data(),fileName[iprod].Data());
1098  printf("iprod %d, file %p\n",iprod,file[iprod]);
1099  }
1100 
1101  if ( !file[iprod] ) continue;
1102 
1103  for(Int_t ism = firstSM; ism < nSM; ism++)
1104  {
1105  for(Int_t im02 = 0; im02 < nShShMaxCut; im02++)
1106  {
1107  gFrac[ism][im02][iprod] = (TGraphErrors*) file[iprod]->Get(Form("gFractionFrom%1.2f_SM%d_PerE",maxShShCut[im02],ism));
1108 
1109  if ( debug )
1110  {
1111  printf("\t gFractionFrom%1.2f_SM%d_PerE \n",maxShShCut[im02],ism);
1112  printf("\t ism %d, im02 %d graph %p \n",ism,im02,gFrac[ism][im02]);
1113  }
1114 
1115  if ( !gFrac[ism][im02][iprod] ) continue;
1116 
1117  // Set the plotting ranges, removing the unwanted points
1118  // Depending on the production name
1119  if(fileName[iprod].Contains("MB") || fileName[iprod].Contains("INT7"))
1120  RemovePointsOutOfRangeOrLargeErrorFromGraph(gFrac[ism][im02][iprod], 2., 10.); // PlotUtils.C
1121  else if( fileName[iprod].Contains("LHC11") && fileName[iprod].Contains("EMC7") )
1122  RemovePointsOutOfRangeOrLargeErrorFromGraph(gFrac[ism][im02][iprod], 6., 60.); // PlotUtils.C
1123 
1124  if(fileName[iprod].Contains("JJ"))
1125  {
1126  if(fileName[iprod].Contains("High"))
1127  RemovePointsOutOfRangeOrLargeErrorFromGraph(gFrac[ism][im02][iprod], 14., 60.); // PlotUtils.C
1128  else if(fileName[iprod].Contains("Low"))
1129  RemovePointsOutOfRangeOrLargeErrorFromGraph(gFrac[ism][im02][iprod], 7., 20.); // PlotUtils.C
1130  }
1131 
1132  } // im02
1133 
1134  } // ism
1135 
1136  } // iprod
1137 
1138 
1139  //
1140  // PLOTS
1141  //
1142  Float_t colorProd[] = { 1, 1,kBlue, kBlue,kBlue,kRed,kRed,kRed, 8, 8, 8};
1143  Float_t styleProd[] = {20,22, 26, 24, 25, 26, 24, 25, 26, 24, 25};
1144 
1145  gStyle->SetPadRightMargin(0.01);
1146  gStyle->SetPadLeftMargin(0.12);
1147  gStyle->SetTitleFontSize(0.06);
1148  // gStyle->SetStatFontSize(0.06);
1149  gStyle->SetOptStat(0);
1150  gStyle->SetOptTitle(1);
1151  gStyle->SetPadTopMargin(0.1);
1152 
1153  //
1154  // Plot: Each frame is a M02 cut, each file is per one SM, SM10 is the sum of all.
1155  //
1156  Int_t ncol = 2;
1157  Int_t nrow = 2;
1158 
1159  if(nShShMaxCut==1)
1160  {
1161  ncol=1;
1162  nrow=1;
1163  }
1164 
1165  for(Int_t ism = firstSM; ism < nSM; ism++)
1166  {
1167  TCanvas * cGraphSM = new TCanvas(Form("cGraphSM_iana%d_SM%d",iana,ism),
1168  Form("iana %d, SM%d",iana, ism),
1169  ncol*2000,nrow*2000);
1170 
1171  cGraphSM->Divide(ncol,nrow);
1172 
1173  gPad->GetLogx();
1174  gPad->GetLogy();
1175 
1176  TLegend *lE = 0;
1177  if(nShShMaxCut == 3)
1178  {
1179  lE = new TLegend(0.1,0.1,1,1);
1180  lE->SetTextSize(0.05);
1181  }
1182  else
1183  {
1184  lE = new TLegend(0.4,0.55,0.95,0.9);
1185  lE->SetTextSize(0.035);
1186  }
1187  lE->SetFillColor(0);
1188  lE->SetFillStyle(0);
1189  lE->SetLineColor(0);
1190  lE->SetBorderSize(0);
1191  //lE->SetHeader(Form("0.10<#sigma^{2}_{long}<%1.2f, %s", maxShShCut[im02], anaTitle[iana].Data()));
1192 
1193  if ( nShShMaxCut > 1 )
1194  {
1195  if (ism==10) lE->SetHeader(Form(" All SM, %s", anaTitle[iana].Data()));
1196  else lE->SetHeader(Form(" SM %d, %s",ism, anaTitle[iana].Data()));
1197  }
1198 
1199  for(Int_t im02 = 0; im02 < nShShMaxCut; im02++)
1200  {
1201  cGraphSM->cd(im02+1);
1202 
1203  gPad->GetLogx();
1204  gPad->GetLogy();
1205 
1206  TString title = Form("0.10<#sigma^{2}_{long}<%1.2f", maxShShCut[im02]);
1207  if ( nShShMaxCut == 1 )
1208  {
1209  if (ism == 10 ) title = Form("0.10<#sigma^{2}_{long}<%1.2f, all SM, %s", maxShShCut[im02],anaTitle[iana].Data());
1210  else title = Form("0.10<#sigma^{2}_{long}<%1.2f, SM %d, %s", maxShShCut[im02],ism, anaTitle[iana].Data());
1211  }
1212 
1213  hAxis[im02][ism] = new TH1F(Form("hAxisBis_SM%d_M02%d",ism,im02),title.Data(),1000,2,50);
1214 
1215  hAxis[im02][ism]->GetYaxis()->SetTitle("#it{R}_{#sigma_{long}}^{#it{n}^{w}_{cell}} (%)");
1216  if ( iana < 2 ) hAxis[im02][ism]->GetXaxis()->SetTitle("#it{E} (GeV)");
1217  else hAxis[im02][ism]->GetXaxis()->SetTitle("#it{E}_{T} (GeV/#it{c})");
1218 
1219  if(ism!=3 && ism !=7 ) hAxis[im02][ism]->SetMaximum(45);
1220  else hAxis[im02][ism]->SetMaximum(65);
1221 
1222  hAxis[im02][ism]->SetMinimum(0);
1223  if( (ism == 3 || ism ==7) && im02 == 2 ) hAxis[im02][ism]->SetMaximum(45);
1224 
1225  hAxis[im02][ism]->Draw("");
1226 
1227  for(Int_t iprod = 0; iprod < nProd; iprod++)
1228  {
1229  if(!file[iprod]) continue;
1230  if(! gFrac[ism][im02][iprod] ) continue;
1231 
1232  gFrac[ism][im02][iprod]->SetMarkerColor(colorProd[iprod]);
1233  gFrac[ism][im02][iprod]->SetLineColor (colorProd[iprod]);
1234  gFrac[ism][im02][iprod]->SetMarkerStyle(styleProd[iprod]);
1235  gFrac[ism][im02][iprod]->SetMarkerSize (5);
1236  if(im02==0) lE->AddEntry(gFrac[ism][im02][iprod],legName[iprod],"PL");
1237 
1238  gFrac[ism][im02][iprod]->Draw("P");
1239  }
1240 
1241  gFrac[ism][im02][0]->Draw("P");
1242  }
1243 
1244  if(nShShMaxCut == 3) cGraphSM->cd(4);
1245  lE->Draw();
1246 
1247  TString fileName = Form("figures/Comparison_Fraction_M02PhotonToAll_NCellCut_XaxisE_PerM02Cut_PerProd_%s_",
1248  anaName[iana].Data());
1249 
1250  if(ism == 10 ) fileName+="AllSM";
1251  else fileName+=Form("SM%d",ism);
1252 
1253  fileName+=fileFormat;
1254  cGraphSM->Print(fileName);
1255  }
1256 
1257  //
1258  // Plot: Each frame is a SM, each file is per one M02 cut
1259  //
1260  for(Int_t im02 = 0; im02 < nShShMaxCut; im02++)
1261  {
1262  TCanvas * cGraphE = new TCanvas(Form("cGraphE_iana%d_%d",iana,im02),
1263  Form("iana %d, %2.2f",iana, maxShShCut[im02]),
1264  ncolSM*2000,nrowSM*2000);
1265 
1266  gPad->GetLogx();
1267  gPad->GetLogy();
1268 
1269  cGraphE->Divide(ncolSM,nrowSM);
1270 
1271  TLegend *lSM = new TLegend(-0.04,0.1,1,1);
1272  lSM->SetFillColor(0);
1273  lSM->SetFillStyle(0);
1274  lSM->SetLineColor(0);
1275  lSM->SetBorderSize(0);
1276  lSM->SetTextSize(0.06);
1277  lSM->SetHeader(Form(" 0.10<#sigma^{2}_{long}<%1.2f, %s", maxShShCut[im02], anaTitle[iana].Data()));
1278 
1279  for(Int_t ism = 0; ism < nSM; ism++)
1280  {
1281  cGraphE->cd(ism+1);
1282 
1283  gPad->GetLogx();
1284  gPad->GetLogy();
1285 
1286  TString title = Form("SM %d",ism);
1287  if(ism == 10)
1288  title = "All SM";
1289 
1290  hAxis[im02][ism] = new TH1F(Form("hAxis_SM%d_M02%d",ism,im02),title.Data(),1000,2,50);
1291 
1292  hAxis[im02][ism]->GetYaxis()->SetTitle("#it{R}_{#sigma_{long}}^{#it{n}^{w}_{cell}} (%)");
1293  if ( iana < 2 ) hAxis[im02][ism]->GetXaxis()->SetTitle("#it{E} (GeV)");
1294  else hAxis[im02][ism]->GetXaxis()->SetTitle("#it{E}_{T} (GeV/#it{c})");
1295 
1296  hAxis[im02][ism]->SetMaximum(35);
1297  hAxis[im02][ism]->SetMinimum(0);
1298 
1299  hAxis[im02][ism]->Draw("");
1300 
1301  for(Int_t iprod = 0; iprod < nProd; iprod++)
1302  {
1303  if(!file[iprod]) continue;
1304  if(! gFrac[ism][im02][iprod] ) continue;
1305 
1306  gFrac[ism][im02][iprod]->SetMarkerColor(colorProd[iprod]);
1307  gFrac[ism][im02][iprod]->SetLineColor (colorProd[iprod]);
1308  gFrac[ism][im02][iprod]->SetMarkerStyle(styleProd[iprod]);
1309  gFrac[ism][im02][iprod]->SetMarkerSize (5);
1310  if(ism==0) lSM->AddEntry(gFrac[ism][im02][iprod],legName[iprod],"PL");
1311 
1312  gFrac[ism][im02][iprod]->Draw("P");
1313  }
1314 
1315  gFrac[ism][im02][0]->Draw("P");
1316  }
1317 
1318  cGraphE->cd(nrowSM*ncolSM);
1319  lSM->Draw();
1320 
1321  TString fileName = Form("figures/Comparison_Fraction_M02PhotonToAll_NCellCut_XaxisE_PerSM_PerProd_%s_M02Cut%d",
1322  anaName[iana].Data(),im02);
1323 
1324  fileName+=fileFormat;
1325  cGraphE->Print(fileName);
1326  }
1327 }
1328 
1329 //------------------------------------------------------------------------------
1341 //------------------------------------------------------------------------------
1342 void CalculateShowerShapeFractionPerNCellCut
1343 (
1344  Bool_t doCalc = kTRUE,
1345  Bool_t doComp = kTRUE,
1346  TString tm = "_TMDep",
1347  Bool_t plotRat = kFALSE,
1348  Bool_t debug = kFALSE
1349  )
1350 {
1351  const Int_t nProd = 8; // Make sure number is equal to entries below, if not it crashes
1352  TString filePath[nProd];
1353  TString title [nProd];
1354  TString legend [nProd];
1355 
1356  if ( debug ) printf("N prod %d\n",nProd);
1357 
1358  Int_t nproditer = 0;
1359 
1360  filePath[nproditer++] = "data/module/TCardChannel3/LHC11cd_EMC7";
1361  title [nproditer-1] = "LHC11cd_EMC7";
1362  legend [nproditer-1] = "data, LHC11c+d, EMC7";
1363  filePath[nproditer++] = "data/module/TCardChannel3/LHC11cd_INT7";
1364  title [nproditer-1] = "LHC11cd_INT7";
1365  legend [nproditer-1] = "data, LHC11c+d, INT7";
1366 
1367  filePath[nproditer++] = "simu/module/pp_7TeV_MB/TCardChannel_Mimic0_Scaled2_v3/ScaledMerged";
1368  title [nproditer-1] = "default_MC_MB";
1369  legend [nproditer-1] = "MC default, MB";
1370  filePath[nproditer++] = "simu/module/pp_7TeV_MB/TCardChannel_Mimic10c_EcellCut_Scaled2_v3/ScaledMerged";
1371  legend [nproditer-1] = "MC mimic, MB";
1372 
1373  filePath[nproditer++] = "simu/module/pp_7TeV_JJ_Dec_GJ/TCardChannel_Mimic0_Scaled2_v3/ScaledMerged";
1374  title [nproditer-1] = "Default_MC_JJ_DecLow";
1375  legend [nproditer-1] = "MC default, #gammaJ-JJ, #it{p}^{EMC}_{T}>3.5 GeV/#it{c}";
1376  filePath[nproditer++] = "simu/module/pp_7TeV_JJ_Dec_GJ/TCardChannel_Mimic10c_EcellCut_Scaled2_v3/ScaledMerged";
1377  title [nproditer-1] = "Mimic_MC_JJ_DecLow";
1378  legend [nproditer-1] = "MC mimic, #gammaJ-JJ, #it{p}^{EMC}_{T}>3.5 GeV/#it{c}";
1379 
1380  filePath[nproditer++] = "simu/module/pp_7TeV_JJ_Dec_High_GJ/TCardChannel_Mimic0_Scaled2_v3/ScaledMerged";
1381  title [nproditer-1] = "Default_MC_JJ_DecHigh";
1382  legend [nproditer-1] = "MC default, #gammaJ-JJ, #it{p}^{EMC}_{T}>7 GeV/#it{c}";
1383  filePath[nproditer++] = "simu/module/pp_7TeV_JJ_Dec_High_GJ/TCardChannel_Mimic10c_EcellCut_Scaled2_v3/ScaledMerged";
1384  title [nproditer-1] = "Mimic_MC_JJ_DecHigh";
1385  legend [nproditer-1] = "MC mimic, #gammaJ-JJ, #it{p}^{EMC}_{T}>7 GeV/#it{c}";
1386 
1387  if ( doCalc )
1388  {
1389  for(Int_t iprod = 0; iprod < nProd; iprod++)
1390  {
1391  CalculateAndPlot(title[iprod],filePath[iprod],0,tm, plotRat, debug); // Not isolated
1392  CalculateAndPlot(title[iprod],filePath[iprod],1,tm, plotRat, debug); // Isolated
1393  CalculateAndPlot(title[iprod],filePath[iprod],2,tm, plotRat, debug); // Inclusive
1394  }
1395  }
1396 
1397  if(doComp)
1398  {
1399  CompareProductions(nProd,title,legend,0,debug);
1400  CompareProductions(nProd,title,legend,1,debug);
1401  CompareProductions(nProd,title,legend,2,debug);
1402  }
1403 }
histoLeg
TString histoLeg[]
Definition: CalculateShowerShapeFractionPerNCellCut.C:49
Double_t
double Double_t
Definition: External.C:58
TH3F
Definition: External.C:260
CompareProductions
void CompareProductions(Int_t nProd, TString *fileName, TString *legName, Int_t iana=0, Bool_t debug=kFALSE)
Definition: CalculateShowerShapeFractionPerNCellCut.C:1065
TH2F
Definition: External.C:236
title
const char * title
Definition: MakeQAPdf.C:27
fileFormat
TString fileFormat
Definition: CalculateShowerShapeFractionPerNCellCut.C:74
fileName
TString fileName
Definition: ProjectCombinHFAndFit.C:25
GetSMNumber
void GetSMNumber(TString titleName, Int_t &totalSM, Int_t &firstSM, Int_t &ncol, Int_t &nrow)
Definition: CalculateShowerShapeFractionPerNCellCut.C:85
c
TCanvas * c
Definition: TestFitELoss.C:172
PlotUtils.C
Plotting utilities.
CalculateAndPlot
void CalculateAndPlot(TString titleName, TString filePath, Int_t iana=0, TString tm="_TMDep", Bool_t plotRat=kFALSE, Bool_t debug=kFALSE)
Definition: CalculateShowerShapeFractionPerNCellCut.C:122
GetFractionError
static Double_t GetFractionError(Double_t num, Double_t den, Double_t numErr, Double_t denErr)
Definition: PlotUtils.C:132
RemovePointsOutOfRangeOrLargeErrorFromGraph
void RemovePointsOutOfRangeOrLargeErrorFromGraph(TGraphErrors *graph, Double_t min, Double_t max, Float_t errFraction=0.5, Float_t value=-1, Float_t valueErr=0)
Definition: PlotUtils.C:257
Int_t
int Int_t
Definition: External.C:63
nProd
const Int_t nProd
Definition: DrawProductionComparison.C:102
Float_t
float Float_t
Definition: External.C:68
filePath
TString filePath[]
Definition: CompareShowerShapeLongPerAnaPerSM.C:61
nShShMaxCut
const Int_t nShShMaxCut
Definition: CalculateShowerShapeFractionPerNCellCut.C:57
marker
Int_t marker[]
Definition: CalculateShowerShapeFractionPerNCellCut.C:71
TH1D
Definition: External.C:212
GetRangeIntegralAndError
static void GetRangeIntegralAndError(TH1D *h, Int_t binMin, Int_t binMax, Double_t &integral, Double_t &integralErr)
Definition: PlotUtils.C:208
nEBins
const Int_t nEBins
Definition: CalculateShowerShapeFractionPerNCellCut.C:52
CalculateShowerShapeFractionPerNCellCut
void CalculateShowerShapeFractionPerNCellCut(Bool_t doCalc=kTRUE, Bool_t doComp=kTRUE, TString tm="_TMDep", Bool_t plotRat=kFALSE, Bool_t debug=kFALSE)
Definition: CalculateShowerShapeFractionPerNCellCut.C:1343
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
shshRef
Int_t shshRef
Definition: CalculateShowerShapeFractionPerNCellCut.C:63
maxShShCut
Float_t maxShShCut[]
Definition: CalculateShowerShapeFractionPerNCellCut.C:58
anaName
TString anaName[]
Definition: CalculateShowerShapeFractionPerNCellCut.C:44
titleName
TString titleName
Definition: CompareShowerShapeLongPerAnaPerSM.C:60
file
TFile * file
TList with histograms for a given trigger.
Definition: DrawAnaCaloTrackQA.C:89
rebin
Int_t rebin
Definition: charmCutsOptimization.C:32
binE
Double_t binE[]
Definition: CalculateShowerShapeFractionPerNCellCut.C:53
GetCanvasColRowNumber
void GetCanvasColRowNumber(Int_t npad, Int_t &ncol, Int_t &nrow, Bool_t square=kFALSE)
Definition: PlotUtils.C:304
Bool_t
bool Bool_t
Definition: External.C:53
fout
TFile * fout
input train file
Definition: DrawAnaCaloTrackQA.C:90
nHisto
const Int_t nHisto
Definition: CalculateShowerShapeFractionPerNCellCut.C:48
binEErr
Double_t binEErr[]
Definition: CalculateShowerShapeFractionPerNCellCut.C:54
lineStyle
Int_t lineStyle[]
Definition: CalculateShowerShapeFractionPerNCellCut.C:68
anaTitle
TString anaTitle[]
Definition: CalculateShowerShapeFractionPerNCellCut.C:45
color
Int_t color[]
Definition: CalculateShowerShapeFractionPerNCellCut.C:66