AliPhysics  8d00e07 (8d00e07)
AliHFInvMassMultiTrialFit.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 2008-2019, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 #include <TMath.h>
17 #include <TPad.h>
18 #include <TCanvas.h>
19 #include <TH1F.h>
20 #include <TH2F.h>
21 #include <TNtuple.h>
22 #include <TF1.h>
23 #include <TLatex.h>
24 #include <TFile.h>
25 #include "AliHFInvMassFitter.h"
27 #include "AliVertexingHFUtils.h"
28 
32 
33 
34 //_________________________________________________________________________
36  TNamed(),
37  fNumOfRebinSteps(4),
38  fRebinSteps(0x0),
39  fNumOfFirstBinSteps(1),
40  fNumOfLowLimFitSteps(6),
41  fLowLimFitSteps(0x0),
42  fNumOfUpLimFitSteps(6),
43  fUpLimFitSteps(0x0),
44  fNumOfnSigmaBinCSteps(11),
45  fnSigmaBinCSteps(0x0),
46  fnSigmaForBkgEval(3),
47  fSigmaGausMC(0.010),
48  fSigmaMCVariation(0.15),
49  fMassD(1.86484),
50  fSuffix(""),
51  fFitOption(0),
52  fUseExpoBkg(kTRUE),
53  fUseLinBkg(kTRUE),
54  fUsePol2Bkg(kTRUE),
55  fUsePol3Bkg(kTRUE),
56  fUsePol4Bkg(kTRUE),
57  fUsePol5Bkg(kFALSE),
58  fUsePowLawBkg(kFALSE),
59  fUsePowLawTimesExpoBkg(kFALSE),
60  fUseFixSigUpFreeMean(kTRUE),
61  fUseFixSigDownFreeMean(kTRUE),
62  fUseFreeS(kTRUE),
63  fUseFixedMeanFreeS(kTRUE),
64  fUseFixSigFreeMean(kTRUE),
65  fUseFixSigFixMean(kTRUE),
66  fUseSecondPeak(kFALSE),
67  fMassSecondPeak(1.86958),
68  fSigmaSecondPeak(0.01),
69  fFixMassSecondPeak(kFALSE),
70  fFixSigmaSecondPeak(kFALSE),
71  fSaveBkgVal(kFALSE),
72  fDrawIndividualFits(kFALSE),
73  fHistoRawYieldDistAll(0x0),
74  fHistoRawYieldTrialAll(0x0),
75  fHistoSigmaTrialAll(0x0),
76  fHistoMeanTrialAll(0x0),
77  fHistoChi2TrialAll(0x0),
78  fHistoSignifTrialAll(0x0),
79  fHistoBkgTrialAll(0x0),
80  fHistoBkgInBinEdgesTrialAll(0x0),
81  fHistoRawYieldDistBinC0All(0x0),
82  fHistoRawYieldTrialBinC0All(0x0),
83  fHistoRawYieldDistBinC1All(0x0),
84  fHistoRawYieldTrialBinC1All(0x0),
85  fHistoRawYieldDist(0x0),
86  fHistoRawYieldTrial(0x0),
87  fHistoSigmaTrial(0x0),
88  fHistoMeanTrial(0x0),
89  fHistoChi2Trial(0x0),
90  fHistoSignifTrial(0x0),
91  fHistoBkgTrial(0x0),
92  fHistoBkgInBinEdgesTrial(0x0),
93  fHistoRawYieldDistBinC0(0x0),
94  fHistoRawYieldTrialBinC0(0x0),
95  fHistoRawYieldDistBinC1(0x0),
96  fHistoRawYieldTrialBinC1(0x0),
97  fhTemplRefl(0x0),
98  fhTemplSign(0x0),
99  fFixRefloS(1.),
100  fNtupleMultiTrials(0x0),
101  fNtupleBinCount(0x0),
102  fMinYieldGlob(0),
103  fMaxYieldGlob(0),
104  fMassFitters()
105 {
106  // constructor
107  Int_t rebinStep[4]={3,4,5,6};
108  Double_t minMassStep[6]={1.68,1.70,1.72,1.74,1.76,1.78};
109  Double_t maxMassStep[6]={2.06,2.04,2.02,2.00,1.98,1.96};
110  Double_t nSigmasBC[11]={2.0,2.5,3.0,3.5,4.0,4.5,5.0,5.5,6.0,6.5,7.0};
111  ConfigureRebinSteps(4,rebinStep);
112  ConfigureLowLimFitSteps(6,minMassStep);
113  ConfigureUpLimFitSteps(6,maxMassStep);
114  ConfigurenSigmaBinCSteps(11,nSigmasBC);
115 }
116 
117 //________________________________________________________________________
119  // destructor
120  delete [] fRebinSteps;
121  delete [] fLowLimFitSteps;
122  delete [] fUpLimFitSteps;
123  if(fhTemplRefl) delete fhTemplRefl;
124  if(fhTemplSign) delete fhTemplSign;
125  for (auto fitter : fMassFitters) delete fitter;
126  delete fHistoRawYieldDistAll;
127  delete fHistoRawYieldTrialAll;
128  delete fHistoSigmaTrialAll;
129  delete fHistoMeanTrialAll;
130  delete fHistoChi2TrialAll;
131  delete fHistoSignifTrialAll;
132  delete fHistoBkgTrialAll;
138  for(Int_t ib=0; ib<kNBkgFuncCases; ib++){
139  for(Int_t igs=0; igs<kNFitConfCases; igs++){
140  Int_t theCase=igs*kNBkgFuncCases+ib;
141  delete fHistoRawYieldDist[theCase];
142  delete fHistoRawYieldDistBinC0[theCase];
143  delete fHistoRawYieldDistBinC1[theCase];
144  delete fHistoRawYieldTrial[theCase];
145  delete fHistoRawYieldTrialBinC0[theCase];
146  delete fHistoRawYieldTrialBinC1[theCase];
147  delete fHistoSigmaTrial[theCase];
148  delete fHistoMeanTrial[theCase];
149  delete fHistoChi2Trial[theCase];
150  delete fHistoSignifTrial[theCase];
151  if(fHistoBkgTrial) delete fHistoBkgTrial[theCase];
153  }
154  }
155  delete [] fHistoRawYieldDist;
156  delete [] fHistoRawYieldDistBinC0;
157  delete [] fHistoRawYieldDistBinC1;
158  delete [] fHistoRawYieldTrial;
159  delete [] fHistoRawYieldTrialBinC0;
160  delete [] fHistoRawYieldTrialBinC1;
161  delete [] fHistoSigmaTrial;
162  delete [] fHistoMeanTrial;
163  delete [] fHistoChi2Trial;
164  delete [] fHistoSignifTrial;
165  delete [] fHistoBkgTrial;
166  delete [] fHistoBkgInBinEdgesTrial;
167 
168  delete fNtupleMultiTrials;
169  delete fNtupleBinCount;
170 }
171 
172 //________________________________________________________________________
174  // creates output histograms
175 
176  const Int_t nCases=kNBkgFuncCases*kNFitConfCases;
177 
178  TString funcBkg[kNBkgFuncCases]={"Expo","Lin","Pol2","Pol3","Pol4","Pol5","PowLaw","PowLawExpo"};
179  TString gausSig[kNFitConfCases]={"FixedS","FixedSp20","FixedSm20","FreeS","FixedMeanFixedS","FixedMeanFreeS"};
180 
182 
183  fHistoRawYieldDistAll = new TH1F(Form("hRawYieldDistAll%s",fSuffix.Data())," ; Raw Yield",5000,0.,50000.);
184  fHistoRawYieldTrialAll = new TH1F(Form("hRawYieldTrialAll%s",fSuffix.Data())," ; Trial # ; Raw Yield",nCases*totTrials,-0.5,nCases*totTrials-0.5);
185  fHistoSigmaTrialAll = new TH1F(Form("hSigmaTrialAll%s",fSuffix.Data())," ; Trial # ; Sigma (GeV/c^{2})",nCases*totTrials,-0.5,nCases*totTrials-0.5);
186  fHistoMeanTrialAll = new TH1F(Form("hMeanTrialAll%s",fSuffix.Data())," ; Trial # ; Mean (GeV/c^{2})",nCases*totTrials,-0.5,nCases*totTrials-0.5);
187  fHistoChi2TrialAll = new TH1F(Form("hChi2TrialAll%s",fSuffix.Data())," ; Trial # ; #chi^{2}",nCases*totTrials,-0.5,nCases*totTrials-0.5);
188  fHistoSignifTrialAll = new TH1F(Form("hSignifTrialAll%s",fSuffix.Data())," ; Trial # ; Significance",nCases*totTrials,-0.5,nCases*totTrials-0.5);
189  if(fSaveBkgVal) {
190  fHistoBkgTrialAll = new TH1F(Form("hBkgTrialAll%s",fSuffix.Data())," ; Background",nCases*totTrials,-0.5,nCases*totTrials-0.5);
191  fHistoBkgInBinEdgesTrialAll = new TH1F(Form("hBkgInBinEdgesTrialAll%s",fSuffix.Data())," ; Background in bin edges",nCases*totTrials,-0.5,nCases*totTrials-0.5);
192  }
193 
194 
195  fHistoRawYieldDistBinC0All = new TH1F(Form("hRawYieldDistBinC0All%s",fSuffix.Data())," ; Raw Yield (bin count)",5000,0.,50000.);
196  fHistoRawYieldTrialBinC0All = new TH2F(Form("hRawYieldTrialBinC0All%s",fSuffix.Data())," ; Trial # ; Range for count ; Raw Yield (bin count)",totTrials,-0.5,totTrials-0.5,fNumOfnSigmaBinCSteps,-0.5,fNumOfnSigmaBinCSteps-0.5);
197  fHistoRawYieldDistBinC1All = new TH1F(Form("hRawYieldDistBinC1All%s",fSuffix.Data())," ; Raw Yield (bin count)",5000,0.,50000.);
198  fHistoRawYieldTrialBinC1All = new TH2F(Form("hRawYieldTrialBinC1All%s",fSuffix.Data())," ; Trial # ; Range for count ; Raw Yield (bin count)",totTrials,-0.5,totTrials-0.5,fNumOfnSigmaBinCSteps,-0.5,fNumOfnSigmaBinCSteps-0.5);
199 
200  fHistoRawYieldDist = new TH1F*[nCases];
201  fHistoRawYieldTrial = new TH1F*[nCases];
202  fHistoSigmaTrial = new TH1F*[nCases];
203  fHistoMeanTrial = new TH1F*[nCases];
204  fHistoChi2Trial = new TH1F*[nCases];
205  fHistoSignifTrial = new TH1F*[nCases];
206  if(fSaveBkgVal) {
207  fHistoBkgTrial = new TH1F*[nCases];
208  fHistoBkgInBinEdgesTrial = new TH1F*[nCases];
209  }
210 
211  fHistoRawYieldDistBinC0 = new TH1F*[nCases];
212  fHistoRawYieldTrialBinC0 = new TH2F*[nCases];
213  fHistoRawYieldDistBinC1 = new TH1F*[nCases];
214  fHistoRawYieldTrialBinC1 = new TH2F*[nCases];
215 
216  for(Int_t ib=0; ib<kNBkgFuncCases; ib++){
217  for(Int_t igs=0; igs<kNFitConfCases; igs++){
218  Int_t theCase=igs*kNBkgFuncCases+ib;
219  fHistoRawYieldDist[theCase]=new TH1F(Form("hRawYieldDist%s%s%s",funcBkg[ib].Data(),gausSig[igs].Data(),fSuffix.Data())," ; Raw Yield",5000,0.,50000.);
220  fHistoRawYieldDistBinC0[theCase]=new TH1F(Form("hRawYieldDistBinC0%s%s%s",funcBkg[ib].Data(),gausSig[igs].Data(),fSuffix.Data())," ; Raw Yield (bin count)",5000,0.,50000.);
221  fHistoRawYieldDistBinC1[theCase]=new TH1F(Form("hRawYieldDistBinC1%s%s%s",funcBkg[ib].Data(),gausSig[igs].Data(),fSuffix.Data())," ; Raw Yield (bin count)",5000,0.,50000.);
222  fHistoRawYieldTrial[theCase]=new TH1F(Form("hRawYieldTrial%s%s%s",funcBkg[ib].Data(),gausSig[igs].Data(),fSuffix.Data())," ; Trial # ; Raw Yield",totTrials,-0.5,totTrials-0.5);
223  fHistoRawYieldTrialBinC0[theCase]=new TH2F(Form("hRawYieldTrialBinC0%s%s%s",funcBkg[ib].Data(),gausSig[igs].Data(),fSuffix.Data())," ; Trial # ; Range for count ; Raw Yield (bin count)",totTrials,-0.5,totTrials-0.5,fNumOfnSigmaBinCSteps,-0.5,fNumOfnSigmaBinCSteps-0.5);
224  fHistoRawYieldTrialBinC1[theCase]=new TH2F(Form("hRawYieldTrialBinC1%s%s%s",funcBkg[ib].Data(),gausSig[igs].Data(),fSuffix.Data())," ; Trial # ; Range for count ; Raw Yield (bin count)",totTrials,-0.5,totTrials-0.5,fNumOfnSigmaBinCSteps,-0.5,fNumOfnSigmaBinCSteps-0.5);
225  fHistoSigmaTrial[theCase]=new TH1F(Form("hSigmaTrial%s%s%s",funcBkg[ib].Data(),gausSig[igs].Data(),fSuffix.Data())," ; Trial # ; Sigma (GeV/c^{2})",totTrials,-0.5,totTrials-0.5);
226  fHistoMeanTrial[theCase]=new TH1F(Form("hMeanTrial%s%s%s",funcBkg[ib].Data(),gausSig[igs].Data(),fSuffix.Data())," ; Trial # ; Mean (GeV/c^{2})",totTrials,-0.5,totTrials-0.5);
227  fHistoChi2Trial[theCase]=new TH1F(Form("hChi2Trial%s%s%s",funcBkg[ib].Data(),gausSig[igs].Data(),fSuffix.Data())," ; Trial # ; #chi^{2}",totTrials,-0.5,totTrials-0.5);
228  fHistoSignifTrial[theCase]=new TH1F(Form("hSignifTrial%s%s%s",funcBkg[ib].Data(),gausSig[igs].Data(),fSuffix.Data())," ; Trial # ; Significance",totTrials,-0.5,totTrials-0.5);
229  if(fSaveBkgVal) {
230  fHistoBkgTrial[theCase] = new TH1F(Form("hBkgTrial%s%s%s",funcBkg[ib].Data(),gausSig[igs].Data(),fSuffix.Data())," ; Background",totTrials,-0.5,totTrials-0.5);
231  fHistoBkgInBinEdgesTrial[theCase] = new TH1F(Form("hBkgInBinEdgesTrial%s%s%s",funcBkg[ib].Data(),gausSig[igs].Data(),fSuffix.Data())," ; Background in bin edges",totTrials,-0.5,totTrials-0.5);
232  }
233 
234  fHistoChi2Trial[theCase]->SetMarkerStyle(7);
235  fHistoSignifTrial[theCase]->SetMarkerStyle(7);
236  fHistoSigmaTrial[theCase]->SetMarkerStyle(7);
237  fHistoMeanTrial[theCase]->SetMarkerStyle(7);
238  if(fSaveBkgVal) {
239  fHistoBkgTrial[theCase]->SetMarkerStyle(7);
240  fHistoBkgInBinEdgesTrial[theCase]->SetMarkerStyle(7);
241  }
242 
243  }
244  }
245  fNtupleMultiTrials = new TNtuple(Form("ntuMultiTrial%s",fSuffix.Data()),Form("ntuMultiTrial%s",fSuffix.Data()),"rebin:firstb:minfit:maxfit:bkgfunc:confsig:confmean:chi2:signif:mean:emean:sigma:esigma:rawy:erawy",128000);
246  fNtupleMultiTrials->SetDirectory(nullptr);
247  fNtupleBinCount = new TNtuple(Form("ntuBinCount%s",fSuffix.Data()),Form("ntuBinCount%s",fSuffix.Data()),"rebin:firstb:minfit:maxfit:bkgfunc:confsig:confmean:chi2:nSigmaBC:rawyBC0:erawyBC0:rawyBC1:erawyBC1",128000);
248  fNtupleBinCount->SetDirectory(nullptr);
249  return kTRUE;
250 
251 }
252 
253 //________________________________________________________________________
255  // perform the multiple fits
256 
257  Bool_t hOK=CreateHistos();
258  if(!hOK) return kFALSE;
259 
260  Int_t itrial=0;
261  Int_t types=0;
262  Int_t itrialBC=0;
264 
265  fMinYieldGlob=999999.;
266  fMaxYieldGlob=0.;
267  Float_t xnt[15];
268  Float_t xntBC[13];
269 
270  for(Int_t ir=0; ir<fNumOfRebinSteps; ir++){
271  Int_t rebin=fRebinSteps[ir];
272  for(Int_t iFirstBin=1; iFirstBin<=fNumOfFirstBinSteps; iFirstBin++) {
273  TH1F* hRebinned=0x0;
274  if(fNumOfFirstBinSteps==1) hRebinned=(TH1F*)AliVertexingHFUtils::RebinHisto(hInvMassHisto,rebin,-1);
275  else hRebinned=(TH1F*)AliVertexingHFUtils::RebinHisto(hInvMassHisto,rebin,iFirstBin);
276  for(Int_t iMinMass=0; iMinMass<fNumOfLowLimFitSteps; iMinMass++){
278  Double_t hmin=TMath::Max(minMassForFit,hRebinned->GetBinLowEdge(2));
279  for(Int_t iMaxMass=0; iMaxMass<fNumOfUpLimFitSteps; iMaxMass++){
281  Double_t hmax=TMath::Min(maxMassForFit,hRebinned->GetBinLowEdge(hRebinned->GetNbinsX()));
282  ++itrial;
283  for(Int_t typeb=0; typeb<kNBkgFuncCases; typeb++){
284  if(typeb==kExpoBkg && !fUseExpoBkg) continue;
285  if(typeb==kLinBkg && !fUseLinBkg) continue;
286  if(typeb==kPol2Bkg && !fUsePol2Bkg) continue;
287  if(typeb==kPol3Bkg && !fUsePol3Bkg) continue;
288  if(typeb==kPol4Bkg && !fUsePol4Bkg) continue;
289  if(typeb==kPol5Bkg && !fUsePol5Bkg) continue;
290  if(typeb==kPowBkg && !fUsePowLawBkg) continue;
292  for(Int_t igs=0; igs<kNFitConfCases; igs++){
293  if (igs==kFixSigUpFreeMean && !fUseFixSigUpFreeMean) continue;
294  if (igs==kFixSigDownFreeMean && !fUseFixSigDownFreeMean) continue;
295  if (igs==kFreeSigFixMean && !fUseFixedMeanFreeS) continue;
296  if (igs==kFreeSigFreeMean && !fUseFreeS) continue;
297  if (igs==kFixSigFreeMean && !fUseFixSigFreeMean) continue;
298  if (igs==kFixSigFixMean && !fUseFixSigFixMean) continue;
299  Int_t theCase=igs*kNBkgFuncCases+typeb;
300  Int_t globBin=itrial+theCase*totTrials;
301  for(Int_t j=0; j<15; j++) xnt[j]=0.;
302 
303  Bool_t mustDeleteFitter = kTRUE;
305  if(typeb==kExpoBkg){
306  fitter=new AliHFInvMassFitter(hRebinned, hmin, hmax, AliHFInvMassFitter::kExpo, types);
307  }else if(typeb==kLinBkg){
308  fitter=new AliHFInvMassFitter(hRebinned, hmin, hmax, AliHFInvMassFitter::kLin, types);
309  }else if(typeb==kPol2Bkg){
310  fitter=new AliHFInvMassFitter(hRebinned, hmin, hmax, AliHFInvMassFitter::kPol2, types);
311  }else if(typeb==kPowBkg){
312  fitter=new AliHFInvMassFitter(hRebinned, hmin, hmax, AliHFInvMassFitter::kPow, types);
313  }else if(typeb==kPowTimesExpoBkg){
314  fitter=new AliHFInvMassFitter(hRebinned, hmin, hmax, AliHFInvMassFitter::kPowEx, types);
315  }else{
316  fitter=new AliHFInvMassFitter(hRebinned, hmin, hmax, 6, types);
317  if(typeb==kPol3Bkg) fitter->SetPolDegreeForBackgroundFit(3);
318  if(typeb==kPol4Bkg) fitter->SetPolDegreeForBackgroundFit(4);
319  if(typeb==kPol5Bkg) fitter->SetPolDegreeForBackgroundFit(5);
320  }
321  // D0 Reflection
322  if(fhTemplRefl && fhTemplSign){
323  TH1F *hReflModif=(TH1F*)AliVertexingHFUtils::AdaptTemplateRangeAndBinning(fhTemplRefl,hRebinned,minMassForFit,maxMassForFit);
324  TH1F *hSigModif=(TH1F*)AliVertexingHFUtils::AdaptTemplateRangeAndBinning(fhTemplSign,hRebinned,minMassForFit,maxMassForFit);
325  TH1F* hrfl=fitter->SetTemplateReflections(hReflModif,"2gaus",minMassForFit,maxMassForFit);
326  if(fFixRefloS>0){
327  Double_t fixSoverRefAt=fFixRefloS*(hReflModif->Integral(hReflModif->FindBin(minMassForFit*1.0001),hReflModif->FindBin(maxMassForFit*0.999))/hSigModif->Integral(hSigModif->FindBin(minMassForFit*1.0001),hSigModif->FindBin(maxMassForFit*0.999)));
328  fitter->SetFixReflOverS(fixSoverRefAt);
329  }
330  delete hReflModif;
331  delete hSigModif;
332  }
333  if(fUseSecondPeak){
335  }
336  if(fFitOption==1) fitter->SetUseChi2Fit();
339  xnt[0]=rebin;
340  xnt[1]=iFirstBin;
341  xnt[2]=minMassForFit;
342  xnt[3]=maxMassForFit;
343  xnt[4]=typeb;
344  xnt[6]=0;
345  if(igs==kFixSigFreeMean){
347  xnt[5]=1;
348  }else if(igs==kFixSigUpFreeMean){
350  xnt[5]=2;
351  }else if(igs==kFixSigDownFreeMean){
353  xnt[5]=3;
354  }else if(igs==kFreeSigFreeMean){
355  xnt[5]=0;
356  }else if(igs==kFixSigFixMean){
358  fitter->SetFixGaussianMean(fMassD);
359  xnt[5]=1;
360  xnt[6]=1;
361  }else if(igs==kFreeSigFixMean){
362  fitter->SetFixGaussianMean(fMassD);
363  xnt[5]=0;
364  xnt[6]=1;
365  }
366  Bool_t out=kFALSE;
367  Double_t chisq=-1.;
368  Double_t sigma=0.;
369  Double_t esigma=0.;
370  Double_t pos=.0;
371  Double_t epos=.0;
372  Double_t ry=.0;
373  Double_t ery=.0;
374  Double_t significance=0.;
375  Double_t erSignif=0.;
376  Double_t bkg=0.;
377  Double_t erbkg=0.;
378  Double_t bkgBEdge=0;
379  Double_t erbkgBEdge=0;
380  if(typeb<kNBkgFuncCases){
381  printf("****** START FIT OF HISTO %s WITH REBIN %d FIRST BIN %d MASS RANGE %f-%f BACKGROUND FIT FUNCTION=%d CONFIG SIGMA/MEAN=%d\n",hInvMassHisto->GetName(),rebin,iFirstBin,minMassForFit,maxMassForFit,typeb,igs);
382  out=fitter->MassFitter(0);
383  chisq=fitter->GetReducedChiSquare();
384  fitter->Significance(fnSigmaForBkgEval,significance,erSignif);
385  sigma=fitter->GetSigma();
386  pos=fitter->GetMean();
387  esigma=fitter->GetSigmaUncertainty();
388  if(esigma<0.00001) esigma=0.0001;
389  epos=fitter->GetMeanUncertainty();
390  if(epos<0.00001) epos=0.0001;
391  ry=fitter->GetRawYield();
392  ery=fitter->GetRawYieldError();
393  fitter->Background(fnSigmaForBkgEval,bkg,erbkg);
394  Double_t minval = hInvMassHisto->GetXaxis()->GetBinLowEdge(hInvMassHisto->FindBin(pos-fnSigmaForBkgEval*sigma));
395  Double_t maxval = hInvMassHisto->GetXaxis()->GetBinUpEdge(hInvMassHisto->FindBin(pos+fnSigmaForBkgEval*sigma));
396  fitter->Background(minval,maxval,bkgBEdge,erbkgBEdge);
397  if(out && fDrawIndividualFits && thePad){
398  thePad->Clear();
399  fitter->DrawHere(thePad, fnSigmaForBkgEval);
400  fMassFitters.push_back(fitter);
401  mustDeleteFitter = kFALSE;
402  for (auto format : fInvMassFitSaveAsFormats) {
403  thePad->SaveAs(Form("FitOutput_%s_Trial%d.%s",hInvMassHisto->GetName(),globBin, format.c_str()));
404  }
405  }
406  }
407  // else{
408  // out=DoFitWithPol3Bkg(hRebinned,hmin,hmax,igs);
409  // if(out && thePad){
410  // thePad->Clear();
411  // hRebinned->Draw();
412  // TF1* fSB=(TF1*)hRebinned->GetListOfFunctions()->FindObject("fSB");
413  // fB1=new TF1("fB1","[0]+[1]*x+[2]*x*x+[3]*x*x*x",hmin,hmax);
414  // for(Int_t j=0; j<4; j++) fB1->SetParameter(j,fSB->GetParameter(3+j));
415  // fB1->SetLineColor(2);
416  // fB1->Draw("same");
417  // fSB->SetLineColor(4);
418  // fSB->Draw("same");
419  // thePad->Update();
420  // chisq=fSB->GetChisquare()/fSB->GetNDF();;
421  // sigma=fSB->GetParameter(2);
422  // esigma=fSB->GetParError(2);
423  // if(esigma<0.00001) esigma=0.0001;
424  // pos=fSB->GetParameter(1);
425  // epos=fSB->GetParError(1);
426  // if(epos<0.00001) epos=0.0001;
427  // ry=fSB->GetParameter(0)/hRebinned->GetBinWidth(1);
428  // ery=fSB->GetParError(0)/hRebinned->GetBinWidth(1);
429  // }
430  // }
431  xnt[7]=chisq;
432  if(out && chisq>0. && sigma>0.5*fSigmaGausMC && sigma<2.0*fSigmaGausMC){
433  xnt[8]=significance;
434  xnt[9]=pos;
435  xnt[10]=epos;
436  xnt[11]=sigma;
437  xnt[12]=esigma;
438  xnt[13]=ry;
439  xnt[14]=ery;
440  fHistoRawYieldDistAll->Fill(ry);
441  fHistoRawYieldTrialAll->SetBinContent(globBin,ry);
442  fHistoRawYieldTrialAll->SetBinError(globBin,ery);
443  fHistoSigmaTrialAll->SetBinContent(globBin,sigma);
444  fHistoSigmaTrialAll->SetBinError(globBin,esigma);
445  fHistoMeanTrialAll->SetBinContent(globBin,pos);
446  fHistoMeanTrialAll->SetBinError(globBin,epos);
447  fHistoChi2TrialAll->SetBinContent(globBin,chisq);
448  fHistoChi2TrialAll->SetBinError(globBin,0.00001);
449  fHistoSignifTrialAll->SetBinContent(globBin,significance);
450  fHistoSignifTrialAll->SetBinError(globBin,erSignif);
451  if(fSaveBkgVal) {
452  fHistoBkgTrialAll->SetBinContent(globBin,bkg);
453  fHistoBkgTrialAll->SetBinError(globBin,erbkg);
454  fHistoBkgInBinEdgesTrialAll->SetBinContent(globBin,bkgBEdge);
455  fHistoBkgInBinEdgesTrialAll->SetBinError(globBin,erbkgBEdge);
456  }
457 
458  if(ry<fMinYieldGlob) fMinYieldGlob=ry;
459  if(ry>fMaxYieldGlob) fMaxYieldGlob=ry;
460  fHistoRawYieldDist[theCase]->Fill(ry);
461  fHistoRawYieldTrial[theCase]->SetBinContent(itrial,ry);
462  fHistoRawYieldTrial[theCase]->SetBinError(itrial,ery);
463  fHistoSigmaTrial[theCase]->SetBinContent(itrial,sigma);
464  fHistoSigmaTrial[theCase]->SetBinError(itrial,esigma);
465  fHistoMeanTrial[theCase]->SetBinContent(itrial,pos);
466  fHistoMeanTrial[theCase]->SetBinError(itrial,epos);
467  fHistoChi2Trial[theCase]->SetBinContent(itrial,chisq);
468  fHistoChi2Trial[theCase]->SetBinError(itrial,0.00001);
469  fHistoSignifTrial[theCase]->SetBinContent(itrial,significance);
470  fHistoSignifTrial[theCase]->SetBinError(itrial,erSignif);
471  if(fSaveBkgVal) {
472  fHistoBkgTrial[theCase]->SetBinContent(itrial,bkg);
473  fHistoBkgTrial[theCase]->SetBinError(itrial,erbkg);
474  fHistoBkgInBinEdgesTrial[theCase]->SetBinContent(itrial,bkgBEdge);
475  fHistoBkgInBinEdgesTrial[theCase]->SetBinError(itrial,erbkgBEdge);
476  }
477  fNtupleMultiTrials->Fill(xnt);
478 
479  for(Int_t j=0; j<8; j++) xntBC[j]=xnt[j];
480 
481 
482  for(Int_t iStepBC=0; iStepBC<fNumOfnSigmaBinCSteps; iStepBC++){
483  Double_t minMassBC=fMassD-fnSigmaBinCSteps[iStepBC]*sigma;
484  Double_t maxMassBC=fMassD+fnSigmaBinCSteps[iStepBC]*sigma;
485  if(minMassBC>minMassForFit &&
486  maxMassBC<maxMassForFit &&
487  minMassBC>(hRebinned->GetXaxis()->GetXmin()) &&
488  maxMassBC<(hRebinned->GetXaxis()->GetXmax())){
489  Double_t cnts0,ecnts0;
490  Double_t cnts1,ecnts1;
491  cnts0=fitter->GetRawYieldBinCounting(ecnts0,minMassBC,maxMassBC,0);
492  cnts1=fitter->GetRawYieldBinCounting(ecnts1,minMassBC,maxMassBC,1);
493  xntBC[8]=fnSigmaBinCSteps[iStepBC];
494  xntBC[9]=cnts0;
495  xntBC[10]=ecnts0;
496  xntBC[11]=cnts1;
497  xntBC[12]=ecnts1;
498  ++itrialBC;
499  fHistoRawYieldDistBinC0All->Fill(cnts0);
500  fHistoRawYieldTrialBinC0All->SetBinContent(globBin,iStepBC+1,cnts0);
501  fHistoRawYieldTrialBinC0All->SetBinError(globBin,iStepBC+1,ecnts0);
502  fHistoRawYieldTrialBinC0[theCase]->SetBinContent(itrial,iStepBC+1,cnts0);
503  fHistoRawYieldTrialBinC0[theCase]->SetBinError(itrial,iStepBC+1,ecnts0);
504  fHistoRawYieldDistBinC0[theCase]->Fill(cnts0);
505  fHistoRawYieldDistBinC1All->Fill(cnts1);
506  fHistoRawYieldTrialBinC1All->SetBinContent(globBin,iStepBC+1,cnts1);
507  fHistoRawYieldTrialBinC1All->SetBinError(globBin,iStepBC+1,ecnts1);
508  fHistoRawYieldTrialBinC1[theCase]->SetBinContent(itrial,iStepBC+1,cnts1);
509  fHistoRawYieldTrialBinC1[theCase]->SetBinError(itrial,iStepBC+1,ecnts1);
510  fHistoRawYieldDistBinC1[theCase]->Fill(cnts1);
511  fNtupleBinCount->Fill(xntBC);
512  }
513  }
514  }
515  if (mustDeleteFitter) delete fitter;
516  }
517  }
518  }
519  }
520  delete hRebinned;
521  }
522  }
523  return kTRUE;
524 }
525 
526 //________________________________________________________________________
528  // save histos in a root file for further analysis
529  const Int_t nCases=kNBkgFuncCases*kNFitConfCases;
530  TFile outHistos(fileName.Data(),option.Data());
531  if (outHistos.IsZombie()) {
532  Printf("Could not open file '%s'!", fileName.Data());
533  return;
534  }
535  outHistos.cd();
536  fHistoRawYieldTrialAll->Write();
537  fHistoSigmaTrialAll->Write();
538  fHistoMeanTrialAll->Write();
539  fHistoChi2TrialAll->Write();
540  fHistoSignifTrialAll->Write();
541  if(fSaveBkgVal) {
542  fHistoBkgTrialAll->Write();
544  }
549  for(Int_t ic=0; ic<nCases; ic++){
550  fHistoRawYieldTrial[ic]->Write();
551  fHistoSigmaTrial[ic]->Write();
552  fHistoMeanTrial[ic]->Write();
553  fHistoChi2Trial[ic]->Write();
554  fHistoSignifTrial[ic]->Write();
555  if(fSaveBkgVal) {
556  fHistoBkgTrial[ic]->Write();
557  fHistoBkgInBinEdgesTrial[ic]->Write();
558  }
559  fHistoRawYieldTrialBinC0[ic]->Write();
560  fHistoRawYieldDistBinC0[ic]->Write();
561  fHistoRawYieldTrialBinC1[ic]->Write();
562  fHistoRawYieldDistBinC1[ic]->Write();
563  }
564  fNtupleMultiTrials->SetDirectory(&outHistos);
565  fNtupleMultiTrials->Write();
566  fNtupleBinCount->SetDirectory(&outHistos);
567  fNtupleBinCount->Write();
568  outHistos.Close();
569 }
570 
571 //________________________________________________________________________
572 void AliHFInvMassMultiTrialFit::DrawHistos(TCanvas* cry) const{
573  // draw histos
574  cry->Divide(2,2);
575  cry->cd(1);
576  fHistoSigmaTrialAll->Draw();
577  cry->cd(3);
578  fHistoChi2TrialAll->Draw();
579  cry->cd(2);
580  fHistoRawYieldTrialAll->Draw();
581  cry->cd(4);
582  fHistoRawYieldDistAll->Draw();
583  TLatex* tmean=new TLatex(0.15,0.8,Form("mean=%.1f",fHistoRawYieldDistAll->GetMean()));
584  tmean->SetNDC();
585  tmean->Draw();
586  TLatex* thrms=new TLatex(0.15,0.72,Form("rms=%.1f",fHistoRawYieldDistAll->GetRMS()));
587  thrms->SetNDC();
588  thrms->Draw();
589  TLatex* tmax=new TLatex(0.6,0.8,Form("max=%.1f",fMaxYieldGlob));
590  tmax->SetNDC();
591  tmax->Draw();
592  TLatex* tmin=new TLatex(0.6,0.72,Form("min=%.1f",fMinYieldGlob));
593  tmin->SetNDC();
594  tmin->Draw();
595  TLatex* trms=new TLatex(0.6,0.64,Form("(max-min)/sqrt(12)=%.1f",(fMaxYieldGlob-fMinYieldGlob)/sqrt(12)));
596  trms->SetNDC();
597  trms->Draw();
598 
599 }
600 //________________________________________________________________________
602  Int_t iCase){
603  //
604 
605  TH1F *hCutTmp=(TH1F*)histoToFit->Clone("hCutTmp");
606  for(Int_t ib=1; ib<=hCutTmp->GetNbinsX(); ib++){
607  Double_t xc=hCutTmp->GetBinCenter(ib);
608  if(xc>(fMassD-5.*fSigmaGausMC) && xc<(fMassD+5.*fSigmaGausMC)){
609  hCutTmp->SetBinContent(ib,0.);
610  hCutTmp->SetBinError(ib,0.);
611  }
612  }
613 
614  hCutTmp->Fit("pol2","E0","",hmin,hmax);
615  TF1* f2=(TF1*)hCutTmp->GetListOfFunctions()->FindObject("pol2");
616  TF1* f3=new TF1("myPol3","pol3");
617  for(Int_t i=0; i<3;i++) f3->SetParameter(i,f2->GetParameter(i));
618  hCutTmp->Fit(f3,"E0","",hmin,hmax);
619  Double_t quickCount=0.;
620  for(Int_t ib=1; ib<=histoToFit->GetNbinsX(); ib++){
621  Double_t xc=hCutTmp->GetBinCenter(ib);
622  if(xc>(fMassD-3.*fSigmaGausMC) && xc<(fMassD+3.*fSigmaGausMC)){
623  quickCount+=(histoToFit->GetBinContent(ib)-f3->Eval(xc));
624  }
625  }
626  TF1* fSB=new TF1("fSB","[0]*1./(TMath::Sqrt(2.*TMath::Pi())*[2])*TMath::Exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))+[3]+[4]*x+[5]*x*x+[6]*x*x*x");
627  fSB->SetParameter(0,quickCount);
628  fSB->SetParameter(1,fMassD);
629  fSB->SetParameter(2,fSigmaGausMC);
630  for(Int_t j=0; j<4; j++) fSB->SetParameter(j+3,f3->GetParameter(j));
631  if(iCase==0) fSB->FixParameter(2,fSigmaGausMC);
632  else if(iCase==1) fSB->FixParameter(2,fSigmaGausMC*(1.+fSigmaMCVariation));
633  else if(iCase==2) fSB->FixParameter(2,fSigmaGausMC*(1.-fSigmaMCVariation));
634  else if(iCase==4){
635  fSB->FixParameter(1,fMassD);
636  fSB->FixParameter(2,fSigmaGausMC);
637  } else if(iCase==5){
638  fSB->FixParameter(1,fMassD);
639  }
640  histoToFit->Fit(fSB,"ME0","",hmin,hmax);
641  // quality cuts
642  if(fSB->GetParError(0)<0.01*fSB->GetParameter(0)) return kFALSE;
643  if(fSB->GetParError(0)>0.6*fSB->GetParameter(0)) return kFALSE;
644 
645  delete hCutTmp;
646  return kTRUE;
647 }
648 //__________________________________________________________________________________
649 void AliHFInvMassMultiTrialFit::SetTemplatesForReflections(const TH1F *hr, const TH1F *hs) {
651  if(fhTemplSign) delete fhTemplSign;
652  if(fhTemplRefl) delete fhTemplRefl;
653  fhTemplRefl=(TH1F*)hr->Clone("hTemplRefl");
654  fhTemplSign=(TH1F*)hs->Clone("hTemplSign");
655  return;
656 }
657 
Bool_t fUsePol4Bkg
switch for pol3 background
Double_t GetMeanUncertainty() const
Int_t * fRebinSteps
number of rebin steps
static TH1D * RebinHisto(TH1 *hOrig, Int_t reb, Int_t firstUse=-1)
Rebinning of invariant mass histograms.
Bool_t fDrawIndividualFits
switch for saving bkg values in nsigma
TH1F ** fHistoChi2Trial
histo with gauss mean from subsamples of trials
TH1F * fHistoBkgInBinEdgesTrialAll
histo with bkg from all trials
Bool_t fUsePowLawTimesExpoBkg
switch for power law background
double Double_t
Definition: External.C:58
void SetFixReflOverS(Double_t rovers)
Bool_t fUseFixSigFixMean
switch for FixSigFreeMean
void DrawHere(TVirtualPad *c, Double_t nsigma=3, Int_t writeFitInfo=1)
void IncludeSecondGausPeak(Double_t mass, Bool_t fixm, Double_t width, Bool_t fixw)
Definition: External.C:236
Double_t fSigmaSecondPeak
position of the 2nd peak
TString fileName
Float_t fFixRefloS
template of signal contribution
TString format
file names tag, basically the trigger and calorimeter combination
Double_t fSigmaMCVariation
sigma of D meson peak from MC
TH1F ** fHistoRawYieldDist
histo with bin counts from all trials
TH2F ** fHistoRawYieldTrialBinC1
histo with bin counts from subsamples of trials
TH1F ** fHistoBkgInBinEdgesTrial
histo with bkg from subsamples of trials
Double_t * fnSigmaBinCSteps
number of steps on the bin counting
Int_t MassFitter(Bool_t draw=kTRUE)
TH1F * fHistoBkgTrialAll
histo with chi2 from all trials
AliHFInvMassFitter class for the fit of invariant mass distribution of charm hadrons.
void ConfigureLowLimFitSteps(Int_t nSteps, Double_t *values)
std::vector< AliHFInvMassFitter * > fMassFitters
maximum yield
Int_t types
static TH1 * AdaptTemplateRangeAndBinning(const TH1 *hRef, TH1 *hData, Double_t minFit, Double_t maxFit)
TH1F * fHistoChi2TrialAll
histo with gauss mean from all trials
Double_t GetSigma() const
void SetPolDegreeForBackgroundFit(Int_t deg)
void ConfigureUpLimFitSteps(Int_t nSteps, Double_t *values)
TH1F ** fHistoSignifTrial
histo with chi2 from subsamples of trials
Double_t fMaxYieldGlob
minimum yield
Bool_t fUseFixSigFreeMean
switch for FixedMeanFreeS
TH1F ** fHistoRawYieldDistBinC0
histo with bkg in mass bin edges from subsamples of trials
Bool_t fUsePol3Bkg
switch for pol2 background
void SetFixGaussianSigma(Double_t sigma)
Bool_t fUseFixSigUpFreeMean
switch for power law background
TH1F * fHistoMeanTrialAll
histo with gauss sigma from all trials
Double_t fMassD
relative variation of the sigma
Int_t typeb
Bool_t fFixSigmaSecondPeak
flag to fix the position of the 2nd peak
Double_t * sigma
TH1F * fHistoRawYieldDistBinC0All
histo with bkg in mass bin edges from all trials
TH1F * fHistoSignifTrialAll
histo with chi2 from all trials
Bool_t fUseSecondPeak
switch for FixSigFixMean
Bool_t fFixMassSecondPeak
width of the 2nd peak
void SetFixGaussianMean(Double_t mean)
int Int_t
Definition: External.C:63
TH1F ** fHistoRawYieldTrial
histo with yield from subsamples of trials
Double_t GetMean() const
void SetTemplatesForReflections(const TH1F *hTemplRefl, const TH1F *hTemplSig)
Int_t fNumOfLowLimFitSteps
number of steps in the first bin for rebin
float Float_t
Definition: External.C:68
TH1F * fHistoRawYieldDistBinC1All
histo with bin counts from all trials
Bool_t DoFitWithPol3Bkg(TH1F *histoToFit, Double_t hmin, Double_t hmax, Int_t theCase)
TH1F * fHistoRawYieldTrialAll
histo with yield from all trials
Double_t * fUpLimFitSteps
number of steps on the max. mass for fit
Double_t GetReducedChiSquare() const
Double_t fMassSecondPeak
switch off/on second peak (for D+->KKpi in Ds)
Bool_t DoMultiTrials(TH1D *hInvMassHisto, TPad *thePad=0x0)
Definition: External.C:212
Int_t totTrials
TH1F ** fHistoRawYieldDistBinC1
histo with bin counts from subsamples of trials
Double_t GetRawYieldBinCounting(Double_t &errRyBC, Double_t nSigma=3., Int_t option=0, Int_t pdgCode=0) const
void ConfigureRebinSteps(Int_t nSteps, Int_t *values)
TH1F * fHistoSigmaTrialAll
histo with yield from all trials
void SaveToRoot(TString fileName, TString option="recreate") const
TH1F * fhTemplRefl
histo with bin counts from subsamples of trials
TH1F * SetTemplateReflections(const TH1 *h, TString opt, Double_t minRange, Double_t maxRange)
TH1F * fHistoRawYieldDistAll
flag for drawing fits
TH2F ** fHistoRawYieldTrialBinC0
histo with bin counts from subsamples of trials
TH1F ** fHistoSigmaTrial
histo with yield from subsamples of trials
TList * fitter
Definition: DrawAnaELoss.C:26
Bool_t fUsePowLawBkg
switch for pol5 background
Bool_t fUseFreeS
switch for FixSigDownFreeMean
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)
Int_t fFitOption
name to characterize analysis case
Bool_t fUsePol5Bkg
switch for pol4 background
Bool_t fSaveBkgVal
flag to fix the width of the 2nd peak
Bool_t fUseFixedMeanFreeS
switch for FreeSigma
Double_t GetSigmaUncertainty() const
void Significance(Double_t nOfSigma, Double_t &significance, Double_t &errsignificance) const
TH1F * fhTemplSign
template of reflection contribution
Bool_t fUseLinBkg
switch for exponential background
Double_t GetRawYieldError() const
Bool_t fUseFixSigDownFreeMean
switch for FixSigUpFreeMean
void Background(Double_t nOfSigma, Double_t &background, Double_t &errbackground) const
void SetInitialGaussianSigma(Double_t sigma)
void SetInitialGaussianMean(Double_t mean)
Bool_t fUseExpoBkg
LL or chi2 fit.
Int_t fNumOfRebinSteps
saves the invariant mass fit canvases in the file formats listed in this vector (if empty...
Bool_t fUsePol2Bkg
switch for linear background
Int_t rebin
Double_t * fLowLimFitSteps
number of steps on the min. mass for fit
void DrawHistos(TCanvas *cry) const
TH1F ** fHistoMeanTrial
histo with gauss sigma from subsamples of trials
bool Bool_t
Definition: External.C:53
std::set< std::string > fInvMassFitSaveAsFormats
void ConfigurenSigmaBinCSteps(Int_t nSteps, Double_t *values)
TH1F ** fHistoBkgTrial
histo with chi2 from subsamples of trials
TH2F * fHistoRawYieldTrialBinC0All
histo with bin counts from all trials
Double_t minMassForFit
Double_t GetRawYield() const
TH2F * fHistoRawYieldTrialBinC1All
histo with bin counts from all trials
Double_t maxMassForFit