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