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