AliPhysics  defa9f6 (defa9f6)
InvMassFit.C
Go to the documentation of this file.
1 
17 #if !defined(__CINT__) || defined(__MAKECINT__)
18 
19 #include <TString.h>
20 #include <TH2F.h>
21 #include <TH1F.h>
22 #include <TH3F.h>
23 #include <TH1D.h>
24 #include <TF1.h>
25 #include <TMath.h>
26 #include <TCanvas.h>
27 #include <TStyle.h>
28 #include <TPad.h>
29 #include <TFile.h>
30 #include <TLegend.h>
31 #include <TObject.h>
32 #include <TDirectoryFile.h>
33 #include <TGraphErrors.h>
34 #include <TList.h>
35 #include <TArrayD.h>
36 #include <TGaxis.h>
37 #include <TROOT.h>
38 #include "PlotUtils.C"
39 
40 #endif
41 
42 // Settings
43 Bool_t kMix = kFALSE;
44 Int_t kPolN = 1;
45 Bool_t kSumw2 = kTRUE;
46 Float_t kNPairCut = 20;
47 Float_t kFirstTRDSM = -1;
48 TString kHistoStartName = "AnaPi0";
49 TString kProdName = "LHC18c3_NystromOn";
50 TString kFileName = "AnalysisResults";
51 TString kPlotFormat = "eps";
52 TString kCalorimeter= "EMCAL";
53 TString kParticle = "Pi0";
54 Bool_t kTrunMixFunc= kTRUE;
55 Float_t kChi2NDFMax = 1000;
56 Int_t kRebin = 1;
57 
58 // Initializations
59 Double_t nEvt = 0;
60 TFile * fil = 0;
61 TFile * fout = 0;
62 TList * lis = 0;
63 TF1 * fitfun= 0;
64 TDirectoryFile * direc =0;
65 // Mixing
66 TF1 *tPolPi0 = 0;
67 TF1 *tPolEta = 0;
68 
69 Int_t modColorIndex[]={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};
70 Int_t modStyleIndex[]={24,25,25,24,25,24,25,24,25,24,25,21,21,21,21,21,22,26,22,26,22,26,22,26};
71 
75 //-----------------------------------------------------------------------------
76 Bool_t GetFileAndEvents( TString dirName , TString listName )
77 {
78  fil = new TFile(Form("%s/%s.root",kProdName.Data(),kFileName.Data()),"read");
79 
80  printf("Open Input File: %s/%s.root\n",kProdName.Data(),kFileName.Data());
81 
82  if ( !fil ) return kFALSE;
83 
84  direc = (TDirectoryFile*) fil->Get(dirName);
85 
86  if ( !direc && dirName != "" ) return kFALSE;
87 
88  //printf("dir %p, list %s\n",dir,listName.Data());
89 
90  if(direc)
91  lis = (TList*) direc ->Get(listName);
92  else
93  lis = (TList*) fil->Get(listName);
94 
95  if ( !lis && listName != "") return kFALSE;
96 
97  if(!lis)
98  nEvt = ((TH1F*) fil->Get("hNEvents"))->GetBinContent(1);
99  else
100  nEvt = ((TH1F*) lis->FindObject("hNEvents"))->GetBinContent(1);
101 
102  printf("nEvt = %2.3e\n",nEvt);
103  //nEvt = 1;
104 
105  return kTRUE;
106 }
107 
112 //-----------------------------------------------------------------------------
113 void SetFitFun()
114 {
115  //Fitting function
116  //if(mix) kPolN = 0;
117 
118  if(kParticle == "Pi0")
119  {
120  if ( kPolN == 0)
121  fitfun = new TF1("fitfun",FunctionGaussPol0,0.100,0.250,4);
122  else if (kPolN == 1)
123  fitfun = new TF1("fitfun",FunctionGaussPol1,0.100,0.250,5);
124  else if (kPolN == 2)
125  fitfun = new TF1("fitfun",FunctionGaussPol2,0.100,0.250,6);
126  else if (kPolN == 3)
127  fitfun = new TF1("fitfun",FunctionGaussPol3,0.100,0.250,7);
128  else
129  {
130  printf("*** <<< Set Crystal Ball!!! >>> ***\n");
131  fitfun = new TF1("fitfun",FunctionCrystalBall,0.100,0.250,6);
132  }
133 
134  if(kPolN < 4)
135  {
136  fitfun->SetParLimits(0, kNPairCut/5,kNPairCut*1.e4);
137  fitfun->SetParLimits(1, 0.105,0.185);
138  fitfun->SetParLimits(2, 0.001,0.040);
139  }
140  else
141  {
142  fitfun->SetParLimits(0, kNPairCut/5,kNPairCut*1.e6);
143  fitfun->SetParLimits(2, 0.105,0.185);
144  fitfun->SetParLimits(1, 0.001,0.040);
145  fitfun->SetParLimits(3, 0.001,0.040);
146  fitfun->SetParLimits(4, 0,10);
147  fitfun->SetParLimits(5, 1,1.e+6);
148  }
149 
150  } // Pi0
151  else // Eta
152  {
153  if ( kPolN == 0)
154  fitfun = new TF1("fitfun",FunctionGaussPol0,0.400,0.650,4);
155  else if (kPolN == 1)
156  fitfun = new TF1("fitfun",FunctionGaussPol1,0.400,0.650,5);
157  else if (kPolN == 2)
158  fitfun = new TF1("fitfun",FunctionGaussPol2,0.400,0.650,6);
159  else if (kPolN == 3)
160  fitfun = new TF1("fitfun",FunctionGaussPol3,0.400,0.650,7);
161  if(kPolN < 4)
162  {
163  fitfun->SetParLimits(0, kNPairCut/10,1.e+6);
164  fitfun->SetParLimits(1, 0.20,0.80);
165  fitfun->SetParLimits(2, 0.001,0.06);
166  }
167  }
168 
169  fitfun->SetLineColor(kRed);
170  fitfun->SetLineWidth(2);
171 
172  fitfun->SetParName(0,"A");
173  fitfun->SetParName(1,"m_{0}");
174  fitfun->SetParName(2,"#sigma");
175  // fitfun->SetParName(3,"a_{0}");
176  // fitfun->SetParName(4,"a_{1}");
177 
178  if (kPolN > 1)
179  {
180  fitfun->SetParName(5,"a_{2}");
181  if (kPolN > 2) fitfun->SetParName(6,"a_{3}");
182  }
183 
184  //Mix fit func
185  tPolPi0 = new TF1("truncatedPolPi0" , FunctionTruncatedPolPi0 , 0.02, 0.25, 2);
186  tPolEta = new TF1("truncatedPolEta" , FunctionTruncatedPolEta , 0.2, 1, 2);
187 }
188 
189 //-----------------------
191 //-----------------------
192 void Efficiency
193 (Int_t nPt,
194  TArrayD xPt , TArrayD exPt ,
195  TArrayD mesonPt , TArrayD emesonPt ,
196  TArrayD mesonPtBC, TArrayD emesonPtBC,
197  TString hname)
198 {
199  gStyle->SetOptTitle(1);
200  gStyle->SetTitleOffset(2,"Y");
201  gStyle->SetOptStat(0);
202  gStyle->SetOptFit(000000);
203 
204  TGraphErrors * gPrim = (TGraphErrors*) fout->Get("Primary");
205  TGraphErrors * gPrimAcc = (TGraphErrors*) fout->Get("PrimaryInAcceptance");
206  if ( !gPrim || !gPrimAcc ) return ;
207 
208  TArrayD effPt ; effPt .Set(nPt);
209  TArrayD effBCPt ; effBCPt .Set(nPt);
210  TArrayD effPtAcc ; effPtAcc .Set(nPt);
211  TArrayD effBCPtAcc ; effBCPtAcc .Set(nPt);
212  TArrayD effPtErr ; effPtErr .Set(nPt);
213  TArrayD effBCPtErr ; effBCPtErr .Set(nPt);
214  TArrayD effPtAccErr ; effPtAccErr .Set(nPt);
215  TArrayD effBCPtAccErr; effBCPtAccErr.Set(nPt);
216 
217  for(Int_t ibin = 0; ibin < nPt; ibin++ )
218  {
219  //printf("Bin %d \n",ibin);
220  //printf("- Fit Reco %2.3e, Prim %2.3e, PrimAcc %2.3e\n",
221  // mesonPt[ibin],gPrim->GetY()[ibin],gPrimAcc->GetY()[ibin]);
222 
223  if ( gPrim->GetY()[ibin] > 0 && mesonPt[ibin] > 0 )
224  {
225  effPt [ibin] = 100 * mesonPt[ibin] / gPrim->GetY()[ibin];
226  effPtErr[ibin] = 100 * GetFractionError( mesonPt[ibin], gPrim->GetY ()[ibin],
227  emesonPt[ibin], gPrim->GetEY()[ibin]); // PlotUtils.C
228 
229  //printf("\t EffxAcc %f, err %f\n",effPt[ibin],effPtErr[ibin]);
230  }
231  else { effPt[ibin] = 0; effPtErr[ibin] = 0; }
232 
233  if ( gPrimAcc->GetY()[ibin] > 0 && mesonPt[ibin] > 0 )
234  {
235  effPtAcc [ibin] = 100 * mesonPt[ibin] / gPrimAcc->GetY()[ibin];
236  effPtAccErr[ibin] = 100 * GetFractionError( mesonPt[ibin], gPrimAcc->GetY ()[ibin],
237  emesonPt[ibin], gPrimAcc->GetEY()[ibin]); // PlotUtils.C
238  //printf("\t Eff %f, err %f\n",effPtAcc[ibin],effPtAccErr[ibin]);
239  }
240  else { effPtAcc[ibin] = 0; effPtAccErr[ibin] = 0; }
241 
242  if ( kMix || hname.Contains("MCpT") )
243  {
244  //printf("- BC Reco %2.3e, Prim %2.3e, PrimAcc %2.3e\n",
245  // mesonPtBC[ibin],gPrim->GetY()[ibin],gPrimAcc->GetY()[ibin]);
246 
247  if ( gPrim->GetY()[ibin] > 0 && mesonPtBC[ibin] > 0 )
248  {
249  effBCPt [ibin] = 100 * mesonPtBC[ibin] / gPrim->GetY()[ibin];
250  effBCPtErr[ibin] = 100 * GetFractionError( mesonPtBC[ibin], gPrim->GetY ()[ibin],
251  emesonPtBC[ibin], gPrim->GetEY()[ibin]); // PlotUtils.C
252  //printf("\t EffxAcc %f, err %f\n",effBCPt[ibin],effBCPtErr[ibin]);
253  }
254  else { effBCPt[ibin] = 0; effBCPtErr[ibin] = 0; }
255 
256  if ( gPrimAcc->GetY()[ibin] > 0 && mesonPtBC[ibin] > 0 )
257  {
258  effBCPtAcc [ibin] = 100 * mesonPtBC[ibin] / gPrimAcc->GetY()[ibin];
259  effBCPtAccErr[ibin] = 100 * GetFractionError( mesonPtBC[ibin], gPrimAcc->GetY ()[ibin],
260  emesonPtBC[ibin], gPrimAcc->GetEY()[ibin]); // PlotUtils.C
261  //printf("\t EffxAcc %f, err %f\n",effBCPtAcc[ibin],effBCPtAccErr[ibin]);
262  }
263  else { effBCPtAcc[ibin] = 0; effBCPtAccErr[ibin] = 0; }
264 
265  } // Bin counting
266 
267  } // pt bin loop
268 
269  TGraphErrors *gEff = new TGraphErrors(nPt,xPt.GetArray(),effPt.GetArray(),exPt.GetArray(),effPtErr.GetArray());
270  gEff->SetName(Form("EfficiencyxAcceptance_%s",hname.Data()));
271  gEff->GetHistogram()->SetYTitle("#epsilon_{Reco #times PID #times Acc} (%)");
272  gEff->GetHistogram()->SetTitleOffset(1.4,"Y");
273  gEff->GetHistogram()->SetXTitle("#it{p}_{T} (GeV/#it{c})");
274  gEff->GetHistogram()->SetTitleOffset(1.2,"X");
275  gEff->GetHistogram()->SetTitle("Reconstruction efficiency #times acceptance");
276  gEff->SetMarkerColor(1);
277  gEff->SetLineColor(1);
278  gEff->SetMarkerStyle(20);
279  gEff->Write();
280 
281  TGraphErrors *gEffAcc = new TGraphErrors(nPt,xPt.GetArray(),effPtAcc.GetArray(),exPt.GetArray(),effPtAccErr.GetArray());
282  gEffAcc->SetName(Form("Efficiency_%s",hname.Data()));
283  gEffAcc->GetHistogram()->SetYTitle("#epsilon_{Reco #times PID} (%)");
284  gEffAcc->GetHistogram()->SetTitleOffset(1.5,"Y");
285  gEffAcc->GetHistogram()->SetXTitle("#it{p}_{T} (GeV/#it{c})");
286  gEffAcc->GetHistogram()->SetTitleOffset(1.2,"X");
287  gEffAcc->SetMarkerColor(1);
288  gEffAcc->SetLineColor(1);
289  gEffAcc->SetMarkerStyle(20);
290  gEffAcc->GetHistogram()->SetTitle("Reconstruction efficiency");
291 
292  gEffAcc->Write();
293 
294  TGraphErrors *gBCEff = 0;
295  TGraphErrors *gBCEffAcc = 0;
296  if(kMix)
297  {
298  gBCEff = new TGraphErrors(nPt,xPt.GetArray(),effBCPt.GetArray(),exPt.GetArray(),effBCPtErr.GetArray());
299  gBCEff->SetName(Form("EfficiencyxAcceptance_BC_%s",hname.Data()));
300  gBCEff->GetHistogram()->SetYTitle("#epsilon_{Reco #times PID #times Acc}");
301  gBCEff->GetHistogram()->SetTitleOffset(1.5,"Y");
302  gBCEff->GetHistogram()->SetXTitle("#it{p}_{T} (GeV/#it{c})");
303  gBCEff->GetHistogram()->SetTitleOffset(1.2,"X");
304  gBCEff->SetMarkerColor(4);
305  gBCEff->SetLineColor(4);
306  gBCEff->SetMarkerStyle(24);
307  gBCEff->Write();
308 
309  gBCEffAcc = new TGraphErrors(nPt,xPt.GetArray(),effBCPtAcc.GetArray(),exPt.GetArray(),effBCPtAccErr.GetArray());
310  gBCEffAcc->SetName(Form("Efficiency_BC_%s",hname.Data()));
311  gBCEffAcc->GetHistogram()->SetYTitle("#epsilon_{Reco #times PID}");
312  gBCEffAcc->GetHistogram()->SetTitleOffset(1.4,"Y");
313  gBCEffAcc->GetHistogram()->SetXTitle("#it{p}_{T} (GeV/#it{c})");
314  gBCEffAcc->GetHistogram()->SetTitleOffset(1.2,"X");
315  gBCEffAcc->SetMarkerColor(4);
316  gBCEffAcc->SetLineColor(4);
317  gBCEffAcc->SetMarkerStyle(24);
318  gBCEffAcc->Write();
319  }
320 
321  // Plot efficiencies
322  //
323  TCanvas *cEff = new TCanvas(Form("cEff_%s",hname.Data()),
324  Form("Efficiency Graphs for %s",hname.Data()),2*600,600);
325  cEff->Divide(2, 1);
326 
327  cEff->cd(1);
328  gPad->SetGridx();
329  gPad->SetGridy();
330  //gPad->SetLogy();
331 
332  //gEff->SetMaximum(1);
333  //gEff->SetMinimum(1e-8);
334 
335  gEff->Draw("AP");
336 
337  if(kMix)
338  {
339  gBCEff ->Draw("P");
340 
341  TLegend * legend = new TLegend(0.5,0.7,0.9,0.9);
342  legend->AddEntry(gEff ,"From fit","P");
343  legend->AddEntry(gBCEff,"From bin counting","P");
344  legend->Draw();
345  }
346 
347  cEff->cd(2);
348  gPad->SetGridx();
349  gPad->SetGridy();
350  //gPad->SetLogy();
351 
352  //gEffAcc->SetMaximum(1);
353  //gEffAcc->SetMinimum(1e-8);
354 
355  gEffAcc->Draw("AP");
356 
357  if(kMix)
358  {
359  gBCEffAcc->Draw("P");
360 
361  TLegend * legend = new TLegend(0.5,0.7,0.9,0.9);
362  legend->AddEntry(gEffAcc ,"From fit","P");
363  legend->AddEntry(gBCEffAcc,"From bin counting","P");
364  legend->Draw();
365  }
366 
367  cEff->Print(Form("IMfigures/%s_%s_Efficiency_%s_%s.%s",
368  kProdName.Data(),kCalorimeter.Data(),kParticle.Data(),hname.Data(),
369  /*kFileName.Data(),*/ kPlotFormat.Data()));
370 }
371 
372 //-------------------------------------
375 //-------------------------------------
376 void ProjectAndFit
377 (
378  const Int_t nPt, TString comment, TString leg, TString hname,
379  TArrayD xPtLimits, TArrayD xPt, TArrayD exPt,
380  TH2F * hRe, TH2F * hMi
381  )
382 {
383  if ( !hRe )
384  {
385  printf("ProjectAndFit() - Null inv mass 2D histo for %s, skip!\n",hname.Data());
386  return;
387  }
388 
389  gStyle->SetOptTitle(1);
390  gStyle->SetTitleOffset(2,"Y");
391  gStyle->SetOptStat(0);
392  gStyle->SetOptFit(000000);
393 
394  Double_t mmin = 0;
395  Double_t mmax = 1;
396  TString particleN = "";
397 
398  if(kParticle=="Pi0")
399  {
400  mmin = 0.04;
401  mmax = 0.44;
402  particleN = " #pi^{0}";
403  }
404  else // eta
405  {
406  mmin = 0.25;
407  mmax = 0.95;
408  particleN = " #eta";
409  }
410 
411  TLegend * pLegendIM[nPt];
412  TH1D* hIM [nPt];
413  TH1D* hMix [nPt];
414  TH1D* hMixCorrected [nPt];
415  TH1D* hRatio [nPt];
416  TH1D* hSignal [nPt];
417  for(Int_t ipt = 0; ipt< nPt; ipt++)
418  {
419  hIM [ipt] = 0 ;
420  hMix [ipt] = 0 ;
421  hMixCorrected[ipt] = 0 ;
422  hRatio [ipt] = 0 ;
423  hSignal [ipt] = 0 ;
424  }
425 
426  TArrayD mesonPt ; mesonPt .Set(nPt);
427  TArrayD mesonPtBC ; mesonPtBC .Set(nPt);
428  TArrayD mesonMass ; mesonMass .Set(nPt);
429  TArrayD mesonWidth; mesonWidth.Set(nPt);
430  TArrayD emesonPt ; emesonPt .Set(nPt);
431  TArrayD emesonPtBC ; emesonPtBC .Set(nPt);
432  TArrayD emesonMass ; emesonMass .Set(nPt);
433  TArrayD emesonWidth; emesonWidth.Set(nPt);
434 
435  Int_t rebin2 = 2*kRebin;
436 
437  Int_t col = TMath::Ceil(TMath::Sqrt(nPt));
438 
439  //hRe->Scale(1./nEvt);
440 
441  // Mix correction factor
442  TF1 *line3[nPt];// = new TF1("line3", "[0]+[1]*x+[2]*x*x+[3]*x*x*x", 0.0, 1);
443 
444  TF1 * fitFunction = 0;
445  Bool_t ok = kFALSE;
446 
447  TCanvas * cIMModi = new TCanvas(Form("c%s", hname.Data()), Form("%s", leg.Data()), 1200, 1200) ;
448  cIMModi->Divide(col, col);
449 
450  for(Int_t i = 0; i < nPt; i++)
451  {
452  cIMModi->cd(i+1) ;
453  //gPad->SetLogy();
454 
455  Double_t ptMin = xPtLimits[i];
456  Double_t ptMax = xPtLimits[i+1];
457  //printf("Bin %d (%2.1f, %2.1f)\n",i, ptMin,ptMax);
458 
459  hIM[i] = hRe->ProjectionY(Form("IM_%s_PtBin%d",hname.Data(),i),
460  hRe->GetXaxis()->FindBin(ptMin),
461  hRe->GetXaxis()->FindBin(ptMax));
462  hIM[i]->SetTitle(Form("%2.1f < #it{p}_{T, #gamma#gamma} < %2.1f GeV/#it{c}",ptMin,ptMax));
463  if(i < nPt-3)
464  hIM[i]->Rebin(kRebin);
465  else
466  hIM[i]->Rebin(rebin2);
467 
468  if ( kSumw2 )
469  hIM[i]->Sumw2();
470 
471  hIM[i]->SetXTitle("M_{#gamma,#gamma} (GeV/#it{c}^{2})");
472  //printf("mmin %f, mmax %f\n",mmin,mmax);
473  hIM[i]->SetAxisRange(mmin,mmax,"X");
474  hIM[i]->SetLineWidth(2);
475  hIM[i]->SetLineColor(4);
476 
477  Double_t mStep = hIM[i]->GetBinWidth(1);
478 
479  hSignal[i] = (TH1D*) hIM[i]->Clone();
480 
481  //--------------------------------------------------
482  // Mix: Project, scale and subract to real pairs
483  //--------------------------------------------------
484  if ( kMix && hMi )
485  {
486  hMix[i] = (TH1D*) hMi->ProjectionY(Form("MiMass_PtBin%d_%s",i,hname.Data()),
487  hMi->GetXaxis()->FindBin(ptMin),
488  hMi->GetXaxis()->FindBin(ptMax));
489  hMix[i]->SetLineColor(1);
490  hMix[i]->SetTitle(Form("%2.2f < #it{p}_{T, #gamma#gamma} < %2.2f GeV/#it{c}",ptMin,ptMax));
491  if(i < nPt-3)
492  hMix[i]->Rebin(kRebin);
493  else
494  hMix[i]->Rebin(rebin2);
495 
496  // Double_t sptmin = 0.2;
497  // Double_t sptmax = 0.45;
498 
499  // Double_t scalereal = hIM->Integral(hIM->FindBin(sptmin),
500  // hIM->FindBin(sptmax));
501  // Double_t scalemix = hMix->Integral(hMix->FindBin(sptmin),
502  // hMix->FindBin(sptmax));
503  // if(scalemix > 0)
504  // hMix->Scale(scalereal/scalemix);
505  // else{
506  // //printf("could not scale: re %f mi %f\n",scalereal,scalemix);
507  // continue;
508  // }
509 
510  // hMix[i]->SetLineWidth(1);
511 
512  if ( kSumw2 )
513  hMix[i]->Sumw2();
514 
515  // --- Ratio ---
516  hRatio[i] = (TH1D*)hIM[i]->Clone(Form("RatioRealMix_PtBin%d_%s",i,hname.Data()));
517  hRatio[i]->SetAxisRange(mmin,mmax,"X");
518  hRatio[i]->Divide(hMix[i]);
519  Double_t p0 =0;
520  Double_t p1 =0;
521  Double_t p2 =0;
522  Double_t p3 =0;
523 
524  // --- Subtract from signal ---
525  hSignal[i] ->SetName(Form("Signal_PtBin%d_%s",i,hname.Data()));
526 
527  if ( hMix[i]->GetEntries() > kNPairCut*3 )
528  {
529  hSignal[i] ->SetLineColor(kViolet);
530 
531  if ( kTrunMixFunc )
532  {
533  if ( kParticle == "Pi0" )
534  {
535  hRatio[i]->Fit("truncatedPolPi0", "NQRL","",0.11,0.4);
536  p0= tPolPi0->GetParameter(0);
537  p1= tPolPi0->GetParameter(1);
538  //p2= tPolPi0->GetParameter(2);
539  //p3= tPolPi0->GetParameter(3);
540  } // pi0
541  else // eta
542  {
543  hRatio[i]->SetAxisRange(mmin,mmax);
544  hRatio[i]->Fit("truncatedPolEta", "NQRL","",0.35,0.95);
545  p0 = tPolEta->GetParameter(0);
546  p1 = tPolEta->GetParameter(1);
547  //p2 = tPolEta->GetParameter(2);
548  //p3 = tPolEta->GetParameter(3);
549  }
550 
551  //line3[i] = new TF1("line3", "[0]+[1]*x+[2]*x*x+[3]*x*x*x", mmin, mmax);
552  line3[i] = new TF1("line3", "[0]+[1]*x", mmin, mmax);
553  line3[i]->SetParameters(p0,p1,p2,p3);
554 
555  //printf("Correct\n");
556  hMixCorrected[i] = (TH1D*) hMix[i]->Clone(Form("MixCorrected_PtBin%d_%s",i,hname.Data()));
557  for(Int_t j = 0; j< hMix[i]->GetNbinsX(); j++)
558  {
559  Double_t x = hMix[i]->GetBinCenter(j);
560  Double_t correction = line3[i]->Eval(x);
561  Double_t corrected = hMix[i]->GetBinContent(j)*correction;
562  Double_t ecorrected = hMix[i]->GetBinError(j) *correction;
563  //printf("bin %d, center %f, mix %f, corrected %f\n",i,x,hMix[i]->GetBinContent(j),corrected);
564  hMixCorrected[i]->SetBinContent(j, corrected);
565  hMixCorrected[i]->SetBinError (j,ecorrected);
566  }
567 
568  // Subtract
569  hSignal[i] ->Add(hMixCorrected[i],-1);
570  }
571  else
572  {
573  Float_t scale = hRatio[i]->GetBinContent(hRatio[i]->FindBin(0.7));
574  hMix[i]->Scale(scale);
575  // Subtract
576  hSignal[i] ->Add(hMix[i],-1);
577  }
578  } // enough mixed entries
579  } // mixed event histogtam exists
580 
581  //----------------------------
582  // ---- Fit subtracted signal
583  //----------------------------
584  Double_t nMax = 0;
585  if(kParticle=="Pi0")
586  {
587  nMax= hSignal[i]->Integral(hIM[i]->FindBin(0.1),
588  hIM[i]->FindBin(0.2));
589  }
590  else
591  {
592  nMax = hSignal[i]->Integral(hIM[i]->FindBin(0.5),
593  hIM[i]->FindBin(0.7));
594  }
595 
596  if(nMax > kNPairCut)
597  {
598  fitfun->SetParLimits(0,nMax/100,nMax*100);
599 
600  if(kParticle=="Pi0")
601  {
602  if(kPolN < 4 )fitfun->SetParameters(nMax/5,0.135,20,0);
603  else fitfun->SetParameters(nMax/5,20,0.135,20,20,nMax/5);
604 
605  if(i < nPt-4)
606  hSignal[i]->Fit("fitfun","QR","",0.11,0.3);
607  else
608  hSignal[i]->Fit("fitfun","QR","",0.11,0.3);
609  }// pi0
610  else // eta
611  {
612 
613  if ( kPolN < 4 ) fitfun->SetParameters(nMax/5,0.54,40,0);
614  else fitfun->SetParameters(nMax/5,40,0.54,40,40,nMax/5);
615  //if(i<4)
616  hSignal[i]->Fit("fitfun","QR","",0.4,0.7);
617  //else if(i < 15)
618  // hSignal[i]->Fit("fitfun","QRL","",0.3,0.8);
619  //else
620  //hSignal[i]->Fit("fitfun","QRL","",0.3,0.8);
621  }
622  }
623  else
624  printf("Skip bin %d: n max %2.3e < cut %2.3e\n",i, nMax, kNPairCut);
625 
626 
627  //----------------------------
628  // ---- Put fit results in arrays
629  //----------------------------
630  // Init results arrays
631  mesonPt .SetAt(-1,i);
632  emesonPt .SetAt(-1,i);
633  mesonPtBC .SetAt(-1,i);
634  emesonPtBC .SetAt(-1,i);
635  mesonMass .SetAt(-1,i);
636  emesonMass .SetAt(-1,i);
637  mesonWidth.SetAt(-1,i);
638  emesonWidth.SetAt(-1,i);
639 
640  fitFunction = (TF1*) hSignal[i]->GetFunction("fitfun");
641 
642  Float_t chi2ndf = 1e6;
643  if ( fitFunction )
644  {
645  // printf("ipt %d: Chi/NDF %f, Chi %f, NDF %d\n",i,
646  // fitFunction->GetChisquare()/fitFunction->GetNDF(),
647  // fitFunction->GetChisquare(),
648  // fitFunction->GetNDF());
649 
650  Float_t chi2 = fitFunction->GetChisquare();
651  Int_t ndf = fitFunction->GetNDF();
652  if( ndf > 0 ) chi2ndf = chi2 / ndf;
653 
654  if ( chi2ndf < kChi2NDFMax )
655  {
656  Double_t A = fitFunction->GetParameter(0);
657  Double_t mean = fitFunction->GetParameter(1);
658  Double_t sigm = fitFunction->GetParameter(2);
659  // Double_t a0 = fitFunction->GetParameter(3);
660  // Double_t a1 = fitFunction->GetParameter(4);
661  // Double_t a2 = 0;
662  // if (kPolN == 2)
663  // a2 = fitFunction->GetParameter(5);
664 
665  Double_t eA = fitFunction->GetParError(0);
666  Double_t emean = fitFunction->GetParError(1);
667  Double_t esigm = fitFunction->GetParError(2);
668  // Double_t ea0 = fitFunction->GetParError(3);
669  // Double_t ea1 = fitFunction->GetParError(4);
670  // Double_t ea2 = 0;
671  // if (kPolN == 2)
672  // ea2 = fitFunction->GetParError(5);
673  pLegendIM[i] = new TLegend(0.48,0.65,0.95,0.93);
674  pLegendIM[i]->SetTextSize(0.035);
675  pLegendIM[i]->SetFillColor(10);
676  pLegendIM[i]->SetBorderSize(1);
677  pLegendIM[i]->SetHeader(Form(" %s - %s",particleN.Data(), leg.Data()));
678  pLegendIM[i]->AddEntry("",Form("A = %2.1e#pm%2.1e ",A,eA),"");
679  pLegendIM[i]->AddEntry("",Form("#mu = %3.1f #pm %3.1f GeV/#it{c}^{2}",mean*1000,emean*1000),"");
680  pLegendIM[i]->AddEntry("",Form("#sigma = %3.1f #pm %3.1f GeV/#it{c}^{2}",sigm*1000,esigm*1000),"");
681 
682  //pLegendIM[i]->AddEntry("",Form("p_{0} = %2.1f #pm %2.1f ",a0,ea0),"");
683  //pLegendIM[i]->AddEntry("",Form("p_{1} = %2.1f #pm %2.1f ",a1,ea1),"");
684  //if(kPolN==2)pLegendIM[i]->AddEntry("",Form("p_{2} = %2.1f#pm %2.1f ",a2,ea2),"");
685 
686  Double_t counts = A*sigm / mStep * TMath::Sqrt(TMath::TwoPi());
687  Double_t eCounts =
688  TMath::Power(eA/A,2) +
689  TMath::Power(esigm/sigm,2);
690  eCounts = TMath::Sqrt(eCounts) * counts;
691  //eCounts = TMath::Min(eCounts, TMath::Sqrt(counts));
692 
693  counts/=(nEvt*(exPt.At(i)*2));
694  eCounts/=(nEvt*(exPt.At(i)*2));
695 
696  mesonPt.SetAt( counts,i);
697  emesonPt.SetAt(eCounts,i);
698 
699  //printf("A %e Aerr %e; mu %f muErr %f; sig %f sigErr %f",
700  // A,eA,mean,emean,sigm,esigm);
701  //cout << "N(pi0) fit = " << counts << "+-"<< eCounts << " mstep " << mStep<<endl;
702 
703  mesonMass .SetAt( mean*1000.,i);
704  emesonMass .SetAt(emean*1000.,i);
705 
706  mesonWidth.SetAt( sigm*1000.,i);
707  emesonWidth.SetAt(esigm*1000.,i);
708  } //Good fit
709  else
710  {
711  printf("Bin %d, Bad fit, Chi2 %f ndf %d, ratio %f!\n",
712  i,chi2,ndf,chi2ndf);
713  }
714 
715  }
716  else printf("Bin %d, NO fit available!\n",i);
717 
718  // Set the integration window, depending on mass mean and width
719  // 2 sigma
720  Double_t mass = mesonMass [i]/1000.;
721  Double_t width = mesonWidth[i]/1000.;
722  Double_t mMinBin = mass - 2 * width;
723  Double_t mMaxBin = mass + 2 * width;
724  if(mass <= 0 || width <= 0)
725  {
726  if(kParticle=="Pi0")
727  {
728  mMinBin = 0.115; mMaxBin = 0.3 ;
729  }
730  else // eta
731  {
732  mMinBin = 0.4; mMaxBin = 0.9 ;
733  }
734  }
735 
736  //
737  // Bin counting instead of fit
738  //
739  //printf("Signal %p, Mix %p, hname %s\n",hSignal[i], hMix[i], hname.Data());
740  if ( hSignal[i] && ( hMix[i] || hname.Contains("MCpT") ) )
741  {
742  Double_t countsBin = 0.;//hSignal[i]->Integral(hSignal[i]->FindBin(mMinBin), hSignal[i]->FindBin(mMaxBin)) ;
743  Double_t eCountsBin = 0.;//TMath::Sqrt(countsBin);
744 
745  // PlotUtils.C
746  GetRangeIntegralAndError(hSignal[i], hSignal[i]->FindBin(mMinBin), hSignal[i]->FindBin(mMaxBin), countsBin, eCountsBin );
747 
748  countsBin/=(nEvt*(exPt.At(i)*2));
749  eCountsBin/=(nEvt*(exPt.At(i)*2));
750 
751  mesonPtBC.SetAt( countsBin,i);
752  emesonPtBC.SetAt(eCountsBin,i);
753 
754  //printf("pt bin %d: [%2.1f,%2.1f] - Mass window [%2.2f, %2.2f] - N(pi0) BC = %2.3e +- %2.4e\n",
755  // i, xPtLimits[i], xPtLimits[i+1], mMinBin, mMaxBin, countsBin, eCountsBin);
756  }
757 
758  if ( nMax > kNPairCut && chi2ndf < kChi2NDFMax )
759  {
760  if( kMix && hMix[i] )
761  {
762  // if ( !kMix ) hIM[i]->SetMinimum( 0.1);
763  // else hIM[i]->SetMinimum(-0.1);
764 
765  hIM[i]->SetMaximum(hIM[i]->GetMaximum()*1.2);
766  hIM[i]->SetMinimum(hSignal[i]->GetMinimum()*0.2);
767 
768  hIM[i]->Draw("HE");
769  pLegendIM[i]->AddEntry(hIM[i],"Raw pairs" ,"L");
770 
771  if ( hMix[i]->GetEntries() > kNPairCut*3 )
772  {
773  if ( kTrunMixFunc )
774  {
775  hMixCorrected[i]->Draw("same");
776  pLegendIM[i]->AddEntry(hMixCorrected[i],"Mixed pairs" ,"L");
777  }
778  else
779  {
780  hMix[i]->Draw("same");
781  pLegendIM[i]->AddEntry(hMix[i],"Mixed pairs" ,"L");
782  }
783  }
784 
785  ok = kTRUE;
786 
787  hSignal[i] -> Draw("same");
788  pLegendIM[i]->AddEntry(hSignal[i],"Signal pairs","L");
789  }
790  else
791  {
792  ok = kTRUE;
793  hSignal[i]->SetMaximum(hSignal[i]->GetMaximum()*1.2);
794  hSignal[i]->SetMinimum(hSignal[i]->GetMinimum()*0.8);
795 
796  pLegendIM[i]->AddEntry(hSignal[i],"Raw pairs","L");
797  hSignal[i] -> Draw("HE");
798  }
799 
800  // Plot mass from pairs originated from pi0/eta
801  if ( hname.Contains("All") )
802  {
803  TH1F* hMesonPairOrigin = (TH1F*) fout->Get(Form("IM_MCpTReco_PtBin%d",i));
804  if ( hMesonPairOrigin )
805  {
806  hMesonPairOrigin->SetLineColor(8);
807  hMesonPairOrigin->Draw("same");
808  pLegendIM[i]->AddEntry(hMesonPairOrigin,Form("%s origin",particleN.Data()),"L");
809  }
810  } // MC origin
811 
812 // if ( fitFunction )
813 // pLegendIM[i]->AddEntry(fitFunction,"Fit","L");
814 
815  pLegendIM[i]->Draw();
816  }
817  else
818  {
819  // Plot just raw pairs
820  hIM[i]->Draw("HE");
821  }
822  } //pT bins
823 
824  if ( ok )
825  cIMModi->Print(Form("IMfigures/%s_%s_Mgg_%s_%s.%s",
826  kProdName.Data(),kCalorimeter.Data(),kParticle.Data(),/*kFileName.Data(),*/
827  hname.Data(),kPlotFormat.Data()));
828 
829  //xxxx Real / Mixxxxx
830  if(kMix)
831  {
832  //printf("Do real/mix\n");
833  Bool_t okR = kFALSE;
834  TCanvas * cRat = new TCanvas(Form("Ratio_%s\n",hname.Data()), Form("Ratio %s\n", leg.Data()), 1200, 1200) ;
835  cRat->Divide(col, col);
836 
837  for(Int_t i = 0; i < nPt; i++)
838  {
839  cRat->cd(i+1) ;
840  //gPad->SetGridy();
841  //gPad->SetLog();
842  if(!hRatio[i])
843  {
844  //printf("No ratio in pt bin %d continue\n",i);
845  continue;
846  }
847 
848  Double_t nMax =0;
849  if(kParticle=="Pi0")
850  nMax = hRatio[i]->Integral(hRatio[i]->FindBin(0.005),
851  hRatio[i]->FindBin(0.20));
852  else
853  nMax = hRatio[i]->Integral(hRatio[i]->FindBin(0.4),
854  hRatio[i]->FindBin(0.6));
855 
856  if ( nMax <= 0 ) continue;
857 
858  //printf("Ratio nMax = %e \n",nMax);
859 
860  hRatio[i]->SetMaximum(nMax/4);
861  hRatio[i]->SetMinimum(1e-6);
862  hRatio[i]->SetAxisRange(mmin,mmax,"X");
863 
864  okR = kTRUE;
865 
866  hRatio[i]->SetYTitle("Real pairs / Mixed pairs");
867  hRatio[i]->SetTitleOffset(1.5,"Y");
868 
869  hRatio[i]->Draw();
870 
871  //if(line3[i]) line3[i]->Draw("same");
872  }
873 
874  if ( okR )
875  cRat->Print(Form("IMfigures/%s_%s_MggRatio_%s_%s.%s",
876  kProdName.Data(),kCalorimeter.Data(),kParticle.Data(),/*kFileName.Data(),*/
877  hname.Data(),kPlotFormat.Data()));
878  } // Mix
879 
880  if( !ok ) return ;
881 
882  //------------------------------
883  // Fit parameters
884  //------------------------------
885  // Put fit results in TGraphErrors
886 
887  TGraphErrors* gPt = new TGraphErrors(nPt,xPt.GetArray(),mesonPt.GetArray(),exPt.GetArray(),emesonPt.GetArray());
888  gPt->SetName(Form("gPt_%s",hname.Data()));
889  gPt->GetHistogram()->SetTitle(Form("#it{p}_{T} of reconstructed %s, %s ",particleN.Data(), comment.Data()));
890  gPt->GetHistogram()->SetXTitle("#it{p}_{T} (GeV/#it{c}) ");
891  gPt->GetHistogram()->SetYTitle(Form("dN_{%s}/d#it{p}_{T} (GeV/#it{c})^{-1} / N_{events} ",particleN.Data()));
892  //gPt->GetHistogram()->SetAxisRange(0.,30);
893  gPt->GetHistogram()->SetTitleOffset(1.5,"Y");
894  gPt->SetMarkerStyle(20);
895  //gPt->SetMarkerSize(1);
896  gPt->SetMarkerColor(1);
897  // gPt->GetHistogram()->SetMaximum(1e8);
898  // gPt->GetHistogram()->SetMinimum(1e-8);
899 
900  TGraphErrors* gPtBC = new TGraphErrors(nPt,xPt.GetArray(),mesonPtBC.GetArray(),exPt.GetArray(),emesonPtBC.GetArray());
901  gPtBC->SetName(Form("gPtBC_%s",hname.Data()));
902  gPtBC->GetHistogram()->SetTitle(Form("#it{p}_{T} of reconstructed %s, %s ",particleN.Data(), comment.Data()));
903  gPtBC->GetHistogram()->SetXTitle("#it{p}_{T} (GeV/#it{c}) ");
904  gPtBC->GetHistogram()->SetYTitle(Form("d#it{N}_{%s}/d#it{p}_{T} (GeV/#it{c})^{-1} / #it{N}_{events} ",particleN.Data()));
905  //gPtBC->GetHistogram()->SetAxisRange(0.,30);
906  gPtBC->GetHistogram()->SetTitleOffset(1.5,"Y");
907  gPtBC->SetMarkerStyle(24);
908  //gPtBC->SetMarkerSize(1);
909  gPtBC->SetMarkerColor(4);
910  // gPtBC->GetHistogram()->SetMaximum(1e8);
911  // gPtBC->GetHistogram()->SetMinimum(1e-8);
912 
913  TGraphErrors* gMass = new TGraphErrors(nPt,xPt.GetArray(),mesonMass.GetArray(),exPt.GetArray(),emesonMass.GetArray());
914  gMass->SetName(Form("gMass_%s",hname.Data()));
915  gMass->GetHistogram()->SetTitle(Form("mass vs #it{p}_{T} of reconstructed %s, %s ",particleN.Data(), comment.Data()));
916  gMass->GetHistogram()->SetXTitle("#it{p}_{T} (GeV/#it{c}) ");
917  gMass->GetHistogram()->SetYTitle(Form("mass_{%s} (MeV/#it{c}^{2}) ",particleN.Data()));
918  //gMass->GetHistogram()->SetAxisRange(0.,30);
919  gMass->GetHistogram()->SetTitleOffset(1.5,"Y");
920  gMass->SetMarkerStyle(20);
921  //gMass->SetMarkerSize(1);
922  gMass->SetMarkerColor(1);
923 
924  TGraphErrors* gWidth= new TGraphErrors(nPt,xPt.GetArray(),mesonWidth.GetArray(),exPt.GetArray(),emesonWidth.GetArray());
925  gWidth->SetName(Form("gWidth_%s",hname.Data()));
926  gWidth->GetHistogram()->SetTitle(Form("Width vs #it{p}_{T} of reconstructed %s, %s ",particleN.Data(), comment.Data()));
927  gWidth->GetHistogram()->SetXTitle("#it{p}_{T} (GeV/#it{c}) ");
928  gWidth->GetHistogram()->SetYTitle(Form("#sigma_{%s} (MeV/#it{c}^{2}) ",particleN.Data()));
929  //gWidth->GetHistogram()->SetAxisRange(0.,30);
930  gWidth->GetHistogram()->SetTitleOffset(1.5,"Y");
931  gWidth->SetMarkerStyle(20);
932  //gWidth->SetMarkerSize(1);
933  gWidth->SetMarkerColor(1);
934 
935  // Plot the fitted results
936 
937  TCanvas *cFitGraph = new TCanvas(Form("cFitGraph_%s",hname.Data()),
938  Form("Fit Graphs for %s",hname.Data()),600,600);
939  cFitGraph->Divide(2, 2);
940 
941  // Mass
942  cFitGraph->cd(1);
943  gPad->SetGridx();
944  gPad->SetGridy();
945 
946  if ( kParticle=="Pi0" )
947  {
948  gMass->SetMaximum(160);
949  gMass->SetMinimum(100);
950  }
951  else
952  {
953  gMass->SetMaximum(660);
954  gMass->SetMinimum(420);
955  }
956 
957  gMass->Draw("AP");
958 
959  cFitGraph->cd(2);
960  gPad->SetGridx();
961  gPad->SetGridy();
962 
963  if(kParticle=="Pi0")
964  {
965  gWidth->SetMaximum(16);
966  gWidth->SetMinimum(6);
967  }
968  else
969  {
970  gWidth->SetMaximum(60);
971  gWidth->SetMinimum(0);
972  }
973 
974  gWidth->Draw("AP");
975 
976  cFitGraph->cd(3);
977  gPad->SetGridx();
978  gPad->SetGridy();
979  gPad->SetLogy();
980 
981  gPt->Draw("AP");
982 
983  Double_t maximum = gPt->GetHistogram()->GetMaximum();
984  Double_t minimum = gPt->GetHistogram()->GetMinimum();
985  //printf("A-maximum %e, minimum %e\n",maximum,minimum);
986 
987 // Double_t maximum = gPrim->GetHistogram()->GetMaximum();
988 // Double_t minimum = gPrim->GetHistogram()->GetMinimum();
989 //
990 // if ( gPrimAcc->GetMaximum() > maximum ) maximum = gPrimAcc->GetMaximum() ;
991 // if ( gPrimAcc->GetMinimum() > minimum ) minimum = gPrimAcc->GetMinimum() ;
992 //
993 // gPrim->SetMaximum(maximum*10);
994 // gPrim->SetMinimum(minimum/10);
995 
996  if(kMix)
997  {
998  gPtBC ->Draw("P");
999 
1000  if(maximum < gPtBC->GetHistogram()->GetMaximum()) maximum = gPtBC->GetHistogram()->GetMaximum();
1001  if(minimum > gPtBC->GetHistogram()->GetMinimum()) minimum = gPtBC->GetHistogram()->GetMaximum();
1002  //printf("B-maximum %e, minimum %e\n",maximum,minimum);
1003 
1004  if(minimum < 0 ) minimum = 2.e-9 ;
1005  gPtBC->SetMaximum(maximum*2.);
1006  gPtBC->SetMinimum(minimum/2.);
1007 
1008  TLegend * legend = new TLegend(0.5,0.7,0.9,0.9);
1009  legend->AddEntry(gPt ,"From fit","P");
1010  legend->AddEntry(gPtBC,"From bin counting","P");
1011  legend->Draw();
1012  }
1013 
1014  gPt->SetMaximum(maximum*2.);
1015  gPt->SetMinimum(minimum/2.);
1016 
1017  cFitGraph->Print(Form("IMfigures/%s_%s_MassWidthPtSpectra_%s_%s.%s",
1018  kProdName.Data(),kCalorimeter.Data(),kParticle.Data(),hname.Data(),
1019  /*kFileName.Data(),*/ kPlotFormat.Data()));
1020 
1021  //-----------------------
1022  // Write results to file
1023  //-----------------------
1024  hRe->Write();
1025  if ( hMi ) hMi->Write();
1026 
1027  for(Int_t ipt = 0; ipt < nPt; ipt++)
1028  {
1029  if ( hIM[ipt] ) hIM[ipt]->Write();
1030 
1031  if ( kMix )
1032  {
1033  if (hMix [ipt] ) hMix [ipt]->Write();
1034  if (hMixCorrected[ipt] ) hMixCorrected[ipt]->Write();
1035  if (hRatio [ipt] ) hRatio [ipt]->Write();
1036  if (hSignal [ipt] ) hSignal [ipt]->Write();
1037  }
1038  }
1039 
1040  gPt ->Write();
1041  gMass ->Write();
1042  gWidth->Write();
1043  if ( kMix ) gPtBC->Write();
1044 
1045  // Do some final efficiency calculations if MC
1046  Efficiency(nPt, xPt, exPt, mesonPt, emesonPt, mesonPtBC, emesonPtBC, hname);
1047 }
1048 
1053 //-----------------------------------------------------------------------------
1054 void PlotGraphs(TString opt, Int_t first, Int_t last)
1055 {
1056  gStyle->SetOptTitle(1);
1057  gStyle->SetTitleOffset(2,"Y");
1058  gStyle->SetOptStat(0);
1059  gStyle->SetOptFit(000000);
1060 
1061  const Int_t nCombi = last-first+1;
1062  Int_t nCombiActive = 0;
1063 
1064  Float_t xmin = 0.7;
1065  Float_t ymin = 0.3;
1066  Float_t xmax = 0.9;
1067  Float_t ymax = 0.9;
1068 
1069  if(opt.Contains("Group"))
1070  {
1071  xmin = 0.3;
1072  ymin = 0.7;
1073  }
1074 
1075  if(opt.Contains("Sector"))
1076  {
1077  xmin = 0.65;
1078  ymin = 0.5;
1079  }
1080 
1081  TLegend * legend = new TLegend(xmin,ymin,xmax,ymax);
1082  legend->SetTextSize(0.05);
1083 
1084  TString particleN = " #pi^{0}";
1085  if(kParticle=="Eta") particleN = " #eta";
1086 
1087  // recover the graphs from output file
1088  //printf("%s\n",opt.Data());
1089 
1090 
1091  // Recover SUM
1092  TString opt2 = opt;
1093  if(opt.Contains("SMGroup")) opt2 = "SM";
1094 
1095  TGraphErrors* gSumPt = (TGraphErrors*) fout->Get(Form("gPt_Same%s" ,opt2.Data()));
1096  TGraphErrors* gSumPtBC = (TGraphErrors*) fout->Get(Form("gPtBC_Same%s" ,opt2.Data()));
1097  TGraphErrors* gSumMass = (TGraphErrors*) fout->Get(Form("gMass_Same%s" ,opt2.Data()));
1098  TGraphErrors* gSumWidth = (TGraphErrors*) fout->Get(Form("gWidth_Same%s",opt2.Data()));
1099  //printf("\t Sum %p %p %p %p\n",gSumPt,gSumPtBC,gSumMass,gSumWidth);
1100 
1101  if ( !gSumMass ) return ;
1102 
1103  gSumPt ->SetMarkerStyle(20);
1104  gSumPt ->SetMarkerColor(1);
1105  gSumPt ->SetLineColor (1);
1106 
1107  gSumMass->SetMarkerStyle(20);
1108  gSumMass->SetMarkerColor(1);
1109  gSumMass->SetLineColor (1);
1110 
1111  gSumWidth->SetMarkerStyle(20);
1112  gSumWidth->SetMarkerColor(1);
1113  gSumWidth->SetLineColor (1);
1114 
1115  if ( kMix )
1116  {
1117  gSumPtBC ->SetMarkerStyle(20);
1118  gSumPtBC ->SetMarkerColor(1);
1119  gSumPtBC ->SetLineColor (1);
1120  }
1121 
1122  legend->AddEntry(gSumPt,"Sum","P");
1123 
1124  // Normalize to total number of SMs/Sectors/Sides
1125  if ( opt != "SMGroup" )
1126  {
1127  for (Int_t ibin = 0; ibin < gSumPt->GetN(); ibin++)
1128  {
1129  gSumPt->GetY ()[ibin] /= nCombi;
1130  gSumPt->GetEY()[ibin] /= nCombi;
1131  if ( kMix )
1132  {
1133  gSumPtBC->GetY ()[ibin] /= nCombi;
1134  gSumPtBC->GetEY()[ibin] /= nCombi;
1135  }
1136  }
1137  }
1138 
1139  TGraphErrors* gSumTRDnotPt = 0;
1140  TGraphErrors* gSumTRDnotPtBC = 0;
1141  TGraphErrors* gSumTRDnotMass = 0;
1142  TGraphErrors* gSumTRDnotWidth = 0;
1143  TGraphErrors* gSumTRDyesPt = 0;
1144  TGraphErrors* gSumTRDyesPtBC = 0;
1145  TGraphErrors* gSumTRDyesMass = 0;
1146  TGraphErrors* gSumTRDyesWidth = 0;
1147 
1148  //
1150  //
1151  if ( kFirstTRDSM > 0 && opt == "SM" ) // it could be extended to sector/side
1152  {
1153  gSumTRDnotPt = (TGraphErrors*) fout->Get(Form("gPt_Same%sTRDNot" ,opt2.Data()));
1154  gSumTRDnotPtBC = (TGraphErrors*) fout->Get(Form("gPtBC_Same%sTRDNot" ,opt2.Data()));
1155  gSumTRDnotMass = (TGraphErrors*) fout->Get(Form("gMass_Same%sTRDNot" ,opt2.Data()));
1156  gSumTRDnotWidth = (TGraphErrors*) fout->Get(Form("gWidth_Same%sTRDNot",opt2.Data()));
1157  printf("\t Sum Not TRD %p %p %p %p\n",gSumTRDnotPt,gSumTRDnotPtBC,gSumTRDnotMass,gSumTRDnotWidth);
1158 
1159  gSumTRDnotPt ->SetMarkerStyle(20);
1160  gSumTRDnotPt ->SetMarkerColor(kGray);
1161  gSumTRDnotPt ->SetLineColor (kGray);
1162 
1163  gSumTRDnotMass->SetMarkerStyle(20);
1164  gSumTRDnotMass->SetMarkerColor(kGray);
1165  gSumTRDnotMass->SetLineColor (kGray);
1166 
1167  gSumTRDnotWidth->SetMarkerStyle(20);
1168  gSumTRDnotWidth->SetMarkerColor(kGray);
1169  gSumTRDnotWidth->SetLineColor (kGray);
1170 
1171  if ( kMix )
1172  {
1173  gSumTRDnotPtBC ->SetMarkerStyle(20);
1174  gSumTRDnotPtBC ->SetMarkerColor(kGray);
1175  gSumTRDnotPtBC ->SetLineColor (kGray);
1176  }
1177 
1178  legend->AddEntry(gSumTRDnotPt,"w/o TRD","P");
1179 
1180  for (Int_t ibin = 0; ibin < gSumPt->GetN(); ibin++)
1181  {
1182  gSumTRDnotPt->GetY ()[ibin] /= kFirstTRDSM;
1183  gSumTRDnotPt->GetEY()[ibin] /= kFirstTRDSM;
1184  if ( kMix )
1185  {
1186  gSumTRDnotPtBC->GetY ()[ibin] /= kFirstTRDSM;
1187  gSumTRDnotPtBC->GetEY()[ibin] /= kFirstTRDSM;
1188  }
1189  }
1190 
1191  gSumTRDyesPt = (TGraphErrors*) fout->Get(Form("gPt_Same%sTRDYes" ,opt2.Data()));
1192  gSumTRDyesPtBC = (TGraphErrors*) fout->Get(Form("gPtBC_Same%sTRDYes" ,opt2.Data()));
1193  gSumTRDyesMass = (TGraphErrors*) fout->Get(Form("gMass_Same%sTRDYes" ,opt2.Data()));
1194  gSumTRDyesWidth = (TGraphErrors*) fout->Get(Form("gWidth_Same%sTRDYes",opt2.Data()));
1195  printf("\t Sum Yes TRD %p %p %p %p\n",gSumTRDyesPt,gSumTRDyesPtBC,gSumTRDyesMass,gSumTRDyesWidth);
1196 
1197  gSumTRDyesPt ->SetMarkerStyle(20);
1198  gSumTRDyesPt ->SetMarkerColor(kCyan);
1199  gSumTRDyesPt ->SetLineColor (kCyan);
1200 
1201  gSumTRDyesMass->SetMarkerStyle(20);
1202  gSumTRDyesMass->SetMarkerColor(kCyan);
1203  gSumTRDyesMass->SetLineColor (kCyan);
1204 
1205  gSumTRDyesWidth->SetMarkerStyle(20);
1206  gSumTRDyesWidth->SetMarkerColor(kCyan);
1207  gSumTRDyesWidth->SetLineColor (kCyan);
1208 
1209  if ( kMix )
1210  {
1211  gSumTRDyesPtBC ->SetMarkerStyle(20);
1212  gSumTRDyesPtBC ->SetMarkerColor(kCyan);
1213  gSumTRDyesPtBC ->SetLineColor (kCyan);
1214  }
1215 
1216  legend->AddEntry(gSumTRDyesPt,"w/ TRD","P");
1217 
1218  Int_t nCovered = 10-kFirstTRDSM;
1219  for (Int_t ibin = 0; ibin < gSumPt->GetN(); ibin++)
1220  {
1221  gSumTRDyesPt->GetY ()[ibin] /= nCovered;
1222  gSumTRDyesPt->GetEY()[ibin] /= nCovered;
1223  if ( kMix )
1224  {
1225  gSumTRDyesPtBC->GetY ()[ibin] /= nCovered;
1226  gSumTRDyesPtBC->GetEY()[ibin] /= nCovered;
1227  }
1228  }
1229 
1230  } // TRD
1231 
1232  //
1233  // Get different combinations
1234  //
1235  TGraphErrors* gPt [nCombi];
1236  TGraphErrors* gPtBC [nCombi];
1237  TGraphErrors* gMass [nCombi];
1238  TGraphErrors* gWidth[nCombi];
1239 
1240  for(Int_t icomb = 0; icomb < nCombi; icomb ++)
1241  {
1242  gPt [icomb] = (TGraphErrors*) fout->Get(Form("gPt_%s%d" ,opt.Data(),icomb+first));
1243  gPtBC [icomb] = (TGraphErrors*) fout->Get(Form("gPtBC_%s%d" ,opt.Data(),icomb+first));
1244  gMass [icomb] = (TGraphErrors*) fout->Get(Form("gMass_%s%d" ,opt.Data(),icomb+first));
1245  gWidth[icomb] = (TGraphErrors*) fout->Get(Form("gWidth_%s%d",opt.Data(),icomb+first));
1246  //printf("\t %d %p %p %p %p\n",icomb,gPt[icomb],gPtBC[icomb],gMass[icomb],gWidth[icomb]);
1247 
1248  if ( !gPt[icomb] ) continue ;
1249 
1250  nCombiActive++;
1251 
1252  gPt [icomb]->SetMarkerStyle(modStyleIndex[icomb]);
1253  gPt [icomb]->SetMarkerColor(modColorIndex[icomb]);
1254  gPt [icomb]->SetLineColor (modColorIndex[icomb]);
1255 
1256  gMass[icomb]->SetMarkerStyle(modStyleIndex[icomb]);
1257  gMass[icomb]->SetMarkerColor(modColorIndex[icomb]);
1258  gMass[icomb]->SetLineColor (modColorIndex[icomb]);
1259 
1260  gWidth[icomb]->SetMarkerStyle(modStyleIndex[icomb]);
1261  gWidth[icomb]->SetMarkerColor(modColorIndex[icomb]);
1262  gWidth[icomb]->SetLineColor (modColorIndex[icomb]);
1263 
1264  gPt [icomb]->GetHistogram()->SetTitle(Form("#it{p}_{T} of %s from fit, %s " ,particleN.Data(), opt.Data()));
1265  gWidth[icomb]->GetHistogram()->SetTitle(Form("%s mass #sigma vs #it{p}_{T}, %s ",particleN.Data(), opt.Data()));
1266  gMass [icomb]->GetHistogram()->SetTitle(Form("%s mass #mu vs #it{p}_{T}, %s " ,particleN.Data(), opt.Data()));
1267 
1268  if ( kMix )
1269  {
1270  gPtBC [icomb]->SetMarkerStyle(modStyleIndex[icomb]);
1271  gPtBC [icomb]->SetMarkerColor(modColorIndex[icomb]);
1272  gPtBC [icomb]->SetLineColor (modColorIndex[icomb]);
1273  gPtBC [icomb]->GetHistogram()->SetTitle(Form("#it{p}_{T} of %s bin counted, %s ",particleN.Data(), opt.Data()));
1274  }
1275 
1276  if ( !opt.Contains("Group") )
1277  legend->AddEntry(gPt[icomb],Form("%s %d",opt.Data(),icomb),"P");
1278  else
1279  {
1280  if(icomb == 0) legend->AddEntry(gPt[icomb],"SM0+4+5+6+7+8+9","P");
1281  if(icomb == 1) legend->AddEntry(gPt[icomb],"SM1+2","P");
1282  if(icomb == 2) legend->AddEntry(gPt[icomb],"SM3+7","P");
1283  }
1284  }
1285 
1286  if ( !gMass[0] ) return ;
1287 
1288  //
1289  // PLOT
1290  //
1291  TCanvas *cFitGraph = new TCanvas(Form("cFitGraph_%sCombinations",opt.Data()),
1292  Form("Fit Graphs for %s combinations",opt.Data()),1000,1000);
1293  cFitGraph->Divide(2, 2);
1294 
1295  // Mass
1296  cFitGraph->cd(1);
1297  gPad->SetGridx();
1298  gPad->SetGridy();
1299 
1300  if ( kParticle == "Pi0" )
1301  {
1302  gMass[0]->SetMaximum(141);
1303  gMass[0]->SetMinimum(116);
1304  }
1305  else
1306  {
1307  gMass[0]->SetMaximum(600);
1308  gMass[0]->SetMinimum(450);
1309  }
1310  gMass[0]->Draw("AP");
1311  for(Int_t icomb = 1; icomb < nCombi; icomb ++)
1312  {
1313  if ( gMass[icomb] ) gMass[icomb]->Draw("P");
1314  }
1315 
1316  gSumMass->Draw("P");
1317  if ( kFirstTRDSM && opt == "SM" )
1318  {
1319  gSumTRDnotMass->Draw("P");
1320  gSumTRDyesMass->Draw("P");
1321  }
1322 
1323  legend->Draw();
1324 
1325  cFitGraph->cd(2);
1326  gPad->SetGridx();
1327  gPad->SetGridy();
1328 
1329  if ( kParticle == "Pi0" )
1330  {
1331  gWidth[0]->SetMaximum(17);
1332  gWidth[0]->SetMinimum(7);
1333  if(opt=="Side") gWidth[0]->SetMinimum(2);
1334  }
1335  else
1336  {
1337  gWidth[0]->SetMaximum(50);
1338  gWidth[0]->SetMinimum(10);
1339  }
1340 
1341  gWidth[0]->Draw("AP");
1342  for(Int_t icomb = 1; icomb < nCombi; icomb ++)
1343  {
1344  if ( gWidth[icomb] ) gWidth[icomb]->Draw("P");
1345  }
1346 
1347  gSumWidth->Draw("P");
1348  if ( kFirstTRDSM && opt == "SM" )
1349  {
1350  gSumTRDnotWidth->Draw("P");
1351  gSumTRDyesWidth->Draw("P");
1352  }
1353 
1354  legend->Draw();
1355 
1356  cFitGraph->cd(3);
1357 
1358  gPad->SetGridx();
1359  gPad->SetGridy();
1360  gPad->SetLogy();
1361 
1362  gPt[0]->SetMaximum(1);
1363  gPt[0]->SetMinimum(1e-8);
1364 
1365  gPt[0]->Draw("AP");
1366  for(Int_t icomb = 1; icomb < nCombi; icomb ++)
1367  {
1368  if ( gPt[icomb] ) gPt[icomb]->Draw("P");
1369  }
1370 
1371  gSumPt->Draw("P");
1372  if ( kFirstTRDSM && opt == "SM" )
1373  {
1374  gSumTRDnotPt->Draw("P");
1375  gSumTRDyesPt->Draw("P");
1376  }
1377 
1378  legend->Draw();
1379 
1380  if ( kMix )
1381  {
1382  cFitGraph->cd(4);
1383 
1384  gPad->SetGridx();
1385  gPad->SetGridy();
1386  gPad->SetLogy();
1387 
1388  gPtBC[0]->SetMaximum(1);
1389  gPtBC[0]->SetMinimum(1e-8);
1390 
1391  gPtBC[0]->Draw("AP");
1392 
1393  for(Int_t icomb = 1; icomb < nCombi; icomb ++)
1394  {
1395  if ( gPtBC[icomb] ) gPtBC[icomb]->Draw("P");
1396  }
1397 
1398  gSumPtBC->Draw("P");
1399  if ( kFirstTRDSM && opt == "SM" )
1400  {
1401  gSumTRDnotPtBC->Draw("P");
1402  gSumTRDyesPtBC->Draw("P");
1403  }
1404  legend->Draw();
1405  }
1406 
1407  cFitGraph->Print(Form("IMfigures/%s_%s_MassWidthPtSpectra_%s_%sCombinations.%s",
1408  kProdName.Data(),kCalorimeter.Data(),kParticle.Data(),opt.Data(),
1409  /*kFileName.Data(),*/kPlotFormat.Data()));
1410  //
1411  // Calculate the ratio to the sum
1412  //
1413  TGraphErrors* gRatPt [nCombi];
1414  TGraphErrors* gRatPtBC [nCombi];
1415  TGraphErrors* gRatMass [nCombi];
1416  TGraphErrors* gRatWidth[nCombi];
1417 
1418  TLegend * legendR = new TLegend(xmin,ymin,xmax,ymax);
1419  legendR->SetTextSize(0.05);
1420 
1421  for(Int_t icomb = 0; icomb < nCombi; icomb ++)
1422  {
1423  //printf("icomb %d\n",icomb);
1424  //printf("\t pt %p %p\n",gPt[icomb],gSumPt);
1425  gRatPt [icomb] = DivideGraphs(gPt [icomb],gSumPt ); // PlotUtils.C
1426  //printf("\t ptBC %p %p\n",gPtBC[icomb],gSumPtBC);
1427  gRatPtBC [icomb] = DivideGraphs(gPtBC [icomb],gSumPtBC ); // PlotUtils.C
1428  //printf("\t pt %p %p\n",gMass[icomb],gSumMass);
1429  gRatMass [icomb] = DivideGraphs(gMass [icomb],gSumMass ); // PlotUtils.C
1430  //printf("\t pt %p %p\n",gWidth[icomb],gSumWidth);
1431  gRatWidth[icomb] = DivideGraphs(gWidth[icomb],gSumWidth); // PlotUtils.C
1432 
1433  if(!gRatMass[icomb]) continue;
1434 
1435  gRatPt [icomb]->SetMarkerStyle(modStyleIndex[icomb]);
1436  gRatPt [icomb]->SetMarkerColor(modColorIndex[icomb]);
1437  gRatPt [icomb]->SetLineColor (modColorIndex[icomb]);
1438 
1439  gRatMass[icomb]->SetMarkerStyle(modStyleIndex[icomb]);
1440  gRatMass[icomb]->SetMarkerColor(modColorIndex[icomb]);
1441  gRatMass[icomb]->SetLineColor (modColorIndex[icomb]);
1442 
1443  gRatWidth[icomb]->SetMarkerStyle(modStyleIndex[icomb]);
1444  gRatWidth[icomb]->SetMarkerColor(modColorIndex[icomb]);
1445  gRatWidth[icomb]->SetLineColor (modColorIndex[icomb]);
1446 
1447  gRatPt [icomb]->GetHistogram()->SetTitle(Form("#it{p}_{T} of %s from fit, %s " ,particleN.Data(), opt.Data()));
1448  gRatWidth[icomb]->GetHistogram()->SetTitle(Form("%s mass #sigma vs #it{p}_{T}, %s ",particleN.Data(), opt.Data()));
1449  gRatMass [icomb]->GetHistogram()->SetTitle(Form("%s mass #mu vs #it{p}_{T}, %s " ,particleN.Data(), opt.Data()));
1450 
1451  gRatPt [icomb]->GetHistogram()->SetYTitle("Ratio to sum");
1452  gRatWidth[icomb]->GetHistogram()->SetYTitle("Ratio to sum");
1453  gRatMass [icomb]->GetHistogram()->SetYTitle("Ratio to sum");
1454 
1455  gRatPt [icomb]->GetHistogram()->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1456  gRatWidth[icomb]->GetHistogram()->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1457  gRatMass [icomb]->GetHistogram()->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1458 
1459  gRatPt [icomb]->GetHistogram()->SetTitleOffset(1.4,"Y");
1460  gRatWidth[icomb]->GetHistogram()->SetTitleOffset(1.4,"Y");
1461  gRatMass [icomb]->GetHistogram()->SetTitleOffset(1.4,"Y");
1462 
1463  if ( kMix )
1464  {
1465  gRatPtBC [icomb]->SetMarkerStyle(modStyleIndex[icomb]);
1466  gRatPtBC [icomb]->SetMarkerColor(modColorIndex[icomb]);
1467  gRatPtBC [icomb]->SetLineColor (modColorIndex[icomb]);
1468  gRatPtBC [icomb]->GetHistogram()->SetTitle(Form("#it{p}_{T} of %s bin counted, %s ",particleN.Data(), opt.Data()));
1469  gRatPtBC [icomb]->GetHistogram()->SetYTitle("Ratio to sum");
1470  gRatPtBC [icomb]->GetHistogram()->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1471  gRatPtBC [icomb]->GetHistogram()->SetTitleOffset(1.5,"Y");
1472  }
1473 
1474  if ( !opt.Contains("Group") )
1475  legendR->AddEntry(gPt[icomb],Form("%s %d",opt.Data(),icomb),"P");
1476  else
1477  {
1478  if(icomb == 0) legendR->AddEntry(gPt[icomb],"SM0+4+5+6+7+8+9","P");
1479  if(icomb == 1) legendR->AddEntry(gPt[icomb],"SM1+2","P");
1480  if(icomb == 2) legendR->AddEntry(gPt[icomb],"SM3+7","P");
1481  }
1482  } // combi loop
1483 
1484  //
1485  // PLOT
1486  //
1487  TCanvas *cFitGraphRatio = new TCanvas(Form("cFitGraphSumRatio_%sCombinations",opt.Data()),
1488  Form("Fit Graphs ratio to sum for %s combinations",opt.Data()),1000,1000);
1489  cFitGraphRatio->Divide(2, 2);
1490 
1491  // Mass
1492  cFitGraphRatio->cd(1);
1493  gPad->SetGridx();
1494  gPad->SetGridy();
1495 
1496  gRatMass[0]->SetMaximum(1.03);
1497  gRatMass[0]->SetMinimum(0.97);
1498 
1499  gRatMass[0]->Draw("AP");
1500  for(Int_t icomb = 1; icomb < nCombi; icomb ++)
1501  {
1502  if ( gRatMass[icomb] ) gRatMass[icomb]->Draw("P");
1503  }
1504 
1505  legendR->Draw();
1506 
1507  cFitGraphRatio->cd(2);
1508  gPad->SetGridx();
1509  gPad->SetGridy();
1510 
1511  gRatWidth[0]->SetMaximum(1.5);
1512  gRatWidth[0]->SetMinimum(0.7);
1513 
1514  gRatWidth[0]->Draw("AP");
1515  for(Int_t icomb = 1; icomb < nCombi; icomb ++)
1516  {
1517  if ( gRatWidth[icomb] ) gRatWidth[icomb]->Draw("P");
1518  }
1519 
1520  legendR->Draw();
1521 
1522  cFitGraphRatio->cd(3);
1523 
1524  gPad->SetGridx();
1525  gPad->SetGridy();
1526 // gPad->SetLogy();
1527 
1528  gRatPt[0]->SetMaximum(2);
1529  gRatPt[0]->SetMinimum(0);
1530 
1531  gRatPt[0]->Draw("AP");
1532  for(Int_t icomb = 1; icomb < nCombi; icomb ++)
1533  {
1534  if ( gRatPt[icomb] ) gRatPt[icomb]->Draw("P");
1535  }
1536 
1537  legendR->Draw();
1538 
1539  if ( kMix )
1540  {
1541  cFitGraphRatio->cd(4);
1542 
1543  gPad->SetGridx();
1544  gPad->SetGridy();
1545 // gPad->SetLogy();
1546 
1547  gRatPtBC[0]->SetMaximum(2);
1548  gRatPtBC[0]->SetMinimum(0);
1549 
1550  gRatPtBC[0]->Draw("AP");
1551  for(Int_t icomb = 1; icomb < nCombi; icomb ++)
1552  {
1553  if ( gRatPtBC[icomb] ) gRatPtBC[icomb]->Draw("P");
1554  }
1555 
1556  legendR->Draw();
1557  }
1558 
1559  cFitGraphRatio->Print(Form("IMfigures/%s_%s_MassWidthPtSpectra_%s_%sCombinations_RatioToSum.%s",
1560  kProdName.Data(),kCalorimeter.Data(),kParticle.Data(),opt.Data(),
1561  /*kFileName.Data(),*/kPlotFormat.Data()));
1562 
1563 
1564  // Ratio to corresponding TRD case
1565  if ( kFirstTRDSM && opt == "SM" ) // it could be extended for Sector/Side
1566  {
1567  //
1568  // Calculate the ratio to the sum
1569  //
1570  TGraphErrors* gRatTRDPt [nCombi];
1571  TGraphErrors* gRatTRDPtBC [nCombi];
1572  TGraphErrors* gRatTRDMass [nCombi];
1573  TGraphErrors* gRatTRDWidth[nCombi];
1574 
1575  TLegend * legendR = new TLegend(xmin,ymin,xmax,ymax);
1576  legendR->SetTextSize(0.05);
1577 
1578  for(Int_t icomb = 0; icomb < nCombi; icomb ++)
1579  {
1580  //printf("icomb %d\n",icomb);
1581  //printf("\t pt %p %p\n",gPt[icomb],gSumPt);
1582 
1583 
1584  if(icomb < kFirstTRDSM)
1585  {
1586  gRatTRDPt [icomb] = DivideGraphs(gPt [icomb],gSumTRDnotPt ); // PlotUtils.C
1587  gRatTRDPtBC [icomb] = DivideGraphs(gPtBC [icomb],gSumTRDnotPtBC ); // PlotUtils.C
1588  gRatTRDMass [icomb] = DivideGraphs(gMass [icomb],gSumTRDnotMass ); // PlotUtils.C
1589  gRatTRDWidth[icomb] = DivideGraphs(gWidth[icomb],gSumTRDnotWidth); // PlotUtils.C
1590  }
1591  else
1592  {
1593  gRatTRDPt [icomb] = DivideGraphs(gPt [icomb],gSumTRDyesPt ); // PlotUtils.C
1594  gRatTRDPtBC [icomb] = DivideGraphs(gPtBC [icomb],gSumTRDyesPtBC ); // PlotUtils.C
1595  gRatTRDMass [icomb] = DivideGraphs(gMass [icomb],gSumTRDyesMass ); // PlotUtils.C
1596  gRatTRDWidth[icomb] = DivideGraphs(gWidth[icomb],gSumTRDyesWidth); // PlotUtils.C
1597  }
1598 
1599 
1600  if(!gRatTRDMass[icomb]) continue;
1601 
1602  gRatTRDPt [icomb]->SetMarkerStyle(modStyleIndex[icomb]);
1603  gRatTRDPt [icomb]->SetMarkerColor(modColorIndex[icomb]);
1604  gRatTRDPt [icomb]->SetLineColor (modColorIndex[icomb]);
1605 
1606  gRatTRDMass[icomb]->SetMarkerStyle(modStyleIndex[icomb]);
1607  gRatTRDMass[icomb]->SetMarkerColor(modColorIndex[icomb]);
1608  gRatTRDMass[icomb]->SetLineColor (modColorIndex[icomb]);
1609 
1610  gRatTRDWidth[icomb]->SetMarkerStyle(modStyleIndex[icomb]);
1611  gRatTRDWidth[icomb]->SetMarkerColor(modColorIndex[icomb]);
1612  gRatTRDWidth[icomb]->SetLineColor (modColorIndex[icomb]);
1613 
1614  gRatTRDPt [icomb]->GetHistogram()->SetTitle(Form("#it{p}_{T} of %s from fit, %s " ,particleN.Data(), opt.Data()));
1615  gRatTRDWidth[icomb]->GetHistogram()->SetTitle(Form("%s mass #sigma vs #it{p}_{T}, %s ",particleN.Data(), opt.Data()));
1616  gRatTRDMass [icomb]->GetHistogram()->SetTitle(Form("%s mass #mu vs #it{p}_{T}, %s " ,particleN.Data(), opt.Data()));
1617 
1618  gRatTRDPt [icomb]->GetHistogram()->SetYTitle("Ratio to sum");
1619  gRatTRDWidth[icomb]->GetHistogram()->SetYTitle("Ratio to sum");
1620  gRatTRDMass [icomb]->GetHistogram()->SetYTitle("Ratio to sum");
1621 
1622  gRatTRDPt [icomb]->GetHistogram()->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1623  gRatTRDWidth[icomb]->GetHistogram()->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1624  gRatTRDMass [icomb]->GetHistogram()->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1625 
1626  gRatTRDPt [icomb]->GetHistogram()->SetTitleOffset(1.4,"Y");
1627  gRatTRDWidth[icomb]->GetHistogram()->SetTitleOffset(1.4,"Y");
1628  gRatTRDMass [icomb]->GetHistogram()->SetTitleOffset(1.4,"Y");
1629 
1630  if ( kMix )
1631  {
1632  gRatTRDPtBC [icomb]->SetMarkerStyle(modStyleIndex[icomb]);
1633  gRatTRDPtBC [icomb]->SetMarkerColor(modColorIndex[icomb]);
1634  gRatTRDPtBC [icomb]->SetLineColor (modColorIndex[icomb]);
1635  gRatTRDPtBC [icomb]->GetHistogram()->SetTitle(Form("#it{p}_{T} of %s bin counted, %s ",particleN.Data(), opt.Data()));
1636  gRatTRDPtBC [icomb]->GetHistogram()->SetYTitle("Ratio to sum");
1637  gRatTRDPtBC [icomb]->GetHistogram()->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1638  gRatTRDPtBC [icomb]->GetHistogram()->SetTitleOffset(1.5,"Y");
1639  }
1640 
1641  if ( !opt.Contains("Group") )
1642  legendR->AddEntry(gPt[icomb],Form("%s %d",opt.Data(),icomb),"P");
1643  else
1644  {
1645  if(icomb == 0) legendR->AddEntry(gPt[icomb],"SM0+4+5+6+7+8+9","P");
1646  if(icomb == 1) legendR->AddEntry(gPt[icomb],"SM1+2","P");
1647  if(icomb == 2) legendR->AddEntry(gPt[icomb],"SM3+7","P");
1648  }
1649  } // combi loop
1650 
1651  //
1652  // PLOT
1653  //
1654  TCanvas *cFitGraphRatioTRD = new TCanvas(Form("cFitGraphSumRatioTRD_%sCombinations",opt.Data()),
1655  Form("Fit Graphs ratio to sum for %s combinations for TRD cases",opt.Data()),1000,1000);
1656  cFitGraphRatioTRD->Divide(2, 2);
1657 
1658  // Mass
1659  cFitGraphRatioTRD->cd(1);
1660  gPad->SetGridx();
1661  gPad->SetGridy();
1662 
1663  gRatTRDMass[0]->SetMaximum(1.03);
1664  gRatTRDMass[0]->SetMinimum(0.97);
1665 
1666  gRatTRDMass[0]->Draw("AP");
1667  for(Int_t icomb = 1; icomb < nCombi; icomb ++)
1668  {
1669  if ( gRatTRDMass[icomb] ) gRatTRDMass[icomb]->Draw("P");
1670  }
1671 
1672  legendR->Draw();
1673 
1674  cFitGraphRatioTRD->cd(2);
1675  gPad->SetGridx();
1676  gPad->SetGridy();
1677 
1678  gRatTRDWidth[0]->SetMaximum(1.5);
1679  gRatTRDWidth[0]->SetMinimum(0.7);
1680 
1681  gRatTRDWidth[0]->Draw("AP");
1682  for(Int_t icomb = 1; icomb < nCombi; icomb ++)
1683  {
1684  if ( gRatTRDWidth[icomb] ) gRatTRDWidth[icomb]->Draw("P");
1685  }
1686 
1687  legendR->Draw();
1688 
1689  cFitGraphRatioTRD->cd(3);
1690 
1691  gPad->SetGridx();
1692  gPad->SetGridy();
1693  // gPad->SetLogy();
1694 
1695  gRatTRDPt[0]->SetMaximum(2);
1696  gRatTRDPt[0]->SetMinimum(0);
1697 
1698  gRatTRDPt[0]->Draw("AP");
1699  for(Int_t icomb = 1; icomb < nCombi; icomb ++)
1700  {
1701  if ( gRatTRDPt[icomb] ) gRatTRDPt[icomb]->Draw("P");
1702  }
1703 
1704  legendR->Draw();
1705 
1706  if ( kMix )
1707  {
1708  cFitGraphRatioTRD->cd(4);
1709 
1710  gPad->SetGridx();
1711  gPad->SetGridy();
1712  // gPad->SetLogy();
1713 
1714  gRatTRDPtBC[0]->SetMaximum(2);
1715  gRatTRDPtBC[0]->SetMinimum(0);
1716 
1717  gRatTRDPtBC[0]->Draw("AP");
1718  for(Int_t icomb = 1; icomb < nCombi; icomb ++)
1719  {
1720  if ( gRatTRDPtBC[icomb] ) gRatTRDPtBC[icomb]->Draw("P");
1721  }
1722 
1723  legendR->Draw();
1724  }
1725 
1726  cFitGraphRatioTRD->Print(Form("IMfigures/%s_%s_MassWidthPtSpectra_%s_%sCombinations_RatioToSumTRD.%s",
1727  kProdName.Data(),kCalorimeter.Data(),kParticle.Data(),opt.Data(),
1728  /*kFileName.Data(),*/kPlotFormat.Data()));
1729 
1730 
1731  TLegend * legendSums = new TLegend(xmin,ymin,xmax,ymax);
1732 
1733 
1734  //
1735  // PLOT
1736  //
1737  TCanvas *cFitGraphSums = new TCanvas(Form("cFitGraph_%sSums",opt.Data()),
1738  Form("Fit Graphs for %s Sums",opt.Data()),1000,1000);
1739  cFitGraphSums->Divide(2, 2);
1740 
1741  // Mass
1742  cFitGraphSums->cd(1);
1743  gPad->SetGridx();
1744  gPad->SetGridy();
1745 
1746  if ( kParticle == "Pi0" )
1747  {
1748  gSumMass->SetMaximum(135);
1749  gSumMass->SetMinimum(125);
1750  }
1751  else
1752  {
1753  gSumMass->SetMaximum(560);
1754  gSumMass->SetMinimum(520);
1755  }
1756 
1757  gSumMass->Draw("AP");
1758 
1759  gSumTRDnotMass->Draw("P");
1760  gSumTRDyesMass->Draw("P");
1761 
1762  legendSums->AddEntry(gSumMass,"Sum","P");
1763  legendSums->AddEntry(gSumTRDnotMass,"w/o TRD","P");
1764  legendSums->AddEntry(gSumTRDyesMass,"w/ TRD","P");
1765 
1766  legendSums->Draw();
1767 
1768  cFitGraphSums->cd(2);
1769  gPad->SetGridx();
1770  gPad->SetGridy();
1771 
1772  if ( kParticle == "Pi0" )
1773  {
1774  gSumWidth->SetMaximum(17);
1775  gSumWidth->SetMinimum(7);
1776  if(opt=="Side") gSumWidth->SetMinimum(2);
1777  }
1778  else
1779  {
1780  gSumWidth->SetMaximum(50);
1781  gSumWidth->SetMinimum(10);
1782  }
1783 
1784  gSumWidth->Draw("AP");
1785  gSumTRDnotWidth->Draw("P");
1786  gSumTRDyesWidth->Draw("P");
1787 
1788  legendSums->Draw();
1789 
1790  cFitGraphSums->cd(3);
1791 
1792  gPad->SetGridx();
1793  gPad->SetGridy();
1794  gPad->SetLogy();
1795 
1796  gSumPt->SetMaximum(1);
1797  gSumPt->SetMinimum(1e-8);
1798 
1799  gSumPt->Draw("AP");
1800  gSumTRDnotPt->Draw("P");
1801  gSumTRDyesPt->Draw("P");
1802 
1803  legendSums->Draw();
1804 
1805  if ( kMix )
1806  {
1807  cFitGraphSums->cd(4);
1808 
1809  gPad->SetGridx();
1810  gPad->SetGridy();
1811  gPad->SetLogy();
1812 
1813  gSumPtBC->SetMaximum(1);
1814  gSumPtBC->SetMinimum(1e-8);
1815 
1816  gSumPtBC->Draw("AP");
1817  gSumTRDnotPtBC->Draw("P");
1818  gSumTRDyesPtBC->Draw("P");
1819 
1820  legend->Draw();
1821  }
1822 
1823  cFitGraphSums->Print(Form("IMfigures/%s_%s_MassWidthPtSpectra_%s_%sSums.%s",
1824  kProdName.Data(),kCalorimeter.Data(),kParticle.Data(),opt.Data(),
1825  /*kFileName.Data(),*/kPlotFormat.Data()));
1826  }
1827 }
1828 
1833 //-----------------------------------------------------------------------------
1834 void PrimarySpectra(Int_t nPt, TArrayD xPtLimits, TArrayD xPt, TArrayD exPt)
1835 {
1836  gStyle->SetOptTitle(1);
1837  gStyle->SetTitleOffset(2,"Y");
1838  gStyle->SetOptStat(0);
1839  gStyle->SetOptFit(000000);
1840 
1841  TArrayD primMesonPt ; primMesonPt .Set(nPt);
1842  TArrayD eprimMesonPt ; eprimMesonPt .Set(nPt);
1843  TArrayD primMesonPtAcc; primMesonPtAcc.Set(nPt);
1844  TArrayD eprimMesonPtAcc; eprimMesonPtAcc.Set(nPt);
1845  TArrayD primMesonAcc ; primMesonAcc.Set(nPt);
1846  TArrayD eprimMesonAcc ; eprimMesonAcc.Set(nPt);
1847 
1848  TH2D* h2PrimMesonPtY = (TH2D*) fil->Get( Form("%s_hPrim%sRapidity" ,kHistoStartName.Data(), kParticle.Data() ) );
1849  TH2D* h2PrimMesonPtYAcc = (TH2D*) fil->Get( Form("%s_hPrim%sAccRapidity",kHistoStartName.Data(), kParticle.Data() ) );
1850 
1851  if ( !h2PrimMesonPtY ) return ;
1852 
1853  if(kSumw2)
1854  {
1855  h2PrimMesonPtY ->Sumw2();
1856  h2PrimMesonPtYAcc->Sumw2();
1857  }
1858 
1859 // if ( nEvt < 1 )
1860 // {
1861 // h2PrimMesonPtY ->Scale(1e10);
1862 // h2PrimMesonPtYAcc->Scale(1e10);
1863 // }
1864 // else
1865  {
1866  h2PrimMesonPtY ->Scale(1./nEvt);
1867  h2PrimMesonPtYAcc->Scale(1./nEvt);
1868  }
1869 
1870  h2PrimMesonPtY ->Write();
1871  h2PrimMesonPtYAcc->Write();
1872 
1873  TH1D* hY = h2PrimMesonPtY->ProjectionY("Rapidity",-1,-1);
1874  Int_t binYmin = hY->FindBin(-0.65);
1875  Int_t binYmax = hY->FindBin(+0.65);
1876  //printf("Y bin min %d, max %d\n",binYmin,binYmax);
1877  TH1D* hPrimMesonPt = (TH1D*) h2PrimMesonPtY->ProjectionX("PrimaryPt" ,binYmin ,binYmax );
1878  hPrimMesonPt->Write();
1879 
1880  TH1D* hYAcc = h2PrimMesonPtYAcc->ProjectionY("RapidityAcc",-1,-1);
1881  binYmin = hYAcc->FindBin(-0.65);
1882  binYmax = hYAcc->FindBin(+0.65);
1883  //printf("Y bin min %d, max %d\n",binYmin,binYmax);
1884  TH1D* hPrimMesonPtAcc = (TH1D*) h2PrimMesonPtYAcc->ProjectionX("PrimaryPtAccepted" ,binYmin ,binYmax );
1885  hPrimMesonPtAcc->Write();
1886 
1887  // PlotUtils.C
1888  ScaleBinBySize(hPrimMesonPt);
1889  ScaleBinBySize(hPrimMesonPtAcc);
1890 
1891  Double_t integralA = 0;
1892  Double_t integralAErr = 0;
1893  Double_t integralB = 0;
1894  Double_t integralBErr = 0;
1895  Double_t ptMin = -1;
1896  Double_t ptMax = -1;
1897  Int_t binMin= -1;
1898  Int_t binMax= -1;
1899  for(Int_t ibin = 0; ibin < nPt; ibin++ )
1900  {
1901  ptMin = xPtLimits[ibin];
1902  ptMax = xPtLimits[ibin+1];
1903 
1904  binMin = hPrimMesonPt->FindBin(ptMin);
1905  binMax = hPrimMesonPt->FindBin(ptMax)-1;
1906 
1907  // PlotUtils.C
1908  GetRangeIntegralAndError(hPrimMesonPt , binMin, binMax, integralA, integralAErr );
1909  GetRangeIntegralAndError(hPrimMesonPtAcc, binMin, binMax, integralB, integralBErr );
1910 
1911  primMesonPt [ibin] = integralA;
1912  eprimMesonPt [ibin] = integralAErr;
1913  primMesonPtAcc[ibin] = integralB;
1914  eprimMesonPtAcc[ibin] = integralBErr;
1915 
1916  if ( integralA > 0 && integralB > 0 )
1917  {
1918  primMesonAcc[ibin] = integralB / integralA ;
1919  eprimMesonAcc[ibin] = GetFractionError(integralB,integralA,integralBErr,integralAErr); // PlotUtils.C
1920  }
1921  else
1922  {
1923  primMesonAcc[ibin] = 0;
1924  eprimMesonAcc[ibin] = 0;
1925  }
1926 
1927 // printf("Bin %d, [%1.1f,%1.1f] bin size %1.1f, content num %2.2e, content den %2.2e: frac %2.3f+%2.3f\n",
1928 // ibin,xPtLimits[ibin],xPtLimits[ibin+1],2*exPt[ibin],primMesonPtAcc[ibin],primMesonPt[ibin],primMesonAcc[ibin],eprimMesonAcc[ibin] );
1929 
1930  primMesonPt [ibin]/=((exPt[ibin]*4)); // It should be 2 not 4???
1931  primMesonPtAcc[ibin]/=((exPt[ibin]*4));
1932 
1933  // // Scale biased production
1934  // if(kFileName.Contains("pp_7TeV_Pi0"))
1935  // {
1936  // primMeson[ibin]*=1.e6;
1937  // primMeson[ibin]*=360./100.;
1938  //
1939  // //eprimMeson[ibin]*=1.e6*1.e6;
1940  // eprimMeson[ibin]*=360./100.;
1941  //
1942  // mesonPtBC [0][ibin] /= (0.62 +xPt[ibin]*0.0017 );
1943  // emesonPtBC[0][ibin] /= (0.62 +xPt[ibin]*0.0017 );
1944  //
1945  // primMeson [ibin] /= (0.56 +xPt[ibin]*0.0096 );
1946  // eprimMeson[ibin] /= (0.56 +xPt[ibin]*0.0096 );
1947  // }
1948  } // pT bin loop
1949 
1950  TGraphErrors * gPrim = new TGraphErrors(nPt,xPt.GetArray(),primMesonPt.GetArray(),exPt.GetArray(),eprimMesonPt.GetArray());
1951  gPrim->SetName("Primary");
1952  gPrim->GetHistogram()->SetYTitle("d #it{N}/d #it{p}_{T}");
1953  gPrim->GetHistogram()->SetTitleOffset(1.4,"Y");
1954  gPrim->GetHistogram()->SetXTitle("#it{p}_{T,gen} (GeV/#it{c})");
1955  gPrim->GetHistogram()->SetTitleOffset(1.2,"X");
1956  gPrim->GetHistogram()->SetTitle("#it{p}_{T} spectra for |#it{Y}| < 0.65");
1957  gPrim->Write();
1958 
1959  TGraphErrors * gPrimAcc = new TGraphErrors(nPt,xPt.GetArray(),primMesonPtAcc.GetArray(),exPt.GetArray(),eprimMesonPtAcc.GetArray());
1960  gPrimAcc->SetName("PrimaryInAcceptance");
1961  gPrimAcc->GetHistogram()->SetYTitle("d#it{N}/d #it{p}_{T}");
1962  gPrimAcc->GetHistogram()->SetTitle(Form("|#it{Y}| < 0.65 in %s",kCalorimeter.Data()));
1963  gPrimAcc->Write();
1964 
1965  // Acceptance
1966  TH1D* hAcc = (TH1D*) hPrimMesonPtAcc->Clone("hAcceptance");
1967  hAcc->Divide(hPrimMesonPt);
1968 
1969  TGraphErrors * gAcc = new TGraphErrors(nPt,xPt.GetArray(),primMesonAcc.GetArray(),exPt.GetArray(),eprimMesonAcc.GetArray());
1970  gAcc->SetName("Acceptance");
1971  gAcc->GetHistogram()->SetYTitle("Acceptance");
1972  gAcc->GetHistogram()->SetTitleOffset(1.4,"Y");
1973  gAcc->GetHistogram()->SetXTitle("#it{p}_{T,gen} (GeV/#it{c})");
1974  gAcc->GetHistogram()->SetTitleOffset(1.2,"X");
1975  gAcc->GetHistogram()->SetTitle(Form("Acceptance for |#it{Y}| < 0.65 in %s",kCalorimeter.Data()));
1976  gAcc->Write();
1977 
1978  // Plot spectra and acceptance
1979  //
1980  TCanvas *cAcc = new TCanvas(Form("cAcceptance"),
1981  Form("Primary generated p_{T} spectra and acceptance"),2*600,600);
1982  cAcc->Divide(2, 1);
1983 
1984  cAcc->cd(1);
1985  //gPad->SetGridx();
1986  //gPad->SetGridy();
1987  gPad->SetLogy();
1988 
1989  Double_t maximum = gPrim->GetHistogram()->GetMaximum();
1990  Double_t minimum = gPrim->GetHistogram()->GetMinimum();
1991 
1992  if ( gPrimAcc->GetMaximum() > maximum ) maximum = gPrimAcc->GetMaximum() ;
1993  if ( gPrimAcc->GetMinimum() > minimum ) minimum = gPrimAcc->GetMinimum() ;
1994 
1995  gPrim->SetMaximum(maximum*10);
1996  gPrim->SetMinimum(minimum/10);
1997 
1998  gPrim ->Draw("AP");
1999  gPrimAcc->Draw("P");
2000 
2001  gPrim ->SetMarkerColor(1);
2002  gPrimAcc->SetMarkerColor(4);
2003  gPrim ->SetLineColor(1);
2004  gPrimAcc->SetLineColor(4);
2005  gPrim ->SetMarkerStyle(24);
2006  gPrimAcc->SetMarkerStyle(24);
2007 
2008  hPrimMesonPt ->Draw("same");
2009  hPrimMesonPtAcc->Draw("same");
2010 
2011  hPrimMesonPt ->Draw("Hsame");
2012  hPrimMesonPtAcc->Draw("Hsame");
2013  hPrimMesonPt ->SetLineColor(1);
2014  hPrimMesonPtAcc->SetLineColor(4);
2015 
2016  TLegend * legendS = new TLegend(0.4,0.7,0.9,0.9);
2017  legendS->AddEntry(gPrim,"|Y| < 0.65","P");
2018  legendS->AddEntry(gPrimAcc,Form("Both in %s",kCalorimeter.Data()),"P");
2019 
2020  legendS->AddEntry(hPrimMesonPt,"Histo |Y| < 0.65","L");
2021  legendS->AddEntry(hPrimMesonPtAcc,Form("Histo Both in %s",kCalorimeter.Data()),"L");
2022  legendS->Draw();
2023 
2024  cAcc->cd(2);
2025  gPad->SetGridx();
2026  gPad->SetGridy();
2027  //gPad->SetLogy();
2028 
2029  //gAcc->SetMaximum(1);
2030  gAcc->SetMaximum(gAcc->GetHistogram()->GetMaximum()*1.3);
2031  gAcc->SetMinimum(0);
2032 
2033  gAcc->Draw("AP");
2034  gAcc->SetMarkerColor(1);
2035  gAcc->SetMarkerStyle(24);
2036 
2037  hAcc->Draw("Hsame");
2038  hAcc->SetLineColor(4);
2039 
2040  TLegend * legendA = new TLegend(0.7,0.75,0.9,0.9);
2041  legendA->AddEntry(gAcc,"Graph","P");
2042  legendA->AddEntry(hAcc,"Histo","L");
2043 
2044  legendA->Draw();
2045 
2046  cAcc->Print(Form("IMfigures/%s_%s_PrimarySpectraAcceptance_%s.%s",
2047  kProdName.Data(),kCalorimeter.Data(),kParticle.Data(),
2048  /*kFileName.Data(),*/ kPlotFormat.Data()));
2049 }
2050 
2076 //-----------------------------------------------------------------------------
2077 void InvMassFit
2078 (TString prodname = "LHC18c3_NystromOn",//"LHC17l3b_fast",
2079  TString filename = "AnalysisResults",
2080  TString histoDir = "Pi0IM_GammaTrackCorr_EMCAL",
2081  TString histoList = "default",
2082  TString histoStart = "",
2083  TString calorimeter = "EMCAL",
2084  TString particle = "Pi0",
2085  Float_t nPairMin = 20,
2086  Bool_t sumw2 = kTRUE,
2087  Int_t rebin = 1,
2088  Int_t firstTRD = 6,
2089  Bool_t mixed = kTRUE,
2090  Bool_t truncmix = kFALSE,
2091  Int_t pol = 1,
2092  Bool_t drawAll = kFALSE,
2093  Bool_t drawMCAll = kFALSE,
2094  Bool_t drawPerSM = kTRUE,
2095  Bool_t drawPerSMGr = kFALSE,
2096  Bool_t drawPerSector= kFALSE,
2097  Bool_t drawPerSide = kFALSE,
2098  TString plotFormat = "eps")
2099 {
2100  // Load plotting utilities
2101  //gROOT->Macro("$ALICE_PHYSICS/PWGGA/CaloTrackCorrelations/macros/PlotUtils.C"); // Already in include
2102  //gROOT->Macro("$MACROS/style.C");//Set different root style parameters
2103 
2104  kSumw2 = sumw2;
2105  kPolN = pol;
2106  kMix = mixed;
2107  kTrunMixFunc= truncmix;
2108  kParticle = particle;
2109  kProdName = prodname;
2110  kFileName = filename;
2111  kPlotFormat = plotFormat;
2112  kCalorimeter= calorimeter;
2113  kNPairCut = nPairMin;
2114  kRebin = rebin;
2115  kFirstTRDSM = firstTRD;
2116 
2117  if(histoStart !="")
2118  kHistoStartName = histoStart;
2119  else
2120  {
2121  kHistoStartName = "AnaPi0_Calo0";
2122  if(kCalorimeter=="DCAL")
2123  kHistoStartName = "AnaPi0_Calo1";
2124  }
2125 
2126  if(drawPerSMGr) drawPerSM = kTRUE;
2127 
2128  //---------------------------------------
2129  // Get input file
2130  //---------------------------------------
2131  Int_t ok = GetFileAndEvents(histoDir, histoList);
2132 
2133 // kProdName.ReplaceAll("module/","");
2134 // kProdName.ReplaceAll("TCardChannel","");
2135  kProdName.ReplaceAll("/","_");
2136 // kProdName.ReplaceAll("2_","_");
2137 // kProdName.ReplaceAll("__","_");
2138  if( kFileName !="AnalysisResults" && !kFileName.Contains("Scale") )
2139  kProdName+=kFileName;
2140 
2141  printf("Settings: prodname <%s>, filename <%s>, histoDir <%s>, histoList <%s>,\n"
2142  " \t histo start name <%s>, particle <%s>, calorimeter <%s>, \n"
2143  " \t mix %d, kPolN %d, sumw2 %d, n pairs cut %2.2e\n",
2144  kProdName.Data(), kFileName.Data(),histoDir.Data(),histoList.Data(), kHistoStartName.Data(),
2145  kParticle.Data(), kCalorimeter.Data(), kMix, kPolN, kSumw2, kNPairCut);
2146 
2147  if ( !ok )
2148  {
2149  printf("Could not recover file <%p> or dir <%p> or list <%p>\n",fil,direc,lis);
2150  return ;
2151  }
2152 
2153  // kNPairCut/=nEvt;
2154 
2155  //---------------------------------------
2156  // Set-up output file with processed histograms
2157  //---------------------------------------
2158 
2159  // Open output file
2160  fout = new TFile(Form("IMfigures/%s_%s_MassWidthPtHistograms_%s.root",
2161  kProdName.Data(),kCalorimeter.Data(),kParticle.Data() /*,kFileName.Data()*/), "recreate");
2162  printf("Open Output File: IMfigures/%s_%s_MassWidthPtHistograms_%s.root\n",
2163  kProdName.Data(),kCalorimeter.Data(),kParticle.Data());//,kFileName.Data());
2164 
2165  //---------------------------------------
2166  // Set the pt bins and total range
2167  //---------------------------------------
2168 
2169  const Int_t nPtEta = 13;
2170  Double_t xPtLimitsEta[] = {2.,3.,4.,5.,6.,8.,10.,14.,18.,24.,28.,34.,40.,50.};
2171 
2172 // const Int_t nPtPi0 = 7;
2173 // Double_t xPtLimitsPi0[] = {1.5,2,2.5,3.,3.5,4.,5.,6.};
2174 
2175  const Int_t nPtPi0 = 23;
2176  Double_t xPtLimitsPi0[] = {1.6,1.8,2.,2.2,2.4,2.6,2.8,3.,3.4,3.8,4.2,4.6,5.,5.5,6.,6.5,7.,7.5,8.,9.,10.,12.,15.,17.,20.};
2177 
2178  Int_t nPt = nPtPi0;
2179  if(kParticle == "Eta") nPt = nPtEta;
2180 
2181  TArrayD xPtLimits; xPtLimits.Set(nPt+1);
2182  TArrayD xPt ; xPt .Set(nPt);
2183  TArrayD exPt ; exPt .Set(nPt);
2184 
2185  for(Int_t i = 0; i < nPt; i++)
2186  {
2187  if(kParticle == "Pi0")
2188  {
2189  xPtLimits.SetAt(xPtLimitsPi0[i],i);
2190  exPt .SetAt((xPtLimitsPi0[i+1]-xPtLimitsPi0[i])/2,i);
2191  xPt .SetAt(xPtLimitsPi0[i]+exPt[i],i);
2192  }
2193  else
2194  {
2195  xPtLimits.SetAt(xPtLimitsEta[i],i);
2196  exPt .SetAt((xPtLimitsEta[i+1]-xPtLimitsEta[i])/2,i);
2197  xPt .SetAt(xPtLimitsEta[i]+exPt[i],i);
2198  }
2199 
2200  //printf("%s pT bin %d, pT %2.2f, dpT %2.2f, limit %2.2f\n",kParticle.Data(),i,xPt[i],exPt[i],xPtLimits[i]);
2201  }
2202  // Extra entry
2203  if ( kParticle == "Pi0" ) xPtLimits.SetAt(xPtLimitsPi0[nPt],nPt);
2204  else xPtLimits.SetAt(xPtLimitsEta[nPt],nPt);
2205 
2206  TString comment = "";
2207  TString leg = "";
2208  TString hname = "";
2209 
2210  //---------------------------------------
2211  // Fitting function and histograms initialization
2212  //---------------------------------------
2213  SetFitFun();
2214 
2215  TH2F * hRe ;
2216  TH2F * hMi ;
2217 
2218  //=======================================
2219  // Get histograms, project per pT bin and fit
2220  // Do it for different kind of histograms:
2221  // Total pairs, MC pairs, pairs per SM combinations
2222  //=======================================
2223 
2224  //---------------------------------------
2225  // MC input
2226  //---------------------------------------
2227 
2228  if ( drawMCAll )
2229  {
2230  // Get the generated primary pi0 spectra, for efficiency calculation
2231  PrimarySpectra(nPt,xPtLimits, xPt, exPt);
2232 
2233  // Reconstructed pT of real meson pairs
2234  //hRe = (TH2F *) fil->Get( Form("%s_hMCOrgMass_2", kHistoStartName.Data()) ); // Pi0
2235  //hRe = (TH2F *) fil->Get( Form("%s_hMCOrgMass_3", kHistoStartName.Data()) ); // Eta
2236  hRe = (TH2F *) fil->Get( Form("%s_hMC%sMassPtRec", kHistoStartName.Data(), kParticle.Data()) );
2237  comment = "MC pT reconstructed";
2238  leg = "MC pT reconstructed";
2239  hname = "MCpTReco" ;
2240  ProjectAndFit(nPt,comment,leg,hname,xPtLimits, xPt, exPt, hRe, 0x0 );
2241 
2242 
2243  // Reconstructed pT of real eta pairs
2244  hRe = (TH2F *) fil->Get( Form("%s_hMC%sMassPtTrue", kHistoStartName.Data(), kParticle.Data()) );
2245  comment = "MC pT generated";
2246  leg = "MC pT generated";
2247  hname = "MCpTGener" ;
2248  ProjectAndFit(nPt,comment,leg,hname,xPtLimits, xPt, exPt, hRe, 0x0 );
2249  } // MC treatment
2250 
2251  //---------------------------------------
2252  // All SM
2253  //---------------------------------------
2254 
2255  if ( drawAll )
2256  {
2257  if ( !lis )
2258  {
2259  hRe = (TH2F *) fil->Get(Form("%s_hRe_cen0_pidbit0_asy0_dist1", kHistoStartName.Data()));
2260  hMi = (TH2F *) fil->Get(Form("%s_hMi_cen0_pidbit0_asy0_dist1", kHistoStartName.Data()));
2261  }
2262  else
2263  {
2264  hRe = (TH2F *) lis->FindObject(Form("%s_hRe_cen0_pidbit0_asy0_dist1", kHistoStartName.Data()));
2265  hMi = (TH2F *) lis->FindObject(Form("%s_hMi_cen0_pidbit0_asy0_dist1", kHistoStartName.Data()));
2266  }
2267 
2268  printf("histo Re %p - Mix %p All SM\n",hRe,hMi);
2269 
2270  if( !hMi ) kMix = kFALSE;
2271 
2272  comment = "All SM";
2273  leg = "All SM";
2274  hname = "AllSM";
2275  ProjectAndFit(nPt,comment,leg,hname,xPtLimits, xPt, exPt, hRe, hMi );
2276  }
2277 
2278  //-------------------------------------------
2279  // Histograms for cluster pairs in different
2280  // combinations of SM
2281  //-------------------------------------------
2282 
2283  //
2284  // Pairs in same SM
2285  //
2286  if ( drawPerSM )
2287  {
2288  Int_t imod;
2289  Int_t firstMod = 0;
2290  Int_t lastMod = 11;
2291 
2292  if ( kCalorimeter=="DCAL" )
2293  {
2294  firstMod = 12;
2295  lastMod = 19;
2296  }
2297 
2298  // Initialize SM grouping histograms
2299  // 0 : SM 0+4+5+6+8+9
2300  // 1 : SM 1+2
2301  // 2 : SM 3+7
2302  TH2F * hReSMGroups[3] ;
2303  TH2F * hMiSMGroups[3] ;
2304  TH2F * hReSM = 0;
2305  TH2F * hMiSM = 0;
2306  TH2F * hReSMTRDYes = 0;
2307  TH2F * hMiSMTRDYes = 0;
2308  TH2F * hReSMTRDNot = 0;
2309  TH2F * hMiSMTRDNot = 0;
2310 
2311  for(Int_t imodgroup = 0; imodgroup<=2; imodgroup++)
2312  {
2313  hReSMGroups[imodgroup] = hMiSMGroups[imodgroup] = 0x0;
2314  }
2315 
2316  for(imod = firstMod; imod<=lastMod; imod++)
2317  {
2318  if(!lis)
2319  {
2320  hRe = (TH2F *) fil->Get(Form("%s_hReMod_%d", kHistoStartName.Data(), imod));
2321  hMi = (TH2F *) fil->Get(Form("%s_hMiMod_%d", kHistoStartName.Data(), imod));
2322  }
2323  else
2324  {
2325  hRe = (TH2F *) lis->FindObject(Form("%s_hReMod_%d", kHistoStartName.Data(), imod));
2326  hMi = (TH2F *) lis->FindObject(Form("%s_hMiMod_%d", kHistoStartName.Data(), imod));
2327  }
2328 
2329  //printf("histo mod %d Re %p - Mix %p\n",imod, hRe,hMi);
2330 
2331  comment = Form("both clusters in SM %d",imod);
2332  leg = Form("SM %d",imod);
2333  hname = Form("SM%d" ,imod);
2334 
2335  ProjectAndFit(nPt,comment,leg,hname,xPtLimits, xPt, exPt, hRe, hMi );
2336 
2337  // Add all SM
2338  if( imod == firstMod )
2339  {
2340  hReSM = (TH2F*) hRe->Clone("h2D_Re_IM_SM");
2341  if(kMix) hMiSM = (TH2F*) hMi->Clone("h2D_Mi_IM_SM");
2342  }
2343  else
2344  {
2345  hReSM->Add( hRe );
2346  if(kMix) hMiSM ->Add( hMi );
2347  }
2348 
2349  // Add TRD covered or not SM
2350  if( kFirstTRDSM > 0 )
2351  {
2352  if ( imod >= kFirstTRDSM && imod < 10 )
2353  {
2354  if( imod == kFirstTRDSM )
2355  {
2356  hReSMTRDYes = (TH2F*) hRe->Clone("h2D_Re_IM_SM_TRD_Yes");
2357  if(kMix) hMiSMTRDYes = (TH2F*) hMi->Clone("h2D_Mi_IM_SM_TRD_Yes");
2358  }
2359  else
2360  {
2361  hReSMTRDYes->Add( hRe );
2362  if(kMix) hMiSMTRDYes ->Add( hMi );
2363  }
2364  }
2365  }
2366 
2367  // Group SMs in with particular behaviors
2368  if ( drawPerSMGr )
2369  {
2370  if ( imod == 0 )
2371  {
2372  hReSMGroups[0] = (TH2F*) hRe->Clone("h2D_Re_IM_SM045689");
2373  if(kMix)
2374  hMiSMGroups[0] = (TH2F*) hMi->Clone("h2D_Mi_IM_SM045689");
2375  }
2376  else if ( imod == 1 )
2377  {
2378  hReSMGroups[1] = (TH2F*) hRe->Clone("h2D_Re_IM_SM12");
2379  if(kMix) hMiSMGroups[1] = (TH2F*) hMi->Clone("h2D_Mi_IM_SM12");
2380  }
2381  else if ( imod == 3 )
2382  {
2383  hReSMGroups[2] = (TH2F*) hRe->Clone("h2D_Re_IM_SM37");
2384  if(kMix) hMiSMGroups[2] = (TH2F*) hMi->Clone("h2D_Mi_IM_SM37");
2385  }
2386  else if ( imod == 2 )
2387  {
2388  hReSMGroups[1]->Add(hRe);
2389  if(kMix) hMiSMGroups[1]->Add(hMi);
2390  }
2391  else if ( imod == 7 )
2392  {
2393  hReSMGroups[2]->Add(hRe);
2394  if(kMix) hMiSMGroups[2]->Add(hMi);
2395  }
2396  else if ( imod < 10 )
2397  {
2398  hReSMGroups[0]->Add(hRe);
2399  if(kMix) hMiSMGroups[0]->Add(hMi);
2400  }
2401  } // do SM groups addition
2402  } // SM loop
2403 
2404  // Sum of pairs in same SM
2405  comment = "both clusters in same SM";
2406  leg = "Same SM";
2407  hname = "SameSM";
2408  ProjectAndFit(nPt,comment,leg,hname,xPtLimits, xPt, exPt, hReSM, hMiSM );
2409 
2410  if ( kFirstTRDSM > 0 )
2411  {
2412  // Sum of pairs in same SM with TRD
2413  comment = "both clusters in same SM, TRD covered";
2414  leg = "Same SM, TRD Yes";
2415  hname = "SameSMTRDYes";
2416  ProjectAndFit(nPt,comment,leg,hname,xPtLimits, xPt, exPt, hReSMTRDYes, hMiSMTRDYes );
2417 
2418  // Sum of pairs in same SM without TRD
2419  hReSMTRDNot = (TH2F*) hReSM->Clone("h2D_Re_IM_SM_TRD_Not");
2420  hReSMTRDNot->Add(hReSMTRDYes,-1);
2421  if(kMix)
2422  {
2423  hMiSMTRDNot = (TH2F*) hMiSM->Clone("h2D_Mi_IM_SM_TRD_Not");
2424  hMiSMTRDNot->Add(hMiSMTRDYes,-1);
2425  }
2426  comment = "both clusters in same SM, NOT TRD covered";
2427  leg = "Same SM, TRD No";
2428  hname = "SameSMTRDNot";
2429  ProjectAndFit(nPt,comment,leg,hname,xPtLimits, xPt, exPt, hReSMTRDNot, hMiSMTRDNot );
2430  }
2431 
2432  // Fit parameters final plotting
2433  PlotGraphs("SM" ,firstMod ,lastMod );
2434 
2435  // Group SMs in with particular behaviors
2436  if( drawPerSMGr )
2437  {
2438  for(Int_t imodgroup = 0; imodgroup<=2; imodgroup++)
2439  {
2440  //printf("histo modgroup %d Re %p - Mix %p\n",
2441  // imod, hReSMGroups[imodgroup], hMiSMGroups[imodgroup]);
2442 
2443  comment = Form("both clusters in same SM, group %d",imodgroup);
2444  leg = Form("SMGroup %d",imodgroup);
2445  hname = Form("SMGroup%d" ,imodgroup);
2446 
2447  if ( imodgroup != 0 )
2448  {
2449  hReSMGroups[imodgroup]->Scale(1./2.);
2450  if ( kMix ) hMiSMGroups[imodgroup]->Scale(1./2.);
2451  }
2452  else
2453  {
2454  hReSMGroups[imodgroup]->Scale(1./6.);
2455  if ( kMix ) hMiSMGroups[imodgroup]->Scale(1./6.);
2456  }
2457 
2458  ProjectAndFit(nPt, comment, leg, hname, xPtLimits, xPt, exPt,
2459  hReSMGroups[imodgroup], hMiSMGroups[imodgroup] );
2460  }
2461 
2462  // Fit parameters final plotting
2463  PlotGraphs("SMGroup",0 ,2 );
2464  } // Per SM groups
2465 
2466  } // Per SM
2467 
2468  //
2469  // Pairs in same sector
2470  //
2471  if(drawPerSector)
2472  {
2473  Int_t firstSector = 0;
2474  Int_t lastSector = 5;
2475  if ( kCalorimeter=="DCAL" )
2476  {
2477  firstSector = 6;
2478  lastSector = 9;
2479  }
2480 
2481  TH2F * hReSector = 0;
2482  TH2F * hMiSector = 0;
2483  Int_t isector;
2484  for(isector = firstSector; isector<=lastSector; isector++)
2485  {
2486  if(!lis)
2487  {
2488  hRe = (TH2F *) fil->Get(Form("%s_hReSameSectorEMCALMod_%d", kHistoStartName.Data(), isector));
2489  hMi = (TH2F *) fil->Get(Form("%s_hMiSameSectorEMCALMod_%d", kHistoStartName.Data(), isector));
2490  }
2491  else
2492  {
2493  hRe = (TH2F *) lis->FindObject(Form("%s_hReSameSectorEMCALMod_%d", kHistoStartName.Data(), isector));
2494  hMi = (TH2F *) lis->FindObject(Form("%s_hMiSameSectorEMCALMod_%d", kHistoStartName.Data(), isector));
2495  }
2496 
2497  // Add all Sectors
2498  if( isector == firstSector )
2499  {
2500  hReSector = (TH2F*) hRe->Clone("h2D_Re_IM_Sector");
2501  if(kMix) hMiSector = (TH2F*) hMi->Clone("h2D_Mi_IM_Sector");
2502  }
2503  else
2504  {
2505  hReSector->Add( hRe );
2506  if(kMix) hMiSector ->Add( hMi );
2507  }
2508 
2509  //printf("histo sector %d Re %p - Mix %p\n",isector, hRe,hMi);
2510 
2511  comment = Form("both clusters in Sector %d",isector);
2512  leg = Form("Sector %d",isector);
2513  hname = Form("Sector%d" ,isector);
2514 
2515  ProjectAndFit(nPt,comment,leg,hname,xPtLimits, xPt, exPt, hRe, hMi );
2516  }
2517 
2518  // Sum of pairs in same Sector
2519  comment = "both clusters in same Sector";
2520  leg = "Same Sector";
2521  hname = "SameSector";
2522  ProjectAndFit(nPt,comment,leg,hname,xPtLimits, xPt, exPt, hReSector, hMiSector );
2523 
2524  // Fit parameters final plotting
2525  PlotGraphs("Sector" ,firstSector,lastSector);
2526  } // Per sector
2527 
2528 
2529  // Pairs in same side
2530  //
2531  if(drawPerSide)
2532  {
2533  Int_t firstSide = 0;
2534  Int_t lastSide = 9;
2535 
2536  if ( kCalorimeter=="DCAL" )
2537  {
2538  firstSide = 10;
2539  lastSide = 15;
2540  }
2541 
2542  TH2F * hReSide = 0;
2543  TH2F * hMiSide = 0;
2544  Int_t iside;
2545  for(iside = firstSide; iside <= lastSide; iside++)
2546  {
2547  if(!lis)
2548  {
2549  hRe = (TH2F *) fil->Get(Form("%s_hReSameSideEMCALMod_%d", kHistoStartName.Data(), iside));
2550  hMi = (TH2F *) fil->Get(Form("%s_hMiSameSideEMCALMod_%d", kHistoStartName.Data(), iside));
2551  }
2552  else
2553  {
2554  hRe = (TH2F *) lis->FindObject(Form("%s_hReSameSideEMCALMod_%d", kHistoStartName.Data(), iside));
2555  hMi = (TH2F *) lis->FindObject(Form("%s_hMiSameSideEMCALMod_%d", kHistoStartName.Data(), iside));
2556  }
2557 
2558  // Add all Sides
2559  if( iside == firstSide )
2560  {
2561  hReSide = (TH2F*) hRe->Clone("h2D_Re_IM_Side");
2562  if(kMix) hMiSide = (TH2F*) hMi->Clone("h2D_Mi_IM_Side");
2563  }
2564  else
2565  {
2566  hReSide->Add( hRe );
2567  if(kMix) hMiSide ->Add( hMi );
2568  }
2569 
2570  //printf("histo side %d Re %p - Mix %p\n",iside, hRe,hMi);
2571 
2572  comment = Form("both clusters in Side %d",iside);
2573  leg = Form("Side %d",iside);
2574  hname = Form("Side%d" ,iside);
2575 
2576  ProjectAndFit(nPt,comment,leg,hname,xPtLimits, xPt, exPt, hRe, hMi );
2577  }
2578 
2579  // Sum of pairs in same Side
2580  comment = "both clusters in same Side";
2581  leg = "Same Side";
2582  hname = "SameSide";
2583  ProjectAndFit(nPt,comment,leg,hname,xPtLimits, xPt, exPt, hReSide, hMiSide );
2584 
2585  // Fit parameters final plotting
2586  PlotGraphs("Side" ,firstSide ,lastSide );
2587 
2588  } // per Side
2589 
2590  fout->Close();
2591 }
hSignal
TH1D * hSignal
Definition: CalculateAveragePt.C:87
kFirstTRDSM
Float_t kFirstTRDSM
Minimum number of cluster pairs in pi0 or eta window.
Definition: InvMassFit.C:47
filename
const char * filename
Definition: TestFCM.C:1
FunctionCrystalBall
static Double_t FunctionCrystalBall(Double_t *x, Double_t *par)
Definition: PlotUtils.C:25
FunctionTruncatedPolPi0
static Double_t FunctionTruncatedPolPi0(Double_t *x, Double_t *par)
Definition: PlotUtils.C:112
ymax
const Double_t ymax
Definition: AddTaskCFDStar.C:7
Double_t
double Double_t
Definition: External.C:58
Draw
void Draw(const char *filename, const char *title="", const char *others="ALL", const char *options="DEFAULT", const char *outFlg="ALL", UShort_t rebin=5, Float_t eff=0, const char *base="")
Definition: DrawdNdeta.C:3603
TH2F
Definition: External.C:236
ScaleBinBySize
static void ScaleBinBySize(TH1D *h)
Scale histogram bins by its size.
Definition: PlotUtils.C:285
kSumw2
Bool_t kSumw2
polinomyal type for residual background under the peak
Definition: InvMassFit.C:45
FunctionGaussPol0
static Double_t FunctionGaussPol0(Double_t *x, Double_t *par)
Definition: PlotUtils.C:84
mass
Double_t mass
Definition: charmCutsOptimization.C:35
kPolN
Int_t kPolN
use mixed event to constrain combinatorial background
Definition: InvMassFit.C:44
FunctionGaussPol1
static Double_t FunctionGaussPol1(Double_t *x, Double_t *par)
Definition: PlotUtils.C:72
kFileName
TString kFileName
Input production directory name where input file is located.
Definition: InvMassFit.C:50
modStyleIndex
Int_t modStyleIndex[]
Definition: InvMassFit.C:70
PlotUtils.C
Fit invariant mass distributions.
ProjectAndFit
void ProjectAndFit(const Int_t nPt, TString comment, TString leg, TString hname, TArrayD xPtLimits, TArrayD xPt, TArrayD exPt, TH2F *hRe, TH2F *hMi)
Definition: InvMassFit.C:377
kParticle
TString kParticle
Calorimeter, EMCAL, DCAL (PHOS)
Definition: InvMassFit.C:53
ptMin
Double_t ptMin
Definition: MuonTrackingEfficiency.C:52
kNPairCut
Float_t kNPairCut
Apply Root method Sumw2(), off for pT hard prod where already done before.
Definition: InvMassFit.C:46
mesonMass
Double_t mesonMass
Definition: CalculateAveragePt.C:60
Efficiency
void Efficiency(Int_t nPt, TArrayD xPt, TArrayD exPt, TArrayD mesonPt, TArrayD emesonPt, TArrayD mesonPtBC, TArrayD emesonPtBC, TString hname)
Efficiency calculation called in ProjectAndFit()
Definition: InvMassFit.C:193
kCalorimeter
TString kCalorimeter
Automatic plots format.
Definition: InvMassFit.C:52
calorimeter
const TString calorimeter
Definition: anaM.C:36
tPolPi0
TF1 * tPolPi0
Definition: InvMassFit.C:66
DivideGraphs
static TGraphErrors * DivideGraphs(TGraphErrors *gNum, TGraphErrors *gDen)
Definition: PlotUtils.C:151
fil
TFile * fil
Definition: InvMassFit.C:60
kMix
Bool_t kMix
Definition: InvMassFit.C:43
GetFractionError
static Double_t GetFractionError(Double_t num, Double_t den, Double_t numErr, Double_t denErr)
Definition: PlotUtils.C:132
kHistoStartName
TString kHistoStartName
First Calo SM covered by a TRD SM, 6 in 2011 and 4 in 2012-13.
Definition: InvMassFit.C:48
lis
TList * lis
Definition: InvMassFit.C:62
kRebin
Int_t kRebin
Maximum value of chi2/ndf to define a fit as good (set lower than 1e6)
Definition: InvMassFit.C:56
Int_t
int Int_t
Definition: External.C:63
GetFileAndEvents
Bool_t GetFileAndEvents(TString dirName, TString listName)
Definition: InvMassFit.C:76
SetFitFun
void SetFitFun()
Definition: InvMassFit.C:113
Float_t
float Float_t
Definition: External.C:68
direc
TDirectoryFile * direc
Definition: InvMassFit.C:64
kTrunMixFunc
Bool_t kTrunMixFunc
Particle searched: "Pi0", "Eta".
Definition: InvMassFit.C:54
TH2D
Definition: External.C:228
InvMassFit
void InvMassFit(TString prodname="LHC18c3_NystromOn", TString filename="AnalysisResults", TString histoDir="Pi0IM_GammaTrackCorr_EMCAL", TString histoList="default", TString histoStart="", TString calorimeter="EMCAL", TString particle="Pi0", Float_t nPairMin=20, Bool_t sumw2=kTRUE, Int_t rebin=1, Int_t firstTRD=6, Bool_t mixed=kTRUE, Bool_t truncmix=kFALSE, Int_t pol=1, Bool_t drawAll=kFALSE, Bool_t drawMCAll=kFALSE, Bool_t drawPerSM=kTRUE, Bool_t drawPerSMGr=kFALSE, Bool_t drawPerSector=kFALSE, Bool_t drawPerSide=kFALSE, TString plotFormat="eps")
Definition: InvMassFit.C:2078
TH1D
Definition: External.C:212
fout
TFile * fout
Definition: InvMassFit.C:61
FunctionGaussPol3
static Double_t FunctionGaussPol3(Double_t *x, Double_t *par)
Definition: PlotUtils.C:48
GetRangeIntegralAndError
static void GetRangeIntegralAndError(TH1D *h, Int_t binMin, Int_t binMax, Double_t &integral, Double_t &integralErr)
Definition: PlotUtils.C:208
PlotGraphs
void PlotGraphs(TString opt, Int_t first, Int_t last)
Definition: InvMassFit.C:1054
PrimarySpectra
void PrimarySpectra(Int_t nPt, TArrayD xPtLimits, TArrayD xPt, TArrayD exPt)
Definition: InvMassFit.C:1834
nEvt
Double_t nEvt
Rebin invariant mass histograms default binning with this value.
Definition: InvMassFit.C:59
FunctionGaussPol2
static Double_t FunctionGaussPol2(Double_t *x, Double_t *par)
Definition: PlotUtils.C:60
FunctionTruncatedPolEta
static Double_t FunctionTruncatedPolEta(Double_t *x, Double_t *par)
Definition: PlotUtils.C:96
kProdName
TString kProdName
Common starting name in histograms.
Definition: InvMassFit.C:49
kChi2NDFMax
Float_t kChi2NDFMax
Use a truncated function to get the mixed event scale.
Definition: InvMassFit.C:55
modColorIndex
Int_t modColorIndex[]
Definition: InvMassFit.C:69
rebin
Int_t rebin
Definition: charmCutsOptimization.C:32
fitfun
TF1 * fitfun
Definition: InvMassFit.C:63
ymin
const Double_t ymin
Definition: AddTaskCFDStar.C:6
Bool_t
bool Bool_t
Definition: External.C:53
tPolEta
TF1 * tPolEta
Definition: InvMassFit.C:67
ptMax
Double_t ptMax
Definition: MuonTrackingEfficiency.C:53
kPlotFormat
TString kPlotFormat
Name of file with input histograms.
Definition: InvMassFit.C:51