AliPhysics  2c6b7ad (2c6b7ad)
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  fUse2GausSignal(kFALSE),
61  fUse2GausSigmaRatioSignal(kFALSE),
62  fFixSecondGausSig(-1.),
63  fFixSecondGausFrac(-1.),
64  fFixSecondGausSigRat(-1.),
65  fUseFixSigUpFreeMean(kTRUE),
66  fUseFixSigDownFreeMean(kTRUE),
67  fUseFreeS(kTRUE),
68  fUseFixedMeanFreeS(kTRUE),
69  fUseFixSigFreeMean(kTRUE),
70  fUseFixSigFixMean(kTRUE),
71  fUseSecondPeak(kFALSE),
72  fMassSecondPeak(1.86958),
73  fSigmaSecondPeak(0.01),
74  fFixMassSecondPeak(kFALSE),
75  fFixSigmaSecondPeak(kFALSE),
76  fSaveBkgVal(kFALSE),
77  fDrawIndividualFits(kFALSE),
78  fHistoRawYieldDistAll(0x0),
79  fHistoRawYieldTrialAll(0x0),
80  fHistoSigmaTrialAll(0x0),
81  fHistoMeanTrialAll(0x0),
82  fHistoChi2TrialAll(0x0),
83  fHistoSignifTrialAll(0x0),
84  fHistoBkgTrialAll(0x0),
85  fHistoBkgInBinEdgesTrialAll(0x0),
86  fHistoRawYieldDistBinC0All(0x0),
87  fHistoRawYieldTrialBinC0All(0x0),
88  fHistoRawYieldDistBinC1All(0x0),
89  fHistoRawYieldTrialBinC1All(0x0),
90  fHistoRawYieldDist(0x0),
91  fHistoRawYieldTrial(0x0),
92  fHistoSigmaTrial(0x0),
93  fHistoMeanTrial(0x0),
94  fHistoChi2Trial(0x0),
95  fHistoSignifTrial(0x0),
96  fHistoBkgTrial(0x0),
97  fHistoBkgInBinEdgesTrial(0x0),
98  fHistoRawYieldDistBinC0(0x0),
99  fHistoRawYieldTrialBinC0(0x0),
100  fHistoRawYieldDistBinC1(0x0),
101  fHistoRawYieldTrialBinC1(0x0),
102  fhTemplRefl(0x0),
103  fhTemplSign(0x0),
104  fFixRefloS(1.),
105  fNtupleMultiTrials(0x0),
106  fNtupleBinCount(0x0),
107  fMinYieldGlob(0),
108  fMaxYieldGlob(0),
109  fMassFitters()
110 {
111  // constructor
112  Int_t rebinStep[4]={3,4,5,6};
113  Double_t minMassStep[6]={1.68,1.70,1.72,1.74,1.76,1.78};
114  Double_t maxMassStep[6]={2.06,2.04,2.02,2.00,1.98,1.96};
115  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};
116  ConfigureRebinSteps(4,rebinStep);
117  ConfigureLowLimFitSteps(6,minMassStep);
118  ConfigureUpLimFitSteps(6,maxMassStep);
119  ConfigurenSigmaBinCSteps(11,nSigmasBC);
120 }
121 
122 //________________________________________________________________________
124  // destructor
125  delete [] fRebinSteps;
126  delete [] fLowLimFitSteps;
127  delete [] fUpLimFitSteps;
128  if(fhTemplRefl) delete fhTemplRefl;
129  if(fhTemplSign) delete fhTemplSign;
130  for (auto fitter : fMassFitters) delete fitter;
131  delete fHistoRawYieldDistAll;
132  delete fHistoRawYieldTrialAll;
133  delete fHistoSigmaTrialAll;
134  delete fHistoMeanTrialAll;
135  delete fHistoChi2TrialAll;
136  delete fHistoSignifTrialAll;
137  delete fHistoBkgTrialAll;
143  for(Int_t typeb=0; typeb<kNBkgFuncCases; typeb++){
144  for(Int_t types=0; types<kNSigFuncCases; types++){
145  for(Int_t igs=0; igs<kNFitConfCases; igs++){
146  Int_t theCase=igs*kNBkgFuncCases*kNSigFuncCases+types*kNBkgFuncCases+typeb;
147  delete fHistoRawYieldDist[theCase];
148  delete fHistoRawYieldDistBinC0[theCase];
149  delete fHistoRawYieldDistBinC1[theCase];
150  delete fHistoRawYieldTrial[theCase];
151  delete fHistoRawYieldTrialBinC0[theCase];
152  delete fHistoRawYieldTrialBinC1[theCase];
153  delete fHistoSigmaTrial[theCase];
154  delete fHistoMeanTrial[theCase];
155  delete fHistoChi2Trial[theCase];
156  delete fHistoSignifTrial[theCase];
157  if(fHistoBkgTrial) delete fHistoBkgTrial[theCase];
159  }
160  }
161  }
162  delete [] fHistoRawYieldDist;
163  delete [] fHistoRawYieldDistBinC0;
164  delete [] fHistoRawYieldDistBinC1;
165  delete [] fHistoRawYieldTrial;
166  delete [] fHistoRawYieldTrialBinC0;
167  delete [] fHistoRawYieldTrialBinC1;
168  delete [] fHistoSigmaTrial;
169  delete [] fHistoMeanTrial;
170  delete [] fHistoChi2Trial;
171  delete [] fHistoSignifTrial;
172  delete [] fHistoBkgTrial;
173  delete [] fHistoBkgInBinEdgesTrial;
174 
175  delete fNtupleMultiTrials;
176  delete fNtupleBinCount;
177 }
178 
179 //________________________________________________________________________
181  // creates output histograms
182 
184 
185  TString funcBkg[kNBkgFuncCases]={"Expo","Lin","Pol2","Pol3","Pol4","Pol5","PowLaw","PowLawExpo"};
186  TString funcSig[kNSigFuncCases]={"","2Gaus","2GausSigmaRatio"};
187  TString gausSig[kNFitConfCases]={"FixedS","FixedSp20","FixedSm20","FreeS","FixedMeanFixedS","FixedMeanFreeS"};
188 
190 
191  fHistoRawYieldDistAll = new TH1F(Form("hRawYieldDistAll%s",fSuffix.Data())," ; Raw Yield",5000,0.,50000.);
192  fHistoRawYieldTrialAll = new TH1F(Form("hRawYieldTrialAll%s",fSuffix.Data())," ; Trial # ; Raw Yield",nCases*totTrials,-0.5,nCases*totTrials-0.5);
193  fHistoSigmaTrialAll = new TH1F(Form("hSigmaTrialAll%s",fSuffix.Data())," ; Trial # ; Sigma (GeV/c^{2})",nCases*totTrials,-0.5,nCases*totTrials-0.5);
194  fHistoMeanTrialAll = new TH1F(Form("hMeanTrialAll%s",fSuffix.Data())," ; Trial # ; Mean (GeV/c^{2})",nCases*totTrials,-0.5,nCases*totTrials-0.5);
195  fHistoChi2TrialAll = new TH1F(Form("hChi2TrialAll%s",fSuffix.Data())," ; Trial # ; #chi^{2}",nCases*totTrials,-0.5,nCases*totTrials-0.5);
196  fHistoSignifTrialAll = new TH1F(Form("hSignifTrialAll%s",fSuffix.Data())," ; Trial # ; Significance",nCases*totTrials,-0.5,nCases*totTrials-0.5);
197  if(fSaveBkgVal) {
198  fHistoBkgTrialAll = new TH1F(Form("hBkgTrialAll%s",fSuffix.Data())," ; Background",nCases*totTrials,-0.5,nCases*totTrials-0.5);
199  fHistoBkgInBinEdgesTrialAll = new TH1F(Form("hBkgInBinEdgesTrialAll%s",fSuffix.Data())," ; Background in bin edges",nCases*totTrials,-0.5,nCases*totTrials-0.5);
200  }
201 
202 
203  fHistoRawYieldDistBinC0All = new TH1F(Form("hRawYieldDistBinC0All%s",fSuffix.Data())," ; Raw Yield (bin count)",5000,0.,50000.);
204  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);
205  fHistoRawYieldDistBinC1All = new TH1F(Form("hRawYieldDistBinC1All%s",fSuffix.Data())," ; Raw Yield (bin count)",5000,0.,50000.);
206  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);
207 
208  fHistoRawYieldDist = new TH1F*[nCases];
209  fHistoRawYieldTrial = new TH1F*[nCases];
210  fHistoSigmaTrial = new TH1F*[nCases];
211  fHistoMeanTrial = new TH1F*[nCases];
212  fHistoChi2Trial = new TH1F*[nCases];
213  fHistoSignifTrial = new TH1F*[nCases];
214  if(fSaveBkgVal) {
215  fHistoBkgTrial = new TH1F*[nCases];
216  fHistoBkgInBinEdgesTrial = new TH1F*[nCases];
217  }
218 
219  fHistoRawYieldDistBinC0 = new TH1F*[nCases];
220  fHistoRawYieldTrialBinC0 = new TH2F*[nCases];
221  fHistoRawYieldDistBinC1 = new TH1F*[nCases];
222  fHistoRawYieldTrialBinC1 = new TH2F*[nCases];
223 
224  for(Int_t typeb=0; typeb<kNBkgFuncCases; typeb++){
225  for(Int_t types=0; types<kNSigFuncCases; types++){
226  for(Int_t igs=0; igs<kNFitConfCases; igs++){
227  Int_t theCase=igs*kNBkgFuncCases*kNSigFuncCases+types*kNBkgFuncCases+typeb;
228  printf("typeb= %d types = %d igs = %d theCase=%d out of %d\n",typeb,types,igs,theCase,nCases);
229  fHistoRawYieldDist[theCase]=new TH1F(Form("hRawYieldDist%s%s%s%s",funcBkg[typeb].Data(),funcSig[types].Data(),gausSig[igs].Data(),fSuffix.Data())," ; Raw Yield",5000,0.,50000.);
230  fHistoRawYieldDistBinC0[theCase]=new TH1F(Form("hRawYieldDistBinC0%s%s%s%s",funcBkg[typeb].Data(),funcSig[types].Data(),gausSig[igs].Data(),fSuffix.Data())," ; Raw Yield (bin count)",5000,0.,50000.);
231  fHistoRawYieldDistBinC1[theCase]=new TH1F(Form("hRawYieldDistBinC1%s%s%s%s",funcBkg[typeb].Data(),funcSig[types].Data(),gausSig[igs].Data(),fSuffix.Data())," ; Raw Yield (bin count)",5000,0.,50000.);
232  fHistoRawYieldTrial[theCase]=new TH1F(Form("hRawYieldTrial%s%s%s%s",funcBkg[typeb].Data(),funcSig[types].Data(),gausSig[igs].Data(),fSuffix.Data())," ; Trial # ; Raw Yield",totTrials,-0.5,totTrials-0.5);
233  fHistoRawYieldTrialBinC0[theCase]=new TH2F(Form("hRawYieldTrialBinC0%s%s%s%s",funcBkg[typeb].Data(),funcSig[types].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);
234  fHistoRawYieldTrialBinC1[theCase]=new TH2F(Form("hRawYieldTrialBinC1%s%s%s%s",funcBkg[typeb].Data(),funcSig[types].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);
235  fHistoSigmaTrial[theCase]=new TH1F(Form("hSigmaTrial%s%s%s%s",funcBkg[typeb].Data(),funcSig[types].Data(),gausSig[igs].Data(),fSuffix.Data())," ; Trial # ; Sigma (GeV/c^{2})",totTrials,-0.5,totTrials-0.5);
236  fHistoMeanTrial[theCase]=new TH1F(Form("hMeanTrial%s%s%s%s",funcBkg[typeb].Data(),funcSig[types].Data(),gausSig[igs].Data(),fSuffix.Data())," ; Trial # ; Mean (GeV/c^{2})",totTrials,-0.5,totTrials-0.5);
237  fHistoChi2Trial[theCase]=new TH1F(Form("hChi2Trial%s%s%s%s",funcBkg[typeb].Data(),funcSig[types].Data(),gausSig[igs].Data(),fSuffix.Data())," ; Trial # ; #chi^{2}",totTrials,-0.5,totTrials-0.5);
238  fHistoSignifTrial[theCase]=new TH1F(Form("hSignifTrial%s%s%s%s",funcBkg[typeb].Data(),funcSig[types].Data(),gausSig[igs].Data(),fSuffix.Data())," ; Trial # ; Significance",totTrials,-0.5,totTrials-0.5);
239  if(fSaveBkgVal) {
240  fHistoBkgTrial[theCase] = new TH1F(Form("hBkgTrial%s%s%s%s",funcBkg[typeb].Data(),funcSig[types].Data(),gausSig[igs].Data(),fSuffix.Data())," ; Background",totTrials,-0.5,totTrials-0.5);
241  fHistoBkgInBinEdgesTrial[theCase] = new TH1F(Form("hBkgInBinEdgesTrial%s%s%s%s",funcBkg[typeb].Data(),funcSig[types].Data(),gausSig[igs].Data(),fSuffix.Data())," ; Background in bin edges",totTrials,-0.5,totTrials-0.5);
242  }
243 
244  fHistoChi2Trial[theCase]->SetMarkerStyle(7);
245  fHistoSignifTrial[theCase]->SetMarkerStyle(7);
246  fHistoSigmaTrial[theCase]->SetMarkerStyle(7);
247  fHistoMeanTrial[theCase]->SetMarkerStyle(7);
248  if(fSaveBkgVal) {
249  fHistoBkgTrial[theCase]->SetMarkerStyle(7);
250  fHistoBkgInBinEdgesTrial[theCase]->SetMarkerStyle(7);
251  }
252  }
253  }
254  }
255  fNtupleMultiTrials = new TNtuple(Form("ntuMultiTrial%s",fSuffix.Data()),Form("ntuMultiTrial%s",fSuffix.Data()),"rebin:firstb:minfit:maxfit:bkgfunc:sigfunc:confsig:confmean:chi2:signif:mean:emean:sigma:esigma:rawy:erawy",128000);
256  fNtupleMultiTrials->SetDirectory(nullptr);
257  fNtupleBinCount = new TNtuple(Form("ntuBinCount%s",fSuffix.Data()),Form("ntuBinCount%s",fSuffix.Data()),"rebin:firstb:minfit:maxfit:bkgfunc:sigfunc:confsig:confmean:chi2:nSigmaBC:rawyBC0:erawyBC0:rawyBC1:erawyBC1",128000);
258  fNtupleBinCount->SetDirectory(nullptr);
259  return kTRUE;
260 
261 }
262 
263 //________________________________________________________________________
265  // perform the multiple fits
266 
267  Bool_t hOK=CreateHistos();
268  if(!hOK) return kFALSE;
269 
270  Int_t itrial=0;
271  Int_t itrialBC=0;
273 
274  fMinYieldGlob=999999.;
275  fMaxYieldGlob=0.;
276  Float_t xnt[16];
277  Float_t xntBC[14];
278 
279  for(Int_t ir=0; ir<fNumOfRebinSteps; ir++){
280  Int_t rebin=fRebinSteps[ir];
281  for(Int_t iFirstBin=1; iFirstBin<=fNumOfFirstBinSteps; iFirstBin++) {
282  TH1F* hRebinned=0x0;
283  if(fNumOfFirstBinSteps==1) hRebinned=(TH1F*)AliVertexingHFUtils::RebinHisto(hInvMassHisto,rebin,-1);
284  else hRebinned=(TH1F*)AliVertexingHFUtils::RebinHisto(hInvMassHisto,rebin,iFirstBin);
285  for(Int_t iMinMass=0; iMinMass<fNumOfLowLimFitSteps; iMinMass++){
287  Double_t hmin=TMath::Max(minMassForFit,hRebinned->GetBinLowEdge(2));
288  for(Int_t iMaxMass=0; iMaxMass<fNumOfUpLimFitSteps; iMaxMass++){
290  Double_t hmax=TMath::Min(maxMassForFit,hRebinned->GetBinLowEdge(hRebinned->GetNbinsX()));
291  ++itrial;
292  for(Int_t typeb=0; typeb<kNBkgFuncCases; typeb++){
293  if(typeb==kExpoBkg && !fUseExpoBkg) continue;
294  if(typeb==kLinBkg && !fUseLinBkg) continue;
295  if(typeb==kPol2Bkg && !fUsePol2Bkg) continue;
296  if(typeb==kPol3Bkg && !fUsePol3Bkg) continue;
297  if(typeb==kPol4Bkg && !fUsePol4Bkg) continue;
298  if(typeb==kPol5Bkg && !fUsePol5Bkg) continue;
299  if(typeb==kPowBkg && !fUsePowLawBkg) continue;
300  if(typeb==kPowTimesExpoBkg && !fUsePowLawTimesExpoBkg) continue;
301  for(Int_t types=0; types<kNSigFuncCases; types++){
302  if(types==k2Gaus && !fUse2GausSignal) continue;
304  for(Int_t igs=0; igs<kNFitConfCases; igs++){
305  if (igs==kFixSigUpFreeMean && !fUseFixSigUpFreeMean) continue;
306  if (igs==kFixSigDownFreeMean && !fUseFixSigDownFreeMean) continue;
307  if (igs==kFreeSigFixMean && !fUseFixedMeanFreeS) continue;
308  if (igs==kFreeSigFreeMean && !fUseFreeS) continue;
309  if (igs==kFixSigFreeMean && !fUseFixSigFreeMean) continue;
310  if (igs==kFixSigFixMean && !fUseFixSigFixMean) continue;
311  Int_t theCase=igs*kNBkgFuncCases*kNSigFuncCases+types*kNBkgFuncCases+typeb;
312  Int_t globBin=itrial+theCase*totTrials;
313  for(Int_t j=0; j<16; j++) xnt[j]=0.;
314 
315  Bool_t mustDeleteFitter = kTRUE;
317  if(typeb==kExpoBkg){
318  fitter=new AliHFInvMassFitter(hRebinned, hmin, hmax, AliHFInvMassFitter::kExpo, types);
319  }else if(typeb==kLinBkg){
320  fitter=new AliHFInvMassFitter(hRebinned, hmin, hmax, AliHFInvMassFitter::kLin, types);
321  }else if(typeb==kPol2Bkg){
322  fitter=new AliHFInvMassFitter(hRebinned, hmin, hmax, AliHFInvMassFitter::kPol2, types);
323  }else if(typeb==kPowBkg){
324  fitter=new AliHFInvMassFitter(hRebinned, hmin, hmax, AliHFInvMassFitter::kPow, types);
325  }else if(typeb==kPowTimesExpoBkg){
326  fitter=new AliHFInvMassFitter(hRebinned, hmin, hmax, AliHFInvMassFitter::kPowEx, types);
327  }else{
328  fitter=new AliHFInvMassFitter(hRebinned, hmin, hmax, 6, types);
329  if(typeb==kPol3Bkg) fitter->SetPolDegreeForBackgroundFit(3);
330  if(typeb==kPol4Bkg) fitter->SetPolDegreeForBackgroundFit(4);
331  if(typeb==kPol5Bkg) fitter->SetPolDegreeForBackgroundFit(5);
332  }
333  if(types==k2Gaus){
336  }else if(types==k2GausSigmaRatioPar){
339  }
340  // D0 Reflection
341  if(fhTemplRefl && fhTemplSign){
342  TH1F *hReflModif=(TH1F*)AliVertexingHFUtils::AdaptTemplateRangeAndBinning(fhTemplRefl,hRebinned,minMassForFit,maxMassForFit);
343  TH1F *hSigModif=(TH1F*)AliVertexingHFUtils::AdaptTemplateRangeAndBinning(fhTemplSign,hRebinned,minMassForFit,maxMassForFit);
344  TH1F* hrfl=fitter->SetTemplateReflections(hReflModif,"2gaus",minMassForFit,maxMassForFit);
345  if(!hrfl) printf("ERROR in SetTemplateReflections\n");
346  if(fFixRefloS>0){
347  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)));
348  fitter->SetFixReflOverS(fixSoverRefAt);
349  }
350  delete hReflModif;
351  delete hSigModif;
352  }
353  if(fUseSecondPeak){
355  }
356  if(fFitOption==1) fitter->SetUseChi2Fit();
359  xnt[0]=rebin;
360  xnt[1]=iFirstBin;
361  xnt[2]=minMassForFit;
362  xnt[3]=maxMassForFit;
363  xnt[4]=typeb;
364  xnt[5]=types;
365  xnt[7]=0;
366  if(igs==kFixSigFreeMean){
368  xnt[6]=1;
369  }else if(igs==kFixSigUpFreeMean){
371  xnt[6]=2;
372  }else if(igs==kFixSigDownFreeMean){
374  xnt[6]=3;
375  }else if(igs==kFreeSigFreeMean){
376  xnt[6]=0;
377  }else if(igs==kFixSigFixMean){
379  fitter->SetFixGaussianMean(fMassD);
380  xnt[6]=1;
381  xnt[7]=1;
382  }else if(igs==kFreeSigFixMean){
383  fitter->SetFixGaussianMean(fMassD);
384  xnt[6]=0;
385  xnt[7]=1;
386  }
387  Bool_t out=kFALSE;
388  Double_t chisq=-1.;
389  Double_t sigma=0.;
390  Double_t esigma=0.;
391  Double_t pos=.0;
392  Double_t epos=.0;
393  Double_t ry=.0;
394  Double_t ery=.0;
395  Double_t significance=0.;
396  Double_t erSignif=0.;
397  Double_t bkg=0.;
398  Double_t erbkg=0.;
399  Double_t bkgBEdge=0;
400  Double_t erbkgBEdge=0;
401  if(typeb<kNBkgFuncCases){
402  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);
403  out=fitter->MassFitter(0);
404  chisq=fitter->GetReducedChiSquare();
405  fitter->Significance(fnSigmaForBkgEval,significance,erSignif);
406  sigma=fitter->GetSigma();
407  pos=fitter->GetMean();
408  esigma=fitter->GetSigmaUncertainty();
409  if(esigma<0.00001) esigma=0.0001;
410  epos=fitter->GetMeanUncertainty();
411  if(epos<0.00001) epos=0.0001;
412  ry=fitter->GetRawYield();
413  ery=fitter->GetRawYieldError();
414  fitter->Background(fnSigmaForBkgEval,bkg,erbkg);
415  Double_t minval = hInvMassHisto->GetXaxis()->GetBinLowEdge(hInvMassHisto->FindBin(pos-fnSigmaForBkgEval*sigma));
416  Double_t maxval = hInvMassHisto->GetXaxis()->GetBinUpEdge(hInvMassHisto->FindBin(pos+fnSigmaForBkgEval*sigma));
417  fitter->Background(minval,maxval,bkgBEdge,erbkgBEdge);
418  if(out && fDrawIndividualFits && thePad){
419  thePad->Clear();
420  fitter->DrawHere(thePad, fnSigmaForBkgEval);
421  fMassFitters.push_back(fitter);
422  mustDeleteFitter = kFALSE;
423  for (auto format : fInvMassFitSaveAsFormats) {
424  thePad->SaveAs(Form("FitOutput_%s_Trial%d.%s",hInvMassHisto->GetName(),globBin, format.c_str()));
425  }
426  }
427  }
428  xnt[8]=chisq;
429  if(out && chisq>0. && sigma>0.5*fSigmaGausMC && sigma<2.0*fSigmaGausMC){
430  xnt[9]=significance;
431  xnt[10]=pos;
432  xnt[11]=epos;
433  xnt[12]=sigma;
434  xnt[13]=esigma;
435  xnt[14]=ry;
436  xnt[15]=ery;
437  fHistoRawYieldDistAll->Fill(ry);
438  fHistoRawYieldTrialAll->SetBinContent(globBin,ry);
439  fHistoRawYieldTrialAll->SetBinError(globBin,ery);
440  fHistoSigmaTrialAll->SetBinContent(globBin,sigma);
441  fHistoSigmaTrialAll->SetBinError(globBin,esigma);
442  fHistoMeanTrialAll->SetBinContent(globBin,pos);
443  fHistoMeanTrialAll->SetBinError(globBin,epos);
444  fHistoChi2TrialAll->SetBinContent(globBin,chisq);
445  fHistoChi2TrialAll->SetBinError(globBin,0.00001);
446  fHistoSignifTrialAll->SetBinContent(globBin,significance);
447  fHistoSignifTrialAll->SetBinError(globBin,erSignif);
448  if(fSaveBkgVal) {
449  fHistoBkgTrialAll->SetBinContent(globBin,bkg);
450  fHistoBkgTrialAll->SetBinError(globBin,erbkg);
451  fHistoBkgInBinEdgesTrialAll->SetBinContent(globBin,bkgBEdge);
452  fHistoBkgInBinEdgesTrialAll->SetBinError(globBin,erbkgBEdge);
453  }
454 
455  if(ry<fMinYieldGlob) fMinYieldGlob=ry;
456  if(ry>fMaxYieldGlob) fMaxYieldGlob=ry;
457  fHistoRawYieldDist[theCase]->Fill(ry);
458  fHistoRawYieldTrial[theCase]->SetBinContent(itrial,ry);
459  fHistoRawYieldTrial[theCase]->SetBinError(itrial,ery);
460  fHistoSigmaTrial[theCase]->SetBinContent(itrial,sigma);
461  fHistoSigmaTrial[theCase]->SetBinError(itrial,esigma);
462  fHistoMeanTrial[theCase]->SetBinContent(itrial,pos);
463  fHistoMeanTrial[theCase]->SetBinError(itrial,epos);
464  fHistoChi2Trial[theCase]->SetBinContent(itrial,chisq);
465  fHistoChi2Trial[theCase]->SetBinError(itrial,0.00001);
466  fHistoSignifTrial[theCase]->SetBinContent(itrial,significance);
467  fHistoSignifTrial[theCase]->SetBinError(itrial,erSignif);
468  if(fSaveBkgVal) {
469  fHistoBkgTrial[theCase]->SetBinContent(itrial,bkg);
470  fHistoBkgTrial[theCase]->SetBinError(itrial,erbkg);
471  fHistoBkgInBinEdgesTrial[theCase]->SetBinContent(itrial,bkgBEdge);
472  fHistoBkgInBinEdgesTrial[theCase]->SetBinError(itrial,erbkgBEdge);
473  }
474  fNtupleMultiTrials->Fill(xnt);
475  if(types==0){
476  // bin counting done only for 1 case of signal line shape
477  for(Int_t j=0; j<9; j++) xntBC[j]=xnt[j];
478 
479 
480  for(Int_t iStepBC=0; iStepBC<fNumOfnSigmaBinCSteps; iStepBC++){
481  Double_t minMassBC=fMassD-fnSigmaBinCSteps[iStepBC]*sigma;
482  Double_t maxMassBC=fMassD+fnSigmaBinCSteps[iStepBC]*sigma;
483  if(minMassBC>minMassForFit &&
484  maxMassBC<maxMassForFit &&
485  minMassBC>(hRebinned->GetXaxis()->GetXmin()) &&
486  maxMassBC<(hRebinned->GetXaxis()->GetXmax())){
487  Double_t cnts0,ecnts0;
488  Double_t cnts1,ecnts1;
489  cnts0=fitter->GetRawYieldBinCounting(ecnts0,fnSigmaBinCSteps[iStepBC],0,0);
490  cnts1=fitter->GetRawYieldBinCounting(ecnts1,fnSigmaBinCSteps[iStepBC],1,0);
491  xntBC[9]=fnSigmaBinCSteps[iStepBC];
492  xntBC[10]=cnts0;
493  xntBC[11]=ecnts0;
494  xntBC[12]=cnts1;
495  xntBC[13]=ecnts1;
496  ++itrialBC;
497  fHistoRawYieldDistBinC0All->Fill(cnts0);
498  fHistoRawYieldTrialBinC0All->SetBinContent(globBin,iStepBC+1,cnts0);
499  fHistoRawYieldTrialBinC0All->SetBinError(globBin,iStepBC+1,ecnts0);
500  fHistoRawYieldTrialBinC0[theCase]->SetBinContent(itrial,iStepBC+1,cnts0);
501  fHistoRawYieldTrialBinC0[theCase]->SetBinError(itrial,iStepBC+1,ecnts0);
502  fHistoRawYieldDistBinC0[theCase]->Fill(cnts0);
503  fHistoRawYieldDistBinC1All->Fill(cnts1);
504  fHistoRawYieldTrialBinC1All->SetBinContent(globBin,iStepBC+1,cnts1);
505  fHistoRawYieldTrialBinC1All->SetBinError(globBin,iStepBC+1,ecnts1);
506  fHistoRawYieldTrialBinC1[theCase]->SetBinContent(itrial,iStepBC+1,cnts1);
507  fHistoRawYieldTrialBinC1[theCase]->SetBinError(itrial,iStepBC+1,ecnts1);
508  fHistoRawYieldDistBinC1[theCase]->Fill(cnts1);
509  fNtupleBinCount->Fill(xntBC);
510  }
511  }
512  }
513  }
514  if (mustDeleteFitter) delete fitter;
515  }
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
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
Bool_t fUse2GausSignal
switch for power law background
AliHFInvMassFitter class for the fit of invariant mass distribution of charm hadrons.
void SetFixSecondGaussianSigma(Double_t sigma)
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
void SetFixFrac2Gaus(Double_t frac)
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)
Double_t fFixSecondGausSigRat
value to fix 2nd gaus area
Bool_t fUseFixSigUpFreeMean
value to fix ratio os sigmas
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
void SetFixRatio2GausSigma(Double_t sigmafrac)
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
Double_t fFixSecondGausFrac
value to fix 2nd gaus sigma
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 fUse2GausSigmaRatioSignal
swicth for 2 gaus line shape for S
Bool_t fUseFixedMeanFreeS
switch for FreeSigma
Double_t GetSigmaUncertainty() const
void Significance(Double_t nOfSigma, Double_t &significance, Double_t &errsignificance) const
Double_t fFixSecondGausSig
swicth for 2 gaus line shape for S
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