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
1338   Int_t index=fParsSize/2 - 3;
1339   valuewitherror[0]=fFitPars[index];
1340   index=fParsSize - 3;
1341   valuewitherror[1]=fFitPars[index];
1342   }
1343
1344
1345 //_________________________________________________________________________
1346 void AliHFMassFitter::AddFunctionsToHisto(){
1347
1348   //Add the background function in the complete range to the list of functions attached to the histogram
1349
1350   //cout<<"AddFunctionsToHisto called"<<endl;
1351   TString bkgname = "funcbkg";
1352
1353   Bool_t done1=kFALSE,done2=kFALSE;
1354
1355   TString bkgnamesave=bkgname;
1356   TString testname=bkgname;
1357   testname += "FullRange";
1358   TF1 *testfunc=(TF1*)fhistoInvMass->FindObject(testname.Data());
1359   if(testfunc){
1360     done1=kTRUE;
1361     testfunc=0x0;
1362   }
1363   testname="funcbkgonly";
1364   testfunc=(TF1*)fhistoInvMass->FindObject(testname.Data());
1365   if(testfunc){
1366     done2=kTRUE;
1367     testfunc=0x0;
1368   }
1369
1370   if(done1 && done2){
1371     cout<<"AddFunctionsToHisto already used: exiting...."<<endl;
1372     return;
1373   }
1374
1375   TList *hlist=fhistoInvMass->GetListOfFunctions();
1376   hlist->ls();
1377
1378   if(!done2){
1379     TF1 *bonly=(TF1*)hlist->FindObject(testname.Data());
1380     if(!bonly){
1381       cout<<testname.Data()<<" not found looking for complete fit"<<endl;
1382     }else{
1383       bonly->SetLineColor(kBlue+3);
1384       hlist->Add((TF1*)bonly->Clone());
1385       if(bonly) {
1386         delete bonly;
1387         bonly=NULL;
1388       }
1389     }
1390
1391   }
1392
1393   if(!done1){
1394     TF1 *b=(TF1*)hlist->FindObject(bkgname.Data());
1395     if(!b){
1396       cout<<bkgname<<" not found, cannot produce "<<bkgname<<"FullRange and "<<bkgname<<"Recalc"<<endl;
1397       return;
1398     }
1399
1400     bkgname += "FullRange";
1401     TF1 *bfullrange=new TF1(bkgname.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,fNFinalPars-3,"AliHFMassFitter","FitFunction4Bkg");
1402     //cout<<bfullrange->GetName()<<endl;
1403     for(Int_t i=0;i<fNFinalPars-3;i++){
1404       bfullrange->SetParName(i,b->GetParName(i));
1405       bfullrange->SetParameter(i,b->GetParameter(i));
1406       bfullrange->SetParError(i,b->GetParError(i));
1407     }
1408     bfullrange->SetLineStyle(4);
1409     bfullrange->SetLineColor(14);
1410
1411     bkgnamesave += "Recalc";
1412
1413     TF1 *blastpar=new TF1(bkgnamesave.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,fNFinalPars-3,"AliHFMassFitter","FitFunction4Bkg");
1414
1415     TF1 *mass=fhistoInvMass->GetFunction("funcmass");
1416
1417     if (!mass){
1418       cout<<"funcmass doesn't exist "<<endl;
1419       return;
1420     }
1421
1422     //intBkg=intTot-intS
1423     blastpar->SetParameter(0,mass->GetParameter(0)-mass->GetParameter(fNFinalPars-3));
1424     blastpar->SetParError(0,mass->GetParError(fNFinalPars-3));
1425     if (fNFinalPars>=5) {
1426       blastpar->SetParameter(1,mass->GetParameter(1));
1427       blastpar->SetParError(1,mass->GetParError(1));
1428     }
1429     if (fNFinalPars==6) {
1430       blastpar->SetParameter(2,mass->GetParameter(2));
1431       blastpar->SetParError(2,mass->GetParError(2));
1432     }
1433
1434     blastpar->SetLineStyle(1);
1435     blastpar->SetLineColor(2);
1436
1437     hlist->Add((TF1*)bfullrange->Clone());
1438     hlist->Add((TF1*)blastpar->Clone());
1439     hlist->ls();
1440
1441     if(bfullrange) {
1442       delete bfullrange;
1443       bfullrange=NULL;
1444     }
1445     if(blastpar) {
1446       delete blastpar;
1447       blastpar=NULL;
1448     }
1449   }
1450
1451
1452 }
1453
1454 //_________________________________________________________________________
1455
1456 TH1F* AliHFMassFitter::GetHistoClone() const{
1457
1458   TH1F* hout=(TH1F*)fhistoInvMass->Clone(fhistoInvMass->GetName());
1459   return hout;
1460 }
1461 //_________________________________________________________________________
1462
1463 void AliHFMassFitter::WriteHisto(TString path) const {
1464
1465   //Write the histogram in the default file HFMassFitterOutput.root
1466
1467   if (fcounter == 0) {
1468     cout<<"Use MassFitter method before WriteHisto"<<endl;
1469     return;
1470   }
1471   TH1F* hget=(TH1F*)fhistoInvMass->Clone();
1472
1473   path += "HFMassFitterOutput.root";
1474   TFile *output;
1475
1476   if (fcounter == 1) output = new TFile(path.Data(),"recreate");
1477   else output = new TFile(path.Data(),"update");
1478   output->cd();
1479   hget->Write();
1480   fContourGraph->Write();
1481
1482
1483   output->Close();
1484
1485   cout<<fcounter<<" "<<hget->GetName()<<" written in "<<path<<endl;
1486
1487   if(output) {
1488     delete output;
1489     output=NULL;
1490   }
1491
1492 }
1493
1494 //_________________________________________________________________________
1495
1496 void AliHFMassFitter::WriteNtuple(TString path) const{
1497   //TNtuple* nget=(TNtuple*)fntuParam->Clone();
1498   path += "HFMassFitterOutput.root";
1499   TFile *output = new TFile(path.Data(),"update");
1500   output->cd();
1501   fntuParam->Write();
1502   //nget->Write();
1503   output->Close();
1504   //cout<<nget->GetName()<<" written in "<<path<<endl;
1505   cout<<fntuParam->GetName()<<" written in "<<path<<endl;
1506   /*
1507   if(nget) {
1508     //delete nget;
1509     nget=NULL;
1510   }
1511   */
1512   if(output) {
1513     delete output;
1514     output=NULL;
1515   }
1516 }
1517
1518 //_________________________________________________________________________
1519 void AliHFMassFitter::WriteCanvas(TString userIDstring,TString path,Double_t nsigma,Int_t writeFitInfo,Bool_t draw) const{
1520
1521  //write the canvas in a root file
1522
1523   gStyle->SetOptStat(0);
1524   gStyle->SetCanvasColor(0);
1525   gStyle->SetFrameFillColor(0);
1526
1527   TString type="";
1528
1529   switch (ftypeOfFit4Bkg){
1530   case 0:
1531     type="Exp"; //3+2
1532     break;
1533   case 1:
1534     type="Lin"; //3+2
1535     break;
1536   case 2:
1537     type="Pl2"; //3+3
1538     break;
1539   case 3:
1540     type="noB"; //3+1
1541     break;
1542   }
1543
1544   TString filename=Form("%sMassFit.root",type.Data());
1545   filename.Prepend(userIDstring);
1546   path.Append(filename);
1547
1548   TFile* outputcv=new TFile(path.Data(),"update");
1549
1550   TCanvas* c=(TCanvas*)GetPad(nsigma,writeFitInfo);
1551   c->SetName(Form("%s%s%s",c->GetName(),userIDstring.Data(),type.Data()));
1552   if(draw)c->DrawClone();
1553   outputcv->cd();
1554   c->Write();
1555   outputcv->Close();
1556 }
1557
1558 //_________________________________________________________________________
1559
1560 TVirtualPad* AliHFMassFitter::GetPad(Double_t nsigma,Int_t writeFitInfo)const{
1561   //return a TVirtualPad with the fitted histograms and info
1562
1563   TString cvtitle="fit of ";
1564   cvtitle+=fhistoInvMass->GetName();
1565   TString cvname="c";
1566   cvname+=fcounter;
1567
1568   TCanvas *c=new TCanvas(cvname,cvtitle);
1569   PlotFit(c->cd(),nsigma,writeFitInfo);
1570   return c->cd();
1571 }
1572 //_________________________________________________________________________
1573
1574 void AliHFMassFitter::PlotFit(TVirtualPad* pd,Double_t nsigma,Int_t writeFitInfo)const{
1575   //plot histogram, fit functions and write parameters according to verbosity level (0,1,>1)
1576   gStyle->SetOptStat(0);
1577   gStyle->SetCanvasColor(0);
1578   gStyle->SetFrameFillColor(0);
1579
1580   cout<<"nsigma = "<<nsigma<<endl;
1581   cout<<"Verbosity = "<<writeFitInfo<<endl;
1582
1583   TH1F* hdraw=GetHistoClone();
1584
1585   if(!hdraw->GetFunction("funcmass") && !hdraw->GetFunction("funcbkgFullRange") && !hdraw->GetFunction("funcbkgRecalc")&& !hdraw->GetFunction("funcbkgonly")){
1586     cout<<"Probably fit failed and you didn't try to refit with background only, there's no function to be drawn"<<endl;
1587     return;
1588   }
1589
1590   if(hdraw->GetFunction("funcbkgonly")){ //Warning! if this function is present, no chance to draw the other!
1591     cout<<"Drawing background fit only"<<endl;
1592     hdraw->SetMinimum(0);
1593     hdraw->GetXaxis()->SetRangeUser(fminMass,fmaxMass);
1594     pd->cd();
1595     hdraw->SetMarkerStyle(20);
1596     hdraw->DrawClone("PE");
1597     hdraw->GetFunction("funcbkgonly")->DrawClone("sames");
1598
1599     if(writeFitInfo > 0){
1600       TPaveText *pinfo=new TPaveText(0.6,0.86,1.,1.,"NDC");
1601       pinfo->SetBorderSize(0);
1602       pinfo->SetFillStyle(0);
1603       TF1* f=hdraw->GetFunction("funcbkgonly");
1604       for (Int_t i=0;i<fNFinalPars-3;i++){
1605         pinfo->SetTextColor(kBlue+3);
1606         TString str=Form("%s = %f #pm %f",f->GetParName(i),f->GetParameter(i),f->GetParError(i));
1607         pinfo->AddText(str);
1608       }
1609
1610       pinfo->AddText(Form("Reduced #chi^{2} = %.3f",f->GetChisquare()/f->GetNDF()));
1611       pd->cd();
1612       pinfo->DrawClone();
1613
1614
1615     }
1616
1617     return;
1618   }
1619
1620   hdraw->SetMinimum(0);
1621   hdraw->GetXaxis()->SetRangeUser(fminMass,fmaxMass);
1622   pd->cd();
1623   hdraw->SetMarkerStyle(20);
1624   hdraw->DrawClone("PE");
1625 //   if(hdraw->GetFunction("funcbkgFullRange")) hdraw->GetFunction("funcbkgFullRange")->DrawClone("same");
1626 //   if(hdraw->GetFunction("funcbkgRecalc")) hdraw->GetFunction("funcbkgRecalc")->DrawClone("same");
1627   if(hdraw->GetFunction("funcmass")) hdraw->GetFunction("funcmass")->DrawClone("same");
1628
1629   if(writeFitInfo > 0){
1630     TPaveText *pinfob=new TPaveText(0.6,0.86,1.,1.,"NDC");
1631     TPaveText *pinfom=new TPaveText(0.6,0.7,1.,.87,"NDC");
1632     pinfob->SetBorderSize(0);
1633     pinfob->SetFillStyle(0);
1634     pinfom->SetBorderSize(0);
1635     pinfom->SetFillStyle(0);
1636     TF1* ff=fhistoInvMass->GetFunction("funcmass");
1637
1638     for (Int_t i=fNFinalPars-3;i<fNFinalPars;i++){
1639       pinfom->SetTextColor(kBlue);
1640       TString str=Form("%s = %f #pm %f",ff->GetParName(i),ff->GetParameter(i),ff->GetParError(i));
1641       if(!(writeFitInfo==1 && i==fNFinalPars-3)) pinfom->AddText(str);
1642     }
1643     pd->cd();
1644     pinfom->DrawClone();
1645
1646     TPaveText *pinfo2=new TPaveText(0.1,0.1,0.6,0.4,"NDC");
1647     pinfo2->SetBorderSize(0);
1648     pinfo2->SetFillStyle(0);
1649
1650     Double_t signif, signal, bkg, errsignif, errsignal, errbkg;
1651
1652     Significance(nsigma,signif,errsignif);
1653     Signal(nsigma,signal,errsignal);
1654     Background(nsigma,bkg, errbkg);
1655     /*
1656       Significance(1.828,1.892,signif,errsignif);
1657       Signal(1.828,1.892,signal,errsignal);
1658       Background(1.828,1.892,bkg, errbkg);
1659     */
1660     TString str=Form("Significance (%.0f#sigma) %.1f #pm %.1f ",nsigma,signif,errsignif);
1661     pinfo2->AddText(str);
1662     str=Form("S (%.0f#sigma) %.0f #pm %.0f ",nsigma,signal,errsignal);
1663     pinfo2->AddText(str);
1664     str=Form("B (%.0f#sigma) %.0f #pm %.0f",nsigma,bkg,errbkg);
1665     pinfo2->AddText(str);
1666
1667     pd->cd();
1668     pinfo2->Draw();
1669
1670     if(writeFitInfo > 1){
1671       for (Int_t i=0;i<fNFinalPars-3;i++){
1672         pinfob->SetTextColor(kRed);
1673         str=Form("%s = %f #pm %f",ff->GetParName(i),ff->GetParameter(i),ff->GetParError(i));
1674         pinfob->AddText(str);
1675       }
1676       pd->cd();
1677       pinfob->DrawClone();
1678     }
1679
1680
1681   }
1682   return;
1683 }
1684
1685 //_________________________________________________________________________
1686
1687 void AliHFMassFitter::DrawHere(TVirtualPad* pd,Double_t nsigma,Int_t writeFitInfo) const {
1688   //draws histogram together with fit functions with default nice colors in user canvas
1689   PlotFit(pd,nsigma,writeFitInfo);
1690
1691   pd->Draw();
1692
1693 }
1694 //_________________________________________________________________________
1695 void AliHFMassFitter::DrawFit(Double_t nsigma) const{
1696
1697   //draws histogram together with fit functions with default nice colors
1698   gStyle->SetOptStat(0);
1699   gStyle->SetCanvasColor(0);
1700   gStyle->SetFrameFillColor(0);
1701
1702
1703   TCanvas* c=(TCanvas*)GetPad(nsigma,1);
1704   c->Draw();
1705
1706 }
1707
1708
1709 //_________________________________________________________________________
1710
1711 void AliHFMassFitter::PrintParTitles() const{
1712
1713   //prints to screen the parameters names
1714
1715   TF1 *f=fhistoInvMass->GetFunction("funcmass");
1716   if(!f) {
1717     cout<<"Fit function not found"<<endl;
1718     return;
1719   }
1720
1721   cout<<"Parameter Titles \n";
1722   for(Int_t i=0;i<fNFinalPars;i++){
1723     cout<<"Par "<<i<<": "<<f->GetParName(i)<<endl;
1724   }
1725   cout<<endl;
1726
1727 }
1728
1729 //************ significance
1730
1731 //_________________________________________________________________________
1732
1733 void AliHFMassFitter::Signal(Double_t nOfSigma,Double_t &signal,Double_t &errsignal) const {
1734   // Return signal integral in mean+- n sigma
1735
1736   if(fcounter==0) {
1737     cout<<"Use MassFitter method before Signal"<<endl;
1738     return;
1739   }
1740
1741   Double_t min=fMass-nOfSigma*fSigmaSgn;
1742   Double_t max=fMass+nOfSigma*fSigmaSgn;
1743
1744   Signal(min,max,signal,errsignal);
1745
1746
1747   return;
1748
1749 }
1750
1751 //_________________________________________________________________________
1752
1753 void AliHFMassFitter::Signal(Double_t min, Double_t max, Double_t &signal,Double_t &errsignal) const {
1754
1755   // Return signal integral in a range
1756
1757   if(fcounter==0) {
1758     cout<<"Use MassFitter method before Signal"<<endl;
1759     return;
1760   }
1761
1762   //functions names
1763   TString bkgname="funcbkgRecalc";
1764   TString bkg1name="funcbkg1Recalc";
1765   TString massname="funcmass";
1766
1767
1768   TF1 *funcbkg=0;
1769   TF1 *funcmass=fhistoInvMass->GetFunction(massname.Data());
1770   if(!funcmass){
1771     cout<<"AliHFMassFitter::Signal() ERROR -> Mass distr function not found!"<<endl;
1772     return;
1773   }
1774
1775   if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
1776   else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
1777
1778   if(!funcbkg){
1779     cout<<"AliHFMassFitter::Signal() ERROR -> Bkg function not found!"<<endl;
1780     return;
1781   }
1782
1783   Int_t np=fNFinalPars-3;
1784
1785   Double_t intS,intSerr;
1786
1787  //relative error evaluation
1788   intS=funcmass->GetParameter(np);
1789   intSerr=funcmass->GetParError(np);
1790
1791   cout<<"Sgn relative error evaluation from fit: "<<intSerr/intS<<endl;
1792   Double_t background,errbackground;
1793   Background(min,max,background,errbackground);
1794
1795   //signal +/- error in the range
1796
1797   Double_t mass=funcmass->Integral(min, max)/fhistoInvMass->GetBinWidth(4);
1798   signal=mass - background;
1799   errsignal=(intSerr/intS)*signal;/*assume relative error is the same as for total integral*/
1800
1801 }
1802
1803 //_________________________________________________________________________
1804
1805 void AliHFMassFitter::Background(Double_t nOfSigma,Double_t &background,Double_t &errbackground) const {
1806   // Return background integral in mean+- n sigma
1807
1808   if(fcounter==0) {
1809     cout<<"Use MassFitter method before Background"<<endl;
1810     return;
1811   }
1812   Double_t min=fMass-nOfSigma*fSigmaSgn;
1813   Double_t max=fMass+nOfSigma*fSigmaSgn;
1814
1815   Background(min,max,background,errbackground);
1816
1817   return;
1818
1819 }
1820 //___________________________________________________________________________
1821
1822 void AliHFMassFitter::Background(Double_t min, Double_t max, Double_t &background,Double_t &errbackground) const {
1823   // Return background integral in a range
1824
1825   if(fcounter==0) {
1826     cout<<"Use MassFitter method before Background"<<endl;
1827     return;
1828   }
1829
1830   //functions names
1831   TString bkgname="funcbkgRecalc";
1832   TString bkg1name="funcbkg1Recalc";
1833
1834   TF1 *funcbkg=0;
1835   if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
1836   else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
1837   if(!funcbkg){
1838     cout<<"AliHFMassFitter::Background() ERROR -> Bkg function not found!"<<endl;
1839     return;
1840   }
1841
1842
1843   Double_t intB,intBerr;
1844
1845   //relative error evaluation: from final parameters of the fit
1846   if(ftypeOfFit4Bkg==3 && ftypeOfFit4Sgn == 0) cout<<"No background fit: Bkg relative error evaluation put to zero"<<endl;
1847   else{
1848     intB=funcbkg->GetParameter(0);
1849     intBerr=funcbkg->GetParError(0);
1850     cout<<"Bkg relative error evaluation: from final parameters of the fit: "<<intBerr/intB<<endl;
1851   }
1852
1853   //relative error evaluation: from histo
1854
1855   intB=fhistoInvMass->Integral(1,fSideBandl)+fhistoInvMass->Integral(fSideBandr,fNbin);
1856   Double_t sum2=0;
1857   for(Int_t i=1;i<=fSideBandl;i++){
1858     sum2+=fhistoInvMass->GetBinError(i)*fhistoInvMass->GetBinError(i);
1859   }
1860   for(Int_t i=fSideBandr;i<=fNbin;i++){
1861     sum2+=fhistoInvMass->GetBinError(i)*fhistoInvMass->GetBinError(i);
1862   }
1863
1864   intBerr=TMath::Sqrt(sum2);
1865   cout<<"Bkg relative error evaluation: from histo: "<<intBerr/intB<<endl;
1866
1867   cout<<"Last estimation of bkg error is used"<<endl;
1868
1869   //backround +/- error in the range
1870   if (ftypeOfFit4Bkg == 3 && ftypeOfFit4Sgn == 0) {
1871     background = 0;
1872     errbackground = 0;
1873   }
1874   else{
1875     background=funcbkg->Integral(min,max)/(Double_t)fhistoInvMass->GetBinWidth(2);
1876     errbackground=intBerr/intB*background; // assume relative error is the same as for total integral
1877     //cout<<"integral = "<<funcbkg->Integral(min, max)<<"\tbinW = "<<fhistoInvMass->GetBinWidth(2)<<endl;
1878   }
1879   return;
1880
1881 }
1882
1883
1884 //__________________________________________________________________________
1885
1886 void AliHFMassFitter::Significance(Double_t nOfSigma,Double_t &significance,Double_t &errsignificance) const  {
1887   // Return significance in mean+- n sigma
1888
1889   Double_t min=fMass-nOfSigma*fSigmaSgn;
1890   Double_t max=fMass+nOfSigma*fSigmaSgn;
1891   Significance(min, max, significance, errsignificance);
1892
1893   return;
1894 }
1895
1896 //__________________________________________________________________________
1897
1898 void AliHFMassFitter::Significance(Double_t min, Double_t max, Double_t &significance,Double_t &errsignificance) const {
1899   // Return significance integral in a range
1900
1901   Double_t signal,errsignal,background,errbackground;
1902   Signal(min, max,signal,errsignal);
1903   Background(min, max,background,errbackground);
1904
1905   if (signal+background <= 0.){
1906     cout<<"Cannot calculate significance because of div by 0!"<<endl;
1907     significance=-1;
1908     errsignificance=0;
1909     return;
1910   }
1911
1912   significance =  signal/TMath::Sqrt(signal+background);
1913
1914   errsignificance = TMath::Sqrt(significance*significance/(signal+background)/(signal+background)*(1/4.*errsignal*errsignal+errbackground*errbackground)+significance*significance/signal/signal*errsignal*errsignal);
1915
1916   return;
1917 }
1918
1919