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 <TFormula.h>
  17 #include <TF1.h>
  18 #include <TCanvas.h>
  19 #include <TMath.h>
  20 #include <TFile.h>
  21
  22 #include "AliDielectronBtoJPSItoEleCDFfitFCN.h"
  23
  24 //_________________________________________________________________________
  25 //                        Class AliDielectronBtoJPSItoEleCDFfitFCN
  26 //                   Definition of main function used in
  27 //                     unbinned log-likelihood fit for
  28 //                 the channel B -> JPsi + X -> e+e- + X
  29 //
  30 //                           Origin: C.Di Giglio
  31 //       Contact: Carmelo.Digiglio@ba.infn.it , Giuseppe.Bruno@ba.infn.it
  32 //_________________________________________________________________________
  33
  34 ClassImp(AliDielectronBtoJPSItoEleCDFfitFCN)
  35
  36         //_________________________________________________________________________________________________
  37         AliDielectronBtoJPSItoEleCDFfitFCN::AliDielectronBtoJPSItoEleCDFfitFCN() :
  38                 fFPlus(0.),
  39                 fFMinus(0.),
  40                 fFSym(0.),
  41                 fintmMassSig(1.),
  42                 fintmMassBkg(1.),
  43                 fhCsiMC(0x0),
  44                 fShiftTemplate(0.),
  45                 fMassWndHigh(0.),
  46                 fMassWndLow(0.),
  47                 fCrystalBallParam(kFALSE),
  48                 fChangeResolution(1.),
  49                 fChangeMass(1.),
  50                 fWeights(0),
  51                 fLoadFunctions(kFALSE),
  52                 fMultivariate(kFALSE),
  53                 fFunBSaved(0x0),
  54                 fFunBkgSaved(0x0),
  55                 fResParams(0x0),
  56                 fBkgParams(0x0),
  57                 fMassWindows(0x0),
  58                 fPtWindows(0x0),
  59                 fExponentialParam(kTRUE),
  60                 fSignalBinForExtrapolation(0)
  61 {
  62         //
  63         // constructor
  64         //
  65         SetMassWndHigh(0.2);
  66         SetMassWndLow(0.5);
  67         fWeightType[0] = 1.; fWeightType[1] = 1.; fWeightType[2] = 1.;
  68         for(Int_t iPar = 0; iPar < 49; iPar++) fParameters[iPar] = 0.;
  69         fParameters[9] = 1.;fParameters[11] = 1.;fParameters[12] = 1.;
  70
  71
  72         AliInfo("Instance of AliDielectronBtoJPSItoEleCDFfitFCN-class created");
  73 }
  74 //_________________________________________________________________________________________________
  75 AliDielectronBtoJPSItoEleCDFfitFCN::AliDielectronBtoJPSItoEleCDFfitFCN(const AliDielectronBtoJPSItoEleCDFfitFCN& source) :
  76         TNamed(source),
  77         fFPlus(source.fFPlus),
  78         fFMinus(source.fFMinus),
  79         fFSym(source.fFSym),
  80         fintmMassSig(source.fintmMassSig),
  81         fintmMassBkg(source.fintmMassBkg),
  82         fhCsiMC(source.fhCsiMC),
  83         fShiftTemplate(source.fShiftTemplate),
  84         fMassWndHigh(source.fMassWndHigh),
  85         fMassWndLow(source.fMassWndLow),
  86         fCrystalBallParam(source.fCrystalBallParam),
  87         fChangeResolution(source.fChangeResolution),
  88         fChangeMass(source.fChangeMass),
  89         fWeights(source.fWeights),
  90         fLoadFunctions(source.fLoadFunctions),
  91         fMultivariate(source.fMultivariate),
  92         fFunBSaved(source.fFunBSaved),
  93         fFunBkgSaved(source.fFunBkgSaved),
  94         fResParams(source.fResParams),
  95         fBkgParams(source.fBkgParams),
  96         fMassWindows(source.fMassWindows),
  97         fPtWindows(source.fPtWindows),
  98         fExponentialParam(source.fExponentialParam),
  99         fSignalBinForExtrapolation(source.fSignalBinForExtrapolation)
 100 {
 101         //
 102         // Copy constructor
 103         //
 104         for(Int_t iPar = 0; iPar < 49; iPar++) fParameters[iPar] = source.fParameters[iPar];
 105         for(Int_t iW=0; iW<2; iW++) fWeightType[iW] = source.fWeightType[iW];
 106 }
 107 //_________________________________________________________________________________________________
 108 AliDielectronBtoJPSItoEleCDFfitFCN& AliDielectronBtoJPSItoEleCDFfitFCN::operator=(const AliDielectronBtoJPSItoEleCDFfitFCN& source)
 109 {
 110         //
 111         // Assignment operator
 112         //
 113         if(&source == this) return *this;
 114         fFPlus = source.fFPlus;
 115         fFMinus = source.fFMinus;
 116         fFSym = source.fFSym;
 117         fintmMassSig = source.fintmMassSig;
 118         fintmMassBkg = source.fintmMassBkg;
 119         fhCsiMC = source.fhCsiMC;
 120         fLoadFunctions = source.fLoadFunctions;
 121         fMultivariate = source.fMultivariate;
 122         fFunBSaved = source.fFunBSaved;
 123         fFunBkgSaved = source.fFunBkgSaved;
 124         fResParams = source.fResParams;
 125         fBkgParams = source.fBkgParams;
 126         fMassWindows = source.fMassWindows;
 127         fPtWindows = source.fPtWindows;
 128         fShiftTemplate = source.fShiftTemplate;
 129         fCrystalBallParam = source.fCrystalBallParam;
 130         fExponentialParam = source.fExponentialParam;
 131
 132         for(Int_t iPar = 0; iPar < 49; iPar++) fParameters[iPar] = source.fParameters[iPar];
 133         return *this;
 134 }
 135 //_________________________________________________________________________________________________
 136 AliDielectronBtoJPSItoEleCDFfitFCN::~AliDielectronBtoJPSItoEleCDFfitFCN()
 137 {
 138         //
 139         // Default destructor
 140         //
 141
 142         delete fhCsiMC;
 143         delete fFunBSaved;
 144         delete fFunBkgSaved;
 145         for(Int_t iPar = 0; iPar < 49; iPar++) fParameters[iPar] = 0.;
 146 }
 147 //_________________________________________________________________________________________________
 148 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateLikelihood(const Double_t* pseudoproperdecaytime,
 149                 const Double_t* invariantmass, const Double_t *pt, const Int_t *type, Int_t ncand) const
 150 {
 151         //
 152         // This function evaluates the Likelihood fnction
 153         // It returns the -Log(of the likelihood function)
 154         //
 155         Double_t f = 0.;
 156         Double_t ret = 0.;
 157
 158         for(Int_t i=0; i < ncand; i++) {
 159                 f = EvaluateCDFfuncNorm(pseudoproperdecaytime[i],invariantmass[i],pt[i],type[i]);
 160                 if(f <= 0.) continue;
 161                 ret += -2.*TMath::Log(f);
 162         }
 163         return ret;
 164 }
 165 //_________________________________________________________________________________________________
 166 void AliDielectronBtoJPSItoEleCDFfitFCN::SetAllParameters(const Double_t* parameters)
 167 {
 168         //
 169         // Sets array of FCN parameters
 170         //
 171         for(Int_t index = 0; index < 49; index++) fParameters[index] = parameters[index];
 172 }
 173 //_________________________________________________________________________________________________
 174 void AliDielectronBtoJPSItoEleCDFfitFCN::ComputeMassIntegral()
 175 {
 176         //
 177         // this function compute the integral of the likelihood function
 178         // (theoretical function) in order to normalize it to unity
 179         //
 180         Double_t npm = 20000.;
 181         Double_t stepm;
 182         Double_t mx=0.;
 183         stepm = (fMassWndHigh-fMassWndLow)/npm;
 184         // compute integrals for  invariant mass terms
 185
 186         Double_t iMassSig;
 187         Double_t intmMassSig = 0.0;
 188         Double_t summMassSig = 0.0;
 189         for(iMassSig = 1.0;  iMassSig<= npm/2.; iMassSig++)  {
 190                 mx = fMassWndLow + (iMassSig - .5)*stepm;
 191                 summMassSig += EvaluateCDFInvMassSigDistr(mx);
 192                 mx = fMassWndHigh - (iMassSig - .5)*stepm;
 193                 summMassSig += EvaluateCDFInvMassSigDistr(mx);
 194         }
 195         intmMassSig = summMassSig*stepm;
 196         SetIntegralMassSig(intmMassSig);
 197         //
 198         Double_t iMassBkg;
 199         Double_t intmMassBkg = 0.0;
 200         Double_t summMassBkg = 0.0;
 201         for(iMassBkg = 1.0; iMassBkg <= npm/2.; iMassBkg++)  {
 202                 mx = fMassWndLow + (iMassBkg - .5)*stepm;
 203                 summMassBkg += EvaluateCDFInvMassBkgDistr(mx);
 204                 mx = fMassWndHigh - (iMassBkg - .5)*stepm;
 205                 summMassBkg += EvaluateCDFInvMassBkgDistr(mx);
 206         }
 207         intmMassBkg = summMassBkg*stepm;
 208         if(intmMassBkg < 1.e-05) intmMassBkg = 1.;
 209         SetIntegralMassBkg(intmMassBkg);
 210 }
 211 //_________________________________________________________________________________________________
 212 void AliDielectronBtoJPSItoEleCDFfitFCN::PrintStatus()
 213 {
 214         //
 215         //  Print the parameters of the fits
 216         //
 217         printf("\n");
 218         // background param
 219         printf("actual value of fWeightRes------------------------------------->> | %f \n", GetResWeight());
 220         printf("actual value of fPos ------------------------------------------>> | %f \n", GetFPlus());
 221         printf("actual value of fNeg ------------------------------------------>> | %f \n", GetFMinus());
 222         printf("actual value of fSym ------------------------------------------>> | %f \n", GetFSym());
 223         printf("actual value of fSym1 ----------------------------------------->> | %f \n", GetFSym1());
 224         printf("actual value of fOneOvLamPlus --------------------------------->> | %f \n", GetLamPlus());
 225         printf("actual value of fOneOvLamMinus -------------------------------->> | %f \n", GetLamMinus());
 226         printf("actual value of fOneOvLamSym ---------------------------------->> | %f \n", GetLamSym());
 227         printf("actual value of fOneOvLamSym1 --------------------------------->> | %f \n", GetLamSym1());
 228         printf("actual value of fFractionJpsiFromBeauty ----------------------->> | %f \n", GetFractionJpsiFromBeauty());
 229         printf("actual value of fFsig ----------------------------------------->> | %f \n", GetFsig());
 230
 231         if(fCrystalBallParam){
 232                 printf("actual value of fCrystalBallMmean ----------------------------->> | %f \n", GetCrystalBallMmean());
 233                 printf("actual value of fCrystalBallNexp ------------------------------>> | %f \n", GetCrystalBallNexp());
 234                 printf("actual value of fCrystalBallSigma ----------------------------->> | %f \n", GetCrystalBallSigma());
 235                 printf("actual value of fCrystalBallAlpha ----------------------------->> | %f \n", GetCrystalBallAlpha());
 236                 printf("actual value of fCrystalBallNorm  ----------------------------->> | %f \n", GetCrystalBallNorm());
 237         }else{
 238                 printf("actual value of fMpv ------------------------------------------>> | %f \n", GetCrystalBallMmean());
 239                 printf("actual value of fConstRovL ------------------------------------>> | %f \n", GetCrystalBallNexp());
 240                 printf("actual value of fSigmaL --------------------------------------->> | %f \n", GetCrystalBallSigma());
 241                 printf("actual value of fSigmaR --------------------------------------->> | %f \n", GetCrystalBallAlpha());
 242         }
 243
 244         // back Mass func
 245         if(fExponentialParam){
 246                 printf("actual value of normBkg ----------------------------------------->> | %f \n", GetBkgInvMassNorm());
 247                 printf("actual value of meanBkg ----------------------------------------->> | %f \n", GetBkgInvMassMean());
 248                 printf("actual value of slopeBkg ---------------------------------------->> | %f \n", GetBkgInvMassSlope());
 249                 printf("actual value of constBkg ---------------------------------------->> | %f \n", GetBkgInvMassConst());
 250         }else{
 251                 printf("actual value of m^{0} ------------------------------------------->> | %f \n", GetBkgInvMassNorm());
 252                 printf("actual value of m^{1} ------------------------------------------->> | %f \n", GetBkgInvMassMean());
 253                 printf("actual value of m^{2} ------------------------------------------->> | %f \n", GetBkgInvMassSlope());
 254                 printf("actual value of m^{3} ------------------------------------------->> | %f \n", GetBkgInvMassConst());
 255                 printf("actual value of m^{4} ------------------------------------------->> | %f \n", GetPolyn4());
 256                 printf("actual value of m^{5} ------------------------------------------->> | %f \n", GetPolyn5());
 257         }
 258
 259         // resolution func (First-First)
 260         printf("actual value of norm1Gauss (FF) --------------------------------->> | %f \n", GetNormGaus1ResFunc(2));
 261         printf("actual value of norm2Gauss (FF) --------------------------------->> | %f \n", GetNormGaus2ResFunc(2));
 262         printf("actual value of fMean1Res (FF) ---------------------------------->> | %f \n", GetResMean1(2));
 263         printf("actual value of fSigma1Res (FF) --------------------------------->> | %f \n", GetResSigma1(2));
 264         printf("actual value of fMean2Res (FF) ---------------------------------->> | %f \n", GetResMean2(2));
 265         printf("actual value of fSigma2Res (FF) --------------------------------->> | %f \n", GetResSigma2(2));
 266
 267         printf("actual value of alfaRes (FF) ------------------------------------>> | %f \n", GetResAlfa(2));
 268         printf("actual value of lambdaRes (FF) ---------------------------------->> | %f \n", GetResLambda(2));
 269         printf("actual value of normExpRes (FF) --------------------------------->> | %f \n", GetResNormExp(2));
 270
 271         // resolution func (First-Second)
 272         printf("actual value of norm1Gauss (FS) --------------------------------->> | %f \n", GetNormGaus1ResFunc(1));
 273         printf("actual value of norm2Gauss (FS) --------------------------------->> | %f \n", GetNormGaus2ResFunc(1));
 274         printf("actual value of fMean1Res (FS) ---------------------------------->> | %f \n", GetResMean1(1));
 275         printf("actual value of fSigma1Res (FS) --------------------------------->> | %f \n", GetResSigma1(1));
 276         printf("actual value of fMean2Res (FS) ---------------------------------->> | %f \n", GetResMean2(1));
 277         printf("actual value of fSigma2Res (FS) --------------------------------->> | %f \n", GetResSigma2(1));
 278
 279         printf("actual value of alfaRes (FS) ------------------------------------>> | %f \n", GetResAlfa(1));
 280         printf("actual value of lambdaRes (FS) ---------------------------------->> | %f \n", GetResLambda(1));
 281         printf("actual value of normExpRes (FS) --------------------------------->> | %f \n", GetResNormExp(1));
 282
 283         // resolution func (Second-Second)
 284         printf("actual value of norm1Gauss (SS) --------------------------------->> | %f \n", GetNormGaus1ResFunc(0));
 285         printf("actual value of norm2Gauss (SS) --------------------------------->> | %f \n", GetNormGaus2ResFunc(0));
 286         printf("actual value of fMean1Res (SS) ---------------------------------->> | %f \n", GetResMean1(0));
 287         printf("actual value of fSigma1Res (SS) --------------------------------->> | %f \n", GetResSigma1(0));
 288         printf("actual value of fMean2Res (SS) ---------------------------------->> | %f \n", GetResMean2(0));
 289         printf("actual value of fSigma2Res (SS) --------------------------------->> | %f \n", GetResSigma2(0));
 290
 291         printf("actual value of alfaRes (SS) ------------------------------------>> | %f \n", GetResAlfa(0));
 292         printf("actual value of lambdaRes (SS) ---------------------------------->> | %f \n", GetResLambda(0));
 293         printf("actual value of normExpRes (SS) --------------------------------->> | %f \n", GetResNormExp(0));
 294
 295         printf("\n");
 296         // integrals constants
 297         printf("Actual value of normalization integral for MassSig ---------------->> | %f \n", GetIntegralMassSig());
 298         printf("Actual value of normalization integral for MassBkg ---------------->> | %f \n", GetIntegralMassBkg());
 299
 300         printf("\n");
 301 }
 302 //_________________________________________________________________________________________________
 303 void AliDielectronBtoJPSItoEleCDFfitFCN::SetResolutionConstants(const Double_t* resolutionConst, Int_t type)
 304 {
 305         //
 306         // Resolution function is parametrized as the sum of two gaussian
 307         //
 308         Int_t index = (2-type)*9;
 309         fParameters[20+index]=resolutionConst[1]; //mean1
 310         fParameters[22+index]=resolutionConst[4]; //mean2
 311         fParameters[18+index]=resolutionConst[0]; //norm1
 312         fParameters[21+index]=resolutionConst[2]; //sigma1
 313         fParameters[23+index]=resolutionConst[5]; //sigma2
 314         fParameters[19+index]=resolutionConst[3]; //norm2
 315
 316         //exp values
 317         fParameters[24+index]=resolutionConst[6]; //alfa
 318         fParameters[25+index]=resolutionConst[7]; //lambda
 319         fParameters[26+index]=resolutionConst[8]; //norm3
 320         return;
 321 }
 322
 323 //________________________________________________________________________________________________
 324 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFfunc(Double_t x, Double_t m, Double_t pt, Int_t type) const
 325 {
 326         // evaluate likelihood function
 327         //printf("CDF func ---> x = %f m = %f pt = %f type = %d \n",x,m,pt,type);
 328         return fParameters[8]*EvaluateCDFfuncSignalPart(x,m, pt, type) + (1. - fParameters[8])*EvaluateCDFfuncBkgPart(x,m, pt, type);
 329 }
 330
 331 //_________________________________________________________________________________________________
 332 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFfuncNorm(Double_t x, Double_t m, Double_t pt, Int_t type) const
 333 {
 334         // evaluate likelihood function
 335         return EvaluateCDFfunc(x,m,pt, type);
 336 }
 337
 338 //_________________________________________________________________________________________________
 339 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFfuncSignalPart(Double_t x, Double_t m, Double_t pt, Int_t type) const
 340 {
 341   // evaluate psproper signal
 342   return EvaluateCDFDecayTimeSigDistr(x,pt, type)*(EvaluateCDFInvMassSigDistr(m)/fintmMassSig);
 343 }
 344
 345 //_________________________________________________________________________________________________
 346 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeSigDistr(Double_t x, Double_t pt, Int_t type) const
 347 {
 348         //
 349         // Implementation of the Background part of the Likelyhood function
 350         //
 351         Double_t retvalue = 0.;
 352         Double_t funBnorm =  fMultivariate ? FunBsaved(x, pt, type) : FunB(x,pt, type)  ;
 353         Double_t funPnorm = ResolutionFunc(x, pt, type);
 354         retvalue = fParameters[7]*funBnorm + (1. - fParameters[7])*funPnorm;
 355         return retvalue;
 356 }
 357
 358 //_________________________________________________________________________________________________
 359 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFInvMassSigDistr(Double_t m) const
 360 {
 361         //
 362         // Parametrization of signal part invariant mass distribution
 363         // It can be either Crystal Ball function or sum of two Landau
 364         //
 365         Double_t fitval = 0.;
 366         // change inv Mass RMS fChangeMass
 367         if(fCrystalBallParam){
 368                 Double_t t = ((m-fParameters[9])/fChangeMass)/fParameters[11]; ;
 369                 if (fParameters[12] < 0) t = -t;
 370                 Double_t absAlpha = TMath::Abs((Double_t)fParameters[12]);
 371
 372                 if (t >= -absAlpha) {
 373                         return fParameters[13]*TMath::Exp(-0.5*t*t);
 374                 }
 375                 else {
 376                         Double_t a =  TMath::Power(fParameters[10]/absAlpha,fParameters[10])* TMath::Exp(-0.5*absAlpha*absAlpha);
 377                         Double_t b= fParameters[10]/absAlpha - absAlpha;
 378                         fitval = (fParameters[13]*a/TMath::Power(b - t, fParameters[10]));
 379                         return fitval;
 380                 }
 381         }else{
 382                 Double_t t=-1*m;
 383                 Double_t tmpv=-1*fParameters[9];
 384                 fitval=TMath::Sqrt(TMath::Landau(t,tmpv,fParameters[11]));
 385                 fitval += fParameters[10]*(TMath::Landau(m,fParameters[9],fParameters[12]));
 386                 return fitval;
 387         }
 388 }
 389 //_________________________________________________________________________________________________
 390 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::FunB(Double_t x, Double_t pt, Int_t type) const
 391 {
 392         //
 393         // Parameterisation of the fit function for the x part of the Background
 394         //
 395         Double_t np = 1000.0;
 396         Double_t sc = 10.;
 397         Double_t sigma3 = 1000.;
 398         Double_t xprime;
 399         Double_t sum = 0.0;
 400         Double_t xlow,xupp;
 401         Double_t step;
 402         Double_t i;
 403         xlow = x - sc * sigma3 ;
 404         xupp = x + sc * sigma3 ;
 405         step = (xupp-xlow) / np;
 406         Double_t csiMCxprime = 0.;
 407         Double_t resolutionxdiff = 0.;
 408         Double_t xdiff = 0.;
 409
 410         for(i=1.0; i<=np; i++){
 411                 xprime = xlow + (i-.5) * step;
 412                 csiMCxprime = CsiMC(xprime);
 413                 xdiff = xprime - x;
 414                 resolutionxdiff = ResolutionFunc(xdiff, pt, type); // normalized value
 415                 sum += csiMCxprime * resolutionxdiff;
 416         }
 417
 418         return step * sum ;
 419 }
 420
 421 //_________________________________________________________________________________________________
 422 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::FunP(Double_t x, Double_t pt, Int_t type) const
 423 {
 424         //
 425         //  Parameterisation of the Prompt part for the x distribution
 426         //
 427         return ResolutionFunc(x, pt, type);
 428 }
 429
 430
 431 //_________________________________________________________________________________________________
 432 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::CsiMC(Double_t x) const
 433 {
 434         //
 435         //  Distribution (template) of the x distribution for the x variable
 436         //  for the J/psi coming from Beauty hadrons
 437         //
 438         Double_t returnvalue = 0.;
 439
 440         if((fhCsiMC->FindBin(x-fShiftTemplate) > 0) && (fhCsiMC->FindBin(x-fShiftTemplate) < fhCsiMC->GetNbinsX()+1))
 441         returnvalue = fhCsiMC->Interpolate(x-fShiftTemplate);
 442
 443
 444         return returnvalue;
 445 }
 446
 447 //_________________________________________________________________________________________________
 448 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFfuncBkgPart(Double_t x,Double_t m, Double_t pt, Int_t type) const
 449 {
 450         //
 451         // Return the part of the likelihood function for the background hypothesis
 452         //
 453          Double_t bkgValx = fMultivariate ? EvaluateCDFDecayTimeBkgDistrSaved(x,type,m,pt) : EvaluateCDFDecayTimeBkgDistr(x,type,m,pt);
 454         return bkgValx*(EvaluateCDFInvMassBkgDistr(m)/fintmMassBkg);
 455 }
 456
 457
 458 //_________________________________________________________________________________________________
 459 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeBkgDistr(Double_t x, Int_t type, Double_t m, Double_t pt) const
 460 {
 461         //
 462         // it returns the value of the probability to have a given x for the background
 463         // in the pt, m , type correspondent range
 464         //
 465         Double_t ret = 0.;
 466         if(fMultivariate)
 467         ret = EvaluateCDFDecayTimeBkgDistrDifferential(x,type,m,pt);
 468         else ret = fParameters[0]/(fParameters[0]+fParameters[1]+fParameters[2]+fParameters[3]+fParameters[46])*ResolutionFunc(x, pt, type) + fParameters[1]/(fParameters[0]+fParameters[1]+fParameters[2]+fParameters[3]+fParameters[46])*FunBkgPos(x, pt,type) +  fParameters[2]/(fParameters[0]+fParameters[1]+fParameters[2]+fParameters[3]+fParameters[46])*FunBkgNeg(x,pt,type) + fParameters[3]/(fParameters[0]+fParameters[1]+fParameters[2]+fParameters[3]+fParameters[46])*FunBkgSym(x, pt,type) + fParameters[46]/(fParameters[0]+fParameters[1]+fParameters[2]+fParameters[3]+fParameters[46])*FunBkgSym1(x,pt,type);
 469         return ret;
 470 }
 471
 472 //_________________________________________________________________________________________________
 473 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFInvMassBkgDistr(Double_t m) const
 474 {
 475         //
 476         // it returns the value of the probability to have a given mass for the background
 477         //
 478         Double_t value = 0.;
 479         if(fExponentialParam)
 480         value = fParameters[14]*TMath::Exp(-1*(m-fParameters[15])/fParameters[16]) + fParameters[17];
 481         else value = fParameters[14] + fParameters[15]*m + fParameters[16]*m*m + fParameters[17]*m*m*m + fParameters[47]*m*m*m*m + fParameters[48]*m*m*m*m*m;
 482         return value;
 483 }
 484 //_________________________________________________________________________________________________
 485 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::FunBkgPos(Double_t x, Double_t pt, Int_t type) const
 486 {
 487         //
 488         // exponential with positive slopes for the background part (x)
 489         //
 490         Double_t np, sc, sigma3;
 491         sc = 10.;
 492         if(fMultivariate){ np = 10000.0; sigma3 = 5000.;}
 493         else{ np = 1000.0; sigma3 = 1000.;}
 494
 495         Double_t xprime;
 496         Double_t sum = 0.0;
 497         Double_t xlow,xupp;
 498         Double_t step;
 499         Double_t i;
 500         xlow = x - sc * sigma3 ;
 501         xupp = x + sc * sigma3 ;
 502         step = (xupp-xlow) / np;
 503
 504         for(i=1.0; i<=np/2; i++) {
 505                 xprime = xlow + (i-.5) * step;
 506                 if (xprime > 0) {sum += fParameters[4] * TMath::Exp(-1*xprime*fParameters[4])*(ResolutionFunc(xprime-x,pt,type));}
 507                 xprime = xupp - (i-.5) * step;
 508                 if (xprime > 0) {sum += fParameters[4] * TMath::Exp(-1*xprime*fParameters[4])*(ResolutionFunc(xprime-x,pt,type));}
 509                 }
 510
 511         return step * sum ;
 512 }
 513 //_________________________________________________________________________________________________
 514 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::FunBkgNeg(Double_t x, Double_t pt, Int_t type) const
 515 {
 516         //
 517         // exponential with negative slopes for the background part (x)
 518         //
 519         Double_t np, sc, sigma3;
 520         sc = 10.;
 521         if(fMultivariate){ np = 10000.0;  sigma3 = 5000.;}
 522         else{ np = 1000.0; sigma3 = 1000.;}
 523
 524         Double_t xprime;
 525         Double_t sum = 0.0;
 526         Double_t xlow,xupp;
 527         Double_t step;
 528         Double_t i;
 529         xlow = x - sc * sigma3 ;
 530         xupp = x + sc * sigma3 ;
 531         step = (xupp-xlow) / np;
 532
 533         for(i=1.0; i<=np/2; i++) {
 534
 535                 xprime = xlow + (i-.5) * step;
 536                 if (xprime < 0) {sum += fParameters[5] * TMath::Exp(xprime*fParameters[5]) * (ResolutionFunc(xprime-x,pt,type));}
 537
 538                 xprime = xupp - (i-.5) * step;
 539                 if (xprime < 0) {sum += fParameters[5] * TMath::Exp(xprime*fParameters[5]) * (ResolutionFunc(xprime-x,pt,type));}
 540         }
 541
 542         return step * sum ;
 543 }
 544 //_________________________________________________________________________________________________
 545 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::FunBkgSym(Double_t x, Double_t pt, Int_t type) const
 546 {
 547         //
 548         // exponential with both positive and negative slopes for the background part (x)
 549         //
 550         Double_t np, sc, sigma3;
 551         sc = 10.;
 552         if(fMultivariate){ np = 10000.0; sigma3 = 5000.;}
 553         else{ np = 1000.0; sigma3 = 1000.;}
 554
 555         Double_t xprime;
 556         Double_t sum1 = 0.0;
 557         Double_t sum2 = 0.0;
 558         Double_t xlow,xupp;
 559         Double_t step;
 560         Double_t i;
 561         xlow = x - sc * sigma3 ;
 562         xupp = x + sc * sigma3 ;
 563         step = (xupp-xlow) / np;
 564
 565         for(i=1.0; i<=np/2; i++) {
 566
 567                 xprime = xlow + (i-.5) * step;
 568                 if (xprime > 0) {sum1 += 0.5 * fParameters[6]*TMath::Exp(-1*xprime*fParameters[6]) * (ResolutionFunc(xprime-x,pt,type));}
 569                 if (xprime < 0) {sum2 += 0.5 * fParameters[6]*TMath::Exp(xprime*fParameters[6]) * (ResolutionFunc(xprime-x,pt,type));}
 570
 571                 xprime = xupp - (i-.5) * step;
 572                 if (xprime > 0) {sum1 += 0.5 * fParameters[6]*TMath::Exp(-1*xprime*fParameters[6]) * (ResolutionFunc(xprime-x,pt,type));}
 573                 if (xprime < 0) {sum2 += 0.5 * fParameters[6]*TMath::Exp(xprime*fParameters[6]) * (ResolutionFunc(xprime-x,pt,type));}
 574         }
 575
 576         return step*(sum1 + sum2) ;
 577 }
 578
 579 //_________________________________________________________________________________________________
 580 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::FunBkgSym1(Double_t x, Double_t pt, Int_t type) const
 581 {
 582         //
 583         // exponential with both positive and negative slopes for the background part (x)
 584         //
 585         Double_t np, sc, sigma3;
 586         sc = 10.;
 587         if(fMultivariate){ np = 10000.0; sigma3 = 5000.;}
 588         else{ np = 1000.0; sigma3 = 1000.;}
 589
 590         Double_t xprime;
 591         Double_t sum1 = 0.0;
 592         Double_t sum2 = 0.0;
 593         Double_t xlow,xupp;
 594         Double_t step;
 595         Double_t i;
 596         xlow = x - sc * sigma3 ;
 597         xupp = x + sc * sigma3 ;
 598         step = (xupp-xlow) / np;
 599
 600         for(i=1.0; i<=np/2; i++) {
 601
 602                 xprime = xlow + (i-.5) * step;
 603                 if (xprime > 0) {sum1 += 0.5 * fParameters[45]*TMath::Exp(-1*xprime*fParameters[45]) * (ResolutionFunc(xprime-x,pt,type));}
 604                 if (xprime < 0) {sum2 += 0.5 * fParameters[45]*TMath::Exp(xprime*fParameters[45]) * (ResolutionFunc(xprime-x,pt,type));}
 605
 606                 xprime = xupp - (i-.5) * step;
 607                 if (xprime > 0) {sum1 += 0.5 * fParameters[45]*TMath::Exp(-1*xprime*fParameters[45]) * (ResolutionFunc(xprime-x,pt,type));}
 608                 if (xprime < 0) {sum2 += 0.5 * fParameters[45]*TMath::Exp(xprime*fParameters[45]) * (ResolutionFunc(xprime-x,pt,type));}
 609         }
 610
 611         return step*(sum1 + sum2) ;
 612 }
 613
 614
 615 //_________________________________________________________________________________________________
 616 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::ResolutionFunc(Double_t x, Double_t pt, Int_t type) const
 617 {
 618         //
 619         // parametrization with 2 gaus
 620         //
 621         x = x/fChangeResolution;
 622         Int_t index = (2-type)*9;
 623         Double_t mean1 = fParameters[20+index];
 624         Double_t mean2 = fParameters[22+index];
 625         Double_t norm1 = fParameters[18+index];
 626         Double_t sigma1 = fParameters[21+index];
 627         Double_t sigma2 = fParameters[23+index];
 628         Double_t norm2 = fParameters[19+index];
 629         //exp values
 630         Double_t alfa = fParameters[24+index];
 631         Double_t lambda = fParameters[25+index];
 632         Double_t norm3 = fParameters[26+index];
 633
 634         if(fMultivariate) // set parameters from matrix
 635         {
 636           //pt dependence
 637           Int_t binPt = -1.;
 638           for(int j=0; j<fPtWindows->GetSize()-1; j++) {if(fPtWindows->At(j)<pt && pt<fPtWindows->At(j+1)) binPt = j;}
 639           norm1 = fResParams[binPt][type][0];
 640           mean1 = fResParams[binPt][type][1];
 641           sigma1 = fResParams[binPt][type][2];
 642           norm2 = fResParams[binPt][type][3];
 643           mean2 = fResParams[binPt][type][4];
 644           sigma2 = fResParams[binPt][type][5];
 645           alfa = fResParams[binPt][type][6];
 646           lambda = fResParams[binPt][type][7];
 647           norm3 = fResParams[binPt][type][8];
 648         }
 649
 650         Double_t ret = 0.; Double_t fitval = 0;
 651         if(TMath::Abs(x)<=alfa) fitval = (lambda-1)/(2*alfa*lambda);
 652         else  fitval = ((lambda-1)/(2*alfa*lambda))*TMath::Power(alfa,lambda)*(TMath::Power(TMath::Abs(x),-1*lambda));
 653
 654         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;
 655
 656         return ret/fChangeResolution;
 657 }
 658
 659
 660 //_________________________________________________________________________________________________
 661 TF1* AliDielectronBtoJPSItoEleCDFfitFCN::GetCsiMC(Double_t xmin, Double_t xmax, Double_t normalization)
 662 {
 663         // return the pointer to the templateMC function
 664         TF1* templateMC = new TF1("MCtemplate",this,&AliDielectronBtoJPSItoEleCDFfitFCN::CsiMCfunc,xmin,xmax,1);
 665         templateMC->SetParameter(0,normalization);
 666         templateMC->SetNpx(5000);
 667         return (TF1*)templateMC->Clone();
 668 }
 669
 670 //__________________________________________________________________________________________________
 671 TF1* AliDielectronBtoJPSItoEleCDFfitFCN::GetResolutionFunc(Double_t xmin, Double_t xmax, Double_t normalization, Double_t pt, Int_t type){
 672         // return the pointer to the resolution function
 673         TF1* resFunc = new TF1(Form("resolutionFunc_%d",type),this,&AliDielectronBtoJPSItoEleCDFfitFCN::ResolutionFuncf,xmin,xmax,3);
 674         resFunc->SetParameter(0,normalization);
 675         resFunc->SetParameter(1,pt);
 676         resFunc->SetParameter(2,(Double_t)type);
 677         resFunc->SetNpx(5000);
 678         return (TF1*)resFunc->Clone();
 679 }
 680
 681 //__________________________________________________________________________________________________
 682 TF1* AliDielectronBtoJPSItoEleCDFfitFCN::GetResolutionFuncAllTypes(Double_t xmin, Double_t xmax, Double_t normalization){
 683         // return the pointer to the resolution function
 684         TF1* resFunc = new TF1(Form("resolutionFunc"),this,&AliDielectronBtoJPSItoEleCDFfitFCN::ResolutionFuncAllTypes,xmin,xmax,1);
 685         resFunc->SetParameter(0,normalization);
 686         resFunc->SetNpx(5000);
 687         return (TF1*)resFunc->Clone();
 688 }
 689
 690 //___________________________________________________________________________________________________
 691 TF1* AliDielectronBtoJPSItoEleCDFfitFCN::GetEvaluateCDFDecayTimeBkgDistr(Double_t xmin, Double_t xmax, Double_t normalization, Int_t type, Double_t mass, Double_t pt, Int_t npx){
 692         // return the pointer to the background x distribution function
 693         TF1 *backFunc = new TF1(Form("backFunc_%d_%1.2f_%1.2f",type,mass,pt),this,&AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeBkgDistrFunc,xmin,xmax,4);
 694         backFunc->SetParameter(0,normalization);
 695         backFunc->SetParameter(1,(Double_t)type);
 696         backFunc->SetParameter(2,mass);
 697         backFunc->SetParameter(3,pt);
 698         backFunc->SetNpx(npx);
 699         return (TF1*)backFunc->Clone();
 700 }
 701
 702 //___________________________________________________________________________________________________
 703 TF1* AliDielectronBtoJPSItoEleCDFfitFCN::GetEvaluateCDFDecayTimeBkgDistrAllTypes(Double_t xmin, Double_t xmax, Double_t normalization){
 704         // return the pointer to the background x distribution function
 705         TF1 *backFuncNew = new TF1(Form("backFunc_%f",normalization),this,&AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeBkgDistrFuncAllTypes,xmin,xmax,1);
 706         backFuncNew->SetParameter(0,normalization);
 707         backFuncNew->SetNpx(5000);
 708         return (TF1*)backFuncNew->Clone();
 709 }
 710
 711 //__________________________________________________________________________________________________
 712 TF1* AliDielectronBtoJPSItoEleCDFfitFCN::GetEvaluateCDFDecayTimeSigDistr(Double_t xmin, Double_t xmax, Double_t normalization, Double_t type){
 713         // return the pointer to the signal x distribution function
 714         TF1 *signFunc = new TF1("signalFunc",this,&AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeSigDistrFunc,xmin,xmax,2);
 715         signFunc->SetParameter(0,normalization);
 716         signFunc->SetParameter(1,type);
 717         signFunc->SetNpx(5000);
 718         return (TF1*)signFunc->Clone();
 719 }
 720
 721 //____________________________________________________________________________________________________
 722 TF1* AliDielectronBtoJPSItoEleCDFfitFCN::GetEvaluateCDFInvMassBkgDistr(Double_t mMin, Double_t mMax, Double_t normalization){
 723   // return the pointer to the invariant mass distribution for the background
 724   TF1 * invMassBkg = new TF1("invMassBkg",this,&AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFInvMassBkgDistrFunc,mMin,mMax,1);
 725   invMassBkg->SetParameter(0,normalization);
 726   invMassBkg->SetNpx(5000);
 727   return (TF1*)invMassBkg->Clone();
 728 }
 729
 730
 731 //____________________________________________________________________________________________________
 732 TF1* AliDielectronBtoJPSItoEleCDFfitFCN::GetEvaluateCDFInvMassSigDistr(Double_t mMin, Double_t mMax, Double_t normalization){
 733   // return the pointer to the invariant mass distribution for the signal
 734   TF1 * invMassSig = new TF1("invMassSig",this,&AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFInvMassSigDistrFunc,mMin,mMax,1);
 735   invMassSig->SetParameter(0,normalization);
 736   invMassSig->SetNpx(5000);
 737   return (TF1*)invMassSig->Clone();
 738 }
 739
 740 //____________________________________________________________________________________________________
 741 TF1 *AliDielectronBtoJPSItoEleCDFfitFCN::GetEvaluateCDFInvMassTotalDistr(Double_t mMin, Double_t mMax, Double_t normalization){
 742   // return the pointer to the invariant mass distribution
 743   TF1 * invMassTot = new TF1("invMassTot",this,&AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFInvMassTotalDistr,mMin,mMax,1);
 744   invMassTot->SetParameter(0,normalization);
 745   invMassTot->SetNpx(5000);
 746   return (TF1*)invMassTot->Clone();
 747 }
 748
 749 //____________________________________________________________________________________________________
 750 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFInvMassTotalDistr(const Double_t* x, const Double_t *par) const
 751 {
 752   // evaluate invariant mass total distribution
 753   Double_t value = 0;
 754   Double_t xx = x[0];
 755   value = ((EvaluateCDFInvMassSigDistr(xx)/fintmMassSig)*fParameters[8] + (1-fParameters[8])*(EvaluateCDFInvMassBkgDistr(xx)/fintmMassBkg))*par[0];
 756   return value;
 757 }
 758
 759 //____________________________________________________________________________________________________
 760 TF1 *AliDielectronBtoJPSItoEleCDFfitFCN::GetEvaluateCDFDecayTimeTotalDistr(Double_t xMin, Double_t xMax, Double_t normalization,Double_t pt, Int_t type){
 761  // return the pointer to the pseudoproper distribution for the background
 762  TF1 *decayTimeTot = new TF1(Form("decayTimeTot_%d",type),this,&AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeTotalDistr,xMin,xMax,3);
 763  decayTimeTot->SetParameter(0,normalization);
 764  decayTimeTot->SetParameter(1,pt);
 765  decayTimeTot->SetParameter(2,(Double_t)type);
 766  decayTimeTot->SetNpx(5000);
 767  return (TF1*)decayTimeTot->Clone();
 768 }
 769
 770 //____________________________________________________________________________________________________
 771 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeTotalDistr(const Double_t* x, const Double_t *par) const
 772 {
 773  // evaluate the total pseudoproper distribution for a given candidate's type. par[1] should be the candidate's type.
 774  Double_t value = 0;
 775  Double_t xx = x[0];
 776  value = (fParameters[8]*EvaluateCDFDecayTimeSigDistr(xx,par[1],(Int_t)par[2]) + (1-fParameters[8])*EvaluateCDFDecayTimeBkgDistr(xx,(Int_t)par[2],par[1]))*par[0];
 777  return value;
 778 }
 779
 780 //____________________________________________________________________________________________________
 781 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeTotalDistrAllTypes(const Double_t* x, const Double_t *par) const
 782 {
 783  // evaluate the total pseudoproper distribution considering all candidate's types
 784  Double_t value = 0;
 785  Double_t xx = x[0];
 786
 787  value = (fParameters[8]*(fWeightType[2]*EvaluateCDFDecayTimeSigDistr(xx,200.,2)+fWeightType[1]*EvaluateCDFDecayTimeSigDistr(xx,200.,1)+fWeightType[0]*EvaluateCDFDecayTimeSigDistr(xx,200.,0)))+((1-fParameters[8])*(fWeightType[2]*EvaluateCDFDecayTimeBkgDistr(xx,2,3.09,200.) + fWeightType[1]*EvaluateCDFDecayTimeBkgDistr(xx,1,3.09,200.)+fWeightType[0]*EvaluateCDFDecayTimeBkgDistr(xx,0,3.09,200.)));
 788
 789  return value*par[0];
 790 }
 791
 792 //____________________________________________________________________________________________________
 793 TF1 *AliDielectronBtoJPSItoEleCDFfitFCN::GetEvaluateCDFDecayTimeTotalDistrAllTypes(Double_t xMin, Double_t xMax, Double_t normalization){
 794  // return the pointer to the pseudoproper function for the background considering all candidate's types
 795  TF1 *decayTimeTot = new TF1(Form("decayTimeTot"),this,&AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeTotalDistrAllTypes,xMin,xMax,1);
 796  decayTimeTot->SetParameter(0,normalization);
 797  decayTimeTot->SetNpx(5000);
 798  return (TF1*)decayTimeTot->Clone();
 799 }
 800
 801
 802 //____________________________________________________________________________________________________
 803 TF1 * AliDielectronBtoJPSItoEleCDFfitFCN::GetFunB(Double_t xmin, Double_t xmax, Double_t normalization, Double_t pt, Int_t type, Int_t npx){
 804  // return the pointer to the function that describe secondary jpsi pseudoproper distribution for a given candidate's type
 805  TF1* funb = new TF1(Form("secondaryJpsiConvolution_%d_%1.2f",type,pt),this,&AliDielectronBtoJPSItoEleCDFfitFCN::FunBfunc,xmin,xmax,3);
 806  funb->SetParameter(0,normalization);
 807  funb->SetParameter(1,pt);
 808  funb->SetParameter(2,(Double_t)type);
 809  funb->SetNpx(npx);
 810  return (TF1*)funb->Clone();
 811  }
 812
 813 //____________________________________________________________________________________________________
 814 TF1 * AliDielectronBtoJPSItoEleCDFfitFCN::GetFunBAllTypes(Double_t xmin, Double_t xmax, Double_t normalization){
 815  // return the pointer to the function that describe secondary jpsi pseudoproper distribution for all types
 816  TF1* funb = new TF1(Form("secondaryJpsiConvolution"),this,&AliDielectronBtoJPSItoEleCDFfitFCN::FunBfuncAllTypes,xmin,xmax,1);
 817  funb->SetParameter(0,normalization);
 818  funb->SetNpx(5000);
 819  return (TF1*)funb->Clone();
 820  }
 821
 822 //
 823 // methods below are needed to perform the multivariate fit (pt,mass,type); this can be enabled
 824 // by the boolean fMultivariate; if functions to describe pseudoproper
 825 // decay lenght in pt,m, type bins have been
 826 // already computed, they can be loaded from the file
 827 // switching-on the option fLoadFunction (this to avoid the
 828 // calculation of convolutions and speed-up the likelihood fit)
 829 //
 830
 831 //____________________________________________________________________________________________________
 832 void AliDielectronBtoJPSItoEleCDFfitFCN::SetBackgroundSpecificParameters(Int_t pt, Int_t mb, Int_t tp){
 833   //
 834   // methods to set specific background parameters in bins(pt, inv. mass, type)
 835   //
 836   for(int j=0; j<4;j++) fParameters[j]=fBkgParams[pt][mb][tp][j];
 837   for(int k=5; k<8;k++) fParameters[k-1]=fBkgParams[pt][mb][tp][k];
 838   fParameters[45] = fBkgParams[pt][mb][tp][8];
 839   fParameters[46] = fBkgParams[pt][mb][tp][4];
 840   return;
 841 }
 842
 843 //_______________________________________________________________________________________________________
 844 void AliDielectronBtoJPSItoEleCDFfitFCN::InitializeFunctions(Int_t ptSize, Int_t massSize){
 845   //
 846   // initialize pointers to save functions for the  multivariate fit
 847   //
 848   fFunBSaved = new TF1**[ptSize];
 849   for(int kpt=0; kpt<ptSize; kpt++) fFunBSaved[kpt] = new TF1*[3]; // type
 850   fFunBkgSaved = new TF1***[ptSize];
 851   for(int kpt=0; kpt<ptSize; kpt++){ fFunBkgSaved[kpt] = new TF1**[massSize];
 852   for(int ks = 0; ks<massSize; ks++) fFunBkgSaved[kpt][ks] = new TF1*[3];
 853   for(int kk=0; kk<3;kk++) {
 854      fFunBSaved[kpt][kk] = new TF1();
 855      for(int kk1=0; kk1<massSize;kk1++){
 856      fFunBkgSaved[kpt][kk1][kk] = new TF1();
 857      }
 858     }
 859   }
 860
 861   // to extrapolate the function under the signal region
 862   fWeights = new Double_t**[massSize - 1]; // mass
 863   for(int km =0; km < (massSize - 1); km++) {fWeights[km] = new Double_t*[ptSize];
 864         for(int kpt =0; kpt<ptSize; kpt++) fWeights[km][kpt] = new Double_t[3];
 865   } // pt
 866  return;
 867 }
 868
 869 //_________________________________________________________________________________________________
 870 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::FunBsaved(Double_t x, Double_t pt, Int_t type) const
 871 {
 872         //
 873         //   functions to describe non-prompt J/psi x distributions
 874         //
 875         Double_t returnvalue = 0.;
 876         Int_t binPt = -1;
 877         for(int j=0; j<fPtWindows->GetSize()-1; j++) {if(fPtWindows->At(j)<pt && pt<fPtWindows->At(j+1)) binPt = j;}
 878         returnvalue = fFunBSaved[binPt][type]->Eval(x);
 879         return returnvalue;
 880 }
 881
 882 //________________________________________________________________________________________________
 883 void AliDielectronBtoJPSItoEleCDFfitFCN::SetFunctionsSaved(Int_t npxFunB, Int_t npxFunBkg, Double_t funBLimits, Double_t funBkgLimits, Int_t signalRegion){
 884         //
 885         // save functions for the multivariate fit
 886         //
 887         if(!fMultivariate)
 888         {AliInfo("Warning: fMultivariate is kFALSE! Functions are not saved! \n"); return;}
 889         SetExtrapolationRegion(signalRegion);
 890
 891                      for(int tp=0;tp<3;tp++)  // type
 892                       {
 893                         // pt
 894                         for(int pt=0; pt<fPtWindows->GetSize()-1;pt++){
 895                         if(fResParams) SetResolutionConstants(fResParams[pt][tp],tp);
 896                         SetFunBFunction(tp,pt,GetFunB(-1.*funBLimits,funBLimits,1.,(fPtWindows->At(pt) + (fPtWindows->At(pt+1)-fPtWindows->At(pt))/2.),tp,npxFunB));
 897                         }
 898                       }
 899
 900         AliInfo("+++++++  Pseudoproper-decay-length function for secondary J/psi saved  ++++++ \n");
 901
 902         if(!fLoadFunctions){
 903         for(int ij = 0; ij<fMassWindows->GetSize()-1;ij++){
 904         if(ij == signalRegion) continue;
 905
 906         Int_t mbin = (ij > signalRegion) ?  ij-1 : ij;
 907         for(int tp=0;tp<3;tp++)  {
 908         for(int pt =0; pt<fPtWindows->GetSize()-1; pt++){
 909          if(fBkgParams) SetBackgroundSpecificParameters(pt,mbin,tp);
 910          SetBkgFunction(ij, tp, pt, GetEvaluateCDFDecayTimeBkgDistr(-1.*funBkgLimits,funBkgLimits,1.,tp,(fMassWindows->At(ij) + (fMassWindows->At(ij+1)-fMassWindows->At(ij))/2.),(fPtWindows->At(pt) + (fPtWindows->At(pt+1)-fPtWindows->At(pt))/2.),npxFunBkg));
 911
 912                 }
 913
 914         }
 915
 916         }
 917         AliInfo("+++++++  Pseudoproper-decay-length function for background saved  +++++++++++ \n");
 918
 919         } // loadFunctions
 920         // evaluate under signal
 921         for(int tp=0;tp<3;tp++)
 922                 {
 923                 for(int pt =0; pt<fPtWindows->GetSize()-1; pt++){
 924                 SetBkgFunction(signalRegion, tp, pt, GetEvaluateCDFDecayTimeBkgDistr(-1.*funBkgLimits,funBkgLimits,1.,tp,(fMassWindows->At(signalRegion) + (fMassWindows->At(signalRegion+1)-fMassWindows->At(signalRegion))/2.),(fPtWindows->At(pt) + (fPtWindows->At(pt+1)-fPtWindows->At(pt))/2.),npxFunBkg));
 925                 }
 926
 927         }
 928         // save functions
 929         TFile func("functions.root","RECREATE");
 930         for(int kpt =0; kpt<fPtWindows->GetSize()-1; kpt++){
 931         for(int ss=0; ss<3;ss++) {fFunBSaved[kpt][ss]->Write();
 932         for(int kk=0; kk<fMassWindows->GetSize()-1; kk++) fFunBkgSaved[kpt][kk][ss]->Write();}}
 933         return;
 934 }
 935
 936 //_________________________________________________________________________________________________
 937 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeBkgDistrDifferential(Double_t x, Int_t type, Double_t m, Double_t pt) const
 938 {
 939         //
 940         // it returns the value of the probability to have a given x for the background
 941         // in the pt, m , type correspondent range
 942         //
 943         Int_t binPt = -1;
 944         for(int j=0; j<fPtWindows->GetSize()-1; j++)
 945         {if(fPtWindows->At(j)<pt && pt<fPtWindows->At(j+1)) binPt = j;}
 946         Bool_t isSignal = (fMassWindows->At(fSignalBinForExtrapolation)<m && m<fMassWindows->At(fSignalBinForExtrapolation+1));
 947         Double_t ret = 0.;
 948         if(!isSignal)
 949         ret = fParameters[0]/(fParameters[0]+fParameters[1]+fParameters[2]+fParameters[3]+fParameters[46])*ResolutionFunc(x, pt, type) + fParameters[1]/(fParameters[0]+fParameters[1]+fParameters[2]+fParameters[3]+fParameters[46])*FunBkgPos(x, pt,type) +  fParameters[2]/(fParameters[0]+fParameters[1]+fParameters[2]+fParameters[3]+fParameters[46])*FunBkgNeg(x,pt,type) + fParameters[3]/(fParameters[0]+fParameters[1]+fParameters[2]+fParameters[3]+fParameters[46])*FunBkgSym(x, pt,type) + fParameters[46]/(fParameters[0]+fParameters[1]+fParameters[2]+fParameters[3]+fParameters[46])*FunBkgSym1(x,pt,type);
 950        else{
 951        for(int k=0; k<fMassWindows->GetSize()-2;k++) {
 952         Int_t mbin = (k > (fSignalBinForExtrapolation-1)) ?  k+1 : k;
 953         ret +=  fWeights[k][binPt][type]*EvaluateCDFDecayTimeBkgDistrSaved(x,type,(fMassWindows->At(mbin) + (fMassWindows->At(mbin+1)-fMassWindows->At(mbin))/2.),pt);}
 954
 955        }
 956         return ret;
 957 }
 958
 959 //_________________________________________________________________________________________________
 960 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeBkgDistrSaved(Double_t x, Int_t type, Double_t m, Double_t pt) const
 961 {
 962         //
 963         // it returns the value of the probability to have a given x for the background
 964         // in the pt, m , type correspondent range
 965         //
 966         Double_t returnvalue = 0.;
 967         Int_t binM = -1.;
 968         for(int j=0; j<fMassWindows->GetSize()-1; j++) {if(fMassWindows->At(j)<m && m<fMassWindows->At(j+1)) binM = j;}
 969         Int_t binPt = -1;
 970         for(int j=0; j<fPtWindows->GetSize()-1; j++) {if(fPtWindows->At(j)<pt && pt<fPtWindows->At(j+1)) binPt = j;}
 971         returnvalue = fFunBkgSaved[binPt][binM][type]->Eval(x);
 972         return returnvalue;
 973 }
 974