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