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