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