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