AliPhysics  9b6b435 (9b6b435)
AliHFInvMassFitter.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2019, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 #include <TH1F.h>
17 #include <TF1.h>
18 #include <TMath.h>
19 #include <TVirtualFitter.h>
20 #include <TDatabasePDG.h>
21 #include <TCanvas.h>
22 #include <TStyle.h>
23 #include <TPaveText.h>
24 
25 #include "AliVertexingHFUtils.h"
26 #include "AliHFInvMassFitter.h"
27 
29 ClassImp(AliHFInvMassFitter);
31 
47 
48 //__________________________________________________________________________
50  TNamed(),
51  fHistoInvMass(0x0),
52  fMinMass(0),
53  fMaxMass(5),
54  fTypeOfFit4Bkg(kExpo),
55  fPolDegreeBkg(4),
56  fCurPolDegreeBkg(-1),
57  fMassParticle(1.864),
58  fTypeOfFit4Sgn(kGaus),
59  fMass(1.865),
60  fMassErr(0.),
61  fSigmaSgn(0.012),
62  fSigmaSgnErr(0.),
63  fSigmaSgn2Gaus(0.012),
64  fFixedMean(kFALSE),
65  fFixedSigma(kFALSE),
66  fFixedSigma2Gaus(kFALSE),
67  fFixedRawYield(-1.),
68  fFrac2Gaus(0.2),
69  fFixedFrac2Gaus(kFALSE),
70  fRatio2GausSigma(0.),
71  fFixedRatio2GausSigma(kFALSE),
72  fNParsSig(3),
73  fNParsBkg(2),
74  fOnlySideBands(kFALSE),
75  fNSigma4SideBands(4.),
76  fFitOption("L,E"),
77  fRawYield(0.),
78  fRawYieldErr(0.),
79  fSigFunc(0x0),
80  fBkgFuncSb(0x0),
81  fBkgFunc(0x0),
82  fBkgFuncRefit(0x0),
83  fReflections(kFALSE),
84  fNParsRfl(0),
85  fRflOverSig(0),
86  fFixRflOverSig(kFALSE),
87  fHistoTemplRfl(0x0),
88  fSmoothRfl(kFALSE),
89  fRawYieldHelp(0),
90  fRflFunc(0x0),
91  fBkRFunc(0x0),
92  fSecondPeak(kFALSE),
93  fNParsSec(0),
94  fSecMass(-999.),
95  fSecWidth(9999.),
96  fFixSecMass(kFALSE),
97  fFixSecWidth(kFALSE),
98  fSecFunc(0x0),
99  fTotFunc(0x0)
100 {
102 }
103 
104 //__________________________________________________________________________
105 AliHFInvMassFitter::AliHFInvMassFitter(const TH1F *histoToFit, Double_t minvalue, Double_t maxvalue, Int_t fittypeb, Int_t fittypes):
106  TNamed(),
107  fHistoInvMass(0x0),
108  fMinMass(minvalue),
109  fMaxMass(maxvalue),
110  fTypeOfFit4Bkg(fittypeb),
111  fPolDegreeBkg(4),
112  fCurPolDegreeBkg(-1),
113  fMassParticle(1.864),
114  fTypeOfFit4Sgn(fittypes),
115  fMass(1.865),
116  fMassErr(0.),
117  fSigmaSgn(0.012),
118  fSigmaSgnErr(0.),
119  fSigmaSgn2Gaus(0.012),
120  fFixedMean(kFALSE),
121  fFixedSigma(kFALSE),
122  fFixedSigma2Gaus(kFALSE),
123  fFixedRawYield(-1.),
124  fFrac2Gaus(0.2),
125  fFixedFrac2Gaus(kFALSE),
126  fRatio2GausSigma(0.),
127  fFixedRatio2GausSigma(kFALSE),
128  fNParsSig(3),
129  fNParsBkg(2),
130  fOnlySideBands(kFALSE),
131  fNSigma4SideBands(4.),
132  fFitOption("L,E"),
133  fRawYield(0.),
134  fRawYieldErr(0.),
135  fSigFunc(0x0),
136  fBkgFuncSb(0x0),
137  fBkgFunc(0x0),
138  fBkgFuncRefit(0x0),
139  fReflections(kFALSE),
140  fNParsRfl(0),
141  fRflOverSig(0),
142  fFixRflOverSig(kFALSE),
143  fHistoTemplRfl(0x0),
144  fSmoothRfl(kFALSE),
145  fRawYieldHelp(0),
146  fRflFunc(0x0),
147  fBkRFunc(0x0),
148  fSecondPeak(kFALSE),
149  fNParsSec(0),
150  fSecMass(-999.),
151  fSecWidth(9999.),
152  fFixSecMass(kFALSE),
153  fFixSecWidth(kFALSE),
154  fSecFunc(0x0),
155  fTotFunc(0x0)
156 {
158  fHistoInvMass=(TH1F*)histoToFit->Clone("fHistoInvMass");
159  fHistoInvMass->SetDirectory(0);
161 }
162 //_________________________________________________________________________
164 
166 
167  delete fSigFunc;
168  delete fBkgFunc;
169  delete fBkgFuncSb;
170  delete fBkgFuncRefit;
171  delete fHistoInvMass;
172  delete fRflFunc;
173  delete fBkRFunc;
174  delete fSecFunc;
175  delete fTotFunc;
176  delete fHistoTemplRfl;
177 
178 }
179 //__________________________________________________________________________
183 
184  switch (fTypeOfFit4Bkg) {
185  case 0:
186  fNParsBkg=2;
187  break;
188  case 1:
189  fNParsBkg=2;
190  break;
191  case 2:
192  fNParsBkg=3;
193  break;
194  case 3:
195  fNParsBkg=1;
196  break;
197  case 4:
198  fNParsBkg=2;
199  break;
200  case 5:
201  fNParsBkg=2;
202  break;
203  case 6:
205  break;
206  default:
207  AliError("Error in computing fNParsBkg: check fTypeOfFit4Bkg");
208  break;
209  }
210 
211  switch (fTypeOfFit4Sgn) {
212  case 0:
213  fNParsSig=3;
214  break;
215  case 1:
216  fNParsSig=5;
217  break;
218  case 2:
219  fNParsSig=5;
220  break;
221  default:
222  AliError("Error in computing fNParsSig: check fTypeOfFit4Sgn");
223  break;
224  }
225 
226  if(fReflections) fNParsRfl=1;
227  else fNParsRfl=0;
228 
229  if(fSecondPeak) fNParsSec=3;
230  else fNParsSec=0;
231 
232 }
233 //__________________________________________________________________________
239 
240  TVirtualFitter::SetDefaultFitter("Minuit");
241 
242  Double_t integralHisto=fHistoInvMass->Integral(fHistoInvMass->FindBin(fMinMass),fHistoInvMass->FindBin(fMaxMass),"width");
243 
244  fOnlySideBands = kTRUE;
245  fBkgFuncSb = CreateBackgroundFitFunction("funcbkgsb",integralHisto);
246  Int_t status=-1;
247  printf("\n--- First fit with only background on the side bands - Exclusion region = %.2f sigma ---\n",fNSigma4SideBands);
248  if(fTypeOfFit4Bkg==6){
250  fHistoInvMass->GetListOfFunctions()->Add(fBkgFuncSb);
251  fHistoInvMass->GetFunction(fBkgFuncSb->GetName())->SetBit(1<<9,kTRUE);
252  status=0;
253  }
254  }
255  else status=fHistoInvMass->Fit("funcbkgsb",Form("R,%s,+,0",fFitOption.Data()));
256  fBkgFuncSb->SetLineColor(kGray+1);
257  if (status != 0){
258  printf(" ---> Failed first fit with only background, minuit status = %d\n",status);
259  return 0;
260  }
261 
262  fOnlySideBands = kFALSE;
263  if(!fBkgFunc){
264  fBkgFunc = CreateBackgroundFitFunction("funcbkg",integralHisto);
265  for(Int_t ipar=0; ipar<fNParsBkg; ipar++) fBkgFunc->SetParameter(ipar,fBkgFuncSb->GetParameter(ipar));
266  }
267  fBkgFunc->SetLineColor(kGray+1);
268 
269  printf("\n--- Estimate signal counts in the peak region ---\n");
270  Double_t estimSignal=CheckForSignal(fMass,fSigmaSgn);
271  Bool_t doFinalFit=kTRUE;
272  if(estimSignal<0.){
273  if(draw) DrawFit();
274  estimSignal=0.;
275  doFinalFit=kFALSE;
276  }
277 
278  fRawYieldHelp=estimSignal; // needed for reflection normalization
279  if(!fBkgFuncRefit){
280  fBkgFuncRefit = CreateBackgroundFitFunction("funcbkgrefit",integralHisto);
281  for(Int_t ipar=0; ipar<fNParsBkg; ipar++) fBkgFuncRefit->SetParameter(ipar,fBkgFunc->GetParameter(ipar));
282  }
283  fBkgFuncRefit->SetLineColor(2);
284  fSigFunc = CreateSignalFitFunction("fsigfit",estimSignal);
285  if(fSecondPeak){
286  printf(" ---> Final fit includes a second inv. mass peak\n");
288  fSecFunc = CreateSecondPeakFunction("fsecpeak",estimSec);
289  }
290  if(fReflections){
291  printf(" ---> Final fit includes reflections\n");
292  fRflFunc = CreateReflectionFunction("freflect");
294  }
295  fTotFunc = CreateTotalFitFunction("funcmass");
296 
297  if(doFinalFit){
298  printf("\n--- Final fit with signal+background on the full range ---\n");
299  status=fHistoInvMass->Fit("funcmass",Form("R,%s,+,0",fFitOption.Data()));
300  if (status != 0){
301  printf(" ---> Failed fit with signal+background, minuit status = %d\n",status);
302  return 0;
303  }
304  }
305 
306  for(Int_t ipar=0; ipar<fNParsBkg; ipar++){
307  fBkgFuncRefit->SetParameter(ipar,fTotFunc->GetParameter(ipar));
308  fBkgFuncRefit->SetParError(ipar,fTotFunc->GetParError(ipar));
309  }
310  for(Int_t ipar=0; ipar<fNParsSig; ipar++){
311  fSigFunc->SetParameter(ipar,fTotFunc->GetParameter(ipar+fNParsBkg));
312  fSigFunc->SetParError(ipar,fTotFunc->GetParError(ipar+fNParsBkg));
313  }
314  if(fSecondPeak){
315  for(Int_t ipar=0; ipar<fNParsSec; ipar++){
316  fSecFunc->SetParameter(ipar,fTotFunc->GetParameter(ipar+fNParsBkg+fNParsSig));
317  fSecFunc->SetParError(ipar,fTotFunc->GetParError(ipar+fNParsBkg+fNParsSig));
318  }
319  fSecFunc->SetLineColor(kMagenta+1);
320  fSecFunc->SetLineStyle(3);
321  }
322  if(fReflections){
323  for(Int_t ipar=0; ipar<fNParsRfl; ipar++){
324  fRflFunc->SetParameter(ipar,fTotFunc->GetParameter(ipar+fNParsBkg+fNParsSig+fNParsSec));
325  fRflFunc->SetParError(ipar,fTotFunc->GetParError(ipar+fNParsBkg+fNParsSig+fNParsSec));
326  fBkRFunc->SetParameter(ipar+fNParsBkg,fTotFunc->GetParameter(ipar+fNParsBkg+fNParsSig+fNParsSec));
327  fBkRFunc->SetParError(ipar+fNParsBkg,fTotFunc->GetParError(ipar+fNParsBkg+fNParsSig+fNParsSec));
328  }
329  for(Int_t ipar=0; ipar<fNParsBkg; ipar++){
330  fBkRFunc->SetParameter(ipar,fTotFunc->GetParameter(ipar));
331  fBkRFunc->SetParError(ipar,fTotFunc->GetParError(ipar));
332  }
333  fRflFunc->SetLineColor(kGreen+1);
334  fBkRFunc->SetLineColor(kRed+1);
335  fBkRFunc->SetLineStyle(7);
336  }
337  fMass=fSigFunc->GetParameter(1);
338  fMassErr=fSigFunc->GetParError(1);
339  fSigmaSgn=fSigFunc->GetParameter(2);
340  fSigmaSgnErr=fSigFunc->GetParError(2);
341  fTotFunc->SetLineColor(4);
342  fRawYield=fTotFunc->GetParameter(fNParsBkg)/fHistoInvMass->GetBinWidth(1);
343  fRawYieldErr=fTotFunc->GetParError(fNParsBkg)/fHistoInvMass->GetBinWidth(1);
345  if(draw) DrawFit();
346  if(doFinalFit) return 1;
347  else return 2;
348 }
349 //______________________________________________________________________________
353 
354  TCanvas* c0=new TCanvas("c0");
355  DrawHere(c0);
356 }
357 //______________________________________________________________________________
358 void AliHFInvMassFitter::DrawHere(TVirtualPad* c, Double_t nsigma,Int_t writeFitInfo){
361 
362  gStyle->SetOptStat(0);
363  gStyle->SetCanvasColor(0);
364  gStyle->SetFrameFillColor(0);
365  c->cd();
366  fHistoInvMass->GetXaxis()->SetRangeUser(fMinMass,fMaxMass);
367  fHistoInvMass->SetMarkerStyle(20);
368  fHistoInvMass->SetMinimum(0.);
369  fHistoInvMass->Draw("PE");
370  if(fBkgFunc) fBkgFunc->Draw("same");
371  if(fBkgFuncRefit) fBkgFuncRefit->Draw("same");
372  if(fBkRFunc) fBkRFunc->Draw("same");
373  if(fRflFunc) fRflFunc->Draw("same");
374  if(fSecFunc) fSecFunc->Draw("same");
375  if(fTotFunc) fTotFunc->Draw("same");
376 
377  if(writeFitInfo > 0){
378  TPaveText *pinfos=new TPaveText(0.12,0.65,0.47,0.89,"NDC");
379  TPaveText *pinfom=new TPaveText(0.6,0.7,1.,.87,"NDC");
380  pinfos->SetBorderSize(0);
381  pinfos->SetFillStyle(0);
382  pinfom->SetBorderSize(0);
383  pinfom->SetFillStyle(0);
384  if(fTotFunc){
385  pinfom->SetTextColor(kBlue);
386  for(Int_t ipar=1; ipar<fNParsSig; ipar++){
387  pinfom->AddText(Form("%s = %.3f #pm %.3f",fTotFunc->GetParName(ipar+fNParsBkg),fTotFunc->GetParameter(ipar+fNParsBkg),fTotFunc->GetParError(ipar+fNParsBkg)));
388  }
389  pinfom->Draw();
390 
391  Double_t bkg,errbkg;
392  Background(nsigma,bkg,errbkg);
393  Double_t signif,errsignif;
394  Significance(nsigma,signif,errsignif);
395 
396  pinfos->AddText(Form("S = %.0f #pm %.0f ",fRawYield,fRawYieldErr));
397  pinfos->AddText(Form("B (%.0f#sigma) = %.0f #pm %.0f",nsigma,bkg,errbkg));
398  pinfos->AddText(Form("S/B (%.0f#sigma) = %.4f ",nsigma,fRawYield/bkg));
399  if(fRflFunc) pinfos->AddText(Form("Refl/Sig = %.3f #pm %.3f ",fRflFunc->GetParameter(0),fRflFunc->GetParError(0)));
400  pinfos->AddText(Form("Signif (%.0f#sigma) = %.1f #pm %.1f ",nsigma,signif,errsignif));
401  pinfos->Draw();
402  }
403  }
404  c->Update();
405  return;
406 }
407 //______________________________________________________________________________
411 
412  Double_t minForSig=mean-4.*sigma;
413  Double_t maxForSig=mean+4.*sigma;
414  Int_t binForMinSig=fHistoInvMass->FindBin(minForSig);
415  Int_t binForMaxSig=fHistoInvMass->FindBin(maxForSig);
416  Double_t sum=0.;
417  Double_t sumback=0.;
418  fBkgFunc->Print();
419  for(Int_t ibin=binForMinSig; ibin<=binForMaxSig; ibin++){
420  sum+=fHistoInvMass->GetBinContent(ibin);
421  sumback+=fBkgFunc->Eval(fHistoInvMass->GetBinCenter(ibin));
422  }
423  Double_t diffUnderPeak=(sum-sumback);
424  printf(" ---> IntegralUnderFitFunc=%f IntegralUnderHisto=%f EstimatedSignal=%f\n",sum,sumback,diffUnderPeak);
425  if(diffUnderPeak/TMath::Sqrt(sum)<1.){
426  printf(" ---> (Tot-Bkg)/sqrt(Tot)=%f ---> Likely no signal/\n",diffUnderPeak/TMath::Sqrt(sum));
427  return -1;
428  }
429  return diffUnderPeak*fHistoInvMass->GetBinWidth(1);
430 }
431 
432 //______________________________________________________________________________
436 
438  TF1* funcbkg = new TF1(fname.Data(),this,&AliHFInvMassFitter::FitFunction4Bkg,fMinMass,fMaxMass,fNParsBkg,"AliHFInvMassFitter","FitFunction4Bkg");
439  switch (fTypeOfFit4Bkg) {
440  case 0: //gaus+expo
441  funcbkg->SetParNames("BkgInt","Slope");
442  funcbkg->SetParameters(integral,-2.);
443  break;
444  case 1:
445  funcbkg->SetParNames("BkgInt","Slope");
446  funcbkg->SetParameters(integral,-100.);
447  break;
448  case 2:
449  funcbkg->SetParNames("BkgInt","Coef1","Coef2");
450  funcbkg->SetParameters(integral,-10.,5);
451  break;
452  case 3:
453  funcbkg->SetParNames("Const");
454  funcbkg->SetParameter(0,0.);
455  funcbkg->FixParameter(0,0.);
456  break;
457  case 4:
458  funcbkg->SetParNames("BkgInt","Coef1");
459  funcbkg->SetParameters(integral,0.5);
460  break;
461  case 5:
462  funcbkg->SetParNames("Coef1","Coef2");
463  funcbkg->SetParameters(-10.,5.);
464  break;
465  case 6:
466  for(Int_t j=0;j<fNParsBkg; j++){
467  funcbkg->SetParName(j,Form("Coef%d",j));
468  funcbkg->SetParameter(j,0);
469  }
470  funcbkg->SetParameter(0,integral);
471  break;
472  default:
473  AliError(Form("Wrong choice of fTypeOfFit4Bkg (%d)",fTypeOfFit4Bkg));
474  delete funcbkg;
475  return 0x0;
476  break;
477  }
478  // if(fFixToHistoIntegral) funcbkg->FixParameter(0,integral);
479  funcbkg->SetLineColor(kBlue+3);
480  return funcbkg;
481 }
482 //______________________________________________________________________________
486 
487  TF1* funcsec = new TF1(fname.Data(),this,&AliHFInvMassFitter::FitFunction4SecPeak,fMinMass,fMaxMass,3,"AliHFInvMassFitter","FitFunction4SecPeak");
488  funcsec->SetParameter(0,integsig);
489  funcsec->SetParameter(1,fSecMass);
490  if(fFixSecMass) funcsec->FixParameter(1,fSecMass);
491  funcsec->SetParameter(2,fSecWidth);
492  if(fFixSecWidth) funcsec->FixParameter(2,fSecWidth);
493  funcsec->SetParNames("SecPeakInt","SecMean","SecSigma");
494  return funcsec;
495 }
496 //______________________________________________________________________________
500  TF1* funcrfl = new TF1(fname.Data(),this,&AliHFInvMassFitter::FitFunction4Refl,fMinMass,fMaxMass,1,"AliHFInvMassFitter","FitFunction4Refl");
501  funcrfl->SetParameter(0,fRflOverSig);
502  funcrfl->SetParLimits(0,0.,1.);
503  if(fFixRflOverSig) funcrfl->FixParameter(0,fRflOverSig);
504  funcrfl->SetParNames("ReflOverS");
505  return funcrfl;
506 }
507 //______________________________________________________________________________
512  Int_t totParams=fNParsBkg+fNParsRfl;
513  TF1* fbr=new TF1(fname.Data(),this,&AliHFInvMassFitter::FitFunction4BkgAndRefl,fMinMass,fMaxMass,totParams,"AliHFInvMassFitter","FitFunction4BkgAndRefl");
514  for(Int_t ipar=0; ipar<fNParsBkg; ipar++){
515  fbr->SetParameter(ipar,fBkgFunc->GetParameter(ipar));
516  fbr->SetParName(ipar,fBkgFunc->GetParName(ipar));
517  }
518  for(Int_t ipar=0; ipar<fNParsRfl; ipar++){
519  fbr->SetParameter(ipar+fNParsBkg,fRflFunc->GetParameter(ipar));
520  fbr->SetParName(ipar+fNParsBkg,fRflFunc->GetParName(ipar));
521  // par limits not set because this function is not used for fitting but only for drawing
522  }
523  return fbr;
524 }
525 //______________________________________________________________________________
529 
531  TF1* funcsig = new TF1(fname.Data(),this,&AliHFInvMassFitter::FitFunction4Sgn,fMinMass,fMaxMass,fNParsSig,"AliHFInvMassFitter","FitFunction4Sgn");
532  if(fTypeOfFit4Sgn==kGaus){
533  funcsig->SetParameter(0,integsig);
534  if(fFixedRawYield>-0.1) funcsig->FixParameter(0,fFixedRawYield);
535  funcsig->SetParameter(1,fMass);
536  if(fFixedMean) funcsig->FixParameter(1,fMass);
537  funcsig->SetParameter(2,fSigmaSgn);
538  if(fFixedSigma) funcsig->FixParameter(2,fSigmaSgn);
539  funcsig->SetParNames("SgnInt","Mean","Sigma");
540  }
541  if(fTypeOfFit4Sgn==k2Gaus){
542  funcsig->SetParameter(0,integsig);
543  if(fFixedRawYield>-0.1) funcsig->FixParameter(0,fFixedRawYield);
544  funcsig->SetParameter(1,fMass);
545  if(fFixedMean) funcsig->FixParameter(1,fMass);
546  funcsig->SetParameter(2,fSigmaSgn);
547  funcsig->SetParLimits(2,0.004,0.05);
548  if(fFixedSigma) funcsig->FixParameter(2,fSigmaSgn);
549  funcsig->SetParameter(3,fFrac2Gaus);
550  if(fFixedFrac2Gaus) funcsig->FixParameter(3,fFrac2Gaus);
551  else funcsig->SetParLimits(3,0.,1.);
552  funcsig->SetParameter(4,fSigmaSgn2Gaus);
553  if(fFixedSigma2Gaus) funcsig->FixParameter(4,fSigmaSgn2Gaus);
554  else funcsig->SetParLimits(4,0.004,0.05);
555  funcsig->SetParNames("SgnInt","Mean","Sigma1","Frac","Sigma2");
556  }
558  funcsig->SetParameter(0,integsig);
559  if(fFixedRawYield>-0.1) funcsig->FixParameter(0,fFixedRawYield);
560  funcsig->SetParameter(1,fMass);
561  if(fFixedMean) funcsig->FixParameter(1,fMass);
562  funcsig->SetParameter(2,fSigmaSgn);
563  funcsig->SetParLimits(2,0.004,0.05);
564  if(fFixedSigma) funcsig->FixParameter(2,fSigmaSgn);
565  funcsig->SetParameter(3,fFrac2Gaus);
566  if(fFixedFrac2Gaus) funcsig->FixParameter(3,fFrac2Gaus);
567  else funcsig->SetParLimits(3,0.,1.);
568  funcsig->SetParameter(4,fSigmaSgn2Gaus);
569  if(fFixedRatio2GausSigma) funcsig->FixParameter(4,fRatio2GausSigma);
570  else funcsig->SetParLimits(4,0.,20.);
571  funcsig->SetParNames("SgnInt","Mean","Sigma1","Frac","RatioSigma12");
572  }
573  return funcsig;
574 }
575 
576 //______________________________________________________________________________
580 
583  TF1* ftot=new TF1(fname.Data(),this,&AliHFInvMassFitter::FitFunction4Mass,fMinMass,fMaxMass,totParams,"AliHFInvMassFitter","FitFunction4Mass");
584  for(Int_t ipar=0; ipar<fNParsBkg; ipar++){
585  ftot->SetParameter(ipar,fBkgFunc->GetParameter(ipar));
586  ftot->SetParName(ipar,fBkgFunc->GetParName(ipar));
587  // Double_t parmin,parmax;
588  // fBkgFunc->GetParLimits(ipar,parmin,parmax);
589  // ftot->SetParLimits(ipar,parmin,parmax);
590  }
591  for(Int_t ipar=0; ipar<fNParsSig; ipar++){
592  ftot->SetParameter(ipar+fNParsBkg,fSigFunc->GetParameter(ipar));
593  ftot->SetParName(ipar+fNParsBkg,fSigFunc->GetParName(ipar));
594  Double_t parmin,parmax;
595  fSigFunc->GetParLimits(ipar,parmin,parmax);
596  ftot->SetParLimits(ipar+fNParsBkg,parmin,parmax);
597  }
598  if(fSecondPeak && fSecFunc){
599  for(Int_t ipar=0; ipar<fNParsSec; ipar++){
600  ftot->SetParameter(ipar+fNParsBkg+fNParsSig,fSecFunc->GetParameter(ipar));
601  ftot->SetParName(ipar+fNParsBkg+fNParsSig,fSecFunc->GetParName(ipar));
602  Double_t parmin,parmax;
603  fSecFunc->GetParLimits(ipar,parmin,parmax);
604  ftot->SetParLimits(ipar+fNParsBkg+fNParsSig,parmin,parmax);
605  }
606  }
607  if(fReflections && fRflFunc){
608  for(Int_t ipar=0; ipar<fNParsRfl; ipar++){
609  ftot->SetParameter(ipar+fNParsBkg+fNParsSig+fNParsSec,fRflFunc->GetParameter(ipar));
610  ftot->SetParName(ipar+fNParsBkg+fNParsSig+fNParsSec,fRflFunc->GetParName(ipar));
611  Double_t parmin,parmax;
612  fRflFunc->GetParLimits(ipar,parmin,parmax);
613  ftot->SetParLimits(ipar+fNParsBkg+fNParsSig+fNParsSec,parmin,parmax);
614  }
615  }
616  return ftot;
617 }
618 //__________________________________________________________________________
622 
624  if(fOnlySideBands && TMath::Abs(x[0]-fMass) < maxDeltaM) {
625  TF1::RejectPoint();
626  return 0;
627  }
628  if(fOnlySideBands && fSecondPeak && TMath::Abs(x[0]-fSecMass) < (fNSigma4SideBands*fSecWidth)){
629  TF1::RejectPoint();
630  return 0;
631  }
632  Double_t total=0;
633 
634  switch (fTypeOfFit4Bkg){
635  case 0:
636  //exponential
637  //exponential = A*exp(B*x) -> integral(exponential)=A/B*exp(B*x)](min,max)
638  //-> A = B*integral/(exp(B*max)-exp(B*min)) where integral can be written
639  //as integralTot- integralGaus (=par [2])
640  //Par:
641  // * [0] = integralBkg;
642  // * [1] = B;
643  //exponential = [1]*[0]/(exp([1]*max)-exp([1]*min))*exp([1]*x)
644  total = par[0]*par[1]/(TMath::Exp(par[1]*fMaxMass)-TMath::Exp(par[1]*fMinMass))*TMath::Exp(par[1]*x[0]);
645  // AliInfo("Background function set to: exponential");
646  break;
647  case 1:
648  //linear
649  //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)
650  // * [0] = integralBkg;
651  // * [1] = b;
652  total= par[0]/(fMaxMass-fMinMass)+par[1]*(x[0]-0.5*(fMaxMass+fMinMass));
653  // AliInfo("Background function set to: linear");
654  break;
655  case 2:
656  //parabola
657  //y=a+b*x+c*x**2 -> integral = a(max-min) + 1/2*b*(max^2-min^2) +
658  //+ 1/3*c*(max^3-min^3) ->
659  //a = (integral-1/2*b*(max^2-min^2)-1/3*c*(max^3-min^3))/(max-min)
660  // * [0] = integralBkg;
661  // * [1] = b;
662  // * [2] = c;
663  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));
664  // AliInfo("Background function set to: polynomial");
665  break;
666  case 3:
667  total=par[0];
668  break;
669  case 4:
670  //power function
671  //y=a(x-m_pi)^b -> integral = a/(b+1)*((max-m_pi)^(b+1)-(min-m_pi)^(b+1))
672  //
673  //a = integral*(b+1)/((max-m_pi)^(b+1)-(min-m_pi)^(b+1))
674  // * [0] = integralBkg;
675  // * [1] = b;
676  // a(power function) = [0]*([1]+1)/((max-m_pi)^([1]+1)-(min-m_pi)^([1]+1))*(x-m_pi)^[1]
677  {
678  Double_t mpi = TDatabasePDG::Instance()->GetParticle(211)->Mass();
679 
680  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]);
681  // AliInfo("Background function set to: powerlaw");
682  }
683  break;
684  case 5:
685  //power function wit exponential
686  //y=a*Sqrt(x-m_pi)*exp(-b*(x-m_pi))
687  {
688  Double_t mpi = TDatabasePDG::Instance()->GetParticle(211)->Mass();
689 
690  total = par[0]*TMath::Sqrt(x[0] - mpi)*TMath::Exp(-1.*par[1]*(x[0]-mpi));
691  // AliInfo("Background function set to: wit exponential");
692  }
693  break;
694  case 6:
695  // the following comment must be removed
696  // // pol 3, following convention for pol 2
697  // //y=a+b*x+c*x**2+d*x**3 -> integral = a(max-min) + 1/2*b*(max^2-min^2) +
698  // //+ 1/3*c*(max^3-min^3) + 1/4 d * (max^4-min^4) ->
699  // //a = (integral-1/2*b*(max^2-min^2)-1/3*c*(max^3-min^3) - 1/4 d * (max^4-min^4) )/(max-min)
700  // // * [0] = integralBkg;
701  // // * [1] = b;
702  // // * [2] = c;
703  // // * [3] = d;
704  {
705  total=par[0];
706  for(Int_t it=1;it<=fPolDegreeBkg;it++){
707  total+=par[it]*TMath::Power(x[0]-fMassParticle,it)/TMath::Factorial(it);
708  }
709  }
710  break;
711  }
712  return total;
713 }
714 //_________________________________________________________________________
718 
719  // AliInfo("Signal function set to: Gaussian");
720  Double_t sigval=0;
721  Double_t g1=0;
722  Double_t g2=0;
723  switch (fTypeOfFit4Sgn){
724  case 0:
725  //gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
726  //Par:
727  // * [0] = integralSgn
728  // * [1] = mean
729  // * [2] = sigma
730  //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
731  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]);
732  break;
733  case 1:
734  //double gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
735  //Par:
736  // * [0] = integralSgn
737  // * [1] = mean
738  // * [2] = sigma1
739  // * [3] = 2nd gaussian ratio
740  // * [4] = deltaSigma
741  //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
742  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]);
743  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]);
744  sigval=par[0]*(g1+g2);
745  break;
746  case 2:
747  //double gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
748  //Par:
749  // * [0] = integralSgn
750  // * [1] = mean
751  // * [2] = sigma1
752  // * [3] = 2nd gaussian ratio
753  // * [4] = ratio sigma12
754  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]);
755  g2=par[3]/TMath::Sqrt(2.*TMath::Pi())/(par[4]*par[2])*TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/2./(par[4]*par[2])/(par[4]*par[2]));
756  sigval=par[0]*(g1+g2);
757  break;
758  }
759  fRawYieldHelp=par[0]/fHistoInvMass->GetBinWidth(1);
760  return sigval;
761 }
762 //_________________________________________________________________________
766  if(!fHistoTemplRfl) return 0;
767 
768  Int_t bin =fHistoTemplRfl->FindBin(x[0]);
769  Double_t value=fHistoTemplRfl->GetBinContent(bin);
770  Int_t binmin=fHistoTemplRfl->FindBin(fMinMass*1.00001);
771  Int_t binmax=fHistoTemplRfl->FindBin(fMaxMass*0.99999);
772  Double_t norm=fHistoTemplRfl->Integral(binmin,binmax)*fHistoTemplRfl->GetBinWidth(bin);
773  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
774  value+=fHistoTemplRfl->GetBinContent(bin-1)+fHistoTemplRfl->GetBinContent(bin+1);
775  value/=3.;
776  }
777  return par[0]*value/norm*fRawYieldHelp*fHistoInvMass->GetBinWidth(1);
778 }
779 //_________________________________________________________________________
783  Double_t bkg=FitFunction4Bkg(x,par);
784  Double_t refl=0;
785  if(fReflections) refl=FitFunction4Refl(x,&par[fNParsBkg]);
786  return bkg+refl;
787 }
788 //_________________________________________________________________________
792 
793  //gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
794  //Par:
795  // * [0] = integralSgn
796  // * [1] = mean
797  // * [2] = sigma
798  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]);
799  return secgaval;
800 }
801 //_________________________________________________________________________
805 
806  Double_t bkg=FitFunction4Bkg(x,par);
807  Double_t sig=FitFunction4Sgn(x,&par[fNParsBkg]);
808  Double_t sec=0.;
809  if(fSecondPeak) sec=FitFunction4SecPeak(x,&par[fNParsBkg+fNParsSig]);
810  Double_t refl=0;
811  if(fReflections) refl=FitFunction4Refl(x,&par[fNParsBkg+fNParsSig+fNParsSec]);
812  return bkg+sig+sec+refl;
813 }
814 
815 //_________________________________________________________________________
816 void AliHFInvMassFitter::Signal(Double_t nOfSigma,Double_t &signal,Double_t &errsignal) const {
819 
820  Double_t minMass=fMass-nOfSigma*fSigmaSgn;
821  Double_t maxMass=fMass+nOfSigma*fSigmaSgn;
822  Signal(minMass,maxMass,signal,errsignal);
823  return;
824 }
825 
826 //_________________________________________________________________________
827 void AliHFInvMassFitter::Signal(Double_t min, Double_t max, Double_t &signal,Double_t &errsignal) const {
830 
831  signal=fSigFunc->Integral(min, max)/(Double_t)fHistoInvMass->GetBinWidth(1);
832  errsignal=(fRawYieldErr/fRawYield)*signal;/*assume relative error is the same as for total integral*/
833  return;
834 }
835 
836 //___________________________________________________________________________
837 void AliHFInvMassFitter::Background(Double_t nOfSigma,Double_t &background,Double_t &errbackground) const {
840 
841  Double_t minMass=fMass-nOfSigma*fSigmaSgn;
842  Double_t maxMass=fMass+nOfSigma*fSigmaSgn;
843  Background(minMass,maxMass,background,errbackground);
844  return;
845 
846 }
847 //___________________________________________________________________________
848 void AliHFInvMassFitter::Background(Double_t min, Double_t max, Double_t &background,Double_t &errbackground) const {
851 
852  TF1 *funcbkg=0x0;
853  if(fBkgFuncRefit) funcbkg=fBkgFuncRefit;
854  else if(fBkgFunc) funcbkg=fBkgFunc;
855  if(!funcbkg){
856  AliError("Bkg function not found!");
857  return;
858  }
859 
860  Double_t intB=funcbkg->GetParameter(0);
861  Double_t intBerr=funcbkg->GetParError(0);
862 
863  //relative error evaluation: from histo
864 
866  Int_t rightBand=fHistoInvMass->FindBin(fMass+fNSigma4SideBands*fSigmaSgn);
867  intB=fHistoInvMass->Integral(1,leftBand)+fHistoInvMass->Integral(rightBand,fHistoInvMass->GetNbinsX());
868  Double_t sum2=0;
869  for(Int_t i=1;i<=leftBand;i++){
870  sum2+=fHistoInvMass->GetBinError(i)*fHistoInvMass->GetBinError(i);
871  }
872  for(Int_t i=rightBand; i<=fHistoInvMass->GetNbinsX();i++){
873  sum2+=fHistoInvMass->GetBinError(i)*fHistoInvMass->GetBinError(i);
874  }
875 
876  intBerr=TMath::Sqrt(sum2);
877 
878  background=funcbkg->Integral(min,max)/(Double_t)fHistoInvMass->GetBinWidth(1);
879  errbackground=intBerr/intB*background;
880 
881  return;
882 
883 }
884 //__________________________________________________________________________
885 
886 void AliHFInvMassFitter::Significance(Double_t nOfSigma,Double_t &significance,Double_t &errsignificance) const {
889 
890  Double_t minMass=fMass-nOfSigma*fSigmaSgn;
891  Double_t maxMass=fMass+nOfSigma*fSigmaSgn;
892  Significance(minMass, maxMass, significance, errsignificance);
893 
894  return;
895 }
896 
897 //__________________________________________________________________________
898 
899 void AliHFInvMassFitter::Significance(Double_t min, Double_t max, Double_t &significance,Double_t &errsignificance) const {
902 
903  Double_t background,errbackground;
904  Background(min,max,background,errbackground);
905 
906  if (fRawYield+background <= 0.){
907  significance=-1;
908  errsignificance=0;
909  return;
910  }
911 
912  AliVertexingHFUtils::ComputeSignificance(fRawYield,fRawYieldErr,background,errbackground,significance,errsignificance);
913 
914  return;
915 }
916 //________________________________________________________________________
920 
921  TH1F *hCp=(TH1F*)fHistoInvMass->Clone("htemp");
922  Double_t estimatecent=0.5*(hCp->GetBinContent(hCp->FindBin(fMass-3.5*fSigmaSgn))+hCp->GetBinContent(hCp->FindBin(fMass+3.5*fSigmaSgn)));// just a first rough estimate
923  Double_t estimateslope=(hCp->GetBinContent(hCp->FindBin(fMass+3.5*fSigmaSgn))-hCp->GetBinContent(hCp->FindBin(fMass-3.5*fSigmaSgn)))/(7*fSigmaSgn);// first rough estimate
924 
926  TF1 *funcbkg,*funcPrev=0x0;
928  funcbkg = new TF1(Form("temp%d",fCurPolDegreeBkg),this,&AliHFInvMassFitter::BackFitFuncPolHelper,fMinMass,fMaxMass,fCurPolDegreeBkg+1,"AliHFInvMassFitter","BackFitFuncPolHelper");
929  if(funcPrev){
930  for(Int_t j=0;j<fCurPolDegreeBkg;j++){// now is +1 degree w.r.t. previous fit funct
931  funcbkg->SetParameter(j,funcPrev->GetParameter(j));
932  }
933  delete funcPrev;
934  }
935  else{
936  funcbkg->SetParameter(0,estimatecent);
937  funcbkg->SetParameter(1,estimateslope);
938  }
939  printf(" ---> Pre-fit of background with pol degree %d ---\n",fCurPolDegreeBkg);
940  hCp->Fit(funcbkg,"REMN","");
941  funcPrev=(TF1*)funcbkg->Clone("ftemp");
942  delete funcbkg;
944  }
945 
946  for(Int_t j=0;j<=fPolDegreeBkg;j++){
947  fback->SetParameter(j,funcPrev->GetParameter(j));
948  fback->SetParError(j,funcPrev->GetParError(j));
949  }
950  printf(" ---> Final background fit with pol degree %d ---\n",fPolDegreeBkg);
951  hCp->Fit(fback,"REMN","");// THIS IS JUST TO SET NOT ONLY THE PARAMETERS BUT ALSO chi2, etc...
952 
953 
954  // The following lines might be useful for debugging
955  // TCanvas *cDebug=new TCanvas();
956  // cDebug->cd();
957  // hCp->Draw();
958  // TString strout=Form("Test%d.root",(Int_t)fhistoInvMass->GetBinContent(fhistoInvMass->FindBin(fMass)));
959  // cDebug->Print(strout.Data());
960  // delete cDebug;
961 
962  delete funcPrev;
963  delete hCp;
964  return kTRUE;
965 
966 }
967 // _______________________________________________________________________
971 
973  if(fOnlySideBands && TMath::Abs(x[0]-fMass) < maxDeltaM) {
974  TF1::RejectPoint();
975  return 0;
976  }
977  Double_t back=par[0];
978  for(Int_t it=1;it<=fCurPolDegreeBkg;it++){
979  back+=par[it]*TMath::Power(x[0]-fMassParticle,it)/TMath::Factorial(it);
980  }
981  return back;
982 }
983 // _______________________________________________________________________
992 
993  if(!h){
994  fReflections=kFALSE;
995  return 0x0;
996  }
997  fHistoTemplRfl=(TH1F*)h->Clone("hTemplRfl");
998  opt.ToLower();
999  fReflections=kTRUE;
1000  printf("\n--- Reflection templates from simulation ---\n");
1001  if(opt.Contains("templ")){
1002  printf(" ---> Reflection contribution using directly the histogram from simulation\n");
1003  return fHistoTemplRfl;
1004  }
1005 
1006  TF1 *f=0x0;
1007  Bool_t isPoissErr=kTRUE;
1008  Double_t xMinForFit=h->GetBinLowEdge(1);
1009  Double_t xMaxForFit=h->GetXaxis()->GetBinUpEdge(h->GetNbinsX());
1010  if(minRange>=0 && maxRange>=0){
1011  xMinForFit=TMath::Max(minRange,h->GetBinLowEdge(1));
1012  xMaxForFit=TMath::Min(maxRange,h->GetXaxis()->GetBinUpEdge(h->GetNbinsX()));
1013  }
1014  if(opt.EqualTo("1gaus") || opt.EqualTo("singlegaus")){
1015  printf(" ---> Reflection contribution from single-Gaussian fit to histogram from simulation\n");
1016  f=new TF1("mygaus","gaus",xMinForFit,xMaxForFit);
1017  f->SetParameter(0,h->GetMaximum());
1018  // f->SetParLimits(0,0,100.*h->Integral());
1019  f->SetParameter(1,1.865);
1020  f->SetParameter(2,0.050);
1021  fHistoTemplRfl->Fit(f,"REM0","");//,h->GetBinLowEdge(1),h->GetXaxis()->GetBinUpEdge(h->GetNbinsX()));
1022  }
1023  else if(opt.EqualTo("2gaus") || opt.EqualTo("doublegaus")){
1024  printf(" ---> Reflection contribution from double-Gaussian fit to histogram from simulation\n");
1025  f=new TF1("my2gaus","[0]*([3]/( TMath::Sqrt(2.*TMath::Pi())*[2])*TMath::Exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))+(1.-[3])/( TMath::Sqrt(2.*TMath::Pi())*[5])*TMath::Exp(-(x-[4])*(x-[4])/(2.*[5]*[5])))",xMinForFit,xMaxForFit);
1026  f->SetParameter(0,h->GetMaximum());
1027  // f->SetParLimits(0,0,100.*h->Integral());
1028  f->SetParLimits(3,0.,1.);
1029  f->SetParameter(3,0.5);
1030  f->SetParameter(1,1.84);
1031  f->SetParameter(2,0.050);
1032  f->SetParameter(4,1.88);
1033  f->SetParameter(5,0.050);
1034  fHistoTemplRfl->Fit(f,"REM0","");//,h->GetBinLowEdge(1),h->GetXaxis()->GetBinUpEdge(h->GetNbinsX()));
1035  }
1036  else if(opt.EqualTo("pol3")){
1037  printf(" ---> Reflection contribution from pol3 fit to histogram from simulation\n");
1038  f=new TF1("mypol3","pol3",xMinForFit,xMaxForFit);
1039  f->SetParameter(0,h->GetMaximum());
1040  // f->SetParLimits(0,0,100.*h->Integral());
1041  // Hard to initialize the other parameters...
1042  fHistoTemplRfl->Fit(f,"REM0","");
1043  // Printf("We USED %d POINTS in the Fit",f->GetNumberFitPoints());
1044  }
1045  else if(opt.EqualTo("pol6")){
1046  printf(" ---> Reflection contribution from pol6 fit to histogram from simulation\n");
1047  f=new TF1("mypol6","pol6",xMinForFit,xMaxForFit);
1048  f->SetParameter(0,h->GetMaximum());
1049  // f->SetParLimits(0,0,100.*h->Integral());
1050  // Hard to initialize the other parameters...
1051  fHistoTemplRfl->Fit(f,"RLEMI0","");//,h->GetBinLowEdge(1),h->GetXaxis()->GetBinUpEdge(h->GetNbinsX()));
1052  }
1053  else{
1054  // no good option passed
1055  printf(" ---> Bad option for reflection configuration -> reflections will not be included in the fit\n");
1056  fReflections=kFALSE;
1057  delete fHistoTemplRfl;
1058  fHistoTemplRfl=0x0;
1059  return 0x0;
1060  }
1061 
1062  // Fill fHistoTemplRfl with values of fit function
1063  if(f){
1064  for(Int_t j=1;j<=fHistoTemplRfl->GetNbinsX();j++){
1065  fHistoTemplRfl->SetBinContent(j,f->Integral(fHistoTemplRfl->GetBinLowEdge(j),fHistoTemplRfl->GetXaxis()->GetBinUpEdge(j))/fHistoTemplRfl->GetBinWidth(j));
1066  if(fHistoTemplRfl->GetBinContent(j)>=0.&&TMath::Abs(h->GetBinError(j)*h->GetBinError(j)-h->GetBinContent(j))>0.1*h->GetBinContent(j))isPoissErr=kFALSE;
1067  }
1068  for(Int_t j=1;j<=fHistoTemplRfl->GetNbinsX();j++){
1069  if(isPoissErr){
1070  fHistoTemplRfl->SetBinError(j,TMath::Sqrt(fHistoTemplRfl->GetBinContent(j)));
1071  }
1072  else fHistoTemplRfl->SetBinError(j,0.001*fHistoTemplRfl->GetBinContent(j));
1073  }
1074  fReflections=kTRUE;
1075  return fHistoTemplRfl;
1076  }else{
1077  printf(" ---> Fit to MC template for reflection failed -> reflections will not be included in the fit\n");
1078  fReflections=kFALSE;
1079  delete fHistoTemplRfl;
1080  fHistoTemplRfl=0x0;
1081  return 0x0;
1082  }
1083  return 0x0;
1084 }
1085 // _______________________________________________________________________
1090  // else (default) mean of gaussian fit
1091 
1092  Double_t massVal=fMass;
1093  switch (pdgCode) {
1094  case 411:
1095  case 421:
1096  case 431:
1097  case 4122:
1098  massVal=TDatabasePDG::Instance()->GetParticle(pdgCode)->Mass();
1099  break;
1100  case 413:
1101  massVal=TDatabasePDG::Instance()->GetParticle(pdgCode)->Mass();
1102  massVal-=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1103  break;
1104  default:
1105  massVal=fMass;
1106  break;
1107  }
1108 
1109  Double_t minMass=massVal-nOfSigma*fSigmaSgn;
1110  Double_t maxMass=massVal+nOfSigma*fSigmaSgn;
1111  return GetRawYieldBinCounting(errRyBC,minMass,maxMass,option);
1112 }
1113 // _______________________________________________________________________
1119 
1120  Int_t minBinSum=fHistoInvMass->FindBin(minMass);
1121  Int_t maxBinSum=fHistoInvMass->FindBin(maxMass);
1122  if(minBinSum<1){
1123  printf("Left range for bin counting smaller than allowed by histogram axis, setting it to the lower edge of the first histo bin\n");
1124  minBinSum=1;
1125  }
1126  if(maxBinSum>fHistoInvMass->GetNbinsX()){
1127  printf("Right range for bin counting larger than allowed by histogram axis, setting it to the upper edge of the last histo bin\n");
1128  maxBinSum=fHistoInvMass->GetNbinsX();
1129  }
1130  Double_t cntSig=0.;
1131  Double_t cntErr=0.;
1132  errRyBC=0;
1133  TF1* fbackground=fBkgFunc;
1134  if(option==1) fbackground=fBkgFuncRefit;
1135  if(!fbackground) return 0.;
1136 
1137  for(Int_t jb=minBinSum; jb<=maxBinSum; jb++){
1138  Double_t cntTot=fHistoInvMass->GetBinContent(jb);
1139  Double_t cntBkg=fbackground->Integral(fHistoInvMass->GetBinLowEdge(jb),fHistoInvMass->GetBinLowEdge(jb)+fHistoInvMass->GetBinWidth(jb))/fHistoInvMass->GetBinWidth(jb);
1140  Double_t cntRefl=0;
1141  if(option==1 && fRflFunc) cntRefl=fRflFunc->Integral(fHistoInvMass->GetBinLowEdge(jb),fHistoInvMass->GetBinLowEdge(jb)+fHistoInvMass->GetBinWidth(jb))/fHistoInvMass->GetBinWidth(jb);
1142  //Double_t cntBkg=fbackground->Eval(fHistoInvMass->GetBinCenter(jb));
1143  Double_t cntSecPeak=0;
1144  if(option==1 && fSecondPeak && fSecFunc) cntSecPeak=fSecFunc->Integral(fHistoInvMass->GetBinLowEdge(jb),fHistoInvMass->GetBinLowEdge(jb)+fHistoInvMass->GetBinWidth(jb))/fHistoInvMass->GetBinWidth(jb);
1145  cntSig+=(cntTot-cntBkg-cntRefl-cntSecPeak);
1146  cntErr+=(fHistoInvMass->GetBinError(jb)*fHistoInvMass->GetBinError(jb));
1147  }
1148  errRyBC=TMath::Sqrt(cntErr);
1149  return cntSig;
1150 }
1151 
1152 
1153 // _______________________________________________________________________
1154 TH1F* AliHFInvMassFitter::GetResidualsAndPulls(TH1 *hPulls,TH1 *hResidualTrend,TH1 *hPullsTrend, Double_t minrange,Double_t maxrange, Int_t option){
1155 
1157 
1158  Int_t binmi=fHistoInvMass->FindBin(fMinMass*1.001);
1159  Int_t binma=fHistoInvMass->FindBin(fMaxMass*0.9999);
1160  if(maxrange>minrange){
1161  binmi=fHistoInvMass->FindBin(minrange*1.001);
1162  binma=fHistoInvMass->FindBin(maxrange*0.9999);
1163  }
1164  if(hResidualTrend){
1165  //fHistoInvMass->Copy(hResidualTrend);
1166  hResidualTrend->SetBins(fHistoInvMass->GetNbinsX(),fHistoInvMass->GetXaxis()->GetXmin(),fHistoInvMass->GetXaxis()->GetXmax());
1167  hResidualTrend->SetName(Form("%s_residualTrend",fHistoInvMass->GetName()));
1168  hResidualTrend->SetTitle(Form("%s (Residuals)",fHistoInvMass->GetTitle()));
1169  hResidualTrend->SetMarkerStyle(20);
1170  hResidualTrend->SetMarkerSize(1.0);
1171  hResidualTrend->Reset();
1172  }
1173  if(hPullsTrend){
1174  hPullsTrend->SetBins(fHistoInvMass->GetNbinsX(),fHistoInvMass->GetXaxis()->GetXmin(),fHistoInvMass->GetXaxis()->GetXmax());
1175  hPullsTrend->Reset();
1176  hPullsTrend->SetName(Form("%s_pullTrend",fHistoInvMass->GetName()));
1177  hPullsTrend->SetTitle(Form("%s (Pulls)",fHistoInvMass->GetTitle()));
1178  hPullsTrend->SetMarkerStyle(20);
1179  hPullsTrend->SetMarkerSize(1.0);
1180  }
1181  if(hPulls){
1182  hPulls->SetName(Form("%s_pulls",fHistoInvMass->GetName()));
1183  hPulls->SetTitle(Form("%s ; Pulls",fHistoInvMass->GetTitle()));
1184  hPulls->SetBins(40,-10,10);
1185  hPulls->Reset();
1186  }
1187 
1188  Double_t res=-1.e-6,min=1.e+12,max=-1.e+12;
1189  TArrayD *arval=new TArrayD(binma-binmi+1);
1190  for(Int_t jst=1;jst<=fHistoInvMass->GetNbinsX();jst++){
1191  Double_t integFit=0;
1192  if(option==0) integFit=fTotFunc->Integral(fHistoInvMass->GetBinLowEdge(jst),fHistoInvMass->GetBinLowEdge(jst)+fHistoInvMass->GetBinWidth(jst));
1193  else{
1194  integFit=fBkgFuncRefit->Integral(fHistoInvMass->GetBinLowEdge(jst),fHistoInvMass->GetBinLowEdge(jst)+fHistoInvMass->GetBinWidth(jst));
1195  if(option==2) integFit+=fRflFunc->Integral(fHistoInvMass->GetBinLowEdge(jst),fHistoInvMass->GetBinLowEdge(jst)+fHistoInvMass->GetBinWidth(jst));
1196  }
1197  res=fHistoInvMass->GetBinContent(jst)-integFit/fHistoInvMass->GetBinWidth(jst);
1198  if(jst>=binmi&&jst<=binma){
1199  arval->AddAt(res,jst-binmi);
1200  if(res<min)min=res;
1201  if(res>max)max=res;
1202  }
1203  // Printf("Res = %f from %f - %f",res,fHistoInvMass->GetBinContent(jst),fTotFunc->Integral(fHistoInvMass->GetBinLowEdge(jst),fHistoInvMass->GetBinLowEdge(jst)+fHistoInvMass->GetBinWidth(jst))/fHistoInvMass->GetBinWidth(jst));
1204  if(hResidualTrend){
1205  hResidualTrend->SetBinContent(jst,res);
1206  hResidualTrend->SetBinError(jst,fHistoInvMass->GetBinError(jst));
1207  }
1208  if(hPulls){
1209  if(jst>=binmi&&jst<=binma)hPulls->Fill(res/fHistoInvMass->GetBinError(jst));
1210  }
1211  if(hPullsTrend){
1212  hPullsTrend->SetBinContent(jst,res/fHistoInvMass->GetBinError(jst));
1213  hPullsTrend->SetBinError(jst,0.0001);
1214  }
1215  }
1216  if(hResidualTrend){
1217  hResidualTrend->GetXaxis()->SetRange(binmi,binma);
1218  if(option!=0){
1219  TF1 *fgauss=new TF1("signalTermForRes","[0]/TMath::Sqrt(2.*TMath::Pi())/[2]*TMath::Exp(-(x-[1])*(x-[1])/2./[2]/[2])",fHistoInvMass->GetBinLowEdge(1),fHistoInvMass->GetBinLowEdge(fHistoInvMass->GetNbinsX()+1));
1220  fgauss->SetParameter(0,fRawYield*fHistoInvMass->GetBinWidth(1));
1221  fgauss->SetParameter(1,fMass);
1222  fgauss->SetParameter(2,fSigmaSgn);
1223  fgauss->SetLineColor(kBlue);
1224  hResidualTrend->GetListOfFunctions()->Add(fgauss);
1225  }
1226  }
1227  if(hPullsTrend){
1228  hPullsTrend->GetXaxis()->SetRange(binmi,binma);
1229  hPullsTrend->SetMinimum(-7);
1230  hPullsTrend->SetMaximum(+7);
1231  }
1232  if(TMath::Abs(min)>TMath::Abs(max))max=min;
1233 
1234  TH1F *hout=new TH1F(Form("%s_residuals",fHistoInvMass->GetName()),Form("%s ; residuals",fHistoInvMass->GetTitle()),25,-TMath::Abs(max)*1.5,TMath::Abs(max)*1.5);
1235  for(Int_t j=0;j<binma-binmi+1;j++){
1236  hout->Fill(arval->At(j));
1237  }
1238  hout->Sumw2();
1239  hout->Fit("gaus","LEM","",-TMath::Abs(max)*1.2,TMath::Abs(max)*1.2);
1240 
1241  if(hPulls){
1242  hPulls->Sumw2();
1243  hPulls->Fit("gaus","LEM","",-3,3);
1244  }
1245  delete arval;
1246  return hout;
1247 }
1248 // _______________________________________________________________________
1249 TH1F* AliHFInvMassFitter::GetOverBackgroundResidualsAndPulls(TH1 *hPulls,TH1 *hResidualTrend,TH1 *hPullsTrend, Double_t minrange,Double_t maxrange){
1251  return GetResidualsAndPulls(hPulls,hResidualTrend,hPullsTrend,minrange,maxrange,1);
1252 }
1253 // _______________________________________________________________________
1254 TH1F* AliHFInvMassFitter::GetOverBackgroundPlusReflResidualsAndPulls(TH1 *hPulls,TH1 *hResidualTrend,TH1 *hPullsTrend, Double_t minrange,Double_t maxrange){
1256  return GetResidualsAndPulls(hPulls,hResidualTrend,hPullsTrend,minrange,maxrange,2);
1257 }
1258 
1259 // _______________________________________________________________________
1263  if(fBkgFunc){
1264  printf("--- Background function in 1st step fit to side bands ---\n");
1265  fBkgFunc->Print();
1266  }
1267  if(fTotFunc){
1268  printf("--- Total fit function from 2nd step fit ---\n");
1269  fTotFunc->Print();
1270  }
1271  if(fBkgFuncRefit){
1272  printf("--- Background function in 2nd step fit ---\n");
1273  fBkgFuncRefit->Print();
1274  }
1275  if(fRflFunc){
1276  printf("--- Reflections ---\n");
1277  fRflFunc->Print();
1278  }
1279  if(fBkRFunc){
1280  printf("--- Background + reflections ---\n");
1281  fBkRFunc->Print();
1282  }
1283  if(fSecFunc){
1284  printf("--- Additional Gaussian peak ---\n");
1285  fSecFunc->Print();
1286  }
1287 }
Int_t fTypeOfFit4Sgn
pdg value of particle mass
Double_t fSigmaSgnErr
signal gaussian sigma
TH1F * GetResidualsAndPulls(TH1 *hPulls=0x0, TH1 *hResidualTrend=0x0, TH1 *hPullsTrend=0x0, Double_t minrange=0, Double_t maxrange=-1, Int_t option=0)
Bool_t fFixSecWidth
flag to fix the position of the 2nd peak
Double_t fRawYieldHelp
switch for smoothing of reflection template
double Double_t
Definition: External.C:58
Double_t fRawYield
L, LW or Chi2.
void DrawHere(TVirtualPad *c, Double_t nsigma=3, Int_t writeFitInfo=1)
TF1 * CreateReflectionFunction(TString fname)
Int_t fCurPolDegreeBkg
degree of polynomial expansion for back fit (option 6 for back)
Bool_t PrepareHighPolFit(TF1 *fback)
Bool_t fOnlySideBands
fit parameters in background fit function
virtual void Signal(Double_t nOfSigma, Double_t &signal, Double_t &errsignal) const
Bool_t fFixedRatio2GausSigma
initialization for ratio between two gaussian sigmas in case of k2GausSigmaRatioPar ...
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
Int_t MassFitter(Bool_t draw=kTRUE)
Int_t fNParsBkg
fit parameters in signal fit function
AliHFInvMassFitter class for the fit of invariant mass distribution of charm hadrons.
Int_t fNParsRfl
flag use/not use reflections
TCanvas * c
Definition: TestFitELoss.C:172
TF1 * CreateSecondPeakFunction(TString fname, Double_t integral)
Bool_t fFixedFrac2Gaus
initialization for fraction of 2nd gaussian in case of k2Gaus or k2GausSigmaRatioPar ...
TH1F * GetOverBackgroundPlusReflResidualsAndPulls(TH1 *hPulls=0x0, TH1 *hResidualTrend=0x0, TH1 *hPullsTrend=0x0, Double_t minrange=0, Double_t maxrange=-1)
TF1 * fTotFunc
fit function for second peak
Double_t fFixedRawYield
switch for fix Sigma of second gaussian in case of k2Gaus
Double_t FitFunction4Bkg(Double_t *x, Double_t *par)
Double_t fSigmaSgn
unc on signal gaussian mean value
Double_t fSecMass
fit parameters in 2nd peak fit function
TH1F * GetOverBackgroundResidualsAndPulls(TH1 *hPulls=0x0, TH1 *hResidualTrend=0x0, TH1 *hPullsTrend=0x0, Double_t minrange=0, Double_t maxrange=-1)
Bool_t fFixSecMass
width of the 2nd peak
TF1 * fSecFunc
flag to fix the width of the 2nd peak
Double_t * sigma
int Int_t
Definition: External.C:63
Bool_t fFixedSigma2Gaus
switch for fix Sigma of gaussian
Double_t fSecWidth
position of the 2nd peak
TString fFitOption
number of sigmas to veto the signal peak
TF1 * fBkRFunc
fit function for reflections
Double_t fFrac2Gaus
initialization for wa yield
Int_t fNParsSec
switch off/on second peak (for D+->KKpi in Ds)
TF1 * CreateSignalFitFunction(TString fname, Double_t integral)
Bool_t fFixedMean
signal second gaussian sigma in case of k2Gaus
Double_t fMassParticle
help variable
Double_t FitFunction4Refl(Double_t *x, Double_t *par)
Double_t GetRawYieldBinCounting(Double_t &errRyBC, Double_t nSigma=3., Int_t option=0, Int_t pdgCode=0) const
TH1F * SetTemplateReflections(const TH1 *h, TString opt, Double_t minRange, Double_t maxRange)
Double_t nsigma
Int_t fTypeOfFit4Bkg
upper mass limit
TH1F * fHistoTemplRfl
switch for fix refl/signal
TF1 * fBkgFunc
background fit function (1st step, side bands only)
Double_t fMassErr
signal gaussian mean value
Double_t fRflOverSig
fit parameters in reflection fit function
Double_t fMinMass
histogram to fit
Bool_t fFixRflOverSig
reflection/signal
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
Double_t BackFitFuncPolHelper(Double_t *x, Double_t *par)
Bool_t draw[nPtBins]
Double_t fMaxMass
lower mass limit
TF1 * CreateBackgroundFitFunction(TString fname, Double_t integral)
TF1 * fRflFunc
internal variable for fit with reflections
Double_t fNSigma4SideBands
kTRUE = only side bands considered
void Background(Double_t nOfSigma, Double_t &background, Double_t &errbackground) const
Double_t minMass
Double_t fRatio2GausSigma
switch for fixed fraction of 2nd gaussian in case of k2Gaus or k2GausSigmaRatioPar ...
Double_t FitFunction4Sgn(Double_t *x, Double_t *par)
Int_t fPolDegreeBkg
background fit func
Double_t fSigmaSgn2Gaus
unc on signal gaussian sigma
Double_t maxMass
bool Bool_t
Definition: External.C:53
TF1 * fBkgFuncRefit
background fit function (1st step, extended in peak region)
Double_t CheckForSignal(Double_t mean, Double_t sigma)
TF1 * fSigFunc
err on signal gaussian integral
TF1 * fBkgFuncSb
Signal fit function.
Bool_t fSecondPeak
fit function for reflections
Bool_t fSmoothRfl
histogram with reflection template
Definition: External.C:196
TF1 * CreateBackgroundPlusReflectionFunction(TString fname)
Double_t FitFunction4BkgAndRefl(Double_t *x, Double_t *par)
Int_t fNParsSig
switch for fixed ratio between two gaussian sigmas in case of k2GausSigmaRatioPar ...
Bool_t fReflections
background fit function (2nd step)
Bool_t fFixedSigma
switch for fix mean of gaussian