PWGDQ/dielectron/AliDielectronBtoJPSItoEleCDFfitFCN.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 #include "AliLog.h"
  16 #include "AliDielectronBtoJPSItoEleCDFfitFCN.h"
  17 #include "TFormula.h"
  18 #include "TF1.h"
  19 #include "TCanvas.h"
  20 #include "TMath.h"
  21
  22 //_________________________________________________________________________
  23 //                        Class AliDielectronBtoJPSItoEleCDFfitFCN
  24 //                   Definition of main function used in
  25 //                     unbinned log-likelihood fit for
  26 //                 the channel B -> JPsi + X -> e+e- + X
  27 //
  28 //                           Origin: C.Di Giglio
  29 //       Contact: Carmelo.Digiglio@ba.infn.it , Giuseppe.Bruno@ba.infn.it
  30 //_________________________________________________________________________
  31
  32 ClassImp(AliDielectronBtoJPSItoEleCDFfitFCN)
  33
  34         //_________________________________________________________________________________________________
  35         AliDielectronBtoJPSItoEleCDFfitFCN::AliDielectronBtoJPSItoEleCDFfitFCN() :
  36                 fFPlus(0.),
  37                 fFMinus(0.),
  38                 fFSym(0.),
  39                 fintmMassSig(1.),
  40                 fintmMassBkg(1.),
  41                 fhCsiMC(0x0),
  42                 fMassWndHigh(0.),
  43                 fMassWndLow(0.),
  44                 fCrystalBallParam(kFALSE)
  45 {
  46         //
  47         // constructor
  48         //
  49         SetCrystalBallFunction(kFALSE);
  50         SetMassWndHigh(0.2);
  51         SetMassWndLow(0.5);
  52         fWeightType[0] = 1.; fWeightType[1] = 1.; fWeightType[2] = 1.;
  53         for(Int_t iPar = 0; iPar < 45; iPar++) fParameters[iPar] = 0.;
  54         fParameters[9] = 1.;fParameters[11] = 1.;fParameters[12] = 1.;
  55         AliInfo("Instance of AliDielectronBtoJPSItoEleCDFfitFCN-class created");
  56 }
  57 //_________________________________________________________________________________________________
  58 AliDielectronBtoJPSItoEleCDFfitFCN::AliDielectronBtoJPSItoEleCDFfitFCN(const AliDielectronBtoJPSItoEleCDFfitFCN& source) :
  59         TNamed(source),
  60         fFPlus(source.fFPlus),
  61         fFMinus(source.fFMinus),
  62         fFSym(source.fFSym),
  63         fintmMassSig(source.fintmMassSig),
  64         fintmMassBkg(source.fintmMassBkg),
  65         fhCsiMC(source.fhCsiMC),
  66         fMassWndHigh(source.fMassWndHigh),
  67         fMassWndLow(source.fMassWndLow),
  68         fCrystalBallParam(source.fCrystalBallParam)
  69 {
  70         //
  71         // Copy constructor
  72         //
  73         for(Int_t iPar = 0; iPar < 45; iPar++) fParameters[iPar] = source.fParameters[iPar];
  74         for(Int_t iW=0; iW<2; iW++) fWeightType[iW] = source.fWeightType[iW];
  75 }
  76 //_________________________________________________________________________________________________
  77 AliDielectronBtoJPSItoEleCDFfitFCN& AliDielectronBtoJPSItoEleCDFfitFCN::operator=(const AliDielectronBtoJPSItoEleCDFfitFCN& source)
  78 {
  79         //
  80         // Assignment operator
  81         //
  82         if(&source == this) return *this;
  83         fFPlus = source.fFPlus;
  84         fFMinus = source.fFMinus;
  85         fFSym = source.fFSym;
  86         fintmMassSig = source.fintmMassSig;
  87         fintmMassBkg = source.fintmMassBkg;
  88         fhCsiMC = source.fhCsiMC;
  89         fCrystalBallParam = source.fCrystalBallParam;
  90
  91         for(Int_t iPar = 0; iPar < 45; iPar++) fParameters[iPar] = source.fParameters[iPar];
  92         return *this;
  93 }
  94 //_________________________________________________________________________________________________
  95 AliDielectronBtoJPSItoEleCDFfitFCN::~AliDielectronBtoJPSItoEleCDFfitFCN()
  96 {
  97         //
  98         // Default destructor
  99         //
 100
 101         delete fhCsiMC;
 102         for(Int_t iPar = 0; iPar < 45; iPar++) fParameters[iPar] = 0.;
 103 }
 104 //_________________________________________________________________________________________________
 105 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateLikelihood(const Double_t* pseudoproperdecaytime,
 106                 const Double_t* invariantmass, const Int_t *type, const Int_t ncand) const
 107 {
 108         //
 109         // This function evaluates the Likelihood fnction
 110         // It returns the -Log(of the likelihood function)
 111         //
 112         Double_t f = 0.;
 113         Double_t ret = 0.;
 114
 115         for(Int_t i=0; i < ncand; i++) {
 116                 f = EvaluateCDFfuncNorm(pseudoproperdecaytime[i],invariantmass[i],type[i]);
 117                 if(f <= 0.) continue;
 118                 ret += -1*TMath::Log(f);
 119         }
 120         return ret;
 121 }
 122 //_________________________________________________________________________________________________
 123 void AliDielectronBtoJPSItoEleCDFfitFCN::SetAllParameters(const Double_t* parameters)
 124 {
 125         //
 126         // Sets array of FCN parameters
 127         //
 128         for(Int_t index = 0; index < 45; index++) fParameters[index] = parameters[index];
 129 }
 130 //_________________________________________________________________________________________________
 131 void AliDielectronBtoJPSItoEleCDFfitFCN::ComputeMassIntegral()
 132 {
 133         //
 134         // this function compute the integral of the likelihood function
 135         // (theoretical function) in order to normalize it to unity
 136         //
 137         Double_t npm = 20000.;
 138         Double_t stepm;
 139         Double_t mx=0.;
 140         stepm = (fMassWndHigh-fMassWndLow)/npm;
 141         // compute integrals for  invariant mass terms
 142
 143         Double_t iMassSig;
 144         Double_t intmMassSig = 0.0;
 145         Double_t summMassSig = 0.0;
 146         for(iMassSig = 1.0;  iMassSig<= npm/2.; iMassSig++)  {
 147                 mx = fMassWndLow + (iMassSig - .5)*stepm;
 148                 summMassSig += EvaluateCDFInvMassSigDistr(mx);
 149                 mx = fMassWndHigh - (iMassSig - .5)*stepm;
 150                 summMassSig += EvaluateCDFInvMassSigDistr(mx);
 151         }
 152         intmMassSig = summMassSig*stepm;
 153         SetIntegralMassSig(intmMassSig);
 154         //
 155         Double_t iMassBkg;
 156         Double_t intmMassBkg = 0.0;
 157         Double_t summMassBkg = 0.0;
 158         for(iMassBkg = 1.0; iMassBkg <= npm/2.; iMassBkg++)  {
 159                 mx = fMassWndLow + (iMassBkg - .5)*stepm;
 160                 summMassBkg += EvaluateCDFInvMassBkgDistr(mx);
 161                 mx = fMassWndHigh - (iMassBkg - .5)*stepm;
 162                 summMassBkg += EvaluateCDFInvMassBkgDistr(mx);
 163         }
 164         intmMassBkg = summMassBkg*stepm;
 165         if(intmMassBkg < 1.e-05) intmMassBkg = 1.;
 166         SetIntegralMassBkg(intmMassBkg);
 167 }
 168 //_________________________________________________________________________________________________
 169 void AliDielectronBtoJPSItoEleCDFfitFCN::PrintStatus()
 170 {
 171         //
 172         //  Print the parameters of the fits
 173         //
 174         printf("\n");
 175         // background param
 176         printf("actual value of fWeightRes------------------------------------->> | %f \n", GetResWeight());
 177         printf("actual value of fPos ------------------------------------------>> | %f \n", GetFPlus());
 178         printf("actual value of fNeg ------------------------------------------>> | %f \n", GetFMinus());
 179         printf("actual value of fSym ------------------------------------------>> | %f \n", GetFSym());
 180         printf("actual value of fOneOvLamPlus --------------------------------->> | %f \n", GetLamPlus());
 181         printf("actual value of fOneOvLamMinus -------------------------------->> | %f \n", GetLamMinus());
 182         printf("actual value of fOneOvLamSym ---------------------------------->> | %f \n", GetLamSym());
 183         printf("actual value of fFractionJpsiFromBeauty ----------------------->> | %f \n", GetFractionJpsiFromBeauty());
 184         printf("actual value of fFsig ----------------------------------------->> | %f \n", GetFsig());
 185
 186         if(fCrystalBallParam){
 187                 printf("actual value of fCrystalBallMmean ----------------------------->> | %f \n", GetCrystalBallMmean());
 188                 printf("actual value of fCrystalBallNexp ------------------------------>> | %f \n", GetCrystalBallNexp());
 189                 printf("actual value of fCrystalBallSigma ----------------------------->> | %f \n", GetCrystalBallSigma());
 190                 printf("actual value of fCrystalBallAlpha ----------------------------->> | %f \n", GetCrystalBallAlpha());
 191                 printf("actual value of fCrystalBallNorm  ----------------------------->> | %f \n", GetCrystalBallNorm());
 192         }else{
 193                 printf("actual value of fMpv ------------------------------------------>> | %f \n", GetCrystalBallMmean());
 194                 printf("actual value of fConstRovL ------------------------------------>> | %f \n", GetCrystalBallNexp());
 195                 printf("actual value of fSigmaL --------------------------------------->> | %f \n", GetCrystalBallSigma());
 196                 printf("actual value of fSigmaR --------------------------------------->> | %f \n", GetCrystalBallAlpha());
 197         }
 198
 199         // back Mass func
 200         printf("actual value of normBkg ----------------------------------------->> | %f \n", GetBkgInvMassNorm());
 201         printf("actual value of meanBkg ----------------------------------------->> | %f \n", GetBkgInvMassMean());
 202         printf("actual value of slopeBkg ---------------------------------------->> | %f \n", GetBkgInvMassSlope());
 203         printf("actual value of constBkg ---------------------------------------->> | %f \n", GetBkgInvMassConst());
 204         // resolution func (First-First)
 205         printf("actual value of norm1Gauss (FF) --------------------------------->> | %f \n", GetNormGaus1ResFunc(2));
 206         printf("actual value of norm2Gauss (FF) --------------------------------->> | %f \n", GetNormGaus2ResFunc(2));
 207         printf("actual value of fMean1Res (FF) ---------------------------------->> | %f \n", GetResMean1(2));
 208         printf("actual value of fSigma1Res (FF) --------------------------------->> | %f \n", GetResSigma1(2));
 209         printf("actual value of fMean2Res (FF) ---------------------------------->> | %f \n", GetResMean2(2));
 210         printf("actual value of fSigma2Res (FF) --------------------------------->> | %f \n", GetResSigma2(2));
 211
 212         printf("actual value of alfaRes (FF) ------------------------------------>> | %f \n", GetResAlfa(2));
 213         printf("actual value of lambdaRes (FF) ---------------------------------->> | %f \n", GetResLambda(2));
 214         printf("actual value of normExpRes (FF) --------------------------------->> | %f \n", GetResNormExp(2));
 215
 216         // resolution func (First-Second)
 217         printf("actual value of norm1Gauss (FS) --------------------------------->> | %f \n", GetNormGaus1ResFunc(1));
 218         printf("actual value of norm2Gauss (FS) --------------------------------->> | %f \n", GetNormGaus2ResFunc(1));
 219         printf("actual value of fMean1Res (FS) ---------------------------------->> | %f \n", GetResMean1(1));
 220         printf("actual value of fSigma1Res (FS) --------------------------------->> | %f \n", GetResSigma1(1));
 221         printf("actual value of fMean2Res (FS) ---------------------------------->> | %f \n", GetResMean2(1));
 222         printf("actual value of fSigma2Res (FS) --------------------------------->> | %f \n", GetResSigma2(1));
 223
 224         printf("actual value of alfaRes (FS) ------------------------------------>> | %f \n", GetResAlfa(1));
 225         printf("actual value of lambdaRes (FS) ---------------------------------->> | %f \n", GetResLambda(1));
 226         printf("actual value of normExpRes (FS) --------------------------------->> | %f \n", GetResNormExp(1));
 227
 228         // resolution func (Second-Second)
 229         printf("actual value of norm1Gauss (SS) --------------------------------->> | %f \n", GetNormGaus1ResFunc(0));
 230         printf("actual value of norm2Gauss (SS) --------------------------------->> | %f \n", GetNormGaus2ResFunc(0));
 231         printf("actual value of fMean1Res (SS) ---------------------------------->> | %f \n", GetResMean1(0));
 232         printf("actual value of fSigma1Res (SS) --------------------------------->> | %f \n", GetResSigma1(0));
 233         printf("actual value of fMean2Res (SS) ---------------------------------->> | %f \n", GetResMean2(0));
 234         printf("actual value of fSigma2Res (SS) --------------------------------->> | %f \n", GetResSigma2(0));
 235
 236         printf("actual value of alfaRes (SS) ------------------------------------>> | %f \n", GetResAlfa(0));
 237         printf("actual value of lambdaRes (SS) ---------------------------------->> | %f \n", GetResLambda(0));
 238         printf("actual value of normExpRes (SS) --------------------------------->> | %f \n", GetResNormExp(0));
 239
 240         printf("\n");
 241         // integrals constants
 242         printf("Actual value of normalization integral for MassSig ---------------->> | %f \n", GetIntegralMassSig());
 243         printf("Actual value of normalization integral for MassBkg ---------------->> | %f \n", GetIntegralMassBkg());
 244
 245         printf("\n");
 246 }
 247 //_________________________________________________________________________________________________
 248 void AliDielectronBtoJPSItoEleCDFfitFCN::SetResolutionConstants(const Double_t* resolutionConst, Int_t type)
 249 {
 250         //
 251         // Resolution function is parametrized as the sum of two gaussian
 252         //
 253         Int_t index = (2-type)*9;
 254         fParameters[20+index]=resolutionConst[1]; //mean1
 255         fParameters[22+index]=resolutionConst[4]; //mean2
 256         fParameters[18+index]=resolutionConst[0]; //norm1
 257         fParameters[21+index]=resolutionConst[2]; //sigma1
 258         fParameters[23+index]=resolutionConst[5]; //sigma2
 259         fParameters[19+index]=resolutionConst[3]; //norm2
 260
 261         //exp values
 262         fParameters[24+index]=resolutionConst[6]; //alfa
 263         fParameters[25+index]=resolutionConst[7]; //lambda
 264         fParameters[26+index]=resolutionConst[8]; //norm3
 265         return;
 266 }
 267
 268 //________________________________________________________________________________________________
 269 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFfunc(Double_t x, Double_t m, Int_t type) const
 270 {
 271         // evaluate likelihood function
 272         return fParameters[8]*EvaluateCDFfuncSignalPart(x,m,type) + (1. - fParameters[8])*EvaluateCDFfuncBkgPart(x,m,type);
 273 }
 274
 275 //_________________________________________________________________________________________________
 276 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFfuncNorm(Double_t x, Double_t m, Int_t type) const
 277 {
 278         // evaluate likelihood function
 279         return EvaluateCDFfunc(x,m,type);
 280 }
 281
 282 //_________________________________________________________________________________________________
 283 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFfuncSignalPart(Double_t x, Double_t m, Int_t type) const
 284 {
 285   // evaluate psproper signal
 286   return EvaluateCDFDecayTimeSigDistr(x,type)*(EvaluateCDFInvMassSigDistr(m)/fintmMassSig);
 287 }
 288
 289 //_________________________________________________________________________________________________
 290 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeSigDistr(Double_t x, Int_t type) const
 291 {
 292         //
 293         // Implementation of the Background part of the Likelyhood function
 294         //
 295         Double_t retvalue = 0.;
 296         Double_t funBnorm = FunB(x,type);
 297         Double_t funPnorm = ResolutionFunc(x,type);
 298         retvalue = fParameters[7]*funBnorm + (1. - fParameters[7])*funPnorm;
 299         return retvalue;
 300 }
 301
 302 //_________________________________________________________________________________________________
 303 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFInvMassSigDistr(Double_t m) const
 304 {
 305         //
 306         // Parametrization of signal part invariant mass distribution
 307         // It can be either Crystal Ball function or sum of two Landau
 308         //
 309         Double_t fitval = 0.;
 310
 311         if(fCrystalBallParam){
 312                 Double_t t = (m-fParameters[9])/fParameters[11]; ;
 313                 if (fParameters[12] < 0) t = -t;
 314                 Double_t absAlpha = TMath::Abs((Double_t)fParameters[12]);
 315
 316                 if (t >= -absAlpha) {
 317                         return fParameters[13]*TMath::Exp(-0.5*t*t);
 318                 }
 319                 else {
 320                         Double_t a =  TMath::Power(fParameters[10]/absAlpha,fParameters[10])* TMath::Exp(-0.5*absAlpha*absAlpha);
 321                         Double_t b= fParameters[10]/absAlpha - absAlpha;
 322                         fitval = (fParameters[13]*a/TMath::Power(b - t, fParameters[10]));
 323                         return fitval;
 324                 }
 325         }else{
 326                 Double_t t=-1*m;
 327                 Double_t tmpv=-1*fParameters[9];
 328                 fitval=TMath::Sqrt(TMath::Landau(t,tmpv,fParameters[11]));
 329                 fitval += fParameters[10]*(TMath::Landau(m,fParameters[9],fParameters[12]));
 330                 return fitval;
 331         }
 332 }
 333 //_________________________________________________________________________________________________
 334 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::FunB(Double_t x, Int_t type) const
 335 {
 336         //
 337         // Parameterisation of the fit function for the x part of the Background
 338         //
 339         Double_t np = 1000.0;
 340         Double_t sc = 10.;
 341         Double_t sigma3 = 1000.; // valore usato nella macro
 342         Double_t xprime;
 343         Double_t sum = 0.0;
 344         Double_t xlow,xupp;
 345         Double_t step;
 346         Double_t i;
 347         xlow = x - sc * sigma3 ;
 348         xupp = x + sc * sigma3 ;
 349         step = (xupp-xlow) / np;
 350         Double_t csiMCxprime = 0.;
 351         Double_t resolutionxdiff = 0.;
 352         Double_t xdiff = 0.;
 353
 354         for(i=1.0; i<=np; i++){
 355                 xprime = xlow + (i-.5) * step;
 356                 csiMCxprime = CsiMC(xprime);
 357                 xdiff = xprime - x;
 358                 resolutionxdiff = ResolutionFunc(xdiff, type); // normalized value
 359                 sum += csiMCxprime * resolutionxdiff;
 360         }
 361
 362         return step * sum ;
 363 }
 364 //_________________________________________________________________________________________________
 365 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::FunP(Double_t x, Int_t type) const
 366 {
 367         //
 368         //  Parameterisation of the Prompt part for the x distribution
 369         //
 370         return ResolutionFunc(x,type);
 371 }
 372
 373
 374 //_________________________________________________________________________________________________
 375 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::CsiMC(Double_t x) const
 376 {
 377         //
 378         //  Distribution (template) of the x distribution for the x variable
 379         //  for the J/psi coming from Beauty hadrons
 380         //
 381         Double_t returnvalue = 0.;
 382
 383         if((fhCsiMC->FindBin(x) > 0) && (fhCsiMC->FindBin(x) < fhCsiMC->GetNbinsX()+1))
 384         returnvalue = fhCsiMC->Interpolate(x);
 385
 386         return returnvalue;
 387 }
 388
 389 //_________________________________________________________________________________________________
 390 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFfuncBkgPart(Double_t x,Double_t m, Int_t type) const
 391 {
 392         //
 393         // Return the part of the likelihood function for the background hypothesis
 394         //
 395         return EvaluateCDFDecayTimeBkgDistr(x,type)*(EvaluateCDFInvMassBkgDistr(m)/fintmMassBkg);
 396 }
 397
 398 //_________________________________________________________________________________________________
 399 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeBkgDistr(Double_t x, Int_t type) const
 400 {
 401         //
 402         // it returns the value of the probability to have a given x for the background
 403         //
 404
 405         Double_t ret = fParameters[0]/(fParameters[0]+fParameters[1]+fParameters[2]+fParameters[3])*ResolutionFunc(x,type) + fParameters[1]/(fParameters[0]+fParameters[1]+fParameters[2]+fParameters[3])*FunBkgPos(x,type) +  fParameters[2]/(fParameters[0]+fParameters[1]+fParameters[2]+fParameters[3])*FunBkgNeg(x,type) + fParameters[3]/(fParameters[0]+fParameters[1]+fParameters[2]+fParameters[3])*FunBkgSym(x,type);
 406         return ret;
 407 }
 408
 409 //_________________________________________________________________________________________________
 410 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFInvMassBkgDistr(Double_t m) const
 411 {
 412         //
 413         // it returns the value of the probability to have a given mass for the background
 414         //
 415         Double_t value = 0.;
 416         value = fParameters[14]*TMath::Exp(-1*(m-fParameters[15])/fParameters[16]) + fParameters[17];
 417         return value;
 418 }
 419 //_________________________________________________________________________________________________
 420 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::FunBkgPos(Double_t x,Int_t type) const
 421 {
 422         //
 423         // exponential with positive slopes for the background part (x)
 424         //
 425         Double_t np = 1000.0;
 426         Double_t sc = 10.;
 427         Double_t sigma3 = 1000.; // valore usato nella macro
 428         Double_t xprime;
 429         Double_t sum = 0.0;
 430         Double_t xlow,xupp;
 431         Double_t step;
 432         Double_t i;
 433         xlow = x - sc * sigma3 ;
 434         xupp = x + sc * sigma3 ;
 435         step = (xupp-xlow) / np;
 436
 437         for(i=1.0; i<=np/2; i++) {
 438                 xprime = xlow + (i-.5) * step;
 439                 if (xprime > 0) {sum += fParameters[4] * TMath::Exp(-1*xprime*fParameters[4])*(ResolutionFunc(xprime-x,type));}
 440                 xprime = xupp - (i-.5) * step;
 441                 if (xprime > 0) {sum += fParameters[4] * TMath::Exp(-1*xprime*fParameters[4])*(ResolutionFunc(xprime-x,type));}
 442         }
 443
 444         return step * sum ;
 445 }
 446 //_________________________________________________________________________________________________
 447 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::FunBkgNeg(Double_t x, Int_t type) const
 448 {
 449         //
 450         // exponential with negative slopes for the background part (x)
 451         //
 452         Double_t np = 1000.0;
 453         Double_t sc = 10.;
 454         Double_t sigma3 = 1000.;
 455         Double_t xprime;
 456         Double_t sum = 0.0;
 457         Double_t xlow,xupp;
 458         Double_t step;
 459         Double_t i;
 460         xlow = x - sc * sigma3 ;
 461         xupp = x + sc * sigma3 ;
 462         step = (xupp-xlow) / np;
 463
 464         for(i=1.0; i<=np/2; i++) {
 465
 466                 xprime = xlow + (i-.5) * step;
 467                 if (xprime < 0) {sum += fParameters[5] * TMath::Exp(xprime*fParameters[5]) * (ResolutionFunc(xprime-x,type));}
 468
 469                 xprime = xupp - (i-.5) * step;
 470                 if (xprime < 0) {sum += fParameters[5] * TMath::Exp(xprime*fParameters[5]) * (ResolutionFunc(xprime-x,type));}
 471         }
 472
 473         return step * sum ;
 474 }
 475 //_________________________________________________________________________________________________
 476 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::FunBkgSym(Double_t x, Int_t type) const
 477 {
 478         //
 479         // exponential with both positive and negative slopes for the background part (x)
 480         //
 481         Double_t np = 1000.0;
 482         Double_t sc = 10.;
 483         Double_t sigma3 = 1000.;
 484         Double_t xprime;
 485         Double_t sum1 = 0.0;
 486         Double_t sum2 = 0.0;
 487         Double_t xlow,xupp;
 488         Double_t step;
 489         Double_t i;
 490         xlow = x - sc * sigma3 ;
 491         xupp = x + sc * sigma3 ;
 492         step = (xupp-xlow) / np;
 493
 494         for(i=1.0; i<=np/2; i++) {
 495
 496                 xprime = xlow + (i-.5) * step;
 497                 if (xprime > 0) {sum1 += 0.5 * fParameters[6]*TMath::Exp(-1*xprime*fParameters[6]) * (ResolutionFunc(xprime-x,type));}
 498                 if (xprime < 0) {sum2 += 0.5 * fParameters[6]*TMath::Exp(xprime*fParameters[6]) * (ResolutionFunc(xprime-x,type));}
 499
 500                 xprime = xupp - (i-.5) * step;
 501                 if (xprime > 0) {sum1 += 0.5 * fParameters[6]*TMath::Exp(-1*xprime*fParameters[6]) * (ResolutionFunc(xprime-x,type));}
 502                 if (xprime < 0) {sum2 += 0.5 * fParameters[6]*TMath::Exp(xprime*fParameters[6]) * (ResolutionFunc(xprime-x,type));}
 503         }
 504
 505         return step*(sum1 + sum2) ;
 506 }
 507 //_________________________________________________________________________________________________
 508 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::ResolutionFunc(Double_t x, Int_t type) const
 509 {
 510         //
 511         // parametrization with 2 gaus
 512         //
 513         Int_t index = (2-type)*9;
 514         Double_t mean1 = fParameters[20+index];
 515         Double_t mean2 = fParameters[22+index];
 516         Double_t norm1 = fParameters[18+index];
 517         Double_t sigma1 = fParameters[21+index];
 518         Double_t sigma2 = fParameters[23+index];
 519         Double_t norm2 = fParameters[19+index];
 520         //exp values
 521         Double_t alfa = fParameters[24+index];
 522         Double_t lambda = fParameters[25+index];
 523         Double_t norm3 = fParameters[26+index];
 524
 525         Double_t ret = 0.; Double_t fitval = 0;
 526         if(TMath::Abs(x)<=alfa) fitval = (lambda-1)/(2*alfa*lambda);
 527         else  fitval = ((lambda-1)/(2*alfa*lambda))*TMath::Power(alfa,lambda)*(TMath::Power(TMath::Abs(x),-1*lambda));
 528
 529         ret = (norm1/(norm1+norm2+norm3))*((1/(sigma1*TMath::Sqrt(2*TMath::Pi())))*TMath::Exp(-0.5*((x-mean1)/sigma1)*((x-mean1)/sigma1))) + (norm2/(norm1+norm2+norm3))*((1/(sigma2*TMath::Sqrt(2*TMath::Pi())))*TMath::Exp(-0.5*((x-mean2)/sigma2)*((x-mean2)/sigma2))) + (norm3/(norm1+norm2+norm3))*fitval;
 530
 531         return ret;
 532 }
 533
 534 //_________________________________________________________________________________________________
 535 TF1* AliDielectronBtoJPSItoEleCDFfitFCN::GetCsiMC(Double_t xmin, Double_t xmax, Double_t normalization)
 536 {
 537         // return the pointer to the templateMC function
 538         TF1* templateMC = new TF1("MCtemplate",this,&AliDielectronBtoJPSItoEleCDFfitFCN::CsiMCfunc,xmin,xmax,1);
 539         templateMC->SetParameter(0,normalization);
 540         templateMC->SetNpx(5000);
 541         return (TF1*)templateMC->Clone();
 542 }
 543
 544 //__________________________________________________________________________________________________
 545 TF1* AliDielectronBtoJPSItoEleCDFfitFCN::GetResolutionFunc(Double_t xmin, Double_t xmax, Double_t normalization, Double_t type){
 546         // return the pointer to the resolution function
 547         TF1* resFunc = new TF1(Form("resolutionFunc_%f",type),this,&AliDielectronBtoJPSItoEleCDFfitFCN::ResolutionFuncf,xmin,xmax,2);
 548         resFunc->SetParameter(0,normalization);
 549         resFunc->SetParameter(1,type);
 550         resFunc->SetNpx(5000);
 551         return (TF1*)resFunc->Clone();
 552 }
 553
 554 //__________________________________________________________________________________________________
 555 TF1* AliDielectronBtoJPSItoEleCDFfitFCN::GetResolutionFuncAllTypes(Double_t xmin, Double_t xmax, Double_t normalization){
 556         // return the pointer to the resolution function
 557         TF1* resFunc = new TF1(Form("resolutionFunc"),this,&AliDielectronBtoJPSItoEleCDFfitFCN::ResolutionFuncAllTypes,xmin,xmax,1);
 558         resFunc->SetParameter(0,normalization);
 559         resFunc->SetNpx(5000);
 560         return (TF1*)resFunc->Clone();
 561 }
 562
 563 //___________________________________________________________________________________________________
 564 TF1* AliDielectronBtoJPSItoEleCDFfitFCN::GetEvaluateCDFDecayTimeBkgDistr(Double_t xmin, Double_t xmax, Double_t normalization, Double_t type){
 565         // return the pointer to the background x distribution function
 566         TF1 *backFunc = new TF1(Form("backFunc_%f",type),this,&AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeBkgDistrFunc,xmin,xmax,2);
 567         backFunc->SetParameter(0,normalization);
 568         backFunc->SetParameter(1,type);
 569         backFunc->SetNpx(5000);
 570         return (TF1*)backFunc->Clone();
 571 }
 572
 573 //___________________________________________________________________________________________________
 574 TF1* AliDielectronBtoJPSItoEleCDFfitFCN::GetEvaluateCDFDecayTimeBkgDistrAllTypes(Double_t xmin, Double_t xmax, Double_t normalization){
 575         // return the pointer to the background x distribution function
 576         TF1 *backFunc = new TF1(Form("backFunc"),this,&AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeBkgDistrFuncAllTypes,xmin,xmax,1);
 577         backFunc->SetParameter(0,normalization);
 578         backFunc->SetNpx(5000);
 579         return (TF1*)backFunc->Clone();
 580 }
 581
 582 //__________________________________________________________________________________________________
 583 TF1* AliDielectronBtoJPSItoEleCDFfitFCN::GetEvaluateCDFDecayTimeSigDistr(Double_t xmin, Double_t xmax, Double_t normalization, Double_t type){
 584         // return the pointer to the signal x distribution function
 585         TF1 *signFunc = new TF1("signalFunc",this,&AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeSigDistrFunc,xmin,xmax,2);
 586         signFunc->SetParameter(0,normalization);
 587         signFunc->SetParameter(1,type);
 588         signFunc->SetNpx(5000);
 589         return (TF1*)signFunc->Clone();
 590 }
 591
 592 //____________________________________________________________________________________________________
 593 TF1* AliDielectronBtoJPSItoEleCDFfitFCN::GetEvaluateCDFInvMassBkgDistr(Double_t mMin, Double_t mMax, Double_t normalization){
 594   // return the pointer to the invariant mass distribution for the background
 595   TF1 * invMassBkg = new TF1("invMassBkg",this,&AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFInvMassBkgDistrFunc,mMin,mMax,1);
 596   invMassBkg->SetParameter(0,normalization);
 597   invMassBkg->SetNpx(5000);
 598   return (TF1*)invMassBkg->Clone();
 599 }
 600
 601
 602 //____________________________________________________________________________________________________
 603 TF1* AliDielectronBtoJPSItoEleCDFfitFCN::GetEvaluateCDFInvMassSigDistr(Double_t mMin, Double_t mMax, Double_t normalization){
 604   // return the pointer to the invariant mass distribution for the signal
 605   TF1 * invMassSig = new TF1("invMassSig",this,&AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFInvMassSigDistrFunc,mMin,mMax,1);
 606   invMassSig->SetParameter(0,normalization);
 607   invMassSig->SetNpx(5000);
 608   return (TF1*)invMassSig->Clone();
 609 }
 610
 611 //____________________________________________________________________________________________________
 612 TF1 *AliDielectronBtoJPSItoEleCDFfitFCN::GetEvaluateCDFInvMassTotalDistr(Double_t mMin, Double_t mMax, Double_t normalization){
 613   // return the pointer to the invariant mass distribution
 614   TF1 * invMassTot = new TF1("invMassTot",this,&AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFInvMassTotalDistr,mMin,mMax,1);
 615   invMassTot->SetParameter(0,normalization);
 616   invMassTot->SetNpx(5000);
 617   return (TF1*)invMassTot->Clone();
 618 }
 619
 620 //____________________________________________________________________________________________________
 621 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFInvMassTotalDistr(const Double_t* x, const Double_t *par) const
 622 {
 623   // evaluate invariant mass total distribution
 624   Double_t value = 0;
 625   Double_t xx = x[0];
 626   value = ((EvaluateCDFInvMassSigDistr(xx)/fintmMassSig)*fParameters[8] + (1-fParameters[8])*(EvaluateCDFInvMassBkgDistr(xx)/fintmMassBkg))*par[0];
 627   return value;
 628 }
 629
 630 //____________________________________________________________________________________________________
 631 TF1 *AliDielectronBtoJPSItoEleCDFfitFCN::GetEvaluateCDFDecayTimeTotalDistr(Double_t xMin, Double_t xMax, Double_t normalization,Double_t type){
 632  // return the pointer to the pseudoproper distribution for the background
 633  TF1 *decayTimeTot = new TF1(Form("decayTimeTot_%f",type),this,&AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeTotalDistr,xMin,xMax,2);
 634  decayTimeTot->SetParameter(0,normalization);
 635  decayTimeTot->SetParameter(1,type);
 636  decayTimeTot->SetNpx(5000);
 637  return (TF1*)decayTimeTot->Clone();
 638 }
 639
 640 //____________________________________________________________________________________________________
 641 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeTotalDistr(const Double_t* x, const Double_t *par) const
 642 {
 643  // evaluate the total pseudoproper distribution for a given candidate's type. par[1] should be the candidate's type.
 644  Double_t value = 0;
 645  Double_t xx = x[0];
 646  value = (fParameters[8]*EvaluateCDFDecayTimeSigDistr(xx,(Int_t)par[1]) + (1-fParameters[8])*EvaluateCDFDecayTimeBkgDistr(xx,(Int_t)par[1]))*par[0];
 647  return value;
 648 }
 649
 650 //____________________________________________________________________________________________________
 651 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeTotalDistrAllTypes(const Double_t* x, const Double_t *par) const
 652 {
 653  // evaluate the total pseudoproper distribution considering all candidate's types
 654  Double_t value = 0;
 655  Double_t xx = x[0];
 656
 657  value = (fParameters[8]*(fWeightType[2]*EvaluateCDFDecayTimeSigDistr(xx,2)+fWeightType[1]*EvaluateCDFDecayTimeSigDistr(xx,1)+fWeightType[0]*EvaluateCDFDecayTimeSigDistr(xx,0)))+((1-fParameters[8])*(fWeightType[2]*EvaluateCDFDecayTimeBkgDistr(xx,2) + fWeightType[1]*EvaluateCDFDecayTimeBkgDistr(xx,1)+fWeightType[0]*EvaluateCDFDecayTimeBkgDistr(xx,0)));
 658
 659  return value*par[0];
 660 }
 661
 662 //____________________________________________________________________________________________________
 663 TF1 *AliDielectronBtoJPSItoEleCDFfitFCN::GetEvaluateCDFDecayTimeTotalDistrAllTypes(Double_t xMin, Double_t xMax, Double_t normalization){
 664  // return the pointer to the pseudoproper function for the background considering all candidate's types
 665  TF1 *decayTimeTot = new TF1(Form("decayTimeTot"),this,&AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeTotalDistrAllTypes,xMin,xMax,1);
 666  decayTimeTot->SetParameter(0,normalization);
 667  decayTimeTot->SetNpx(5000);
 668  return (TF1*)decayTimeTot->Clone();
 669 }
 670
 671
 672 //____________________________________________________________________________________________________
 673 TF1 * AliDielectronBtoJPSItoEleCDFfitFCN::GetFunB(Double_t xmin, Double_t xmax, Double_t normalization, Double_t type){
 674  // return the pointer to the function that describe secondary jpsi pseudoproper distribution for a given candidate's type
 675  TF1* funb = new TF1(Form("secondaryJpsiConvolution_%f",type),this,&AliDielectronBtoJPSItoEleCDFfitFCN::FunBfunc,xmin,xmax,2);
 676  funb->SetParameter(0,normalization);
 677  funb->SetParameter(1,type);
 678  funb->SetNpx(5000);
 679  return (TF1*)funb->Clone();
 680  }
 681
 682 //____________________________________________________________________________________________________
 683 TF1 * AliDielectronBtoJPSItoEleCDFfitFCN::GetFunBAllTypes(Double_t xmin, Double_t xmax, Double_t normalization){
 684 // return the pointer to the function that describe secondary jpsi pseudoproper distribution for all types
 685  TF1* funb = new TF1(Form("secondaryJpsiConvolution"),this,&AliDielectronBtoJPSItoEleCDFfitFCN::FunBfuncAllTypes,xmin,xmax,1);
 686  funb->SetParameter(0,normalization);
 687  funb->SetNpx(5000);
 688  return (TF1*)funb->Clone();
 689  }
 690
 691