AliPhysics  ec7afe5 (ec7afe5)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
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  fFixedMean(kFALSE),
64  fFixedSigma(kFALSE),
65  fFixedRawYield(-1.),
66  fNParsSig(3),
67  fNParsBkg(2),
68  fOnlySideBands(kFALSE),
69  fNSigma4SideBands(4.),
70  fFitOption("L,E"),
71  fRawYield(0.),
72  fRawYieldErr(0.),
73  fSigFunc(0x0),
74  fBkgFuncSb(0x0),
75  fBkgFunc(0x0),
76  fBkgFuncRefit(0x0),
77  fReflections(kFALSE),
78  fNParsRfl(0),
79  fRflOverSig(0),
80  fFixRflOverSig(kFALSE),
81  fHistoTemplRfl(0x0),
82  fSmoothRfl(kFALSE),
83  fRawYieldHelp(0),
84  fRflFunc(0x0),
85  fBkRFunc(0x0),
86  fSecondPeak(kFALSE),
87  fNParsSec(0),
88  fSecMass(-999.),
89  fSecWidth(9999.),
90  fFixSecMass(kFALSE),
91  fFixSecWidth(kFALSE),
92  fSecFunc(0x0),
93  fTotFunc(0x0)
94 {
96 }
97 
98 //__________________________________________________________________________
99 AliHFInvMassFitter::AliHFInvMassFitter(const TH1F *histoToFit, Double_t minvalue, Double_t maxvalue, Int_t fittypeb, Int_t fittypes):
100  TNamed(),
101  fHistoInvMass(0x0),
102  fMinMass(minvalue),
103  fMaxMass(maxvalue),
104  fTypeOfFit4Bkg(fittypeb),
105  fPolDegreeBkg(4),
106  fCurPolDegreeBkg(-1),
107  fMassParticle(1.864),
108  fTypeOfFit4Sgn(fittypes),
109  fMass(1.865),
110  fMassErr(0.),
111  fSigmaSgn(0.012),
112  fSigmaSgnErr(0.),
113  fFixedMean(kFALSE),
114  fFixedSigma(kFALSE),
115  fFixedRawYield(-1.),
116  fNParsSig(3),
117  fNParsBkg(2),
118  fOnlySideBands(kFALSE),
119  fNSigma4SideBands(4.),
120  fFitOption("L,E"),
121  fRawYield(0.),
122  fRawYieldErr(0.),
123  fSigFunc(0x0),
124  fBkgFuncSb(0x0),
125  fBkgFunc(0x0),
126  fBkgFuncRefit(0x0),
127  fReflections(kFALSE),
128  fNParsRfl(0),
129  fRflOverSig(0),
130  fFixRflOverSig(kFALSE),
131  fHistoTemplRfl(0x0),
132  fSmoothRfl(kFALSE),
133  fRawYieldHelp(0),
134  fRflFunc(0x0),
135  fBkRFunc(0x0),
136  fSecondPeak(kFALSE),
137  fNParsSec(0),
138  fSecMass(-999.),
139  fSecWidth(9999.),
140  fFixSecMass(kFALSE),
141  fFixSecWidth(kFALSE),
142  fSecFunc(0x0),
143  fTotFunc(0x0)
144 {
146  fHistoInvMass=(TH1F*)histoToFit->Clone("fHistoInvMass");
147  fHistoInvMass->SetDirectory(0);
149 }
150 //_________________________________________________________________________
152 
154 
155  delete fSigFunc;
156  delete fBkgFunc;
157  delete fBkgFuncSb;
158  delete fBkgFuncRefit;
159  delete fHistoInvMass;
160  delete fRflFunc;
161  delete fBkRFunc;
162  delete fSecFunc;
163  delete fTotFunc;
164  delete fHistoTemplRfl;
165 
166 }
167 //__________________________________________________________________________
171 
172  switch (fTypeOfFit4Bkg) {
173  case 0:
174  fNParsBkg=2;
175  break;
176  case 1:
177  fNParsBkg=2;
178  break;
179  case 2:
180  fNParsBkg=3;
181  break;
182  case 3:
183  fNParsBkg=1;
184  break;
185  case 4:
186  fNParsBkg=2;
187  break;
188  case 5:
189  fNParsBkg=3;
190  break;
191  case 6:
193  break;
194  default:
195  AliError("Error in computing fNParsBkg: check fTypeOfFit4Bkg");
196  break;
197  }
198 
199  switch (fTypeOfFit4Sgn) {
200  case 0:
201  fNParsSig=3;
202  break;
203  case 1:
204  fNParsSig=5;
205  break;
206  default:
207  AliError("Error in computing fNParsSig: check fTypeOfFit4Sgn");
208  break;
209  }
210 
211  if(fReflections) fNParsRfl=1;
212  else fNParsRfl=0;
213 
214  if(fSecondPeak) fNParsSec=3;
215  else fNParsSec=0;
216 
217 }
218 //__________________________________________________________________________
222 
223  TVirtualFitter::SetDefaultFitter("Minuit");
224 
225  Double_t integralHisto=fHistoInvMass->Integral(fHistoInvMass->FindBin(fMinMass),fHistoInvMass->FindBin(fMaxMass),"width");
226 
227  fOnlySideBands = kTRUE;
228  fBkgFuncSb = CreateBackgroundFitFunction("funcbkgsb",integralHisto);
229  Int_t status=-1;
230  printf("\n--- First fit with only background on the side bands - Exclusion region = %.2f sigma ---\n",fNSigma4SideBands);
231  if(fTypeOfFit4Bkg==6){
233  fHistoInvMass->GetListOfFunctions()->Add(fBkgFuncSb);
234  fHistoInvMass->GetFunction(fBkgFuncSb->GetName())->SetBit(1<<9,kTRUE);
235  status=0;
236  }
237  }
238  else status=fHistoInvMass->Fit("funcbkgsb",Form("R,%s,+,0",fFitOption.Data()));
239  fBkgFuncSb->SetLineColor(kGray+1);
240  if (status != 0){
241  printf(" ---> Failed first fit with only background, minuit status = %d\n",status);
242  return kFALSE;
243  }
244 
245  fOnlySideBands = kFALSE;
246  if(!fBkgFunc){
247  fBkgFunc = CreateBackgroundFitFunction("funcbkg",integralHisto);
248  for(Int_t ipar=0; ipar<fNParsBkg; ipar++) fBkgFunc->SetParameter(ipar,fBkgFuncSb->GetParameter(ipar));
249  }
250  fBkgFunc->SetLineColor(kGray+1);
251 
252  printf("\n--- Estimate signal counts in the peak region ---\n");
253  Double_t estimSignal=CheckForSignal(fMass,fSigmaSgn);
254  if(estimSignal<0.){
255  if(draw) DrawFit();
256  return kTRUE;
257  }
258 
259  printf("\n--- Final fit with signal+background on the full range ---\n");
260  fRawYieldHelp=estimSignal; // needed for reflection normalization
261  if(!fBkgFuncRefit){
262  fBkgFuncRefit = CreateBackgroundFitFunction("funcbkgrefit",integralHisto);
263  for(Int_t ipar=0; ipar<fNParsBkg; ipar++) fBkgFuncRefit->SetParameter(ipar,fBkgFunc->GetParameter(ipar));
264  }
265  fBkgFuncRefit->SetLineColor(2);
266  fSigFunc = CreateSignalFitFunction("fsigfit",estimSignal);
267  if(fSecondPeak){
268  printf(" ---> Final fit includes a second inv. mass peak\n");
270  fSecFunc = CreateSecondPeakFunction("fsecpeak",estimSec);
271  }
272  if(fReflections){
273  printf(" ---> Final fit includes reflections\n");
274  fRflFunc = CreateReflectionFunction("freflect");
276  }
277  fTotFunc = CreateTotalFitFunction("funcmass");
278  status=fHistoInvMass->Fit("funcmass",Form("R,%s,+,0",fFitOption.Data()));
279  if (status != 0){
280  printf(" ---> Failed fit with signal+background, minuit status = %d\n",status);
281  return kFALSE;
282  }
283  for(Int_t ipar=0; ipar<fNParsBkg; ipar++){
284  fBkgFuncRefit->SetParameter(ipar,fTotFunc->GetParameter(ipar));
285  fBkgFuncRefit->SetParError(ipar,fTotFunc->GetParError(ipar));
286  }
287  for(Int_t ipar=0; ipar<fNParsSig; ipar++){
288  fSigFunc->SetParameter(ipar,fTotFunc->GetParameter(ipar+fNParsBkg));
289  fSigFunc->SetParError(ipar,fTotFunc->GetParError(ipar+fNParsBkg));
290  }
291  if(fSecondPeak){
292  for(Int_t ipar=0; ipar<fNParsSec; ipar++){
293  fSecFunc->SetParameter(ipar,fTotFunc->GetParameter(ipar+fNParsBkg+fNParsSig));
294  fSecFunc->SetParError(ipar,fTotFunc->GetParError(ipar+fNParsBkg+fNParsSig));
295  }
296  fSecFunc->SetLineColor(kMagenta+1);
297  fSecFunc->SetLineStyle(3);
298  }
299  if(fReflections){
300  for(Int_t ipar=0; ipar<fNParsRfl; ipar++){
301  fRflFunc->SetParameter(ipar,fTotFunc->GetParameter(ipar+fNParsBkg+fNParsSig+fNParsSec));
302  fRflFunc->SetParError(ipar,fTotFunc->GetParError(ipar+fNParsBkg+fNParsSig+fNParsSec));
303  fBkRFunc->SetParameter(ipar+fNParsBkg,fTotFunc->GetParameter(ipar+fNParsBkg+fNParsSig+fNParsSec));
304  fBkRFunc->SetParError(ipar+fNParsBkg,fTotFunc->GetParError(ipar+fNParsBkg+fNParsSig+fNParsSec));
305  }
306  for(Int_t ipar=0; ipar<fNParsBkg; ipar++){
307  fBkRFunc->SetParameter(ipar,fTotFunc->GetParameter(ipar));
308  fBkRFunc->SetParError(ipar,fTotFunc->GetParError(ipar));
309  }
310  fRflFunc->SetLineColor(kGreen+1);
311  fBkRFunc->SetLineColor(kRed+1);
312  fBkRFunc->SetLineStyle(7);
313  }
314  fMass=fSigFunc->GetParameter(1);
315  fMassErr=fSigFunc->GetParError(1);
316  fSigmaSgn=fSigFunc->GetParameter(2);
317  fSigmaSgnErr=fSigFunc->GetParError(2);
318  fTotFunc->SetLineColor(4);
319  fRawYield=fTotFunc->GetParameter(fNParsBkg)/fHistoInvMass->GetBinWidth(1);
320  fRawYieldErr=fTotFunc->GetParError(fNParsBkg)/fHistoInvMass->GetBinWidth(1);
322  if(draw) DrawFit();
323  return kTRUE;
324 }
325 //______________________________________________________________________________
329 
330  TCanvas* c0=new TCanvas("c0");
331  DrawHere(c0);
332 }
333 //______________________________________________________________________________
334 void AliHFInvMassFitter::DrawHere(TVirtualPad* c, Double_t nsigma,Int_t writeFitInfo){
337 
338  gStyle->SetOptStat(0);
339  gStyle->SetCanvasColor(0);
340  gStyle->SetFrameFillColor(0);
341  c->cd();
342  fHistoInvMass->GetXaxis()->SetRangeUser(fMinMass,fMaxMass);
343  fHistoInvMass->SetMarkerStyle(20);
344  fHistoInvMass->SetMinimum(0.);
345  fHistoInvMass->Draw("PE");
346  if(fBkgFunc) fBkgFunc->Draw("same");
347  if(fBkgFuncRefit) fBkgFuncRefit->Draw("same");
348  if(fBkRFunc) fBkRFunc->Draw("same");
349  if(fRflFunc) fRflFunc->Draw("same");
350  if(fSecFunc) fSecFunc->Draw("same");
351  if(fTotFunc) fTotFunc->Draw("same");
352 
353  if(writeFitInfo > 0){
354  TPaveText *pinfos=new TPaveText(0.12,0.65,0.47,0.89,"NDC");
355  TPaveText *pinfom=new TPaveText(0.6,0.7,1.,.87,"NDC");
356  pinfos->SetBorderSize(0);
357  pinfos->SetFillStyle(0);
358  pinfom->SetBorderSize(0);
359  pinfom->SetFillStyle(0);
360  if(fTotFunc){
361  pinfom->SetTextColor(kBlue);
362  for(Int_t ipar=1; ipar<fNParsSig; ipar++){
363  pinfom->AddText(Form("%s = %.3f #pm %.3f",fTotFunc->GetParName(ipar+fNParsBkg),fTotFunc->GetParameter(ipar+fNParsBkg),fTotFunc->GetParError(ipar+fNParsBkg)));
364  }
365  pinfom->Draw();
366 
367  Double_t bkg,errbkg;
368  Background(nsigma,bkg,errbkg);
369  Double_t signif,errsignif;
370  Significance(nsigma,signif,errsignif);
371 
372  pinfos->AddText(Form("S = %.0f #pm %.0f ",fRawYield,fRawYieldErr));
373  pinfos->AddText(Form("B (%.0f#sigma) = %.0f #pm %.0f",nsigma,bkg,errbkg));
374  pinfos->AddText(Form("S/B = (%.0f#sigma) %.4f ",nsigma,fRawYield/bkg));
375  if(fRflFunc) pinfos->AddText(Form("Refl/Sig = %.3f #pm %.3f ",fRflFunc->GetParameter(0),fRflFunc->GetParError(0)));
376  pinfos->AddText(Form("Signif (%.0f#sigma) = %.1f #pm %.1f ",nsigma,signif,errsignif));
377  pinfos->Draw();
378  }
379  }
380  c->Update();
381  return;
382 }
383 //______________________________________________________________________________
387 
388  Double_t minForSig=mean-4.*sigma;
389  Double_t maxForSig=mean+4.*sigma;
390  Int_t binForMinSig=fHistoInvMass->FindBin(minForSig);
391  Int_t binForMaxSig=fHistoInvMass->FindBin(maxForSig);
392  Double_t sum=0.;
393  Double_t sumback=0.;
394  fBkgFunc->Print();
395  for(Int_t ibin=binForMinSig; ibin<=binForMaxSig; ibin++){
396  sum+=fHistoInvMass->GetBinContent(ibin);
397  sumback+=fBkgFunc->Eval(fHistoInvMass->GetBinCenter(ibin));
398  }
399  Double_t diffUnderPeak=(sum-sumback);
400  printf(" ---> IntegralUnderFitFunc=%f IntegralUnderHisto=%f EstimatedSignal=%f\n",sum,sumback,diffUnderPeak);
401  if(diffUnderPeak/TMath::Sqrt(sum)<1.){
402  printf(" ---> (Tot-Bkg)/sqrt(Tot)=%f ---> Likely no signal/\n",diffUnderPeak/TMath::Sqrt(sum));
403  return -1;
404  }
405  return diffUnderPeak*fHistoInvMass->GetBinWidth(1);
406 }
407 
408 //______________________________________________________________________________
412 
414  TF1* funcbkg = new TF1(fname.Data(),this,&AliHFInvMassFitter::FitFunction4Bkg,fMinMass,fMaxMass,fNParsBkg,"AliHFInvMassFitter","FitFunction4Bkg");
415  switch (fTypeOfFit4Bkg) {
416  case 0: //gaus+expo
417  funcbkg->SetParNames("BkgInt","Slope");
418  funcbkg->SetParameters(integral,-2.);
419  break;
420  case 1:
421  funcbkg->SetParNames("BkgInt","Slope");
422  funcbkg->SetParameters(integral,-100.);
423  break;
424  case 2:
425  funcbkg->SetParNames("BkgInt","Coef1","Coef2");
426  funcbkg->SetParameters(integral,-10.,5);
427  break;
428  case 3:
429  funcbkg->SetParNames("Const");
430  funcbkg->SetParameter(0,0.);
431  funcbkg->FixParameter(0,0.);
432  break;
433  case 4:
434  funcbkg->SetParNames("BkgInt","Coef1");
435  funcbkg->SetParameters(integral,0.5);
436  break;
437  case 5:
438  funcbkg->SetParNames("BkgInt","Coef1","Coef2");
439  funcbkg->SetParameters(integral,-10.,5.);
440  break;
441  case 6:
442  for(Int_t j=0;j<fNParsBkg; j++){
443  funcbkg->SetParName(j,Form("Coef%d",j));
444  funcbkg->SetParameter(j,0);
445  }
446  funcbkg->SetParameter(0,integral);
447  break;
448  default:
449  AliError(Form("Wrong choice of fTypeOfFit4Bkg (%d)",fTypeOfFit4Bkg));
450  delete funcbkg;
451  return 0x0;
452  break;
453  }
454  funcbkg->SetLineColor(kBlue+3);
455  return funcbkg;
456 }
457 //______________________________________________________________________________
461 
462  TF1* funcsec = new TF1(fname.Data(),this,&AliHFInvMassFitter::FitFunction4SecPeak,fMinMass,fMaxMass,3,"AliHFInvMassFitter","FitFunction4SecPeak");
463  funcsec->SetParameter(0,integsig);
464  funcsec->SetParameter(1,fSecMass);
465  if(fFixSecMass) funcsec->FixParameter(1,fSecMass);
466  funcsec->SetParameter(2,fSecWidth);
467  if(fFixSecWidth) funcsec->FixParameter(2,fSecWidth);
468  funcsec->SetParNames("SecPeakInt","SecMean","SecSigma");
469  return funcsec;
470 }
471 //______________________________________________________________________________
475  TF1* funcrfl = new TF1(fname.Data(),this,&AliHFInvMassFitter::FitFunction4Refl,fMinMass,fMaxMass,1,"AliHFInvMassFitter","FitFunction4Refl");
476  funcrfl->SetParameter(0,fRflOverSig);
477  funcrfl->SetParLimits(0,0.,1.);
478  if(fFixRflOverSig) funcrfl->FixParameter(0,fRflOverSig);
479  funcrfl->SetParNames("ReflOverS");
480  return funcrfl;
481 }
482 //______________________________________________________________________________
487  Int_t totParams=fNParsBkg+fNParsRfl;
488  TF1* fbr=new TF1(fname.Data(),this,&AliHFInvMassFitter::FitFunction4BkgAndRefl,fMinMass,fMaxMass,totParams,"AliHFInvMassFitter","FitFunction4BkgAndRefl");
489  for(Int_t ipar=0; ipar<fNParsBkg; ipar++){
490  fbr->SetParameter(ipar,fBkgFunc->GetParameter(ipar));
491  fbr->SetParName(ipar,fBkgFunc->GetParName(ipar));
492  }
493  for(Int_t ipar=0; ipar<fNParsRfl; ipar++){
494  fbr->SetParameter(ipar+fNParsBkg,fRflFunc->GetParameter(ipar));
495  fbr->SetParName(ipar+fNParsBkg,fRflFunc->GetParName(ipar));
496  // par limits not set because this function is not used for fitting but only for drawing
497  }
498  return fbr;
499 }
500 //______________________________________________________________________________
504 
506  TF1* funcsig = new TF1(fname.Data(),this,&AliHFInvMassFitter::FitFunction4Sgn,fMinMass,fMaxMass,fNParsSig,"AliHFInvMassFitter","FitFunction4Sgn");
507  if(fTypeOfFit4Sgn==kGaus){
508  funcsig->SetParameter(0,integsig);
509  if(fFixedRawYield>-0.1) funcsig->FixParameter(0,fFixedRawYield);
510  funcsig->SetParameter(1,fMass);
511  if(fFixedMean) funcsig->FixParameter(1,fMass);
512  funcsig->SetParameter(2,fSigmaSgn);
513  if(fFixedSigma) funcsig->FixParameter(2,fSigmaSgn);
514  funcsig->SetParNames("SgnInt","Mean","Sigma");
515  }
516  if(fTypeOfFit4Sgn==k2Gaus){
517  funcsig->SetParameter(0,integsig);
518  if(fFixedRawYield>-0.1) funcsig->FixParameter(0,fFixedRawYield);
519  funcsig->SetParameter(1,fMass);
520  if(fFixedMean) funcsig->FixParameter(1,fMass);
521  funcsig->SetParameter(2,fSigmaSgn);
522  funcsig->SetParLimits(2,0.004,0.05);
523  if(fFixedSigma) funcsig->FixParameter(2,fSigmaSgn);
524  funcsig->SetParameter(3,0.2);
525  funcsig->SetParLimits(3,0.,1.);
526  funcsig->SetParameter(4,fSigmaSgn);
527  funcsig->SetParLimits(4,0.004,0.05);
528  funcsig->SetParNames("SgnInt","Mean","Sigma1","Frac","Sigma2");
529  }
530  return funcsig;
531 }
532 
533 //______________________________________________________________________________
537 
540  TF1* ftot=ftot=new TF1(fname.Data(),this,&AliHFInvMassFitter::FitFunction4Mass,fMinMass,fMaxMass,totParams,"AliHFInvMassFitter","FitFunction4Mass");
541  for(Int_t ipar=0; ipar<fNParsBkg; ipar++){
542  ftot->SetParameter(ipar,fBkgFunc->GetParameter(ipar));
543  ftot->SetParName(ipar,fBkgFunc->GetParName(ipar));
544  }
545  for(Int_t ipar=0; ipar<fNParsSig; ipar++){
546  ftot->SetParameter(ipar+fNParsBkg,fSigFunc->GetParameter(ipar));
547  ftot->SetParName(ipar+fNParsBkg,fSigFunc->GetParName(ipar));
548  Double_t parmin,parmax;
549  fSigFunc->GetParLimits(ipar,parmin,parmax);
550  ftot->SetParLimits(ipar+fNParsBkg,parmin,parmax);
551  }
552  if(fSecondPeak && fSecFunc){
553  for(Int_t ipar=0; ipar<fNParsSec; ipar++){
554  ftot->SetParameter(ipar+fNParsBkg+fNParsSig,fSecFunc->GetParameter(ipar));
555  ftot->SetParName(ipar+fNParsBkg+fNParsSig,fSecFunc->GetParName(ipar));
556  Double_t parmin,parmax;
557  fSecFunc->GetParLimits(ipar,parmin,parmax);
558  ftot->SetParLimits(ipar+fNParsBkg+fNParsSig,parmin,parmax);
559  }
560  }
561  if(fReflections && fRflFunc){
562  for(Int_t ipar=0; ipar<fNParsRfl; ipar++){
563  ftot->SetParameter(ipar+fNParsBkg+fNParsSig+fNParsSec,fRflFunc->GetParameter(ipar));
564  ftot->SetParName(ipar+fNParsBkg+fNParsSig+fNParsSec,fRflFunc->GetParName(ipar));
565  Double_t parmin,parmax;
566  fRflFunc->GetParLimits(ipar,parmin,parmax);
567  ftot->SetParLimits(ipar+fNParsBkg+fNParsSig+fNParsSec,parmin,parmax);
568  }
569  }
570  return ftot;
571 }
572 //__________________________________________________________________________
576 
578  if(fOnlySideBands && TMath::Abs(x[0]-fMass) < maxDeltaM) {
579  TF1::RejectPoint();
580  return 0;
581  }
582  if(fOnlySideBands && fSecondPeak && TMath::Abs(x[0]-fSecMass) < (fNSigma4SideBands*fSecWidth)){
583  TF1::RejectPoint();
584  return 0;
585  }
586  Double_t total=0;
587 
588  switch (fTypeOfFit4Bkg){
589  case 0:
590  //exponential
591  //exponential = A*exp(B*x) -> integral(exponential)=A/B*exp(B*x)](min,max)
592  //-> A = B*integral/(exp(B*max)-exp(B*min)) where integral can be written
593  //as integralTot- integralGaus (=par [2])
594  //Par:
595  // * [0] = integralBkg;
596  // * [1] = B;
597  //exponential = [1]*[0]/(exp([1]*max)-exp([1]*min))*exp([1]*x)
598  total = par[0]*par[1]/(TMath::Exp(par[1]*fMaxMass)-TMath::Exp(par[1]*fMinMass))*TMath::Exp(par[1]*x[0]);
599  // AliInfo("Background function set to: exponential");
600  break;
601  case 1:
602  //linear
603  //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)
604  // * [0] = integralBkg;
605  // * [1] = b;
606  total= par[0]/(fMaxMass-fMinMass)+par[1]*(x[0]-0.5*(fMaxMass+fMinMass));
607  // AliInfo("Background function set to: linear");
608  break;
609  case 2:
610  //parabola
611  //y=a+b*x+c*x**2 -> integral = a(max-min) + 1/2*b*(max^2-min^2) +
612  //+ 1/3*c*(max^3-min^3) ->
613  //a = (integral-1/2*b*(max^2-min^2)-1/3*c*(max^3-min^3))/(max-min)
614  // * [0] = integralBkg;
615  // * [1] = b;
616  // * [2] = c;
617  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));
618  // AliInfo("Background function set to: polynomial");
619  break;
620  case 3:
621  total=par[0];
622  break;
623  case 4:
624  //power function
625  //y=a(x-m_pi)^b -> integral = a/(b+1)*((max-m_pi)^(b+1)-(min-m_pi)^(b+1))
626  //
627  //a = integral*(b+1)/((max-m_pi)^(b+1)-(min-m_pi)^(b+1))
628  // * [0] = integralBkg;
629  // * [1] = b;
630  // a(power function) = [0]*([1]+1)/((max-m_pi)^([1]+1)-(min-m_pi)^([1]+1))*(x-m_pi)^[1]
631  {
632  Double_t mpi = TDatabasePDG::Instance()->GetParticle(211)->Mass();
633 
634  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]);
635  // AliInfo("Background function set to: powerlaw");
636  }
637  break;
638  case 5:
639  //power function wit exponential
640  //y=a*Sqrt(x-m_pi)*exp(-b*(x-m_pi))
641  {
642  Double_t mpi = TDatabasePDG::Instance()->GetParticle(211)->Mass();
643 
644  total = par[1]*TMath::Sqrt(x[0] - mpi)*TMath::Exp(-1.*par[2]*(x[0]-mpi));
645  // AliInfo("Background function set to: wit exponential");
646  }
647  break;
648  case 6:
649  // the following comment must be removed
650  // // pol 3, following convention for pol 2
651  // //y=a+b*x+c*x**2+d*x**3 -> integral = a(max-min) + 1/2*b*(max^2-min^2) +
652  // //+ 1/3*c*(max^3-min^3) + 1/4 d * (max^4-min^4) ->
653  // //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)
654  // // * [0] = integralBkg;
655  // // * [1] = b;
656  // // * [2] = c;
657  // // * [3] = d;
658  {
659  total=par[0];
660  for(Int_t it=1;it<=fPolDegreeBkg;it++){
661  total+=par[it]*TMath::Power(x[0]-fMassParticle,it)/TMath::Factorial(it);
662  }
663  }
664  break;
665  }
666  return total;
667 }
668 //_________________________________________________________________________
672 
673  // AliInfo("Signal function set to: Gaussian");
674  Double_t sigval=0;
675  switch (fTypeOfFit4Sgn){
676  case 0:
677  //gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
678  //Par:
679  // * [0] = integralSgn
680  // * [1] = mean
681  // * [2] = sigma
682  //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
683  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]);
684  break;
685  case 1:
686  //double gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
687  //Par:
688  // * [0] = integralSgn
689  // * [1] = mean
690  // * [2] = sigma1
691  // * [3] = 2nd gaussian ratio
692  // * [4] = deltaSigma
693  //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
694  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]);
695  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]);
696  sigval=par[0]*(g1+g2);
697  break;
698  }
699  fRawYieldHelp=par[0]/fHistoInvMass->GetBinWidth(1);
700  return sigval;
701 }
702 //_________________________________________________________________________
706  if(!fHistoTemplRfl) return 0;
707 
708  Int_t bin =fHistoTemplRfl->FindBin(x[0]);
709  Double_t value=fHistoTemplRfl->GetBinContent(bin);
710  Int_t binmin=fHistoTemplRfl->FindBin(fMinMass*1.00001);
711  Int_t binmax=fHistoTemplRfl->FindBin(fMaxMass*0.99999);
712  Double_t norm=fHistoTemplRfl->Integral(binmin,binmax)*fHistoTemplRfl->GetBinWidth(bin);
713  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
714  value+=fHistoTemplRfl->GetBinContent(bin-1)+fHistoTemplRfl->GetBinContent(bin+1);
715  value/=3.;
716  }
717  return par[0]*value/norm*fRawYieldHelp*fHistoInvMass->GetBinWidth(1);
718 }
719 //_________________________________________________________________________
723  Double_t bkg=FitFunction4Bkg(x,par);
724  Double_t refl=0;
725  if(fReflections) refl=FitFunction4Refl(x,&par[fNParsBkg]);
726  return bkg+refl;
727 }
728 //_________________________________________________________________________
732 
733  //gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
734  //Par:
735  // * [0] = integralSgn
736  // * [1] = mean
737  // * [2] = sigma
738  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]);
739  return secgaval;
740 }
741 //_________________________________________________________________________
745 
746  Double_t bkg=FitFunction4Bkg(x,par);
747  Double_t sig=FitFunction4Sgn(x,&par[fNParsBkg]);
748  Double_t sec=0.;
749  if(fSecondPeak) sec=FitFunction4SecPeak(x,&par[fNParsBkg+fNParsSig]);
750  Double_t refl=0;
751  if(fReflections) refl=FitFunction4Refl(x,&par[fNParsBkg+fNParsSig+fNParsSec]);
752  return bkg+sig+sec+refl;
753 }
754 
755 //_________________________________________________________________________
756 void AliHFInvMassFitter::Signal(Double_t nOfSigma,Double_t &signal,Double_t &errsignal) const {
759 
760  Double_t minMass=fMass-nOfSigma*fSigmaSgn;
761  Double_t maxMass=fMass+nOfSigma*fSigmaSgn;
762  Signal(minMass,maxMass,signal,errsignal);
763  return;
764 }
765 
766 //_________________________________________________________________________
767 void AliHFInvMassFitter::Signal(Double_t min, Double_t max, Double_t &signal,Double_t &errsignal) const {
770 
771  signal=fSigFunc->Integral(min, max)/(Double_t)fHistoInvMass->GetBinWidth(1);
772  errsignal=(fRawYieldErr/fRawYield)*signal;/*assume relative error is the same as for total integral*/
773  return;
774 }
775 
776 //___________________________________________________________________________
777 void AliHFInvMassFitter::Background(Double_t nOfSigma,Double_t &background,Double_t &errbackground) const {
780 
781  Double_t minMass=fMass-nOfSigma*fSigmaSgn;
782  Double_t maxMass=fMass+nOfSigma*fSigmaSgn;
783  Background(minMass,maxMass,background,errbackground);
784  return;
785 
786 }
787 //___________________________________________________________________________
788 void AliHFInvMassFitter::Background(Double_t min, Double_t max, Double_t &background,Double_t &errbackground) const {
791 
792  TF1 *funcbkg=0x0;
793  if(fBkgFuncRefit) funcbkg=fBkgFuncRefit;
794  else if(fBkgFunc) funcbkg=fBkgFunc;
795  if(!funcbkg){
796  AliError("Bkg function not found!");
797  return;
798  }
799 
800  Double_t intB=funcbkg->GetParameter(0);
801  Double_t intBerr=funcbkg->GetParError(0);
802 
803  //relative error evaluation: from histo
804 
806  Int_t rightBand=fHistoInvMass->FindBin(fMass+fNSigma4SideBands*fSigmaSgn);
807  intB=fHistoInvMass->Integral(1,leftBand)+fHistoInvMass->Integral(rightBand,fHistoInvMass->GetNbinsX());
808  Double_t sum2=0;
809  for(Int_t i=1;i<=leftBand;i++){
810  sum2+=fHistoInvMass->GetBinError(i)*fHistoInvMass->GetBinError(i);
811  }
812  for(Int_t i=rightBand; i<=fHistoInvMass->GetNbinsX();i++){
813  sum2+=fHistoInvMass->GetBinError(i)*fHistoInvMass->GetBinError(i);
814  }
815 
816  intBerr=TMath::Sqrt(sum2);
817 
818  background=funcbkg->Integral(min,max)/(Double_t)fHistoInvMass->GetBinWidth(1);
819  errbackground=intBerr/intB*background;
820 
821  return;
822 
823 }
824 //__________________________________________________________________________
825 
826 void AliHFInvMassFitter::Significance(Double_t nOfSigma,Double_t &significance,Double_t &errsignificance) const {
829 
830  Double_t minMass=fMass-nOfSigma*fSigmaSgn;
831  Double_t maxMass=fMass+nOfSigma*fSigmaSgn;
832  Significance(minMass, maxMass, significance, errsignificance);
833 
834  return;
835 }
836 
837 //__________________________________________________________________________
838 
839 void AliHFInvMassFitter::Significance(Double_t min, Double_t max, Double_t &significance,Double_t &errsignificance) const {
842 
843  Double_t background,errbackground;
844  Background(min,max,background,errbackground);
845 
846  if (fRawYield+background <= 0.){
847  significance=-1;
848  errsignificance=0;
849  return;
850  }
851 
852  AliVertexingHFUtils::ComputeSignificance(fRawYield,fRawYieldErr,background,errbackground,significance,errsignificance);
853 
854  return;
855 }
856 //________________________________________________________________________
860 
861  TH1F *hCp=(TH1F*)fHistoInvMass->Clone("htemp");
862  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
863  Double_t estimateslope=(hCp->GetBinContent(hCp->FindBin(fMass+3.5*fSigmaSgn))-hCp->GetBinContent(hCp->FindBin(fMass-3.5*fSigmaSgn)))/(7*fSigmaSgn);// first rough estimate
864 
866  TF1 *funcbkg,*funcPrev=0x0;
868  funcbkg = new TF1(Form("temp%d",fCurPolDegreeBkg),this,&AliHFInvMassFitter::BackFitFuncPolHelper,fMinMass,fMaxMass,fCurPolDegreeBkg+1,"AliHFInvMassFitter","BackFitFuncPolHelper");
869  if(funcPrev){
870  for(Int_t j=0;j<fCurPolDegreeBkg;j++){// now is +1 degree w.r.t. previous fit funct
871  funcbkg->SetParameter(j,funcPrev->GetParameter(j));
872  }
873  delete funcPrev;
874  }
875  else{
876  funcbkg->SetParameter(0,estimatecent);
877  funcbkg->SetParameter(1,estimateslope);
878  }
879  printf(" ---> Pre-fit of background with pol degree %d ---\n",fCurPolDegreeBkg);
880  hCp->Fit(funcbkg,"REMN","");
881  funcPrev=(TF1*)funcbkg->Clone("ftemp");
882  delete funcbkg;
884  }
885 
886  for(Int_t j=0;j<=fPolDegreeBkg;j++){
887  fback->SetParameter(j,funcPrev->GetParameter(j));
888  fback->SetParError(j,funcPrev->GetParError(j));
889  }
890  printf(" ---> Final background fit with pol degree %d ---\n",fPolDegreeBkg);
891  hCp->Fit(fback,"REMN","");// THIS IS JUST TO SET NOT ONLY THE PARAMETERS BUT ALSO chi2, etc...
892 
893 
894  // The following lines might be useful for debugging
895  // TCanvas *cDebug=new TCanvas();
896  // cDebug->cd();
897  // hCp->Draw();
898  // TString strout=Form("Test%d.root",(Int_t)fhistoInvMass->GetBinContent(fhistoInvMass->FindBin(fMass)));
899  // cDebug->Print(strout.Data());
900  // delete cDebug;
901 
902  delete funcPrev;
903  delete hCp;
904  return kTRUE;
905 
906 }
907 // _______________________________________________________________________
911 
913  if(fOnlySideBands && TMath::Abs(x[0]-fMass) < maxDeltaM) {
914  TF1::RejectPoint();
915  return 0;
916  }
917  Double_t back=par[0];
918  for(Int_t it=1;it<=fCurPolDegreeBkg;it++){
919  back+=par[it]*TMath::Power(x[0]-fMassParticle,it)/TMath::Factorial(it);
920  }
921  return back;
922 }
923 // _______________________________________________________________________
932 
933  if(!h){
934  fReflections=kFALSE;
935  return 0x0;
936  }
937  fHistoTemplRfl=(TH1F*)h->Clone("hTemplRfl");
938  opt.ToLower();
939  fReflections=kTRUE;
940  printf("\n--- Reflection templates from simulation ---\n");
941  if(opt.Contains("templ")){
942  printf(" ---> Reflection contribution using directly the histogram from simulation\n");
943  return fHistoTemplRfl;
944  }
945 
946  TF1 *f=0x0;
947  Bool_t isPoissErr=kTRUE;
948  Double_t xMinForFit=h->GetBinLowEdge(1);
949  Double_t xMaxForFit=h->GetXaxis()->GetBinUpEdge(h->GetNbinsX());
950  if(minRange>=0 && maxRange>=0){
951  xMinForFit=TMath::Max(minRange,h->GetBinLowEdge(1));
952  xMaxForFit=TMath::Min(maxRange,h->GetXaxis()->GetBinUpEdge(h->GetNbinsX()));
953  }
954  if(opt.EqualTo("1gaus") || opt.EqualTo("singlegaus")){
955  printf(" ---> Reflection contribution from single-Gaussian fit to histogram from simulation\n");
956  f=new TF1("mygaus","gaus",xMinForFit,xMaxForFit);
957  f->SetParameter(0,h->GetMaximum());
958  // f->SetParLimits(0,0,100.*h->Integral());
959  f->SetParameter(1,1.865);
960  f->SetParameter(2,0.050);
961  fHistoTemplRfl->Fit(f,"REM","");//,h->GetBinLowEdge(1),h->GetXaxis()->GetBinUpEdge(h->GetNbinsX()));
962  }
963  else if(opt.EqualTo("2gaus") || opt.EqualTo("doublegaus")){
964  printf(" ---> Reflection contribution from double-Gaussian fit to histogram from simulation\n");
965  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);
966  f->SetParameter(0,h->GetMaximum());
967  // f->SetParLimits(0,0,100.*h->Integral());
968  f->SetParLimits(3,0.,1.);
969  f->SetParameter(3,0.5);
970  f->SetParameter(1,1.84);
971  f->SetParameter(2,0.050);
972  f->SetParameter(4,1.88);
973  f->SetParameter(5,0.050);
974  fHistoTemplRfl->Fit(f,"REM","");//,h->GetBinLowEdge(1),h->GetXaxis()->GetBinUpEdge(h->GetNbinsX()));
975  }
976  else if(opt.EqualTo("pol3")){
977  printf(" ---> Reflection contribution from pol3 fit to histogram from simulation\n");
978  f=new TF1("mypol3","pol3",xMinForFit,xMaxForFit);
979  f->SetParameter(0,h->GetMaximum());
980  // f->SetParLimits(0,0,100.*h->Integral());
981  // Hard to initialize the other parameters...
982  fHistoTemplRfl->Fit(f,"REM","");
983  // Printf("We USED %d POINTS in the Fit",f->GetNumberFitPoints());
984  }
985  else if(opt.EqualTo("pol6")){
986  printf(" ---> Reflection contribution from pol6 fit to histogram from simulation\n");
987  f=new TF1("mypol6","pol6",xMinForFit,xMaxForFit);
988  f->SetParameter(0,h->GetMaximum());
989  // f->SetParLimits(0,0,100.*h->Integral());
990  // Hard to initialize the other parameters...
991  fHistoTemplRfl->Fit(f,"RLEMI","");//,h->GetBinLowEdge(1),h->GetXaxis()->GetBinUpEdge(h->GetNbinsX()));
992  }
993  else{
994  // no good option passed
995  printf(" ---> Bad option for reflection configuration -> reflections will not be included in the fit\n");
996  fReflections=kFALSE;
997  delete fHistoTemplRfl;
998  fHistoTemplRfl=0x0;
999  return 0x0;
1000  }
1001 
1002  // Fill fHistoTemplRfl with values of fit function
1003  if(f){
1004  for(Int_t j=1;j<=fHistoTemplRfl->GetNbinsX();j++){
1005  fHistoTemplRfl->SetBinContent(j,f->Integral(fHistoTemplRfl->GetBinLowEdge(j),fHistoTemplRfl->GetXaxis()->GetBinUpEdge(j))/fHistoTemplRfl->GetBinWidth(j));
1006  if(fHistoTemplRfl->GetBinContent(j)>=0.&&TMath::Abs(h->GetBinError(j)*h->GetBinError(j)-h->GetBinContent(j))>0.1*h->GetBinContent(j))isPoissErr=kFALSE;
1007  }
1008  for(Int_t j=1;j<=fHistoTemplRfl->GetNbinsX();j++){
1009  if(isPoissErr){
1010  fHistoTemplRfl->SetBinError(j,TMath::Sqrt(fHistoTemplRfl->GetBinContent(j)));
1011  }
1012  else fHistoTemplRfl->SetBinError(j,0.001*fHistoTemplRfl->GetBinContent(j));
1013  }
1014  fReflections=kTRUE;
1015  return fHistoTemplRfl;
1016  }else{
1017  printf(" ---> Fit to MC template for reflection failed -> reflections will not be included in the fit\n");
1018  fReflections=kFALSE;
1019  delete fHistoTemplRfl;
1020  fHistoTemplRfl=0x0;
1021  return 0x0;
1022  }
1023  return 0x0;
1024 }
1025 // _______________________________________________________________________
1030  // else (default) mean of gaussian fit
1031 
1032  Double_t massVal=fMass;
1033  switch (pdgCode) {
1034  case 411:
1035  case 421:
1036  case 431:
1037  case 4122:
1038  massVal=TDatabasePDG::Instance()->GetParticle(pdgCode)->Mass();
1039  break;
1040  case 413:
1041  massVal=TDatabasePDG::Instance()->GetParticle(pdgCode)->Mass();
1042  massVal-=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1043  break;
1044  default:
1045  massVal=fMass;
1046  break;
1047  }
1048 
1049  Double_t minMass=massVal-nOfSigma*fSigmaSgn;
1050  Double_t maxMass=massVal+nOfSigma*fSigmaSgn;
1051  return GetRawYieldBinCounting(errRyBC,minMass,maxMass,option);
1052 }
1053 // _______________________________________________________________________
1059 
1060  Int_t minBinSum=fHistoInvMass->FindBin(minMass);
1061  Int_t maxBinSum=fHistoInvMass->FindBin(maxMass);
1062  if(minBinSum<1){
1063  printf("Left range for bin counting smaller than allowed by histogram axis, setting it to the lower edge of the first histo bin\n");
1064  minBinSum=1;
1065  }
1066  if(maxBinSum>fHistoInvMass->GetNbinsX()){
1067  printf("Right range for bin counting larger than allowed by histogram axis, setting it to the upper edge of the last histo bin\n");
1068  maxBinSum=fHistoInvMass->GetNbinsX();
1069  }
1070  Double_t cntSig=0.;
1071  Double_t cntErr=0.;
1072  errRyBC=0;
1073  TF1* fbackground=fBkgFunc;
1074  if(option==1) fbackground=fBkgFuncRefit;
1075  if(!fbackground) return 0.;
1076 
1077  for(Int_t jb=minBinSum; jb<=maxBinSum; jb++){
1078  Double_t cntTot=fHistoInvMass->GetBinContent(jb);
1079  Double_t cntBkg=fbackground->Integral(fHistoInvMass->GetBinLowEdge(jb),fHistoInvMass->GetBinLowEdge(jb)+fHistoInvMass->GetBinWidth(jb))/fHistoInvMass->GetBinWidth(jb);
1080  //Double_t cntBkg=fbackground->Eval(fHistoInvMass->GetBinCenter(jb));
1081  cntSig+=(cntTot-cntBkg);
1082  cntErr+=(fHistoInvMass->GetBinError(jb)*fHistoInvMass->GetBinError(jb));
1083  }
1084  errRyBC=TMath::Sqrt(cntErr);
1085  return cntSig;
1086 }
1087 
1088 
1089 // _______________________________________________________________________
1090 TH1F* AliHFInvMassFitter::GetResidualsAndPulls(TH1 *hPulls,TH1 *hResidualTrend,TH1 *hPullsTrend, Double_t minrange,Double_t maxrange){
1091 
1093 
1094  Int_t binmi=fHistoInvMass->FindBin(fMinMass*1.001);
1095  Int_t binma=fHistoInvMass->FindBin(fMaxMass*0.9999);
1096  if(maxrange>minrange){
1097  binmi=fHistoInvMass->FindBin(minrange*1.001);
1098  binma=fHistoInvMass->FindBin(maxrange*0.9999);
1099  }
1100  if(hResidualTrend){
1101  //fHistoInvMass->Copy(hResidualTrend);
1102  hResidualTrend->SetBins(fHistoInvMass->GetNbinsX(),fHistoInvMass->GetXaxis()->GetXmin(),fHistoInvMass->GetXaxis()->GetXmax());
1103  hResidualTrend->SetName(Form("%s_residualTrend",fHistoInvMass->GetName()));
1104  hResidualTrend->SetTitle(Form("%s (Residuals)",fHistoInvMass->GetTitle()));
1105  hResidualTrend->SetMarkerStyle(20);
1106  hResidualTrend->SetMarkerSize(1.0);
1107  hResidualTrend->Reset();
1108  }
1109  if(hPullsTrend){
1110  hPullsTrend->SetBins(fHistoInvMass->GetNbinsX(),fHistoInvMass->GetXaxis()->GetXmin(),fHistoInvMass->GetXaxis()->GetXmax());
1111  hPullsTrend->Reset();
1112  hPullsTrend->SetName(Form("%s_pullTrend",fHistoInvMass->GetName()));
1113  hPullsTrend->SetTitle(Form("%s (Pulls)",fHistoInvMass->GetTitle()));
1114  hPullsTrend->SetMarkerStyle(20);
1115  hPullsTrend->SetMarkerSize(1.0);
1116  }
1117  if(hPulls){
1118  hPulls->SetName(Form("%s_pulls",fHistoInvMass->GetName()));
1119  hPulls->SetTitle(Form("%s ; Pulls",fHistoInvMass->GetTitle()));
1120  hPulls->SetBins(40,-10,10);
1121  hPulls->Reset();
1122  }
1123 
1124  Double_t res=-1.e-6,min=1.e+12,max=-1.e+12;
1125  TArrayD *arval=new TArrayD(binma-binmi+1);
1126  for(Int_t jst=1;jst<=fHistoInvMass->GetNbinsX();jst++){
1127 
1128  res=fHistoInvMass->GetBinContent(jst)-fTotFunc->Integral(fHistoInvMass->GetBinLowEdge(jst),fHistoInvMass->GetBinLowEdge(jst)+fHistoInvMass->GetBinWidth(jst))/fHistoInvMass->GetBinWidth(jst);
1129  if(jst>=binmi&&jst<=binma){
1130  arval->AddAt(res,jst-binmi);
1131  if(res<min)min=res;
1132  if(res>max)max=res;
1133  }
1134  // Printf("Res = %f from %f - %f",res,fHistoInvMass->GetBinContent(jst),fTotFunc->Integral(fHistoInvMass->GetBinLowEdge(jst),fHistoInvMass->GetBinLowEdge(jst)+fHistoInvMass->GetBinWidth(jst))/fHistoInvMass->GetBinWidth(jst));
1135  if(hResidualTrend){
1136  hResidualTrend->SetBinContent(jst,res);
1137  hResidualTrend->SetBinError(jst,fHistoInvMass->GetBinError(jst));
1138  }
1139  if(hPulls){
1140  if(jst>=binmi&&jst<=binma)hPulls->Fill(res/fHistoInvMass->GetBinError(jst));
1141  }
1142  if(hPullsTrend){
1143  hPullsTrend->SetBinContent(jst,res/fHistoInvMass->GetBinError(jst));
1144  hPullsTrend->SetBinError(jst,0.0001);
1145  }
1146  }
1147  if(hResidualTrend)hResidualTrend->GetXaxis()->SetRange(binmi,binma);
1148  if(hPullsTrend){
1149  hPullsTrend->GetXaxis()->SetRange(binmi,binma);
1150  hPullsTrend->SetMinimum(-7);
1151  hPullsTrend->SetMaximum(+7);
1152  }
1153  if(TMath::Abs(min)>TMath::Abs(max))max=min;
1154 
1155  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);
1156  for(Int_t j=0;j<binma-binmi+1;j++){
1157  hout->Fill(arval->At(j));
1158  }
1159  hout->Sumw2();
1160  hout->Fit("gaus","LEM","",-TMath::Abs(max)*1.2,TMath::Abs(max)*1.2);
1161 
1162  if(hPulls){
1163  hPulls->Sumw2();
1164  hPulls->Fit("gaus","LEM","",-3,3);
1165  }
1166  delete arval;
1167  return hout;
1168 }
1169 // _______________________________________________________________________
1173  if(fBkgFunc){
1174  printf("--- Background function in 1st step fit to side bands ---\n");
1175  fBkgFunc->Print();
1176  }
1177  if(fTotFunc){
1178  printf("--- Total fit function from 2nd step fit ---\n");
1179  fTotFunc->Print();
1180  }
1181  if(fBkgFuncRefit){
1182  printf("--- Background function in 2nd step fit ---\n");
1183  fBkgFuncRefit->Print();
1184  }
1185  if(fRflFunc){
1186  printf("--- Reflections ---\n");
1187  fRflFunc->Print();
1188  }
1189  if(fBkRFunc){
1190  printf("--- Background + reflections ---\n");
1191  fBkRFunc->Print();
1192  }
1193  if(fSecFunc){
1194  printf("--- Additional Gaussian peak ---\n");
1195  fSecFunc->Print();
1196  }
1197 }
Int_t fTypeOfFit4Sgn
pdg value of particle mass
Double_t fSigmaSgnErr
signal gaussian sigma
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
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 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)
TF1 * fTotFunc
fit function for second peak
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
fit parameters in 2nd peak fit function
Bool_t fFixSecMass
width of the 2nd peak
TF1 * fSecFunc
flag to fix the width of the 2nd peak
Double_t * sigma
TH1F * GetResidualsAndPulls(TH1 *hPulls=0x0, TH1 *hResidualTrend=0x0, TH1 *hPullsTrend=0x0, Double_t minrange=0, Double_t maxrange=-1)
int Int_t
Definition: External.C:63
Double_t fSecWidth
position of the 2nd peak
TString fFitOption
number of sigmas to veto the signal peak
TF1 * fBkRFunc
fit function for reflections
Int_t fNParsSec
switch off/on second peak (for D+->KKpi in Ds)
TF1 * CreateSignalFitFunction(TString fname, Double_t integral)
Bool_t fFixedMean
unc on signal gaussian sigma
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 MassFitter(Bool_t draw=kTRUE)
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
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
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 FitFunction4Sgn(Double_t *x, Double_t *par)
Int_t fPolDegreeBkg
background fit func
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
initialization for wa yield
Bool_t fReflections
background fit function (2nd step)
Bool_t fFixedSigma
switch for fix mean of gaussian