AliPhysics  608b256 (608b256)
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  if(!vnprefit) {AliError("Impossible to perform the bkg vn prefit"); return kFALSE;}
228 
229  std::vector<Double_t> initpars;
230  for(Int_t iBkgPar=0; iBkgPar<fNParsMassBkg; iBkgPar++) {
231  initpars.push_back(fMassFuncFromPrefit->GetParameter(iBkgPar));
232  }
233  for(Int_t iSgnPar=0; iSgnPar<fNParsMassSgn; iSgnPar++) {
234  initpars.push_back(fMassFuncFromPrefit->GetParameter(iSgnPar+fNParsMassBkg));
235  }
236  for(Int_t iSecPeakPar=0; iSecPeakPar<fNParsSec; iSecPeakPar++) {
237  initpars.push_back(fMassFuncFromPrefit->GetParameter(iSecPeakPar+fNParsMassBkg+fNParsMassSgn));
238  }
239  for(Int_t iReflPar=0; iReflPar<fNParsRfl; iReflPar++) {
240  initpars.push_back(fMassFuncFromPrefit->GetParameter(iReflPar+fNParsMassBkg+fNParsMassSgn+fNParsSec));
241  }
242  for(Int_t iVnBkgPar=0; iVnBkgPar<fNParsVnBkg; iVnBkgPar++) {
243  if(vnprefit) {initpars.push_back(fVnBkgFuncSb->GetParameter(iVnBkgPar));}
244  else {initpars.push_back(0.05);}
245  }
246  initpars.push_back(0.10); //initial parameter for signal vn
247  if(fSecondPeak && fDoSecondPeakVn) {initpars.push_back(0.10);} //initial parameter for second peak vn
248 
249  fMassTotFunc = new TF1("fMassTotFunc",this,&AliHFVnVsMassFitter::MassFunc,fMassMin,fMassMax,nparsmass,"AliHFVnVsMassFitter","MassFunc");
250  fVnTotFunc = new TF1("fVnTotFunc",this,&AliHFVnVsMassFitter::vnFunc,fMassMin,fMassMax,nparsvn,"AliHFVnVsMassFitter","vnFunc");
251  SetParNames();
252 
253  ROOT::Math::WrappedMultiTF1 wfTotMass(*fMassTotFunc,1);
254  ROOT::Math::WrappedMultiTF1 wfTotVn(*fVnTotFunc,1);
255 
256  // set data options and ranges
257  ROOT::Fit::DataOptions opt;
258  ROOT::Fit::DataRange rangeMass; //same range for two functions
259 
260  rangeMass.SetRange(fMassMin,fMassMax);
261  ROOT::Fit::BinData dataMass(opt,rangeMass);
262  ROOT::Fit::FillData(dataMass, fMassHisto);
263  ROOT::Fit::BinData dataVn(opt,rangeMass);
264  ROOT::Fit::FillData(dataVn, fVnVsMassHisto);
265 
266  //define the 2 chi squares
267  ROOT::Fit::Chi2Function chi2Mass(dataMass, wfTotMass);
268  ROOT::Fit::Chi2Function chi2Vn(dataVn, wfTotVn);
269 
270  //define the global chi square
271  AliHFGlobalChi2 globalChi2(chi2Mass, chi2Vn);
272 
273  //define fitter
274  ROOT::Fit::Fitter fitter;
275  // create before the parameter settings in order to fix or set range on them
276  fitter.Config().SetParamsSettings(nparsvn,initpars.data()); //set initial parameters from prefits
277  if(fMeanFixed==2 || fMeanFixedFromMassFit) {fitter.Config().ParSettings(fNParsMassBkg+1).Fix();}
278  if(fSigmaFixed==2 || fSigmaFixedFromMassFit) {fitter.Config().ParSettings(fNParsMassBkg+2).Fix();}
279  if(fMassSgnFuncType==k2Gaus) {
280  if(fFrac2GausFixed==2 || fFrac2GausFixedFromMassFit) {fitter.Config().ParSettings(fNParsMassBkg+3).Fix();}
281  if(fSigma2GausFixed==2 || fSigma2GausFixedFromMassFit) {fitter.Config().ParSettings(fNParsMassBkg+4).Fix();}
282  }
283  if(fSecondPeak) {
284  if(fFixSecMass) {fitter.Config().ParSettings(fNParsMassBkg+fNParsMassSgn+1).Fix();}
285  if(fFixSecWidth) {fitter.Config().ParSettings(fNParsMassBkg+fNParsMassSgn+2).Fix();}
286  }
287  if(fReflections) {
288  if(fFixRflOverSig) fitter.Config().ParSettings(fNParsMassBkg+fNParsMassSgn+fNParsSec).Fix();
289  if(fVnRflLimited) fitter.Config().ParSettings(nparsmass+fNParsVnBkg+NvnParsSgn-1).SetLimits(fVnRflMin,fVnRflMax);
290  }
291 
292  fitter.Config().MinimizerOptions().SetPrintLevel(0);
293  fitter.Config().SetMinimizer("Minuit2","Migrad");
294  for(Int_t iPar=0; iPar<nparsvn; iPar++) {fitter.Config().ParSettings(iPar).SetName(fVnTotFunc->GetParName(iPar));}
295  // fit FCN function directly
296  // (specify optionally data size and flag to indicate that is a chi2 fit
297  Bool_t isFitOk = fitter.FitFCN(nparsvn,globalChi2,0,dataMass.Size()+dataVn.Size(),kFALSE);
298  if(!isFitOk) return kFALSE;
299 
300  ROOT::Fit::FitResult result = fitter.Result();
301  result.Print(std::cout);
302 
303  //set parameters in every function
304  fVnBkgFunc = new TF1("fVnBkgFunc",this,&AliHFVnVsMassFitter::vnBkgFunc,fMassMin,fMassMax,fNParsVnBkg,"AliHFVnVsMassFitter","vnBkgFunc");
305  fMassBkgFunc = new TF1("fMassBkgFunc",this,&AliHFVnVsMassFitter::MassBkg,fMassMin,fMassMax,fNParsMassBkg,"AliHFVnVsMassFitter","MassBkg");
306  fMassSgnFunc = new TF1("fMassSgnFunc",this,&AliHFVnVsMassFitter::MassSignal,fMassMin,fMassMax,fNParsMassSgn,"AliHFVnVsMassFitter","MassSignal");
307  if(fReflections) {fMassRflFunc = new TF1("fMassRflFunc",this,&AliHFVnVsMassFitter::MassRfl,fMassMin,fMassMax,fNParsRfl,"AliHFVnVsMassFitter","MassRfl");}
308  if(fReflections) {fMassBkgRflFunc = new TF1("fMassBkgRflFunc",this,&AliHFVnVsMassFitter::MassBkgRfl,fMassMin,fMassMax,fNParsMassBkg+fNParsRfl,"AliHFVnVsMassFitter","MassBkgRfl");}
309  if(fSecondPeak) {fMassSecPeakFunc = new TF1("fMassSecPeakFunc",this,&AliHFVnVsMassFitter::MassSecondPeak,fMassMin,fMassMax,fNParsSec,"AliHFVnVsMassFitter","MassSecondPeak");}
310  for(Int_t iPar=0; iPar<nparsvn; iPar++) {
311  fVnTotFunc->SetParameter(iPar,result.Parameter(iPar));
312  fVnTotFunc->SetParError(iPar,result.ParError(iPar));
313  if(iPar<nparsmass) {
314  fMassTotFunc->SetParameter(iPar,result.Parameter(iPar));
315  fMassTotFunc->SetParError(iPar,result.ParError(iPar));
316  }
317  if(iPar>=nparsmass && iPar<nparsvn-NvnParsSgn) {
318  fVnBkgFunc->SetParameter(iPar-nparsmass,result.Parameter(iPar));
319  fVnBkgFunc->SetParError(iPar-nparsmass,result.ParError(iPar));
320  }
321  if(iPar>=fNParsMassBkg && iPar<fNParsMassBkg+fNParsMassSgn) {
322  fMassSgnFunc->SetParameter(iPar-fNParsMassBkg,result.Parameter(iPar));
323  }
324  if(iPar<fNParsMassBkg) {
325  fMassBkgFunc->SetParameter(iPar,result.Parameter(iPar));
326  fMassBkgFunc->SetParError(iPar,result.ParError(iPar));
327  if(fReflections) {
328  fMassBkgRflFunc->SetParameter(iPar,result.Parameter(iPar));
329  fMassBkgRflFunc->SetParError(iPar,result.ParError(iPar));
330  }
331  }
332  if(fReflections && (iPar>=fNParsMassBkg+fNParsMassSgn+fNParsSec && iPar<fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl)) {
333  fMassRflFunc->SetParameter(iPar-(fNParsMassBkg+fNParsMassSgn+fNParsSec),result.Parameter(iPar));
334  fMassRflFunc->SetParError(iPar-(fNParsMassBkg+fNParsMassSgn+fNParsSec),result.ParError(iPar));
335  fMassBkgRflFunc->SetParameter(iPar-(fNParsMassSgn+fNParsSec),result.Parameter(iPar));
336  fMassBkgRflFunc->SetParError(iPar-(fNParsMassSgn+fNParsSec),result.ParError(iPar));
337  }
338  if(fSecondPeak && (iPar>=fNParsMassBkg+fNParsMassSgn && iPar<fNParsMassBkg+fNParsMassSgn+fNParsSec)) {
339  fMassSecPeakFunc->SetParameter(iPar-(fNParsMassBkg+fNParsMassSgn),result.Parameter(iPar));
340  fMassSecPeakFunc->SetParError(iPar-(fNParsMassBkg+fNParsMassSgn),result.ParError(iPar));
341  }
342  }
343  if(drawFit) {DrawFit();}
344 
345  fVn = fVnTotFunc->GetParameter(fVnTotFunc->GetNpar()-NvnParsSgn);
346  fVnUncertainty = fVnTotFunc->GetParError(fVnTotFunc->GetNpar()-NvnParsSgn);
347  if(fDoSecondPeakVn) {
348  fVnSecPeak = fVnTotFunc->GetParameter(fVnTotFunc->GetNpar()-1);
349  fVnSecPeakUncertainty = fVnTotFunc->GetParError(fVnTotFunc->GetNpar()-1);
350  }
351  fRawYield = fVnTotFunc->GetParameter(fNParsMassBkg)/fMassHisto->GetBinWidth(10);
352  fRawYieldUncertainty = fVnTotFunc->GetParError(fNParsMassBkg)/fMassHisto->GetBinWidth(10);
353  fMean = fVnTotFunc->GetParameter(fNParsMassBkg+1);
354  fMeanUncertainty = fVnTotFunc->GetParError(fNParsMassBkg+1);
355  fSigma = fVnTotFunc->GetParameter(fNParsMassBkg+2);
356  fSigmaUncertainty = fVnTotFunc->GetParError(fNParsMassBkg+2);
357  fChiSquare = result.MinFcnValue();
358  fNDF = result.Ndf();
359  fProb = result.Prob();
360 
361  return kTRUE;
362 }
363 
364 //______________________________________________________________________________
365 void AliHFVnVsMassFitter::DrawHere(TVirtualPad* c){
367 
368  gStyle->SetOptStat(0);
369  gStyle->SetCanvasColor(0);
370  gStyle->SetFrameFillColor(0);
371  c->Divide(1,2);
372 
373  c->cd(1);
374  fMassHisto->SetTitle("");
375  fMassHisto->SetMarkerStyle(20);
376  fMassHisto->SetMarkerSize(1);
377  fMassHisto->SetMarkerColor(kBlack);
378  fMassHisto->SetLineColor(kBlack);
379  fMassHisto->GetYaxis()->SetRangeUser(0.,fMassHisto->GetMaximum()*1.2);
380  fMassHisto->GetXaxis()->SetRangeUser(fMassMin,fMassMax);
381  fMassHisto->Draw("E");
382  if(fMassFuncFromPrefit) {
383  fMassFuncFromPrefit->SetLineColor(kGray+1);
384  fMassFuncFromPrefit->SetLineStyle(7);
386  fMassFuncFromPrefit->Draw("same");
387  }
388  if(fMassBkgFunc) {
389  fMassBkgFunc->SetLineColor(kRed);
390  fMassBkgFunc->SetRange(fMassMin,fMassMax);
391  fMassBkgFunc->Draw("same");
392  }
393  if(fMassRflFunc) {
394  fMassRflFunc->SetLineColor(kGreen+2);
395  fMassRflFunc->SetRange(fMassMin,fMassMax);
396  fMassRflFunc->Draw("same");
397  }
398  if(fMassBkgRflFunc) {
399  fMassBkgRflFunc->SetLineColor(kOrange+7);
400  fMassBkgRflFunc->SetLineStyle(9);
401  fMassBkgRflFunc->SetRange(fMassMin,fMassMax);
402  fMassBkgRflFunc->Draw("same");
403  }
404  if(fMassSecPeakFunc) {
405  fMassSecPeakFunc->SetLineColor(kMagenta+1);
406  fMassSecPeakFunc->SetLineStyle(7);
408  fMassSecPeakFunc->Draw("same");
409  }
410  if(fMassTotFunc) {
411  fMassTotFunc->SetLineColor(kBlue);
412  fMassTotFunc->SetRange(fMassMin,fMassMax);
413  fMassTotFunc->Draw("same");
414  }
415  TPaveText* massinfo = new TPaveText(0.45,0.7,1.,0.87,"NDC");
416  massinfo->SetTextFont(42);
417  massinfo->SetTextSize(0.05);
418  massinfo->SetBorderSize(0);
419  massinfo->SetFillStyle(0);
420  massinfo->SetTextColor(kBlue);
421  massinfo->AddText(Form("mean = %.3f #pm %.3f",fVnTotFunc->GetParameter(fNParsMassBkg+1),fVnTotFunc->GetParError(fNParsMassBkg+1)));
422  massinfo->AddText(Form("sigma = %.3f #pm %.3f",fVnTotFunc->GetParameter(fNParsMassBkg+2),fVnTotFunc->GetParError(fNParsMassBkg+2)));
423  if(fMassSgnFuncType==k2Gaus) {
424  massinfo->AddText(Form("sigma2 = %.3f #pm %.3f",fVnTotFunc->GetParameter(fNParsMassBkg+3),fVnTotFunc->GetParError(fNParsMassBkg+3)));
425  }
426  massinfo->Draw("same");
427 
428  c->cd(2);
429  fVnVsMassHisto->SetTitle("");
430  fVnVsMassHisto->SetMarkerStyle(20);
431  fVnVsMassHisto->SetMarkerSize(1);
432  fVnVsMassHisto->SetMarkerColor(kBlack);
433  fVnVsMassHisto->SetLineColor(kBlack);
434  fVnVsMassHisto->GetYaxis()->SetRangeUser(fVnVsMassHisto->GetMinimum()-0.15,fVnVsMassHisto->GetMaximum()+0.20);
435  fVnVsMassHisto->GetXaxis()->SetRangeUser(fMassMin,fMassMax);
436  fVnBkgFuncSb->SetRange(fMassMin,fMassMax);
437  fVnVsMassHisto->Draw("E");
438  if(fVnBkgFuncSb) {
439  fVnBkgFuncSb->SetLineColor(kGray+1);
440  fVnBkgFuncSb->SetLineStyle(7);
441  fVnBkgFuncSb->Draw("same");
442  }
443  if(fVnBkgFunc) {
444  fVnBkgFunc->SetLineColor(kRed);
445  fVnBkgFunc->SetRange(fMassMin,fMassMax);
446  fVnBkgFunc->Draw("same");
447  }
448  if(fVnTotFunc) {
449  fVnTotFunc->SetLineColor(kBlue);
450  fVnTotFunc->SetRange(fMassMin,fMassMax);
451  fVnTotFunc->Draw("same");
452  }
453 
454  TPaveText* vninfo = new TPaveText(-0.45,0.7,1.,0.87,"NDC");
455  vninfo->SetTextFont(42);
456  vninfo->SetTextSize(0.05);
457  vninfo->SetBorderSize(0);
458  vninfo->SetFillStyle(0);
459  Int_t NvnParsSgn = 1;
460  if(fSecondPeak && fDoSecondPeakVn) {NvnParsSgn+=1;}
461  if(fReflections && fVnRflOpt==kFreePar) {NvnParsSgn+=1;}
462  vninfo->AddText(Form("#it{v}_{%d}^{sgn} = %.3f #pm %.3f",fHarmonic,fVnTotFunc->GetParameter(fVnTotFunc->GetNpar()-NvnParsSgn),fVnTotFunc->GetParError(fVnTotFunc->GetNpar()-NvnParsSgn)));
463  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)));}
464  vninfo->AddText(Form("#chi^{2}/#it{ndf} = %.2f/%d",fChiSquare,fNDF));
465  vninfo->Draw("same");
466 
467  c->Update();
468 }
469 
470 //________________________________________________________________
472 
473  //define proper maxs and mins from histos
474  Double_t tmpmin = TMath::Max(fMassHisto->GetBinLowEdge(1),fVnVsMassHisto->GetBinLowEdge(1));
475  fMassMin=TMath::Max(fMassMin,tmpmin);
476  Double_t tmpmax = TMath::Min(fMassHisto->GetBinLowEdge(fMassHisto->GetNbinsX())+fMassHisto->GetBinWidth(fMassHisto->GetNbinsX()),fVnVsMassHisto->GetBinLowEdge(fVnVsMassHisto->GetNbinsX())+fVnVsMassHisto->GetBinWidth(fVnVsMassHisto->GetNbinsX()));
477  fMassMax=TMath::Min(fMassMax,tmpmax);
478 
484  if(fMassSgnFuncType==k2Gaus) {
489  }
493  if(fReflections) {
497  }
498  Bool_t status = fMassFitter->MassFitter(kFALSE);
499 
500  if(status) {
502  fMassFuncFromPrefit->SetName("fMassFuncFromPrefit");
503  }
505 
506  return status;
507 }
508 
509 //________________________________________________________________
511 
512  Double_t mean = fMassFitter->GetMean();
514  const Int_t nMassBins = fVnVsMassHisto->GetNbinsX();
515  Double_t SBbins[nMassBins];
516  Int_t nSBbins=0;
517  for(Int_t iBin=0; iBin<nMassBins; iBin++) {
518  Double_t min = fVnVsMassHisto->GetBinLowEdge(iBin+1);
519  Double_t max = fVnVsMassHisto->GetBinLowEdge(iBin+1)+fVnVsMassHisto->GetBinWidth(iBin+1);
520  if(max<mean){
521  if(max<(mean-fNSigmaForSB*sigma)) {SBbins[iBin]=1; nSBbins++;}
522  else {SBbins[iBin]=0;}
523  }
524  if(min>=mean){
525  if(min>(mean+fNSigmaForSB*sigma)) {SBbins[iBin]=1; nSBbins++;}
526  else {SBbins[iBin]=0;}
527  }
528  }
529  TGraphErrors* gVnVsMassSB = new TGraphErrors(nSBbins);
530  for(Int_t iBin=0; iBin<nMassBins; iBin++) {
531  if(SBbins[iBin]==1) {
532  gVnVsMassSB->SetPoint(iBin,fVnVsMassHisto->GetBinCenter(iBin+1),fVnVsMassHisto->GetBinContent(iBin+1));
533  gVnVsMassSB->SetPointError(iBin,fVnVsMassHisto->GetBinWidth(iBin+1)/2,fVnVsMassHisto->GetBinError(iBin+1));
534  }
535  }
536  fVnBkgFuncSb = new TF1("fVnBkgFuncSb",this,&AliHFVnVsMassFitter::vnBkgFunc,fMassMin,fMassMax,fNParsVnBkg,"AliHFVnVsMassFitter","vnBkgFunc");
537  switch(fVnBkgFuncType) {
538  case 1:
539  fVnBkgFuncSb->SetParName(0,"ConstVnBkg");
540  fVnBkgFuncSb->SetParName(1,"SlopeVnBkg");
541  break;
542  case 2:
543  fVnBkgFuncSb->SetParName(0,"ConstVnBkg");
544  fVnBkgFuncSb->SetParName(1,"Coef1VnBkg");
545  fVnBkgFuncSb->SetParName(2,"Coef2VnBkg");
546  break;
547  default:
548  AliError("Error in setting signal par names: check fVnBkgFuncType");
549  break;
550  }
551  Int_t status = gVnVsMassSB->Fit(fVnBkgFuncSb,"","",fMassMin,fMassMax);
552  delete gVnVsMassSB;
553 
554  if(status==0) return kTRUE;
555  return kFALSE;
556 }
557 
558 //________________________________________________________________
560 
561  switch(fMassSgnFuncType) {
562  case 0: //single gaus
563  fNParsMassSgn=3;
564  break;
565  case 1: //double gaus
566  fNParsMassSgn=5;
567  break;
568  default:
569  AliError("Error in computing fMassSgnFuncType: check fMassSgnFuncType");
570  break;
571  }
572 
573  switch(fMassBkgFuncType) {
574  case 0: //expo
575  fNParsMassBkg=2;
576  break;
577  case 1: //lin
578  fNParsMassBkg=2;
579  break;
580  case 2: //pol2
581  fNParsMassBkg=3;
582  break;
583  case 3: //no bkg
584  fNParsMassBkg=1;
585  break;
586  case 4: //power law
587  fNParsMassBkg=2;
588  break;
589  case 5: //power expo
590  fNParsMassBkg=2;
591  break;
592  case 6: //high degree pol
594  break;
595  default:
596  AliError("Error in computing fNParsMassBkg: check fMassBkgFuncType");
597  break;
598  }
599 
600  switch(fVnBkgFuncType) {
601  case 0: //expo
602  fNParsVnBkg=2;
603  break;
604  case 1: //lin
605  fNParsVnBkg=2;
606  break;
607  case 2: //pol2
608  fNParsVnBkg=3;
609  break;
610  case 6: //high degree pol
612  break;
613  default:
614  AliError("Error in computing fNParsVnBkg: check fVnBkgFuncType");
615  break;
616  }
617 
618  fNParsVnSgn=1;
619 
620  if(fReflections) {
621  fNParsRfl=1;
622  if(fVnRflOpt==3) fNParsVnRfl=1;
623  else fNParsVnRfl=0;
624  }
625  else {
626  fNParsRfl=0;
627  fNParsVnRfl=0;
628  }
629 
630  if(fSecondPeak) {
631  fNParsSec=3;
632  fNParsVnSecPeak=1;
633  }
634  else {
635  fNParsSec=0;
636  fNParsVnSecPeak=0;
637  }
638 }
639 
640 //________________________________________________________________
642 
643  switch(fMassSgnFuncType) {
644  case 0: //single gaus
645  fVnTotFunc->SetParName(fNParsMassBkg,"SgnInt");
646  fVnTotFunc->SetParName(fNParsMassBkg+1,"Mean");
647  fVnTotFunc->SetParName(fNParsMassBkg+2,"Sigma");
648  break;
649  case 1: //double gaus
650  fVnTotFunc->SetParName(fNParsMassBkg,"SgnInt");
651  fVnTotFunc->SetParName(fNParsMassBkg+1,"Mean");
652  fVnTotFunc->SetParName(fNParsMassBkg+2,"Sigma1");
653  fVnTotFunc->SetParName(fNParsMassBkg+3,"Frac");
654  fVnTotFunc->SetParName(fNParsMassBkg+4,"Sigma2");
655  break;
656  default:
657  AliError("Error in setting signal par names: check fMassSgnFuncType");
658  break;
659  }
660 
661  switch(fMassBkgFuncType) {
662  case 0: //expo
663  fVnTotFunc->SetParName(0,"BkgInt");
664  fVnTotFunc->SetParName(1,"Slope");
665  break;
666  case 1: //lin
667  fVnTotFunc->SetParName(0,"BkgInt");
668  fVnTotFunc->SetParName(1,"Slope");
669  break;
670  case 2: //pol2
671  fVnTotFunc->SetParName(0,"BkgInt");
672  fVnTotFunc->SetParName(1,"Coef1");
673  fVnTotFunc->SetParName(2,"Coef2");
674  break;
675  case 3: //no bkg
676  fVnTotFunc->SetParName(0,"Const");
677  break;
678  case 4: //power law
679  fVnTotFunc->SetParName(0,"BkgInt");
680  fVnTotFunc->SetParName(1,"Coef1");
681  break;
682  case 5: //power expo
683  fVnTotFunc->SetParName(0,"Coef1");
684  fVnTotFunc->SetParName(1,"Coef2");
685  break;
686  case 6: //high degree pol
687  fVnTotFunc->SetParName(0,"BkgInt");
688  for(Int_t iPar=1; iPar<fNParsMassBkg; iPar++) {fVnTotFunc->SetParName(iPar,Form("Coef%d",iPar));}
689  break;
690  default:
691  AliError("Error in setting signal par names: check fMassBkgFuncType");
692  break;
693  }
694 
695  for(Int_t iPar=0; iPar<fNParsVnBkg; iPar++) {fVnTotFunc->SetParName(fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl+iPar,fVnBkgFuncSb->GetParName(iPar));}
696 
697  if(fReflections) {fVnTotFunc->SetParName(fNParsMassBkg+fNParsMassSgn+fNParsSec,"ReflOverS");}
698 
699  if(fSecondPeak) {
700  fVnTotFunc->SetParName(fNParsMassBkg+fNParsMassSgn,"SecPeakInt");
701  fVnTotFunc->SetParName(fNParsMassBkg+fNParsMassSgn+1,"SecPeakMean");
702  fVnTotFunc->SetParName(fNParsMassBkg+fNParsMassSgn+2,"SecPeakSigma");
703  }
704  fVnTotFunc->SetParName(fNParsMassBkg+fNParsMassSgn+fNParsRfl+fNParsSec+fNParsVnBkg,Form("v%dSgn",fHarmonic));
705 
706  if(fSecondPeak && fDoSecondPeakVn) {fVnTotFunc->SetParName(fNParsMassBkg+fNParsMassSgn+fNParsRfl+fNParsSec+fNParsVnBkg+1,Form("v%dSecPeak",fHarmonic));}
708 }
709 
710 //_________________________________________________________________________
711 void AliHFVnVsMassFitter::Signal(Double_t nOfSigma,Double_t &signal,Double_t &errsignal) const {
714 
715  Double_t minMass=fMean-nOfSigma*fSigma;
716  Double_t maxMass=fMean+nOfSigma*fSigma;
717  Signal(minMass,maxMass,signal,errsignal);
718  return;
719 }
720 
721 //_________________________________________________________________________
722 void AliHFVnVsMassFitter::Signal(Double_t min, Double_t max, Double_t &signal,Double_t &errsignal) const {
725  if(!fMassSgnFunc) {signal=-1; errsignal=0; return;}
726 
727  signal=fMassSgnFunc->Integral(min, max)/(Double_t)fMassHisto->GetBinWidth(1);
728  errsignal=(fRawYieldUncertainty/fRawYield)*signal;/*assume relative error is the same as for total integral*/
729 
730  return;
731 }
732 
733 //___________________________________________________________________________
734 void AliHFVnVsMassFitter::Background(Double_t nOfSigma,Double_t &background,Double_t &errbackground) const {
737 
738  Double_t minMass=fMean-nOfSigma*fSigma;
739  Double_t maxMass=fMean+nOfSigma*fSigma;
740  Background(minMass,maxMass,background,errbackground);
741 
742  return;
743 }
744 
745 //___________________________________________________________________________
746 void AliHFVnVsMassFitter::Background(Double_t min, Double_t max, Double_t &background,Double_t &errbackground) const {
749 
750  if(!fMassBkgFunc) {background=-1; errbackground=0; return;}
751 
752  Double_t intB=fMassBkgFunc->GetParameter(0);
753  Double_t intBerr=fMassBkgFunc->GetParError(0);
754  //relative error evaluation: from histo
755 
756  Int_t leftBand=fMassHisto->FindBin(fMean-4*fSigma);
757  Int_t rightBand=fMassHisto->FindBin(fMean+4*fSigma);
758  intB=fMassHisto->Integral(1,leftBand)+fMassHisto->Integral(rightBand,fMassHisto->GetNbinsX());
759  Double_t sum2=0;
760  for(Int_t iBin=1; iBin<=leftBand; iBin++){
761  sum2+=fMassHisto->GetBinError(iBin)*fMassHisto->GetBinError(iBin);
762  }
763  for(Int_t iBin=rightBand; iBin<=fMassHisto->GetNbinsX(); iBin++){
764  sum2+=fMassHisto->GetBinError(iBin)*fMassHisto->GetBinError(iBin);
765  }
766 
767  intBerr=TMath::Sqrt(sum2);
768 
769  background=fMassBkgFunc->Integral(min,max)/(Double_t)fMassHisto->GetBinWidth(1);
770  errbackground=intBerr/intB*background;
771 
772  return;
773 }
774 
775 //__________________________________________________________________________
776 
777 void AliHFVnVsMassFitter::Significance(Double_t nOfSigma,Double_t &significance,Double_t &errsignificance) const {
780 
781  Double_t minMass=fMean-nOfSigma*fSigma;
782  Double_t maxMass=fMean+nOfSigma*fSigma;
783  Significance(minMass, maxMass, significance, errsignificance);
784 
785  return;
786 }
787 
788 //__________________________________________________________________________
789 
790 void AliHFVnVsMassFitter::Significance(Double_t min, Double_t max, Double_t &significance,Double_t &errsignificance) const {
793 
794  Double_t background,errbackground;
795  Background(min,max,background,errbackground);
796 
797  if (fRawYield+background <= 0.){
798  significance=-1;
799  errsignificance=0;
800  return;
801  }
802 
803  AliVertexingHFUtils::ComputeSignificance(fRawYield,fRawYieldUncertainty,background,errbackground,significance,errsignificance);
804 
805  return;
806 }
807 
808 //________________________________________________________________
810 
811  return TMath::Gaus(x,mean,sigma,kTRUE);
812 }
813 
814 //________________________________________________________________
816 
817  if(isnorm) {return TMath::Exp(x/coeff)/(coeff*(TMath::Exp(fMassMax/coeff)-TMath::Exp(fMassMin/coeff)));}
818  else return TMath::Exp(x/coeff);
819 }
820 
821 //________________________________________________________________
823 
824  switch(order) {
825  case 0:
826  if(isnorm) {return 1./(fMassMax-fMassMin);}
827  else {return pars[0];}
828  break;
829  case 1:
830  if(isnorm) {return (pars[0]-pars[1]/2*(fMassMax*fMassMax-fMassMin*fMassMin))/(fMassMax-fMassMin)+pars[1]*x;}
831  else {return pars[0]+pars[1]*x;}
832  break;
833  case 2:
834  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;}
835  else {return pars[0]+pars[1]*x+pars[2]*x*x;}
836  }
837  return 0;
838 }
839 
840 //________________________________________________________________
842 
843  Double_t mpi = TDatabasePDG::Instance()->GetParticle(211)->Mass();
844  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]);
845 }
846 
847 //________________________________________________________________
849 
850  Double_t mpi = TDatabasePDG::Instance()->GetParticle(211)->Mass();
851  return pars[0]*TMath::Sqrt(x - mpi)*TMath::Exp(-1.*pars[1]*(x-mpi));
852 }
853 
854 //________________________________________________________________
856 
857  Double_t total=pars[0];
858  for(Int_t iT=1; iT<=Ndeg; iT++){
859  if(isnorm) total+=pars[iT]*TMath::Power(x-fMassParticle,iT)/TMath::Factorial(iT);
860  else total+=pars[iT]*TMath::Power(x,iT);
861  }
862  return total;
863 }
864 
865 //________________________________________________________________
867 
868  switch(fMassSgnFuncType) {
869  case 0:
870  return pars[0]*GetGausPDF(m[0],pars[1],pars[2]);
871  break;
872  case 1:
873  return pars[0]*(pars[3]*GetGausPDF(m[0],pars[1],pars[2])+(1-pars[3])*GetGausPDF(m[0],pars[1],pars[4]));
874  break;
875  }
876 
877  return 0;
878 }
879 
880 //________________________________________________________________
882 
883  switch(fMassBkgFuncType) {
884  case 0: //exponential
885  return pars[0]*GetExpoPDF(m[0],pars[1],kTRUE);
886  break;
887  case 1: //linear
888  return GetPolPDF(m[0],pars,1,kTRUE);
889  break;
890  case 2: //parabolic
891  return GetPolPDF(m[0],pars,2,kTRUE);
892  break;
893  case 3: //constant
894  return GetPolPDF(m[0],pars,0,kTRUE);
895  break;
896  case 4: //power law
897  return GetPowerFuncPDF(m[0],pars);
898  break;
899  case 5: //power law expo
900  return GetPowerExpoPDF(m[0],pars);
901  break;
902  case 6: //higher order (>=3) polinomial
903  return GetHigherPolFuncPDF(m[0],pars,fPolDegreeBkg,kTRUE);
904  break;
905  }
906  return 0;
907 }
908 
909 //_________________________________________________________________________
913  if(!fHistoTemplRfl) return 0;
914 
915  Int_t bin =fHistoTemplRfl->FindBin(m[0]);
916  Double_t value=fHistoTemplRfl->GetBinContent(bin);
917  Int_t binmin=fHistoTemplRfl->FindBin(fMassMin*1.00001);
918  Int_t binmax=fHistoTemplRfl->FindBin(fMassMax*0.99999);
919  Double_t norm=fHistoTemplRfl->Integral(binmin,binmax)*fHistoTemplRfl->GetBinWidth(bin);
920  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
921  value+=fHistoTemplRfl->GetBinContent(bin-1)+fHistoTemplRfl->GetBinContent(bin+1);
922  value/=3.;
923  }
924 
925  return pars[0]*value/norm*fRawYieldHelp*fMassHisto->GetBinWidth(1);
926 }
927 
928 //_________________________________________________________________________
930 
931  if(!fHistoTemplRfl) {return MassBkg(m,pars);}
932  else {
933  //bkg mass parameters
934  const Int_t nBkgPars = fNParsMassBkg;
935  Double_t bkgpars[nBkgPars];
936  for(Int_t iPar=0; iPar<fNParsMassBkg; iPar++) {bkgpars[iPar] = pars[iPar];}
937  //reflection parameters
938  Double_t rflpars[1]; //maximum number of parameters for rfl = 1 for the implemented functions
939  for(Int_t iPar=0; iPar<fNParsRfl; iPar++) {rflpars[iPar] = pars[iPar+fNParsMassBkg];}
940  return MassBkg(m,bkgpars)+MassRfl(m,rflpars);
941  }
942 }
943 
944 //_________________________________________________________________________
948 
949  return pars[0]*GetGausPDF(m[0],pars[1],pars[2]);
950 }
951 
952 //________________________________________________________________
954 
955  switch(fVnBkgFuncType) {
956  case 0: //expo
957  return pars[0]*GetExpoPDF(m[0],pars[1],kFALSE);
958  break;
959  case 1: //linear
960  return GetPolPDF(m[0],pars,1,kFALSE);
961  break;
962  case 2: //parabolic
963  return GetPolPDF(m[0],pars,2,kFALSE);
964  break;
965  case 6: //higher order (>=3) polinomial
966  return GetHigherPolFuncPDF(m[0],pars,fPolDegreeVnBkg,kFALSE);
967  break;
968  }
969  return 0;
970 }
971 
972 //________________________________________________________________
974 
975  //bkg mass parameters
976  const Int_t nBkgPars = fNParsMassBkg;
977  Double_t bkgpars[nBkgPars];
978  for(Int_t iPar=0; iPar<fNParsMassBkg; iPar++) {bkgpars[iPar] = pars[iPar];}
979  //signal mass parameters
980  Double_t sgnpars[5]; //maximum number of parameters for sgn = 5 for the implemented functions
981  for(Int_t iPar=0; iPar<fNParsMassSgn; iPar++) {sgnpars[iPar] = pars[iPar+fNParsMassBkg];}
982  //second peak parameters
983  Double_t secpeakpars[3]; //maximum number of parameters for second peak = 3 for the implemented functions
984  for(Int_t iPar=0; iPar<fNParsSec; iPar++) {secpeakpars[iPar] = pars[iPar+fNParsMassBkg+fNParsMassSgn];}
985  //reflection parameters
986  Double_t rflpars[1]; //maximum number of parameters for rfl = 1 for the implemented functions
987  for(Int_t iPar=0; iPar<fNParsRfl; iPar++) {rflpars[iPar] = pars[iPar+fNParsMassBkg+fNParsMassSgn+fNParsSec];}
988 
989  Double_t total = MassSignal(m,sgnpars)+MassBkg(m,bkgpars);
990  if(fSecondPeak) {total += MassSecondPeak(m,secpeakpars);}
991  if(fReflections) {total += MassRfl(m,rflpars);}
992 
993  return total;
994 }
995 
996 //________________________________________________________________
998 
999  //bkg mass parameters
1000  const Int_t nBkgPars = fNParsMassBkg;
1001  Double_t massbkgpars[nBkgPars];
1002  for(Int_t iPar=0; iPar<fNParsMassBkg; iPar++) {massbkgpars[iPar] = pars[iPar];}
1003  //signal mass parameters
1004  Double_t masssgnpars[5]; //maximum number of parameters for mass sgn = 5 for the implemented functions
1005  for(Int_t iPar=0; iPar<fNParsMassSgn; iPar++) {masssgnpars[iPar] = pars[iPar+fNParsMassBkg];}
1006  //second peak parameters
1007  Double_t secpeakpars[3]; //maximum number of parameters for second peak = 3 for the implemented functions
1008  for(Int_t iPar=0; iPar<fNParsSec; iPar++) {secpeakpars[iPar] = pars[iPar+fNParsMassBkg+fNParsMassSgn];}
1009  //reflection parameters
1010  Double_t rflpars[1]; //maximum number of parameters for rfl = 1 for the implemented functions
1011  for(Int_t iPar=0; iPar<fNParsRfl; iPar++) {rflpars[iPar] = pars[iPar+fNParsMassBkg+fNParsMassSgn+fNParsSec];}
1012  //bkg vn parameters
1013  const Int_t nVnBkgPars = fNParsVnBkg;
1014  Double_t vnbkgpars[nVnBkgPars];
1015  for(Int_t iPar=0; iPar<fNParsVnBkg; iPar++) {vnbkgpars[iPar] = pars[iPar+fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl];}
1016  //signal vn parameter
1017  Double_t vnSgn = pars[fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl+fNParsVnBkg];
1018  //second peak vn parameter
1019  Double_t vnSecPeak = 0;
1020  if(fSecondPeak && fDoSecondPeakVn) {vnSecPeak = pars[fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl+fNParsVnBkg+fNParsVnSgn];}
1021  //refl vn parameter
1022  Double_t vnRefl = 0;
1023  if(fReflections) {
1024  switch(fVnRflOpt) {
1025  case 0:
1026  vnRefl = pars[fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl+fNParsVnBkg];
1027  break;
1028  case 1:
1029  vnRefl = -pars[fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl+fNParsVnBkg];
1030  break;
1031  case 2:
1032  vnRefl = 0; //not used
1033  break;
1034  case 3:
1035  vnRefl = pars[fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl+fNParsVnBkg+fNParsVnSgn+fNParsVnSecPeak];
1036  break;
1037  default:
1038  AliError("Error in setting reflection vn option: check fVnRflOpt");
1039  break;
1040  }
1041  }
1042 
1043  Double_t vnBkg = vnBkgFunc(m,vnbkgpars);
1044  Double_t Sgn = MassSignal(m,masssgnpars);
1045  Double_t Bkg = MassBkg(m,massbkgpars);
1046  Double_t SecPeak = 0;
1047  if(fSecondPeak) {
1048  if(fDoSecondPeakVn) SecPeak += MassSecondPeak(m,secpeakpars);
1049  else Bkg += MassSecondPeak(m,secpeakpars);
1050  }
1051  Double_t Refl=0;
1052  if(fReflections) {
1053  if(fVnRflOpt==kSameVnBkg) Bkg += MassRfl(m,rflpars);
1054  else Refl += MassRfl(m,rflpars);
1055  }
1056 
1057  return (vnSgn*Sgn+vnBkg*Bkg+vnSecPeak*SecPeak+vnRefl*Refl)/(Sgn+Bkg+SecPeak+Refl);
1058 }
1059 
1060 //______________________________________________________________________________
1064  TCanvas* c0=new TCanvas("c0","",600,800);
1065  DrawHere(c0);
1066 }
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