AliPhysics  ff1d528 (ff1d528)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
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==1) fitter->SetUseChi2Fit();
264  xnt[0]=rebin;
265  xnt[1]=iFirstBin;
266  xnt[2]=minMassForFit;
267  xnt[3]=maxMassForFit;
268  xnt[4]=typeb;
269  xnt[6]=0;
270  if(igs==kFixSigFreeMean){
271  fitter->SetFixGaussianSigma(fSigmaGausMC,kTRUE);
272  xnt[5]=1;
273  }else if(igs==kFixSigUpFreeMean){
275  xnt[5]=2;
276  }else if(igs==kFixSigDownFreeMean){
278  xnt[5]=3;
279  }else if(igs==kFreeSigFreeMean){
280  xnt[5]=0;
281  }else if(igs==kFixSigFixMean){
282  fitter->SetFixGaussianSigma(fSigmaGausMC,kTRUE);
283  fitter->SetFixGaussianMean(fMassD,kTRUE);
284  xnt[5]=1;
285  xnt[6]=1;
286  }else if(igs==kFreeSigFixMean){
287  fitter->SetFixGaussianMean(fMassD,kTRUE);
288  xnt[5]=0;
289  xnt[6]=1;
290  }
291  Bool_t out=kFALSE;
292  Double_t chisq=-1.;
293  Double_t sigma=0.;
294  Double_t esigma=0.;
295  Double_t pos=.0;
296  Double_t epos=.0;
297  Double_t ry=.0;
298  Double_t ery=.0;
299  Double_t significance=0.;
300  Double_t erSignif=0.;
301  Double_t bkg=0.;
302  Double_t erbkg=0.;
303  Double_t bkgBEdge=0;
304  Double_t erbkgBEdge=0;
305  TF1* fB1=0x0;
306  if(typeb<kNBkgFuncCases){
307  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);
308  out=fitter->MassFitter(0);
309  chisq=fitter->GetReducedChiSquare();
310  fitter->Significance(fnSigmaForBkgEval,significance,erSignif);
311  sigma=fitter->GetSigma();
312  pos=fitter->GetMean();
313  esigma=fitter->GetSigmaUncertainty();
314  if(esigma<0.00001) esigma=0.0001;
315  epos=fitter->GetMeanUncertainty();
316  if(epos<0.00001) epos=0.0001;
317  ry=fitter->GetRawYield();
318  ery=fitter->GetRawYieldError();
319  fB1=fitter->GetBackgroundFullRangeFunc();
320  fitter->Background(fnSigmaForBkgEval,bkg,erbkg);
321  Double_t minval = hInvMassHisto->GetXaxis()->GetBinLowEdge(hInvMassHisto->FindBin(pos-fnSigmaForBkgEval*sigma));
322  Double_t maxval = hInvMassHisto->GetXaxis()->GetBinUpEdge(hInvMassHisto->FindBin(pos+fnSigmaForBkgEval*sigma));
323  fitter->Background(minval,maxval,bkgBEdge,erbkgBEdge);
324  if(out && fDrawIndividualFits && thePad){
325  thePad->Clear();
326  fitter->DrawHere(thePad, fnSigmaForBkgEval);
327  fMassFitters.push_back(fitter);
328  mustDeleteFitter = kFALSE;
329  for (auto format : fInvMassFitSaveAsFormats) {
330  thePad->SaveAs(Form("FitOutput_%s_Trial%d.%s",hInvMassHisto->GetName(),globBin, format.c_str()));
331  }
332  }
333  }
334  // else{
335  // out=DoFitWithPol3Bkg(hRebinned,hmin,hmax,igs);
336  // if(out && thePad){
337  // thePad->Clear();
338  // hRebinned->Draw();
339  // TF1* fSB=(TF1*)hRebinned->GetListOfFunctions()->FindObject("fSB");
340  // fB1=new TF1("fB1","[0]+[1]*x+[2]*x*x+[3]*x*x*x",hmin,hmax);
341  // for(Int_t j=0; j<4; j++) fB1->SetParameter(j,fSB->GetParameter(3+j));
342  // fB1->SetLineColor(2);
343  // fB1->Draw("same");
344  // fSB->SetLineColor(4);
345  // fSB->Draw("same");
346  // thePad->Update();
347  // chisq=fSB->GetChisquare()/fSB->GetNDF();;
348  // sigma=fSB->GetParameter(2);
349  // esigma=fSB->GetParError(2);
350  // if(esigma<0.00001) esigma=0.0001;
351  // pos=fSB->GetParameter(1);
352  // epos=fSB->GetParError(1);
353  // if(epos<0.00001) epos=0.0001;
354  // ry=fSB->GetParameter(0)/hRebinned->GetBinWidth(1);
355  // ery=fSB->GetParError(0)/hRebinned->GetBinWidth(1);
356  // }
357  // }
358  xnt[7]=chisq;
359  if(out && chisq>0. && sigma>0.5*fSigmaGausMC && sigma<2.0*fSigmaGausMC){
360  xnt[8]=significance;
361  xnt[9]=pos;
362  xnt[10]=epos;
363  xnt[11]=sigma;
364  xnt[12]=esigma;
365  xnt[13]=ry;
366  xnt[14]=ery;
367  fHistoRawYieldDistAll->Fill(ry);
368  fHistoRawYieldTrialAll->SetBinContent(globBin,ry);
369  fHistoRawYieldTrialAll->SetBinError(globBin,ery);
370  fHistoSigmaTrialAll->SetBinContent(globBin,sigma);
371  fHistoSigmaTrialAll->SetBinError(globBin,esigma);
372  fHistoMeanTrialAll->SetBinContent(globBin,pos);
373  fHistoMeanTrialAll->SetBinError(globBin,epos);
374  fHistoChi2TrialAll->SetBinContent(globBin,chisq);
375  fHistoChi2TrialAll->SetBinError(globBin,0.00001);
376  fHistoSignifTrialAll->SetBinContent(globBin,significance);
377  fHistoSignifTrialAll->SetBinError(globBin,erSignif);
378  if(fSaveBkgVal) {
379  fHistoBkgTrialAll->SetBinContent(globBin,bkg);
380  fHistoBkgTrialAll->SetBinError(globBin,erbkg);
381  fHistoBkgInBinEdgesTrialAll->SetBinContent(globBin,bkgBEdge);
382  fHistoBkgInBinEdgesTrialAll->SetBinError(globBin,erbkgBEdge);
383  }
384 
385  if(ry<fMinYieldGlob) fMinYieldGlob=ry;
386  if(ry>fMaxYieldGlob) fMaxYieldGlob=ry;
387  fHistoRawYieldDist[theCase]->Fill(ry);
388  fHistoRawYieldTrial[theCase]->SetBinContent(itrial,ry);
389  fHistoRawYieldTrial[theCase]->SetBinError(itrial,ery);
390  fHistoSigmaTrial[theCase]->SetBinContent(itrial,sigma);
391  fHistoSigmaTrial[theCase]->SetBinError(itrial,esigma);
392  fHistoMeanTrial[theCase]->SetBinContent(itrial,pos);
393  fHistoMeanTrial[theCase]->SetBinError(itrial,epos);
394  fHistoChi2Trial[theCase]->SetBinContent(itrial,chisq);
395  fHistoChi2Trial[theCase]->SetBinError(itrial,0.00001);
396  fHistoSignifTrial[theCase]->SetBinContent(itrial,significance);
397  fHistoSignifTrial[theCase]->SetBinError(itrial,erSignif);
398  if(fSaveBkgVal) {
399  fHistoBkgTrial[theCase]->SetBinContent(itrial,bkg);
400  fHistoBkgTrial[theCase]->SetBinError(itrial,erbkg);
401  fHistoBkgInBinEdgesTrial[theCase]->SetBinContent(itrial,bkgBEdge);
402  fHistoBkgInBinEdgesTrial[theCase]->SetBinError(itrial,erbkgBEdge);
403  }
404 
405  for(Int_t iStepBC=0; iStepBC<fNumOfnSigmaBinCSteps; iStepBC++){
406  Double_t minMassBC=fMassD-fnSigmaBinCSteps[iStepBC]*sigma;
407  Double_t maxMassBC=fMassD+fnSigmaBinCSteps[iStepBC]*sigma;
408  if(minMassBC>minMassForFit &&
409  maxMassBC<maxMassForFit &&
410  minMassBC>(hRebinned->GetXaxis()->GetXmin()) &&
411  maxMassBC<(hRebinned->GetXaxis()->GetXmax())){
412  Double_t cnts,ecnts;
413  BinCount(hRebinned,fB1,1,minMassBC,maxMassBC,cnts,ecnts);
414  ++itrialBC;
415  fHistoRawYieldDistBinCAll->Fill(cnts);
416  fHistoRawYieldTrialBinCAll->SetBinContent(globBin,iStepBC+1,cnts);
417  fHistoRawYieldTrialBinCAll->SetBinError(globBin,iStepBC+1,ecnts);
418  fHistoRawYieldTrialBinC[theCase]->SetBinContent(itrial,iStepBC+1,cnts);
419  fHistoRawYieldTrialBinC[theCase]->SetBinError(itrial,iStepBC+1,ecnts);
420  fHistoRawYieldDistBinC[theCase]->Fill(cnts);
421  }
422  }
423  }
424  if (mustDeleteFitter) delete fitter;
425  fNtupleMultiTrials->Fill(xnt);
426  }
427  }
428  }
429  }
430  delete hRebinned;
431  }
432  }
433  return kTRUE;
434 }
435 
436 //________________________________________________________________________
438  // save histos in a root file for further analysis
439  const Int_t nCases=kNBkgFuncCases*kNFitConfCases;
440  TFile outHistos(fileName.Data(),option.Data());
441  if (outHistos.IsZombie()) {
442  Printf("Could not open file '%s'!", fileName.Data());
443  return;
444  }
445  outHistos.cd();
446  fHistoRawYieldTrialAll->Write();
447  fHistoSigmaTrialAll->Write();
448  fHistoMeanTrialAll->Write();
449  fHistoChi2TrialAll->Write();
450  fHistoSignifTrialAll->Write();
451  if(fSaveBkgVal) {
452  fHistoBkgTrialAll->Write();
453  fHistoBkgInBinEdgesTrialAll->Write();
454  }
455  fHistoRawYieldDistBinCAll->Write();
456  fHistoRawYieldTrialBinCAll->Write();
457  for(Int_t ic=0; ic<nCases; ic++){
458  fHistoRawYieldTrial[ic]->Write();
459  fHistoSigmaTrial[ic]->Write();
460  fHistoMeanTrial[ic]->Write();
461  fHistoChi2Trial[ic]->Write();
462  fHistoSignifTrial[ic]->Write();
463  if(fSaveBkgVal) {
464  fHistoBkgTrial[ic]->Write();
465  fHistoBkgInBinEdgesTrial[ic]->Write();
466  }
467  fHistoRawYieldTrialBinC[ic]->Write();
468  fHistoRawYieldDistBinC[ic]->Write();
469  }
470  fNtupleMultiTrials->SetDirectory(&outHistos);
471  fNtupleMultiTrials->Write();
472  outHistos.Close();
473 }
474 
475 //________________________________________________________________________
476 void AliHFMultiTrials::DrawHistos(TCanvas* cry) const{
477  // draw histos
478  cry->Divide(2,2);
479  cry->cd(1);
480  fHistoSigmaTrialAll->Draw();
481  cry->cd(3);
482  fHistoChi2TrialAll->Draw();
483  cry->cd(2);
484  fHistoRawYieldTrialAll->Draw();
485  cry->cd(4);
486  fHistoRawYieldDistAll->Draw();
487  TLatex* tmean=new TLatex(0.15,0.8,Form("mean=%.1f",fHistoRawYieldDistAll->GetMean()));
488  tmean->SetNDC();
489  tmean->Draw();
490  TLatex* thrms=new TLatex(0.15,0.72,Form("rms=%.1f",fHistoRawYieldDistAll->GetRMS()));
491  thrms->SetNDC();
492  thrms->Draw();
493  TLatex* tmax=new TLatex(0.6,0.8,Form("max=%.1f",fMaxYieldGlob));
494  tmax->SetNDC();
495  tmax->Draw();
496  TLatex* tmin=new TLatex(0.6,0.72,Form("min=%.1f",fMinYieldGlob));
497  tmin->SetNDC();
498  tmin->Draw();
499  TLatex* trms=new TLatex(0.6,0.64,Form("(max-min)/sqrt(12)=%.1f",(fMaxYieldGlob-fMinYieldGlob)/sqrt(12)));
500  trms->SetNDC();
501  trms->Draw();
502 
503 }
504 //________________________________________________________________________
505 TH1F* AliHFMultiTrials::RebinHisto(TH1D* hOrig, Int_t reb, Int_t firstUse) const{
506  // Rebin histogram, from bin firstUse to lastUse
507  // Use all bins if firstUse=-1
508 
509  Int_t nBinOrig=hOrig->GetNbinsX();
510  Int_t firstBinOrig=1;
511  Int_t lastBinOrig=nBinOrig;
512  Int_t nBinOrigUsed=nBinOrig;
513  Int_t nBinFinal=nBinOrig/reb;
514  if(firstUse>=1){
515  firstBinOrig=firstUse;
516  nBinFinal=(nBinOrig-firstUse+1)/reb;
517  nBinOrigUsed=nBinFinal*reb;
518  lastBinOrig=firstBinOrig+nBinOrigUsed-1;
519  }else{
520  Int_t exc=nBinOrigUsed%reb;
521  if(exc!=0){
522  nBinOrigUsed-=exc;
523  firstBinOrig+=exc/2;
524  lastBinOrig=firstBinOrig+nBinOrigUsed-1;
525  }
526  }
527 
528  printf("Rebin from %d bins to %d bins -- Used bins=%d in range %d-%d\n",nBinOrig,nBinFinal,nBinOrigUsed,firstBinOrig,lastBinOrig);
529  Float_t lowLim=hOrig->GetXaxis()->GetBinLowEdge(firstBinOrig);
530  Float_t hiLim=hOrig->GetXaxis()->GetBinUpEdge(lastBinOrig);
531  TH1F* hRebin=new TH1F(Form("%s-rebin%d_%d",hOrig->GetName(),reb,firstUse),hOrig->GetTitle(),nBinFinal,lowLim,hiLim);
532  Int_t lastSummed=firstBinOrig-1;
533  for(Int_t iBin=1;iBin<=nBinFinal; iBin++){
534  Float_t sum=0.;
535  Float_t sum2=0.;
536  for(Int_t iOrigBin=0;iOrigBin<reb;iOrigBin++){
537  sum+=hOrig->GetBinContent(lastSummed+1);
538  sum2+=hOrig->GetBinError(lastSummed+1)*hOrig->GetBinError(lastSummed+1);
539  lastSummed++;
540  }
541  hRebin->SetBinContent(iBin,sum);
542  hRebin->SetBinError(iBin,TMath::Sqrt(sum2));
543  }
544  return hRebin;
545 }
546 
547 //________________________________________________________________________
548 void AliHFMultiTrials::BinCount(TH1F* h, TF1* fB, Int_t rebin, Double_t minMass, Double_t maxMass, Double_t& count, Double_t& ecount) const{
549  // compute yield with bin couting
550  Int_t minBinSum=h->FindBin(minMass);
551  Int_t maxBinSum=h->FindBin(maxMass);
552  Double_t cntSig=0.;
553  Double_t cntErr=0.;
554  for(Int_t iMB=minBinSum; iMB<=maxBinSum; iMB++){
555  Double_t bkg=fB ? fB->Eval(h->GetBinCenter(iMB))/(Double_t)rebin : 0;
556  cntSig+=(h->GetBinContent(iMB)-bkg);
557  cntErr+=(h->GetBinError(iMB)*h->GetBinError(iMB));
558  }
559  count=cntSig;
560  ecount=TMath::Sqrt(cntErr);
561 }
562 //________________________________________________________________________
564  Int_t iCase){
565  //
566 
567  TH1F *hCutTmp=(TH1F*)histoToFit->Clone("hCutTmp");
568  for(Int_t ib=1; ib<=hCutTmp->GetNbinsX(); ib++){
569  Double_t xc=hCutTmp->GetBinCenter(ib);
570  if(xc>(fMassD-5.*fSigmaGausMC) && xc<(fMassD+5.*fSigmaGausMC)){
571  hCutTmp->SetBinContent(ib,0.);
572  hCutTmp->SetBinError(ib,0.);
573  }
574  }
575 
576  hCutTmp->Fit("pol2","E0","",hmin,hmax);
577  TF1* f2=(TF1*)hCutTmp->GetListOfFunctions()->FindObject("pol2");
578  TF1* f3=new TF1("myPol3","pol3");
579  for(Int_t i=0; i<3;i++) f3->SetParameter(i,f2->GetParameter(i));
580  hCutTmp->Fit(f3,"E0","",hmin,hmax);
581  Double_t quickCount=0.;
582  for(Int_t ib=1; ib<=histoToFit->GetNbinsX(); ib++){
583  Double_t xc=hCutTmp->GetBinCenter(ib);
584  if(xc>(fMassD-3.*fSigmaGausMC) && xc<(fMassD+3.*fSigmaGausMC)){
585  quickCount+=(histoToFit->GetBinContent(ib)-f3->Eval(xc));
586  }
587  }
588  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");
589  fSB->SetParameter(0,quickCount);
590  fSB->SetParameter(1,fMassD);
591  fSB->SetParameter(2,fSigmaGausMC);
592  for(Int_t j=0; j<4; j++) fSB->SetParameter(j+3,f3->GetParameter(j));
593  if(iCase==0) fSB->FixParameter(2,fSigmaGausMC);
594  else if(iCase==1) fSB->FixParameter(2,fSigmaGausMC*(1.+fSigmaMCVariation));
595  else if(iCase==2) fSB->FixParameter(2,fSigmaGausMC*(1.-fSigmaMCVariation));
596  else if(iCase==4){
597  fSB->FixParameter(1,fMassD);
598  fSB->FixParameter(2,fSigmaGausMC);
599  } else if(iCase==5){
600  fSB->FixParameter(1,fMassD);
601  }
602  histoToFit->Fit(fSB,"ME0","",hmin,hmax);
603  // quality cuts
604  if(fSB->GetParError(0)<0.01*fSB->GetParameter(0)) return kFALSE;
605  if(fSB->GetParError(0)>0.6*fSB->GetParameter(0)) return kFALSE;
606 
607  delete hCutTmp;
608  return kTRUE;
609 }
610 //__________________________________________________________________________________
611 TH1F* AliHFMultiTrials::SetTemplateRefl(const TH1F *h) {
612  fhTemplRefl=(TH1F*)h->Clone("hTemplRefl");
613  return fhTemplRefl;
614 }
615 
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
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.
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
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
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
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