AliPhysics  63e47e1 (63e47e1)
FitMassSpectra.C
Go to the documentation of this file.
1 #if !defined(__CINT__) || defined(__MAKECINT__)
2 #include <TInterpreter.h>
3 #include <TString.h>
4 #include <TObjString.h>
5 #include <TObjArray.h>
6 #include <TMath.h>
7 #include <TFile.h>
8 #include <TCanvas.h>
9 #include <TH1F.h>
10 #include <TH2F.h>
11 #include <TH1D.h>
12 #include <TF1.h>
13 #include <TStyle.h>
14 #include <TLegend.h>
15 #include <TLegendEntry.h>
16 #include <TDatabasePDG.h>
17 #include <TGraph.h>
18 
19 #include "AliAODRecoDecayHF.h"
20 #include "AliRDHFCuts.h"
23 #include "AliRDHFCutsDstoKKpi.h"
24 #include "AliRDHFCutsD0toKpi.h"
25 #include "AliHFMassFitter.h"
27 #endif
28 
29 
30 // MACRO to perform fits to D meson invariant mass spectra
31 // and store raw yields and cut object into a root output file
32 // Origin: F. Prino (prino@to.infn.it)
33 // D0: C. Bianchin (cbianchi@pd.infn.it)
34 
35 
36 
37 //
40 enum {kExpo=0, kLinear, kPol2, kThrExpo=5};
41 enum {kGaus=0, kDoubleGaus};
42 
43 
44 // Common variables: to be configured by the user
45 const Int_t nPtBins=6;
46 Double_t ptlims[nPtBins+1]={2.,3.,4.,5.,6.,8.,12.};
47 Int_t rebin[nPtBins]={2,4,4,4,4,4};
48 Int_t firstUsedBin[nPtBins]={-1,-1,-1,-1,-1,-1}; // -1 uses all bins, >=1 to set the lower bin to be accepted from original histo
49 
50 TString suffix="Loose_SecondSet1236_ForCF08";
51 
52 
53 //const Int_t nPtBins=7;//6;
54 //Double_t ptlims[nPtBins+1]={1.,2.,3.,4.,5.,6.,8.,12.};
55 //Int_t rebin[nPtBins+1]={8,6,10,10,10,10,10,10}; //for looser cuts
56 //Int_t rebin[nPtBins+1]={10,10,10,10,10,10,10,10}; //for Chiara's cuts
66 
67 //for D0only
69 const Int_t nsamples=2;//3;
70 //Int_t nevents[nsamples]={8.5635859e+07 /*LHC10b+cpart*/,8.9700624e+07/*LHC10dpart*/};
71 //Int_t nevents[nsamples]={9.0374946e+07 /*LHC10b+c*/,9.7593785e+07/*LHC10d*/};
72 //Int_t nevents[nsamples]={1.1777812e+08 /*LHC10dnewTPCpid*/,9.7593785e+07/*LHC10d*/};
73 //Int_t nevents[nsamples]={1.1777812e+08 /*LHC10dnewTPCpid*/,9.7593785e+07/*LHC10d*/,9.0374946e+07 /*LHC10b+c*/};
74 Int_t nevents[nsamples]={1.18860695e+08 /*LHC10dnewTPCpid*/,9.0374946e+07 /*LHC10b+c*/};
75 //Int_t nevents[nsamples]={1. /*LHC10dnewTPCpid*/,1 /*LHC10dnewTPCpidrescale*/};
76 // Functions
77 
78 Bool_t LoadDplusHistos(TObjArray* listFiles, TH1F** hMass);
79 Bool_t LoadDsHistos(TObjArray* listFiles, TH1F** hMass);
80 Bool_t LoadD0toKpiHistos(TObjArray* listFiles, TH1F** hMass);
81 Bool_t LoadDstarD0piHistos(TObjArray* listFiles, TH1F** hMass);
82 
83 
84 
85 
86 void FitMassSpectra(Int_t analysisType=kDplusKpipi,
87  TString fileNameb="",
88  TString fileNamec="",
89  TString fileNamed="",
90  TString fileNamee=""
91  ){
92  //
93 
94  gInterpreter->ExecuteMacro("$ALICE_ROOT/PWGHF/vertexingHF/macros/LoadLibraries.C");
95  gStyle->SetOptTitle(0);
96 
97  TObjArray* listFiles=new TObjArray();
98  if(fileNameb!="") listFiles->AddLast(new TObjString(fileNameb.Data()));
99  if(fileNamec!="") listFiles->AddLast(new TObjString(fileNamec.Data()));
100  if(fileNamed!="") listFiles->AddLast(new TObjString(fileNamed.Data()));
101  if(fileNamee!="") listFiles->AddLast(new TObjString(fileNamee.Data()));
102  if(listFiles->GetEntries()==0){
103  printf("Missing file names in input\n");
104  return;
105  }
106 
107  TH1F** hmass=new TH1F*[nPtBins];
108  for(Int_t i=0;i<nPtBins;i++) hmass[i]=0x0;
109 
110  Float_t massD, massD0_fDstar;
111  Bool_t retCode;
112  if(analysisType==kD0toKpi){
113  retCode=LoadD0toKpiHistos(listFiles,hmass);
114  massD=TDatabasePDG::Instance()->GetParticle(421)->Mass();
115  }
116  else if(analysisType==kDplusKpipi){
117  retCode=LoadDplusHistos(listFiles,hmass);
118  massD=TDatabasePDG::Instance()->GetParticle(411)->Mass();
119  }
120  else if(analysisType==kDStarD0pi){
121  retCode=LoadDstarD0piHistos(listFiles,hmass);
122  massD=TDatabasePDG::Instance()->GetParticle(413)->Mass();
123  massD0_fDstar=TDatabasePDG::Instance()->GetParticle(421)->Mass();
124  massD =massD-massD0_fDstar;
125  }
126  else if(analysisType==kDsKKpi){
127  retCode=LoadDsHistos(listFiles,hmass);
128  massD=TDatabasePDG::Instance()->GetParticle(431)->Mass();
129  }
130  else{
131  printf("Wrong analysis type parameter\n");
132  return;
133  }
134  if(!retCode){
135  printf("ERROR in reading input files\n");
136  return;
137  }
138 
139 
140 
141  TH1D* hCntSig1=new TH1D("hCntSig1","hCntSig1",nPtBins,ptlims);
142  TH1D* hCntSig2=new TH1D("hCntSig2","hCntSig2",nPtBins,ptlims);
143  TH1D* hNDiffCntSig1=new TH1D("hNDiffCntSig1","hNDiffCntSig1",nPtBins,ptlims);
144  TH1D* hNDiffCntSig2=new TH1D("hNDiffCntSig2","hNDiffCntSig2",nPtBins,ptlims);
145  TH1D* hSignal=new TH1D("hSignal","hSignal",nPtBins,ptlims);
146  TH1D* hRelErrSig=new TH1D("hRelErrSig","hRelErrSig",nPtBins,ptlims);
147  TH1D* hInvSignif=new TH1D("hInvSignif","hInvSignif",nPtBins,ptlims);
148  TH1D* hBackground=new TH1D("hBackground","hBackground",nPtBins,ptlims);
149  TH1D* hBackgroundNormSigma=new TH1D("hBackgroundNormSigma","hBackgroundNormSigma",nPtBins,ptlims);
150  TH1D* hSignificance=new TH1D("hSignificance","hSignificance",nPtBins,ptlims);
151  TH1D* hMass=new TH1D("hMass","hMass",nPtBins,ptlims);
152  TH1D* hSigma=new TH1D("hSigma","hSigma",nPtBins,ptlims);
153 
154 
155  Int_t nMassBins=hmass[0]->GetNbinsX();
156  Double_t hmin=TMath::Max(minMassForFit,hmass[0]->GetBinLowEdge(2));
157  Double_t hmax=TMath::Min(maxMassForFit,hmass[0]->GetBinLowEdge(nMassBins-2)+hmass[0]->GetBinWidth(nMassBins-2));
158  Float_t minBinSum=hmass[0]->FindBin(massD-massRangeForCounting);
159  Float_t maxBinSum=hmass[0]->FindBin(massD+massRangeForCounting);
160  Int_t iPad=1;
161 
162  TF1* funBckStore1=0x0;
163  TF1* funBckStore2=0x0;
164  TF1* funBckStore3=0x0;
165 
167  Double_t arrchisquare[nPtBins];
168  TCanvas* c1= new TCanvas("c1","MassSpectra",1000,800);
169  Int_t nx=3, ny=2;
170  if(nPtBins>6){
171  nx=4;
172  ny=3;
173  }
174  if(nPtBins>12){
175  nx=5;
176  ny=4;
177  }
178  c1->Divide(nx,ny);
179 
180  Double_t sig,errsig,s,errs,b,errb;
181  for(Int_t iBin=0; iBin<nPtBins; iBin++){
182  hmass[iBin]->SetName(Form("hInvMass_PtBin%d",iBin));
183  hmass[iBin]->SetTitle(Form("%.1f<pt<%.1f",ptlims[iBin],ptlims[iBin+1]));
184  c1->cd(iPad++);
185  Int_t origNbins=hmass[iBin]->GetNbinsX();
186  TH1F* hRebinned=(TH1F*)AliVertexingHFUtils::RebinHisto(hmass[iBin],rebin[iBin],firstUsedBin[iBin]);
187  hRebinned->GetXaxis()->SetTitle("Invariant Mass (GeV/c^{2})");
188  hRebinned->GetYaxis()->SetTitle(Form("Entries/(%.0f MeV/c^{2})",(hRebinned->GetBinWidth(1)*1000)));
189  hRebinned->GetYaxis()->SetTitleOffset(1.1);
190  hmin=TMath::Max(minMassForFit,hRebinned->GetBinLowEdge(2));
191  hmax=TMath::Min(maxMassForFit,hRebinned->GetBinLowEdge(hRebinned->GetNbinsX()));
192  fitter[iBin]=new AliHFMassFitter(hRebinned,hmin, hmax,1,typeb,types);
193  rebin[iBin]=origNbins/fitter[iBin]->GetBinN();
194  fitter[iBin]->SetReflectionSigmaFactor(factor4refl);
195  fitter[iBin]->SetInitialGaussianMean(massD);
196  if(analysisType==kDStarD0pi) fitter[iBin]->SetInitialGaussianSigma(0.00065);
197  Bool_t out=fitter[iBin]->MassFitter(0);
198  if(!out) {
199  fitter[iBin]->GetHistoClone()->Draw();
200  continue;
201  }
202  Double_t mass=fitter[iBin]->GetMean();
203  Double_t sigma=fitter[iBin]->GetSigma();
204  arrchisquare[iBin]=fitter[iBin]->GetReducedChiSquare();
205  TF1* fB1=fitter[iBin]->GetBackgroundFullRangeFunc();
206  TF1* fB2=fitter[iBin]->GetBackgroundRecalcFunc();
207  TF1* fM=fitter[iBin]->GetMassFunc();
208  if(iBin==0 && fB1) funBckStore1=new TF1(*fB1);
209  if(iBin==0 && fB2) funBckStore2=new TF1(*fB2);
210  if(iBin==0 && fM) funBckStore3=new TF1(*fM);
211 
212  fitter[iBin]->DrawHere(gPad);
213  fitter[iBin]->Signal(3,s,errs);
214  fitter[iBin]->Background(3,b,errb);
215  fitter[iBin]->Significance(3,sig,errsig);
216  Double_t ry=fitter[iBin]->GetRawYield();
217  Double_t ery=fitter[iBin]->GetRawYieldError();
218  Float_t cntSig1=0.;
219  Float_t cntSig2=0.;
220  Float_t cntErr=0.;
221  for(Int_t iMB=minBinSum; iMB<=maxBinSum; iMB++){
222  Float_t bkg1=fB1 ? fB1->Eval(hmass[iBin]->GetBinCenter(iMB))/rebin[iBin] : 0;
223  Float_t bkg2=fB2 ? fB2->Eval(hmass[iBin]->GetBinCenter(iMB))/rebin[iBin] : 0;
224  cntSig1+=(hmass[iBin]->GetBinContent(iMB)-bkg1);
225  cntSig2+=(hmass[iBin]->GetBinContent(iMB)-bkg2);
226  cntErr+=(hmass[iBin]->GetBinContent(iMB));
227  }
228  hCntSig1->SetBinContent(iBin+1,cntSig1);
229  hCntSig1->SetBinError(iBin+1,TMath::Sqrt(cntErr));
230  hNDiffCntSig1->SetBinContent(iBin+1,(s-cntSig1)/s);
231  hNDiffCntSig1->SetBinError(iBin+1,TMath::Sqrt(cntErr)/s);
232  hCntSig2->SetBinContent(iBin+1,cntSig2);
233  hNDiffCntSig2->SetBinContent(iBin+1,(s-cntSig2)/s);
234  hNDiffCntSig2->SetBinError(iBin+1,TMath::Sqrt(cntErr)/s);
235  hCntSig2->SetBinError(iBin+1,TMath::Sqrt(cntErr));
236  hSignal->SetBinContent(iBin+1,ry);
237  hSignal->SetBinError(iBin+1,ery);
238  hRelErrSig->SetBinContent(iBin+1,errs/s);
239  hInvSignif->SetBinContent(iBin+1,1/sig);
240  hInvSignif->SetBinError(iBin+1,errsig/(sig*sig));
241  hBackground->SetBinContent(iBin+1,b); //consider sigma
242  hBackground->SetBinError(iBin+1,errb);
243  hBackgroundNormSigma->SetBinContent(iBin+1,b/(3*fitter[iBin]->GetSigma())*(3*0.012)); //consider sigma
244  hBackgroundNormSigma->SetBinError(iBin+1,errb);
245  hSignificance->SetBinContent(iBin+1,sig);
246  hSignificance->SetBinError(iBin+1,errsig);
247  hMass->SetBinContent(iBin+1,mass);
248  hMass->SetBinError(iBin+1,0.0001);
249  hSigma->SetBinContent(iBin+1,sigma);
250  hSigma->SetBinError(iBin+1,fitter[iBin]->GetSigmaUncertainty());
251 
252  }
253 
254  /*
255  c1->cd(1); // is some cases the fitting function of 1st bin get lost
256  funBckStore1->Draw("same");
257  funBckStore2->Draw("same");
258  funBckStore3->Draw("same");
259  */
260 
261  TCanvas *cpar=new TCanvas("cpar","Fit params",1200,600);
262  cpar->Divide(2,1);
263  cpar->cd(1);
264  hMass->SetMarkerStyle(20);
265  hMass->Draw("PE");
266  hMass->GetXaxis()->SetTitle("Pt (GeV/c)");
267  hMass->GetXaxis()->SetTitle("Mass (GeV/c^{2})");
268  cpar->cd(2);
269  hSigma->SetMarkerStyle(20);
270  hSigma->Draw("PE");
271  hSigma->GetXaxis()->SetTitle("Pt (GeV/c)");
272  hSigma->GetXaxis()->SetTitle("Sigma (GeV/c^{2})");
273 
274  TCanvas* csig=new TCanvas("csig","Results",1200,600);
275  csig->Divide(3,1);
276  csig->cd(1);
277  hSignal->SetMarkerStyle(20);
278  hSignal->SetMarkerColor(4);
279  hSignal->SetLineColor(4);
280  hSignal->GetXaxis()->SetTitle("Pt (GeV/c)");
281  hSignal->GetYaxis()->SetTitle("Signal");
282  hSignal->Draw("P");
283  hCntSig1->SetMarkerStyle(26);
284  hCntSig1->SetMarkerColor(2);
285  hCntSig1->SetLineColor(2);
286  hCntSig1->Draw("PSAME");
287  hCntSig2->SetMarkerStyle(29);
288  hCntSig2->SetMarkerColor(kGray+1);
289  hCntSig2->SetLineColor(kGray+1);
290  hCntSig2->Draw("PSAME");
291  TLegend* leg=new TLegend(0.4,0.7,0.89,0.89);
292  leg->SetFillColor(0);
293  TLegendEntry* ent=leg->AddEntry(hSignal,"From Fit","PL");
294  ent->SetTextColor(hSignal->GetMarkerColor());
295  ent=leg->AddEntry(hCntSig1,"From Counting1","PL");
296  ent->SetTextColor(hCntSig1->GetMarkerColor());
297  ent=leg->AddEntry(hCntSig2,"From Counting2","PL");
298  ent->SetTextColor(hCntSig2->GetMarkerColor());
299  leg->Draw();
300  csig->cd(2);
301  hBackground->SetMarkerStyle(20);
302  hBackground->Draw("P");
303  hBackground->GetXaxis()->SetTitle("Pt (GeV/c)");
304  hBackground->GetYaxis()->SetTitle("Background");
305  csig->cd(3);
306  hSignificance->SetMarkerStyle(20);
307  hSignificance->Draw("P");
308  hSignificance->GetXaxis()->SetTitle("Pt (GeV/c)");
309  hSignificance->GetYaxis()->SetTitle("Significance");
310 
311  TCanvas* cDiffS=new TCanvas("cDiffS","Difference",1200,600);
312  cDiffS->Divide(2,1);
313  cDiffS->cd(1);
314  hRelErrSig->SetMarkerStyle(20); //fullcircle
315  hRelErrSig->SetTitleOffset(1.2);
316  hRelErrSig->SetTitle("Rel Error from Fit;p_{t} (GeV/c);Signal Relative Error");
317  hRelErrSig->Draw("P");
318  hInvSignif->SetMarkerStyle(21); //fullsquare
319  hInvSignif->SetMarkerColor(kMagenta+1);
320  hInvSignif->SetLineColor(kMagenta+1);
321  hInvSignif->Draw("PSAMES");
322  TLegend* leg2=new TLegend(0.4,0.7,0.89,0.89);
323  leg2->SetFillColor(0);
324  TLegendEntry* ent2=leg2->AddEntry(hRelErrSig,"From Fit","P");
325  ent2->SetTextColor(hRelErrSig->GetMarkerColor());
326  ent2=leg2->AddEntry(hInvSignif,"1/Significance","PL");
327  ent2->SetTextColor(hInvSignif->GetMarkerColor());
328  leg2->Draw();
329 
330  cDiffS->cd(2);
331  hNDiffCntSig1->SetMarkerStyle(26);
332  hNDiffCntSig1->SetMarkerColor(2);
333  hNDiffCntSig1->SetLineColor(2);
334  hNDiffCntSig1->SetTitle("Cmp Fit-Count;p_{t} (GeV/c);(S_{fit}-S_{count})/S_{fit}");
335  hNDiffCntSig1->Draw("P");
336  hNDiffCntSig2->SetMarkerStyle(29);
337  hNDiffCntSig2->SetMarkerColor(kGray+1);
338  hNDiffCntSig2->SetLineColor(kGray+1);
339  hNDiffCntSig2->Draw("PSAME");
340  TLegend* leg1=new TLegend(0.4,0.7,0.89,0.89);
341  leg1->SetFillColor(0);
342  TLegendEntry* ent1=leg1->AddEntry(hNDiffCntSig1,"From Counting1","PL");
343  ent1->SetTextColor(hNDiffCntSig1->GetMarkerColor());
344  ent1=leg1->AddEntry(hNDiffCntSig2,"From Counting2","PL");
345  ent1->SetTextColor(hNDiffCntSig2->GetMarkerColor());
346  leg1->Draw();
347 
348  TGraph* grReducedChiSquare=new TGraph(nPtBins,ptlims,arrchisquare);
349  grReducedChiSquare->SetName("grReducedChiSquare");
350  grReducedChiSquare->SetTitle("Reduced Chi2;p_{t} (GeV/c);#tilde{#chi}^{2}");
351  TCanvas *cChi2=new TCanvas("cChi2","reduced chi square",600,600);
352  cChi2->cd();
353  grReducedChiSquare->SetMarkerStyle(21);
354  grReducedChiSquare->Draw("AP");
355 
356  TCanvas* cbkgNormSigma=new TCanvas("cbkgNormSigma","Background normalized to sigma",400,600);
357  cbkgNormSigma->cd();
358  hBackgroundNormSigma->SetMarkerStyle(20);
359  hBackgroundNormSigma->Draw("P");
360  hBackgroundNormSigma->GetXaxis()->SetTitle("Pt (GeV/c)");
361  hBackgroundNormSigma->GetYaxis()->SetTitle("Background #times 3 #times 0.012/ (3 #times #sigma)");
362  hBackgroundNormSigma->Draw();
363 
364 
365  TString partname="Both";
367  if(analysisType==kD0toKpi) partname="D0";
368  if(analysisType==kDplusKpipi) partname="Dplus";
369  if(analysisType==kDsKKpi) partname="Dsplus";
370  }
372  if(analysisType==kD0toKpi) partname="D0bar";
373  if(analysisType==kDplusKpipi) partname="Dminus";
374  if(analysisType==kDsKKpi) partname="Dsminus";
375  }
376 
377  printf("Events for norm = %f\n",nEventsForNorm);
378  TH1F* hNEvents=new TH1F("hNEvents","",1,0.,1.);
379  hNEvents->SetBinContent(1,nEventsForNorm);
380 
381  TFile* outf=new TFile(Form("RawYield%s.root",partname.Data()),"update");
382  outf->cd();
383  hNEvents->Write();
384  hMass->Write();
385  hSigma->Write();
386  hCntSig1->Write();
387  hCntSig2->Write();
388  hNDiffCntSig1->Write();
389  hNDiffCntSig2->Write();
390  hSignal->Write();
391  hRelErrSig->Write();
392  hInvSignif->Write();
393  hBackground->Write();
394  hBackgroundNormSigma->Write();
395  hSignificance->Write();
396  grReducedChiSquare->Write();
397  hPtMass->Write();
398  outf->Close();
399 }
400 
401 
402 Bool_t LoadDplusHistos(TObjArray* listFiles, TH1F** hMass){
403 
404  Int_t nFiles=listFiles->GetEntries();
405  TList **hlist=new TList*[nFiles];
407 
408  Int_t nReadFiles=0;
409  for(Int_t iFile=0; iFile<nFiles; iFile++){
410  TString fName=((TObjString*)listFiles->At(iFile))->GetString();
411  TFile *f=TFile::Open(fName.Data());
412  if(!f){
413  printf("ERROR: file %s does not exist\n",fName.Data());
414  continue;
415  }
416  printf("Open File %s\n",f->GetName());
417  TDirectory *dir = (TDirectory*)f->Get("PWG3_D2H_InvMassDplus");
418  if(!dir){
419  printf("ERROR: directory PWG3_D2H_InvMassDplus not found in %s\n",fName.Data());
420  continue;
421  }
422  hlist[nReadFiles]=(TList*)dir->Get(Form("coutputDplus%s",suffix.Data()));
423  TList *listcut = (TList*)dir->Get(Form("coutputDplusCuts%s",suffix.Data()));
424  dcuts[nReadFiles]=(AliRDHFCutsDplustoKpipi*)listcut->At(0);
425  if(nReadFiles>0){
426  Bool_t sameCuts=dcuts[nReadFiles]->CompareCuts(dcuts[0]);
427  if(!sameCuts){
428  printf("ERROR: Cut objects do not match\n");
429  return kFALSE;
430  }
431  }
432  AliNormalizationCounter* c=(AliNormalizationCounter*)dir->Get(Form("coutputDplusNorm%s",suffix.Data()));
433  printf("Events for normalization = %f\n",c->GetNEventsForNorm());
435  nReadFiles++;
436  }
437  if(nReadFiles<nFiles){
438  printf("WARNING: not all requested files have been found\n");
439  }
440 
441  Int_t nPtBinsCuts=dcuts[0]->GetNPtBins();
442  printf("Number of pt bins for cut object = %d\n",nPtBins);
443  Float_t *ptlimsCuts=dcuts[0]->GetPtBinLimits();
444  ptlimsCuts[nPtBinsCuts]=ptlimsCuts[nPtBinsCuts-1]+4.;
445 
446  Int_t iFinBin=0;
447  for(Int_t i=0;i<nPtBinsCuts;i++){
448  if(ptlimsCuts[i]>=ptlims[iFinBin+1]) iFinBin+=1;
449  if(iFinBin>nPtBins) break;
450  if(ptlimsCuts[i]>=ptlims[iFinBin] &&
451  ptlimsCuts[i+1]<=ptlims[iFinBin+1]){
452  for(Int_t iFile=0; iFile<nReadFiles; iFile++){
453  TString histoName;
454  if(optPartAntiPart==kBoth) histoName.Form("hMassPt%dTC",i);
455  else if(optPartAntiPart==kParticleOnly) histoName.Form("hMassPt%dTCPlus",i);
456  else if(optPartAntiPart==kAntiParticleOnly) histoName.Form("hMassPt%dTCMinus",i);
457  TH1F* htemp=(TH1F*)hlist[iFile]->FindObject(histoName.Data());
458  if(!htemp){
459  printf("ERROR: Histogram %s not found\n",histoName.Data());
460  return kFALSE;
461  }
462  if(!hMass[iFinBin]){
463  hMass[iFinBin]=new TH1F(*htemp);
464  }else{
465  hMass[iFinBin]->Add(htemp);
466  }
467  }
468  }
469  }
470  TString partname="Both";
471  if(optPartAntiPart==kParticleOnly) partname="Dplus";
472  if(optPartAntiPart==kAntiParticleOnly) partname="Dminus";
473 
474  for(Int_t iFile=0; iFile<nReadFiles; iFile++){
475  TH2F* htemp2=(TH2F*)hlist[iFile]->FindObject("hPtVsMassTC");
476  if(iFile==0){
477  hPtMass=new TH2F(*htemp2);
478  }else{
479  hPtMass->Add(htemp2);
480  }
481  }
482 
483  TFile* outf=new TFile(Form("RawYield%s.root",partname.Data()),"recreate");
484  outf->cd();
485  dcuts[0]->Write();
486  outf->Close();
487 
488  return kTRUE;
489 
490 }
491 
492 
493 Bool_t LoadDsHistos(TObjArray* listFiles, TH1F** hMass){
494 
495  Int_t nFiles=listFiles->GetEntries();
496  TList **hlist=new TList*[nFiles];
497  AliRDHFCutsDstoKKpi** dcuts=new AliRDHFCutsDstoKKpi*[nFiles];
498 
499  Int_t nReadFiles=0;
500  for(Int_t iFile=0; iFile<nFiles; iFile++){
501  TString fName=((TObjString*)listFiles->At(iFile))->GetString();
502  TFile *f=TFile::Open(fName.Data());
503  if(!f){
504  printf("ERROR: file %s does not exist\n",fName.Data());
505  continue;
506  }
507  printf("Open File %s\n",f->GetName());
508  TDirectory *dir = (TDirectory*)f->Get("PWG3_D2H_InvMassDs");
509  if(!dir){
510  printf("ERROR: directory PWG3_D2H_InvMassDs not found in %s\n",fName.Data());
511  continue;
512  }
513  hlist[nReadFiles]=(TList*)dir->Get("coutputDs");
514  TList *listcut = (TList*)dir->Get("coutputDsCuts");
515  dcuts[nReadFiles]=(AliRDHFCutsDstoKKpi*)listcut->At(0);
516  cout<< dcuts[nReadFiles]<<endl;
517  if(nReadFiles>0){
518  Bool_t sameCuts=dcuts[nReadFiles]->CompareCuts(dcuts[0]);
519  if(!sameCuts){
520  printf("ERROR: Cut objects do not match\n");
521  return kFALSE;
522  }
523  }
524  nReadFiles++;
525  }
526  if(nReadFiles<nFiles){
527  printf("WARNING: not all requested files have been found\n");
528  }
529 
530  Int_t nPtBinsCuts=dcuts[0]->GetNPtBins();
531  printf("Number of pt bins for cut object = %d\n",nPtBins);
532  Float_t *ptlimsCuts=dcuts[0]->GetPtBinLimits();
533  ptlimsCuts[nPtBinsCuts]=ptlimsCuts[nPtBinsCuts-1]+4.;
534 
535  Int_t iFinBin=0;
536  for(Int_t i=0;i<nPtBinsCuts;i++){
537  if(ptlimsCuts[i]>=ptlims[iFinBin+1]) iFinBin+=1;
538  if(iFinBin>nPtBins) break;
539  if(ptlimsCuts[i]>=ptlims[iFinBin] &&
540  ptlimsCuts[i+1]<=ptlims[iFinBin+1]){
541  for(Int_t iFile=0; iFile<nReadFiles; iFile++){
542  TString histoName;
543  if(optPartAntiPart==kBoth) histoName.Form("hMassAllPt%dphi",i);
544  else if(optPartAntiPart==kParticleOnly){
545  printf("Particle/Antiparticle not yet enabled for Ds");
546  histoName.Form("hMassAllPt%dphi",i);
547  }
549  printf("Particle/Antiparticle not yet enabled for Ds");
550  histoName.Form("hMassAllPt%dphi",i);
551  }
552  TH1F* htemp=(TH1F*)hlist[iFile]->FindObject(histoName.Data());
553  if(!htemp){
554  printf("ERROR: Histogram %s not found\n",histoName.Data());
555  return kFALSE;
556  }
557  if(!hMass[iFinBin]){
558  hMass[iFinBin]=new TH1F(*htemp);
559  }else{
560  hMass[iFinBin]->Add(htemp);
561  }
562  }
563  }
564  }
565  TString partname="Both";
566  if(optPartAntiPart==kParticleOnly) partname="Both";
567  if(optPartAntiPart==kAntiParticleOnly) partname="Both";
568 
569  for(Int_t iFile=0; iFile<nReadFiles; iFile++){
570  TH2F* htemp2=(TH2F*)hlist[iFile]->FindObject("hPtVsMassPhi");
571  if(iFile==0){
572  hPtMass=new TH2F(*htemp2);
573  }else{
574  hPtMass->Add(htemp2);
575  }
576  }
577 
578  TFile* outf=new TFile(Form("RawYield%s.root",partname.Data()),"recreate");
579  outf->cd();
580  dcuts[0]->Write();
581  outf->Close();
582 
583  return kTRUE;
584 
585 }
586 
587 
588 Bool_t LoadD0toKpiHistos(TObjArray* listFiles, TH1F** hMass){
589  //
590  Int_t nFiles=listFiles->GetEntries();
591  TList **hlist=new TList*[nFiles];
592  AliRDHFCutsD0toKpi** dcuts=new AliRDHFCutsD0toKpi*[nFiles];
593 
594  Int_t nReadFiles=0;
595  for(Int_t iFile=0; iFile<nFiles; iFile++){
596  TString fName=((TObjString*)listFiles->At(iFile))->GetString();
597  TFile *f=TFile::Open(fName.Data());
598  if(!f){
599  printf("ERROR: file %s does not exist\n",fName.Data());
600  continue;
601  }
602  printf("Open File %s\n",f->GetName());
603 
604  TString dirname="PWG3_D2H_D0InvMass";
605  if(optPartAntiPart==kParticleOnly) dirname+="D0";
606  else if(optPartAntiPart==kAntiParticleOnly) dirname+="D0bar";
607  if(cutsappliedondistr) dirname+="C";
608  TDirectory *dir = (TDirectory*)f->Get(dirname);
609  if(!dir){
610  printf("ERROR: directory %s not found in %s\n",dirname.Data(),fName.Data());
611  continue;
612  }
613  TString listmassname="coutputmassD0Mass";
614  if(optPartAntiPart==kParticleOnly) listmassname+="D0";
615  else if(optPartAntiPart==kAntiParticleOnly) listmassname+="D0bar";
616  if(cutsappliedondistr) listmassname+="C";
617 
618  hlist[nReadFiles]=(TList*)dir->Get(listmassname);
619 
620  TString cutsobjname="cutsD0";
621  if(optPartAntiPart==kParticleOnly) cutsobjname+="D0";
622  else if(optPartAntiPart==kAntiParticleOnly) cutsobjname+="D0bar";
623  if(cutsappliedondistr) cutsobjname+="C";
624 
625  dcuts[nReadFiles]=(AliRDHFCutsD0toKpi*)dir->Get(cutsobjname);
626  if(nReadFiles>0){
627  Bool_t sameCuts=dcuts[nReadFiles]->CompareCuts(dcuts[0]);
628  if(!sameCuts){
629  printf("ERROR: Cut objects do not match\n");
630  return kFALSE;
631  }
632  }
633  nReadFiles++;
634  }
635  if(nReadFiles<nFiles){
636  printf("WARNING: not all requested files have been found\n");
637  if (nReadFiles==0) {
638  printf("ERROR: Any file/dir found\n");
639  return kFALSE;
640  }
641  }
642 
643  Int_t nPtBinsCuts=dcuts[0]->GetNPtBins();
644  printf("Number of pt bins for cut object = %d\n",nPtBins);
645  Float_t *ptlimsCuts=dcuts[0]->GetPtBinLimits();
646  ptlimsCuts[nPtBinsCuts]=ptlimsCuts[nPtBinsCuts-1]+4.;
647 
648  Int_t iFinBin=0;
649  for(Int_t i=0;i<nPtBinsCuts;i++){
650  if(ptlimsCuts[i]>=ptlims[iFinBin+1]) iFinBin+=1;
651  if(iFinBin>nPtBins) break;
652  if(ptlimsCuts[i]>=ptlims[iFinBin] &&
653  ptlimsCuts[i+1]<=ptlims[iFinBin+1]){
654  for(Int_t iFile=0; iFile<nReadFiles; iFile++){
655  TH1F* htemp=(TH1F*)hlist[iFile]->FindObject(Form("histMass_%d",i));
656  if(!hMass[iFinBin]){
657  hMass[iFinBin]=new TH1F(*htemp);
658  }else{
659  hMass[iFinBin]->Add(htemp);
660  }
661  }
662  }
663  }
664  TString partname="Both";
665  if(optPartAntiPart==kParticleOnly) partname="D0";
666  if(optPartAntiPart==kAntiParticleOnly) partname="D0bar";
667 
668  TFile* outf=new TFile(Form("RawYield%s.root",partname.Data()),"recreate");
669  outf->cd();
670  dcuts[0]->Write();
671  outf->Close();
672  return kTRUE;
673 }
674 
675 Bool_t LoadDstarD0piHistos(TObjArray* listFiles, TH1F** hMass){
676  //
677  Int_t nFiles=listFiles->GetEntries();
678  TList **hlist=new TList*[nFiles];
680  Int_t nReadFiles=0;
681  for(Int_t iFile=0; iFile<nFiles; iFile++){
682  TString fName=((TObjString*)listFiles->At(iFile))->GetString();
683  TFile *f=TFile::Open(fName.Data());
684  if(!f){
685  printf("ERROR: file %s does not exist\n",fName.Data());
686  continue;
687  }
688  printf("Open File %s\n",f->GetName());
689  TString dirname="PWG3_D2H_DStarSpectra";
690  TDirectory *dir = (TDirectory*)f->Get(dirname);
691  if(!dir){
692  printf("ERROR: directory %s not found in %s\n",dirname.Data(),fName.Data());
693  continue;
694  }
695  TString listmassname="DStarPID00";
696 
697  hlist[nReadFiles]=(TList*)dir->Get(listmassname);
698  TString cutsobjname="DStartoKpipiCuts";
699  dcuts[nReadFiles]=(AliRDHFCutsDStartoKpipi*)dir->Get(cutsobjname);
700  if(nReadFiles>0){
701  Bool_t sameCuts=dcuts[nReadFiles]->CompareCuts(dcuts[0]);
702  if(!sameCuts){
703  printf("ERROR: Cut objects do not match\n");
704  return kFALSE;
705  }
706  }
707  nReadFiles++;
708  }
709  if(nReadFiles<nFiles){
710  printf("WARNING: not all requested files have been found\n");
711  if (nReadFiles==0) {
712  printf("ERROR: Any file/dir found\n");
713  return kFALSE;
714  }
715  }
716  Int_t nPtBinsCuts=dcuts[0]->GetNPtBins();
717  printf("Number of pt bins for cut object = %d\n",nPtBins);
718  Float_t *ptlimsCuts=dcuts[0]->GetPtBinLimits();
719  ptlimsCuts[nPtBinsCuts]=ptlimsCuts[nPtBinsCuts-1]+4.;
720 
721 
722 
723 
724  Int_t iFinBin=0;
725  for(Int_t i=0;i<nPtBinsCuts;i++){
726  if(ptlimsCuts[i]>=ptlims[iFinBin+1]) iFinBin+=1;
727  if(iFinBin>nPtBins) break;
728  if(ptlimsCuts[i]>=ptlims[iFinBin] &&
729  ptlimsCuts[i+1]<=ptlims[iFinBin+1]){
730  for(Int_t iFile=0; iFile<nReadFiles; iFile++){
731  TH1F* htemp=(TH1F*)hlist[iFile]->FindObject(Form("histDeltaMass_%d",i));
732  if(!hMass[iFinBin]){
733  hMass[iFinBin]=new TH1F(*htemp);
734  }else{
735  hMass[iFinBin]->Add(htemp);
736  }
737  }
738  }
739  }
740  TString partname="Both";
741  if(optPartAntiPart==kParticleOnly) partname="DStar";
742  if(optPartAntiPart==kAntiParticleOnly) partname="DStarbar";
743 
744  TFile* outf=new TFile(Form("RawYield%s.root",partname.Data()),"recreate");
745  outf->cd();
746  dcuts[0]->Write();
747  outf->Close();
748  return kTRUE;
749 }
750 
751 void CompareFitTypes(TString* paths, TString* legtext,Int_t ncmp=3,TString* filenameYield=0x0){
752  //read ncmp RawYield.roots and draw them together
753  //set the global variable nevents before running
754  //arguments:
755  // - paths= vector of ncmp dimension with the paths of the file RawYield.root
756  // - legtext= vector of ncmp dimension with the label for the legend
757  // - ncmp= number of files to compare (default is 3)
758  // - filenameYield= optional ncmp-dimensional array with the difference between the names of the files to be compared (useful if the 2 files are in the same directory but have different names)
759 
760  gStyle->SetOptStat(0);
761  gStyle->SetOptTitle(0);
762  gStyle->SetCanvasColor(0);
763  gStyle->SetFrameFillColor(0);
764  gStyle->SetTitleFillColor(0);
765  gStyle->SetFrameBorderMode(0);
766 
767  if(!filenameYield) filenameYield=new TString[ncmp];
768 
769  for(Int_t k=0;k<ncmp;k++){
770  if(!filenameYield) filenameYield[k]="RawYield.root";
771  filenameYield[k].Prepend(paths[k]);
772  }
773 
774  TCanvas* cSig=new TCanvas("cSig","Results",1200,600);
775  cSig->Divide(3,1);
776  TCanvas* cBkgN=new TCanvas("cBkgN","Background normalized to sigma",400,600);
777  TCanvas* cDiffS=new TCanvas("cDiffS","Difference",1200,600);
778  cDiffS->Divide(2,1);
779  TCanvas *cChi2=new TCanvas("cChi2","reduced chi square",600,600);
780 
781  TLegend* leg1=new TLegend(0.4,0.7,0.89,0.89);
782  leg1->SetFillColor(0);
783  TLegend* leg2=(TLegend*)leg1->Clone();
784  TLegend* leg3=(TLegend*)leg1->Clone();
785  TLegend* leg4=new TLegend(0.4,0.6,0.8,0.89);
786  leg4->SetFillColor(0);
787 
788  for(Int_t iTypes=0;iTypes<ncmp;iTypes++){
789  TFile* fin=new TFile(filenameYield[iTypes]);
790  if(!fin){
791  printf("WARNING: %s not found",filenameYield[iTypes].Data());
792  continue;
793  }
794 
795  TH1F* hSignal=(TH1F*)fin->Get("hSignal");
796  TH1F* hBackground=(TH1F*)fin->Get("hBackground");
797  TH1F* hBackgroundNormSigma=(TH1F*)fin->Get("hBackgroundNormSigma");
798  TH1F* hSignificance=(TH1F*)fin->Get("hSignificance");
799  hSignal->SetName(Form("%s%d",hSignal->GetName(),iTypes));
800  hBackground->SetName(Form("%s%d",hBackground->GetName(),iTypes));
801  hBackgroundNormSigma->SetName(Form("%s%d",hBackgroundNormSigma->GetName(),iTypes));
802  hSignificance->SetName(Form("%s%d",hSignificance->GetName(),iTypes));
803 
804  hSignal->SetMarkerColor(iTypes+2);
805  hSignal->SetLineColor(iTypes+2);
806  hBackground->SetMarkerColor(iTypes+2);
807  hBackground->SetLineColor(iTypes+2);
808  hBackgroundNormSigma->SetMarkerColor(iTypes+2);
809  hBackgroundNormSigma->SetLineColor(iTypes+2);
810  hSignificance->SetMarkerColor(iTypes+2);
811  hSignificance->SetLineColor(iTypes+2);
812 
813  TLegendEntry* ent4=leg4->AddEntry(hSignal,Form("%s",legtext[iTypes].Data()),"PL");
814  ent4->SetTextColor(hSignal->GetMarkerColor());
815 
816  if(ncmp==nsamples){
817  printf("Info: Normalizing signal, background and significance to the number of events\n");
818  hSignal->Scale(1./nevents[iTypes]);
819  hBackground->Scale(1./nevents[iTypes]);
820  hBackgroundNormSigma->Scale(1./nevents[iTypes]);
821  hSignificance->Scale(1./TMath::Sqrt(nevents[iTypes]));
822  }
823 
824  if(iTypes==0){
825  cSig->cd(1);
826  hSignal->DrawClone("P");
827  cSig->cd(2);
828  hBackground->DrawClone("P");
829  cSig->cd(3);
830  hSignificance->DrawClone("P");
831  cBkgN->cd();
832  hBackgroundNormSigma->DrawClone("P");
833  } else{
834  cSig->cd(1);
835  hSignal->DrawClone("Psames");
836  cSig->cd(2);
837  hBackground->DrawClone("Psames");
838  cSig->cd(3);
839  hSignificance->DrawClone("Psames");
840  cBkgN->cd();
841  hBackgroundNormSigma->DrawClone("Psames");
842  }
843 
844  TH1F* hRelErrSig=(TH1F*)fin->Get("hRelErrSig");
845  TH1F* hInvSignif=(TH1F*)fin->Get("hInvSignif");
846  hRelErrSig->SetName(Form("%s%d",hRelErrSig->GetName(),iTypes));
847  hInvSignif->SetName(Form("%s%d",hInvSignif->GetName(),iTypes));
848 
849  hRelErrSig->SetMarkerColor(iTypes+2);
850  hRelErrSig->SetLineColor(iTypes+2);
851  hInvSignif->SetMarkerColor(iTypes+2);
852  hInvSignif->SetLineColor(iTypes+2);
853 
854  TLegendEntry* ent1=leg1->AddEntry(hRelErrSig,Form("From Fit (%s)",legtext[iTypes].Data()),"P");
855  ent1->SetTextColor(hRelErrSig->GetMarkerColor());
856  ent1=leg1->AddEntry(hInvSignif,Form("1/Significance (%s)",legtext[iTypes].Data()),"PL");
857  ent1->SetTextColor(hInvSignif->GetMarkerColor());
858 
859  cDiffS->cd(1);
860  if(iTypes==0){
861  hRelErrSig->DrawClone("P");
862  hInvSignif->DrawClone();
863  } else{
864  hRelErrSig->DrawClone("Psames");
865  hInvSignif->DrawClone("sames");
866  }
867 
868  TH1F* hNDiffCntSig1=(TH1F*)fin->Get("hNDiffCntSig1");
869  TH1F* hNDiffCntSig2=(TH1F*)fin->Get("hNDiffCntSig2");
870  hNDiffCntSig1->SetName(Form("%s%d",hNDiffCntSig1->GetName(),iTypes));
871  hNDiffCntSig2->SetName(Form("%s%d",hNDiffCntSig2->GetName(),iTypes));
872 
873  hNDiffCntSig1->SetMarkerColor(iTypes+2);
874  hNDiffCntSig1->SetLineColor(iTypes+2);
875  hNDiffCntSig2->SetMarkerColor(iTypes+2);
876  hNDiffCntSig2->SetLineColor(iTypes+2);
877  TLegendEntry* ent2=leg2->AddEntry(hNDiffCntSig1,Form("From Counting1 (%s)",legtext[iTypes].Data()),"PL");
878  ent2->SetTextColor(hNDiffCntSig1->GetMarkerColor());
879  ent2=leg2->AddEntry(hNDiffCntSig2,Form("From Counting2 (%s)",legtext[iTypes].Data()),"PL");
880  ent2->SetTextColor(hNDiffCntSig2->GetMarkerColor());
881 
882  cDiffS->cd(2);
883  if(iTypes==0){
884  hNDiffCntSig1->DrawClone();
885  hNDiffCntSig2->DrawClone();
886  }else{
887  hNDiffCntSig1->DrawClone("sames");
888  hNDiffCntSig2->DrawClone("sames");
889  }
890 
891  TGraph* grReducedChiSquare=(TGraph*)fin->Get("grReducedChiSquare");
892  grReducedChiSquare->SetName(Form("%s%d",grReducedChiSquare->GetName(),iTypes));
893 
894  grReducedChiSquare->SetMarkerColor(iTypes+2);
895  grReducedChiSquare->SetLineColor(iTypes+2);
896  TLegendEntry* ent3=leg3->AddEntry(grReducedChiSquare,Form("%s",legtext[iTypes].Data()),"PL");
897  ent3->SetTextColor(grReducedChiSquare->GetMarkerColor());
898 
899  cChi2->cd();
900  if(iTypes==0) grReducedChiSquare->DrawClone("AP");
901  else grReducedChiSquare->DrawClone("P");
902  }
903 
904  cSig->cd(1);
905  leg4->Draw();
906 
907  cDiffS->cd(1);
908  leg1->Draw();
909 
910  cDiffS->cd(2);
911  leg2->Draw();
912 
913  cChi2->cd();
914  leg3->Draw();
915 
916  TFile* fout=new TFile("ComparisonRawYield.root","RECREATE");
917  fout->cd();
918  cDiffS->Write();
919  cChi2->Write();
920  fout->Close();
921 }
922 
923 
TH1D * hSignal
TH1F * GetHistoClone() const
void DrawHere(TVirtualPad *pd, Double_t nsigma=3, Int_t writeFitInfo=1) const
write the canvas in a root file
const Int_t nsamples
static TH1D * RebinHisto(TH1 *hOrig, Int_t reb, Int_t firstUse=-1)
Rebinning of invariant mass histograms.
double Double_t
Definition: External.C:58
Definition: External.C:236
void FitMassSpectra(Int_t analysisType=kDplusKpipi, TString fileNameb="", TString fileNamec="", TString fileNamed="", TString fileNamee="")
TH1D * hSigma
Bool_t LoadDstarD0piHistos(TObjArray *listFiles, TH1F **hMass)
Bool_t LoadDplusHistos(TObjArray *listFiles, TH1F **hMass)
TString suffix
void SetInitialGaussianMean(Double_t mean)
Double_t mass
Bool_t LoadD0toKpiHistos(TObjArray *listFiles, TH1F **hMass)
Double_t ptlims[nPtBins+1]
void CompareFitTypes(TString *paths, TString *legtext, Int_t ncmp=3, TString *filenameYield=0x0)
TCanvas * c
Definition: TestFitELoss.C:172
Int_t types
Int_t firstUsedBin[nPtBins]
TH2F * hPtMass
TF1 * GetBackgroundFullRangeFunc()
TF1 * GetBackgroundRecalcFunc()
Int_t optPartAntiPart
Double_t GetRawYieldError() const
Int_t typeb
Double_t * sigma
virtual void Signal(Double_t nOfSigma, Double_t &signal, Double_t &errsignal) const
return total integral of the histogram
Class for cuts on AOD reconstructed D+->Kpipi.
Double_t massD
int Int_t
Definition: External.C:63
Bool_t LoadDsHistos(TObjArray *listFiles, TH1F **hMass)
float Float_t
Definition: External.C:68
Double_t GetMean() const
Double_t GetReducedChiSquare() const
Int_t factor4refl
void SetReflectionSigmaFactor(Int_t constant)
Int_t GetBinN() const
Definition: External.C:212
virtual void Background(Double_t nOfSigma, Double_t &background, Double_t &errbackground) const
signal in (min, max) with error
Double_t GetRawYield() const
TList * fitter
Definition: DrawAnaELoss.C:26
void SetInitialGaussianSigma(Double_t sigma)
change the default value of the mean
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
TObject * FindObject(int bin, const char *nameH, const TList *lst, Bool_t normPerEvent=kTRUE)
const Int_t nPtBins
Bool_t CompareCuts(const AliRDHFCuts *obj) const
virtual Bool_t MassFitter(Bool_t draw=kTRUE)
Bool_t cutsappliedondistr
Float_t * GetPtBinLimits() const
Definition: AliRDHFCuts.h:248
Int_t nevents[nsamples]
Int_t GetNPtBins() const
Definition: AliRDHFCuts.h:249
Double_t GetSigma() const
bool Bool_t
Definition: External.C:53
TFile * fout
input train file
Float_t massRangeForCounting
void Significance(Double_t nOfSigma, Double_t &significance, Double_t &errsignificance) const
backgournd in (min, max) with error
Double_t nEventsForNorm
AliHFMassFitter for the fit of invariant mass distribution of charmed mesons.
Int_t rebin[nPtBins]
Double_t minMassForFit
Double_t maxMassForFit
TDirectoryFile * dir