PWG3/vertexingHF/AliHFMassFitter.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2009, 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 /////////////////////////////////////////////////////////////
  17 //
  18 // AliHFMassFitter for the fit of invariant mass distribution
  19 // of charmed mesons
  20 //
  21 // Author: C.Bianchin, chiara.bianchin@pd.infn.it
  22 /////////////////////////////////////////////////////////////
  23
  24 #include <TCanvas.h>
  25
  26 #include "AliHFMassFitter.h"
  27
  28
  29 ClassImp(AliHFMassFitter)
  30
  31  //************** constructors
  32 AliHFMassFitter::AliHFMassFitter() :
  33   TNamed(),
  34   fhistoInvMass(0),
  35   fminMass(0),
  36   fmaxMass(0),
  37   fNbin(0),
  38   fParsSize(1),
  39   fFitPars(0),
  40   fWithBkg(0),
  41   ftypeOfFit4Bkg(0),
  42   ftypeOfFit4Sgn(0),
  43   ffactor(1),
  44   fntuParam(0),
  45   fMass(1.85),
  46   fSigmaSgn(0.012),
  47   fSideBands(0),
  48   fcounter(0)
  49 {
  50   // default constructor
  51
  52   cout<<"Default constructor"<<endl;
  53 }
  54
  55 //___________________________________________________________________________
  56
  57 AliHFMassFitter::AliHFMassFitter (TH1F *histoToFit, Double_t minvalue, Double_t maxvalue, Int_t rebin,Int_t fittypeb,Int_t fittypes):
  58  TNamed(),
  59  fhistoInvMass(0),
  60  fminMass(0),
  61  fmaxMass(0),
  62  fNbin(0),
  63  fParsSize(1),
  64  fFitPars(0),
  65  fWithBkg(0),
  66  ftypeOfFit4Bkg(0),
  67  ftypeOfFit4Sgn(0),
  68  ffactor(1),
  69  fntuParam(0),
  70  fMass(1.85),
  71  fSigmaSgn(0.012),
  72  fSideBands(0),
  73  fcounter(0)
  74 {
  75   // standard constructor
  76
  77   fhistoInvMass= (TH1F*)histoToFit->Clone("fhistoInvMass");
  78   fminMass=minvalue;
  79   fmaxMass=maxvalue;
  80   if(rebin!=1) RebinMass(rebin);
  81   else fNbin=(Int_t)fhistoInvMass->GetNbinsX();
  82   ftypeOfFit4Bkg=fittypeb;
  83   ftypeOfFit4Sgn=fittypes;
  84   if(ftypeOfFit4Bkg!=0 && ftypeOfFit4Bkg!=1 && ftypeOfFit4Bkg!=2) fWithBkg=kFALSE;
  85   else fWithBkg=kTRUE;
  86   if (!fWithBkg) cout<<"Fit Histogram of Signal only"<<endl;
  87   else  cout<<"Type of fit For Background = "<<ftypeOfFit4Bkg<<endl;
  88
  89   ComputeParSize();
  90   fFitPars=new Float_t[fParsSize];
  91 }
  92
  93
  94 AliHFMassFitter::AliHFMassFitter(const AliHFMassFitter &mfit):
  95   TNamed(),
  96   fhistoInvMass(mfit.fhistoInvMass),
  97   fminMass(mfit.fminMass),
  98   fmaxMass(mfit.fmaxMass),
  99   fNbin(mfit.fNbin),
 100   fParsSize(mfit.fParsSize),
 101   fFitPars(0),
 102   fWithBkg(mfit.fWithBkg),
 103   ftypeOfFit4Bkg(mfit.ftypeOfFit4Bkg),
 104   ftypeOfFit4Sgn(mfit.ftypeOfFit4Sgn),
 105   ffactor(mfit.ffactor),
 106   fntuParam(mfit.fntuParam),
 107   fMass(mfit.fMass),
 108   fSigmaSgn(mfit.fSigmaSgn),
 109   fSideBands(mfit.fSideBands),
 110   fcounter(mfit.fcounter)
 111
 112 {
 113   //copy constructor
 114
 115   if(mfit.fParsSize > 0){
 116     fFitPars=new Float_t[fParsSize];
 117     memcpy(fFitPars,mfit.fFitPars,mfit.fParsSize*sizeof(Float_t));
 118   }
 119   //for(Int_t i=0;i<fParsSize;i++) fFitPars[i]=mfit.fFitPars[i];
 120 }
 121
 122 //_________________________________________________________________________
 123
 124 AliHFMassFitter::~AliHFMassFitter() {
 125   cout<<"AliHFMassFitter destructor called"<<endl;
 126   if(fhistoInvMass) delete   fhistoInvMass;
 127   if(fntuParam)     delete   fntuParam;
 128   if(fFitPars)      delete[] fFitPars;
 129   fcounter = 0;
 130 }
 131
 132 //_________________________________________________________________________
 133
 134 AliHFMassFitter& AliHFMassFitter::operator=(const AliHFMassFitter &mfit){
 135
 136   //assignment operator
 137
 138   if(&mfit == this) return *this;
 139   fhistoInvMass= mfit.fhistoInvMass;
 140   fminMass= mfit.fminMass;
 141   fmaxMass= mfit.fmaxMass;
 142   fNbin= mfit.fNbin;
 143   fParsSize= mfit.fParsSize;
 144   fWithBkg= mfit.fWithBkg;
 145   ftypeOfFit4Bkg= mfit.ftypeOfFit4Bkg;
 146   ftypeOfFit4Sgn= mfit.ftypeOfFit4Sgn;
 147   ffactor= mfit.ffactor;
 148   fntuParam= mfit.fntuParam;
 149   fMass= mfit.fMass;
 150   fSigmaSgn= mfit.fSigmaSgn;
 151   fSideBands= mfit.fSideBands;
 152   fcounter= mfit.fcounter;
 153   if(mfit.fParsSize > 0){
 154     if(fFitPars) delete[] fFitPars;
 155     fFitPars=new Float_t[fParsSize];
 156     memcpy(fFitPars,mfit.fFitPars,mfit.fParsSize*sizeof(Float_t));
 157   }
 158 // fFitPars=new Float_t[fParsSize];
 159 //   for(Int_t i=0;i<fParsSize;i++) fFitPars[i]=mfit.fFitPars[i];
 160
 161   return *this;
 162 }
 163
 164 //************ tools & settings
 165
 166 //__________________________________________________________________________
 167
 168 void AliHFMassFitter::ComputeParSize() {
 169   switch (ftypeOfFit4Bkg) {//npar backround func
 170   case 0:
 171     fParsSize = 2*3;
 172     break;
 173   case 1:
 174     fParsSize = 2*3;
 175     break;
 176   case 2:
 177     fParsSize = 3*3;
 178     break;
 179   case 3:
 180     fParsSize = 1*3;
 181     break;
 182   default:
 183     cout<<"Error in computing fParsSize: check ftypeOfFit4Bkg"<<endl;
 184     break;
 185   }
 186
 187   fParsSize += 3; // npar refl
 188   fParsSize += 3; // npar signal gaus
 189
 190   fParsSize*=2;   // add errors
 191   cout<<"Parameters array size "<<fParsSize<<endl;
 192 }
 193
 194 //___________________________________________________________________________
 195 void AliHFMassFitter::SetHisto(TH1F *histoToFit){
 196   fhistoInvMass = (TH1F*)histoToFit->Clone();
 197 }
 198
 199 //___________________________________________________________________________
 200
 201   void AliHFMassFitter::SetType(Int_t fittypeb, Int_t fittypes) {
 202     ftypeOfFit4Bkg = fittypeb;
 203     ftypeOfFit4Sgn = fittypes;
 204     ComputeParSize();
 205     fFitPars = new Float_t[fParsSize];
 206
 207 }
 208
 209 //___________________________________________________________________________
 210
 211 void AliHFMassFitter::Reset() {
 212   fMass=1.85;
 213   fSigmaSgn=0.012;
 214   delete fhistoInvMass;
 215 }
 216
 217 //_________________________________________________________________________
 218
 219 void AliHFMassFitter::InitNtuParam(char *ntuname) {
 220   // Create ntuple to keep fit parameters
 221   fntuParam=0;
 222   fntuParam=new TNtuple(ntuname,"Contains fit parameters","intbkg1:slope1:conc1:intGB:meanGB:sigmaGB:intbkg2:slope2:conc2:inttot:slope3:conc3:intsgn:meansgn:sigmasgn:intbkg1Err:slope1Err:conc1Err:intGBErr:meanGBErr:sigmaGBErr:intbkg2Err:slope2Err:conc2Err:inttotErr:slope3Err:conc3Err:intsgnErr:meansgnErr:sigmasgnErr");
 223
 224 }
 225
 226 //_________________________________________________________________________
 227
 228 void AliHFMassFitter::FillNtuParam() {
 229   // Fill ntuple with fit parameters
 230
 231   Float_t nothing=0.;
 232
 233   if (ftypeOfFit4Bkg==2) {
 234       fntuParam->SetBranchAddress("intbkg1",&fFitPars[0]);
 235       fntuParam->SetBranchAddress("slope1",&fFitPars[1]);
 236       fntuParam->SetBranchAddress("conc1",&fFitPars[2]);
 237       fntuParam->SetBranchAddress("intGB",&fFitPars[3]);
 238       fntuParam->SetBranchAddress("meanGB",&fFitPars[4]);
 239       fntuParam->SetBranchAddress("sigmaGB",&fFitPars[5]);
 240       fntuParam->SetBranchAddress("intbkg2",&fFitPars[6]);
 241       fntuParam->SetBranchAddress("slope2",&fFitPars[7]);
 242       fntuParam->SetBranchAddress("conc2",&fFitPars[8]);
 243       fntuParam->SetBranchAddress("inttot",&fFitPars[9]);
 244       fntuParam->SetBranchAddress("slope3",&fFitPars[10]);
 245       fntuParam->SetBranchAddress("conc3",&fFitPars[11]);
 246       fntuParam->SetBranchAddress("intsgn",&fFitPars[12]);
 247       fntuParam->SetBranchAddress("meansgn",&fFitPars[13]);
 248       fntuParam->SetBranchAddress("sigmasgn",&fFitPars[14]);
 249
 250       fntuParam->SetBranchAddress("intbkg1Err",&fFitPars[15]);
 251       fntuParam->SetBranchAddress("slope1Err",&fFitPars[16]);
 252       fntuParam->SetBranchAddress("conc1Err",&fFitPars[17]);
 253       fntuParam->SetBranchAddress("intGBErr",&fFitPars[18]);
 254       fntuParam->SetBranchAddress("meanGBErr",&fFitPars[19]);
 255       fntuParam->SetBranchAddress("sigmaGBErr",&fFitPars[20]);
 256       fntuParam->SetBranchAddress("intbkg2Err",&fFitPars[21]);
 257       fntuParam->SetBranchAddress("slope2Err",&fFitPars[22]);
 258       fntuParam->SetBranchAddress("conc2Err",&fFitPars[23]);
 259       fntuParam->SetBranchAddress("inttotErr",&fFitPars[24]);
 260       fntuParam->SetBranchAddress("slope3Err",&fFitPars[25]);
 261       fntuParam->SetBranchAddress("conc3Err",&fFitPars[26]);
 262       fntuParam->SetBranchAddress("intsgnErr",&fFitPars[27]);
 263       fntuParam->SetBranchAddress("meansgnErr",&fFitPars[28]);
 264       fntuParam->SetBranchAddress("sigmasgnErr",&fFitPars[29]);
 265
 266   } else {
 267
 268     if(ftypeOfFit4Bkg==3){
 269       fntuParam->SetBranchAddress("intbkg1",&fFitPars[0]);
 270       fntuParam->SetBranchAddress("slope1",&nothing);
 271       fntuParam->SetBranchAddress("conc1",&nothing);
 272       fntuParam->SetBranchAddress("intGB",&fFitPars[1]);
 273       fntuParam->SetBranchAddress("meanGB",&fFitPars[2]);
 274       fntuParam->SetBranchAddress("sigmaGB",&fFitPars[3]);
 275       fntuParam->SetBranchAddress("intbkg2",&fFitPars[4]);
 276       fntuParam->SetBranchAddress("slope2",&nothing);
 277       fntuParam->SetBranchAddress("conc2",&nothing);
 278       fntuParam->SetBranchAddress("inttot",&fFitPars[6]);
 279       fntuParam->SetBranchAddress("slope3",&nothing);
 280       fntuParam->SetBranchAddress("conc3",&nothing);
 281       fntuParam->SetBranchAddress("intsgn",&fFitPars[6]);
 282       fntuParam->SetBranchAddress("meansgn",&fFitPars[7]);
 283       fntuParam->SetBranchAddress("sigmasgn",&fFitPars[8]);
 284
 285       fntuParam->SetBranchAddress("intbkg1Err",&fFitPars[9]);
 286       fntuParam->SetBranchAddress("slope1Err",&nothing);
 287       fntuParam->SetBranchAddress("conc1Err",&nothing);
 288       fntuParam->SetBranchAddress("intGBErr",&fFitPars[10]);
 289       fntuParam->SetBranchAddress("meanGBErr",&fFitPars[11]);
 290       fntuParam->SetBranchAddress("sigmaGBErr",&fFitPars[12]);
 291       fntuParam->SetBranchAddress("intbkg2Err",&fFitPars[13]);
 292       fntuParam->SetBranchAddress("slope2Err",&nothing);
 293       fntuParam->SetBranchAddress("conc2Err",&nothing);
 294       fntuParam->SetBranchAddress("inttotErr",&fFitPars[15]);
 295       fntuParam->SetBranchAddress("slope3Err",&nothing);
 296       fntuParam->SetBranchAddress("conc3Err",&nothing);
 297       fntuParam->SetBranchAddress("intsgnErr",&fFitPars[15]);
 298       fntuParam->SetBranchAddress("meansgnErr",&fFitPars[16]);
 299       fntuParam->SetBranchAddress("sigmasgnErr",&fFitPars[17]);
 300
 301     }
 302     else{
 303       fntuParam->SetBranchAddress("intbkg1",&fFitPars[0]);
 304       fntuParam->SetBranchAddress("slope1",&fFitPars[1]);
 305       fntuParam->SetBranchAddress("conc1",&nothing);
 306       fntuParam->SetBranchAddress("intGB",&fFitPars[2]);
 307       fntuParam->SetBranchAddress("meanGB",&fFitPars[3]);
 308       fntuParam->SetBranchAddress("sigmaGB",&fFitPars[4]);
 309       fntuParam->SetBranchAddress("intbkg2",&fFitPars[5]);
 310       fntuParam->SetBranchAddress("slope2",&fFitPars[6]);
 311       fntuParam->SetBranchAddress("conc2",&nothing);
 312       fntuParam->SetBranchAddress("inttot",&fFitPars[7]);
 313       fntuParam->SetBranchAddress("slope3",&fFitPars[8]);
 314       fntuParam->SetBranchAddress("conc3",&nothing);
 315       fntuParam->SetBranchAddress("intsgn",&fFitPars[9]);
 316       fntuParam->SetBranchAddress("meansgn",&fFitPars[10]);
 317       fntuParam->SetBranchAddress("sigmasgn",&fFitPars[11]);
 318
 319       fntuParam->SetBranchAddress("intbkg1Err",&fFitPars[12]);
 320       fntuParam->SetBranchAddress("slope1Err",&fFitPars[13]);
 321       fntuParam->SetBranchAddress("conc1Err",&nothing);
 322       fntuParam->SetBranchAddress("intGBErr",&fFitPars[14]);
 323       fntuParam->SetBranchAddress("meanGBErr",&fFitPars[15]);
 324       fntuParam->SetBranchAddress("sigmaGBErr",&fFitPars[16]);
 325       fntuParam->SetBranchAddress("intbkg2Err",&fFitPars[17]);
 326       fntuParam->SetBranchAddress("slope2Err",&fFitPars[18]);
 327       fntuParam->SetBranchAddress("conc2Err",&nothing);
 328       fntuParam->SetBranchAddress("inttotErr",&fFitPars[19]);
 329       fntuParam->SetBranchAddress("slope3Err",&fFitPars[20]);
 330       fntuParam->SetBranchAddress("conc3Err",&nothing);
 331       fntuParam->SetBranchAddress("intsgnErr",&fFitPars[21]);
 332       fntuParam->SetBranchAddress("meansgnErr",&fFitPars[22]);
 333       fntuParam->SetBranchAddress("sigmasgnErr",&fFitPars[23]);
 334     }
 335
 336   }
 337   fntuParam->TTree::Fill();
 338 }
 339
 340 //_________________________________________________________________________
 341
 342 TNtuple* AliHFMassFitter::NtuParamOneShot(char *ntuname){
 343   // Create, fill and return ntuple with fit parameters
 344
 345   InitNtuParam(ntuname);
 346   FillNtuParam();
 347   return fntuParam;
 348 }
 349 //_________________________________________________________________________
 350
 351 void AliHFMassFitter::RebinMass(Int_t bingroup){
 352   // Rebin invariant mass histogram
 353
 354   if(bingroup<1){
 355     cout<<"Error! Cannot group "<<bingroup<<" bins\n";
 356     fNbin=fhistoInvMass->GetNbinsX();
 357     cout<<"Kept original number of bins: "<<fNbin<<endl;
 358   } else{
 359     fhistoInvMass->Rebin(bingroup);
 360     fNbin = fhistoInvMass->GetNbinsX();
 361     cout<<"New number of bins: "<<fNbin<<endl;
 362   }
 363
 364
 365 }
 366
 367 //************* fit
 368
 369 //___________________________________________________________________________
 370
 371 Double_t AliHFMassFitter::FitFunction4MassDistr (Double_t *x, Double_t *par){
 372   // Fit function for signal+background
 373
 374
 375   //exponential or linear fit
 376   //
 377   // par[0] = tot integral
 378   // par[1] = slope
 379   // par[2] = gaussian integral
 380   // par[3] = gaussian mean
 381   // par[4] = gaussian sigma
 382
 383   Double_t total,bkg=0,sgn=0;
 384
 385   if (ftypeOfFit4Bkg==0 || ftypeOfFit4Bkg==1) {
 386     if(ftypeOfFit4Sgn == 0) {
 387
 388       Double_t parbkg[2] = {par[0]-par[2], par[1]};
 389       bkg = FitFunction4Bkg(x,parbkg);
 390     }
 391     if(ftypeOfFit4Sgn == 1) {
 392       Double_t parbkg[5] = {par[2],par[3],ffactor*par[4],par[0]-2*par[2], par[1]};
 393       bkg = FitFunction4Bkg(x,parbkg);
 394     }
 395
 396     sgn = FitFunction4Sgn(x,&par[2]);
 397
 398   }
 399
 400   //polynomial fit
 401
 402     // par[0] = tot integral
 403     // par[1] = coef1
 404     // par[2] = coef2
 405     // par[3] = gaussian integral
 406     // par[4] = gaussian mean
 407     // par[5] = gaussian sigma
 408
 409   if (ftypeOfFit4Bkg==2) {
 410
 411     if(ftypeOfFit4Sgn == 0) {
 412
 413       Double_t parbkg[3] = {par[0]-par[3], par[1], par[2]};
 414       bkg = FitFunction4Bkg(x,parbkg);
 415     }
 416     if(ftypeOfFit4Sgn == 1) {
 417
 418       Double_t parbkg[6] = {par[3],par[4],ffactor*par[5],par[0]-2*par[3], par[1], par[2]};
 419       bkg = FitFunction4Bkg(x,parbkg);
 420     }
 421
 422     sgn = FitFunction4Sgn(x,&par[3]);
 423   }
 424
 425   if (ftypeOfFit4Bkg==3) {
 426
 427     if(ftypeOfFit4Sgn == 0) {
 428         bkg=FitFunction4Bkg(x,par);
 429         sgn=FitFunction4Sgn(x,&par[1]);
 430     }
 431     if(ftypeOfFit4Sgn == 1) {
 432       Double_t parbkg[4]={par[1],par[2],ffactor*par[3],par[0]};
 433       bkg=FitFunction4Bkg(x,parbkg);
 434       sgn=FitFunction4Sgn(x,&par[1]);
 435     }
 436   }
 437
 438   total = bkg + sgn;
 439
 440   return  total;
 441 }
 442
 443 //_________________________________________________________________________
 444 Double_t AliHFMassFitter::FitFunction4Sgn (Double_t *x, Double_t *par){
 445   // Fit function for the signal
 446
 447   //gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
 448   //Par:
 449   // * [0] = integralSgn
 450   // * [1] = mean
 451   // * [2] = sigma
 452   //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
 453
 454   return par[0]/TMath::Sqrt(2.*TMath::Pi())/par[2]*TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/2./par[2]/par[2]);
 455
 456 }
 457
 458 //__________________________________________________________________________
 459
 460 Double_t AliHFMassFitter::FitFunction4Bkg (Double_t *x, Double_t *par){
 461   // Fit function for the background
 462
 463   Double_t maxDeltaM = 4.*fSigmaSgn;
 464   if(fSideBands && TMath::Abs(x[0]-fMass) < maxDeltaM) {
 465     TF1::RejectPoint();
 466     return 0;
 467   }
 468   Int_t firstPar=0;
 469   Double_t gaus2=0,total=-1;
 470   if(ftypeOfFit4Sgn == 1){
 471     firstPar=3;
 472     //gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
 473     //Par:
 474     // * [0] = integralSgn
 475     // * [1] = mean
 476     // * [2] = sigma
 477     //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
 478     gaus2 = FitFunction4Sgn(x,par);
 479   }
 480
 481   switch (ftypeOfFit4Bkg){
 482   case 0:
 483     //exponential
 484     //exponential = A*exp(B*x) -> integral(exponential)=A/B*exp(B*x)](min,max)
 485     //-> A = B*integral/(exp(B*max)-exp(B*min)) where integral can be written
 486     //as integralTot- integralGaus (=par [2])
 487     //Par:
 488     // * [0] = integralBkg;
 489     // * [1] = B;
 490     //exponential = [1]*[0]/(exp([1]*max)-exp([1]*min))*exp([1]*x)
 491     total = par[0+firstPar]*par[1+firstPar]/(TMath::Exp(par[1+firstPar]*fmaxMass)-TMath::Exp(par[1+firstPar]*fminMass))*TMath::Exp(par[1+firstPar]*x[0]);
 492     break;
 493   case 1:
 494     //linear
 495     //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)
 496     // * [0] = integralBkg;
 497     // * [1] = b;
 498     total= par[0+firstPar]/(fmaxMass-fminMass)+par[1+firstPar]*(x[0]-0.5*(fmaxMass+fminMass));
 499     break;
 500   case 2:
 501     //polynomial
 502     //y=a+b*x+c*x**2 -> integral = a(max-min) + 1/2*b*(max^2-min^2) +
 503     //+ 1/3*c*(max^3-min^3) ->
 504     //a = (integral-1/2*b*(max^2-min^2)-1/3*c*(max^3-min^3))/(max-min)
 505     // * [0] = integralBkg;
 506     // * [1] = b;
 507     // * [2] = c;
 508     total = par[0+firstPar]/(fmaxMass-fminMass)+par[1]*(x[0]-0.5*(fmaxMass+fminMass))+par[2+firstPar]*(x[0]*x[0]-1/3.*(fmaxMass*fmaxMass*fmaxMass-fminMass*fminMass*fminMass)/(fmaxMass-fminMass));
 509     break;
 510   case 3:
 511     total=par[0+firstPar];
 512     break;
 513 //   default:
 514 //     Types of Fit Functions for Background:
 515 //     * 0 = exponential;
 516 //     * 1 = linear;
 517 //     * 2 = polynomial 2nd order
 518 //     * 3 = no background"<<endl;
 519
 520   }
 521   return total+gaus2;
 522 }
 523
 524 //__________________________________________________________________________
 525
 526 Bool_t AliHFMassFitter::MassFitter(Bool_t draw){
 527   // Main method of the class: performs the fit of the histogram
 528
 529   Int_t nFitPars=0; //total function's number of parameters
 530   switch (ftypeOfFit4Bkg){
 531   case 0:
 532     nFitPars=5; //3+2
 533     break;
 534   case 1:
 535     nFitPars=5; //3+2
 536     break;
 537   case 2:
 538     nFitPars=6; //3+3
 539     break;
 540   case 3:
 541     nFitPars=4; //3+1
 542     break;
 543   }
 544
 545   Int_t bkgPar = nFitPars-3; //background function's number of parameters
 546
 547   cout<<"nFitPars = "<<nFitPars<<"\tbkgPar = "<<bkgPar<<endl;
 548
 549   //function names
 550   TString bkgname="funcbkg";
 551   TString bkg1name="funcbkg1";
 552   TString massname="funcmass";
 553
 554   //Total integral
 555   Double_t totInt = fhistoInvMass->Integral("width");
 556
 557   fSideBands = kTRUE;
 558   Double_t width=fhistoInvMass->GetBinWidth(8);
 559   Int_t binleft,binright;
 560   if (fNbin==0) fNbin=fhistoInvMass->GetNbinsX();
 561   Double_t minHisto=fhistoInvMass->GetBinLowEdge(1);
 562   Double_t maxHisto=fhistoInvMass->GetBinLowEdge(fNbin)+width;
 563
 564   if(fMass-fminMass < 0) {
 565     cout<<"Left limit of range > mean: change left limit or initial mean value"<<endl;
 566     return kFALSE;
 567   }
 568
 569   binleft=(Int_t)((fMass-4.*fSigmaSgn-minHisto)/width);
 570   binright=(Int_t)((fMass+4.*fSigmaSgn-minHisto)/width);
 571
 572
 573   if((fminMass!=minHisto || fmaxMass!= maxHisto) && (fMass-4.*fSigmaSgn-fminMass) < 0){
 574       Double_t coeff = (fMass-fminMass)/fSigmaSgn;
 575       binleft=(Int_t)((fMass-0.5*coeff*fSigmaSgn-fminMass)/width);
 576       binright=(Int_t)((fMass+0.5*coeff*fSigmaSgn-fminMass)/width);
 577       cout<<"Changed number of sigma from 4 to "<<0.5*coeff<<" for the estimation of the side bands"<<endl;
 578     }
 579
 580   if (binleft==0) {
 581     cout<<"Error! Range too little";
 582     return kFALSE;
 583   }
 584   Int_t fbin=(Int_t)((fminMass-minHisto)/width)+1, lbin=(Int_t)((fmaxMass-minHisto)/width)+1;
 585   cout<<"Range: "<<fminMass<<", "<<fmaxMass<<endl;
 586   cout<<"Range in bin = "<<fbin<<" --> "<<lbin<<endl;
 587
 588   //sidebands integral - first approx (from histo)
 589   Double_t sideBandsInt=(Double_t)fhistoInvMass->Integral(1,binleft,"width") + (Double_t)fhistoInvMass->Integral(binright,fNbin,"width");
 590   cout<<"------nbin = "<<fNbin<<"\twidth = "<<width<<"\tbinleft = "<<binleft<<"\tbinright = "<<binright<<endl;
 591   cout<<"------sideBandsInt - first approx = "<<sideBandsInt<<endl;
 592   /*
 593   //side bands usando il range utente
 594 Double_t sideBandsInt=(Double_t)fhistoInvMass->Integral(fbin,binleft,"width") + (Double_t)fhistoInvMass->Integral(binright,lbin,"width");
 595   cout<<"------nbin = "<<fNbin<<"\twidth = "<<width<<"\tbinleft = "<<binleft<<"\tbinright = "<<binright<<endl;
 596   cout<<"------sideBandsInt - first approx = "<<sideBandsInt<<endl;
 597   */
 598   /*Fit Bkg*/
 599
 600   //TF1 *funcbkg = new TF1(bkgname.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,bkgPar,"AliHFMassFitter","FitFunction4Bkg");
 601   TF1 *funcbkg = new TF1(bkgname.Data(),this,&AliHFMassFitter::FitFunction4Bkg,minHisto,maxHisto,bkgPar,"AliHFMassFitter","FitFunction4Bkg");
 602   cout<<"Function name = "<<funcbkg->GetName()<<endl<<endl;
 603
 604   funcbkg->SetLineColor(2); //red
 605
 606   //first fit for bkg: approx bkgint
 607
 608   switch (ftypeOfFit4Bkg) {
 609   case 0: //gaus+expo
 610     funcbkg->SetParNames("BkgInt","Slope");
 611     funcbkg->SetParameters(sideBandsInt,-2.);
 612     break;
 613   case 1:
 614     funcbkg->SetParNames("BkgInt","Slope");
 615     funcbkg->SetParameters(sideBandsInt,-100.);
 616     break;
 617   case 2:
 618     funcbkg->SetParNames("BkgInt","Coef1","Coef2");
 619     funcbkg->SetParameters(sideBandsInt,-10.,5);
 620     break;
 621   case 3:
 622     if(ftypeOfFit4Sgn==0){
 623       funcbkg->SetParNames("Const");
 624       funcbkg->SetParameter(0,0.);
 625       funcbkg->FixParameter(0,0.);
 626     }
 627     break;
 628   default:
 629     cout<<"Wrong choise of ftypeOfFit4Bkg ("<<ftypeOfFit4Bkg<<")"<<endl;
 630     return kFALSE;
 631     break;
 632   }
 633   cout<<"\nBACKGROUND FIT - only combinatorial"<<endl;
 634   Int_t ftypeOfFit4SgnBkp=ftypeOfFit4Sgn;
 635
 636   Double_t intbkg1=0,slope1=0,conc1=0;
 637   //if only signal and reflection: skip
 638   if (!(ftypeOfFit4Bkg==3 && ftypeOfFit4Sgn==1)) {
 639     ftypeOfFit4Sgn=0;
 640     fhistoInvMass->Fit(bkgname.Data(),"R,L,E,0");
 641
 642     for(Int_t i=0;i<bkgPar;i++){
 643       fFitPars[i]=funcbkg->GetParameter(i);
 644       //cout<<i<<"\t"<<funcbkg->GetParameter(i)<<"\t";
 645       fFitPars[nFitPars+2*bkgPar+3+i]= funcbkg->GetParError(i);
 646       //cout<<nFitPars+2*bkgPar+3+i<<"\t"<< funcbkg->GetParError(i)<<endl;
 647     }
 648
 649     //intbkg1 = funcbkg->GetParameter(0);
 650     funcbkg->SetRange(fminMass,fmaxMass);
 651     intbkg1 = funcbkg->Integral(fminMass,fmaxMass);
 652     if(ftypeOfFit4Bkg!=3) slope1 = funcbkg->GetParameter(1);
 653     if(ftypeOfFit4Bkg==2) conc1 = funcbkg->GetParameter(2);
 654     cout<<"Primo fit: \nintbkg1 = "<<intbkg1<<"\t(Compare with par0 = "<<funcbkg->GetParameter(0)<<")\nslope1= "<<slope1<<"\nconc1 = "<<conc1<<endl;
 655   }
 656   else cout<<"\t\t//"<<endl;
 657
 658   ftypeOfFit4Sgn=ftypeOfFit4SgnBkp;
 659   TF1 *funcbkg1=0;
 660   if (ftypeOfFit4Sgn == 1) {
 661     cout<<"\nBACKGROUND FIT WITH REFLECTION"<<endl;
 662     bkgPar+=3;
 663
 664     cout<<"nFitPars = "<<nFitPars<<"\tbkgPar = "<<bkgPar<<endl;
 665
 666     funcbkg1 = new TF1(bkg1name.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,bkgPar,"AliHFMassFitter","FitFunction4Bkg");
 667     cout<<"Function name = "<<funcbkg1->GetName()<<endl;
 668
 669     funcbkg1->SetLineColor(2); //red
 670
 671     if(ftypeOfFit4Bkg==2){
 672       cout<<"*** Polynomial Fit ***"<<endl;
 673       funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Coef1","Coef2");
 674       funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1,conc1);
 675 //     cout<<"Parameters set to: "<<0.5*(totInt-intbkg1)<<"\t"<<fMass<<"\t"<<ffactor*fSigmaSgn<<"\t"<<intbkg1<<"\t"<<slope1<<"\t"<<conc1<<"\t"<<endl;
 676 //     cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<bkgPar<<"\tgsidebands = "<<fSideBands<<endl;
 677     } else{
 678       if(ftypeOfFit4Bkg==3) //no background: gaus sign+ gaus broadened
 679         {
 680           cout<<"*** No background Fit ***"<<endl;
 681           funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","Const");
 682           funcbkg1->SetParameters(0.5*totInt,fMass,ffactor*fSigmaSgn,0.);
 683           funcbkg1->FixParameter(3,0.);
 684         } else{ //expo or linear
 685           if(ftypeOfFit4Bkg==0) cout<<"*** Exponential Fit ***"<<endl;
 686           if(ftypeOfFit4Bkg==1) cout<<"*** Linear Fit ***"<<endl;
 687           funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Slope");
 688           funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1);
 689         }
 690     }
 691     fhistoInvMass->Fit(bkg1name.Data(),"R,L,E,+,0");
 692
 693     for(Int_t i=0;i<bkgPar;i++){
 694       fFitPars[bkgPar-3+i]=funcbkg1->GetParameter(i);
 695       //cout<<bkgPar-3+i<<"\t"<<funcbkg1->GetParameter(i);
 696       fFitPars[nFitPars+3*bkgPar-6+i]= funcbkg1->GetParError(i);
 697       //cout<<"\t"<<nFitPars+3*bkgPar-6+i<<"\t"<<funcbkg1->GetParError(i)<<endl;
 698     }
 699
 700     intbkg1=funcbkg1->GetParameter(3);
 701     if(ftypeOfFit4Bkg!=3) slope1 = funcbkg1->GetParameter(4);
 702     if(ftypeOfFit4Bkg==2) conc1 = funcbkg1->GetParameter(5);
 703
 704   } else {
 705     bkgPar+=3;
 706
 707     for(Int_t i=0;i<3;i++){
 708       fFitPars[bkgPar-3+i]=0.;
 709       cout<<bkgPar-3+i<<"\t"<<0.<<"\t";
 710       fFitPars[nFitPars+3*bkgPar-6+i]= 0.;
 711       cout<<nFitPars+3*bkgPar-6+i<<"\t"<<0.<<endl;
 712     }
 713
 714     for(Int_t i=0;i<bkgPar-3;i++){
 715       fFitPars[bkgPar+i]=funcbkg->GetParameter(i);
 716       cout<<bkgPar+i<<"\t"<<funcbkg->GetParameter(i)<<"\t";
 717       fFitPars[nFitPars+3*bkgPar-3+i]= funcbkg->GetParError(i);
 718       cout<<nFitPars+3*bkgPar-3+i<<"\t"<< funcbkg->GetParError(i)<<endl;
 719     }
 720
 721
 722   }
 723
 724   //sidebands integral - second approx (from fit)
 725   fSideBands = kFALSE;
 726   Double_t bkgInt;
 727
 728   if(ftypeOfFit4Sgn == 1) bkgInt=funcbkg1->Integral(fminMass,fmaxMass);
 729   else bkgInt=funcbkg->Integral(fminMass,fmaxMass);
 730   cout<<"------BkgInt(Fit) = "<<bkgInt<<endl;
 731
 732   //Signal integral - first approx
 733   Double_t sgnInt;
 734   sgnInt = totInt-bkgInt;
 735   cout<<"------TotInt = "<<totInt<<"\tsgnInt = "<<sgnInt<<endl;
 736
 737   /*Fit All Mass distribution with exponential + gaussian (+gaussiam braodened) */
 738   TF1 *funcmass = new TF1(massname.Data(),this,&AliHFMassFitter::FitFunction4MassDistr,fminMass,fmaxMass,nFitPars,"AliHFMassFitter","FitFunction4MassDistr");
 739   cout<<"Function name = "<<funcmass->GetName()<<endl<<endl;
 740   funcmass->SetLineColor(4); //blue
 741
 742   //Set parameters
 743   cout<<"\nTOTAL FIT"<<endl;
 744
 745   if(nFitPars==5){
 746     funcmass->SetParNames("TotInt","Slope","SgnInt","Mean","Sigma");
 747     funcmass->SetParameters(totInt,slope1,sgnInt,fMass,fSigmaSgn);
 748     //cout<<"Parameters set to: "<<totInt<<"\t"<<slope1<<"\t"<<sgnInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl;
 749 //     cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<nFitPars<<"\tgsidebands = "<<fSideBands<<endl;
 750     funcmass->FixParameter(0,totInt);
 751   }
 752   if (nFitPars==6){
 753     funcmass->SetParNames("TotInt","Coef1","Coef2","SgnInt","Mean","Sigma");
 754     funcmass->SetParameters(totInt,slope1,conc1,sgnInt,fMass,fSigmaSgn);
 755 //     cout<<"Parameters set to: "<<totInt<<"\t"<<slope1<<"\t"<<conc1<<"\t"<<sgnInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl;
 756 //     cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<nFitPars<<"\tgsidebands = "<<fSideBands<<endl;
 757     funcmass->FixParameter(0,totInt);
 758   }
 759   if(nFitPars==4){
 760     funcmass->SetParNames("Const","SgnInt","Mean","Sigma");
 761     if(ftypeOfFit4Sgn == 1) funcmass->SetParameters(0.,0.5*totInt,fMass,fSigmaSgn);
 762     else funcmass->SetParameters(0.,totInt,fMass,fSigmaSgn);
 763     funcmass->FixParameter(0,0.);
 764     //cout<<"Parameters set to: "<<0.5*totInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl;
 765     //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<nFitPars<<"\tgsidebands = "<<fSideBands<<endl;
 766
 767   }
 768
 769   fhistoInvMass->Fit(massname.Data(),"R,L,E,+,0");
 770   cout<<"fit done"<<endl;
 771   //reset value of fMass and fSigmaSgn to those found from fit
 772   fMass=funcmass->GetParameter(nFitPars-2);
 773   fSigmaSgn=funcmass->GetParameter(nFitPars-1);
 774
 775   for(Int_t i=0;i<nFitPars;i++){
 776     fFitPars[i+2*bkgPar-3]=funcmass->GetParameter(i);
 777     fFitPars[nFitPars+4*bkgPar-6+i]= funcmass->GetParError(i);
 778     //cout<<i+2*bkgPar-3<<"\t"<<funcmass->GetParameter(i)<<"\t\t"<<nFitPars+4*bkgPar-6+i<<"\t"<<funcmass->GetParError(i)<<endl;
 779   }
 780   /*
 781   //check: cout parameters
 782   for(Int_t i=0;i<2*(nFitPars+2*bkgPar-3);i++){
 783     cout<<i<<"\t"<<fFitPars[i]<<endl;
 784     }
 785   */
 786
 787   if(draw){
 788     TCanvas *canvas=new TCanvas(fhistoInvMass->GetName(),fhistoInvMass->GetName());
 789     TH1F *fhistocopy=new TH1F(*fhistoInvMass);
 790     canvas->cd();
 791     fhistocopy->DrawClone();
 792     if(ftypeOfFit4Sgn == 1) funcbkg1->DrawClone("sames");
 793     else funcbkg->DrawClone("sames");
 794     funcmass->DrawClone("sames");
 795     cout<<"Drawn"<<endl;
 796   }
 797
 798
 799   //increase counter of number of fits done
 800   fcounter++;
 801
 802   if (ftypeOfFit4Sgn == 1) delete funcbkg1;
 803   delete funcbkg;
 804   delete funcmass;
 805   return kTRUE;
 806 }
 807 //_________________________________________________________________________
 808 Double_t AliHFMassFitter::GetChiSquare() const{
 809   TF1 *funcmass=(TF1*)fhistoInvMass->GetFunction("funcmass");
 810   return funcmass->GetChisquare();
 811 }
 812
 813 //*********output
 814
 815 //_________________________________________________________________________
 816 void  AliHFMassFitter::GetFitPars(Float_t *vector) const {
 817   // Return fit parameters
 818
 819   for(Int_t i=0;i<fParsSize;i++){
 820     vector[i]=fFitPars[i];
 821   }
 822 }
 823 //_________________________________________________________________________
 824
 825 TH1F* AliHFMassFitter::GetHistoClone() const{
 826   TH1F* hout=(TH1F*)fhistoInvMass->Clone(fhistoInvMass->GetName());
 827   return hout;
 828 }
 829 //_________________________________________________________________________
 830
 831 void AliHFMassFitter::WriteHisto(TString path) {
 832   if (fcounter == 0) {
 833     cout<<"Use MassFitter method before WriteHisto"<<endl;
 834     return;
 835   }
 836
 837   TH1F* hget=(TH1F*)fhistoInvMass->Clone();
 838
 839   TString bkgname = "funcbkg";
 840   Int_t np=-99;
 841   switch (ftypeOfFit4Bkg){
 842   case 0:
 843     np=2;
 844     break;
 845   case 1:
 846     np=2;
 847     break;
 848   case 2:
 849     np=3;
 850     break;
 851   case 3:
 852     np=1;
 853     break;
 854   }
 855   if (ftypeOfFit4Sgn == 1) {
 856     bkgname += 1;
 857     np+=3;
 858   }
 859   TList *hlist=hget->GetListOfFunctions();
 860   hlist->ls();
 861   TF1 *b=(TF1*)hlist->FindObject(bkgname.Data());
 862   cout<< b->GetParameter(0)<<endl;
 863   //cout<<"WRITEHISTO np = "<<np<<"\n"<<bkgname<<endl;
 864   bkgname += "FullRange";
 865   TF1 *bwrite=new TF1(bkgname.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,np,"AliHFMassFitter","FitFunction4Bkg");
 866   cout<<bwrite->GetName()<<endl;
 867   for(Int_t i=0;i<np;i++){
 868     //cout<<i<<" di "<<np<<endl;
 869     bwrite->SetParName(i,b->GetParName(i));
 870     bwrite->SetParameter(i,b->GetParameter(i));
 871   }
 872
 873   hlist->Add(bwrite);
 874   path += "HFMassFitterOutput.root";
 875   TFile *output;
 876
 877   if (fcounter == 1) output = new TFile(path.Data(),"recreate");
 878   else output = new TFile(path.Data(),"update");
 879   output->cd();
 880   hget->Write();
 881   output->Close();
 882
 883   cout<<fcounter<<" "<<hget->GetName()<<" written in "<<path<<endl;
 884
 885   delete bwrite;
 886   delete output;
 887
 888 }
 889
 890 //_________________________________________________________________________
 891
 892 void AliHFMassFitter::WriteNtuple(TString path) const{
 893   TNtuple* nget=(TNtuple*)fntuParam->Clone();
 894   path += "HFMassFitterOutput.root";
 895   TFile *output = new TFile(path.Data(),"update");
 896   output->cd();
 897   nget->Write();
 898   output->Close();
 899   cout<<nget->GetName()<<" written in "<<path<<endl;
 900   delete nget;
 901   delete output;
 902 }
 903
 904
 905 //************ significance
 906
 907 //_________________________________________________________________________
 908
 909 void AliHFMassFitter::Signal(Double_t nOfSigma,Double_t &signal,Double_t &errsignal) const {
 910   // Return signal integral in mean+- n sigma
 911
 912
 913   //functions names
 914   TString bkgname="funcbkg";
 915   TString bkg1name="funcbkg1";
 916   TString massname="funcmass";
 917
 918   TF1 *funcmass=fhistoInvMass->GetFunction(massname.Data());
 919   Double_t mean=0,sigma=1;
 920   Double_t intS,intSerr;
 921
 922   if(ftypeOfFit4Bkg == 2) { //polynomial
 923     mean=funcmass->GetParameter(4); //mean
 924     sigma=funcmass->GetParameter(5); //sigma
 925     intS=fFitPars[12];
 926     intSerr=fFitPars[27];
 927   } else if(ftypeOfFit4Bkg == 3){ //no background
 928     mean=funcmass->GetParameter(2); //mean
 929     sigma=funcmass->GetParameter(3); //sigma
 930     intS=fFitPars[6];
 931     intSerr=fFitPars[15];
 932   } else { //expo or linear
 933     mean=funcmass->GetParameter(3); //mean
 934     sigma=funcmass->GetParameter(4); //sigma
 935     intS=fFitPars[9];
 936     intSerr=fFitPars[21];
 937   }
 938
 939   TF1 *funcbkg=0;
 940   if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
 941   else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
 942   Double_t min=mean-nOfSigma*sigma;
 943   Double_t max=mean+nOfSigma*sigma;
 944   if(ftypeOfFit4Bkg == 3 && ftypeOfFit4Sgn == 0) signal=funcmass->Integral(min,max)/(Double_t)fhistoInvMass->GetBinWidth(2);
 945   else signal=(funcmass->Integral(min,max)-funcbkg->Integral(min,max))/(Double_t)fhistoInvMass->GetBinWidth(2);
 946   errsignal=intSerr/intS*signal; // assume relative error is the same as for total integral
 947
 948   return;
 949 }
 950 //_________________________________________________________________________
 951
 952 void AliHFMassFitter::Background(Double_t nOfSigma,Double_t &background,Double_t &errbackground) const {
 953   // Return background integral in mean+- n sigma
 954
 955   //functions names
 956   TString bkgname="funcbkg";
 957   TString bkg1name="funcbkg1";
 958   TString massname="funcmass";
 959
 960   TF1 *funcmass=fhistoInvMass->GetFunction(massname.Data());
 961   Double_t mean=0,sigma=1;
 962   Double_t intB,intBerr,err;
 963
 964   if(ftypeOfFit4Bkg == 2) { //polynomial
 965     mean=funcmass->GetParameter(4); //mean
 966     sigma=funcmass->GetParameter(5); //sigma
 967     intB=fFitPars[9]-fFitPars[12];
 968     intBerr=fFitPars[27];
 969   } else if(ftypeOfFit4Bkg == 3){ //no background
 970     mean=funcmass->GetParameter(2); //mean
 971     sigma=funcmass->GetParameter(3); //sigma
 972     if (ftypeOfFit4Sgn == 1){ //reflection
 973       intB=fFitPars[1];
 974       intBerr=fFitPars[9];
 975     } else {
 976       intB=-1;intBerr=-1;
 977       err=0;
 978     }
 979   } else { //expo or linear
 980     mean=funcmass->GetParameter(3); //mean
 981     sigma=funcmass->GetParameter(4); //sigma
 982     intB=fFitPars[7]-fFitPars[9];
 983     intBerr=fFitPars[21];
 984   }
 985
 986   TF1 *funcbkg=0;
 987
 988   if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
 989   else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
 990   Double_t min=mean-nOfSigma*sigma;
 991   Double_t max=mean+nOfSigma*sigma;
 992   if (ftypeOfFit4Bkg == 3 && ftypeOfFit4Sgn == 0) {
 993     background = 0;
 994     errbackground = 0;
 995   }
 996   else{
 997     background=funcbkg->Integral(min,max)/(Double_t)fhistoInvMass->GetBinWidth(2);
 998     errbackground=intBerr/intB*background;
 999   }
1000   return;
1001
1002 }
1003
1004 //__________________________________________________________________________
1005
1006 void AliHFMassFitter::Significance(Double_t nOfSigma,Double_t &significance,Double_t &errsignificance)  const {
1007   // Return significance in mean+- n sigma
1008
1009   Double_t signal,errsignal,background,errbackground;
1010   Signal(nOfSigma,signal,errsignal);
1011   Background(nOfSigma,background,errbackground);
1012
1013   significance =  signal/TMath::Sqrt(signal+background);
1014   errsignificance = TMath::Sqrt(significance*significance/(signal+background)/(signal+background)*(1/4.*errsignal*errsignal+errbackground*errbackground)+significance*significance/signal/signal*errsignal*errsignal);
1015
1016   return;
1017 }
1018
1019 //__________________________________________________________________________
1020