72ff24cbe89d729242f679688b564a9f3aba5b42
[u/mrichter/AliRoot.git] / 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 "AliDielectronBtoJPSItoEleCDFfitFCN.h"
17 #include "TFormula.h"
18 #include "TF1.h"
19 #include "TCanvas.h"
20 #include "TMath.h"
21
22 //_________________________________________________________________________
23 //                        Class AliDielectronBtoJPSItoEleCDFfitFCN
24 //                   Definition of main function used in 
25 //                     unbinned log-likelihood fit for
26 //                 the channel B -> JPsi + X -> e+e- + X
27 //      
28 //                           Origin: C.Di Giglio
29 //       Contact: Carmelo.Digiglio@ba.infn.it , Giuseppe.Bruno@ba.infn.it
30 //_________________________________________________________________________
31
32 ClassImp(AliDielectronBtoJPSItoEleCDFfitFCN)
33
34         //_________________________________________________________________________________________________
35         AliDielectronBtoJPSItoEleCDFfitFCN::AliDielectronBtoJPSItoEleCDFfitFCN() :
36                 fFPlus(0.),
37                 fFMinus(0.),
38                 fFSym(0.),
39                 fintmMassSig(1.),
40                 fintmMassBkg(1.),
41                 fhCsiMC(0x0),
42                 fMassWndHigh(0.),
43                 fMassWndLow(0.),
44                 fCrystalBallParam(kFALSE)
45 {
46         //
47         // constructor
48         //
49         SetCrystalBallFunction(kFALSE);
50         SetMassWndHigh(0.2);
51         SetMassWndLow(0.5);
52         fWeightType[0] = 1.; fWeightType[1] = 1.; fWeightType[2] = 1.;
53         for(Int_t iPar = 0; iPar < 45; iPar++) fParameters[iPar] = 0.;
54         fParameters[9] = 1.;fParameters[11] = 1.;fParameters[12] = 1.;
55         AliInfo("Instance of AliDielectronBtoJPSItoEleCDFfitFCN-class created");
56 }
57 //_________________________________________________________________________________________________
58 AliDielectronBtoJPSItoEleCDFfitFCN::AliDielectronBtoJPSItoEleCDFfitFCN(const AliDielectronBtoJPSItoEleCDFfitFCN& source) :
59         TNamed(source),
60         fFPlus(source.fFPlus),
61         fFMinus(source.fFMinus),
62         fFSym(source.fFSym),
63         fintmMassSig(source.fintmMassSig),
64         fintmMassBkg(source.fintmMassBkg),
65         fhCsiMC(source.fhCsiMC),
66         fMassWndHigh(source.fMassWndHigh),
67         fMassWndLow(source.fMassWndLow),
68         fCrystalBallParam(source.fCrystalBallParam) 
69 {
70         //
71         // Copy constructor
72         //
73         for(Int_t iPar = 0; iPar < 45; iPar++) fParameters[iPar] = source.fParameters[iPar];
74         for(Int_t iW=0; iW<2; iW++) fWeightType[iW] = source.fWeightType[iW]; 
75 }
76 //_________________________________________________________________________________________________
77 AliDielectronBtoJPSItoEleCDFfitFCN& AliDielectronBtoJPSItoEleCDFfitFCN::operator=(const AliDielectronBtoJPSItoEleCDFfitFCN& source) 
78 {
79         //
80         // Assignment operator
81         //
82         if(&source == this) return *this;
83         fFPlus = source.fFPlus;
84         fFMinus = source.fFMinus;
85         fFSym = source.fFSym;
86         fintmMassSig = source.fintmMassSig;
87         fintmMassBkg = source.fintmMassBkg;
88         fhCsiMC = source.fhCsiMC;
89         fCrystalBallParam = source.fCrystalBallParam;
90
91         for(Int_t iPar = 0; iPar < 45; iPar++) fParameters[iPar] = source.fParameters[iPar];
92         return *this;
93 }  
94 //_________________________________________________________________________________________________
95 AliDielectronBtoJPSItoEleCDFfitFCN::~AliDielectronBtoJPSItoEleCDFfitFCN()
96 {
97         //
98         // Default destructor
99         //
100
101         delete fhCsiMC;
102         for(Int_t iPar = 0; iPar < 45; iPar++) fParameters[iPar] = 0.;
103 }
104 //_________________________________________________________________________________________________
105 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateLikelihood(const Double_t* pseudoproperdecaytime,
106                 const Double_t* invariantmass, const Int_t *type, const Int_t ncand) const
107 {
108         //
109         // This function evaluates the Likelihood fnction
110         // It returns the -Log(of the likelihood function)
111         //
112         Double_t f = 0.;
113         Double_t ret = 0.;
114
115         for(Int_t i=0; i < ncand; i++) {
116                 f = EvaluateCDFfuncNorm(pseudoproperdecaytime[i],invariantmass[i],type[i]);
117                 if(f <= 0.) continue;  
118                 ret += -1*TMath::Log(f);  
119         }
120         return ret;
121 }
122 //_________________________________________________________________________________________________
123 void AliDielectronBtoJPSItoEleCDFfitFCN::SetAllParameters(const Double_t* parameters)
124
125         //
126         // Sets array of FCN parameters
127         //
128         for(Int_t index = 0; index < 45; index++) fParameters[index] = parameters[index];
129 }
130 //_________________________________________________________________________________________________
131 void AliDielectronBtoJPSItoEleCDFfitFCN::ComputeMassIntegral() 
132
133         //
134         // this function compute the integral of the likelihood function 
135         // (theoretical function) in order to normalize it to unity
136         //
137         Double_t npm = 20000.;
138         Double_t stepm;
139         Double_t mx=0.;
140         stepm = (fMassWndHigh-fMassWndLow)/npm; 
141         // compute integrals for  invariant mass terms        
142
143         Double_t iMassSig;
144         Double_t intmMassSig = 0.0;
145         Double_t summMassSig = 0.0;
146         for(iMassSig = 1.0;  iMassSig<= npm/2.; iMassSig++)  {
147                 mx = fMassWndLow + (iMassSig - .5)*stepm;
148                 summMassSig += EvaluateCDFInvMassSigDistr(mx);
149                 mx = fMassWndHigh - (iMassSig - .5)*stepm;
150                 summMassSig += EvaluateCDFInvMassSigDistr(mx);
151         }
152         intmMassSig = summMassSig*stepm;
153         SetIntegralMassSig(intmMassSig);
154         //
155         Double_t iMassBkg;
156         Double_t intmMassBkg = 0.0;
157         Double_t summMassBkg = 0.0;
158         for(iMassBkg = 1.0; iMassBkg <= npm/2.; iMassBkg++)  {
159                 mx = fMassWndLow + (iMassBkg - .5)*stepm;
160                 summMassBkg += EvaluateCDFInvMassBkgDistr(mx);
161                 mx = fMassWndHigh - (iMassBkg - .5)*stepm;
162                 summMassBkg += EvaluateCDFInvMassBkgDistr(mx);
163         }
164         intmMassBkg = summMassBkg*stepm;
165         if(intmMassBkg < 1.e-05) intmMassBkg = 1.;
166         SetIntegralMassBkg(intmMassBkg);
167 }
168 //_________________________________________________________________________________________________
169 void AliDielectronBtoJPSItoEleCDFfitFCN::PrintStatus()
170 {
171         //
172         //  Print the parameters of the fits 
173         //
174         printf("\n");
175         // background param
176         printf("actual value of fWeightRes------------------------------------->> | %f \n", GetResWeight());
177         printf("actual value of fPos ------------------------------------------>> | %f \n", GetFPlus());
178         printf("actual value of fNeg ------------------------------------------>> | %f \n", GetFMinus());
179         printf("actual value of fSym ------------------------------------------>> | %f \n", GetFSym()); 
180         printf("actual value of fOneOvLamPlus --------------------------------->> | %f \n", GetLamPlus());
181         printf("actual value of fOneOvLamMinus -------------------------------->> | %f \n", GetLamMinus());
182         printf("actual value of fOneOvLamSym ---------------------------------->> | %f \n", GetLamSym());
183         printf("actual value of fFractionJpsiFromBeauty ----------------------->> | %f \n", GetFractionJpsiFromBeauty());
184         printf("actual value of fFsig ----------------------------------------->> | %f \n", GetFsig());
185
186         if(fCrystalBallParam){
187                 printf("actual value of fCrystalBallMmean ----------------------------->> | %f \n", GetCrystalBallMmean());
188                 printf("actual value of fCrystalBallNexp ------------------------------>> | %f \n", GetCrystalBallNexp());
189                 printf("actual value of fCrystalBallSigma ----------------------------->> | %f \n", GetCrystalBallSigma());
190                 printf("actual value of fCrystalBallAlpha ----------------------------->> | %f \n", GetCrystalBallAlpha());
191                 printf("actual value of fCrystalBallNorm  ----------------------------->> | %f \n", GetCrystalBallNorm());
192         }else{
193                 printf("actual value of fMpv ------------------------------------------>> | %f \n", GetCrystalBallMmean());
194                 printf("actual value of fConstRovL ------------------------------------>> | %f \n", GetCrystalBallNexp());
195                 printf("actual value of fSigmaL --------------------------------------->> | %f \n", GetCrystalBallSigma());
196                 printf("actual value of fSigmaR --------------------------------------->> | %f \n", GetCrystalBallAlpha());
197         }
198
199         // back Mass func
200         printf("actual value of normBkg ----------------------------------------->> | %f \n", GetBkgInvMassNorm());
201         printf("actual value of meanBkg ----------------------------------------->> | %f \n", GetBkgInvMassMean());
202         printf("actual value of slopeBkg ---------------------------------------->> | %f \n", GetBkgInvMassSlope());
203         printf("actual value of constBkg ---------------------------------------->> | %f \n", GetBkgInvMassConst());
204         // resolution func (First-First)
205         printf("actual value of norm1Gauss (FF) --------------------------------->> | %f \n", GetNormGaus1ResFunc(2));
206         printf("actual value of norm2Gauss (FF) --------------------------------->> | %f \n", GetNormGaus2ResFunc(2));
207         printf("actual value of fMean1Res (FF) ---------------------------------->> | %f \n", GetResMean1(2));
208         printf("actual value of fSigma1Res (FF) --------------------------------->> | %f \n", GetResSigma1(2));
209         printf("actual value of fMean2Res (FF) ---------------------------------->> | %f \n", GetResMean2(2));
210         printf("actual value of fSigma2Res (FF) --------------------------------->> | %f \n", GetResSigma2(2));
211         
212         printf("actual value of alfaRes (FF) ------------------------------------>> | %f \n", GetResAlfa(2));
213         printf("actual value of lambdaRes (FF) ---------------------------------->> | %f \n", GetResLambda(2)); 
214         printf("actual value of normExpRes (FF) --------------------------------->> | %f \n", GetResNormExp(2)); 
215        
216         // resolution func (First-Second)
217         printf("actual value of norm1Gauss (FS) --------------------------------->> | %f \n", GetNormGaus1ResFunc(1));
218         printf("actual value of norm2Gauss (FS) --------------------------------->> | %f \n", GetNormGaus2ResFunc(1));
219         printf("actual value of fMean1Res (FS) ---------------------------------->> | %f \n", GetResMean1(1));
220         printf("actual value of fSigma1Res (FS) --------------------------------->> | %f \n", GetResSigma1(1));
221         printf("actual value of fMean2Res (FS) ---------------------------------->> | %f \n", GetResMean2(1));
222         printf("actual value of fSigma2Res (FS) --------------------------------->> | %f \n", GetResSigma2(1));
223         
224         printf("actual value of alfaRes (FS) ------------------------------------>> | %f \n", GetResAlfa(1));
225         printf("actual value of lambdaRes (FS) ---------------------------------->> | %f \n", GetResLambda(1));    
226         printf("actual value of normExpRes (FS) --------------------------------->> | %f \n", GetResNormExp(1));    
227         
228         // resolution func (Second-Second) 
229         printf("actual value of norm1Gauss (SS) --------------------------------->> | %f \n", GetNormGaus1ResFunc(0));
230         printf("actual value of norm2Gauss (SS) --------------------------------->> | %f \n", GetNormGaus2ResFunc(0));
231         printf("actual value of fMean1Res (SS) ---------------------------------->> | %f \n", GetResMean1(0));
232         printf("actual value of fSigma1Res (SS) --------------------------------->> | %f \n", GetResSigma1(0));
233         printf("actual value of fMean2Res (SS) ---------------------------------->> | %f \n", GetResMean2(0));
234         printf("actual value of fSigma2Res (SS) --------------------------------->> | %f \n", GetResSigma2(0));
235         
236         printf("actual value of alfaRes (SS) ------------------------------------>> | %f \n", GetResAlfa(0));
237         printf("actual value of lambdaRes (SS) ---------------------------------->> | %f \n", GetResLambda(0));    
238         printf("actual value of normExpRes (SS) --------------------------------->> | %f \n", GetResNormExp(0));    
239
240         printf("\n");
241         // integrals constants
242         printf("Actual value of normalization integral for MassSig ---------------->> | %f \n", GetIntegralMassSig());
243         printf("Actual value of normalization integral for MassBkg ---------------->> | %f \n", GetIntegralMassBkg());
244
245         printf("\n");
246 }
247 //_________________________________________________________________________________________________
248 void AliDielectronBtoJPSItoEleCDFfitFCN::SetResolutionConstants(const Double_t* resolutionConst, Int_t type) 
249 {
250         //
251         // Resolution function is parametrized as the sum of two gaussian
252         //
253         Int_t index = (2-type)*9; 
254         fParameters[20+index]=resolutionConst[1]; //mean1
255         fParameters[22+index]=resolutionConst[4]; //mean2
256         fParameters[18+index]=resolutionConst[0]; //norm1
257         fParameters[21+index]=resolutionConst[2]; //sigma1
258         fParameters[23+index]=resolutionConst[5]; //sigma2
259         fParameters[19+index]=resolutionConst[3]; //norm2
260
261         //exp values
262         fParameters[24+index]=resolutionConst[6]; //alfa
263         fParameters[25+index]=resolutionConst[7]; //lambda
264         fParameters[26+index]=resolutionConst[8]; //norm3
265         return;
266 }
267
268 //________________________________________________________________________________________________
269 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFfunc(Double_t x, Double_t m, Int_t type) const 
270 {
271         // evaluate likelihood function
272         return fParameters[8]*EvaluateCDFfuncSignalPart(x,m,type) + (1. - fParameters[8])*EvaluateCDFfuncBkgPart(x,m,type);
273 }
274
275 //_________________________________________________________________________________________________
276 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFfuncNorm(Double_t x, Double_t m, Int_t type) const
277 {
278         // evaluate likelihood function
279         return EvaluateCDFfunc(x,m,type);
280 }
281
282 //_________________________________________________________________________________________________
283 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFfuncSignalPart(Double_t x, Double_t m, Int_t type) const 
284 {
285   // evaluate psproper signal   
286   return EvaluateCDFDecayTimeSigDistr(x,type)*(EvaluateCDFInvMassSigDistr(m)/fintmMassSig); 
287 }
288
289 //_________________________________________________________________________________________________
290 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeSigDistr(Double_t x, Int_t type) const
291 {
292         //
293         // Implementation of the Background part of the Likelyhood function
294         // 
295         Double_t retvalue = 0.;
296         Double_t funBnorm = FunB(x,type);
297         Double_t funPnorm = ResolutionFunc(x,type);
298         retvalue = fParameters[7]*funBnorm + (1. - fParameters[7])*funPnorm;
299         return retvalue;
300 }
301
302 //_________________________________________________________________________________________________
303 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFInvMassSigDistr(Double_t m) const
304
305         //
306         // Parametrization of signal part invariant mass distribution
307         // It can be either Crystal Ball function or sum of two Landau
308         //
309         Double_t fitval = 0.;
310
311         if(fCrystalBallParam){
312                 Double_t t = (m-fParameters[9])/fParameters[11]; ;
313                 if (fParameters[12] < 0) t = -t;
314                 Double_t absAlpha = TMath::Abs((Double_t)fParameters[12]);
315
316                 if (t >= -absAlpha) {
317                         return fParameters[13]*TMath::Exp(-0.5*t*t);
318                 }
319                 else {
320                         Double_t a =  TMath::Power(fParameters[10]/absAlpha,fParameters[10])* TMath::Exp(-0.5*absAlpha*absAlpha);
321                         Double_t b= fParameters[10]/absAlpha - absAlpha;
322                         fitval = (fParameters[13]*a/TMath::Power(b - t, fParameters[10]));
323                         return fitval;
324                 }
325         }else{
326                 Double_t t=-1*m;
327                 Double_t tmpv=-1*fParameters[9];
328                 fitval=TMath::Sqrt(TMath::Landau(t,tmpv,fParameters[11]));
329                 fitval += fParameters[10]*(TMath::Landau(m,fParameters[9],fParameters[12]));
330                 return fitval;
331         }
332 }
333 //_________________________________________________________________________________________________
334 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::FunB(Double_t x, Int_t type) const  
335 {
336         //  
337         // Parameterisation of the fit function for the x part of the Background
338         //
339         Double_t np = 1000.0;
340         Double_t sc = 10.;
341         Double_t sigma3 = 1000.; // valore usato nella macro
342         Double_t xprime;
343         Double_t sum = 0.0;
344         Double_t xlow,xupp;
345         Double_t step;
346         Double_t i;
347         xlow = x - sc * sigma3 ;
348         xupp = x + sc * sigma3 ;
349         step = (xupp-xlow) / np;
350         Double_t csiMCxprime = 0.;
351         Double_t resolutionxdiff = 0.;
352         Double_t xdiff = 0.;
353
354         for(i=1.0; i<=np; i++){
355                 xprime = xlow + (i-.5) * step;
356                 csiMCxprime = CsiMC(xprime);
357                 xdiff = xprime - x;
358                 resolutionxdiff = ResolutionFunc(xdiff, type); // normalized value
359                 sum += csiMCxprime * resolutionxdiff;
360         }
361
362         return step * sum ;
363 }
364 //_________________________________________________________________________________________________
365 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::FunP(Double_t x, Int_t type) const 
366 {
367         //
368         //  Parameterisation of the Prompt part for the x distribution
369         //
370         return ResolutionFunc(x,type);
371 }
372
373
374 //_________________________________________________________________________________________________
375 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::CsiMC(Double_t x) const 
376 {
377         //
378         //  Distribution (template) of the x distribution for the x variable 
379         //  for the J/psi coming from Beauty hadrons
380         //
381         Double_t returnvalue = 0.;
382
383         if((fhCsiMC->FindBin(x) > 0) && (fhCsiMC->FindBin(x) < fhCsiMC->GetNbinsX()+1))  
384         returnvalue = fhCsiMC->Interpolate(x);
385
386         return returnvalue;
387 }
388
389 //_________________________________________________________________________________________________
390 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFfuncBkgPart(Double_t x,Double_t m, Int_t type) const 
391 {
392         //
393         // Return the part of the likelihood function for the background hypothesis
394         //
395         return EvaluateCDFDecayTimeBkgDistr(x,type)*(EvaluateCDFInvMassBkgDistr(m)/fintmMassBkg);
396 }  
397
398 //_________________________________________________________________________________________________
399 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeBkgDistr(Double_t x, Int_t type) const 
400 {
401         //
402         // it returns the value of the probability to have a given x for the background 
403         //
404
405         Double_t ret = fParameters[0]/(fParameters[0]+fParameters[1]+fParameters[2]+fParameters[3])*ResolutionFunc(x,type) + fParameters[1]/(fParameters[0]+fParameters[1]+fParameters[2]+fParameters[3])*FunBkgPos(x,type) +  fParameters[2]/(fParameters[0]+fParameters[1]+fParameters[2]+fParameters[3])*FunBkgNeg(x,type) + fParameters[3]/(fParameters[0]+fParameters[1]+fParameters[2]+fParameters[3])*FunBkgSym(x,type);
406         return ret;
407 }
408
409 //_________________________________________________________________________________________________
410 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFInvMassBkgDistr(Double_t m) const 
411 {
412         //
413         // it returns the value of the probability to have a given mass for the background
414         //
415         Double_t value = 0.;
416         value = fParameters[14]*TMath::Exp(-1*(m-fParameters[15])/fParameters[16]) + fParameters[17];  
417         return value;
418 }
419 //_________________________________________________________________________________________________
420 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::FunBkgPos(Double_t x,Int_t type) const 
421 {
422         //
423         // exponential with positive slopes for the background part (x)
424         //
425         Double_t np = 1000.0;
426         Double_t sc = 10.;      
427         Double_t sigma3 = 1000.; // valore usato nella macro 
428         Double_t xprime;
429         Double_t sum = 0.0;
430         Double_t xlow,xupp;
431         Double_t step;
432         Double_t i;
433         xlow = x - sc * sigma3 ;
434         xupp = x + sc * sigma3 ;
435         step = (xupp-xlow) / np;
436
437         for(i=1.0; i<=np/2; i++) {
438                 xprime = xlow + (i-.5) * step;
439                 if (xprime > 0) {sum += fParameters[4] * TMath::Exp(-1*xprime*fParameters[4])*(ResolutionFunc(xprime-x,type));}
440                 xprime = xupp - (i-.5) * step;
441                 if (xprime > 0) {sum += fParameters[4] * TMath::Exp(-1*xprime*fParameters[4])*(ResolutionFunc(xprime-x,type));}
442         }
443
444         return step * sum ;
445 }
446 //_________________________________________________________________________________________________
447 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::FunBkgNeg(Double_t x, Int_t type) const 
448 {
449         //
450         // exponential with negative slopes for the background part (x)
451         //
452         Double_t np = 1000.0;
453         Double_t sc = 10.;      
454         Double_t sigma3 = 1000.; 
455         Double_t xprime;
456         Double_t sum = 0.0;
457         Double_t xlow,xupp;
458         Double_t step;
459         Double_t i;
460         xlow = x - sc * sigma3 ;
461         xupp = x + sc * sigma3 ;
462         step = (xupp-xlow) / np;
463
464         for(i=1.0; i<=np/2; i++) {
465
466                 xprime = xlow + (i-.5) * step;
467                 if (xprime < 0) {sum += fParameters[5] * TMath::Exp(xprime*fParameters[5]) * (ResolutionFunc(xprime-x,type));}
468
469                 xprime = xupp - (i-.5) * step;
470                 if (xprime < 0) {sum += fParameters[5] * TMath::Exp(xprime*fParameters[5]) * (ResolutionFunc(xprime-x,type));}
471         }
472
473         return step * sum ;
474 }
475 //_________________________________________________________________________________________________
476 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::FunBkgSym(Double_t x, Int_t type) const 
477 {
478         //
479         // exponential with both positive and negative slopes for the background part (x)
480         //
481         Double_t np = 1000.0;
482         Double_t sc = 10.;      
483         Double_t sigma3 = 1000.; 
484         Double_t xprime;
485         Double_t sum1 = 0.0;
486         Double_t sum2 = 0.0;
487         Double_t xlow,xupp;
488         Double_t step;
489         Double_t i;
490         xlow = x - sc * sigma3 ;
491         xupp = x + sc * sigma3 ;
492         step = (xupp-xlow) / np;
493
494         for(i=1.0; i<=np/2; i++) {
495
496                 xprime = xlow + (i-.5) * step;
497                 if (xprime > 0) {sum1 += 0.5 * fParameters[6]*TMath::Exp(-1*xprime*fParameters[6]) * (ResolutionFunc(xprime-x,type));}
498                 if (xprime < 0) {sum2 += 0.5 * fParameters[6]*TMath::Exp(xprime*fParameters[6]) * (ResolutionFunc(xprime-x,type));}
499
500                 xprime = xupp - (i-.5) * step;
501                 if (xprime > 0) {sum1 += 0.5 * fParameters[6]*TMath::Exp(-1*xprime*fParameters[6]) * (ResolutionFunc(xprime-x,type));} 
502                 if (xprime < 0) {sum2 += 0.5 * fParameters[6]*TMath::Exp(xprime*fParameters[6]) * (ResolutionFunc(xprime-x,type));}
503         }
504
505         return step*(sum1 + sum2) ;
506 }
507 //_________________________________________________________________________________________________
508 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::ResolutionFunc(Double_t x, Int_t type) const  
509 {
510         //
511         // parametrization with 2 gaus
512         //
513         Int_t index = (2-type)*9;
514         Double_t mean1 = fParameters[20+index];
515         Double_t mean2 = fParameters[22+index];
516         Double_t norm1 = fParameters[18+index];
517         Double_t sigma1 = fParameters[21+index];
518         Double_t sigma2 = fParameters[23+index];
519         Double_t norm2 = fParameters[19+index];
520         //exp values
521         Double_t alfa = fParameters[24+index];
522         Double_t lambda = fParameters[25+index];
523         Double_t norm3 = fParameters[26+index];
524  
525         Double_t ret = 0.; Double_t fitval = 0; 
526         if(TMath::Abs(x)<=alfa) fitval = (lambda-1)/(2*alfa*lambda);
527         else  fitval = ((lambda-1)/(2*alfa*lambda))*TMath::Power(alfa,lambda)*(TMath::Power(TMath::Abs(x),-1*lambda));
528
529         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;
530
531         return ret;
532 }     
533
534 //_________________________________________________________________________________________________
535 TF1* AliDielectronBtoJPSItoEleCDFfitFCN::GetCsiMC(Double_t xmin, Double_t xmax, Double_t normalization) 
536 {
537         // return the pointer to the templateMC function 
538         TF1* templateMC = new TF1("MCtemplate",this,&AliDielectronBtoJPSItoEleCDFfitFCN::CsiMCfunc,xmin,xmax,1);
539         templateMC->SetParameter(0,normalization);
540         templateMC->SetNpx(5000);
541         return (TF1*)templateMC->Clone();
542 }
543
544 //__________________________________________________________________________________________________
545 TF1* AliDielectronBtoJPSItoEleCDFfitFCN::GetResolutionFunc(Double_t xmin, Double_t xmax, Double_t normalization, Double_t type){
546         // return the pointer to the resolution function
547         TF1* resFunc = new TF1(Form("resolutionFunc_%f",type),this,&AliDielectronBtoJPSItoEleCDFfitFCN::ResolutionFuncf,xmin,xmax,2);
548         resFunc->SetParameter(0,normalization);
549         resFunc->SetParameter(1,type);
550         resFunc->SetNpx(5000);
551         return (TF1*)resFunc->Clone();
552 }
553
554 //__________________________________________________________________________________________________
555 TF1* AliDielectronBtoJPSItoEleCDFfitFCN::GetResolutionFuncAllTypes(Double_t xmin, Double_t xmax, Double_t normalization){
556         // return the pointer to the resolution function
557         TF1* resFunc = new TF1(Form("resolutionFunc"),this,&AliDielectronBtoJPSItoEleCDFfitFCN::ResolutionFuncAllTypes,xmin,xmax,1);
558         resFunc->SetParameter(0,normalization);
559         resFunc->SetNpx(5000);
560         return (TF1*)resFunc->Clone();
561 }
562
563 //___________________________________________________________________________________________________
564 TF1* AliDielectronBtoJPSItoEleCDFfitFCN::GetEvaluateCDFDecayTimeBkgDistr(Double_t xmin, Double_t xmax, Double_t normalization, Double_t type){
565         // return the pointer to the background x distribution function
566         TF1 *backFunc = new TF1(Form("backFunc_%f",type),this,&AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeBkgDistrFunc,xmin,xmax,2);
567         backFunc->SetParameter(0,normalization);
568         backFunc->SetParameter(1,type);
569         backFunc->SetNpx(5000);
570         return (TF1*)backFunc->Clone();
571 }
572
573 //___________________________________________________________________________________________________
574 TF1* AliDielectronBtoJPSItoEleCDFfitFCN::GetEvaluateCDFDecayTimeBkgDistrAllTypes(Double_t xmin, Double_t xmax, Double_t normalization){
575         // return the pointer to the background x distribution function
576         TF1 *backFunc = new TF1(Form("backFunc"),this,&AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeBkgDistrFuncAllTypes,xmin,xmax,1);
577         backFunc->SetParameter(0,normalization);
578         backFunc->SetNpx(5000);
579         return (TF1*)backFunc->Clone();
580 }
581
582 //__________________________________________________________________________________________________
583 TF1* AliDielectronBtoJPSItoEleCDFfitFCN::GetEvaluateCDFDecayTimeSigDistr(Double_t xmin, Double_t xmax, Double_t normalization, Double_t type){
584         // return the pointer to the signal x distribution function
585         TF1 *signFunc = new TF1("signalFunc",this,&AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeSigDistrFunc,xmin,xmax,2);
586         signFunc->SetParameter(0,normalization);
587         signFunc->SetParameter(1,type); 
588         signFunc->SetNpx(5000);
589         return (TF1*)signFunc->Clone();
590 }
591
592 //____________________________________________________________________________________________________
593 TF1* AliDielectronBtoJPSItoEleCDFfitFCN::GetEvaluateCDFInvMassBkgDistr(Double_t mMin, Double_t mMax, Double_t normalization){
594   // return the pointer to the invariant mass distribution for the background 
595   TF1 * invMassBkg = new TF1("invMassBkg",this,&AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFInvMassBkgDistrFunc,mMin,mMax,1);
596   invMassBkg->SetParameter(0,normalization);
597   invMassBkg->SetNpx(5000);
598   return (TF1*)invMassBkg->Clone();
599 }
600
601
602 //____________________________________________________________________________________________________
603 TF1* AliDielectronBtoJPSItoEleCDFfitFCN::GetEvaluateCDFInvMassSigDistr(Double_t mMin, Double_t mMax, Double_t normalization){
604   // return the pointer to the invariant mass distribution for the signal
605   TF1 * invMassSig = new TF1("invMassSig",this,&AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFInvMassSigDistrFunc,mMin,mMax,1);
606   invMassSig->SetParameter(0,normalization);
607   invMassSig->SetNpx(5000);
608   return (TF1*)invMassSig->Clone();
609 }
610
611 //____________________________________________________________________________________________________
612 TF1 *AliDielectronBtoJPSItoEleCDFfitFCN::GetEvaluateCDFInvMassTotalDistr(Double_t mMin, Double_t mMax, Double_t normalization){
613   // return the pointer to the invariant mass distribution
614   TF1 * invMassTot = new TF1("invMassTot",this,&AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFInvMassTotalDistr,mMin,mMax,1);
615   invMassTot->SetParameter(0,normalization);
616   invMassTot->SetNpx(5000);
617   return (TF1*)invMassTot->Clone();
618 }
619
620 //____________________________________________________________________________________________________
621 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFInvMassTotalDistr(const Double_t* x, const Double_t *par) const
622 {
623   // evaluate invariant mass total distribution
624   Double_t value = 0;
625   Double_t xx = x[0];
626   value = ((EvaluateCDFInvMassSigDistr(xx)/fintmMassSig)*fParameters[8] + (1-fParameters[8])*(EvaluateCDFInvMassBkgDistr(xx)/fintmMassBkg))*par[0];
627   return value;  
628 }
629
630 //____________________________________________________________________________________________________
631 TF1 *AliDielectronBtoJPSItoEleCDFfitFCN::GetEvaluateCDFDecayTimeTotalDistr(Double_t xMin, Double_t xMax, Double_t normalization,Double_t type){
632  // return the pointer to the pseudoproper distribution for the background
633  TF1 *decayTimeTot = new TF1(Form("decayTimeTot_%f",type),this,&AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeTotalDistr,xMin,xMax,2);
634  decayTimeTot->SetParameter(0,normalization);
635  decayTimeTot->SetParameter(1,type);
636  decayTimeTot->SetNpx(5000);
637  return (TF1*)decayTimeTot->Clone();
638 }
639
640 //____________________________________________________________________________________________________
641 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeTotalDistr(const Double_t* x, const Double_t *par) const
642 {
643  // evaluate the total pseudoproper distribution for a given candidate's type. par[1] should be the candidate's type.
644  Double_t value = 0;
645  Double_t xx = x[0];
646  value = (fParameters[8]*EvaluateCDFDecayTimeSigDistr(xx,(Int_t)par[1]) + (1-fParameters[8])*EvaluateCDFDecayTimeBkgDistr(xx,(Int_t)par[1]))*par[0];
647  return value;
648 }
649
650 //____________________________________________________________________________________________________
651 Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeTotalDistrAllTypes(const Double_t* x, const Double_t *par) const
652 {
653  // evaluate the total pseudoproper distribution considering all candidate's types
654  Double_t value = 0;
655  Double_t xx = x[0];
656
657  value = (fParameters[8]*(fWeightType[2]*EvaluateCDFDecayTimeSigDistr(xx,2)+fWeightType[1]*EvaluateCDFDecayTimeSigDistr(xx,1)+fWeightType[0]*EvaluateCDFDecayTimeSigDistr(xx,0)))+((1-fParameters[8])*(fWeightType[2]*EvaluateCDFDecayTimeBkgDistr(xx,2) + fWeightType[1]*EvaluateCDFDecayTimeBkgDistr(xx,1)+fWeightType[0]*EvaluateCDFDecayTimeBkgDistr(xx,0))); 
658
659  return value*par[0];
660 }
661
662 //____________________________________________________________________________________________________
663 TF1 *AliDielectronBtoJPSItoEleCDFfitFCN::GetEvaluateCDFDecayTimeTotalDistrAllTypes(Double_t xMin, Double_t xMax, Double_t normalization){
664  // return the pointer to the pseudoproper function for the background considering all candidate's types
665  TF1 *decayTimeTot = new TF1(Form("decayTimeTot"),this,&AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeTotalDistrAllTypes,xMin,xMax,1);
666  decayTimeTot->SetParameter(0,normalization);
667  decayTimeTot->SetNpx(5000);
668  return (TF1*)decayTimeTot->Clone();
669 }
670
671
672 //____________________________________________________________________________________________________
673 TF1 * AliDielectronBtoJPSItoEleCDFfitFCN::GetFunB(Double_t xmin, Double_t xmax, Double_t normalization, Double_t type){
674  // return the pointer to the function that describe secondary jpsi pseudoproper distribution for a given candidate's type
675  TF1* funb = new TF1(Form("secondaryJpsiConvolution_%f",type),this,&AliDielectronBtoJPSItoEleCDFfitFCN::FunBfunc,xmin,xmax,2);
676  funb->SetParameter(0,normalization);
677  funb->SetParameter(1,type);
678  funb->SetNpx(5000);
679  return (TF1*)funb->Clone();
680  }
681
682 //____________________________________________________________________________________________________
683 TF1 * AliDielectronBtoJPSItoEleCDFfitFCN::GetFunBAllTypes(Double_t xmin, Double_t xmax, Double_t normalization){
684 // return the pointer to the function that describe secondary jpsi pseudoproper distribution for all types
685  TF1* funb = new TF1(Form("secondaryJpsiConvolution"),this,&AliDielectronBtoJPSItoEleCDFfitFCN::FunBfuncAllTypes,xmin,xmax,1);
686  funb->SetParameter(0,normalization);
687  funb->SetNpx(5000);
688  return (TF1*)funb->Clone();
689  }
690
691