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   fSideBandl(0),
  49   fSideBandr(0),
  50   fcounter(0)
  51 {
  52   // default constructor
  53
  54   cout<<"Default constructor"<<endl;
  55 }
  56
  57 //___________________________________________________________________________
  58
  59 AliHFMassFitter::AliHFMassFitter (TH1F *histoToFit, Double_t minvalue, Double_t maxvalue, Int_t rebin,Int_t fittypeb,Int_t fittypes):
  60  TNamed(),
  61  fhistoInvMass(0),
  62  fminMass(0),
  63  fmaxMass(0),
  64  fNbin(0),
  65  fParsSize(1),
  66  fFitPars(0),
  67  fWithBkg(0),
  68  ftypeOfFit4Bkg(0),
  69  ftypeOfFit4Sgn(0),
  70  ffactor(1),
  71  fntuParam(0),
  72  fMass(1.85),
  73  fSigmaSgn(0.012),
  74  fSideBands(0),
  75  fSideBandl(0),
  76  fSideBandr(0),
  77  fcounter(0)
  78 {
  79   // standard constructor
  80
  81   fhistoInvMass= (TH1F*)histoToFit->Clone("fhistoInvMass");
  82   fminMass=minvalue;
  83   fmaxMass=maxvalue;
  84   if(rebin!=1) RebinMass(rebin);
  85   else fNbin=(Int_t)fhistoInvMass->GetNbinsX();
  86   ftypeOfFit4Bkg=fittypeb;
  87   ftypeOfFit4Sgn=fittypes;
  88   if(ftypeOfFit4Bkg!=0 && ftypeOfFit4Bkg!=1 && ftypeOfFit4Bkg!=2) fWithBkg=kFALSE;
  89   else fWithBkg=kTRUE;
  90   if (!fWithBkg) cout<<"Fit Histogram of Signal only"<<endl;
  91   else  cout<<"Type of fit For Background = "<<ftypeOfFit4Bkg<<endl;
  92
  93   ComputeParSize();
  94   fFitPars=new Float_t[fParsSize];
  95 }
  96
  97
  98 AliHFMassFitter::AliHFMassFitter(const AliHFMassFitter &mfit):
  99   TNamed(),
 100   fhistoInvMass(mfit.fhistoInvMass),
 101   fminMass(mfit.fminMass),
 102   fmaxMass(mfit.fmaxMass),
 103   fNbin(mfit.fNbin),
 104   fParsSize(mfit.fParsSize),
 105   fFitPars(0),
 106   fWithBkg(mfit.fWithBkg),
 107   ftypeOfFit4Bkg(mfit.ftypeOfFit4Bkg),
 108   ftypeOfFit4Sgn(mfit.ftypeOfFit4Sgn),
 109   ffactor(mfit.ffactor),
 110   fntuParam(mfit.fntuParam),
 111   fMass(mfit.fMass),
 112   fSigmaSgn(mfit.fSigmaSgn),
 113   fSideBands(mfit.fSideBands),
 114   fSideBandl(mfit.fSideBandl),
 115   fSideBandr(mfit.fSideBandr),
 116   fcounter(mfit.fcounter)
 117
 118 {
 119   //copy constructor
 120
 121   if(mfit.fParsSize > 0){
 122     fFitPars=new Float_t[fParsSize];
 123     memcpy(fFitPars,mfit.fFitPars,mfit.fParsSize*sizeof(Float_t));
 124   }
 125   //for(Int_t i=0;i<fParsSize;i++) fFitPars[i]=mfit.fFitPars[i];
 126 }
 127
 128 //_________________________________________________________________________
 129
 130 AliHFMassFitter::~AliHFMassFitter() {
 131   cout<<"AliHFMassFitter destructor called"<<endl;
 132   if(fhistoInvMass) {
 133     cout<<"deleting histogram..."<<endl;
 134     delete fhistoInvMass;
 135     fhistoInvMass=NULL;
 136   }
 137   if(fntuParam){
 138     cout<<"deleting ntuple..."<<endl;
 139     delete fntuParam;
 140
 141     fntuParam=NULL;
 142   }
 143
 144   if(fFitPars) {
 145     delete[] fFitPars;
 146     cout<<"deleting parameter array..."<<endl;
 147     fFitPars=NULL;
 148   }
 149
 150   fcounter = 0;
 151 }
 152
 153 //_________________________________________________________________________
 154
 155 AliHFMassFitter& AliHFMassFitter::operator=(const AliHFMassFitter &mfit){
 156
 157   //assignment operator
 158
 159   if(&mfit == this) return *this;
 160   fhistoInvMass= mfit.fhistoInvMass;
 161   fminMass= mfit.fminMass;
 162   fmaxMass= mfit.fmaxMass;
 163   fNbin= mfit.fNbin;
 164   fParsSize= mfit.fParsSize;
 165   fWithBkg= mfit.fWithBkg;
 166   ftypeOfFit4Bkg= mfit.ftypeOfFit4Bkg;
 167   ftypeOfFit4Sgn= mfit.ftypeOfFit4Sgn;
 168   ffactor= mfit.ffactor;
 169   fntuParam= mfit.fntuParam;
 170   fMass= mfit.fMass;
 171   fSigmaSgn= mfit.fSigmaSgn;
 172   fSideBands= mfit.fSideBands;
 173   fSideBandl= mfit.fSideBandl;
 174   fSideBandr= mfit.fSideBandr;
 175   fcounter= mfit.fcounter;
 176   if(mfit.fParsSize > 0){
 177     if(fFitPars) {
 178       delete[] fFitPars;
 179       fFitPars=NULL;
 180     }
 181     fFitPars=new Float_t[fParsSize];
 182     memcpy(fFitPars,mfit.fFitPars,mfit.fParsSize*sizeof(Float_t));
 183   }
 184 // fFitPars=new Float_t[fParsSize];
 185 //   for(Int_t i=0;i<fParsSize;i++) fFitPars[i]=mfit.fFitPars[i];
 186
 187   return *this;
 188 }
 189
 190 //************ tools & settings
 191
 192 //__________________________________________________________________________
 193
 194 void AliHFMassFitter::ComputeParSize() {
 195   switch (ftypeOfFit4Bkg) {//npar backround func
 196   case 0:
 197     fParsSize = 2*3;
 198     break;
 199   case 1:
 200     fParsSize = 2*3;
 201     break;
 202   case 2:
 203     fParsSize = 3*3;
 204     break;
 205   case 3:
 206     fParsSize = 1*3;
 207     break;
 208   default:
 209     cout<<"Error in computing fParsSize: check ftypeOfFit4Bkg"<<endl;
 210     break;
 211   }
 212
 213   fParsSize += 3; // npar refl
 214   fParsSize += 3; // npar signal gaus
 215
 216   fParsSize*=2;   // add errors
 217   cout<<"Parameters array size "<<fParsSize<<endl;
 218 }
 219
 220 //___________________________________________________________________________
 221 void AliHFMassFitter::SetHisto(TH1F *histoToFit){
 222   //fhistoInvMass = (TH1F*)histoToFit->Clone();
 223   fhistoInvMass = new TH1F(*histoToFit);
 224   cout<<"SetHisto pointer "<<fhistoInvMass<<endl;
 225 }
 226
 227 //___________________________________________________________________________
 228
 229   void AliHFMassFitter::SetType(Int_t fittypeb, Int_t fittypes) {
 230     ftypeOfFit4Bkg = fittypeb;
 231     ftypeOfFit4Sgn = fittypes;
 232     ComputeParSize();
 233     fFitPars = new Float_t[fParsSize];
 234
 235 }
 236
 237 //___________________________________________________________________________
 238
 239 void AliHFMassFitter::Reset() {
 240   cout<<"Reset called: delete histogram, set mean value to 1.85 and sigma to 0.012"<<endl;
 241   fMass=1.85;
 242   fSigmaSgn=0.012;
 243   cout<<"Reset "<<fhistoInvMass<<endl;
 244   if(fhistoInvMass) {
 245     cout<<"esiste"<<endl;
 246     //delete fhistoInvMass;
 247     fhistoInvMass=NULL;
 248     cout<<fhistoInvMass<<endl;
 249   }
 250   else cout<<"histogram doesn't exist, do not delete"<<endl;
 251
 252 }
 253
 254 //_________________________________________________________________________
 255
 256 void AliHFMassFitter::InitNtuParam(char *ntuname) {
 257   // Create ntuple to keep fit parameters
 258   fntuParam=0;
 259   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");
 260
 261 }
 262
 263 //_________________________________________________________________________
 264
 265 void AliHFMassFitter::FillNtuParam() {
 266   // Fill ntuple with fit parameters
 267
 268   Float_t nothing=0.;
 269
 270   if (ftypeOfFit4Bkg==2) {
 271       fntuParam->SetBranchAddress("intbkg1",&fFitPars[0]);
 272       fntuParam->SetBranchAddress("slope1",&fFitPars[1]);
 273       fntuParam->SetBranchAddress("conc1",&fFitPars[2]);
 274       fntuParam->SetBranchAddress("intGB",&fFitPars[3]);
 275       fntuParam->SetBranchAddress("meanGB",&fFitPars[4]);
 276       fntuParam->SetBranchAddress("sigmaGB",&fFitPars[5]);
 277       fntuParam->SetBranchAddress("intbkg2",&fFitPars[6]);
 278       fntuParam->SetBranchAddress("slope2",&fFitPars[7]);
 279       fntuParam->SetBranchAddress("conc2",&fFitPars[8]);
 280       fntuParam->SetBranchAddress("inttot",&fFitPars[9]);
 281       fntuParam->SetBranchAddress("slope3",&fFitPars[10]);
 282       fntuParam->SetBranchAddress("conc3",&fFitPars[11]);
 283       fntuParam->SetBranchAddress("intsgn",&fFitPars[12]);
 284       fntuParam->SetBranchAddress("meansgn",&fFitPars[13]);
 285       fntuParam->SetBranchAddress("sigmasgn",&fFitPars[14]);
 286
 287       fntuParam->SetBranchAddress("intbkg1Err",&fFitPars[15]);
 288       fntuParam->SetBranchAddress("slope1Err",&fFitPars[16]);
 289       fntuParam->SetBranchAddress("conc1Err",&fFitPars[17]);
 290       fntuParam->SetBranchAddress("intGBErr",&fFitPars[18]);
 291       fntuParam->SetBranchAddress("meanGBErr",&fFitPars[19]);
 292       fntuParam->SetBranchAddress("sigmaGBErr",&fFitPars[20]);
 293       fntuParam->SetBranchAddress("intbkg2Err",&fFitPars[21]);
 294       fntuParam->SetBranchAddress("slope2Err",&fFitPars[22]);
 295       fntuParam->SetBranchAddress("conc2Err",&fFitPars[23]);
 296       fntuParam->SetBranchAddress("inttotErr",&fFitPars[24]);
 297       fntuParam->SetBranchAddress("slope3Err",&fFitPars[25]);
 298       fntuParam->SetBranchAddress("conc3Err",&fFitPars[26]);
 299       fntuParam->SetBranchAddress("intsgnErr",&fFitPars[27]);
 300       fntuParam->SetBranchAddress("meansgnErr",&fFitPars[28]);
 301       fntuParam->SetBranchAddress("sigmasgnErr",&fFitPars[29]);
 302
 303   } else {
 304
 305     if(ftypeOfFit4Bkg==3){
 306       fntuParam->SetBranchAddress("intbkg1",&fFitPars[0]);
 307       fntuParam->SetBranchAddress("slope1",&nothing);
 308       fntuParam->SetBranchAddress("conc1",&nothing);
 309       fntuParam->SetBranchAddress("intGB",&fFitPars[1]);
 310       fntuParam->SetBranchAddress("meanGB",&fFitPars[2]);
 311       fntuParam->SetBranchAddress("sigmaGB",&fFitPars[3]);
 312       fntuParam->SetBranchAddress("intbkg2",&fFitPars[4]);
 313       fntuParam->SetBranchAddress("slope2",&nothing);
 314       fntuParam->SetBranchAddress("conc2",&nothing);
 315       fntuParam->SetBranchAddress("inttot",&fFitPars[6]);
 316       fntuParam->SetBranchAddress("slope3",&nothing);
 317       fntuParam->SetBranchAddress("conc3",&nothing);
 318       fntuParam->SetBranchAddress("intsgn",&fFitPars[6]);
 319       fntuParam->SetBranchAddress("meansgn",&fFitPars[7]);
 320       fntuParam->SetBranchAddress("sigmasgn",&fFitPars[8]);
 321
 322       fntuParam->SetBranchAddress("intbkg1Err",&fFitPars[9]);
 323       fntuParam->SetBranchAddress("slope1Err",&nothing);
 324       fntuParam->SetBranchAddress("conc1Err",&nothing);
 325       fntuParam->SetBranchAddress("intGBErr",&fFitPars[10]);
 326       fntuParam->SetBranchAddress("meanGBErr",&fFitPars[11]);
 327       fntuParam->SetBranchAddress("sigmaGBErr",&fFitPars[12]);
 328       fntuParam->SetBranchAddress("intbkg2Err",&fFitPars[13]);
 329       fntuParam->SetBranchAddress("slope2Err",&nothing);
 330       fntuParam->SetBranchAddress("conc2Err",&nothing);
 331       fntuParam->SetBranchAddress("inttotErr",&fFitPars[15]);
 332       fntuParam->SetBranchAddress("slope3Err",&nothing);
 333       fntuParam->SetBranchAddress("conc3Err",&nothing);
 334       fntuParam->SetBranchAddress("intsgnErr",&fFitPars[15]);
 335       fntuParam->SetBranchAddress("meansgnErr",&fFitPars[16]);
 336       fntuParam->SetBranchAddress("sigmasgnErr",&fFitPars[17]);
 337
 338     }
 339     else{
 340       fntuParam->SetBranchAddress("intbkg1",&fFitPars[0]);
 341       fntuParam->SetBranchAddress("slope1",&fFitPars[1]);
 342       fntuParam->SetBranchAddress("conc1",&nothing);
 343       fntuParam->SetBranchAddress("intGB",&fFitPars[2]);
 344       fntuParam->SetBranchAddress("meanGB",&fFitPars[3]);
 345       fntuParam->SetBranchAddress("sigmaGB",&fFitPars[4]);
 346       fntuParam->SetBranchAddress("intbkg2",&fFitPars[5]);
 347       fntuParam->SetBranchAddress("slope2",&fFitPars[6]);
 348       fntuParam->SetBranchAddress("conc2",&nothing);
 349       fntuParam->SetBranchAddress("inttot",&fFitPars[7]);
 350       fntuParam->SetBranchAddress("slope3",&fFitPars[8]);
 351       fntuParam->SetBranchAddress("conc3",&nothing);
 352       fntuParam->SetBranchAddress("intsgn",&fFitPars[9]);
 353       fntuParam->SetBranchAddress("meansgn",&fFitPars[10]);
 354       fntuParam->SetBranchAddress("sigmasgn",&fFitPars[11]);
 355
 356       fntuParam->SetBranchAddress("intbkg1Err",&fFitPars[12]);
 357       fntuParam->SetBranchAddress("slope1Err",&fFitPars[13]);
 358       fntuParam->SetBranchAddress("conc1Err",&nothing);
 359       fntuParam->SetBranchAddress("intGBErr",&fFitPars[14]);
 360       fntuParam->SetBranchAddress("meanGBErr",&fFitPars[15]);
 361       fntuParam->SetBranchAddress("sigmaGBErr",&fFitPars[16]);
 362       fntuParam->SetBranchAddress("intbkg2Err",&fFitPars[17]);
 363       fntuParam->SetBranchAddress("slope2Err",&fFitPars[18]);
 364       fntuParam->SetBranchAddress("conc2Err",&nothing);
 365       fntuParam->SetBranchAddress("inttotErr",&fFitPars[19]);
 366       fntuParam->SetBranchAddress("slope3Err",&fFitPars[20]);
 367       fntuParam->SetBranchAddress("conc3Err",&nothing);
 368       fntuParam->SetBranchAddress("intsgnErr",&fFitPars[21]);
 369       fntuParam->SetBranchAddress("meansgnErr",&fFitPars[22]);
 370       fntuParam->SetBranchAddress("sigmasgnErr",&fFitPars[23]);
 371     }
 372
 373   }
 374   fntuParam->TTree::Fill();
 375 }
 376
 377 //_________________________________________________________________________
 378
 379 TNtuple* AliHFMassFitter::NtuParamOneShot(char *ntuname){
 380   // Create, fill and return ntuple with fit parameters
 381
 382   InitNtuParam(ntuname);
 383   FillNtuParam();
 384   return fntuParam;
 385 }
 386 //_________________________________________________________________________
 387
 388 void AliHFMassFitter::RebinMass(Int_t bingroup){
 389   // Rebin invariant mass histogram
 390
 391   if(bingroup<1){
 392     cout<<"Error! Cannot group "<<bingroup<<" bins\n";
 393     fNbin=fhistoInvMass->GetNbinsX();
 394     cout<<"Kept original number of bins: "<<fNbin<<endl;
 395   } else{
 396     fhistoInvMass->Rebin(bingroup);
 397     fNbin = fhistoInvMass->GetNbinsX();
 398     cout<<"New number of bins: "<<fNbin<<endl;
 399   }
 400
 401
 402 }
 403
 404 //************* fit
 405
 406 //___________________________________________________________________________
 407
 408 Double_t AliHFMassFitter::FitFunction4MassDistr (Double_t *x, Double_t *par){
 409   // Fit function for signal+background
 410
 411
 412   //exponential or linear fit
 413   //
 414   // par[0] = tot integral
 415   // par[1] = slope
 416   // par[2] = gaussian integral
 417   // par[3] = gaussian mean
 418   // par[4] = gaussian sigma
 419
 420   Double_t total,bkg=0,sgn=0;
 421
 422   if (ftypeOfFit4Bkg==0 || ftypeOfFit4Bkg==1) {
 423     if(ftypeOfFit4Sgn == 0) {
 424
 425       Double_t parbkg[2] = {par[0]-par[2], par[1]};
 426       bkg = FitFunction4Bkg(x,parbkg);
 427     }
 428     if(ftypeOfFit4Sgn == 1) {
 429       Double_t parbkg[5] = {par[2],par[3],ffactor*par[4],par[0]-2*par[2], par[1]};
 430       bkg = FitFunction4Bkg(x,parbkg);
 431     }
 432
 433     sgn = FitFunction4Sgn(x,&par[2]);
 434
 435   }
 436
 437   //polynomial fit
 438
 439     // par[0] = tot integral
 440     // par[1] = coef1
 441     // par[2] = coef2
 442     // par[3] = gaussian integral
 443     // par[4] = gaussian mean
 444     // par[5] = gaussian sigma
 445
 446   if (ftypeOfFit4Bkg==2) {
 447
 448     if(ftypeOfFit4Sgn == 0) {
 449
 450       Double_t parbkg[3] = {par[0]-par[3], par[1], par[2]};
 451       bkg = FitFunction4Bkg(x,parbkg);
 452     }
 453     if(ftypeOfFit4Sgn == 1) {
 454
 455       Double_t parbkg[6] = {par[3],par[4],ffactor*par[5],par[0]-2*par[3], par[1], par[2]};
 456      bkg = FitFunction4Bkg(x,parbkg);
 457     }
 458
 459     sgn = FitFunction4Sgn(x,&par[3]);
 460   }
 461
 462   if (ftypeOfFit4Bkg==3) {
 463
 464     if(ftypeOfFit4Sgn == 0) {
 465         bkg=FitFunction4Bkg(x,par);
 466         sgn=FitFunction4Sgn(x,&par[1]);
 467     }
 468     if(ftypeOfFit4Sgn == 1) {
 469       Double_t parbkg[4]={par[1],par[2],ffactor*par[3],par[0]};
 470       bkg=FitFunction4Bkg(x,parbkg);
 471       sgn=FitFunction4Sgn(x,&par[1]);
 472     }
 473   }
 474
 475   total = bkg + sgn;
 476
 477   return  total;
 478 }
 479
 480 //_________________________________________________________________________
 481 Double_t AliHFMassFitter::FitFunction4Sgn (Double_t *x, Double_t *par){
 482   // Fit function for the signal
 483
 484   //gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
 485   //Par:
 486   // * [0] = integralSgn
 487   // * [1] = mean
 488   // * [2] = sigma
 489   //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
 490
 491   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]);
 492
 493 }
 494
 495 //__________________________________________________________________________
 496
 497 Double_t AliHFMassFitter::FitFunction4Bkg (Double_t *x, Double_t *par){
 498   // Fit function for the background
 499
 500   Double_t maxDeltaM = 4.*fSigmaSgn;
 501   if(fSideBands && TMath::Abs(x[0]-fMass) < maxDeltaM) {
 502     TF1::RejectPoint();
 503     return 0;
 504   }
 505   Int_t firstPar=0;
 506   Double_t gaus2=0,total=-1;
 507   if(ftypeOfFit4Sgn == 1){
 508     firstPar=3;
 509     //gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
 510     //Par:
 511     // * [0] = integralSgn
 512     // * [1] = mean
 513     // * [2] = sigma
 514     //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
 515     gaus2 = FitFunction4Sgn(x,par);
 516   }
 517
 518   switch (ftypeOfFit4Bkg){
 519   case 0:
 520     //exponential
 521     //exponential = A*exp(B*x) -> integral(exponential)=A/B*exp(B*x)](min,max)
 522     //-> A = B*integral/(exp(B*max)-exp(B*min)) where integral can be written
 523     //as integralTot- integralGaus (=par [2])
 524     //Par:
 525     // * [0] = integralBkg;
 526     // * [1] = B;
 527     //exponential = [1]*[0]/(exp([1]*max)-exp([1]*min))*exp([1]*x)
 528     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]);
 529     break;
 530   case 1:
 531     //linear
 532     //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)
 533     // * [0] = integralBkg;
 534     // * [1] = b;
 535     total= par[0+firstPar]/(fmaxMass-fminMass)+par[1+firstPar]*(x[0]-0.5*(fmaxMass+fminMass));
 536     break;
 537   case 2:
 538     //polynomial
 539     //y=a+b*x+c*x**2 -> integral = a(max-min) + 1/2*b*(max^2-min^2) +
 540     //+ 1/3*c*(max^3-min^3) ->
 541     //a = (integral-1/2*b*(max^2-min^2)-1/3*c*(max^3-min^3))/(max-min)
 542     // * [0] = integralBkg;
 543     // * [1] = b;
 544     // * [2] = c;
 545     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));
 546     break;
 547   case 3:
 548     total=par[0+firstPar];
 549     break;
 550 //   default:
 551 //     Types of Fit Functions for Background:
 552 //     * 0 = exponential;
 553 //     * 1 = linear;
 554 //     * 2 = polynomial 2nd order
 555 //     * 3 = no background"<<endl;
 556
 557   }
 558   return total+gaus2;
 559 }
 560
 561 //__________________________________________________________________________
 562 Bool_t AliHFMassFitter::SideBandsBounds(){
 563
 564   Double_t width=fhistoInvMass->GetBinWidth(8);
 565   if (fNbin==0) fNbin=fhistoInvMass->GetNbinsX();
 566   Double_t minHisto=fhistoInvMass->GetBinLowEdge(1);
 567   Double_t maxHisto=fhistoInvMass->GetBinLowEdge(fNbin)+width;
 568
 569   if(fMass-fminMass < 0 || fmaxMass-fMass <0) {
 570     cout<<"Left limit of range > mean or right limit of range < mean: change left/right limit or initial mean value"<<endl;
 571     return kFALSE;
 572   }
 573
 574   if((fminMass!=minHisto || fmaxMass!= maxHisto) && (fMass-4.*fSigmaSgn-fminMass) <= 0){
 575     Double_t coeff = (fMass-fminMass)/fSigmaSgn;
 576
 577     fSideBandl=(Int_t)((fMass-0.5*coeff*fSigmaSgn-fminMass)/width);
 578     fSideBandr=(Int_t)((fMass+0.5*coeff*fSigmaSgn-fminMass)/width);
 579     cout<<"Changed number of sigma from 4 to "<<0.5*coeff<<" for the estimation of the side bands"<<endl;
 580     if (coeff<3) cout<<"Side bands inside 3 sigma, may be better use ftypeOfFit4Bkg = 3 (only signal)"<<endl;
 581
 582     if (coeff<2) {
 583       cout<<"Side bands inside 2 sigma. Change mode: ftypeOfFit4Bkg = 3"<<endl;
 584       ftypeOfFit4Bkg=3;
 585       //set binleft and right without considering SetRangeFit- anyway no bkg!
 586       fSideBandl=(Int_t)((fMass-4.*fSigmaSgn-minHisto)/width);
 587       fSideBandr=(Int_t)((fMass+4.*fSigmaSgn-minHisto)/width);
 588     }
 589   }
 590   else {
 591   fSideBandl=(Int_t)((fMass-4.*fSigmaSgn-minHisto)/width);
 592   fSideBandr=(Int_t)((fMass+4.*fSigmaSgn-minHisto)/width);
 593 //   cout<<"\tfMass = "<<fMass<<"\tfSigmaSgn = "<<fSigmaSgn<<"\tminHisto = "<<minHisto<<endl;
 594 //   cout<<"\tbinleft = "<<fSideBandl<<"\tbinright = "<<fSideBandr<<endl;
 595   }
 596
 597   if (fSideBandl==0) {
 598     cout<<"Error! Range too little";
 599     return kFALSE;
 600   }
 601   return kTRUE;
 602 }
 603
 604 //__________________________________________________________________________
 605
 606 void AliHFMassFitter::GetSideBandsBounds(Int_t &left, Int_t &right) const{
 607   if (fSideBandl==0 && fSideBandr==0){
 608     cout<<"Use MassFitter method first"<<endl;
 609     return;
 610   }
 611   left=fSideBandl;
 612   right=fSideBandr;
 613 }
 614
 615 //__________________________________________________________________________
 616
 617 Bool_t AliHFMassFitter::MassFitter(Bool_t draw){
 618   // Main method of the class: performs the fit of the histogram
 619
 620   Int_t nFitPars=0; //total function's number of parameters
 621   switch (ftypeOfFit4Bkg){
 622   case 0:
 623     nFitPars=5; //3+2
 624     break;
 625   case 1:
 626     nFitPars=5; //3+2
 627     break;
 628   case 2:
 629     nFitPars=6; //3+3
 630     break;
 631   case 3:
 632     nFitPars=4; //3+1
 633     break;
 634   }
 635
 636   Int_t bkgPar = nFitPars-3; //background function's number of parameters
 637
 638   cout<<"nFitPars = "<<nFitPars<<"\tbkgPar = "<<bkgPar<<endl;
 639
 640   //function names
 641   TString bkgname="funcbkg";
 642   TString bkg1name="funcbkg1";
 643   TString massname="funcmass";
 644
 645   //Total integral
 646   Double_t totInt = fhistoInvMass->Integral("width");
 647
 648   fSideBands = kTRUE;
 649   Double_t width=fhistoInvMass->GetBinWidth(8);
 650
 651   if (fNbin==0) fNbin=fhistoInvMass->GetNbinsX();
 652   Double_t minHisto=fhistoInvMass->GetBinLowEdge(1);
 653   Double_t maxHisto=fhistoInvMass->GetBinLowEdge(fNbin)+width;
 654   Bool_t ok=SideBandsBounds();
 655   if(!ok) return kFALSE;
 656
 657   //sidebands integral - first approx (from histo)
 658   Double_t sideBandsInt=(Double_t)fhistoInvMass->Integral(1,fSideBandl,"width") + (Double_t)fhistoInvMass->Integral(fSideBandr,fNbin,"width");
 659   cout<<"------nbin = "<<fNbin<<"\twidth = "<<width<<"\tbinleft = "<<fSideBandl<<"\tbinright = "<<fSideBandr<<endl;
 660   cout<<"------sideBandsInt - first approx = "<<sideBandsInt<<endl;
 661   if (sideBandsInt<=0) {
 662     cout<<"! sideBandsInt <=0. There's a problem, cannot start the fit"<<endl;
 663     return kFALSE;
 664   }
 665
 666   /*Fit Bkg*/
 667
 668
 669   TF1 *funcbkg = new TF1(bkgname.Data(),this,&AliHFMassFitter::FitFunction4Bkg,minHisto,maxHisto,bkgPar,"AliHFMassFitter","FitFunction4Bkg");
 670   cout<<"Function name = "<<funcbkg->GetName()<<endl<<endl;
 671
 672   funcbkg->SetLineColor(2); //red
 673
 674   //first fit for bkg: approx bkgint
 675
 676   switch (ftypeOfFit4Bkg) {
 677   case 0: //gaus+expo
 678     funcbkg->SetParNames("BkgInt","Slope");
 679     funcbkg->SetParameters(sideBandsInt,-2.);
 680     break;
 681   case 1:
 682     funcbkg->SetParNames("BkgInt","Slope");
 683     funcbkg->SetParameters(sideBandsInt,-100.);
 684     break;
 685   case 2:
 686     funcbkg->SetParNames("BkgInt","Coef1","Coef2");
 687     funcbkg->SetParameters(sideBandsInt,-10.,5);
 688     break;
 689   case 3:
 690     if(ftypeOfFit4Sgn==0){
 691       funcbkg->SetParNames("Const");
 692       funcbkg->SetParameter(0,0.);
 693       funcbkg->FixParameter(0,0.);
 694     }
 695     break;
 696   default:
 697     cout<<"Wrong choise of ftypeOfFit4Bkg ("<<ftypeOfFit4Bkg<<")"<<endl;
 698     return kFALSE;
 699     break;
 700   }
 701   cout<<"\nBACKGROUND FIT - only combinatorial"<<endl;
 702   Int_t ftypeOfFit4SgnBkp=ftypeOfFit4Sgn;
 703
 704   Double_t intbkg1=0,slope1=0,conc1=0;
 705   //if only signal and reflection: skip
 706   if (!(ftypeOfFit4Bkg==3 && ftypeOfFit4Sgn==1)) {
 707     ftypeOfFit4Sgn=0;
 708     fhistoInvMass->Fit(bkgname.Data(),"R,L,E,0");
 709
 710     for(Int_t i=0;i<bkgPar;i++){
 711       fFitPars[i]=funcbkg->GetParameter(i);
 712       //cout<<i<<"\t"<<funcbkg->GetParameter(i)<<"\t";
 713       fFitPars[nFitPars+2*bkgPar+3+i]= funcbkg->GetParError(i);
 714       //cout<<nFitPars+2*bkgPar+3+i<<"\t"<< funcbkg->GetParError(i)<<endl;
 715     }
 716     fSideBands = kFALSE;
 717     //intbkg1 = funcbkg->GetParameter(0);
 718     funcbkg->SetRange(fminMass,fmaxMass);
 719     intbkg1 = funcbkg->Integral(fminMass,fmaxMass);
 720     if(ftypeOfFit4Bkg!=3) slope1 = funcbkg->GetParameter(1);
 721     if(ftypeOfFit4Bkg==2) conc1 = funcbkg->GetParameter(2);
 722     cout<<"Primo fit: \nintbkg1 = "<<intbkg1<<"\t(Compare with par0 = "<<funcbkg->GetParameter(0)<<")\nslope1= "<<slope1<<"\nconc1 = "<<conc1<<endl;
 723   }
 724   else cout<<"\t\t//"<<endl;
 725
 726   ftypeOfFit4Sgn=ftypeOfFit4SgnBkp;
 727   TF1 *funcbkg1=0;
 728   if (ftypeOfFit4Sgn == 1) {
 729     cout<<"\nBACKGROUND FIT WITH REFLECTION"<<endl;
 730     bkgPar+=3;
 731
 732     cout<<"nFitPars = "<<nFitPars<<"\tbkgPar = "<<bkgPar<<endl;
 733
 734     funcbkg1 = new TF1(bkg1name.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,bkgPar,"AliHFMassFitter","FitFunction4Bkg");
 735     cout<<"Function name = "<<funcbkg1->GetName()<<endl;
 736
 737     funcbkg1->SetLineColor(2); //red
 738
 739     if(ftypeOfFit4Bkg==2){
 740       cout<<"*** Polynomial Fit ***"<<endl;
 741       funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Coef1","Coef2");
 742       funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1,conc1);
 743 //     cout<<"Parameters set to: "<<0.5*(totInt-intbkg1)<<"\t"<<fMass<<"\t"<<ffactor*fSigmaSgn<<"\t"<<intbkg1<<"\t"<<slope1<<"\t"<<conc1<<"\t"<<endl;
 744 //     cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<bkgPar<<"\tgsidebands = "<<fSideBands<<endl;
 745     } else{
 746       if(ftypeOfFit4Bkg==3) //no background: gaus sign+ gaus broadened
 747         {
 748           cout<<"*** No background Fit ***"<<endl;
 749           funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","Const");
 750           funcbkg1->SetParameters(0.5*totInt,fMass,ffactor*fSigmaSgn,0.);
 751           funcbkg1->FixParameter(3,0.);
 752         } else{ //expo or linear
 753           if(ftypeOfFit4Bkg==0) cout<<"*** Exponential Fit ***"<<endl;
 754           if(ftypeOfFit4Bkg==1) cout<<"*** Linear Fit ***"<<endl;
 755           funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Slope");
 756           funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1);
 757         }
 758     }
 759     fhistoInvMass->Fit(bkg1name.Data(),"R,L,E,+,0");
 760
 761     for(Int_t i=0;i<bkgPar;i++){
 762       fFitPars[bkgPar-3+i]=funcbkg1->GetParameter(i);
 763       //cout<<bkgPar-3+i<<"\t"<<funcbkg1->GetParameter(i);
 764       fFitPars[nFitPars+3*bkgPar-6+i]= funcbkg1->GetParError(i);
 765       //cout<<"\t"<<nFitPars+3*bkgPar-6+i<<"\t"<<funcbkg1->GetParError(i)<<endl;
 766     }
 767
 768     intbkg1=funcbkg1->GetParameter(3);
 769     if(ftypeOfFit4Bkg!=3) slope1 = funcbkg1->GetParameter(4);
 770     if(ftypeOfFit4Bkg==2) conc1 = funcbkg1->GetParameter(5);
 771
 772   } else {
 773     bkgPar+=3;
 774
 775     for(Int_t i=0;i<3;i++){
 776       fFitPars[bkgPar-3+i]=0.;
 777       cout<<bkgPar-3+i<<"\t"<<0.<<"\t";
 778       fFitPars[nFitPars+3*bkgPar-6+i]= 0.;
 779       cout<<nFitPars+3*bkgPar-6+i<<"\t"<<0.<<endl;
 780     }
 781
 782     for(Int_t i=0;i<bkgPar-3;i++){
 783       fFitPars[bkgPar+i]=funcbkg->GetParameter(i);
 784       cout<<bkgPar+i<<"\t"<<funcbkg->GetParameter(i)<<"\t";
 785       fFitPars[nFitPars+3*bkgPar-3+i]= funcbkg->GetParError(i);
 786       cout<<nFitPars+3*bkgPar-3+i<<"\t"<< funcbkg->GetParError(i)<<endl;
 787     }
 788
 789
 790   }
 791
 792   //sidebands integral - second approx (from fit)
 793   fSideBands = kFALSE;
 794   Double_t bkgInt;
 795   cout<<"Compare intbkg1 = "<<intbkg1<<" and integral = ";
 796   if(ftypeOfFit4Sgn == 1) bkgInt=funcbkg1->Integral(fminMass,fmaxMass);
 797   else bkgInt=funcbkg->Integral(fminMass,fmaxMass);
 798   cout<</*"------BkgInt(Fit) = "<<*/bkgInt<<endl;
 799
 800   //Signal integral - first approx
 801   Double_t sgnInt;
 802   sgnInt = totInt-bkgInt;
 803   cout<<"------TotInt = "<<totInt<<"\tsgnInt = "<<sgnInt<<endl;
 804
 805   /*Fit All Mass distribution with exponential + gaussian (+gaussiam braodened) */
 806   TF1 *funcmass = new TF1(massname.Data(),this,&AliHFMassFitter::FitFunction4MassDistr,fminMass,fmaxMass,nFitPars,"AliHFMassFitter","FitFunction4MassDistr");
 807   cout<<"Function name = "<<funcmass->GetName()<<endl<<endl;
 808   funcmass->SetLineColor(4); //blue
 809
 810   //Set parameters
 811   cout<<"\nTOTAL FIT"<<endl;
 812
 813   if(nFitPars==5){
 814     funcmass->SetParNames("TotInt","Slope","SgnInt","Mean","Sigma");
 815     funcmass->SetParameters(totInt,slope1,sgnInt,fMass,fSigmaSgn);
 816     //cout<<"Parameters set to: "<<totInt<<"\t"<<slope1<<"\t"<<sgnInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl;
 817 //     cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<nFitPars<<"\tgsidebands = "<<fSideBands<<endl;
 818     funcmass->FixParameter(0,totInt);
 819   }
 820   if (nFitPars==6){
 821     funcmass->SetParNames("TotInt","Coef1","Coef2","SgnInt","Mean","Sigma");
 822     funcmass->SetParameters(totInt,slope1,conc1,sgnInt,fMass,fSigmaSgn);
 823 //     cout<<"Parameters set to: "<<totInt<<"\t"<<slope1<<"\t"<<conc1<<"\t"<<sgnInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl;
 824 //     cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<nFitPars<<"\tgsidebands = "<<fSideBands<<endl;
 825     funcmass->FixParameter(0,totInt);
 826   }
 827   if(nFitPars==4){
 828     funcmass->SetParNames("Const","SgnInt","Mean","Sigma");
 829     if(ftypeOfFit4Sgn == 1) funcmass->SetParameters(0.,0.5*totInt,fMass,fSigmaSgn);
 830     else funcmass->SetParameters(0.,totInt,fMass,fSigmaSgn);
 831     funcmass->FixParameter(0,0.);
 832     //cout<<"Parameters set to: "<<0.5*totInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl;
 833     //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<nFitPars<<"\tgsidebands = "<<fSideBands<<endl;
 834
 835   }
 836
 837   fhistoInvMass->Fit(massname.Data(),"R,L,E,+,0");
 838   cout<<"fit done"<<endl;
 839   //reset value of fMass and fSigmaSgn to those found from fit
 840   fMass=funcmass->GetParameter(nFitPars-2);
 841   fSigmaSgn=funcmass->GetParameter(nFitPars-1);
 842
 843   for(Int_t i=0;i<nFitPars;i++){
 844     fFitPars[i+2*bkgPar-3]=funcmass->GetParameter(i);
 845     fFitPars[nFitPars+4*bkgPar-6+i]= funcmass->GetParError(i);
 846     //cout<<i+2*bkgPar-3<<"\t"<<funcmass->GetParameter(i)<<"\t\t"<<nFitPars+4*bkgPar-6+i<<"\t"<<funcmass->GetParError(i)<<endl;
 847   }
 848   /*
 849   //check: cout parameters
 850   for(Int_t i=0;i<2*(nFitPars+2*bkgPar-3);i++){
 851     cout<<i<<"\t"<<fFitPars[i]<<endl;
 852     }
 853   */
 854
 855 //   if(draw){
 856 //     TCanvas *canvas=new TCanvas(fhistoInvMass->GetName(),fhistoInvMass->GetName());
 857 //     TH1F *fhistocopy=new TH1F(*fhistoInvMass);
 858 //     canvas->cd();
 859 //     fhistocopy->DrawClone();
 860 //     if(ftypeOfFit4Sgn == 1) funcbkg1->DrawClone("sames");
 861 //     else funcbkg->DrawClone("sames");
 862 //     funcmass->DrawClone("sames");
 863 //     cout<<"Drawn"<<endl;
 864 //   }
 865
 866   if(funcmass->GetParameter(nFitPars-1) <0 || funcmass->GetParameter(nFitPars-2) <0 || funcmass->GetParameter(nFitPars-3) <0 ) {
 867     cout<<"IntS or mean or sigma negative. You may tray to SetInitialGaussianSigma(..) and SetInitialGaussianMean(..)"<<endl;
 868     return kFALSE;
 869   }
 870
 871   //increase counter of number of fits done
 872   fcounter++;
 873
 874   if (ftypeOfFit4Sgn == 1 && funcbkg1) {
 875     delete funcbkg1;
 876     funcbkg1=NULL;
 877   }
 878   if (funcbkg) {
 879     delete funcbkg;
 880     funcbkg=NULL;
 881   }
 882   if (funcmass) {
 883     delete funcmass;
 884     funcmass=NULL;
 885   }
 886
 887   AddFunctionsToHisto();
 888   if (draw) DrawFit();
 889
 890
 891   return kTRUE;
 892 }
 893 //_________________________________________________________________________
 894 Double_t AliHFMassFitter::GetChiSquare() const{
 895   TF1 *funcmass=(TF1*)fhistoInvMass->GetFunction("funcmass");
 896   return funcmass->GetChisquare();
 897 }
 898
 899 //_________________________________________________________________________
 900 Double_t AliHFMassFitter::GetReducedChiSquare() const{
 901   TF1 *funcmass=(TF1*)fhistoInvMass->GetFunction("funcmass");
 902   return funcmass->GetChisquare()/funcmass->GetNDF();
 903 }
 904
 905 //*********output
 906
 907 //_________________________________________________________________________
 908 void  AliHFMassFitter::GetFitPars(Float_t *vector) const {
 909   // Return fit parameters
 910
 911   for(Int_t i=0;i<fParsSize;i++){
 912     vector[i]=fFitPars[i];
 913   }
 914 }
 915
 916
 917 //_________________________________________________________________________
 918 void AliHFMassFitter::AddFunctionsToHisto(){
 919
 920   cout<<"AddFunctionsToHisto called"<<endl;
 921   TString bkgname = "funcbkg";
 922   Int_t np=-99;
 923   switch (ftypeOfFit4Bkg){
 924   case 0: //expo
 925     np=2;
 926     break;
 927   case 1: //linear
 928     np=2;
 929     break;
 930   case 2: //pol2
 931     np=3;
 932     break;
 933   case 3: //no bkg
 934     np=1;
 935     break;
 936   }
 937   if (ftypeOfFit4Sgn == 1) {
 938     bkgname += 1;
 939     np+=3;
 940   }
 941
 942   TString bkgnamesave=bkgname;
 943   TString testname=bkgname;
 944   testname += "FullRange";
 945   TF1 *testfunc=(TF1*)fhistoInvMass->FindObject(testname.Data());
 946   if(testfunc){
 947     cout<<"AddFunctionsToHisto already used: exiting...."<<endl;
 948     return;
 949   }
 950
 951   TList *hlist=fhistoInvMass->GetListOfFunctions();
 952   hlist->ls();
 953
 954   TF1 *b=(TF1*)hlist->FindObject(bkgname.Data());
 955   if(!b){
 956     cout<<bkgname<<" not found, cannot produce "<<bkgname<<"FullRange and "<<bkgname<<"Recalc"<<endl;
 957     return;
 958   }
 959
 960   bkgname += "FullRange";
 961   TF1 *bfullrange=new TF1(bkgname.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,np,"AliHFMassFitter","FitFunction4Bkg");
 962   //cout<<bfullrange->GetName()<<endl;
 963   for(Int_t i=0;i<np;i++){
 964     //cout<<i<<" di "<<np<<endl;
 965     bfullrange->SetParName(i,b->GetParName(i));
 966     bfullrange->SetParameter(i,b->GetParameter(i));
 967     bfullrange->SetParError(i,b->GetParError(i));
 968   }
 969   bfullrange->SetLineStyle(4);
 970   bfullrange->SetLineColor(14);
 971
 972   bkgnamesave += "Recalc";
 973
 974   TF1 *blastpar=new TF1(bkgnamesave.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,np,"AliHFMassFitter","FitFunction4Bkg");
 975
 976   TF1 *mass=fhistoInvMass->GetFunction("funcmass");
 977
 978   if (!mass){
 979     cout<<"funcmass doesn't exist "<<endl;
 980     return;
 981   }
 982
 983   blastpar->SetParameter(0,mass->GetParameter(0)-mass->GetParameter(np));
 984   blastpar->SetParError(0,mass->GetParError(np));
 985   if (np>=2) {
 986     blastpar->SetParameter(1,mass->GetParameter(1));
 987     blastpar->SetParError(1,mass->GetParError(1));
 988   }
 989   if (np==3) {
 990     blastpar->SetParameter(2,mass->GetParameter(2));
 991     blastpar->SetParError(2,mass->GetParError(2));
 992   }
 993
 994   blastpar->SetLineStyle(1);
 995   blastpar->SetLineColor(2);
 996
 997   hlist->Add(bfullrange->Clone());
 998   hlist->Add(blastpar->Clone());
 999   hlist->ls();
1000
1001   if(bfullrange) {
1002     delete bfullrange;
1003     bfullrange=NULL;
1004   }
1005   if(blastpar) {
1006     delete blastpar;
1007     blastpar=NULL;
1008   }
1009 }
1010
1011 //_________________________________________________________________________
1012
1013 TH1F* AliHFMassFitter::GetHistoClone() const{
1014
1015   TH1F* hout=(TH1F*)fhistoInvMass->Clone(fhistoInvMass->GetName());
1016   return hout;
1017 }
1018 //_________________________________________________________________________
1019
1020 void AliHFMassFitter::WriteHisto(TString path) {
1021   if (fcounter == 0) {
1022     cout<<"Use MassFitter method before WriteHisto"<<endl;
1023     return;
1024   }
1025   TH1F* hget=(TH1F*)fhistoInvMass->Clone();
1026
1027   path += "HFMassFitterOutput.root";
1028   TFile *output;
1029
1030   if (fcounter == 1) output = new TFile(path.Data(),"recreate");
1031   else output = new TFile(path.Data(),"update");
1032   output->cd();
1033   hget->Write();
1034   output->Close();
1035
1036   cout<<fcounter<<" "<<hget->GetName()<<" written in "<<path<<endl;
1037
1038   if(output) {
1039     delete output;
1040     output=NULL;
1041   }
1042
1043 }
1044
1045 //_________________________________________________________________________
1046
1047 void AliHFMassFitter::WriteNtuple(TString path) const{
1048   //TNtuple* nget=(TNtuple*)fntuParam->Clone();
1049   path += "HFMassFitterOutput.root";
1050   TFile *output = new TFile(path.Data(),"update");
1051   output->cd();
1052   fntuParam->Write();
1053   //nget->Write();
1054   output->Close();
1055   //cout<<nget->GetName()<<" written in "<<path<<endl;
1056   cout<<fntuParam->GetName()<<" written in "<<path<<endl;
1057   /*
1058   if(nget) {
1059     //delete nget;
1060     nget=NULL;
1061   }
1062   */
1063   if(output) {
1064     delete output;
1065     output=NULL;
1066   }
1067 }
1068
1069 //_________________________________________________________________________
1070
1071 void AliHFMassFitter::DrawFit() const{
1072
1073   TString cvtitle="fit of ";
1074   cvtitle+=fhistoInvMass->GetName();
1075   TString cvname="c";
1076   cvname+=fcounter;
1077
1078   TCanvas *c=new TCanvas(cvname,cvtitle);
1079   c->cd();
1080   GetHistoClone()->Draw("hist");
1081   GetHistoClone()->GetFunction("funcbkgFullRange")->Draw("same");
1082   GetHistoClone()->GetFunction("funcbkgRecalc")->Draw("same");
1083   GetHistoClone()->GetFunction("funcmass")->Draw("same");
1084
1085
1086 }
1087
1088
1089 //************ significance
1090
1091 //_________________________________________________________________________
1092
1093 void AliHFMassFitter::Signal(Double_t nOfSigma,Double_t &signal,Double_t &errsignal) const {
1094   // Return signal integral in mean+- n sigma
1095
1096   if(fcounter==0) {
1097     cout<<"Use MassFitter method before Signal"<<endl;
1098     return;
1099   }
1100
1101   //functions names
1102   TString bkgname= "funcbkgRecalc";
1103   TString bkg1name="funcbkg1Recalc";
1104   TString massname="funcmass";
1105
1106
1107   TF1 *funcbkg=0;
1108   TF1 *funcmass=fhistoInvMass->GetFunction(massname.Data());
1109   if(!funcmass){
1110     cout<<"AliHFMassFitter::Signal() ERROR -> Mass distr function not found!"<<endl;
1111     return;
1112   }
1113
1114   if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
1115   else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
1116
1117   if(!funcbkg){
1118     cout<<"AliHFMassFitter::Signal() ERROR -> Bkg function not found!"<<endl;
1119     return;
1120   }
1121
1122   Int_t np=-99;
1123   switch (ftypeOfFit4Bkg){
1124   case 0: //expo
1125     np=2;
1126     break;
1127   case 1: //linear
1128     np=2;
1129     break;
1130   case 2: //pol2
1131     np=3;
1132     break;
1133   case 3: //no bkg
1134     np=1;
1135     break;
1136   }
1137
1138   Double_t intS,intSerr;
1139
1140  //relative error evaluation
1141   intS=funcmass->GetParameter(np);
1142   intSerr=funcmass->GetParError(np);
1143   cout<<"Sgn relative error evaluation from fit: "<<intSerr/intS<<endl;
1144   Double_t background,errbackground;
1145   Background(nOfSigma,background,errbackground);
1146
1147   //signal +/- error in nsigma
1148   Double_t min=fMass-nOfSigma*fSigmaSgn;
1149   Double_t max=fMass+nOfSigma*fSigmaSgn;
1150
1151   Double_t mass=funcmass->Integral(min, max)/fhistoInvMass->GetBinWidth(4);
1152   signal=mass - background;
1153   errsignal=TMath::Sqrt((intSerr/intS*mass)*(intSerr/intS*mass)/*assume relative error is the same as for total integral*/ + errbackground*errbackground);
1154   return;
1155 }
1156
1157 //_________________________________________________________________________
1158
1159 void AliHFMassFitter::Signal(Double_t min, Double_t max, Double_t &signal,Double_t &errsignal) const {
1160
1161   // Return signal integral in a range
1162
1163   if(fcounter==0) {
1164     cout<<"Use MassFitter method before Signal"<<endl;
1165     return;
1166   }
1167
1168   //functions names
1169   TString bkgname="funcbkgRecalc";
1170   TString bkg1name="funcbkg1Recalc";
1171   TString massname="funcmass";
1172
1173
1174   TF1 *funcbkg=0;
1175   TF1 *funcmass=fhistoInvMass->GetFunction(massname.Data());
1176   if(!funcmass){
1177     cout<<"AliHFMassFitter::Signal() ERROR -> Mass distr function not found!"<<endl;
1178     return;
1179   }
1180
1181   if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
1182   else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
1183
1184   if(!funcbkg){
1185     cout<<"AliHFMassFitter::Signal() ERROR -> Bkg function not found!"<<endl;
1186     return;
1187   }
1188
1189   Int_t np=-99;
1190   switch (ftypeOfFit4Bkg){
1191   case 0: //expo
1192     np=2;
1193     break;
1194   case 1: //linear
1195     np=2;
1196     break;
1197   case 2: //pol2
1198     np=3;
1199     break;
1200   case 3: //no bkg
1201     np=1;
1202     break;
1203   }
1204
1205   Double_t intS,intSerr;
1206
1207  //relative error evaluation
1208   intS=funcmass->GetParameter(np);
1209   intSerr=funcmass->GetParError(np);
1210
1211   cout<<"Sgn relative error evaluation from fit: "<<intSerr/intS<<endl;
1212   Double_t background,errbackground;
1213   Background(min,max,background,errbackground);
1214
1215   //signal +/- error in the range
1216
1217   Double_t mass=funcmass->Integral(min, max)/fhistoInvMass->GetBinWidth(4);
1218   signal=mass - background;
1219   errsignal=(intSerr/intS)*signal;/*assume relative error is the same as for total integral*/
1220
1221 }
1222
1223 //_________________________________________________________________________
1224
1225 void AliHFMassFitter::Background(Double_t nOfSigma,Double_t &background,Double_t &errbackground) const {
1226   // Return background integral in mean+- n sigma
1227
1228   if(fcounter==0) {
1229     cout<<"Use MassFitter method before Background"<<endl;
1230     return;
1231   }
1232
1233   //functions names
1234   TString bkgname="funcbkgRecalc";
1235   TString bkg1name="funcbkg1Recalc";
1236
1237   TF1 *funcbkg=0;
1238   if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
1239   else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
1240   if(!funcbkg){
1241     cout<<"AliHFMassFitter::Background() ERROR -> Bkg function not found!"<<endl;
1242     return;
1243   }
1244
1245   Double_t intB,intBerr;
1246
1247   //relative error evaluation: from final parameters of the fit
1248   if(ftypeOfFit4Bkg==3 && ftypeOfFit4Sgn == 0) cout<<"No background fit: Bkg relative error evaluation put to zero"<<endl;
1249   else{
1250     intB=funcbkg->GetParameter(0);
1251     intBerr=funcbkg->GetParError(0);
1252     cout<<"Bkg relative error evaluation: from final parameters of the fit: "<<intBerr/intB<<endl;
1253   }
1254
1255   Double_t min=fMass-nOfSigma*fSigmaSgn;
1256   Double_t max=fMass+nOfSigma*fSigmaSgn;
1257
1258   //relative error evaluation: from histo
1259
1260   intB=fhistoInvMass->Integral(1,fSideBandl)+fhistoInvMass->Integral(fSideBandr,fNbin);
1261   Double_t sum2=0;
1262   for(Int_t i=1;i<=fSideBandl;i++){
1263     sum2+=fhistoInvMass->GetBinError(i)*fhistoInvMass->GetBinError(i);
1264   }
1265   for(Int_t i=fSideBandr;i<=fNbin;i++){
1266     sum2+=fhistoInvMass->GetBinError(i)*fhistoInvMass->GetBinError(i);
1267   }
1268
1269   intBerr=TMath::Sqrt(sum2);
1270   cout<<"Bkg relative error evaluation: from histo: "<<intBerr/intB<<endl;
1271
1272   cout<<"Last estimation of bkg error is used"<<endl;
1273
1274   //backround +/- error in nsigma
1275   if (ftypeOfFit4Bkg == 3 && ftypeOfFit4Sgn == 0) {
1276     background = 0;
1277     errbackground = 0;
1278   }
1279   else{
1280     background=funcbkg->Integral(min,max)/(Double_t)fhistoInvMass->GetBinWidth(2);
1281     errbackground=intBerr/intB*background; // assume relative error is the same as for total integral
1282     //cout<<"integral = "<<funcbkg->Integral(min, max)<<"\tbinW = "<<fhistoInvMass->GetBinWidth(2)<<endl;
1283   }
1284   return;
1285
1286 }
1287 //___________________________________________________________________________
1288
1289 void AliHFMassFitter::Background(Double_t min, Double_t max, Double_t &background,Double_t &errbackground) const {
1290   // Return background integral in a range
1291
1292   if(fcounter==0) {
1293     cout<<"Use MassFitter method before Background"<<endl;
1294     return;
1295   }
1296
1297   //functions names
1298   TString bkgname="funcbkgRecalc";
1299   TString bkg1name="funcbkg1Recalc";
1300
1301   TF1 *funcbkg=0;
1302   if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
1303   else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
1304   if(!funcbkg){
1305     cout<<"AliHFMassFitter::Background() ERROR -> Bkg function not found!"<<endl;
1306     return;
1307   }
1308
1309
1310   Double_t intB,intBerr;
1311
1312   //relative error evaluation: from final parameters of the fit
1313   if(ftypeOfFit4Bkg==3 && ftypeOfFit4Sgn == 0) cout<<"No background fit: Bkg relative error evaluation put to zero"<<endl;
1314   else{
1315     intB=funcbkg->GetParameter(0);
1316     intBerr=funcbkg->GetParError(0);
1317     cout<<"Bkg relative error evaluation: from final parameters of the fit: "<<intBerr/intB<<endl;
1318   }
1319
1320   //relative error evaluation: from histo
1321
1322   intB=fhistoInvMass->Integral(1,fSideBandl)+fhistoInvMass->Integral(fSideBandr,fNbin);
1323   Double_t sum2=0;
1324   for(Int_t i=1;i<=fSideBandl;i++){
1325     sum2+=fhistoInvMass->GetBinError(i)*fhistoInvMass->GetBinError(i);
1326   }
1327   for(Int_t i=fSideBandr;i<=fNbin;i++){
1328     sum2+=fhistoInvMass->GetBinError(i)*fhistoInvMass->GetBinError(i);
1329   }
1330
1331   intBerr=TMath::Sqrt(sum2);
1332   cout<<"Bkg relative error evaluation: from histo: "<<intBerr/intB<<endl;
1333
1334   cout<<"Last estimation of bkg error is used"<<endl;
1335
1336   //backround +/- error in the range
1337   if (ftypeOfFit4Bkg == 3 && ftypeOfFit4Sgn == 0) {
1338     background = 0;
1339     errbackground = 0;
1340   }
1341   else{
1342     background=funcbkg->Integral(min,max)/(Double_t)fhistoInvMass->GetBinWidth(2);
1343     errbackground=intBerr/intB*background; // assume relative error is the same as for total integral
1344     //cout<<"integral = "<<funcbkg->Integral(min, max)<<"\tbinW = "<<fhistoInvMass->GetBinWidth(2)<<endl;
1345   }
1346   return;
1347
1348 }
1349
1350
1351 //__________________________________________________________________________
1352
1353 void AliHFMassFitter::Significance(Double_t nOfSigma,Double_t &significance,Double_t &errsignificance) const  {
1354   // Return significance in mean+- n sigma
1355
1356   Double_t signal,errsignal,background,errbackground;
1357   Signal(nOfSigma,signal,errsignal);
1358   Background(nOfSigma,background,errbackground);
1359
1360   significance =  signal/TMath::Sqrt(signal+background);
1361
1362   errsignificance = TMath::Sqrt(significance*significance/(signal+background)/(signal+background)*(1/4.*errsignal*errsignal+errbackground*errbackground)+significance*significance/signal/signal*errsignal*errsignal);
1363
1364   return;
1365 }
1366
1367 //__________________________________________________________________________
1368
1369 void AliHFMassFitter::Significance(Double_t min, Double_t max, Double_t &significance,Double_t &errsignificance) const {
1370   // Return significance integral in a range
1371
1372   Double_t signal,errsignal,background,errbackground;
1373   Signal(min, max,signal,errsignal);
1374   Background(min, max,background,errbackground);
1375
1376   significance =  signal/TMath::Sqrt(signal+background);
1377
1378   errsignificance = TMath::Sqrt(significance*significance/(signal+background)/(signal+background)*(1/4.*errsignal*errsignal+errbackground*errbackground)+significance*significance/signal/signal*errsignal*errsignal);
1379
1380   return;
1381 }