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(0.),
  36 fhCsiMC(0x0),
  37 fMassWndHigh(0.),
  38 fMassWndLow(0.),
  39 fCrystalBallParam(kFALSE)
  40 {
  41   //
  42   // constructor
  43   //
  44   SetCrystalBallParam(kFALSE);
  45   SetMassWndHigh(0.2);
  46   SetMassWndLow(0.5);
  47   for(Int_t iPar = 0; iPar < 13; iPar++) fParameters[iPar] = 0.;
  48   fParameters[9] = 1.;fParameters[11] = 1.;fParameters[12] = 1.;
  49   for(Int_t index=0; index<7; index++) fResolutionConstants[index] = 0.;
  50   AliInfo("Instance of AliBtoJPSItoEleCDFfitFCN-class created");
  51 }
  52 //_________________________________________________________________________________________________
  53 AliBtoJPSItoEleCDFfitFCN::AliBtoJPSItoEleCDFfitFCN(const AliBtoJPSItoEleCDFfitFCN& source) :
  54 TNamed(source),
  55 fFPlus(source.fFPlus),
  56 fFMinus(source.fFMinus),
  57 fFSym(source.fFSym),
  58 fIntegral(source.fIntegral),
  59 fhCsiMC(source.fhCsiMC),
  60 fMassWndHigh(source.fMassWndHigh),
  61 fMassWndLow(source.fMassWndLow),
  62 fCrystalBallParam(source.fCrystalBallParam)
  63 {
  64   //
  65   // Copy constructor
  66   //
  67   for(Int_t iPar = 0; iPar < 13; iPar++) fParameters[iPar] = source.fParameters[iPar];
  68   for(Int_t index=0; index<7; index++) fResolutionConstants[index] = source.fResolutionConstants[index];
  69 }
  70 //_________________________________________________________________________________________________
  71 AliBtoJPSItoEleCDFfitFCN& AliBtoJPSItoEleCDFfitFCN::operator=(const AliBtoJPSItoEleCDFfitFCN& source)
  72 {
  73   //
  74   // Assignment operator
  75   //
  76   if(&source == this) return *this;
  77   fFPlus = source.fFPlus;
  78   fFMinus = source.fFMinus;
  79   fFSym = source.fFSym;
  80   fIntegral = source.fIntegral;
  81   fhCsiMC = source.fhCsiMC;
  82   fCrystalBallParam = source.fCrystalBallParam;
  83
  84   for(Int_t iPar = 0; iPar < 13; iPar++) fParameters[iPar] = source.fParameters[iPar];
  85   for(Int_t index=0; index<7; index++) fResolutionConstants[index] = source.fResolutionConstants[index];
  86
  87   return *this;
  88 }
  89 //_________________________________________________________________________________________________
  90 AliBtoJPSItoEleCDFfitFCN::~AliBtoJPSItoEleCDFfitFCN()
  91 {
  92   //
  93  // Default destructor
  94   //
  95
  96   delete fhCsiMC;
  97   for(Int_t iPar = 0; iPar < 13; iPar++) fParameters[iPar] = 0.;
  98   for(Int_t index=0; index<7; index++) fResolutionConstants[index] = 0.;
  99 }
 100 //_________________________________________________________________________________________________
 101 Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateLikelihood(const Double_t* pseudoproperdecaytime,
 102            const Double_t* invariantmass, const Int_t ncand)
 103 {
 104 //
 105 // This function evaluates the Likelihood fnction
 106 // It returns the -Log(of the likelihood function)
 107 //
 108   Double_t f = 0.;
 109   Double_t ret = 0.;
 110
 111   for(Int_t i=0; i < ncand; i++) {
 112       f = EvaluateCDFfuncNorm(pseudoproperdecaytime[i],invariantmass[i]);
 113       if(f < 0.) {
 114         //AliWarning("One negative contributors in the Log(Likely) ! ");
 115         continue;
 116       }
 117       ret+=-1.*TMath::Log(f);
 118       }
 119   return ret;
 120 }
 121 //_________________________________________________________________________________________________
 122 void AliBtoJPSItoEleCDFfitFCN::SetAllParameters(const Double_t* parameters)
 123 {
 124   //
 125   // Sets array of FCN parameters
 126   //
 127   for(Int_t index = 0; index < 13; index++) fParameters[index] = parameters[index];
 128 }
 129 //_________________________________________________________________________________________________
 130 void AliBtoJPSItoEleCDFfitFCN::ComputeIntegral()
 131 {
 132 //
 133 // this function compute the integral of the likelihood function
 134 // (theoretical function) in order to normalize it to unity
 135 //
 136   Double_t np = 50.0;                  //number of integration steps
 137   Double_t stepm;Double_t stepx;       //integration step width in variable m,x
 138   Double_t mx;Double_t xx;
 139   Double_t xlow = -4000.; Double_t xup = 4000.;
 140   Double_t i; Double_t j;
 141   Double_t sumx = 0.0;Double_t intx = 0.0;Double_t intm = 0.0;
 142   stepm = (fMassWndHigh-fMassWndLow)/np;
 143   stepx = (xup-xlow)/np;
 144
 145         for(i = 1.0; i <= np; i++)  {
 146             Double_t summ = 0.0;
 147             xx = xlow + (i - .5)*stepx;
 148           for(j = 1.0; j <= np/2; j++)  {
 149               mx = fMassWndLow + (j - .5)*stepm;
 150               summ += EvaluateCDFfunc(xx,mx);
 151               mx = fMassWndHigh - (j - .5)*stepm;
 152               summ += EvaluateCDFfunc(xx,mx);
 153               }
 154             intm = summ*stepm;
 155             sumx += intm;
 156             }
 157         intx = sumx*stepx;
 158         SetIntegral(intx);
 159 }
 160 //_________________________________________________________________________________________________
 161 void AliBtoJPSItoEleCDFfitFCN::PrintStatus()
 162 {
 163 //
 164 //  Print the parameters of the fits
 165 //
 166   printf("\n");
 167   printf("actual value of fRadius---------------------------------------->> | %f \n", GetRadius());
 168   printf("actual value of fTheta ---------------------------------------->> | %f \n", GetTheta());
 169   printf("actual value of fPhi ------------------------------------------>> | %f \n", GetPhi());
 170   printf("actual value of fFPlus ---------------------------------------->> | %f \n", GetFPlus());
 171   printf("actual value of fFMinus --------------------------------------->> | %f \n", GetFMinus());
 172   printf("actual value of fFSym ----------------------------------------->> | %f \n", GetFSym());
 173   printf("actual value of fOneOvLamPlus --------------------------------->> | %f \n", GetLamPlus());
 174   printf("actual value of fOneOvLamMinus -------------------------------->> | %f \n", GetLamMinus());
 175   printf("actual value of fOneOvLamSym ---------------------------------->> | %f \n", GetLamSym());
 176   printf("actual value of fMassBkgSlope --------------------------------->> | %f \n", GetMassSlope());
 177   printf("actual value of fFractionJpsiFromBeauty ----------------------->> | %f \n", GetFractionJpsiFromBeauty());
 178   printf("actual value of fFsig ----------------------------------------->> | %f \n", GetFsig());
 179   if(fCrystalBallParam){
 180   printf("actual value of fCrystalBallMmean ----------------------------->> | %f \n", GetCrystalBallMmean());
 181   printf("actual value of fCrystalBallNexp ------------------------------>> | %f \n", GetCrystalBallNexp());
 182   printf("actual value of fCrystalBallSigma ----------------------------->> | %f \n", GetCrystalBallSigma());
 183   printf("actual value of fCrystalBallAlpha ----------------------------->> | %f \n", GetCrystalBallAlpha());
 184   }else{
 185    printf("actual value of fMpv ------------------------------------------>> | %f \n", GetCrystalBallMmean());
 186    printf("actual value of fConstRovL ------------------------------------>> | %f \n", GetCrystalBallNexp());
 187    printf("actual value of fSigmaL --------------------------------------->> | %f \n", GetCrystalBallSigma());
 188    printf("actual value of fSigmaR --------------------------------------->> | %f \n", GetCrystalBallAlpha());
 189   }
 190   printf("\n");
 191   printf("Actual value of normalization integral for FCN ---------------->> | %f \n", GetIntegral());
 192   printf("\n");
 193 }
 194 //_________________________________________________________________________________________________
 195 void AliBtoJPSItoEleCDFfitFCN::SetResolutionConstants(Int_t BinNum)
 196 {
 197 //  This method must be update.
 198 //  For the time beeing the values are hard-wired.
 199 //  Implementations have to be done to set the values from outside (e.g. from a ConfigHF file)
 200 //
 201   switch(BinNum){
 202
 203    case(kallpt):
 204     fResolutionConstants[0]  = 0.326     ; fResolutionConstants[2] = 0.3622*535.9 ; fResolutionConstants[4] = 1.171*5535.9  ;//from fit integrated in pt
 205     fResolutionConstants[1]  = 0.0998*535.9; fResolutionConstants[3] = 0.1072   ; fResolutionConstants[5] = 0.04115    ; fResolutionConstants[6] = 1e-04;
 206     break;
 207    case(kptbin1):
 208     fResolutionConstants[0]  = 0.326     ; fResolutionConstants[2] = 0.3622*1087  ; fResolutionConstants[4] = 1.171*1087    ;//from fit integrated in pt
 209     fResolutionConstants[1]  = 0.04253*1087 ; fResolutionConstants[3] = 0.1482  ; fResolutionConstants[5] = 0.09778    ; fResolutionConstants[6] = 3.773e-04;
 210     break;
 211    case(kptbin2):
 212     fResolutionConstants[0]  = 0.326     ; fResolutionConstants[2] = 0.3622*661.5 ; fResolutionConstants[4] = 1.171*661.5   ;//from fit integrated in pt
 213     fResolutionConstants[1]  = 0.1*661.5    ; fResolutionConstants[3] = 0.2809  ; fResolutionConstants[5] =  0.09771   ; fResolutionConstants[6] = 1.916e-04;
 214     break;
 215    case(kptbin3):
 216     fResolutionConstants[0]  = 0.326     ; fResolutionConstants[2] = 0.3622*502.8 ; fResolutionConstants[4] = 1.171*502.8   ;//from fit integrated in pt
 217     fResolutionConstants[1]  = 0.1578*502.8 ; fResolutionConstants[3] = 0.3547  ; fResolutionConstants[5] = 0.09896    ; fResolutionConstants[6] = 5.241e-04;
 218     break;
 219    case(kptbin4):
 220     fResolutionConstants[0]  = 0.326     ; fResolutionConstants[2] = 0.3622*502.8 ; fResolutionConstants[4] = 1.171*502.8   ;//from fit integrated in pt
 221     fResolutionConstants[1] = 0.2048*415.9  ; fResolutionConstants[3] = 0.4265  ; fResolutionConstants[5] = 0.09597    ; fResolutionConstants[6] = 6.469e-04;
 222     break;
 223    case(kptbin5):
 224     fResolutionConstants[0]  = 0.326     ; fResolutionConstants[2] = 0.3622*502.8 ; fResolutionConstants[4] = 1.171*502.8   ;//from fit integrated in pt
 225     fResolutionConstants[1] = 0.2219*379.7  ; fResolutionConstants[3] = 0.5414 ;  fResolutionConstants[5] = 0.07506    ; fResolutionConstants[6] = 7.465e-04;
 226     break;
 227    case(kptbin6):
 228     fResolutionConstants[0]  = 0.326     ; fResolutionConstants[2] = 0.3622*502.8 ; fResolutionConstants[4] = 1.171*502.8   ;//from fit integrated in pt
 229     fResolutionConstants[1] = 0.2481*307.   ; fResolutionConstants[3] = 0.8073 ;  fResolutionConstants[5] = 0.09664    ; fResolutionConstants[6] = 8.481e-04;
 230     break;
 231    case(kptbin7):
 232     fResolutionConstants[0]  = 0.326     ; fResolutionConstants[2] = 0.3622*502.8 ; fResolutionConstants[4] = 1.171*502.8   ;//from fit integrated in pt
 233     fResolutionConstants[1] = 0.262*283.5   ; fResolutionConstants[3] = 0.9639 ;  fResolutionConstants[5] = 0.07943    ; fResolutionConstants[6] = 6.873e-04;
 234     break;
 235    case(kptbin8):
 236     fResolutionConstants[0]  = 0.326     ; fResolutionConstants[2] = 0.3622*502.8 ; fResolutionConstants[4] = 1.171*502.8   ;//from fit integrated in pt
 237     fResolutionConstants[1] =  0.4514*204.8 ; fResolutionConstants[3] = 0.98  ;   fResolutionConstants[5] = 0.1192     ; fResolutionConstants[6] = 8.646e-04;
 238     break;
 239    case(kptbin9):
 240     fResolutionConstants[0]  = 0.326     ; fResolutionConstants[2] = 0.3622*502.8 ; fResolutionConstants[4] = 1.171*502.8   ;//from fit integrated in pt
 241     fResolutionConstants[1] =  0.525*181.   ; fResolutionConstants[3] = 0.99  ;   fResolutionConstants[5] = 0.1097     ; fResolutionConstants[6] = 9.637e-04;
 242     break;
 243    }
 244 }
 245 //_________________________________________________________________________________________________
 246 Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFfunc(Double_t x, Double_t m) const
 247 {
 248   return fParameters[8]*EvaluateCDFfuncSignalPart(x,m) + (1. - fParameters[8])*EvaluateCDFfuncBkgPart(x,m);
 249 }
 250 //_________________________________________________________________________________________________
 251 Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFfuncNorm(Double_t x, Double_t m) const
 252 {
 253   return EvaluateCDFfunc(x,m)/fIntegral;
 254 }
 255 //_________________________________________________________________________________________________
 256 Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFfuncSignalPart(Double_t x, Double_t m) const
 257 {
 258   return EvaluateCDFDecayTimeSigDistr(x)*EvaluateCDFInvMassSigDistr(m);
 259 }
 260 //_________________________________________________________________________________________________
 261 Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeSigDistr(Double_t x) const
 262 {
 263 //
 264 // Implementation of the Background part of the Likelyhood function
 265 //
 266 //
 267   Double_t retvalue = 0.;
 268   retvalue = fParameters[7]*FunB(x) + (1. - fParameters[7])*FunP(x);
 269   return retvalue;
 270 }
 271 //_________________________________________________________________________________________________
 272 Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFInvMassSigDistr(Double_t m) const
 273 {
 274   //
 275   // Parametrization of signal part invariant mass distribution
 276   // It can be either Crystal Ball function or sum of two Landau
 277   //
 278   Double_t fitval = 0.;
 279
 280   if(fCrystalBallParam){
 281    Double_t fitvalCB = 0.;
 282    Double_t normFactorCB = 1.;
 283    Double_t arg = (m - fParameters[9])/fParameters[11];
 284    Double_t numfactor = fParameters[10];
 285    Double_t denomfactor = numfactor - TMath::Abs(fParameters[12]) - arg;
 286
 287    if(arg <= -1*TMath::Abs(fParameters[12])){
 288       Double_t exponent = fParameters[10]*TMath::Abs(fParameters[12]);
 289       Double_t numer = TMath::Exp(-0.5*fParameters[12]*fParameters[12])*TMath::Power(numfactor,exponent);
 290       Double_t denom = TMath::Power(denomfactor,exponent);
 291       fitvalCB += numer/denom;
 292       }
 293    if(arg > -1*TMath::Abs(fParameters[12])){
 294       fitvalCB += TMath::Exp(-0.5*arg*arg);
 295       }
 296    fitval = normFactorCB*fitvalCB;
 297    return fitval;
 298    }else{
 299       Double_t t=-1*m;
 300       Double_t tmpv=-1*fParameters[9];
 301       fitval=TMath::Sqrt(TMath::Landau(t,tmpv,fParameters[11]));
 302       fitval += fParameters[10]*(TMath::Landau(m,fParameters[9],fParameters[12]));
 303       return fitval;
 304      }
 305 }
 306 //_________________________________________________________________________________________________
 307 Double_t AliBtoJPSItoEleCDFfitFCN::FunB(Double_t x) const
 308 {
 309 //
 310 // Parameterisation of the fit function for the x part of the Background
 311 //
 312   Double_t np = 50.0;
 313   Double_t sc = 100.;
 314   Double_t sigma3 =  fResolutionConstants[2];
 315   Double_t xx;
 316   Double_t sum = 0.0;
 317   Double_t xlow,xupp;
 318   Double_t step;
 319   Double_t i;
 320   xlow = x - sc * sigma3 ;
 321   xupp = x + sc * sigma3 ;
 322   step = (xupp-xlow) / np;
 323   for(i=1.0; i<=np/2; i++) {
 324       xx = xlow + (i-.5) * step;
 325       sum += CsiMC(xx) * ResolutionFunc(xx-x);
 326
 327       xx = xupp - (i-.5) * step;
 328       sum += CsiMC(xx) * ResolutionFunc(xx-x);
 329       }
 330   return (step * sum) ;
 331 }
 332 //_________________________________________________________________________________________________
 333 Double_t AliBtoJPSItoEleCDFfitFCN::FunP(Double_t x) const
 334 //
 335 //  Parameterisation of the Prompt part for the x distribution
 336 //
 337 {
 338   return ResolutionFunc(x);
 339 }
 340 //_________________________________________________________________________________________________
 341 Double_t AliBtoJPSItoEleCDFfitFCN::CsiMC(Double_t x) const
 342 {
 343 //  Distribution (template) of the x distribution for the x variable
 344 //  for the J/psi coming from Beauty hadrons
 345 //
 346   Double_t returnvalue = 0.;
 347   returnvalue = fhCsiMC->GetBinContent(fhCsiMC->FindBin(x));
 348   return returnvalue;
 349 }
 350 //_________________________________________________________________________________________________
 351 Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFfuncBkgPart(Double_t x,Double_t m) const
 352 {
 353 //
 354 // Return the part of the likelihood function for the background hypothesis
 355 //
 356   return EvaluateCDFDecayTimeBkgDistr(x)*EvaluateCDFInvMassBkgDistr(m);
 357 }
 358 //_________________________________________________________________________________________________
 359 Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeBkgDistr(Double_t x) const
 360 {
 361 // it returns the value of the probability to have a given x for the background
 362 //
 363 //
 364   Double_t ret = (1 - fParameters[0]*fParameters[0])*ResolutionFunc(x);
 365   ret += FunBkgPos(x);
 366   ret += FunBkgNeg(x);
 367   ret += FunBkgSym(x);
 368   return ret;
 369 }
 370 //_________________________________________________________________________________________________
 371 Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFInvMassBkgDistr(Double_t m) const
 372 {
 373 //
 374 // it returns the value of the probability to have a given mass for the background
 375 //
 376   return 1/(fMassWndHigh-fMassWndLow)+fParameters[6]*(m-(fMassWndHigh+fMassWndLow)/2);
 377 }
 378 //_________________________________________________________________________________________________
 379 Double_t AliBtoJPSItoEleCDFfitFCN::FunBkgPos(Double_t x) const
 380 {
 381 //
 382 // exponential with positive slopes for the background part (x)
 383 //
 384   Double_t np = 50.0;
 385   Double_t sc = 100.;
 386   Double_t sigma3 = fResolutionConstants[2];
 387   Double_t xx;
 388   Double_t sum = 0.0;
 389   Double_t xlow,xupp;
 390   Double_t step;
 391   Double_t i;
 392   xlow = x - sc * sigma3 ;
 393   xupp = x + sc * sigma3 ;
 394   step = (xupp-xlow) / np;
 395   for(i=1.0; i<=np/2; i++) {
 396       xx = xlow + (i-.5) * step;
 397       if (xx > 0) sum += (fParameters[0]*TMath::Cos(fParameters[1]))*(fParameters[0]*TMath::Cos(fParameters[1]))*fParameters[3]*TMath::Exp(-1*xx*fParameters[3]) * ResolutionFunc(xx-x);
 398
 399       xx = xupp - (i-.5) * step;
 400       if (xx > 0) sum += fParameters[0]*TMath::Cos(fParameters[1])*(fParameters[0]*TMath::Cos(fParameters[1]))*fParameters[3]*TMath::Exp(-1*xx*fParameters[3]) * ResolutionFunc(xx-x);
 401       }
 402   return (step * sum) ;
 403 }
 404 //_________________________________________________________________________________________________
 405 Double_t AliBtoJPSItoEleCDFfitFCN::FunBkgNeg(Double_t x) const
 406 {
 407 //
 408 // exponential with negative slopes for the background part (x)
 409 //
 410   Double_t np = 50.0;
 411   Double_t sc = 100.;
 412   Double_t sigma3 =  fResolutionConstants[2];
 413   Double_t xx;
 414   Double_t sum = 0.0;
 415   Double_t xlow,xupp;
 416   Double_t step;
 417   Double_t i;
 418   xlow = x - sc * sigma3 ;
 419   xupp = x + sc * sigma3 ;
 420   step = (xupp-xlow) / np;
 421   for(i=1.0; i<=np/2; i++) {
 422       xx = xlow + (i-.5) * step;
 423       if (xx < 0) sum += (fParameters[0]*TMath::Sin(fParameters[1])*TMath::Sin(fParameters[2]))*(fParameters[0]*TMath::Sin(fParameters[1])*TMath::Sin(fParameters[2]))*
 424                          fParameters[4]*TMath::Exp(xx*fParameters[4]) * ResolutionFunc(xx-x);
 425       xx = xupp - (i-.5) * step;
 426       if (xx < 0) sum += (fParameters[0]*TMath::Sin(fParameters[1])*TMath::Sin(fParameters[2]))*(fParameters[0]*TMath::Sin(fParameters[1])*TMath::Sin(fParameters[2]))*
 427                          fParameters[4]*TMath::Exp(xx*fParameters[4]) * ResolutionFunc(xx-x);
 428       }
 429   return (step * sum) ;
 430 }
 431 //_________________________________________________________________________________________________
 432 Double_t AliBtoJPSItoEleCDFfitFCN::FunBkgSym(Double_t x) const
 433 {
 434 //
 435 // exponential with both positive and negative slopes for the background part (x)
 436 //
 437   Double_t np = 50.0;
 438   Double_t sc = 100.;
 439   Double_t sigma3 =  fResolutionConstants[2];
 440   Double_t xx;
 441   Double_t sum1 = 0.0;
 442   Double_t sum2 = 0.0;
 443   Double_t xlow,xupp;
 444   Double_t step;
 445   Double_t i;
 446   xlow = x - sc * sigma3 ;
 447   xupp = x + sc * sigma3 ;
 448   step = (xupp-xlow) / np;
 449   for(i=1.0; i<=np/2; i++) {
 450       xx = xlow + (i-.5) * step;
 451       if (xx > 0) sum1 += 0.5*(fParameters[0]*TMath::Sin(fParameters[1])*TMath::Cos(fParameters[2]))*(fParameters[0]*TMath::Sin(fParameters[1])*TMath::Cos(fParameters[2]))*
 452                               fParameters[5]*TMath::Exp(-1*xx*fParameters[5]) * ResolutionFunc(xx-x);
 453
 454       xx = xupp - (i-.5) * step;
 455       if (xx > 0) sum1 += 0.5*(fParameters[0]*TMath::Sin(fParameters[1])*TMath::Cos(fParameters[2]))*(fParameters[0]*TMath::Sin(fParameters[1])*TMath::Cos(fParameters[2]))*
 456                               fParameters[5]*TMath::Exp(-1*xx*fParameters[5]) * ResolutionFunc(xx-x);
 457       }
 458   for(i=1.0; i<=np/2; i++) {
 459       xx = xlow + (i-.5) * step;
 460       if (xx < 0) sum2 += 0.5*(fParameters[0]*TMath::Sin(fParameters[1])*TMath::Cos(fParameters[2]))*(fParameters[0]*TMath::Sin(fParameters[1])*TMath::Cos(fParameters[2]))*
 461                               fParameters[5]*TMath::Exp(xx*fParameters[5]) * ResolutionFunc(xx-x);
 462
 463       xx = xupp - (i-.5) * step;
 464       if (xx < 0) sum2 += 0.5*(fParameters[0]*TMath::Sin(fParameters[1])*TMath::Cos(fParameters[2]))*(fParameters[0]*TMath::Sin(fParameters[1])*TMath::Cos(fParameters[2]))*
 465                               fParameters[5]*TMath::Exp(xx*fParameters[5]) * ResolutionFunc(xx-x);
 466       }
 467   return step*(sum1 + sum2) ;
 468 }
 469 //_________________________________________________________________________________________________
 470 Double_t AliBtoJPSItoEleCDFfitFCN::ResolutionFunc(Double_t x) const
 471 {
 472   //
 473   //parametrization with 2 gaus + 1 exponential + 1 constant
 474   //
 475   Double_t arg=0;
 476   arg=x/fResolutionConstants[1];
 477   Double_t ret=TMath::Exp(-0.5*arg*arg);
 478   arg=x/fResolutionConstants[2];
 479   ret+=fResolutionConstants[3]*TMath::Exp(-0.5*arg*arg);
 480   arg=x/fResolutionConstants[4];
 481   if(x > 0)  { ret+=fResolutionConstants[5]*TMath::Exp(-arg); }
 482   if(x <= 0) { ret+=fResolutionConstants[5]*TMath::Exp(arg); }
 483   return fResolutionConstants[0]*(ret + fResolutionConstants[6]);
 484 }