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