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() :
37 fintxDecayTimeBkgPos(1.),
38 fintxDecayTimeBkgNeg(1.),
39 fintxDecayTimeBkgSym(1.),
46 fCrystalBallParam(kFALSE)
51 SetCrystalBallFunction(kFALSE);
54 for(Int_t iPar = 0; iPar < 16; iPar++) fParameters[iPar] = 0.;
55 fParameters[9] = 1.;fParameters[11] = 1.;fParameters[12] = 1.;
56 for(Int_t index=0; index<6; index++) fResolutionConstants[index] = 0.;
57 AliInfo("Instance of AliBtoJPSItoEleCDFfitFCN-class created");
59 //_________________________________________________________________________________________________
60 AliBtoJPSItoEleCDFfitFCN::AliBtoJPSItoEleCDFfitFCN(const AliBtoJPSItoEleCDFfitFCN& source) :
62 fFPlus(source.fFPlus),
63 fFMinus(source.fFMinus),
65 fIntegral(source.fIntegral),
66 fintxFunB(source.fintxFunB),
67 fintxDecayTimeBkgPos(source.fintxDecayTimeBkgPos),
68 fintxDecayTimeBkgNeg(source.fintxDecayTimeBkgNeg),
69 fintxDecayTimeBkgSym(source.fintxDecayTimeBkgSym),
70 fintmMassSig(source.fintmMassSig),
71 fintxRes(source.fintxRes),
72 fintmMassBkg(source.fintmMassBkg),
73 fhCsiMC(source.fhCsiMC),
74 fMassWndHigh(source.fMassWndHigh),
75 fMassWndLow(source.fMassWndLow),
76 fCrystalBallParam(source.fCrystalBallParam)
81 for(Int_t iPar = 0; iPar < 16; iPar++) fParameters[iPar] = source.fParameters[iPar];
82 for(Int_t index=0; index<6; index++) fResolutionConstants[index] = source.fResolutionConstants[index];
84 //_________________________________________________________________________________________________
85 AliBtoJPSItoEleCDFfitFCN& AliBtoJPSItoEleCDFfitFCN::operator=(const AliBtoJPSItoEleCDFfitFCN& source)
88 // Assignment operator
90 if(&source == this) return *this;
91 fFPlus = source.fFPlus;
92 fFMinus = source.fFMinus;
94 fIntegral = source.fIntegral;
95 fintxFunB = source.fintxFunB;
96 fintxDecayTimeBkgPos = source.fintxDecayTimeBkgPos;
97 fintxDecayTimeBkgNeg = source.fintxDecayTimeBkgNeg;
98 fintxDecayTimeBkgSym = source.fintxDecayTimeBkgSym;
99 fintmMassSig = source.fintmMassSig;
100 fintxRes = source.fintxRes;
101 fintmMassBkg = source.fintmMassBkg;
102 fhCsiMC = source.fhCsiMC;
103 fCrystalBallParam = source.fCrystalBallParam;
105 for(Int_t iPar = 0; iPar < 16; iPar++) fParameters[iPar] = source.fParameters[iPar];
106 for(Int_t index=0; index<6; index++) fResolutionConstants[index] = source.fResolutionConstants[index];
110 //_________________________________________________________________________________________________
111 AliBtoJPSItoEleCDFfitFCN::~AliBtoJPSItoEleCDFfitFCN()
114 // Default destructor
118 for(Int_t iPar = 0; iPar < 16; iPar++) fParameters[iPar] = 0.;
119 for(Int_t index=0; index<6; index++) fResolutionConstants[index] = 0.;
121 //_________________________________________________________________________________________________
122 Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateLikelihood(const Double_t* pseudoproperdecaytime,
123 const Double_t* invariantmass, const Int_t ncand)
126 // This function evaluates the Likelihood fnction
127 // It returns the -Log(of the likelihood function)
132 for(Int_t i=0; i < ncand; i++) {
133 f = EvaluateCDFfuncNorm(pseudoproperdecaytime[i],invariantmass[i]);
137 ret+=-1.*TMath::Log(f);
141 //_________________________________________________________________________________________________
142 void AliBtoJPSItoEleCDFfitFCN::SetAllParameters(const Double_t* parameters)
145 // Sets array of FCN parameters
147 for(Int_t index = 0; index < 16; index++) fParameters[index] = parameters[index];
149 //_________________________________________________________________________________________________
150 void AliBtoJPSItoEleCDFfitFCN::ComputeIntegral()
153 // this function compute the integral of the likelihood function
154 // (theoretical function) in order to normalize it to unity
156 Double_t np = 100.0; //number of integration steps
157 Double_t npres = 200.0; //number of integration steps for the resolution function
159 Double_t stepm;Double_t stepx;Double_t stepxres; //integration step width in variable m,x
160 Double_t mx=0.;Double_t xprime=0.;
161 Double_t xlow = -4000.; Double_t xup = 4000.;
162 Double_t i; Double_t j;
163 Double_t sumx = 0.0;Double_t intx = 0.0;Double_t intm = 0.0;
164 stepm = (fMassWndHigh-fMassWndLow)/npm;
165 stepx = (xup-xlow)/np;
166 stepxres = (xup-xlow)/npres;
168 // compute integrals for all the terms
171 Double_t intxRes = 0.0;
172 Double_t sumxRes = 0.0;
173 for(iRes = 1.0; iRes <= npres/2.; iRes++) {
174 xprime = xlow + (iRes - .5)*stepxres;
175 sumxRes += ResolutionFunc(xprime);
176 xprime = xup - (iRes - .5)*stepxres;
177 sumxRes += ResolutionFunc(xprime);
179 intxRes = sumxRes*stepxres;
180 SetIntegralRes(intxRes);
184 Double_t intxFunB = 0.0;
185 Double_t sumxFunB = 0.0;
186 for(iFunB = 1.0; iFunB <= np/2; iFunB++) {
187 xprime = xlow + (iFunB - .5)*stepx;
188 sumxFunB += FunB(xprime);
189 xprime = xup - (iFunB - .5)*stepx;
190 sumxFunB += FunB(xprime);
192 intxFunB = sumxFunB*stepx;
193 SetIntegralFunB(intxFunB);
196 Double_t iDecayTimeBkgPos;
197 Double_t intxDecayTimeBkgPos = 0.0;
198 Double_t sumxDecayTimeBkgPos = 0.0;
199 for(iDecayTimeBkgPos = 1.0; iDecayTimeBkgPos <= np/2; iDecayTimeBkgPos++) {
200 xprime = xlow + (iDecayTimeBkgPos - .5)*stepx;
201 sumxDecayTimeBkgPos += FunBkgPos(xprime);
202 xprime = xup - (iDecayTimeBkgPos - .5)*stepx;
203 sumxDecayTimeBkgPos += FunBkgPos(xprime);
205 intxDecayTimeBkgPos = sumxDecayTimeBkgPos*stepx;
206 SetIntegralBkgPos(intxDecayTimeBkgPos);
209 Double_t iDecayTimeBkgNeg;
210 Double_t intxDecayTimeBkgNeg = 0.0;
211 Double_t sumxDecayTimeBkgNeg = 0.0;
212 for(iDecayTimeBkgNeg = 1.0; iDecayTimeBkgNeg<= np/2; iDecayTimeBkgNeg++) {
213 xprime = xlow + (iDecayTimeBkgNeg - .5)*stepx;
214 sumxDecayTimeBkgNeg += FunBkgNeg(xprime);
215 xprime = xup - (iDecayTimeBkgNeg - .5)*stepx;
216 sumxDecayTimeBkgNeg += FunBkgNeg(xprime);
218 intxDecayTimeBkgNeg = sumxDecayTimeBkgNeg*stepx;
219 SetIntegralBkgNeg(intxDecayTimeBkgNeg);
221 Double_t iDecayTimeBkgSym;
222 Double_t intxDecayTimeBkgSym = 0.0;
223 Double_t sumxDecayTimeBkgSym = 0.0;
224 for(iDecayTimeBkgSym = 1.0; intxDecayTimeBkgSym <= np/2; intxDecayTimeBkgSym++) {
225 xprime = xlow + (intxDecayTimeBkgSym - .5)*stepx;
226 sumxDecayTimeBkgSym += FunBkgSym(xprime);
227 xprime = xup - (intxDecayTimeBkgSym - .5)*stepx;
228 sumxDecayTimeBkgSym += FunBkgSym(xprime);
230 intxDecayTimeBkgSym = sumxDecayTimeBkgSym*stepx;
231 SetIntegralBkgSym(intxDecayTimeBkgSym);
235 Double_t intmMassSig = 0.0;
236 Double_t summMassSig = 0.0;
237 for(iMassSig = 1.0; iMassSig<= npm/2.; iMassSig++) {
238 mx = fMassWndLow + (iMassSig - .5)*stepm;
239 summMassSig += EvaluateCDFInvMassSigDistr(mx);
240 mx = fMassWndHigh - (iMassSig - .5)*stepm;
241 summMassSig += EvaluateCDFInvMassSigDistr(mx);
243 intmMassSig = summMassSig*stepm;
244 SetIntegralMassSig(intmMassSig);
248 Double_t intmMassBkg = 0.0;
249 Double_t summMassBkg = 0.0;
250 for(iMassBkg = 1.0; iMassBkg <= npm/2.; iMassBkg++) {
251 mx = fMassWndLow + (iMassBkg - .5)*stepm;
252 summMassBkg += EvaluateCDFInvMassBkgDistr(mx);
253 mx = fMassWndHigh - (iMassBkg - .5)*stepm;
254 summMassBkg += EvaluateCDFInvMassBkgDistr(mx);
256 intmMassBkg = summMassBkg*stepm;
257 SetIntegralMassBkg(intmMassBkg);
260 // Compute integral of the whole distribution function
262 for(i = 1.0; i <= np; i++) {
264 xprime = xlow + (i - .5)*stepx;
265 for(j = 1.0; j <= npm/2; j++) {
266 mx = fMassWndLow + (j - .5)*stepm;
267 summ += EvaluateCDFfunc(xprime,mx);
268 mx = fMassWndHigh - (j - .5)*stepm;
269 summ += EvaluateCDFfunc(xprime,mx);
278 //_________________________________________________________________________________________________
279 void AliBtoJPSItoEleCDFfitFCN::PrintStatus()
282 // Print the parameters of the fits
285 printf("actual value of fRadius---------------------------------------->> | %f \n", GetRadius());
286 printf("actual value of fTheta ---------------------------------------->> | %f \n", GetTheta());
287 printf("actual value of fPhi ------------------------------------------>> | %f \n", GetPhi());
288 printf("actual value of fFPlus ---------------------------------------->> | %f \n", GetFPlus());
289 printf("actual value of fFMinus --------------------------------------->> | %f \n", GetFMinus());
290 printf("actual value of fFSym ----------------------------------------->> | %f \n", GetFSym());
291 printf("actual value of fOneOvLamPlus --------------------------------->> | %f \n", GetLamPlus());
292 printf("actual value of fOneOvLamMinus -------------------------------->> | %f \n", GetLamMinus());
293 printf("actual value of fOneOvLamSym ---------------------------------->> | %f \n", GetLamSym());
294 printf("actual value of fMassBkgSlope --------------------------------->> | %f \n", GetMassSlope());
295 printf("actual value of fFractionJpsiFromBeauty ----------------------->> | %f \n", GetFractionJpsiFromBeauty());
296 printf("actual value of fFsig ----------------------------------------->> | %f \n", GetFsig());
297 if(fCrystalBallParam){
298 printf("actual value of fCrystalBallMmean ----------------------------->> | %f \n", GetCrystalBallMmean());
299 printf("actual value of fCrystalBallNexp ------------------------------>> | %f \n", GetCrystalBallNexp());
300 printf("actual value of fCrystalBallSigma ----------------------------->> | %f \n", GetCrystalBallSigma());
301 printf("actual value of fCrystalBallAlpha ----------------------------->> | %f \n", GetCrystalBallAlpha());
302 printf("actual value of fCrystalBallNorm ----------------------------->> | %f \n", GetCrystalBallNorm());
304 printf("actual value of fMpv ------------------------------------------>> | %f \n", GetCrystalBallMmean());
305 printf("actual value of fConstRovL ------------------------------------>> | %f \n", GetCrystalBallNexp());
306 printf("actual value of fSigmaL --------------------------------------->> | %f \n", GetCrystalBallSigma());
307 printf("actual value of fSigmaR --------------------------------------->> | %f \n", GetCrystalBallAlpha());
309 printf("actual value of fSigmaResol ----------------------------------->> | %f \n", GetSigmaResol());
310 printf("actual value of fNResol --------------------------------------->> | %f \n", GetNResol());
312 printf("Actual value of normalization integral for FunB ------------------->> | %f \n", GetIntegralFunB());
313 printf("Actual value of normalization integral for BkgPos ----------------->> | %f \n", GetIntegralBkgPos());
314 printf("Actual value of normalization integral for BkgNeg ----------------->> | %f \n", GetIntegralBkgNeg());
315 printf("Actual value of normalization integral for BkgSym ----------------->> | %f \n", GetIntegralBkgSym());
316 printf("Actual value of normalization integral for MassSig ---------------->> | %f \n", GetIntegralMassSig());
317 printf("Actual value of normalization integral for MassBkg ---------------->> | %f \n", GetIntegralMassBkg());
318 printf("Actual value of normalization integral for Resolution ------------->> | %f \n", GetIntegralRes());
319 printf("Actual value of normalization integral for FCN -------------------->> | %f \n", GetIntegral());
323 //_________________________________________________________________________________________________
324 void AliBtoJPSItoEleCDFfitFCN::SetResolutionConstants()
327 // This method must be update:
328 // for the time beeing the values are hard-wired.
329 // Implementations have to be done to set the values from outside
330 // (e.g. from a ConfigHF file) starting from an indipendent fit
331 // of primary JPSI distribution.
334 fResolutionConstants[0] = 8.; // mean sigma2/sigma1
335 fResolutionConstants[1] = 0.1675; // mean Integral2/Integral1
336 fResolutionConstants[2] = 1374.; // sigma2
337 fResolutionConstants[3] = 0.001022; // N2
338 fResolutionConstants[4] = 686.6; // mu2
340 //_________________________________________________________________________________________________
341 Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFfunc(Double_t x, Double_t m) const
343 return fParameters[8]*EvaluateCDFfuncSignalPart(x,m) + (1. - fParameters[8])*EvaluateCDFfuncBkgPart(x,m);
345 //_________________________________________________________________________________________________
346 Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFfuncNorm(Double_t x, Double_t m) const
348 return EvaluateCDFfunc(x,m)/fIntegral;
350 //_________________________________________________________________________________________________
351 Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFfuncSignalPart(Double_t x, Double_t m) const
353 return EvaluateCDFDecayTimeSigDistr(x)*(EvaluateCDFInvMassSigDistr(m)/fintmMassSig);
355 //_________________________________________________________________________________________________
356 Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeSigDistr(Double_t x) const
359 // Implementation of the Background part of the Likelyhood function
362 Double_t retvalue = 0.;
363 Double_t FunBnorm = FunB(x)/fintxFunB;
364 Double_t FunPnorm = ResolutionFunc(x)/fintxRes;
365 retvalue = fParameters[7]*FunBnorm + (1. - fParameters[7])*FunPnorm;
368 //_________________________________________________________________________________________________
369 Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFInvMassSigDistr(Double_t m) const
372 // Parametrization of signal part invariant mass distribution
373 // It can be either Crystal Ball function or sum of two Landau
376 Double_t fitval = 0.;
378 if(fCrystalBallParam){
379 Double_t t = (m-fParameters[9])/fParameters[11]; ;
380 if (fParameters[12] < 0) t = -t;
382 Double_t absAlpha = TMath::Abs((Double_t)fParameters[12]);
384 if (t >= -absAlpha) {
385 return fParameters[13]*TMath::Exp(-0.5*t*t);
388 Double_t a = TMath::Power(fParameters[10]/absAlpha,fParameters[10])* TMath::Exp(-0.5*absAlpha*absAlpha);
389 Double_t b= fParameters[10]/absAlpha - absAlpha;
390 fitval = (fParameters[13]*a/TMath::Power(b - t, fParameters[10]));
395 Double_t tmpv=-1*fParameters[9];
396 fitval=TMath::Sqrt(TMath::Landau(t,tmpv,fParameters[11]));
397 fitval += fParameters[10]*(TMath::Landau(m,fParameters[9],fParameters[12]));
401 //_________________________________________________________________________________________________
402 Double_t AliBtoJPSItoEleCDFfitFCN::FunB(Double_t x) const
405 // Parameterisation of the fit function for the x part of the Background
410 Double_t sigma3 = fResolutionConstants[2];
416 xlow = x - sc * sigma3 ;
417 xupp = x + sc * sigma3 ;
418 step = (xupp-xlow) / np;
419 Double_t CsiMCxprime = 0.;
420 Double_t Resolutionxdiff = 0.;
423 for(i=1.0; i<=np; i++) {
425 xprime = xlow + (i-.5) * step;
426 CsiMCxprime = CsiMC(xprime);
428 Resolutionxdiff = ResolutionFunc(xdiff)/fintxRes; // normalized value
429 sum += CsiMCxprime * Resolutionxdiff;
435 //_________________________________________________________________________________________________
436 Double_t AliBtoJPSItoEleCDFfitFCN::FunP(Double_t x) const
439 // Parameterisation of the Prompt part for the x distribution
442 return ResolutionFunc(x)/fintxRes; // normalized value
444 //_________________________________________________________________________________________________
445 Double_t AliBtoJPSItoEleCDFfitFCN::CsiMC(Double_t x) const
448 // Distribution (template) of the x distribution for the x variable
449 // for the J/psi coming from Beauty hadrons
452 Double_t returnvalue = 0.;
453 returnvalue = fhCsiMC->GetBinContent(fhCsiMC->FindBin(x));
456 //_________________________________________________________________________________________________
457 Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFfuncBkgPart(Double_t x,Double_t m) const
460 // Return the part of the likelihood function for the background hypothesis
463 return EvaluateCDFDecayTimeBkgDistr(x)*(EvaluateCDFInvMassBkgDistr(m)/fintmMassBkg);
465 //_________________________________________________________________________________________________
466 Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeBkgDistr(Double_t x) const
469 // it returns the value of the probability to have a given x for the background
472 Double_t ret = (1. - TMath::Power(fParameters[0],2.))*(ResolutionFunc(x)/fintxRes)
473 + TMath::Power(fParameters[0]*TMath::Cos(fParameters[1]),2.)*
474 (FunBkgPos(x)/fintxDecayTimeBkgPos)
475 + TMath::Power(fParameters[0]*TMath::Sin(fParameters[1])*TMath::Sin(fParameters[2]),2.)*
476 (FunBkgNeg(x)/fintxDecayTimeBkgNeg)
477 + TMath::Power(fParameters[0]*TMath::Sin(fParameters[1])*TMath::Cos(fParameters[2]),2.)*
478 (FunBkgSym(x)/fintxDecayTimeBkgSym);
481 //_________________________________________________________________________________________________
482 Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFInvMassBkgDistr(Double_t m) const
485 // it returns the value of the probability to have a given mass for the background
488 return 1/(fMassWndHigh-fMassWndLow) +
490 fParameters[6] * ((fMassWndHigh+fMassWndLow)/2);
492 //_________________________________________________________________________________________________
493 Double_t AliBtoJPSItoEleCDFfitFCN::FunBkgPos(Double_t x) const
496 // exponential with positive slopes for the background part (x)
501 Double_t sigma3 = fResolutionConstants[2];
507 xlow = x - sc * sigma3 ;
508 xupp = x + sc * sigma3 ;
509 step = (xupp-xlow) / np;
511 for(i=1.0; i<=np/2; i++) {
513 xprime = xlow + (i-.5) * step;
514 if (xprime > 0) {sum += fParameters[3] * TMath::Exp(-1*xprime*fParameters[3]) * (ResolutionFunc(xprime-x)/fintxRes);} else {sum = 0.;}
516 xprime = xupp - (i-.5) * step;
517 if (xprime > 0) {sum += fParameters[3] * TMath::Exp(-1*xprime*fParameters[3]) * (ResolutionFunc(xprime-x)/fintxRes);} else {sum = 0.;}
523 //_________________________________________________________________________________________________
524 Double_t AliBtoJPSItoEleCDFfitFCN::FunBkgNeg(Double_t x) const
527 // exponential with negative slopes for the background part (x)
532 Double_t sigma3 = fResolutionConstants[2];
538 xlow = x - sc * sigma3 ;
539 xupp = x + sc * sigma3 ;
540 step = (xupp-xlow) / np;
542 for(i=1.0; i<=np/2; i++) {
544 xprime = xlow + (i-.5) * step;
545 if (xprime < 0) {sum += fParameters[4] * TMath::Exp(xprime*fParameters[4]) * (ResolutionFunc(xprime-x)/fintxRes);} else {sum = 0.;}
547 xprime = xupp - (i-.5) * step;
548 if (xprime < 0) {sum += fParameters[4] * TMath::Exp(xprime*fParameters[4]) * (ResolutionFunc(xprime-x)/fintxRes);} else {sum = 0.;}
553 //_________________________________________________________________________________________________
554 Double_t AliBtoJPSItoEleCDFfitFCN::FunBkgSym(Double_t x) const
557 // exponential with both positive and negative slopes for the background part (x)
562 Double_t sigma3 = fResolutionConstants[2];
569 xlow = x - sc * sigma3 ;
570 xupp = x + sc * sigma3 ;
571 step = (xupp-xlow) / np;
573 for(i=1.0; i<=np/2; i++) {
575 xprime = xlow + (i-.5) * step;
576 if (xprime > 0) {sum1 += 0.5 * fParameters[5]*TMath::Exp(-1*xprime*fParameters[5]) * (ResolutionFunc(xprime-x)/fintxRes);} else {sum1 = 0.;}
578 xprime = xupp - (i-.5) * step;
579 if (xprime > 0) {sum1 += 0.5 * fParameters[5]*TMath::Exp(-1*xprime*fParameters[5]) * (ResolutionFunc(xprime-x)/fintxRes);} else {sum1 = 0.;}
582 for(i=1.0; i<=np/2; i++) {
584 xprime = xlow + (i-.5) * step;
585 if (xprime < 0) {sum2 += 0.5 * fParameters[5]*TMath::Exp(xprime*fParameters[5]) * (ResolutionFunc(xprime-x)/fintxRes);} else {sum2 = 0.;}
587 xprime = xupp - (i-.5) * step;
588 if (xprime < 0) {sum2 += 0.5 * fParameters[5]*TMath::Exp(xprime*fParameters[5]) * (ResolutionFunc(xprime-x)/fintxRes);} else {sum2 = 0.;}
591 return step*(sum1 + sum2) ;
593 //_________________________________________________________________________________________________
594 Double_t AliBtoJPSItoEleCDFfitFCN::ResolutionFunc(Double_t x) const
597 // parametrization with 2 gaus
603 //Double_t mean1 = 0.;
604 Double_t mean2 = fResolutionConstants[4];
605 Double_t sigma1 = fParameters[14];
606 Double_t sigma2 = fResolutionConstants[2];
607 Double_t n1 = fParameters[15];
608 Double_t n2 = fResolutionConstants[3];
609 Double_t arg1 = x1/sigma1;
610 Double_t arg2 = (x2-mean2)/sigma2;
611 Double_t sqrt2Pi = TMath::Sqrt(2*TMath::Pi());
613 ret = n2*((n1/n2)*TMath::Exp(-0.5*arg1*arg1) + TMath::Exp(-0.5*arg2*arg2));
615 return ret/(sqrt2Pi*(n1*sigma1+n2*sigma2));//return value is normalized