]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGDQ/dielectron/AliDielectronBtoJPSItoEleCDFfitFCN.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
add task pp filter
[u/mrichter/AliRoot.git] / PWGDQ / dielectron / AliDielectronBtoJPSItoEleCDFfitFCN.cxx
CommitLineData
ba15fdfb 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 **************************************************************************/
bc7a3220 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
ba15fdfb 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
34ClassImp(AliDielectronBtoJPSItoEleCDFfitFCN)
35
36 //_________________________________________________________________________________________________
37 AliDielectronBtoJPSItoEleCDFfitFCN::AliDielectronBtoJPSItoEleCDFfitFCN() :
38 fFPlus(0.),
39 fFMinus(0.),
40 fFSym(0.),
41 fintmMassSig(1.),
42 fintmMassBkg(1.),
43 fhCsiMC(0x0),
bc7a3220 44 fShiftTemplate(0.),
ba15fdfb 45 fMassWndHigh(0.),
46 fMassWndLow(0.),
bc7a3220 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)
ba15fdfb 61{
62 //
63 // constructor
64 //
ba15fdfb 65 SetMassWndHigh(0.2);
66 SetMassWndLow(0.5);
5720c765 67 fWeightType[0] = 1.; fWeightType[1] = 1.; fWeightType[2] = 1.;
bc7a3220 68 for(Int_t iPar = 0; iPar < 49; iPar++) fParameters[iPar] = 0.;
5720c765 69 fParameters[9] = 1.;fParameters[11] = 1.;fParameters[12] = 1.;
bc7a3220 70
71
72 AliInfo("Instance of AliDielectronBtoJPSItoEleCDFfitFCN-class created");
ba15fdfb 73}
74//_________________________________________________________________________________________________
75AliDielectronBtoJPSItoEleCDFfitFCN::AliDielectronBtoJPSItoEleCDFfitFCN(const AliDielectronBtoJPSItoEleCDFfitFCN& source) :
76 TNamed(source),
77 fFPlus(source.fFPlus),
bc7a3220 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)
ba15fdfb 100{
101 //
102 // Copy constructor
103 //
bc7a3220 104 for(Int_t iPar = 0; iPar < 49; iPar++) fParameters[iPar] = source.fParameters[iPar];
5720c765 105 for(Int_t iW=0; iW<2; iW++) fWeightType[iW] = source.fWeightType[iW];
ba15fdfb 106}
107//_________________________________________________________________________________________________
108AliDielectronBtoJPSItoEleCDFfitFCN& AliDielectronBtoJPSItoEleCDFfitFCN::operator=(const AliDielectronBtoJPSItoEleCDFfitFCN& source)
109{
110 //
111 // Assignment operator
112 //
113 if(&source == this) return *this;
bc7a3220 114 fFPlus = source.fFPlus;
ba15fdfb 115 fFMinus = source.fFMinus;
116 fFSym = source.fFSym;
117 fintmMassSig = source.fintmMassSig;
118 fintmMassBkg = source.fintmMassBkg;
119 fhCsiMC = source.fhCsiMC;
bc7a3220 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;
ba15fdfb 129 fCrystalBallParam = source.fCrystalBallParam;
bc7a3220 130 fExponentialParam = source.fExponentialParam;
ba15fdfb 131
bc7a3220 132 for(Int_t iPar = 0; iPar < 49; iPar++) fParameters[iPar] = source.fParameters[iPar];
ba15fdfb 133 return *this;
134}
135//_________________________________________________________________________________________________
136AliDielectronBtoJPSItoEleCDFfitFCN::~AliDielectronBtoJPSItoEleCDFfitFCN()
137{
138 //
139 // Default destructor
140 //
141
142 delete fhCsiMC;
bc7a3220 143 delete fFunBSaved;
144 delete fFunBkgSaved;
145 for(Int_t iPar = 0; iPar < 49; iPar++) fParameters[iPar] = 0.;
ba15fdfb 146}
147//_________________________________________________________________________________________________
148Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateLikelihood(const Double_t* pseudoproperdecaytime,
34e8fe1d 149 const Double_t* invariantmass, const Double_t *pt, const Int_t *type, Int_t ncand) const
ba15fdfb 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++) {
bc7a3220 159 f = EvaluateCDFfuncNorm(pseudoproperdecaytime[i],invariantmass[i],pt[i],type[i]);
160 if(f <= 0.) continue;
161 ret += -2.*TMath::Log(f);
ba15fdfb 162 }
5720c765 163 return ret;
ba15fdfb 164}
165//_________________________________________________________________________________________________
166void AliDielectronBtoJPSItoEleCDFfitFCN::SetAllParameters(const Double_t* parameters)
167{
168 //
169 // Sets array of FCN parameters
170 //
bc7a3220 171 for(Int_t index = 0; index < 49; index++) fParameters[index] = parameters[index];
ba15fdfb 172}
173//_________________________________________________________________________________________________
174void 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;
5720c765 208 if(intmMassBkg < 1.e-05) intmMassBkg = 1.;
209 SetIntegralMassBkg(intmMassBkg);
ba15fdfb 210}
211//_________________________________________________________________________________________________
212void 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());
bc7a3220 223 printf("actual value of fSym1 ----------------------------------------->> | %f \n", GetFSym1());
ba15fdfb 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());
bc7a3220 227 printf("actual value of fOneOvLamSym1 --------------------------------->> | %f \n", GetLamSym1());
ba15fdfb 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
bc7a3220 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
5720c765 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");
ba15fdfb 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//_________________________________________________________________________________________________
d07ad639 303void AliDielectronBtoJPSItoEleCDFfitFCN::SetResolutionConstants(const Double_t* resolutionConst, Int_t type)
ba15fdfb 304{
305 //
306 // Resolution function is parametrized as the sum of two gaussian
307 //
5720c765 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;
ba15fdfb 321}
322
5720c765 323//________________________________________________________________________________________________
bc7a3220 324Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFfunc(Double_t x, Double_t m, Double_t pt, Int_t type) const
ba15fdfb 325{
5720c765 326 // evaluate likelihood function
bc7a3220 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);
ba15fdfb 329}
330
331//_________________________________________________________________________________________________
bc7a3220 332Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFfuncNorm(Double_t x, Double_t m, Double_t pt, Int_t type) const
ba15fdfb 333{
5720c765 334 // evaluate likelihood function
bc7a3220 335 return EvaluateCDFfunc(x,m,pt, type);
ba15fdfb 336}
337
338//_________________________________________________________________________________________________
bc7a3220 339Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFfuncSignalPart(Double_t x, Double_t m, Double_t pt, Int_t type) const
ba15fdfb 340{
5720c765 341 // evaluate psproper signal
bc7a3220 342 return EvaluateCDFDecayTimeSigDistr(x,pt, type)*(EvaluateCDFInvMassSigDistr(m)/fintmMassSig);
ba15fdfb 343}
344
345//_________________________________________________________________________________________________
bc7a3220 346Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeSigDistr(Double_t x, Double_t pt, Int_t type) const
ba15fdfb 347{
348 //
349 // Implementation of the Background part of the Likelyhood function
350 //
ba15fdfb 351 Double_t retvalue = 0.;
bc7a3220 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;
ba15fdfb 355 return retvalue;
356}
357
358//_________________________________________________________________________________________________
359Double_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 //
ba15fdfb 365 Double_t fitval = 0.;
bc7a3220 366 // change inv Mass RMS fChangeMass
ba15fdfb 367 if(fCrystalBallParam){
bc7a3220 368 Double_t t = ((m-fParameters[9])/fChangeMass)/fParameters[11]; ;
ba15fdfb 369 if (fParameters[12] < 0) t = -t;
ba15fdfb 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//_________________________________________________________________________________________________
bc7a3220 390Double_t AliDielectronBtoJPSItoEleCDFfitFCN::FunB(Double_t x, Double_t pt, Int_t type) const
ba15fdfb 391{
392 //
393 // Parameterisation of the fit function for the x part of the Background
394 //
bc7a3220 395 Double_t np = 1000.0;
ba15fdfb 396 Double_t sc = 10.;
bc7a3220 397 Double_t sigma3 = 1000.;
ba15fdfb 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;
bc7a3220 414 resolutionxdiff = ResolutionFunc(xdiff, pt, type); // normalized value
ba15fdfb 415 sum += csiMCxprime * resolutionxdiff;
416 }
bc7a3220 417
ba15fdfb 418 return step * sum ;
419}
bc7a3220 420
ba15fdfb 421//_________________________________________________________________________________________________
bc7a3220 422Double_t AliDielectronBtoJPSItoEleCDFfitFCN::FunP(Double_t x, Double_t pt, Int_t type) const
ba15fdfb 423{
424 //
425 // Parameterisation of the Prompt part for the x distribution
426 //
bc7a3220 427 return ResolutionFunc(x, pt, type);
ba15fdfb 428}
429
430
431//_________________________________________________________________________________________________
432Double_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.;
bc7a3220 439
440 if((fhCsiMC->FindBin(x-fShiftTemplate) > 0) && (fhCsiMC->FindBin(x-fShiftTemplate) < fhCsiMC->GetNbinsX()+1))
441 returnvalue = fhCsiMC->Interpolate(x-fShiftTemplate);
ba15fdfb 442
ba15fdfb 443
444 return returnvalue;
445}
446
447//_________________________________________________________________________________________________
bc7a3220 448Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFfuncBkgPart(Double_t x,Double_t m, Double_t pt, Int_t type) const
ba15fdfb 449{
450 //
451 // Return the part of the likelihood function for the background hypothesis
452 //
bc7a3220 453 Double_t bkgValx = fMultivariate ? EvaluateCDFDecayTimeBkgDistrSaved(x,type,m,pt) : EvaluateCDFDecayTimeBkgDistr(x,type,m,pt);
454 return bkgValx*(EvaluateCDFInvMassBkgDistr(m)/fintmMassBkg);
455}
456
ba15fdfb 457
458//_________________________________________________________________________________________________
bc7a3220 459Double_t AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeBkgDistr(Double_t x, Int_t type, Double_t m, Double_t pt) const
ba15fdfb 460{
bc7a3220 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);
ba15fdfb 469 return ret;
470}
471
472//_________________________________________________________________________________________________
473Double_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.;
bc7a3220 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;
ba15fdfb 483}
484//_________________________________________________________________________________________________
bc7a3220 485Double_t AliDielectronBtoJPSItoEleCDFfitFCN::FunBkgPos(Double_t x, Double_t pt, Int_t type) const
ba15fdfb 486{
bc7a3220 487 //
ba15fdfb 488 // exponential with positive slopes for the background part (x)
489 //
bc7a3220 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
ba15fdfb 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;
bc7a3220 506 if (xprime > 0) {sum += fParameters[4] * TMath::Exp(-1*xprime*fParameters[4])*(ResolutionFunc(xprime-x,pt,type));}
ba15fdfb 507 xprime = xupp - (i-.5) * step;
bc7a3220 508 if (xprime > 0) {sum += fParameters[4] * TMath::Exp(-1*xprime*fParameters[4])*(ResolutionFunc(xprime-x,pt,type));}
509 }
ba15fdfb 510
511 return step * sum ;
512}
513//_________________________________________________________________________________________________
bc7a3220 514Double_t AliDielectronBtoJPSItoEleCDFfitFCN::FunBkgNeg(Double_t x, Double_t pt, Int_t type) const
515{
ba15fdfb 516 //
517 // exponential with negative slopes for the background part (x)
518 //
bc7a3220 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
ba15fdfb 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;
bc7a3220 536 if (xprime < 0) {sum += fParameters[5] * TMath::Exp(xprime*fParameters[5]) * (ResolutionFunc(xprime-x,pt,type));}
ba15fdfb 537
538 xprime = xupp - (i-.5) * step;
bc7a3220 539 if (xprime < 0) {sum += fParameters[5] * TMath::Exp(xprime*fParameters[5]) * (ResolutionFunc(xprime-x,pt,type));}
ba15fdfb 540 }
541
542 return step * sum ;
543}
544//_________________________________________________________________________________________________
bc7a3220 545Double_t AliDielectronBtoJPSItoEleCDFfitFCN::FunBkgSym(Double_t x, Double_t pt, Int_t type) const
ba15fdfb 546{
547 //
548 // exponential with both positive and negative slopes for the background part (x)
549 //
bc7a3220 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
ba15fdfb 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;
bc7a3220 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));}
ba15fdfb 570
571 xprime = xupp - (i-.5) * step;
bc7a3220 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));}
ba15fdfb 574 }
575
576 return step*(sum1 + sum2) ;
577}
bc7a3220 578
579//_________________________________________________________________________________________________
580Double_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
ba15fdfb 615//_________________________________________________________________________________________________
bc7a3220 616Double_t AliDielectronBtoJPSItoEleCDFfitFCN::ResolutionFunc(Double_t x, Double_t pt, Int_t type) const
ba15fdfb 617{
618 //
619 // parametrization with 2 gaus
620 //
bc7a3220 621 x = x/fChangeResolution;
5720c765 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
bc7a3220 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
5720c765 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
bc7a3220 656 return ret/fChangeResolution;
5720c765 657}
ba15fdfb 658
bc7a3220 659
ba15fdfb 660//_________________________________________________________________________________________________
5720c765 661TF1* AliDielectronBtoJPSItoEleCDFfitFCN::GetCsiMC(Double_t xmin, Double_t xmax, Double_t normalization)
ba15fdfb 662{
663 // return the pointer to the templateMC function
5720c765 664 TF1* templateMC = new TF1("MCtemplate",this,&AliDielectronBtoJPSItoEleCDFfitFCN::CsiMCfunc,xmin,xmax,1);
665 templateMC->SetParameter(0,normalization);
666 templateMC->SetNpx(5000);
ba15fdfb 667 return (TF1*)templateMC->Clone();
668}
669
670//__________________________________________________________________________________________________
bc7a3220 671TF1* AliDielectronBtoJPSItoEleCDFfitFCN::GetResolutionFunc(Double_t xmin, Double_t xmax, Double_t normalization, Double_t pt, Int_t type){
ba15fdfb 672 // return the pointer to the resolution function
bc7a3220 673 TF1* resFunc = new TF1(Form("resolutionFunc_%d",type),this,&AliDielectronBtoJPSItoEleCDFfitFCN::ResolutionFuncf,xmin,xmax,3);
5720c765 674 resFunc->SetParameter(0,normalization);
bc7a3220 675 resFunc->SetParameter(1,pt);
676 resFunc->SetParameter(2,(Double_t)type);
ba15fdfb 677 resFunc->SetNpx(5000);
678 return (TF1*)resFunc->Clone();
679}
680
5720c765 681//__________________________________________________________________________________________________
682TF1* 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
ba15fdfb 690//___________________________________________________________________________________________________
bc7a3220 691TF1* AliDielectronBtoJPSItoEleCDFfitFCN::GetEvaluateCDFDecayTimeBkgDistr(Double_t xmin, Double_t xmax, Double_t normalization, Int_t type, Double_t mass, Double_t pt, Int_t npx){
ba15fdfb 692 // return the pointer to the background x distribution function
bc7a3220 693 TF1 *backFunc = new TF1(Form("backFunc_%d_%1.2f_%1.2f",type,mass,pt),this,&AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeBkgDistrFunc,xmin,xmax,4);
5720c765 694 backFunc->SetParameter(0,normalization);
bc7a3220 695 backFunc->SetParameter(1,(Double_t)type);
696 backFunc->SetParameter(2,mass);
697 backFunc->SetParameter(3,pt);
698 backFunc->SetNpx(npx);
ba15fdfb 699 return (TF1*)backFunc->Clone();
700}
701
5720c765 702//___________________________________________________________________________________________________
703TF1* AliDielectronBtoJPSItoEleCDFfitFCN::GetEvaluateCDFDecayTimeBkgDistrAllTypes(Double_t xmin, Double_t xmax, Double_t normalization){
704 // return the pointer to the background x distribution function
bc7a3220 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();
5720c765 709}
710
ba15fdfb 711//__________________________________________________________________________________________________
5720c765 712TF1* AliDielectronBtoJPSItoEleCDFfitFCN::GetEvaluateCDFDecayTimeSigDistr(Double_t xmin, Double_t xmax, Double_t normalization, Double_t type){
ba15fdfb 713 // return the pointer to the signal x distribution function
5720c765 714 TF1 *signFunc = new TF1("signalFunc",this,&AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeSigDistrFunc,xmin,xmax,2);
715 signFunc->SetParameter(0,normalization);
716 signFunc->SetParameter(1,type);
ba15fdfb 717 signFunc->SetNpx(5000);
718 return (TF1*)signFunc->Clone();
719}
5720c765 720
721//____________________________________________________________________________________________________
722TF1* 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//____________________________________________________________________________________________________
732TF1* 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//____________________________________________________________________________________________________
741TF1 *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//____________________________________________________________________________________________________
750Double_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//____________________________________________________________________________________________________
bc7a3220 760TF1 *AliDielectronBtoJPSItoEleCDFfitFCN::GetEvaluateCDFDecayTimeTotalDistr(Double_t xMin, Double_t xMax, Double_t normalization,Double_t pt, Int_t type){
5720c765 761 // return the pointer to the pseudoproper distribution for the background
bc7a3220 762 TF1 *decayTimeTot = new TF1(Form("decayTimeTot_%d",type),this,&AliDielectronBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeTotalDistr,xMin,xMax,3);
5720c765 763 decayTimeTot->SetParameter(0,normalization);
bc7a3220 764 decayTimeTot->SetParameter(1,pt);
765 decayTimeTot->SetParameter(2,(Double_t)type);
5720c765 766 decayTimeTot->SetNpx(5000);
767 return (TF1*)decayTimeTot->Clone();
768}
769
770//____________________________________________________________________________________________________
771Double_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];
bc7a3220 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];
5720c765 777 return value;
778}
779
780//____________________________________________________________________________________________________
781Double_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
bc7a3220 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.)));
5720c765 788
789 return value*par[0];
790}
791
792//____________________________________________________________________________________________________
793TF1 *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//____________________________________________________________________________________________________
bc7a3220 803TF1 * AliDielectronBtoJPSItoEleCDFfitFCN::GetFunB(Double_t xmin, Double_t xmax, Double_t normalization, Double_t pt, Int_t type, Int_t npx){
5720c765 804 // return the pointer to the function that describe secondary jpsi pseudoproper distribution for a given candidate's type
bc7a3220 805 TF1* funb = new TF1(Form("secondaryJpsiConvolution_%d_%1.2f",type,pt),this,&AliDielectronBtoJPSItoEleCDFfitFCN::FunBfunc,xmin,xmax,3);
5720c765 806 funb->SetParameter(0,normalization);
bc7a3220 807 funb->SetParameter(1,pt);
808 funb->SetParameter(2,(Double_t)type);
809 funb->SetNpx(npx);
5720c765 810 return (TF1*)funb->Clone();
811 }
812
813//____________________________________________________________________________________________________
814TF1 * AliDielectronBtoJPSItoEleCDFfitFCN::GetFunBAllTypes(Double_t xmin, Double_t xmax, Double_t normalization){
bc7a3220 815 // return the pointer to the function that describe secondary jpsi pseudoproper distribution for all types
5720c765 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
bc7a3220 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//____________________________________________________________________________________________________
832void 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//_______________________________________________________________________________________________________
844void 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//_________________________________________________________________________________________________
870Double_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//________________________________________________________________________________________________
883void 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//_________________________________________________________________________________________________
937Double_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//_________________________________________________________________________________________________
960Double_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}
5720c765 974