AliPhysics  449db5a (449db5a)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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  fMinYieldGlob(0),
102  fMaxYieldGlob(0),
103  fMassFitters()
104 {
105  // constructor
106  Int_t rebinStep[4]={3,4,5,6};
107  Double_t minMassStep[6]={1.68,1.70,1.72,1.74,1.76,1.78};
108  Double_t maxMassStep[6]={2.06,2.04,2.02,2.00,1.98,1.96};
109  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};
110  ConfigureRebinSteps(4,rebinStep);
111  ConfigureLowLimFitSteps(6,minMassStep);
112  ConfigureUpLimFitSteps(6,maxMassStep);
113  ConfigurenSigmaBinCSteps(11,nSigmasBC);
114 }
115 
116 //________________________________________________________________________
118  // destructor
119  delete [] fRebinSteps;
120  delete [] fLowLimFitSteps;
121  delete [] fUpLimFitSteps;
122  if(fhTemplRefl) delete fhTemplRefl;
123  if(fhTemplSign) delete fhTemplSign;
124  for (auto fitter : fMassFitters) delete fitter;
125 }
126 
127 //________________________________________________________________________
129  // creates output histograms
130 
131  const Int_t nCases=kNBkgFuncCases*kNFitConfCases;
132 
133  TString funcBkg[kNBkgFuncCases]={"Expo","Lin","Pol2","Pol3","Pol4","Pol5","PowLaw","PowLawExpo"};
134  TString gausSig[kNFitConfCases]={"FixedS","FixedSp20","FixedSm20","FreeS","FixedMeanFixedS","FixedMeanFreeS"};
135 
137 
138  fHistoRawYieldDistAll = new TH1F(Form("hRawYieldDistAll%s",fSuffix.Data())," ; Raw Yield",5000,0.,50000.);
139  fHistoRawYieldTrialAll = new TH1F(Form("hRawYieldTrialAll%s",fSuffix.Data())," ; Trial # ; Raw Yield",nCases*totTrials,-0.5,nCases*totTrials-0.5);
140  fHistoSigmaTrialAll = new TH1F(Form("hSigmaTrialAll%s",fSuffix.Data())," ; Trial # ; Sigma (GeV/c^{2})",nCases*totTrials,-0.5,nCases*totTrials-0.5);
141  fHistoMeanTrialAll = new TH1F(Form("hMeanTrialAll%s",fSuffix.Data())," ; Trial # ; Mean (GeV/c^{2})",nCases*totTrials,-0.5,nCases*totTrials-0.5);
142  fHistoChi2TrialAll = new TH1F(Form("hChi2TrialAll%s",fSuffix.Data())," ; Trial # ; #chi^{2}",nCases*totTrials,-0.5,nCases*totTrials-0.5);
143  fHistoSignifTrialAll = new TH1F(Form("hSignifTrialAll%s",fSuffix.Data())," ; Trial # ; Significance",nCases*totTrials,-0.5,nCases*totTrials-0.5);
144  if(fSaveBkgVal) {
145  fHistoBkgTrialAll = new TH1F(Form("hBkgTrialAll%s",fSuffix.Data())," ; Background",nCases*totTrials,-0.5,nCases*totTrials-0.5);
146  fHistoBkgInBinEdgesTrialAll = new TH1F(Form("hBkgInBinEdgesTrialAll%s",fSuffix.Data())," ; Background in bin edges",nCases*totTrials,-0.5,nCases*totTrials-0.5);
147  }
148 
149 
150  fHistoRawYieldDistBinC0All = new TH1F(Form("hRawYieldDistBinC0All%s",fSuffix.Data())," ; Raw Yield (bin count)",5000,0.,50000.);
151  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);
152  fHistoRawYieldDistBinC1All = new TH1F(Form("hRawYieldDistBinC1All%s",fSuffix.Data())," ; Raw Yield (bin count)",5000,0.,50000.);
153  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);
154 
155  fHistoRawYieldDist = new TH1F*[nCases];
156  fHistoRawYieldTrial = new TH1F*[nCases];
157  fHistoSigmaTrial = new TH1F*[nCases];
158  fHistoMeanTrial = new TH1F*[nCases];
159  fHistoChi2Trial = new TH1F*[nCases];
160  fHistoSignifTrial = new TH1F*[nCases];
161  if(fSaveBkgVal) {
162  fHistoBkgTrial = new TH1F*[nCases];
163  fHistoBkgInBinEdgesTrial = new TH1F*[nCases];
164  }
165 
166  fHistoRawYieldDistBinC0 = new TH1F*[nCases];
167  fHistoRawYieldTrialBinC0 = new TH2F*[nCases];
168  fHistoRawYieldDistBinC1 = new TH1F*[nCases];
169  fHistoRawYieldTrialBinC1 = new TH2F*[nCases];
170 
171  for(Int_t ib=0; ib<kNBkgFuncCases; ib++){
172  for(Int_t igs=0; igs<kNFitConfCases; igs++){
173  Int_t theCase=igs*kNBkgFuncCases+ib;
174  fHistoRawYieldDist[theCase]=new TH1F(Form("hRawYieldDist%s%s%s",funcBkg[ib].Data(),gausSig[igs].Data(),fSuffix.Data())," ; Raw Yield",5000,0.,50000.);
175  fHistoRawYieldDistBinC0[theCase]=new TH1F(Form("hRawYieldDistBinC0%s%s%s",funcBkg[ib].Data(),gausSig[igs].Data(),fSuffix.Data())," ; Raw Yield (bin count)",5000,0.,50000.);
176  fHistoRawYieldDistBinC1[theCase]=new TH1F(Form("hRawYieldDistBinC1%s%s%s",funcBkg[ib].Data(),gausSig[igs].Data(),fSuffix.Data())," ; Raw Yield (bin count)",5000,0.,50000.);
177  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);
178  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);
179  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);
180  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);
181  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);
182  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);
183  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);
184  if(fSaveBkgVal) {
185  fHistoBkgTrial[theCase] = new TH1F(Form("hBkgTrial%s%s%s",funcBkg[ib].Data(),gausSig[igs].Data(),fSuffix.Data())," ; Background",totTrials,-0.5,totTrials-0.5);
186  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);
187  }
188 
189  fHistoChi2Trial[theCase]->SetMarkerStyle(7);
190  fHistoSignifTrial[theCase]->SetMarkerStyle(7);
191  fHistoSigmaTrial[theCase]->SetMarkerStyle(7);
192  fHistoMeanTrial[theCase]->SetMarkerStyle(7);
193  if(fSaveBkgVal) {
194  fHistoBkgTrial[theCase]->SetMarkerStyle(7);
195  fHistoBkgInBinEdgesTrial[theCase]->SetMarkerStyle(7);
196  }
197 
198  }
199  }
200  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);
201  fNtupleMultiTrials->SetDirectory(nullptr);
202  return kTRUE;
203 
204 }
205 
206 //________________________________________________________________________
208  // perform the multiple fits
209 
210  Bool_t hOK=CreateHistos();
211  if(!hOK) return kFALSE;
212 
213  Int_t itrial=0;
214  Int_t types=0;
215  Int_t itrialBC=0;
217 
218  fMinYieldGlob=999999.;
219  fMaxYieldGlob=0.;
220  Float_t xnt[15];
221 
222  for(Int_t ir=0; ir<fNumOfRebinSteps; ir++){
223  Int_t rebin=fRebinSteps[ir];
224  for(Int_t iFirstBin=1; iFirstBin<=fNumOfFirstBinSteps; iFirstBin++) {
225  TH1F* hRebinned=0x0;
226  if(fNumOfFirstBinSteps==1) hRebinned=(TH1F*)AliVertexingHFUtils::RebinHisto(hInvMassHisto,rebin,-1);
227  else hRebinned=(TH1F*)AliVertexingHFUtils::RebinHisto(hInvMassHisto,rebin,iFirstBin);
228  for(Int_t iMinMass=0; iMinMass<fNumOfLowLimFitSteps; iMinMass++){
230  Double_t hmin=TMath::Max(minMassForFit,hRebinned->GetBinLowEdge(2));
231  for(Int_t iMaxMass=0; iMaxMass<fNumOfUpLimFitSteps; iMaxMass++){
233  Double_t hmax=TMath::Min(maxMassForFit,hRebinned->GetBinLowEdge(hRebinned->GetNbinsX()));
234  ++itrial;
235  for(Int_t typeb=0; typeb<kNBkgFuncCases; typeb++){
236  if(typeb==kExpoBkg && !fUseExpoBkg) continue;
237  if(typeb==kLinBkg && !fUseLinBkg) continue;
238  if(typeb==kPol2Bkg && !fUsePol2Bkg) continue;
239  if(typeb==kPol3Bkg && !fUsePol3Bkg) continue;
240  if(typeb==kPol4Bkg && !fUsePol4Bkg) continue;
241  if(typeb==kPol5Bkg && !fUsePol5Bkg) continue;
242  if(typeb==kPowBkg && !fUsePowLawBkg) continue;
244  for(Int_t igs=0; igs<kNFitConfCases; igs++){
245  if (igs==kFixSigUpFreeMean && !fUseFixSigUpFreeMean) continue;
246  if (igs==kFixSigDownFreeMean && !fUseFixSigDownFreeMean) continue;
247  if (igs==kFreeSigFixMean && !fUseFixedMeanFreeS) continue;
248  if (igs==kFreeSigFreeMean && !fUseFreeS) continue;
249  if (igs==kFixSigFreeMean && !fUseFixSigFreeMean) continue;
250  if (igs==kFixSigFixMean && !fUseFixSigFixMean) continue;
251  Int_t theCase=igs*kNBkgFuncCases+typeb;
252  Int_t globBin=itrial+theCase*totTrials;
253  for(Int_t j=0; j<15; j++) xnt[j]=0.;
254 
255  Bool_t mustDeleteFitter = kTRUE;
257  if(typeb==kExpoBkg){
258  fitter=new AliHFInvMassFitter(hRebinned, hmin, hmax, AliHFInvMassFitter::kExpo, types);
259  }else if(typeb==kLinBkg){
260  fitter=new AliHFInvMassFitter(hRebinned, hmin, hmax, AliHFInvMassFitter::kLin, types);
261  }else if(typeb==kPol2Bkg){
262  fitter=new AliHFInvMassFitter(hRebinned, hmin, hmax, AliHFInvMassFitter::kPol2, types);
263  }else if(typeb==kPowBkg){
264  fitter=new AliHFInvMassFitter(hRebinned, hmin, hmax, AliHFInvMassFitter::kPow, types);
265  }else if(typeb==kPowTimesExpoBkg){
266  fitter=new AliHFInvMassFitter(hRebinned, hmin, hmax, AliHFInvMassFitter::kPowEx, types);
267  }else{
268  fitter=new AliHFInvMassFitter(hRebinned, hmin, hmax, 6, types);
269  if(typeb==kPol3Bkg) fitter->SetPolDegreeForBackgroundFit(3);
270  if(typeb==kPol4Bkg) fitter->SetPolDegreeForBackgroundFit(4);
271  if(typeb==kPol5Bkg) fitter->SetPolDegreeForBackgroundFit(5);
272  }
273  // D0 Reflection
274  if(fhTemplRefl && fhTemplSign){
275  TH1F *hReflModif=(TH1F*)AliVertexingHFUtils::AdaptTemplateRangeAndBinning(fhTemplRefl,hRebinned,minMassForFit,maxMassForFit);
276  TH1F *hSigModif=(TH1F*)AliVertexingHFUtils::AdaptTemplateRangeAndBinning(fhTemplSign,hRebinned,minMassForFit,maxMassForFit);
277  TH1F* hrfl=fitter->SetTemplateReflections(hReflModif,"2gaus",minMassForFit,maxMassForFit);
278  if(fFixRefloS>0){
279  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)));
280  fitter->SetFixReflOverS(fixSoverRefAt);
281  }
282  delete hReflModif;
283  delete hSigModif;
284  }
285  if(fUseSecondPeak){
287  }
288  if(fFitOption==1) fitter->SetUseChi2Fit();
291  xnt[0]=rebin;
292  xnt[1]=iFirstBin;
293  xnt[2]=minMassForFit;
294  xnt[3]=maxMassForFit;
295  xnt[4]=typeb;
296  xnt[6]=0;
297  if(igs==kFixSigFreeMean){
299  xnt[5]=1;
300  }else if(igs==kFixSigUpFreeMean){
302  xnt[5]=2;
303  }else if(igs==kFixSigDownFreeMean){
305  xnt[5]=3;
306  }else if(igs==kFreeSigFreeMean){
307  xnt[5]=0;
308  }else if(igs==kFixSigFixMean){
310  fitter->SetFixGaussianMean(fMassD);
311  xnt[5]=1;
312  xnt[6]=1;
313  }else if(igs==kFreeSigFixMean){
314  fitter->SetFixGaussianMean(fMassD);
315  xnt[5]=0;
316  xnt[6]=1;
317  }
318  Bool_t out=kFALSE;
319  Double_t chisq=-1.;
320  Double_t sigma=0.;
321  Double_t esigma=0.;
322  Double_t pos=.0;
323  Double_t epos=.0;
324  Double_t ry=.0;
325  Double_t ery=.0;
326  Double_t significance=0.;
327  Double_t erSignif=0.;
328  Double_t bkg=0.;
329  Double_t erbkg=0.;
330  Double_t bkgBEdge=0;
331  Double_t erbkgBEdge=0;
332  TF1* fB1=0x0;
333  if(typeb<kNBkgFuncCases){
334  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);
335  out=fitter->MassFitter(0);
336  chisq=fitter->GetReducedChiSquare();
337  fitter->Significance(fnSigmaForBkgEval,significance,erSignif);
338  sigma=fitter->GetSigma();
339  pos=fitter->GetMean();
340  esigma=fitter->GetSigmaUncertainty();
341  if(esigma<0.00001) esigma=0.0001;
342  epos=fitter->GetMeanUncertainty();
343  if(epos<0.00001) epos=0.0001;
344  ry=fitter->GetRawYield();
345  ery=fitter->GetRawYieldError();
346  fB1=fitter->GetBackgroundFullRangeFunc();
347  fitter->Background(fnSigmaForBkgEval,bkg,erbkg);
348  Double_t minval = hInvMassHisto->GetXaxis()->GetBinLowEdge(hInvMassHisto->FindBin(pos-fnSigmaForBkgEval*sigma));
349  Double_t maxval = hInvMassHisto->GetXaxis()->GetBinUpEdge(hInvMassHisto->FindBin(pos+fnSigmaForBkgEval*sigma));
350  fitter->Background(minval,maxval,bkgBEdge,erbkgBEdge);
351  if(out && fDrawIndividualFits && thePad){
352  thePad->Clear();
353  fitter->DrawHere(thePad, fnSigmaForBkgEval);
354  fMassFitters.push_back(fitter);
355  mustDeleteFitter = kFALSE;
356  for (auto format : fInvMassFitSaveAsFormats) {
357  thePad->SaveAs(Form("FitOutput_%s_Trial%d.%s",hInvMassHisto->GetName(),globBin, format.c_str()));
358  }
359  }
360  }
361  // else{
362  // out=DoFitWithPol3Bkg(hRebinned,hmin,hmax,igs);
363  // if(out && thePad){
364  // thePad->Clear();
365  // hRebinned->Draw();
366  // TF1* fSB=(TF1*)hRebinned->GetListOfFunctions()->FindObject("fSB");
367  // fB1=new TF1("fB1","[0]+[1]*x+[2]*x*x+[3]*x*x*x",hmin,hmax);
368  // for(Int_t j=0; j<4; j++) fB1->SetParameter(j,fSB->GetParameter(3+j));
369  // fB1->SetLineColor(2);
370  // fB1->Draw("same");
371  // fSB->SetLineColor(4);
372  // fSB->Draw("same");
373  // thePad->Update();
374  // chisq=fSB->GetChisquare()/fSB->GetNDF();;
375  // sigma=fSB->GetParameter(2);
376  // esigma=fSB->GetParError(2);
377  // if(esigma<0.00001) esigma=0.0001;
378  // pos=fSB->GetParameter(1);
379  // epos=fSB->GetParError(1);
380  // if(epos<0.00001) epos=0.0001;
381  // ry=fSB->GetParameter(0)/hRebinned->GetBinWidth(1);
382  // ery=fSB->GetParError(0)/hRebinned->GetBinWidth(1);
383  // }
384  // }
385  xnt[7]=chisq;
386  if(out && chisq>0. && sigma>0.5*fSigmaGausMC && sigma<2.0*fSigmaGausMC){
387  xnt[8]=significance;
388  xnt[9]=pos;
389  xnt[10]=epos;
390  xnt[11]=sigma;
391  xnt[12]=esigma;
392  xnt[13]=ry;
393  xnt[14]=ery;
394  fHistoRawYieldDistAll->Fill(ry);
395  fHistoRawYieldTrialAll->SetBinContent(globBin,ry);
396  fHistoRawYieldTrialAll->SetBinError(globBin,ery);
397  fHistoSigmaTrialAll->SetBinContent(globBin,sigma);
398  fHistoSigmaTrialAll->SetBinError(globBin,esigma);
399  fHistoMeanTrialAll->SetBinContent(globBin,pos);
400  fHistoMeanTrialAll->SetBinError(globBin,epos);
401  fHistoChi2TrialAll->SetBinContent(globBin,chisq);
402  fHistoChi2TrialAll->SetBinError(globBin,0.00001);
403  fHistoSignifTrialAll->SetBinContent(globBin,significance);
404  fHistoSignifTrialAll->SetBinError(globBin,erSignif);
405  if(fSaveBkgVal) {
406  fHistoBkgTrialAll->SetBinContent(globBin,bkg);
407  fHistoBkgTrialAll->SetBinError(globBin,erbkg);
408  fHistoBkgInBinEdgesTrialAll->SetBinContent(globBin,bkgBEdge);
409  fHistoBkgInBinEdgesTrialAll->SetBinError(globBin,erbkgBEdge);
410  }
411 
412  if(ry<fMinYieldGlob) fMinYieldGlob=ry;
413  if(ry>fMaxYieldGlob) fMaxYieldGlob=ry;
414  fHistoRawYieldDist[theCase]->Fill(ry);
415  fHistoRawYieldTrial[theCase]->SetBinContent(itrial,ry);
416  fHistoRawYieldTrial[theCase]->SetBinError(itrial,ery);
417  fHistoSigmaTrial[theCase]->SetBinContent(itrial,sigma);
418  fHistoSigmaTrial[theCase]->SetBinError(itrial,esigma);
419  fHistoMeanTrial[theCase]->SetBinContent(itrial,pos);
420  fHistoMeanTrial[theCase]->SetBinError(itrial,epos);
421  fHistoChi2Trial[theCase]->SetBinContent(itrial,chisq);
422  fHistoChi2Trial[theCase]->SetBinError(itrial,0.00001);
423  fHistoSignifTrial[theCase]->SetBinContent(itrial,significance);
424  fHistoSignifTrial[theCase]->SetBinError(itrial,erSignif);
425  if(fSaveBkgVal) {
426  fHistoBkgTrial[theCase]->SetBinContent(itrial,bkg);
427  fHistoBkgTrial[theCase]->SetBinError(itrial,erbkg);
428  fHistoBkgInBinEdgesTrial[theCase]->SetBinContent(itrial,bkgBEdge);
429  fHistoBkgInBinEdgesTrial[theCase]->SetBinError(itrial,erbkgBEdge);
430  }
431 
432  for(Int_t iStepBC=0; iStepBC<fNumOfnSigmaBinCSteps; iStepBC++){
433  Double_t minMassBC=fMassD-fnSigmaBinCSteps[iStepBC]*sigma;
434  Double_t maxMassBC=fMassD+fnSigmaBinCSteps[iStepBC]*sigma;
435  if(minMassBC>minMassForFit &&
436  maxMassBC<maxMassForFit &&
437  minMassBC>(hRebinned->GetXaxis()->GetXmin()) &&
438  maxMassBC<(hRebinned->GetXaxis()->GetXmax())){
439  Double_t cnts0,ecnts0;
440  Double_t cnts1,ecnts1;
441  cnts0=fitter->GetRawYieldBinCounting(ecnts0,minMassBC,maxMassBC,0);
442  cnts1=fitter->GetRawYieldBinCounting(ecnts1,minMassBC,maxMassBC,1);
443  ++itrialBC;
444  fHistoRawYieldDistBinC0All->Fill(cnts0);
445  fHistoRawYieldTrialBinC0All->SetBinContent(globBin,iStepBC+1,cnts0);
446  fHistoRawYieldTrialBinC0All->SetBinError(globBin,iStepBC+1,ecnts0);
447  fHistoRawYieldTrialBinC0[theCase]->SetBinContent(itrial,iStepBC+1,cnts0);
448  fHistoRawYieldTrialBinC0[theCase]->SetBinError(itrial,iStepBC+1,ecnts0);
449  fHistoRawYieldDistBinC0[theCase]->Fill(cnts0);
450  fHistoRawYieldDistBinC1All->Fill(cnts1);
451  fHistoRawYieldTrialBinC1All->SetBinContent(globBin,iStepBC+1,cnts1);
452  fHistoRawYieldTrialBinC1All->SetBinError(globBin,iStepBC+1,ecnts1);
453  fHistoRawYieldTrialBinC1[theCase]->SetBinContent(itrial,iStepBC+1,cnts1);
454  fHistoRawYieldTrialBinC1[theCase]->SetBinError(itrial,iStepBC+1,ecnts1);
455  fHistoRawYieldDistBinC1[theCase]->Fill(cnts1);
456  }
457  }
458  }
459  if (mustDeleteFitter) delete fitter;
460  fNtupleMultiTrials->Fill(xnt);
461  }
462  }
463  }
464  }
465  delete hRebinned;
466  }
467  }
468  return kTRUE;
469 }
470 
471 //________________________________________________________________________
473  // save histos in a root file for further analysis
474  const Int_t nCases=kNBkgFuncCases*kNFitConfCases;
475  TFile outHistos(fileName.Data(),option.Data());
476  if (outHistos.IsZombie()) {
477  Printf("Could not open file '%s'!", fileName.Data());
478  return;
479  }
480  outHistos.cd();
481  fHistoRawYieldTrialAll->Write();
482  fHistoSigmaTrialAll->Write();
483  fHistoMeanTrialAll->Write();
484  fHistoChi2TrialAll->Write();
485  fHistoSignifTrialAll->Write();
486  if(fSaveBkgVal) {
487  fHistoBkgTrialAll->Write();
489  }
494  for(Int_t ic=0; ic<nCases; ic++){
495  fHistoRawYieldTrial[ic]->Write();
496  fHistoSigmaTrial[ic]->Write();
497  fHistoMeanTrial[ic]->Write();
498  fHistoChi2Trial[ic]->Write();
499  fHistoSignifTrial[ic]->Write();
500  if(fSaveBkgVal) {
501  fHistoBkgTrial[ic]->Write();
502  fHistoBkgInBinEdgesTrial[ic]->Write();
503  }
504  fHistoRawYieldTrialBinC0[ic]->Write();
505  fHistoRawYieldDistBinC0[ic]->Write();
506  fHistoRawYieldTrialBinC1[ic]->Write();
507  fHistoRawYieldDistBinC1[ic]->Write();
508  }
509  fNtupleMultiTrials->SetDirectory(&outHistos);
510  fNtupleMultiTrials->Write();
511  outHistos.Close();
512 }
513 
514 //________________________________________________________________________
515 void AliHFInvMassMultiTrialFit::DrawHistos(TCanvas* cry) const{
516  // draw histos
517  cry->Divide(2,2);
518  cry->cd(1);
519  fHistoSigmaTrialAll->Draw();
520  cry->cd(3);
521  fHistoChi2TrialAll->Draw();
522  cry->cd(2);
523  fHistoRawYieldTrialAll->Draw();
524  cry->cd(4);
525  fHistoRawYieldDistAll->Draw();
526  TLatex* tmean=new TLatex(0.15,0.8,Form("mean=%.1f",fHistoRawYieldDistAll->GetMean()));
527  tmean->SetNDC();
528  tmean->Draw();
529  TLatex* thrms=new TLatex(0.15,0.72,Form("rms=%.1f",fHistoRawYieldDistAll->GetRMS()));
530  thrms->SetNDC();
531  thrms->Draw();
532  TLatex* tmax=new TLatex(0.6,0.8,Form("max=%.1f",fMaxYieldGlob));
533  tmax->SetNDC();
534  tmax->Draw();
535  TLatex* tmin=new TLatex(0.6,0.72,Form("min=%.1f",fMinYieldGlob));
536  tmin->SetNDC();
537  tmin->Draw();
538  TLatex* trms=new TLatex(0.6,0.64,Form("(max-min)/sqrt(12)=%.1f",(fMaxYieldGlob-fMinYieldGlob)/sqrt(12)));
539  trms->SetNDC();
540  trms->Draw();
541 
542 }
543 //________________________________________________________________________
545  Int_t iCase){
546  //
547 
548  TH1F *hCutTmp=(TH1F*)histoToFit->Clone("hCutTmp");
549  for(Int_t ib=1; ib<=hCutTmp->GetNbinsX(); ib++){
550  Double_t xc=hCutTmp->GetBinCenter(ib);
551  if(xc>(fMassD-5.*fSigmaGausMC) && xc<(fMassD+5.*fSigmaGausMC)){
552  hCutTmp->SetBinContent(ib,0.);
553  hCutTmp->SetBinError(ib,0.);
554  }
555  }
556 
557  hCutTmp->Fit("pol2","E0","",hmin,hmax);
558  TF1* f2=(TF1*)hCutTmp->GetListOfFunctions()->FindObject("pol2");
559  TF1* f3=new TF1("myPol3","pol3");
560  for(Int_t i=0; i<3;i++) f3->SetParameter(i,f2->GetParameter(i));
561  hCutTmp->Fit(f3,"E0","",hmin,hmax);
562  Double_t quickCount=0.;
563  for(Int_t ib=1; ib<=histoToFit->GetNbinsX(); ib++){
564  Double_t xc=hCutTmp->GetBinCenter(ib);
565  if(xc>(fMassD-3.*fSigmaGausMC) && xc<(fMassD+3.*fSigmaGausMC)){
566  quickCount+=(histoToFit->GetBinContent(ib)-f3->Eval(xc));
567  }
568  }
569  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");
570  fSB->SetParameter(0,quickCount);
571  fSB->SetParameter(1,fMassD);
572  fSB->SetParameter(2,fSigmaGausMC);
573  for(Int_t j=0; j<4; j++) fSB->SetParameter(j+3,f3->GetParameter(j));
574  if(iCase==0) fSB->FixParameter(2,fSigmaGausMC);
575  else if(iCase==1) fSB->FixParameter(2,fSigmaGausMC*(1.+fSigmaMCVariation));
576  else if(iCase==2) fSB->FixParameter(2,fSigmaGausMC*(1.-fSigmaMCVariation));
577  else if(iCase==4){
578  fSB->FixParameter(1,fMassD);
579  fSB->FixParameter(2,fSigmaGausMC);
580  } else if(iCase==5){
581  fSB->FixParameter(1,fMassD);
582  }
583  histoToFit->Fit(fSB,"ME0","",hmin,hmax);
584  // quality cuts
585  if(fSB->GetParError(0)<0.01*fSB->GetParameter(0)) return kFALSE;
586  if(fSB->GetParError(0)>0.6*fSB->GetParameter(0)) return kFALSE;
587 
588  delete hCutTmp;
589  return kTRUE;
590 }
591 //__________________________________________________________________________________
592 void AliHFInvMassMultiTrialFit::SetTemplatesForReflections(const TH1F *hr, const TH1F *hs) {
594  if(fhTemplSign) delete fhTemplSign;
595  if(fhTemplRefl) delete fhTemplRefl;
596  fhTemplRefl=(TH1F*)hr->Clone("hTemplRefl");
597  fhTemplSign=(TH1F*)hs->Clone("hTemplSign");
598  return;
599 }
600 
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)
ClassImp(AliAnalysisTaskDeltaPt) AliAnalysisTaskDeltaPt
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
TF1 * GetBackgroundFullRangeFunc()
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