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