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=- 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       cout<<"fix1"<<endl;
1097       funcmass->FixParameter(0,totInt);
1098     }
1099     if(fFixPar[1]){
1100       cout<<"fix2"<<endl;
1101       funcmass->FixParameter(1,slope1);
1102     }
1103     if(fFixPar[2]){
1104       cout<<"fix3"<<endl;
1105       funcmass->FixParameter(2,sgnInt);
1106     }
1107     if(fFixPar[3]){
1108       cout<<"fix4"<<endl;
1109       funcmass->FixParameter(3,fMass);
1110     }
1111     if(fFixPar[4]){
1112       cout<<"fix5"<<endl;
1113       funcmass->FixParameter(4,fSigmaSgn);
1114     }
1115   }
1116   if (fNFinalPars==6){
1117     funcmass->SetParNames("TotInt","Coef1","Coef2","SgnInt","Mean","Sigma");
1118     funcmass->SetParameters(totInt,slope1,conc1,sgnInt,fMass,fSigmaSgn);
1119
1120     //cout<<"Parameters set to: "<<totInt<<"\t"<<slope1<<"\t"<<conc1<<"\t"<<sgnInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl;
1121     //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<fNFinalPars<<"\tgsidebands = "<<fSideBands<<endl;
1122     if(fFixPar[0])funcmass->FixParameter(0,totInt);
1123     if(fFixPar[1])funcmass->FixParameter(1,slope1);
1124     if(fFixPar[2])funcmass->FixParameter(2,conc1);
1125     if(fFixPar[3])funcmass->FixParameter(3,sgnInt);
1126     if(fFixPar[4])funcmass->FixParameter(4,fMass);
1127     if(fFixPar[5])funcmass->FixParameter(5,fSigmaSgn);
1128     //
1129     //funcmass->FixParameter(2,sgnInt);
1130   }
1131   if(fNFinalPars==4){
1132     funcmass->SetParNames("Const","SgnInt","Mean","Sigma");
1133     if(ftypeOfFit4Sgn == 1) funcmass->SetParameters(0.,0.5*totInt,fMass,fSigmaSgn);
1134     else funcmass->SetParameters(0.,totInt,fMass,fSigmaSgn);
1135     if(fFixPar[0]) funcmass->FixParameter(0,0.);
1136     //cout<<"Parameters set to: "<<0.5*totInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl;
1137     //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<fNFinalPars<<"\tgsidebands = "<<fSideBands<<endl;
1138
1139   }
1140
1141   Int_t status;
1142
1143   status = fhistoInvMass->Fit(massname.Data(),"R,L,E,+,0");
1144   if (status != 0){
1145     cout<<"Minuit returned "<<status<<endl;
1146     return kFALSE;
1147   }
1148
1149   cout<<"fit done"<<endl;
1150   //reset value of fMass and fSigmaSgn to those found from fit
1151   fMass=funcmass->GetParameter(fNFinalPars-2);
1152   fSigmaSgn=funcmass->GetParameter(fNFinalPars-1);
1153
1154   for(Int_t i=0;i<fNFinalPars;i++){
1155     fFitPars[i+2*bkgPar-3]=funcmass->GetParameter(i);
1156     fFitPars[fNFinalPars+4*bkgPar-6+i]= funcmass->GetParError(i);
1157     //cout<<i+2*bkgPar-3<<"\t"<<funcmass->GetParameter(i)<<"\t\t"<<fNFinalPars+4*bkgPar-6+i<<"\t"<<funcmass->GetParError(i)<<endl;
1158   }
1159   /*
1160   //check: cout parameters
1161   for(Int_t i=0;i<2*(fNFinalPars+2*bkgPar-3);i++){
1162     cout<<i<<"\t"<<fFitPars[i]<<endl;
1163     }
1164   */
1165
1166   if(funcmass->GetParameter(fNFinalPars-1) <0 || funcmass->GetParameter(fNFinalPars-2) <0 || funcmass->GetParameter(fNFinalPars-3) <0 ) {
1167     cout<<"IntS or mean or sigma negative. You may tray to SetInitialGaussianSigma(..) and SetInitialGaussianMean(..)"<<endl;
1168     return kFALSE;
1169   }
1170
1171   //increase counter of number of fits done
1172   fcounter++;
1173
1174   //contour plots
1175   if(draw){
1176
1177     for (Int_t kpar=1; kpar<fNFinalPars;kpar++){
1178
1179       for(Int_t jpar=kpar+1;jpar<fNFinalPars;jpar++){
1180         cout<<"Par "<<kpar<<" and "<<jpar<<endl;
1181
1182         // produce 2 contours per couple of parameters
1183         TGraph* cont[2] = {0x0, 0x0};
1184         const Double_t errDef[2] = {1., 4.};
1185         for (Int_t i=0; i<2; i++) {
1186           gMinuit->SetErrorDef(errDef[i]);
1187           cont[i] = (TGraph*)gMinuit->Contour(80,kpar,jpar);
1188           cout<<"Minuit Status = "<<gMinuit->GetStatus()<<endl;
1189         }
1190
1191         if(!cont[0] || !cont[1]){
1192           cout<<"Skipping par "<<kpar<<" vs par "<<jpar<<endl;
1193           continue;
1194         }
1195
1196         // set graph titles and add them to the list
1197         TString title = "Contour plot";
1198         TString titleX = funcmass->GetParName(kpar);
1199         TString titleY = funcmass->GetParName(jpar);
1200         for (Int_t i=0; i<2; i++) {
1201           cont[i]->SetName( Form("cperr%d_%d%d", i, kpar, jpar) );
1202           cont[i]->SetTitle(title);
1203           cont[i]->GetXaxis()->SetTitle(titleX);
1204           cont[i]->GetYaxis()->SetTitle(titleY);
1205           cont[i]->GetYaxis()->SetLabelSize(0.033);
1206           cont[i]->GetYaxis()->SetTitleSize(0.033);
1207           cont[i]->GetYaxis()->SetTitleOffset(1.67);
1208
1209           fContourGraph->Add(cont[i]);
1210         }
1211
1212         // plot them
1213         TString cvname = Form("c%d%d", kpar, jpar);
1214         TCanvas *c4=new TCanvas(cvname,cvname,600,600);
1215         c4->cd();
1216         cont[1]->SetFillColor(38);
1217         cont[1]->Draw("alf");
1218         cont[0]->SetFillColor(9);
1219         cont[0]->Draw("lf");
1220
1221       }
1222
1223     }
1224
1225   }
1226
1227   if (ftypeOfFit4Sgn == 1 && funcbkg1) {
1228     delete funcbkg1;
1229     funcbkg1=NULL;
1230   }
1231   if (funcbkg) {
1232     delete funcbkg;
1233     funcbkg=NULL;
1234   }
1235   if (funcmass) {
1236     delete funcmass;
1237     funcmass=NULL;
1238   }
1239
1240   AddFunctionsToHisto();
1241   if (draw) DrawFit();
1242
1243
1244   return kTRUE;
1245 }
1246
1247 //______________________________________________________________________________
1248
1249 Bool_t AliHFMassFitter::RefitWithBkgOnly(Bool_t draw){
1250
1251   //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
1252   //If you want to change the backgroud function or range use SetType or SetRangeFit before
1253
1254   TString bkgname="funcbkgonly";
1255   fSideBands = kFALSE;
1256
1257   TF1* funcbkg = new TF1(bkgname.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,fNFinalPars-3,"AliHFMassFitter","FitFunction4Bkg");
1258
1259   funcbkg->SetLineColor(kBlue+3); //dark blue
1260
1261   Double_t integral=fhistoInvMass->Integral(fhistoInvMass->FindBin(fminMass),fhistoInvMass->FindBin(fmaxMass),"width");
1262
1263   switch (ftypeOfFit4Bkg) {
1264   case 0: //gaus+expo
1265     funcbkg->SetParNames("BkgInt","Slope");
1266     funcbkg->SetParameters(integral,-2.);
1267     break;
1268   case 1:
1269     funcbkg->SetParNames("BkgInt","Slope");
1270     funcbkg->SetParameters(integral,-100.);
1271     break;
1272   case 2:
1273     funcbkg->SetParNames("BkgInt","Coef1","Coef2");
1274     funcbkg->SetParameters(integral,-10.,5);
1275     break;
1276   case 3:
1277     cout<<"Warning! This choice does not have a lot of sense..."<<endl;
1278     if(ftypeOfFit4Sgn==0){
1279       funcbkg->SetParNames("Const");
1280       funcbkg->SetParameter(0,0.);
1281       funcbkg->FixParameter(0,0.);
1282     }
1283     break;
1284   default:
1285     cout<<"Wrong choise of ftypeOfFit4Bkg ("<<ftypeOfFit4Bkg<<")"<<endl;
1286     return kFALSE;
1287     break;
1288   }
1289
1290
1291   Int_t status=fhistoInvMass->Fit(bkgname.Data(),"R,L,E,+,0");
1292   if (status != 0){
1293     cout<<"Minuit returned "<<status<<endl;
1294     return kFALSE;
1295   }
1296   AddFunctionsToHisto();
1297
1298   if(draw) DrawFit();
1299
1300   return kTRUE;
1301
1302 }
1303 //_________________________________________________________________________
1304 Double_t AliHFMassFitter::GetChiSquare() const{
1305   TF1 *funcmass=(TF1*)fhistoInvMass->GetFunction("funcmass");
1306   if(!funcmass) {
1307     cout<<"funcmass not found"<<endl;
1308     return -1;
1309   }
1310   return funcmass->GetChisquare();
1311 }
1312
1313 //_________________________________________________________________________
1314 Double_t AliHFMassFitter::GetReducedChiSquare() const{
1315   TF1 *funcmass=(TF1*)fhistoInvMass->GetFunction("funcmass");
1316   if(!funcmass) {
1317     cout<<"funcmass not found"<<endl;
1318     return -1;
1319   }
1320   return funcmass->GetChisquare()/funcmass->GetNDF();
1321 }
1322
1323 //*********output
1324
1325 //_________________________________________________________________________
1326 void  AliHFMassFitter::GetFitPars(Float_t *vector) const {
1327   // Return fit parameters
1328
1329   for(Int_t i=0;i<fParsSize;i++){
1330     vector[i]=fFitPars[i];
1331   }
1332 }
1333
1334
1335 //_________________________________________________________________________
1336 void AliHFMassFitter::IntS(Float_t *valuewitherror) const {
1337
1338   //gives the integral of signal obtained from fit parameters
1339   if(!valuewitherror)valuewitherror=new Float_t[2];
1340
1341   Int_t index=fParsSize/2 - 3;
1342   valuewitherror[0]=fFitPars[index];
1343   index=fParsSize - 3;
1344   valuewitherror[1]=fFitPars[index];
1345   }
1346
1347
1348 //_________________________________________________________________________
1349 void AliHFMassFitter::AddFunctionsToHisto(){
1350
1351   //Add the background function in the complete range to the list of functions attached to the histogram
1352
1353   cout<<"AddFunctionsToHisto called"<<endl;
1354   TString bkgname = "funcbkg";
1355
1356   Bool_t done1=kFALSE,done2=kFALSE;
1357
1358   TString bkgnamesave=bkgname;
1359   TString testname=bkgname;
1360   testname += "FullRange";
1361   TF1 *testfunc=(TF1*)fhistoInvMass->FindObject(testname.Data());
1362   if(testfunc){
1363     done1=kTRUE;
1364     testfunc=0x0;
1365   }
1366   testname="funcbkgonly";
1367   testfunc=(TF1*)fhistoInvMass->FindObject(testname.Data());
1368   if(testfunc){
1369     done2=kTRUE;
1370     testfunc=0x0;
1371   }
1372
1373   if(done1 && done2){
1374     cout<<"AddFunctionsToHisto already used: exiting...."<<endl;
1375     return;
1376   }
1377
1378   TList *hlist=fhistoInvMass->GetListOfFunctions();
1379   hlist->ls();
1380
1381   if(!done2){
1382     TF1 *bonly=(TF1*)hlist->FindObject(testname.Data());
1383     if(!bonly){
1384       cout<<testname.Data()<<" not found looking for complete fit"<<endl;
1385     }else{
1386       bonly->SetLineColor(kBlue+3);
1387       hlist->Add((TF1*)bonly->Clone());
1388       if(bonly) {
1389         delete bonly;
1390         bonly=NULL;
1391       }
1392     }
1393
1394   }
1395
1396   if(!done1){
1397     TF1 *b=(TF1*)hlist->FindObject(bkgname.Data());
1398     if(!b){
1399       cout<<bkgname<<" not found, cannot produce "<<bkgname<<"FullRange and "<<bkgname<<"Recalc"<<endl;
1400       return;
1401     }
1402
1403     bkgname += "FullRange";
1404     TF1 *bfullrange=new TF1(bkgname.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,fNFinalPars,"AliHFMassFitter","FitFunction4Bkg");
1405     //cout<<bfullrange->GetName()<<endl;
1406     for(Int_t i=0;i<fNFinalPars;i++){
1407       //cout<<i<<" di "<<fNFinalPars<<endl;
1408       bfullrange->SetParName(i,b->GetParName(i));
1409       bfullrange->SetParameter(i,b->GetParameter(i));
1410       bfullrange->SetParError(i,b->GetParError(i));
1411     }
1412     bfullrange->SetLineStyle(4);
1413     bfullrange->SetLineColor(14);
1414
1415     bkgnamesave += "Recalc";
1416
1417     TF1 *blastpar=new TF1(bkgnamesave.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,fNFinalPars,"AliHFMassFitter","FitFunction4Bkg");
1418
1419     TF1 *mass=fhistoInvMass->GetFunction("funcmass");
1420
1421     if (!mass){
1422       cout<<"funcmass doesn't exist "<<endl;
1423       return;
1424     }
1425
1426     //intBkg=intTot-intS
1427     blastpar->SetParameter(0,mass->GetParameter(0)-mass->GetParameter(fNFinalPars-3));
1428     blastpar->SetParError(0,mass->GetParError(fNFinalPars-3));
1429     if (fNFinalPars>=2) {
1430       blastpar->SetParameter(1,mass->GetParameter(1));
1431       blastpar->SetParError(1,mass->GetParError(1));
1432     }
1433     if (fNFinalPars==3) {
1434       blastpar->SetParameter(2,mass->GetParameter(2));
1435       blastpar->SetParError(2,mass->GetParError(2));
1436     }
1437
1438     blastpar->SetLineStyle(1);
1439     blastpar->SetLineColor(2);
1440
1441     hlist->Add((TF1*)bfullrange->Clone());
1442     hlist->Add((TF1*)blastpar->Clone());
1443     hlist->ls();
1444
1445     if(bfullrange) {
1446       delete bfullrange;
1447       bfullrange=NULL;
1448     }
1449     if(blastpar) {
1450       delete blastpar;
1451       blastpar=NULL;
1452     }
1453   }
1454
1455
1456 }
1457
1458 //_________________________________________________________________________
1459
1460 TH1F* AliHFMassFitter::GetHistoClone() const{
1461
1462   TH1F* hout=(TH1F*)fhistoInvMass->Clone(fhistoInvMass->GetName());
1463   return hout;
1464 }
1465 //_________________________________________________________________________
1466
1467 void AliHFMassFitter::WriteHisto(TString path) const {
1468
1469   //Write the histogram in the default file HFMassFitterOutput.root
1470
1471   if (fcounter == 0) {
1472     cout<<"Use MassFitter method before WriteHisto"<<endl;
1473     return;
1474   }
1475   TH1F* hget=(TH1F*)fhistoInvMass->Clone();
1476
1477   path += "HFMassFitterOutput.root";
1478   TFile *output;
1479
1480   if (fcounter == 1) output = new TFile(path.Data(),"recreate");
1481   else output = new TFile(path.Data(),"update");
1482   output->cd();
1483   hget->Write();
1484   fContourGraph->Write();
1485
1486
1487   output->Close();
1488
1489   cout<<fcounter<<" "<<hget->GetName()<<" written in "<<path<<endl;
1490
1491   if(output) {
1492     delete output;
1493     output=NULL;
1494   }
1495
1496 }
1497
1498 //_________________________________________________________________________
1499
1500 void AliHFMassFitter::WriteNtuple(TString path) const{
1501   //TNtuple* nget=(TNtuple*)fntuParam->Clone();
1502   path += "HFMassFitterOutput.root";
1503   TFile *output = new TFile(path.Data(),"update");
1504   output->cd();
1505   fntuParam->Write();
1506   //nget->Write();
1507   output->Close();
1508   //cout<<nget->GetName()<<" written in "<<path<<endl;
1509   cout<<fntuParam->GetName()<<" written in "<<path<<endl;
1510   /*
1511   if(nget) {
1512     //delete nget;
1513     nget=NULL;
1514   }
1515   */
1516   if(output) {
1517     delete output;
1518     output=NULL;
1519   }
1520 }
1521
1522 //_________________________________________________________________________
1523 void AliHFMassFitter::WriteCanvas(TString userIDstring,TString path,Double_t nsigma,Int_t writeFitInfo,Bool_t draw) const{
1524
1525  //write the canvas in a root file
1526
1527   gStyle->SetOptStat(0);
1528   gStyle->SetCanvasColor(0);
1529   gStyle->SetFrameFillColor(0);
1530
1531   TString type="";
1532
1533   switch (ftypeOfFit4Bkg){
1534   case 0:
1535     type="Exp"; //3+2
1536     break;
1537   case 1:
1538     type="Lin"; //3+2
1539     break;
1540   case 2:
1541     type="Pl2"; //3+3
1542     break;
1543   case 3:
1544     type="noB"; //3+1
1545     break;
1546   }
1547
1548   TString filename=Form("%sMassFit.root",type.Data());
1549   filename.Prepend(userIDstring);
1550   path.Append(filename);
1551
1552   TFile* outputcv=new TFile(path.Data(),"update");
1553
1554   TCanvas* c=(TCanvas*)GetPad(nsigma,writeFitInfo);
1555   c->SetName(Form("%s%s%s",c->GetName(),userIDstring.Data(),type.Data()));
1556   if(draw)c->DrawClone();
1557   outputcv->cd();
1558   c->Write();
1559   outputcv->Close();
1560 }
1561
1562 //_________________________________________________________________________
1563
1564 TVirtualPad* AliHFMassFitter::GetPad(Double_t nsigma,Int_t writeFitInfo)const{
1565   //return a TVirtualPad with the fitted histograms and info
1566
1567   TString cvtitle="fit of ";
1568   cvtitle+=fhistoInvMass->GetName();
1569   TString cvname="c";
1570   cvname+=fcounter;
1571
1572   TCanvas *c=new TCanvas(cvname,cvtitle);
1573   PlotFit(c->cd(),nsigma,writeFitInfo);
1574   return c->cd();
1575 }
1576 //_________________________________________________________________________
1577
1578 void AliHFMassFitter::PlotFit(TVirtualPad* pd,Double_t nsigma,Int_t writeFitInfo)const{
1579   //plot histogram, fit functions and write parameters according to verbosity level (0,1,>1)
1580   gStyle->SetOptStat(0);
1581   gStyle->SetCanvasColor(0);
1582   gStyle->SetFrameFillColor(0);
1583
1584   cout<<"nsigma = "<<nsigma<<endl;
1585   cout<<"Verbosity = "<<writeFitInfo<<endl;
1586
1587   TH1F* hdraw=GetHistoClone();
1588
1589   if(!hdraw->GetFunction("funcmass") && !hdraw->GetFunction("funcbkgFullRange") && !hdraw->GetFunction("funcbkgRecalc")&& !hdraw->GetFunction("funcbkgonly")){
1590     cout<<"Probably fit failed and you didn't try to refit with background only, there's no function to be drawn"<<endl;
1591     return;
1592   }
1593
1594   if(hdraw->GetFunction("funcbkgonly")){ //Warning! if this function is present, no chance to draw the other!
1595     cout<<"Drawing background fit only"<<endl;
1596     hdraw->SetMinimum(0);
1597     hdraw->GetXaxis()->SetRangeUser(fminMass,fmaxMass);
1598     pd->cd();
1599     hdraw->SetMarkerStyle(20);
1600     hdraw->DrawClone("PE");
1601     hdraw->GetFunction("funcbkgonly")->DrawClone("sames");
1602
1603     if(writeFitInfo > 0){
1604       TPaveText *pinfo=new TPaveText(0.6,0.86,1.,1.,"NDC");
1605       pinfo->SetBorderSize(0);
1606       pinfo->SetFillStyle(0);
1607       TF1* f=hdraw->GetFunction("funcbkgonly");
1608       for (Int_t i=0;i<fNFinalPars-3;i++){
1609         pinfo->SetTextColor(kBlue+3);
1610         TString str=Form("%s = %f #pm %f",f->GetParName(i),f->GetParameter(i),f->GetParError(i));
1611         pinfo->AddText(str);
1612       }
1613
1614       pinfo->AddText(Form("Reduced #chi^{2} = %.3f",f->GetChisquare()/f->GetNDF()));
1615       pd->cd();
1616       pinfo->DrawClone();
1617
1618
1619     }
1620
1621     return;
1622   }
1623
1624   hdraw->SetMinimum(0);
1625   hdraw->GetXaxis()->SetRangeUser(fminMass,fmaxMass);
1626   pd->cd();
1627   hdraw->SetMarkerStyle(20);
1628   hdraw->DrawClone("PE");
1629   if(hdraw->GetFunction("funcbkgFullRange")) hdraw->GetFunction("funcbkgFullRange")->DrawClone("same");
1630   if(hdraw->GetFunction("funcbkgRecalc")) hdraw->GetFunction("funcbkgRecalc")->DrawClone("same");
1631   if(hdraw->GetFunction("funcmass")) hdraw->GetFunction("funcmass")->DrawClone("same");
1632
1633   if(writeFitInfo > 0){
1634     TPaveText *pinfob=new TPaveText(0.6,0.86,1.,1.,"NDC");
1635     TPaveText *pinfom=new TPaveText(0.6,0.7,1.,.87,"NDC");
1636     pinfob->SetBorderSize(0);
1637     pinfob->SetFillStyle(0);
1638     pinfom->SetBorderSize(0);
1639     pinfom->SetFillStyle(0);
1640     TF1* ff=fhistoInvMass->GetFunction("funcmass");
1641
1642     for (Int_t i=fNFinalPars-3;i<fNFinalPars;i++){
1643       pinfom->SetTextColor(kBlue);
1644       TString str=Form("%s = %f #pm %f",ff->GetParName(i),ff->GetParameter(i),ff->GetParError(i));
1645       if(!(writeFitInfo==1 && i==fNFinalPars-3)) pinfom->AddText(str);
1646     }
1647     pd->cd();
1648     pinfom->DrawClone();
1649
1650     TPaveText *pinfo2=new TPaveText(0.1,0.1,0.6,0.4,"NDC");
1651     pinfo2->SetBorderSize(0);
1652     pinfo2->SetFillStyle(0);
1653
1654     Double_t signif, signal, bkg, errsignif, errsignal, errbkg;
1655
1656     Significance(nsigma,signif,errsignif);
1657     Signal(nsigma,signal,errsignal);
1658     Background(nsigma,bkg, errbkg);
1659     /*
1660       Significance(1.828,1.892,signif,errsignif);
1661       Signal(1.828,1.892,signal,errsignal);
1662       Background(1.828,1.892,bkg, errbkg);
1663     */
1664     TString str=Form("Significance (%.0f#sigma) %.1f #pm %.1f ",nsigma,signif,errsignif);
1665     pinfo2->AddText(str);
1666     str=Form("S (%.0f#sigma) %.0f #pm %.0f ",nsigma,signal,errsignal);
1667     pinfo2->AddText(str);
1668     str=Form("B (%.0f#sigma) %.0f #pm %.0f",nsigma,bkg,errbkg);
1669     pinfo2->AddText(str);
1670
1671     pd->cd();
1672     pinfo2->Draw();
1673
1674     if(writeFitInfo > 1){
1675       for (Int_t i=0;i<fNFinalPars-3;i++){
1676         pinfob->SetTextColor(kRed);
1677         str=Form("%s = %f #pm %f",ff->GetParName(i),ff->GetParameter(i),ff->GetParError(i));
1678         pinfob->AddText(str);
1679       }
1680       pd->cd();
1681       pinfob->DrawClone();
1682     }
1683
1684
1685   }
1686   return;
1687 }
1688
1689 //_________________________________________________________________________
1690
1691 void AliHFMassFitter::DrawHere(TVirtualPad* pd,Double_t nsigma,Int_t writeFitInfo) const {
1692   //draws histogram together with fit functions with default nice colors in user canvas
1693   PlotFit(pd,nsigma,writeFitInfo);
1694
1695   pd->Draw();
1696
1697 }
1698 //_________________________________________________________________________
1699 void AliHFMassFitter::DrawFit(Double_t nsigma) const{
1700
1701   //draws histogram together with fit functions with default nice colors
1702   gStyle->SetOptStat(0);
1703   gStyle->SetCanvasColor(0);
1704   gStyle->SetFrameFillColor(0);
1705
1706
1707   TCanvas* c=(TCanvas*)GetPad(nsigma,1);
1708   c->Draw();
1709
1710 }
1711
1712
1713 //_________________________________________________________________________
1714
1715 void AliHFMassFitter::PrintParTitles() const{
1716
1717   //prints to screen the parameters names
1718
1719   TF1 *f=fhistoInvMass->GetFunction("funcmass");
1720   if(!f) {
1721     cout<<"Fit function not found"<<endl;
1722     return;
1723   }
1724
1725   cout<<"Parameter Titles \n";
1726   for(Int_t i=0;i<fNFinalPars;i++){
1727     cout<<"Par "<<i<<": "<<f->GetParName(i)<<endl;
1728   }
1729   cout<<endl;
1730
1731 }
1732
1733 //************ significance
1734
1735 //_________________________________________________________________________
1736
1737 void AliHFMassFitter::Signal(Double_t nOfSigma,Double_t &signal,Double_t &errsignal) const {
1738   // Return signal integral in mean+- n sigma
1739
1740   if(fcounter==0) {
1741     cout<<"Use MassFitter method before Signal"<<endl;
1742     return;
1743   }
1744
1745   Double_t min=fMass-nOfSigma*fSigmaSgn;
1746   Double_t max=fMass+nOfSigma*fSigmaSgn;
1747
1748   Signal(min,max,signal,errsignal);
1749
1750
1751   return;
1752
1753 }
1754
1755 //_________________________________________________________________________
1756
1757 void AliHFMassFitter::Signal(Double_t min, Double_t max, Double_t &signal,Double_t &errsignal) const {
1758
1759   // Return signal integral in a range
1760
1761   if(fcounter==0) {
1762     cout<<"Use MassFitter method before Signal"<<endl;
1763     return;
1764   }
1765
1766   //functions names
1767   TString bkgname="funcbkgRecalc";
1768   TString bkg1name="funcbkg1Recalc";
1769   TString massname="funcmass";
1770
1771
1772   TF1 *funcbkg=0;
1773   TF1 *funcmass=fhistoInvMass->GetFunction(massname.Data());
1774   if(!funcmass){
1775     cout<<"AliHFMassFitter::Signal() ERROR -> Mass distr function not found!"<<endl;
1776     return;
1777   }
1778
1779   if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
1780   else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
1781
1782   if(!funcbkg){
1783     cout<<"AliHFMassFitter::Signal() ERROR -> Bkg function not found!"<<endl;
1784     return;
1785   }
1786
1787   Int_t np=fNFinalPars-3;
1788
1789   Double_t intS,intSerr;
1790
1791  //relative error evaluation
1792   intS=funcmass->GetParameter(np);
1793   intSerr=funcmass->GetParError(np);
1794
1795   cout<<"Sgn relative error evaluation from fit: "<<intSerr/intS<<endl;
1796   Double_t background,errbackground;
1797   Background(min,max,background,errbackground);
1798
1799   //signal +/- error in the range
1800
1801   Double_t mass=funcmass->Integral(min, max)/fhistoInvMass->GetBinWidth(4);
1802   signal=mass - background;
1803   errsignal=(intSerr/intS)*signal;/*assume relative error is the same as for total integral*/
1804
1805 }
1806
1807 //_________________________________________________________________________
1808
1809 void AliHFMassFitter::Background(Double_t nOfSigma,Double_t &background,Double_t &errbackground) const {
1810   // Return background integral in mean+- n sigma
1811
1812   if(fcounter==0) {
1813     cout<<"Use MassFitter method before Background"<<endl;
1814     return;
1815   }
1816   Double_t min=fMass-nOfSigma*fSigmaSgn;
1817   Double_t max=fMass+nOfSigma*fSigmaSgn;
1818
1819   Background(min,max,background,errbackground);
1820
1821   return;
1822
1823 }
1824 //___________________________________________________________________________
1825
1826 void AliHFMassFitter::Background(Double_t min, Double_t max, Double_t &background,Double_t &errbackground) const {
1827   // Return background integral in a range
1828
1829   if(fcounter==0) {
1830     cout<<"Use MassFitter method before Background"<<endl;
1831     return;
1832   }
1833
1834   //functions names
1835   TString bkgname="funcbkgRecalc";
1836   TString bkg1name="funcbkg1Recalc";
1837
1838   TF1 *funcbkg=0;
1839   if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
1840   else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
1841   if(!funcbkg){
1842     cout<<"AliHFMassFitter::Background() ERROR -> Bkg function not found!"<<endl;
1843     return;
1844   }
1845
1846
1847   Double_t intB,intBerr;
1848
1849   //relative error evaluation: from final parameters of the fit
1850   if(ftypeOfFit4Bkg==3 && ftypeOfFit4Sgn == 0) cout<<"No background fit: Bkg relative error evaluation put to zero"<<endl;
1851   else{
1852     intB=funcbkg->GetParameter(0);
1853     intBerr=funcbkg->GetParError(0);
1854     cout<<"Bkg relative error evaluation: from final parameters of the fit: "<<intBerr/intB<<endl;
1855   }
1856
1857   //relative error evaluation: from histo
1858
1859   intB=fhistoInvMass->Integral(1,fSideBandl)+fhistoInvMass->Integral(fSideBandr,fNbin);
1860   Double_t sum2=0;
1861   for(Int_t i=1;i<=fSideBandl;i++){
1862     sum2+=fhistoInvMass->GetBinError(i)*fhistoInvMass->GetBinError(i);
1863   }
1864   for(Int_t i=fSideBandr;i<=fNbin;i++){
1865     sum2+=fhistoInvMass->GetBinError(i)*fhistoInvMass->GetBinError(i);
1866   }
1867
1868   intBerr=TMath::Sqrt(sum2);
1869   cout<<"Bkg relative error evaluation: from histo: "<<intBerr/intB<<endl;
1870
1871   cout<<"Last estimation of bkg error is used"<<endl;
1872
1873   //backround +/- error in the range
1874   if (ftypeOfFit4Bkg == 3 && ftypeOfFit4Sgn == 0) {
1875     background = 0;
1876     errbackground = 0;
1877   }
1878   else{
1879     background=funcbkg->Integral(min,max)/(Double_t)fhistoInvMass->GetBinWidth(2);
1880     errbackground=intBerr/intB*background; // assume relative error is the same as for total integral
1881     //cout<<"integral = "<<funcbkg->Integral(min, max)<<"\tbinW = "<<fhistoInvMass->GetBinWidth(2)<<endl;
1882   }
1883   return;
1884
1885 }
1886
1887
1888 //__________________________________________________________________________
1889
1890 void AliHFMassFitter::Significance(Double_t nOfSigma,Double_t &significance,Double_t &errsignificance) const  {
1891   // Return significance in mean+- n sigma
1892
1893   Double_t min=fMass-nOfSigma*fSigmaSgn;
1894   Double_t max=fMass+nOfSigma*fSigmaSgn;
1895   Significance(min, max, significance, errsignificance);
1896
1897   return;
1898 }
1899
1900 //__________________________________________________________________________
1901
1902 void AliHFMassFitter::Significance(Double_t min, Double_t max, Double_t &significance,Double_t &errsignificance) const {
1903   // Return significance integral in a range
1904
1905   Double_t signal,errsignal,background,errbackground;
1906   Signal(min, max,signal,errsignal);
1907   Background(min, max,background,errbackground);
1908
1909   if (signal+background <= 0.){
1910     cout<<"Cannot calculate significance because of div by 0!"<<endl;
1911     significance=-1;
1912     errsignificance=0;
1913   }
1914
1915   significance =  signal/TMath::Sqrt(signal+background);
1916
1917   errsignificance = TMath::Sqrt(significance*significance/(signal+background)/(signal+background)*(1/4.*errsignal*errsignal+errbackground*errbackground)+significance*significance/signal/signal*errsignal*errsignal);
1918
1919   return;
1920 }
1921
1922