AliPhysics  fde8a9f (fde8a9f)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliHFInvMassFitter.cxx
Go to the documentation of this file.
1 #include <TH1F.h>
2 #include <TF1.h>
3 #include <TMath.h>
4 #include <TVirtualFitter.h>
5 #include <TDatabasePDG.h>
6 #include <TCanvas.h>
7 #include <TStyle.h>
8 #include <TPaveText.h>
9 
10 #include "AliLog.h"
11 #include "AliVertexingHFUtils.h"
12 #include "AliHFInvMassFitter.h"
13 
17 
18 //__________________________________________________________________________
20  TNamed(),
21  fHistoInvMass(0x0),
22  fMinMass(0),
23  fMaxMass(5),
24  fTypeOfFit4Bkg(kExpo),
25  fTypeOfFit4Sgn(kGaus),
26  fMass(1.865),
27  fMassErr(0.),
28  fSigmaSgn(0.012),
29  fSigmaSgnErr(0.),
30  fFixedMean(kFALSE),
31  fFixedSigma(kFALSE),
32  fFixedRawYield(-1.),
33  fNSigPars(3),
34  fNBkgPars(2),
35  fOnlySideBands(kFALSE),
36  fFitOption("L,E"),
37  fRawYield(0.),
38  fRawYieldErr(0.),
39  fSigFunc(0x0),
40  fBkgFuncSb(0x0),
41  fBkgFunc(0x0),
42  fBkgFuncRef(0x0),
43  fSecondPeak(kFALSE),
44  fSecMass(-999.),
45  fSecWidth(9999.),
46  fFixSecMass(kFALSE),
47  fFixSecWidth(kFALSE),
48  fSecFunc(0x0),
49  fFuncTot(0x0)
50 {
52 }
53 
54 //__________________________________________________________________________
55 AliHFInvMassFitter::AliHFInvMassFitter(const TH1F *histoToFit, Double_t minvalue, Double_t maxvalue, Int_t fittypeb, Int_t fittypes):
56  TNamed(),
57  fHistoInvMass(0x0),
58  fMinMass(minvalue),
59  fMaxMass(maxvalue),
60  fTypeOfFit4Bkg(fittypeb),
61  fTypeOfFit4Sgn(fittypes),
62  fMass(1.865),
63  fMassErr(0.),
64  fSigmaSgn(0.012),
65  fSigmaSgnErr(0.),
66  fFixedMean(kFALSE),
67  fFixedSigma(kFALSE),
68  fFixedRawYield(-1.),
69  fNSigPars(3),
70  fNBkgPars(2),
71  fOnlySideBands(kFALSE),
72  fFitOption("L,E"),
73  fRawYield(0.),
74  fRawYieldErr(0.),
75  fSigFunc(0x0),
76  fBkgFuncSb(0x0),
77  fBkgFunc(0x0),
78  fBkgFuncRef(0x0),
79  fSecondPeak(kFALSE),
80  fSecMass(-999.),
81  fSecWidth(9999.),
82  fFixSecMass(kFALSE),
83  fFixSecWidth(kFALSE),
84  fSecFunc(0x0),
85  fFuncTot(0x0)
86 {
88  fHistoInvMass=(TH1F*)histoToFit->Clone("fHistoInvMass");
89  fHistoInvMass->SetDirectory(0);
91 }
92 //_________________________________________________________________________
94 
96 
97  delete fHistoInvMass;
98  delete fSigFunc;
99  delete fBkgFunc;
100  delete fBkgFuncRef;
101  delete fSecFunc;
102  delete fFuncTot;
103 
104 }
105 //__________________________________________________________________________
107  switch (fTypeOfFit4Bkg) {
108  case 0:
109  fNBkgPars=2;
110  break;
111  case 1:
112  fNBkgPars=2;
113  break;
114  case 2:
115  fNBkgPars=3;
116  break;
117  case 3:
118  fNBkgPars=1;
119  break;
120  case 4:
121  fNBkgPars=2;
122  break;
123  case 5:
124  fNBkgPars=3;
125  break;
126  default:
127  AliError("Error in computing fNBkgPars: check fTypeOfFit4Bkg");
128  break;
129  }
130  switch (fTypeOfFit4Sgn) {
131  case 0:
132  fNSigPars=3;
133  break;
134  case 1:
135  fNSigPars=5;
136  break;
137  default:
138  AliError("Error in computing fNSigPars: check fTypeOfFit4Sgn");
139  break;
140  }
141 
142 }
143 //__________________________________________________________________________
145  TVirtualFitter::SetDefaultFitter("Minuit");
146 
147  Double_t integralHisto=fHistoInvMass->Integral(fHistoInvMass->FindBin(fMinMass),fHistoInvMass->FindBin(fMaxMass),"width");
148 
149  fOnlySideBands = kTRUE;
150  fBkgFuncSb = CreateBackgroundFitFunction("funcbkgsb",integralHisto);
151  Int_t status=fHistoInvMass->Fit("funcbkgsb",Form("R,%s,+,0",fFitOption.Data()));
152  fBkgFuncSb->SetLineColor(kGray+1);
153  if (status != 0){
154  printf("Failed first fit with only background, minuit status = %d\n",status);
155  return kFALSE;
156  }
157 
158  fOnlySideBands = kFALSE;
159  if(!fBkgFunc){
160  fBkgFunc = CreateBackgroundFitFunction("funcbkg",integralHisto);
161  for(Int_t ipar=0; ipar<fNBkgPars; ipar++) fBkgFunc->SetParameter(ipar,fBkgFuncSb->GetParameter(ipar));
162  }
163  fBkgFunc->SetLineColor(kGray+1);
164  Double_t estimSignal=CheckForSignal(fMass,fSigmaSgn);
165  if(estimSignal<0.){
166  if(draw) DrawFit();
167  return kTRUE;
168  }
169  if(!fBkgFuncRef){
170  fBkgFuncRef = CreateBackgroundFitFunction("funcbkgrefit",integralHisto);
171  for(Int_t ipar=0; ipar<fNBkgPars; ipar++) fBkgFuncRef->SetParameter(ipar,fBkgFunc->GetParameter(ipar));
172  }
173  fBkgFuncRef->SetLineColor(2);
174  fSigFunc = CreateSignalFitFunction("fsigfit",estimSignal);
175  if(fSecondPeak){
177  fSecFunc = CreateSecondPeakFunction("fsecpeak",estimSec);
178  }
179  fFuncTot = CreateTotalFitFunction("funcmass");
180  status=fHistoInvMass->Fit("funcmass",Form("R,%s,+,0",fFitOption.Data()));
181  if (status != 0){
182  printf("Failed fit with signal+background, minuit status = %d\n",status);
183  return kFALSE;
184  }
185  for(Int_t ipar=0; ipar<fNBkgPars; ipar++) fBkgFuncRef->SetParameter(ipar,fFuncTot->GetParameter(ipar));
186  for(Int_t ipar=0; ipar<fNSigPars; ipar++) fSigFunc->SetParameter(ipar,fFuncTot->GetParameter(ipar+fNBkgPars));
187  fMass=fSigFunc->GetParameter(1);
188  fMassErr=fSigFunc->GetParError(1);
189  fSigmaSgn=fSigFunc->GetParameter(2);
190  fSigmaSgnErr=fSigFunc->GetParError(2);
191  fFuncTot->SetLineColor(4);
192  fRawYield=fFuncTot->GetParameter(fNBkgPars)/fHistoInvMass->GetBinWidth(1);
193  fRawYieldErr=fFuncTot->GetParError(fNBkgPars)/fHistoInvMass->GetBinWidth(1);
194  if(draw) DrawFit();
195  return kTRUE;
196 }
197 //______________________________________________________________________________
199 
200  TCanvas* c0=new TCanvas("c0");
201  DrawHere(c0);
202 }
203 //______________________________________________________________________________
204 void AliHFInvMassFitter::DrawHere(TVirtualPad* c){
205  gStyle->SetOptStat(0);
206  gStyle->SetCanvasColor(0);
207  gStyle->SetFrameFillColor(0);
208  c->cd();
209  fHistoInvMass->GetXaxis()->SetRangeUser(fMinMass,fMaxMass);
210  fHistoInvMass->SetMarkerStyle(20);
211  fHistoInvMass->SetMinimum(0.);
212  fHistoInvMass->Draw("PE");
213  if(fBkgFunc) fBkgFunc->Draw("same");
214  if(fBkgFuncRef) fBkgFuncRef->Draw("same");
215  if(fFuncTot) fFuncTot->Draw("same");
216  TPaveText *pinfos=new TPaveText(0.12,0.65,0.47,0.89,"NDC");
217  TPaveText *pinfom=new TPaveText(0.6,0.7,1.,.87,"NDC");
218  pinfos->SetBorderSize(0);
219  pinfos->SetFillStyle(0);
220  pinfom->SetBorderSize(0);
221  pinfom->SetFillStyle(0);
222  if(fFuncTot){
223  pinfom->SetTextColor(kBlue);
224  for(Int_t ipar=1; ipar<fNSigPars; ipar++){
225  pinfom->AddText(Form("%s = %.3f #pm %.3f",fFuncTot->GetParName(ipar+fNBkgPars),fFuncTot->GetParameter(ipar+fNBkgPars),fFuncTot->GetParError(ipar+fNBkgPars)));
226  }
227  pinfom->Draw();
228 
229  Double_t nsigma=3;
230  Double_t bkg,errbkg;
231  Background(nsigma,bkg,errbkg);
232  Double_t signif,errsignif;
233  Significance(nsigma,signif,errsignif);
234 
235  pinfos->AddText(Form("S = %.0f #pm %.0f ",fRawYield,fRawYieldErr));
236  pinfos->AddText(Form("B (%.0f#sigma) = %.0f #pm %.0f",nsigma,bkg,errbkg));
237  pinfos->AddText(Form("S/B = (%.0f#sigma) %.4f ",nsigma,fRawYield/bkg));
238  pinfos->AddText(Form("Signif (%.0f#sigma) = %.1f #pm %.1f ",nsigma,signif,errsignif));
239  pinfos->Draw();
240  }
241  c->Update();
242  return;
243 }
244 //______________________________________________________________________________
246  Double_t minForSig=mean-4.*sigma;
247  Double_t maxForSig=mean+4.*sigma;
248  Int_t binForMinSig=fHistoInvMass->FindBin(minForSig);
249  Int_t binForMaxSig=fHistoInvMass->FindBin(maxForSig);
250  Double_t sum=0.;
251  Double_t sumback=0.;
252  fBkgFunc->Print();
253  for(Int_t ibin=binForMinSig; ibin<=binForMaxSig; ibin++){
254  sum+=fHistoInvMass->GetBinContent(ibin);
255  sumback+=fBkgFunc->Eval(fHistoInvMass->GetBinCenter(ibin));
256  }
257  Double_t diffUnderPeak=(sum-sumback);
258  printf("intUnderFunc=%f intUnderHisto=%f sestimSig=%f\n",sum,sumback,diffUnderPeak);
259  if(diffUnderPeak/TMath::Sqrt(sum)<1.){
260  printf("(Tot-Bkg)/sqrt(Tot)=%f ---> Likely no signal/\n",diffUnderPeak/TMath::Sqrt(sum));
261  return -1;
262  }
263  return diffUnderPeak*fHistoInvMass->GetBinWidth(1);
264 }
265 
266 //______________________________________________________________________________
269  TF1* funcbkg = new TF1(fname.Data(),this,&AliHFInvMassFitter::FitFunction4Bkg,fMinMass,fMaxMass,fNBkgPars,"AliHFInvMassFitter","FitFunction4Bkg");
270  switch (fTypeOfFit4Bkg) {
271  case 0: //gaus+expo
272  funcbkg->SetParNames("BkgInt","Slope");
273  funcbkg->SetParameters(integral,-2.);
274  break;
275  case 1:
276  funcbkg->SetParNames("BkgInt","Slope");
277  funcbkg->SetParameters(integral,-100.);
278  break;
279  case 2:
280  funcbkg->SetParNames("BkgInt","Coef1","Coef2");
281  funcbkg->SetParameters(integral,-10.,5);
282  break;
283  case 3:
284  funcbkg->SetParNames("Const");
285  funcbkg->SetParameter(0,0.);
286  funcbkg->FixParameter(0,0.);
287  break;
288  case 4:
289  funcbkg->SetParNames("BkgInt","Coef1");
290  funcbkg->SetParameters(integral,0.5);
291  break;
292  case 5:
293  funcbkg->SetParNames("BkgInt","Coef1","Coef2");
294  funcbkg->SetParameters(integral,-10.,5.);
295  break;
296  default:
297  AliError(Form("Wrong choice of fTypeOfFit4Bkg (%d)",fTypeOfFit4Bkg));
298  delete funcbkg;
299  return 0x0;
300  break;
301  }
302  funcbkg->SetLineColor(kBlue+3);
303  return funcbkg;
304 }
305 //______________________________________________________________________________
307  TF1* funcsec = new TF1(fname.Data(),this,&AliHFInvMassFitter::FitFunction4SecPeak,fMinMass,fMaxMass,3,"AliHFInvMassFitter","FitFunction4SecPeak");
308  funcsec->SetParameter(0,integsig);
309  funcsec->SetParameter(1,fSecMass);
310  if(fFixSecMass) funcsec->FixParameter(1,fSecMass);
311  funcsec->SetParameter(2,fSecWidth);
312  if(fFixSecWidth) funcsec->FixParameter(2,fSecWidth);
313  funcsec->SetParNames("SecPeakInt","SecMean","SecSigma");
314  return funcsec;
315 }
316 
317 //______________________________________________________________________________
320  TF1* funcsig = new TF1(fname.Data(),this,&AliHFInvMassFitter::FitFunction4Sgn,fMinMass,fMaxMass,fNSigPars,"AliHFInvMassFitter","FitFunction4Sgn");
321  if(fTypeOfFit4Sgn==kGaus){
322  funcsig->SetParameter(0,integsig);
323  if(fFixedRawYield>-0.1) funcsig->FixParameter(0,fFixedRawYield);
324  funcsig->SetParameter(1,fMass);
325  if(fFixedMean) funcsig->FixParameter(1,fMass);
326  funcsig->SetParameter(2,fSigmaSgn);
327  if(fFixedSigma) funcsig->FixParameter(2,fSigmaSgn);
328  funcsig->SetParNames("SgnInt","Mean","Sigma");
329  }
330  if(fTypeOfFit4Sgn==k2Gaus){
331  funcsig->SetParameter(0,integsig);
332  if(fFixedRawYield>-0.1) funcsig->FixParameter(0,fFixedRawYield);
333  funcsig->SetParameter(1,fMass);
334  if(fFixedMean) funcsig->FixParameter(1,fMass);
335  funcsig->SetParameter(2,fSigmaSgn);
336  funcsig->SetParLimits(2,0.004,0.05);
337  if(fFixedSigma) funcsig->FixParameter(2,fSigmaSgn);
338  funcsig->SetParameter(3,0.2);
339  funcsig->SetParLimits(3,0.,1.);
340  funcsig->SetParameter(4,fSigmaSgn);
341  funcsig->SetParLimits(4,0.004,0.05);
342  funcsig->SetParNames("SgnInt","Mean","Sigma1","Frac","Sigma2");
343  }
344  return funcsig;
345 }
346 
347 //______________________________________________________________________________
350  Int_t nParSecPeak=0;
351  if(fSecondPeak && fSecFunc) nParSecPeak=3;
352  TF1* ftot=ftot=new TF1(fname.Data(),this,&AliHFInvMassFitter::FitFunction4Mass,fMinMass,fMaxMass,(fNBkgPars+nParSecPeak+fNSigPars),"AliHFInvMassFitter","FitFunction4Mass");
353  for(Int_t ipar=0; ipar<fNBkgPars; ipar++){
354  ftot->SetParameter(ipar,fBkgFunc->GetParameter(ipar));
355  ftot->SetParName(ipar,fBkgFunc->GetParName(ipar));
356  }
357  printf("Copy parameters for signal = %d\n",fNSigPars);
358  for(Int_t ipar=0; ipar<fNSigPars; ipar++){
359  ftot->SetParameter(ipar+fNBkgPars,fSigFunc->GetParameter(ipar));
360  ftot->SetParName(ipar+fNBkgPars,fSigFunc->GetParName(ipar));
361  Double_t parmin,parmax;
362  fSigFunc->GetParLimits(ipar,parmin,parmax);
363  ftot->SetParLimits(ipar+fNBkgPars,parmin,parmax);
364  printf("Par %d val %f vs %f\n",ipar,fSigFunc->GetParameter(ipar),ftot->GetParameter(ipar+fNBkgPars));
365  }
366  if(fSecondPeak && fSecFunc){
367  for(Int_t ipar=0; ipar<nParSecPeak; ipar++){
368  ftot->SetParameter(ipar+fNBkgPars+fNSigPars,fSecFunc->GetParameter(ipar));
369  ftot->SetParName(ipar+fNBkgPars+fNSigPars,fSecFunc->GetParName(ipar));
370  Double_t parmin,parmax;
371  fSecFunc->GetParLimits(ipar,parmin,parmax);
372  ftot->SetParLimits(ipar+fNBkgPars+fNSigPars,parmin,parmax);
373  }
374  }
375  return ftot;
376 }
377 //__________________________________________________________________________
379  // Fit function for the background
380 
381  Double_t maxDeltaM = 4.*fSigmaSgn;
382  if(fOnlySideBands && TMath::Abs(x[0]-fMass) < maxDeltaM) {
383  TF1::RejectPoint();
384  return 0;
385  }
386  if(fOnlySideBands && fSecondPeak && TMath::Abs(x[0]-fSecMass) < (4.*fSecWidth)){
387  TF1::RejectPoint();
388  return 0;
389  }
390  Double_t total=0;
391 
392  switch (fTypeOfFit4Bkg){
393  case 0:
394  //exponential
395  //exponential = A*exp(B*x) -> integral(exponential)=A/B*exp(B*x)](min,max)
396  //-> A = B*integral/(exp(B*max)-exp(B*min)) where integral can be written
397  //as integralTot- integralGaus (=par [2])
398  //Par:
399  // * [0] = integralBkg;
400  // * [1] = B;
401  //exponential = [1]*[0]/(exp([1]*max)-exp([1]*min))*exp([1]*x)
402  total = par[0]*par[1]/(TMath::Exp(par[1]*fMaxMass)-TMath::Exp(par[1]*fMinMass))*TMath::Exp(par[1]*x[0]);
403  // AliInfo("Background function set to: exponential");
404  break;
405  case 1:
406  //linear
407  //y=a+b*x -> integral = a(max-min)+1/2*b*(max^2-min^2) -> a = (integral-1/2*b*(max^2-min^2))/(max-min)=integral/(max-min)-1/2*b*(max+min)
408  // * [0] = integralBkg;
409  // * [1] = b;
410  total= par[0]/(fMaxMass-fMinMass)+par[1]*(x[0]-0.5*(fMaxMass+fMinMass));
411  // AliInfo("Background function set to: linear");
412  break;
413  case 2:
414  //parabola
415  //y=a+b*x+c*x**2 -> integral = a(max-min) + 1/2*b*(max^2-min^2) +
416  //+ 1/3*c*(max^3-min^3) ->
417  //a = (integral-1/2*b*(max^2-min^2)-1/3*c*(max^3-min^3))/(max-min)
418  // * [0] = integralBkg;
419  // * [1] = b;
420  // * [2] = c;
421  total = par[0]/(fMaxMass-fMinMass)+par[1]*(x[0]-0.5*(fMaxMass+fMinMass))+par[2]*(x[0]*x[0]-1/3.*(fMaxMass*fMaxMass*fMaxMass-fMinMass*fMinMass*fMinMass)/(fMaxMass-fMinMass));
422  // AliInfo("Background function set to: polynomial");
423  break;
424  case 3:
425  total=par[0];
426  break;
427  case 4:
428  //power function
429  //y=a(x-m_pi)^b -> integral = a/(b+1)*((max-m_pi)^(b+1)-(min-m_pi)^(b+1))
430  //
431  //a = integral*(b+1)/((max-m_pi)^(b+1)-(min-m_pi)^(b+1))
432  // * [0] = integralBkg;
433  // * [1] = b;
434  // a(power function) = [0]*([1]+1)/((max-m_pi)^([1]+1)-(min-m_pi)^([1]+1))*(x-m_pi)^[1]
435  {
436  Double_t mpi = TDatabasePDG::Instance()->GetParticle(211)->Mass();
437 
438  total = par[0]*(par[1]+1.)/(TMath::Power(fMaxMass-mpi,par[1]+1.)-TMath::Power(fMinMass-mpi,par[1]+1.))*TMath::Power(x[0]-mpi,par[1]);
439  // AliInfo("Background function set to: powerlaw");
440  }
441  break;
442  case 5:
443  //power function wit exponential
444  //y=a*Sqrt(x-m_pi)*exp(-b*(x-m_pi))
445  {
446  Double_t mpi = TDatabasePDG::Instance()->GetParticle(211)->Mass();
447 
448  total = par[1]*TMath::Sqrt(x[0] - mpi)*TMath::Exp(-1.*par[2]*(x[0]-mpi));
449  // AliInfo("Background function set to: wit exponential");
450  }
451  break;
452  }
453  return total;
454 }
455 //_________________________________________________________________________
457  // Fit function for the signal
458 
459 
460  // AliInfo("Signal function set to: Gaussian");
461  Double_t sigval=0;
462  switch (fTypeOfFit4Sgn){
463  case 0:
464  //gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
465  //Par:
466  // * [0] = integralSgn
467  // * [1] = mean
468  // * [2] = sigma
469  //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
470  sigval=par[0]/TMath::Sqrt(2.*TMath::Pi())/par[2]*TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/2./par[2]/par[2]);
471  break;
472  case 1:
473  //double gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
474  //Par:
475  // * [0] = integralSgn
476  // * [1] = mean
477  // * [2] = sigma1
478  // * [3] = 2nd gaussian ratio
479  // * [4] = deltaSigma
480  //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
481  Double_t g1=(1.-par[3])/TMath::Sqrt(2.*TMath::Pi())/par[2]*TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/2./par[2]/par[2]);
482  Double_t g2=par[3]/TMath::Sqrt(2.*TMath::Pi())/par[4]*TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/2./par[4]/par[4]);
483  sigval=par[0]*(g1+g2);
484  break;
485  }
486  return sigval;
487 }
488 //_________________________________________________________________________
490  // Fit function for a second gaussian peak (for D+->KKpi in teh Ds mass spectrum)
491 
492  //gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
493  //Par:
494  // * [0] = integralSgn
495  // * [1] = mean
496  // * [2] = sigma
497  Double_t secgaval=par[0]/TMath::Sqrt(2.*TMath::Pi())/par[2]*TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/2./par[2]/par[2]);
498  return secgaval;
499 }
500 //_________________________________________________________________________
502  // Total
503 
504  Double_t bkg=FitFunction4Bkg(x,par);
505  Double_t sig=FitFunction4Sgn(x,&par[fNBkgPars]);
506  Double_t sec=0.;
507  if(fSecondPeak) sec=FitFunction4SecPeak(x,&par[fNBkgPars+fNSigPars]);
508  return bkg+sig+sec;
509 }
510 
511 //_________________________________________________________________________
512 void AliHFInvMassFitter::Signal(Double_t nOfSigma,Double_t &signal,Double_t &errsignal) const {
513  // Return signal integral in mean+- n sigma
514 
515  Double_t min=fMass-nOfSigma*fSigmaSgn;
516  Double_t max=fMass+nOfSigma*fSigmaSgn;
517  Signal(min,max,signal,errsignal);
518  return;
519 }
520 
521 //_________________________________________________________________________
522 void AliHFInvMassFitter::Signal(Double_t min, Double_t max, Double_t &signal,Double_t &errsignal) const {
523 
524  // Return signal integral in a range
525  signal=fSigFunc->Integral(min, max)/(Double_t)fHistoInvMass->GetBinWidth(1);
526  errsignal=(fRawYieldErr/fRawYield)*signal;/*assume relative error is the same as for total integral*/
527  return;
528 }
529 
530 //___________________________________________________________________________
531 void AliHFInvMassFitter::Background(Double_t nOfSigma,Double_t &background,Double_t &errbackground) const {
532  // Return background integral in mean+- n sigma
533 
534  Double_t min=fMass-nOfSigma*fSigmaSgn;
535  Double_t max=fMass+nOfSigma*fSigmaSgn;
536  Background(min,max,background,errbackground);
537  return;
538 
539 }
540 //___________________________________________________________________________
541 void AliHFInvMassFitter::Background(Double_t min, Double_t max, Double_t &background,Double_t &errbackground) const {
542  // Return background integral in a range
543 
544  TF1 *funcbkg=0x0;
545  if(fBkgFuncRef) funcbkg=fBkgFuncRef;
546  else if(fBkgFunc) funcbkg=fBkgFunc;
547  if(!funcbkg){
548  AliError("Bkg function not found!");
549  return;
550  }
551 
552  Double_t intB=funcbkg->GetParameter(0);
553  Double_t intBerr=funcbkg->GetParError(0);
554 
555  //relative error evaluation: from histo
556 
557  Int_t leftBand=fHistoInvMass->FindBin(fMass-4*fSigmaSgn);
558  Int_t rightBand=fHistoInvMass->FindBin(fMass+4*fSigmaSgn);
559  intB=fHistoInvMass->Integral(1,leftBand)+fHistoInvMass->Integral(rightBand,fHistoInvMass->GetNbinsX());
560  Double_t sum2=0;
561  for(Int_t i=1;i<=leftBand;i++){
562  sum2+=fHistoInvMass->GetBinError(i)*fHistoInvMass->GetBinError(i);
563  }
564  for(Int_t i=rightBand; i<=fHistoInvMass->GetNbinsX();i++){
565  sum2+=fHistoInvMass->GetBinError(i)*fHistoInvMass->GetBinError(i);
566  }
567 
568  intBerr=TMath::Sqrt(sum2);
569 
570  background=funcbkg->Integral(min,max)/(Double_t)fHistoInvMass->GetBinWidth(1);
571  errbackground=intBerr/intB*background;
572 
573  return;
574 
575 }
576 //__________________________________________________________________________
577 
578 void AliHFInvMassFitter::Significance(Double_t nOfSigma,Double_t &significance,Double_t &errsignificance) const {
579  // Return significance in mean+- n sigma
580 
581  Double_t min=fMass-nOfSigma*fSigmaSgn;
582  Double_t max=fMass+nOfSigma*fSigmaSgn;
583  Significance(min, max, significance, errsignificance);
584 
585  return;
586 }
587 
588 //__________________________________________________________________________
589 
590 void AliHFInvMassFitter::Significance(Double_t min, Double_t max, Double_t &significance,Double_t &errsignificance) const {
591  // Return significance integral in a range
592 
593  Double_t background,errbackground;
594  Background(min,max,background,errbackground);
595 
596  if (fRawYield+background <= 0.){
597  significance=-1;
598  errsignificance=0;
599  return;
600  }
601 
602  AliVertexingHFUtils::ComputeSignificance(fRawYield,fRawYieldErr,background,errbackground,significance,errsignificance);
603 
604  return;
605 }
Int_t fTypeOfFit4Sgn
background fit func
Double_t fSigmaSgnErr
signal gaussian sigma
Bool_t fFixSecWidth
flag to fix the position of the 2nd peak
double Double_t
Definition: External.C:58
Double_t fRawYield
L, LW or Chi2.
Bool_t fOnlySideBands
fit parameters in background fit function
virtual void Signal(Double_t nOfSigma, Double_t &signal, Double_t &errsignal) const
Double_t FitFunction4SecPeak(Double_t *x, Double_t *par)
Double_t fMass
signal fit func
TF1 * CreateTotalFitFunction(TString fname)
Double_t fRawYieldErr
signal gaussian integral
TCanvas * c
Definition: TestFitELoss.C:172
TF1 * CreateSecondPeakFunction(TString fname, Double_t integral)
Double_t fFixedRawYield
switch for fix Sigma of gaussian
Double_t FitFunction4Bkg(Double_t *x, Double_t *par)
Double_t fSigmaSgn
unc on signal gaussian mean value
Double_t fSecMass
swicth off/on second peak (for D+->KKpi in Ds)
Bool_t fFixSecMass
width of the 2nd peak
TF1 * fSecFunc
flag to fix the width of the 2nd peak
Double_t * sigma
TF1 * fBkgFuncRef
background fit function (1st step)
int Int_t
Definition: External.C:63
Double_t fSecWidth
position of the 2nd peak
Int_t fNSigPars
initialization for wa yield
TString fFitOption
kTRUE = only side bands considered
TF1 * CreateSignalFitFunction(TString fname, Double_t integral)
Bool_t fFixedMean
unc on signal gaussian sigma
Double_t nsigma
Int_t fTypeOfFit4Bkg
upper mass limit
TF1 * fBkgFunc
background fit function (1st step)
Double_t fMassErr
signal gaussian mean value
Double_t fMinMass
histogram to fit
TF1 * fFuncTot
fit function for second peak
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 FitFunction4Mass(Double_t *x, Double_t *par)
void Significance(Double_t nOfSigma, Double_t &significance, Double_t &errsignificance) const
Bool_t MassFitter(Bool_t draw=kTRUE)
Bool_t draw[nPtBins]
Double_t fMaxMass
lower mass limit
TF1 * CreateBackgroundFitFunction(TString fname, Double_t integral)
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
void Background(Double_t nOfSigma, Double_t &background, Double_t &errbackground) const
Double_t FitFunction4Sgn(Double_t *x, Double_t *par)
bool Bool_t
Definition: External.C:53
Double_t CheckForSignal(Double_t mean, Double_t sigma)
TF1 * fSigFunc
err on signal gaussian integral
TF1 * fBkgFuncSb
Signal fit function.
Bool_t fSecondPeak
background fit function (2nd step)
void DrawHere(TVirtualPad *c)
Int_t fNBkgPars
fit parameters in signal fit function
Bool_t fFixedSigma
switch for fix mean of gaussian