AliPhysics  71e3bc7 (71e3bc7)
 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 
31 
32 
33 //_________________________________________________________________________
35  TNamed(),
36  fNumOfRebinSteps(4),
37  fRebinSteps(0x0),
38  fNumOfFirstBinSteps(1),
39  fNumOfLowLimFitSteps(6),
40  fLowLimFitSteps(0x0),
41  fNumOfUpLimFitSteps(6),
42  fUpLimFitSteps(0x0),
43  fNumOfnSigmaBinCSteps(11),
44  fnSigmaBinCSteps(0x0),
45  fnSigmaForBkgEval(3),
46  fSigmaGausMC(0.010),
47  fSigmaMCVariation(0.15),
48  fMassD(1.86484),
49  fSuffix(""),
50  fFitOption(0),
51  fUseExpoBkg(kTRUE),
52  fUseLinBkg(kTRUE),
53  fUsePol2Bkg(kTRUE),
54  fUsePol3Bkg(kTRUE),
55  fUsePol4Bkg(kTRUE),
56  fUsePol5Bkg(kFALSE),
57  fUsePowLawBkg(kFALSE),
58  fUsePowLawTimesExpoBkg(kFALSE),
59  fUseFixSigUpFreeMean(kTRUE),
60  fUseFixSigDownFreeMean(kTRUE),
61  fUseFreeS(kTRUE),
62  fUseFixedMeanFreeS(kTRUE),
63  fUseFixSigFreeMean(kTRUE),
64  fUseFixSigFixMean(kTRUE),
65  fUseSecondPeak(kFALSE),
66  fMassSecondPeak(1.86958),
67  fSigmaSecondPeak(0.01),
68  fFixMassSecondPeak(kFALSE),
69  fFixSigmaSecondPeak(kFALSE),
70  fSaveBkgVal(kFALSE),
71  fDrawIndividualFits(kFALSE),
72  fHistoRawYieldDistAll(0x0),
73  fHistoRawYieldTrialAll(0x0),
74  fHistoSigmaTrialAll(0x0),
75  fHistoMeanTrialAll(0x0),
76  fHistoChi2TrialAll(0x0),
77  fHistoSignifTrialAll(0x0),
78  fHistoBkgTrialAll(0x0),
79  fHistoBkgInBinEdgesTrialAll(0x0),
80  fHistoRawYieldDistBinCAll(0x0),
81  fHistoRawYieldTrialBinCAll(0x0),
82  fHistoRawYieldDist(0x0),
83  fHistoRawYieldTrial(0x0),
84  fHistoSigmaTrial(0x0),
85  fHistoMeanTrial(0x0),
86  fHistoChi2Trial(0x0),
87  fHistoSignifTrial(0x0),
88  fHistoBkgTrial(0x0),
89  fHistoBkgInBinEdgesTrial(0x0),
90  fHistoRawYieldDistBinC(0x0),
91  fHistoRawYieldTrialBinC(0x0),
92  fhTemplRefl(0x0),
93  fFixRefloS(0),
94  fNtupleMultiTrials(0x0),
95  fMinYieldGlob(0),
96  fMaxYieldGlob(0),
97  fMassFitters()
98 {
99  // constructor
100  Int_t rebinStep[4]={3,4,5,6};
101  Double_t minMassStep[6]={1.68,1.70,1.72,1.74,1.76,1.78};
102  Double_t maxMassStep[6]={2.06,2.04,2.02,2.00,1.98,1.96};
103  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};
104  ConfigureRebinSteps(4,rebinStep);
105  ConfigureLowLimFitSteps(6,minMassStep);
106  ConfigureUpLimFitSteps(6,maxMassStep);
107  ConfigurenSigmaBinCSteps(11,nSigmasBC);
108 }
109 
110 //________________________________________________________________________
112  // destructor
113  delete [] fRebinSteps;
114  delete [] fLowLimFitSteps;
115  delete [] fUpLimFitSteps;
116  if(fhTemplRefl) delete fhTemplRefl;
117  for (auto fitter : fMassFitters) delete fitter;
118 }
119 
120 //________________________________________________________________________
122  // creates output histograms
123 
124  const Int_t nCases=kNBkgFuncCases*kNFitConfCases;
125 
126  TString funcBkg[kNBkgFuncCases]={"Expo","Lin","Pol2","Pol3","Pol4","Pol5","PowLaw","PowLawExpo"};
127  TString gausSig[kNFitConfCases]={"FixedS","FixedSp20","FixedSm20","FreeS","FixedMeanFixedS","FixedMeanFreeS"};
128 
130 
131  fHistoRawYieldDistAll = new TH1F(Form("hRawYieldDistAll%s",fSuffix.Data())," ; Raw Yield",5000,0.,50000.);
132  fHistoRawYieldTrialAll = new TH1F(Form("hRawYieldTrialAll%s",fSuffix.Data())," ; Trial # ; Raw Yield",nCases*totTrials,-0.5,nCases*totTrials-0.5);
133  fHistoSigmaTrialAll = new TH1F(Form("hSigmaTrialAll%s",fSuffix.Data())," ; Trial # ; Sigma (GeV/c^{2})",nCases*totTrials,-0.5,nCases*totTrials-0.5);
134  fHistoMeanTrialAll = new TH1F(Form("hMeanTrialAll%s",fSuffix.Data())," ; Trial # ; Mean (GeV/c^{2})",nCases*totTrials,-0.5,nCases*totTrials-0.5);
135  fHistoChi2TrialAll = new TH1F(Form("hChi2TrialAll%s",fSuffix.Data())," ; Trial # ; #chi^{2}",nCases*totTrials,-0.5,nCases*totTrials-0.5);
136  fHistoSignifTrialAll = new TH1F(Form("hSignifTrialAll%s",fSuffix.Data())," ; Trial # ; Significance",nCases*totTrials,-0.5,nCases*totTrials-0.5);
137  if(fSaveBkgVal) {
138  fHistoBkgTrialAll = new TH1F(Form("hBkgTrialAll%s",fSuffix.Data())," ; Background",nCases*totTrials,-0.5,nCases*totTrials-0.5);
139  fHistoBkgInBinEdgesTrialAll = new TH1F(Form("hBkgInBinEdgesTrialAll%s",fSuffix.Data())," ; Background in bin edges",nCases*totTrials,-0.5,nCases*totTrials-0.5);
140  }
141 
142 
143  fHistoRawYieldDistBinCAll = new TH1F(Form("hRawYieldDistBinCAll%s",fSuffix.Data())," ; Raw Yield (bin count)",5000,0.,50000.);
144  fHistoRawYieldTrialBinCAll = new TH2F(Form("hRawYieldTrialBinCAll%s",fSuffix.Data())," ; Trial # ; Range for count ; Raw Yield (bin count)",totTrials,-0.5,totTrials-0.5,fNumOfnSigmaBinCSteps,-0.5,fNumOfnSigmaBinCSteps-0.5);
145 
146  fHistoRawYieldDist = new TH1F*[nCases];
147  fHistoRawYieldTrial = new TH1F*[nCases];
148  fHistoSigmaTrial = new TH1F*[nCases];
149  fHistoMeanTrial = new TH1F*[nCases];
150  fHistoChi2Trial = new TH1F*[nCases];
151  fHistoSignifTrial = new TH1F*[nCases];
152  if(fSaveBkgVal) {
153  fHistoBkgTrial = new TH1F*[nCases];
154  fHistoBkgInBinEdgesTrial = new TH1F*[nCases];
155  }
156 
157  fHistoRawYieldDistBinC = new TH1F*[nCases];
158  fHistoRawYieldTrialBinC = new TH2F*[nCases];
159 
160  for(Int_t ib=0; ib<kNBkgFuncCases; ib++){
161  for(Int_t igs=0; igs<kNFitConfCases; igs++){
162  Int_t theCase=igs*kNBkgFuncCases+ib;
163  fHistoRawYieldDist[theCase]=new TH1F(Form("hRawYieldDist%s%s%s",funcBkg[ib].Data(),gausSig[igs].Data(),fSuffix.Data())," ; Raw Yield",5000,0.,50000.);
164  fHistoRawYieldDistBinC[theCase]=new TH1F(Form("hRawYieldDistBinC%s%s%s",funcBkg[ib].Data(),gausSig[igs].Data(),fSuffix.Data())," ; Raw Yield (bin count)",5000,0.,50000.);
165  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);
166  fHistoRawYieldTrialBinC[theCase]=new TH2F(Form("hRawYieldTrialBinC%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);
167  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);
168  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);
169  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);
170  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);
171  if(fSaveBkgVal) {
172  fHistoBkgTrial[theCase] = new TH1F(Form("hBkgTrial%s%s%s",funcBkg[ib].Data(),gausSig[igs].Data(),fSuffix.Data())," ; Background",totTrials,-0.5,totTrials-0.5);
173  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);
174  }
175 
176  fHistoChi2Trial[theCase]->SetMarkerStyle(7);
177  fHistoSignifTrial[theCase]->SetMarkerStyle(7);
178  fHistoSigmaTrial[theCase]->SetMarkerStyle(7);
179  fHistoMeanTrial[theCase]->SetMarkerStyle(7);
180  if(fSaveBkgVal) {
181  fHistoBkgTrial[theCase]->SetMarkerStyle(7);
182  fHistoBkgInBinEdgesTrial[theCase]->SetMarkerStyle(7);
183  }
184 
185  }
186  }
187  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);
188  fNtupleMultiTrials->SetDirectory(nullptr);
189  return kTRUE;
190 
191 }
192 
193 //________________________________________________________________________
195  // perform the multiple fits
196 
197  Bool_t hOK=CreateHistos();
198  if(!hOK) return kFALSE;
199 
200  Int_t itrial=0;
201  Int_t types=0;
202  Int_t itrialBC=0;
204 
205  fMinYieldGlob=999999.;
206  fMaxYieldGlob=0.;
207  Float_t xnt[15];
208 
209  for(Int_t ir=0; ir<fNumOfRebinSteps; ir++){
210  Int_t rebin=fRebinSteps[ir];
211  for(Int_t iFirstBin=1; iFirstBin<=fNumOfFirstBinSteps; iFirstBin++) {
212  TH1F* hRebinned=0x0;
213  if(fNumOfFirstBinSteps==1) hRebinned=RebinHisto(hInvMassHisto,rebin,-1);
214  else hRebinned=RebinHisto(hInvMassHisto,rebin,iFirstBin);
215  for(Int_t iMinMass=0; iMinMass<fNumOfLowLimFitSteps; iMinMass++){
217  Double_t hmin=TMath::Max(minMassForFit,hRebinned->GetBinLowEdge(2));
218  for(Int_t iMaxMass=0; iMaxMass<fNumOfUpLimFitSteps; iMaxMass++){
220  Double_t hmax=TMath::Min(maxMassForFit,hRebinned->GetBinLowEdge(hRebinned->GetNbinsX()));
221  ++itrial;
222  for(Int_t typeb=0; typeb<kNBkgFuncCases; typeb++){
223  if(typeb==kExpoBkg && !fUseExpoBkg) continue;
224  if(typeb==kLinBkg && !fUseLinBkg) continue;
225  if(typeb==kPol2Bkg && !fUsePol2Bkg) continue;
226  if(typeb==kPol3Bkg && !fUsePol3Bkg) continue;
227  if(typeb==kPol4Bkg && !fUsePol4Bkg) continue;
228  if(typeb==kPol5Bkg && !fUsePol5Bkg) continue;
229  if(typeb==kPowBkg && !fUsePowLawBkg) continue;
231  for(Int_t igs=0; igs<kNFitConfCases; igs++){
232  if (igs==kFixSigUpFreeMean && !fUseFixSigUpFreeMean) continue;
233  if (igs==kFixSigDownFreeMean && !fUseFixSigDownFreeMean) continue;
234  if (igs==kFreeSigFixMean && !fUseFixedMeanFreeS) continue;
235  if (igs==kFreeSigFreeMean && !fUseFreeS) continue;
236  if (igs==kFixSigFreeMean && !fUseFixSigFreeMean) continue;
237  if (igs==kFixSigFixMean && !fUseFixSigFixMean) continue;
238  Int_t theCase=igs*kNBkgFuncCases+typeb;
239  Int_t globBin=itrial+theCase*totTrials;
240  for(Int_t j=0; j<15; j++) xnt[j]=0.;
241 
242  Bool_t mustDeleteFitter = kTRUE;
244  if(typeb==kExpoBkg){
245  fitter=new AliHFInvMassFitter(hRebinned, hmin, hmax, AliHFInvMassFitter::kExpo, types);
246  }else if(typeb==kLinBkg){
247  fitter=new AliHFInvMassFitter(hRebinned, hmin, hmax, AliHFInvMassFitter::kLin, types);
248  }else if(typeb==kPol2Bkg){
249  fitter=new AliHFInvMassFitter(hRebinned, hmin, hmax, AliHFInvMassFitter::kPol2, types);
250  }else if(typeb==kPowBkg){
251  fitter=new AliHFInvMassFitter(hRebinned, hmin, hmax, AliHFInvMassFitter::kPow, types);
252  }else if(typeb==kPowTimesExpoBkg){
253  fitter=new AliHFInvMassFitter(hRebinned, hmin, hmax, AliHFInvMassFitter::kPowEx, types);
254  }else{
255  fitter=new AliHFInvMassFitter(hRebinned, hmin, hmax, 6, types);
256  if(typeb==kPol3Bkg) fitter->SetPolDegreeForBackgroundFit(3);
257  if(typeb==kPol4Bkg) fitter->SetPolDegreeForBackgroundFit(4);
258  if(typeb==kPol5Bkg) fitter->SetPolDegreeForBackgroundFit(5);
259  }
260  // D0 Reflection
261  if(fhTemplRefl){
262  TH1F* hrfl=fitter->SetTemplateReflections(fhTemplRefl,"2gaus",minMassForFit,maxMassForFit);
263  if(hrfl){
264  fitter->SetFixReflOverS(fFixRefloS);
265  }
266  }
267  if(fUseSecondPeak){
269  }
270  if(fFitOption==1) fitter->SetUseChi2Fit();
273  xnt[0]=rebin;
274  xnt[1]=iFirstBin;
275  xnt[2]=minMassForFit;
276  xnt[3]=maxMassForFit;
277  xnt[4]=typeb;
278  xnt[6]=0;
279  if(igs==kFixSigFreeMean){
281  xnt[5]=1;
282  }else if(igs==kFixSigUpFreeMean){
284  xnt[5]=2;
285  }else if(igs==kFixSigDownFreeMean){
287  xnt[5]=3;
288  }else if(igs==kFreeSigFreeMean){
289  xnt[5]=0;
290  }else if(igs==kFixSigFixMean){
292  fitter->SetFixGaussianMean(fMassD);
293  xnt[5]=1;
294  xnt[6]=1;
295  }else if(igs==kFreeSigFixMean){
296  fitter->SetFixGaussianMean(fMassD);
297  xnt[5]=0;
298  xnt[6]=1;
299  }
300  Bool_t out=kFALSE;
301  Double_t chisq=-1.;
302  Double_t sigma=0.;
303  Double_t esigma=0.;
304  Double_t pos=.0;
305  Double_t epos=.0;
306  Double_t ry=.0;
307  Double_t ery=.0;
308  Double_t significance=0.;
309  Double_t erSignif=0.;
310  Double_t bkg=0.;
311  Double_t erbkg=0.;
312  Double_t bkgBEdge=0;
313  Double_t erbkgBEdge=0;
314  TF1* fB1=0x0;
315  if(typeb<kNBkgFuncCases){
316  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);
317  out=fitter->MassFitter(0);
318  chisq=fitter->GetReducedChiSquare();
319  fitter->Significance(fnSigmaForBkgEval,significance,erSignif);
320  sigma=fitter->GetSigma();
321  pos=fitter->GetMean();
322  esigma=fitter->GetSigmaUncertainty();
323  if(esigma<0.00001) esigma=0.0001;
324  epos=fitter->GetMeanUncertainty();
325  if(epos<0.00001) epos=0.0001;
326  ry=fitter->GetRawYield();
327  ery=fitter->GetRawYieldError();
328  fB1=fitter->GetBackgroundFullRangeFunc();
329  fitter->Background(fnSigmaForBkgEval,bkg,erbkg);
330  Double_t minval = hInvMassHisto->GetXaxis()->GetBinLowEdge(hInvMassHisto->FindBin(pos-fnSigmaForBkgEval*sigma));
331  Double_t maxval = hInvMassHisto->GetXaxis()->GetBinUpEdge(hInvMassHisto->FindBin(pos+fnSigmaForBkgEval*sigma));
332  fitter->Background(minval,maxval,bkgBEdge,erbkgBEdge);
333  if(out && fDrawIndividualFits && thePad){
334  thePad->Clear();
335  fitter->DrawHere(thePad, fnSigmaForBkgEval);
336  fMassFitters.push_back(fitter);
337  mustDeleteFitter = kFALSE;
338  for (auto format : fInvMassFitSaveAsFormats) {
339  thePad->SaveAs(Form("FitOutput_%s_Trial%d.%s",hInvMassHisto->GetName(),globBin, format.c_str()));
340  }
341  }
342  }
343  // else{
344  // out=DoFitWithPol3Bkg(hRebinned,hmin,hmax,igs);
345  // if(out && thePad){
346  // thePad->Clear();
347  // hRebinned->Draw();
348  // TF1* fSB=(TF1*)hRebinned->GetListOfFunctions()->FindObject("fSB");
349  // fB1=new TF1("fB1","[0]+[1]*x+[2]*x*x+[3]*x*x*x",hmin,hmax);
350  // for(Int_t j=0; j<4; j++) fB1->SetParameter(j,fSB->GetParameter(3+j));
351  // fB1->SetLineColor(2);
352  // fB1->Draw("same");
353  // fSB->SetLineColor(4);
354  // fSB->Draw("same");
355  // thePad->Update();
356  // chisq=fSB->GetChisquare()/fSB->GetNDF();;
357  // sigma=fSB->GetParameter(2);
358  // esigma=fSB->GetParError(2);
359  // if(esigma<0.00001) esigma=0.0001;
360  // pos=fSB->GetParameter(1);
361  // epos=fSB->GetParError(1);
362  // if(epos<0.00001) epos=0.0001;
363  // ry=fSB->GetParameter(0)/hRebinned->GetBinWidth(1);
364  // ery=fSB->GetParError(0)/hRebinned->GetBinWidth(1);
365  // }
366  // }
367  xnt[7]=chisq;
368  if(out && chisq>0. && sigma>0.5*fSigmaGausMC && sigma<2.0*fSigmaGausMC){
369  xnt[8]=significance;
370  xnt[9]=pos;
371  xnt[10]=epos;
372  xnt[11]=sigma;
373  xnt[12]=esigma;
374  xnt[13]=ry;
375  xnt[14]=ery;
376  fHistoRawYieldDistAll->Fill(ry);
377  fHistoRawYieldTrialAll->SetBinContent(globBin,ry);
378  fHistoRawYieldTrialAll->SetBinError(globBin,ery);
379  fHistoSigmaTrialAll->SetBinContent(globBin,sigma);
380  fHistoSigmaTrialAll->SetBinError(globBin,esigma);
381  fHistoMeanTrialAll->SetBinContent(globBin,pos);
382  fHistoMeanTrialAll->SetBinError(globBin,epos);
383  fHistoChi2TrialAll->SetBinContent(globBin,chisq);
384  fHistoChi2TrialAll->SetBinError(globBin,0.00001);
385  fHistoSignifTrialAll->SetBinContent(globBin,significance);
386  fHistoSignifTrialAll->SetBinError(globBin,erSignif);
387  if(fSaveBkgVal) {
388  fHistoBkgTrialAll->SetBinContent(globBin,bkg);
389  fHistoBkgTrialAll->SetBinError(globBin,erbkg);
390  fHistoBkgInBinEdgesTrialAll->SetBinContent(globBin,bkgBEdge);
391  fHistoBkgInBinEdgesTrialAll->SetBinError(globBin,erbkgBEdge);
392  }
393 
394  if(ry<fMinYieldGlob) fMinYieldGlob=ry;
395  if(ry>fMaxYieldGlob) fMaxYieldGlob=ry;
396  fHistoRawYieldDist[theCase]->Fill(ry);
397  fHistoRawYieldTrial[theCase]->SetBinContent(itrial,ry);
398  fHistoRawYieldTrial[theCase]->SetBinError(itrial,ery);
399  fHistoSigmaTrial[theCase]->SetBinContent(itrial,sigma);
400  fHistoSigmaTrial[theCase]->SetBinError(itrial,esigma);
401  fHistoMeanTrial[theCase]->SetBinContent(itrial,pos);
402  fHistoMeanTrial[theCase]->SetBinError(itrial,epos);
403  fHistoChi2Trial[theCase]->SetBinContent(itrial,chisq);
404  fHistoChi2Trial[theCase]->SetBinError(itrial,0.00001);
405  fHistoSignifTrial[theCase]->SetBinContent(itrial,significance);
406  fHistoSignifTrial[theCase]->SetBinError(itrial,erSignif);
407  if(fSaveBkgVal) {
408  fHistoBkgTrial[theCase]->SetBinContent(itrial,bkg);
409  fHistoBkgTrial[theCase]->SetBinError(itrial,erbkg);
410  fHistoBkgInBinEdgesTrial[theCase]->SetBinContent(itrial,bkgBEdge);
411  fHistoBkgInBinEdgesTrial[theCase]->SetBinError(itrial,erbkgBEdge);
412  }
413 
414  for(Int_t iStepBC=0; iStepBC<fNumOfnSigmaBinCSteps; iStepBC++){
415  Double_t minMassBC=fMassD-fnSigmaBinCSteps[iStepBC]*sigma;
416  Double_t maxMassBC=fMassD+fnSigmaBinCSteps[iStepBC]*sigma;
417  if(minMassBC>minMassForFit &&
418  maxMassBC<maxMassForFit &&
419  minMassBC>(hRebinned->GetXaxis()->GetXmin()) &&
420  maxMassBC<(hRebinned->GetXaxis()->GetXmax())){
421  Double_t cnts,ecnts;
422  BinCount(hRebinned,fB1,1,minMassBC,maxMassBC,cnts,ecnts);
423  ++itrialBC;
424  fHistoRawYieldDistBinCAll->Fill(cnts);
425  fHistoRawYieldTrialBinCAll->SetBinContent(globBin,iStepBC+1,cnts);
426  fHistoRawYieldTrialBinCAll->SetBinError(globBin,iStepBC+1,ecnts);
427  fHistoRawYieldTrialBinC[theCase]->SetBinContent(itrial,iStepBC+1,cnts);
428  fHistoRawYieldTrialBinC[theCase]->SetBinError(itrial,iStepBC+1,ecnts);
429  fHistoRawYieldDistBinC[theCase]->Fill(cnts);
430  }
431  }
432  }
433  if (mustDeleteFitter) delete fitter;
434  fNtupleMultiTrials->Fill(xnt);
435  }
436  }
437  }
438  }
439  delete hRebinned;
440  }
441  }
442  return kTRUE;
443 }
444 
445 //________________________________________________________________________
447  // save histos in a root file for further analysis
448  const Int_t nCases=kNBkgFuncCases*kNFitConfCases;
449  TFile outHistos(fileName.Data(),option.Data());
450  if (outHistos.IsZombie()) {
451  Printf("Could not open file '%s'!", fileName.Data());
452  return;
453  }
454  outHistos.cd();
455  fHistoRawYieldTrialAll->Write();
456  fHistoSigmaTrialAll->Write();
457  fHistoMeanTrialAll->Write();
458  fHistoChi2TrialAll->Write();
459  fHistoSignifTrialAll->Write();
460  if(fSaveBkgVal) {
461  fHistoBkgTrialAll->Write();
462  fHistoBkgInBinEdgesTrialAll->Write();
463  }
464  fHistoRawYieldDistBinCAll->Write();
465  fHistoRawYieldTrialBinCAll->Write();
466  for(Int_t ic=0; ic<nCases; ic++){
467  fHistoRawYieldTrial[ic]->Write();
468  fHistoSigmaTrial[ic]->Write();
469  fHistoMeanTrial[ic]->Write();
470  fHistoChi2Trial[ic]->Write();
471  fHistoSignifTrial[ic]->Write();
472  if(fSaveBkgVal) {
473  fHistoBkgTrial[ic]->Write();
474  fHistoBkgInBinEdgesTrial[ic]->Write();
475  }
476  fHistoRawYieldTrialBinC[ic]->Write();
477  fHistoRawYieldDistBinC[ic]->Write();
478  }
479  fNtupleMultiTrials->SetDirectory(&outHistos);
480  fNtupleMultiTrials->Write();
481  outHistos.Close();
482 }
483 
484 //________________________________________________________________________
485 void AliHFInvMassMultiTrialFit::DrawHistos(TCanvas* cry) const{
486  // draw histos
487  cry->Divide(2,2);
488  cry->cd(1);
489  fHistoSigmaTrialAll->Draw();
490  cry->cd(3);
491  fHistoChi2TrialAll->Draw();
492  cry->cd(2);
493  fHistoRawYieldTrialAll->Draw();
494  cry->cd(4);
495  fHistoRawYieldDistAll->Draw();
496  TLatex* tmean=new TLatex(0.15,0.8,Form("mean=%.1f",fHistoRawYieldDistAll->GetMean()));
497  tmean->SetNDC();
498  tmean->Draw();
499  TLatex* thrms=new TLatex(0.15,0.72,Form("rms=%.1f",fHistoRawYieldDistAll->GetRMS()));
500  thrms->SetNDC();
501  thrms->Draw();
502  TLatex* tmax=new TLatex(0.6,0.8,Form("max=%.1f",fMaxYieldGlob));
503  tmax->SetNDC();
504  tmax->Draw();
505  TLatex* tmin=new TLatex(0.6,0.72,Form("min=%.1f",fMinYieldGlob));
506  tmin->SetNDC();
507  tmin->Draw();
508  TLatex* trms=new TLatex(0.6,0.64,Form("(max-min)/sqrt(12)=%.1f",(fMaxYieldGlob-fMinYieldGlob)/sqrt(12)));
509  trms->SetNDC();
510  trms->Draw();
511 
512 }
513 //________________________________________________________________________
514 TH1F* AliHFInvMassMultiTrialFit::RebinHisto(TH1D* hOrig, Int_t reb, Int_t firstUse) const{
515  // Rebin histogram, from bin firstUse to lastUse
516  // Use all bins if firstUse=-1
517 
518  Int_t nBinOrig=hOrig->GetNbinsX();
519  Int_t firstBinOrig=1;
520  Int_t lastBinOrig=nBinOrig;
521  Int_t nBinOrigUsed=nBinOrig;
522  Int_t nBinFinal=nBinOrig/reb;
523  if(firstUse>=1){
524  firstBinOrig=firstUse;
525  nBinFinal=(nBinOrig-firstUse+1)/reb;
526  nBinOrigUsed=nBinFinal*reb;
527  lastBinOrig=firstBinOrig+nBinOrigUsed-1;
528  }else{
529  Int_t exc=nBinOrigUsed%reb;
530  if(exc!=0){
531  nBinOrigUsed-=exc;
532  firstBinOrig+=exc/2;
533  lastBinOrig=firstBinOrig+nBinOrigUsed-1;
534  }
535  }
536 
537  printf("Rebin from %d bins to %d bins -- Used bins=%d in range %d-%d\n",nBinOrig,nBinFinal,nBinOrigUsed,firstBinOrig,lastBinOrig);
538  Float_t lowLim=hOrig->GetXaxis()->GetBinLowEdge(firstBinOrig);
539  Float_t hiLim=hOrig->GetXaxis()->GetBinUpEdge(lastBinOrig);
540  TH1F* hRebin=new TH1F(Form("%s-rebin%d_%d",hOrig->GetName(),reb,firstUse),hOrig->GetTitle(),nBinFinal,lowLim,hiLim);
541  Int_t lastSummed=firstBinOrig-1;
542  for(Int_t iBin=1;iBin<=nBinFinal; iBin++){
543  Float_t sum=0.;
544  Float_t sum2=0.;
545  for(Int_t iOrigBin=0;iOrigBin<reb;iOrigBin++){
546  sum+=hOrig->GetBinContent(lastSummed+1);
547  sum2+=hOrig->GetBinError(lastSummed+1)*hOrig->GetBinError(lastSummed+1);
548  lastSummed++;
549  }
550  hRebin->SetBinContent(iBin,sum);
551  hRebin->SetBinError(iBin,TMath::Sqrt(sum2));
552  }
553  return hRebin;
554 }
555 
556 //________________________________________________________________________
558  // compute yield with bin couting
559  Int_t minBinSum=h->FindBin(minMass);
560  Int_t maxBinSum=h->FindBin(maxMass);
561  Double_t cntSig=0.;
562  Double_t cntErr=0.;
563  for(Int_t iMB=minBinSum; iMB<=maxBinSum; iMB++){
564  Double_t bkg=fB ? fB->Eval(h->GetBinCenter(iMB))/(Double_t)rebin : 0;
565  cntSig+=(h->GetBinContent(iMB)-bkg);
566  cntErr+=(h->GetBinError(iMB)*h->GetBinError(iMB));
567  }
568  count=cntSig;
569  ecount=TMath::Sqrt(cntErr);
570 }
571 //________________________________________________________________________
573  Int_t iCase){
574  //
575 
576  TH1F *hCutTmp=(TH1F*)histoToFit->Clone("hCutTmp");
577  for(Int_t ib=1; ib<=hCutTmp->GetNbinsX(); ib++){
578  Double_t xc=hCutTmp->GetBinCenter(ib);
579  if(xc>(fMassD-5.*fSigmaGausMC) && xc<(fMassD+5.*fSigmaGausMC)){
580  hCutTmp->SetBinContent(ib,0.);
581  hCutTmp->SetBinError(ib,0.);
582  }
583  }
584 
585  hCutTmp->Fit("pol2","E0","",hmin,hmax);
586  TF1* f2=(TF1*)hCutTmp->GetListOfFunctions()->FindObject("pol2");
587  TF1* f3=new TF1("myPol3","pol3");
588  for(Int_t i=0; i<3;i++) f3->SetParameter(i,f2->GetParameter(i));
589  hCutTmp->Fit(f3,"E0","",hmin,hmax);
590  Double_t quickCount=0.;
591  for(Int_t ib=1; ib<=histoToFit->GetNbinsX(); ib++){
592  Double_t xc=hCutTmp->GetBinCenter(ib);
593  if(xc>(fMassD-3.*fSigmaGausMC) && xc<(fMassD+3.*fSigmaGausMC)){
594  quickCount+=(histoToFit->GetBinContent(ib)-f3->Eval(xc));
595  }
596  }
597  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");
598  fSB->SetParameter(0,quickCount);
599  fSB->SetParameter(1,fMassD);
600  fSB->SetParameter(2,fSigmaGausMC);
601  for(Int_t j=0; j<4; j++) fSB->SetParameter(j+3,f3->GetParameter(j));
602  if(iCase==0) fSB->FixParameter(2,fSigmaGausMC);
603  else if(iCase==1) fSB->FixParameter(2,fSigmaGausMC*(1.+fSigmaMCVariation));
604  else if(iCase==2) fSB->FixParameter(2,fSigmaGausMC*(1.-fSigmaMCVariation));
605  else if(iCase==4){
606  fSB->FixParameter(1,fMassD);
607  fSB->FixParameter(2,fSigmaGausMC);
608  } else if(iCase==5){
609  fSB->FixParameter(1,fMassD);
610  }
611  histoToFit->Fit(fSB,"ME0","",hmin,hmax);
612  // quality cuts
613  if(fSB->GetParError(0)<0.01*fSB->GetParameter(0)) return kFALSE;
614  if(fSB->GetParError(0)>0.6*fSB->GetParameter(0)) return kFALSE;
615 
616  delete hCutTmp;
617  return kTRUE;
618 }
619 //__________________________________________________________________________________
621  fhTemplRefl=(TH1F*)h->Clone("hTemplRefl");
622  return fhTemplRefl;
623 }
624 
Bool_t fUsePol4Bkg
switch for pol3 background
Double_t GetMeanUncertainty() const
Int_t * fRebinSteps
number of rebin steps
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
TH1F * RebinHisto(TH1D *hOrig, Int_t reb, Int_t firstUse) const
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
TH2F * fHistoRawYieldTrialBinCAll
histo with bin counts from all trials
Double_t fSigmaSecondPeak
position of the 2nd peak
TString fileName
Float_t fFixRefloS
template of reflection contribution
TString format
file names tag, basically the trigger and calorimeter combination
TH2F ** fHistoRawYieldTrialBinC
histo with bin counts from subsamples of trials
Double_t fSigmaMCVariation
sigma of D meson peak from MC
TH1F ** fHistoRawYieldDist
histo with bin counts from all 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
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
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 * 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
TH1F ** fHistoRawYieldDistBinC
histo with bkg in mass bin edges from subsamples of trials
Int_t fNumOfLowLimFitSteps
number of steps in the first bin for rebin
float Float_t
Definition: External.C:68
Bool_t DoFitWithPol3Bkg(TH1F *histoToFit, Double_t hmin, Double_t hmax, Int_t theCase)
TH1F * fHistoRawYieldTrialAll
histo with yield from all trials
TH1F * SetTemplateRefl(const TH1F *hTemplRefl)
void BinCount(TH1F *h, TF1 *fB, Int_t rebin, Double_t minMass, Double_t maxMass, Double_t &count, Double_t &ecount) const
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
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
TH1F ** fHistoSigmaTrial
histo with yield from subsamples of trials
TList * fitter
Definition: DrawAnaELoss.C:26
Bool_t fUsePowLawBkg
switch for pol5 background
Bool_t fUseFreeS
switch for FixSigDownFreeMean
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
Int_t fFitOption
name to characterize analysis case
Bool_t fUsePol5Bkg
switch for pol4 background
Bool_t fSaveBkgVal
flag to fix the width of the 2nd peak
Bool_t fUseFixedMeanFreeS
switch for FreeSigma
Double_t GetSigmaUncertainty() const
void Significance(Double_t nOfSigma, Double_t &significance, Double_t &errsignificance) const
Bool_t fUseLinBkg
switch for exponential background
Double_t GetRawYieldError() const
Bool_t fUseFixSigDownFreeMean
switch for FixSigUpFreeMean
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
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.
Double_t minMass
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 * fHistoRawYieldDistBinCAll
histo with bkg in mass bin edges from all trials
TH1F ** fHistoMeanTrial
histo with gauss sigma from subsamples of trials
Double_t maxMass
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
Double_t minMassForFit
Double_t GetRawYield() const
Double_t maxMassForFit