PWG3/vertexingHF/AliBtoJPSItoEleCDFfitFCN.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 "AliBtoJPSItoEleCDFfitFCN.h"
  17
  18 //_________________________________________________________________________
  19 //                        Class AliBtoJPSItoEleCDFfitFCN
  20 //                   Definition of main function used in
  21 //                     unbinned log-likelihood fit for
  22 //                 the channel B -> JPsi + X -> e+e- + X
  23 //
  24 //                           Origin: C.Di Giglio
  25 //       Contact: Carmelo.Digiglio@ba.infn.it , Giuseppe.Bruno@ba.infn.it
  26 //_________________________________________________________________________
  27
  28 ClassImp(AliBtoJPSItoEleCDFfitFCN)
  29
  30 //_________________________________________________________________________________________________
  31 AliBtoJPSItoEleCDFfitFCN::AliBtoJPSItoEleCDFfitFCN() :
  32 fFPlus(0.),
  33 fFMinus(0.),
  34 fFSym(0.),
  35 fIntegral(1.),
  36 fintxFunB(1.),
  37 fintxDecayTimeBkgPos(1.),
  38 fintxDecayTimeBkgNeg(1.),
  39 fintxDecayTimeBkgSym(1.),
  40 fintmMassSig(1.),
  41 fintxRes(1.),
  42 fintmMassBkg(1.),
  43 fhCsiMC(0x0),
  44 fMassWndHigh(0.),
  45 fMassWndLow(0.),
  46 fCrystalBallParam(kFALSE)
  47 {
  48   //
  49   // constructor
  50   //
  51   SetCrystalBallFunction(kFALSE);
  52   SetMassWndHigh(0.2);
  53   SetMassWndLow(0.5);
  54   for(Int_t iPar = 0; iPar < 16; iPar++) fParameters[iPar] = 0.;
  55   fParameters[9] = 1.;fParameters[11] = 1.;fParameters[12] = 1.;
  56   for(Int_t index=0; index<6; index++) fResolutionConstants[index] = 0.;
  57   AliInfo("Instance of AliBtoJPSItoEleCDFfitFCN-class created");
  58 }
  59 //_________________________________________________________________________________________________
  60 AliBtoJPSItoEleCDFfitFCN::AliBtoJPSItoEleCDFfitFCN(const AliBtoJPSItoEleCDFfitFCN& source) :
  61 TNamed(source),
  62 fFPlus(source.fFPlus),
  63 fFMinus(source.fFMinus),
  64 fFSym(source.fFSym),
  65 fIntegral(source.fIntegral),
  66 fintxFunB(source.fintxFunB),
  67 fintxDecayTimeBkgPos(source.fintxDecayTimeBkgPos),
  68 fintxDecayTimeBkgNeg(source.fintxDecayTimeBkgNeg),
  69 fintxDecayTimeBkgSym(source.fintxDecayTimeBkgSym),
  70 fintmMassSig(source.fintmMassSig),
  71 fintxRes(source.fintxRes),
  72 fintmMassBkg(source.fintmMassBkg),
  73 fhCsiMC(source.fhCsiMC),
  74 fMassWndHigh(source.fMassWndHigh),
  75 fMassWndLow(source.fMassWndLow),
  76 fCrystalBallParam(source.fCrystalBallParam)
  77 {
  78   //
  79   // Copy constructor
  80   //
  81   for(Int_t iPar = 0; iPar < 16; iPar++) fParameters[iPar] = source.fParameters[iPar];
  82   for(Int_t index=0; index<6; index++) fResolutionConstants[index] = source.fResolutionConstants[index];
  83 }
  84 //_________________________________________________________________________________________________
  85 AliBtoJPSItoEleCDFfitFCN& AliBtoJPSItoEleCDFfitFCN::operator=(const AliBtoJPSItoEleCDFfitFCN& source)
  86 {
  87   //
  88   // Assignment operator
  89   //
  90   if(&source == this) return *this;
  91   fFPlus = source.fFPlus;
  92   fFMinus = source.fFMinus;
  93   fFSym = source.fFSym;
  94   fIntegral = source.fIntegral;
  95   fintxFunB = source.fintxFunB;
  96   fintxDecayTimeBkgPos = source.fintxDecayTimeBkgPos;
  97   fintxDecayTimeBkgNeg = source.fintxDecayTimeBkgNeg;
  98   fintxDecayTimeBkgSym = source.fintxDecayTimeBkgSym;
  99   fintmMassSig = source.fintmMassSig;
 100   fintxRes = source.fintxRes;
 101   fintmMassBkg = source.fintmMassBkg;
 102   fhCsiMC = source.fhCsiMC;
 103   fCrystalBallParam = source.fCrystalBallParam;
 104
 105   for(Int_t iPar = 0; iPar < 16; iPar++) fParameters[iPar] = source.fParameters[iPar];
 106   for(Int_t index=0; index<6; index++) fResolutionConstants[index] = source.fResolutionConstants[index];
 107
 108   return *this;
 109 }
 110 //_________________________________________________________________________________________________
 111 AliBtoJPSItoEleCDFfitFCN::~AliBtoJPSItoEleCDFfitFCN()
 112 {
 113   //
 114   // Default destructor
 115   //
 116
 117   delete fhCsiMC;
 118   for(Int_t iPar = 0; iPar < 16; iPar++) fParameters[iPar] = 0.;
 119   for(Int_t index=0; index<6; index++) fResolutionConstants[index] = 0.;
 120 }
 121 //_________________________________________________________________________________________________
 122 Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateLikelihood(const Double_t* pseudoproperdecaytime,
 123            const Double_t* invariantmass, const Int_t ncand)
 124 {
 125   //
 126   // This function evaluates the Likelihood fnction
 127   // It returns the -Log(of the likelihood function)
 128   //
 129   Double_t f = 0.;
 130   Double_t ret = 0.;
 131
 132   for(Int_t i=0; i < ncand; i++) {
 133       f = EvaluateCDFfuncNorm(pseudoproperdecaytime[i],invariantmass[i]);
 134       if(f < 0.) {
 135         continue;
 136       }
 137       ret+=-1.*TMath::Log(f);
 138       }
 139   return ret;
 140 }
 141 //_________________________________________________________________________________________________
 142 void AliBtoJPSItoEleCDFfitFCN::SetAllParameters(const Double_t* parameters)
 143 {
 144   //
 145   // Sets array of FCN parameters
 146   //
 147   for(Int_t index = 0; index < 16; index++) fParameters[index] = parameters[index];
 148 }
 149 //_________________________________________________________________________________________________
 150 void AliBtoJPSItoEleCDFfitFCN::ComputeIntegral()
 151 {
 152   //
 153   // this function compute the integral of the likelihood function
 154   // (theoretical function) in order to normalize it to unity
 155   //
 156   Double_t np = 100.0;                 //number of integration steps
 157   Double_t npres = 200.0;              //number of integration steps for the resolution function
 158   Double_t npm = 200.;
 159   Double_t stepm;Double_t stepx;Double_t stepxres;       //integration step width in variable m,x
 160   Double_t mx=0.;Double_t xprime=0.;
 161   Double_t xlow = -4000.; Double_t xup = 4000.;
 162   Double_t i; Double_t j;
 163   Double_t sumx = 0.0;Double_t intx = 0.0;Double_t intm = 0.0;
 164   stepm = (fMassWndHigh-fMassWndLow)/npm;
 165   stepx = (xup-xlow)/np;
 166   stepxres = (xup-xlow)/npres;
 167
 168 // compute integrals for all the terms
 169
 170   Double_t iRes;
 171   Double_t intxRes = 0.0;
 172   Double_t sumxRes = 0.0;
 173        for(iRes = 1.0; iRes <= npres/2.; iRes++)  {
 174            xprime = xlow + (iRes - .5)*stepxres;
 175            sumxRes += ResolutionFunc(xprime);
 176            xprime = xup - (iRes - .5)*stepxres;
 177            sumxRes += ResolutionFunc(xprime);
 178            }
 179        intxRes = sumxRes*stepxres;
 180        SetIntegralRes(intxRes);
 181
 182 //
 183   Double_t iFunB;
 184   Double_t intxFunB = 0.0;
 185   Double_t sumxFunB = 0.0;
 186        for(iFunB = 1.0; iFunB <= np/2; iFunB++)  {
 187            xprime = xlow + (iFunB - .5)*stepx;
 188            sumxFunB += FunB(xprime);
 189            xprime = xup - (iFunB - .5)*stepx;
 190            sumxFunB += FunB(xprime);
 191            }
 192        intxFunB = sumxFunB*stepx;
 193        SetIntegralFunB(intxFunB);
 194
 195 //
 196   Double_t iDecayTimeBkgPos;
 197   Double_t intxDecayTimeBkgPos = 0.0;
 198   Double_t sumxDecayTimeBkgPos = 0.0;
 199        for(iDecayTimeBkgPos = 1.0; iDecayTimeBkgPos <= np/2; iDecayTimeBkgPos++)  {
 200            xprime = xlow + (iDecayTimeBkgPos - .5)*stepx;
 201            sumxDecayTimeBkgPos += FunBkgPos(xprime);
 202            xprime = xup - (iDecayTimeBkgPos - .5)*stepx;
 203            sumxDecayTimeBkgPos += FunBkgPos(xprime);
 204            }
 205        intxDecayTimeBkgPos = sumxDecayTimeBkgPos*stepx;
 206        SetIntegralBkgPos(intxDecayTimeBkgPos);
 207
 208 //
 209   Double_t iDecayTimeBkgNeg;
 210   Double_t intxDecayTimeBkgNeg = 0.0;
 211   Double_t sumxDecayTimeBkgNeg = 0.0;
 212        for(iDecayTimeBkgNeg = 1.0;  iDecayTimeBkgNeg<= np/2; iDecayTimeBkgNeg++)  {
 213            xprime = xlow + (iDecayTimeBkgNeg - .5)*stepx;
 214            sumxDecayTimeBkgNeg += FunBkgNeg(xprime);
 215            xprime = xup - (iDecayTimeBkgNeg - .5)*stepx;
 216            sumxDecayTimeBkgNeg += FunBkgNeg(xprime);
 217            }
 218        intxDecayTimeBkgNeg = sumxDecayTimeBkgNeg*stepx;
 219        SetIntegralBkgNeg(intxDecayTimeBkgNeg);
 220 //
 221   Double_t iDecayTimeBkgSym;
 222   Double_t intxDecayTimeBkgSym = 0.0;
 223   Double_t sumxDecayTimeBkgSym = 0.0;
 224        for(iDecayTimeBkgSym = 1.0; intxDecayTimeBkgSym <= np/2; intxDecayTimeBkgSym++)  {
 225            xprime = xlow + (intxDecayTimeBkgSym - .5)*stepx;
 226            sumxDecayTimeBkgSym += FunBkgSym(xprime);
 227            xprime = xup - (intxDecayTimeBkgSym - .5)*stepx;
 228            sumxDecayTimeBkgSym += FunBkgSym(xprime);
 229            }
 230        intxDecayTimeBkgSym = sumxDecayTimeBkgSym*stepx;
 231        SetIntegralBkgSym(intxDecayTimeBkgSym);
 232
 233 //
 234   Double_t iMassSig;
 235   Double_t intmMassSig = 0.0;
 236   Double_t summMassSig = 0.0;
 237        for(iMassSig = 1.0;  iMassSig<= npm/2.; iMassSig++)  {
 238            mx = fMassWndLow + (iMassSig - .5)*stepm;
 239            summMassSig += EvaluateCDFInvMassSigDistr(mx);
 240            mx = fMassWndHigh - (iMassSig - .5)*stepm;
 241            summMassSig += EvaluateCDFInvMassSigDistr(mx);
 242            }
 243        intmMassSig = summMassSig*stepm;
 244        SetIntegralMassSig(intmMassSig);
 245
 246 //
 247   Double_t iMassBkg;
 248   Double_t intmMassBkg = 0.0;
 249   Double_t summMassBkg = 0.0;
 250        for(iMassBkg = 1.0; iMassBkg <= npm/2.; iMassBkg++)  {
 251            mx = fMassWndLow + (iMassBkg - .5)*stepm;
 252            summMassBkg += EvaluateCDFInvMassBkgDistr(mx);
 253            mx = fMassWndHigh - (iMassBkg - .5)*stepm;
 254            summMassBkg += EvaluateCDFInvMassBkgDistr(mx);
 255            }
 256        intmMassBkg = summMassBkg*stepm;
 257        SetIntegralMassBkg(intmMassBkg);
 258
 259 //
 260 // Compute integral of the whole distribution function
 261 //
 262        for(i = 1.0; i <= np; i++)  {
 263            Double_t summ = 0.0;
 264            xprime = xlow + (i - .5)*stepx;
 265          for(j = 1.0; j <= npm/2; j++)  {
 266              mx = fMassWndLow + (j - .5)*stepm;
 267              summ += EvaluateCDFfunc(xprime,mx);
 268              mx = fMassWndHigh - (j - .5)*stepm;
 269              summ += EvaluateCDFfunc(xprime,mx);
 270              }
 271            intm = summ*stepm;
 272            sumx += intm;
 273            }
 274        intx = sumx*stepx;
 275        SetIntegral(intx);
 276
 277 }
 278 //_________________________________________________________________________________________________
 279 void AliBtoJPSItoEleCDFfitFCN::PrintStatus()
 280 {
 281   //
 282   //  Print the parameters of the fits
 283   //
 284   printf("\n");
 285   printf("actual value of fRadius---------------------------------------->> | %f \n", GetRadius());
 286   printf("actual value of fTheta ---------------------------------------->> | %f \n", GetTheta());
 287   printf("actual value of fPhi ------------------------------------------>> | %f \n", GetPhi());
 288   printf("actual value of fFPlus ---------------------------------------->> | %f \n", GetFPlus());
 289   printf("actual value of fFMinus --------------------------------------->> | %f \n", GetFMinus());
 290   printf("actual value of fFSym ----------------------------------------->> | %f \n", GetFSym());
 291   printf("actual value of fOneOvLamPlus --------------------------------->> | %f \n", GetLamPlus());
 292   printf("actual value of fOneOvLamMinus -------------------------------->> | %f \n", GetLamMinus());
 293   printf("actual value of fOneOvLamSym ---------------------------------->> | %f \n", GetLamSym());
 294   printf("actual value of fMassBkgSlope --------------------------------->> | %f \n", GetMassSlope());
 295   printf("actual value of fFractionJpsiFromBeauty ----------------------->> | %f \n", GetFractionJpsiFromBeauty());
 296   printf("actual value of fFsig ----------------------------------------->> | %f \n", GetFsig());
 297   if(fCrystalBallParam){
 298   printf("actual value of fCrystalBallMmean ----------------------------->> | %f \n", GetCrystalBallMmean());
 299   printf("actual value of fCrystalBallNexp ------------------------------>> | %f \n", GetCrystalBallNexp());
 300   printf("actual value of fCrystalBallSigma ----------------------------->> | %f \n", GetCrystalBallSigma());
 301   printf("actual value of fCrystalBallAlpha ----------------------------->> | %f \n", GetCrystalBallAlpha());
 302   printf("actual value of fCrystalBallNorm  ----------------------------->> | %f \n", GetCrystalBallNorm());
 303   }else{
 304    printf("actual value of fMpv ------------------------------------------>> | %f \n", GetCrystalBallMmean());
 305    printf("actual value of fConstRovL ------------------------------------>> | %f \n", GetCrystalBallNexp());
 306    printf("actual value of fSigmaL --------------------------------------->> | %f \n", GetCrystalBallSigma());
 307    printf("actual value of fSigmaR --------------------------------------->> | %f \n", GetCrystalBallAlpha());
 308   }
 309   printf("actual value of fSigmaResol ----------------------------------->> | %f \n", GetSigmaResol());
 310   printf("actual value of fNResol --------------------------------------->> | %f \n", GetNResol());
 311   printf("\n");
 312   printf("Actual value of normalization integral for FunB ------------------->> | %f \n", GetIntegralFunB());
 313   printf("Actual value of normalization integral for BkgPos ----------------->> | %f \n", GetIntegralBkgPos());
 314   printf("Actual value of normalization integral for BkgNeg ----------------->> | %f \n", GetIntegralBkgNeg());
 315   printf("Actual value of normalization integral for BkgSym ----------------->> | %f \n", GetIntegralBkgSym());
 316   printf("Actual value of normalization integral for MassSig ---------------->> | %f \n", GetIntegralMassSig());
 317   printf("Actual value of normalization integral for MassBkg ---------------->> | %f \n", GetIntegralMassBkg());
 318   printf("Actual value of normalization integral for Resolution ------------->> | %f \n", GetIntegralRes());
 319   printf("Actual value of normalization integral for FCN -------------------->> | %f \n", GetIntegral());
 320
 321   printf("\n");
 322 }
 323 //_________________________________________________________________________________________________
 324 void AliBtoJPSItoEleCDFfitFCN::SetResolutionConstants()
 325 {
 326   //
 327   //  This method must be update:
 328   //  for the time beeing the values are hard-wired.
 329   //  Implementations have to be done to set the values from outside
 330   //  (e.g. from a ConfigHF file) starting from an indipendent fit
 331   //  of primary JPSI distribution.
 332   //
 333
 334   fResolutionConstants[0]  = 8.; // mean sigma2/sigma1
 335   fResolutionConstants[1]  = 0.1675; // mean Integral2/Integral1
 336   fResolutionConstants[2]   = 1374.; // sigma2
 337   fResolutionConstants[3]  = 0.001022; // N2
 338   fResolutionConstants[4]  = 686.6; // mu2
 339 }
 340 //_________________________________________________________________________________________________
 341 Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFfunc(Double_t x, Double_t m) const
 342 {
 343   return fParameters[8]*EvaluateCDFfuncSignalPart(x,m) + (1. - fParameters[8])*EvaluateCDFfuncBkgPart(x,m);
 344 }
 345 //_________________________________________________________________________________________________
 346 Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFfuncNorm(Double_t x, Double_t m) const
 347 {
 348   return EvaluateCDFfunc(x,m)/fIntegral;
 349 }
 350 //_________________________________________________________________________________________________
 351 Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFfuncSignalPart(Double_t x, Double_t m) const
 352 {
 353   return EvaluateCDFDecayTimeSigDistr(x)*(EvaluateCDFInvMassSigDistr(m)/fintmMassSig);
 354 }
 355 //_________________________________________________________________________________________________
 356 Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeSigDistr(Double_t x) const
 357 {
 358   //
 359   // Implementation of the Background part of the Likelyhood function
 360   //
 361
 362   Double_t retvalue = 0.;
 363   Double_t FunBnorm = FunB(x)/fintxFunB;
 364   Double_t FunPnorm = ResolutionFunc(x)/fintxRes;
 365   retvalue = fParameters[7]*FunBnorm + (1. - fParameters[7])*FunPnorm;
 366   return retvalue;
 367 }
 368 //_________________________________________________________________________________________________
 369 Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFInvMassSigDistr(Double_t m) const
 370 {
 371   //
 372   // Parametrization of signal part invariant mass distribution
 373   // It can be either Crystal Ball function or sum of two Landau
 374   //
 375
 376   Double_t fitval = 0.;
 377
 378   if(fCrystalBallParam){
 379    Double_t t = (m-fParameters[9])/fParameters[11]; ;
 380    if (fParameters[12] < 0) t = -t;
 381
 382    Double_t absAlpha = TMath::Abs((Double_t)fParameters[12]);
 383
 384    if (t >= -absAlpha) {
 385       return fParameters[13]*TMath::Exp(-0.5*t*t);
 386       }
 387    else {
 388      Double_t a =  TMath::Power(fParameters[10]/absAlpha,fParameters[10])* TMath::Exp(-0.5*absAlpha*absAlpha);
 389      Double_t b= fParameters[10]/absAlpha - absAlpha;
 390     fitval = (fParameters[13]*a/TMath::Power(b - t, fParameters[10]));
 391     return fitval;
 392     }
 393    }else{
 394       Double_t t=-1*m;
 395       Double_t tmpv=-1*fParameters[9];
 396       fitval=TMath::Sqrt(TMath::Landau(t,tmpv,fParameters[11]));
 397       fitval += fParameters[10]*(TMath::Landau(m,fParameters[9],fParameters[12]));
 398       return fitval;
 399      }
 400 }
 401 //_________________________________________________________________________________________________
 402 Double_t AliBtoJPSItoEleCDFfitFCN::FunB(Double_t x) const
 403 {
 404   //
 405   // Parameterisation of the fit function for the x part of the Background
 406   //
 407
 408   Double_t np = 100.0;
 409   Double_t sc = 10.;
 410   Double_t sigma3 =  fResolutionConstants[2];
 411   Double_t xprime;
 412   Double_t sum = 0.0;
 413   Double_t xlow,xupp;
 414   Double_t step;
 415   Double_t i;
 416   xlow = x - sc * sigma3 ;
 417   xupp = x + sc * sigma3 ;
 418   step = (xupp-xlow) / np;
 419   Double_t CsiMCxprime = 0.;
 420   Double_t Resolutionxdiff = 0.;
 421   Double_t xdiff = 0.;
 422
 423   for(i=1.0; i<=np; i++) {
 424
 425       xprime = xlow + (i-.5) * step;
 426       CsiMCxprime = CsiMC(xprime);
 427       xdiff = xprime - x;
 428       Resolutionxdiff = ResolutionFunc(xdiff)/fintxRes; // normalized value
 429       sum += CsiMCxprime * Resolutionxdiff;
 430
 431       }
 432
 433   return step * sum ;
 434 }
 435 //_________________________________________________________________________________________________
 436 Double_t AliBtoJPSItoEleCDFfitFCN::FunP(Double_t x) const
 437 {
 438   //
 439   //  Parameterisation of the Prompt part for the x distribution
 440   //
 441
 442   return ResolutionFunc(x)/fintxRes; // normalized value
 443 }
 444 //_________________________________________________________________________________________________
 445 Double_t AliBtoJPSItoEleCDFfitFCN::CsiMC(Double_t x) const
 446 {
 447   //
 448   //  Distribution (template) of the x distribution for the x variable
 449   //  for the J/psi coming from Beauty hadrons
 450   //
 451
 452   Double_t returnvalue = 0.;
 453   returnvalue = fhCsiMC->GetBinContent(fhCsiMC->FindBin(x));
 454   return returnvalue;
 455 }
 456 //_________________________________________________________________________________________________
 457 Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFfuncBkgPart(Double_t x,Double_t m) const
 458 {
 459   //
 460   // Return the part of the likelihood function for the background hypothesis
 461   //
 462
 463   return EvaluateCDFDecayTimeBkgDistr(x)*(EvaluateCDFInvMassBkgDistr(m)/fintmMassBkg);
 464 }
 465 //_________________________________________________________________________________________________
 466 Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeBkgDistr(Double_t x) const
 467 {
 468   //
 469   // it returns the value of the probability to have a given x for the background
 470   //
 471
 472   Double_t ret = (1. - TMath::Power(fParameters[0],2.))*(ResolutionFunc(x)/fintxRes)
 473                      + TMath::Power(fParameters[0]*TMath::Cos(fParameters[1]),2.)*
 474                                               (FunBkgPos(x)/fintxDecayTimeBkgPos)
 475                      + TMath::Power(fParameters[0]*TMath::Sin(fParameters[1])*TMath::Sin(fParameters[2]),2.)*
 476                                               (FunBkgNeg(x)/fintxDecayTimeBkgNeg)
 477                      + TMath::Power(fParameters[0]*TMath::Sin(fParameters[1])*TMath::Cos(fParameters[2]),2.)*
 478                                               (FunBkgSym(x)/fintxDecayTimeBkgSym);
 479   return ret;
 480 }
 481 //_________________________________________________________________________________________________
 482 Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFInvMassBkgDistr(Double_t m) const
 483 {
 484   //
 485   // it returns the value of the probability to have a given mass for the background
 486   //
 487
 488   return 1/(fMassWndHigh-fMassWndLow) +
 489          fParameters[6] * m -
 490          fParameters[6] * ((fMassWndHigh+fMassWndLow)/2);
 491 }
 492 //_________________________________________________________________________________________________
 493 Double_t AliBtoJPSItoEleCDFfitFCN::FunBkgPos(Double_t x) const
 494 {
 495   //
 496   // exponential with positive slopes for the background part (x)
 497   //
 498
 499   Double_t np = 100.0;
 500   Double_t sc = 10.;
 501   Double_t sigma3 = fResolutionConstants[2];
 502   Double_t xprime;
 503   Double_t sum = 0.0;
 504   Double_t xlow,xupp;
 505   Double_t step;
 506   Double_t i;
 507   xlow = x - sc * sigma3 ;
 508   xupp = x + sc * sigma3 ;
 509   step = (xupp-xlow) / np;
 510
 511   for(i=1.0; i<=np/2; i++) {
 512
 513       xprime = xlow + (i-.5) * step;
 514        if (xprime > 0) {sum += fParameters[3] * TMath::Exp(-1*xprime*fParameters[3]) * (ResolutionFunc(xprime-x)/fintxRes);} else {sum = 0.;}
 515
 516       xprime = xupp - (i-.5) * step;
 517       if (xprime > 0) {sum += fParameters[3] * TMath::Exp(-1*xprime*fParameters[3]) * (ResolutionFunc(xprime-x)/fintxRes);} else {sum = 0.;}
 518
 519       }
 520
 521   return step * sum ;
 522 }
 523 //_________________________________________________________________________________________________
 524 Double_t AliBtoJPSItoEleCDFfitFCN::FunBkgNeg(Double_t x) const
 525 {
 526   //
 527   // exponential with negative slopes for the background part (x)
 528   //
 529
 530   Double_t np = 100.0;
 531   Double_t sc = 10.;
 532   Double_t sigma3 =  fResolutionConstants[2];
 533   Double_t xprime;
 534   Double_t sum = 0.0;
 535   Double_t xlow,xupp;
 536   Double_t step;
 537   Double_t i;
 538   xlow = x - sc * sigma3 ;
 539   xupp = x + sc * sigma3 ;
 540   step = (xupp-xlow) / np;
 541
 542   for(i=1.0; i<=np/2; i++) {
 543
 544       xprime = xlow + (i-.5) * step;
 545       if (xprime < 0) {sum += fParameters[4] * TMath::Exp(xprime*fParameters[4]) * (ResolutionFunc(xprime-x)/fintxRes);} else {sum = 0.;}
 546
 547       xprime = xupp - (i-.5) * step;
 548       if (xprime < 0) {sum += fParameters[4] * TMath::Exp(xprime*fParameters[4]) * (ResolutionFunc(xprime-x)/fintxRes);} else {sum = 0.;}
 549       }
 550
 551   return step * sum ;
 552 }
 553 //_________________________________________________________________________________________________
 554 Double_t AliBtoJPSItoEleCDFfitFCN::FunBkgSym(Double_t x) const
 555 {
 556   //
 557   // exponential with both positive and negative slopes for the background part (x)
 558   //
 559
 560   Double_t np = 100.0;
 561   Double_t sc = 10.;
 562   Double_t sigma3 =  fResolutionConstants[2];
 563   Double_t xprime;
 564   Double_t sum1 = 0.0;
 565   Double_t sum2 = 0.0;
 566   Double_t xlow,xupp;
 567   Double_t step;
 568   Double_t i;
 569   xlow = x - sc * sigma3 ;
 570   xupp = x + sc * sigma3 ;
 571   step = (xupp-xlow) / np;
 572
 573   for(i=1.0; i<=np/2; i++) {
 574
 575       xprime = xlow + (i-.5) * step;
 576       if (xprime > 0) {sum1 += 0.5 * fParameters[5]*TMath::Exp(-1*xprime*fParameters[5]) * (ResolutionFunc(xprime-x)/fintxRes);} else {sum1 = 0.;}
 577
 578       xprime = xupp - (i-.5) * step;
 579       if (xprime > 0) {sum1 += 0.5 * fParameters[5]*TMath::Exp(-1*xprime*fParameters[5]) * (ResolutionFunc(xprime-x)/fintxRes);} else {sum1 = 0.;}
 580       }
 581
 582   for(i=1.0; i<=np/2; i++) {
 583
 584       xprime = xlow + (i-.5) * step;
 585       if (xprime < 0) {sum2 += 0.5 * fParameters[5]*TMath::Exp(xprime*fParameters[5]) * (ResolutionFunc(xprime-x)/fintxRes);} else {sum2 = 0.;}
 586
 587       xprime = xupp - (i-.5) * step;
 588       if (xprime < 0) {sum2 += 0.5 * fParameters[5]*TMath::Exp(xprime*fParameters[5]) * (ResolutionFunc(xprime-x)/fintxRes);} else {sum2 = 0.;}
 589       }
 590
 591   return step*(sum1 + sum2) ;
 592 }
 593 //_________________________________________________________________________________________________
 594 Double_t AliBtoJPSItoEleCDFfitFCN::ResolutionFunc(Double_t x) const
 595 {
 596   //
 597   // parametrization with 2 gaus
 598   //
 599
 600   Double_t ret = 0.;
 601   Double_t x1 = x;
 602   Double_t x2 = x;
 603   //Double_t mean1 = 0.;
 604   Double_t mean2 = fResolutionConstants[4];
 605   Double_t sigma1 = fParameters[14];
 606   Double_t sigma2 = fResolutionConstants[2];
 607   Double_t n1 = fParameters[15];
 608   Double_t n2 = fResolutionConstants[3];
 609   Double_t arg1 = x1/sigma1;
 610   Double_t arg2 = (x2-mean2)/sigma2;
 611   Double_t sqrt2Pi = TMath::Sqrt(2*TMath::Pi());
 612
 613   ret = n2*((n1/n2)*TMath::Exp(-0.5*arg1*arg1) + TMath::Exp(-0.5*arg2*arg2));
 614
 615   return ret/(sqrt2Pi*(n1*sigma1+n2*sigma2));//return value is normalized
 616
 617 }
 618