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