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