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