PWG3/vertexingHF/AliHFMassFitter.cxx

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