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