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