1 /**************************************************************************
2 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 #include "AliBtoJPSItoEleCDFfitFCN.h"
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
24 // Origin: C.Di Giglio
25 // Contact: Carmelo.Digiglio@ba.infn.it , Giuseppe.Bruno@ba.infn.it
26 //_________________________________________________________________________
28 ClassImp(AliBtoJPSItoEleCDFfitFCN)
30 //_________________________________________________________________________________________________
31 AliBtoJPSItoEleCDFfitFCN::AliBtoJPSItoEleCDFfitFCN() :
39 fCrystalBallParam(kFALSE)
44 SetCrystalBallParam(kFALSE);
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");
52 //_________________________________________________________________________________________________
53 AliBtoJPSItoEleCDFfitFCN::AliBtoJPSItoEleCDFfitFCN(const AliBtoJPSItoEleCDFfitFCN& source) :
55 fFPlus(source.fFPlus),
56 fFMinus(source.fFMinus),
58 fIntegral(source.fIntegral),
59 fhCsiMC(source.fhCsiMC),
60 fMassWndHigh(source.fMassWndHigh),
61 fMassWndLow(source.fMassWndLow),
62 fCrystalBallParam(source.fCrystalBallParam)
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];
70 //_________________________________________________________________________________________________
71 AliBtoJPSItoEleCDFfitFCN& AliBtoJPSItoEleCDFfitFCN::operator=(const AliBtoJPSItoEleCDFfitFCN& source)
74 // Assignment operator
76 if(&source == this) return *this;
77 fFPlus = source.fFPlus;
78 fFMinus = source.fFMinus;
80 fIntegral = source.fIntegral;
81 fhCsiMC = source.fhCsiMC;
82 fCrystalBallParam = source.fCrystalBallParam;
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];
89 //_________________________________________________________________________________________________
90 AliBtoJPSItoEleCDFfitFCN::~AliBtoJPSItoEleCDFfitFCN()
97 for(Int_t iPar = 0; iPar < 13; iPar++) fParameters[iPar] = 0.;
98 for(Int_t index=0; index<7; index++) fResolutionConstants[index] = 0.;
100 //_________________________________________________________________________________________________
101 Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateLikelihood(const Double_t* pseudoproperdecaytime,
102 const Double_t* invariantmass, const Int_t ncand)
105 // This function evaluates the Likelihood fnction
106 // It returns the -Log(of the likelihood function)
111 for(Int_t i=0; i < ncand; i++) {
112 f = EvaluateCDFfuncNorm(pseudoproperdecaytime[i],invariantmass[i]);
114 //AliWarning("One negative contributors in the Log(Likely) ! ");
117 ret+=-1.*TMath::Log(f);
121 //_________________________________________________________________________________________________
122 void AliBtoJPSItoEleCDFfitFCN::SetAllParameters(const Double_t* parameters)
125 // Sets array of FCN parameters
127 for(Int_t index = 0; index < 13; index++) fParameters[index] = parameters[index];
129 //_________________________________________________________________________________________________
130 void AliBtoJPSItoEleCDFfitFCN::ComputeIntegral()
133 // this function compute the integral of the likelihood function
134 // (theoretical function) in order to normalize it to unity
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;
145 for(i = 1.0; i <= np; i++) {
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);
160 //_________________________________________________________________________________________________
161 void AliBtoJPSItoEleCDFfitFCN::PrintStatus()
164 // Print the parameters of the fits
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());
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());
191 printf("Actual value of normalization integral for FCN ---------------->> | %f \n", GetIntegral());
194 //_________________________________________________________________________________________________
195 void AliBtoJPSItoEleCDFfitFCN::SetResolutionConstants(Int_t BinNum)
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)
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;
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;
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;
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;
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;
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;
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;
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;
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;
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;
245 //_________________________________________________________________________________________________
246 Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFfunc(Double_t x, Double_t m) const
248 return fParameters[8]*EvaluateCDFfuncSignalPart(x,m) + (1. - fParameters[8])*EvaluateCDFfuncBkgPart(x,m);
250 //_________________________________________________________________________________________________
251 Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFfuncNorm(Double_t x, Double_t m) const
253 return EvaluateCDFfunc(x,m)/fIntegral;
255 //_________________________________________________________________________________________________
256 Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFfuncSignalPart(Double_t x, Double_t m) const
258 return EvaluateCDFDecayTimeSigDistr(x)*EvaluateCDFInvMassSigDistr(m);
260 //_________________________________________________________________________________________________
261 Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeSigDistr(Double_t x) const
264 // Implementation of the Background part of the Likelyhood function
267 Double_t retvalue = 0.;
268 retvalue = fParameters[7]*FunB(x) + (1. - fParameters[7])*FunP(x);
271 //_________________________________________________________________________________________________
272 Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFInvMassSigDistr(Double_t m) const
275 // Parametrization of signal part invariant mass distribution
276 // It can be either Crystal Ball function or sum of two Landau
278 Double_t fitval = 0.;
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;
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;
293 if(arg > -1*TMath::Abs(fParameters[12])){
294 fitvalCB += TMath::Exp(-0.5*arg*arg);
296 fitval = normFactorCB*fitvalCB;
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]));
306 //_________________________________________________________________________________________________
307 Double_t AliBtoJPSItoEleCDFfitFCN::FunB(Double_t x) const
310 // Parameterisation of the fit function for the x part of the Background
314 Double_t sigma3 = fResolutionConstants[2];
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);
327 xx = xupp - (i-.5) * step;
328 sum += CsiMC(xx) * ResolutionFunc(xx-x);
330 return (step * sum) ;
332 //_________________________________________________________________________________________________
333 Double_t AliBtoJPSItoEleCDFfitFCN::FunP(Double_t x) const
335 // Parameterisation of the Prompt part for the x distribution
338 return ResolutionFunc(x);
340 //_________________________________________________________________________________________________
341 Double_t AliBtoJPSItoEleCDFfitFCN::CsiMC(Double_t x) const
343 // Distribution (template) of the x distribution for the x variable
344 // for the J/psi coming from Beauty hadrons
346 Double_t returnvalue = 0.;
347 returnvalue = fhCsiMC->GetBinContent(fhCsiMC->FindBin(x));
350 //_________________________________________________________________________________________________
351 Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFfuncBkgPart(Double_t x,Double_t m) const
354 // Return the part of the likelihood function for the background hypothesis
356 return EvaluateCDFDecayTimeBkgDistr(x)*EvaluateCDFInvMassBkgDistr(m);
358 //_________________________________________________________________________________________________
359 Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeBkgDistr(Double_t x) const
361 // it returns the value of the probability to have a given x for the background
364 Double_t ret = (1 - fParameters[0]*fParameters[0])*ResolutionFunc(x);
370 //_________________________________________________________________________________________________
371 Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFInvMassBkgDistr(Double_t m) const
374 // it returns the value of the probability to have a given mass for the background
376 return 1/(fMassWndHigh-fMassWndLow)+fParameters[6]*(m-(fMassWndHigh+fMassWndLow)/2);
378 //_________________________________________________________________________________________________
379 Double_t AliBtoJPSItoEleCDFfitFCN::FunBkgPos(Double_t x) const
382 // exponential with positive slopes for the background part (x)
386 Double_t sigma3 = fResolutionConstants[2];
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);
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);
402 return (step * sum) ;
404 //_________________________________________________________________________________________________
405 Double_t AliBtoJPSItoEleCDFfitFCN::FunBkgNeg(Double_t x) const
408 // exponential with negative slopes for the background part (x)
412 Double_t sigma3 = fResolutionConstants[2];
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);
429 return (step * sum) ;
431 //_________________________________________________________________________________________________
432 Double_t AliBtoJPSItoEleCDFfitFCN::FunBkgSym(Double_t x) const
435 // exponential with both positive and negative slopes for the background part (x)
439 Double_t sigma3 = fResolutionConstants[2];
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);
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);
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);
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);
467 return step*(sum1 + sum2) ;
469 //_________________________________________________________________________________________________
470 Double_t AliBtoJPSItoEleCDFfitFCN::ResolutionFunc(Double_t x) const
473 //parametrization with 2 gaus + 1 exponential + 1 constant
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]);