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