AliPhysics  d497547 (d497547)
 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 
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  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 //__________________________________________________________________________
224 
225  TVirtualFitter::SetDefaultFitter("Minuit");
226 
227  Double_t integralHisto=fHistoInvMass->Integral(fHistoInvMass->FindBin(fMinMass),fHistoInvMass->FindBin(fMaxMass),"width");
228 
229  fOnlySideBands = kTRUE;
230  fBkgFuncSb = CreateBackgroundFitFunction("funcbkgsb",integralHisto);
231  Int_t status=-1;
232  printf("\n--- First fit with only background on the side bands - Exclusion region = %.2f sigma ---\n",fNSigma4SideBands);
233  if(fTypeOfFit4Bkg==6){
235  fHistoInvMass->GetListOfFunctions()->Add(fBkgFuncSb);
236  fHistoInvMass->GetFunction(fBkgFuncSb->GetName())->SetBit(1<<9,kTRUE);
237  status=0;
238  }
239  }
240  else status=fHistoInvMass->Fit("funcbkgsb",Form("R,%s,+,0",fFitOption.Data()));
241  fBkgFuncSb->SetLineColor(kGray+1);
242  if (status != 0){
243  printf(" ---> Failed first fit with only background, minuit status = %d\n",status);
244  return 0;
245  }
246 
247  fOnlySideBands = kFALSE;
248  if(!fBkgFunc){
249  fBkgFunc = CreateBackgroundFitFunction("funcbkg",integralHisto);
250  for(Int_t ipar=0; ipar<fNParsBkg; ipar++) fBkgFunc->SetParameter(ipar,fBkgFuncSb->GetParameter(ipar));
251  }
252  fBkgFunc->SetLineColor(kGray+1);
253 
254  printf("\n--- Estimate signal counts in the peak region ---\n");
255  Double_t estimSignal=CheckForSignal(fMass,fSigmaSgn);
256  Bool_t doFinalFit=kTRUE;
257  if(estimSignal<0.){
258  if(draw) DrawFit();
259  estimSignal=0.;
260  doFinalFit=kFALSE;
261  }
262 
263  fRawYieldHelp=estimSignal; // needed for reflection normalization
264  if(!fBkgFuncRefit){
265  fBkgFuncRefit = CreateBackgroundFitFunction("funcbkgrefit",integralHisto);
266  for(Int_t ipar=0; ipar<fNParsBkg; ipar++) fBkgFuncRefit->SetParameter(ipar,fBkgFunc->GetParameter(ipar));
267  }
268  fBkgFuncRefit->SetLineColor(2);
269  fSigFunc = CreateSignalFitFunction("fsigfit",estimSignal);
270  if(fSecondPeak){
271  printf(" ---> Final fit includes a second inv. mass peak\n");
273  fSecFunc = CreateSecondPeakFunction("fsecpeak",estimSec);
274  }
275  if(fReflections){
276  printf(" ---> Final fit includes reflections\n");
277  fRflFunc = CreateReflectionFunction("freflect");
279  }
280  fTotFunc = CreateTotalFitFunction("funcmass");
281 
282  if(doFinalFit){
283  printf("\n--- Final fit with signal+background on the full range ---\n");
284  status=fHistoInvMass->Fit("funcmass",Form("R,%s,+,0",fFitOption.Data()));
285  if (status != 0){
286  printf(" ---> Failed fit with signal+background, minuit status = %d\n",status);
287  return 0;
288  }
289  }
290 
291  for(Int_t ipar=0; ipar<fNParsBkg; ipar++){
292  fBkgFuncRefit->SetParameter(ipar,fTotFunc->GetParameter(ipar));
293  fBkgFuncRefit->SetParError(ipar,fTotFunc->GetParError(ipar));
294  }
295  for(Int_t ipar=0; ipar<fNParsSig; ipar++){
296  fSigFunc->SetParameter(ipar,fTotFunc->GetParameter(ipar+fNParsBkg));
297  fSigFunc->SetParError(ipar,fTotFunc->GetParError(ipar+fNParsBkg));
298  }
299  if(fSecondPeak){
300  for(Int_t ipar=0; ipar<fNParsSec; ipar++){
301  fSecFunc->SetParameter(ipar,fTotFunc->GetParameter(ipar+fNParsBkg+fNParsSig));
302  fSecFunc->SetParError(ipar,fTotFunc->GetParError(ipar+fNParsBkg+fNParsSig));
303  }
304  fSecFunc->SetLineColor(kMagenta+1);
305  fSecFunc->SetLineStyle(3);
306  }
307  if(fReflections){
308  for(Int_t ipar=0; ipar<fNParsRfl; ipar++){
309  fRflFunc->SetParameter(ipar,fTotFunc->GetParameter(ipar+fNParsBkg+fNParsSig+fNParsSec));
310  fRflFunc->SetParError(ipar,fTotFunc->GetParError(ipar+fNParsBkg+fNParsSig+fNParsSec));
311  fBkRFunc->SetParameter(ipar+fNParsBkg,fTotFunc->GetParameter(ipar+fNParsBkg+fNParsSig+fNParsSec));
312  fBkRFunc->SetParError(ipar+fNParsBkg,fTotFunc->GetParError(ipar+fNParsBkg+fNParsSig+fNParsSec));
313  }
314  for(Int_t ipar=0; ipar<fNParsBkg; ipar++){
315  fBkRFunc->SetParameter(ipar,fTotFunc->GetParameter(ipar));
316  fBkRFunc->SetParError(ipar,fTotFunc->GetParError(ipar));
317  }
318  fRflFunc->SetLineColor(kGreen+1);
319  fBkRFunc->SetLineColor(kRed+1);
320  fBkRFunc->SetLineStyle(7);
321  }
322  fMass=fSigFunc->GetParameter(1);
323  fMassErr=fSigFunc->GetParError(1);
324  fSigmaSgn=fSigFunc->GetParameter(2);
325  fSigmaSgnErr=fSigFunc->GetParError(2);
326  fTotFunc->SetLineColor(4);
327  fRawYield=fTotFunc->GetParameter(fNParsBkg)/fHistoInvMass->GetBinWidth(1);
328  fRawYieldErr=fTotFunc->GetParError(fNParsBkg)/fHistoInvMass->GetBinWidth(1);
330  if(draw) DrawFit();
331  if(doFinalFit) return 1;
332  else return 2;
333 }
334 //______________________________________________________________________________
338 
339  TCanvas* c0=new TCanvas("c0");
340  DrawHere(c0);
341 }
342 //______________________________________________________________________________
343 void AliHFInvMassFitter::DrawHere(TVirtualPad* c, Double_t nsigma,Int_t writeFitInfo){
346 
347  gStyle->SetOptStat(0);
348  gStyle->SetCanvasColor(0);
349  gStyle->SetFrameFillColor(0);
350  c->cd();
351  fHistoInvMass->GetXaxis()->SetRangeUser(fMinMass,fMaxMass);
352  fHistoInvMass->SetMarkerStyle(20);
353  fHistoInvMass->SetMinimum(0.);
354  fHistoInvMass->Draw("PE");
355  if(fBkgFunc) fBkgFunc->Draw("same");
356  if(fBkgFuncRefit) fBkgFuncRefit->Draw("same");
357  if(fBkRFunc) fBkRFunc->Draw("same");
358  if(fRflFunc) fRflFunc->Draw("same");
359  if(fSecFunc) fSecFunc->Draw("same");
360  if(fTotFunc) fTotFunc->Draw("same");
361 
362  if(writeFitInfo > 0){
363  TPaveText *pinfos=new TPaveText(0.12,0.65,0.47,0.89,"NDC");
364  TPaveText *pinfom=new TPaveText(0.6,0.7,1.,.87,"NDC");
365  pinfos->SetBorderSize(0);
366  pinfos->SetFillStyle(0);
367  pinfom->SetBorderSize(0);
368  pinfom->SetFillStyle(0);
369  if(fTotFunc){
370  pinfom->SetTextColor(kBlue);
371  for(Int_t ipar=1; ipar<fNParsSig; ipar++){
372  pinfom->AddText(Form("%s = %.3f #pm %.3f",fTotFunc->GetParName(ipar+fNParsBkg),fTotFunc->GetParameter(ipar+fNParsBkg),fTotFunc->GetParError(ipar+fNParsBkg)));
373  }
374  pinfom->Draw();
375 
376  Double_t bkg,errbkg;
377  Background(nsigma,bkg,errbkg);
378  Double_t signif,errsignif;
379  Significance(nsigma,signif,errsignif);
380 
381  pinfos->AddText(Form("S = %.0f #pm %.0f ",fRawYield,fRawYieldErr));
382  pinfos->AddText(Form("B (%.0f#sigma) = %.0f #pm %.0f",nsigma,bkg,errbkg));
383  pinfos->AddText(Form("S/B = (%.0f#sigma) %.4f ",nsigma,fRawYield/bkg));
384  if(fRflFunc) pinfos->AddText(Form("Refl/Sig = %.3f #pm %.3f ",fRflFunc->GetParameter(0),fRflFunc->GetParError(0)));
385  pinfos->AddText(Form("Signif (%.0f#sigma) = %.1f #pm %.1f ",nsigma,signif,errsignif));
386  pinfos->Draw();
387  }
388  }
389  c->Update();
390  return;
391 }
392 //______________________________________________________________________________
396 
397  Double_t minForSig=mean-4.*sigma;
398  Double_t maxForSig=mean+4.*sigma;
399  Int_t binForMinSig=fHistoInvMass->FindBin(minForSig);
400  Int_t binForMaxSig=fHistoInvMass->FindBin(maxForSig);
401  Double_t sum=0.;
402  Double_t sumback=0.;
403  fBkgFunc->Print();
404  for(Int_t ibin=binForMinSig; ibin<=binForMaxSig; ibin++){
405  sum+=fHistoInvMass->GetBinContent(ibin);
406  sumback+=fBkgFunc->Eval(fHistoInvMass->GetBinCenter(ibin));
407  }
408  Double_t diffUnderPeak=(sum-sumback);
409  printf(" ---> IntegralUnderFitFunc=%f IntegralUnderHisto=%f EstimatedSignal=%f\n",sum,sumback,diffUnderPeak);
410  if(diffUnderPeak/TMath::Sqrt(sum)<1.){
411  printf(" ---> (Tot-Bkg)/sqrt(Tot)=%f ---> Likely no signal/\n",diffUnderPeak/TMath::Sqrt(sum));
412  return -1;
413  }
414  return diffUnderPeak*fHistoInvMass->GetBinWidth(1);
415 }
416 
417 //______________________________________________________________________________
421 
423  TF1* funcbkg = new TF1(fname.Data(),this,&AliHFInvMassFitter::FitFunction4Bkg,fMinMass,fMaxMass,fNParsBkg,"AliHFInvMassFitter","FitFunction4Bkg");
424  switch (fTypeOfFit4Bkg) {
425  case 0: //gaus+expo
426  funcbkg->SetParNames("BkgInt","Slope");
427  funcbkg->SetParameters(integral,-2.);
428  break;
429  case 1:
430  funcbkg->SetParNames("BkgInt","Slope");
431  funcbkg->SetParameters(integral,-100.);
432  break;
433  case 2:
434  funcbkg->SetParNames("BkgInt","Coef1","Coef2");
435  funcbkg->SetParameters(integral,-10.,5);
436  break;
437  case 3:
438  funcbkg->SetParNames("Const");
439  funcbkg->SetParameter(0,0.);
440  funcbkg->FixParameter(0,0.);
441  break;
442  case 4:
443  funcbkg->SetParNames("BkgInt","Coef1");
444  funcbkg->SetParameters(integral,0.5);
445  break;
446  case 5:
447  funcbkg->SetParNames("BkgInt","Coef1","Coef2");
448  funcbkg->SetParameters(integral,-10.,5.);
449  break;
450  case 6:
451  for(Int_t j=0;j<fNParsBkg; j++){
452  funcbkg->SetParName(j,Form("Coef%d",j));
453  funcbkg->SetParameter(j,0);
454  }
455  funcbkg->SetParameter(0,integral);
456  break;
457  default:
458  AliError(Form("Wrong choice of fTypeOfFit4Bkg (%d)",fTypeOfFit4Bkg));
459  delete funcbkg;
460  return 0x0;
461  break;
462  }
463  funcbkg->SetLineColor(kBlue+3);
464  return funcbkg;
465 }
466 //______________________________________________________________________________
470 
471  TF1* funcsec = new TF1(fname.Data(),this,&AliHFInvMassFitter::FitFunction4SecPeak,fMinMass,fMaxMass,3,"AliHFInvMassFitter","FitFunction4SecPeak");
472  funcsec->SetParameter(0,integsig);
473  funcsec->SetParameter(1,fSecMass);
474  if(fFixSecMass) funcsec->FixParameter(1,fSecMass);
475  funcsec->SetParameter(2,fSecWidth);
476  if(fFixSecWidth) funcsec->FixParameter(2,fSecWidth);
477  funcsec->SetParNames("SecPeakInt","SecMean","SecSigma");
478  return funcsec;
479 }
480 //______________________________________________________________________________
484  TF1* funcrfl = new TF1(fname.Data(),this,&AliHFInvMassFitter::FitFunction4Refl,fMinMass,fMaxMass,1,"AliHFInvMassFitter","FitFunction4Refl");
485  funcrfl->SetParameter(0,fRflOverSig);
486  funcrfl->SetParLimits(0,0.,1.);
487  if(fFixRflOverSig) funcrfl->FixParameter(0,fRflOverSig);
488  funcrfl->SetParNames("ReflOverS");
489  return funcrfl;
490 }
491 //______________________________________________________________________________
496  Int_t totParams=fNParsBkg+fNParsRfl;
497  TF1* fbr=new TF1(fname.Data(),this,&AliHFInvMassFitter::FitFunction4BkgAndRefl,fMinMass,fMaxMass,totParams,"AliHFInvMassFitter","FitFunction4BkgAndRefl");
498  for(Int_t ipar=0; ipar<fNParsBkg; ipar++){
499  fbr->SetParameter(ipar,fBkgFunc->GetParameter(ipar));
500  fbr->SetParName(ipar,fBkgFunc->GetParName(ipar));
501  }
502  for(Int_t ipar=0; ipar<fNParsRfl; ipar++){
503  fbr->SetParameter(ipar+fNParsBkg,fRflFunc->GetParameter(ipar));
504  fbr->SetParName(ipar+fNParsBkg,fRflFunc->GetParName(ipar));
505  // par limits not set because this function is not used for fitting but only for drawing
506  }
507  return fbr;
508 }
509 //______________________________________________________________________________
513 
515  TF1* funcsig = new TF1(fname.Data(),this,&AliHFInvMassFitter::FitFunction4Sgn,fMinMass,fMaxMass,fNParsSig,"AliHFInvMassFitter","FitFunction4Sgn");
516  if(fTypeOfFit4Sgn==kGaus){
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  if(fFixedSigma) funcsig->FixParameter(2,fSigmaSgn);
523  funcsig->SetParNames("SgnInt","Mean","Sigma");
524  }
525  if(fTypeOfFit4Sgn==k2Gaus){
526  funcsig->SetParameter(0,integsig);
527  if(fFixedRawYield>-0.1) funcsig->FixParameter(0,fFixedRawYield);
528  funcsig->SetParameter(1,fMass);
529  if(fFixedMean) funcsig->FixParameter(1,fMass);
530  funcsig->SetParameter(2,fSigmaSgn);
531  funcsig->SetParLimits(2,0.004,0.05);
532  if(fFixedSigma) funcsig->FixParameter(2,fSigmaSgn);
533  funcsig->SetParameter(3,0.2);
534  funcsig->SetParLimits(3,0.,1.);
535  funcsig->SetParameter(4,fSigmaSgn);
536  funcsig->SetParLimits(4,0.004,0.05);
537  funcsig->SetParNames("SgnInt","Mean","Sigma1","Frac","Sigma2");
538  }
539  return funcsig;
540 }
541 
542 //______________________________________________________________________________
546 
549  TF1* ftot=ftot=new TF1(fname.Data(),this,&AliHFInvMassFitter::FitFunction4Mass,fMinMass,fMaxMass,totParams,"AliHFInvMassFitter","FitFunction4Mass");
550  for(Int_t ipar=0; ipar<fNParsBkg; ipar++){
551  ftot->SetParameter(ipar,fBkgFunc->GetParameter(ipar));
552  ftot->SetParName(ipar,fBkgFunc->GetParName(ipar));
553  }
554  for(Int_t ipar=0; ipar<fNParsSig; ipar++){
555  ftot->SetParameter(ipar+fNParsBkg,fSigFunc->GetParameter(ipar));
556  ftot->SetParName(ipar+fNParsBkg,fSigFunc->GetParName(ipar));
557  Double_t parmin,parmax;
558  fSigFunc->GetParLimits(ipar,parmin,parmax);
559  ftot->SetParLimits(ipar+fNParsBkg,parmin,parmax);
560  }
561  if(fSecondPeak && fSecFunc){
562  for(Int_t ipar=0; ipar<fNParsSec; ipar++){
563  ftot->SetParameter(ipar+fNParsBkg+fNParsSig,fSecFunc->GetParameter(ipar));
564  ftot->SetParName(ipar+fNParsBkg+fNParsSig,fSecFunc->GetParName(ipar));
565  Double_t parmin,parmax;
566  fSecFunc->GetParLimits(ipar,parmin,parmax);
567  ftot->SetParLimits(ipar+fNParsBkg+fNParsSig,parmin,parmax);
568  }
569  }
570  if(fReflections && fRflFunc){
571  for(Int_t ipar=0; ipar<fNParsRfl; ipar++){
572  ftot->SetParameter(ipar+fNParsBkg+fNParsSig+fNParsSec,fRflFunc->GetParameter(ipar));
573  ftot->SetParName(ipar+fNParsBkg+fNParsSig+fNParsSec,fRflFunc->GetParName(ipar));
574  Double_t parmin,parmax;
575  fRflFunc->GetParLimits(ipar,parmin,parmax);
576  ftot->SetParLimits(ipar+fNParsBkg+fNParsSig+fNParsSec,parmin,parmax);
577  }
578  }
579  return ftot;
580 }
581 //__________________________________________________________________________
585 
587  if(fOnlySideBands && TMath::Abs(x[0]-fMass) < maxDeltaM) {
588  TF1::RejectPoint();
589  return 0;
590  }
591  if(fOnlySideBands && fSecondPeak && TMath::Abs(x[0]-fSecMass) < (fNSigma4SideBands*fSecWidth)){
592  TF1::RejectPoint();
593  return 0;
594  }
595  Double_t total=0;
596 
597  switch (fTypeOfFit4Bkg){
598  case 0:
599  //exponential
600  //exponential = A*exp(B*x) -> integral(exponential)=A/B*exp(B*x)](min,max)
601  //-> A = B*integral/(exp(B*max)-exp(B*min)) where integral can be written
602  //as integralTot- integralGaus (=par [2])
603  //Par:
604  // * [0] = integralBkg;
605  // * [1] = B;
606  //exponential = [1]*[0]/(exp([1]*max)-exp([1]*min))*exp([1]*x)
607  total = par[0]*par[1]/(TMath::Exp(par[1]*fMaxMass)-TMath::Exp(par[1]*fMinMass))*TMath::Exp(par[1]*x[0]);
608  // AliInfo("Background function set to: exponential");
609  break;
610  case 1:
611  //linear
612  //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)
613  // * [0] = integralBkg;
614  // * [1] = b;
615  total= par[0]/(fMaxMass-fMinMass)+par[1]*(x[0]-0.5*(fMaxMass+fMinMass));
616  // AliInfo("Background function set to: linear");
617  break;
618  case 2:
619  //parabola
620  //y=a+b*x+c*x**2 -> integral = a(max-min) + 1/2*b*(max^2-min^2) +
621  //+ 1/3*c*(max^3-min^3) ->
622  //a = (integral-1/2*b*(max^2-min^2)-1/3*c*(max^3-min^3))/(max-min)
623  // * [0] = integralBkg;
624  // * [1] = b;
625  // * [2] = c;
626  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));
627  // AliInfo("Background function set to: polynomial");
628  break;
629  case 3:
630  total=par[0];
631  break;
632  case 4:
633  //power function
634  //y=a(x-m_pi)^b -> integral = a/(b+1)*((max-m_pi)^(b+1)-(min-m_pi)^(b+1))
635  //
636  //a = integral*(b+1)/((max-m_pi)^(b+1)-(min-m_pi)^(b+1))
637  // * [0] = integralBkg;
638  // * [1] = b;
639  // a(power function) = [0]*([1]+1)/((max-m_pi)^([1]+1)-(min-m_pi)^([1]+1))*(x-m_pi)^[1]
640  {
641  Double_t mpi = TDatabasePDG::Instance()->GetParticle(211)->Mass();
642 
643  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]);
644  // AliInfo("Background function set to: powerlaw");
645  }
646  break;
647  case 5:
648  //power function wit exponential
649  //y=a*Sqrt(x-m_pi)*exp(-b*(x-m_pi))
650  {
651  Double_t mpi = TDatabasePDG::Instance()->GetParticle(211)->Mass();
652 
653  total = par[1]*TMath::Sqrt(x[0] - mpi)*TMath::Exp(-1.*par[2]*(x[0]-mpi));
654  // AliInfo("Background function set to: wit exponential");
655  }
656  break;
657  case 6:
658  // the following comment must be removed
659  // // pol 3, following convention for pol 2
660  // //y=a+b*x+c*x**2+d*x**3 -> integral = a(max-min) + 1/2*b*(max^2-min^2) +
661  // //+ 1/3*c*(max^3-min^3) + 1/4 d * (max^4-min^4) ->
662  // //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)
663  // // * [0] = integralBkg;
664  // // * [1] = b;
665  // // * [2] = c;
666  // // * [3] = d;
667  {
668  total=par[0];
669  for(Int_t it=1;it<=fPolDegreeBkg;it++){
670  total+=par[it]*TMath::Power(x[0]-fMassParticle,it)/TMath::Factorial(it);
671  }
672  }
673  break;
674  }
675  return total;
676 }
677 //_________________________________________________________________________
681 
682  // AliInfo("Signal function set to: Gaussian");
683  Double_t sigval=0;
684  switch (fTypeOfFit4Sgn){
685  case 0:
686  //gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
687  //Par:
688  // * [0] = integralSgn
689  // * [1] = mean
690  // * [2] = sigma
691  //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
692  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]);
693  break;
694  case 1:
695  //double gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
696  //Par:
697  // * [0] = integralSgn
698  // * [1] = mean
699  // * [2] = sigma1
700  // * [3] = 2nd gaussian ratio
701  // * [4] = deltaSigma
702  //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
703  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]);
704  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]);
705  sigval=par[0]*(g1+g2);
706  break;
707  }
708  fRawYieldHelp=par[0]/fHistoInvMass->GetBinWidth(1);
709  return sigval;
710 }
711 //_________________________________________________________________________
715  if(!fHistoTemplRfl) return 0;
716 
717  Int_t bin =fHistoTemplRfl->FindBin(x[0]);
718  Double_t value=fHistoTemplRfl->GetBinContent(bin);
719  Int_t binmin=fHistoTemplRfl->FindBin(fMinMass*1.00001);
720  Int_t binmax=fHistoTemplRfl->FindBin(fMaxMass*0.99999);
721  Double_t norm=fHistoTemplRfl->Integral(binmin,binmax)*fHistoTemplRfl->GetBinWidth(bin);
722  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
723  value+=fHistoTemplRfl->GetBinContent(bin-1)+fHistoTemplRfl->GetBinContent(bin+1);
724  value/=3.;
725  }
726  return par[0]*value/norm*fRawYieldHelp*fHistoInvMass->GetBinWidth(1);
727 }
728 //_________________________________________________________________________
732  Double_t bkg=FitFunction4Bkg(x,par);
733  Double_t refl=0;
734  if(fReflections) refl=FitFunction4Refl(x,&par[fNParsBkg]);
735  return bkg+refl;
736 }
737 //_________________________________________________________________________
741 
742  //gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
743  //Par:
744  // * [0] = integralSgn
745  // * [1] = mean
746  // * [2] = sigma
747  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]);
748  return secgaval;
749 }
750 //_________________________________________________________________________
754 
755  Double_t bkg=FitFunction4Bkg(x,par);
756  Double_t sig=FitFunction4Sgn(x,&par[fNParsBkg]);
757  Double_t sec=0.;
758  if(fSecondPeak) sec=FitFunction4SecPeak(x,&par[fNParsBkg+fNParsSig]);
759  Double_t refl=0;
760  if(fReflections) refl=FitFunction4Refl(x,&par[fNParsBkg+fNParsSig+fNParsSec]);
761  return bkg+sig+sec+refl;
762 }
763 
764 //_________________________________________________________________________
765 void AliHFInvMassFitter::Signal(Double_t nOfSigma,Double_t &signal,Double_t &errsignal) const {
768 
769  Double_t minMass=fMass-nOfSigma*fSigmaSgn;
770  Double_t maxMass=fMass+nOfSigma*fSigmaSgn;
771  Signal(minMass,maxMass,signal,errsignal);
772  return;
773 }
774 
775 //_________________________________________________________________________
776 void AliHFInvMassFitter::Signal(Double_t min, Double_t max, Double_t &signal,Double_t &errsignal) const {
779 
780  signal=fSigFunc->Integral(min, max)/(Double_t)fHistoInvMass->GetBinWidth(1);
781  errsignal=(fRawYieldErr/fRawYield)*signal;/*assume relative error is the same as for total integral*/
782  return;
783 }
784 
785 //___________________________________________________________________________
786 void AliHFInvMassFitter::Background(Double_t nOfSigma,Double_t &background,Double_t &errbackground) const {
789 
790  Double_t minMass=fMass-nOfSigma*fSigmaSgn;
791  Double_t maxMass=fMass+nOfSigma*fSigmaSgn;
792  Background(minMass,maxMass,background,errbackground);
793  return;
794 
795 }
796 //___________________________________________________________________________
797 void AliHFInvMassFitter::Background(Double_t min, Double_t max, Double_t &background,Double_t &errbackground) const {
800 
801  TF1 *funcbkg=0x0;
802  if(fBkgFuncRefit) funcbkg=fBkgFuncRefit;
803  else if(fBkgFunc) funcbkg=fBkgFunc;
804  if(!funcbkg){
805  AliError("Bkg function not found!");
806  return;
807  }
808 
809  Double_t intB=funcbkg->GetParameter(0);
810  Double_t intBerr=funcbkg->GetParError(0);
811 
812  //relative error evaluation: from histo
813 
815  Int_t rightBand=fHistoInvMass->FindBin(fMass+fNSigma4SideBands*fSigmaSgn);
816  intB=fHistoInvMass->Integral(1,leftBand)+fHistoInvMass->Integral(rightBand,fHistoInvMass->GetNbinsX());
817  Double_t sum2=0;
818  for(Int_t i=1;i<=leftBand;i++){
819  sum2+=fHistoInvMass->GetBinError(i)*fHistoInvMass->GetBinError(i);
820  }
821  for(Int_t i=rightBand; i<=fHistoInvMass->GetNbinsX();i++){
822  sum2+=fHistoInvMass->GetBinError(i)*fHistoInvMass->GetBinError(i);
823  }
824 
825  intBerr=TMath::Sqrt(sum2);
826 
827  background=funcbkg->Integral(min,max)/(Double_t)fHistoInvMass->GetBinWidth(1);
828  errbackground=intBerr/intB*background;
829 
830  return;
831 
832 }
833 //__________________________________________________________________________
834 
835 void AliHFInvMassFitter::Significance(Double_t nOfSigma,Double_t &significance,Double_t &errsignificance) const {
838 
839  Double_t minMass=fMass-nOfSigma*fSigmaSgn;
840  Double_t maxMass=fMass+nOfSigma*fSigmaSgn;
841  Significance(minMass, maxMass, significance, errsignificance);
842 
843  return;
844 }
845 
846 //__________________________________________________________________________
847 
848 void AliHFInvMassFitter::Significance(Double_t min, Double_t max, Double_t &significance,Double_t &errsignificance) const {
851 
852  Double_t background,errbackground;
853  Background(min,max,background,errbackground);
854 
855  if (fRawYield+background <= 0.){
856  significance=-1;
857  errsignificance=0;
858  return;
859  }
860 
861  AliVertexingHFUtils::ComputeSignificance(fRawYield,fRawYieldErr,background,errbackground,significance,errsignificance);
862 
863  return;
864 }
865 //________________________________________________________________________
869 
870  TH1F *hCp=(TH1F*)fHistoInvMass->Clone("htemp");
871  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
872  Double_t estimateslope=(hCp->GetBinContent(hCp->FindBin(fMass+3.5*fSigmaSgn))-hCp->GetBinContent(hCp->FindBin(fMass-3.5*fSigmaSgn)))/(7*fSigmaSgn);// first rough estimate
873 
875  TF1 *funcbkg,*funcPrev=0x0;
877  funcbkg = new TF1(Form("temp%d",fCurPolDegreeBkg),this,&AliHFInvMassFitter::BackFitFuncPolHelper,fMinMass,fMaxMass,fCurPolDegreeBkg+1,"AliHFInvMassFitter","BackFitFuncPolHelper");
878  if(funcPrev){
879  for(Int_t j=0;j<fCurPolDegreeBkg;j++){// now is +1 degree w.r.t. previous fit funct
880  funcbkg->SetParameter(j,funcPrev->GetParameter(j));
881  }
882  delete funcPrev;
883  }
884  else{
885  funcbkg->SetParameter(0,estimatecent);
886  funcbkg->SetParameter(1,estimateslope);
887  }
888  printf(" ---> Pre-fit of background with pol degree %d ---\n",fCurPolDegreeBkg);
889  hCp->Fit(funcbkg,"REMN","");
890  funcPrev=(TF1*)funcbkg->Clone("ftemp");
891  delete funcbkg;
893  }
894 
895  for(Int_t j=0;j<=fPolDegreeBkg;j++){
896  fback->SetParameter(j,funcPrev->GetParameter(j));
897  fback->SetParError(j,funcPrev->GetParError(j));
898  }
899  printf(" ---> Final background fit with pol degree %d ---\n",fPolDegreeBkg);
900  hCp->Fit(fback,"REMN","");// THIS IS JUST TO SET NOT ONLY THE PARAMETERS BUT ALSO chi2, etc...
901 
902 
903  // The following lines might be useful for debugging
904  // TCanvas *cDebug=new TCanvas();
905  // cDebug->cd();
906  // hCp->Draw();
907  // TString strout=Form("Test%d.root",(Int_t)fhistoInvMass->GetBinContent(fhistoInvMass->FindBin(fMass)));
908  // cDebug->Print(strout.Data());
909  // delete cDebug;
910 
911  delete funcPrev;
912  delete hCp;
913  return kTRUE;
914 
915 }
916 // _______________________________________________________________________
920 
922  if(fOnlySideBands && TMath::Abs(x[0]-fMass) < maxDeltaM) {
923  TF1::RejectPoint();
924  return 0;
925  }
926  Double_t back=par[0];
927  for(Int_t it=1;it<=fCurPolDegreeBkg;it++){
928  back+=par[it]*TMath::Power(x[0]-fMassParticle,it)/TMath::Factorial(it);
929  }
930  return back;
931 }
932 // _______________________________________________________________________
941 
942  if(!h){
943  fReflections=kFALSE;
944  return 0x0;
945  }
946  fHistoTemplRfl=(TH1F*)h->Clone("hTemplRfl");
947  opt.ToLower();
948  fReflections=kTRUE;
949  printf("\n--- Reflection templates from simulation ---\n");
950  if(opt.Contains("templ")){
951  printf(" ---> Reflection contribution using directly the histogram from simulation\n");
952  return fHistoTemplRfl;
953  }
954 
955  TF1 *f=0x0;
956  Bool_t isPoissErr=kTRUE;
957  Double_t xMinForFit=h->GetBinLowEdge(1);
958  Double_t xMaxForFit=h->GetXaxis()->GetBinUpEdge(h->GetNbinsX());
959  if(minRange>=0 && maxRange>=0){
960  xMinForFit=TMath::Max(minRange,h->GetBinLowEdge(1));
961  xMaxForFit=TMath::Min(maxRange,h->GetXaxis()->GetBinUpEdge(h->GetNbinsX()));
962  }
963  if(opt.EqualTo("1gaus") || opt.EqualTo("singlegaus")){
964  printf(" ---> Reflection contribution from single-Gaussian fit to histogram from simulation\n");
965  f=new TF1("mygaus","gaus",xMinForFit,xMaxForFit);
966  f->SetParameter(0,h->GetMaximum());
967  // f->SetParLimits(0,0,100.*h->Integral());
968  f->SetParameter(1,1.865);
969  f->SetParameter(2,0.050);
970  fHistoTemplRfl->Fit(f,"REM","");//,h->GetBinLowEdge(1),h->GetXaxis()->GetBinUpEdge(h->GetNbinsX()));
971  }
972  else if(opt.EqualTo("2gaus") || opt.EqualTo("doublegaus")){
973  printf(" ---> Reflection contribution from double-Gaussian fit to histogram from simulation\n");
974  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);
975  f->SetParameter(0,h->GetMaximum());
976  // f->SetParLimits(0,0,100.*h->Integral());
977  f->SetParLimits(3,0.,1.);
978  f->SetParameter(3,0.5);
979  f->SetParameter(1,1.84);
980  f->SetParameter(2,0.050);
981  f->SetParameter(4,1.88);
982  f->SetParameter(5,0.050);
983  fHistoTemplRfl->Fit(f,"REM","");//,h->GetBinLowEdge(1),h->GetXaxis()->GetBinUpEdge(h->GetNbinsX()));
984  }
985  else if(opt.EqualTo("pol3")){
986  printf(" ---> Reflection contribution from pol3 fit to histogram from simulation\n");
987  f=new TF1("mypol3","pol3",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,"REM","");
992  // Printf("We USED %d POINTS in the Fit",f->GetNumberFitPoints());
993  }
994  else if(opt.EqualTo("pol6")){
995  printf(" ---> Reflection contribution from pol6 fit to histogram from simulation\n");
996  f=new TF1("mypol6","pol6",xMinForFit,xMaxForFit);
997  f->SetParameter(0,h->GetMaximum());
998  // f->SetParLimits(0,0,100.*h->Integral());
999  // Hard to initialize the other parameters...
1000  fHistoTemplRfl->Fit(f,"RLEMI","");//,h->GetBinLowEdge(1),h->GetXaxis()->GetBinUpEdge(h->GetNbinsX()));
1001  }
1002  else{
1003  // no good option passed
1004  printf(" ---> Bad option for reflection configuration -> reflections will not be included in the fit\n");
1005  fReflections=kFALSE;
1006  delete fHistoTemplRfl;
1007  fHistoTemplRfl=0x0;
1008  return 0x0;
1009  }
1010 
1011  // Fill fHistoTemplRfl with values of fit function
1012  if(f){
1013  for(Int_t j=1;j<=fHistoTemplRfl->GetNbinsX();j++){
1014  fHistoTemplRfl->SetBinContent(j,f->Integral(fHistoTemplRfl->GetBinLowEdge(j),fHistoTemplRfl->GetXaxis()->GetBinUpEdge(j))/fHistoTemplRfl->GetBinWidth(j));
1015  if(fHistoTemplRfl->GetBinContent(j)>=0.&&TMath::Abs(h->GetBinError(j)*h->GetBinError(j)-h->GetBinContent(j))>0.1*h->GetBinContent(j))isPoissErr=kFALSE;
1016  }
1017  for(Int_t j=1;j<=fHistoTemplRfl->GetNbinsX();j++){
1018  if(isPoissErr){
1019  fHistoTemplRfl->SetBinError(j,TMath::Sqrt(fHistoTemplRfl->GetBinContent(j)));
1020  }
1021  else fHistoTemplRfl->SetBinError(j,0.001*fHistoTemplRfl->GetBinContent(j));
1022  }
1023  fReflections=kTRUE;
1024  return fHistoTemplRfl;
1025  }else{
1026  printf(" ---> Fit to MC template for reflection failed -> reflections will not be included in the fit\n");
1027  fReflections=kFALSE;
1028  delete fHistoTemplRfl;
1029  fHistoTemplRfl=0x0;
1030  return 0x0;
1031  }
1032  return 0x0;
1033 }
1034 // _______________________________________________________________________
1039  // else (default) mean of gaussian fit
1040 
1041  Double_t massVal=fMass;
1042  switch (pdgCode) {
1043  case 411:
1044  case 421:
1045  case 431:
1046  case 4122:
1047  massVal=TDatabasePDG::Instance()->GetParticle(pdgCode)->Mass();
1048  break;
1049  case 413:
1050  massVal=TDatabasePDG::Instance()->GetParticle(pdgCode)->Mass();
1051  massVal-=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1052  break;
1053  default:
1054  massVal=fMass;
1055  break;
1056  }
1057 
1058  Double_t minMass=massVal-nOfSigma*fSigmaSgn;
1059  Double_t maxMass=massVal+nOfSigma*fSigmaSgn;
1060  return GetRawYieldBinCounting(errRyBC,minMass,maxMass,option);
1061 }
1062 // _______________________________________________________________________
1068 
1069  Int_t minBinSum=fHistoInvMass->FindBin(minMass);
1070  Int_t maxBinSum=fHistoInvMass->FindBin(maxMass);
1071  if(minBinSum<1){
1072  printf("Left range for bin counting smaller than allowed by histogram axis, setting it to the lower edge of the first histo bin\n");
1073  minBinSum=1;
1074  }
1075  if(maxBinSum>fHistoInvMass->GetNbinsX()){
1076  printf("Right range for bin counting larger than allowed by histogram axis, setting it to the upper edge of the last histo bin\n");
1077  maxBinSum=fHistoInvMass->GetNbinsX();
1078  }
1079  Double_t cntSig=0.;
1080  Double_t cntErr=0.;
1081  errRyBC=0;
1082  TF1* fbackground=fBkgFunc;
1083  if(option==1) fbackground=fBkgFuncRefit;
1084  if(!fbackground) return 0.;
1085 
1086  for(Int_t jb=minBinSum; jb<=maxBinSum; jb++){
1087  Double_t cntTot=fHistoInvMass->GetBinContent(jb);
1088  Double_t cntBkg=fbackground->Integral(fHistoInvMass->GetBinLowEdge(jb),fHistoInvMass->GetBinLowEdge(jb)+fHistoInvMass->GetBinWidth(jb))/fHistoInvMass->GetBinWidth(jb);
1089  Double_t cntRefl=0;
1090  if(option==1 && fRflFunc) cntRefl=fRflFunc->Integral(fHistoInvMass->GetBinLowEdge(jb),fHistoInvMass->GetBinLowEdge(jb)+fHistoInvMass->GetBinWidth(jb))/fHistoInvMass->GetBinWidth(jb);
1091  //Double_t cntBkg=fbackground->Eval(fHistoInvMass->GetBinCenter(jb));
1092  cntSig+=(cntTot-cntBkg-cntRefl);
1093  cntErr+=(fHistoInvMass->GetBinError(jb)*fHistoInvMass->GetBinError(jb));
1094  }
1095  errRyBC=TMath::Sqrt(cntErr);
1096  return cntSig;
1097 }
1098 
1099 
1100 // _______________________________________________________________________
1101 TH1F* AliHFInvMassFitter::GetResidualsAndPulls(TH1 *hPulls,TH1 *hResidualTrend,TH1 *hPullsTrend, Double_t minrange,Double_t maxrange){
1102 
1104 
1105  Int_t binmi=fHistoInvMass->FindBin(fMinMass*1.001);
1106  Int_t binma=fHistoInvMass->FindBin(fMaxMass*0.9999);
1107  if(maxrange>minrange){
1108  binmi=fHistoInvMass->FindBin(minrange*1.001);
1109  binma=fHistoInvMass->FindBin(maxrange*0.9999);
1110  }
1111  if(hResidualTrend){
1112  //fHistoInvMass->Copy(hResidualTrend);
1113  hResidualTrend->SetBins(fHistoInvMass->GetNbinsX(),fHistoInvMass->GetXaxis()->GetXmin(),fHistoInvMass->GetXaxis()->GetXmax());
1114  hResidualTrend->SetName(Form("%s_residualTrend",fHistoInvMass->GetName()));
1115  hResidualTrend->SetTitle(Form("%s (Residuals)",fHistoInvMass->GetTitle()));
1116  hResidualTrend->SetMarkerStyle(20);
1117  hResidualTrend->SetMarkerSize(1.0);
1118  hResidualTrend->Reset();
1119  }
1120  if(hPullsTrend){
1121  hPullsTrend->SetBins(fHistoInvMass->GetNbinsX(),fHistoInvMass->GetXaxis()->GetXmin(),fHistoInvMass->GetXaxis()->GetXmax());
1122  hPullsTrend->Reset();
1123  hPullsTrend->SetName(Form("%s_pullTrend",fHistoInvMass->GetName()));
1124  hPullsTrend->SetTitle(Form("%s (Pulls)",fHistoInvMass->GetTitle()));
1125  hPullsTrend->SetMarkerStyle(20);
1126  hPullsTrend->SetMarkerSize(1.0);
1127  }
1128  if(hPulls){
1129  hPulls->SetName(Form("%s_pulls",fHistoInvMass->GetName()));
1130  hPulls->SetTitle(Form("%s ; Pulls",fHistoInvMass->GetTitle()));
1131  hPulls->SetBins(40,-10,10);
1132  hPulls->Reset();
1133  }
1134 
1135  Double_t res=-1.e-6,min=1.e+12,max=-1.e+12;
1136  TArrayD *arval=new TArrayD(binma-binmi+1);
1137  for(Int_t jst=1;jst<=fHistoInvMass->GetNbinsX();jst++){
1138 
1139  res=fHistoInvMass->GetBinContent(jst)-fTotFunc->Integral(fHistoInvMass->GetBinLowEdge(jst),fHistoInvMass->GetBinLowEdge(jst)+fHistoInvMass->GetBinWidth(jst))/fHistoInvMass->GetBinWidth(jst);
1140  if(jst>=binmi&&jst<=binma){
1141  arval->AddAt(res,jst-binmi);
1142  if(res<min)min=res;
1143  if(res>max)max=res;
1144  }
1145  // Printf("Res = %f from %f - %f",res,fHistoInvMass->GetBinContent(jst),fTotFunc->Integral(fHistoInvMass->GetBinLowEdge(jst),fHistoInvMass->GetBinLowEdge(jst)+fHistoInvMass->GetBinWidth(jst))/fHistoInvMass->GetBinWidth(jst));
1146  if(hResidualTrend){
1147  hResidualTrend->SetBinContent(jst,res);
1148  hResidualTrend->SetBinError(jst,fHistoInvMass->GetBinError(jst));
1149  }
1150  if(hPulls){
1151  if(jst>=binmi&&jst<=binma)hPulls->Fill(res/fHistoInvMass->GetBinError(jst));
1152  }
1153  if(hPullsTrend){
1154  hPullsTrend->SetBinContent(jst,res/fHistoInvMass->GetBinError(jst));
1155  hPullsTrend->SetBinError(jst,0.0001);
1156  }
1157  }
1158  if(hResidualTrend)hResidualTrend->GetXaxis()->SetRange(binmi,binma);
1159  if(hPullsTrend){
1160  hPullsTrend->GetXaxis()->SetRange(binmi,binma);
1161  hPullsTrend->SetMinimum(-7);
1162  hPullsTrend->SetMaximum(+7);
1163  }
1164  if(TMath::Abs(min)>TMath::Abs(max))max=min;
1165 
1166  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);
1167  for(Int_t j=0;j<binma-binmi+1;j++){
1168  hout->Fill(arval->At(j));
1169  }
1170  hout->Sumw2();
1171  hout->Fit("gaus","LEM","",-TMath::Abs(max)*1.2,TMath::Abs(max)*1.2);
1172 
1173  if(hPulls){
1174  hPulls->Sumw2();
1175  hPulls->Fit("gaus","LEM","",-3,3);
1176  }
1177  delete arval;
1178  return hout;
1179 }
1180 // _______________________________________________________________________
1184  if(fBkgFunc){
1185  printf("--- Background function in 1st step fit to side bands ---\n");
1186  fBkgFunc->Print();
1187  }
1188  if(fTotFunc){
1189  printf("--- Total fit function from 2nd step fit ---\n");
1190  fTotFunc->Print();
1191  }
1192  if(fBkgFuncRefit){
1193  printf("--- Background function in 2nd step fit ---\n");
1194  fBkgFuncRefit->Print();
1195  }
1196  if(fRflFunc){
1197  printf("--- Reflections ---\n");
1198  fRflFunc->Print();
1199  }
1200  if(fBkRFunc){
1201  printf("--- Background + reflections ---\n");
1202  fBkRFunc->Print();
1203  }
1204  if(fSecFunc){
1205  printf("--- Additional Gaussian peak ---\n");
1206  fSecFunc->Print();
1207  }
1208 }
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 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)
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 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 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