AliPhysics  e6d2b2b (e6d2b2b)
AliHFVnVsMassFitter.cxx
Go to the documentation of this file.
1 #include "AliHFVnVsMassFitter.h"
2 
3 #include <TROOT.h>
4 #include <TMath.h>
5 #include <TF1.h>
6 #include <TGraphErrors.h>
7 #include <TColor.h>
8 #include <TLegend.h>
9 #include <TList.h>
10 #include <TCanvas.h>
11 #include <TStyle.h>
12 #include <TVirtualPad.h>
13 #include <TDatabasePDG.h>
14 #include <TPaveText.h>
15 #include "Fit/BinData.h"
16 #include "HFitInterface.h"
17 #include <vector>
18 
19 #include "AliVertexingHFUtils.h"
20 
22 ClassImp(AliHFVnVsMassFitter);
24 
25 //________________________________________________________________
27  :TObject()
28  ,fMassHisto(0x0)
29  ,fVnVsMassHisto(0x0)
30  ,fMassSgnFuncType(kGaus)
31  ,fMassBkgFuncType(kExpo)
32  ,fVnBkgFuncType(kLin)
33  ,fMassFuncFromPrefit(0x0)
34  ,fMassBkgFunc(0x0)
35  ,fMassSgnFunc(0x0)
36  ,fMassTotFunc(0x0)
37  ,fVnBkgFuncSb(0x0)
38  ,fVnBkgFunc(0x0)
39  ,fVnTotFunc(0x0)
40  ,fMassFitter(0x0)
41  ,fMassMin(1.69)
42  ,fMassMax(2.05)
43  ,fVn(0.)
44  ,fVnUncertainty(0.)
45  ,fSigma(0.)
46  ,fSigmaUncertainty(0.)
47  ,fMean(0.)
48  ,fMeanUncertainty(0.)
49  ,fRawYield(0.)
50  ,fRawYieldUncertainty(0.)
51  ,fChiSquare(0.)
52  ,fNDF(0)
53  ,fProb(0.)
54  ,fNSigmaForSB(3.)
55  ,fSigmaInit(0.012)
56  ,fMeanInit(1.870)
57  ,fSigma2GausInit(0.012)
58  ,fFrac2GausInit(0.2)
59  ,fMeanFixedFromMassFit(kFALSE)
60  ,fSigmaFixedFromMassFit(kFALSE)
61  ,fSigma2GausFixedFromMassFit(kFALSE)
62  ,fFrac2GausFixedFromMassFit(kFALSE)
63  ,fMassParticle(1.870)
64  ,fNParsMassSgn(3)
65  ,fNParsMassBkg(2)
66  ,fNParsVnBkg(2)
67  ,fNParsVnSgn(1)
68  ,fNParsVnSecPeak(0)
69  ,fNParsVnRfl(0)
70  ,fSigmaFixed(0)
71  ,fMeanFixed(0)
72  ,fSigma2GausFixed(0)
73  ,fFrac2GausFixed(0)
74  ,fPolDegreeBkg(3)
75  ,fPolDegreeVnBkg(3)
76  ,fReflections(kFALSE)
77  ,fNParsRfl(0)
78  ,fRflOverSig(0.)
79  ,fFixRflOverSig(kFALSE)
80  ,fHistoTemplRfl(0x0)
81  ,fHistoTemplRflInit(0x0)
82  ,fMassRflFunc(0x0)
83  ,fMassBkgRflFunc(0x0)
84  ,fRflOpt("1gaus")
85  ,fMinRefl(0.)
86  ,fMaxRefl(0.)
87  ,fSmoothRfl(kFALSE)
88  ,fRawYieldHelp(0.)
89  ,fVnRflOpt(0)
90  ,fVnRflLimited(kFALSE)
91  ,fVnRflMin(-1.)
92  ,fVnRflMax(1.)
93  ,fSecondPeak(kFALSE)
94  ,fMassSecPeakFunc(0x0)
95  ,fNParsSec(0)
96  ,fSecMass(-999.)
97  ,fSecWidth(9999.)
98  ,fFixSecMass(kFALSE)
99  ,fFixSecWidth(kFALSE)
100  ,fDoSecondPeakVn(kFALSE)
101  ,fHarmonic(2) {
102 
103  //default constructor
104 }
105 
106 //________________________________________________________________
107 AliHFVnVsMassFitter::AliHFVnVsMassFitter(TH1F* hMass, TH1F* hvn, Double_t min, Double_t max, Int_t funcMassBkg, Int_t funcMassSgn, Int_t funcvnBkg)
108  :TObject()
109  ,fMassSgnFuncType(funcMassSgn)
110  ,fMassBkgFuncType(funcMassBkg)
111  ,fVnBkgFuncType(funcvnBkg)
112  ,fMassFuncFromPrefit(0x0)
113  ,fMassBkgFunc(0x0)
114  ,fMassSgnFunc(0x0)
115  ,fMassTotFunc(0x0)
116  ,fVnBkgFuncSb(0x0)
117  ,fVnBkgFunc(0x0)
118  ,fVnTotFunc(0x0)
119  ,fMassFitter(0x0)
120  ,fMassMin(min)
121  ,fMassMax(max)
122  ,fVn(0.)
123  ,fVnUncertainty(0.)
124  ,fSigma(0.)
125  ,fSigmaUncertainty(0.)
126  ,fMean(0.)
127  ,fMeanUncertainty(0.)
128  ,fRawYield(0.)
130  ,fChiSquare(0.)
131  ,fNDF(0)
132  ,fProb(0)
133  ,fNSigmaForSB(3.)
134  ,fSigmaInit(0.012)
135  ,fMeanInit(1.870)
136  ,fSigma2GausInit(0.012)
137  ,fFrac2GausInit(0.2)
138  ,fMeanFixedFromMassFit(kFALSE)
139  ,fSigmaFixedFromMassFit(kFALSE)
142  ,fMassParticle(1.870)
143  ,fNParsMassSgn(3)
144  ,fNParsMassBkg(2)
145  ,fNParsVnBkg(2)
146  ,fNParsVnSgn(1)
147  ,fNParsVnSecPeak(0)
148  ,fNParsVnRfl(0)
149  ,fSigmaFixed(0)
150  ,fMeanFixed(0)
151  ,fSigma2GausFixed(0)
152  ,fFrac2GausFixed(0)
153  ,fPolDegreeBkg(3)
154  ,fPolDegreeVnBkg(3)
155  ,fReflections(kFALSE)
156  ,fNParsRfl(0)
157  ,fRflOverSig(0.)
158  ,fFixRflOverSig(kFALSE)
159  ,fHistoTemplRfl(0x0)
160  ,fHistoTemplRflInit(0x0)
161  ,fMassRflFunc(0x0)
162  ,fMassBkgRflFunc(0x0)
163  ,fRflOpt("1gaus")
164  ,fMinRefl(0.)
165  ,fMaxRefl(0.)
166  ,fSmoothRfl(kFALSE)
167  ,fRawYieldHelp(0.)
168  ,fVnRflOpt(0)
169  ,fVnRflLimited(kFALSE)
170  ,fVnRflMin(-1.)
171  ,fVnRflMax(1.)
172  ,fSecondPeak(kFALSE)
173  ,fMassSecPeakFunc(0x0)
174  ,fNParsSec(0)
175  ,fSecMass(-999.)
176  ,fSecWidth(9999.)
177  ,fFixSecMass(kFALSE)
178  ,fFixSecWidth(kFALSE)
179  ,fDoSecondPeakVn(kFALSE)
180  ,fHarmonic(2) {
181 
182  //standard constructor
183  fMassHisto = (TH1F*)hMass->Clone("fHistoInvMass");
184  fMassHisto->SetDirectory(0);
185  fVnVsMassHisto = (TH1F*)hvn->Clone(Form("fHistoV%dVsMass",fHarmonic));
186  fVnVsMassHisto->SetDirectory(0);
187 
189 }
190 
191 //________________________________________________________________
193 
194  //destructor
195  if(fMassHisto) delete fMassHisto;
196  if(fVnVsMassHisto) delete fVnVsMassHisto;
198  if(fMassBkgFunc) delete fMassBkgFunc;
199  if(fMassBkgRflFunc) delete fMassBkgRflFunc;
200  if(fMassSgnFunc) delete fMassSgnFunc;
201  if(fMassTotFunc) delete fMassTotFunc;
202  if(fVnBkgFuncSb) delete fVnBkgFuncSb;
203  if(fVnBkgFunc) delete fVnBkgFunc;
204  if(fVnTotFunc) delete fVnTotFunc;
205  if(fMassFitter) delete fMassFitter;
206  if(fHistoTemplRfl) delete fHistoTemplRfl;
208  if(fMassRflFunc) delete fMassRflFunc;
210 }
211 
212 //________________________________________________________________
214 
215  if(!fMassHisto || !fVnVsMassHisto) {AliError("Histograms not set! Exit."); return kFALSE;}
217 
219  Int_t NvnParsSgn = 1;
220  if(fSecondPeak && fDoSecondPeakVn) {NvnParsSgn+=1;}
221  if(fReflections && fVnRflOpt==kFreePar) {NvnParsSgn+=1;}
222  const Int_t nparsvn = nparsmass+fNParsVnBkg+NvnParsSgn;
223 
224  Bool_t massprefit=MassPrefit();
225  if(!massprefit) {AliError("Impossible to perform the mass prefit"); return kFALSE;}
226  Bool_t vnprefit=VnSBPrefit();
227 
228  std::vector<Double_t> initpars;
229  for(Int_t iBkgPar=0; iBkgPar<fNParsMassBkg; iBkgPar++) {
230  initpars.push_back(fMassFuncFromPrefit->GetParameter(iBkgPar));
231  }
232  for(Int_t iSgnPar=0; iSgnPar<fNParsMassSgn; iSgnPar++) {
233  initpars.push_back(fMassFuncFromPrefit->GetParameter(iSgnPar+fNParsMassBkg));
234  }
235  for(Int_t iSecPeakPar=0; iSecPeakPar<fNParsSec; iSecPeakPar++) {
236  initpars.push_back(fMassFuncFromPrefit->GetParameter(iSecPeakPar+fNParsMassBkg+fNParsMassSgn));
237  }
238  for(Int_t iReflPar=0; iReflPar<fNParsRfl; iReflPar++) {
239  initpars.push_back(fMassFuncFromPrefit->GetParameter(iReflPar+fNParsMassBkg+fNParsMassSgn+fNParsSec));
240  }
241  for(Int_t iVnBkgPar=0; iVnBkgPar<fNParsVnBkg; iVnBkgPar++) {
242  if(vnprefit) {initpars.push_back(fVnBkgFuncSb->GetParameter(iVnBkgPar));}
243  else {initpars.push_back(0.05);}
244  }
245  initpars.push_back(0.10); //initial parameter for signal vn
246  if(fSecondPeak && fDoSecondPeakVn) {initpars.push_back(0.10);} //initial parameter for second peak vn
247 
248  fMassTotFunc = new TF1("fMassTotFunc",this,&AliHFVnVsMassFitter::MassFunc,fMassMin,fMassMax,nparsmass,"AliHFVnVsMassFitter","MassFunc");
249  fVnTotFunc = new TF1("fVnTotFunc",this,&AliHFVnVsMassFitter::vnFunc,fMassMin,fMassMax,nparsvn,"AliHFVnVsMassFitter","vnFunc");
250  SetParNames();
251 
252  ROOT::Math::WrappedMultiTF1 wfTotMass(*fMassTotFunc,1);
253  ROOT::Math::WrappedMultiTF1 wfTotVn(*fVnTotFunc,1);
254 
255  // set data options and ranges
256  ROOT::Fit::DataOptions opt;
257  ROOT::Fit::DataRange rangeMass; //same range for two functions
258 
259  rangeMass.SetRange(fMassMin,fMassMax);
260  ROOT::Fit::BinData dataMass(opt,rangeMass);
261  ROOT::Fit::FillData(dataMass, fMassHisto);
262  ROOT::Fit::BinData dataVn(opt,rangeMass);
263  ROOT::Fit::FillData(dataVn, fVnVsMassHisto);
264 
265  //define the 2 chi squares
266  ROOT::Fit::Chi2Function chi2Mass(dataMass, wfTotMass);
267  ROOT::Fit::Chi2Function chi2Vn(dataVn, wfTotVn);
268 
269  //define the global chi square
270  AliHFGlobalChi2 globalChi2(chi2Mass, chi2Vn);
271 
272  //define fitter
273  ROOT::Fit::Fitter fitter;
274  // create before the parameter settings in order to fix or set range on them
275  fitter.Config().SetParamsSettings(nparsvn,initpars.data()); //set initial parameters from prefits
276  if(fMeanFixed==2 || fMeanFixedFromMassFit) {fitter.Config().ParSettings(fNParsMassBkg+1).Fix();}
277  if(fSigmaFixed==2 || fSigmaFixedFromMassFit) {fitter.Config().ParSettings(fNParsMassBkg+2).Fix();}
278  if(fMassSgnFuncType==k2Gaus) {
279  if(fFrac2GausFixed==2 || fFrac2GausFixedFromMassFit) {fitter.Config().ParSettings(fNParsMassBkg+3).Fix();}
280  if(fSigma2GausFixed==2 || fSigma2GausFixedFromMassFit) {fitter.Config().ParSettings(fNParsMassBkg+4).Fix();}
281  }
282  if(fSecondPeak) {
283  if(fFixSecMass) {fitter.Config().ParSettings(fNParsMassBkg+fNParsMassSgn+1).Fix();}
284  if(fFixSecWidth) {fitter.Config().ParSettings(fNParsMassBkg+fNParsMassSgn+2).Fix();}
285  }
286  if(fReflections) {
287  if(fFixRflOverSig) fitter.Config().ParSettings(fNParsMassBkg+fNParsMassSgn+fNParsSec).Fix();
288  if(fVnRflLimited) fitter.Config().ParSettings(nparsmass+fNParsVnBkg+NvnParsSgn-1).SetLimits(fVnRflMin,fVnRflMax);
289  }
290 
291  fitter.Config().MinimizerOptions().SetPrintLevel(0);
292  fitter.Config().SetMinimizer("Minuit2","Migrad");
293  for(Int_t iPar=0; iPar<nparsvn; iPar++) {fitter.Config().ParSettings(iPar).SetName(fVnTotFunc->GetParName(iPar));}
294  // fit FCN function directly
295  // (specify optionally data size and flag to indicate that is a chi2 fit
296  fitter.FitFCN(nparsvn,globalChi2,0,dataMass.Size()+dataVn.Size(),kFALSE);
297  ROOT::Fit::FitResult result = fitter.Result();
298  result.Print(std::cout);
299 
300  //set parameters in every function
301  fVnBkgFunc = new TF1("fVnBkgFunc",this,&AliHFVnVsMassFitter::vnBkgFunc,fMassMin,fMassMax,fNParsVnBkg,"AliHFVnVsMassFitter","vnBkgFunc");
302  fMassBkgFunc = new TF1("fMassBkgFunc",this,&AliHFVnVsMassFitter::MassBkg,fMassMin,fMassMax,fNParsMassBkg,"AliHFVnVsMassFitter","MassBkg");
303  fMassSgnFunc = new TF1("fMassSgnFunc",this,&AliHFVnVsMassFitter::MassSignal,fMassMin,fMassMax,fNParsMassSgn,"AliHFVnVsMassFitter","MassSignal");
304  if(fReflections) {fMassRflFunc = new TF1("fMassRflFunc",this,&AliHFVnVsMassFitter::MassRfl,fMassMin,fMassMax,fNParsRfl,"AliHFVnVsMassFitter","MassRfl");}
305  if(fReflections) {fMassBkgRflFunc = new TF1("fMassBkgRflFunc",this,&AliHFVnVsMassFitter::MassBkgRfl,fMassMin,fMassMax,fNParsMassBkg+fNParsRfl,"AliHFVnVsMassFitter","MassBkgRfl");}
306  if(fSecondPeak) {fMassSecPeakFunc = new TF1("fMassSecPeakFunc",this,&AliHFVnVsMassFitter::MassSecondPeak,fMassMin,fMassMax,fNParsSec,"AliHFVnVsMassFitter","MassSecondPeak");}
307  for(Int_t iPar=0; iPar<nparsvn; iPar++) {
308  fVnTotFunc->SetParameter(iPar,result.Parameter(iPar));
309  fVnTotFunc->SetParError(iPar,result.ParError(iPar));
310  if(iPar<nparsmass) {
311  fMassTotFunc->SetParameter(iPar,result.Parameter(iPar));
312  fMassTotFunc->SetParError(iPar,result.ParError(iPar));
313  }
314  if(iPar>=nparsmass && iPar<nparsvn-NvnParsSgn) {
315  fVnBkgFunc->SetParameter(iPar-nparsmass,result.Parameter(iPar));
316  fVnBkgFunc->SetParError(iPar-nparsmass,result.ParError(iPar));
317  }
318  if(iPar>=fNParsMassBkg && iPar<fNParsMassBkg+fNParsMassSgn) {
319  fMassSgnFunc->SetParameter(iPar-fNParsMassBkg,result.Parameter(iPar));
320  }
321  if(iPar<fNParsMassBkg) {
322  fMassBkgFunc->SetParameter(iPar,result.Parameter(iPar));
323  fMassBkgFunc->SetParError(iPar,result.ParError(iPar));
324  if(fReflections) {
325  fMassBkgRflFunc->SetParameter(iPar,result.Parameter(iPar));
326  fMassBkgRflFunc->SetParError(iPar,result.ParError(iPar));
327  }
328  }
329  if(fReflections && (iPar>=fNParsMassBkg+fNParsMassSgn+fNParsSec && iPar<fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl)) {
330  fMassRflFunc->SetParameter(iPar-(fNParsMassBkg+fNParsMassSgn+fNParsSec),result.Parameter(iPar));
331  fMassRflFunc->SetParError(iPar-(fNParsMassBkg+fNParsMassSgn+fNParsSec),result.ParError(iPar));
332  fMassBkgRflFunc->SetParameter(iPar-(fNParsMassSgn+fNParsSec),result.Parameter(iPar));
333  fMassBkgRflFunc->SetParError(iPar-(fNParsMassSgn+fNParsSec),result.ParError(iPar));
334  }
335  if(fSecondPeak && (iPar>=fNParsMassBkg+fNParsMassSgn && iPar<fNParsMassBkg+fNParsMassSgn+fNParsSec)) {
336  fMassSecPeakFunc->SetParameter(iPar-(fNParsMassBkg+fNParsMassSgn),result.Parameter(iPar));
337  fMassSecPeakFunc->SetParError(iPar-(fNParsMassBkg+fNParsMassSgn),result.ParError(iPar));
338  }
339  }
340  if(drawFit) {DrawFit();}
341 
342  fVn = fVnTotFunc->GetParameter(fVnTotFunc->GetNpar()-NvnParsSgn);
343  fVnUncertainty = fVnTotFunc->GetParError(fVnTotFunc->GetNpar()-NvnParsSgn);
344  if(fDoSecondPeakVn) {
345  fVnSecPeak = fVnTotFunc->GetParameter(fVnTotFunc->GetNpar()-1);
346  fVnSecPeakUncertainty = fVnTotFunc->GetParError(fVnTotFunc->GetNpar()-1);
347  }
348  fRawYield = fVnTotFunc->GetParameter(fNParsMassBkg)/fMassHisto->GetBinWidth(10);
349  fRawYieldUncertainty = fVnTotFunc->GetParError(fNParsMassBkg)/fMassHisto->GetBinWidth(10);
350  fMean = fVnTotFunc->GetParameter(fNParsMassBkg+1);
351  fMeanUncertainty = fVnTotFunc->GetParError(fNParsMassBkg+1);
352  fSigma = fVnTotFunc->GetParameter(fNParsMassBkg+2);
353  fSigmaUncertainty = fVnTotFunc->GetParError(fNParsMassBkg+2);
354  fChiSquare = result.MinFcnValue();
355  fNDF = result.Ndf();
356  fProb = result.Prob();
357 
358  return kTRUE;
359 }
360 
361 //______________________________________________________________________________
362 void AliHFVnVsMassFitter::DrawHere(TVirtualPad* c){
364 
365  gStyle->SetOptStat(0);
366  gStyle->SetCanvasColor(0);
367  gStyle->SetFrameFillColor(0);
368  c->Divide(1,2);
369 
370  c->cd(1);
371  fMassHisto->SetTitle("");
372  fMassHisto->SetMarkerStyle(20);
373  fMassHisto->SetMarkerSize(1);
374  fMassHisto->SetMarkerColor(kBlack);
375  fMassHisto->SetLineColor(kBlack);
376  fMassHisto->GetYaxis()->SetRangeUser(0.,fMassHisto->GetMaximum()*1.2);
377  fMassHisto->GetXaxis()->SetRangeUser(fMassMin,fMassMax);
378  fMassHisto->Draw("E");
379  if(fMassFuncFromPrefit) {
380  fMassFuncFromPrefit->SetLineColor(kGray+1);
381  fMassFuncFromPrefit->SetLineStyle(7);
383  fMassFuncFromPrefit->Draw("same");
384  }
385  if(fMassBkgFunc) {
386  fMassBkgFunc->SetLineColor(kRed);
387  fMassBkgFunc->SetRange(fMassMin,fMassMax);
388  fMassBkgFunc->Draw("same");
389  }
390  if(fMassRflFunc) {
391  fMassRflFunc->SetLineColor(kGreen+2);
392  fMassRflFunc->SetRange(fMassMin,fMassMax);
393  fMassRflFunc->Draw("same");
394  }
395  if(fMassBkgRflFunc) {
396  fMassBkgRflFunc->SetLineColor(kOrange+7);
397  fMassBkgRflFunc->SetLineStyle(9);
398  fMassBkgRflFunc->SetRange(fMassMin,fMassMax);
399  fMassBkgRflFunc->Draw("same");
400  }
401  if(fMassSecPeakFunc) {
402  fMassSecPeakFunc->SetLineColor(kMagenta+1);
403  fMassSecPeakFunc->SetLineStyle(7);
405  fMassSecPeakFunc->Draw("same");
406  }
407  if(fMassTotFunc) {
408  fMassTotFunc->SetLineColor(kBlue);
409  fMassTotFunc->SetRange(fMassMin,fMassMax);
410  fMassTotFunc->Draw("same");
411  }
412  TPaveText* massinfo = new TPaveText(0.45,0.7,1.,0.87,"NDC");
413  massinfo->SetTextFont(42);
414  massinfo->SetTextSize(0.05);
415  massinfo->SetBorderSize(0);
416  massinfo->SetFillStyle(0);
417  massinfo->SetTextColor(kBlue);
418  massinfo->AddText(Form("mean = %.3f #pm %.3f",fVnTotFunc->GetParameter(fNParsMassBkg+1),fVnTotFunc->GetParError(fNParsMassBkg+1)));
419  massinfo->AddText(Form("sigma = %.3f #pm %.3f",fVnTotFunc->GetParameter(fNParsMassBkg+2),fVnTotFunc->GetParError(fNParsMassBkg+2)));
420  if(fMassSgnFuncType==k2Gaus) {
421  massinfo->AddText(Form("sigma2 = %.3f #pm %.3f",fVnTotFunc->GetParameter(fNParsMassBkg+3),fVnTotFunc->GetParError(fNParsMassBkg+3)));
422  }
423  massinfo->Draw("same");
424 
425  c->cd(2);
426  fVnVsMassHisto->SetTitle("");
427  fVnVsMassHisto->SetMarkerStyle(20);
428  fVnVsMassHisto->SetMarkerSize(1);
429  fVnVsMassHisto->SetMarkerColor(kBlack);
430  fVnVsMassHisto->SetLineColor(kBlack);
431  fVnVsMassHisto->GetYaxis()->SetRangeUser(fVnVsMassHisto->GetMinimum()-0.15,fVnVsMassHisto->GetMaximum()+0.20);
432  fVnVsMassHisto->GetXaxis()->SetRangeUser(fMassMin,fMassMax);
433  fVnBkgFuncSb->SetRange(fMassMin,fMassMax);
434  fVnVsMassHisto->Draw("E");
435  if(fVnBkgFuncSb) {
436  fVnBkgFuncSb->SetLineColor(kGray+1);
437  fVnBkgFuncSb->SetLineStyle(7);
438  fVnBkgFuncSb->Draw("same");
439  }
440  if(fVnBkgFunc) {
441  fVnBkgFunc->SetLineColor(kRed);
442  fVnBkgFunc->SetRange(fMassMin,fMassMax);
443  fVnBkgFunc->Draw("same");
444  }
445  if(fVnTotFunc) {
446  fVnTotFunc->SetLineColor(kBlue);
447  fVnTotFunc->SetRange(fMassMin,fMassMax);
448  fVnTotFunc->Draw("same");
449  }
450 
451  TPaveText* vninfo = new TPaveText(-0.45,0.7,1.,0.87,"NDC");
452  vninfo->SetTextFont(42);
453  vninfo->SetTextSize(0.05);
454  vninfo->SetBorderSize(0);
455  vninfo->SetFillStyle(0);
456  Int_t NvnParsSgn = 1;
457  if(fSecondPeak && fDoSecondPeakVn) {NvnParsSgn+=1;}
458  if(fReflections && fVnRflOpt==kFreePar) {NvnParsSgn+=1;}
459  vninfo->AddText(Form("#it{v}_{%d}^{sgn} = %.3f #pm %.3f",fHarmonic,fVnTotFunc->GetParameter(fVnTotFunc->GetNpar()-NvnParsSgn),fVnTotFunc->GetParError(fVnTotFunc->GetNpar()-NvnParsSgn)));
460  if(fSecondPeak && fDoSecondPeakVn) {vninfo->AddText(Form("#it{v}_{%d}^{sec peak} = %.3f #pm %.3f",fHarmonic,fVnTotFunc->GetParameter(fVnTotFunc->GetNpar()-1),fVnTotFunc->GetParError(fVnTotFunc->GetNpar()-1)));}
461  vninfo->AddText(Form("#chi^{2}/#it{ndf} = %.2f/%d",fChiSquare,fNDF));
462  vninfo->Draw("same");
463 
464  c->Update();
465 }
466 
467 //________________________________________________________________
469 
470  //define proper maxs and mins from histos
471  Double_t tmpmin = TMath::Max(fMassHisto->GetBinLowEdge(1),fVnVsMassHisto->GetBinLowEdge(1));
472  fMassMin=TMath::Max(fMassMin,tmpmin);
473  Double_t tmpmax = TMath::Min(fMassHisto->GetBinLowEdge(fMassHisto->GetNbinsX())+fMassHisto->GetBinWidth(fMassHisto->GetNbinsX()),fVnVsMassHisto->GetBinLowEdge(fVnVsMassHisto->GetNbinsX())+fVnVsMassHisto->GetBinWidth(fVnVsMassHisto->GetNbinsX()));
474  fMassMax=TMath::Min(fMassMax,tmpmax);
475 
481  if(fMassSgnFuncType==k2Gaus) {
486  }
490  if(fReflections) {
494  }
495  Bool_t status = fMassFitter->MassFitter(kFALSE);
496 
497  if(status) {
499  fMassFuncFromPrefit->SetName("fMassFuncFromPrefit");
500  }
502 
503  return status;
504 }
505 
506 //________________________________________________________________
508 
509  Double_t mean = fMassFitter->GetMean();
511  const Int_t nMassBins = fVnVsMassHisto->GetNbinsX();
512  Double_t SBbins[nMassBins];
513  Int_t nSBbins=0;
514  for(Int_t iBin=0; iBin<nMassBins; iBin++) {
515  Double_t min = fVnVsMassHisto->GetBinLowEdge(iBin+1);
516  Double_t max = fVnVsMassHisto->GetBinLowEdge(iBin+1)+fVnVsMassHisto->GetBinWidth(iBin+1);
517  if(max<mean){
518  if(max<(mean-fNSigmaForSB*sigma)) {SBbins[iBin]=1; nSBbins++;}
519  else {SBbins[iBin]=0;}
520  }
521  if(min>=mean){
522  if(min>(mean+fNSigmaForSB*sigma)) {SBbins[iBin]=1; nSBbins++;}
523  else {SBbins[iBin]=0;}
524  }
525  }
526  TGraphErrors* gVnVsMassSB = new TGraphErrors(nSBbins);
527  for(Int_t iBin=0; iBin<nMassBins; iBin++) {
528  if(SBbins[iBin]==1) {
529  gVnVsMassSB->SetPoint(iBin,fVnVsMassHisto->GetBinCenter(iBin+1),fVnVsMassHisto->GetBinContent(iBin+1));
530  gVnVsMassSB->SetPointError(iBin,fVnVsMassHisto->GetBinWidth(iBin+1)/2,fVnVsMassHisto->GetBinError(iBin+1));
531  }
532  }
533  fVnBkgFuncSb = new TF1("fVnBkgFuncSb",this,&AliHFVnVsMassFitter::vnBkgFunc,fMassMin,fMassMax,fNParsVnBkg,"AliHFVnVsMassFitter","vnBkgFunc");
534  switch(fVnBkgFuncType) {
535  case 1:
536  fVnBkgFuncSb->SetParName(0,"ConstVnBkg");
537  fVnBkgFuncSb->SetParName(1,"SlopeVnBkg");
538  break;
539  case 2:
540  fVnBkgFuncSb->SetParName(0,"ConstVnBkg");
541  fVnBkgFuncSb->SetParName(1,"Coef1VnBkg");
542  fVnBkgFuncSb->SetParName(2,"Coef2VnBkg");
543  break;
544  default:
545  AliError("Error in setting signal par names: check fVnBkgFuncType");
546  break;
547  }
548  gVnVsMassSB->Fit(fVnBkgFuncSb,"","",fMassMin,fMassMax);
549  Bool_t status=kFALSE;
550  if(fVnBkgFuncSb->GetChisquare()<1000) status=kTRUE;
551 
552  delete gVnVsMassSB;
553  return status;
554 }
555 
556 //________________________________________________________________
558 
559  switch(fMassSgnFuncType) {
560  case 0: //single gaus
561  fNParsMassSgn=3;
562  break;
563  case 1: //double gaus
564  fNParsMassSgn=5;
565  break;
566  default:
567  AliError("Error in computing fMassSgnFuncType: check fMassSgnFuncType");
568  break;
569  }
570 
571  switch(fMassBkgFuncType) {
572  case 0: //expo
573  fNParsMassBkg=2;
574  break;
575  case 1: //lin
576  fNParsMassBkg=2;
577  break;
578  case 2: //pol2
579  fNParsMassBkg=3;
580  break;
581  case 3: //no bkg
582  fNParsMassBkg=1;
583  break;
584  case 4: //power law
585  fNParsMassBkg=2;
586  break;
587  case 5: //power expo
588  fNParsMassBkg=2;
589  break;
590  case 6: //high degree pol
592  break;
593  default:
594  AliError("Error in computing fNParsMassBkg: check fMassBkgFuncType");
595  break;
596  }
597 
598  switch(fVnBkgFuncType) {
599  case 0: //expo
600  fNParsVnBkg=2;
601  break;
602  case 1: //lin
603  fNParsVnBkg=2;
604  break;
605  case 2: //pol2
606  fNParsVnBkg=3;
607  break;
608  case 6: //high degree pol
610  break;
611  default:
612  AliError("Error in computing fNParsVnBkg: check fVnBkgFuncType");
613  break;
614  }
615 
616  fNParsVnSgn=1;
617 
618  if(fReflections) {
619  fNParsRfl=1;
620  if(fVnRflOpt==3) fNParsVnRfl=1;
621  else fNParsVnRfl=0;
622  }
623  else {
624  fNParsRfl=0;
625  fNParsVnRfl=0;
626  }
627 
628  if(fSecondPeak) {
629  fNParsSec=3;
630  fNParsVnSecPeak=1;
631  }
632  else {
633  fNParsSec=0;
634  fNParsVnSecPeak=0;
635  }
636 }
637 
638 //________________________________________________________________
640 
641  switch(fMassSgnFuncType) {
642  case 0: //single gaus
643  fVnTotFunc->SetParName(fNParsMassBkg,"SgnInt");
644  fVnTotFunc->SetParName(fNParsMassBkg+1,"Mean");
645  fVnTotFunc->SetParName(fNParsMassBkg+2,"Sigma");
646  break;
647  case 1: //double gaus
648  fVnTotFunc->SetParName(fNParsMassBkg,"SgnInt");
649  fVnTotFunc->SetParName(fNParsMassBkg+1,"Mean");
650  fVnTotFunc->SetParName(fNParsMassBkg+2,"Sigma1");
651  fVnTotFunc->SetParName(fNParsMassBkg+3,"Frac");
652  fVnTotFunc->SetParName(fNParsMassBkg+4,"Sigma2");
653  break;
654  default:
655  AliError("Error in setting signal par names: check fMassSgnFuncType");
656  break;
657  }
658 
659  switch(fMassBkgFuncType) {
660  case 0: //expo
661  fVnTotFunc->SetParName(0,"BkgInt");
662  fVnTotFunc->SetParName(1,"Slope");
663  break;
664  case 1: //lin
665  fVnTotFunc->SetParName(0,"BkgInt");
666  fVnTotFunc->SetParName(1,"Slope");
667  break;
668  case 2: //pol2
669  fVnTotFunc->SetParName(0,"BkgInt");
670  fVnTotFunc->SetParName(1,"Coef1");
671  fVnTotFunc->SetParName(2,"Coef1");
672  break;
673  case 3: //no bkg
674  fVnTotFunc->SetParName(0,"Const");
675  break;
676  case 4: //power law
677  fVnTotFunc->SetParName(0,"BkgInt");
678  fVnTotFunc->SetParName(1,"Coef1");
679  break;
680  case 5: //power expo
681  fVnTotFunc->SetParName(0,"Coef1");
682  fVnTotFunc->SetParName(1,"Coef2");
683  break;
684  case 6: //high degree pol
685  fVnTotFunc->SetParName(0,"BkgInt");
686  for(Int_t iPar=1; iPar<fNParsMassBkg; iPar++) {fVnTotFunc->SetParName(iPar,Form("Coef%d",iPar));}
687  break;
688  default:
689  AliError("Error in setting signal par names: check fMassBkgFuncType");
690  break;
691  }
692 
693  for(Int_t iPar=0; iPar<fNParsVnBkg; iPar++) {fVnTotFunc->SetParName(fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl+iPar,fVnBkgFuncSb->GetParName(iPar));}
694 
695  if(fReflections) {fVnTotFunc->SetParName(fNParsMassBkg+fNParsMassSgn+fNParsSec,"ReflOverS");}
696 
697  if(fSecondPeak) {
698  fVnTotFunc->SetParName(fNParsMassBkg+fNParsMassSgn,"SecPeakInt");
699  fVnTotFunc->SetParName(fNParsMassBkg+fNParsMassSgn+1,"SecPeakMean");
700  fVnTotFunc->SetParName(fNParsMassBkg+fNParsMassSgn+2,"SecPeakSigma");
701  }
702  fVnTotFunc->SetParName(fNParsMassBkg+fNParsMassSgn+fNParsRfl+fNParsSec+fNParsVnBkg,Form("v%dSgn",fHarmonic));
703 
704  if(fSecondPeak && fDoSecondPeakVn) {fVnTotFunc->SetParName(fNParsMassBkg+fNParsMassSgn+fNParsRfl+fNParsSec+fNParsVnBkg+1,Form("v%dSecPeak",fHarmonic));}
706 }
707 
708 //_________________________________________________________________________
709 void AliHFVnVsMassFitter::Signal(Double_t nOfSigma,Double_t &signal,Double_t &errsignal) const {
712 
713  Double_t minMass=fMean-nOfSigma*fSigma;
714  Double_t maxMass=fMean+nOfSigma*fSigma;
715  Signal(minMass,maxMass,signal,errsignal);
716  return;
717 }
718 
719 //_________________________________________________________________________
720 void AliHFVnVsMassFitter::Signal(Double_t min, Double_t max, Double_t &signal,Double_t &errsignal) const {
723  if(!fMassSgnFunc) {signal=-1; errsignal=0; return;}
724 
725  signal=fMassSgnFunc->Integral(min, max)/(Double_t)fMassHisto->GetBinWidth(1);
726  errsignal=(fRawYieldUncertainty/fRawYield)*signal;/*assume relative error is the same as for total integral*/
727 
728  return;
729 }
730 
731 //___________________________________________________________________________
732 void AliHFVnVsMassFitter::Background(Double_t nOfSigma,Double_t &background,Double_t &errbackground) const {
735 
736  Double_t minMass=fMean-nOfSigma*fSigma;
737  Double_t maxMass=fMean+nOfSigma*fSigma;
738  Background(minMass,maxMass,background,errbackground);
739 
740  return;
741 }
742 
743 //___________________________________________________________________________
744 void AliHFVnVsMassFitter::Background(Double_t min, Double_t max, Double_t &background,Double_t &errbackground) const {
747 
748  if(!fMassBkgFunc) {background=-1; errbackground=0; return;}
749 
750  Double_t intB=fMassBkgFunc->GetParameter(0);
751  Double_t intBerr=fMassBkgFunc->GetParError(0);
752  //relative error evaluation: from histo
753 
754  Int_t leftBand=fMassHisto->FindBin(fMean-4*fSigma);
755  Int_t rightBand=fMassHisto->FindBin(fMean+4*fSigma);
756  intB=fMassHisto->Integral(1,leftBand)+fMassHisto->Integral(rightBand,fMassHisto->GetNbinsX());
757  Double_t sum2=0;
758  for(Int_t iBin=1; iBin<=leftBand; iBin++){
759  sum2+=fMassHisto->GetBinError(iBin)*fMassHisto->GetBinError(iBin);
760  }
761  for(Int_t iBin=rightBand; iBin<=fMassHisto->GetNbinsX(); iBin++){
762  sum2+=fMassHisto->GetBinError(iBin)*fMassHisto->GetBinError(iBin);
763  }
764 
765  intBerr=TMath::Sqrt(sum2);
766 
767  background=fMassBkgFunc->Integral(min,max)/(Double_t)fMassHisto->GetBinWidth(1);
768  errbackground=intBerr/intB*background;
769 
770  return;
771 }
772 
773 //__________________________________________________________________________
774 
775 void AliHFVnVsMassFitter::Significance(Double_t nOfSigma,Double_t &significance,Double_t &errsignificance) const {
778 
779  Double_t minMass=fMean-nOfSigma*fSigma;
780  Double_t maxMass=fMean+nOfSigma*fSigma;
781  Significance(minMass, maxMass, significance, errsignificance);
782 
783  return;
784 }
785 
786 //__________________________________________________________________________
787 
788 void AliHFVnVsMassFitter::Significance(Double_t min, Double_t max, Double_t &significance,Double_t &errsignificance) const {
791 
792  Double_t background,errbackground;
793  Background(min,max,background,errbackground);
794 
795  if (fRawYield+background <= 0.){
796  significance=-1;
797  errsignificance=0;
798  return;
799  }
800 
801  AliVertexingHFUtils::ComputeSignificance(fRawYield,fRawYieldUncertainty,background,errbackground,significance,errsignificance);
802 
803  return;
804 }
805 
806 //________________________________________________________________
808 
809  return TMath::Gaus(x,mean,sigma,kTRUE);
810 }
811 
812 //________________________________________________________________
814 
815  if(isnorm) {return TMath::Exp(x/coeff)/(coeff*(TMath::Exp(fMassMax/coeff)-TMath::Exp(fMassMin/coeff)));}
816  else TMath::Exp(x/coeff);
817 }
818 
819 //________________________________________________________________
821 
822  switch(order) {
823  case 0:
824  if(isnorm) {return 1./(fMassMax-fMassMin);}
825  else {return pars[0];}
826  break;
827  case 1:
828  if(isnorm) {return (pars[0]-pars[1]/2*(fMassMax*fMassMax-fMassMin*fMassMin))/(fMassMax-fMassMin)+pars[1]*x;}
829  else {return pars[0]+pars[1]*x;}
830  break;
831  case 2:
832  if(isnorm) {return (pars[0]-pars[1]/2*(fMassMax*fMassMax-fMassMin*fMassMin)-pars[2]/3*(fMassMax*fMassMax*fMassMax-fMassMin*fMassMin*fMassMin))/(fMassMax-fMassMin)+pars[1]*x;}
833  else {return pars[0]+pars[1]*x+pars[2]*x*x;}
834  }
835  return 0;
836 }
837 
838 //________________________________________________________________
840 
841  Double_t mpi = TDatabasePDG::Instance()->GetParticle(211)->Mass();
842  return pars[0]*(pars[1]+1.)/(TMath::Power(fMassMax-mpi,pars[1]+1.)-TMath::Power(fMassMin-mpi,pars[1]+1.))*TMath::Power(x-mpi,pars[1]);
843 }
844 
845 //________________________________________________________________
847 
848  Double_t mpi = TDatabasePDG::Instance()->GetParticle(211)->Mass();
849  return pars[0]*TMath::Sqrt(x - mpi)*TMath::Exp(-1.*pars[1]*(x-mpi));
850 }
851 
852 //________________________________________________________________
854 
855  Double_t total=pars[0];
856  for(Int_t iT=1; iT<=Ndeg; iT++){
857  if(isnorm) total+=pars[iT]*TMath::Power(x-fMassParticle,iT)/TMath::Factorial(iT);
858  else total+=pars[iT]*TMath::Power(x,iT);
859  }
860  return total;
861 }
862 
863 //________________________________________________________________
865 
866  switch(fMassSgnFuncType) {
867  case 0:
868  return pars[0]*GetGausPDF(m[0],pars[1],pars[2]);
869  break;
870  case 1:
871  return pars[0]*(pars[3]*GetGausPDF(m[0],pars[1],pars[2])+(1-pars[3])*GetGausPDF(m[0],pars[1],pars[4]));
872  break;
873  }
874 
875  return 0;
876 }
877 
878 //________________________________________________________________
880 
881  switch(fMassBkgFuncType) {
882  case 0: //exponential
883  return pars[0]*GetExpoPDF(m[0],pars[1],kTRUE);
884  break;
885  case 1: //linear
886  return GetPolPDF(m[0],pars,1,kTRUE);
887  break;
888  case 2: //parabolic
889  return GetPolPDF(m[0],pars,2,kTRUE);
890  break;
891  case 3: //constant
892  return GetPolPDF(m[0],pars,0,kTRUE);
893  break;
894  case 4: //power law
895  return GetPowerFuncPDF(m[0],pars);
896  break;
897  case 5: //power law expo
898  return GetPowerExpoPDF(m[0],pars);
899  break;
900  case 6: //higher order (>=3) polinomial
901  return GetHigherPolFuncPDF(m[0],pars,fPolDegreeBkg,kTRUE);
902  break;
903  }
904  return 0;
905 }
906 
907 //_________________________________________________________________________
911  if(!fHistoTemplRfl) return 0;
912 
913  Int_t bin =fHistoTemplRfl->FindBin(m[0]);
914  Double_t value=fHistoTemplRfl->GetBinContent(bin);
915  Int_t binmin=fHistoTemplRfl->FindBin(fMassMin*1.00001);
916  Int_t binmax=fHistoTemplRfl->FindBin(fMassMax*0.99999);
917  Double_t norm=fHistoTemplRfl->Integral(binmin,binmax)*fHistoTemplRfl->GetBinWidth(bin);
918  if(TMath::Abs(value)<1.e-14 && fSmoothRfl){// very rough, assume a constant trend, much better would be a pol1 or pol2 over a broader range
919  value+=fHistoTemplRfl->GetBinContent(bin-1)+fHistoTemplRfl->GetBinContent(bin+1);
920  value/=3.;
921  }
922 
923  return pars[0]*value/norm*fRawYieldHelp*fMassHisto->GetBinWidth(1);
924 }
925 
926 //_________________________________________________________________________
928 
929  if(!fHistoTemplRfl) {return MassBkg(m,pars);}
930  else {
931  //bkg mass parameters
932  const Int_t nBkgPars = fNParsMassBkg;
933  Double_t bkgpars[nBkgPars];
934  for(Int_t iPar=0; iPar<fNParsMassBkg; iPar++) {bkgpars[iPar] = pars[iPar];}
935  //reflection parameters
936  Double_t rflpars[1]; //maximum number of parameters for rfl = 1 for the implemented functions
937  for(Int_t iPar=0; iPar<fNParsRfl; iPar++) {rflpars[iPar] = pars[iPar+fNParsMassBkg];}
938  return MassBkg(m,bkgpars)+MassRfl(m,rflpars);
939  }
940 }
941 
942 //_________________________________________________________________________
946 
947  return pars[0]*GetGausPDF(m[0],pars[1],pars[2]);
948 }
949 
950 //________________________________________________________________
952 
953  switch(fVnBkgFuncType) {
954  case 0: //expo
955  return pars[0]*GetExpoPDF(m[0],pars[1],kFALSE);
956  break;
957  case 1: //linear
958  return GetPolPDF(m[0],pars,1,kFALSE);
959  break;
960  case 2: //parabolic
961  return GetPolPDF(m[0],pars,2,kFALSE);
962  break;
963  case 6: //higher order (>=3) polinomial
964  return GetHigherPolFuncPDF(m[0],pars,fPolDegreeVnBkg,kFALSE);
965  break;
966  }
967  return 0;
968 }
969 
970 //________________________________________________________________
972 
973  //bkg mass parameters
974  const Int_t nBkgPars = fNParsMassBkg;
975  Double_t bkgpars[nBkgPars];
976  for(Int_t iPar=0; iPar<fNParsMassBkg; iPar++) {bkgpars[iPar] = pars[iPar];}
977  //signal mass parameters
978  Double_t sgnpars[5]; //maximum number of parameters for sgn = 5 for the implemented functions
979  for(Int_t iPar=0; iPar<fNParsMassSgn; iPar++) {sgnpars[iPar] = pars[iPar+fNParsMassBkg];}
980  //second peak parameters
981  Double_t secpeakpars[3]; //maximum number of parameters for second peak = 3 for the implemented functions
982  for(Int_t iPar=0; iPar<fNParsSec; iPar++) {secpeakpars[iPar] = pars[iPar+fNParsMassBkg+fNParsMassSgn];}
983  //reflection parameters
984  Double_t rflpars[1]; //maximum number of parameters for rfl = 1 for the implemented functions
985  for(Int_t iPar=0; iPar<fNParsRfl; iPar++) {rflpars[iPar] = pars[iPar+fNParsMassBkg+fNParsMassSgn+fNParsSec];}
986 
987  Double_t total = MassSignal(m,sgnpars)+MassBkg(m,bkgpars);
988  if(fSecondPeak) {total += MassSecondPeak(m,secpeakpars);}
989  if(fReflections) {total += MassRfl(m,rflpars);}
990 
991  return total;
992 }
993 
994 //________________________________________________________________
996 
997  //bkg mass parameters
998  const Int_t nBkgPars = fNParsMassBkg;
999  Double_t massbkgpars[nBkgPars];
1000  for(Int_t iPar=0; iPar<fNParsMassBkg; iPar++) {massbkgpars[iPar] = pars[iPar];}
1001  //signal mass parameters
1002  Double_t masssgnpars[5]; //maximum number of parameters for mass sgn = 5 for the implemented functions
1003  for(Int_t iPar=0; iPar<fNParsMassSgn; iPar++) {masssgnpars[iPar] = pars[iPar+fNParsMassBkg];}
1004  //second peak parameters
1005  Double_t secpeakpars[3]; //maximum number of parameters for second peak = 3 for the implemented functions
1006  for(Int_t iPar=0; iPar<fNParsSec; iPar++) {secpeakpars[iPar] = pars[iPar+fNParsMassBkg+fNParsMassSgn];}
1007  //reflection parameters
1008  Double_t rflpars[1]; //maximum number of parameters for rfl = 1 for the implemented functions
1009  for(Int_t iPar=0; iPar<fNParsRfl; iPar++) {rflpars[iPar] = pars[iPar+fNParsMassBkg+fNParsMassSgn+fNParsSec];}
1010  //bkg vn parameters
1011  const Int_t nVnBkgPars = fNParsVnBkg;
1012  Double_t vnbkgpars[nVnBkgPars];
1013  for(Int_t iPar=0; iPar<fNParsVnBkg; iPar++) {vnbkgpars[iPar] = pars[iPar+fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl];}
1014  //signal vn parameter
1015  Double_t vnSgn = pars[fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl+fNParsVnBkg];
1016  //second peak vn parameter
1017  Double_t vnSecPeak = 0;
1018  if(fSecondPeak && fDoSecondPeakVn) {vnSecPeak = pars[fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl+fNParsVnBkg+fNParsVnSgn];}
1019  //refl vn parameter
1020  Double_t vnRefl = 0;
1021  if(fReflections) {
1022  switch(fVnRflOpt) {
1023  case 0:
1024  vnRefl = pars[fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl+fNParsVnBkg];
1025  break;
1026  case 1:
1027  vnRefl = -pars[fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl+fNParsVnBkg];
1028  break;
1029  case 2:
1030  vnRefl = 0; //not used
1031  break;
1032  case 3:
1033  vnRefl = pars[fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl+fNParsVnBkg+fNParsVnSgn+fNParsVnSecPeak];
1034  break;
1035  default:
1036  AliError("Error in setting reflection vn option: check fVnRflOpt");
1037  break;
1038  }
1039  }
1040 
1041  Double_t vnBkg = vnBkgFunc(m,vnbkgpars);
1042  Double_t Sgn = MassSignal(m,masssgnpars);
1043  Double_t Bkg = MassBkg(m,massbkgpars);
1044  Double_t SecPeak = 0;
1045  if(fSecondPeak) {
1046  if(fDoSecondPeakVn) SecPeak += MassSecondPeak(m,secpeakpars);
1047  else Bkg += MassSecondPeak(m,secpeakpars);
1048  }
1049  Double_t Refl=0;
1050  if(fReflections) {
1051  if(fVnRflOpt==kSameVnBkg) Bkg += MassRfl(m,rflpars);
1052  else Refl += MassRfl(m,rflpars);
1053  }
1054 
1055  return (vnSgn*Sgn+vnBkg*Bkg+vnSecPeak*SecPeak+vnRefl*Refl)/(Sgn+Bkg+SecPeak+Refl);
1056 }
1057 
1058 //______________________________________________________________________________
1062  TCanvas* c0=new TCanvas("c0","",600,800);
1063  DrawHere(c0);
1064 }
Int_t fVnRflOpt
internal variable for fit with reflections
Int_t fFrac2GausFixed
flag to fix second peak width in case of k2Gaus
Int_t fPolDegreeVnBkg
degree of polynomial expansion for back fit (option 6 for back)
Double_t MassSecondPeak(Double_t *m, Double_t *par)
Double_t GetExpoPDF(Double_t x, Double_t slope, Bool_t isnorm=kTRUE)
Int_t fNParsVnBkg
number of parameters in mass bkg fit function
Bool_t fFixSecWidth
flag to fix the position of the 2nd peak
double Double_t
Definition: External.C:58
void SetFixReflOverS(Double_t rovers)
Int_t fNParsVnSecPeak
number of parameters in vn sgn fit function (1)
Int_t fNParsRfl
flag use/not use reflections
Bool_t fReflections
degree of polynomial expansion for vn back fit (option 6 for back)
void IncludeSecondGausPeak(Double_t mass, Bool_t fixm, Double_t width, Bool_t fixw)
Double_t MassRfl(Double_t *m, Double_t *par)
Bool_t SimultaneusFit(Bool_t drawFit=kTRUE)
Double_t GetPowerFuncPDF(Double_t x, Double_t *pars)
Bool_t fMeanFixedFromMassFit
initialization for fraction of second gaussian in case of k2Gaus
Double_t fMean
uncertainty on mass peak width from simultaneus fit
TH1F * fHistoTemplRfl
switch for fix refl/signal
Double_t fMassParticle
flag to fix fraction of second gaussian in case of k2Gaus
Bool_t fFixSecMass
width of the 2nd peak
Int_t fNParsSec
fit function for second peak
Bool_t fSecondPeak
maximum vn of reflections
Double_t fSigma2GausInit
initialization for peak position
Double_t GetHigherPolFuncPDF(Double_t x, Double_t *pars, Int_t Ndeg, Bool_t isnorm=kTRUE)
Bool_t fSigmaFixedFromMassFit
flag to fix peak position from mass prefit
void Signal(Double_t nOfSigma, Double_t &signal, Double_t &errsignal) const
Int_t MassFitter(Bool_t draw=kTRUE)
Double_t fSecMass
number of parameters in second peak fit function
TF1 * fMassTotFunc
mass signal fit function (final, after simultaneus fit)
Double_t vnBkgFunc(Double_t *m, Double_t *pars)
Double_t fRflOverSig
fit parameters in reflection fit function
Double_t fMassMax
upper mass limit
AliHFInvMassFitter class for the fit of invariant mass distribution of charm hadrons.
void SetFixSecondGaussianSigma(Double_t sigma)
TCanvas * c
Definition: TestFitELoss.C:172
Int_t fSigma2GausFixed
flag to fix peak position
Int_t fMeanFixed
flag to fix peak width
Double_t fVn
lower mass limit
Int_t fMassSgnFuncType
vn vs. mass histogram to fit
void DefineNumberOfParameters()
private methods
void SetFixFrac2Gaus(Double_t frac)
Double_t GetSigma() const
void SetPolDegreeForBackgroundFit(Int_t deg)
Double_t fChiSquare
uncertainty raw yield from simultaneus fit
Bool_t fSigma2GausFixedFromMassFit
flag to fix peak width from mass prefit
void SetFixGaussianSigma(Double_t sigma)
Double_t fRawYieldUncertainty
raw yield from simultaneus fit
Double_t GetPowerExpoPDF(Double_t x, Double_t *pars)
Int_t fNSigmaForSB
simultaneus fit probability
Double_t * sigma
TF1 * fMassBkgFunc
mass fit function (1st step, from prefit)
Double_t fMaxRefl
minimum for refelction histo
Double_t vnFunc(Double_t *m, Double_t *pars)
Double_t fVnUncertainty
vn of the signal from fit
Double_t fProb
simultaneus fit number of degree of freedom
Double_t MassSignal(Double_t *m, Double_t *pars)
void SetFixGaussianMean(Double_t mean)
Bool_t fFrac2GausFixedFromMassFit
flag to fix second peak width from mass prefit in case of k2Gaus
int Int_t
Definition: External.C:63
TF1 * fMassSecPeakFunc
switch off/on second peak (for D+->KKpi in Ds)
Double_t GetMean() const
TF1 * fVnBkgFunc
vn bkg fit function (1st step from SB prefit)
Int_t fPolDegreeBkg
flag to fix fraction of second gaussian in case of k2Gaus
Double_t fSecWidth
position of the 2nd peak
void SetInitialSecondGaussianSigma(Double_t sigma)
Bool_t fVnRflLimited
option for reflection vn type
TF1 * fVnTotFunc
vn bkg fit function (final, after simultaneus fit)
Int_t fHarmonic
vn uncertainty of second peak from fit
Double_t fSigmaInit
number of sigma for sidebands region (vn bkg prefit)
Int_t fNParsMassBkg
number of parameters in mass signal fit function
Double_t MassBkgRfl(Double_t *m, Double_t *par)
Double_t fMeanUncertainty
mass peak position from simultaneus fit
Double_t GetPolPDF(Double_t x, Double_t *pars, Int_t order, Bool_t isnorm=kTRUE)
TH1F * SetTemplateReflections(const TH1 *h, TString opt, Double_t minRange, Double_t maxRange)
Double_t fSigma
uncertainty on vn of the signal from simultaneus fit
TF1 * fVnBkgFuncSb
mass fit function (final, after simultaneus fit)
Double_t fMeanInit
initialization for peak width
TList * fitter
Definition: DrawAnaELoss.C:26
TH1F * fMassHisto
data members
void DrawHere(TVirtualPad *c)
TH1F * fHistoTemplRflInit
histogram with reflection template
Double_t GetGausPDF(Double_t x, Double_t mean, Double_t sigma)
fit functions
Int_t fNDF
simultaneus fit chi square
static void ComputeSignificance(Double_t signal, Double_t errsignal, Double_t background, Double_t errbackground, Double_t &significance, Double_t &errsignificance)
Significance calculator.
Bool_t fDoSecondPeakVn
vn of second peak from fit
TF1 * fMassRflFunc
initial histogram with reflection template
Double_t MassFunc(Double_t *m, Double_t *pars)
void Background(Double_t nOfSigma, Double_t &background, Double_t &errbackground) const
AliHFInvMassFitter * fMassFitter
vn fit function (final, after simultaneus fit)
TF1 * fMassBkgRflFunc
fit function for reflections
void SetInitialReflOverS(Double_t rovers)
Int_t fNParsVnRfl
number of parameters in vn sec peak fit function (1 if included, 0 otherwise)
TF1 * fMassFuncFromPrefit
type of vn bkg fit function
TH1F * fVnVsMassHisto
mass histogram to fit
Int_t fNParsMassSgn
mass of selected particle
void SetInitialGaussianSigma(Double_t sigma)
void SetInitialGaussianMean(Double_t mean)
Double_t minMass
Int_t fNParsVnSgn
number of parameters in vn bkg fit function
Double_t fSigmaUncertainty
mass peak width from simultaneus fit
Double_t fFrac2GausInit
initialization for second peak width in case of k2Gaus
Double_t fRawYield
uncertainty on mass peak position from simultaneus fit
Double_t fVnRflMax
minimum vn of reflections
Int_t fSigmaFixed
number of parameters in vn refl fit function (1 if included, 0 otherwise)
Double_t fMinRefl
refelction option
Int_t fVnBkgFuncType
type of mass bkg fit function
void Significance(Double_t nOfSigma, Double_t &significance, Double_t &errsignificance) const
Int_t fMassBkgFuncType
type of mass signal fit function
Double_t fVnRflMin
flag to limit or not the vn of reflections
TF1 * fMassSgnFunc
mass bkg fit function (final, after simultaneus fit)
void SetInitialFrac2Gaus(Double_t frac)
Double_t maxMass
bool Bool_t
Definition: External.C:53
Double_t fRawYieldHelp
switch for smoothing of reflection template
Bool_t fFixRflOverSig
reflection/signal
TString fRflOpt
mass bkg fit function plus reflections (final, after simultaneus fit)
Double_t fMassMin
mass fitter for mass prefit
Double_t MassBkg(Double_t *m, Double_t *pars)
Double_t fVnSecPeakUncertainty
flag to introduce second peak vn in the vn vs. mass fit
Bool_t fSmoothRfl
maximum for refelction histo
Double_t GetRawYield() const
Double_t fVnSecPeak
flag to fix the width of the 2nd peak