PWGHF/vertexingHF/AliHFMassFitter.cxx

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