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   cout<<"Parameter "<<thispar<<" is now fixed"<<endl;
 336   return kTRUE;
 337 }
 338
 339 //___________________________________________________________________________
 340 Bool_t AliHFMassFitter::GetFixThisParam(Int_t thispar)const{
 341   //return the value of fFixPar[thispar]
 342   if(thispar>=fNFinalPars) {
 343     cout<<"Error! Parameter out of bounds! Max is "<<fNFinalPars-1<<endl;
 344     return kFALSE;
 345   }
 346   if(!fFixPar) {
 347     cout<<"Error! Parameters to be fixed still not set"<<endl;
 348     return kFALSE;
 349
 350   }
 351   return fFixPar[thispar];
 352
 353 }
 354
 355 //___________________________________________________________________________
 356 void AliHFMassFitter::SetHisto(const TH1F *histoToFit){
 357   //fhistoInvMass = (TH1F*)histoToFit->Clone();
 358   fhistoInvMass = new TH1F(*histoToFit);
 359   fhistoInvMass->SetDirectory(0);
 360   cout<<"SetHisto pointer "<<fhistoInvMass<<endl;
 361 }
 362
 363 //___________________________________________________________________________
 364
 365   void AliHFMassFitter::SetType(Int_t fittypeb, Int_t fittypes) {
 366
 367     //set the type of fit to perform for signal and background
 368
 369     ftypeOfFit4Bkg = fittypeb;
 370     ftypeOfFit4Sgn = fittypes;
 371
 372     /*
 373     if(fFitPars) {
 374       delete[] fFitPars;
 375       fFitPars=NULL;
 376     }
 377     */
 378     ComputeParSize();
 379     fFitPars = new Float_t[fParsSize];
 380     /*
 381     if(fFixPar){
 382       delete[] fFixPar;
 383       fFixPar=NULL;
 384     }
 385     */
 386     SetDefaultFixParam();
 387
 388
 389 }
 390
 391 //___________________________________________________________________________
 392
 393 void AliHFMassFitter::Reset() {
 394
 395   //delete the histogram and reset the mean and sigma to default
 396
 397   cout<<"Reset called: delete histogram, set mean value to 1.85 and sigma to 0.012"<<endl;
 398   fMass=1.85;
 399   fSigmaSgn=0.012;
 400   cout<<"Reset "<<fhistoInvMass<<endl;
 401   if(fhistoInvMass) {
 402     //cout<<"esiste"<<endl;
 403     delete fhistoInvMass;
 404     fhistoInvMass=NULL;
 405     cout<<fhistoInvMass<<endl;
 406   }
 407   else cout<<"histogram doesn't exist, do not delete"<<endl;
 408
 409
 410 }
 411
 412 //_________________________________________________________________________
 413
 414 void AliHFMassFitter::InitNtuParam(char *ntuname) {
 415
 416   // Create ntuple to keep fit parameters
 417
 418   fntuParam=0;
 419   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");
 420
 421 }
 422
 423 //_________________________________________________________________________
 424
 425 void AliHFMassFitter::FillNtuParam() {
 426   // Fill ntuple with fit parameters
 427
 428   Float_t nothing=0.;
 429
 430   if (ftypeOfFit4Bkg==2) {
 431       fntuParam->SetBranchAddress("intbkg1",&fFitPars[0]);
 432       fntuParam->SetBranchAddress("slope1",&fFitPars[1]);
 433       fntuParam->SetBranchAddress("conc1",&fFitPars[2]);
 434       fntuParam->SetBranchAddress("intGB",&fFitPars[3]);
 435       fntuParam->SetBranchAddress("meanGB",&fFitPars[4]);
 436       fntuParam->SetBranchAddress("sigmaGB",&fFitPars[5]);
 437       fntuParam->SetBranchAddress("intbkg2",&fFitPars[6]);
 438       fntuParam->SetBranchAddress("slope2",&fFitPars[7]);
 439       fntuParam->SetBranchAddress("conc2",&fFitPars[8]);
 440       fntuParam->SetBranchAddress("inttot",&fFitPars[9]);
 441       fntuParam->SetBranchAddress("slope3",&fFitPars[10]);
 442       fntuParam->SetBranchAddress("conc3",&fFitPars[11]);
 443       fntuParam->SetBranchAddress("intsgn",&fFitPars[12]);
 444       fntuParam->SetBranchAddress("meansgn",&fFitPars[13]);
 445       fntuParam->SetBranchAddress("sigmasgn",&fFitPars[14]);
 446
 447       fntuParam->SetBranchAddress("intbkg1Err",&fFitPars[15]);
 448       fntuParam->SetBranchAddress("slope1Err",&fFitPars[16]);
 449       fntuParam->SetBranchAddress("conc1Err",&fFitPars[17]);
 450       fntuParam->SetBranchAddress("intGBErr",&fFitPars[18]);
 451       fntuParam->SetBranchAddress("meanGBErr",&fFitPars[19]);
 452       fntuParam->SetBranchAddress("sigmaGBErr",&fFitPars[20]);
 453       fntuParam->SetBranchAddress("intbkg2Err",&fFitPars[21]);
 454       fntuParam->SetBranchAddress("slope2Err",&fFitPars[22]);
 455       fntuParam->SetBranchAddress("conc2Err",&fFitPars[23]);
 456       fntuParam->SetBranchAddress("inttotErr",&fFitPars[24]);
 457       fntuParam->SetBranchAddress("slope3Err",&fFitPars[25]);
 458       fntuParam->SetBranchAddress("conc3Err",&fFitPars[26]);
 459       fntuParam->SetBranchAddress("intsgnErr",&fFitPars[27]);
 460       fntuParam->SetBranchAddress("meansgnErr",&fFitPars[28]);
 461       fntuParam->SetBranchAddress("sigmasgnErr",&fFitPars[29]);
 462
 463   } else {
 464
 465     if(ftypeOfFit4Bkg==3){
 466       fntuParam->SetBranchAddress("intbkg1",&fFitPars[0]);
 467       fntuParam->SetBranchAddress("slope1",&nothing);
 468       fntuParam->SetBranchAddress("conc1",&nothing);
 469       fntuParam->SetBranchAddress("intGB",&fFitPars[1]);
 470       fntuParam->SetBranchAddress("meanGB",&fFitPars[2]);
 471       fntuParam->SetBranchAddress("sigmaGB",&fFitPars[3]);
 472       fntuParam->SetBranchAddress("intbkg2",&fFitPars[4]);
 473       fntuParam->SetBranchAddress("slope2",&nothing);
 474       fntuParam->SetBranchAddress("conc2",&nothing);
 475       fntuParam->SetBranchAddress("inttot",&fFitPars[6]);
 476       fntuParam->SetBranchAddress("slope3",&nothing);
 477       fntuParam->SetBranchAddress("conc3",&nothing);
 478       fntuParam->SetBranchAddress("intsgn",&fFitPars[6]);
 479       fntuParam->SetBranchAddress("meansgn",&fFitPars[7]);
 480       fntuParam->SetBranchAddress("sigmasgn",&fFitPars[8]);
 481
 482       fntuParam->SetBranchAddress("intbkg1Err",&fFitPars[9]);
 483       fntuParam->SetBranchAddress("slope1Err",&nothing);
 484       fntuParam->SetBranchAddress("conc1Err",&nothing);
 485       fntuParam->SetBranchAddress("intGBErr",&fFitPars[10]);
 486       fntuParam->SetBranchAddress("meanGBErr",&fFitPars[11]);
 487       fntuParam->SetBranchAddress("sigmaGBErr",&fFitPars[12]);
 488       fntuParam->SetBranchAddress("intbkg2Err",&fFitPars[13]);
 489       fntuParam->SetBranchAddress("slope2Err",&nothing);
 490       fntuParam->SetBranchAddress("conc2Err",&nothing);
 491       fntuParam->SetBranchAddress("inttotErr",&fFitPars[15]);
 492       fntuParam->SetBranchAddress("slope3Err",&nothing);
 493       fntuParam->SetBranchAddress("conc3Err",&nothing);
 494       fntuParam->SetBranchAddress("intsgnErr",&fFitPars[15]);
 495       fntuParam->SetBranchAddress("meansgnErr",&fFitPars[16]);
 496       fntuParam->SetBranchAddress("sigmasgnErr",&fFitPars[17]);
 497
 498     }
 499     else{
 500       fntuParam->SetBranchAddress("intbkg1",&fFitPars[0]);
 501       fntuParam->SetBranchAddress("slope1",&fFitPars[1]);
 502       fntuParam->SetBranchAddress("conc1",&nothing);
 503       fntuParam->SetBranchAddress("intGB",&fFitPars[2]);
 504       fntuParam->SetBranchAddress("meanGB",&fFitPars[3]);
 505       fntuParam->SetBranchAddress("sigmaGB",&fFitPars[4]);
 506       fntuParam->SetBranchAddress("intbkg2",&fFitPars[5]);
 507       fntuParam->SetBranchAddress("slope2",&fFitPars[6]);
 508       fntuParam->SetBranchAddress("conc2",&nothing);
 509       fntuParam->SetBranchAddress("inttot",&fFitPars[7]);
 510       fntuParam->SetBranchAddress("slope3",&fFitPars[8]);
 511       fntuParam->SetBranchAddress("conc3",&nothing);
 512       fntuParam->SetBranchAddress("intsgn",&fFitPars[9]);
 513       fntuParam->SetBranchAddress("meansgn",&fFitPars[10]);
 514       fntuParam->SetBranchAddress("sigmasgn",&fFitPars[11]);
 515
 516       fntuParam->SetBranchAddress("intbkg1Err",&fFitPars[12]);
 517       fntuParam->SetBranchAddress("slope1Err",&fFitPars[13]);
 518       fntuParam->SetBranchAddress("conc1Err",&nothing);
 519       fntuParam->SetBranchAddress("intGBErr",&fFitPars[14]);
 520       fntuParam->SetBranchAddress("meanGBErr",&fFitPars[15]);
 521       fntuParam->SetBranchAddress("sigmaGBErr",&fFitPars[16]);
 522       fntuParam->SetBranchAddress("intbkg2Err",&fFitPars[17]);
 523       fntuParam->SetBranchAddress("slope2Err",&fFitPars[18]);
 524       fntuParam->SetBranchAddress("conc2Err",&nothing);
 525       fntuParam->SetBranchAddress("inttotErr",&fFitPars[19]);
 526       fntuParam->SetBranchAddress("slope3Err",&fFitPars[20]);
 527       fntuParam->SetBranchAddress("conc3Err",&nothing);
 528       fntuParam->SetBranchAddress("intsgnErr",&fFitPars[21]);
 529       fntuParam->SetBranchAddress("meansgnErr",&fFitPars[22]);
 530       fntuParam->SetBranchAddress("sigmasgnErr",&fFitPars[23]);
 531     }
 532
 533   }
 534   fntuParam->TTree::Fill();
 535 }
 536
 537 //_________________________________________________________________________
 538
 539 TNtuple* AliHFMassFitter::NtuParamOneShot(char *ntuname){
 540   // Create, fill and return ntuple with fit parameters
 541
 542   InitNtuParam(ntuname);
 543   FillNtuParam();
 544   return fntuParam;
 545 }
 546 //_________________________________________________________________________
 547
 548 void AliHFMassFitter::RebinMass(Int_t bingroup){
 549   // Rebin invariant mass histogram
 550
 551   if(!fhistoInvMass){
 552     cout<<"Histogram not set"<<endl;
 553     return;
 554   }
 555   Int_t nbinshisto=fhistoInvMass->GetNbinsX();
 556   if(bingroup<1){
 557     cout<<"Error! Cannot group "<<bingroup<<" bins\n";
 558     fNbin=nbinshisto;
 559     cout<<"Kept original number of bins: "<<fNbin<<endl;
 560   } else{
 561     fhistoInvMass->Rebin(bingroup);
 562     if(nbinshisto%bingroup != 0) CheckRangeFit();
 563     fNbin = fhistoInvMass->GetNbinsX();
 564     cout<<"New number of bins: "<<fNbin<<endl;
 565   }
 566
 567 }
 568
 569 //************* fit
 570
 571 //___________________________________________________________________________
 572
 573 Double_t AliHFMassFitter::FitFunction4MassDistr (Double_t *x, Double_t *par){
 574   // Fit function for signal+background
 575
 576
 577   //exponential or linear fit
 578   //
 579   // par[0] = tot integral
 580   // par[1] = slope
 581   // par[2] = gaussian integral
 582   // par[3] = gaussian mean
 583   // par[4] = gaussian sigma
 584
 585   Double_t total,bkg=0,sgn=0;
 586
 587   if (ftypeOfFit4Bkg==0 || ftypeOfFit4Bkg==1) {
 588     if(ftypeOfFit4Sgn == 0) {
 589
 590       Double_t parbkg[2] = {par[0]-par[2], par[1]};
 591       bkg = FitFunction4Bkg(x,parbkg);
 592     }
 593     if(ftypeOfFit4Sgn == 1) {
 594       Double_t parbkg[5] = {par[2],par[3],ffactor*par[4],par[0]-2*par[2], par[1]};
 595       bkg = FitFunction4Bkg(x,parbkg);
 596     }
 597
 598     sgn = FitFunction4Sgn(x,&par[2]);
 599
 600   }
 601
 602   //polynomial fit
 603
 604     // par[0] = tot integral
 605     // par[1] = coef1
 606     // par[2] = coef2
 607     // par[3] = gaussian integral
 608     // par[4] = gaussian mean
 609     // par[5] = gaussian sigma
 610
 611   if (ftypeOfFit4Bkg==2) {
 612
 613     if(ftypeOfFit4Sgn == 0) {
 614
 615       Double_t parbkg[3] = {par[0]-par[3], par[1], par[2]};
 616       bkg = FitFunction4Bkg(x,parbkg);
 617     }
 618     if(ftypeOfFit4Sgn == 1) {
 619
 620       Double_t parbkg[6] = {par[3],par[4],ffactor*par[5],par[0]-2*par[3], par[1], par[2]};
 621      bkg = FitFunction4Bkg(x,parbkg);
 622     }
 623
 624     sgn = FitFunction4Sgn(x,&par[3]);
 625   }
 626
 627   if (ftypeOfFit4Bkg==3) {
 628
 629     if(ftypeOfFit4Sgn == 0) {
 630         bkg=FitFunction4Bkg(x,par);
 631         sgn=FitFunction4Sgn(x,&par[1]);
 632     }
 633     if(ftypeOfFit4Sgn == 1) {
 634       Double_t parbkg[4]={par[1],par[2],ffactor*par[3],par[0]};
 635       bkg=FitFunction4Bkg(x,parbkg);
 636       sgn=FitFunction4Sgn(x,&par[1]);
 637     }
 638   }
 639
 640   total = bkg + sgn;
 641
 642   return  total;
 643 }
 644
 645 //_________________________________________________________________________
 646 Double_t AliHFMassFitter::FitFunction4Sgn (Double_t *x, Double_t *par){
 647   // Fit function for the signal
 648
 649   //gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
 650   //Par:
 651   // * [0] = integralSgn
 652   // * [1] = mean
 653   // * [2] = sigma
 654   //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
 655
 656   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]);
 657
 658 }
 659
 660 //__________________________________________________________________________
 661
 662 Double_t AliHFMassFitter::FitFunction4Bkg (Double_t *x, Double_t *par){
 663   // Fit function for the background
 664
 665   Double_t maxDeltaM = 4.*fSigmaSgn;
 666   if(fSideBands && TMath::Abs(x[0]-fMass) < maxDeltaM) {
 667     TF1::RejectPoint();
 668     return 0;
 669   }
 670   Int_t firstPar=0;
 671   Double_t gaus2=0,total=-1;
 672   if(ftypeOfFit4Sgn == 1){
 673     firstPar=3;
 674     //gaussian = A/(sigma*sqrt(2*pi))*exp(-(x-mean)^2/2/sigma^2)
 675     //Par:
 676     // * [0] = integralSgn
 677     // * [1] = mean
 678     // * [2] = sigma
 679     //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
 680     gaus2 = FitFunction4Sgn(x,par);
 681   }
 682
 683   switch (ftypeOfFit4Bkg){
 684   case 0:
 685     //exponential
 686     //exponential = A*exp(B*x) -> integral(exponential)=A/B*exp(B*x)](min,max)
 687     //-> A = B*integral/(exp(B*max)-exp(B*min)) where integral can be written
 688     //as integralTot- integralGaus (=par [2])
 689     //Par:
 690     // * [0] = integralBkg;
 691     // * [1] = B;
 692     //exponential = [1]*[0]/(exp([1]*max)-exp([1]*min))*exp([1]*x)
 693     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]);
 694     break;
 695   case 1:
 696     //linear
 697     //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)
 698     // * [0] = integralBkg;
 699     // * [1] = b;
 700     total= par[0+firstPar]/(fmaxMass-fminMass)+par[1+firstPar]*(x[0]-0.5*(fmaxMass+fminMass));
 701     break;
 702   case 2:
 703     //polynomial
 704     //y=a+b*x+c*x**2 -> integral = a(max-min) + 1/2*b*(max^2-min^2) +
 705     //+ 1/3*c*(max^3-min^3) ->
 706     //a = (integral-1/2*b*(max^2-min^2)-1/3*c*(max^3-min^3))/(max-min)
 707     // * [0] = integralBkg;
 708     // * [1] = b;
 709     // * [2] = c;
 710     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));
 711     break;
 712   case 3:
 713     total=par[0+firstPar];
 714     break;
 715 //   default:
 716 //     Types of Fit Functions for Background:
 717 //     * 0 = exponential;
 718 //     * 1 = linear;
 719 //     * 2 = polynomial 2nd order
 720 //     * 3 = no background"<<endl;
 721
 722   }
 723   return total+gaus2;
 724 }
 725
 726 //__________________________________________________________________________
 727 Bool_t AliHFMassFitter::SideBandsBounds(){
 728
 729   //determines the ranges of the side bands
 730
 731   Double_t width=fhistoInvMass->GetBinWidth(8);
 732   if (fNbin==0) fNbin=fhistoInvMass->GetNbinsX();
 733   Double_t minHisto=fhistoInvMass->GetBinLowEdge(1);
 734   Double_t maxHisto=fhistoInvMass->GetBinLowEdge(fNbin)+width;
 735
 736   if(fMass-fminMass < 0 || fmaxMass-fMass <0) {
 737     cout<<"Left limit of range > mean or right limit of range < mean: change left/right limit or initial mean value"<<endl;
 738     return kFALSE;
 739   }
 740
 741   if((TMath::Abs(fminMass-minHisto) < 10e6 || TMath::Abs(fmaxMass - maxHisto) < 10e6) && (fMass-4.*fSigmaSgn-fminMass) < 10e6){
 742     Double_t coeff = (fMass-fminMass)/fSigmaSgn;
 743
 744     fSideBandl=(Int_t)((fMass-0.5*coeff*fSigmaSgn-fminMass)/width);
 745     fSideBandr=(Int_t)((fMass+0.5*coeff*fSigmaSgn-fminMass)/width);
 746     cout<<"Changed number of sigma from 4 to "<<0.5*coeff<<" for the estimation of the side bands"<<endl;
 747     if (coeff<3) cout<<"Side bands inside 3 sigma, may be better use ftypeOfFit4Bkg = 3 (only signal)"<<endl;
 748
 749     if (coeff<2) {
 750       cout<<"Side bands inside 2 sigma. Change mode: ftypeOfFit4Bkg = 3"<<endl;
 751       ftypeOfFit4Bkg=3;
 752       //set binleft and right without considering SetRangeFit- anyway no bkg!
 753       fSideBandl=(Int_t)((fMass-4.*fSigmaSgn-minHisto)/width);
 754       fSideBandr=(Int_t)((fMass+4.*fSigmaSgn-minHisto)/width);
 755     }
 756   }
 757   else {
 758   fSideBandl=(Int_t)((fMass-4.*fSigmaSgn-minHisto)/width);
 759   fSideBandr=(Int_t)((fMass+4.*fSigmaSgn-minHisto)/width);
 760 //   cout<<"\tfMass = "<<fMass<<"\tfSigmaSgn = "<<fSigmaSgn<<"\tminHisto = "<<minHisto<<endl;
 761 //   cout<<"\tbinleft = "<<fSideBandl<<"\tbinright = "<<fSideBandr<<endl;
 762   }
 763
 764   if (fSideBandl==0) {
 765     cout<<"Error! Range too little";
 766     return kFALSE;
 767   }
 768   return kTRUE;
 769 }
 770
 771 //__________________________________________________________________________
 772
 773 void AliHFMassFitter::GetSideBandsBounds(Int_t &left, Int_t &right) const{
 774
 775   // get the range of the side bands
 776
 777   if (fSideBandl==0 && fSideBandr==0){
 778     cout<<"Use MassFitter method first"<<endl;
 779     return;
 780   }
 781   left=fSideBandl;
 782   right=fSideBandr;
 783 }
 784
 785 //__________________________________________________________________________
 786 Bool_t AliHFMassFitter::CheckRangeFit(){
 787   //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
 788
 789   if (!fhistoInvMass) {
 790     cout<<"No histogram to fit! SetHisto(TH1F* h) before! "<<endl;
 791     return kFALSE;
 792   }
 793   Bool_t leftok=kFALSE, rightok=kFALSE;
 794   Int_t nbins=fhistoInvMass->GetNbinsX();
 795   Double_t width=fhistoInvMass->GetBinWidth(1);
 796   Double_t minhisto=fhistoInvMass->GetBinLowEdge(1), maxhisto=fhistoInvMass->GetBinLowEdge(nbins)+width;
 797
 798   //check if limits are inside histogram range
 799   if( fminMass-minhisto < 0. ) {
 800     cout<<"Out of histogram left bound!"<<endl;
 801     fminMass=minhisto;
 802   }
 803   if( fmaxMass-maxhisto > 0. ) {
 804     cout<<"Out of histogram right bound!"<<endl;
 805     fmaxMass=maxhisto;
 806   }
 807
 808   Int_t binl,binr;
 809   Double_t tmp=0.;
 810   tmp=fminMass;
 811   //calculate bin corresponding to fminMass
 812   binl=fhistoInvMass->FindBin(fminMass);
 813   if (fminMass > fhistoInvMass->GetBinCenter(binl)) binl++;
 814   fminMass=fhistoInvMass->GetBinLowEdge(binl);
 815   if(TMath::Abs(tmp-fminMass) > 1e-6){
 816     cout<<"Left bound "<<tmp<<" is not allowed, changing it to the nearest value allowed: "<<fminMass<<endl;
 817     leftok=kTRUE;
 818   }
 819
 820   tmp=fmaxMass;
 821   //calculate bin corresponding to fmaxMass
 822   binr=fhistoInvMass->FindBin(fmaxMass);
 823   if (fmaxMass < fhistoInvMass->GetBinCenter(binr)) binr--;
 824   fmaxMass=fhistoInvMass->GetBinLowEdge(binr)+width;
 825   if(TMath::Abs(tmp-fmaxMass) > 1e-6){
 826     cout<<"Right bound "<<tmp<<" is not allowed, changing it to the nearest value allowed: "<<fmaxMass<<endl;
 827     rightok=kTRUE;
 828   }
 829
 830   return (leftok && rightok);
 831
 832 }
 833
 834 //__________________________________________________________________________
 835
 836 Bool_t AliHFMassFitter::MassFitter(Bool_t draw){
 837   // Main method of the class: performs the fit of the histogram
 838
 839   //Set default fitter Minuit in order to use gMinuit in the contour plots
 840   TVirtualFitter::SetDefaultFitter("Minuit");
 841
 842   Int_t nFitPars=0; //total function's number of parameters
 843   switch (ftypeOfFit4Bkg){
 844   case 0:
 845     nFitPars=5; //3+2
 846     break;
 847   case 1:
 848     nFitPars=5; //3+2
 849     break;
 850   case 2:
 851     nFitPars=6; //3+3
 852     break;
 853   case 3:
 854     nFitPars=4; //3+1
 855     break;
 856   }
 857
 858   Int_t bkgPar = nFitPars-3; //background function's number of parameters
 859
 860   cout<<"nFitPars = "<<nFitPars<<"\tbkgPar = "<<bkgPar<<endl;
 861
 862
 863   TString listname="contourplot";
 864   listname+=fcounter;
 865   if(!fContourGraph){
 866     fContourGraph=new TList();
 867     fContourGraph->SetOwner();
 868   }
 869
 870   fContourGraph->SetName(listname);
 871
 872
 873   //function names
 874   TString bkgname="funcbkg";
 875   TString bkg1name="funcbkg1";
 876   TString massname="funcmass";
 877
 878   //Total integral
 879   Double_t totInt = fhistoInvMass->Integral(fhistoInvMass->FindBin(fminMass), fhistoInvMass->FindBin(fmaxMass), "width");
 880
 881   fSideBands = kTRUE;
 882   Double_t width=fhistoInvMass->GetBinWidth(8);
 883   cout<<"fNbin"<<fNbin<<endl;
 884   if (fNbin==0) fNbin=fhistoInvMass->GetNbinsX();
 885   //Double_t minHisto=fhistoInvMass->GetBinLowEdge(1);
 886   //Double_t maxHisto=fhistoInvMass->GetBinLowEdge(fNbin)+width;
 887   Bool_t ok=SideBandsBounds();
 888   if(!ok) return kFALSE;
 889
 890   //sidebands integral - first approx (from histo)
 891   Double_t sideBandsInt=(Double_t)fhistoInvMass->Integral(1,fSideBandl,"width") + (Double_t)fhistoInvMass->Integral(fSideBandr,fNbin,"width");
 892   cout<<"------nbin = "<<fNbin<<"\twidth = "<<width<<"\tbinleft = "<<fSideBandl<<"\tbinright = "<<fSideBandr<<endl;
 893   cout<<"------sideBandsInt - first approx = "<<sideBandsInt<<endl;
 894   if (sideBandsInt<=0) {
 895     cout<<"! sideBandsInt <=0. There's a problem, cannot start the fit"<<endl;
 896     return kFALSE;
 897   }
 898
 899   /*Fit Bkg*/
 900
 901
 902   TF1 *funcbkg = new TF1(bkgname.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,bkgPar,"AliHFMassFitter","FitFunction4Bkg");
 903   cout<<"Function name = "<<funcbkg->GetName()<<endl<<endl;
 904
 905   funcbkg->SetLineColor(2); //red
 906
 907   //first fit for bkg: approx bkgint
 908
 909   switch (ftypeOfFit4Bkg) {
 910   case 0: //gaus+expo
 911     funcbkg->SetParNames("BkgInt","Slope");
 912     funcbkg->SetParameters(sideBandsInt,-2.);
 913     break;
 914   case 1:
 915     funcbkg->SetParNames("BkgInt","Slope");
 916     funcbkg->SetParameters(sideBandsInt,-100.);
 917     break;
 918   case 2:
 919     funcbkg->SetParNames("BkgInt","Coef1","Coef2");
 920     funcbkg->SetParameters(sideBandsInt,-10.,5);
 921     break;
 922   case 3:
 923     if(ftypeOfFit4Sgn==0){
 924       funcbkg->SetParNames("Const");
 925       funcbkg->SetParameter(0,0.);
 926       funcbkg->FixParameter(0,0.);
 927     }
 928     break;
 929   default:
 930     cout<<"Wrong choise of ftypeOfFit4Bkg ("<<ftypeOfFit4Bkg<<")"<<endl;
 931     return kFALSE;
 932     break;
 933   }
 934   cout<<"\nBACKGROUND FIT - only combinatorial"<<endl;
 935   Int_t ftypeOfFit4SgnBkp=ftypeOfFit4Sgn;
 936
 937   Double_t intbkg1=0,slope1=0,conc1=0;
 938   //if only signal and reflection: skip
 939   if (!(ftypeOfFit4Bkg==3 && ftypeOfFit4Sgn==1)) {
 940     ftypeOfFit4Sgn=0;
 941     fhistoInvMass->Fit(bkgname.Data(),"R,L,E,0");
 942
 943     for(Int_t i=0;i<bkgPar;i++){
 944       fFitPars[i]=funcbkg->GetParameter(i);
 945       //cout<<i<<"\t"<<funcbkg->GetParameter(i)<<"\t";
 946       fFitPars[nFitPars+2*bkgPar+3+i]= funcbkg->GetParError(i);
 947       //cout<<nFitPars+2*bkgPar+3+i<<"\t"<< funcbkg->GetParError(i)<<endl;
 948     }
 949     fSideBands = kFALSE;
 950     //intbkg1 = funcbkg->GetParameter(0);
 951     funcbkg->SetRange(fminMass,fmaxMass);
 952     intbkg1 = funcbkg->Integral(fminMass,fmaxMass);
 953     if(ftypeOfFit4Bkg!=3) slope1 = funcbkg->GetParameter(1);
 954     if(ftypeOfFit4Bkg==2) conc1 = funcbkg->GetParameter(2);
 955     cout<<"First fit: \nintbkg1 = "<<intbkg1<<"\t(Compare with par0 = "<<funcbkg->GetParameter(0)<<")\nslope1= "<<slope1<<"\nconc1 = "<<conc1<<endl;
 956   }
 957   else cout<<"\t\t//"<<endl;
 958
 959   ftypeOfFit4Sgn=ftypeOfFit4SgnBkp;
 960   TF1 *funcbkg1=0;
 961   if (ftypeOfFit4Sgn == 1) {
 962     cout<<"\nBACKGROUND FIT WITH REFLECTION"<<endl;
 963     bkgPar+=3;
 964
 965     cout<<"nFitPars = "<<nFitPars<<"\tbkgPar = "<<bkgPar<<endl;
 966
 967     funcbkg1 = new TF1(bkg1name.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,bkgPar,"AliHFMassFitter","FitFunction4Bkg");
 968     cout<<"Function name = "<<funcbkg1->GetName()<<endl;
 969
 970     funcbkg1->SetLineColor(2); //red
 971
 972     if(ftypeOfFit4Bkg==2){
 973       cout<<"*** Polynomial Fit ***"<<endl;
 974       funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Coef1","Coef2");
 975       funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1,conc1);
 976
 977       //cout<<"Parameters set to: "<<0.5*(totInt-intbkg1)<<"\t"<<fMass<<"\t"<<ffactor*fSigmaSgn<<"\t"<<intbkg1<<"\t"<<slope1<<"\t"<<conc1<<"\t"<<endl;
 978       //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<bkgPar<<"\tgsidebands = "<<fSideBands<<endl;
 979     } else{
 980       if(ftypeOfFit4Bkg==3) //no background: gaus sign+ gaus broadened
 981         {
 982           cout<<"*** No background Fit ***"<<endl;
 983           funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","Const");
 984           funcbkg1->SetParameters(0.5*totInt,fMass,ffactor*fSigmaSgn,0.);
 985           funcbkg1->FixParameter(3,0.);
 986         } else{ //expo or linear
 987           if(ftypeOfFit4Bkg==0) cout<<"*** Exponential Fit ***"<<endl;
 988           if(ftypeOfFit4Bkg==1) cout<<"*** Linear Fit ***"<<endl;
 989           funcbkg1->SetParNames("IntGB","MeanGB","SigmaGB","BkgInt","Slope");
 990           funcbkg1->SetParameters(0.5*(totInt-intbkg1),fMass,ffactor*fSigmaSgn,intbkg1,slope1);
 991         }
 992     }
 993     Int_t status=fhistoInvMass->Fit(bkg1name.Data(),"R,L,E,+,0");
 994     if (status != 0){
 995       cout<<"Minuit returned "<<status<<endl;
 996       return kFALSE;
 997     }
 998
 999     for(Int_t i=0;i<bkgPar;i++){
1000       fFitPars[bkgPar-3+i]=funcbkg1->GetParameter(i);
1001       //cout<<bkgPar-3+i<<"\t"<<funcbkg1->GetParameter(i);
1002       fFitPars[nFitPars+3*bkgPar-6+i]= funcbkg1->GetParError(i);
1003       //cout<<"\t"<<nFitPars+3*bkgPar-6+i<<"\t"<<funcbkg1->GetParError(i)<<endl;
1004     }
1005
1006     intbkg1=funcbkg1->GetParameter(3);
1007     if(ftypeOfFit4Bkg!=3) slope1 = funcbkg1->GetParameter(4);
1008     if(ftypeOfFit4Bkg==2) conc1 = funcbkg1->GetParameter(5);
1009
1010   } else {
1011     bkgPar+=3;
1012
1013     for(Int_t i=0;i<3;i++){
1014       fFitPars[bkgPar-3+i]=0.;
1015       cout<<bkgPar-3+i<<"\t"<<0.<<"\t";
1016       fFitPars[nFitPars+3*bkgPar-6+i]= 0.;
1017       cout<<nFitPars+3*bkgPar-6+i<<"\t"<<0.<<endl;
1018     }
1019
1020     for(Int_t i=0;i<bkgPar-3;i++){
1021       fFitPars[bkgPar+i]=funcbkg->GetParameter(i);
1022       cout<<bkgPar+i<<"\t"<<funcbkg->GetParameter(i)<<"\t";
1023       fFitPars[nFitPars+3*bkgPar-3+i]= funcbkg->GetParError(i);
1024       cout<<nFitPars+3*bkgPar-3+i<<"\t"<< funcbkg->GetParError(i)<<endl;
1025     }
1026
1027
1028   }
1029
1030   //sidebands integral - second approx (from fit)
1031   fSideBands = kFALSE;
1032   Double_t bkgInt;
1033   cout<<"Compare intbkg1 = "<<intbkg1<<" and integral = ";
1034   if(ftypeOfFit4Sgn == 1) bkgInt=funcbkg1->Integral(fminMass,fmaxMass);
1035   else bkgInt=funcbkg->Integral(fminMass,fmaxMass);
1036   cout<</*"------BkgInt(Fit) = "<<*/bkgInt<<endl;
1037
1038   //Signal integral - first approx
1039   Double_t sgnInt;
1040   sgnInt = totInt-bkgInt;
1041   cout<<"------TotInt = "<<totInt<<"\tsgnInt = "<<sgnInt<<endl;
1042   if (sgnInt <= 0){
1043     cout<<"Setting sgnInt = - sgnInt"<<endl;
1044     sgnInt=- sgnInt;
1045   }
1046   /*Fit All Mass distribution with exponential + gaussian (+gaussiam braodened) */
1047   TF1 *funcmass = new TF1(massname.Data(),this,&AliHFMassFitter::FitFunction4MassDistr,fminMass,fmaxMass,nFitPars,"AliHFMassFitter","FitFunction4MassDistr");
1048   cout<<"Function name = "<<funcmass->GetName()<<endl<<endl;
1049   funcmass->SetLineColor(4); //blue
1050
1051   //Set parameters
1052   cout<<"\nTOTAL FIT"<<endl;
1053
1054   if(nFitPars==5){
1055     funcmass->SetParNames("TotInt","Slope","SgnInt","Mean","Sigma");
1056     funcmass->SetParameters(totInt,slope1,sgnInt,fMass,fSigmaSgn);
1057
1058     //cout<<"Parameters set to: "<<totInt<<"\t"<<slope1<<"\t"<<sgnInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl;
1059     //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<nFitPars<<"\tgsidebands = "<<fSideBands<<endl;
1060     if(fFixPar[0]){
1061       cout<<"fix1"<<endl;
1062       funcmass->FixParameter(0,totInt);
1063     }
1064     if(fFixPar[1]){
1065       cout<<"fix2"<<endl;
1066       funcmass->FixParameter(1,slope1);
1067     }
1068     if(fFixPar[2]){
1069       cout<<"fix3"<<endl;
1070       funcmass->FixParameter(2,sgnInt);
1071     }
1072     if(fFixPar[3]){
1073       cout<<"fix4"<<endl;
1074       funcmass->FixParameter(3,fMass);
1075     }
1076     if(fFixPar[4]){
1077       cout<<"fix5"<<endl;
1078       funcmass->FixParameter(4,fSigmaSgn);
1079     }
1080   }
1081   if (nFitPars==6){
1082     funcmass->SetParNames("TotInt","Coef1","Coef2","SgnInt","Mean","Sigma");
1083     funcmass->SetParameters(totInt,slope1,conc1,sgnInt,fMass,fSigmaSgn);
1084
1085     //cout<<"Parameters set to: "<<totInt<<"\t"<<slope1<<"\t"<<conc1<<"\t"<<sgnInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl;
1086     //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<nFitPars<<"\tgsidebands = "<<fSideBands<<endl;
1087     if(fFixPar[0])funcmass->FixParameter(0,totInt);
1088     if(fFixPar[1])funcmass->FixParameter(1,slope1);
1089     if(fFixPar[2])funcmass->FixParameter(2,conc1);
1090     if(fFixPar[3])funcmass->FixParameter(3,sgnInt);
1091     if(fFixPar[4])funcmass->FixParameter(4,fMass);
1092     if(fFixPar[5])funcmass->FixParameter(5,fSigmaSgn);
1093     //
1094     //funcmass->FixParameter(2,sgnInt);
1095   }
1096   if(nFitPars==4){
1097     funcmass->SetParNames("Const","SgnInt","Mean","Sigma");
1098     if(ftypeOfFit4Sgn == 1) funcmass->SetParameters(0.,0.5*totInt,fMass,fSigmaSgn);
1099     else funcmass->SetParameters(0.,totInt,fMass,fSigmaSgn);
1100     if(fFixPar[0]) funcmass->FixParameter(0,0.);
1101     //cout<<"Parameters set to: "<<0.5*totInt<<"\t"<<fMass<<"\t"<<fSigmaSgn<<"\t"<<endl;
1102     //cout<<"Limits: ("<<fminMass<<","<<fmaxMass<<")\tnPar = "<<nFitPars<<"\tgsidebands = "<<fSideBands<<endl;
1103
1104   }
1105   //funcmass->FixParameter(nFitPars-2,fMass);
1106   //funcmass->SetParLimits(nFitPars-1,0.007,0.05);
1107     //funcmass->SetParLimits(nFitPars-2,fMass-0.01,fMass+0.01);
1108   Int_t status;
1109
1110   status = fhistoInvMass->Fit(massname.Data(),"R,L,E,+,0");
1111   if (status != 0){
1112     cout<<"Minuit returned "<<status<<endl;
1113     return kFALSE;
1114   }
1115
1116   cout<<"fit done"<<endl;
1117   //reset value of fMass and fSigmaSgn to those found from fit
1118   fMass=funcmass->GetParameter(nFitPars-2);
1119   fSigmaSgn=funcmass->GetParameter(nFitPars-1);
1120
1121   for(Int_t i=0;i<nFitPars;i++){
1122     fFitPars[i+2*bkgPar-3]=funcmass->GetParameter(i);
1123     fFitPars[nFitPars+4*bkgPar-6+i]= funcmass->GetParError(i);
1124     //cout<<i+2*bkgPar-3<<"\t"<<funcmass->GetParameter(i)<<"\t\t"<<nFitPars+4*bkgPar-6+i<<"\t"<<funcmass->GetParError(i)<<endl;
1125   }
1126   /*
1127   //check: cout parameters
1128   for(Int_t i=0;i<2*(nFitPars+2*bkgPar-3);i++){
1129     cout<<i<<"\t"<<fFitPars[i]<<endl;
1130     }
1131   */
1132
1133 //   if(draw){
1134 //     TCanvas *canvas=new TCanvas(fhistoInvMass->GetName(),fhistoInvMass->GetName());
1135 //     TH1F *fhistocopy=new TH1F(*fhistoInvMass);
1136 //     canvas->cd();
1137 //     fhistocopy->DrawClone();
1138 //     if(ftypeOfFit4Sgn == 1) funcbkg1->DrawClone("sames");
1139 //     else funcbkg->DrawClone("sames");
1140 //     funcmass->DrawClone("sames");
1141 //     cout<<"Drawn"<<endl;
1142 //   }
1143
1144   if(funcmass->GetParameter(nFitPars-1) <0 || funcmass->GetParameter(nFitPars-2) <0 || funcmass->GetParameter(nFitPars-3) <0 ) {
1145     cout<<"IntS or mean or sigma negative. You may tray to SetInitialGaussianSigma(..) and SetInitialGaussianMean(..)"<<endl;
1146     return kFALSE;
1147   }
1148
1149   //increase counter of number of fits done
1150   fcounter++;
1151
1152   //contour plots
1153   if(draw){
1154
1155     for (Int_t kpar=1; kpar<nFitPars;kpar++){
1156
1157       for(Int_t jpar=kpar+1;jpar<nFitPars;jpar++){
1158         cout<<"Par "<<kpar<<" and "<<jpar<<endl;
1159
1160         // produce 2 contours per couple of parameters
1161         TGraph* cont[2] = {0x0, 0x0};
1162         const Double_t errDef[2] = {1., 4.};
1163         for (Int_t i=0; i<2; i++) {
1164           gMinuit->SetErrorDef(errDef[i]);
1165           cont[i] = (TGraph*)gMinuit->Contour(80,kpar,jpar);
1166           cout<<"Minuit Status = "<<gMinuit->GetStatus()<<endl;
1167         }
1168
1169         if(!cont[0] || !cont[1]){
1170           cout<<"Skipping par "<<kpar<<" vs par "<<jpar<<endl;
1171           continue;
1172         }
1173
1174         // set graph titles and add them to the list
1175         TString title = "Contour plot";
1176         TString titleX = funcmass->GetParName(kpar);
1177         TString titleY = funcmass->GetParName(jpar);
1178         for (Int_t i=0; i<2; i++) {
1179           cont[i]->SetName( Form("cperr%d_%d%d", i, kpar, jpar) );
1180           cont[i]->SetTitle(title);
1181           cont[i]->GetXaxis()->SetTitle(titleX);
1182           cont[i]->GetYaxis()->SetTitle(titleY);
1183           cont[i]->GetYaxis()->SetLabelSize(0.033);
1184           cont[i]->GetYaxis()->SetTitleSize(0.033);
1185           cont[i]->GetYaxis()->SetTitleOffset(1.67);
1186
1187           fContourGraph->Add(cont[i]);
1188         }
1189
1190         // plot them
1191         TString cvname = Form("c%d%d", kpar, jpar);
1192         TCanvas *c4=new TCanvas(cvname,cvname,600,600);
1193         c4->cd();
1194         cont[1]->SetFillColor(38);
1195         cont[1]->Draw("alf");
1196         cont[0]->SetFillColor(9);
1197         cont[0]->Draw("lf");
1198
1199       }
1200
1201     }
1202
1203   }
1204
1205   if (ftypeOfFit4Sgn == 1 && funcbkg1) {
1206     delete funcbkg1;
1207     funcbkg1=NULL;
1208   }
1209   if (funcbkg) {
1210     delete funcbkg;
1211     funcbkg=NULL;
1212   }
1213   if (funcmass) {
1214     delete funcmass;
1215     funcmass=NULL;
1216   }
1217
1218   AddFunctionsToHisto();
1219   if (draw) DrawFit();
1220
1221
1222   return kTRUE;
1223 }
1224 //_________________________________________________________________________
1225 Double_t AliHFMassFitter::GetChiSquare() const{
1226   TF1 *funcmass=(TF1*)fhistoInvMass->GetFunction("funcmass");
1227   return funcmass->GetChisquare();
1228 }
1229
1230 //_________________________________________________________________________
1231 Double_t AliHFMassFitter::GetReducedChiSquare() const{
1232   TF1 *funcmass=(TF1*)fhistoInvMass->GetFunction("funcmass");
1233   return funcmass->GetChisquare()/funcmass->GetNDF();
1234 }
1235
1236 //*********output
1237
1238 //_________________________________________________________________________
1239 void  AliHFMassFitter::GetFitPars(Float_t *vector) const {
1240   // Return fit parameters
1241
1242   for(Int_t i=0;i<fParsSize;i++){
1243     vector[i]=fFitPars[i];
1244   }
1245 }
1246
1247
1248 //_________________________________________________________________________
1249 void AliHFMassFitter::IntS(Float_t *valuewitherror) const {
1250
1251   //gives the integral of signal obtained from fit parameters
1252   if(!valuewitherror)valuewitherror=new Float_t[2];
1253
1254   Int_t index=fParsSize/2 - 3;
1255   valuewitherror[0]=fFitPars[index];
1256   index=fParsSize - 3;
1257   valuewitherror[1]=fFitPars[index];
1258   }
1259
1260
1261 //_________________________________________________________________________
1262 void AliHFMassFitter::AddFunctionsToHisto(){
1263
1264   //Add the background function in the complete range to the list of functions attached to the histogram
1265
1266   cout<<"AddFunctionsToHisto called"<<endl;
1267   TString bkgname = "funcbkg";
1268   Int_t np=-99;
1269   switch (ftypeOfFit4Bkg){
1270   case 0: //expo
1271     np=2;
1272     break;
1273   case 1: //linear
1274     np=2;
1275     break;
1276   case 2: //pol2
1277     np=3;
1278     break;
1279   case 3: //no bkg
1280     np=1;
1281     break;
1282   }
1283   if (ftypeOfFit4Sgn == 1) {
1284     bkgname += 1;
1285     np+=3;
1286   }
1287
1288   TString bkgnamesave=bkgname;
1289   TString testname=bkgname;
1290   testname += "FullRange";
1291   TF1 *testfunc=(TF1*)fhistoInvMass->FindObject(testname.Data());
1292   if(testfunc){
1293     cout<<"AddFunctionsToHisto already used: exiting...."<<endl;
1294     return;
1295   }
1296
1297   TList *hlist=fhistoInvMass->GetListOfFunctions();
1298   hlist->ls();
1299
1300   TF1 *b=(TF1*)hlist->FindObject(bkgname.Data());
1301   if(!b){
1302     cout<<bkgname<<" not found, cannot produce "<<bkgname<<"FullRange and "<<bkgname<<"Recalc"<<endl;
1303     return;
1304   }
1305
1306   bkgname += "FullRange";
1307   TF1 *bfullrange=new TF1(bkgname.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,np,"AliHFMassFitter","FitFunction4Bkg");
1308   //cout<<bfullrange->GetName()<<endl;
1309   for(Int_t i=0;i<np;i++){
1310     //cout<<i<<" di "<<np<<endl;
1311     bfullrange->SetParName(i,b->GetParName(i));
1312     bfullrange->SetParameter(i,b->GetParameter(i));
1313     bfullrange->SetParError(i,b->GetParError(i));
1314   }
1315   bfullrange->SetLineStyle(4);
1316   bfullrange->SetLineColor(14);
1317
1318   bkgnamesave += "Recalc";
1319
1320   TF1 *blastpar=new TF1(bkgnamesave.Data(),this,&AliHFMassFitter::FitFunction4Bkg,fminMass,fmaxMass,np,"AliHFMassFitter","FitFunction4Bkg");
1321
1322   TF1 *mass=fhistoInvMass->GetFunction("funcmass");
1323
1324   if (!mass){
1325     cout<<"funcmass doesn't exist "<<endl;
1326     return;
1327   }
1328
1329   blastpar->SetParameter(0,mass->GetParameter(0)-mass->GetParameter(np));
1330   blastpar->SetParError(0,mass->GetParError(np));
1331   if (np>=2) {
1332     blastpar->SetParameter(1,mass->GetParameter(1));
1333     blastpar->SetParError(1,mass->GetParError(1));
1334   }
1335   if (np==3) {
1336     blastpar->SetParameter(2,mass->GetParameter(2));
1337     blastpar->SetParError(2,mass->GetParError(2));
1338   }
1339
1340   blastpar->SetLineStyle(1);
1341   blastpar->SetLineColor(2);
1342
1343   hlist->Add((TF1*)bfullrange->Clone());
1344   hlist->Add((TF1*)blastpar->Clone());
1345   hlist->ls();
1346
1347   if(bfullrange) {
1348     delete bfullrange;
1349     bfullrange=NULL;
1350   }
1351   if(blastpar) {
1352     delete blastpar;
1353     blastpar=NULL;
1354   }
1355 }
1356
1357 //_________________________________________________________________________
1358
1359 TH1F* AliHFMassFitter::GetHistoClone() const{
1360
1361   TH1F* hout=(TH1F*)fhistoInvMass->Clone(fhistoInvMass->GetName());
1362   return hout;
1363 }
1364 //_________________________________________________________________________
1365
1366 void AliHFMassFitter::WriteHisto(TString path) const {
1367
1368   //Write the histogram in the default file HFMassFitterOutput.root
1369
1370   if (fcounter == 0) {
1371     cout<<"Use MassFitter method before WriteHisto"<<endl;
1372     return;
1373   }
1374   TH1F* hget=(TH1F*)fhistoInvMass->Clone();
1375
1376   path += "HFMassFitterOutput.root";
1377   TFile *output;
1378
1379   if (fcounter == 1) output = new TFile(path.Data(),"recreate");
1380   else output = new TFile(path.Data(),"update");
1381   output->cd();
1382   hget->Write();
1383   fContourGraph->Write();
1384
1385
1386   output->Close();
1387
1388   cout<<fcounter<<" "<<hget->GetName()<<" written in "<<path<<endl;
1389
1390   if(output) {
1391     delete output;
1392     output=NULL;
1393   }
1394
1395 }
1396
1397 //_________________________________________________________________________
1398
1399 void AliHFMassFitter::WriteNtuple(TString path) const{
1400   //TNtuple* nget=(TNtuple*)fntuParam->Clone();
1401   path += "HFMassFitterOutput.root";
1402   TFile *output = new TFile(path.Data(),"update");
1403   output->cd();
1404   fntuParam->Write();
1405   //nget->Write();
1406   output->Close();
1407   //cout<<nget->GetName()<<" written in "<<path<<endl;
1408   cout<<fntuParam->GetName()<<" written in "<<path<<endl;
1409   /*
1410   if(nget) {
1411     //delete nget;
1412     nget=NULL;
1413   }
1414   */
1415   if(output) {
1416     delete output;
1417     output=NULL;
1418   }
1419 }
1420
1421 //_________________________________________________________________________
1422 void AliHFMassFitter::WriteCanvas(TString userIDstring,TString path,Double_t nsigma,Int_t writeFitInfo,Bool_t draw) const{
1423
1424  //write the canvas in a root file
1425
1426   gStyle->SetOptStat(0);
1427   gStyle->SetCanvasColor(0);
1428   gStyle->SetFrameFillColor(0);
1429
1430   TString type="";
1431
1432   switch (ftypeOfFit4Bkg){
1433   case 0:
1434     type="Exp"; //3+2
1435     break;
1436   case 1:
1437     type="Lin"; //3+2
1438     break;
1439   case 2:
1440     type="Pl2"; //3+3
1441     break;
1442   case 3:
1443     type="noB"; //3+1
1444     break;
1445   }
1446
1447   TString filename=Form("%sMassFit.root",type.Data());
1448   filename.Prepend(userIDstring);
1449   path.Append(filename);
1450
1451   TFile* outputcv=new TFile(path.Data(),"update");
1452
1453   TCanvas* c=(TCanvas*)GetPad(nsigma,writeFitInfo);
1454   c->SetName(Form("%s%s%s",c->GetName(),userIDstring.Data(),type.Data()));
1455   if(draw)c->DrawClone();
1456   outputcv->cd();
1457   c->Write();
1458   outputcv->Close();
1459 }
1460
1461 //_________________________________________________________________________
1462
1463 TVirtualPad* AliHFMassFitter::GetPad(Double_t nsigma,Int_t writeFitInfo)const{
1464   //return a TVirtualPad with the fitted histograms and info
1465
1466   TString cvtitle="fit of ";
1467   cvtitle+=fhistoInvMass->GetName();
1468   TString cvname="c";
1469   cvname+=fcounter;
1470
1471   TCanvas *c=new TCanvas(cvname,cvtitle);
1472   PlotFit(c->cd(),nsigma,writeFitInfo);
1473   return c->cd();
1474 }
1475 //_________________________________________________________________________
1476
1477 void AliHFMassFitter::PlotFit(TVirtualPad* pd,Double_t nsigma,Int_t writeFitInfo)const{
1478   //plot histogram, fit functions and write parameters according to verbosity level (0,1,>1)
1479   gStyle->SetOptStat(0);
1480   gStyle->SetCanvasColor(0);
1481   gStyle->SetFrameFillColor(0);
1482
1483   cout<<"nsigma = "<<nsigma<<endl;
1484   cout<<"Verbosity = "<<writeFitInfo<<endl;
1485
1486   TH1F* hdraw=GetHistoClone();
1487   hdraw->SetMinimum(0);
1488   hdraw->GetXaxis()->SetRangeUser(fminMass,fmaxMass);
1489   pd->cd();
1490   hdraw->SetMarkerStyle(20);
1491   hdraw->DrawClone("PE");
1492   hdraw->GetFunction("funcbkgFullRange")->DrawClone("same");
1493   hdraw->GetFunction("funcbkgRecalc")->DrawClone("same");
1494   hdraw->GetFunction("funcmass")->DrawClone("same");
1495
1496   if(writeFitInfo > 0){
1497     TPaveText *pinfob=new TPaveText(0.6,0.86,1.,1.,"NDC");
1498     TPaveText *pinfom=new TPaveText(0.6,0.7,1.,.87,"NDC");
1499     pinfob->SetBorderSize(0);
1500     pinfob->SetFillStyle(0);
1501     pinfom->SetBorderSize(0);
1502     pinfom->SetFillStyle(0);
1503     TF1* ff=fhistoInvMass->GetFunction("funcmass");
1504
1505     for (Int_t i=fNFinalPars-3;i<fNFinalPars;i++){
1506       pinfom->SetTextColor(kBlue);
1507       TString str=Form("%s = %f #pm %f",ff->GetParName(i),ff->GetParameter(i),ff->GetParError(i));
1508       if(!(writeFitInfo==1 && i==fNFinalPars-3)) pinfom->AddText(str);
1509     }
1510     pd->cd();
1511     pinfom->DrawClone();
1512
1513     TPaveText *pinfo2=new TPaveText(0.1,0.1,0.6,0.4,"NDC");
1514     pinfo2->SetBorderSize(0);
1515     pinfo2->SetFillStyle(0);
1516
1517     Double_t signif, signal, bkg, errsignif, errsignal, errbkg;
1518
1519     Significance(nsigma,signif,errsignif);
1520     Signal(nsigma,signal,errsignal);
1521     Background(nsigma,bkg, errbkg);
1522     /*
1523       Significance(1.828,1.892,signif,errsignif);
1524       Signal(1.828,1.892,signal,errsignal);
1525       Background(1.828,1.892,bkg, errbkg);
1526     */
1527     TString str=Form("Significance (%.0f#sigma) %.1f #pm %.1f ",nsigma,signif,errsignif);
1528     pinfo2->AddText(str);
1529     str=Form("S (%.0f#sigma) %.0f #pm %.0f ",nsigma,signal,errsignal);
1530     pinfo2->AddText(str);
1531     str=Form("B (%.0f#sigma) %.0f #pm %.0f",nsigma,bkg,errbkg);
1532     pinfo2->AddText(str);
1533
1534     pd->cd();
1535     pinfo2->Draw();
1536
1537     if(writeFitInfo > 1){
1538       for (Int_t i=0;i<fNFinalPars-3;i++){
1539         pinfob->SetTextColor(kRed);
1540         str=Form("%s = %f #pm %f",ff->GetParName(i),ff->GetParameter(i),ff->GetParError(i));
1541         pinfob->AddText(str);
1542       }
1543       pd->cd();
1544       pinfob->DrawClone();
1545     }
1546
1547
1548   }
1549   return;
1550 }
1551
1552 //_________________________________________________________________________
1553
1554 void AliHFMassFitter::DrawHere(TVirtualPad* pd,Double_t nsigma,Int_t writeFitInfo) const {
1555   //draws histogram together with fit functions with default nice colors in user canvas
1556   PlotFit(pd,nsigma,writeFitInfo);
1557
1558   pd->Draw();
1559
1560 }
1561 //_________________________________________________________________________
1562 void AliHFMassFitter::DrawFit(Double_t nsigma) const{
1563
1564   //draws histogram together with fit functions with default nice colors
1565   gStyle->SetOptStat(0);
1566   gStyle->SetCanvasColor(0);
1567   gStyle->SetFrameFillColor(0);
1568
1569
1570   TCanvas* c=(TCanvas*)GetPad(nsigma,1);
1571   c->Draw();
1572
1573 }
1574
1575
1576 //_________________________________________________________________________
1577
1578 void AliHFMassFitter::PrintParTitles() const{
1579
1580   //prints to screen the parameters names
1581
1582   TF1 *f=fhistoInvMass->GetFunction("funcmass");
1583   if(!f) {
1584     cout<<"Fit function not found"<<endl;
1585     return;
1586   }
1587
1588   Int_t np=0;
1589   switch (ftypeOfFit4Bkg){
1590   case 0: //expo
1591     np=2;
1592     break;
1593   case 1: //linear
1594     np=2;
1595     break;
1596   case 2: //pol2
1597     np=3;
1598     break;
1599   case 3: //no bkg
1600     np=1;
1601     break;
1602   }
1603
1604   np+=3; //3 parameter for signal
1605   if (ftypeOfFit4Sgn == 1) np+=3;
1606
1607   cout<<"Parameter Titles \n";
1608   for(Int_t i=0;i<np;i++){
1609     cout<<"Par "<<i<<": "<<f->GetParName(i)<<endl;
1610   }
1611   cout<<endl;
1612
1613 }
1614
1615 //************ significance
1616
1617 //_________________________________________________________________________
1618
1619 void AliHFMassFitter::Signal(Double_t nOfSigma,Double_t &signal,Double_t &errsignal) const {
1620   // Return signal integral in mean+- n sigma
1621
1622   if(fcounter==0) {
1623     cout<<"Use MassFitter method before Signal"<<endl;
1624     return;
1625   }
1626
1627   Double_t min=fMass-nOfSigma*fSigmaSgn;
1628   Double_t max=fMass+nOfSigma*fSigmaSgn;
1629
1630   Signal(min,max,signal,errsignal);
1631
1632 //   //functions names
1633 //   TString bkgname= "funcbkgRecalc";
1634 //   TString bkg1name="funcbkg1Recalc";
1635 //   TString massname="funcmass";
1636
1637
1638 //   TF1 *funcbkg=0;
1639 //   TF1 *funcmass=fhistoInvMass->GetFunction(massname.Data());
1640 //   if(!funcmass){
1641 //     cout<<"AliHFMassFitter::Signal() ERROR -> Mass distr function not found!"<<endl;
1642 //     return;
1643 //   }
1644
1645 //   if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
1646 //   else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
1647
1648 //   if(!funcbkg){
1649 //     cout<<"AliHFMassFitter::Signal() ERROR -> Bkg function not found!"<<endl;
1650 //     return;
1651 //   }
1652
1653 //   Int_t np=-99;
1654 //   switch (ftypeOfFit4Bkg){
1655 //   case 0: //expo
1656 //     np=2;
1657 //     break;
1658 //   case 1: //linear
1659 //     np=2;
1660 //     break;
1661 //   case 2: //pol2
1662 //     np=3;
1663 //     break;
1664 //   case 3: //no bkg
1665 //     np=1;
1666 //     break;
1667 //   }
1668
1669 //   Float_t intS,intSerr;
1670
1671 //  //relative error evaluation
1672
1673 //   intS=funcmass->GetParameter(np);
1674 //   intSerr=funcmass->GetParError(np);
1675
1676 //   cout<<"Sgn relative error evaluation from fit: "<<intSerr/intS<<endl;
1677 //   Double_t background,errbackground;
1678 //   Background(nOfSigma,background,errbackground);
1679
1680 //   //signal +/- error in nsigma
1681 //   Double_t min=fMass-nOfSigma*fSigmaSgn;
1682 //   Double_t max=fMass+nOfSigma*fSigmaSgn;
1683
1684 //   Double_t mass=funcmass->Integral(min, max)/fhistoInvMass->GetBinWidth(4);
1685 //   signal=mass - background;
1686 //   errsignal=TMath::Sqrt((intSerr/intS*mass)*(intSerr/intS*mass)/*assume relative error is the same as for total integral*/ + errbackground*errbackground);
1687   return;
1688
1689 }
1690
1691 //_________________________________________________________________________
1692
1693 void AliHFMassFitter::Signal(Double_t min, Double_t max, Double_t &signal,Double_t &errsignal) const {
1694
1695   // Return signal integral in a range
1696
1697   if(fcounter==0) {
1698     cout<<"Use MassFitter method before Signal"<<endl;
1699     return;
1700   }
1701
1702   //functions names
1703   TString bkgname="funcbkgRecalc";
1704   TString bkg1name="funcbkg1Recalc";
1705   TString massname="funcmass";
1706
1707
1708   TF1 *funcbkg=0;
1709   TF1 *funcmass=fhistoInvMass->GetFunction(massname.Data());
1710   if(!funcmass){
1711     cout<<"AliHFMassFitter::Signal() ERROR -> Mass distr function not found!"<<endl;
1712     return;
1713   }
1714
1715   if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
1716   else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
1717
1718   if(!funcbkg){
1719     cout<<"AliHFMassFitter::Signal() ERROR -> Bkg function not found!"<<endl;
1720     return;
1721   }
1722
1723   Int_t np=-99;
1724   switch (ftypeOfFit4Bkg){
1725   case 0: //expo
1726     np=2;
1727     break;
1728   case 1: //linear
1729     np=2;
1730     break;
1731   case 2: //pol2
1732     np=3;
1733     break;
1734   case 3: //no bkg
1735     np=1;
1736     break;
1737   }
1738
1739   Double_t intS,intSerr;
1740
1741  //relative error evaluation
1742   intS=funcmass->GetParameter(np);
1743   intSerr=funcmass->GetParError(np);
1744
1745   cout<<"Sgn relative error evaluation from fit: "<<intSerr/intS<<endl;
1746   Double_t background,errbackground;
1747   Background(min,max,background,errbackground);
1748
1749   //signal +/- error in the range
1750
1751   Double_t mass=funcmass->Integral(min, max)/fhistoInvMass->GetBinWidth(4);
1752   signal=mass - background;
1753   errsignal=(intSerr/intS)*signal;/*assume relative error is the same as for total integral*/
1754
1755 }
1756
1757 //_________________________________________________________________________
1758
1759 void AliHFMassFitter::Background(Double_t nOfSigma,Double_t &background,Double_t &errbackground) const {
1760   // Return background integral in mean+- n sigma
1761
1762   if(fcounter==0) {
1763     cout<<"Use MassFitter method before Background"<<endl;
1764     return;
1765   }
1766   Double_t min=fMass-nOfSigma*fSigmaSgn;
1767   Double_t max=fMass+nOfSigma*fSigmaSgn;
1768
1769   Background(min,max,background,errbackground);
1770
1771 //   //functions names
1772 //   TString bkgname="funcbkgRecalc";
1773 //   TString bkg1name="funcbkg1Recalc";
1774
1775 //   TF1 *funcbkg=0;
1776 //   if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
1777 //   else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
1778 //   if(!funcbkg){
1779 //     cout<<"AliHFMassFitter::Background() ERROR -> Bkg function not found!"<<endl;
1780 //     return;
1781 //   }
1782
1783 //   Float_t intB,intBerr;//, intT,intTerr,intS,intSerr;
1784
1785 //   //relative error evaluation: from final parameters of the fit
1786 //   if(ftypeOfFit4Bkg==3 && ftypeOfFit4Sgn == 0) cout<<"No background fit: Bkg relative error evaluation put to zero"<<endl;
1787 //   else{
1788
1789 //     intB=funcbkg->GetParameter(0);
1790 //     intBerr=funcbkg->GetParError(0);
1791
1792 //     cout<<"Bkg relative error evaluation: from final parameters of the fit: "<<intBerr/intB<<endl;
1793 //   }
1794
1795 //   Double_t min=fMass-nOfSigma*fSigmaSgn;
1796 //   Double_t max=fMass+nOfSigma*fSigmaSgn;
1797
1798 //   //relative error evaluation: from histo
1799
1800 //   intB=fhistoInvMass->Integral(1,fSideBandl)+fhistoInvMass->Integral(fSideBandr,fNbin);
1801 //   Double_t sum2=0;
1802
1803 //   for(Int_t i=1;i<=fSideBandl;i++){
1804 //     sum2+=fhistoInvMass->GetBinError(i)*fhistoInvMass->GetBinError(i);
1805 //   }
1806 //   for(Int_t i=fSideBandr;i<=fNbin;i++){
1807 //     sum2+=fhistoInvMass->GetBinError(i)*fhistoInvMass->GetBinError(i);
1808 //   }
1809
1810 //   intBerr=TMath::Sqrt(sum2);
1811 //   cout<<"Bkg relative error evaluation: from histo: "<<intBerr/intB<<endl;
1812
1813 //   cout<<"Last estimation of bkg error is used"<<endl;
1814
1815 //   //backround +/- error in nsigma
1816 //   if (ftypeOfFit4Bkg == 3 && ftypeOfFit4Sgn == 0) {
1817 //     background = 0;
1818 //     errbackground = 0;
1819 //   }
1820 //   else{
1821 //     background=funcbkg->Integral(min,max)/(Double_t)fhistoInvMass->GetBinWidth(2);
1822 //     errbackground=intBerr/intB*background; // assume relative error is the same as for total integral
1823 //     //cout<<"integral = "<<funcbkg->Integral(min, max)<<"\tbinW = "<<fhistoInvMass->GetBinWidth(2)<<endl;
1824 //   }
1825   return;
1826
1827 }
1828 //___________________________________________________________________________
1829
1830 void AliHFMassFitter::Background(Double_t min, Double_t max, Double_t &background,Double_t &errbackground) const {
1831   // Return background integral in a range
1832
1833   if(fcounter==0) {
1834     cout<<"Use MassFitter method before Background"<<endl;
1835     return;
1836   }
1837
1838   //functions names
1839   TString bkgname="funcbkgRecalc";
1840   TString bkg1name="funcbkg1Recalc";
1841
1842   TF1 *funcbkg=0;
1843   if(ftypeOfFit4Sgn == 0) funcbkg=fhistoInvMass->GetFunction(bkgname.Data());
1844   else funcbkg=fhistoInvMass->GetFunction(bkg1name.Data());
1845   if(!funcbkg){
1846     cout<<"AliHFMassFitter::Background() ERROR -> Bkg function not found!"<<endl;
1847     return;
1848   }
1849
1850
1851   Double_t intB,intBerr;
1852
1853   //relative error evaluation: from final parameters of the fit
1854   if(ftypeOfFit4Bkg==3 && ftypeOfFit4Sgn == 0) cout<<"No background fit: Bkg relative error evaluation put to zero"<<endl;
1855   else{
1856     intB=funcbkg->GetParameter(0);
1857     intBerr=funcbkg->GetParError(0);
1858     cout<<"Bkg relative error evaluation: from final parameters of the fit: "<<intBerr/intB<<endl;
1859   }
1860
1861   //relative error evaluation: from histo
1862
1863   intB=fhistoInvMass->Integral(1,fSideBandl)+fhistoInvMass->Integral(fSideBandr,fNbin);
1864   Double_t sum2=0;
1865   for(Int_t i=1;i<=fSideBandl;i++){
1866     sum2+=fhistoInvMass->GetBinError(i)*fhistoInvMass->GetBinError(i);
1867   }
1868   for(Int_t i=fSideBandr;i<=fNbin;i++){
1869     sum2+=fhistoInvMass->GetBinError(i)*fhistoInvMass->GetBinError(i);
1870   }
1871
1872   intBerr=TMath::Sqrt(sum2);
1873   cout<<"Bkg relative error evaluation: from histo: "<<intBerr/intB<<endl;
1874
1875   cout<<"Last estimation of bkg error is used"<<endl;
1876
1877   //backround +/- error in the range
1878   if (ftypeOfFit4Bkg == 3 && ftypeOfFit4Sgn == 0) {
1879     background = 0;
1880     errbackground = 0;
1881   }
1882   else{
1883     background=funcbkg->Integral(min,max)/(Double_t)fhistoInvMass->GetBinWidth(2);
1884     errbackground=intBerr/intB*background; // assume relative error is the same as for total integral
1885     //cout<<"integral = "<<funcbkg->Integral(min, max)<<"\tbinW = "<<fhistoInvMass->GetBinWidth(2)<<endl;
1886   }
1887   return;
1888
1889 }
1890
1891
1892 //__________________________________________________________________________
1893
1894 void AliHFMassFitter::Significance(Double_t nOfSigma,Double_t &significance,Double_t &errsignificance) const  {
1895   // Return significance in mean+- n sigma
1896
1897   Double_t signal,errsignal,background,errbackground;
1898   Signal(nOfSigma,signal,errsignal);
1899   Background(nOfSigma,background,errbackground);
1900
1901   significance =  signal/TMath::Sqrt(signal+background);
1902
1903   errsignificance = TMath::Sqrt(significance*significance/(signal+background)/(signal+background)*(1/4.*errsignal*errsignal+errbackground*errbackground)+significance*significance/signal/signal*errsignal*errsignal);
1904
1905   return;
1906 }
1907
1908 //__________________________________________________________________________
1909
1910 void AliHFMassFitter::Significance(Double_t min, Double_t max, Double_t &significance,Double_t &errsignificance) const {
1911   // Return significance integral in a range
1912
1913   Double_t signal,errsignal,background,errbackground;
1914   Signal(min, max,signal,errsignal);
1915   Background(min, max,background,errbackground);
1916
1917   significance =  signal/TMath::Sqrt(signal+background);
1918
1919   errsignificance = TMath::Sqrt(significance*significance/(signal+background)/(signal+background)*(1/4.*errsignal*errsignal+errbackground*errbackground)+significance*significance/signal/signal*errsignal*errsignal);
1920
1921   return;
1922 }
1923
1924