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