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