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 <TVirtualPad.h>
  33 #include <TGraph.h>
  34 #include <TVirtualFitter.h>
  35 #include <TMinuit.h>
  36 #include <TStyle.h>
  37 #include <TPaveText.h>
  38
  39 #include "AliHFMassFitter.h"
  40
  41
  42
  43 ClassImp(AliHFMassFitter)
  44
  45  //************** constructors
  46 AliHFMassFitter::AliHFMassFitter() :
  47   TNamed(),
  48   fhistoInvMass(0),
  49   fminMass(0),
  50   fmaxMass(0),
  51   fminBinMass(0),
  52   fmaxBinMass(0),
  53   fNbin(0),
  54   fParsSize(1),
  55   fNFinalPars(1),
  56   fFitPars(0),
  57   fWithBkg(0),
  58   ftypeOfFit4Bkg(0),
  59   ftypeOfFit4Sgn(0),
  60   ffactor(1),
  61   fntuParam(0),
  62   fMass(1.865),
  63   fSigmaSgn(0.012),
  64   fSideBands(0),
  65   fFixPar(0),
  66   fSideBandl(0),
  67   fSideBandr(0),
  68   fcounter(0),
  69   fContourGraph(0)
  70 {
  71   // default constructor
  72
  73   cout<<"Default constructor:"<<endl;
  74   cout<<"Remember to set the Histo, the Type, the FixPar"<<endl;
  75
  76 }
  77
  78 //___________________________________________________________________________
  79
  80 AliHFMassFitter::AliHFMassFitter (const TH1F *histoToFit, Double_t minvalue, Double_t maxvalue, Int_t rebin,Int_t fittypeb,Int_t fittypes):
  81  TNamed(),
  82  fhistoInvMass(0),
  83  fminMass(0),
  84  fmaxMass(0),
  85  fminBinMass(0),
  86  fmaxBinMass(0),
  87  fNbin(0),
  88  fParsSize(1),
  89  fNFinalPars(1),
  90  fFitPars(0),
  91  fWithBkg(0),
  92  ftypeOfFit4Bkg(0),
  93  ftypeOfFit4Sgn(0),
  94  ffactor(1),
  95  fntuParam(0),
  96  fMass(1.865),
  97  fSigmaSgn(0.012),
  98  fSideBands(0),
  99  fFixPar(0),
 100  fSideBandl(0),
 101  fSideBandr(0),
 102  fcounter(0),
 103  fContourGraph(0)
 104 {
 105   // standard constructor
 106
 107   fhistoInvMass= (TH1F*)histoToFit->Clone("fhistoInvMass");
 108   fhistoInvMass->SetDirectory(0);
 109   fminMass=minvalue;
 110   fmaxMass=maxvalue;
 111   if(rebin!=1) RebinMass(rebin);
 112   else fNbin=(Int_t)fhistoInvMass->GetNbinsX();
 113   CheckRangeFit();
 114   ftypeOfFit4Bkg=fittypeb;
 115   ftypeOfFit4Sgn=fittypes;
 116   if(ftypeOfFit4Bkg!=0 && ftypeOfFit4Bkg!=1 && ftypeOfFit4Bkg!=2) fWithBkg=kFALSE;
 117   else fWithBkg=kTRUE;
 118   if (!fWithBkg) cout<<"Fit Histogram of Signal only"<<endl;
 119   else  cout<<"Type of fit For Background = "<<ftypeOfFit4Bkg<<endl;
 120
 121   ComputeParSize();
 122   fFitPars=new Float_t[fParsSize];
 123
 124   SetDefaultFixParam();
 125
 126   fContourGraph=new TList();
 127   fContourGraph->SetOwner();
 128
 129 }
 130
 131
 132 AliHFMassFitter::AliHFMassFitter(const AliHFMassFitter &mfit):
 133   TNamed(),
 134   fhistoInvMass(mfit.fhistoInvMass),
 135   fminMass(mfit.fminMass),
 136   fmaxMass(mfit.fmaxMass),
 137   fminBinMass(mfit.fminBinMass),
 138   fmaxBinMass(mfit.fmaxBinMass),
 139   fNbin(mfit.fNbin),
 140   fParsSize(mfit.fParsSize),
 141   fNFinalPars(mfit.fNFinalPars),
 142   fFitPars(0),
 143   fWithBkg(mfit.fWithBkg),
 144   ftypeOfFit4Bkg(mfit.ftypeOfFit4Bkg),
 145   ftypeOfFit4Sgn(mfit.ftypeOfFit4Sgn),
 146   ffactor(mfit.ffactor),
 147   fntuParam(mfit.fntuParam),
 148   fMass(mfit.fMass),
 149   fSigmaSgn(mfit.fSigmaSgn),
 150   fSideBands(mfit.fSideBands),
 151   fFixPar(0),
 152   fSideBandl(mfit.fSideBandl),
 153   fSideBandr(mfit.fSideBandr),
 154   fcounter(mfit.fcounter),
 155   fContourGraph(mfit.fContourGraph)
 156 {
 157   //copy constructor
 158
 159   if(mfit.fParsSize > 0){
 160     fFitPars=new Float_t[fParsSize];
 161     fFixPar=new Bool_t[fNFinalPars];
 162     memcpy(fFitPars,mfit.fFitPars,mfit.fParsSize*sizeof(Float_t));
 163     memcpy(fFixPar,mfit.fFixPar,mfit.fNFinalPars*sizeof(Bool_t));
 164   }
 165
 166 }
 167
 168 //_________________________________________________________________________
 169
 170 AliHFMassFitter::~AliHFMassFitter() {
 171
 172   //destructor
 173
 174   cout<<"AliHFMassFitter destructor called"<<endl;
 175   if(fhistoInvMass) {
 176     cout<<"deleting histogram..."<<endl;
 177     delete fhistoInvMass;
 178     fhistoInvMass=NULL;
 179   }
 180   if(fntuParam){
 181     cout<<"deleting ntuple..."<<endl;
 182     delete fntuParam;
 183
 184     fntuParam=NULL;
 185   }
 186
 187   if(fFitPars) {
 188     delete[] fFitPars;
 189     cout<<"deleting parameter array..."<<endl;
 190     fFitPars=NULL;
 191   }
 192
 193   if(fFixPar) {
 194     delete[] fFixPar;
 195     cout<<"deleting bool array..."<<endl;
 196     fFixPar=NULL;
 197   }
 198
 199   fcounter = 0;
 200 }
 201
 202 //_________________________________________________________________________
 203
 204 AliHFMassFitter& AliHFMassFitter::operator=(const AliHFMassFitter &mfit){
 205
 206   //assignment operator
 207
 208   if(&mfit == this) return *this;
 209   fhistoInvMass= mfit.fhistoInvMass;
 210   fminMass= mfit.fminMass;
 211   fmaxMass= mfit.fmaxMass;
 212   fNbin= mfit.fNbin;
 213   fParsSize= mfit.fParsSize;
 214   fWithBkg= mfit.fWithBkg;
 215   ftypeOfFit4Bkg= mfit.ftypeOfFit4Bkg;
 216   ftypeOfFit4Sgn= mfit.ftypeOfFit4Sgn;
 217   ffactor= mfit.ffactor;
 218   fntuParam= mfit.fntuParam;
 219   fMass= mfit.fMass;
 220   fSigmaSgn= mfit.fSigmaSgn;
 221   fSideBands= mfit.fSideBands;
 222   fSideBandl= mfit.fSideBandl;
 223   fSideBandr= mfit.fSideBandr;
 224   fcounter= mfit.fcounter;
 225   fContourGraph= mfit.fContourGraph;
 226
 227   if(mfit.fParsSize > 0){
 228     if(fFitPars) {
 229       delete[] fFitPars;
 230       fFitPars=NULL;
 231     }
 232     fFitPars=new Float_t[fParsSize];
 233     memcpy(fFitPars,mfit.fFitPars,mfit.fParsSize*sizeof(Float_t));
 234
 235     if(fFixPar) {
 236       delete[] fFixPar;
 237       fFixPar=NULL;
 238     }
 239     fFixPar=new Bool_t[fNFinalPars];
 240     memcpy(fFixPar,mfit.fFixPar,mfit.fNFinalPars*sizeof(Float_t));
 241   }
 242
 243   return *this;
 244 }
 245
 246 //************ tools & settings
 247
 248 //__________________________________________________________________________
 249
 250 void AliHFMassFitter::ComputeNFinalPars() {
 251
 252   //compute the number of parameters of the total (signal+bgk) function
 253   cout<<"Info:ComputeNFinalPars... ";
 254   switch (ftypeOfFit4Bkg) {//npar background func
 255   case 0:
 256     fNFinalPars=2;
 257     break;
 258   case 1:
 259     fNFinalPars=2;
 260     break;
 261   case 2:
 262     fNFinalPars=3;
 263     break;
 264   case 3:
 265     fNFinalPars=1;
 266     break;
 267   default:
 268     cout<<"Error in computing fNFinalPars: check ftypeOfFit4Bkg"<<endl;
 269     break;
 270   }
 271
 272   fNFinalPars+=3; //gaussian signal
 273   cout<<": "<<fNFinalPars<<endl;
 274 }
 275 //__________________________________________________________________________
 276
 277 void AliHFMassFitter::ComputeParSize() {
 278
 279   //compute the size of the parameter array and set the data member
 280
 281   switch (ftypeOfFit4Bkg) {//npar background func
 282   case 0:
 283     fParsSize = 2*3;
 284     break;
 285   case 1:
 286     fParsSize = 2*3;
 287     break;
 288   case 2:
 289     fParsSize = 3*3;
 290     break;
 291   case 3:
 292     fParsSize = 1*3;
 293     break;
 294   default:
 295     cout<<"Error in computing fParsSize: check ftypeOfFit4Bkg"<<endl;
 296     break;
 297   }
 298
 299   fParsSize += 3; // npar refl
 300   fParsSize += 3; // npar signal gaus
 301
 302   fParsSize*=2;   // add errors
 303   cout<<"Parameters array size "<<fParsSize<<endl;
 304 }
 305
 306 //___________________________________________________________________________
 307 void AliHFMassFitter::SetDefaultFixParam(){
 308
 309   //Set default values for fFixPar (only total integral fixed)
 310
 311   ComputeNFinalPars();
 312   fFixPar=new Bool_t[fNFinalPars];
 313
 314   fFixPar[0]=kTRUE; //default: IntTot fixed
 315   cout<<"Parameter 0 is fixed"<<endl;
 316   for(Int_t i=1;i<fNFinalPars;i++){
 317     fFixPar[i]=kFALSE;
 318   }
 319
 320 }
 321
 322 //___________________________________________________________________________
 323 Bool_t AliHFMassFitter::SetFixThisParam(Int_t thispar,Bool_t fixpar){
 324
 325   //set the value (kFALSE or kTRUE) of one element of fFixPar
 326   //return kFALSE if something wrong
 327
 328   if(thispar>=fNFinalPars) {
 329     cout<<"Error! Parameter out of bounds! Max is "<<fNFinalPars-1<<endl;
 330     return kFALSE;
 331   }
 332   if(!fFixPar){
 333     cout<<"Initializing fFixPar...";
 334     SetDefaultFixParam();
 335     cout<<" done."<<endl;
 336   }
 337
 338   fFixPar[thispar]=fixpar;
 339   if(fixpar)cout<<"Parameter "<<thispar<<" is now fixed"<<endl;
 340   else cout<<"Parameter "<<thispar<<" is now free"<<endl;
 341   return kTRUE;
 342 }
 343
 344 //___________________________________________________________________________
 345 Bool_t AliHFMassFitter::GetFixThisParam(Int_t thispar)const{
 346   //return the value of fFixPar[thispar]
 347   if(thispar>=fNFinalPars) {
 348     cout<<"Error! Parameter out of bounds! Max is "<<fNFinalPars-1<<endl;
 349     return kFALSE;
 350   }
 351   if(!fFixPar) {
 352     cout<<"Error! Parameters to be fixed still not set"<<endl;
 353     return kFALSE;
 354
 355   }
 356   return fFixPar[thispar];
 357
 358 }
 359
 360 //___________________________________________________________________________
 361 void AliHFMassFitter::SetHisto(const TH1F *histoToFit){
 362
 363   fhistoInvMass = new TH1F(*histoToFit);
 364   fhistoInvMass->SetDirectory(0);
 365   //cout<<"SetHisto pointer "<<fhistoInvMass<<endl;
 366 }
 367
 368 //___________________________________________________________________________
 369
 370 void AliHFMassFitter::SetType(Int_t fittypeb, Int_t fittypes) {
 371
 372   //set the type of fit to perform for signal and background
 373
 374   ftypeOfFit4Bkg = fittypeb;
 375   ftypeOfFit4Sgn = fittypes;
 376
 377   ComputeParSize();
 378   fFitPars = new Float_t[fParsSize];
 379
 380   SetDefaultFixParam();
 381 }
 382
 383 //___________________________________________________________________________
 384
 385 void AliHFMassFitter::Reset() {
 386
 387   //delete the histogram and reset the mean and sigma to default
 388
 389   cout<<"Reset called: delete histogram, set mean value to 1.85 and sigma to 0.012"<<endl;
 390   fMass=1.85;
 391   fSigmaSgn=0.012;
 392   cout<<"Reset "<<fhistoInvMass<<endl;
 393   if(fhistoInvMass) {
 394     delete fhistoInvMass;
 395     fhistoInvMass=NULL;
 396     cout<<fhistoInvMass<<endl;
 397   }
 398   else cout<<"histogram doesn't exist, do not delete"<<endl;
 399
 400
 401 }
 402
 403 //_________________________________________________________________________
 404
 405 void AliHFMassFitter::InitNtuParam(char *ntuname) {
 406
 407   // Create ntuple to keep fit parameters
 408
 409   fntuParam=0;
 410   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");
 411
 412 }
 413
 414 //_________________________________________________________________________
 415
 416 void AliHFMassFitter::FillNtuParam() {
 417   // Fill ntuple with fit parameters
 418
 419   Float_t nothing=0.;
 420
 421   if (ftypeOfFit4Bkg==2) {
 422       fntuParam->SetBranchAddress("intbkg1",&fFitPars[0]);
 423       fntuParam->SetBranchAddress("slope1",&fFitPars[1]);
 424       fntuParam->SetBranchAddress("conc1",&fFitPars[2]);
 425       fntuParam->SetBranchAddress("intGB",&fFitPars[3]);
 426       fntuParam->SetBranchAddress("meanGB",&fFitPars[4]);
 427       fntuParam->SetBranchAddress("sigmaGB",&fFitPars[5]);
 428       fntuParam->SetBranchAddress("intbkg2",&fFitPars[6]);
 429       fntuParam->SetBranchAddress("slope2",&fFitPars[7]);
 430       fntuParam->SetBranchAddress("conc2",&fFitPars[8]);
 431       fntuParam->SetBranchAddress("inttot",&fFitPars[9]);
 432       fntuParam->SetBranchAddress("slope3",&fFitPars[10]);
 433       fntuParam->SetBranchAddress("conc3",&fFitPars[11]);
 434       fntuParam->SetBranchAddress("intsgn",&fFitPars[12]);
 435       fntuParam->SetBranchAddress("meansgn",&fFitPars[13]);
 436       fntuParam->SetBranchAddress("sigmasgn",&fFitPars[14]);
 437
 438       fntuParam->SetBranchAddress("intbkg1Err",&fFitPars[15]);
 439       fntuParam->SetBranchAddress("slope1Err",&fFitPars[16]);
 440       fntuParam->SetBranchAddress("conc1Err",&fFitPars[17]);
 441       fntuParam->SetBranchAddress("intGBErr",&fFitPars[18]);
 442       fntuParam->SetBranchAddress("meanGBErr",&fFitPars[19]);
 443       fntuParam->SetBranchAddress("sigmaGBErr",&fFitPars[20]);
 444       fntuParam->SetBranchAddress("intbkg2Err",&fFitPars[21]);
 445       fntuParam->SetBranchAddress("slope2Err",&fFitPars[22]);
 446       fntuParam->SetBranchAddress("conc2Err",&fFitPars[23]);
 447       fntuParam->SetBranchAddress("inttotErr",&fFitPars[24]);
 448       fntuParam->SetBranchAddress("slope3Err",&fFitPars[25]);
 449       fntuParam->SetBranchAddress("conc3Err",&fFitPars[26]);
 450       fntuParam->SetBranchAddress("intsgnErr",&fFitPars[27]);
 451       fntuParam->SetBranchAddress("meansgnErr",&fFitPars[28]);
 452       fntuParam->SetBranchAddress("sigmasgnErr",&fFitPars[29]);
 453
 454   } else {
 455
 456     if(ftypeOfFit4Bkg==3){
 457       fntuParam->SetBranchAddress("intbkg1",&fFitPars[0]);
 458       fntuParam->SetBranchAddress("slope1",&nothing);
 459       fntuParam->SetBranchAddress("conc1",&nothing);
 460       fntuParam->SetBranchAddress("intGB",&fFitPars[1]);
 461       fntuParam->SetBranchAddress("meanGB",&fFitPars[2]);
 462       fntuParam->SetBranchAddress("sigmaGB",&fFitPars[3]);
 463       fntuParam->SetBranchAddress("intbkg2",&fFitPars[4]);
 464       fntuParam->SetBranchAddress("slope2",&nothing);
 465       fntuParam->SetBranchAddress("conc2",&nothing);
 466       fntuParam->SetBranchAddress("inttot",&fFitPars[6]);
 467       fntuParam->SetBranchAddress("slope3",&nothing);
 468       fntuParam->SetBranchAddress("conc3",&nothing);
 469       fntuParam->SetBranchAddress("intsgn",&fFitPars[6]);
 470       fntuParam->SetBranchAddress("meansgn",&fFitPars[7]);
 471       fntuParam->SetBranchAddress("sigmasgn",&fFitPars[8]);
 472
 473       fntuParam->SetBranchAddress("intbkg1Err",&fFitPars[9]);
 474       fntuParam->SetBranchAddress("slope1Err",&nothing);
 475       fntuParam->SetBranchAddress("conc1Err",&nothing);
 476       fntuParam->SetBranchAddress("intGBErr",&fFitPars[10]);
 477       fntuParam->SetBranchAddress("meanGBErr",&fFitPars[11]);
 478       fntuParam->SetBranchAddress("sigmaGBErr",&fFitPars[12]);
 479       fntuParam->SetBranchAddress("intbkg2Err",&fFitPars[13]);
 480       fntuParam->SetBranchAddress("slope2Err",&nothing);
 481       fntuParam->SetBranchAddress("conc2Err",&nothing);
 482       fntuParam->SetBranchAddress("inttotErr",&fFitPars[15]);
 483       fntuParam->SetBranchAddress("slope3Err",&nothing);
 484       fntuParam->SetBranchAddress("conc3Err",&nothing);
 485       fntuParam->SetBranchAddress("intsgnErr",&fFitPars[15]);
 486       fntuParam->SetBranchAddress("meansgnErr",&fFitPars[16]);
 487       fntuParam->SetBranchAddress("sigmasgnErr",&fFitPars[17]);
 488
 489     }
 490     else{
 491       fntuParam->SetBranchAddress("intbkg1",&fFitPars[0]);
 492       fntuParam->SetBranchAddress("slope1",&fFitPars[1]);
 493       fntuParam->SetBranchAddress("conc1",&nothing);
 494       fntuParam->SetBranchAddress("intGB",&fFitPars[2]);
 495       fntuParam->SetBranchAddress("meanGB",&fFitPars[3]);
 496       fntuParam->SetBranchAddress("sigmaGB",&fFitPars[4]);
 497       fntuParam->SetBranchAddress("intbkg2",&fFitPars[5]);
 498       fntuParam->SetBranchAddress("slope2",&fFitPars[6]);
 499       fntuParam->SetBranchAddress("conc2",&nothing);
 500       fntuParam->SetBranchAddress("inttot",&fFitPars[7]);
 501       fntuParam->SetBranchAddress("slope3",&fFitPars[8]);
 502       fntuParam->SetBranchAddress("conc3",&nothing);
 503       fntuParam->SetBranchAddress("intsgn",&fFitPars[9]);
 504       fntuParam->SetBranchAddress("meansgn",&fFitPars[10]);
 505       fntuParam->SetBranchAddress("sigmasgn",&fFitPars[11]);
 506
 507       fntuParam->SetBranchAddress("intbkg1Err",&fFitPars[12]);
 508       fntuParam->SetBranchAddress("slope1Err",&fFitPars[13]);
 509       fntuParam->SetBranchAddress("conc1Err",&nothing);
 510       fntuParam->SetBranchAddress("intGBErr",&fFitPars[14]);
 511       fntuParam->SetBranchAddress("meanGBErr",&fFitPars[15]);
 512       fntuParam->SetBranchAddress("sigmaGBErr",&fFitPars[16]);
 513       fntuParam->SetBranchAddress("intbkg2Err",&fFitPars[17]);
 514       fntuParam->SetBranchAddress("slope2Err",&fFitPars[18]);
 515       fntuParam->SetBranchAddress("conc2Err",&nothing);
 516       fntuParam->SetBranchAddress("inttotErr",&fFitPars[19]);
 517       fntuParam->SetBranchAddress("slope3Err",&fFitPars[20]);
 518       fntuParam->SetBranchAddress("conc3Err",&nothing);
 519       fntuParam->SetBranchAddress("intsgnErr",&fFitPars[21]);
 520       fntuParam->SetBranchAddress("meansgnErr",&fFitPars[22]);
 521       fntuParam->SetBranchAddress("sigmasgnErr",&fFitPars[23]);
 522     }
 523
 524   }
 525   fntuParam->TTree::Fill();
 526 }
 527
 528 //_________________________________________________________________________
 529
 530 TNtuple* AliHFMassFitter::NtuParamOneShot(char *ntuname){
 531   // Create, fill and return ntuple with fit parameters
 532
 533   InitNtuParam(ntuname);
 534   FillNtuParam();
 535   return fntuParam;
 536 }
 537 //_________________________________________________________________________
 538
 539 void AliHFMassFitter::RebinMass(Int_t bingroup){
 540   // Rebin invariant mass histogram
 541
 542   if(!fhistoInvMass){
 543     cout<<"Histogram not set"<<endl;
 544     return;
 545   }
 546   Int_t nbinshisto=fhistoInvMass->GetNbinsX();
 547   if(bingroup<1){
 548     cout<<"Error! Cannot group "<<bingroup<<" bins\n";
 549     fNbin=nbinshisto;
 550     cout<<"Kept original number of bins: "<<fNbin<<endl;
 551   } else{
 552
 553     while(nbinshisto%bingroup != 0) {
 554       bingroup--;
 555     }
 556     cout<<"Group "<<bingroup<<" bins"<<endl;
 557     fhistoInvMass->Rebin(bingroup);
 558     fNbin = fhistoInvMass->GetNbinsX();
 559     cout<<"New number of bins: "<<fNbin<<endl;
 560   }
 561
 562 }
 563
 564 //************* fit
 565
 566 //___________________________________________________________________________
 567
 568 Double_t AliHFMassFitter::FitFunction4MassDistr (Double_t *x, Double_t *par){
 569   // Fit function for signal+background
 570
 571
 572   //exponential or linear fit
 573   //
 574   // par[0] = tot integral
 575   // par[1] = slope
 576   // par[2] = gaussian integral
 577   // par[3] = gaussian mean
 578   // par[4] = gaussian sigma
 579
 580   Double_t total,bkg=0,sgn=0;
 581
 582   if (ftypeOfFit4Bkg==0 || ftypeOfFit4Bkg==1) {
 583     if(ftypeOfFit4Sgn == 0) {
 584
 585       Double_t parbkg[2] = {par[0]-par[2], par[1]};
 586       bkg = FitFunction4Bkg(x,parbkg);
 587     }
 588     if(ftypeOfFit4Sgn == 1) {
 589       Double_t parbkg[5] = {par[2],par[3],ffactor*par[4],par[0]-2*par[2], par[1]};
 590       bkg = FitFunction4Bkg(x,parbkg);
 591     }
 592
 593     sgn = FitFunction4Sgn(x,&par[2]);
 594
 595   }
 596
 597   //polynomial fit
 598
 599     // par[0] = tot integral
 600     // par[1] = coef1
 601     // par[2] = coef2
 602     // par[3] = gaussian integral
 603     // par[4] = gaussian mean
 604     // par[5] = gaussian sigma
 605
 606   if (ftypeOfFit4Bkg==2) {
 607
 608     if(ftypeOfFit4Sgn == 0) {
 609
 610       Double_t parbkg[3] = {par[0]-par[3], par[1], par[2]};
 611       bkg = FitFunction4Bkg(x,parbkg);
 612     }
 613     if(ftypeOfFit4Sgn == 1) {
 614
 615       Double_t parbkg[6] = {par[3],par[4],ffactor*par[5],par[0]-2*par[3], par[1], par[2]};
 616      bkg = FitFunction4Bkg(x,parbkg);
 617     }
 618
 619     sgn = FitFunction4Sgn(x,&par[3]);
 620   }
 621
 622   if (ftypeOfFit4Bkg==3) {
 623
 624     if(ftypeOfFit4Sgn == 0) {
 625         bkg=FitFunction4Bkg(x,par);
 626         sgn=FitFunction4Sgn(x,&par[1]);
 627     }
 628     if(ftypeOfFit4Sgn == 1) {
 629       Double_t parbkg[4]={par[1],par[2],ffactor*par[3],par[0]};
 630       bkg=FitFunction4Bkg(x,parbkg);
 631       sgn=FitFunction4Sgn(x,&par[1]);
 632     }
 633   }
 634
 635   total = bkg + sgn;
 636
 637   return  total;
 638 }
 639
 640 //_________________________________________________________________________
 641 Double_t AliHFMassFitter::FitFunction4Sgn (Double_t *x, Double_t *par){
 642   // Fit function for the signal
 643
 644   //gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
 645   //Par:
 646   // * [0] = integralSgn
 647   // * [1] = mean
 648   // * [2] = sigma
 649   //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
 650
 651   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]);
 652
 653 }
 654
 655 //__________________________________________________________________________
 656
 657 Double_t AliHFMassFitter::FitFunction4Bkg (Double_t *x, Double_t *par){
 658   // Fit function for the background
 659
 660   Double_t maxDeltaM = 4.*fSigmaSgn;
 661   if(fSideBands && TMath::Abs(x[0]-fMass) < maxDeltaM) {
 662     TF1::RejectPoint();
 663     return 0;
 664   }
 665   Int_t firstPar=0;
 666   Double_t gaus2=0,total=-1;
 667   if(ftypeOfFit4Sgn == 1){
 668     firstPar=3;
 669     //gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
 670     //Par:
 671     // * [0] = integralSgn
 672     // * [1] = mean
 673     // * [2] = sigma
 674     //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
 675     gaus2 = FitFunction4Sgn(x,par);
 676   }
 677
 678   switch (ftypeOfFit4Bkg){
 679   case 0:
 680     //exponential
 681     //exponential = A*exp(B*x) -> integral(exponential)=A/B*exp(B*x)](min,max)
 682     //-> A = B*integral/(exp(B*max)-exp(B*min)) where integral can be written
 683     //as integralTot- integralGaus (=par [2])
 684     //Par:
 685     // * [0] = integralBkg;
 686     // * [1] = B;
 687     //exponential = [1]*[0]/(exp([1]*max)-exp([1]*min))*exp([1]*x)
 688     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]);
 689     break;
 690   case 1:
 691     //linear
 692     //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)
 693     // * [0] = integralBkg;
 694     // * [1] = b;
 695     total= par[0+firstPar]/(fmaxMass-fminMass)+par[1+firstPar]*(x[0]-0.5*(fmaxMass+fminMass));
 696     break;
 697   case 2:
 698     //polynomial
 699     //y=a+b*x+c*x**2 -> integral = a(max-min) + 1/2*b*(max^2-min^2) +
 700     //+ 1/3*c*(max^3-min^3) ->
 701     //a = (integral-1/2*b*(max^2-min^2)-1/3*c*(max^3-min^3))/(max-min)
 702     // * [0] = integralBkg;
 703     // * [1] = b;
 704     // * [2] = c;
 705     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));
 706     break;
 707   case 3:
 708     total=par[0+firstPar];
 709     break;
 710 //   default:
 711 //     Types of Fit Functions for Background:
 712 //     * 0 = exponential;
 713 //     * 1 = linear;
 714 //     * 2 = polynomial 2nd order
 715 //     * 3 = no background"<<endl;
 716
 717   }
 718   return total+gaus2;
 719 }
 720
 721 //__________________________________________________________________________
 722 Bool_t AliHFMassFitter::SideBandsBounds(){
 723
 724   //determines the ranges of the side bands
 725
 726   if (fNbin==0) fNbin=fhistoInvMass->GetNbinsX();
 727   Double_t minHisto=fhistoInvMass->GetBinLowEdge(1);
 728   Double_t maxHisto=fhistoInvMass->GetBinLowEdge(fNbin+1);
 729
 730   Double_t sidebandldouble,sidebandrdouble;
 731   Bool_t leftok=kFALSE, rightok=kFALSE;
 732
 733   if(fMass-fminMass < 0 || fmaxMass-fMass <0) {
 734     cout<<"Left limit of range > mean or right limit of range < mean: change left/right limit or initial mean value"<<endl;
 735     return kFALSE;
 736   }
 737
 738   //histo limit = fit function limit
 739   if((TMath::Abs(fminMass-minHisto) < 1e6 || TMath::Abs(fmaxMass - maxHisto) < 1e6) && (fMass-4.*fSigmaSgn-fminMass) < 1e6){
 740     Double_t coeff = (fMass-fminMass)/fSigmaSgn;
 741     sidebandldouble=(fMass-0.5*coeff*fSigmaSgn);
 742     sidebandrdouble=(fMass+0.5*coeff*fSigmaSgn);
 743     cout<<"Changed number of sigma from 4 to "<<0.5*coeff<<" for the estimation of the side bands"<<endl;
 744     if (coeff<3) cout<<"Side bands inside 3 sigma, may be better use ftypeOfFit4Bkg = 3 (only signal)"<<endl;
 745     if (coeff<2) {
 746       cout<<"Side bands inside 2 sigma. Change mode: ftypeOfFit4Bkg = 3"<<endl;
 747       ftypeOfFit4Bkg=3;
 748       //set binleft and right without considering SetRangeFit- anyway no bkg!
 749       sidebandldouble=(fMass-4.*fSigmaSgn);
 750       sidebandrdouble=(fMass+4.*fSigmaSgn);
 751     }
 752   }
 753   else {
 754     sidebandldouble=(fMass-4.*fSigmaSgn);
 755     sidebandrdouble=(fMass+4.*fSigmaSgn);
 756   }
 757
 758   cout<<"Left side band ";
 759   Double_t tmp=0.;
 760   tmp=sidebandldouble;
 761   //calculate bin corresponding to fSideBandl
 762   fSideBandl=fhistoInvMass->FindBin(sidebandldouble);
 763   if (sidebandldouble >= fhistoInvMass->GetBinCenter(fSideBandl)) fSideBandl++;
 764   sidebandldouble=fhistoInvMass->GetBinLowEdge(fSideBandl);
 765
 766   if(TMath::Abs(tmp-sidebandldouble) > 1e-6){
 767     cout<<tmp<<" is not allowed, changing it to the nearest value allowed: ";
 768     leftok=kTRUE;
 769   }
 770   cout<<sidebandldouble<<" (bin "<<fSideBandl<<")"<<endl;
 771
 772   cout<<"Right side band ";
 773   tmp=sidebandrdouble;
 774   //calculate bin corresponding to fSideBandr
 775   fSideBandr=fhistoInvMass->FindBin(sidebandrdouble);
 776   if (sidebandrdouble < fhistoInvMass->GetBinCenter(fSideBandr)) fSideBandr--;
 777   sidebandrdouble=fhistoInvMass->GetBinLowEdge(fSideBandr+1);
 778
 779   if(TMath::Abs(tmp-sidebandrdouble) > 1e-6){
 780     cout<<tmp<<" is not allowed, changing it to the nearest value allowed: ";
 781     rightok=kTRUE;
 782   }
 783   cout<<sidebandrdouble<<" (bin "<<fSideBandr<<")"<<endl;
 784   if (fSideBandl==0 || fSideBandr==fNbin) {
 785     cout<<"Error! Range too little";
 786     return kFALSE;
 787   }
 788   return kTRUE;
 789 }
 790
 791 //__________________________________________________________________________
 792
 793 void AliHFMassFitter::GetSideBandsBounds(Int_t &left, Int_t &right) const{
 794
 795   // get the range of the side bands
 796
 797   if (fSideBandl==0 && fSideBandr==0){
 798     cout<<"Use MassFitter method first"<<endl;
 799     return;
 800   }
 801   left=fSideBandl;
 802   right=fSideBandr;
 803 }
 804
 805 //__________________________________________________________________________
 806 Bool_t AliHFMassFitter::CheckRangeFit(){
 807   //check if the limit of the range correspond to the limit of bins. If not reset the limit to the nearer value which satisfy this condition
 808
 809   if (!fhistoInvMass) {
 810     cout<<"No histogram to fit! SetHisto(TH1F* h) before! "<<endl;
 811     return kFALSE;
 812   }
 813   Bool_t leftok=kFALSE, rightok=kFALSE;
 814   Int_t nbins=fhistoInvMass->GetNbinsX();
 815   Double_t minhisto=fhistoInvMass->GetBinLowEdge(1), maxhisto=fhistoInvMass->GetBinLowEdge(nbins+1);
 816
 817   //check if limits are inside histogram range
 818
 819   if( fminMass-minhisto < 0. ) {
 820     cout<<"Out of histogram left bound! Setting to "<<minhisto<<endl;
 821     fminMass=minhisto;
 822   }
 823   if( fmaxMass-maxhisto > 0. ) {
 824     cout<<"Out of histogram right bound! Setting to"<<maxhisto<<endl;
 825     fmaxMass=maxhisto;
 826   }
 827
 828   Double_t tmp=0.;
 829   tmp=fminMass;
 830   //calculate bin corresponding to fminMass
 831   fminBinMass=fhistoInvMass->FindBin(fminMass);
 832   if (fminMass >= fhistoInvMass->GetBinCenter(fminBinMass)) fminBinMass++;
 833   fminMass=fhistoInvMass->GetBinLowEdge(fminBinMass);
 834   if(TMath::Abs(tmp-fminMass) > 1e-6){
 835     cout<<"Left bound "<<tmp<<" is not allowed, changing it to the nearest value allowed: "<<fminMass<<endl;
 836     leftok=kTRUE;
 837   }
 838
 839   tmp=fmaxMass;
 840   //calculate bin corresponding to fmaxMass
 841   fmaxBinMass=fhistoInvMass->FindBin(fmaxMass);
 842   if (fmaxMass < fhistoInvMass->GetBinCenter(fmaxBinMass)) fmaxBinMass--;
 843   fmaxMass=fhistoInvMass->GetBinLowEdge(fmaxBinMass+1);
 844   if(TMath::Abs(tmp-fmaxMass) > 1e-6){
 845     cout<<"Right bound "<<tmp<<" is not allowed, changing it to the nearest value allowed: "<<fmaxMass<<endl;
 846     rightok=kTRUE;
 847   }
 848
 849   return (leftok && rightok);
 850
 851 }
 852
 853 //__________________________________________________________________________
 854
 855 Bool_t AliHFMassFitter::MassFitter(Bool_t draw){
 856   // Main method of the class: performs the fit of the histogram
 857
 858   //Set default fitter Minuit in order to use gMinuit in the contour plots
 859   TVirtualFitter::SetDefaultFitter("Minuit");
 860
 861
 862   Bool_t isBkgOnly=kFALSE;
 863
 864   Int_t fit1status=RefitWithBkgOnly(kFALSE);
 865   if(fit1status){
 866     Int_t checkinnsigma=4;
 867     Double_t range[2]={fMass-checkinnsigma*fSigmaSgn,fMass+checkinnsigma*fSigmaSgn};
 868     TF1* func=GetHistoClone()->GetFunction("funcbkgonly");
 869     Double_t intUnderFunc=func->Integral(range[0],range[1]);
 870     Double_t intUnderHisto=fhistoInvMass->Integral(fhistoInvMass->FindBin(range[0]),fhistoInvMass->FindBin(range[1]),"width");
 871     cout<<"Pick zone: IntFunc = "<<intUnderFunc<<"; IntHist = "<<intUnderHisto<<"\tDiff = "<<intUnderHisto-intUnderFunc<<"\tRelDiff = "<<(intUnderHisto-intUnderFunc)/intUnderFunc<<endl;
 872     Double_t diffUnderPick=(intUnderHisto-intUnderFunc);
 873     intUnderFunc=func->Integral(fminMass,fminMass+checkinnsigma*fSigmaSgn);
 874     intUnderHisto=fhistoInvMass->Integral(fhistoInvMass->FindBin(fminMass),fhistoInvMass->FindBin(fminMass+checkinnsigma*fSigmaSgn),"width");
 875     cout<<"Band (l) zone: IntFunc = "<<intUnderFunc<<"; IntHist = "<<intUnderHisto<<"\tDiff = "<<intUnderHisto-intUnderFunc<<"\tRelDiff = "<<(intUnderHisto-intUnderFunc)/intUnderFunc<<endl;
 876     Double_t diffUnderBands=(intUnderHisto-intUnderFunc);
 877     Double_t relDiff=diffUnderPick/diffUnderBands;
 878     cout<<"Relative difference = "<<relDiff<<endl;
 879     if(TMath::Abs(relDiff) < 1) isBkgOnly=kTRUE;
 880     else{
 881       cout<<"Relative difference = "<<relDiff<<": I suppose there is some signal, continue with total fit!"<<endl;
 882     }
 883   }
 884   if(isBkgOnly) {
 885     cout<<"INFO!! The histogram contains only background"<<endl;
 886     if(draw)DrawFit();
 887
 888     //increase counter of number of fits done
 889     fcounter++;
 890
 891     return kTRUE;
 892   }
 893
 894   Int_t bkgPar = fNFinalPars-3; //background function's number of parameters
 895
 896   cout<<"fNFinalPars = "<<fNFinalPars<<"\tbkgPar = "<<bkgPar<<endl;
 897
 898
 899   TString listname="contourplot";
 900   listname+=fcounter;
 901   if(!fContourGraph){
 902     fContourGraph=new TList();
 903     fContourGraph->SetOwner();
 904   }
 905
 906   fContourGraph->SetName(listname);
 907
 908
 909   //function names
 910   TString bkgname="funcbkg";
 911   TString bkg1name="funcbkg1";
 912   TString massname="funcmass";
 913
 914   //Total integral
 915   Double_t totInt = fhistoInvMass->Integral(fminBinMass,fmaxBinMass, "width");
 916   //cout<<"Here tot integral is = "<<totInt<<"; integral in whole range is "<<fhistoInvMass->Integral("width")<<endl;
 917   fSideBands = kTRUE;
 918   Double_t width=fhistoInvMass->GetBinWidth(8);
 919   //cout<<"fNbin = "<<fNbin<<endl;
 920   if (fNbin==0) fNbin=fhistoInvMass->GetNbinsX();
 921
 922   Bool_t ok=SideBandsBounds();
 923   if(!ok) return kFALSE;
 924
 925   //sidebands integral - first approx (from histo)
 926   Double_t sideBandsInt=(Double_t)fhistoInvMass->Integral(1,fSideBandl,"width") + (Double_t)fhistoInvMass->Integral(fSideBandr,fNbin,"width");
 927   cout<<"------nbin = "<<fNbin<<"\twidth = "<<width<<"\tbinleft = "<<fSideBandl<<"\tbinright = "<<fSideBandr<<endl;
 928   cout<<"------sideBandsInt - first approx = "<<sideBandsInt<<endl;
 929   if (sideBandsInt<=0) {
 930     cout<<"! sideBandsInt <=0. There's a problem, cannot start the fit"<<endl;
 931     return kFALSE;
 932   }
 933
 934   /*Fit Bkg*/
 935
 936
 937   TF1 *funcbkg = new TF1(bkgname.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,bkgPar,"AliHFMassFitter","FitFunction4Bkg");
 938   cout<<"Function name = "<<funcbkg->GetName()<<endl<<endl;
 939
 940   funcbkg->SetLineColor(2); //red
 941
 942   //first fit for bkg: approx bkgint
 943
 944   switch (ftypeOfFit4Bkg) {
 945   case 0: //gaus+expo
 946     funcbkg->SetParNames("BkgInt","Slope");
 947     funcbkg->SetParameters(sideBandsInt,-2.);
 948     break;
 949   case 1:
 950     funcbkg->SetParNames("BkgInt","Slope");
 951     funcbkg->SetParameters(sideBandsInt,-100.);
 952     break;
 953   case 2:
 954     funcbkg->SetParNames("BkgInt","Coef1","Coef2");
 955     funcbkg->SetParameters(sideBandsInt,-10.,5);
 956     break;
 957   case 3:
 958     if(ftypeOfFit4Sgn==0){
 959       funcbkg->SetParNames("Const");
 960       funcbkg->SetParameter(0,0.);
 961       funcbkg->FixParameter(0,0.);
 962     }
 963     break;
 964   default:
 965     cout<<"Wrong choise of ftypeOfFit4Bkg ("<<ftypeOfFit4Bkg<<")"<<endl;
 966     return kFALSE;
 967     break;
 968   }
 969   cout<<"\nBACKGROUND FIT - only combinatorial"<<endl;
 970   Int_t ftypeOfFit4SgnBkp=ftypeOfFit4Sgn;
 971
 972   Double_t intbkg1=0,slope1=0,conc1=0;
 973   //if only signal and reflection: skip
 974   if (!(ftypeOfFit4Bkg==3 && ftypeOfFit4Sgn==1)) {
 975     ftypeOfFit4Sgn=0;
 976     fhistoInvMass->Fit(bkgname.Data(),"R,L,E,0");
 977
 978     for(Int_t i=0;i<bkgPar;i++){
 979       fFitPars[i]=funcbkg->GetParameter(i);
 980       //cout<<i<<"\t"<<funcbkg->GetParameter(i)<<"\t";
 981       fFitPars[fNFinalPars+2*bkgPar+3+i]= funcbkg->GetParError(i);
 982       //cout<<fNFinalPars+2*bkgPar+3+i<<"\t"<< funcbkg->GetParError(i)<<endl;
 983     }
 984     fSideBands = kFALSE;
 985     //intbkg1 = funcbkg->GetParameter(0);
 986
 987     intbkg1 = funcbkg->Integral(fminMass,fmaxMass);
 988     if(ftypeOfFit4Bkg!=3) slope1 = funcbkg->GetParameter(1);
 989     if(ftypeOfFit4Bkg==2) conc1 = funcbkg->GetParameter(2);
 990     //cout<<"First fit: \nintbkg1 = "<<intbkg1<<"\t(Compare with par0 = "<<funcbkg->GetParameter(0)<<")\nslope1= "<<slope1<<"\nconc1 = "<<conc1<<endl;
 991   }
 992   else cout<<"\t\t//"<<endl;
 993
 994   ftypeOfFit4Sgn=ftypeOfFit4SgnBkp;
 995   TF1 *funcbkg1=0;
 996   if (ftypeOfFit4Sgn == 1) {
 997     cout<<"\nBACKGROUND FIT WITH REFLECTION"<<endl;
 998     bkgPar+=3;
 999
1000     //cout<<"fNFinalPars = "<<fNFinalPars<<"\tbkgPar = "<<bkgPar<<endl;
1001
1002     funcbkg1 = new TF1(bkg1name.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,bkgPar,"AliHFMassFitter","FitFunction4Bkg");
1003     cout<<"Function name = "<<funcbkg1->GetName()<<endl;
1004
1005     funcbkg1->SetLineColor(2); //red
1006
1007     if(ftypeOfFit4Bkg==2){
1008       cout<<"*** Polynomial Fit ***"<<endl;
1009       funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Coef1","Coef2");
1010       funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1,conc1);
1011
1012       //cout<<"Parameters set to: "<<0.5*(totInt-intbkg1)<<"\t"<<fMass<<"\t"<<ffactor*fSigmaSgn<<"\t"<<intbkg1<<"\t"<<slope1<<"\t"<<conc1<<"\t"<<endl;
1013       //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<bkgPar<<"\tgsidebands = "<<fSideBands<<endl;
1014     } else{
1015       if(ftypeOfFit4Bkg==3) //no background: gaus sign+ gaus broadened
1016         {
1017           cout<<"*** No background Fit ***"<<endl;
1018           funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","Const");
1019           funcbkg1->SetParameters(0.5*totInt,fMass,ffactor*fSigmaSgn,0.);
1020           funcbkg1->FixParameter(3,0.);
1021         } else{ //expo or linear
1022           if(ftypeOfFit4Bkg==0) cout<<"*** Exponential Fit ***"<<endl;
1023           if(ftypeOfFit4Bkg==1) cout<<"*** Linear Fit ***"<<endl;
1024           funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Slope");
1025           funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1);
1026         }
1027     }
1028     Int_t status=fhistoInvMass->Fit(bkg1name.Data(),"R,L,E,+,0");
1029     if (status != 0){
1030       cout<<"Minuit returned "<<status<<endl;
1031       return kFALSE;
1032     }
1033
1034     for(Int_t i=0;i<bkgPar;i++){
1035       fFitPars[bkgPar-3+i]=funcbkg1->GetParameter(i);
1036       //cout<<bkgPar-3+i<<"\t"<<funcbkg1->GetParameter(i);
1037       fFitPars[fNFinalPars+3*bkgPar-6+i]= funcbkg1->GetParError(i);
1038       //cout<<"\t"<<fNFinalPars+3*bkgPar-6+i<<"\t"<<funcbkg1->GetParError(i)<<endl;
1039     }
1040
1041     intbkg1=funcbkg1->GetParameter(3);
1042     if(ftypeOfFit4Bkg!=3) slope1 = funcbkg1->GetParameter(4);
1043     if(ftypeOfFit4Bkg==2) conc1 = funcbkg1->GetParameter(5);
1044
1045   } else {
1046     bkgPar+=3;
1047
1048     for(Int_t i=0;i<3;i++){
1049       fFitPars[bkgPar-3+i]=0.;
1050       cout<<bkgPar-3+i<<"\t"<<0.<<"\t";
1051       fFitPars[fNFinalPars+3*bkgPar-6+i]= 0.;
1052       cout<<fNFinalPars+3*bkgPar-6+i<<"\t"<<0.<<endl;
1053     }
1054
1055     for(Int_t i=0;i<bkgPar-3;i++){
1056       fFitPars[bkgPar+i]=funcbkg->GetParameter(i);
1057       cout<<bkgPar+i<<"\t"<<funcbkg->GetParameter(i)<<"\t";
1058       fFitPars[fNFinalPars+3*bkgPar-3+i]= funcbkg->GetParError(i);
1059       cout<<fNFinalPars+3*bkgPar-3+i<<"\t"<< funcbkg->GetParError(i)<<endl;
1060     }
1061
1062
1063   }
1064
1065   //sidebands integral - second approx (from fit)
1066   fSideBands = kFALSE;
1067   Double_t bkgInt;
1068   //cout<<"Compare intbkg1 = "<<intbkg1<<" and integral = ";
1069   if(ftypeOfFit4Sgn == 1) bkgInt=funcbkg1->Integral(fminMass,fmaxMass);
1070   else bkgInt=funcbkg->Integral(fminMass,fmaxMass);
1071   //cout<</*"------BkgInt(Fit) = "<<*/bkgInt<<endl;
1072
1073   //Signal integral - first approx
1074   Double_t sgnInt;
1075   sgnInt = totInt-bkgInt;
1076   //cout<<"------TotInt = "<<totInt<<"\tsgnInt = "<<sgnInt<<endl;
1077   if (sgnInt <= 0){
1078     cout<<"Setting sgnInt = - sgnInt"<<endl;
1079     sgnInt=(-1)*sgnInt;
1080   }
1081   /*Fit All Mass distribution with exponential + gaussian (+gaussian braodened) */
1082   TF1 *funcmass = new TF1(massname.Data(),this,&AliHFMassFitter::FitFunction4MassDistr,fminMass,fmaxMass,fNFinalPars,"AliHFMassFitter","FitFunction4MassDistr");
1083   cout<<"Function name = "<<funcmass->GetName()<<endl<<endl;
1084   funcmass->SetLineColor(4); //blue
1085
1086   //Set parameters
1087   cout<<"\nTOTAL FIT"<<endl;
1088
1089   if(fNFinalPars==5){
1090     funcmass->SetParNames("TotInt","Slope","SgnInt","Mean","Sigma");
1091     funcmass->SetParameters(totInt,slope1,sgnInt,fMass,fSigmaSgn);
1092
1093     //cout<<"Parameters set to: "<<totInt<<"\t"<<slope1<<"\t"<<sgnInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl;
1094     //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<fNFinalPars<<"\tgsidebands = "<<fSideBands<<endl;
1095     if(fFixPar[0]){
1096       funcmass->FixParameter(0,totInt);
1097     }
1098     if(fFixPar[1]){
1099       funcmass->FixParameter(1,slope1);
1100     }
1101     if(fFixPar[2]){
1102       funcmass->FixParameter(2,sgnInt);
1103     }
1104     if(fFixPar[3]){
1105       funcmass->FixParameter(3,fMass);
1106     }
1107     if(fFixPar[4]){
1108       funcmass->FixParameter(4,fSigmaSgn);
1109     }
1110   }
1111   if (fNFinalPars==6){
1112     funcmass->SetParNames("TotInt","Coef1","Coef2","SgnInt","Mean","Sigma");
1113     funcmass->SetParameters(totInt,slope1,conc1,sgnInt,fMass,fSigmaSgn);
1114
1115     //cout<<"Parameters set to: "<<totInt<<"\t"<<slope1<<"\t"<<conc1<<"\t"<<sgnInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl;
1116     //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<fNFinalPars<<"\tgsidebands = "<<fSideBands<<endl;
1117     if(fFixPar[0])funcmass->FixParameter(0,totInt);
1118     if(fFixPar[1])funcmass->FixParameter(1,slope1);
1119     if(fFixPar[2])funcmass->FixParameter(2,conc1);
1120     if(fFixPar[3])funcmass->FixParameter(3,sgnInt);
1121     if(fFixPar[4])funcmass->FixParameter(4,fMass);
1122     if(fFixPar[5])funcmass->FixParameter(5,fSigmaSgn);
1123     //
1124     //funcmass->FixParameter(2,sgnInt);
1125   }
1126   if(fNFinalPars==4){
1127     funcmass->SetParNames("Const","SgnInt","Mean","Sigma");
1128     if(ftypeOfFit4Sgn == 1) funcmass->SetParameters(0.,0.5*totInt,fMass,fSigmaSgn);
1129     else funcmass->SetParameters(0.,totInt,fMass,fSigmaSgn);
1130     if(fFixPar[0]) funcmass->FixParameter(0,0.);
1131     if(fFixPar[1])funcmass->FixParameter(1,sgnInt);
1132     if(fFixPar[2])funcmass->FixParameter(2,fMass);
1133     if(fFixPar[3])funcmass->FixParameter(3,fSigmaSgn);
1134    //cout<<"Parameters set to: "<<0.5*totInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl;
1135     //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<fNFinalPars<<"\tgsidebands = "<<fSideBands<<endl;
1136
1137   }
1138
1139   Int_t status;
1140
1141   status = fhistoInvMass->Fit(massname.Data(),"R,L,E,+,0");
1142   if (status != 0){
1143     cout<<"Minuit returned "<<status<<endl;
1144     return kFALSE;
1145   }
1146
1147   cout<<"fit done"<<endl;
1148   //reset value of fMass and fSigmaSgn to those found from fit
1149   fMass=funcmass->GetParameter(fNFinalPars-2);
1150   fSigmaSgn=funcmass->GetParameter(fNFinalPars-1);
1151
1152   for(Int_t i=0;i<fNFinalPars;i++){
1153     fFitPars[i+2*bkgPar-3]=funcmass->GetParameter(i);
1154     fFitPars[fNFinalPars+4*bkgPar-6+i]= funcmass->GetParError(i);
1155     //cout<<i+2*bkgPar-3<<"\t"<<funcmass->GetParameter(i)<<"\t\t"<<fNFinalPars+4*bkgPar-6+i<<"\t"<<funcmass->GetParError(i)<<endl;
1156   }
1157   /*
1158   //check: cout parameters
1159   for(Int_t i=0;i<2*(fNFinalPars+2*bkgPar-3);i++){
1160     cout<<i<<"\t"<<fFitPars[i]<<endl;
1161     }
1162   */
1163
1164   if(funcmass->GetParameter(fNFinalPars-1) <0 || funcmass->GetParameter(fNFinalPars-2) <0 || funcmass->GetParameter(fNFinalPars-3) <0 ) {
1165     cout<<"IntS or mean or sigma negative. You may tray to SetInitialGaussianSigma(..) and SetInitialGaussianMean(..)"<<endl;
1166     return kFALSE;
1167   }
1168
1169   //increase counter of number of fits done
1170   fcounter++;
1171
1172   //contour plots
1173   if(draw){
1174
1175     for (Int_t kpar=1; kpar<fNFinalPars;kpar++){
1176
1177       for(Int_t jpar=kpar+1;jpar<fNFinalPars;jpar++){
1178         cout<<"Par "<<kpar<<" and "<<jpar<<endl;
1179
1180         // produce 2 contours per couple of parameters
1181         TGraph* cont[2] = {0x0, 0x0};
1182         const Double_t errDef[2] = {1., 4.};
1183         for (Int_t i=0; i<2; i++) {
1184           gMinuit->SetErrorDef(errDef[i]);
1185           cont[i] = (TGraph*)gMinuit->Contour(80,kpar,jpar);
1186           cout<<"Minuit Status = "<<gMinuit->GetStatus()<<endl;
1187         }
1188
1189         if(!cont[0] || !cont[1]){
1190           cout<<"Skipping par "<<kpar<<" vs par "<<jpar<<endl;
1191           continue;
1192         }
1193
1194         // set graph titles and add them to the list
1195         TString title = "Contour plot";
1196         TString titleX = funcmass->GetParName(kpar);
1197         TString titleY = funcmass->GetParName(jpar);
1198         for (Int_t i=0; i<2; i++) {
1199           cont[i]->SetName( Form("cperr%d_%d%d", i, kpar, jpar) );
1200           cont[i]->SetTitle(title);
1201           cont[i]->GetXaxis()->SetTitle(titleX);
1202           cont[i]->GetYaxis()->SetTitle(titleY);
1203           cont[i]->GetYaxis()->SetLabelSize(0.033);
1204           cont[i]->GetYaxis()->SetTitleSize(0.033);
1205           cont[i]->GetYaxis()->SetTitleOffset(1.67);
1206
1207           fContourGraph->Add(cont[i]);
1208         }
1209
1210         // plot them
1211         TString cvname = Form("c%d%d", kpar, jpar);
1212         TCanvas *c4=new TCanvas(cvname,cvname,600,600);
1213         c4->cd();
1214         cont[1]->SetFillColor(38);
1215         cont[1]->Draw("alf");
1216         cont[0]->SetFillColor(9);
1217         cont[0]->Draw("lf");
1218
1219       }
1220
1221     }
1222
1223   }
1224
1225   if (ftypeOfFit4Sgn == 1 && funcbkg1) {
1226     delete funcbkg1;
1227     funcbkg1=NULL;
1228   }
1229   if (funcbkg) {
1230     delete funcbkg;
1231     funcbkg=NULL;
1232   }
1233   if (funcmass) {
1234     delete funcmass;
1235     funcmass=NULL;
1236   }
1237
1238   AddFunctionsToHisto();
1239   if (draw) DrawFit();
1240
1241
1242   return kTRUE;
1243 }
1244
1245 //______________________________________________________________________________
1246
1247 Bool_t AliHFMassFitter::RefitWithBkgOnly(Bool_t draw){
1248
1249   //perform a fit with background function only. Can be useful to try when fit fails to understand if it is because there's no signal
1250   //If you want to change the backgroud function or range use SetType or SetRangeFit before
1251
1252   TString bkgname="funcbkgonly";
1253   fSideBands = kFALSE;
1254
1255   TF1* funcbkg = new TF1(bkgname.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,fNFinalPars-3,"AliHFMassFitter","FitFunction4Bkg");
1256
1257   funcbkg->SetLineColor(kBlue+3); //dark blue
1258
1259   Double_t integral=fhistoInvMass->Integral(fhistoInvMass->FindBin(fminMass),fhistoInvMass->FindBin(fmaxMass),"width");
1260
1261   switch (ftypeOfFit4Bkg) {
1262   case 0: //gaus+expo
1263     funcbkg->SetParNames("BkgInt","Slope");
1264     funcbkg->SetParameters(integral,-2.);
1265     break;
1266   case 1:
1267     funcbkg->SetParNames("BkgInt","Slope");
1268     funcbkg->SetParameters(integral,-100.);
1269     break;
1270   case 2:
1271     funcbkg->SetParNames("BkgInt","Coef1","Coef2");
1272     funcbkg->SetParameters(integral,-10.,5);
1273     break;
1274   case 3:
1275     cout<<"Warning! This choice does not have a lot of sense..."<<endl;
1276     if(ftypeOfFit4Sgn==0){
1277       funcbkg->SetParNames("Const");
1278       funcbkg->SetParameter(0,0.);
1279       funcbkg->FixParameter(0,0.);
1280     }
1281     break;
1282   default:
1283     cout<<"Wrong choise of ftypeOfFit4Bkg ("<<ftypeOfFit4Bkg<<")"<<endl;
1284     return kFALSE;
1285     break;
1286   }
1287
1288
1289   Int_t status=fhistoInvMass->Fit(bkgname.Data(),"R,L,E,+,0");
1290   if (status != 0){
1291     cout<<"Minuit returned "<<status<<endl;
1292     return kFALSE;
1293   }
1294   AddFunctionsToHisto();
1295
1296   if(draw) DrawFit();
1297
1298   return kTRUE;
1299
1300 }
1301 //_________________________________________________________________________
1302 Double_t AliHFMassFitter::GetChiSquare() const{
1303   TF1 *funcmass=(TF1*)fhistoInvMass->GetFunction("funcmass");
1304   if(!funcmass) {
1305     cout<<"funcmass not found"<<endl;
1306     return -1;
1307   }
1308   return funcmass->GetChisquare();
1309 }
1310
1311 //_________________________________________________________________________
1312 Double_t AliHFMassFitter::GetReducedChiSquare() const{
1313   TF1 *funcmass=(TF1*)fhistoInvMass->GetFunction("funcmass");
1314   if(!funcmass) {
1315     cout<<"funcmass not found"<<endl;
1316     return -1;
1317   }
1318   return funcmass->GetChisquare()/funcmass->GetNDF();
1319 }
1320
1321 //*********output
1322
1323 //_________________________________________________________________________
1324 void  AliHFMassFitter::GetFitPars(Float_t *vector) const {
1325   // Return fit parameters
1326
1327   for(Int_t i=0;i<fParsSize;i++){
1328     vector[i]=fFitPars[i];
1329   }
1330 }
1331
1332
1333 //_________________________________________________________________________
1334 void AliHFMassFitter::IntS(Float_t *valuewitherror) const {
1335
1336   //gives the integral of signal obtained from fit parameters
1337   if(!valuewitherror)valuewitherror=new Float_t[2];
1338
1339   Int_t index=fParsSize/2 - 3;
1340   valuewitherror[0]=fFitPars[index];
1341   index=fParsSize - 3;
1342   valuewitherror[1]=fFitPars[index];
1343   }
1344
1345
1346 //_________________________________________________________________________
1347 void AliHFMassFitter::AddFunctionsToHisto(){
1348
1349   //Add the background function in the complete range to the list of functions attached to the histogram
1350
1351   //cout<<"AddFunctionsToHisto called"<<endl;
1352   TString bkgname = "funcbkg";
1353
1354   Bool_t done1=kFALSE,done2=kFALSE;
1355
1356   TString bkgnamesave=bkgname;
1357   TString testname=bkgname;
1358   testname += "FullRange";
1359   TF1 *testfunc=(TF1*)fhistoInvMass->FindObject(testname.Data());
1360   if(testfunc){
1361     done1=kTRUE;
1362     testfunc=0x0;
1363   }
1364   testname="funcbkgonly";
1365   testfunc=(TF1*)fhistoInvMass->FindObject(testname.Data());
1366   if(testfunc){
1367     done2=kTRUE;
1368     testfunc=0x0;
1369   }
1370
1371   if(done1 && done2){
1372     cout<<"AddFunctionsToHisto already used: exiting...."<<endl;
1373     return;
1374   }
1375
1376   TList *hlist=fhistoInvMass->GetListOfFunctions();
1377   hlist->ls();
1378
1379   if(!done2){
1380     TF1 *bonly=(TF1*)hlist->FindObject(testname.Data());
1381     if(!bonly){
1382       cout<<testname.Data()<<" not found looking for complete fit"<<endl;
1383     }else{
1384       bonly->SetLineColor(kBlue+3);
1385       hlist->Add((TF1*)bonly->Clone());
1386       if(bonly) {
1387         delete bonly;
1388         bonly=NULL;
1389       }
1390     }
1391
1392   }
1393
1394   if(!done1){
1395     TF1 *b=(TF1*)hlist->FindObject(bkgname.Data());
1396     if(!b){
1397       cout<<bkgname<<" not found, cannot produce "<<bkgname<<"FullRange and "<<bkgname<<"Recalc"<<endl;
1398       return;
1399     }
1400
1401     bkgname += "FullRange";
1402     TF1 *bfullrange=new TF1(bkgname.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,fNFinalPars-3,"AliHFMassFitter","FitFunction4Bkg");
1403     //cout<<bfullrange->GetName()<<endl;
1404     for(Int_t i=0;i<fNFinalPars-3;i++){
1405       bfullrange->SetParName(i,b->GetParName(i));
1406       bfullrange->SetParameter(i,b->GetParameter(i));
1407       bfullrange->SetParError(i,b->GetParError(i));
1408     }
1409     bfullrange->SetLineStyle(4);
1410     bfullrange->SetLineColor(14);
1411
1412     bkgnamesave += "Recalc";
1413
1414     TF1 *blastpar=new TF1(bkgnamesave.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,fNFinalPars-3,"AliHFMassFitter","FitFunction4Bkg");
1415
1416     TF1 *mass=fhistoInvMass->GetFunction("funcmass");
1417
1418     if (!mass){
1419       cout<<"funcmass doesn't exist "<<endl;
1420       return;
1421     }
1422
1423     //intBkg=intTot-intS
1424     blastpar->SetParameter(0,mass->GetParameter(0)-mass->GetParameter(fNFinalPars-3));
1425     blastpar->SetParError(0,mass->GetParError(fNFinalPars-3));
1426     if (fNFinalPars>=5) {
1427       blastpar->SetParameter(1,mass->GetParameter(1));
1428       blastpar->SetParError(1,mass->GetParError(1));
1429     }
1430     if (fNFinalPars==6) {
1431       blastpar->SetParameter(2,mass->GetParameter(2));
1432       blastpar->SetParError(2,mass->GetParError(2));
1433     }
1434
1435     blastpar->SetLineStyle(1);
1436     blastpar->SetLineColor(2);
1437
1438     hlist->Add((TF1*)bfullrange->Clone());
1439     hlist->Add((TF1*)blastpar->Clone());
1440     hlist->ls();
1441
1442     if(bfullrange) {
1443       delete bfullrange;
1444       bfullrange=NULL;
1445     }
1446     if(blastpar) {
1447       delete blastpar;
1448       blastpar=NULL;
1449     }
1450   }
1451
1452
1453 }
1454
1455 //_________________________________________________________________________
1456
1457 TH1F* AliHFMassFitter::GetHistoClone() const{
1458
1459   TH1F* hout=(TH1F*)fhistoInvMass->Clone(fhistoInvMass->GetName());
1460   return hout;
1461 }
1462 //_________________________________________________________________________
1463
1464 void AliHFMassFitter::WriteHisto(TString path) const {
1465
1466   //Write the histogram in the default file HFMassFitterOutput.root
1467
1468   if (fcounter == 0) {
1469     cout<<"Use MassFitter method before WriteHisto"<<endl;
1470     return;
1471   }
1472   TH1F* hget=(TH1F*)fhistoInvMass->Clone();
1473
1474   path += "HFMassFitterOutput.root";
1475   TFile *output;
1476
1477   if (fcounter == 1) output = new TFile(path.Data(),"recreate");
1478   else output = new TFile(path.Data(),"update");
1479   output->cd();
1480   hget->Write();
1481   fContourGraph->Write();
1482
1483
1484   output->Close();
1485
1486   cout<<fcounter<<" "<<hget->GetName()<<" written in "<<path<<endl;
1487
1488   if(output) {
1489     delete output;
1490     output=NULL;
1491   }
1492
1493 }
1494
1495 //_________________________________________________________________________
1496
1497 void AliHFMassFitter::WriteNtuple(TString path) const{
1498   //TNtuple* nget=(TNtuple*)fntuParam->Clone();
1499   path += "HFMassFitterOutput.root";
1500   TFile *output = new TFile(path.Data(),"update");
1501   output->cd();
1502   fntuParam->Write();
1503   //nget->Write();
1504   output->Close();
1505   //cout<<nget->GetName()<<" written in "<<path<<endl;
1506   cout<<fntuParam->GetName()<<" written in "<<path<<endl;
1507   /*
1508   if(nget) {
1509     //delete nget;
1510     nget=NULL;
1511   }
1512   */
1513   if(output) {
1514     delete output;
1515     output=NULL;
1516   }
1517 }
1518
1519 //_________________________________________________________________________
1520 void AliHFMassFitter::WriteCanvas(TString userIDstring,TString path,Double_t nsigma,Int_t writeFitInfo,Bool_t draw) const{
1521
1522  //write the canvas in a root file
1523
1524   gStyle->SetOptStat(0);
1525   gStyle->SetCanvasColor(0);
1526   gStyle->SetFrameFillColor(0);
1527
1528   TString type="";
1529
1530   switch (ftypeOfFit4Bkg){
1531   case 0:
1532     type="Exp"; //3+2
1533     break;
1534   case 1:
1535     type="Lin"; //3+2
1536     break;
1537   case 2:
1538     type="Pl2"; //3+3
1539     break;
1540   case 3:
1541     type="noB"; //3+1
1542     break;
1543   }
1544
1545   TString filename=Form("%sMassFit.root",type.Data());
1546   filename.Prepend(userIDstring);
1547   path.Append(filename);
1548
1549   TFile* outputcv=new TFile(path.Data(),"update");
1550
1551   TCanvas* c=(TCanvas*)GetPad(nsigma,writeFitInfo);
1552   c->SetName(Form("%s%s%s",c->GetName(),userIDstring.Data(),type.Data()));
1553   if(draw)c->DrawClone();
1554   outputcv->cd();
1555   c->Write();
1556   outputcv->Close();
1557 }
1558
1559 //_________________________________________________________________________
1560
1561 TVirtualPad* AliHFMassFitter::GetPad(Double_t nsigma,Int_t writeFitInfo)const{
1562   //return a TVirtualPad with the fitted histograms and info
1563
1564   TString cvtitle="fit of ";
1565   cvtitle+=fhistoInvMass->GetName();
1566   TString cvname="c";
1567   cvname+=fcounter;
1568
1569   TCanvas *c=new TCanvas(cvname,cvtitle);
1570   PlotFit(c->cd(),nsigma,writeFitInfo);
1571   return c->cd();
1572 }
1573 //_________________________________________________________________________
1574
1575 void AliHFMassFitter::PlotFit(TVirtualPad* pd,Double_t nsigma,Int_t writeFitInfo)const{
1576   //plot histogram, fit functions and write parameters according to verbosity level (0,1,>1)
1577   gStyle->SetOptStat(0);
1578   gStyle->SetCanvasColor(0);
1579   gStyle->SetFrameFillColor(0);
1580
1581   cout<<"nsigma = "<<nsigma<<endl;
1582   cout<<"Verbosity = "<<writeFitInfo<<endl;
1583
1584   TH1F* hdraw=GetHistoClone();
1585
1586   if(!hdraw->GetFunction("funcmass") && !hdraw->GetFunction("funcbkgFullRange") && !hdraw->GetFunction("funcbkgRecalc")&& !hdraw->GetFunction("funcbkgonly")){
1587     cout<<"Probably fit failed and you didn't try to refit with background only, there's no function to be drawn"<<endl;
1588     return;
1589   }
1590
1591   if(hdraw->GetFunction("funcbkgonly")){ //Warning! if this function is present, no chance to draw the other!
1592     cout<<"Drawing background fit only"<<endl;
1593     hdraw->SetMinimum(0);
1594     hdraw->GetXaxis()->SetRangeUser(fminMass,fmaxMass);
1595     pd->cd();
1596     hdraw->SetMarkerStyle(20);
1597     hdraw->DrawClone("PE");
1598     hdraw->GetFunction("funcbkgonly")->DrawClone("sames");
1599
1600     if(writeFitInfo > 0){
1601       TPaveText *pinfo=new TPaveText(0.6,0.86,1.,1.,"NDC");
1602       pinfo->SetBorderSize(0);
1603       pinfo->SetFillStyle(0);
1604       TF1* f=hdraw->GetFunction("funcbkgonly");
1605       for (Int_t i=0;i<fNFinalPars-3;i++){
1606         pinfo->SetTextColor(kBlue+3);
1607         TString str=Form("%s = %f #pm %f",f->GetParName(i),f->GetParameter(i),f->GetParError(i));
1608         pinfo->AddText(str);
1609       }
1610
1611       pinfo->AddText(Form("Reduced #chi^{2} = %.3f",f->GetChisquare()/f->GetNDF()));
1612       pd->cd();
1613       pinfo->DrawClone();
1614
1615
1616     }
1617
1618     return;
1619   }
1620
1621   hdraw->SetMinimum(0);
1622   hdraw->GetXaxis()->SetRangeUser(fminMass,fmaxMass);
1623   pd->cd();
1624   hdraw->SetMarkerStyle(20);
1625   hdraw->DrawClone("PE");
1626 //   if(hdraw->GetFunction("funcbkgFullRange")) hdraw->GetFunction("funcbkgFullRange")->DrawClone("same");
1627 //   if(hdraw->GetFunction("funcbkgRecalc")) hdraw->GetFunction("funcbkgRecalc")->DrawClone("same");
1628   if(hdraw->GetFunction("funcmass")) hdraw->GetFunction("funcmass")->DrawClone("same");
1629
1630   if(writeFitInfo > 0){
1631     TPaveText *pinfob=new TPaveText(0.6,0.86,1.,1.,"NDC");
1632     TPaveText *pinfom=new TPaveText(0.6,0.7,1.,.87,"NDC");
1633     pinfob->SetBorderSize(0);
1634     pinfob->SetFillStyle(0);
1635     pinfom->SetBorderSize(0);
1636     pinfom->SetFillStyle(0);
1637     TF1* ff=fhistoInvMass->GetFunction("funcmass");
1638
1639     for (Int_t i=fNFinalPars-3;i<fNFinalPars;i++){
1640       pinfom->SetTextColor(kBlue);
1641       TString str=Form("%s = %f #pm %f",ff->GetParName(i),ff->GetParameter(i),ff->GetParError(i));
1642       if(!(writeFitInfo==1 && i==fNFinalPars-3)) pinfom->AddText(str);
1643     }
1644     pd->cd();
1645     pinfom->DrawClone();
1646
1647     TPaveText *pinfo2=new TPaveText(0.1,0.1,0.6,0.4,"NDC");
1648     pinfo2->SetBorderSize(0);
1649     pinfo2->SetFillStyle(0);
1650
1651     Double_t signif, signal, bkg, errsignif, errsignal, errbkg;
1652
1653     Significance(nsigma,signif,errsignif);
1654     Signal(nsigma,signal,errsignal);
1655     Background(nsigma,bkg, errbkg);
1656     /*
1657       Significance(1.828,1.892,signif,errsignif);
1658       Signal(1.828,1.892,signal,errsignal);
1659       Background(1.828,1.892,bkg, errbkg);
1660     */
1661     TString str=Form("Significance (%.0f#sigma) %.1f #pm %.1f ",nsigma,signif,errsignif);
1662     pinfo2->AddText(str);
1663     str=Form("S (%.0f#sigma) %.0f #pm %.0f ",nsigma,signal,errsignal);
1664     pinfo2->AddText(str);
1665     str=Form("B (%.0f#sigma) %.0f #pm %.0f",nsigma,bkg,errbkg);
1666     pinfo2->AddText(str);
1667
1668     pd->cd();
1669     pinfo2->Draw();
1670
1671     if(writeFitInfo > 1){
1672       for (Int_t i=0;i<fNFinalPars-3;i++){
1673         pinfob->SetTextColor(kRed);
1674         str=Form("%s = %f #pm %f",ff->GetParName(i),ff->GetParameter(i),ff->GetParError(i));
1675         pinfob->AddText(str);
1676       }
1677       pd->cd();
1678       pinfob->DrawClone();
1679     }
1680
1681
1682   }
1683   return;
1684 }
1685
1686 //_________________________________________________________________________
1687
1688 void AliHFMassFitter::DrawHere(TVirtualPad* pd,Double_t nsigma,Int_t writeFitInfo) const {
1689   //draws histogram together with fit functions with default nice colors in user canvas
1690   PlotFit(pd,nsigma,writeFitInfo);
1691
1692   pd->Draw();
1693
1694 }
1695 //_________________________________________________________________________
1696 void AliHFMassFitter::DrawFit(Double_t nsigma) const{
1697
1698   //draws histogram together with fit functions with default nice colors
1699   gStyle->SetOptStat(0);
1700   gStyle->SetCanvasColor(0);
1701   gStyle->SetFrameFillColor(0);
1702
1703
1704   TCanvas* c=(TCanvas*)GetPad(nsigma,1);
1705   c->Draw();
1706
1707 }
1708
1709
1710 //_________________________________________________________________________
1711
1712 void AliHFMassFitter::PrintParTitles() const{
1713
1714   //prints to screen the parameters names
1715
1716   TF1 *f=fhistoInvMass->GetFunction("funcmass");
1717   if(!f) {
1718     cout<<"Fit function not found"<<endl;
1719     return;
1720   }
1721
1722   cout<<"Parameter Titles \n";
1723   for(Int_t i=0;i<fNFinalPars;i++){
1724     cout<<"Par "<<i<<": "<<f->GetParName(i)<<endl;
1725   }
1726   cout<<endl;
1727
1728 }
1729
1730 //************ significance
1731
1732 //_________________________________________________________________________
1733
1734 void AliHFMassFitter::Signal(Double_t nOfSigma,Double_t &signal,Double_t &errsignal) const {
1735   // Return signal integral in mean+- n sigma
1736
1737   if(fcounter==0) {
1738     cout<<"Use MassFitter method before Signal"<<endl;
1739     return;
1740   }
1741
1742   Double_t min=fMass-nOfSigma*fSigmaSgn;
1743   Double_t max=fMass+nOfSigma*fSigmaSgn;
1744
1745   Signal(min,max,signal,errsignal);
1746
1747
1748   return;
1749
1750 }
1751
1752 //_________________________________________________________________________
1753
1754 void AliHFMassFitter::Signal(Double_t min, Double_t max, Double_t &signal,Double_t &errsignal) const {
1755
1756   // Return signal integral in a range
1757
1758   if(fcounter==0) {
1759     cout<<"Use MassFitter method before Signal"<<endl;
1760     return;
1761   }
1762
1763   //functions names
1764   TString bkgname="funcbkgRecalc";
1765   TString bkg1name="funcbkg1Recalc";
1766   TString massname="funcmass";
1767
1768
1769   TF1 *funcbkg=0;
1770   TF1 *funcmass=fhistoInvMass->GetFunction(massname.Data());
1771   if(!funcmass){
1772     cout<<"AliHFMassFitter::Signal() ERROR -> Mass distr function not found!"<<endl;
1773     return;
1774   }
1775
1776   if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
1777   else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
1778
1779   if(!funcbkg){
1780     cout<<"AliHFMassFitter::Signal() ERROR -> Bkg function not found!"<<endl;
1781     return;
1782   }
1783
1784   Int_t np=fNFinalPars-3;
1785
1786   Double_t intS,intSerr;
1787
1788  //relative error evaluation
1789   intS=funcmass->GetParameter(np);
1790   intSerr=funcmass->GetParError(np);
1791
1792   cout<<"Sgn relative error evaluation from fit: "<<intSerr/intS<<endl;
1793   Double_t background,errbackground;
1794   Background(min,max,background,errbackground);
1795
1796   //signal +/- error in the range
1797
1798   Double_t mass=funcmass->Integral(min, max)/fhistoInvMass->GetBinWidth(4);
1799   signal=mass - background;
1800   errsignal=(intSerr/intS)*signal;/*assume relative error is the same as for total integral*/
1801
1802 }
1803
1804 //_________________________________________________________________________
1805
1806 void AliHFMassFitter::Background(Double_t nOfSigma,Double_t &background,Double_t &errbackground) const {
1807   // Return background integral in mean+- n sigma
1808
1809   if(fcounter==0) {
1810     cout<<"Use MassFitter method before Background"<<endl;
1811     return;
1812   }
1813   Double_t min=fMass-nOfSigma*fSigmaSgn;
1814   Double_t max=fMass+nOfSigma*fSigmaSgn;
1815
1816   Background(min,max,background,errbackground);
1817
1818   return;
1819
1820 }
1821 //___________________________________________________________________________
1822
1823 void AliHFMassFitter::Background(Double_t min, Double_t max, Double_t &background,Double_t &errbackground) const {
1824   // Return background integral in a range
1825
1826   if(fcounter==0) {
1827     cout<<"Use MassFitter method before Background"<<endl;
1828     return;
1829   }
1830
1831   //functions names
1832   TString bkgname="funcbkgRecalc";
1833   TString bkg1name="funcbkg1Recalc";
1834
1835   TF1 *funcbkg=0;
1836   if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
1837   else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
1838   if(!funcbkg){
1839     cout<<"AliHFMassFitter::Background() ERROR -> Bkg function not found!"<<endl;
1840     return;
1841   }
1842
1843
1844   Double_t intB,intBerr;
1845
1846   //relative error evaluation: from final parameters of the fit
1847   if(ftypeOfFit4Bkg==3 && ftypeOfFit4Sgn == 0) cout<<"No background fit: Bkg relative error evaluation put to zero"<<endl;
1848   else{
1849     intB=funcbkg->GetParameter(0);
1850     intBerr=funcbkg->GetParError(0);
1851     cout<<"Bkg relative error evaluation: from final parameters of the fit: "<<intBerr/intB<<endl;
1852   }
1853
1854   //relative error evaluation: from histo
1855
1856   intB=fhistoInvMass->Integral(1,fSideBandl)+fhistoInvMass->Integral(fSideBandr,fNbin);
1857   Double_t sum2=0;
1858   for(Int_t i=1;i<=fSideBandl;i++){
1859     sum2+=fhistoInvMass->GetBinError(i)*fhistoInvMass->GetBinError(i);
1860   }
1861   for(Int_t i=fSideBandr;i<=fNbin;i++){
1862     sum2+=fhistoInvMass->GetBinError(i)*fhistoInvMass->GetBinError(i);
1863   }
1864
1865   intBerr=TMath::Sqrt(sum2);
1866   cout<<"Bkg relative error evaluation: from histo: "<<intBerr/intB<<endl;
1867
1868   cout<<"Last estimation of bkg error is used"<<endl;
1869
1870   //backround +/- error in the range
1871   if (ftypeOfFit4Bkg == 3 && ftypeOfFit4Sgn == 0) {
1872     background = 0;
1873     errbackground = 0;
1874   }
1875   else{
1876     background=funcbkg->Integral(min,max)/(Double_t)fhistoInvMass->GetBinWidth(2);
1877     errbackground=intBerr/intB*background; // assume relative error is the same as for total integral
1878     //cout<<"integral = "<<funcbkg->Integral(min, max)<<"\tbinW = "<<fhistoInvMass->GetBinWidth(2)<<endl;
1879   }
1880   return;
1881
1882 }
1883
1884
1885 //__________________________________________________________________________
1886
1887 void AliHFMassFitter::Significance(Double_t nOfSigma,Double_t &significance,Double_t &errsignificance) const  {
1888   // Return significance in mean+- n sigma
1889
1890   Double_t min=fMass-nOfSigma*fSigmaSgn;
1891   Double_t max=fMass+nOfSigma*fSigmaSgn;
1892   Significance(min, max, significance, errsignificance);
1893
1894   return;
1895 }
1896
1897 //__________________________________________________________________________
1898
1899 void AliHFMassFitter::Significance(Double_t min, Double_t max, Double_t &significance,Double_t &errsignificance) const {
1900   // Return significance integral in a range
1901
1902   Double_t signal,errsignal,background,errbackground;
1903   Signal(min, max,signal,errsignal);
1904   Background(min, max,background,errbackground);
1905
1906   if (signal+background <= 0.){
1907     cout<<"Cannot calculate significance because of div by 0!"<<endl;
1908     significance=-1;
1909     errsignificance=0;
1910     return;
1911   }
1912
1913   significance =  signal/TMath::Sqrt(signal+background);
1914
1915   errsignificance = TMath::Sqrt(significance*significance/(signal+background)/(signal+background)*(1/4.*errsignal*errsignal+errbackground*errbackground)+significance*significance/signal/signal*errsignal*errsignal);
1916
1917   return;
1918 }
1919
1920