AliPhysics  ed43440 (ed43440)
ProjectCombinHFAndFit.C
Go to the documentation of this file.
1 #if !defined(__CINT__) || defined(__MAKECINT__)
2 #include <TString.h>
3 #include <TMath.h>
4 #include <TFile.h>
5 #include <TCanvas.h>
6 #include <TH1F.h>
7 #include <TH2F.h>
8 #include <TH3F.h>
9 #include <TH1D.h>
10 #include <TF1.h>
11 #include <TStyle.h>
12 #include <TLatex.h>
13 #include <TPaveText.h>
14 #include <TLegend.h>
15 #include <TLegendEntry.h>
16 #include <TDatabasePDG.h>
17 #include "AliHFInvMassFitter.h"
18 #include "AliVertexingHFUtils.h"
19 #include "AliNormalizationCounter.h"
20 #endif
21 
22 enum Method{kME,kRot,kLS,kSB};
23 
24 // input files and pt binning
25 TString fileName="DataTrains/AnalysisResults_17pq_FAST_wSDD_train2210.root";
26 TString suffix="3SigPID_Pt300_FidY_PilMV_EM1";
27 TString fileNameMC="MCTrains/AnalysisResults_LHC17pq_FAST_CENTwSDD_G3_train1169-1168.root";
28 //TString fileNameMC="MCTrains/AnalysisResults_LHC17pq_FAST_CENTwSDD_G4_train892-891.root";
29 TString suffixMC="_3SigPID_Pt300_FidY_PilMV_EM1";
30 
31 // TString fileName="DataTrains/AnalysisResults_17pq_FAST_wSDD_train2104.root";
32 // TString suffix="2SigPID_Pt300_FidY_PilSPD5_EM1";
33 // //TString fileNameMC="MCTrains/AnalysisResults_LHC17pq_FAST_CENTwSDD_G4_train870-869.root";
34 // TString fileNameMC="MCTrains/AnalysisResults_LHC17pq_FAST_CENTwSDD_G3_train868-867.root";
35 // TString suffixMC="_Prompt_2SigPID_Pt300_FidY_PilSPD5_EM1";
36 
37 TString meson="Dzero";
38 const Int_t nPtBins=8;
39 Double_t binLims[nPtBins+1]={0.,1.,2.,3.,4.,5.,6.,8.,12.};
40 //Double_t binLims[nPtBins+1]={0.,0.5,1.,1.5,2.,2.5,3.,3.5,4.};
41 Double_t sigmas[nPtBins]={0.006,0.008,0.009,0.010,0.011,0.012,0.013,0.013};
42 
43 // outputfiles
44 Bool_t saveCanvasAsRoot=kTRUE;
45 Int_t saveCanvasAsEps=1; //0=none, 1=main ones, 2=all
46 
47 // fit configuration
48 Int_t rebin[nPtBins]={4,6,7,8,9,10,10,12};//8,9,10,12,12,12,12};
49 Int_t fixSigmaConf=0; // 0= all free, 1=all fixed, 2=use values per pt bin
50 Bool_t fixSigma[nPtBins]={kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE};
51 Double_t tuneSigmaOnData=-1.00; // scaling factor for the Gaussian sigma, if -1. use MC out of the box
52 Int_t fixMeanConf=0; // 0= all free, 1=all fixed, 2=use values per pt bin
53 Bool_t fixMean[nPtBins]={kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE};
54 Double_t minMass=1.72;
55 Double_t maxMass=2.04;
56 Int_t optForNorm=1;
57 Double_t rangeForNorm=0.05;
58 TString fitoption="E";
59 Bool_t useEMwithLS=kTRUE;
60 Int_t typeb=2;
61 Double_t nsigmaBinCounting=4.; // defines the mass interval over which the signal is bin counted
62 Int_t optBkgBinCount=1;
63 Double_t massD=-1.; // use negative value to use PDG mass
64 
65 Int_t smoothLS=0;
66 
67 // reflection option
68 TString reflopt="2gaus";
69 Bool_t correctForRefl=kTRUE;
70 Double_t rOverSmodif=1;
71 
72 // objects and options related to side-band fit method
73 Bool_t tryDirectFit=kTRUE;
74 
75 Int_t nparback=0;
76 Double_t fitrangelow[nPtBins]={1.74,1.74,1.74,1.72,1.72,1.72,1.72,1.72};
77 Double_t fitrangeup[nPtBins]={2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.04};
78 Int_t nDegreeBackPol[nPtBins]={4,4,4,2,2,2,2,2}; // degree of polynomial function describing the background
79 
80 TH1D* hMCReflPtBin;
81 TH1D* hMCSigPtBin;
82 
83 void WriteFitInfo(AliHFInvMassFitter *fitter, TH1D* histo);
84 void WriteFitFunctionsToFile(AliHFInvMassFitter *fitter, TString meth, Int_t iPtBin);
85 TH1F* FitMCInvMassSpectra(TList* lMC);
86 
87 AliHFInvMassFitter* ConfigureFitter(TH1D* histo, Int_t iPtBin, Int_t backcase, Double_t minFit, Double_t maxFit){
88  TH1F* histof=(TH1F*)histo->Clone(Form("%s_Fl",histo->GetName()));
89 
90 
91  AliHFInvMassFitter* fitter=new AliHFInvMassFitter(histof,minFit,maxFit,backcase,0);
92  if(backcase==6)fitter->SetPolDegreeForBackgroundFit(nDegreeBackPol[iPtBin]);
93  fitter->SetUseChi2Fit();
94  fitter->SetFitOption(fitoption.Data());
95  fitter->SetInitialGaussianMean(massD);
96  fitter->SetInitialGaussianSigma(sigmas[iPtBin]);
97  if(fixSigmaConf==1 || (fixSigmaConf==2 && fixSigma[iPtBin])) fitter->SetFixGaussianSigma(sigmas[iPtBin]);
98  if(fixMeanConf==1 || (fixMeanConf==2 && fixMean[iPtBin])) fitter->SetFixGaussianMean(massD);
99  if(correctForRefl){
100  TCanvas *cTest=new TCanvas("cTest","cTest",800,800);
101  TH1F *hmasstemp=fitter->GetHistoClone();
102  TH1F *hReflModif=(TH1F*)AliVertexingHFUtils::AdaptTemplateRangeAndBinning(hMCReflPtBin,hmasstemp,minFit,maxFit);
103  TH1F *hSigModif=(TH1F*)AliVertexingHFUtils::AdaptTemplateRangeAndBinning(hMCSigPtBin,hmasstemp,minFit,maxFit);
104  hReflModif->SetLineColor(kRed);
105  hSigModif->SetLineColor(kBlue);
106  hSigModif->Draw();
107  hReflModif->Draw("same");
108  cTest->SaveAs(Form("figures/cTest%d.eps",iPtBin));
109  delete hmasstemp;
110  Double_t fixSoverRefAt=rOverSmodif*(hReflModif->Integral(hReflModif->FindBin(minFit*1.0001),hReflModif->FindBin(maxFit*0.999))/hSigModif->Integral(hSigModif->FindBin(minFit*1.0001),hSigModif->FindBin(maxFit*0.999)));
111  TH1F* hrfl=fitter->SetTemplateReflections(hReflModif,reflopt,minFit,maxFit);
112  if(!hrfl){
113  Printf("SOMETHING WENT WRONG WHILE SETTINGS REFLECTIONS TEMPLATE");
114  delete hReflModif;
115  delete hSigModif;
116  delete cTest;
117  return 0x0;
118  }
119  fitter->SetFixReflOverS(fixSoverRefAt);
120  delete cTest;
121  delete hReflModif;
122  delete hSigModif;
123  }
124  return fitter;
125 }
126 
127 
128 
129 void PrintGausParams(TH1F* hPulls){
130  TF1* fgfit=(TF1*)hPulls->GetListOfFunctions()->FindObject("gaus");
131  TLatex* tg1=new TLatex(0.2,0.8,Form("Gauss mean = %.2f#pm%.2f",fgfit->GetParameter(1),fgfit->GetParError(1)));
132  tg1->SetNDC();
133  tg1->Draw();
134  TLatex* tg2=new TLatex(0.2,0.7,Form("Gauss #sigma = %.2f#pm%.2f",fgfit->GetParameter(2),fgfit->GetParError(2)));
135  tg2->SetNDC();
136  tg2->Draw();
137 
138 }
139 
140 
141 Bool_t QuadraticSmooth(TH1 *h,Int_t ntimes=1){// quadratic fit of 5 points
142  ntimes--;
143  TF1 *fp2=new TF1("fp2","pol2",h->GetXaxis()->GetBinLowEdge(1),h->GetXaxis()->GetBinUpEdge(h->GetNbinsX()));
144  TH1D *htmp=(TH1D*)h->Clone("htmp");
145  for(Int_t i=3;i<=h->GetNbinsX()-3;i++){ // first intermediate bins
146  htmp->Fit(fp2,"RLEMN0","",h->GetXaxis()->GetBinLowEdge(i-2),h->GetXaxis()->GetBinUpEdge(i+2));
147  h->SetBinContent(i,(fp2->Integral(h->GetXaxis()->GetBinLowEdge(i),h->GetXaxis()->GetBinUpEdge(i)))/h->GetBinWidth(i));// change only content, errors unchanged
148  }
149  // now initial and final bins
150  htmp->Fit(fp2,"RLEMN0","",h->GetXaxis()->GetBinLowEdge(1),h->GetXaxis()->GetBinUpEdge(5));// 5 pts, asymmetric fit region
151  h->SetBinContent(2,(fp2->Integral(h->GetXaxis()->GetBinLowEdge(2),h->GetXaxis()->GetBinUpEdge(2)))/h->GetBinWidth(2));// change only content, errors unchanged
152 
153  htmp->Fit(fp2,"RLEMN0","",h->GetXaxis()->GetBinLowEdge(1),h->GetXaxis()->GetBinUpEdge(4));// use only 4 pts
154  h->SetBinContent(1,(fp2->Integral(h->GetXaxis()->GetBinLowEdge(1),h->GetXaxis()->GetBinUpEdge(1)))/h->GetBinWidth(1));// change only content, errors unchanged
155 
156 
157  htmp->Fit(fp2,"RLEMN0","",h->GetXaxis()->GetBinLowEdge(h->GetNbinsX()-4),h->GetXaxis()->GetBinUpEdge(h->GetNbinsX()));// 5 pts, asymmetric fit region
158  h->SetBinContent(h->GetNbinsX()-1,(fp2->Integral(h->GetXaxis()->GetBinLowEdge(h->GetNbinsX()-1),h->GetXaxis()->GetBinUpEdge(h->GetNbinsX()-1)))/h->GetBinWidth(h->GetNbinsX()-1));// change only content, errors unchanged
159 
160  htmp->Fit(fp2,"RLEMN0","",h->GetXaxis()->GetBinLowEdge(h->GetNbinsX()-3),h->GetXaxis()->GetBinUpEdge(h->GetNbinsX()));
161  h->SetBinContent(h->GetNbinsX(),(fp2->Integral(h->GetXaxis()->GetBinLowEdge(h->GetNbinsX()),h->GetXaxis()->GetBinUpEdge(h->GetNbinsX())))/h->GetBinWidth(h->GetNbinsX()));// change only content, errors unchanged
162  if(ntimes>0)QuadraticSmooth(h,ntimes);
163  delete htmp;
164  return kTRUE;
165 }
166 
167 
168 void SetStyleHisto(TH1 *h,Int_t method,Int_t isXpt=-1){
169  if(isXpt==1)h->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");// just for convenience, most of the time it's pt :)
170 
171  if(method==kRot){
172  h->SetMarkerStyle(21);
173  h->GetYaxis()->SetTitleOffset(1.8);
174  h->SetMarkerColor(kBlack);
175  h->SetMarkerSize(1.0);
176  h->SetLineColor(kBlack);
177  return;
178  }
179 
180  if(method==kME){
181  h->SetMarkerStyle(25);
182  h->GetYaxis()->SetTitleOffset(1.8);
183  h->SetMarkerColor(4);
184  h->SetMarkerSize(1.0);
185  h->SetLineColor(4);
186  return;
187  }
188 
189  if(method==kLS){
190  h->SetMarkerStyle(22);
191  h->GetYaxis()->SetTitleOffset(1.8);
192  h->SetMarkerColor(kGreen+2);
193  h->SetMarkerSize(1.0);
194  h->SetLineColor(kGreen+2);
195  return;
196  }
197 
198  if(method==kSB){
199  h->SetMarkerStyle(27);
200  h->GetYaxis()->SetTitleOffset(1.8);
201  h->SetMarkerColor(6);
202  h->SetMarkerSize(1.0);
203  h->SetLineColor(6);
204  return;
205  }
206 
207  return;
208 }
209 
210 void DivideCanvas(TCanvas *c,Int_t ndivisions){
211  if(ndivisions>1){
212  if(ndivisions<3){
213  c->Divide(2,1);
214  }
215  else if(ndivisions<5){
216  c->Divide(2,2);
217  }
218  else if(ndivisions<7){
219  c->Divide(3,2);
220 
221  }
222  else if(ndivisions<=8){
223  c->Divide(4,2);
224 
225  }
226  else if(ndivisions<10){
227  c->Divide(3,3);
228 
229  }
230  else if(ndivisions<13){
231  c->Divide(4,3);
232 
233  }
234  else if(ndivisions<17){
235  c->Divide(4,4);
236 
237  }
238  else {
239  c->Divide(5,5);
240 
241  }
242  }else{
243  c->Divide(1,1);
244  }
245  return;
246 }
247 
248 
249 
250 TF1 *GausPlusLine(Double_t minRange=1.72,Double_t maxRange=2.05){
251  TF1 *f=new TF1("GausPlusLine","[0]*1./(TMath::Sqrt(2.*TMath::Pi())*[2])*TMath::Exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))+[3]+[4]*x",minRange,maxRange);
252 
253  f->SetParameter(1,massD);
254  f->SetParLimits(0,0.,1.e+7);
255  f->SetParLimits(1,1.8,1.92);
256  f->SetParameter(2,0.010);
257  f->SetParLimits(2,0,0.1);
258  f->SetParameter(3,0.);
259  // f->SetParLimits(3,-10000,10000);
260  return f;
261 
262 }
263 
264 
265 
266 Double_t GetBackgroundNormalizationFactor(TH1D* hRatio, Int_t reb=1){
267  //
268  Double_t norm=hRatio->GetMaximum();
269  Double_t massForNorm=0;
270 
271  if(optForNorm==0){
272  Int_t binl=hRatio->GetXaxis()->FindBin(minMass);
273  Int_t binh=hRatio->GetXaxis()->FindBin(maxMass);
274  Double_t norml=hRatio->GetBinContent(1);
275  if(binl>0 && binl<=hRatio->GetNbinsX()) norml=hRatio->GetBinContent(binl);
276  Double_t normh=hRatio->GetBinContent(hRatio->GetNbinsX());
277  if(binh>0 && binh<=hRatio->GetNbinsX()) normh=hRatio->GetBinContent(binh);
278  if(norml>normh){
279  norm=norml;
280  massForNorm=hRatio->GetBinCenter(binl);
281  }else{
282  norm=normh;
283  massForNorm=hRatio->GetBinCenter(binh);
284  }
285  }else if(optForNorm==1){
286  norm=0.0001;
287  for(Int_t iMassBin=1; iMassBin<hRatio->GetNbinsX(); iMassBin++){
288  Double_t bce=hRatio->GetBinCenter(iMassBin);
289  if(bce>minMass && bce<maxMass){
290  Double_t bco=hRatio->GetBinContent(iMassBin);
291  if(bco>norm){
292  norm=bco;
293  massForNorm=bce;
294  }
295  }
296  }
297  }else if(optForNorm==2){
298  norm=0.0001;
299  for(Int_t iMassBin=1; iMassBin<hRatio->GetNbinsX(); iMassBin++){
300  Double_t bce=hRatio->GetBinCenter(iMassBin);
301  Double_t bco=hRatio->GetBinContent(iMassBin);
302  if(bco>norm){
303  norm=bco;
304  massForNorm=bce;
305  }
306  }
307  }else if(optForNorm==3){
308  Int_t binl=hRatio->GetXaxis()->FindBin(minMass);
309  Int_t binh=hRatio->GetXaxis()->FindBin(maxMass);
310  Double_t norml=0;
311  Double_t normReb=0;
312  Int_t iFirstBin=TMath::Max(binl-reb/2,1);
313  Int_t iLastBin=iFirstBin+reb;
314  if(iLastBin>hRatio->GetNbinsX()){
315  iLastBin=hRatio->GetNbinsX();
316  iFirstBin=iLastBin-reb;
317  }
318  for(Int_t jjj=iFirstBin; jjj<=iLastBin; jjj++){
319  norml+=hRatio->GetBinContent(jjj);
320  normReb+=1;
321  }
322  if(normReb>=1) norml/=normReb;
323  Double_t normh=0;
324  normReb=0;
325  iFirstBin=TMath::Max(binh-reb/2,1);
326  iLastBin=iFirstBin+reb;
327  if(iLastBin>hRatio->GetNbinsX()){
328  iLastBin=hRatio->GetNbinsX();
329  iFirstBin=iLastBin-reb;
330  }
331  for(Int_t jjj=iFirstBin; jjj<=iLastBin; jjj++){
332  normh+=hRatio->GetBinContent(jjj);
333  normReb+=1;
334  }
335  if(normReb>=1) normh/=normReb;
336  if(norml>normh){
337  norm=norml;
338  massForNorm=hRatio->GetBinCenter(binl);
339  }else{
340  norm=normh;
341  massForNorm=hRatio->GetBinCenter(binh);
342  }
343  }else if(optForNorm==4){
344  norm=0.0001;
345  for(Int_t iMassBin=1; iMassBin<hRatio->GetNbinsX(); iMassBin++){
346  Double_t bce=hRatio->GetBinCenter(iMassBin);
347  if(bce>minMass && bce<maxMass){
348  Double_t bco=0;
349  Double_t normReb=0;
350  Int_t iFirstBin=TMath::Max(iMassBin-reb/2,1);
351  Int_t iLastBin=iFirstBin+reb;
352  if(iLastBin>hRatio->GetNbinsX()){
353  iLastBin=hRatio->GetNbinsX();
354  iFirstBin=iLastBin-reb;
355  }
356  for(Int_t jjj=iFirstBin; jjj<=iLastBin; jjj++){
357  bco+=hRatio->GetBinContent(jjj);
358  normReb+=1;
359  }
360  if(normReb>=1){
361  bco/=normReb;
362  if(bco>norm){
363  norm=bco;
364  massForNorm=bce;
365  }
366  }
367  }
368  }
369  }else if(optForNorm==5){
370  hRatio->Fit("pol0","","",minMass,minMass+rangeForNorm);
371  TF1* func0=(TF1*)hRatio->GetListOfFunctions()->FindObject("pol0");
372  Double_t norml=func0->GetParameter(0);
373  hRatio->Fit("pol0","","",maxMass-rangeForNorm,maxMass);
374  func0=(TF1*)hRatio->GetListOfFunctions()->FindObject("pol0");
375  Double_t normh=func0->GetParameter(0);
376  if(norml>normh){
377  norm=norml;
378  massForNorm=maxMass-0.5*rangeForNorm;
379  }else{
380  norm=normh;
381  massForNorm=minMass+0.5*rangeForNorm;
382  }
383  }
384 
385  printf("Normalization factor = %f --> taken at mass=%f\n",norm,massForNorm);
386 
387  return norm;
388 }
389 
390 void ProjectCombinHFAndFit(){
391 
392 
393 
394  TString dirName=Form("PWG3_D2H_InvMass%sLowPt%s",meson.Data(),suffix.Data());
395  TString lstName=Form("coutput%s%s",meson.Data(),suffix.Data());
396  TString dirNameMC=Form("PWG3_D2H_InvMass%sLowPt%s",meson.Data(),suffixMC.Data());
397  TString lstNameMC=Form("coutput%s%s",meson.Data(),suffixMC.Data());
398 
399  if(correctForRefl) suffix.Prepend("Refl_");
400  if(fileName.Contains("FAST") && !fileName.Contains("wSDD")){
401  suffix.Prepend("FAST_");
402  }else if(!fileName.Contains("FAST") && fileName.Contains("wSDD")){
403  suffix.Prepend("wSDD_");
404  }
405  if(fileNameMC.Contains("_G3")) suffix.Append("_Geant3MC");
406  else if(fileNameMC.Contains("_G4")) suffix.Append("_Geant4MC");
407  if(smoothLS!=0) suffix.Append(Form("_smoothLS%d",smoothLS));
408 
409  if(massD<0){
410  if(meson=="Dplus") massD=TDatabasePDG::Instance()->GetParticle(411)->Mass();
411  else if(meson=="Dzero") massD=TDatabasePDG::Instance()->GetParticle(421)->Mass();
412  }
413 
414  Int_t nBinsWithFixSig=0;
415  Int_t nBinsWithFixMean=0;
416  Int_t patSig=0;
417  Int_t patMean=0;
418  for(Int_t iPtBin=0; iPtBin<nPtBins; iPtBin++){
419  if(fixSigmaConf==1 || (fixSigmaConf==2 && fixSigma[iPtBin])){
420  nBinsWithFixSig++;
421  patSig+=1<<iPtBin;
422  }
423  if(fixMeanConf==1 || (fixMeanConf==2 && fixMean[iPtBin])){
424  nBinsWithFixMean++;
425  patMean+=1<<iPtBin;
426  }
427  }
428 
429  TFile* fil=new TFile(fileName.Data());
430  TDirectoryFile* df=(TDirectoryFile*)fil->Get(dirName.Data());
431 
432 
433  AliNormalizationCounter *nC=(AliNormalizationCounter*)df->Get("NormalizationCounter");
434  TH1F* hnEv=new TH1F("hEvForNorm","events for normalization",1,0,1);
435  printf("Number of Ev. for norm=%d\n",(Int_t)nC->GetNEventsForNorm());
436  hnEv->SetBinContent(1,nC->GetNEventsForNorm());
437 
438  TH1F* hRebin=new TH1F("hRebin","",nPtBins,binLims);
439  TH1F* hBkgFitFunc=new TH1F("hBkgFitFunc","",nPtBins,binLims);
440  TH1F* hBkgFitFuncSB=new TH1F("hBkgFitFuncSB","",nPtBins,binLims);
441 
442  TList* l=(TList*)df->Get(lstName.Data());
443  l->ls();
444 
445  TH3F* h3d=(TH3F*)l->FindObject("hMassVsPtVsY");
446  TH3F* h3dr=(TH3F*)l->FindObject("hMassVsPtVsYRot");
447  TH3F* h3dme=(TH3F*)l->FindObject("hMassVsPtVsYME");
448  TH3F* h3dmepp=(TH3F*)l->FindObject("hMassVsPtVsYMELSpp");
449  TH3F* h3dmemm=(TH3F*)l->FindObject("hMassVsPtVsYMELSmm");
450  TH2F* hpoolMix=(TH2F*)l->FindObject("hMixingsPerPool");
451  TH2F* hpoolEv=(TH2F*)l->FindObject("hEventsPerPool");
452  TH3F* h3dlsp=(TH3F*)l->FindObject("hMassVsPtVsYLSpp");
453  TH3F* h3dlsm=(TH3F*)l->FindObject("hMassVsPtVsYLSmm");
454 
455  TH3F* h3drefl=0x0;
456  TH3F* h3dmcsig=0x0;
457  TH1F* hSigmaMC=0x0;
458  TFile* filMC=new TFile(fileNameMC.Data());
459  if(filMC && filMC->IsOpen()){
460  TDirectoryFile* dfMC=(TDirectoryFile*)filMC->Get(dirNameMC.Data());
461  TList* lMC=(TList*)dfMC->Get(lstNameMC.Data());
462  hSigmaMC=FitMCInvMassSpectra(lMC);
463  if(nBinsWithFixSig>0 && !hSigmaMC){
464  printf("Fit to MC inv. mass spectra failed\n");
465  return;
466  }
467  if(correctForRefl){
468  h3drefl=(TH3F*)lMC->FindObject("hMassVsPtVsYRefl");
469  h3dmcsig=(TH3F*)lMC->FindObject("hMassVsPtVsYSig");
470  }
471  }
472 
473  TString sigConf="FixedSigma";
474  if(nBinsWithFixSig==0) sigConf="FreeSigma";
475  else if(nBinsWithFixSig==nPtBins) sigConf="FixedSigmaAll";
476  else sigConf=Form("FixedSigma%d",patSig);
477  if(nBinsWithFixSig>0 && tuneSigmaOnData>0.) sigConf+=Form("%d",TMath::Nint(tuneSigmaOnData*100.));
478  if(nBinsWithFixMean==nPtBins) sigConf+="FixedMeanAll";
479  else if(nBinsWithFixMean>0) sigConf+=Form("FixedMean%d",patMean);
480 
481  TCanvas* cem=new TCanvas("cem","Pools",1200,600);
482  cem->Divide(2,1);
483  cem->cd(1);
484  gPad->SetRightMargin(0.12);
485  gPad->SetLeftMargin(0.12);
486  hpoolMix->GetXaxis()->SetTitle("zVertex");
487  hpoolMix->GetYaxis()->SetTitle("Multiplicity");
488  hpoolMix->GetYaxis()->SetTitleOffset(1.4);
489  hpoolMix->Draw("colztext");
490  cem->cd(2);
491  gPad->SetRightMargin(0.12);
492  gPad->SetLeftMargin(0.12);
493  hpoolEv->Draw("colztext");
494  hpoolEv->GetXaxis()->SetTitle("zVertex");
495  hpoolEv->GetYaxis()->SetTitle("Multiplicity");
496  hpoolEv->GetYaxis()->SetTitleOffset(1.4);
497 
498 
499  TCanvas* c1=new TCanvas("c1","Mass",1200,800);
500  DivideCanvas(c1,nPtBins);
501 
502  TCanvas* c2=new TCanvas("c2","Mass-Bkg Rot",1200,800);
503  DivideCanvas(c2,nPtBins);
504  TCanvas* c2pulls=new TCanvas("c2pulls","Mass-Bkg Rot pulls",1200,800);
505  DivideCanvas(c2pulls,nPtBins);
506  TCanvas* c2residuals=new TCanvas("c2residuals","Mass-Bkg Rot residuals",1200,800);
507  DivideCanvas(c2residuals,nPtBins);
508  TCanvas* c2residualTrend=new TCanvas("c2residualTrend","Mass-Bkg Rot residual trend vs. mass",1200,800);
509  DivideCanvas(c2residualTrend,nPtBins);
510  TCanvas* c2pullTrend=new TCanvas("c2pullTrend","Mass-Bkg Rot pull trend vs. mass",1200,800);
511  DivideCanvas(c2pullTrend,nPtBins);
512 
513  TCanvas* c3=new TCanvas("c3","Mass-Bkg LS",1200,800);
514  DivideCanvas(c3,nPtBins);
515  TCanvas* c3pulls=new TCanvas("c3pulls","Mass-Bkg LS pulls",1200,800);
516  DivideCanvas(c3pulls,nPtBins);
517  TCanvas* c3residuals=new TCanvas("c3residuals","Mass-Bkg LS residuals",1200,800);
518  DivideCanvas(c3residuals,nPtBins);
519  TCanvas* c3residualTrend=new TCanvas("c3residualTrend","Mass-Bkg LS residual trend vs. mass",1200,800);
520  DivideCanvas(c3residualTrend,nPtBins);
521  TCanvas* c3pullTrend=new TCanvas("c3pullTrend","Mass-Bkg LS pull trend vs. mass",1200,800);
522  DivideCanvas(c3pullTrend,nPtBins);
523 
524  TCanvas* c4=new TCanvas("c4","Mass-Bkg EM",1200,800);
525  DivideCanvas(c4,nPtBins);
526  TCanvas* c4pulls=new TCanvas("c4pulls","Mass-Bkg EM pulls",1200,800);
527  DivideCanvas(c4pulls,nPtBins);
528  TCanvas* c4residuals=new TCanvas("c4residuals","Mass-Bkg EM residuals",1200,800);
529  DivideCanvas(c4residuals,nPtBins);
530  TCanvas* c4residualTrend=new TCanvas("c4residualTrend","Mass-Bkg EM residual trend vs. mass",1200,800);
531  DivideCanvas(c4residualTrend,nPtBins);
532  TCanvas* c4pullTrend=new TCanvas("c4pullTrend","Mass-Bkg EM pull trend vs. mass",1200,800);
533  DivideCanvas(c4pullTrend,nPtBins);
534 
535  TCanvas *c5=0x0;
536  TCanvas *c5sub=0x0;
537  TCanvas *c5pulls=0x0;
538  TCanvas *c5residuals=0x0;
539  TCanvas *c5residualTrend=0x0;
540  TCanvas *c5pullTrend=0x0;
541 
542  if(tryDirectFit){
543  c5=new TCanvas("c5","Mass SB Fit",1200,800);
544  DivideCanvas(c5,nPtBins);
545  c5sub=new TCanvas("c5sub","Mass-Bkg SB",1200,800);
546  DivideCanvas(c5sub,nPtBins);
547  c5pulls=new TCanvas("c5pulls","Mass-Bkg SB pulls",1200,800);
548  DivideCanvas(c5pulls,nPtBins);
549  c5residuals=new TCanvas("c5residuals","Mass-Bkg SB residuals",1200,800);
550  DivideCanvas(c5residuals,nPtBins);
551  c5residualTrend=new TCanvas("c5residualTrend","Mass-Bkg SB residual trend vs. mass",1200,800);
552  DivideCanvas(c5residualTrend,nPtBins);
553  c5pullTrend=new TCanvas("c5pullTrend","Mass-Bkg SB pull trend vs. mass",1200,800);
554  DivideCanvas(c5pullTrend,nPtBins);
555  }
556 
557  TCanvas *cCompareResidualTrends=new TCanvas("cCompareResidualTrends","cCompareResidualTrends",1200,800);
558  DivideCanvas(cCompareResidualTrends,nPtBins);
559 
560 
561  AliHFInvMassFitter* fitterRot[nPtBins];
562  AliHFInvMassFitter* fitterLS[nPtBins];
563  AliHFInvMassFitter* fitterME[nPtBins];
564  AliHFInvMassFitter* fitterSB[nPtBins];
565 
566  TH1F* hRawYieldRot=new TH1F("hRawYieldRot"," ; p_{T} (GeV/c) ; Raw yield",nPtBins,binLims);
567  TH1F* hRawYieldLS=new TH1F("hRawYieldLS"," ; p_{T} (GeV/c) ; Raw yield",nPtBins,binLims);
568  TH1F* hRawYieldME=new TH1F("hRawYieldME"," ; p_{T} (GeV/c) ; Raw yield",nPtBins,binLims);
569  TH1F* hRawYieldSB=new TH1F("hRawYieldSB"," ; p_{T} (GeV/c) ; Raw yield",nPtBins,binLims);
570 
571  TH1F* hRelStatRot=new TH1F("hRelStatRot"," ; p_{T} (GeV/c) ; Relative stat. unc.",nPtBins,binLims);
572  TH1F* hRelStatLS=new TH1F("hRelStatLS"," ; p_{T} (GeV/c) ; Relative stat. unc.",nPtBins,binLims);
573  TH1F* hRelStatME=new TH1F("hRelStatME"," ; p_{T} (GeV/c) ; Relative stat. unc.",nPtBins,binLims);
574  TH1F* hRelStatSB=new TH1F("hRelStatSB"," ; p_{T} (GeV/c) ; Relative stat. unc.",nPtBins,binLims);
575 
576  TH1F* hSignifRot=new TH1F("hSignifRot","",nPtBins,binLims);
577  TH1F* hSignifLS=new TH1F("hSignifLS","",nPtBins,binLims);
578  TH1F* hSignifME=new TH1F("hSignifME","",nPtBins,binLims);
579  TH1F* hSignifSB=new TH1F("hSignifSB","",nPtBins,binLims);
580 
581  TH1F* hSoverBRot=new TH1F("hSoverBRot","",nPtBins,binLims);
582  TH1F* hSoverBLS=new TH1F("hSoverBLS","",nPtBins,binLims);
583  TH1F* hSoverBME=new TH1F("hSoverBME","",nPtBins,binLims);
584  TH1F* hSoverBSB=new TH1F("hSoverBSB","",nPtBins,binLims);
585 
586  TH1F* hGausMeanRot=new TH1F("hGausMeanRot","",nPtBins,binLims);
587  TH1F* hGausMeanLS=new TH1F("hGausMeanLS","",nPtBins,binLims);
588  TH1F* hGausMeanME=new TH1F("hGausMeanME","",nPtBins,binLims);
589  TH1F* hGausMeanSB=new TH1F("hGausMeanSB","",nPtBins,binLims);
590 
591  TH1F* hGausSigmaRot=new TH1F("hGausSigmaRot","",nPtBins,binLims);
592  TH1F* hGausSigmaLS=new TH1F("hGausSigmaLS","",nPtBins,binLims);
593  TH1F* hGausSigmaME=new TH1F("hGausSigmaME","",nPtBins,binLims);
594  TH1F* hGausSigmaSB=new TH1F("hGausSigmaSB","",nPtBins,binLims);
595 
596  TH1F* hChiSqRot=new TH1F("hChiSqRot","",nPtBins,binLims);
597  TH1F* hChiSqLS=new TH1F("hChiSqLS","",nPtBins,binLims);
598  TH1F* hChiSqME=new TH1F("hChiSqME","",nPtBins,binLims);
599  TH1F* hChiSqSB=new TH1F("hChiSqSB","",nPtBins,binLims);
600 
601  TH1F* hNdfRot=new TH1F("hNdfRot","",nPtBins,binLims);
602  TH1F* hNdfLS=new TH1F("hNdfLS","",nPtBins,binLims);
603  TH1F* hNdfME=new TH1F("hNdfME","",nPtBins,binLims);
604  TH1F* hNdfSB=new TH1F("hNdfSB","",nPtBins,binLims);
605 
606  TH1F* hRawYieldRotBC=new TH1F("hRawYieldRotBC","BC yield (rotational background)",nPtBins,binLims);
607  TH1F* hRawYieldLSBC=new TH1F("hRawYieldLSBC","BC yield (like-sign background)",nPtBins,binLims);
608  TH1F* hRawYieldMEBC=new TH1F("hRawYieldMEBC","BC yield (mixed-event background)",nPtBins,binLims);
609  TH1F* hRawYieldSBBC=new TH1F("hRawYieldSBBC","BC yield (side-band fit background)",nPtBins,binLims);
610 
611  TH1F* hInvMassHistoBinWidthRot=new TH1F("hInvMassHistoBinWidthRot","",nPtBins,binLims);
612  TH1F* hInvMassHistoBinWidthLS=new TH1F("hInvMassHistoBinWidthLS","",nPtBins,binLims);
613  TH1F* hInvMassHistoBinWidthME=new TH1F("hInvMassHistoBinWidthME","",nPtBins,binLims);
614  TH1F* hInvMassHistoBinWidthSB=new TH1F("hInvMassHistoBinWidthSB","",nPtBins,binLims);
615 
616  TLatex* tME=new TLatex(0.65,0.82,"MixEv +- pairs");
617  tME->SetNDC();
618  TLatex* tMEpp=new TLatex(0.65,0.75,"MixEv ++ pairs");
619  tMEpp->SetNDC();
620  TLatex* tMEmm=new TLatex(0.65,0.68,"MixEv -- pairs");
621  tMEmm->SetNDC();
622 
623  TF1 *fpeak=new TF1("fpeak","[0]*1./(TMath::Sqrt(2.*TMath::Pi())*[2])*TMath::Exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",minMass,maxMass);
624 
625 
626  TDirectory *current = gDirectory;
627  TFile* fout=new TFile(Form("outputMassFits_%s_%s.root",sigConf.Data(),suffix.Data()),"recreate");
628  current->cd();
629 
630 
631  for(Int_t iPtBin=0; iPtBin<nPtBins; iPtBin++){
632 
633  Int_t bin1=h3d->GetYaxis()->FindBin(binLims[iPtBin]);
634  Int_t bin2=h3d->GetYaxis()->FindBin(binLims[iPtBin+1]-0.0001);
635  printf("Bin %d Pt range=%f %f\n",iPtBin,h3d->GetYaxis()->GetBinLowEdge(bin1),h3d->GetYaxis()->GetBinUpEdge(bin2));
636  TH1D* hMassPtBin=h3d->ProjectionX(Form("hMassPtBin%d",iPtBin),bin1,bin2);
637  TH1D* hMassPtBinr=h3dr->ProjectionX(Form("hMassPtBinr%d",iPtBin),bin1,bin2);
638  TH1D* hMassPtBinme=h3dme->ProjectionX(Form("hMassPtBinme%d",iPtBin),bin1,bin2);
639  TH1D* hMassPtBinmeLSpp=h3dmepp->ProjectionX(Form("hMassPtBinmeLSpp%d",iPtBin),bin1,bin2);
640  TH1D* hMassPtBinmeLSmm=h3dmemm->ProjectionX(Form("hMassPtBinmeLSmm%d",iPtBin),bin1,bin2);
641 
642  if(correctForRefl){
643  Int_t bin1MC=h3drefl->GetYaxis()->FindBin(binLims[iPtBin]);
644  Int_t bin2MC=h3drefl->GetYaxis()->FindBin(binLims[iPtBin+1]-0.0001);
645  hMCReflPtBin=h3drefl->ProjectionX(Form("hMCReflPtBin%d",iPtBin),bin1MC,bin2MC);
646  hMCSigPtBin=h3dmcsig->ProjectionX(Form("hMCSigPtBin%d",iPtBin),bin1MC,bin2MC);
647  }
648 
649  TH1D* hMassPtBinlsp=0x0;
650  TH1D* hMassPtBinlsm=0x0;
651  TH1D* hMassPtBinls=0x0;
652  if(h3dlsp){
653  hMassPtBinlsp=h3dlsp->ProjectionX(Form("hMassPtBinlsp%d",iPtBin),bin1,bin2);
654  hMassPtBinlsm=h3dlsm->ProjectionX(Form("hMassPtBinlsm%d",iPtBin),bin1,bin2);
655  hMassPtBinls=(TH1D*)hMassPtBinlsp->Clone(Form("hMassPtBinls%d",iPtBin));
656  hMassPtBinls->Reset();
657  for(Int_t iBin=1; iBin<=hMassPtBinlsp->GetNbinsX(); iBin++){
658  Double_t np=hMassPtBinlsp->GetBinContent(iBin);
659  Double_t nm=hMassPtBinlsm->GetBinContent(iBin);
660  Double_t tt=2*TMath::Sqrt(np*nm);
661  Double_t enp=hMassPtBinlsp->GetBinError(iBin);
662  Double_t enm=hMassPtBinlsm->GetBinError(iBin);
663  Double_t ett=0;
664  if(tt>0) ett=2./tt*TMath::Sqrt(np*np*enm*enm+nm*nm*enp*enp);
665  hMassPtBinls->SetBinContent(iBin,tt);
666  hMassPtBinls->SetBinError(iBin,ett);
667  }
668  if(smoothLS<-0.5)hMassPtBinls->Smooth(-1*smoothLS);
669  else if(smoothLS>0.5)QuadraticSmooth(hMassPtBinls,smoothLS);
670  hMassPtBinls->SetLineColor(kGreen+1);
671  }
672 
673  // hMassPtBin->Sumw2();
674  // hMassPtBinr->Sumw2();
675  // hMassPtBinme->Sumw2();
676  // hMassPtBinlsp->Sumw2();
677  // hMassPtBinls->Sumw2();
678  hMassPtBin->SetLineColor(1);
679  hMassPtBinr->SetLineColor(4);
680  hMassPtBinme->SetLineColor(kOrange+1);
681  hMassPtBinmeLSpp->SetLineColor(kGreen+2);
682  hMassPtBinmeLSmm->SetLineColor(kCyan);
683  hMassPtBin->SetTitle(Form("%.1f<p_{T}<%.1f GeV/c",binLims[iPtBin],binLims[iPtBin+1]));
684  hMassPtBinr->SetTitle(Form("%.1f<p_{T}<%.1f GeV/c",binLims[iPtBin],binLims[iPtBin+1]));
685  hMassPtBinme->SetTitle(Form("%.1f<p_{T}<%.1f GeV/c",binLims[iPtBin],binLims[iPtBin+1]));
686  TH1D* hRatioMEpp=(TH1D*)hMassPtBinmeLSpp->Clone(Form("hRatioPtBinmeLSpp%d",iPtBin));
687  hRatioMEpp->Divide(hMassPtBinme);
688  TH1D* hRatioMEmm=(TH1D*)hMassPtBinmeLSmm->Clone(Form("hRatioPtBinmeLSmm%d",iPtBin));
689  hRatioMEmm->Divide(hMassPtBinme);
690  TH1D* hMassPtBinmeAll=(TH1D*)hMassPtBinme->Clone(Form("hRatioPtBinmeAll%d",iPtBin));
691  hMassPtBinmeAll->Add(hMassPtBinmeLSpp);
692  hMassPtBinmeAll->Add(hMassPtBinmeLSmm);
693  hMassPtBinmeAll->SetLineColor(kRed);
694  TH1D* hRatioME=(TH1D*)hMassPtBinme->Clone(Form("hRatioME%d",iPtBin));
695  hRatioME->Divide(hMassPtBin);
696  TH1D* hRatioMEAll=(TH1D*)hMassPtBinmeAll->Clone(Form("hRatioMEAll%d",iPtBin));
697  hRatioMEAll->Divide(hMassPtBin);
698 
699 
700  TCanvas* c0=new TCanvas(Form("CBin%d",iPtBin),Form("Bin%d norm",iPtBin),1000,700);
701  c0->Divide(2,2);
702  c0->cd(1);
703  TH1D* hRatio=(TH1D*)hMassPtBinr->Clone(Form("hRatioFormNorm%d",iPtBin));
704  hRatio->Divide(hMassPtBin);
705  hRatio->Draw();
706  hRatio->GetYaxis()->SetTitle("Rotational/All");
707  hRatio->GetXaxis()->SetTitle("Invariant mass (GeV/c^{2})");
708  Double_t normRot=GetBackgroundNormalizationFactor(hRatio,rebin[iPtBin]);
709  hMassPtBinr->Scale(1./normRot);
710  TLatex* tnr=new TLatex(0.2,0.2,Form("Normaliz. factor = %f",normRot));
711  tnr->SetNDC();
712  tnr->Draw();
713  c0->cd(2);
714  hMassPtBinme->GetYaxis()->SetTitle("Entries (EvMix)");
715  hMassPtBinme->GetXaxis()->SetTitle("Invariant mass (GeV/c^{2})");
716  hMassPtBinme->DrawCopy();
717  hMassPtBinmeLSpp->Draw("same");
718  hMassPtBinmeLSmm->Draw("same");
719  tME->SetTextColor(hMassPtBinme->GetLineColor());
720  tMEpp->SetTextColor(hMassPtBinmeLSpp->GetLineColor());
721  tMEmm->SetTextColor(hMassPtBinmeLSmm->GetLineColor());
722  tME->Draw();
723  tMEpp->Draw();
724  tMEmm->Draw();
725  c0->cd(4);
726  hRatioMEpp->Draw();
727  hRatioMEpp->SetMinimum(0.4);
728  hRatioMEpp->SetMaximum(0.6);
729  hRatioMEpp->GetYaxis()->SetTitle("ME with LS / ME with OS");
730  hRatioMEpp->GetXaxis()->SetTitle("Invariant mass (GeV/c^{2})");
731  hRatioMEmm->Draw("same");
732  c0->cd(3);
733  hRatioME->Draw();
734  hRatioME->SetMaximum(hRatioMEAll->GetMaximum()*1.05);
735  hRatioMEAll->Draw("same");
736  hRatioME->GetYaxis()->SetTitle("EvMix/All");
737  hRatioME->GetXaxis()->SetTitle("Invariant mass (GeV/c^{2})");
738 
739  Double_t normME=GetBackgroundNormalizationFactor(hRatioME,rebin[iPtBin]);
740  Double_t normMEAll=GetBackgroundNormalizationFactor(hRatioMEAll,rebin[iPtBin]);
741  hMassPtBinme->Scale(1./normME);
742  hMassPtBinmeAll->Scale(1./normMEAll);
743  printf("Background norm bin %d DONE\n",iPtBin);
744 
745  c1->cd(iPtBin+1);
746  hMassPtBin->GetXaxis()->SetRangeUser(minMass,maxMass);
747  hMassPtBin->SetMinimum(0);
748  hMassPtBin->Draw();
749  hMassPtBin->GetYaxis()->SetTitle("Counts");
750  hMassPtBin->GetYaxis()->SetTitleOffset(2.);
751  hMassPtBin->GetXaxis()->SetTitle("Invariant mass (GeV/c^{2})");
752  hMassPtBinr->Draw("same");
753  if(!useEMwithLS) hMassPtBinme->Draw("same");
754  else hMassPtBinmeAll->Draw("same");
755  if(hMassPtBinls) hMassPtBinls->Draw("same");
756  if(iPtBin==0){
757  TLegend* leg=new TLegend(0.5,0.6,0.89,0.89);
758  leg->SetFillStyle(0);
759  leg->AddEntry(hMassPtBin,"All candidates","L")->SetTextColor(hMassPtBin->GetLineColor());
760  leg->AddEntry(hMassPtBinr,"Background (rotations)","L")->SetTextColor(hMassPtBinr->GetLineColor());
761  if(!useEMwithLS) leg->AddEntry(hMassPtBinme,"Background (ME)","L")->SetTextColor(hMassPtBinme->GetLineColor());
762  else leg->AddEntry(hMassPtBinmeAll,"Background (ME)","L")->SetTextColor(hMassPtBinmeAll->GetLineColor());
763  if(hMassPtBinls) leg->AddEntry(hMassPtBinls,"Like-sign","L")->SetTextColor(hMassPtBinls->GetLineColor());
764  leg->Draw();
765  }
766  gPad->Update();
767 
768  TH1D* hMassSubRot=(TH1D*)hMassPtBin->Clone(Form("hMassSubRot_bin%d",iPtBin));
769  hMassSubRot->Add(hMassPtBinr,-1);
770  hMassSubRot->SetTitle(Form("%.1f<p_{T}<%.1f GeV/c -- Rotational",binLims[iPtBin],binLims[iPtBin+1]));
771  TH1D* hMassSubME=(TH1D*)hMassPtBin->Clone(Form("hMassSubME_bin%d",iPtBin));
772  if(!useEMwithLS) hMassSubME->Add(hMassPtBinme,-1);
773  else hMassSubME->Add(hMassPtBinmeAll,-1);
774  hMassSubME->SetTitle(Form("%.1f<p_{T}<%.1f GeV/c -- Mixed Ev",binLims[iPtBin],binLims[iPtBin+1]));
775  TH1D* hMassSubLS=0x0;
776  if(hMassPtBinls){
777  hMassSubLS=(TH1D*)hMassPtBin->Clone(Form("hMassSubLS_bin%d",iPtBin));
778  hMassSubLS->Add(hMassPtBinls,-1);
779  hMassSubLS->SetTitle(Form("%.1f<p_{T}<%.1f GeV/c -- Like Sign",binLims[iPtBin],binLims[iPtBin+1]));
780  }
781 
782  fout->cd();
783  hMassPtBin->Write();
784  hMassSubRot->Write();
785  hMassSubME->Write();
786  if(hMassPtBinls) hMassSubLS->Write();
787  if(hMCReflPtBin) hMCReflPtBin->Write();
788  if(hMCSigPtBin) hMCSigPtBin->Write();
789  current->cd();
790 
791  hMassSubRot=AliVertexingHFUtils::RebinHisto(hMassSubRot,rebin[iPtBin]);
792  hMassSubME=AliVertexingHFUtils::RebinHisto(hMassSubME,rebin[iPtBin]);
793  if(hMassPtBinls) hMassSubLS=AliVertexingHFUtils::RebinHisto(hMassSubLS,rebin[iPtBin]);
794 
795  hRebin->SetBinContent(iPtBin+1,rebin[iPtBin]);
796  Int_t bkgToFill=typeb;
797  if(typeb==6) bkgToFill=typeb+nDegreeBackPol[iPtBin];
798  Int_t bkgToFillSB=6+nDegreeBackPol[iPtBin];
799 
800  hBkgFitFunc->SetBinContent(iPtBin+1,bkgToFill);
801  hBkgFitFuncSB->SetBinContent(iPtBin+1,bkgToFillSB);
802 
803  fitterRot[iPtBin]=ConfigureFitter(hMassSubRot,iPtBin,typeb,minMass,maxMass);
804  if(hMassPtBinls) fitterLS[iPtBin]=ConfigureFitter(hMassSubLS,iPtBin,typeb,minMass,maxMass);
805  fitterME[iPtBin]=ConfigureFitter(hMassSubME,iPtBin,typeb,minMass,maxMass);
806 
807  Bool_t out1=fitterRot[iPtBin]->MassFitter(0);
808  Bool_t out2=kFALSE;
809  if(hMassPtBinls) out2=fitterLS[iPtBin]->MassFitter(0);
810  Bool_t out3=fitterME[iPtBin]->MassFitter(0);
811 
812  Double_t background=999999999.;
813  Bool_t out4=kFALSE;
814  if(tryDirectFit){
815  TH1D *hMassDirectFit=(TH1D*)hMassPtBin->Clone(Form("hMassDirectFit_bin%d",iPtBin));
816  hMassDirectFit=AliVertexingHFUtils::RebinHisto(hMassDirectFit,rebin[iPtBin]);
817  fitterSB[iPtBin]=ConfigureFitter(hMassDirectFit,iPtBin,6,fitrangelow[iPtBin],fitrangeup[iPtBin]);
818  out4=fitterSB[iPtBin]->MassFitter(0);//DirectFit(hMassDirectFit,iPtBin,hRawYieldSB);
819 
820  Double_t background,ebkg;
821 
822  if(out4 && fitterSB[iPtBin]->GetMassFunc()){
823  c5->cd(iPtBin+1);
824  fitterSB[iPtBin]->DrawHere(gPad,3,0);
825  hRawYieldSB->SetBinContent(iPtBin+1,fitterSB[iPtBin]->GetRawYield());
826  hRawYieldSB->SetBinError(iPtBin+1,fitterSB[iPtBin]->GetRawYieldError());
827  if(fitterSB[iPtBin]->GetRawYield()>0){
828  hRelStatSB->SetBinContent(iPtBin+1,fitterSB[iPtBin]->GetRawYieldError()/fitterSB[iPtBin]->GetRawYield());
829  hRelStatSB->SetBinError(iPtBin+1,0.00000001);
830  }
831  fitterSB[iPtBin]->Background(3.,background,ebkg);
832  hSignifSB->SetBinContent(iPtBin+1,fitterSB[iPtBin]->GetRawYield()/TMath::Sqrt(background+fitterSB[iPtBin]->GetRawYield()));
833  hSignifSB->SetBinError(iPtBin+1,0.00000001);
834  hSoverBSB->SetBinContent(iPtBin+1,fitterSB[iPtBin]->GetRawYield()/background);
835  hSoverBSB->SetBinError(iPtBin+1,0.00000001);
836  hInvMassHistoBinWidthSB->SetBinContent(iPtBin+1,hMassDirectFit->GetBinWidth(1));
837  hGausMeanSB->SetBinContent(iPtBin+1,fitterSB[iPtBin]->GetMean());
838  hGausMeanSB->SetBinError(iPtBin+1,fitterSB[iPtBin]->GetMeanUncertainty());
839  hGausSigmaSB->SetBinContent(iPtBin+1,fitterSB[iPtBin]->GetSigma());
840  hGausSigmaSB->SetBinError(iPtBin+1,fitterSB[iPtBin]->GetSigmaUncertainty());
841  hChiSqSB->SetBinContent(iPtBin+1,fitterSB[iPtBin]->GetReducedChiSquare());
842  hChiSqSB->SetBinError(iPtBin+1,0.00001); // very small number, for graphics
843  hNdfSB->SetBinContent(iPtBin+1,fitterSB[iPtBin]->GetMassFunc()->GetNDF());
844  hNdfSB->SetBinError(iPtBin+1,0.00001); // very small number, for graphics
845  // if(!correctForRefl)
846  WriteFitInfo(fitterSB[iPtBin],hMassPtBin);
847  fout->cd();
848  WriteFitFunctionsToFile(fitterSB[iPtBin],"SB",iPtBin);
849  current->cd();
850 
851  c5sub->cd(iPtBin+1);
852  TH1F* hsubTemp=(TH1F*)hMassDirectFit->Clone(Form("%sSubBack%d",hMassDirectFit->GetName(),iPtBin));
853  TH1F* hsubTempAllRange=(TH1F*)hMassDirectFit->Clone(Form("%sSubBackAllRange%d",hMassDirectFit->GetName(),iPtBin));
854  TF1* funcAll=fitterSB[iPtBin]->GetMassFunc();
855  TF1* funcBkg=fitterSB[iPtBin]->GetBackgroundRecalcFunc();
856  // TF1* funcBkg2=fitterSB[iPtBin]->GetBackgroundFullRangeFunc();
857  for(Int_t jst=1;jst<=hsubTemp->GetNbinsX();jst++){
858  Double_t backg=funcBkg->Integral(hsubTemp->GetBinLowEdge(jst),hsubTemp->GetBinLowEdge(jst)+hsubTemp->GetBinWidth(jst))/hsubTemp->GetBinWidth(jst);
859  Double_t tot=funcAll->Integral(hsubTempAllRange->GetBinLowEdge(jst),hsubTempAllRange->GetBinLowEdge(jst)+hsubTempAllRange->GetBinWidth(jst))/hsubTempAllRange->GetBinWidth(jst);
860  hsubTemp->SetBinContent(jst,hsubTemp->GetBinContent(jst)-backg);
861  hsubTempAllRange->SetBinContent(jst,hsubTempAllRange->GetBinContent(jst)-tot);
862  }
863  hsubTemp->SetLineColor(kBlue);
864  hsubTempAllRange->SetLineColor(kGray+2);
865  Double_t ymin=0;
866  Double_t ymax=1;
867  for(Int_t ibs=1; ibs<hsubTemp->GetNbinsX(); ibs++){
868  Double_t binc=hsubTemp->GetBinCenter(ibs);
869  if(binc>fitrangelow[iPtBin] && binc<fitrangeup[iPtBin]){
870  Double_t yl=hsubTemp->GetBinContent(ibs)-hsubTemp->GetBinError(ibs);
871  Double_t yu=hsubTemp->GetBinContent(ibs)+hsubTemp->GetBinError(ibs);
872  if(yl<ymin) ymin=yl;
873  if(yu>ymax) ymax=yu;
874  }
875  }
876  if(ymax>0) ymax*=1.2;
877  else ymax*=0.8;
878  if(ymin<0) ymin*=1.2;
879  else ymin*=0.8;
880  hsubTemp->GetXaxis()->SetRangeUser(fitrangelow[iPtBin],fitrangeup[iPtBin]);
881  hsubTemp->SetMinimum(ymin);
882  hsubTemp->SetMaximum(ymax);
883  hsubTemp->SetMarkerStyle(20);
884  hsubTemp->SetMarkerColor(hsubTemp->GetLineColor());
885  hsubTemp->DrawCopy();
886  hsubTempAllRange->DrawCopy("same");
887  hsubTemp->DrawCopy("same");
888  fpeak->SetRange(fitrangelow[iPtBin],fitrangeup[iPtBin]);
889  fpeak->SetParameter(0,funcAll->GetParameter(nDegreeBackPol[iPtBin]+1));
890  fpeak->SetParameter(1,funcAll->GetParameter(nDegreeBackPol[iPtBin]+2));
891  fpeak->SetParameter(2,funcAll->GetParameter(nDegreeBackPol[iPtBin]+3));
892  fpeak->DrawCopy("same");
893 
894  Double_t errbc;
895  Double_t bc=fitterSB[iPtBin]->GetRawYieldBinCounting(errbc,nsigmaBinCounting,optBkgBinCount);
896  hRawYieldSBBC->SetBinContent(iPtBin+1,bc);
897  hRawYieldSBBC->SetBinError(iPtBin+1,errbc);
898 
899  c5pulls->cd(iPtBin+1);
900  TH1F *hPulls=new TH1F();//the name is changed afterwards, histo must not be deleted
901  TH1F *hPullsTrend=new TH1F();// the name is changed afterwards, histo must not be deleted
902  TH1F *hResidualTrend=new TH1F();// the name is changed afterwards, histo must not be deleted
903  TH1F *hResiduals=fitterSB[iPtBin]->GetResidualsAndPulls(hPulls,hResidualTrend,hPullsTrend);
904  hResiduals->SetName(Form("hResidualsSB_PtBin%d",iPtBin));
905 
906  hPulls->Draw();
907  PrintGausParams(hPulls);
908 
909  c5residuals->cd(iPtBin+1);
910  hResiduals->Draw();
911 
912  SetStyleHisto(hResidualTrend,kSB);
913  cCompareResidualTrends->cd(iPtBin+1);
914  hResidualTrend->Draw();
915 
916  c5residualTrend->cd(iPtBin+1);
917  hResidualTrend->Draw();
918 
919  c5pullTrend->cd(iPtBin+1);
920  hPullsTrend->Draw();
921  //TPaveStats *tp=(TPaveStats*)hPulls->FindObject("stats");
922  // tp->SetOptStat("remnpcev");
923  delete hMassDirectFit;
924  }
925  }
926 
927  c2->cd(iPtBin+1);
928  if(out1 && fitterRot[iPtBin]->GetMassFunc()){
929  fitterRot[iPtBin]->DrawHere(gPad,3,0);
930  gPad->Update();
931  hRawYieldRot->SetBinContent(iPtBin+1,fitterRot[iPtBin]->GetRawYield());
932  hRawYieldRot->SetBinError(iPtBin+1,fitterRot[iPtBin]->GetRawYieldError());
933  Double_t minBinBkg=hMassPtBin->FindBin(fitterRot[iPtBin]->GetMean()-3.*fitterRot[iPtBin]->GetSigma());
934  Double_t maxBinBkg=hMassPtBin->FindBin(fitterRot[iPtBin]->GetMean()+3.*fitterRot[iPtBin]->GetSigma());
935  background=hMassPtBin->Integral(minBinBkg,maxBinBkg);
936  hSoverBRot->SetBinContent(iPtBin+1,fitterRot[iPtBin]->GetRawYield()/background);
937  hSoverBRot->SetBinError(iPtBin+1,0.000001);
938  hRelStatRot->SetBinContent(iPtBin+1,fitterRot[iPtBin]->GetRawYieldError()/fitterRot[iPtBin]->GetRawYield());
939  hRelStatRot->SetBinError(iPtBin+1,0.000001);
940  hSignifRot->SetBinContent(iPtBin+1,fitterRot[iPtBin]->GetRawYield()/TMath::Sqrt(background+fitterRot[iPtBin]->GetRawYield()));
941  hSignifRot->SetBinError(iPtBin+1,0.00000001);
942  hInvMassHistoBinWidthRot->SetBinContent(iPtBin+1,hMassSubRot->GetBinWidth(1));
943  hGausMeanRot->SetBinContent(iPtBin+1,fitterRot[iPtBin]->GetMean());
944  hGausMeanRot->SetBinError(iPtBin+1,fitterRot[iPtBin]->GetMeanUncertainty());
945  hGausSigmaRot->SetBinContent(iPtBin+1,fitterRot[iPtBin]->GetSigma());
946  hGausSigmaRot->SetBinError(iPtBin+1,fitterRot[iPtBin]->GetSigmaUncertainty());
947  hChiSqRot->SetBinContent(iPtBin+1,fitterRot[iPtBin]->GetReducedChiSquare());
948  hChiSqRot->SetBinError(iPtBin+1,0.00001); // very small number, for graphics
949  hNdfRot->SetBinContent(iPtBin+1,fitterRot[iPtBin]->GetMassFunc()->GetNDF());
950  hNdfRot->SetBinError(iPtBin+1,0.00001); // very small number, for graphics
951  // if(!correctForRefl)
952  WriteFitInfo(fitterRot[iPtBin],hMassPtBin);
953  fout->cd();
954  WriteFitFunctionsToFile(fitterRot[iPtBin],"Rot",iPtBin);
955  current->cd();
956 
957  Double_t errbc;
958  Double_t bc=fitterRot[iPtBin]->GetRawYieldBinCounting(errbc,nsigmaBinCounting,optBkgBinCount);
959  hRawYieldRotBC->SetBinContent(iPtBin+1,bc);
960  hRawYieldRotBC->SetBinError(iPtBin+1,errbc);
961 
962  c2pulls->cd(iPtBin+1);
963  TH1F *hPullsTrend=new TH1F();// the name is changed afterwards, histo must not be deleted
964  TH1F *hPulls=new TH1F();//the name is changed afterwards, histo must not be deleted
965  TH1F *hResidualTrend=new TH1F();// the name is changed afterwards, histo must not be deleted
966  TH1F *hResiduals=fitterRot[iPtBin]->GetResidualsAndPulls(hPulls,hResidualTrend,hPullsTrend);
967  hResiduals->SetName(Form("hResidualsRot_PtBin%d",iPtBin));
968 
969  hPulls->Draw();
970  PrintGausParams(hPulls);
971 
972  c2residuals->cd(iPtBin+1);
973  hResiduals->Draw();
974 
975  SetStyleHisto(hResidualTrend,kRot);
976  c2residualTrend->cd(iPtBin+1);
977  hResidualTrend->Draw();
978  cCompareResidualTrends->cd(iPtBin+1);
979  if(out4)hResidualTrend->Draw("same");
980  else hResidualTrend->Draw();
981 
982  c2pullTrend->cd(iPtBin+1);
983  hPullsTrend->Draw();
984 
985  }
986  else hMassSubRot->Draw("");
987 
988 
989  c3->cd(iPtBin+1);
990  if(out2 && fitterLS[iPtBin]->GetMassFunc()){
991  fitterLS[iPtBin]->DrawHere(gPad,3,0);
992  hRawYieldLS->SetBinContent(iPtBin+1,fitterLS[iPtBin]->GetRawYield());
993  hRawYieldLS->SetBinError(iPtBin+1,fitterLS[iPtBin]->GetRawYieldError());
994  Double_t minBinBkg=hMassPtBin->FindBin(fitterLS[iPtBin]->GetMean()-3.*fitterLS[iPtBin]->GetSigma());
995  Double_t maxBinBkg=hMassPtBin->FindBin(fitterLS[iPtBin]->GetMean()+3.*fitterLS[iPtBin]->GetSigma());
996  background=hMassPtBin->Integral(minBinBkg,maxBinBkg);
997  hSoverBLS->SetBinContent(iPtBin+1,fitterLS[iPtBin]->GetRawYield()/background);
998  hSoverBLS->SetBinError(iPtBin+1,0.000001);
999  hRelStatLS->SetBinContent(iPtBin+1,fitterLS[iPtBin]->GetRawYieldError()/fitterLS[iPtBin]->GetRawYield());
1000  hRelStatLS->SetBinError(iPtBin+1,0.0000001);
1001  hSignifLS->SetBinContent(iPtBin+1,fitterLS[iPtBin]->GetRawYield()/TMath::Sqrt(background+fitterLS[iPtBin]->GetRawYield()));
1002  hSignifLS->SetBinError(iPtBin+1,0.00000001);
1003  hInvMassHistoBinWidthLS->SetBinContent(iPtBin+1,hMassSubLS->GetBinWidth(1));
1004  hGausMeanLS->SetBinContent(iPtBin+1,fitterLS[iPtBin]->GetMean());
1005  hGausMeanLS->SetBinError(iPtBin+1,fitterLS[iPtBin]->GetMeanUncertainty());
1006  hGausSigmaLS->SetBinContent(iPtBin+1,fitterLS[iPtBin]->GetSigma());
1007  hGausSigmaLS->SetBinError(iPtBin+1,fitterLS[iPtBin]->GetSigmaUncertainty());
1008  hChiSqLS->SetBinContent(iPtBin+1,fitterLS[iPtBin]->GetReducedChiSquare());
1009  hChiSqLS->SetBinError(iPtBin+1,0.00001); // very small number, for graphics
1010  hNdfLS->SetBinContent(iPtBin+1,fitterLS[iPtBin]->GetMassFunc()->GetNDF());
1011  hNdfLS->SetBinError(iPtBin+1,0.00001); // very small number, for graphics
1012  // if(!correctForRefl)
1013  WriteFitInfo(fitterLS[iPtBin],hMassPtBin);
1014  fout->cd();
1015  WriteFitFunctionsToFile(fitterLS[iPtBin],"LS",iPtBin);
1016  current->cd();
1017 
1018 
1019  Double_t errbc;
1020  Double_t bc=fitterLS[iPtBin]->GetRawYieldBinCounting(errbc,nsigmaBinCounting,optBkgBinCount);
1021  hRawYieldLSBC->SetBinContent(iPtBin+1,bc);
1022  hRawYieldLSBC->SetBinError(iPtBin+1,errbc);
1023 
1024  c3pulls->cd(iPtBin+1);
1025  TH1F *hPullsTrend=new TH1F();// the name is changed afterwards, histo must not be deleted
1026  TH1F *hPulls=new TH1F();//the name is changed afterwards, histo must not be deleted
1027  TH1F *hResidualTrend=new TH1F();// the name is changed afterwards, histo must not be deleted
1028  TH1F *hResiduals=fitterLS[iPtBin]->GetResidualsAndPulls(hPulls,hResidualTrend,hPullsTrend);
1029  hResiduals->SetName(Form("hResidualsLS_PtBin%d",iPtBin));
1030 
1031  hPulls->Draw();
1032  PrintGausParams(hPulls);
1033 
1034  c3residuals->cd(iPtBin+1);
1035  hResiduals->Draw();
1036 
1037  SetStyleHisto(hResidualTrend,kLS);
1038  c3residualTrend->cd(iPtBin+1);
1039  hResidualTrend->Draw();
1040 
1041  cCompareResidualTrends->cd(iPtBin+1);
1042  if(out4||out1)hResidualTrend->Draw("same");
1043  else hResidualTrend->Draw();
1044 
1045  c3pullTrend->cd(iPtBin+1);
1046  hPullsTrend->Draw();
1047  }
1048  else if(hMassPtBinls) hMassSubLS->Draw("");
1049 
1050  c4->cd(iPtBin+1);
1051  if(out3 && fitterME[iPtBin]->GetMassFunc()){
1052  fitterME[iPtBin]->DrawHere(gPad,3,0);
1053  hRawYieldME->SetBinContent(iPtBin+1,fitterME[iPtBin]->GetRawYield());
1054  hRawYieldME->SetBinError(iPtBin+1,fitterME[iPtBin]->GetRawYieldError());
1055  Double_t minBinBkg=hMassPtBin->FindBin(fitterME[iPtBin]->GetMean()-3.*fitterME[iPtBin]->GetSigma());
1056  Double_t maxBinBkg=hMassPtBin->FindBin(fitterME[iPtBin]->GetMean()+3.*fitterME[iPtBin]->GetSigma());
1057  background=hMassPtBin->Integral(minBinBkg,maxBinBkg);
1058  hSoverBME->SetBinContent(iPtBin+1,fitterME[iPtBin]->GetRawYield()/background);
1059  hSoverBME->SetBinError(iPtBin+1,0.000001);
1060  hRelStatME->SetBinContent(iPtBin+1,fitterME[iPtBin]->GetRawYieldError()/fitterME[iPtBin]->GetRawYield());
1061  hRelStatME->SetBinError(iPtBin+1,0.000001);
1062  hSignifME->SetBinContent(iPtBin+1,fitterME[iPtBin]->GetRawYield()/TMath::Sqrt(background+fitterME[iPtBin]->GetRawYield()));
1063  hSignifME->SetBinError(iPtBin+1,0.00000001);
1064  hInvMassHistoBinWidthME->SetBinContent(iPtBin+1,hMassSubME->GetBinWidth(1));
1065  hGausMeanME->SetBinContent(iPtBin+1,fitterME[iPtBin]->GetMean());
1066  hGausMeanME->SetBinError(iPtBin+1,fitterME[iPtBin]->GetMeanUncertainty());
1067  hGausSigmaME->SetBinContent(iPtBin+1,fitterME[iPtBin]->GetSigma());
1068  hGausSigmaME->SetBinError(iPtBin+1,fitterME[iPtBin]->GetSigmaUncertainty());
1069  hChiSqME->SetBinContent(iPtBin+1,fitterME[iPtBin]->GetReducedChiSquare());
1070  hChiSqME->SetBinError(iPtBin+1,0.00001); // very small number, for graphics
1071  hNdfME->SetBinContent(iPtBin+1,fitterME[iPtBin]->GetMassFunc()->GetNDF());
1072  hNdfME->SetBinError(iPtBin+1,0.00001); // very small number, for graphics
1073  // if(!correctForRefl)
1074  WriteFitInfo(fitterME[iPtBin],hMassPtBin);
1075  fout->cd();
1076  WriteFitFunctionsToFile(fitterME[iPtBin],"ME",iPtBin);
1077  current->cd();
1078 
1079  Double_t errbc;
1080  Double_t bc=fitterME[iPtBin]->GetRawYieldBinCounting(errbc,nsigmaBinCounting,optBkgBinCount);
1081  hRawYieldMEBC->SetBinContent(iPtBin+1,bc);
1082  hRawYieldMEBC->SetBinError(iPtBin+1,errbc);
1083 
1084 
1085  c4pulls->cd(iPtBin+1);
1086  TH1F *hPullsTrend=new TH1F();// the name is changed afterwards, histo must not be deleted
1087  TH1F *hPulls=new TH1F();//the name is changed afterwards, histo must not be deleted
1088  TH1F *hResidualTrend=new TH1F();// the name is changed afterwards, histo must not be deleted
1089  TH1F *hResiduals=fitterME[iPtBin]->GetResidualsAndPulls(hPulls,hResidualTrend,hPullsTrend);
1090  hResiduals->SetName(Form("hResidualsME_PtBin%d",iPtBin));
1091  hPulls->Draw();
1092  PrintGausParams(hPulls);
1093 
1094  c4residuals->cd(iPtBin+1);
1095  hResiduals->Draw();
1096 
1097  SetStyleHisto(hResidualTrend,kME);
1098  c4residualTrend->cd(iPtBin+1);
1099  hResidualTrend->Draw();
1100 
1101  cCompareResidualTrends->cd(iPtBin+1);
1102  if(out4||out1||out2)hResidualTrend->Draw("same");
1103  else hResidualTrend->Draw();
1104 
1105  c4pullTrend->cd(iPtBin+1);
1106  hPullsTrend->Draw();
1107  }
1108  else hMassSubME->Draw("");
1109  if(correctForRefl){
1110  delete hMCReflPtBin;
1111  delete hMCSigPtBin;
1112  }
1113 
1114  }
1115  TString path(gSystem->pwd());
1116  path.Append("/figures");
1117  if(gSystem->AccessPathName(path.Data())){
1118  gROOT->ProcessLine(Form(".!mkdir -p %s",path.Data()));
1119  }
1120 
1121  if(saveCanvasAsEps>0){
1122  c1->SaveAs(Form("figures/InvMassSpectra_%s_%s_NoBkgSub.eps",sigConf.Data(),suffix.Data()));
1123  c2->SaveAs(Form("figures/InvMassSpectra_%s_%s_Rot.eps",sigConf.Data(),suffix.Data()));
1124  c3->SaveAs(Form("figures/InvMassSpectra_%s_%s_LS.eps",sigConf.Data(),suffix.Data()));
1125  c4->SaveAs(Form("figures/InvMassSpectra_%s_%s_EM.eps",sigConf.Data(),suffix.Data()));
1126  if(tryDirectFit){
1127  c5->SaveAs(Form("figures/InvMassSpectra_%s_%s_SB.eps",sigConf.Data(),suffix.Data()));
1128  c5sub->SaveAs(Form("figures/InvMassSpectra_%s_%s_SBsub.eps",sigConf.Data(),suffix.Data()));
1129  }
1130 
1131  if(saveCanvasAsEps>1){
1132  c2residuals->SaveAs(Form("figures/ResidualDistribution_%s_%s_Rot.eps",sigConf.Data(),suffix.Data()));
1133  c3residuals->SaveAs(Form("figures/ResidualDistribution_%s_%s_LS.eps",sigConf.Data(),suffix.Data()));
1134  c4residuals->SaveAs(Form("figures/ResidualDistribution_%s_%s_EM.eps",sigConf.Data(),suffix.Data()));
1135 
1136  c2residualTrend->SaveAs(Form("figures/residualTrendvsMass_%s_%s_Rot.eps",sigConf.Data(),suffix.Data()));
1137  c3residualTrend->SaveAs(Form("figures/residualTrendvsMass_%s_%s_LS.eps",sigConf.Data(),suffix.Data()));
1138  c4residualTrend->SaveAs(Form("figures/residualTrendvsMass_%s_%s_EM.eps",sigConf.Data(),suffix.Data()));
1139 
1140  c2pulls->SaveAs(Form("figures/PullDistribution_%s_%s_Rot.eps",sigConf.Data(),suffix.Data()));
1141  c3pulls->SaveAs(Form("figures/PullDistribution_%s_%s_LS.eps",sigConf.Data(),suffix.Data()));
1142  c4pulls->SaveAs(Form("figures/PullDistribution_%s_%s_EM.eps",sigConf.Data(),suffix.Data()));
1143 
1144  c2pullTrend->SaveAs(Form("figures/pullTrendvsMass_%s_%s_Rot.eps",sigConf.Data(),suffix.Data()));
1145  c3pullTrend->SaveAs(Form("figures/pullTrendvsMass_%s_%s_LS.eps",sigConf.Data(),suffix.Data()));
1146  c4pullTrend->SaveAs(Form("figures/pullTrendvsMass_%s_%s_EM.eps",sigConf.Data(),suffix.Data()));
1147 
1148  if(tryDirectFit){
1149  c5residuals->SaveAs(Form("figures/ResidualDistribution_%s_%s_SB.eps",sigConf.Data(),suffix.Data()));
1150  c5residualTrend->SaveAs(Form("figures/residualTrendvsMass_%s_%s_SB.eps",sigConf.Data(),suffix.Data()));
1151  c5pulls->SaveAs(Form("figures/PullDistribution_%s_%s_SB.eps",sigConf.Data(),suffix.Data()));
1152  c5pullTrend->SaveAs(Form("figures/pullTrendvsMass_%s_%s_SB.eps",sigConf.Data(),suffix.Data()));
1153  }
1154  }
1155  }
1156 
1157  // save also .root
1158  if(saveCanvasAsRoot){
1159  c2->SaveAs(Form("figures/InvMassSpectra_%s_%s_Rot.root",sigConf.Data(),suffix.Data()));
1160  c3->SaveAs(Form("figures/InvMassSpectra_%s_%s_LS.root",sigConf.Data(),suffix.Data()));
1161  c4->SaveAs(Form("figures/InvMassSpectra_%s_%s_EM.root",sigConf.Data(),suffix.Data()));
1162 
1163  c2residuals->SaveAs(Form("figures/ResidualDistribution_%s_%s_Rot.root",sigConf.Data(),suffix.Data()));
1164  c3residuals->SaveAs(Form("figures/ResidualDistribution_%s_%s_LS.root",sigConf.Data(),suffix.Data()));
1165  c4residuals->SaveAs(Form("figures/ResidualDistribution_%s_%s_EM.root",sigConf.Data(),suffix.Data()));
1166 
1167  c2residualTrend->SaveAs(Form("figures/residualTrendvsMass_%s_%s_Rot.root",sigConf.Data(),suffix.Data()));
1168  c3residualTrend->SaveAs(Form("figures/residualTrendvsMass_%s_%s_LS.root",sigConf.Data(),suffix.Data()));
1169  c4residualTrend->SaveAs(Form("figures/residualTrendvsMass_%s_%s_EM.root",sigConf.Data(),suffix.Data()));
1170 
1171  c2pulls->SaveAs(Form("figures/PullDistribution_%s_%s_Rot.root",sigConf.Data(),suffix.Data()));
1172  c3pulls->SaveAs(Form("figures/PullDistribution_%s_%s_LS.root",sigConf.Data(),suffix.Data()));
1173  c4pulls->SaveAs(Form("figures/PullDistribution_%s_%s_EM.root",sigConf.Data(),suffix.Data()));
1174 
1175  c2pullTrend->SaveAs(Form("figures/pullTrendvsMass_%s_%s_Rot.root",sigConf.Data(),suffix.Data()));
1176  c3pullTrend->SaveAs(Form("figures/pullTrendvsMass_%s_%s_LS.root",sigConf.Data(),suffix.Data()));
1177  c4pullTrend->SaveAs(Form("figures/pullTrendvsMass_%s_%s_EM.root",sigConf.Data(),suffix.Data()));
1178 
1179  if(tryDirectFit){
1180  c5->SaveAs(Form("figures/InvMassSpectra_%s_%s_SB.root",sigConf.Data(),suffix.Data()));
1181  c5sub->SaveAs(Form("figures/InvMassSpectra_%s_%s_SBsub.root",sigConf.Data(),suffix.Data()));
1182  c5residuals->SaveAs(Form("figures/ResidualDistribution_%s_%s_SB.root",sigConf.Data(),suffix.Data()));
1183  c5residualTrend->SaveAs(Form("figures/residualTrendvsMass_%s_%s_SB.root",sigConf.Data(),suffix.Data()));
1184  c5pulls->SaveAs(Form("figures/PullDistribution_%s_%s_SB.root",sigConf.Data(),suffix.Data()));
1185  c5pullTrend->SaveAs(Form("figures/pullTrendvsMass_%s_%s_SB.root",sigConf.Data(),suffix.Data()));
1186  }
1187  }
1188 
1189 
1190  hRawYieldRot->SetMarkerStyle(21);
1191  hRawYieldLS->SetMarkerStyle(22);
1192  hRawYieldLS->SetMarkerColor(kGreen+2);
1193  hRawYieldLS->SetLineColor(kGreen+2);
1194  hRawYieldME->SetMarkerStyle(25);
1195  hRawYieldME->SetMarkerColor(4);
1196  hRawYieldME->SetLineColor(4);
1197  hRawYieldSB->SetMarkerStyle(27);
1198  hRawYieldSB->SetMarkerColor(6);
1199  hRawYieldSB->SetLineColor(6);
1200 
1201 
1202  TCanvas* cry=new TCanvas("cry","RawYield",800,700);
1203  cry->SetLeftMargin(0.15);
1204  hRawYieldRot->Draw("P");
1205  hRawYieldRot->SetMinimum(0);
1206  hRawYieldRot->GetYaxis()->SetTitleOffset(1.8);
1207  Double_t max=hRawYieldRot->GetMaximum();
1208  if(hRawYieldLS->GetMaximum()>max)max=hRawYieldLS->GetMaximum();
1209  if(hRawYieldME->GetMaximum()>max)max=hRawYieldME->GetMaximum();
1210  if(tryDirectFit){
1211  if(hRawYieldSB->GetMaximum()>max)max=hRawYieldSB->GetMaximum();
1212  }
1213  hRawYieldRot->SetMaximum(max*1.2);
1214  hRawYieldLS->Draw("PZSAME");
1215  hRawYieldME->Draw("PSAME");
1216  if(tryDirectFit) hRawYieldSB->Draw("PSAME");
1217  TLegend* legry=new TLegend(0.7,0.7,0.89,0.89);
1218  legry->SetFillStyle(0);
1219  legry->SetBorderSize(0);
1220  legry->AddEntry(hRawYieldRot,"Rotational","PL")->SetTextColor(1);
1221  legry->AddEntry(hRawYieldLS,"Like Sign","PL")->SetTextColor(kGreen+2);
1222  legry->AddEntry(hRawYieldME,"Ev Mix","PL")->SetTextColor(4);
1223  if(tryDirectFit){
1224  legry->AddEntry(hRawYieldSB,"Side-Band Fit","PL")->SetTextColor(6);
1225  }
1226  legry->Draw();
1227  if(saveCanvasAsEps>0) cry->SaveAs(Form("figures/RawYield_%s_%s.eps",sigConf.Data(),suffix.Data()));
1228  if(saveCanvasAsRoot) cry->SaveAs(Form("figures/RawYield_%s_%s.root",sigConf.Data(),suffix.Data()));
1229 
1230 
1231 
1232  TCanvas* cch2=new TCanvas("cch2","Chi2",800,700);
1233  cch2->SetLeftMargin(0.15);
1234  hChiSqRot->SetMarkerStyle(21);
1235  hChiSqRot->Draw("P");
1236  hChiSqRot->SetMinimum(0);
1237  hChiSqRot->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1238  hChiSqRot->GetYaxis()->SetTitle("#chi^{2}/ndf");
1239  hChiSqRot->GetYaxis()->SetTitleOffset(1.8);
1240  Double_t maxc=hChiSqRot->GetMaximum();
1241  if(hChiSqLS->GetMaximum()>maxc)maxc=hChiSqLS->GetMaximum();
1242  if(hChiSqME->GetMaximum()>maxc)maxc=hChiSqME->GetMaximum();
1243  if(tryDirectFit){
1244  if(hChiSqSB->GetMaximum()>maxc)maxc=hChiSqSB->GetMaximum();
1245  }
1246  hChiSqRot->SetMaximum(maxc*1.2);
1247 
1248  hChiSqLS->SetMarkerStyle(22);
1249  hChiSqLS->SetMarkerColor(kGreen+2);
1250  hChiSqLS->SetLineColor(kGreen+2);
1251  hChiSqLS->Draw("PZSAME");
1252  hChiSqME->SetMarkerStyle(25);
1253  hChiSqME->SetMarkerColor(4);
1254  hChiSqME->SetLineColor(4);
1255  hChiSqME->Draw("PSAME");
1256  if(tryDirectFit){
1257  hChiSqSB->SetMarkerStyle(27);
1258  hChiSqSB->SetMarkerColor(6);
1259  hChiSqSB->SetLineColor(6);
1260  hChiSqSB->Draw("PSAME");
1261  }
1262  legry->Draw();
1263  if(saveCanvasAsEps>0) cch2->SaveAs(Form("figures/ChiSq_%s_%s.eps",sigConf.Data(),suffix.Data()));
1264  if(saveCanvasAsRoot) cch2->SaveAs(Form("figures/ChiSq_%s_%s.root",sigConf.Data(),suffix.Data()));
1265 
1266  TH1F* hRatioLSToME=(TH1F*)hRawYieldLS->Clone("hRatioLStoME");
1267  TH1F* hRatioRotToME=(TH1F*)hRawYieldRot->Clone("hRatioRottoME");
1268  TH1F* hRatioMEToME=(TH1F*)hRawYieldME->Clone("hRatioMEtoME");
1269  TH1F* hRatioSBToME=(TH1F*)hRawYieldSB->Clone("hRatioSBtoME");
1270  for(Int_t ib=1; ib<=hRawYieldME->GetNbinsX(); ib++){
1271  Double_t yme=hRawYieldME->GetBinContent(ib);
1272  if(yme>0.){
1273  hRatioLSToME->SetBinContent(ib,hRawYieldLS->GetBinContent(ib)/yme);
1274  hRatioLSToME->SetBinError(ib,hRawYieldLS->GetBinError(ib)/yme);
1275  hRatioMEToME->SetBinContent(ib,hRawYieldME->GetBinContent(ib)/yme);
1276  hRatioMEToME->SetBinError(ib,hRawYieldME->GetBinError(ib)/yme);
1277  hRatioRotToME->SetBinContent(ib,hRawYieldRot->GetBinContent(ib)/yme);
1278  hRatioRotToME->SetBinError(ib,hRawYieldRot->GetBinError(ib)/yme);
1279  hRatioSBToME->SetBinContent(ib,hRawYieldSB->GetBinContent(ib)/yme);
1280  hRatioSBToME->SetBinError(ib,hRawYieldSB->GetBinError(ib)/yme);
1281  }
1282  }
1283 
1284  TCanvas* cry2=new TCanvas("cry2","RawYield+Ratios",1400,700);
1285  cry2->Divide(2,1);
1286  cry2->cd(1);
1287  gPad->SetLeftMargin(0.15);
1288  gPad->SetRightMargin(0.05);
1289  hRawYieldRot->Draw("P");
1290  hRawYieldLS->Draw("PZSAME");
1291  hRawYieldME->Draw("PSAME");
1292  if(tryDirectFit){
1293  hRawYieldSB->Draw("PSAME");
1294  }
1295  legry->Draw();
1296  cry2->cd(2);
1297  hRatioLSToME->SetStats(0);
1298  hRatioLSToME->SetMinimum(0.3);
1299  hRatioLSToME->SetMaximum(1.7);
1300  hRatioLSToME->GetYaxis()->SetTitle("Ratio To EvMix");
1301  hRatioLSToME->Draw("same");
1302  hRatioRotToME->Draw("same");
1303  hRatioMEToME->Draw("same");
1304  hRatioSBToME->Draw("same");
1305  if(saveCanvasAsEps>0) cry2->SaveAs(Form("figures/RawYieldAndRatios_%s_%s.eps",sigConf.Data(),suffix.Data()));
1306  if(saveCanvasAsRoot) cry2->SaveAs(Form("figures/RawYieldAndRatios_%s_%s.root",sigConf.Data(),suffix.Data()));
1307 
1308  hRelStatRot->GetYaxis()->SetTitleOffset(1.8);
1309  hRelStatRot->SetMarkerStyle(21);
1310  hRelStatLS->SetMarkerStyle(22);
1311  hRelStatLS->SetMarkerColor(kGreen+2);
1312  hRelStatLS->SetLineColor(kGreen+2);
1313  hRelStatME->SetMarkerStyle(25);
1314  hRelStatME->SetMarkerColor(4);
1315  hRelStatME->SetLineColor(4);
1316  hRelStatSB->SetMarkerStyle(27);
1317  hRelStatSB->SetMarkerColor(6);
1318  hRelStatSB->SetLineColor(6);
1319 
1320  TCanvas* cry3=new TCanvas("cry3","RawYield+Ratios+Unc",1800,600);
1321  cry3->Divide(3,1);
1322  cry3->cd(1);
1323  gPad->SetLeftMargin(0.15);
1324  gPad->SetRightMargin(0.05);
1325  hRawYieldRot->Draw("P");
1326  hRawYieldLS->Draw("PZSAME");
1327  hRawYieldME->Draw("PSAME");
1328  if(tryDirectFit){
1329  hRawYieldSB->Draw("PSAME");
1330  }
1331  legry->Draw();
1332  cry3->cd(2);
1333  hRatioLSToME->SetStats(0);
1334  hRatioLSToME->SetMinimum(0.3);
1335  hRatioLSToME->SetMaximum(1.7);
1336  hRatioLSToME->GetYaxis()->SetTitle("Ratio To EvMix");
1337  hRatioLSToME->Draw("same");
1338  hRatioRotToME->Draw("same");
1339  hRatioMEToME->Draw("same");
1340  hRatioSBToME->Draw("same");
1341  cry3->cd(3);
1342  gPad->SetLeftMargin(0.15);
1343  gPad->SetRightMargin(0.05);
1344  hRelStatRot->SetStats(0);
1345  hRelStatRot->SetMinimum(0.04);
1346  hRelStatRot->SetMaximum(0.2);
1347  hRelStatRot->Draw();
1348  hRelStatLS->Draw("same");
1349  hRelStatME->Draw("same");
1350  if(tryDirectFit) hRelStatSB->Draw("same");
1351 
1352  if(saveCanvasAsEps>0) cry3->SaveAs(Form("figures/RawYieldRatiosAndUnc_%s_%s.eps",sigConf.Data(),suffix.Data()));
1353  if(saveCanvasAsRoot) cry3->SaveAs(Form("figures/RawYieldRatiosAndUnc_%s_%s.root",sigConf.Data(),suffix.Data()));
1354 
1355 
1356 
1357 
1358  cCompareResidualTrends->cd(1);
1359  TLegend *legRT=new TLegend(*legry);
1360  legRT->SetX1NDC(0.12);
1361  legRT->SetY1NDC(0.7);
1362  legRT->SetX2NDC(0.3);
1363  legRT->SetY2NDC(0.9);
1364  legRT->Draw();
1365 
1366  // NOW COMPARE YIELDS OBTAINED WITH BC FROM DIFFERENT APPROACHES
1367  TCanvas* cryBC=new TCanvas("cryBC","RawYield with BC",800,700);
1368  cryBC->SetLeftMargin(0.15);
1369  hRawYieldRot->Draw("P");
1370  hRawYieldME->Draw("PSAME");
1371  hRawYieldLS->Draw("PZSAME");
1372  if(tryDirectFit){
1373  hRawYieldSB->Draw("PSAME");
1374  }
1375 
1376  if(hRawYieldRotBC->GetMaximum()>max)max=hRawYieldRotBC->GetMaximum();
1377  if(hRawYieldLSBC->GetMaximum()>max)max=hRawYieldLSBC->GetMaximum();
1378  if(hRawYieldMEBC->GetMaximum()>max)max=hRawYieldMEBC->GetMaximum();
1379  if(tryDirectFit){
1380  if(hRawYieldSBBC->GetMaximum()>max)max=hRawYieldSBBC->GetMaximum();
1381  }
1382  hRawYieldRot->SetMaximum(max*1.2);
1383 
1384  hRawYieldRotBC->SetMarkerStyle(21);
1385  hRawYieldRotBC->SetLineStyle(2);
1386  hRawYieldRotBC->Draw("PSAME");
1387 
1388 
1389  hRawYieldLSBC->SetMarkerStyle(22);
1390  hRawYieldLSBC->SetMarkerColor(kGreen+2);
1391  hRawYieldLSBC->SetLineColor(kGreen+2);
1392  hRawYieldLSBC->SetLineStyle(2);
1393  hRawYieldLSBC->Draw("PZSAME");
1394 
1395  hRawYieldMEBC->SetMarkerStyle(25);
1396  hRawYieldMEBC->SetMarkerColor(4);
1397  hRawYieldMEBC->SetLineColor(4);
1398  hRawYieldMEBC->SetLineStyle(2);
1399  hRawYieldMEBC->Draw("PSAME");
1400  if(tryDirectFit){
1401  hRawYieldSBBC->SetMarkerStyle(27);
1402  hRawYieldSBBC->SetMarkerColor(6);
1403  hRawYieldSBBC->SetLineColor(6);
1404  hRawYieldSBBC->SetLineStyle(2);
1405  hRawYieldSBBC->Draw("PSAME");
1406  }
1407 
1408  TLegend* legryBC=new TLegend(0.7,0.7,0.89,0.89);
1409  legryBC->SetFillStyle(0);
1410  legryBC->SetBorderSize(0);
1411  legryBC->AddEntry(hRawYieldRot,"Rotational","PL")->SetTextColor(1);
1412  legryBC->AddEntry(hRawYieldRotBC,"Rotational BC","PL")->SetTextColor(1);
1413  legryBC->AddEntry(hRawYieldLS,"Like Sign ","PL")->SetTextColor(kGreen+2);
1414  legryBC->AddEntry(hRawYieldLSBC,"Like Sign BC","PL")->SetTextColor(kGreen+2);
1415  legryBC->AddEntry(hRawYieldME,"Ev Mix","PL")->SetTextColor(4);
1416  legryBC->AddEntry(hRawYieldMEBC,"Ev Mix BC","PL")->SetTextColor(4);
1417  if(tryDirectFit){
1418  legryBC->AddEntry(hRawYieldSB,"Side-Band Fit","PL")->SetTextColor(6);
1419  legryBC->AddEntry(hRawYieldSBBC,"Side-Band Fit BC","PL")->SetTextColor(6);
1420  }
1421  legryBC->Draw();
1422  if(saveCanvasAsEps>0) cryBC->SaveAs(Form("figures/RawYieldBC_%s_%s.eps",sigConf.Data(),suffix.Data()));
1423  if(saveCanvasAsRoot) cryBC->SaveAs(Form("figures/RawYieldBC_%s_%s.root",sigConf.Data(),suffix.Data()));
1424 
1425 
1426 
1427 
1428  fout->cd();
1429  hRawYieldRot->Write();
1430  hRawYieldLS->Write();
1431  hRawYieldME->Write();
1432  hRawYieldSB->Write();
1433  hRelStatRot->Write();
1434  hRelStatLS->Write();
1435  hRelStatME->Write();
1436  hRelStatSB->Write();
1437  hSignifRot->Write();
1438  hSignifLS->Write();
1439  hSignifME->Write();
1440  hSignifSB->Write();
1441  hSoverBRot->Write();
1442  hSoverBLS->Write();
1443  hSoverBME->Write();
1444  hSoverBSB->Write();
1445  hGausMeanRot->Write();
1446  hGausMeanLS->Write();
1447  hGausMeanME->Write();
1448  hGausMeanSB->Write();
1449  hGausSigmaRot->Write();
1450  hGausSigmaLS->Write();
1451  hGausSigmaME->Write();
1452  hGausSigmaSB->Write();
1453  hChiSqRot->Write();
1454  hChiSqLS->Write();
1455  hChiSqME->Write();
1456  hChiSqSB->Write();
1457  hNdfRot->Write();
1458  hNdfLS->Write();
1459  hNdfME->Write();
1460  hNdfSB->Write();
1461  hInvMassHistoBinWidthRot->Write();
1462  hInvMassHistoBinWidthLS->Write();
1463  hInvMassHistoBinWidthME->Write();
1464  hInvMassHistoBinWidthSB->Write();
1465  hRebin->Write();
1466  hBkgFitFunc->Write();
1467  hBkgFitFuncSB->Write();
1468  hnEv->Write();
1469  if(hSigmaMC) hSigmaMC->Write();
1470  fout->Close();
1471 
1472  return;
1473 }
1474 
1475 void WriteFitFunctionsToFile(AliHFInvMassFitter *fitter, TString meth, Int_t iPtBin){
1476  TString nameh;
1477  TF1* fTot=fitter->GetMassFunc();
1478  nameh.Form("FitFuncTot_%s_PtBin%d",meth.Data(),iPtBin);
1479  fTot->SetRange(1.6,2.2);
1480  fTot->SetNpx(500);
1481  fTot->Write(nameh.Data());
1482  TF1* fSig=fitter->GetSignalFunc();
1483  nameh.Form("FitFuncSig_%s_PtBin%d",meth.Data(),iPtBin);
1484  fSig->SetRange(1.6,2.2);
1485  fSig->SetNpx(500);
1486  fSig->Write(nameh.Data());
1487  TF1* fBkg=fitter->GetBackgroundRecalcFunc();
1488  nameh.Form("FitFuncBkg_%s_PtBin%d",meth.Data(),iPtBin);
1489  fBkg->SetRange(1.6,2.2);
1490  fBkg->SetNpx(500);
1491  fBkg->Write(nameh.Data());
1492  if(meson=="Dzero"){
1493  TF1* fBkgR=fitter->GetBkgPlusReflFunc();
1494  nameh.Form("FitFuncBkgRefl_%s_PtBin%d",meth.Data(),iPtBin);
1495  fBkgR->SetRange(1.6,2.2);
1496  fBkgR->SetNpx(500);
1497  fBkgR->Write(nameh.Data());
1498  }
1499  return;
1500 }
1501 
1502 void WriteFitInfo(AliHFInvMassFitter *fitter, TH1D* histo){
1503  Double_t sig=fitter->GetRawYield();
1504  Double_t esig=fitter->GetRawYieldError();
1505  Double_t mean=fitter->GetMean();
1506  Double_t emean=fitter->GetMeanUncertainty();
1507  Double_t sigma=fitter->GetSigma();
1508  Double_t esigma=fitter->GetSigmaUncertainty();
1509  Double_t minBin=histo->FindBin(mean-3.*sigma);
1510  Double_t maxBin=histo->FindBin(mean+3.*sigma);
1511  Double_t back=histo->Integral(minBin,maxBin);
1512  TPaveText* tpar=new TPaveText(0.5,0.7,0.89,.87,"NDC");
1513  tpar->SetBorderSize(0);
1514  tpar->SetFillStyle(0);
1515  tpar->AddText(Form("Mean = %.3f #pm %.3f",mean,emean));
1516  tpar->AddText(Form("Sigma = %.3f #pm %.3f",sigma,esigma));
1517  tpar->SetTextColor(4);
1518  tpar->Draw();
1519 
1520  TPaveText* tss=new TPaveText(0.15,0.15,0.5,0.4,"NDC");
1521  tss->SetBorderSize(0);
1522  tss->SetFillStyle(0);
1523  tss->AddText(Form("S = %.0f #pm %.0f",sig,esig));
1524  tss->AddText(Form("B(3#sigma) = %.3g",back));
1525  tss->AddText(Form("S/B (3#sigma) = %.4f",sig/back));
1526  if(correctForRefl) tss->AddText(Form("Refl/Sig = %.3f #pm %.3f ",fitter->GetReflOverSig(),fitter->GetReflOverSigUncertainty()));
1527  tss->AddText(Form("Significance(3#sigma) = %.2f",sig/TMath::Sqrt(back+sig)));
1528  tss->SetTextColor(1);
1529  tss->Draw();
1530 
1531 }
1532 
1533 TH1F* FitMCInvMassSpectra(TList* lMC){
1534  TH3F* h3d=(TH3F*)lMC->FindObject("hMassVsPtVsYSig");
1535  if(!h3d){
1536  printf("hMassVsPtVsYSig not found\n");
1537  return 0x0;
1538  }
1539  TH3F* h3dr=(TH3F*)lMC->FindObject("hMassVsPtVsYRefl");
1540  if(!h3dr){
1541  printf("hMassVsPtVsYRefl not found\n");
1542  }
1543  TH1F* hSigmaMC=new TH1F("hSigmaMC","",nPtBins,binLims);
1544 
1545  TCanvas* cmc1=new TCanvas("InvMassMC","InvMassMC",1200,800);
1546  cmc1->Divide(4,2);
1547 
1548  gStyle->SetOptFit(0);
1549  gStyle->SetOptStat(0);
1550  for(Int_t iPtBin=0; iPtBin<nPtBins; iPtBin++){
1551  Int_t bin1=h3d->GetYaxis()->FindBin(binLims[iPtBin]);
1552  Int_t bin2=h3d->GetYaxis()->FindBin(binLims[iPtBin+1]-0.0001);
1553  Double_t ptmed=(binLims[iPtBin]+binLims[iPtBin+1])*0.5;
1554  Double_t pthalfwid=(binLims[iPtBin+1]-binLims[iPtBin])*0.5;
1555  printf("Bin %d Pt range=%f %f\n",iPtBin,h3d->GetYaxis()->GetBinLowEdge(bin1),h3d->GetYaxis()->GetBinUpEdge(bin2));
1556  TH1D* hMassMCPtBin=h3d->ProjectionX(Form("hMassMCPtBin%d",iPtBin),bin1,bin2);
1557  hMassMCPtBin->SetTitle(Form("%.1f<p_{T}<%.1f GeV/c",binLims[iPtBin],binLims[iPtBin+1]));
1558  hMassMCPtBin->GetXaxis()->SetTitle("Invariant mass (GeV/c^{2})");
1559  hMassMCPtBin->GetXaxis()->SetRangeUser(1.65,2.06);
1560  cmc1->cd(iPtBin+1);
1561  gPad->SetLogy();
1562  hMassMCPtBin->Fit("gaus","");
1563  TF1* fg=(TF1*)hMassMCPtBin->GetListOfFunctions()->FindObject("gaus");
1564  TLatex* tsig=new TLatex(0.65,0.82,"Signal");
1565  tsig->SetNDC();
1566  tsig->SetTextColor(kBlue+1);
1567  tsig->Draw();
1568 
1569  TPaveText* t1=new TPaveText(0.15,0.55,0.4,0.88,"ndc");
1570  t1->SetFillStyle(0);
1571  t1->SetBorderSize(0);
1572  t1->AddText(Form("#mu = %.4f GeV/c^{2}",fg->GetParameter(1)));
1573  t1->AddText(Form("#sigma = %.4f GeV/c^{2}",fg->GetParameter(2)));
1574  t1->AddText(Form("Integral = %.0f",fg->Integral(1.7,2.1)/hMassMCPtBin->GetBinWidth(1)));
1575  t1->AddText(Form("Entries = %.0f",hMassMCPtBin->GetEntries()));
1576 
1577  if(h3dr){
1578  TH1D* hReflPtBin=h3dr->ProjectionX(Form("hReflPtBin%d",iPtBin),bin1,bin2);
1579  hReflPtBin->SetTitle(Form("%.1f<p_{T}<%.1f GeV/c",binLims[iPtBin],binLims[iPtBin+1]));
1580  hReflPtBin->GetXaxis()->SetTitle("Invariant mass (GeV/c^{2})");
1581  hReflPtBin->SetLineColor(2);
1582  Int_t bin1=hReflPtBin->FindBin(fg->GetParameter(1)-3.*fg->GetParameter(2));
1583  Int_t bin2=hReflPtBin->FindBin(fg->GetParameter(1)+3.*fg->GetParameter(2));
1584  TLatex* tref=new TLatex(0.65,0.75,"Reflections");
1585  tref->SetNDC();
1586  tref->SetTextColor(2);
1587  tref->Draw();
1588  t1->AddText(Form("Refl(3#sigma) = %.0f",hReflPtBin->Integral(bin1,bin2)));
1589  hReflPtBin->Draw("same");
1590  t1->Draw();
1591  }
1592  sigmas[iPtBin]=fg->GetParameter(2);
1593  Double_t errMCsigma=fg->GetParError(2);
1594  if(tuneSigmaOnData>0.){
1595  sigmas[iPtBin]*=tuneSigmaOnData;
1596  errMCsigma*=tuneSigmaOnData;
1597  }
1598  hSigmaMC->SetBinContent(iPtBin+1,sigmas[iPtBin]);
1599  hSigmaMC->SetBinError(iPtBin+1,errMCsigma);
1600  }
1601  return hSigmaMC;
1602 }
rOverSmodif
Double_t rOverSmodif
Definition: ProjectCombinHFAndFit.C:70
AliHFInvMassFitter::GetMeanUncertainty
Double_t GetMeanUncertainty() const
Definition: AliHFInvMassFitter.h:106
AliHFInvMassFitter::GetReflOverSig
Double_t GetReflOverSig() const
Definition: AliHFInvMassFitter.h:109
AliHFInvMassFitter::GetResidualsAndPulls
TH1F * GetResidualsAndPulls(TH1 *hPulls=0x0, TH1 *hResidualTrend=0x0, TH1 *hPullsTrend=0x0, Double_t minrange=0, Double_t maxrange=-1, Int_t option=0)
Definition: AliHFInvMassFitter.cxx:1154
AliVertexingHFUtils::RebinHisto
static TH1D * RebinHisto(TH1 *hOrig, Int_t reb, Int_t firstUse=-1)
Rebinning of invariant mass histograms.
Definition: AliVertexingHFUtils.cxx:2577
kRot
Definition: ProjectCombinHFAndFit.C:22
ymax
const Double_t ymax
Definition: AddTaskCFDStar.C:7
Double_t
double Double_t
Definition: External.C:58
TH3F
Definition: External.C:260
AliHFInvMassFitter::SetFixReflOverS
void SetFixReflOverS(Double_t rovers)
Definition: AliHFInvMassFitter.h:93
AliHFInvMassFitter::DrawHere
void DrawHere(TVirtualPad *c, Double_t nsigma=3, Int_t writeFitInfo=1)
Definition: AliHFInvMassFitter.cxx:358
TH2F
Definition: External.C:236
saveCanvasAsEps
Int_t saveCanvasAsEps
Definition: ProjectCombinHFAndFit.C:45
reflopt
TString reflopt
Definition: ProjectCombinHFAndFit.C:68
AliHFInvMassFitter.h
DivideCanvas
void DivideCanvas(TCanvas *c, Int_t ndivisions)
Definition: ProjectCombinHFAndFit.C:210
fileName
TString fileName
Definition: ProjectCombinHFAndFit.C:25
AliVertexingHFUtils.h
gSystem
TSystem * gSystem
AliHFInvMassFitter::MassFitter
Int_t MassFitter(Bool_t draw=kTRUE)
Definition: AliHFInvMassFitter.cxx:234
Method
Method
Definition: ProjectCombinHFAndFit.C:22
AliHFInvMassFitter
AliHFInvMassFitter class for the fit of invariant mass distribution of charm hadrons.
Definition: AliHFInvMassFitter.h:24
PrintGausParams
void PrintGausParams(TH1F *hPulls)
Definition: ProjectCombinHFAndFit.C:129
nPtBins
const Int_t nPtBins
Definition: ProjectCombinHFAndFit.C:38
AliHFInvMassFitter::SetFitOption
void SetFitOption(TString opt)
Definition: AliHFInvMassFitter.h:49
c
TCanvas * c
Definition: TestFitELoss.C:172
nparback
Int_t nparback
Definition: ProjectCombinHFAndFit.C:75
optBkgBinCount
Int_t optBkgBinCount
Definition: ProjectCombinHFAndFit.C:62
AliVertexingHFUtils::AdaptTemplateRangeAndBinning
static TH1 * AdaptTemplateRangeAndBinning(const TH1 *hRef, TH1 *hData, Double_t minFit, Double_t maxFit)
Definition: AliVertexingHFUtils.cxx:2621
smoothLS
Int_t smoothLS
Definition: ProjectCombinHFAndFit.C:65
AliHFInvMassFitter::GetSigma
Double_t GetSigma() const
Definition: AliHFInvMassFitter.h:107
GetBackgroundNormalizationFactor
Double_t GetBackgroundNormalizationFactor(TH1D *hRatio, Int_t reb=1)
Definition: ProjectCombinHFAndFit.C:266
fileNameMC
TString fileNameMC
Definition: ProjectCombinHFAndFit.C:27
AliHFInvMassFitter::SetPolDegreeForBackgroundFit
void SetPolDegreeForBackgroundFit(Int_t deg)
Definition: AliHFInvMassFitter.h:52
fil
TFile * fil
Definition: InvMassFit.C:61
ProjectCombinHFAndFit
void ProjectCombinHFAndFit()
Definition: ProjectCombinHFAndFit.C:390
AliHFInvMassFitter::SetFixGaussianSigma
void SetFixGaussianSigma(Double_t sigma)
Definition: AliHFInvMassFitter.h:66
suffix
TString suffix
Definition: ProjectCombinHFAndFit.C:26
FitMCInvMassSpectra
TH1F * FitMCInvMassSpectra(TList *lMC)
Definition: ProjectCombinHFAndFit.C:1533
sigma
Double_t * sigma
Definition: CalculateAveragePt.C:58
binLims
Double_t binLims[nPtBins+1]
Definition: ProjectCombinHFAndFit.C:39
SetStyleHisto
void SetStyleHisto(TH1 *h, Int_t method, Int_t isXpt=-1)
Definition: ProjectCombinHFAndFit.C:168
massD
Double_t massD
Definition: ProjectCombinHFAndFit.C:63
tryDirectFit
Bool_t tryDirectFit
Definition: ProjectCombinHFAndFit.C:73
kSB
Definition: ProjectCombinHFAndFit.C:22
correctForRefl
Bool_t correctForRefl
Definition: ProjectCombinHFAndFit.C:69
AliHFInvMassFitter::SetFixGaussianMean
void SetFixGaussianMean(Double_t mean)
Definition: AliHFInvMassFitter.h:62
fixMeanConf
Int_t fixMeanConf
Definition: ProjectCombinHFAndFit.C:52
Int_t
int Int_t
Definition: External.C:63
AliHFInvMassFitter::GetMean
Double_t GetMean() const
Definition: AliHFInvMassFitter.h:105
fitrangelow
Double_t fitrangelow[nPtBins]
Definition: ProjectCombinHFAndFit.C:76
fixMean
Bool_t fixMean[nPtBins]
Definition: ProjectCombinHFAndFit.C:53
method
Int_t method
Definition: ComputeDmesonYield.C:32
suffixMC
TString suffixMC
Definition: ProjectCombinHFAndFit.C:29
sigmas
Double_t sigmas[nPtBins]
Definition: ProjectCombinHFAndFit.C:41
WriteFitFunctionsToFile
void WriteFitFunctionsToFile(AliHFInvMassFitter *fitter, TString meth, Int_t iPtBin)
Definition: ProjectCombinHFAndFit.C:1475
TH1D
Definition: External.C:212
QuadraticSmooth
Bool_t QuadraticSmooth(TH1 *h, Int_t ntimes=1)
Definition: ProjectCombinHFAndFit.C:141
nsigmaBinCounting
Double_t nsigmaBinCounting
Definition: ProjectCombinHFAndFit.C:61
hMCReflPtBin
TH1D * hMCReflPtBin
Definition: ProjectCombinHFAndFit.C:80
fitrangeup
Double_t fitrangeup[nPtBins]
Definition: ProjectCombinHFAndFit.C:77
AliHFInvMassFitter::GetRawYieldBinCounting
Double_t GetRawYieldBinCounting(Double_t &errRyBC, Double_t nSigma=3., Int_t option=0, Int_t pdgCode=0) const
Definition: AliHFInvMassFitter.cxx:1086
AliHFInvMassFitter::SetTemplateReflections
TH1F * SetTemplateReflections(const TH1 *h, TString opt, Double_t minRange, Double_t maxRange)
Definition: AliHFInvMassFitter.cxx:984
fitter
TList * fitter
Definition: DrawAnaELoss.C:26
fixSigma
Bool_t fixSigma[nPtBins]
Definition: ProjectCombinHFAndFit.C:50
kME
Definition: ProjectCombinHFAndFit.C:22
useEMwithLS
Bool_t useEMwithLS
Definition: ProjectCombinHFAndFit.C:59
AliHFInvMassFitter::GetReflOverSigUncertainty
Double_t GetReflOverSigUncertainty() const
Definition: AliHFInvMassFitter.h:113
saveCanvasAsRoot
Bool_t saveCanvasAsRoot
Definition: ProjectCombinHFAndFit.C:44
fitoption
TString fitoption
Definition: ProjectCombinHFAndFit.C:58
AliHFInvMassFitter::GetSigmaUncertainty
Double_t GetSigmaUncertainty() const
Definition: AliHFInvMassFitter.h:108
tuneSigmaOnData
Double_t tuneSigmaOnData
Definition: ProjectCombinHFAndFit.C:51
AliHFInvMassFitter::GetRawYieldError
Double_t GetRawYieldError() const
Definition: AliHFInvMassFitter.h:104
hMCSigPtBin
TH1D * hMCSigPtBin
Definition: ProjectCombinHFAndFit.C:81
nDegreeBackPol
Int_t nDegreeBackPol[nPtBins]
Definition: ProjectCombinHFAndFit.C:78
AliHFInvMassFitter::Background
void Background(Double_t nOfSigma, Double_t &background, Double_t &errbackground) const
Definition: AliHFInvMassFitter.cxx:837
AliNormalizationCounter.h
AliHFInvMassFitter::SetInitialGaussianSigma
void SetInitialGaussianSigma(Double_t sigma)
Definition: AliHFInvMassFitter.h:58
AliHFInvMassFitter::SetInitialGaussianMean
void SetInitialGaussianMean(Double_t mean)
Definition: AliHFInvMassFitter.h:57
fixSigmaConf
Int_t fixSigmaConf
Definition: ProjectCombinHFAndFit.C:49
minMass
Double_t minMass
Definition: ProjectCombinHFAndFit.C:54
optForNorm
Int_t optForNorm
Definition: ProjectCombinHFAndFit.C:56
ymin
const Double_t ymin
Definition: AddTaskCFDStar.C:6
rebin
Int_t rebin[nPtBins]
Definition: ProjectCombinHFAndFit.C:48
ConfigureFitter
AliHFInvMassFitter * ConfigureFitter(TH1D *histo, Int_t iPtBin, Int_t backcase, Double_t minFit, Double_t maxFit)
Definition: ProjectCombinHFAndFit.C:87
kLS
Definition: ProjectCombinHFAndFit.C:22
maxMass
Double_t maxMass
Definition: ProjectCombinHFAndFit.C:55
Bool_t
bool Bool_t
Definition: External.C:53
fout
TFile * fout
input train file
Definition: DrawAnaCaloTrackQA.C:90
GausPlusLine
TF1 * GausPlusLine(Double_t minRange=1.72, Double_t maxRange=2.05)
Definition: ProjectCombinHFAndFit.C:250
meson
TString meson
Definition: ProjectCombinHFAndFit.C:37
typeb
Int_t typeb
Definition: ProjectCombinHFAndFit.C:60
rangeForNorm
Double_t rangeForNorm
Definition: ProjectCombinHFAndFit.C:57
AliHFInvMassFitter::GetHistoClone
TH1F * GetHistoClone() const
Definition: AliHFInvMassFitter.h:134
TH1
Definition: External.C:196
AliHFInvMassFitter::GetRawYield
Double_t GetRawYield() const
Definition: AliHFInvMassFitter.h:103
WriteFitInfo
void WriteFitInfo(AliHFInvMassFitter *fitter, TH1D *histo)
Definition: ProjectCombinHFAndFit.C:1502