AliPhysics  9c66e61 (9c66e61)
CompareInvMassGraphs.C
Go to the documentation of this file.
1 
14 #if !defined(__CINT__) || defined(__MAKECINT__)
15 
16 #include <TFile.h>
17 #include <TDirectoryFile.h>
18 #include <TList.h>
19 #include <TString.h>
20 #include <TROOT.h>
21 #include <TStyle.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 "PlotUtils.C"
36 
37 #endif
38 
39 TString fileFormat = ".eps";
40 
41 //------------------------------------------------------------------------------
48 //------------------------------------------------------------------------------
50 (
51  TString particle = "Pi0",
52  TString calorimeter = "EMCAL",
53  TString titleMC = "simu_pp_7TeV_MB",
54  TString titleData = "LHC11cd_INT7"
55 )
56 {
57  // SM
58  const Int_t nSM = 10; // 10 SM + 1 sum all
59  Int_t ncol = 4;
60  Int_t nrow = 3;
61 
62  TString daLeg = "data";
63  if (titleData .Contains("LHC11cd_EMC7"))
64  daLeg = "pp@7 TeV, LHC11c+d EMC7";
65  else if(titleData .Contains("LHC11cd_INT7"))
66  daLeg = "pp@7 TeV, LHC11c+d INT7";
67  else if(titleData == "LHC12_EMC7")
68  daLeg = "pp@8 TeV, LHC12 EMC7";
69  else if(titleData == "LHC12_EMG1")
70  daLeg = "pp@8 TeV, LHC12 EGA";
71  else if(titleData == "LHC12_INT7")
72  daLeg = "pp@8 TeV, LHC12 INT7";
73  else if(titleData == "LHC17pq_EMC7")
74  daLeg = "pp@5 TeV, LHC11p+q EMC7";
75  else if(titleData == "LHC17pq_EMG1")
76  daLeg = "pp@5 TeV, LHC11p+q EGA";
77  else if(titleData == "LHC17pq_INT7")
78  daLeg = "pp@5 TeV, LHC17p+q INT7";
79  else if(titleData == "LHC17pq_DCAL_DMC7")
80  daLeg = "pp@5 TeV, LHC11p+q DMC7";
81  else if(titleData == "LHC17pq_DCAL_DMG1")
82  daLeg = "pp@5 TeV, LHC11p+q DGA";
83  else if(titleData == "LHC17pq_DCAL_INT7")
84  daLeg = "pp@5 TeV, LHC17p+q INT7";
85 
86  //------------------------------------------------
87  // Put the path to the files in prod string array
88  // also the explanation title on prodLeg string array
89  //------------------------------------------------
90 
91  const Int_t nProd = 3;
92 
93  TString prod[] =
94  {
95  Form("data_%s" ,titleData.Data())
96  , Form("%s_Mimic0_Scaled" ,titleMC .Data())
97  , Form("%s_Mimic10c_EcellCut_Scaled",titleMC .Data())
98  };
99 
100  TString prodLeg[] =
101  {
102  "Data"
103  , "MC default"
104  , "MC xTalk + #it{E}_{inc+cell}>100 MeV"
105  };
106 
107  // Init the graph arrays.
108  //
109  TGraphErrors* gMass[nSM][nProd];
110  TGraphErrors* gMassAllSM[nProd];
111  TGraphErrors* gMassTRDNo[nProd];
112  TGraphErrors* gMassTRDYe[nProd];
113  TH1F* hAxis[nSM];
114 
115  TFile* file[nProd];
116 
117  TGraphErrors* gMassRat[nSM][nProd];
118  TGraphErrors* gMassRatAllSM[nProd];
119  TGraphErrors* gMassRatTRDNo[nProd];
120  TGraphErrors* gMassRatTRDYe[nProd];
121 
122  //===============================
123  // Open the input files and recover the graphs
124  // apply some style settings
125  //===============================
126 
127  Float_t colorProd[] = { 1,kBlue,kRed,kViolet,kRed,kOrange-2,kCyan};
128  Float_t styleProd[] = {20, 24, 24, 24, 25, 25, 25};
129  Float_t markerSize = 8;
130 
131  titleData.ReplaceAll("/","_");
132 
133  for(Int_t iprod = 0; iprod < nProd; iprod++)
134  {
135  file[iprod] = new TFile(Form("IMfigures/%s_%s_MassWidthPtHistograms_%s.root",
136  prod[iprod].Data(),calorimeter.Data(),particle.Data()),"read");
137 
138  printf("iprod %d, IMfigures/%s_%s_MassWidthPtHistograms_%s.root %p\n",
139  iprod,prod[iprod].Data(),calorimeter.Data(),particle.Data(),file[iprod]);
140 
141  if(!file[iprod]) continue;
142 
143  gMassAllSM[iprod] = (TGraphErrors*) file[iprod]->Get("gMass_SameSM");
144  //printf("\t gMass_SameSM %p\n",gMassAllSM[iprod]);
145 
146  gMassAllSM[iprod]->SetMarkerColor(colorProd[iprod]);
147  gMassAllSM[iprod]->SetLineColor (colorProd[iprod]);
148  gMassAllSM[iprod]->SetMarkerStyle(styleProd[iprod]);
149  gMassAllSM[iprod]->SetMarkerSize (markerSize);
150 
151  if ( iprod > 0 && gMassAllSM[0] )
152  {
153  gMassRatAllSM[iprod] = DivideGraphs(gMassAllSM[iprod],gMassAllSM[0]);
154  gMassRatAllSM[iprod]->SetMarkerColor(colorProd[iprod]);
155  gMassRatAllSM[iprod]->SetLineColor (colorProd[iprod]);
156  gMassRatAllSM[iprod]->SetMarkerStyle(styleProd[iprod]);
157  gMassRatAllSM[iprod]->SetMarkerSize (markerSize);
158  }
159 
160  gMassTRDNo[iprod] = (TGraphErrors*) file[iprod]->Get("gMass_SameSMTRDNot");
161  //printf("\t gMass_SameSMTRDNot %p\n",gMassTRDNo[iprod]);
162 
163  if ( gMassTRDNo[iprod] )
164  {
165  gMassTRDNo[iprod]->SetMarkerColor(colorProd[iprod]);
166  gMassTRDNo[iprod]->SetLineColor (colorProd[iprod]);
167  gMassTRDNo[iprod]->SetMarkerStyle(styleProd[iprod]);
168  gMassTRDNo[iprod]->SetMarkerSize (markerSize);
169 
170  if ( iprod > 0 && gMassTRDNo[0] )
171  {
172  gMassRatTRDNo[iprod] = DivideGraphs(gMassTRDNo[iprod],gMassTRDNo[0]);
173  gMassRatTRDNo[iprod]->SetMarkerColor(colorProd[iprod]);
174  gMassRatTRDNo[iprod]->SetLineColor (colorProd[iprod]);
175  gMassRatTRDNo[iprod]->SetMarkerStyle(styleProd[iprod]);
176  gMassRatTRDNo[iprod]->SetMarkerSize (markerSize);
177  }
178  }
179 
180  gMassTRDYe[iprod] = (TGraphErrors*) file[iprod]->Get("gMass_SameSMTRDYes");
181  //printf("\t gMass_SameSMTRDYes %p\n",gMassTRDYe[iprod]);
182 
183  if ( gMassTRDYe[iprod] )
184  {
185  gMassTRDYe[iprod]->SetMarkerColor(colorProd[iprod]);
186  gMassTRDYe[iprod]->SetLineColor (colorProd[iprod]);
187  gMassTRDYe[iprod]->SetMarkerStyle(styleProd[iprod]);
188  gMassTRDYe[iprod]->SetMarkerSize (markerSize);
189 
190  if ( iprod > 0 && gMassTRDYe[0] )
191  {
192  gMassRatTRDYe[iprod] = DivideGraphs(gMassTRDYe[iprod],gMassTRDYe[0]);
193  gMassRatTRDYe[iprod]->SetMarkerColor(colorProd[iprod]);
194  gMassRatTRDYe[iprod]->SetLineColor (colorProd[iprod]);
195  gMassRatTRDYe[iprod]->SetMarkerStyle(styleProd[iprod]);
196  gMassRatTRDYe[iprod]->SetMarkerSize (markerSize);
197  }
198  }
199 
200  for(Int_t ism = 0; ism < nSM; ism++)
201  {
202  gMass[ism][iprod] = (TGraphErrors*) file[iprod]->Get(Form("gMass_SM%d",ism));
203  // printf("\t gMass_SM%d %p\n",ism,gMass[ism][iprod]);
204 
205  gMass[ism][iprod]->SetMarkerColor(colorProd[iprod]);
206  gMass[ism][iprod]->SetLineColor (colorProd[iprod]);
207  gMass[ism][iprod]->SetMarkerStyle(styleProd[iprod]);
208  gMass[ism][iprod]->SetMarkerSize (markerSize);
209 
210  if ( iprod > 0 && gMass[ism][0] )
211  {
212  gMassRat[ism][iprod] = DivideGraphs(gMass[ism][iprod],gMass[ism][0]);
213  gMassRat[ism][iprod]->SetMarkerColor(colorProd[iprod]);
214  gMassRat[ism][iprod]->SetLineColor (colorProd[iprod]);
215  gMassRat[ism][iprod]->SetMarkerStyle(styleProd[iprod]);
216  gMassRat[ism][iprod]->SetMarkerSize (markerSize);
217  }
218 
219  } // ism
220 
221  } // iprod
222 
223  //===============================
224  // PLOTS
225  //===============================
226  gStyle->SetPadRightMargin(0.01);
227  gStyle->SetPadLeftMargin(0.12);
228  gStyle->SetTitleFontSize(0.06);
229  // gStyle->SetStatFontSize(0.06);
230  gStyle->SetOptStat(0);
231  gStyle->SetOptTitle(1);
232  gStyle->SetPadTopMargin(0.07);
233 
234  //-----------------------------------------------------
235  // Plot per SM the Mass distributions and their ratios
236  //-----------------------------------------------------
238 
239  TLegend *lE = new TLegend(-0.04,0.1,1,1);
240  lE->SetFillColor(0);
241  lE->SetFillStyle(0);
242  lE->SetLineColor(0);
243  lE->SetBorderSize(0);
244  lE->SetTextSize(0.05);
245  //lE->SetHeader(Form(" SM %d, %s",ism, isoTitle[iso].Data()));
246 
247  TCanvas * cGraphSM = new TCanvas("cGraphSM","Per SM",
248  ncol*2000,nrow*2000);
249 
250  cGraphSM->Divide(ncol,nrow);
251 
252  for(Int_t ism = 0; ism < nSM; ism++)
253  {
254  cGraphSM->cd(ism+1);
255 
256  //gPad->SetGridx();
257  gPad->SetGridy();
258 
259  hAxis[ism] = new TH1F(Form("hAxisBisRat_SM%d",ism),Form("SM %d",ism),1000,1,12);
260 
261  hAxis[ism]->GetYaxis()->SetTitle("Mass (MeV/#it{c}^{2})");
262  hAxis[ism]->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
263  hAxis[ism]->SetMaximum(140);
264  hAxis[ism]->SetMinimum(115);
265 // hAxis[ism]->SetMaximum(137);
266 // hAxis[ism]->SetMinimum(123);
267 
268 // hAxis[ism]->SetMaximum(140);
269 // hAxis[ism]->SetMinimum(128);
270 
271  hAxis[ism]->Draw("");
272 
273  for(Int_t iprod = 0; iprod < nProd; iprod++)
274  {
275  if(!file[iprod]) continue;
276 
277  if ( ism == 0 )
278  lE->AddEntry(gMass[ism][iprod],prodLeg[iprod],"PL");
279 
280  gMass[ism][iprod]->Draw("P");
281  }
282 
283  }
284  cGraphSM->cd(nSM+1);
285  lE->Draw();
286 
287  fileName = Form("IMfigures/Comparison_Mass_%s_%s_%s_%s_PerSM",
288  titleData.Data(),titleMC.Data(),calorimeter.Data(),particle.Data());
289 
290  fileName+=fileFormat;
291  cGraphSM->Print(fileName);
292 
293  // Ratio
294  //
295  TCanvas * cGraphSMRat = new TCanvas("cGraphSMRat","Ratio Per SM",
296  ncol*2000,nrow*2000);
297 
298  cGraphSMRat->Divide(ncol,nrow);
299 
300  for(Int_t ism = 0; ism < nSM; ism++)
301  {
302  cGraphSMRat->cd(ism+1);
303 
304  //gPad->SetGridx();
305  gPad->SetGridy();
306 
307  hAxis[ism] = new TH1F(Form("hAxisBis_SM%d",ism),Form("SM %d",ism),1000,1,12);
308 
309  hAxis[ism]->GetYaxis()->SetTitle("Mass MC / Data");
310  hAxis[ism]->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
311 // hAxis[ism]->SetMaximum(1.002);
312 // hAxis[ism]->SetMinimum(0.932);
313 // hAxis[ism]->SetMaximum(1.05);
314 // hAxis[ism]->SetMinimum(0.90);
315  hAxis[ism]->SetMaximum(0.98);
316  hAxis[ism]->SetMinimum(0.85);
317  hAxis[ism]->Draw("");
318 
319  for(Int_t iprod = 1; iprod < nProd; iprod++)
320  {
321  if ( !gMassRat[ism][iprod] ) continue;
322 
323  //printf("%d %p\n",iprod,gMassRat[ism][iprod]);
324  gMassRat[ism][iprod]->Draw("P");
325  }
326 
327  }
328 
329  cGraphSMRat->cd(nSM+1);
330 
331  lE->Draw();
332 
333  fileName = Form("IMfigures/Comparison_MassRatio_%s_%s_%s_%s_PerSM",
334  titleData.Data(),titleMC.Data(),calorimeter.Data(),particle.Data());
335 
336  fileName+=fileFormat;
337  cGraphSMRat->Print(fileName);
338 
339 
340  //-----------------------------------------------------
341  // Plot sum of all SM or group of SM
342  // the Mass distributions and their ratios
343  //-----------------------------------------------------
344 
345  ncol = 2;
346  nrow = 2;
347 
348  if(gMassTRDYe[0])
349  {
350  ncol = 4;
351  nrow = 2;
352  }
353 
354  TH1F* hAxisAll = new TH1F("hAxisAll","All SM",1000,1,12);
355 
356  hAxisAll->GetYaxis()->SetTitle("Mass (MeV/#it{c}^{2})");
357  hAxisAll->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
358  hAxisAll->SetMaximum(140);
359  hAxisAll->SetMinimum(125);
360 
361 // hAxisAll->SetMaximum(140);
362 // hAxisAll->SetMinimum(115);
363 
364  TH1F* hAxisAllRat = new TH1F("hAxisAllRat","MC/Data ratio, All SM",1000,1,12);
365  hAxisAllRat->GetYaxis()->SetTitle("Mass Data/MC");
366  hAxisAllRat->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
367 // hAxisAllRat->SetMaximum(1.002);
368 // hAxisAllRat->SetMinimum(0.932);
369 
370 // hAxisAllRat->SetMaximum(1.02);
371 // hAxisAllRat->SetMinimum(0.92);
372  hAxisAllRat->SetMaximum(1.05);
373  hAxisAllRat->SetMinimum(0.98);
374 
375 
376  TCanvas * cGraph = new TCanvas("cGraphAllSM","All SM",
377  ncol*2000,nrow*2000);
378 
379  cGraph->Divide(ncol,nrow);
380 
381  cGraph->cd(1);
382 
383  //gPad->SetGridx();
384  gPad->SetGridy();
385 
386  hAxisAll->Draw("");
387 
388  for(Int_t iprod = 0; iprod < nProd; iprod++)
389  {
390  if(!file[iprod]) continue;
391  gMassAllSM[iprod]->Draw("P");
392  }
393 
394  gMassAllSM[1]->Draw("P");
395 
396  if ( gMassTRDYe[0] ) cGraph->cd(5);
397  else cGraph->cd(3);
398 
399  //gPad->SetGridx();
400  gPad->SetGridy();
401 
402 
403  TLegend *lA = new TLegend(0.2,0.8,0.9,0.9);
404  lA->SetFillColor(0);
405  lA->SetFillStyle(0);
406  lA->SetLineColor(0);
407  lA->SetBorderSize(0);
408  lA->SetTextSize(0.05);
409 
410  hAxisAllRat->Draw("");
411 
412  for(Int_t iprod = 1; iprod < nProd; iprod++)
413  {
414  if(!file[iprod]) continue;
415 
416  gMassRatAllSM[iprod]->Draw("P");
417  gMassRatAllSM[iprod]->Fit("pol0","R","",2,8);
418  TF1* func = gMassRatAllSM[iprod]->GetFunction("pol0");
419  lA->AddEntry(func,Form("Fit %2.3f #pm %2.3f",func->GetParameter(0),func->GetParError(0)),"L");
420  }
421  gMassRatAllSM[1]->Draw("P");
422 
423  lA->Draw();
424 
425  if ( gMassTRDNo[0] )
426  {
427  cGraph->cd(2);
428 
429  //gPad->SetGridx();
430  gPad->SetGridy();
431 
432  TH1F* hAxisNot = (TH1F*) hAxisAll->Clone("hAxisNot");
433  hAxisNot->SetTitle("SM without TRD");
434  hAxisNot->Draw("");
435 
436  for(Int_t iprod = 0; iprod < nProd; iprod++)
437  {
438  if(!file[iprod]) continue;
439  gMassTRDNo[iprod]->Draw("P");
440  }
441  gMassTRDNo[1]->Draw("P");
442 
443  cGraph->cd(6);
444 
445  TLegend *lN = new TLegend(0.2,0.8,0.9,0.9);
446  lN->SetFillColor(0);
447  lN->SetFillStyle(0);
448  lN->SetLineColor(0);
449  lN->SetBorderSize(0);
450  lN->SetTextSize(0.05);
451 
452  //gPad->SetGridx();
453  gPad->SetGridy();
454 
455  TH1F* hAxisNotRat = (TH1F*) hAxisAllRat->Clone("hAxisNotRat");
456  hAxisNotRat->SetTitle("MC/Data ratio, SM without TRD");
457  hAxisNotRat->Draw("");
458 
459  for(Int_t iprod = 1; iprod < nProd; iprod++)
460  {
461  if(!file[iprod]) continue;
462 
463  gMassRatTRDNo[iprod]->Draw("P");
464  gMassRatTRDNo[iprod]->Fit("pol0","R","",2,8);
465  TF1* func = gMassRatTRDNo[iprod]->GetFunction("pol0");
466  lN->AddEntry(func,Form("Fit %2.3f #pm %2.3f",func->GetParameter(0),func->GetParError(0)),"L");
467 
468  }
469  gMassRatTRDNo[1]->Draw("P");
470  lN->Draw();
471 
472  }
473 
474  if ( gMassTRDYe[0] )
475  {
476  cGraph->cd(3);
477 
478  //gPad->SetGridx();
479  gPad->SetGridy();
480 
481  TLegend *lY = new TLegend(0.2,0.8,0.9,0.9);
482  lY->SetFillColor(0);
483  lY->SetFillStyle(0);
484  lY->SetLineColor(0);
485  lY->SetBorderSize(0);
486  lY->SetTextSize(0.05);
487 
488  TH1F* hAxisYes = (TH1F*) hAxisAll->Clone("hAxisYes");
489  hAxisYes->SetTitle("SM with TRD");
490  hAxisYes->Draw("");
491 
492  for(Int_t iprod = 0; iprod < nProd; iprod++)
493  {
494  if(!file[iprod]) continue;
495  gMassTRDYe[iprod]->Draw("P");
496  }
497 
498  gMassTRDYe[1]->Draw("P");
499 
500  cGraph->cd(7);
501 
502  //gPad->SetGridx();
503  gPad->SetGridy();
504 
505  TH1F* hAxisYesRat = (TH1F*) hAxisAllRat->Clone("hAxisYesRat");
506  hAxisYesRat->SetTitle("MC/Data ratio, SM with TRD");
507  hAxisYesRat->Draw("");
508 
509  for(Int_t iprod = 1; iprod < nProd; iprod++)
510  {
511  if(!file[iprod]) continue;
512 
513  gMassRatTRDYe[iprod]->Draw("P");
514  gMassRatTRDYe[iprod]->Fit("pol0","R","",2,8);
515  TF1* func = gMassRatTRDYe[iprod]->GetFunction("pol0");
516  lY->AddEntry(func,Form("Fit %2.3f #pm %2.3f",func->GetParameter(0),func->GetParError(0)),"L");
517  }
518  gMassRatTRDYe[1]->Draw("P");
519  lY->Draw();
520  }
521 
522  if ( gMassTRDYe[0] ) cGraph->cd(4);
523  else cGraph->cd(2);
524 
525  lE->Draw();
526 
527  //cGraph->cd(ncol*nrow);
528  //lE->Draw();
529 
530  fileName = Form("IMfigures/Comparison_Mass_%s_%s_%s_%s",
531  titleData.Data(),titleMC.Data(),calorimeter.Data(),particle.Data());
532 
533  fileName+=fileFormat;
534  cGraph->Print(fileName);
535 
536 }
537 
538 
void CompareInvMassGraphs(TString particle="Pi0", TString calorimeter="EMCAL", TString titleMC="simu_pp_7TeV_MB", TString titleData="LHC11cd_INT7")
File format: ".eps" ".pdf" ...
TString fileName
Plotting utilities.
const TString calorimeter
Definition: anaM.C:36
static TGraphErrors * DivideGraphs(TGraphErrors *gNum, TGraphErrors *gDen)
Definition: PlotUtils.C:151
int Int_t
Definition: External.C:63
const Int_t nProd
float Float_t
Definition: External.C:68
TString fileFormat
TString prodLeg[]
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)
TFile * file
TList with histograms for a given trigger.
TString prod[]
productions to be compared, directory name