AliPhysics  9b6b435 (9b6b435)
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  fSigmaMCVariationUp(0.1),
49  fSigmaMCVariationDw(0.1),
50  fMassD(1.86484),
51  fSuffix(""),
52  fFitOption(0),
53  fUseExpoBkg(kTRUE),
54  fUseLinBkg(kTRUE),
55  fUsePol2Bkg(kTRUE),
56  fUsePol3Bkg(kTRUE),
57  fUsePol4Bkg(kTRUE),
58  fUsePol5Bkg(kFALSE),
59  fUsePowLawBkg(kFALSE),
60  fUsePowLawTimesExpoBkg(kFALSE),
61  fUse2GausSignal(kFALSE),
62  fUse2GausSigmaRatioSignal(kFALSE),
63  fFixSecondGausSig(-1.),
64  fFixSecondGausFrac(-1.),
65  fFixSecondGausSigRat(-1.),
66  fUseFixSigUpFreeMean(kTRUE),
67  fUseFixSigDownFreeMean(kTRUE),
68  fUseFreeS(kTRUE),
69  fUseFixedMeanFreeS(kTRUE),
70  fUseFixSigFreeMean(kTRUE),
71  fUseFixSigFixMean(kTRUE),
72  fUseSecondPeak(kFALSE),
73  fMassSecondPeak(1.86958),
74  fSigmaSecondPeak(0.01),
75  fFixMassSecondPeak(kFALSE),
76  fFixSigmaSecondPeak(kFALSE),
77  fSaveBkgVal(kFALSE),
78  fDrawIndividualFits(kFALSE),
79  fHistoRawYieldDistAll(0x0),
80  fHistoRawYieldTrialAll(0x0),
81  fHistoSigmaTrialAll(0x0),
82  fHistoMeanTrialAll(0x0),
83  fHistoChi2TrialAll(0x0),
84  fHistoSignifTrialAll(0x0),
85  fHistoBkgTrialAll(0x0),
86  fHistoBkgInBinEdgesTrialAll(0x0),
87  fHistoRawYieldDistBinC0All(0x0),
88  fHistoRawYieldTrialBinC0All(0x0),
89  fHistoRawYieldDistBinC1All(0x0),
90  fHistoRawYieldTrialBinC1All(0x0),
91  fHistoRawYieldDist(0x0),
92  fHistoRawYieldTrial(0x0),
93  fHistoSigmaTrial(0x0),
94  fHistoMeanTrial(0x0),
95  fHistoChi2Trial(0x0),
96  fHistoSignifTrial(0x0),
97  fHistoBkgTrial(0x0),
98  fHistoBkgInBinEdgesTrial(0x0),
99  fHistoRawYieldDistBinC0(0x0),
100  fHistoRawYieldTrialBinC0(0x0),
101  fHistoRawYieldDistBinC1(0x0),
102  fHistoRawYieldTrialBinC1(0x0),
103  fhTemplRefl(0x0),
104  fhTemplSign(0x0),
105  fFixRefloS(1.),
106  fNtupleMultiTrials(0x0),
107  fNtupleBinCount(0x0),
108  fMinYieldGlob(0),
109  fMaxYieldGlob(0),
110  fMassFitters()
111 {
112  // constructor
113  Int_t rebinStep[4]={3,4,5,6};
114  Double_t minMassStep[6]={1.68,1.70,1.72,1.74,1.76,1.78};
115  Double_t maxMassStep[6]={2.06,2.04,2.02,2.00,1.98,1.96};
116  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};
117  ConfigureRebinSteps(4,rebinStep);
118  ConfigureLowLimFitSteps(6,minMassStep);
119  ConfigureUpLimFitSteps(6,maxMassStep);
120  ConfigurenSigmaBinCSteps(11,nSigmasBC);
121 }
122 
123 //________________________________________________________________________
125  // destructor
126  delete [] fRebinSteps;
127  delete [] fLowLimFitSteps;
128  delete [] fUpLimFitSteps;
129  if(fhTemplRefl) delete fhTemplRefl;
130  if(fhTemplSign) delete fhTemplSign;
131  for (auto fitter : fMassFitters) delete fitter;
132  delete fHistoRawYieldDistAll;
133  delete fHistoRawYieldTrialAll;
134  delete fHistoSigmaTrialAll;
135  delete fHistoMeanTrialAll;
136  delete fHistoChi2TrialAll;
137  delete fHistoSignifTrialAll;
138  delete fHistoBkgTrialAll;
144  for(Int_t typeb=0; typeb<kNBkgFuncCases; typeb++){
145  for(Int_t types=0; types<kNSigFuncCases; types++){
146  for(Int_t igs=0; igs<kNFitConfCases; igs++){
147  Int_t theCase=igs*kNBkgFuncCases*kNSigFuncCases+types*kNBkgFuncCases+typeb;
148  delete fHistoRawYieldDist[theCase];
149  delete fHistoRawYieldDistBinC0[theCase];
150  delete fHistoRawYieldDistBinC1[theCase];
151  delete fHistoRawYieldTrial[theCase];
152  delete fHistoRawYieldTrialBinC0[theCase];
153  delete fHistoRawYieldTrialBinC1[theCase];
154  delete fHistoSigmaTrial[theCase];
155  delete fHistoMeanTrial[theCase];
156  delete fHistoChi2Trial[theCase];
157  delete fHistoSignifTrial[theCase];
158  if(fHistoBkgTrial) delete fHistoBkgTrial[theCase];
160  }
161  }
162  }
163  delete [] fHistoRawYieldDist;
164  delete [] fHistoRawYieldDistBinC0;
165  delete [] fHistoRawYieldDistBinC1;
166  delete [] fHistoRawYieldTrial;
167  delete [] fHistoRawYieldTrialBinC0;
168  delete [] fHistoRawYieldTrialBinC1;
169  delete [] fHistoSigmaTrial;
170  delete [] fHistoMeanTrial;
171  delete [] fHistoChi2Trial;
172  delete [] fHistoSignifTrial;
173  delete [] fHistoBkgTrial;
174  delete [] fHistoBkgInBinEdgesTrial;
175 
176  delete fNtupleMultiTrials;
177  delete fNtupleBinCount;
178 }
179 
180 //________________________________________________________________________
182  // creates output histograms
183 
185 
186  TString funcBkg[kNBkgFuncCases]={"Expo","Lin","Pol2","Pol3","Pol4","Pol5","PowLaw","PowLawExpo"};
187  TString funcSig[kNSigFuncCases]={"","2Gaus","2GausSigmaRatio"};
188  TString gausSig[kNFitConfCases]={"FixedS","FixedSp20","FixedSm20","FreeS","FixedMeanFixedS","FixedMeanFreeS"};
189 
191 
192  fHistoRawYieldDistAll = new TH1F(Form("hRawYieldDistAll%s",fSuffix.Data())," ; Raw Yield",5000,0.,50000.);
193  fHistoRawYieldTrialAll = new TH1F(Form("hRawYieldTrialAll%s",fSuffix.Data())," ; Trial # ; Raw Yield",nCases*totTrials,-0.5,nCases*totTrials-0.5);
194  fHistoSigmaTrialAll = new TH1F(Form("hSigmaTrialAll%s",fSuffix.Data())," ; Trial # ; Sigma (GeV/c^{2})",nCases*totTrials,-0.5,nCases*totTrials-0.5);
195  fHistoMeanTrialAll = new TH1F(Form("hMeanTrialAll%s",fSuffix.Data())," ; Trial # ; Mean (GeV/c^{2})",nCases*totTrials,-0.5,nCases*totTrials-0.5);
196  fHistoChi2TrialAll = new TH1F(Form("hChi2TrialAll%s",fSuffix.Data())," ; Trial # ; #chi^{2}",nCases*totTrials,-0.5,nCases*totTrials-0.5);
197  fHistoSignifTrialAll = new TH1F(Form("hSignifTrialAll%s",fSuffix.Data())," ; Trial # ; Significance",nCases*totTrials,-0.5,nCases*totTrials-0.5);
198  if(fSaveBkgVal) {
199  fHistoBkgTrialAll = new TH1F(Form("hBkgTrialAll%s",fSuffix.Data())," ; Background",nCases*totTrials,-0.5,nCases*totTrials-0.5);
200  fHistoBkgInBinEdgesTrialAll = new TH1F(Form("hBkgInBinEdgesTrialAll%s",fSuffix.Data())," ; Background in bin edges",nCases*totTrials,-0.5,nCases*totTrials-0.5);
201  }
202 
203 
204  fHistoRawYieldDistBinC0All = new TH1F(Form("hRawYieldDistBinC0All%s",fSuffix.Data())," ; Raw Yield (bin count)",5000,0.,50000.);
205  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);
206  fHistoRawYieldDistBinC1All = new TH1F(Form("hRawYieldDistBinC1All%s",fSuffix.Data())," ; Raw Yield (bin count)",5000,0.,50000.);
207  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);
208 
209  fHistoRawYieldDist = new TH1F*[nCases];
210  fHistoRawYieldTrial = new TH1F*[nCases];
211  fHistoSigmaTrial = new TH1F*[nCases];
212  fHistoMeanTrial = new TH1F*[nCases];
213  fHistoChi2Trial = new TH1F*[nCases];
214  fHistoSignifTrial = new TH1F*[nCases];
215  if(fSaveBkgVal) {
216  fHistoBkgTrial = new TH1F*[nCases];
217  fHistoBkgInBinEdgesTrial = new TH1F*[nCases];
218  }
219 
220  fHistoRawYieldDistBinC0 = new TH1F*[nCases];
221  fHistoRawYieldTrialBinC0 = new TH2F*[nCases];
222  fHistoRawYieldDistBinC1 = new TH1F*[nCases];
223  fHistoRawYieldTrialBinC1 = new TH2F*[nCases];
224 
225  for(Int_t typeb=0; typeb<kNBkgFuncCases; typeb++){
226  for(Int_t types=0; types<kNSigFuncCases; types++){
227  for(Int_t igs=0; igs<kNFitConfCases; igs++){
228  Int_t theCase=igs*kNBkgFuncCases*kNSigFuncCases+types*kNBkgFuncCases+typeb;
229  printf("typeb= %d types = %d igs = %d theCase=%d out of %d\n",typeb,types,igs,theCase,nCases);
230  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.);
231  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.);
232  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.);
233  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);
234  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);
235  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);
236  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);
237  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);
238  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);
239  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);
240  if(fSaveBkgVal) {
241  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);
242  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);
243  }
244 
245  fHistoChi2Trial[theCase]->SetMarkerStyle(7);
246  fHistoSignifTrial[theCase]->SetMarkerStyle(7);
247  fHistoSigmaTrial[theCase]->SetMarkerStyle(7);
248  fHistoMeanTrial[theCase]->SetMarkerStyle(7);
249  if(fSaveBkgVal) {
250  fHistoBkgTrial[theCase]->SetMarkerStyle(7);
251  fHistoBkgInBinEdgesTrial[theCase]->SetMarkerStyle(7);
252  }
253  }
254  }
255  }
256  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);
257  fNtupleMultiTrials->SetDirectory(nullptr);
258  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);
259  fNtupleBinCount->SetDirectory(nullptr);
260  return kTRUE;
261 
262 }
263 
264 //________________________________________________________________________
266  // perform the multiple fits
267 
268  Bool_t hOK=CreateHistos();
269  if(!hOK) return kFALSE;
270 
271  Int_t itrial=0;
272  Int_t itrialBC=0;
274 
275  fMinYieldGlob=999999.;
276  fMaxYieldGlob=0.;
277  Float_t xnt[16];
278  Float_t xntBC[14];
279 
280  for(Int_t ir=0; ir<fNumOfRebinSteps; ir++){
281  Int_t rebin=fRebinSteps[ir];
282  for(Int_t iFirstBin=1; iFirstBin<=fNumOfFirstBinSteps; iFirstBin++) {
283  TH1F* hRebinned=0x0;
284  if(fNumOfFirstBinSteps==1) hRebinned=(TH1F*)AliVertexingHFUtils::RebinHisto(hInvMassHisto,rebin,-1);
285  else hRebinned=(TH1F*)AliVertexingHFUtils::RebinHisto(hInvMassHisto,rebin,iFirstBin);
286  for(Int_t iMinMass=0; iMinMass<fNumOfLowLimFitSteps; iMinMass++){
288  Double_t hmin=TMath::Max(minMassForFit,hRebinned->GetBinLowEdge(2));
289  for(Int_t iMaxMass=0; iMaxMass<fNumOfUpLimFitSteps; iMaxMass++){
291  Double_t hmax=TMath::Min(maxMassForFit,hRebinned->GetBinLowEdge(hRebinned->GetNbinsX()));
292  ++itrial;
293  for(Int_t typeb=0; typeb<kNBkgFuncCases; typeb++){
294  if(typeb==kExpoBkg && !fUseExpoBkg) continue;
295  if(typeb==kLinBkg && !fUseLinBkg) continue;
296  if(typeb==kPol2Bkg && !fUsePol2Bkg) continue;
297  if(typeb==kPol3Bkg && !fUsePol3Bkg) continue;
298  if(typeb==kPol4Bkg && !fUsePol4Bkg) continue;
299  if(typeb==kPol5Bkg && !fUsePol5Bkg) continue;
300  if(typeb==kPowBkg && !fUsePowLawBkg) continue;
301  if(typeb==kPowTimesExpoBkg && !fUsePowLawTimesExpoBkg) continue;
302  for(Int_t types=0; types<kNSigFuncCases; types++){
303  if(types==k2Gaus && !fUse2GausSignal) continue;
305  for(Int_t igs=0; igs<kNFitConfCases; igs++){
306  if (igs==kFixSigUpFreeMean && !fUseFixSigUpFreeMean) continue;
307  if (igs==kFixSigDownFreeMean && !fUseFixSigDownFreeMean) continue;
308  if (igs==kFreeSigFixMean && !fUseFixedMeanFreeS) continue;
309  if (igs==kFreeSigFreeMean && !fUseFreeS) continue;
310  if (igs==kFixSigFreeMean && !fUseFixSigFreeMean) continue;
311  if (igs==kFixSigFixMean && !fUseFixSigFixMean) continue;
312  Int_t theCase=igs*kNBkgFuncCases*kNSigFuncCases+types*kNBkgFuncCases+typeb;
313  Int_t globBin=itrial+theCase*totTrials;
314  for(Int_t j=0; j<16; j++) xnt[j]=0.;
315 
316  Bool_t mustDeleteFitter = kTRUE;
318  if(typeb==kExpoBkg){
319  fitter=new AliHFInvMassFitter(hRebinned, hmin, hmax, AliHFInvMassFitter::kExpo, types);
320  }else if(typeb==kLinBkg){
321  fitter=new AliHFInvMassFitter(hRebinned, hmin, hmax, AliHFInvMassFitter::kLin, types);
322  }else if(typeb==kPol2Bkg){
323  fitter=new AliHFInvMassFitter(hRebinned, hmin, hmax, AliHFInvMassFitter::kPol2, types);
324  }else if(typeb==kPowBkg){
325  fitter=new AliHFInvMassFitter(hRebinned, hmin, hmax, AliHFInvMassFitter::kPow, types);
326  }else if(typeb==kPowTimesExpoBkg){
327  fitter=new AliHFInvMassFitter(hRebinned, hmin, hmax, AliHFInvMassFitter::kPowEx, types);
328  }else{
329  fitter=new AliHFInvMassFitter(hRebinned, hmin, hmax, 6, types);
330  if(typeb==kPol3Bkg) fitter->SetPolDegreeForBackgroundFit(3);
331  if(typeb==kPol4Bkg) fitter->SetPolDegreeForBackgroundFit(4);
332  if(typeb==kPol5Bkg) fitter->SetPolDegreeForBackgroundFit(5);
333  }
334  if(types==k2Gaus){
337  }else if(types==k2GausSigmaRatioPar){
340  }
341  // D0 Reflection
342  if(fhTemplRefl && fhTemplSign){
343  TH1F *hReflModif=(TH1F*)AliVertexingHFUtils::AdaptTemplateRangeAndBinning(fhTemplRefl,hRebinned,minMassForFit,maxMassForFit);
344  TH1F *hSigModif=(TH1F*)AliVertexingHFUtils::AdaptTemplateRangeAndBinning(fhTemplSign,hRebinned,minMassForFit,maxMassForFit);
345  TH1F* hrfl=fitter->SetTemplateReflections(hReflModif,"2gaus",minMassForFit,maxMassForFit);
346  if(!hrfl) printf("ERROR in SetTemplateReflections\n");
347  if(fFixRefloS>0){
348  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)));
349  fitter->SetFixReflOverS(fixSoverRefAt);
350  }
351  delete hReflModif;
352  delete hSigModif;
353  }
354  if(fUseSecondPeak){
356  }
357  if(fFitOption==1) fitter->SetUseChi2Fit();
360  xnt[0]=rebin;
361  xnt[1]=iFirstBin;
362  xnt[2]=minMassForFit;
363  xnt[3]=maxMassForFit;
364  xnt[4]=typeb;
365  xnt[5]=types;
366  xnt[7]=0;
367  if(igs==kFixSigFreeMean){
369  xnt[6]=1;
370  }else if(igs==kFixSigUpFreeMean){
372  xnt[6]=2;
373  }else if(igs==kFixSigDownFreeMean){
375  xnt[6]=3;
376  }else if(igs==kFreeSigFreeMean){
377  xnt[6]=0;
378  }else if(igs==kFixSigFixMean){
380  fitter->SetFixGaussianMean(fMassD);
381  xnt[6]=1;
382  xnt[7]=1;
383  }else if(igs==kFreeSigFixMean){
384  fitter->SetFixGaussianMean(fMassD);
385  xnt[6]=0;
386  xnt[7]=1;
387  }
388  Bool_t out=kFALSE;
389  Double_t chisq=-1.;
390  Double_t sigma=0.;
391  Double_t esigma=0.;
392  Double_t pos=.0;
393  Double_t epos=.0;
394  Double_t ry=.0;
395  Double_t ery=.0;
396  Double_t significance=0.;
397  Double_t erSignif=0.;
398  Double_t bkg=0.;
399  Double_t erbkg=0.;
400  Double_t bkgBEdge=0;
401  Double_t erbkgBEdge=0;
402  if(typeb<kNBkgFuncCases){
403  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);
404  out=fitter->MassFitter(0);
405  chisq=fitter->GetReducedChiSquare();
406  fitter->Significance(fnSigmaForBkgEval,significance,erSignif);
407  sigma=fitter->GetSigma();
408  pos=fitter->GetMean();
409  esigma=fitter->GetSigmaUncertainty();
410  if(esigma<0.00001) esigma=0.0001;
411  epos=fitter->GetMeanUncertainty();
412  if(epos<0.00001) epos=0.0001;
413  ry=fitter->GetRawYield();
414  ery=fitter->GetRawYieldError();
415  fitter->Background(fnSigmaForBkgEval,bkg,erbkg);
416  Double_t minval = hInvMassHisto->GetXaxis()->GetBinLowEdge(hInvMassHisto->FindBin(pos-fnSigmaForBkgEval*sigma));
417  Double_t maxval = hInvMassHisto->GetXaxis()->GetBinUpEdge(hInvMassHisto->FindBin(pos+fnSigmaForBkgEval*sigma));
418  fitter->Background(minval,maxval,bkgBEdge,erbkgBEdge);
419  if(out && fDrawIndividualFits && thePad){
420  thePad->Clear();
421  fitter->DrawHere(thePad, fnSigmaForBkgEval);
422  fMassFitters.push_back(fitter);
423  mustDeleteFitter = kFALSE;
424  for (auto format : fInvMassFitSaveAsFormats) {
425  thePad->SaveAs(Form("FitOutput_%s_Trial%d.%s",hInvMassHisto->GetName(),globBin, format.c_str()));
426  }
427  }
428  }
429  xnt[8]=chisq;
430  if(out && chisq>0. && sigma>0.5*fSigmaGausMC && sigma<2.0*fSigmaGausMC){
431  xnt[9]=significance;
432  xnt[10]=pos;
433  xnt[11]=epos;
434  xnt[12]=sigma;
435  xnt[13]=esigma;
436  xnt[14]=ry;
437  xnt[15]=ery;
438  fHistoRawYieldDistAll->Fill(ry);
439  fHistoRawYieldTrialAll->SetBinContent(globBin,ry);
440  fHistoRawYieldTrialAll->SetBinError(globBin,ery);
441  fHistoSigmaTrialAll->SetBinContent(globBin,sigma);
442  fHistoSigmaTrialAll->SetBinError(globBin,esigma);
443  fHistoMeanTrialAll->SetBinContent(globBin,pos);
444  fHistoMeanTrialAll->SetBinError(globBin,epos);
445  fHistoChi2TrialAll->SetBinContent(globBin,chisq);
446  fHistoChi2TrialAll->SetBinError(globBin,0.00001);
447  fHistoSignifTrialAll->SetBinContent(globBin,significance);
448  fHistoSignifTrialAll->SetBinError(globBin,erSignif);
449  if(fSaveBkgVal) {
450  fHistoBkgTrialAll->SetBinContent(globBin,bkg);
451  fHistoBkgTrialAll->SetBinError(globBin,erbkg);
452  fHistoBkgInBinEdgesTrialAll->SetBinContent(globBin,bkgBEdge);
453  fHistoBkgInBinEdgesTrialAll->SetBinError(globBin,erbkgBEdge);
454  }
455 
456  if(ry<fMinYieldGlob) fMinYieldGlob=ry;
457  if(ry>fMaxYieldGlob) fMaxYieldGlob=ry;
458  fHistoRawYieldDist[theCase]->Fill(ry);
459  fHistoRawYieldTrial[theCase]->SetBinContent(itrial,ry);
460  fHistoRawYieldTrial[theCase]->SetBinError(itrial,ery);
461  fHistoSigmaTrial[theCase]->SetBinContent(itrial,sigma);
462  fHistoSigmaTrial[theCase]->SetBinError(itrial,esigma);
463  fHistoMeanTrial[theCase]->SetBinContent(itrial,pos);
464  fHistoMeanTrial[theCase]->SetBinError(itrial,epos);
465  fHistoChi2Trial[theCase]->SetBinContent(itrial,chisq);
466  fHistoChi2Trial[theCase]->SetBinError(itrial,0.00001);
467  fHistoSignifTrial[theCase]->SetBinContent(itrial,significance);
468  fHistoSignifTrial[theCase]->SetBinError(itrial,erSignif);
469  if(fSaveBkgVal) {
470  fHistoBkgTrial[theCase]->SetBinContent(itrial,bkg);
471  fHistoBkgTrial[theCase]->SetBinError(itrial,erbkg);
472  fHistoBkgInBinEdgesTrial[theCase]->SetBinContent(itrial,bkgBEdge);
473  fHistoBkgInBinEdgesTrial[theCase]->SetBinError(itrial,erbkgBEdge);
474  }
475  fNtupleMultiTrials->Fill(xnt);
476  if(types==0){
477  // bin counting done only for 1 case of signal line shape
478  for(Int_t j=0; j<9; j++) xntBC[j]=xnt[j];
479 
480 
481  for(Int_t iStepBC=0; iStepBC<fNumOfnSigmaBinCSteps; iStepBC++){
482  Double_t minMassBC=fMassD-fnSigmaBinCSteps[iStepBC]*sigma;
483  Double_t maxMassBC=fMassD+fnSigmaBinCSteps[iStepBC]*sigma;
484  if(minMassBC>minMassForFit &&
485  maxMassBC<maxMassForFit &&
486  minMassBC>(hRebinned->GetXaxis()->GetXmin()) &&
487  maxMassBC<(hRebinned->GetXaxis()->GetXmax())){
488  Double_t cnts0,ecnts0;
489  Double_t cnts1,ecnts1;
490  cnts0=fitter->GetRawYieldBinCounting(ecnts0,fnSigmaBinCSteps[iStepBC],0,0);
491  cnts1=fitter->GetRawYieldBinCounting(ecnts1,fnSigmaBinCSteps[iStepBC],1,0);
492  xntBC[9]=fnSigmaBinCSteps[iStepBC];
493  xntBC[10]=cnts0;
494  xntBC[11]=ecnts0;
495  xntBC[12]=cnts1;
496  xntBC[13]=ecnts1;
497  ++itrialBC;
498  fHistoRawYieldDistBinC0All->Fill(cnts0);
499  fHistoRawYieldTrialBinC0All->SetBinContent(globBin,iStepBC+1,cnts0);
500  fHistoRawYieldTrialBinC0All->SetBinError(globBin,iStepBC+1,ecnts0);
501  fHistoRawYieldTrialBinC0[theCase]->SetBinContent(itrial,iStepBC+1,cnts0);
502  fHistoRawYieldTrialBinC0[theCase]->SetBinError(itrial,iStepBC+1,ecnts0);
503  fHistoRawYieldDistBinC0[theCase]->Fill(cnts0);
504  fHistoRawYieldDistBinC1All->Fill(cnts1);
505  fHistoRawYieldTrialBinC1All->SetBinContent(globBin,iStepBC+1,cnts1);
506  fHistoRawYieldTrialBinC1All->SetBinError(globBin,iStepBC+1,ecnts1);
507  fHistoRawYieldTrialBinC1[theCase]->SetBinContent(itrial,iStepBC+1,cnts1);
508  fHistoRawYieldTrialBinC1[theCase]->SetBinError(itrial,iStepBC+1,ecnts1);
509  fHistoRawYieldDistBinC1[theCase]->Fill(cnts1);
510  fNtupleBinCount->Fill(xntBC);
511  }
512  }
513  }
514  }
515  if (mustDeleteFitter) delete fitter;
516  }
517  }
518  }
519  }
520  }
521  delete hRebinned;
522  }
523  }
524  return kTRUE;
525 }
526 
527 //________________________________________________________________________
529  // save histos in a root file for further analysis
531  TFile outHistos(fileName.Data(),option.Data());
532  if (outHistos.IsZombie()) {
533  Printf("Could not open file '%s'!", fileName.Data());
534  return;
535  }
536  outHistos.cd();
537  fHistoRawYieldTrialAll->Write();
538  fHistoSigmaTrialAll->Write();
539  fHistoMeanTrialAll->Write();
540  fHistoChi2TrialAll->Write();
541  fHistoSignifTrialAll->Write();
542  if(fSaveBkgVal) {
543  fHistoBkgTrialAll->Write();
545  }
550  for(Int_t ic=0; ic<nCases; ic++){
551  fHistoRawYieldTrial[ic]->Write();
552  fHistoSigmaTrial[ic]->Write();
553  fHistoMeanTrial[ic]->Write();
554  fHistoChi2Trial[ic]->Write();
555  fHistoSignifTrial[ic]->Write();
556  if(fSaveBkgVal) {
557  fHistoBkgTrial[ic]->Write();
558  fHistoBkgInBinEdgesTrial[ic]->Write();
559  }
560  fHistoRawYieldTrialBinC0[ic]->Write();
561  fHistoRawYieldDistBinC0[ic]->Write();
562  fHistoRawYieldTrialBinC1[ic]->Write();
563  fHistoRawYieldDistBinC1[ic]->Write();
564  }
565  fNtupleMultiTrials->SetDirectory(&outHistos);
566  fNtupleMultiTrials->Write();
567  fNtupleBinCount->SetDirectory(&outHistos);
568  fNtupleBinCount->Write();
569  outHistos.Close();
570 }
571 
572 //________________________________________________________________________
573 void AliHFInvMassMultiTrialFit::DrawHistos(TCanvas* cry) const{
574  // draw histos
575  cry->Divide(2,2);
576  cry->cd(1);
577  fHistoSigmaTrialAll->Draw();
578  cry->cd(3);
579  fHistoChi2TrialAll->Draw();
580  cry->cd(2);
581  fHistoRawYieldTrialAll->Draw();
582  cry->cd(4);
583  fHistoRawYieldDistAll->Draw();
584  TLatex* tmean=new TLatex(0.15,0.8,Form("mean=%.1f",fHistoRawYieldDistAll->GetMean()));
585  tmean->SetNDC();
586  tmean->Draw();
587  TLatex* thrms=new TLatex(0.15,0.72,Form("rms=%.1f",fHistoRawYieldDistAll->GetRMS()));
588  thrms->SetNDC();
589  thrms->Draw();
590  TLatex* tmax=new TLatex(0.6,0.8,Form("max=%.1f",fMaxYieldGlob));
591  tmax->SetNDC();
592  tmax->Draw();
593  TLatex* tmin=new TLatex(0.6,0.72,Form("min=%.1f",fMinYieldGlob));
594  tmin->SetNDC();
595  tmin->Draw();
596  TLatex* trms=new TLatex(0.6,0.64,Form("(max-min)/sqrt(12)=%.1f",(fMaxYieldGlob-fMinYieldGlob)/sqrt(12)));
597  trms->SetNDC();
598  trms->Draw();
599 
600 }
601 //________________________________________________________________________
603  Int_t iCase){
604  //
605 
606  TH1F *hCutTmp=(TH1F*)histoToFit->Clone("hCutTmp");
607  for(Int_t ib=1; ib<=hCutTmp->GetNbinsX(); ib++){
608  Double_t xc=hCutTmp->GetBinCenter(ib);
609  if(xc>(fMassD-5.*fSigmaGausMC) && xc<(fMassD+5.*fSigmaGausMC)){
610  hCutTmp->SetBinContent(ib,0.);
611  hCutTmp->SetBinError(ib,0.);
612  }
613  }
614 
615  hCutTmp->Fit("pol2","E0","",hmin,hmax);
616  TF1* f2=(TF1*)hCutTmp->GetListOfFunctions()->FindObject("pol2");
617  TF1* f3=new TF1("myPol3","pol3");
618  for(Int_t i=0; i<3;i++) f3->SetParameter(i,f2->GetParameter(i));
619  hCutTmp->Fit(f3,"E0","",hmin,hmax);
620  Double_t quickCount=0.;
621  for(Int_t ib=1; ib<=histoToFit->GetNbinsX(); ib++){
622  Double_t xc=hCutTmp->GetBinCenter(ib);
623  if(xc>(fMassD-3.*fSigmaGausMC) && xc<(fMassD+3.*fSigmaGausMC)){
624  quickCount+=(histoToFit->GetBinContent(ib)-f3->Eval(xc));
625  }
626  }
627  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");
628  fSB->SetParameter(0,quickCount);
629  fSB->SetParameter(1,fMassD);
630  fSB->SetParameter(2,fSigmaGausMC);
631  for(Int_t j=0; j<4; j++) fSB->SetParameter(j+3,f3->GetParameter(j));
632  if(iCase==0) fSB->FixParameter(2,fSigmaGausMC);
633  else if(iCase==1) fSB->FixParameter(2,fSigmaGausMC*(1.+fSigmaMCVariationUp));
634  else if(iCase==2) fSB->FixParameter(2,fSigmaGausMC*(1.-fSigmaMCVariationDw));
635  else if(iCase==4){
636  fSB->FixParameter(1,fMassD);
637  fSB->FixParameter(2,fSigmaGausMC);
638  } else if(iCase==5){
639  fSB->FixParameter(1,fMassD);
640  }
641  histoToFit->Fit(fSB,"ME0","",hmin,hmax);
642  // quality cuts
643  if(fSB->GetParError(0)<0.01*fSB->GetParameter(0)) return kFALSE;
644  if(fSB->GetParError(0)>0.6*fSB->GetParameter(0)) return kFALSE;
645 
646  delete hCutTmp;
647  return kTRUE;
648 }
649 //__________________________________________________________________________________
650 void AliHFInvMassMultiTrialFit::SetTemplatesForReflections(const TH1F *hr, const TH1F *hs) {
652  if(fhTemplSign) delete fhTemplSign;
653  if(fhTemplRefl) delete fhTemplRefl;
654  fhTemplRefl=(TH1F*)hr->Clone("hTemplRefl");
655  fhTemplSign=(TH1F*)hs->Clone("hTemplSign");
656  return;
657 }
658 
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
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)
Double_t fSigmaMCVariationUp
sigma of D meson peak from MC
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
Double_t fSigmaMCVariationDw
relative variation of the sigma
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