AliPhysics  5e2c166 (5e2c166)
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;
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(meson=="Dplus") massD=TDatabasePDG::Instance()->GetParticle(411)->Mass();
410  else if(meson=="Dzero") massD=TDatabasePDG::Instance()->GetParticle(421)->Mass();
411 
412  Int_t nBinsWithFixSig=0;
413  Int_t nBinsWithFixMean=0;
414  Int_t patSig=0;
415  Int_t patMean=0;
416  for(Int_t iPtBin=0; iPtBin<nPtBins; iPtBin++){
417  if(fixSigmaConf==1 || (fixSigmaConf==2 && fixSigma[iPtBin])){
418  nBinsWithFixSig++;
419  patSig+=1<<iPtBin;
420  }
421  if(fixMeanConf==1 || (fixMeanConf==2 && fixMean[iPtBin])){
422  nBinsWithFixMean++;
423  patMean+=1<<iPtBin;
424  }
425  }
426 
427  TFile* fil=new TFile(fileName.Data());
428  TDirectoryFile* df=(TDirectoryFile*)fil->Get(dirName.Data());
429 
430 
431  AliNormalizationCounter *nC=(AliNormalizationCounter*)df->Get("NormalizationCounter");
432  TH1F* hnEv=new TH1F("hEvForNorm","events for normalization",1,0,1);
433  printf("Number of Ev. for norm=%d\n",(Int_t)nC->GetNEventsForNorm());
434  hnEv->SetBinContent(1,nC->GetNEventsForNorm());
435 
436  TH1F* hRebin=new TH1F("hRebin","",nPtBins,binLims);
437  TH1F* hBkgFitFunc=new TH1F("hBkgFitFunc","",nPtBins,binLims);
438  TH1F* hBkgFitFuncSB=new TH1F("hBkgFitFuncSB","",nPtBins,binLims);
439 
440  TList* l=(TList*)df->Get(lstName.Data());
441  l->ls();
442 
443  TH3F* h3d=(TH3F*)l->FindObject("hMassVsPtVsY");
444  TH3F* h3dr=(TH3F*)l->FindObject("hMassVsPtVsYRot");
445  TH3F* h3dme=(TH3F*)l->FindObject("hMassVsPtVsYME");
446  TH3F* h3dmepp=(TH3F*)l->FindObject("hMassVsPtVsYMELSpp");
447  TH3F* h3dmemm=(TH3F*)l->FindObject("hMassVsPtVsYMELSmm");
448  TH2F* hpoolMix=(TH2F*)l->FindObject("hMixingsPerPool");
449  TH2F* hpoolEv=(TH2F*)l->FindObject("hEventsPerPool");
450  TH3F* h3dlsp=(TH3F*)l->FindObject("hMassVsPtVsYLSpp");
451  TH3F* h3dlsm=(TH3F*)l->FindObject("hMassVsPtVsYLSmm");
452 
453  TH3F* h3drefl=0x0;
454  TH3F* h3dmcsig=0x0;
455  TH1F* hSigmaMC=0x0;
456  TFile* filMC=new TFile(fileNameMC.Data());
457  if(filMC && filMC->IsOpen()){
458  TDirectoryFile* dfMC=(TDirectoryFile*)filMC->Get(dirNameMC.Data());
459  TList* lMC=(TList*)dfMC->Get(lstNameMC.Data());
460  hSigmaMC=FitMCInvMassSpectra(lMC);
461  if(nBinsWithFixSig>0 && !hSigmaMC){
462  printf("Fit to MC inv. mass spectra failed\n");
463  return;
464  }
465  if(correctForRefl){
466  h3drefl=(TH3F*)lMC->FindObject("hMassVsPtVsYRefl");
467  h3dmcsig=(TH3F*)lMC->FindObject("hMassVsPtVsYSig");
468  }
469  }
470 
471  TString sigConf="FixedSigma";
472  if(nBinsWithFixSig==0) sigConf="FreeSigma";
473  else if(nBinsWithFixSig==nPtBins) sigConf="FixedSigmaAll";
474  else sigConf=Form("FixedSigma%d",patSig);
475  if(nBinsWithFixSig>0 && tuneSigmaOnData>0.) sigConf+=Form("%d",TMath::Nint(tuneSigmaOnData*100.));
476  if(nBinsWithFixMean==nPtBins) sigConf+="FixedMeanAll";
477  else if(nBinsWithFixMean>0) sigConf+=Form("FixedMean%d",patMean);
478 
479  TCanvas* cem=new TCanvas("cem","Pools",1200,600);
480  cem->Divide(2,1);
481  cem->cd(1);
482  gPad->SetRightMargin(0.12);
483  gPad->SetLeftMargin(0.12);
484  hpoolMix->GetXaxis()->SetTitle("zVertex");
485  hpoolMix->GetYaxis()->SetTitle("Multiplicity");
486  hpoolMix->GetYaxis()->SetTitleOffset(1.4);
487  hpoolMix->Draw("colztext");
488  cem->cd(2);
489  gPad->SetRightMargin(0.12);
490  gPad->SetLeftMargin(0.12);
491  hpoolEv->Draw("colztext");
492  hpoolEv->GetXaxis()->SetTitle("zVertex");
493  hpoolEv->GetYaxis()->SetTitle("Multiplicity");
494  hpoolEv->GetYaxis()->SetTitleOffset(1.4);
495 
496 
497  TCanvas* c1=new TCanvas("c1","Mass",1200,800);
498  DivideCanvas(c1,nPtBins);
499 
500  TCanvas* c2=new TCanvas("c2","Mass-Bkg Rot",1200,800);
501  DivideCanvas(c2,nPtBins);
502  TCanvas* c2pulls=new TCanvas("c2pulls","Mass-Bkg Rot pulls",1200,800);
503  DivideCanvas(c2pulls,nPtBins);
504  TCanvas* c2residuals=new TCanvas("c2residuals","Mass-Bkg Rot residuals",1200,800);
505  DivideCanvas(c2residuals,nPtBins);
506  TCanvas* c2residualTrend=new TCanvas("c2residualTrend","Mass-Bkg Rot residual trend vs. mass",1200,800);
507  DivideCanvas(c2residualTrend,nPtBins);
508  TCanvas* c2pullTrend=new TCanvas("c2pullTrend","Mass-Bkg Rot pull trend vs. mass",1200,800);
509  DivideCanvas(c2pullTrend,nPtBins);
510 
511  TCanvas* c3=new TCanvas("c3","Mass-Bkg LS",1200,800);
512  DivideCanvas(c3,nPtBins);
513  TCanvas* c3pulls=new TCanvas("c3pulls","Mass-Bkg LS pulls",1200,800);
514  DivideCanvas(c3pulls,nPtBins);
515  TCanvas* c3residuals=new TCanvas("c3residuals","Mass-Bkg LS residuals",1200,800);
516  DivideCanvas(c3residuals,nPtBins);
517  TCanvas* c3residualTrend=new TCanvas("c3residualTrend","Mass-Bkg LS residual trend vs. mass",1200,800);
518  DivideCanvas(c3residualTrend,nPtBins);
519  TCanvas* c3pullTrend=new TCanvas("c3pullTrend","Mass-Bkg LS pull trend vs. mass",1200,800);
520  DivideCanvas(c3pullTrend,nPtBins);
521 
522  TCanvas* c4=new TCanvas("c4","Mass-Bkg EM",1200,800);
523  DivideCanvas(c4,nPtBins);
524  TCanvas* c4pulls=new TCanvas("c4pulls","Mass-Bkg EM pulls",1200,800);
525  DivideCanvas(c4pulls,nPtBins);
526  TCanvas* c4residuals=new TCanvas("c4residuals","Mass-Bkg EM residuals",1200,800);
527  DivideCanvas(c4residuals,nPtBins);
528  TCanvas* c4residualTrend=new TCanvas("c4residualTrend","Mass-Bkg EM residual trend vs. mass",1200,800);
529  DivideCanvas(c4residualTrend,nPtBins);
530  TCanvas* c4pullTrend=new TCanvas("c4pullTrend","Mass-Bkg EM pull trend vs. mass",1200,800);
531  DivideCanvas(c4pullTrend,nPtBins);
532 
533  TCanvas *c5=0x0;
534  TCanvas *c5sub=0x0;
535  TCanvas *c5pulls=0x0;
536  TCanvas *c5residuals=0x0;
537  TCanvas *c5residualTrend=0x0;
538  TCanvas *c5pullTrend=0x0;
539 
540  if(tryDirectFit){
541  c5=new TCanvas("c5","Mass SB Fit",1200,800);
542  DivideCanvas(c5,nPtBins);
543  c5sub=new TCanvas("c5sub","Mass-Bkg SB",1200,800);
544  DivideCanvas(c5sub,nPtBins);
545  c5pulls=new TCanvas("c5pulls","Mass-Bkg SB pulls",1200,800);
546  DivideCanvas(c5pulls,nPtBins);
547  c5residuals=new TCanvas("c5residuals","Mass-Bkg SB residuals",1200,800);
548  DivideCanvas(c5residuals,nPtBins);
549  c5residualTrend=new TCanvas("c5residualTrend","Mass-Bkg SB residual trend vs. mass",1200,800);
550  DivideCanvas(c5residualTrend,nPtBins);
551  c5pullTrend=new TCanvas("c5pullTrend","Mass-Bkg SB pull trend vs. mass",1200,800);
552  DivideCanvas(c5pullTrend,nPtBins);
553  }
554 
555  TCanvas *cCompareResidualTrends=new TCanvas("cCompareResidualTrends","cCompareResidualTrends",1200,800);
556  DivideCanvas(cCompareResidualTrends,nPtBins);
557 
558 
559  AliHFInvMassFitter* fitterRot[nPtBins];
560  AliHFInvMassFitter* fitterLS[nPtBins];
561  AliHFInvMassFitter* fitterME[nPtBins];
562  AliHFInvMassFitter* fitterSB[nPtBins];
563 
564  TH1F* hRawYieldRot=new TH1F("hRawYieldRot"," ; p_{T} (GeV/c) ; Raw yield",nPtBins,binLims);
565  TH1F* hRawYieldLS=new TH1F("hRawYieldLS"," ; p_{T} (GeV/c) ; Raw yield",nPtBins,binLims);
566  TH1F* hRawYieldME=new TH1F("hRawYieldME"," ; p_{T} (GeV/c) ; Raw yield",nPtBins,binLims);
567  TH1F* hRawYieldSB=new TH1F("hRawYieldSB"," ; p_{T} (GeV/c) ; Raw yield",nPtBins,binLims);
568 
569  TH1F* hRelStatRot=new TH1F("hRelStatRot"," ; p_{T} (GeV/c) ; Relative stat. unc.",nPtBins,binLims);
570  TH1F* hRelStatLS=new TH1F("hRelStatLS"," ; p_{T} (GeV/c) ; Relative stat. unc.",nPtBins,binLims);
571  TH1F* hRelStatME=new TH1F("hRelStatME"," ; p_{T} (GeV/c) ; Relative stat. unc.",nPtBins,binLims);
572  TH1F* hRelStatSB=new TH1F("hRelStatSB"," ; p_{T} (GeV/c) ; Relative stat. unc.",nPtBins,binLims);
573 
574  TH1F* hSignifRot=new TH1F("hSignifRot","",nPtBins,binLims);
575  TH1F* hSignifLS=new TH1F("hSignifLS","",nPtBins,binLims);
576  TH1F* hSignifME=new TH1F("hSignifME","",nPtBins,binLims);
577  TH1F* hSignifSB=new TH1F("hSignifSB","",nPtBins,binLims);
578 
579  TH1F* hSoverBRot=new TH1F("hSoverBRot","",nPtBins,binLims);
580  TH1F* hSoverBLS=new TH1F("hSoverBLS","",nPtBins,binLims);
581  TH1F* hSoverBME=new TH1F("hSoverBME","",nPtBins,binLims);
582  TH1F* hSoverBSB=new TH1F("hSoverBSB","",nPtBins,binLims);
583 
584  TH1F* hGausMeanRot=new TH1F("hGausMeanRot","",nPtBins,binLims);
585  TH1F* hGausMeanLS=new TH1F("hGausMeanLS","",nPtBins,binLims);
586  TH1F* hGausMeanME=new TH1F("hGausMeanME","",nPtBins,binLims);
587  TH1F* hGausMeanSB=new TH1F("hGausMeanSB","",nPtBins,binLims);
588 
589  TH1F* hGausSigmaRot=new TH1F("hGausSigmaRot","",nPtBins,binLims);
590  TH1F* hGausSigmaLS=new TH1F("hGausSigmaLS","",nPtBins,binLims);
591  TH1F* hGausSigmaME=new TH1F("hGausSigmaME","",nPtBins,binLims);
592  TH1F* hGausSigmaSB=new TH1F("hGausSigmaSB","",nPtBins,binLims);
593 
594  TH1F* hChiSqRot=new TH1F("hChiSqRot","",nPtBins,binLims);
595  TH1F* hChiSqLS=new TH1F("hChiSqLS","",nPtBins,binLims);
596  TH1F* hChiSqME=new TH1F("hChiSqME","",nPtBins,binLims);
597  TH1F* hChiSqSB=new TH1F("hChiSqSB","",nPtBins,binLims);
598 
599  TH1F* hNdfRot=new TH1F("hNdfRot","",nPtBins,binLims);
600  TH1F* hNdfLS=new TH1F("hNdfLS","",nPtBins,binLims);
601  TH1F* hNdfME=new TH1F("hNdfME","",nPtBins,binLims);
602  TH1F* hNdfSB=new TH1F("hNdfSB","",nPtBins,binLims);
603 
604  TH1F* hRawYieldRotBC=new TH1F("hRawYieldRotBC","BC yield (rotational background)",nPtBins,binLims);
605  TH1F* hRawYieldLSBC=new TH1F("hRawYieldLSBC","BC yield (like-sign background)",nPtBins,binLims);
606  TH1F* hRawYieldMEBC=new TH1F("hRawYieldMEBC","BC yield (mixed-event background)",nPtBins,binLims);
607  TH1F* hRawYieldSBBC=new TH1F("hRawYieldSBBC","BC yield (side-band fit background)",nPtBins,binLims);
608 
609  TH1F* hInvMassHistoBinWidthRot=new TH1F("hInvMassHistoBinWidthRot","",nPtBins,binLims);
610  TH1F* hInvMassHistoBinWidthLS=new TH1F("hInvMassHistoBinWidthLS","",nPtBins,binLims);
611  TH1F* hInvMassHistoBinWidthME=new TH1F("hInvMassHistoBinWidthME","",nPtBins,binLims);
612  TH1F* hInvMassHistoBinWidthSB=new TH1F("hInvMassHistoBinWidthSB","",nPtBins,binLims);
613 
614  TLatex* tME=new TLatex(0.65,0.82,"MixEv +- pairs");
615  tME->SetNDC();
616  TLatex* tMEpp=new TLatex(0.65,0.75,"MixEv ++ pairs");
617  tMEpp->SetNDC();
618  TLatex* tMEmm=new TLatex(0.65,0.68,"MixEv -- pairs");
619  tMEmm->SetNDC();
620 
621  TF1 *fpeak=new TF1("fpeak","[0]*1./(TMath::Sqrt(2.*TMath::Pi())*[2])*TMath::Exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",minMass,maxMass);
622 
623 
624  TDirectory *current = gDirectory;
625  TFile* fout=new TFile(Form("outputMassFits_%s_%s.root",sigConf.Data(),suffix.Data()),"recreate");
626  current->cd();
627 
628 
629  for(Int_t iPtBin=0; iPtBin<nPtBins; iPtBin++){
630 
631  Int_t bin1=h3d->GetYaxis()->FindBin(binLims[iPtBin]);
632  Int_t bin2=h3d->GetYaxis()->FindBin(binLims[iPtBin+1]-0.0001);
633  printf("Bin %d Pt range=%f %f\n",iPtBin,h3d->GetYaxis()->GetBinLowEdge(bin1),h3d->GetYaxis()->GetBinUpEdge(bin2));
634  TH1D* hMassPtBin=h3d->ProjectionX(Form("hMassPtBin%d",iPtBin),bin1,bin2);
635  TH1D* hMassPtBinr=h3dr->ProjectionX(Form("hMassPtBinr%d",iPtBin),bin1,bin2);
636  TH1D* hMassPtBinme=h3dme->ProjectionX(Form("hMassPtBinme%d",iPtBin),bin1,bin2);
637  TH1D* hMassPtBinmeLSpp=h3dmepp->ProjectionX(Form("hMassPtBinmeLSpp%d",iPtBin),bin1,bin2);
638  TH1D* hMassPtBinmeLSmm=h3dmemm->ProjectionX(Form("hMassPtBinmeLSmm%d",iPtBin),bin1,bin2);
639 
640  if(correctForRefl){
641  Int_t bin1MC=h3drefl->GetYaxis()->FindBin(binLims[iPtBin]);
642  Int_t bin2MC=h3drefl->GetYaxis()->FindBin(binLims[iPtBin+1]-0.0001);
643  hMCReflPtBin=h3drefl->ProjectionX(Form("hMCReflPtBin%d",iPtBin),bin1MC,bin2MC);
644  hMCSigPtBin=h3dmcsig->ProjectionX(Form("hMCSigPtBin%d",iPtBin),bin1MC,bin2MC);
645  }
646 
647  TH1D* hMassPtBinlsp=0x0;
648  TH1D* hMassPtBinlsm=0x0;
649  TH1D* hMassPtBinls=0x0;
650  if(h3dlsp){
651  hMassPtBinlsp=h3dlsp->ProjectionX(Form("hMassPtBinlsp%d",iPtBin),bin1,bin2);
652  hMassPtBinlsm=h3dlsm->ProjectionX(Form("hMassPtBinlsm%d",iPtBin),bin1,bin2);
653  hMassPtBinls=(TH1D*)hMassPtBinlsp->Clone(Form("hMassPtBinls%d",iPtBin));
654  hMassPtBinls->Reset();
655  for(Int_t iBin=1; iBin<=hMassPtBinlsp->GetNbinsX(); iBin++){
656  Double_t np=hMassPtBinlsp->GetBinContent(iBin);
657  Double_t nm=hMassPtBinlsm->GetBinContent(iBin);
658  Double_t tt=2*TMath::Sqrt(np*nm);
659  Double_t enp=hMassPtBinlsp->GetBinError(iBin);
660  Double_t enm=hMassPtBinlsm->GetBinError(iBin);
661  Double_t ett=0;
662  if(tt>0) ett=2./tt*TMath::Sqrt(np*np*enm*enm+nm*nm*enp*enp);
663  hMassPtBinls->SetBinContent(iBin,tt);
664  hMassPtBinls->SetBinError(iBin,ett);
665  }
666  if(smoothLS<-0.5)hMassPtBinls->Smooth(-1*smoothLS);
667  else if(smoothLS>0.5)QuadraticSmooth(hMassPtBinls,smoothLS);
668  hMassPtBinls->SetLineColor(kGreen+1);
669  }
670 
671  // hMassPtBin->Sumw2();
672  // hMassPtBinr->Sumw2();
673  // hMassPtBinme->Sumw2();
674  // hMassPtBinlsp->Sumw2();
675  // hMassPtBinls->Sumw2();
676  hMassPtBin->SetLineColor(1);
677  hMassPtBinr->SetLineColor(4);
678  hMassPtBinme->SetLineColor(kOrange+1);
679  hMassPtBinmeLSpp->SetLineColor(kGreen+2);
680  hMassPtBinmeLSmm->SetLineColor(kCyan);
681  hMassPtBin->SetTitle(Form("%.1f<p_{T}<%.1f GeV/c",binLims[iPtBin],binLims[iPtBin+1]));
682  hMassPtBinr->SetTitle(Form("%.1f<p_{T}<%.1f GeV/c",binLims[iPtBin],binLims[iPtBin+1]));
683  hMassPtBinme->SetTitle(Form("%.1f<p_{T}<%.1f GeV/c",binLims[iPtBin],binLims[iPtBin+1]));
684  TH1D* hRatioMEpp=(TH1D*)hMassPtBinmeLSpp->Clone(Form("hRatioPtBinmeLSpp%d",iPtBin));
685  hRatioMEpp->Divide(hMassPtBinme);
686  TH1D* hRatioMEmm=(TH1D*)hMassPtBinmeLSmm->Clone(Form("hRatioPtBinmeLSmm%d",iPtBin));
687  hRatioMEmm->Divide(hMassPtBinme);
688  TH1D* hMassPtBinmeAll=(TH1D*)hMassPtBinme->Clone(Form("hRatioPtBinmeAll%d",iPtBin));
689  hMassPtBinmeAll->Add(hMassPtBinmeLSpp);
690  hMassPtBinmeAll->Add(hMassPtBinmeLSmm);
691  hMassPtBinmeAll->SetLineColor(kRed);
692  TH1D* hRatioME=(TH1D*)hMassPtBinme->Clone(Form("hRatioME%d",iPtBin));
693  hRatioME->Divide(hMassPtBin);
694  TH1D* hRatioMEAll=(TH1D*)hMassPtBinmeAll->Clone(Form("hRatioMEAll%d",iPtBin));
695  hRatioMEAll->Divide(hMassPtBin);
696 
697 
698  TCanvas* c0=new TCanvas(Form("CBin%d",iPtBin),Form("Bin%d norm",iPtBin),1000,700);
699  c0->Divide(2,2);
700  c0->cd(1);
701  TH1D* hRatio=(TH1D*)hMassPtBinr->Clone(Form("hRatioFormNorm%d",iPtBin));
702  hRatio->Divide(hMassPtBin);
703  hRatio->Draw();
704  hRatio->GetYaxis()->SetTitle("Rotational/All");
705  hRatio->GetXaxis()->SetTitle("Invariant mass (GeV/c^{2})");
706  Double_t normRot=GetBackgroundNormalizationFactor(hRatio,rebin[iPtBin]);
707  hMassPtBinr->Scale(1./normRot);
708  TLatex* tnr=new TLatex(0.2,0.2,Form("Normaliz. factor = %f",normRot));
709  tnr->SetNDC();
710  tnr->Draw();
711  c0->cd(2);
712  hMassPtBinme->GetYaxis()->SetTitle("Entries (EvMix)");
713  hMassPtBinme->GetXaxis()->SetTitle("Invariant mass (GeV/c^{2})");
714  hMassPtBinme->DrawCopy();
715  hMassPtBinmeLSpp->Draw("same");
716  hMassPtBinmeLSmm->Draw("same");
717  tME->SetTextColor(hMassPtBinme->GetLineColor());
718  tMEpp->SetTextColor(hMassPtBinmeLSpp->GetLineColor());
719  tMEmm->SetTextColor(hMassPtBinmeLSmm->GetLineColor());
720  tME->Draw();
721  tMEpp->Draw();
722  tMEmm->Draw();
723  c0->cd(4);
724  hRatioMEpp->Draw();
725  hRatioMEpp->SetMinimum(0.4);
726  hRatioMEpp->SetMaximum(0.6);
727  hRatioMEpp->GetYaxis()->SetTitle("ME with LS / ME with OS");
728  hRatioMEpp->GetXaxis()->SetTitle("Invariant mass (GeV/c^{2})");
729  hRatioMEmm->Draw("same");
730  c0->cd(3);
731  hRatioME->Draw();
732  hRatioME->SetMaximum(hRatioMEAll->GetMaximum()*1.05);
733  hRatioMEAll->Draw("same");
734  hRatioME->GetYaxis()->SetTitle("EvMix/All");
735  hRatioME->GetXaxis()->SetTitle("Invariant mass (GeV/c^{2})");
736 
737  Double_t normME=GetBackgroundNormalizationFactor(hRatioME,rebin[iPtBin]);
738  Double_t normMEAll=GetBackgroundNormalizationFactor(hRatioMEAll,rebin[iPtBin]);
739  hMassPtBinme->Scale(1./normME);
740  hMassPtBinmeAll->Scale(1./normMEAll);
741  printf("Background norm bin %d DONE\n",iPtBin);
742 
743  c1->cd(iPtBin+1);
744  hMassPtBin->GetXaxis()->SetRangeUser(minMass,maxMass);
745  hMassPtBin->SetMinimum(0);
746  hMassPtBin->Draw();
747  hMassPtBin->GetYaxis()->SetTitle("Counts");
748  hMassPtBin->GetYaxis()->SetTitleOffset(2.);
749  hMassPtBin->GetXaxis()->SetTitle("Invariant mass (GeV/c^{2})");
750  hMassPtBinr->Draw("same");
751  if(!useEMwithLS) hMassPtBinme->Draw("same");
752  else hMassPtBinmeAll->Draw("same");
753  if(hMassPtBinls) hMassPtBinls->Draw("same");
754  if(iPtBin==0){
755  TLegend* leg=new TLegend(0.5,0.6,0.89,0.89);
756  leg->SetFillStyle(0);
757  leg->AddEntry(hMassPtBin,"All candidates","L")->SetTextColor(hMassPtBin->GetLineColor());
758  leg->AddEntry(hMassPtBinr,"Background (rotations)","L")->SetTextColor(hMassPtBinr->GetLineColor());
759  if(!useEMwithLS) leg->AddEntry(hMassPtBinme,"Background (ME)","L")->SetTextColor(hMassPtBinme->GetLineColor());
760  else leg->AddEntry(hMassPtBinmeAll,"Background (ME)","L")->SetTextColor(hMassPtBinmeAll->GetLineColor());
761  if(hMassPtBinls) leg->AddEntry(hMassPtBinls,"Like-sign","L")->SetTextColor(hMassPtBinls->GetLineColor());
762  leg->Draw();
763  }
764  gPad->Update();
765 
766  TH1D* hMassSubRot=(TH1D*)hMassPtBin->Clone(Form("hMassSubRot_bin%d",iPtBin));
767  hMassSubRot->Add(hMassPtBinr,-1);
768  hMassSubRot->SetTitle(Form("%.1f<p_{T}<%.1f GeV/c -- Rotational",binLims[iPtBin],binLims[iPtBin+1]));
769  TH1D* hMassSubME=(TH1D*)hMassPtBin->Clone(Form("hMassSubME_bin%d",iPtBin));
770  if(!useEMwithLS) hMassSubME->Add(hMassPtBinme,-1);
771  else hMassSubME->Add(hMassPtBinmeAll,-1);
772  hMassSubME->SetTitle(Form("%.1f<p_{T}<%.1f GeV/c -- Mixed Ev",binLims[iPtBin],binLims[iPtBin+1]));
773  TH1D* hMassSubLS=0x0;
774  if(hMassPtBinls){
775  hMassSubLS=(TH1D*)hMassPtBin->Clone(Form("hMassSubLS_bin%d",iPtBin));
776  hMassSubLS->Add(hMassPtBinls,-1);
777  hMassSubLS->SetTitle(Form("%.1f<p_{T}<%.1f GeV/c -- Like Sign",binLims[iPtBin],binLims[iPtBin+1]));
778  }
779 
780  fout->cd();
781  hMassPtBin->Write();
782  hMassSubRot->Write();
783  hMassSubME->Write();
784  if(hMassPtBinls) hMassSubLS->Write();
785  if(hMCReflPtBin) hMCReflPtBin->Write();
786  if(hMCSigPtBin) hMCSigPtBin->Write();
787  current->cd();
788 
789  hMassSubRot=AliVertexingHFUtils::RebinHisto(hMassSubRot,rebin[iPtBin]);
790  hMassSubME=AliVertexingHFUtils::RebinHisto(hMassSubME,rebin[iPtBin]);
791  if(hMassPtBinls) hMassSubLS=AliVertexingHFUtils::RebinHisto(hMassSubLS,rebin[iPtBin]);
792 
793  hRebin->SetBinContent(iPtBin+1,rebin[iPtBin]);
794  Int_t bkgToFill=typeb;
795  if(typeb==6) bkgToFill=typeb+nDegreeBackPol[iPtBin];
796  Int_t bkgToFillSB=6+nDegreeBackPol[iPtBin];
797 
798  hBkgFitFunc->SetBinContent(iPtBin+1,bkgToFill);
799  hBkgFitFuncSB->SetBinContent(iPtBin+1,bkgToFillSB);
800 
801  fitterRot[iPtBin]=ConfigureFitter(hMassSubRot,iPtBin,typeb,minMass,maxMass);
802  if(hMassPtBinls) fitterLS[iPtBin]=ConfigureFitter(hMassSubLS,iPtBin,typeb,minMass,maxMass);
803  fitterME[iPtBin]=ConfigureFitter(hMassSubME,iPtBin,typeb,minMass,maxMass);
804 
805  Bool_t out1=fitterRot[iPtBin]->MassFitter(0);
806  Bool_t out2=kFALSE;
807  if(hMassPtBinls) out2=fitterLS[iPtBin]->MassFitter(0);
808  Bool_t out3=fitterME[iPtBin]->MassFitter(0);
809 
810  Double_t background=999999999.;
811  Bool_t out4=kFALSE;
812  if(tryDirectFit){
813  TH1D *hMassDirectFit=(TH1D*)hMassPtBin->Clone(Form("hMassDirectFit_bin%d",iPtBin));
814  hMassDirectFit=AliVertexingHFUtils::RebinHisto(hMassDirectFit,rebin[iPtBin]);
815  fitterSB[iPtBin]=ConfigureFitter(hMassDirectFit,iPtBin,6,fitrangelow[iPtBin],fitrangeup[iPtBin]);
816  out4=fitterSB[iPtBin]->MassFitter(0);//DirectFit(hMassDirectFit,iPtBin,hRawYieldSB);
817 
818  Double_t background,ebkg;
819 
820  if(out4 && fitterSB[iPtBin]->GetMassFunc()){
821  c5->cd(iPtBin+1);
822  fitterSB[iPtBin]->DrawHere(gPad,3,0);
823  hRawYieldSB->SetBinContent(iPtBin+1,fitterSB[iPtBin]->GetRawYield());
824  hRawYieldSB->SetBinError(iPtBin+1,fitterSB[iPtBin]->GetRawYieldError());
825  if(fitterSB[iPtBin]->GetRawYield()>0){
826  hRelStatSB->SetBinContent(iPtBin+1,fitterSB[iPtBin]->GetRawYieldError()/fitterSB[iPtBin]->GetRawYield());
827  hRelStatSB->SetBinError(iPtBin+1,0.00000001);
828  }
829  fitterSB[iPtBin]->Background(3.,background,ebkg);
830  hSignifSB->SetBinContent(iPtBin+1,fitterSB[iPtBin]->GetRawYield()/TMath::Sqrt(background+fitterSB[iPtBin]->GetRawYield()));
831  hSignifSB->SetBinError(iPtBin+1,0.00000001);
832  hSoverBSB->SetBinContent(iPtBin+1,fitterSB[iPtBin]->GetRawYield()/background);
833  hSoverBSB->SetBinError(iPtBin+1,0.00000001);
834  hInvMassHistoBinWidthSB->SetBinContent(iPtBin+1,hMassDirectFit->GetBinWidth(1));
835  hGausMeanSB->SetBinContent(iPtBin+1,fitterSB[iPtBin]->GetMean());
836  hGausMeanSB->SetBinError(iPtBin+1,fitterSB[iPtBin]->GetMeanUncertainty());
837  hGausSigmaSB->SetBinContent(iPtBin+1,fitterSB[iPtBin]->GetSigma());
838  hGausSigmaSB->SetBinError(iPtBin+1,fitterSB[iPtBin]->GetSigmaUncertainty());
839  hChiSqSB->SetBinContent(iPtBin+1,fitterSB[iPtBin]->GetReducedChiSquare());
840  hChiSqSB->SetBinError(iPtBin+1,0.00001); // very small number, for graphics
841  hNdfSB->SetBinContent(iPtBin+1,fitterSB[iPtBin]->GetMassFunc()->GetNDF());
842  hNdfSB->SetBinError(iPtBin+1,0.00001); // very small number, for graphics
843  // if(!correctForRefl)
844  WriteFitInfo(fitterSB[iPtBin],hMassPtBin);
845  fout->cd();
846  WriteFitFunctionsToFile(fitterSB[iPtBin],"SB",iPtBin);
847  current->cd();
848 
849  c5sub->cd(iPtBin+1);
850  TH1F* hsubTemp=(TH1F*)hMassDirectFit->Clone(Form("%sSubBack%d",hMassDirectFit->GetName(),iPtBin));
851  TH1F* hsubTempAllRange=(TH1F*)hMassDirectFit->Clone(Form("%sSubBackAllRange%d",hMassDirectFit->GetName(),iPtBin));
852  TF1* funcAll=fitterSB[iPtBin]->GetMassFunc();
853  TF1* funcBkg=fitterSB[iPtBin]->GetBackgroundRecalcFunc();
854  // TF1* funcBkg2=fitterSB[iPtBin]->GetBackgroundFullRangeFunc();
855  for(Int_t jst=1;jst<=hsubTemp->GetNbinsX();jst++){
856  Double_t backg=funcBkg->Integral(hsubTemp->GetBinLowEdge(jst),hsubTemp->GetBinLowEdge(jst)+hsubTemp->GetBinWidth(jst))/hsubTemp->GetBinWidth(jst);
857  Double_t tot=funcAll->Integral(hsubTempAllRange->GetBinLowEdge(jst),hsubTempAllRange->GetBinLowEdge(jst)+hsubTempAllRange->GetBinWidth(jst))/hsubTempAllRange->GetBinWidth(jst);
858  hsubTemp->SetBinContent(jst,hsubTemp->GetBinContent(jst)-backg);
859  hsubTempAllRange->SetBinContent(jst,hsubTempAllRange->GetBinContent(jst)-tot);
860  }
861  hsubTemp->SetLineColor(kBlue);
862  hsubTempAllRange->SetLineColor(kGray+2);
863  Double_t ymin=0;
864  Double_t ymax=1;
865  for(Int_t ibs=1; ibs<hsubTemp->GetNbinsX(); ibs++){
866  Double_t binc=hsubTemp->GetBinCenter(ibs);
867  if(binc>fitrangelow[iPtBin] && binc<fitrangeup[iPtBin]){
868  Double_t yl=hsubTemp->GetBinContent(ibs)-hsubTemp->GetBinError(ibs);
869  Double_t yu=hsubTemp->GetBinContent(ibs)+hsubTemp->GetBinError(ibs);
870  if(yl<ymin) ymin=yl;
871  if(yu>ymax) ymax=yu;
872  }
873  }
874  if(ymax>0) ymax*=1.2;
875  else ymax*=0.8;
876  if(ymin<0) ymin*=1.2;
877  else ymin*=0.8;
878  hsubTemp->GetXaxis()->SetRangeUser(fitrangelow[iPtBin],fitrangeup[iPtBin]);
879  hsubTemp->SetMinimum(ymin);
880  hsubTemp->SetMaximum(ymax);
881  hsubTemp->SetMarkerStyle(20);
882  hsubTemp->SetMarkerColor(hsubTemp->GetLineColor());
883  hsubTemp->DrawCopy();
884  hsubTempAllRange->DrawCopy("same");
885  hsubTemp->DrawCopy("same");
886  fpeak->SetRange(fitrangelow[iPtBin],fitrangeup[iPtBin]);
887  fpeak->SetParameter(0,funcAll->GetParameter(nDegreeBackPol[iPtBin]+1));
888  fpeak->SetParameter(1,funcAll->GetParameter(nDegreeBackPol[iPtBin]+2));
889  fpeak->SetParameter(2,funcAll->GetParameter(nDegreeBackPol[iPtBin]+3));
890  fpeak->DrawCopy("same");
891 
892  Double_t errbc;
893  Double_t bc=fitterSB[iPtBin]->GetRawYieldBinCounting(errbc,nsigmaBinCounting,optBkgBinCount);
894  hRawYieldSBBC->SetBinContent(iPtBin+1,bc);
895  hRawYieldSBBC->SetBinError(iPtBin+1,errbc);
896 
897  c5pulls->cd(iPtBin+1);
898  TH1F *hPulls=new TH1F();//the name is changed afterwards, histo must not be deleted
899  TH1F *hPullsTrend=new TH1F();// the name is changed afterwards, histo must not be deleted
900  TH1F *hResidualTrend=new TH1F();// the name is changed afterwards, histo must not be deleted
901  TH1F *hResiduals=fitterSB[iPtBin]->GetResidualsAndPulls(hPulls,hResidualTrend,hPullsTrend);
902  hResiduals->SetName(Form("hResidualsSB_PtBin%d",iPtBin));
903 
904  hPulls->Draw();
905  PrintGausParams(hPulls);
906 
907  c5residuals->cd(iPtBin+1);
908  hResiduals->Draw();
909 
910  SetStyleHisto(hResidualTrend,kSB);
911  cCompareResidualTrends->cd(iPtBin+1);
912  hResidualTrend->Draw();
913 
914  c5residualTrend->cd(iPtBin+1);
915  hResidualTrend->Draw();
916 
917  c5pullTrend->cd(iPtBin+1);
918  hPullsTrend->Draw();
919  //TPaveStats *tp=(TPaveStats*)hPulls->FindObject("stats");
920  // tp->SetOptStat("remnpcev");
921  delete hMassDirectFit;
922  }
923  }
924 
925  c2->cd(iPtBin+1);
926  if(out1 && fitterRot[iPtBin]->GetMassFunc()){
927  fitterRot[iPtBin]->DrawHere(gPad,3,0);
928  gPad->Update();
929  hRawYieldRot->SetBinContent(iPtBin+1,fitterRot[iPtBin]->GetRawYield());
930  hRawYieldRot->SetBinError(iPtBin+1,fitterRot[iPtBin]->GetRawYieldError());
931  Double_t minBinBkg=hMassPtBin->FindBin(fitterRot[iPtBin]->GetMean()-3.*fitterRot[iPtBin]->GetSigma());
932  Double_t maxBinBkg=hMassPtBin->FindBin(fitterRot[iPtBin]->GetMean()+3.*fitterRot[iPtBin]->GetSigma());
933  background=hMassPtBin->Integral(minBinBkg,maxBinBkg);
934  hSoverBRot->SetBinContent(iPtBin+1,fitterRot[iPtBin]->GetRawYield()/background);
935  hSoverBRot->SetBinError(iPtBin+1,0.000001);
936  hRelStatRot->SetBinContent(iPtBin+1,fitterRot[iPtBin]->GetRawYieldError()/fitterRot[iPtBin]->GetRawYield());
937  hRelStatRot->SetBinError(iPtBin+1,0.000001);
938  hSignifRot->SetBinContent(iPtBin+1,fitterRot[iPtBin]->GetRawYield()/TMath::Sqrt(background+fitterRot[iPtBin]->GetRawYield()));
939  hSignifRot->SetBinError(iPtBin+1,0.00000001);
940  hInvMassHistoBinWidthRot->SetBinContent(iPtBin+1,hMassSubRot->GetBinWidth(1));
941  hGausMeanRot->SetBinContent(iPtBin+1,fitterRot[iPtBin]->GetMean());
942  hGausMeanRot->SetBinError(iPtBin+1,fitterRot[iPtBin]->GetMeanUncertainty());
943  hGausSigmaRot->SetBinContent(iPtBin+1,fitterRot[iPtBin]->GetSigma());
944  hGausSigmaRot->SetBinError(iPtBin+1,fitterRot[iPtBin]->GetSigmaUncertainty());
945  hChiSqRot->SetBinContent(iPtBin+1,fitterRot[iPtBin]->GetReducedChiSquare());
946  hChiSqRot->SetBinError(iPtBin+1,0.00001); // very small number, for graphics
947  hNdfRot->SetBinContent(iPtBin+1,fitterRot[iPtBin]->GetMassFunc()->GetNDF());
948  hNdfRot->SetBinError(iPtBin+1,0.00001); // very small number, for graphics
949  // if(!correctForRefl)
950  WriteFitInfo(fitterRot[iPtBin],hMassPtBin);
951  fout->cd();
952  WriteFitFunctionsToFile(fitterRot[iPtBin],"Rot",iPtBin);
953  current->cd();
954 
955  Double_t errbc;
956  Double_t bc=fitterRot[iPtBin]->GetRawYieldBinCounting(errbc,nsigmaBinCounting,optBkgBinCount);
957  hRawYieldRotBC->SetBinContent(iPtBin+1,bc);
958  hRawYieldRotBC->SetBinError(iPtBin+1,errbc);
959 
960  c2pulls->cd(iPtBin+1);
961  TH1F *hPullsTrend=new TH1F();// the name is changed afterwards, histo must not be deleted
962  TH1F *hPulls=new TH1F();//the name is changed afterwards, histo must not be deleted
963  TH1F *hResidualTrend=new TH1F();// the name is changed afterwards, histo must not be deleted
964  TH1F *hResiduals=fitterRot[iPtBin]->GetResidualsAndPulls(hPulls,hResidualTrend,hPullsTrend);
965  hResiduals->SetName(Form("hResidualsRot_PtBin%d",iPtBin));
966 
967  hPulls->Draw();
968  PrintGausParams(hPulls);
969 
970  c2residuals->cd(iPtBin+1);
971  hResiduals->Draw();
972 
973  SetStyleHisto(hResidualTrend,kRot);
974  c2residualTrend->cd(iPtBin+1);
975  hResidualTrend->Draw();
976  cCompareResidualTrends->cd(iPtBin+1);
977  if(out4)hResidualTrend->Draw("same");
978  else hResidualTrend->Draw();
979 
980  c2pullTrend->cd(iPtBin+1);
981  hPullsTrend->Draw();
982 
983  }
984  else hMassSubRot->Draw("");
985 
986 
987  c3->cd(iPtBin+1);
988  if(out2 && fitterLS[iPtBin]->GetMassFunc()){
989  fitterLS[iPtBin]->DrawHere(gPad,3,0);
990  hRawYieldLS->SetBinContent(iPtBin+1,fitterLS[iPtBin]->GetRawYield());
991  hRawYieldLS->SetBinError(iPtBin+1,fitterLS[iPtBin]->GetRawYieldError());
992  Double_t minBinBkg=hMassPtBin->FindBin(fitterLS[iPtBin]->GetMean()-3.*fitterLS[iPtBin]->GetSigma());
993  Double_t maxBinBkg=hMassPtBin->FindBin(fitterLS[iPtBin]->GetMean()+3.*fitterLS[iPtBin]->GetSigma());
994  background=hMassPtBin->Integral(minBinBkg,maxBinBkg);
995  hSoverBLS->SetBinContent(iPtBin+1,fitterLS[iPtBin]->GetRawYield()/background);
996  hSoverBLS->SetBinError(iPtBin+1,0.000001);
997  hRelStatLS->SetBinContent(iPtBin+1,fitterLS[iPtBin]->GetRawYieldError()/fitterLS[iPtBin]->GetRawYield());
998  hRelStatLS->SetBinError(iPtBin+1,0.0000001);
999  hSignifLS->SetBinContent(iPtBin+1,fitterLS[iPtBin]->GetRawYield()/TMath::Sqrt(background+fitterLS[iPtBin]->GetRawYield()));
1000  hSignifLS->SetBinError(iPtBin+1,0.00000001);
1001  hInvMassHistoBinWidthLS->SetBinContent(iPtBin+1,hMassSubLS->GetBinWidth(1));
1002  hGausMeanLS->SetBinContent(iPtBin+1,fitterLS[iPtBin]->GetMean());
1003  hGausMeanLS->SetBinError(iPtBin+1,fitterLS[iPtBin]->GetMeanUncertainty());
1004  hGausSigmaLS->SetBinContent(iPtBin+1,fitterLS[iPtBin]->GetSigma());
1005  hGausSigmaLS->SetBinError(iPtBin+1,fitterLS[iPtBin]->GetSigmaUncertainty());
1006  hChiSqLS->SetBinContent(iPtBin+1,fitterLS[iPtBin]->GetReducedChiSquare());
1007  hChiSqLS->SetBinError(iPtBin+1,0.00001); // very small number, for graphics
1008  hNdfLS->SetBinContent(iPtBin+1,fitterLS[iPtBin]->GetMassFunc()->GetNDF());
1009  hNdfLS->SetBinError(iPtBin+1,0.00001); // very small number, for graphics
1010  // if(!correctForRefl)
1011  WriteFitInfo(fitterLS[iPtBin],hMassPtBin);
1012  fout->cd();
1013  WriteFitFunctionsToFile(fitterLS[iPtBin],"LS",iPtBin);
1014  current->cd();
1015 
1016 
1017  Double_t errbc;
1018  Double_t bc=fitterLS[iPtBin]->GetRawYieldBinCounting(errbc,nsigmaBinCounting,optBkgBinCount);
1019  hRawYieldLSBC->SetBinContent(iPtBin+1,bc);
1020  hRawYieldLSBC->SetBinError(iPtBin+1,errbc);
1021 
1022  c3pulls->cd(iPtBin+1);
1023  TH1F *hPullsTrend=new TH1F();// the name is changed afterwards, histo must not be deleted
1024  TH1F *hPulls=new TH1F();//the name is changed afterwards, histo must not be deleted
1025  TH1F *hResidualTrend=new TH1F();// the name is changed afterwards, histo must not be deleted
1026  TH1F *hResiduals=fitterLS[iPtBin]->GetResidualsAndPulls(hPulls,hResidualTrend,hPullsTrend);
1027  hResiduals->SetName(Form("hResidualsLS_PtBin%d",iPtBin));
1028 
1029  hPulls->Draw();
1030  PrintGausParams(hPulls);
1031 
1032  c3residuals->cd(iPtBin+1);
1033  hResiduals->Draw();
1034 
1035  SetStyleHisto(hResidualTrend,kLS);
1036  c3residualTrend->cd(iPtBin+1);
1037  hResidualTrend->Draw();
1038 
1039  cCompareResidualTrends->cd(iPtBin+1);
1040  if(out4||out1)hResidualTrend->Draw("same");
1041  else hResidualTrend->Draw();
1042 
1043  c3pullTrend->cd(iPtBin+1);
1044  hPullsTrend->Draw();
1045  }
1046  else if(hMassPtBinls) hMassSubLS->Draw("");
1047 
1048  c4->cd(iPtBin+1);
1049  if(out3 && fitterME[iPtBin]->GetMassFunc()){
1050  fitterME[iPtBin]->DrawHere(gPad,3,0);
1051  hRawYieldME->SetBinContent(iPtBin+1,fitterME[iPtBin]->GetRawYield());
1052  hRawYieldME->SetBinError(iPtBin+1,fitterME[iPtBin]->GetRawYieldError());
1053  Double_t minBinBkg=hMassPtBin->FindBin(fitterME[iPtBin]->GetMean()-3.*fitterME[iPtBin]->GetSigma());
1054  Double_t maxBinBkg=hMassPtBin->FindBin(fitterME[iPtBin]->GetMean()+3.*fitterME[iPtBin]->GetSigma());
1055  background=hMassPtBin->Integral(minBinBkg,maxBinBkg);
1056  hSoverBME->SetBinContent(iPtBin+1,fitterME[iPtBin]->GetRawYield()/background);
1057  hSoverBME->SetBinError(iPtBin+1,0.000001);
1058  hRelStatME->SetBinContent(iPtBin+1,fitterME[iPtBin]->GetRawYieldError()/fitterME[iPtBin]->GetRawYield());
1059  hRelStatME->SetBinError(iPtBin+1,0.000001);
1060  hSignifME->SetBinContent(iPtBin+1,fitterME[iPtBin]->GetRawYield()/TMath::Sqrt(background+fitterME[iPtBin]->GetRawYield()));
1061  hSignifME->SetBinError(iPtBin+1,0.00000001);
1062  hInvMassHistoBinWidthME->SetBinContent(iPtBin+1,hMassSubME->GetBinWidth(1));
1063  hGausMeanME->SetBinContent(iPtBin+1,fitterME[iPtBin]->GetMean());
1064  hGausMeanME->SetBinError(iPtBin+1,fitterME[iPtBin]->GetMeanUncertainty());
1065  hGausSigmaME->SetBinContent(iPtBin+1,fitterME[iPtBin]->GetSigma());
1066  hGausSigmaME->SetBinError(iPtBin+1,fitterME[iPtBin]->GetSigmaUncertainty());
1067  hChiSqME->SetBinContent(iPtBin+1,fitterME[iPtBin]->GetReducedChiSquare());
1068  hChiSqME->SetBinError(iPtBin+1,0.00001); // very small number, for graphics
1069  hNdfME->SetBinContent(iPtBin+1,fitterME[iPtBin]->GetMassFunc()->GetNDF());
1070  hNdfME->SetBinError(iPtBin+1,0.00001); // very small number, for graphics
1071  // if(!correctForRefl)
1072  WriteFitInfo(fitterME[iPtBin],hMassPtBin);
1073  fout->cd();
1074  WriteFitFunctionsToFile(fitterME[iPtBin],"ME",iPtBin);
1075  current->cd();
1076 
1077  Double_t errbc;
1078  Double_t bc=fitterME[iPtBin]->GetRawYieldBinCounting(errbc,nsigmaBinCounting,optBkgBinCount);
1079  hRawYieldMEBC->SetBinContent(iPtBin+1,bc);
1080  hRawYieldMEBC->SetBinError(iPtBin+1,errbc);
1081 
1082 
1083  c4pulls->cd(iPtBin+1);
1084  TH1F *hPullsTrend=new TH1F();// the name is changed afterwards, histo must not be deleted
1085  TH1F *hPulls=new TH1F();//the name is changed afterwards, histo must not be deleted
1086  TH1F *hResidualTrend=new TH1F();// the name is changed afterwards, histo must not be deleted
1087  TH1F *hResiduals=fitterME[iPtBin]->GetResidualsAndPulls(hPulls,hResidualTrend,hPullsTrend);
1088  hResiduals->SetName(Form("hResidualsME_PtBin%d",iPtBin));
1089  hPulls->Draw();
1090  PrintGausParams(hPulls);
1091 
1092  c4residuals->cd(iPtBin+1);
1093  hResiduals->Draw();
1094 
1095  SetStyleHisto(hResidualTrend,kME);
1096  c4residualTrend->cd(iPtBin+1);
1097  hResidualTrend->Draw();
1098 
1099  cCompareResidualTrends->cd(iPtBin+1);
1100  if(out4||out1||out2)hResidualTrend->Draw("same");
1101  else hResidualTrend->Draw();
1102 
1103  c4pullTrend->cd(iPtBin+1);
1104  hPullsTrend->Draw();
1105  }
1106  else hMassSubME->Draw("");
1107  if(correctForRefl){
1108  delete hMCReflPtBin;
1109  delete hMCSigPtBin;
1110  }
1111 
1112  }
1113  TString path(gSystem->pwd());
1114  path.Append("/figures");
1115  if(gSystem->AccessPathName(path.Data())){
1116  gROOT->ProcessLine(Form(".!mkdir -p %s",path.Data()));
1117  }
1118 
1119  if(saveCanvasAsEps>0){
1120  c1->SaveAs(Form("figures/InvMassSpectra_%s_%s_NoBkgSub.eps",sigConf.Data(),suffix.Data()));
1121  c2->SaveAs(Form("figures/InvMassSpectra_%s_%s_Rot.eps",sigConf.Data(),suffix.Data()));
1122  c3->SaveAs(Form("figures/InvMassSpectra_%s_%s_LS.eps",sigConf.Data(),suffix.Data()));
1123  c4->SaveAs(Form("figures/InvMassSpectra_%s_%s_EM.eps",sigConf.Data(),suffix.Data()));
1124  if(tryDirectFit){
1125  c5->SaveAs(Form("figures/InvMassSpectra_%s_%s_SB.eps",sigConf.Data(),suffix.Data()));
1126  c5sub->SaveAs(Form("figures/InvMassSpectra_%s_%s_SBsub.eps",sigConf.Data(),suffix.Data()));
1127  }
1128 
1129  if(saveCanvasAsEps>1){
1130  c2residuals->SaveAs(Form("figures/ResidualDistribution_%s_%s_Rot.eps",sigConf.Data(),suffix.Data()));
1131  c3residuals->SaveAs(Form("figures/ResidualDistribution_%s_%s_LS.eps",sigConf.Data(),suffix.Data()));
1132  c4residuals->SaveAs(Form("figures/ResidualDistribution_%s_%s_EM.eps",sigConf.Data(),suffix.Data()));
1133 
1134  c2residualTrend->SaveAs(Form("figures/residualTrendvsMass_%s_%s_Rot.eps",sigConf.Data(),suffix.Data()));
1135  c3residualTrend->SaveAs(Form("figures/residualTrendvsMass_%s_%s_LS.eps",sigConf.Data(),suffix.Data()));
1136  c4residualTrend->SaveAs(Form("figures/residualTrendvsMass_%s_%s_EM.eps",sigConf.Data(),suffix.Data()));
1137 
1138  c2pulls->SaveAs(Form("figures/PullDistribution_%s_%s_Rot.eps",sigConf.Data(),suffix.Data()));
1139  c3pulls->SaveAs(Form("figures/PullDistribution_%s_%s_LS.eps",sigConf.Data(),suffix.Data()));
1140  c4pulls->SaveAs(Form("figures/PullDistribution_%s_%s_EM.eps",sigConf.Data(),suffix.Data()));
1141 
1142  c2pullTrend->SaveAs(Form("figures/pullTrendvsMass_%s_%s_Rot.eps",sigConf.Data(),suffix.Data()));
1143  c3pullTrend->SaveAs(Form("figures/pullTrendvsMass_%s_%s_LS.eps",sigConf.Data(),suffix.Data()));
1144  c4pullTrend->SaveAs(Form("figures/pullTrendvsMass_%s_%s_EM.eps",sigConf.Data(),suffix.Data()));
1145 
1146  if(tryDirectFit){
1147  c5residuals->SaveAs(Form("figures/ResidualDistribution_%s_%s_SB.eps",sigConf.Data(),suffix.Data()));
1148  c5residualTrend->SaveAs(Form("figures/residualTrendvsMass_%s_%s_SB.eps",sigConf.Data(),suffix.Data()));
1149  c5pulls->SaveAs(Form("figures/PullDistribution_%s_%s_SB.eps",sigConf.Data(),suffix.Data()));
1150  c5pullTrend->SaveAs(Form("figures/pullTrendvsMass_%s_%s_SB.eps",sigConf.Data(),suffix.Data()));
1151  }
1152  }
1153  }
1154 
1155  // save also .root
1156  if(saveCanvasAsRoot){
1157  c2->SaveAs(Form("figures/InvMassSpectra_%s_%s_Rot.root",sigConf.Data(),suffix.Data()));
1158  c3->SaveAs(Form("figures/InvMassSpectra_%s_%s_LS.root",sigConf.Data(),suffix.Data()));
1159  c4->SaveAs(Form("figures/InvMassSpectra_%s_%s_EM.root",sigConf.Data(),suffix.Data()));
1160 
1161  c2residuals->SaveAs(Form("figures/ResidualDistribution_%s_%s_Rot.root",sigConf.Data(),suffix.Data()));
1162  c3residuals->SaveAs(Form("figures/ResidualDistribution_%s_%s_LS.root",sigConf.Data(),suffix.Data()));
1163  c4residuals->SaveAs(Form("figures/ResidualDistribution_%s_%s_EM.root",sigConf.Data(),suffix.Data()));
1164 
1165  c2residualTrend->SaveAs(Form("figures/residualTrendvsMass_%s_%s_Rot.root",sigConf.Data(),suffix.Data()));
1166  c3residualTrend->SaveAs(Form("figures/residualTrendvsMass_%s_%s_LS.root",sigConf.Data(),suffix.Data()));
1167  c4residualTrend->SaveAs(Form("figures/residualTrendvsMass_%s_%s_EM.root",sigConf.Data(),suffix.Data()));
1168 
1169  c2pulls->SaveAs(Form("figures/PullDistribution_%s_%s_Rot.root",sigConf.Data(),suffix.Data()));
1170  c3pulls->SaveAs(Form("figures/PullDistribution_%s_%s_LS.root",sigConf.Data(),suffix.Data()));
1171  c4pulls->SaveAs(Form("figures/PullDistribution_%s_%s_EM.root",sigConf.Data(),suffix.Data()));
1172 
1173  c2pullTrend->SaveAs(Form("figures/pullTrendvsMass_%s_%s_Rot.root",sigConf.Data(),suffix.Data()));
1174  c3pullTrend->SaveAs(Form("figures/pullTrendvsMass_%s_%s_LS.root",sigConf.Data(),suffix.Data()));
1175  c4pullTrend->SaveAs(Form("figures/pullTrendvsMass_%s_%s_EM.root",sigConf.Data(),suffix.Data()));
1176 
1177  if(tryDirectFit){
1178  c5->SaveAs(Form("figures/InvMassSpectra_%s_%s_SB.root",sigConf.Data(),suffix.Data()));
1179  c5sub->SaveAs(Form("figures/InvMassSpectra_%s_%s_SBsub.root",sigConf.Data(),suffix.Data()));
1180  c5residuals->SaveAs(Form("figures/ResidualDistribution_%s_%s_SB.root",sigConf.Data(),suffix.Data()));
1181  c5residualTrend->SaveAs(Form("figures/residualTrendvsMass_%s_%s_SB.root",sigConf.Data(),suffix.Data()));
1182  c5pulls->SaveAs(Form("figures/PullDistribution_%s_%s_SB.root",sigConf.Data(),suffix.Data()));
1183  c5pullTrend->SaveAs(Form("figures/pullTrendvsMass_%s_%s_SB.root",sigConf.Data(),suffix.Data()));
1184  }
1185  }
1186 
1187 
1188  hRawYieldRot->SetMarkerStyle(21);
1189  hRawYieldLS->SetMarkerStyle(22);
1190  hRawYieldLS->SetMarkerColor(kGreen+2);
1191  hRawYieldLS->SetLineColor(kGreen+2);
1192  hRawYieldME->SetMarkerStyle(25);
1193  hRawYieldME->SetMarkerColor(4);
1194  hRawYieldME->SetLineColor(4);
1195  hRawYieldSB->SetMarkerStyle(27);
1196  hRawYieldSB->SetMarkerColor(6);
1197  hRawYieldSB->SetLineColor(6);
1198 
1199 
1200  TCanvas* cry=new TCanvas("cry","RawYield",800,700);
1201  cry->SetLeftMargin(0.15);
1202  hRawYieldRot->Draw("P");
1203  hRawYieldRot->SetMinimum(0);
1204  hRawYieldRot->GetYaxis()->SetTitleOffset(1.8);
1205  Double_t max=hRawYieldRot->GetMaximum();
1206  if(hRawYieldLS->GetMaximum()>max)max=hRawYieldLS->GetMaximum();
1207  if(hRawYieldME->GetMaximum()>max)max=hRawYieldME->GetMaximum();
1208  if(tryDirectFit){
1209  if(hRawYieldSB->GetMaximum()>max)max=hRawYieldSB->GetMaximum();
1210  }
1211  hRawYieldRot->SetMaximum(max*1.2);
1212  hRawYieldLS->Draw("PZSAME");
1213  hRawYieldME->Draw("PSAME");
1214  if(tryDirectFit) hRawYieldSB->Draw("PSAME");
1215  TLegend* legry=new TLegend(0.7,0.7,0.89,0.89);
1216  legry->SetFillStyle(0);
1217  legry->SetBorderSize(0);
1218  legry->AddEntry(hRawYieldRot,"Rotational","PL")->SetTextColor(1);
1219  legry->AddEntry(hRawYieldLS,"Like Sign","PL")->SetTextColor(kGreen+2);
1220  legry->AddEntry(hRawYieldME,"Ev Mix","PL")->SetTextColor(4);
1221  if(tryDirectFit){
1222  legry->AddEntry(hRawYieldSB,"Side-Band Fit","PL")->SetTextColor(6);
1223  }
1224  legry->Draw();
1225  if(saveCanvasAsEps>0) cry->SaveAs(Form("figures/RawYield_%s_%s.eps",sigConf.Data(),suffix.Data()));
1226  if(saveCanvasAsRoot) cry->SaveAs(Form("figures/RawYield_%s_%s.root",sigConf.Data(),suffix.Data()));
1227 
1228 
1229 
1230  TCanvas* cch2=new TCanvas("cch2","Chi2",800,700);
1231  cch2->SetLeftMargin(0.15);
1232  hChiSqRot->SetMarkerStyle(21);
1233  hChiSqRot->Draw("P");
1234  hChiSqRot->SetMinimum(0);
1235  hChiSqRot->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1236  hChiSqRot->GetYaxis()->SetTitle("#chi^{2}/ndf");
1237  hChiSqRot->GetYaxis()->SetTitleOffset(1.8);
1238  Double_t maxc=hChiSqRot->GetMaximum();
1239  if(hChiSqLS->GetMaximum()>maxc)maxc=hChiSqLS->GetMaximum();
1240  if(hChiSqME->GetMaximum()>maxc)maxc=hChiSqME->GetMaximum();
1241  if(tryDirectFit){
1242  if(hChiSqSB->GetMaximum()>maxc)maxc=hChiSqSB->GetMaximum();
1243  }
1244  hChiSqRot->SetMaximum(maxc*1.2);
1245 
1246  hChiSqLS->SetMarkerStyle(22);
1247  hChiSqLS->SetMarkerColor(kGreen+2);
1248  hChiSqLS->SetLineColor(kGreen+2);
1249  hChiSqLS->Draw("PZSAME");
1250  hChiSqME->SetMarkerStyle(25);
1251  hChiSqME->SetMarkerColor(4);
1252  hChiSqME->SetLineColor(4);
1253  hChiSqME->Draw("PSAME");
1254  if(tryDirectFit){
1255  hChiSqSB->SetMarkerStyle(27);
1256  hChiSqSB->SetMarkerColor(6);
1257  hChiSqSB->SetLineColor(6);
1258  hChiSqSB->Draw("PSAME");
1259  }
1260  legry->Draw();
1261  if(saveCanvasAsEps>0) cch2->SaveAs(Form("figures/ChiSq_%s_%s.eps",sigConf.Data(),suffix.Data()));
1262  if(saveCanvasAsRoot) cch2->SaveAs(Form("figures/ChiSq_%s_%s.root",sigConf.Data(),suffix.Data()));
1263 
1264  TH1F* hRatioLSToME=(TH1F*)hRawYieldLS->Clone("hRatioLStoME");
1265  TH1F* hRatioRotToME=(TH1F*)hRawYieldRot->Clone("hRatioRottoME");
1266  TH1F* hRatioMEToME=(TH1F*)hRawYieldME->Clone("hRatioMEtoME");
1267  TH1F* hRatioSBToME=(TH1F*)hRawYieldSB->Clone("hRatioSBtoME");
1268  for(Int_t ib=1; ib<=hRawYieldME->GetNbinsX(); ib++){
1269  Double_t yme=hRawYieldME->GetBinContent(ib);
1270  if(yme>0.){
1271  hRatioLSToME->SetBinContent(ib,hRawYieldLS->GetBinContent(ib)/yme);
1272  hRatioLSToME->SetBinError(ib,hRawYieldLS->GetBinError(ib)/yme);
1273  hRatioMEToME->SetBinContent(ib,hRawYieldME->GetBinContent(ib)/yme);
1274  hRatioMEToME->SetBinError(ib,hRawYieldME->GetBinError(ib)/yme);
1275  hRatioRotToME->SetBinContent(ib,hRawYieldRot->GetBinContent(ib)/yme);
1276  hRatioRotToME->SetBinError(ib,hRawYieldRot->GetBinError(ib)/yme);
1277  hRatioSBToME->SetBinContent(ib,hRawYieldSB->GetBinContent(ib)/yme);
1278  hRatioSBToME->SetBinError(ib,hRawYieldSB->GetBinError(ib)/yme);
1279  }
1280  }
1281 
1282  TCanvas* cry2=new TCanvas("cry2","RawYield+Ratios",1400,700);
1283  cry2->Divide(2,1);
1284  cry2->cd(1);
1285  gPad->SetLeftMargin(0.15);
1286  gPad->SetRightMargin(0.05);
1287  hRawYieldRot->Draw("P");
1288  hRawYieldLS->Draw("PZSAME");
1289  hRawYieldME->Draw("PSAME");
1290  if(tryDirectFit){
1291  hRawYieldSB->Draw("PSAME");
1292  }
1293  legry->Draw();
1294  cry2->cd(2);
1295  hRatioLSToME->SetStats(0);
1296  hRatioLSToME->SetMinimum(0.3);
1297  hRatioLSToME->SetMaximum(1.7);
1298  hRatioLSToME->GetYaxis()->SetTitle("Ratio To EvMix");
1299  hRatioLSToME->Draw("same");
1300  hRatioRotToME->Draw("same");
1301  hRatioMEToME->Draw("same");
1302  hRatioSBToME->Draw("same");
1303  if(saveCanvasAsEps>0) cry2->SaveAs(Form("figures/RawYieldAndRatios_%s_%s.eps",sigConf.Data(),suffix.Data()));
1304  if(saveCanvasAsRoot) cry2->SaveAs(Form("figures/RawYieldAndRatios_%s_%s.root",sigConf.Data(),suffix.Data()));
1305 
1306  hRelStatRot->GetYaxis()->SetTitleOffset(1.8);
1307  hRelStatRot->SetMarkerStyle(21);
1308  hRelStatLS->SetMarkerStyle(22);
1309  hRelStatLS->SetMarkerColor(kGreen+2);
1310  hRelStatLS->SetLineColor(kGreen+2);
1311  hRelStatME->SetMarkerStyle(25);
1312  hRelStatME->SetMarkerColor(4);
1313  hRelStatME->SetLineColor(4);
1314  hRelStatSB->SetMarkerStyle(27);
1315  hRelStatSB->SetMarkerColor(6);
1316  hRelStatSB->SetLineColor(6);
1317 
1318  TCanvas* cry3=new TCanvas("cry3","RawYield+Ratios+Unc",1800,600);
1319  cry3->Divide(3,1);
1320  cry3->cd(1);
1321  gPad->SetLeftMargin(0.15);
1322  gPad->SetRightMargin(0.05);
1323  hRawYieldRot->Draw("P");
1324  hRawYieldLS->Draw("PZSAME");
1325  hRawYieldME->Draw("PSAME");
1326  if(tryDirectFit){
1327  hRawYieldSB->Draw("PSAME");
1328  }
1329  legry->Draw();
1330  cry3->cd(2);
1331  hRatioLSToME->SetStats(0);
1332  hRatioLSToME->SetMinimum(0.3);
1333  hRatioLSToME->SetMaximum(1.7);
1334  hRatioLSToME->GetYaxis()->SetTitle("Ratio To EvMix");
1335  hRatioLSToME->Draw("same");
1336  hRatioRotToME->Draw("same");
1337  hRatioMEToME->Draw("same");
1338  hRatioSBToME->Draw("same");
1339  cry3->cd(3);
1340  gPad->SetLeftMargin(0.15);
1341  gPad->SetRightMargin(0.05);
1342  hRelStatRot->SetStats(0);
1343  hRelStatRot->SetMinimum(0.04);
1344  hRelStatRot->SetMaximum(0.2);
1345  hRelStatRot->Draw();
1346  hRelStatLS->Draw("same");
1347  hRelStatME->Draw("same");
1348  if(tryDirectFit) hRelStatSB->Draw("same");
1349 
1350  if(saveCanvasAsEps>0) cry3->SaveAs(Form("figures/RawYieldRatiosAndUnc_%s_%s.eps",sigConf.Data(),suffix.Data()));
1351  if(saveCanvasAsRoot) cry3->SaveAs(Form("figures/RawYieldRatiosAndUnc_%s_%s.root",sigConf.Data(),suffix.Data()));
1352 
1353 
1354 
1355 
1356  cCompareResidualTrends->cd(1);
1357  TLegend *legRT=new TLegend(*legry);
1358  legRT->SetX1NDC(0.12);
1359  legRT->SetY1NDC(0.7);
1360  legRT->SetX2NDC(0.3);
1361  legRT->SetY2NDC(0.9);
1362  legRT->Draw();
1363 
1364  // NOW COMPARE YIELDS OBTAINED WITH BC FROM DIFFERENT APPROACHES
1365  TCanvas* cryBC=new TCanvas("cryBC","RawYield with BC",800,700);
1366  cryBC->SetLeftMargin(0.15);
1367  hRawYieldRot->Draw("P");
1368  hRawYieldME->Draw("PSAME");
1369  hRawYieldLS->Draw("PZSAME");
1370  if(tryDirectFit){
1371  hRawYieldSB->Draw("PSAME");
1372  }
1373 
1374  if(hRawYieldRotBC->GetMaximum()>max)max=hRawYieldRotBC->GetMaximum();
1375  if(hRawYieldLSBC->GetMaximum()>max)max=hRawYieldLSBC->GetMaximum();
1376  if(hRawYieldMEBC->GetMaximum()>max)max=hRawYieldMEBC->GetMaximum();
1377  if(tryDirectFit){
1378  if(hRawYieldSBBC->GetMaximum()>max)max=hRawYieldSBBC->GetMaximum();
1379  }
1380  hRawYieldRot->SetMaximum(max*1.2);
1381 
1382  hRawYieldRotBC->SetMarkerStyle(21);
1383  hRawYieldRotBC->SetLineStyle(2);
1384  hRawYieldRotBC->Draw("PSAME");
1385 
1386 
1387  hRawYieldLSBC->SetMarkerStyle(22);
1388  hRawYieldLSBC->SetMarkerColor(kGreen+2);
1389  hRawYieldLSBC->SetLineColor(kGreen+2);
1390  hRawYieldLSBC->SetLineStyle(2);
1391  hRawYieldLSBC->Draw("PZSAME");
1392 
1393  hRawYieldMEBC->SetMarkerStyle(25);
1394  hRawYieldMEBC->SetMarkerColor(4);
1395  hRawYieldMEBC->SetLineColor(4);
1396  hRawYieldMEBC->SetLineStyle(2);
1397  hRawYieldMEBC->Draw("PSAME");
1398  if(tryDirectFit){
1399  hRawYieldSBBC->SetMarkerStyle(27);
1400  hRawYieldSBBC->SetMarkerColor(6);
1401  hRawYieldSBBC->SetLineColor(6);
1402  hRawYieldSBBC->SetLineStyle(2);
1403  hRawYieldSBBC->Draw("PSAME");
1404  }
1405 
1406  TLegend* legryBC=new TLegend(0.7,0.7,0.89,0.89);
1407  legryBC->SetFillStyle(0);
1408  legryBC->SetBorderSize(0);
1409  legryBC->AddEntry(hRawYieldRot,"Rotational","PL")->SetTextColor(1);
1410  legryBC->AddEntry(hRawYieldRotBC,"Rotational BC","PL")->SetTextColor(1);
1411  legryBC->AddEntry(hRawYieldLS,"Like Sign ","PL")->SetTextColor(kGreen+2);
1412  legryBC->AddEntry(hRawYieldLSBC,"Like Sign BC","PL")->SetTextColor(kGreen+2);
1413  legryBC->AddEntry(hRawYieldME,"Ev Mix","PL")->SetTextColor(4);
1414  legryBC->AddEntry(hRawYieldMEBC,"Ev Mix BC","PL")->SetTextColor(4);
1415  if(tryDirectFit){
1416  legryBC->AddEntry(hRawYieldSB,"Side-Band Fit","PL")->SetTextColor(6);
1417  legryBC->AddEntry(hRawYieldSBBC,"Side-Band Fit BC","PL")->SetTextColor(6);
1418  }
1419  legryBC->Draw();
1420  if(saveCanvasAsEps>0) cryBC->SaveAs(Form("figures/RawYieldBC_%s_%s.eps",sigConf.Data(),suffix.Data()));
1421  if(saveCanvasAsRoot) cryBC->SaveAs(Form("figures/RawYieldBC_%s_%s.root",sigConf.Data(),suffix.Data()));
1422 
1423 
1424 
1425 
1426  fout->cd();
1427  hRawYieldRot->Write();
1428  hRawYieldLS->Write();
1429  hRawYieldME->Write();
1430  hRawYieldSB->Write();
1431  hRelStatRot->Write();
1432  hRelStatLS->Write();
1433  hRelStatME->Write();
1434  hRelStatSB->Write();
1435  hSignifRot->Write();
1436  hSignifLS->Write();
1437  hSignifME->Write();
1438  hSignifSB->Write();
1439  hSoverBRot->Write();
1440  hSoverBLS->Write();
1441  hSoverBME->Write();
1442  hSoverBSB->Write();
1443  hGausMeanRot->Write();
1444  hGausMeanLS->Write();
1445  hGausMeanME->Write();
1446  hGausMeanSB->Write();
1447  hGausSigmaRot->Write();
1448  hGausSigmaLS->Write();
1449  hGausSigmaME->Write();
1450  hGausSigmaSB->Write();
1451  hChiSqRot->Write();
1452  hChiSqLS->Write();
1453  hChiSqME->Write();
1454  hChiSqSB->Write();
1455  hNdfRot->Write();
1456  hNdfLS->Write();
1457  hNdfME->Write();
1458  hNdfSB->Write();
1459  hInvMassHistoBinWidthRot->Write();
1460  hInvMassHistoBinWidthLS->Write();
1461  hInvMassHistoBinWidthME->Write();
1462  hInvMassHistoBinWidthSB->Write();
1463  hRebin->Write();
1464  hBkgFitFunc->Write();
1465  hBkgFitFuncSB->Write();
1466  hnEv->Write();
1467  if(hSigmaMC) hSigmaMC->Write();
1468  fout->Close();
1469 
1470  return;
1471 }
1472 
1473 void WriteFitFunctionsToFile(AliHFInvMassFitter *fitter, TString meth, Int_t iPtBin){
1474  TString nameh;
1475  TF1* fTot=fitter->GetMassFunc();
1476  nameh.Form("FitFuncTot_%s_PtBin%d",meth.Data(),iPtBin);
1477  fTot->SetRange(1.6,2.2);
1478  fTot->SetNpx(500);
1479  fTot->Write(nameh.Data());
1480  TF1* fSig=fitter->GetSignalFunc();
1481  nameh.Form("FitFuncSig_%s_PtBin%d",meth.Data(),iPtBin);
1482  fSig->SetRange(1.6,2.2);
1483  fSig->SetNpx(500);
1484  fSig->Write(nameh.Data());
1485  TF1* fBkg=fitter->GetBackgroundRecalcFunc();
1486  nameh.Form("FitFuncBkg_%s_PtBin%d",meth.Data(),iPtBin);
1487  fBkg->SetRange(1.6,2.2);
1488  fBkg->SetNpx(500);
1489  fBkg->Write(nameh.Data());
1490  if(meson=="Dzero"){
1491  TF1* fBkgR=fitter->GetBkgPlusReflFunc();
1492  nameh.Form("FitFuncBkgRefl_%s_PtBin%d",meth.Data(),iPtBin);
1493  fBkgR->SetRange(1.6,2.2);
1494  fBkgR->SetNpx(500);
1495  fBkgR->Write(nameh.Data());
1496  }
1497  return;
1498 }
1499 
1500 void WriteFitInfo(AliHFInvMassFitter *fitter, TH1D* histo){
1501  Double_t sig=fitter->GetRawYield();
1502  Double_t esig=fitter->GetRawYieldError();
1503  Double_t mean=fitter->GetMean();
1504  Double_t emean=fitter->GetMeanUncertainty();
1505  Double_t sigma=fitter->GetSigma();
1506  Double_t esigma=fitter->GetSigmaUncertainty();
1507  Double_t minBin=histo->FindBin(mean-3.*sigma);
1508  Double_t maxBin=histo->FindBin(mean+3.*sigma);
1509  Double_t back=histo->Integral(minBin,maxBin);
1510  TPaveText* tpar=new TPaveText(0.5,0.7,0.89,.87,"NDC");
1511  tpar->SetBorderSize(0);
1512  tpar->SetFillStyle(0);
1513  tpar->AddText(Form("Mean = %.3f #pm %.3f",mean,emean));
1514  tpar->AddText(Form("Sigma = %.3f #pm %.3f",sigma,esigma));
1515  tpar->SetTextColor(4);
1516  tpar->Draw();
1517 
1518  TPaveText* tss=new TPaveText(0.15,0.15,0.5,0.4,"NDC");
1519  tss->SetBorderSize(0);
1520  tss->SetFillStyle(0);
1521  tss->AddText(Form("S = %.0f #pm %.0f",sig,esig));
1522  tss->AddText(Form("B(3#sigma) = %.3g",back));
1523  tss->AddText(Form("S/B (3#sigma) = %.4f",sig/back));
1524  if(correctForRefl) tss->AddText(Form("Refl/Sig = %.3f #pm %.3f ",fitter->GetReflOverSig(),fitter->GetReflOverSigUncertainty()));
1525  tss->AddText(Form("Significance(3#sigma) = %.2f",sig/TMath::Sqrt(back+sig)));
1526  tss->SetTextColor(1);
1527  tss->Draw();
1528 
1529 }
1530 
1531 TH1F* FitMCInvMassSpectra(TList* lMC){
1532  TH3F* h3d=(TH3F*)lMC->FindObject("hMassVsPtVsYSig");
1533  if(!h3d){
1534  printf("hMassVsPtVsYSig not found\n");
1535  return 0x0;
1536  }
1537  TH3F* h3dr=(TH3F*)lMC->FindObject("hMassVsPtVsYRefl");
1538  if(!h3dr){
1539  printf("hMassVsPtVsYRefl not found\n");
1540  }
1541  TH1F* hSigmaMC=new TH1F("hSigmaMC","",nPtBins,binLims);
1542 
1543  TCanvas* cmc1=new TCanvas("InvMassMC","InvMassMC",1200,800);
1544  cmc1->Divide(4,2);
1545 
1546  gStyle->SetOptFit(0);
1547  gStyle->SetOptStat(0);
1548  for(Int_t iPtBin=0; iPtBin<nPtBins; iPtBin++){
1549  Int_t bin1=h3d->GetYaxis()->FindBin(binLims[iPtBin]);
1550  Int_t bin2=h3d->GetYaxis()->FindBin(binLims[iPtBin+1]-0.0001);
1551  Double_t ptmed=(binLims[iPtBin]+binLims[iPtBin+1])*0.5;
1552  Double_t pthalfwid=(binLims[iPtBin+1]-binLims[iPtBin])*0.5;
1553  printf("Bin %d Pt range=%f %f\n",iPtBin,h3d->GetYaxis()->GetBinLowEdge(bin1),h3d->GetYaxis()->GetBinUpEdge(bin2));
1554  TH1D* hMassMCPtBin=h3d->ProjectionX(Form("hMassMCPtBin%d",iPtBin),bin1,bin2);
1555  hMassMCPtBin->SetTitle(Form("%.1f<p_{T}<%.1f GeV/c",binLims[iPtBin],binLims[iPtBin+1]));
1556  hMassMCPtBin->GetXaxis()->SetTitle("Invariant mass (GeV/c^{2})");
1557  hMassMCPtBin->GetXaxis()->SetRangeUser(1.65,2.06);
1558  cmc1->cd(iPtBin+1);
1559  gPad->SetLogy();
1560  hMassMCPtBin->Fit("gaus","");
1561  TF1* fg=(TF1*)hMassMCPtBin->GetListOfFunctions()->FindObject("gaus");
1562  TLatex* tsig=new TLatex(0.65,0.82,"Signal");
1563  tsig->SetNDC();
1564  tsig->SetTextColor(kBlue+1);
1565  tsig->Draw();
1566 
1567  TPaveText* t1=new TPaveText(0.15,0.55,0.4,0.88,"ndc");
1568  t1->SetFillStyle(0);
1569  t1->SetBorderSize(0);
1570  t1->AddText(Form("#mu = %.4f GeV/c^{2}",fg->GetParameter(1)));
1571  t1->AddText(Form("#sigma = %.4f GeV/c^{2}",fg->GetParameter(2)));
1572  t1->AddText(Form("Integral = %.0f",fg->Integral(1.7,2.1)/hMassMCPtBin->GetBinWidth(1)));
1573  t1->AddText(Form("Entries = %.0f",hMassMCPtBin->GetEntries()));
1574 
1575  if(h3dr){
1576  TH1D* hReflPtBin=h3dr->ProjectionX(Form("hReflPtBin%d",iPtBin),bin1,bin2);
1577  hReflPtBin->SetTitle(Form("%.1f<p_{T}<%.1f GeV/c",binLims[iPtBin],binLims[iPtBin+1]));
1578  hReflPtBin->GetXaxis()->SetTitle("Invariant mass (GeV/c^{2})");
1579  hReflPtBin->SetLineColor(2);
1580  Int_t bin1=hReflPtBin->FindBin(fg->GetParameter(1)-3.*fg->GetParameter(2));
1581  Int_t bin2=hReflPtBin->FindBin(fg->GetParameter(1)+3.*fg->GetParameter(2));
1582  TLatex* tref=new TLatex(0.65,0.75,"Reflections");
1583  tref->SetNDC();
1584  tref->SetTextColor(2);
1585  tref->Draw();
1586  t1->AddText(Form("Refl(3#sigma) = %.0f",hReflPtBin->Integral(bin1,bin2)));
1587  hReflPtBin->Draw("same");
1588  t1->Draw();
1589  }
1590  sigmas[iPtBin]=fg->GetParameter(2);
1591  Double_t errMCsigma=fg->GetParError(2);
1592  if(tuneSigmaOnData>0.){
1593  sigmas[iPtBin]*=tuneSigmaOnData;
1594  errMCsigma*=tuneSigmaOnData;
1595  }
1596  hSigmaMC->SetBinContent(iPtBin+1,sigmas[iPtBin]);
1597  hSigmaMC->SetBinError(iPtBin+1,errMCsigma);
1598  }
1599  return hSigmaMC;
1600 }
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
AliVertexingHFUtils::RebinHisto
static TH1D * RebinHisto(TH1 *hOrig, Int_t reb, Int_t firstUse=-1)
Rebinning of invariant mass histograms.
Definition: AliVertexingHFUtils.cxx:2301
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:2345
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:60
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:1531
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
AliHFInvMassFitter::GetResidualsAndPulls
TH1F * GetResidualsAndPulls(TH1 *hPulls=0x0, TH1 *hResidualTrend=0x0, TH1 *hPullsTrend=0x0, Double_t minrange=0, Double_t maxrange=-1)
Definition: AliHFInvMassFitter.cxx:1154
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:1473
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:1500