AliPhysics  f6e6b3f (f6e6b3f)
InvMassFit.C
Go to the documentation of this file.
1 
16 #if !defined(__CINT__) || defined(__MAKECINT__)
17 
18 #include "TString.h"
19 #include "TH2F.h"
20 #include "TH1F.h"
21 #include "TH3F.h"
22 #include "TH1D.h"
23 #include "TF1.h"
24 #include "TMath.h"
25 #include "TCanvas.h"
26 #include "TStyle.h"
27 #include "TPad.h"
28 #include "TFile.h"
29 #include "TLegend.h"
30 #include "TObject.h"
31 #include "TDirectoryFile.h"
32 #include "TGraphErrors.h"
33 #include "TList.h"
34 #include <TGaxis.h>
35 
36 #endif
37 
38 //-----------------------------------------------------------------------------
39 Double_t pi0massP0(Double_t *x, Double_t *par);
40 Double_t pi0massP1(Double_t *x, Double_t *par);
41 Double_t pi0massP2(Double_t *x, Double_t *par);
42 Double_t pi0massP3(Double_t *x, Double_t *par);
43 Double_t truncatedPolPi0(Double_t *x, Double_t *par);
44 Double_t truncatedPolEta(Double_t *x, Double_t *par);
45 Double_t CrystalBall(Double_t *x, Double_t *par);
46 
47 Bool_t GetFileAndEvents( TString prodname, TString filename, TString dirName, TString listName );
48 void SetFitFun();
49 
50 // Settings
51 Bool_t mix = kFALSE;
52 TString part = "Eta";
53 Int_t polN = 1;
54 Bool_t sumw2 = kTRUE;
55 Bool_t drawAllCombi = 0;
56 Float_t nPairCut = 20;
57 
58 // Initializations
59 Int_t nEvt = 0;
60 TFile * fil = 0;
61 TList * lis = 0;
62 TDirectoryFile * direc =0;
63 TF1 *fitfun = 0;
64 
80 //-----------------------------------------------------------------------------
81 void InvMassFit
82 (TString prodname = "LHC18c3_NystromOn",//"LHC17l3b_fast",
83  TString filename = "AnalysisResults",
84  TString histoDir = "Pi0IM_GammaTrackCorr_EMCAL",
85  TString histoList = "default",
86  TString calorimeter = "DCAL",
87  TString particle = "Pi0",
88  Bool_t mixed = 0,
89  Int_t pol = 1,
90  Int_t ncomb = 1,
91  Float_t nPairMin = 20,
92  TString fileFormat = "eps")
93 {
94  part = particle;
95  polN = pol;
96  mix = mixed;
97 
98  Int_t ok = GetFileAndEvents(prodname,filename,histoDir, histoList);
99 
100  nPairCut = nPairMin;
101 // nPairCut/=nEvt;
102 
103  printf("Settings: prodname %s, filename %s, histoDir %s, histoList %s, particle %s, calorimeter %s, \n"
104  " \t mix %d, polN %d, n combi %d, n pairs cut %2.2e\n",
105  prodname.Data(),filename.Data(),histoDir.Data(),histoList.Data(), part.Data(), calorimeter.Data(),
106  mix, polN, ncomb,nPairCut);
107 
108  if ( !ok )
109  {
110  printf("Could not recover file <%p> or dir <%p> or list <%p>\n",fil,direc,lis);
111  return ;
112  }
113 
114  //---------------------------------------
115  // Set-up fitting and histogram plotting
116  //---------------------------------------
117  const Int_t ncombFix = 100;
118  const Int_t nFixBins = 50; // to initialize arrays size
119 
120  const Int_t nPtEta = 8;
121  Double_t xPtLimitsEta[] = {2.,4.,6.,8.,10.,14.,18.,24.};
122 
123  const Int_t nPtPi0 = 12;
124  Double_t xPtLimitsPi0[] = {2.,3.,4.,5.,6.,7.,8.,9.,10.,12.,14.,16.,18.};
125 
126  Double_t xPtLimits[nFixBins];
127  Double_t xPt [nFixBins] ;
128  Double_t exPt [nFixBins] ;
129 
130  Int_t nPt = nPtPi0;
131  if(part == "Eta") nPt = nPtEta;
132 
133  for(Int_t i = 0; i < nPt+1; i++)
134  {
135  if(part == "Pi0")
136  {
137  xPtLimits[i] = xPtLimitsPi0[i];
138  exPt[i] = (xPtLimitsPi0[i+1]-xPtLimitsPi0[i])/2;
139  xPt [i] = xPtLimitsPi0[i]+exPt[i];
140  }
141  else
142  {
143  xPtLimits[i] = xPtLimitsEta[i];
144  exPt[i] = (xPtLimitsEta[i+1]-xPtLimitsEta[i])/2;
145  xPt [i] = xPtLimitsEta[i]+exPt[i];
146  }
147  //printf("%s pT bin %d, pT %f, dpT %f\n",part.Data(),i,xPt[i],exPt[i]);
148  }
149 
150  Double_t mmin = 0;
151  Double_t mmax = 1;
152 
153  if(part=="Pi0")
154  {
155  mmin = 0.00;
156  mmax = 0.45;
157  }
158  else // eta
159  {
160  mmin = 0.25;
161  mmax = 0.95;
162  }
163 
164  Int_t rebin = 1 ;
165  Int_t rebin2 = 2*rebin;
166 
167  // - SM SM-Sector - Side Side
168  Int_t modColorIndex[]={1 , 2, 2, 3, 3, 4, 4, 7, 7, 6, 6, 2, 3, 4, 7, 6, 2, 2, 3, 3, 4, 4, 6, 6};
169  Int_t modStyleIndex[]={20,24,25,24,25,24,25,24,25,24,25,21,21,21,21,21,22,26,22,26,22,26,22,26};
170 
171  //Mix fit func
172 // TF1 *tPolPi0 = new TF1("truncatedPolPi0", truncatedPolPi0, 0.02, 0.25, 2);
173 // TF1 *tPolEta = new TF1("truncatedPolEta", truncatedPolEta, 0.2, 1, 2);
174 
175  // Fitting function initialization
176  SetFitFun();
177 
178  //---------------------------------------
179  // Get the histograms from file
180  //---------------------------------------
181 
182  TH2F * hRe [ncombFix];
183  TH2F * hPu ;
184  TH2F * hMi [ncombFix];
185 
186  TString comment[ncombFix];
187  TString leg [ncombFix];
188  TString hname [ncombFix];
189 
190  Int_t icalo = 0;
191  if(calorimeter=="DCAL") icalo = 1;
192 
193  if ( !lis )
194  {
195  hRe[0] = (TH2F *) fil->Get(Form("AnaPi0_Calo%d_hRe_cen0_pidbit0_asy0_dist1",icalo));
196  hMi[0] = (TH2F *) fil->Get(Form("AnaPi0_Calo%d_hMi_cen0_pidbit0_asy0_dist1",icalo));
197  }
198  else
199  {
200  hRe[0] = (TH2F *) lis->FindObject(Form("AnaPi0_Calo%d_hRe_cen0_pidbit0_asy0_dist1",icalo));
201  hMi[0] = (TH2F *) lis->FindObject(Form("AnaPi0_Calo%d_hMi_cen0_pidbit0_asy0_dist1",icalo));
202  }
203 
204  //hPu = (TH2F *) fil->Get("AnaPi0_TM1_MCOrgMass_2" );
205 
206  printf("histo Re %p - Mix %p\n",hRe[0],hMi[0]);
207 
208  comment[0] = "All SM";
209  leg [0] = "All SM";
210  hname [0] = "AllSM";
211 
212  //hRe[0]->Scale(1./nEvt);
213 
214  // Histograms for cluster pairs in different combinations of SM
215  //
216 
217  // Pairs in same SM
218  //
219  if(drawAllCombi)
220  {
221  Int_t imod;
222  Int_t firstMod = 0;
223  Int_t lastMod = 11;
224  if(calorimeter=="DCAL")
225  {
226  Int_t firstMod = 12;
227  Int_t lastMod = 19;
228  }
229 
230  for(imod = firstMod; imod<lastMod; imod++)
231  {
232  printf("imod+1 %d\n",imod+1);
233  if(!lis)
234  {
235  hRe[imod+1] = (TH2F *) fil->Get(Form("AnaPi0_Calo%d_hReMod_%d",icalo,imod));
236  hMi[imod+1] = (TH2F *) fil->Get(Form("AnaPi0_Calo%d_hMiMod_%d",icalo,imod));
237  }
238  else
239  {
240  hRe[imod+1] = (TH2F *) lis->FindObject(Form("AnaPi0_Calo%d_hReMod_%d",icalo,imod));
241  hMi[imod+1] = (TH2F *) lis->FindObject(Form("AnaPi0_Calo%d_hMiMod_%d",icalo,imod));
242  }
243 
244  comment[imod+1] = Form("both clusters in SM %d",imod);
245  leg [imod+1] = Form("SM %d",imod);
246  hname [imod+1] = Form("SM%d" ,imod);
247 
248  printf("histo mod %d Re %p - Mix %p\n",imod, hRe[imod],hMi[imod]);
249  }
250 
251  // Pairs in same sector
252  //
253  Int_t isector;
254  for(isector = 0; isector<5; isector++)
255  {
256  printf("imod+1+isector %d\n",imod+1+isector);
257 
258  if(!lis)
259  {
260  hRe[imod+1+isector] = (TH2F *) fil->Get(Form("AnaPi0_Calo%d_hReSameSectorEMCALMod_%d",icalo,isector));
261  hMi[imod+1+isector] = (TH2F *) fil->Get(Form("AnaPi0_Calo%d_hMiSameSectorEMCALMod_%d",icalo,isector));
262  }
263  else
264  {
265  hRe[imod+1+isector] = (TH2F *) lis->FindObject(Form("AnaPi0_Calo%d_hReSameSectorEMCALMod_%d",icalo,isector));
266  hMi[imod+1+isector] = (TH2F *) lis->FindObject(Form("AnaPi0_Calo%d_hMiSameSectorEMCALMod_%d",icalo,isector));
267  }
268 
269  comment[imod+1+isector] = Form("both clusters in Sector %d",isector);
270  leg [imod+1+isector] = Form("Sector %d",isector);
271  hname [imod+1+isector] = Form("Sector%d" ,isector);
272 
273  printf("histo sector %d Re %p - Mix %p\n",isector, hRe[imod+1+isector],hMi[imod+1+isector]);
274  }
275 
276  // Pairs in same side
277  //
278  Int_t iside;
279  for(iside = 0; iside < 8; iside++)
280  {
281  printf("imod+1+isector+iside %d\n",imod+1+isector+iside);
282 
283  if(!lis)
284  {
285  hRe[imod+1+isector+iside] = (TH2F *) fil->Get(Form("AnaPi0_Calo%d_hReSameSideEMCALMod_%d",icalo,iside));
286  hMi[imod+1+isector+iside] = (TH2F *) fil->Get(Form("AnaPi0_Calo%d_hMiSameSideEMCALMod_%d",icalo,iside));
287  }
288  else
289  {
290  hRe[imod+1+isector+iside] = (TH2F *) lis->FindObject(Form("AnaPi0_Calo%d_hReSameSideEMCALMod_%d",icalo,iside));
291  hMi[imod+1+isector+iside] = (TH2F *) lis->FindObject(Form("AnaPi0_Calo%d_hMiSameSideEMCALMod_%d",icalo,iside));
292  }
293  comment[imod+1+isector+iside] = Form("both clusters in Side %d",iside);
294  leg [imod+1+isector+iside] = Form("Side %d",iside);
295  hname [imod+1+isector+iside] = Form("Side%d" ,iside);
296 
297  printf("histo side %d Re %p - Mix %p\n",isector, hRe[imod+1+isector+iside],hMi[imod+1+isector+iside]);
298  }
299 
300  for(Int_t icomb = 0; icomb<ncomb; icomb++)
301  printf("icomb %d p %p ; %p ; %s, %s, %s\n",icomb,hRe[icomb],hMi[icomb],comment[icomb].Data(),leg[icomb].Data(),hname[icomb].Data());
302 
303  }
304 
305  //---------------------------------------
306  // Init some histograms array size
307  //---------------------------------------
308 
309  TLegend * pLegendIM[nFixBins];
310 
311  TH1D* hIM [ncombFix][nFixBins];
312  TH1D* hIMPu [nFixBins];
313  TH1D* hMix [ncombFix][nFixBins];
314  TH1D* hMixCorrected [ncombFix][nFixBins];
315  TH1D* hRatio [ncombFix][nFixBins];
316  TH1D* hSignal [ncombFix][nFixBins];
317 
318  Double_t mesonPt [ncombFix][nFixBins];
319  Double_t mesonPtBC [ncombFix][nFixBins];
320  Double_t mesonPtBCPu [nFixBins];
321  Double_t mesonMass [ncombFix][nFixBins];
322  Double_t mesonWidth [ncombFix][nFixBins];
323  Double_t mesonEffPt [ncombFix][nFixBins];
324  Double_t primMeson[nFixBins];
325 
326  Double_t emesonPt [ncombFix][nFixBins];
327  Double_t emesonPtBC [ncombFix][nFixBins];
328  Double_t emesonPtBCPu [nFixBins];
329  Double_t emesonMass [ncombFix][nFixBins];
330  Double_t emesonWidth[ncombFix][nFixBins];
331  Double_t emesonEffPt[ncombFix][nFixBins];
332  Double_t eprimMeson[nFixBins];
333 
334  // Mix correction factor
335  //TF1 *line3[ncombFix][nFixBins];// = new TF1("line3", "[0]+[1]*x+[2]*x*x+[3]*x*x*x", 0.0, 1);
336 
337 
338  for(Int_t icomb = 0; icomb< ncomb ; icomb++)
339  {
340  for(Int_t ipt = 0; ipt< nPt; ipt++)
341  {
342  hIM [icomb][ipt] = 0 ;
343  hIMPu [ipt] = 0 ;
344  hMix [icomb][ipt] = 0 ;
345  hMixCorrected[icomb][ipt] = 0 ;
346  hRatio [icomb][ipt] = 0 ;
347  hSignal [icomb][ipt] = 0 ;
348 
349  mesonPt [icomb][ipt] = 0 ;
350  mesonPtBC [icomb][ipt] = 0 ;
351  mesonPtBCPu [ipt] = 0 ;
352  mesonMass [icomb][ipt] = 0 ;
353  mesonWidth [icomb][ipt] = 0 ;
354  mesonEffPt [icomb][ipt] = 0 ;
355 
356  emesonPt [icomb][ipt] = 0 ;
357  emesonPtBC [icomb][ipt] = 0 ;
358  emesonPtBCPu [ipt] = 0 ;
359  emesonMass [icomb][ipt] = 0 ;
360  emesonWidth [icomb][ipt] = 0 ;
361  emesonEffPt [icomb][ipt] = 0 ;
362 
363  //line3 [icomb][ipt] = 0;
364  }
365  }
366 
367  //-------------------------------------
368  // Do energy projections, fits and plotting
369  //-------------------------------------
370 
371 
372  // Invariant mass
373  //
374  //gROOT->Macro("$MACROS/style.C");//Set different root style parameters
375  gStyle->SetOptTitle(1);
376  gStyle->SetTitleOffset(2,"Y");
377  gStyle->SetOptStat(0);
378  gStyle->SetOptFit(000000);
379 
380  Int_t col = TMath::Ceil(TMath::Sqrt(nPt));
381 
382  for(Int_t icomb = 0; icomb< ncomb; icomb++)
383  {
384  TCanvas * cIMModi = new TCanvas(Form("Combination_%d\n", icomb), Form("Combination_%d\n", icomb), 1200, 1200) ;
385  cIMModi->Divide(col, col);
386 
387  for(Int_t i = 0; i < nPt; i++)
388  {
389  if(icomb >10 && i > 8) continue;
390 
391  cIMModi->cd(i+1) ;
392  //gPad->SetLogy();
393  Double_t ptMin = xPtLimits[i];
394  Double_t ptMax = xPtLimits[i+1];
395  printf("Bin %d (%f, %f)\n",i, ptMin,ptMax);
396 
397  hIM[icomb][i] = hRe[icomb]->ProjectionY(Form("IM_Comb%d_PtBin%d",icomb,i),
398  hRe[icomb]->GetXaxis()->FindBin(ptMin),
399  hRe[icomb]->GetXaxis()->FindBin(ptMax));
400  hIM[icomb][i]->SetTitle(Form("%2.1f < p_{T, #gamma#gamma} < %2.1f GeV/c",ptMin,ptMax));
401  if(i < nPt-3)
402  hIM[icomb][i]->Rebin(rebin);
403  else
404  hIM[icomb][i]->Rebin(rebin2);
405 
406  if(sumw2)
407  hIM[icomb][i]->Sumw2();
408 
409  hIM[icomb][i]->SetXTitle("M_{#gamma,#gamma} (GeV/c^{2})");
410  hIM[icomb][i]->SetLineColor(modColorIndex[icomb]);
411  //printf("mmin %f, mmax %f\n",mmin,mmax);
412  hIM[icomb][i]->SetAxisRange(mmin,mmax,"X");
413  hIM[icomb][i] ->SetLineWidth(2);
414  hIM[icomb][i] ->SetLineColor(4);
415  //hIM[icomb][i]->Draw();
416 
417 
418 // if(filename.Contains("imu") && icomb==0)
419 // {
420 // hIMPu[i] = hPu->ProjectionY(Form("IMPu_Comb%d_PtBin%d",icomb,i),
421 // hPu->GetXaxis()->FindBin(ptMin),
422 // hPu->GetXaxis()->FindBin(ptMax));
423 // hIMPu[i]->SetTitle(Form("%2.1f < p_{T, #gamma#gamma} < %2.1f GeV/c",ptMin,ptMax));
424 // if(i < nPt-3)
425 // hIMPu[i]->Rebin(rebin);
426 // else
427 // hIMPu[i]->Rebin(rebin2);
428 //
429 // if(sumw2)
430 // hIMPu[icomb][i]->Sumw2();
431 // else
432 // hIMPu[icomb][i]->Scale(1.e8);
433 //
434 // hIMPu[i]->SetXTitle("M_{#gamma,#gamma} (GeV/c^{2})");
435 // hIMPu[i]->SetLineColor(modColorIndex[icomb]);
436 // //printf("mmin %f, mmax %f\n",mmin,mmax);
437 // hIMPu[i]->SetAxisRange(mmin,mmax,"X");
438 // hIMPu[i]->SetLineWidth(2);
439 // hIMPu[i]->SetLineColor(kOrange-3);
440 // }
441 
442  Double_t nMax = 0;
443  if(part=="Pi0")
444  {
445  nMax= hIM[icomb][i]->Integral(hIM[icomb][i]->FindBin(0.09),
446  hIM[icomb][i]->FindBin(0.25));
447  }
448  else
449  {
450  nMax = hIM[icomb][i]->Integral(hIM[icomb][i]->FindBin(0.4),
451  hIM[icomb][i]->FindBin(0.65));
452  }
453 
454  //printf("icombi %d, bin %d, nMax %f \n",icomb,i,nMax);
455 
456  fitfun->SetParLimits(0,nMax/100,nMax*100);
457 
458  if(nMax > nPairCut && !mix)
459  {
460  if(part=="Pi0")
461  {
462  fitfun->SetParameters(nMax/4,0.135,0.02,1.,0.);
463  if(i<9)
464  hIM[icomb][i]->Fit("fitfun","QR","",0.09,0.25);
465  else
466  hIM[icomb][i]->Fit("fitfun","QR","",0.10,0.25);
467  } // pi0
468  else //eta
469  {
470  fitfun->SetParameters(nMax/4,0.54,0.04,1.,0.);
471  if(i<4)
472  hIM[icomb][i]->Fit("fitfun","QR","",0.4,0.7);
473  else if(i < 15)
474  hIM[icomb][i]->Fit("fitfun","QR","",0.3,0.8);
475  else
476  hIM[icomb][i]->Fit("fitfun","QR","",0.3,0.8);
477  }
478  }
479 
480  if(part=="Pi0")
481  {
482  if(i < nPt -3)
483  hIM[icomb][i]->SetMaximum(nMax/5.);
484  else
485  hIM[icomb][i]->SetMaximum(nMax/2.);
486 
487  }
488  else
489  {
490  if(i < nPt -3)
491  hIM[icomb][i]->SetMaximum(nMax/20.);
492  else
493  hIM[icomb][i]->SetMaximum(nMax/8.);
494 
495  }
496 
497  Double_t mStep = hIM[icomb][i]->GetBinWidth(i);
498  hIM[icomb][i]->SetMinimum(0.1);
499  hIM[icomb][i]->SetAxisRange(mmin,mmax,"X");
500 
501  pLegendIM[i] = 0;
502 
503  if ( hIM[icomb][i]->GetFunction("fitfun") && !mix )
504  {
505  if(hIM[icomb][i]->GetFunction("fitfun")->GetChisquare()/hIM[icomb][i]->GetFunction("fitfun")->GetNDF()<100)
506  {
507  Double_t A = hIM[icomb][i]->GetFunction("fitfun")->GetParameter(0);
508  Double_t mean = hIM[icomb][i]->GetFunction("fitfun")->GetParameter(1);
509  Double_t sigm = hIM[icomb][i]->GetFunction("fitfun")->GetParameter(2);
510 // Double_t a0 = hIM[icomb][i]->GetFunction("fitfun")->GetParameter(3);
511 // Double_t a1 = hIM[icomb][i]->GetFunction("fitfun")->GetParameter(4);
512 // Double_t a2 = 0;
513 // if (polN == 2)
514 // a2 = hIM[icomb][i]->GetFunction("fitfun")->GetParameter(5);
515 
516  Double_t eA = hIM[icomb][i]->GetFunction("fitfun")->GetParError(0);
517  Double_t emean = hIM[icomb][i]->GetFunction("fitfun")->GetParError(1);
518  Double_t esigm = hIM[icomb][i]->GetFunction("fitfun")->GetParError(2);
519 // Double_t ea0 = hIM[icomb][i]->GetFunction("fitfun")->GetParError(3);
520 // Double_t ea1 = hIM[icomb][i]->GetFunction("fitfun")->GetParError(4);
521 // Double_t ea2 = 0;
522 // if (polN == 2)
523 // ea2 = hIM[icomb][i]->GetFunction("fitfun")->GetParError(5);
524 
525  pLegendIM[i] = new TLegend(0.55,0.65,0.95,0.93);
526  pLegendIM[i]->SetTextSize(0.035);
527  pLegendIM[i]->SetFillColor(10);
528  pLegendIM[i]->SetBorderSize(1);
529  pLegendIM[i]->SetHeader(Form("#pi %s",leg[icomb].Data()));
530  pLegendIM[i]->AddEntry("",Form("A = %2.2f #pm %2.2f ",A,eA),"");
531  pLegendIM[i]->AddEntry("",Form("#mu = %2.3f #pm %2.3f ",mean,emean),"");
532  pLegendIM[i]->AddEntry("",Form("#sigma = %2.3f #pm %2.3f ",sigm,esigm),"");
533  //pLegendIM[i]->AddEntry("",Form("p_{0} = %2.1f #pm %2.1f ",a0,ea0),"");
534  //pLegendIM[i]->AddEntry("",Form("p_{1} = %2.1f #pm %2.1f ",a1,ea1),"");
535  //if(polN==2)pLegendIM[i]->AddEntry("",Form("p_{2} = %2.1f#pm %2.1f ",a2,ea2),"");
536 
537  mesonPt[icomb] [i] = A*sigm / mStep * TMath::Sqrt(TMath::TwoPi());
538  Double_t ePi0 =
539  TMath::Power(eA/A,2) +
540  TMath::Power(esigm/sigm,2);
541  emesonPt[icomb][i]= TMath::Sqrt(ePi0)* mesonPt[icomb][i];
542 
543  emesonPt[icomb][i] = TMath::Min(emesonPt[icomb][i], TMath::Sqrt(mesonPt[icomb] [i]));
544 
545  mesonPt [icomb] [i]/=(nEvt*(exPt[i]*2));
546  emesonPt[icomb] [i]/=(nEvt*(exPt[i]*2));
547 
548  cout << "N(pi0) fit = " << mesonPt[icomb][i] << "+-"<< emesonPt[icomb][i] << " mstep " << mStep<<endl;
549 
550  mesonMass [icomb][i] = mean*1000.;
551  emesonMass[icomb][i] = emean*1000.;
552 
553  mesonWidth [icomb][i] = sigm*1000.;
554  emesonWidth[icomb][i] = esigm*1000.;
555 
556  } //Good fit
557  else{
558 
559  mesonPt[icomb][i] =-1;
560  emesonPt[icomb][i] = -1;
561 
562  mesonMass[icomb][i] = -1;
563  emesonMass[icomb][i] = -1;
564 
565  mesonWidth[icomb][i] = -1;
566  emesonWidth[icomb][i] = -1;
567  }
568  }
569  else{
570  mesonPt[icomb][i] =-1;
571  emesonPt[icomb][i] = -1;
572 
573  mesonMass[icomb][i] = -1;
574  emesonMass[icomb][i] = -1;
575 
576  mesonWidth[icomb][i] = -1;
577  emesonWidth[icomb][i] = -1;
578  }
579 
580 
581  //--------------------------------------------------
582  // Mix
583  //--------------------------------------------------
584  if(mix)
585  {
586  hMix[icomb][i] = (TH1D*) hMi[icomb]->ProjectionY(Form("MiMass_PtBin%d_Combi%d",i,icomb),
587  hMi[icomb]->GetXaxis()->FindBin(ptMin),
588  hMi[icomb]->GetXaxis()->FindBin(ptMax));
589  hMix[icomb][i]->SetLineColor(1);
590  hMix[icomb][i]->SetTitle(Form("%2.2f < p_{T, #gamma#gamma} < %2.2f GeV/c",ptMin,ptMax));
591  if(i < nPt-3)
592  hMix[icomb][i]->Rebin(rebin);
593  else
594  hMix[icomb][i]->Rebin(rebin2);
595 
596 
597 
598  // Double_t sptmin = 0.2;
599  // Double_t sptmax = 0.45;
600 
601  // Double_t scalereal = hIM[icomb]->Integral(hIM[icomb]->FindBin(sptmin),
602  // hIM[icomb]->FindBin(sptmax));
603  // Double_t scalemix = hMix[icomb]->Integral(hMix[icomb]->FindBin(sptmin),
604  // hMix[icomb]->FindBin(sptmax));
605  // if(scalemix > 0)
606  // hMix[icomb]->Scale(scalereal/scalemix);
607  // else{
608  // //printf("could not scale: re %f mi %f\n",scalereal,scalemix);
609  // continue;
610  // }
611 
612  // hMix[icomb][i]->SetLineWidth(1);
613 
614  if(!sumw2)
615  hMix[icomb][i]->Sumw2();
616 
617  // --- Ratio ---
618  hRatio[icomb][i] = (TH1D*)hIM[icomb][i]->Clone(Form("RatioRealMix_PtBin%d_comb%d",i,icomb));
619  hRatio[icomb][i]->SetAxisRange(mmin,mmax,"X");
620  hRatio[icomb][i]->Divide(hMix[icomb][i]);
621  Double_t p0 =0;
622  Double_t p1 =0;
623  Double_t p2 =0;
624  Double_t p3 =0;
625 
626  if(hRatio[icomb][i]->GetEntries()>100)
627  {
628 // if(part=="Pi0")
629 // {
630 // hRatio[icomb][i]->Fit("truncatedPolPi0", "NQRL","",0.11,0.4);
631 // p0= tPolPi0->GetParameter(0);
632 // p1= tPolPi0->GetParameter(1);
633 // p2= tPolPi0->GetParameter(2);
634 // p3= tPolPi0->GetParameter(3);
635 // } // pi0
636 // else // eta
637 // {
638 // hRatio[icomb][i]->SetAxisRange(mmin,mmax);
639 // hRatio[icomb][i]->Fit("truncatedPolEta", "NQRL","",0.25,0.95);
640 // p0 = tPolEta->GetParameter(0);
641 // p1 = tPolEta->GetParameter(1);
642 // p2 = tPolEta->GetParameter(2);
643 // p3 = tPolEta->GetParameter(3);
644 // }
645 //
646 // //hRatio[icomb][i]->Draw("same");
647 //
648 // //line3[icomb][i]= new TF1("line3", "[0]+[1]*x+[2]*x*x+[3]*x*x*x", mmin, mmax);
649 // line3[icomb][i]= new TF1("line3", "[0]+[1]*x", mmin, mmax);
650 // line3[icomb][i]->SetParameters(p0,p1,p2,p3);
651 //
652 // if(hMix[icomb][i])
653 // {
654 // //printf("Correct\n");
655 // hMixCorrected[icomb][i] = (TH1D*)hMix[icomb][i]->Clone(Form("MixCorrected_PtBin%d_comb%d",i,icomb));
656 // for(Int_t j = 0; j< hMix[icomb][i]->GetNbinsX(); j++){
657 // Double_t x = hMix[icomb][i]->GetBinCenter(j);
658 // Double_t correction = line3[icomb][i]->Eval(x);
659 // Double_t corrected = hMix[icomb][i]->GetBinContent(j)*correction;
660 // Double_t ecorrected = hMix[icomb][i]->GetBinError(j) *correction;
661 // //printf("bin %d, center %f, mix %f, corrected %f\n",i,x,hMix[icomb][i]->GetBinContent(j),corrected);
662 // hMixCorrected[icomb][i]->SetBinContent(j, corrected);
663 // hMixCorrected[icomb][i]->SetBinError (j,ecorrected);
664 // }
665 // }
666 
667  Float_t scale = hRatio[icomb][i]->GetBinContent(hRatio[icomb][i]->FindBin(0.9));
668  hMix[icomb][i]->Scale(scale);
669 
670  // --- Subtracted ---
671  //printf("Signal\n");
672  hSignal[icomb][i] = (TH1D*)hIM[icomb][i]->Clone(Form("Signal_PtBin%d_comb%d",i,icomb));
673  //hSignal[icomb][i] ->Add(hMixCorrected[icomb][i],-1);
674  hSignal[icomb][i] ->Add(hMix[icomb][i],-1);
675  hSignal[icomb][i] ->SetLineColor(kViolet);
676 
677  // ---- Fit subtracted signal
678  if(part=="Pi0")
679  {
680  nMax= hSignal[icomb][i]->Integral(hIM[icomb][i]->FindBin(0.1),
681  hIM[icomb][i]->FindBin(0.2));
682  }
683  else
684  {
685  nMax = hSignal[icomb][i]->Integral(hIM[icomb][i]->FindBin(0.4),
686  hIM[icomb][i]->FindBin(0.65));
687  }
688 
689  if(nMax > nPairCut)
690  {
691  fitfun->SetParLimits(0,nMax/100,nMax*100);
692 
693  if(part=="Pi0"){
694  if(polN < 4 )fitfun->SetParameters(nMax/5,0.135,20,0);
695  else fitfun->SetParameters(nMax/5,20,0.135,20,20,nMax/5);
696 
697  if(i < nPt-4)
698  hSignal[icomb][i]->Fit("fitfun","QR","",0.11,0.3);
699  else
700  hSignal[icomb][i]->Fit("fitfun","QR","",0.11,0.3);
701  }// pi0
702  else // eta
703  {
704 
705  if(polN < 4 )fitfun->SetParameters(nMax/5,0.54,40,0);
706  else fitfun->SetParameters(nMax/5,40,0.54,40,40,nMax/5);
707  //if(i<4)
708  hSignal[icomb][i]->Fit("fitfun","QR","",0.4,0.7);
709  //else if(i < 15)
710  // hSignal[icomb][i]->Fit("fitfun","QRL","",0.3,0.8);
711  //else
712  //hSignal[icomb][i]->Fit("fitfun","QRL","",0.3,0.8);
713  }
714  }
715 
716  if(hSignal[icomb][i]->GetFunction("fitfun"))
717  {
718 // printf("ipt %d: Chi/NDF %f, Chi %f, NDF %d\n",i,
719 // hSignal[icomb][i]->GetFunction("fitfun")->GetChisquare()/hSignal[icomb][i]->GetFunction("fitfun")->GetNDF(),
720 // hSignal[icomb][i]->GetFunction("fitfun")->GetChisquare(),
721 // hSignal[icomb][i]->GetFunction("fitfun")->GetNDF());
722  if(hSignal[icomb][i]->GetFunction("fitfun")->GetChisquare()/hSignal[icomb][i]->GetFunction("fitfun")->GetNDF()<1000){
723  Double_t A = hSignal[icomb][i]->GetFunction("fitfun")->GetParameter(0);
724  Double_t mean = hSignal[icomb][i]->GetFunction("fitfun")->GetParameter(1);
725  Double_t sigm = hSignal[icomb][i]->GetFunction("fitfun")->GetParameter(2);
726 // Double_t a0 = hSignal[icomb][i]->GetFunction("fitfun")->GetParameter(3);
727 // Double_t a1 = hSignal[icomb][i]->GetFunction("fitfun")->GetParameter(4);
728 // Double_t a2 = 0;
729 // if (polN == 2)
730 // a2 = hSignal[icomb][i]->GetFunction("fitfun")->GetParameter(5);
731 
732  Double_t eA = hSignal[icomb][i]->GetFunction("fitfun")->GetParError(0);
733  Double_t emean = hSignal[icomb][i]->GetFunction("fitfun")->GetParError(1);
734  Double_t esigm = hSignal[icomb][i]->GetFunction("fitfun")->GetParError(2);
735 // Double_t ea0 = hSignal[icomb][i]->GetFunction("fitfun")->GetParError(3);
736 // Double_t ea1 = hSignal[icomb][i]->GetFunction("fitfun")->GetParError(4);
737 // Double_t ea2 = 0;
738 // if (polN == 2)
739 // ea2 = hSignal[icomb][i]->GetFunction("fitfun")->GetParError(5);
740  pLegendIM[i] = new TLegend(0.55,0.65,0.95,0.93);
741  pLegendIM[i]->SetTextSize(0.035);
742  pLegendIM[i]->SetFillColor(10);
743  pLegendIM[i]->SetBorderSize(1);
744  pLegendIM[i]->SetHeader(Form("#pi %s",leg[icomb].Data()));
745  pLegendIM[i]->AddEntry("",Form("A = %2.2f #pm %2.2f ",A,eA),"");
746  pLegendIM[i]->AddEntry("",Form("#mu = %2.3f #pm %2.3f ",mean,emean),"");
747  pLegendIM[i]->AddEntry("",Form("#sigma = %2.3f #pm %2.3f ",sigm,esigm),"");
748  //pLegendIM[i]->AddEntry("",Form("p_{0} = %2.1f #pm %2.1f ",a0,ea0),"");
749  //pLegendIM[i]->AddEntry("",Form("p_{1} = %2.1f #pm %2.1f ",a1,ea1),"");
750  //if(polN==2)pLegendIM[i]->AddEntry("",Form("p_{2} = %2.1f#pm %2.1f ",a2,ea2),"");
751 
752 
753  mesonPt[icomb] [i] = A*sigm / mStep * TMath::Sqrt(TMath::TwoPi());
754  Double_t ePi0 =
755  TMath::Power(eA/A,2) +
756  TMath::Power(esigm/sigm,2);
757  emesonPt[icomb][i]= TMath::Sqrt(ePi0)* mesonPt[icomb][i];
758 
759  emesonPt[icomb][i] = TMath::Min(emesonPt[icomb][i], TMath::Sqrt(mesonPt[icomb] [i]));
760 
761  mesonPt [icomb] [i]/=(nEvt*(exPt[i]*2));
762  emesonPt[icomb] [i]/=(nEvt*(exPt[i]*2));
763 
764  //cout << "N(pi0) fit = " << mesonPt[icomb][i] << "+-"<< emesonPt[icomb][i] << " mstep "<<mStep<<endl;
765 
766  mesonMass [icomb][i] = mean *1000.;
767  emesonMass[icomb][i] = emean*1000.;
768 
769  mesonWidth [icomb][i] = sigm *1000.;
770  emesonWidth[icomb][i] = esigm*1000.;
771 
772  } //Good fit
773  else{
774 
775  mesonPt[icomb][i] =-1;
776  emesonPt[icomb][i] = -1;
777 
778  mesonMass[icomb][i] = -1;
779  emesonMass[icomb][i] = -1;
780 
781  mesonWidth[icomb][i] = -1;
782  emesonWidth[icomb][i] = -1;
783  }
784  }
785  else{
786  mesonPt[icomb][i] =-1;
787  emesonPt[icomb][i] = -1;
788 
789  mesonMass[icomb][i] = -1;
790  emesonMass[icomb][i] = -1;
791 
792  mesonWidth[icomb][i] = -1;
793  emesonWidth[icomb][i] = -1;
794  }
795 
796  Double_t mass = mesonMass [icomb][i];
797  Double_t width = mesonWidth[icomb][i];
798  Double_t mMinBin = mass - 2 * width;
799  Double_t mMaxBin = mass + 2 * width;
800  if(mass <= 0 || width <= 0)
801  {
802  if(part=="Pi0")
803  {
804  mMinBin = 0.115; mMaxBin = 0.3 ;
805  }
806  else{
807  mMinBin = 0.4; mMaxBin = 0.9 ;
808  }
809 
810  }
811 
812  if(hSignal[icomb][i])
813  {
814  mesonPtBC[icomb] [i] = hSignal[icomb][i]->Integral(hSignal[icomb][i]->FindBin(mMinBin), hSignal[icomb][i]->FindBin(mMaxBin)) ;
815 
816  emesonPtBC[icomb] [i] = TMath::Sqrt(mesonPtBC[icomb] [i]);
817 
818  mesonPtBC[icomb] [i]/=(nEvt*(exPt[i]*2));
819  emesonPtBC[icomb] [i]/=(nEvt*(exPt[i]*2));
820  //printf("Bin %d, [%1.1f,%1.1f]\n",i,xPt[i],xPt[i+1]);
821  }
822 
823 // if(filename.Contains("imu") && icomb==0)
824 // {
825 // mesonPtBCPu [i] = hIMPu[i]->Integral(hIMPu[i]->FindBin(mMinBin), hIMPu[i]->FindBin(mMaxBin)) ;
826 //
827 // emesonPtBCPu[i] = TMath::Sqrt(mesonPtBCPu [i]);
828 //
829 // mesonPtBCPu [i]/=(nEvt*(exPt[i]*2));
830 // emesonPtBCPu[i]/=(nEvt*(exPt[i]*2));
831 // //printf("Bin %d, [%1.1f,%1.1f]\n",i,xPt[i],xPt[i+1]);
832 //
833 // }
834 
835  printf("N(pi0) BC = %2.3e +- %2.4e\n",mesonPtBC[icomb][i],emesonPtBC[icomb][i]);
836  if(filename.Contains("imu") && icomb==0)
837  {
838  printf("N(pi0) Pu = %2.3e +- %2.4e\n",mesonPtBCPu[i],emesonPtBCPu[i]);
839  if(mesonPtBCPu[i]>0)printf("\t ratio BC / Pu %f \n",mesonPtBC[icomb][i]/mesonPtBCPu[i]);
840  }
841  }
842 
843 
844  if(nMax > nPairCut)
845  {
846  hIM[icomb][i]->Draw("HE");
847 
848  hMix[icomb][i]->Draw("same");
849  //if(hMixCorrected[icomb][i])hMixCorrected[icomb][i] ->Draw("same");
850 
851  if(hSignal[icomb][i]) hSignal[icomb][i] -> Draw("same");
852 
853  if(hSignal[icomb][i] && hSignal[icomb][i]->GetFunction("fitfun") && pLegendIM[i]) pLegendIM[i]->Draw();
854  }
855 
856  }//mix on
857  else if(nMax > nPairCut)
858  {
859  //printf("HERE 1\n");
860 
861  hIM[icomb][i]->Draw("HE");
862  //printf("HERE 2\n");
863 
864  //if(hIM[icomb][i]->GetFunction("fitfun") && pLegendIM[i])pLegendIM[i]->Draw();
865  //printf("HERE 3\n");
866 
867  }
868  } //pT bins
869 
870  cIMModi->Print(Form("IMfigures/%s_%s_Mgg_%s_%s_%s.%s",prodname.Data(),calorimeter.Data(),part.Data(),filename.Data(), hname[icomb].Data(),fileFormat.Data()));
871 
872  }// combinations
873 
874 
875  //xxxx Real / Mixxxxx
876  if(mix)
877  {
878  for(Int_t icomb = 0; icomb< ncomb; icomb++)
879  {
880  printf("Do real/mix\n");
881 
882  TCanvas * cRat = new TCanvas(Form("CombinationRatio_%d\n", icomb), Form("CombinationRatio_%d\n", icomb), 1200, 1200) ;
883  cRat->Divide(col, col);
884 
885  for(Int_t i = 0; i < nPt; i++){
886  cRat->cd(i+1) ;
887  //gPad->SetGridy();
888  //gPad->SetLog();
889  if(!hRatio[icomb][i]) {
890  //printf("No ratio in pt bin %d continue\n",i);
891  continue;
892  }
893 
894  Double_t nMax =0;
895  if(part=="Pi0")
896  nMax = hRatio[icomb][i]->Integral(hRatio[icomb][i]->FindBin(0.005),
897  hRatio[icomb][i]->FindBin(0.20));
898  else
899  nMax = hRatio[icomb][i]->Integral(hRatio[icomb][i]->FindBin(0.4),
900  hRatio[icomb][i]->FindBin(0.6));
901 
902  hRatio[icomb][i]->SetMaximum(nMax/4);
903  hRatio[icomb][i]->SetMinimum(1e-6);
904  hRatio[icomb][i]->SetAxisRange(mmin,mmax,"X");
905 
906  hRatio[icomb][i]->Draw();
907 
908  //if(line3[icomb][i]) line3[icomb][i]->Draw("same");
909  }
910 
911 
912  cRat->Print(Form("IMfigures/%s_%s_MggRatio_%s_%s_%s.%s",prodname.Data(),calorimeter.Data(),part.Data(),filename.Data(), hname[icomb].Data(),fileFormat.Data()));
913 
914  }
915  }
916 
917  //------------------------------
918  // Fit parameters final plotting
919  //------------------------------
920 
921  gStyle->SetOptTitle(0);
922 
923  // Put fit results in TGraphErrors
924  //
925  TGraphErrors* gPt [ncombFix];
926  TGraphErrors* gPtBC [ncombFix];
927  TGraphErrors* gMass [ncombFix];
928  TGraphErrors* gWidth[ncombFix];
929 
930  TString particleN = " #pi^{0}";
931  if(part=="Eta") particleN = " #eta";
932 
933  for(Int_t icomb = 0; icomb<ncomb; icomb++)
934  {
935  //cout<<"icomb "<<icomb<<" " <<gPt[icomb]<<" "<<comment[icomb]<<endl;
936  gPt[icomb] = new TGraphErrors(nPt,xPt,mesonPt[icomb],exPt,emesonPt[icomb]);
937  gPt[icomb] ->SetName(Form("gPt_%s",hname[icomb].Data()));
938  gPt[icomb] ->GetHistogram()->SetTitle(Form("p_{T} of reconstructed %s, %s ",particleN.Data(), comment[icomb].Data()));
939  gPt[icomb] ->GetHistogram()->SetXTitle("p_{T} (GeV/c) ");
940  gPt[icomb] ->GetHistogram()->SetYTitle(Form("dN_{%s}/dp_{T} (GeV/c)^{-1} / N_{events} ",particleN.Data()));
941  gPt[icomb] ->GetHistogram()->SetAxisRange(0.,30);
942  gPt[icomb] ->GetHistogram()->SetTitleOffset(1.5,"Y");
943  gPt[icomb] ->SetMarkerStyle(20);
944  gPt[icomb] ->SetMarkerSize(1);
945  gPt[icomb] ->SetMarkerColor(1);
946  // gPt[icomb] ->GetHistogram()->SetMaximum(1e8);
947  // gPt[icomb] ->GetHistogram()->SetMinimum(1e-8);
948 
949  gPtBC[icomb] = new TGraphErrors(nPt,xPt,mesonPtBC[icomb],exPt,emesonPtBC[icomb]);
950  gPtBC[icomb] ->SetName(Form("gPtBC_%s",hname[icomb].Data()));
951  gPtBC[icomb] ->GetHistogram()->SetTitle(Form("p_{T} of reconstructed %s, %s ",particleN.Data(), comment[icomb].Data()));
952  gPtBC[icomb] ->GetHistogram()->SetXTitle("p_{T} (GeV/c) ");
953  gPtBC[icomb] ->GetHistogram()->SetYTitle(Form("dN_{%s}/dp_{T} (GeV/c)^{-1} / N_{events} ",particleN.Data()));
954  gPtBC[icomb] ->GetHistogram()->SetAxisRange(0.,30);
955  gPtBC[icomb] ->GetHistogram()->SetTitleOffset(1.5,"Y");
956  gPtBC[icomb] ->SetMarkerStyle(20);
957  gPtBC[icomb] ->SetMarkerSize(1);
958  gPtBC[icomb] ->SetMarkerColor(1);
959 
960 
961  // gPtBC[icomb] ->GetHistogram()->SetMaximum(1e8);
962  // gPtBC[icomb] ->GetHistogram()->SetMinimum(1e-8);
963 
964  gMass[icomb] = new TGraphErrors(nPt,xPt,mesonMass[icomb],exPt,emesonMass[icomb]);
965  gMass[icomb] ->SetName(Form("gMass_%s",hname[icomb].Data()));
966  gMass[icomb] ->GetHistogram()->SetTitle(Form("mass vs p_{T} of reconstructed %s, %s ",particleN.Data(), comment[icomb].Data()));
967  gMass[icomb] ->GetHistogram()->SetXTitle("p_{T} (GeV/c) ");
968  gMass[icomb] ->GetHistogram()->SetYTitle(Form("mass_{%s} (MeV/c)^{2} ",particleN.Data()));
969  //gMass[icomb] ->GetHistogram()->SetAxisRange(0.,30);
970  gMass[icomb] ->GetHistogram()->SetTitleOffset(1.5,"Y");
971  gMass[icomb] ->SetMarkerStyle(20);
972  gMass[icomb] ->SetMarkerSize(1);
973  gMass[icomb] ->SetMarkerColor(1);
974 
975 
976  gWidth[icomb] = new TGraphErrors(nPt,xPt,mesonWidth[icomb],exPt,emesonWidth[icomb]);
977  gWidth[icomb] ->SetName(Form("gWidth_%s",hname[icomb].Data()));
978  gWidth[icomb] ->GetHistogram()->SetTitle(Form("Width vs p_{T} of reconstructed %s, %s ",particleN.Data(), comment[icomb].Data()));
979  gWidth[icomb] ->GetHistogram()->SetXTitle("p_{T} (GeV/c) ");
980  gWidth[icomb] ->GetHistogram()->SetYTitle(Form("#sigma_{%s} (MeV/c)^{2} ",particleN.Data()));
981  //gWidth[icomb] ->GetHistogram()->SetAxisRange(0.,30);
982  gWidth[icomb] ->GetHistogram()->SetTitleOffset(1.5,"Y");
983  gWidth[icomb] ->SetMarkerStyle(20);
984  gWidth[icomb] ->SetMarkerSize(1);
985  gWidth[icomb] ->SetMarkerColor(1);
986 
987  }
988 
989  TGraphErrors *gPtBCPu = 0;
990  if(filename.Contains("imu"))
991  {
992  gPtBCPu = new TGraphErrors(nPt,xPt,mesonPtBCPu,exPt,emesonPtBCPu);
993  gPtBCPu ->SetName(Form("gPtBC_%s_Pure",hname[0].Data()));
994  gPtBCPu ->GetHistogram()->SetTitle(Form("p_{T} of reconstructed %s, %s ",particleN.Data(), comment[0].Data()));
995  gPtBCPu ->GetHistogram()->SetXTitle("p_{T} (GeV/c) ");
996  gPtBCPu ->GetHistogram()->SetYTitle(Form("dN_{%s}/dp_{T} (GeV/c)^{-1} / N_{events} ",particleN.Data()));
997  //gPtBCPu ->GetHistogram()->SetAxisRange(0.,30);
998  gPtBCPu ->GetHistogram()->SetTitleOffset(1.6,"Y");
999  gPtBCPu ->SetMarkerStyle(24);
1000  gPtBCPu ->SetMarkerSize(1);
1001  gPtBCPu ->SetMarkerColor(2);
1002  }
1003 
1004 
1005  //PLOTS
1006  //
1007 
1008  //SM
1009  TLegend pLegendFitM(0.15,0.53,0.3,0.93);
1010  pLegendFitM.SetTextSize(0.03);
1011  pLegendFitM.SetFillColor(10);
1012  pLegendFitM.SetBorderSize(1);
1013 
1014  TLegend pLegendFitW(0.4,0.53,0.6,0.93);
1015  pLegendFitW.SetTextSize(0.03);
1016  pLegendFitW.SetFillColor(10);
1017  pLegendFitW.SetBorderSize(1);
1018 
1019  TCanvas *cFitSM = new TCanvas("cFitSM","cFitSM",1200,600);
1020  cFitSM->Divide(2, 1);
1021 
1022  cFitSM->cd(1);
1023  //......................................................................
1024  gPad->SetGridx();
1025  gPad->SetGridy();
1026 
1027  if(part=="Pi0")
1028  {
1029  gMass[0]->SetMaximum(180);
1030  gMass[0]->SetMinimum(120);
1031  }
1032  else
1033  {
1034  gMass[0]->SetMaximum(660);
1035  gMass[0]->SetMinimum(420);
1036  }
1037 
1038  gMass[0]->SetMarkerStyle(modStyleIndex[0]);
1039  gMass[0]->SetMarkerColor(modColorIndex[0]);
1040  gMass[0]->SetLineColor(modColorIndex[0]);
1041 
1042  gMass[0]->Draw("AP");
1043 
1044  pLegendFitM.AddEntry(gMass[0],Form("%s",leg[0].Data()),"P");
1045 
1046  if(drawAllCombi)
1047  {
1048  for(Int_t icomb = 1; icomb <11; icomb++)
1049  {
1050  //printf("option %d\n",icomb);
1051  pLegendFitM.AddEntry(gMass[icomb],Form("%s",leg[icomb].Data()),"P");
1052  gMass[icomb]->SetMarkerStyle(modStyleIndex[icomb]);
1053  gMass[icomb]->SetMarkerColor(modColorIndex[icomb]);
1054  gMass[icomb]->SetLineColor(modColorIndex[icomb]);
1055  gMass[icomb]->Draw("P");
1056  }
1057  }
1058 
1059  gMass[0]->Draw("P");
1060  if(ncomb > 1) pLegendFitM.Draw();
1061 
1062  cFitSM->cd(2);
1063  gPad->SetGridx();
1064  gPad->SetGridy();
1065 
1066  if(part=="Pi0")
1067  {
1068  gWidth[0]->SetMaximum(16);
1069  gWidth[0]->SetMinimum(6);
1070  }
1071  else
1072  {
1073  gWidth[0]->SetMaximum(60);
1074  gWidth[0]->SetMinimum(0);
1075  }
1076 
1077  gWidth[0]->SetMarkerStyle(modStyleIndex[0]);
1078  gWidth[0]->SetMarkerColor(modColorIndex[0]);
1079  gWidth[0]->SetLineColor (modColorIndex[0]);
1080 
1081  gWidth[0]->Draw("AP");
1082  pLegendFitW.AddEntry(gWidth[0],Form("%s",leg[0].Data()),"P");
1083 
1084  if(drawAllCombi)
1085  {
1086  for(Int_t icomb = 1; icomb <11; icomb++)
1087  {
1088  //printf("option %d\n",icomb);
1089  pLegendFitW.AddEntry(gWidth[icomb],Form("%s",leg[icomb].Data()),"P");
1090  gWidth[icomb]->SetMarkerStyle(modStyleIndex[icomb]);
1091  gWidth[icomb]->SetMarkerColor(modColorIndex[icomb]);
1092  gWidth[icomb]->SetLineColor(modColorIndex[icomb]);
1093  gWidth[icomb]->Draw("P");
1094 
1095  }
1096  }
1097 
1098  gWidth[0]->Draw("P");
1099  if(ncomb > 1) pLegendFitW.Draw();
1100 
1101  cFitSM->Print(Form("IMfigures/%s_%s_SingleSM_MassWidth_%s_%s.%s",prodname.Data(),calorimeter.Data(),part.Data(),filename.Data(),fileFormat.Data()));
1102 
1103 
1104 // if(drawAllCombi)
1105 // {
1106 // //Pi0
1107 // //Sector
1108 // TLegend pLegendFitMS(0.75,0.15,0.95,0.45);
1109 // pLegendFitMS.SetTextSize(0.03);
1110 // pLegendFitMS.SetFillColor(10);
1111 // pLegendFitMS.SetBorderSize(1);
1112 //
1113 //
1114 // TLegend pLegendFitWS(0.7,0.6,0.9,0.9);
1115 // pLegendFitWS.SetTextSize(0.03);
1116 // pLegendFitWS.SetFillColor(10);
1117 // pLegendFitWS.SetBorderSize(1);
1118 //
1119 // TCanvas *cFitS = new TCanvas("cFitS","cFitS",1200,600);
1120 // cFitS->Divide(2, 1);
1121 //
1122 // cFitS->cd(1);
1123 // // //......................................................................
1124 // // TCanvas *c5 = new TCanvas("c5","c5",0,0,600,600);
1125 // gPad->SetGridx();
1126 // gPad->SetGridy();
1127 
1128 // gMass[0]->SetMaximum(0.15);
1129 //
1130 // //printf("Fitting Mass plot\n");
1131 // gMass[0]->Draw("AP");
1132 // pLegendFitMS.AddEntry(gMass[0],Form("%s",leg[0].Data()),"P");
1133 //
1134 // for(Int_t icomb = 11; icomb <16; icomb++){
1135 // //printf("option %d\n",icomb);
1136 // pLegendFitMS.AddEntry(gMass[icomb],Form("%s",leg[icomb].Data()),"P");
1137 // gMass[icomb]->SetMarkerStyle(modStyleIndex[icomb]);
1138 // gMass[icomb]->SetMarkerColor(modColorIndex[icomb]);
1139 // gMass[icomb]->SetLineColor(modColorIndex[icomb]);
1140 // gMass[icomb]->Draw("P");
1141 // }
1142 // pLegendFitMS.Draw();
1143 //
1144 // cFitS->cd(2);
1145 // gPad->SetGridx();
1146 // gPad->SetGridy();
1147 // gWidth[0]->Draw("AP");
1148 // //printf("Width plots\n");
1149 // pLegendFitWS.AddEntry(gWidth[0],Form("%s",leg[0].Data()),"P");
1150 //
1151 // for(Int_t icomb = 11; icomb <16; icomb++){
1152 // //printf("option %d\n",icomb);
1153 // pLegendFitWS.AddEntry(gWidth[icomb],Form("%s",leg[icomb].Data()),"P");
1154 //
1155 // gWidth[icomb]->SetMarkerStyle(modStyleIndex[icomb]);
1156 // gWidth[icomb]->SetMarkerColor(modColorIndex[icomb]);
1157 // gWidth[icomb]->SetLineColor(modColorIndex[icomb]);
1158 // gWidth[icomb]->Draw("P");
1159 //
1160 // }
1161 // pLegendFitWS.Draw();
1162 //
1163 // cFitS->Print(Form("IMfigures/%s_%s_Sector_MassWidth_%s_%s.%s",prodname.Data(),calorimeter.Data(),part.Data(),filename.Data(),fileFormat.Data()));
1164 //
1165 // //Side
1166 // TLegend pLegendFitMSi(0.7,0.15,0.9,0.45);
1167 // pLegendFitMSi.SetTextSize(0.03);
1168 // pLegendFitMSi.SetFillColor(10);
1169 // pLegendFitMSi.SetBorderSize(1);
1170 //
1171 //
1172 // TLegend pLegendFitWSi(0.7,0.5,0.9,0.9);
1173 // pLegendFitWSi.SetTextSize(0.03);
1174 // pLegendFitWSi.SetFillColor(10);
1175 // pLegendFitWSi.SetBorderSize(1);
1176 //
1177 // TCanvas *cFitSi = new TCanvas("cFitSi","cFitSi",1200,600);
1178 // cFitSi->Divide(2, 1);
1179 // //printf("Side \n");
1180 // cFitSi->cd(1);
1181 // // //......................................................................
1182 // // TCanvas *c5 = new TCanvas("c5","c5",0,0,600,600);
1183 // gPad->SetGridx();
1184 // gPad->SetGridy();
1185 // // c5->SetLeftMargin(0.12);
1186 // // c5->SetRightMargin(0.08);
1187 // //printf("Fitting Mass plot\n");
1188 // gMass[0]->Draw("AP");
1189 // pLegendFitMSi.AddEntry(gMass[0],Form("%s",leg[0].Data()),"P");
1190 //
1191 // for(Int_t icomb = 16; icomb <24; icomb++){
1192 // //printf("option %d\n",icomb);
1193 // pLegendFitMSi.AddEntry(gMass[icomb],Form("%s",leg[icomb].Data()),"P");
1194 // gMass[icomb]->SetMarkerStyle(modStyleIndex[icomb]);
1195 // gMass[icomb]->SetMarkerColor(modColorIndex[icomb]);
1196 // gMass[icomb]->SetLineColor(modColorIndex[icomb]);
1197 // gMass[icomb]->Draw("P");
1198 // }
1199 // pLegendFitMSi.Draw();
1200 //
1201 // cFitSi->cd(2);
1202 // gPad->SetGridx();
1203 // gPad->SetGridy();
1204 // // c6->SetLeftMargin(0.12);
1205 // // c6->SetRightMargin(0.08);
1206 // gWidth[0]->Draw("AP");
1207 // //printf("Width plots\n");
1208 // pLegendFitWSi.AddEntry(gWidth[0],Form("%s",leg[0].Data()),"P");
1209 //
1210 // for(Int_t icomb = 16; icomb <24; icomb++){
1211 // //printf("option w %d\n",icomb);
1212 // pLegendFitWSi.AddEntry(gWidth[icomb],Form("%s",leg[icomb].Data()),"P");
1213 //
1214 // gWidth[icomb]->SetMarkerStyle(modStyleIndex[icomb]);
1215 // gWidth[icomb]->SetMarkerColor(modColorIndex[icomb]);
1216 // gWidth[icomb]->SetLineColor(modColorIndex[icomb]);
1217 // gWidth[icomb]->Draw("P");
1218 //
1219 // }
1220 // pLegendFitWSi.Draw();
1221 //
1222 // cFitSi->Print(Form("IMfigures/%s_%s_Side_MassWidth_%s_%s.%s",prodname.Data(),calorimeter.Data(),part.Data(),filename.Data(),fileFormat.Data()));
1223 //
1224 // }
1225 
1226  // ---- Pt ----
1227 
1228  TLegend pLegendFitPt(0.6,0.6,0.9,0.9);
1229  pLegendFitPt.SetTextSize(0.03);
1230  pLegendFitPt.SetFillColor(10);
1231  pLegendFitPt.SetBorderSize(1);
1232 
1233  TCanvas *cFitPt = new TCanvas("cFitPt","cFitPt",600,600);
1234  cFitPt->Divide(1, 1);
1235 
1236  cFitPt->cd(1);
1237 
1238  // gPad->SetGridx();
1239  // gPad->SetGridy();
1240  gPad->SetLogy();
1241 
1242  gPt[0]->SetTitle(Form(" Number of %s vs p_{T}, %s ",particle.Data(),comment[0].Data()));
1243  gPt[0]->SetMarkerStyle(modStyleIndex[0]);
1244  gPt[0]->SetMarkerColor(modColorIndex[0]);
1245  gPt[0]->SetLineColor(modColorIndex[0]);
1246 
1247  gPt[0]->Draw("AP");
1248 
1249  // for(Int_t icomb = 1; icomb <ncomb; icomb++){
1250  // gPt[icomb]->SetMarkerStyle(modStyleIndex[icomb]);
1251  // gPt[icomb]->SetMarkerColor(modColorIndex[icomb]);
1252  // gPt[icomb]->SetLineColor(modColorIndex[icomb]);
1253  // gPt[icomb]->Draw("AP");
1254  // //pLegendFitPt.AddEntry(gPt[icomb],Form("%s",leg[icomb].Data()),"P");
1255  // }
1256 
1257  //pLegendFitPt.Draw();
1258  cFitPt->Print(Form("IMfigures/%s_%s_Pt_%s_%s.%s",prodname.Data(),calorimeter.Data(),part.Data(),filename.Data(),fileFormat.Data()));
1259 
1260 
1261  //
1262 
1263  TGraphErrors *gEff = 0;
1264  TGraphErrors *gPrim = 0;
1265 // if(filename.Contains("simu"))
1266 // {
1267 // TH2D* h2Den = (TH2D*) file ->Get("AnaPi0_TM1_PrimPi0Rapidity");
1268 // TH1D* hY = h2Den ->ProjectionY("PtYproj",-1,-1);
1269 // Int_t binP07 = hY->FindBin(+0.65);
1270 // Int_t binM07 = hY->FindBin(-0.65);
1271 // TH1D* hDen = (TH1D*)h2Den->ProjectionX("PrimaryPt" ,binM07 ,binP07 );
1272 //
1273 // for(Int_t ibin = 0; ibin < nPt; ibin++ )
1274 // {
1275 // Double_t ptMin = xPtLimits[ibin];
1276 // Double_t ptMax = xPtLimits[ibin+1];
1277 // primMeson [ibin] = hDen->Integral(hDen->FindBin(ptMin), hDen->FindBin(ptMax));
1278 // eprimMeson[ibin] = TMath::Sqrt(primMeson[ibin]);
1279 // //printf("Bin %d, [%1.1f,%1.1f], content num %f, content den %f times evt %e, bin size %f\n",ibin,xPt[ibin],xPt[ibin+1],mesonPtBC[0][ibin],denBin,nEvt,exPt[ibin]*2);
1280 //
1281 // primMeson [ibin]/=(nEvt*(exPt[ibin]*2));
1282 // eprimMeson[ibin]/=(nEvt*(exPt[ibin]*2));
1283 // if(filename.Contains("pp_7TeV_Pi0"))
1284 // {
1285 // primMeson[ibin]*=1.e6;
1286 // primMeson[ibin]*=360./100.;
1287 //
1288 // //eprimMeson[ibin]*=1.e6*1.e6;
1289 // eprimMeson[ibin]*=360./100.;
1290 //
1291 // mesonPtBC [0][ibin] /= (0.62 +xPt[ibin]*0.0017 );
1292 // emesonPtBC[0][ibin] /= (0.62 +xPt[ibin]*0.0017 );
1293 //
1294 // primMeson [ibin] /= (0.56 +xPt[ibin]*0.0096 );
1295 // eprimMeson[ibin] /= (0.56 +xPt[ibin]*0.0096 );
1296 // }
1297 //
1298 // if(primMeson[ibin] <=0) continue;
1299 //
1300 // mesonEffPt [0][ibin] = mesonPtBC [0][ibin] / primMeson[ibin] ;
1301 //
1302 // if(mesonEffPt[0][ibin] < 0)
1303 // {
1304 // mesonEffPt [0][ibin] = 0;
1305 // emesonEffPt[0][ibin] = 0;
1306 // }
1307 // else emesonEffPt[0][ibin] = mesonEffPt[0][ibin] * TMath::Sqrt(1./(primMeson[ibin]*nEvt)+1./(mesonPtBC[0][ibin]*nEvt)*2);
1308 //
1309 // printf("Bin %d, [%1.1f,%1.1f], Num %2.3e, Den %2.3e, eDen %e, Eff %f, err %f\n",ibin,xPt[ibin],xPt[ibin+1],
1310 // mesonPtBC[0][ibin],primMeson[ibin],eprimMeson[ibin], mesonEffPt[0][ibin],emesonEffPt[0][ibin]);
1311 // }
1312 //
1313 // gEff = new TGraphErrors(nPt,xPt,mesonEffPt[0],exPt,emesonEffPt[0]);
1314 // gEff->SetName("EffxAcc");
1315 // gEff->GetHistogram()->SetYTitle("#epsilon_{Reco #times PID #times Acc}");
1316 //
1317 //
1318 // gPrim = new TGraphErrors(nPt,xPt,primMeson,exPt,eprimMeson);
1319 // gPrim->SetName("Primary");
1320 // gPrim->GetHistogram()->SetYTitle("Primary");
1321 // }
1322 
1323 
1324  //
1325  if(mix)
1326  {
1327  TCanvas *cFitPtBC = new TCanvas("cFitPtBC","cFitPtBC",600,600);
1328  cFitPtBC->Divide(1, 1);
1329 
1330  cFitPtBC->cd(1);
1331 
1332  // gPad->SetGridx();
1333  // gPad->SetGridy();
1334  gPad->SetLogy();
1335  // // c3->SetLeftMargin(0.12);
1336  // // c3->SetRightMargin(0.08);
1337  gPtBC[0]->SetTitle(Form(" Number of %s vs p_{T}, bin counting, %s ",particle.Data(),comment[0].Data()));
1338  gPtBC[0]->SetMarkerStyle(modStyleIndex[0]);
1339  gPtBC[0]->SetMarkerColor(modColorIndex[0]);
1340  gPtBC[0]->SetLineColor(modColorIndex[0]);
1341  // if(part=="Pi0"){
1342  // gPtBC[0]->SetMaximum(1e-1);
1343  // gPtBC[0]->SetMinimum(1e-6);
1344  // }
1345  // else{
1346  // gPtBC[0]->SetMaximum(1e-1);
1347  // gPtBC[0]->SetMinimum(1e-6);
1348  // }
1349 
1350  gPtBC[0]->Draw("AP");
1351 
1352  gPt[0]->SetMarkerColor(2);
1353  gPt[0]->SetMarkerStyle(22);
1354 
1355  gPt[0]->Draw("P");
1356  if(filename.Contains("imu"))
1357  {
1358  gPtBCPu->SetMaximum(1e-1);
1359  gPtBCPu->SetMinimum(1e-6);
1360  gPtBCPu->Draw("P");
1361 
1362  // gPrim->SetMarkerStyle(25);
1363  // gPrim->SetMarkerColor(4);
1364  // gPrim->Draw("P");
1365  }
1366 
1367  // for(Int_t icomb = 1; icomb <ncomb; icomb++){
1368  // gPtBC[icomb]->SetMarkerStyle(modStyleIndex[icomb]);
1369  // gPtBC[icomb]->SetMarkerColor(modColorIndex[icomb]);
1370  // gPtBC[icomb]->SetLineColor(modColorIndex[icomb]);
1371  // gPtBC[icomb]->Draw("AP");
1372  // //pLegendFitPtBC.AddEntry(gPtBC[icomb],Form("%s",leg[icomb].Data()),"P");
1373  // }
1374 
1375  //pLegendFitPtBC.Draw();
1376 
1377  TLegend legBC(0.5,0.5,0.9,0.9);
1378  legBC.AddEntry(gPt[0],"Fit","P");
1379  legBC.AddEntry(gPtBC[0],"Integrated","P");
1380  legBC.AddEntry(gPtBCPu ,"Integrated pure","P");
1381  // legBC.AddEntry(gPrim ,"Primary","P");
1382  legBC.Draw();
1383 
1384  cFitPtBC->Print(Form("IMfigures/%s_%s_PtBC_%s_%s.%s",prodname.Data(),calorimeter.Data(),part.Data(),filename.Data(),fileFormat.Data()));
1385  }
1386 
1387  //--------------------
1388  // Store plots in file
1389  //--------------------
1390 
1391  TFile * fout = new TFile(Form("IMfigures/%s_%s_MassWidthPtHistograms_%s_%s_%s.root",prodname.Data(),calorimeter.Data(),part.Data(),filename.Data(),hname[0].Data()), "recreate");
1392 
1393  for(Int_t icomb = 0; icomb <ncomb; icomb++)
1394  {
1395  for(Int_t ipt = 0; ipt < nPt; ipt++)
1396  {
1397  //if(hIM[icomb][ipt]) { hIM[icomb][ipt]->Scale(1./nEvt); hIM[icomb][ipt]->Write(); }
1398  if(hIM[icomb][ipt]) { hIM[icomb][ipt]->Write(); }
1399  if(mix)
1400  {
1401 // if(hMix [icomb][ipt]) { hMix [icomb][ipt]->Scale(1./nEvt); hMix [icomb][ipt]->Write();}
1402 // if(hMixCorrected[icomb][ipt]) { hMixCorrected[icomb][ipt]->Scale(1./nEvt); hMixCorrected[icomb][ipt]->Write();}
1403 // if(hRatio [icomb][ipt]) { hRatio [icomb][ipt]->Scale(1./nEvt); hRatio [icomb][ipt]->Write();}
1404 // if(hSignal [icomb][ipt]) { hSignal [icomb][ipt]->Scale(1./nEvt); hSignal [icomb][ipt]->Write();}
1405  if(hMix [icomb][ipt]) { hMix [icomb][ipt]->Write();}
1406  if(hMixCorrected[icomb][ipt]) { hMixCorrected[icomb][ipt]->Write();}
1407  if(hRatio [icomb][ipt]) { hRatio [icomb][ipt]->Write();}
1408  if(hSignal [icomb][ipt]) { hSignal [icomb][ipt]->Write();}
1409  }
1410  }
1411 
1412  gPt [icomb]->Write();
1413  gMass [icomb]->Write();
1414  gWidth[icomb]->Write();
1415  if(mix) gPtBC [icomb]->Write();
1416  }
1417 
1418  if(gEff) gEff->Write();
1419  if(gPtBCPu) gPtBCPu->Write();
1420  // if(gPrim) gPrim->Write();
1421  hRe[0]->Write();
1422  if(hMi[0])hMi[0]->Write();
1423 
1424  fout->Close();
1425 
1426 }
1427 
1430 //-----------------------------------------------------------------------------
1431 Bool_t GetFileAndEvents( TString prodname, TString filename,
1432  TString dirName , TString listName )
1433 {
1434  fil = new TFile(Form("%s/%s.root",prodname.Data(),filename.Data()),"read");
1435 
1436  if ( !fil ) return kFALSE;
1437 
1438  direc = (TDirectoryFile*) fil->Get(dirName);
1439 
1440  if ( !direc && dirName != "" ) return kFALSE;
1441 
1442  //printf("dir %p, list %s\n",dir,listName.Data());
1443 
1444  if(direc)
1445  lis = (TList*) direc ->Get(listName);
1446  else
1447  lis = (TList*) fil->Get(listName);
1448 
1449  if ( !lis && listName != "") return kFALSE;
1450 
1451  if(!lis)
1452  nEvt = ((TH1F*) fil->Get("hNEvents"))->GetEntries();
1453  else
1454  nEvt = ((TH1F*) lis->FindObject("hNEvents"))->GetEntries();
1455 
1456  printf("nEvt = %d\n",nEvt);
1457  //nEvt = 1;
1458 
1459  return kTRUE;
1460 }
1461 
1465 //-----------------------------------------------------------------------------
1466 void SetFitFun()
1467 {
1468  //Fitting function
1469  //if(mix) polN = 0;
1470 
1471  if(part == "Pi0")
1472  {
1473  if ( polN == 0)
1474  fitfun = new TF1("fitfun",pi0massP0,0.100,0.250,4);
1475  else if (polN == 1)
1476  fitfun = new TF1("fitfun",pi0massP1,0.100,0.250,5);
1477  else if (polN == 2)
1478  fitfun = new TF1("fitfun",pi0massP2,0.100,0.250,6);
1479  else if (polN == 3)
1480  fitfun = new TF1("fitfun",pi0massP3,0.100,0.250,7);
1481  else
1482  {
1483  printf("*** <<< Set Crystal Ball!!! >>> ***\n");
1484  fitfun = new TF1("fitfun",CrystalBall,0.100,0.250,6);
1485  }
1486 
1487  if(polN < 4)
1488  {
1489  fitfun->SetParLimits(0, nPairCut/5,nPairCut*1.e4);
1490  fitfun->SetParLimits(1, 0.105,0.185);
1491  fitfun->SetParLimits(2, 0.001,0.040);
1492  }
1493  else
1494  {
1495  fitfun->SetParLimits(0, nPairCut/5,nPairCut*1.e6);
1496  fitfun->SetParLimits(2, 0.105,0.185);
1497  fitfun->SetParLimits(1, 0.001,0.040);
1498  fitfun->SetParLimits(3, 0.001,0.040);
1499  fitfun->SetParLimits(4, 0,10);
1500  fitfun->SetParLimits(5, 1,1.e+6);
1501  }
1502 
1503  } // Pi0
1504  else // Eta
1505  {
1506  if ( polN == 0)
1507  fitfun = new TF1("fitfun",pi0massP0,0.400,0.650,4);
1508  else if (polN == 1)
1509  fitfun = new TF1("fitfun",pi0massP1,0.400,0.650,5);
1510  else if (polN == 2)
1511  fitfun = new TF1("fitfun",pi0massP2,0.400,0.650,6);
1512  else if (polN == 3)
1513  fitfun = new TF1("fitfun",pi0massP3,0.400,0.650,7);
1514  if(polN < 4){
1515  fitfun->SetParLimits(0, nPairCut/10,1.e+6);
1516  fitfun->SetParLimits(1, 0.20,0.80);
1517  fitfun->SetParLimits(2, 0.001,0.06);
1518  }
1519  }
1520 
1521  fitfun->SetLineColor(kRed);
1522  fitfun->SetLineWidth(2);
1523 
1524  fitfun->SetParName(0,"A");
1525  fitfun->SetParName(1,"m_{0}");
1526  fitfun->SetParName(2,"#sigma");
1527 // fitfun->SetParName(3,"a_{0}");
1528 // fitfun->SetParName(4,"a_{1}");
1529 
1530  if (polN > 1)
1531  {
1532  fitfun->SetParName(5,"a_{2}");
1533  if (polN > 2) fitfun->SetParName(6,"a_{3}");
1534  }
1535 
1536 }
1537 
1538 //-----------------------------------------------------------------------------
1539 Double_t pi0massP3(Double_t *x, Double_t *par)
1540 {
1541  Double_t gaus = par[0] * TMath::Exp( -(x[0]-par[1])*(x[0]-par[1]) /
1542  (2*par[2]*par[2]) );
1543  Double_t back = par[3] + par[4]*x[0] + par[5]*x[0]*x[0] + par[6]*x[0]*x[0]*x[0];
1544  return gaus+back;
1545 }
1546 
1547 //-----------------------------------------------------------------------------
1548 Double_t pi0massP2(Double_t *x, Double_t *par)
1549 {
1550  Double_t gaus = par[0] * TMath::Exp( -(x[0]-par[1])*(x[0]-par[1]) /
1551  (2*par[2]*par[2]) );
1552  Double_t back = par[3] + par[4]*x[0] + par[5]*x[0]*x[0];
1553  return gaus+back;
1554 }
1555 
1556 //-----------------------------------------------------------------------------
1557 Double_t pi0massP1(Double_t *x, Double_t *par)
1558 {
1559  Double_t gaus = par[0] * TMath::Exp( -(x[0]-par[1])*(x[0]-par[1]) /
1560  (2*par[2]*par[2]) );
1561  Double_t back = par[3] + par[4]*x[0];
1562  return gaus+back;
1563 }
1564 
1565 //-----------------------------------------------------------------------------
1566 Double_t pi0massP0(Double_t *x, Double_t *par)
1567 {
1568  Double_t gaus = par[0] * TMath::Exp( -(x[0]-par[1])*(x[0]-par[1]) /
1569  (2*par[2]*par[2]) );
1570  Double_t back = 0;//par[3];
1571  return gaus+back;
1572 }
1573 
1574 //-----------------------------------------------------------------------------
1575 Double_t truncatedPolEta(Double_t *x, Double_t *par)
1576 {
1577  if ((x[0] > 0.425 && x[0] < 0.65) ) {
1578  TF1::RejectPoint();
1579  return 0;
1580  }
1581 
1582  return par[0] + par[1]*x[0];// + x[0]*x[0]*par[2] + x[0]*x[0]*x[0]*par[3];//+ x[0]*x[0]*x[0]*x[0]*par[4];
1583 }
1584 
1585 //-----------------------------------------------------------------------------
1586 Double_t truncatedPolPi0(Double_t *x, Double_t *par)
1587 {
1588  if ((x[0] > 0.07 && x[0] < 0.2)) {
1589  TF1::RejectPoint();
1590  return 0;
1591  }
1592 
1593  return par[0] + par[1]*x[0];// + x[0]*x[0]*par[2] + x[0]*x[0]*x[0]*par[3] ;//+ x[0]*x[0]*x[0]*x[0]*par[4];
1594 }
1595 
1596 //-----------------------------------------------------------------------------
1597 Double_t CrystalBall(Double_t *x, Double_t *par)
1598 {
1599  Double_t N = par[0];
1600  Double_t width = par[1];
1601 
1602  Double_t mean = par[2];
1603  Double_t sigma = par[3];
1604  Double_t alpha = par[4];
1605  Double_t n = par[5];
1606 
1607  Double_t A = pow(n/fabs(alpha),n)*exp(-pow(alpha,2)/2);
1608  Double_t B = n/fabs(alpha) - fabs(alpha);
1609 
1610  if ((x[0]-mean)/sigma>-alpha)
1611  return N*width*TMath::Gaus(x[0],mean,sigma,1);
1612  else
1613  return N/(sqrt(TMath::TwoPi())*sigma)*width*A*pow(B-(x[0]-mean)/sigma,-n);
1614 }
hSignal
TH1D * hSignal
Definition: CalculateAveragePt.C:87
filename
const char * filename
Definition: TestFCM.C:1
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
nPairCut
Float_t nPairCut
Plot also many SM combinations.
Definition: InvMassFit.C:56
drawAllCombi
Bool_t drawAllCombi
Apply Root method Sumw2()
Definition: InvMassFit.C:55
pi0massP3
Double_t pi0massP3(Double_t *x, Double_t *par)
Definition: InvMassFit.C:1539
mix
Bool_t mix
Definition: InvMassFit.C:51
mass
Double_t mass
Definition: charmCutsOptimization.C:35
InvMassFit
void InvMassFit(TString prodname="LHC18c3_NystromOn", TString filename="AnalysisResults", TString histoDir="Pi0IM_GammaTrackCorr_EMCAL", TString histoList="default", TString calorimeter="DCAL", TString particle="Pi0", Bool_t mixed=0, Int_t pol=1, Int_t ncomb=1, Float_t nPairMin=20, TString fileFormat="eps")
Definition: InvMassFit.C:82
nEvt
Int_t nEvt
Minimum number of cluster pairs in pi0 or eta window.
Definition: InvMassFit.C:59
sumw2
Bool_t sumw2
polinomyal type for residual background under the peak
Definition: InvMassFit.C:54
GetFileAndEvents
Bool_t GetFileAndEvents(TString prodname, TString filename, TString dirName, TString listName)
Definition: InvMassFit.C:1431
ptMin
Double_t ptMin
Definition: MuonTrackingEfficiency.C:52
mesonMass
Double_t mesonMass
Definition: CalculateAveragePt.C:60
calorimeter
const TString calorimeter
Definition: anaM.C:36
part
TString part
use mixed event to constrain combinatorial background
Definition: InvMassFit.C:52
fil
TFile * fil
Definition: InvMassFit.C:60
pi0massP1
Double_t pi0massP1(Double_t *x, Double_t *par)
Definition: InvMassFit.C:1557
sigma
Double_t * sigma
Definition: CalculateAveragePt.C:58
lis
TList * lis
Definition: InvMassFit.C:61
Int_t
int Int_t
Definition: External.C:63
SetFitFun
void SetFitFun()
Definition: InvMassFit.C:1466
Float_t
float Float_t
Definition: External.C:68
CrystalBall
Double_t CrystalBall(Double_t *x, Double_t *par)
Definition: InvMassFit.C:1597
direc
TDirectoryFile * direc
Definition: InvMassFit.C:62
TH1D
Definition: External.C:212
truncatedPolPi0
Double_t truncatedPolPi0(Double_t *x, Double_t *par)
Definition: InvMassFit.C:1586
Data
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
Definition: charmCutsOptimization.C:371
pi0massP0
Double_t pi0massP0(Double_t *x, Double_t *par)
Definition: InvMassFit.C:1566
rebin
Int_t rebin
Definition: charmCutsOptimization.C:32
fitfun
TF1 * fitfun
Definition: InvMassFit.C:63
pi0massP2
Double_t pi0massP2(Double_t *x, Double_t *par)
Definition: InvMassFit.C:1548
polN
Int_t polN
define fitting and plotting ranges for particle
Definition: InvMassFit.C:53
Bool_t
bool Bool_t
Definition: External.C:53
fout
TFile * fout
input train file
Definition: DrawAnaCaloTrackQA.C:90
ptMax
Double_t ptMax
Definition: MuonTrackingEfficiency.C:53
truncatedPolEta
Double_t truncatedPolEta(Double_t *x, Double_t *par)
Definition: InvMassFit.C:1575