]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/AliESDv0.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merging of AliESDV0MI and AliESDv0 classes (Yu.Belikov)
[u/mrichter/AliRoot.git] / STEER / AliESDv0.cxx
CommitLineData
e23730c7 1/**************************************************************************
2 * Copyright(c) 1998-1999, 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
4f679a16 16/* $Id$ */
17
e23730c7 18//-------------------------------------------------------------------------
19// Implementation of the ESD V0 vertex class
4f679a16 20// This class is part of the Event Data Summary
21// set of classes and contains information about
22// V0 kind vertexes generated by a neutral particle
e23730c7 23// Origin: Iouri Belikov, IReS, Strasbourg, Jouri.Belikov@cern.ch
d6a49f20 24// Modified by: Marian Ivanov, CERN, Marian.Ivanov@cern.ch
25// and Boris Hippolyte,IPHC, hippolyt@in2p3.fr
e23730c7 26//-------------------------------------------------------------------------
4f679a16 27
e23730c7 28#include <Riostream.h>
29#include <TMath.h>
90e48c0c 30#include <TDatabasePDG.h>
e23730c7 31#include <TPDGCode.h>
90e48c0c 32#include <TParticlePDG.h>
e23730c7 33
5f7789fc 34#include "AliLog.h"
e23730c7 35#include "AliESDv0.h"
c7bafca9 36#include "AliExternalTrackParam.h"
e23730c7 37
38ClassImp(AliESDv0)
39
d6a49f20 40AliESDV0Params AliESDv0::fgkParams;
41
90e48c0c 42AliESDv0::AliESDv0() :
43 TObject(),
d6a49f20 44 fOnFlyStatus(kFALSE),
90e48c0c 45 fPdgCode(kK0Short),
46 fEffMass(TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass()),
c028b974 47 fDcaV0Daughters(0),
48 fChi2V0(1.e+33),
90e48c0c 49 fNidx(0),
d6a49f20 50 fPidx(0),
51 fParamP(),
52 fParamN(),
53 fID(0),
54 fDist1(-1),
55 fDist2(-1),
56 fRr(-1),
57 fStatus(0),
58 fRow0(-1),
59 fDistNorm(0),
60 fDistSigma(0),
61 fChi2Before(0),
62 fNBefore(0),
63 fChi2After(0),
64 fNAfter(0),
65 fPointAngleFi(0),
66 fPointAngleTh(0),
67 fPointAngle(0)
90e48c0c 68{
e23730c7 69 //--------------------------------------------------------------------
70 // Default constructor (K0s)
71 //--------------------------------------------------------------------
6605de26 72
73 for (Int_t i=0; i<3; i++) {
74 fPos[i] = 0.;
75 fNmom[i] = 0.;
76 fPmom[i] = 0.;
77 }
78
79 for (Int_t i=0; i<6; i++) {
80 fPosCov[i]= 0.;
81 fNmomCov[i] = 0.;
82 fPmomCov[i] = 0.;
83 }
d6a49f20 84
85 for (Int_t i=0;i<5;i++){
86 fRP[i]=fRM[i]=0;
87 }
88 fLab[0]=fLab[1]=-1;
89 fIndex[0]=fIndex[1]=-1;
90 for (Int_t i=0;i<6;i++){fClusters[0][i]=0; fClusters[1][i]=0;}
91 fNormDCAPrim[0]=fNormDCAPrim[1]=0;
92 for (Int_t i=0;i<3;i++){fPP[i]=fPM[i]=fXr[i]=fAngle[i]=0;}
93 for (Int_t i=0;i<3;i++){fOrder[i]=0;}
94 for (Int_t i=0;i<4;i++){fCausality[i]=0;}
e23730c7 95}
96
d6a49f20 97AliESDv0::AliESDv0(const AliESDv0& v0) :
98 TObject(v0),
99 fOnFlyStatus(v0.fOnFlyStatus),
100 fPdgCode(v0.fPdgCode),
101 fEffMass(v0.fEffMass),
102 fDcaV0Daughters(v0.fDcaV0Daughters),
103 fChi2V0(v0.fChi2V0),
104 fNidx(v0.fNidx),
105 fPidx(v0.fPidx),
106 fParamP(v0.fParamP),
107 fParamN(v0.fParamN),
108 fID(v0.fID),
109 fDist1(v0.fDist1),
110 fDist2(v0.fDist2),
111 fRr(v0.fRr),
112 fStatus(v0.fStatus),
113 fRow0(v0.fRow0),
114 fDistNorm(v0.fDistNorm),
115 fDistSigma(v0.fDistSigma),
116 fChi2Before(v0.fChi2Before),
117 fNBefore(v0.fNBefore),
118 fChi2After(v0.fChi2After),
119 fNAfter(v0.fNAfter),
120 fPointAngleFi(v0.fPointAngleFi),
121 fPointAngleTh(v0.fPointAngleTh),
122 fPointAngle(v0.fPointAngle)
c028b974 123{
d6a49f20 124 //--------------------------------------------------------------------
125 // The copy constructor
126 //--------------------------------------------------------------------
c028b974 127
128 for (int i=0; i<3; i++) {
d6a49f20 129 fPos[i] = v0.fPos[i];
130 fNmom[i] = v0.fNmom[i];
131 fPmom[i] = v0.fPmom[i];
c028b974 132 }
133 for (int i=0; i<6; i++) {
d6a49f20 134 fPosCov[i] = v0.fPosCov[i];
135 fNmomCov[i] = v0.fNmomCov[i];
136 fPmomCov[i] = v0.fPmomCov[i];
c028b974 137 }
c028b974 138
d6a49f20 139 for (Int_t i=0;i<5;i++){
140 fRP[i]=v0.fRP[i];
141 fRM[i]=v0.fRM[i];
142 }
143 for (Int_t i=0; i<2; i++) {
144 fLab[i]=v0.fLab[i];
145 fIndex[i]=v0.fIndex[i];
146 fNormDCAPrim[i]=v0.fNormDCAPrim[i];
147 }
148 for (Int_t i=0;i<6;i++){
149 fClusters[0][i]=v0.fClusters[0][i];
150 fClusters[1][i]=v0.fClusters[1][i];
151 }
152 for (Int_t i=0;i<3;i++){
153 fPP[i]=v0.fPP[i];
154 fPM[i]=v0.fPM[i];
155 fXr[i]=v0.fXr[i];
156 fAngle[i]=v0.fAngle[i];
157 fOrder[i]=v0.fOrder[i];
c028b974 158 }
d6a49f20 159 for (Int_t i=0;i<4;i++){fCausality[i]=v0.fCausality[i];}
c028b974 160}
161
c7bafca9 162AliESDv0::AliESDv0(const AliExternalTrackParam &t1, Int_t i1,
163 const AliExternalTrackParam &t2, Int_t i2) :
164 TObject(),
d6a49f20 165 fOnFlyStatus(kFALSE),
c7bafca9 166 fPdgCode(kK0Short),
167 fEffMass(TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass()),
c028b974 168 fDcaV0Daughters(0),
169 fChi2V0(1.e+33),
c7bafca9 170 fNidx(i1),
d6a49f20 171 fPidx(i2),
172 fParamP(),
173 fParamN(),
174 fID(0),
175 fDist1(-1),
176 fDist2(-1),
177 fRr(-1),
178 fStatus(0),
179 fRow0(-1),
180 fDistNorm(0),
181 fDistSigma(0),
182 fChi2Before(0),
183 fNBefore(0),
184 fChi2After(0),
185 fNAfter(0),
186 fPointAngleFi(0),
187 fPointAngleTh(0),
188 fPointAngle(0)
c7bafca9 189{
190 //--------------------------------------------------------------------
191 // Main constructor (K0s)
192 //--------------------------------------------------------------------
193
194 for (Int_t i=0; i<6; i++) {
195 fPosCov[i]= 0.;
196 fNmomCov[i] = 0.;
197 fPmomCov[i] = 0.;
198 }
199
200 //Trivial estimation of the vertex parameters
201 Double_t x=t1.GetX(), alpha=t1.GetAlpha();
202 const Double_t *par=t1.GetParameter();
203 Double_t pt=1./TMath::Abs(par[4]),
204 phi=TMath::ASin(par[2]) + alpha,
205 cs=TMath::Cos(alpha), sn=TMath::Sin(alpha);
206
207 Double_t px1=pt*TMath::Cos(phi), py1=pt*TMath::Sin(phi), pz1=pt*par[3];
208 Double_t x1=x*cs - par[0]*sn;
209 Double_t y1=x*sn + par[0]*cs;
210 Double_t z1=par[1];
37919f0b 211 const Double_t ss=0.0005*0.0005;//a kind of a residual misalignment precision
212 Double_t sx1=sn*sn*t1.GetSigmaY2()+ss, sy1=cs*cs*t1.GetSigmaY2()+ss;
c7bafca9 213
214
215
216 x=t2.GetX(); alpha=t2.GetAlpha(); par=t2.GetParameter();
217 pt=1./TMath::Abs(par[4]);
218 phi=TMath::ASin(par[2]) + alpha;
219 cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
220
221 Double_t px2=pt*TMath::Cos(phi), py2=pt*TMath::Sin(phi), pz2=pt*par[3];
222 Double_t x2=x*cs - par[0]*sn;
223 Double_t y2=x*sn + par[0]*cs;
224 Double_t z2=par[1];
37919f0b 225 Double_t sx2=sn*sn*t2.GetSigmaY2()+ss, sy2=cs*cs*t2.GetSigmaY2()+ss;
c7bafca9 226
227 Double_t sz1=t1.GetSigmaZ2(), sz2=t2.GetSigmaZ2();
37919f0b 228 Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
229 Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
c7bafca9 230 Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
231 fPos[0]=wx1*x1 + wx2*x2; fPos[1]=wy1*y1 + wy2*y2; fPos[2]=wz1*z1 + wz2*z2;
232
233 //fPos[0]=0.5*(x1+x2); fPos[1]=0.5*(y1+y2); fPos[2]=0.5*(z1+z2);
234 fNmom[0]=px1; fNmom[1]=py1; fNmom[2]=pz1;
235 fPmom[0]=px2; fPmom[1]=py2; fPmom[2]=pz2;
236
237 Double_t e1=TMath::Sqrt(0.13957*0.13957 + px1*px1 + py1*py1 + pz1*pz1);
238 Double_t e2=TMath::Sqrt(0.13957*0.13957 + px2*px2 + py2*py2 + pz2*pz2);
239 fEffMass=TMath::Sqrt((e1+e2)*(e1+e2)-
240 (px1+px2)*(px1+px2)-(py1+py2)*(py1+py2)-(pz1+pz2)*(pz1+pz2));
241
c028b974 242 fChi2V0=7.;
243
244}
c7bafca9 245
c028b974 246AliESDv0::~AliESDv0(){
247 //--------------------------------------------------------------------
248 // Empty destructor
249 //--------------------------------------------------------------------
c7bafca9 250}
251
c028b974 252
253
e23730c7 254Double_t AliESDv0::ChangeMassHypothesis(Int_t code) {
255 //--------------------------------------------------------------------
256 // This function changes the mass hypothesis for this V0
257 // and returns the "kinematical quality" of this hypothesis
258 //--------------------------------------------------------------------
259 Double_t nmass=0.13957, pmass=0.13957, mass=0.49767, ps=0.206;
260
261 fPdgCode=code;
262
263 switch (code) {
264 case kLambda0:
265 nmass=0.13957; pmass=0.93827; mass=1.1157; ps=0.101; break;
266 case kLambda0Bar:
267 pmass=0.13957; nmass=0.93827; mass=1.1157; ps=0.101; break;
268 case kK0Short:
269 break;
270 default:
5f7789fc 271 AliError("invalide PDG code ! Assuming K0s...");
e23730c7 272 fPdgCode=kK0Short;
273 break;
274 }
275
276 Double_t pxn=fNmom[0], pyn=fNmom[1], pzn=fNmom[2];
277 Double_t pxp=fPmom[0], pyp=fPmom[1], pzp=fPmom[2];
278
279 Double_t en=TMath::Sqrt(nmass*nmass + pxn*pxn + pyn*pyn + pzn*pzn);
280 Double_t ep=TMath::Sqrt(pmass*pmass + pxp*pxp + pyp*pyp + pzp*pzp);
281 Double_t pxl=pxn+pxp, pyl=pyn+pyp, pzl=pzn+pzp;
282 Double_t pl=TMath::Sqrt(pxl*pxl + pyl*pyl + pzl*pzl);
283
284 fEffMass=TMath::Sqrt((en+ep)*(en+ep)-pl*pl);
285
286 Double_t beta=pl/(en+ep);
287 Double_t pln=(pxn*pxl + pyn*pyl + pzn*pzl)/pl;
288 Double_t plp=(pxp*pxl + pyp*pyl + pzp*pzl)/pl;
289
290 Double_t pt2=pxp*pxp + pyp*pyp + pzp*pzp - plp*plp;
291
292 Double_t a=(plp-pln)/(plp+pln);
293 a -= (pmass*pmass-nmass*nmass)/(mass*mass);
294 a = 0.25*beta*beta*mass*mass*a*a + pt2;
295
296 return (a - ps*ps);
297
298}
299
300void AliESDv0::GetPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const {
301 //--------------------------------------------------------------------
302 // This function returns V0's momentum (global)
303 //--------------------------------------------------------------------
304 px=fNmom[0]+fPmom[0];
305 py=fNmom[1]+fPmom[1];
306 pz=fNmom[2]+fPmom[2];
307}
308
309void AliESDv0::GetXYZ(Double_t &x, Double_t &y, Double_t &z) const {
310 //--------------------------------------------------------------------
311 // This function returns V0's position (global)
312 //--------------------------------------------------------------------
313 x=fPos[0];
314 y=fPos[1];
315 z=fPos[2];
316}
317
318Double_t AliESDv0::GetD(Double_t x0, Double_t y0, Double_t z0) const {
319 //--------------------------------------------------------------------
320 // This function returns V0's impact parameter
321 //--------------------------------------------------------------------
322 Double_t x=fPos[0],y=fPos[1],z=fPos[2];
323 Double_t px=fNmom[0]+fPmom[0];
324 Double_t py=fNmom[1]+fPmom[1];
325 Double_t pz=fNmom[2]+fPmom[2];
326
327 Double_t dx=(y0-y)*pz - (z0-z)*py;
328 Double_t dy=(x0-x)*pz - (z0-z)*px;
329 Double_t dz=(x0-x)*py - (y0-y)*px;
330 Double_t d=TMath::Sqrt((dx*dx+dy*dy+dz*dz)/(px*px+py*py+pz*pz));
331 return d;
332}
c028b974 333
334
335Double_t AliESDv0::GetV0CosineOfPointingAngle(Double_t& refPointX, Double_t& refPointY, Double_t& refPointZ) const {
336 // calculates the pointing angle of the V0 wrt a reference point
337
338 Double_t momV0[3]; //momentum of the V0
339 GetPxPyPz(momV0[0],momV0[1],momV0[2]);
340
341 Double_t deltaPos[3]; //vector between the reference point and the V0 vertex
342 deltaPos[0] = fPos[0] - refPointX;
343 deltaPos[1] = fPos[1] - refPointY;
344 deltaPos[2] = fPos[2] - refPointZ;
345
346 Double_t momV02 = momV0[0]*momV0[0] + momV0[1]*momV0[1] + momV0[2]*momV0[2];
347 Double_t deltaPos2 = deltaPos[0]*deltaPos[0] + deltaPos[1]*deltaPos[1] + deltaPos[2]*deltaPos[2];
348
349 Double_t cosinePointingAngle = (deltaPos[0]*momV0[0] +
350 deltaPos[1]*momV0[1] +
351 deltaPos[2]*momV0[2] ) /
352 TMath::Sqrt(momV02 * deltaPos2);
353
354 return cosinePointingAngle;
355}
d6a49f20 356
357
358// **** The following functions need to be revised
359
360Double_t AliESDv0::GetSigmaY(){
361 //
362 // return sigmay in y at vertex position using covariance matrix
363 //
364 const Double_t * cp = fParamP.GetCovariance();
365 const Double_t * cm = fParamN.GetCovariance();
366 Double_t sigmay = cp[0]+cm[0]+ cp[5]*(fParamP.GetX()-fRr)*(fParamP.GetX()-fRr)+ cm[5]*(fParamN.GetX()-fRr)*(fParamN.GetX()-fRr);
367 return (sigmay>0) ? TMath::Sqrt(sigmay):100;
368}
369
370Double_t AliESDv0::GetSigmaZ(){
371 //
372 // return sigmay in y at vertex position using covariance matrix
373 //
374 const Double_t * cp = fParamP.GetCovariance();
375 const Double_t * cm = fParamN.GetCovariance();
376 Double_t sigmaz = cp[2]+cm[2]+ cp[9]*(fParamP.GetX()-fRr)*(fParamP.GetX()-fRr)+ cm[9]*(fParamN.GetX()-fRr)*(fParamN.GetX()-fRr);
377 return (sigmaz>0) ? TMath::Sqrt(sigmaz):100;
378}
379
380Double_t AliESDv0::GetSigmaD0(){
381 //
382 // Sigma parameterization using covariance matrix
383 //
384 // sigma of distance between two tracks in vertex position
385 // sigma of DCA is proportianal to sigmaD0
386 // factor 2 difference is explained by the fact that the DCA is calculated at the position
387 // where the tracks as closest together ( not exact position of the vertex)
388 //
389 const Double_t * cp = fParamP.GetCovariance();
390 const Double_t * cm = fParamN.GetCovariance();
391 Double_t sigmaD0 = cp[0]+cm[0]+cp[2]+cm[2]+fgkParams.fPSigmaOffsetD0*fgkParams.fPSigmaOffsetD0;
392 sigmaD0 += ((fParamP.GetX()-fRr)*(fParamP.GetX()-fRr))*(cp[5]+cp[9]);
393 sigmaD0 += ((fParamN.GetX()-fRr)*(fParamN.GetX()-fRr))*(cm[5]+cm[9]);
394 return (sigmaD0>0)? TMath::Sqrt(sigmaD0):100;
395}
396
397
398Double_t AliESDv0::GetSigmaAP0(){
399 //
400 //Sigma parameterization using covariance matrices
401 //
402 Double_t prec = TMath::Sqrt((fPM[0]+fPP[0])*(fPM[0]+fPP[0])
403 +(fPM[1]+fPP[1])*(fPM[1]+fPP[1])
404 +(fPM[2]+fPP[2])*(fPM[2]+fPP[2]));
405 Double_t normp = TMath::Sqrt(fPP[0]*fPP[0]+fPP[1]*fPP[1]+fPP[2]*fPP[2])/prec; // fraction of the momenta
406 Double_t normm = TMath::Sqrt(fPM[0]*fPM[0]+fPM[1]*fPM[1]+fPM[2]*fPM[2])/prec;
407 const Double_t * cp = fParamP.GetCovariance();
408 const Double_t * cm = fParamN.GetCovariance();
409 Double_t sigmaAP0 = fgkParams.fPSigmaOffsetAP0*fgkParams.fPSigmaOffsetAP0; // minimal part
410 sigmaAP0 += (cp[5]+cp[9])*(normp*normp)+(cm[5]+cm[9])*(normm*normm); // angular resolution part
411 Double_t sigmaAP1 = GetSigmaD0()/(TMath::Abs(fRr)+0.01); // vertex position part
412 sigmaAP0 += 0.5*sigmaAP1*sigmaAP1;
413 return (sigmaAP0>0)? TMath::Sqrt(sigmaAP0):100;
414}
415
416Double_t AliESDv0::GetEffectiveSigmaD0(){
417 //
418 // minimax - effective Sigma parameterization
419 // p12 effective curvature and v0 radius postion used as parameters
420 //
421 Double_t p12 = TMath::Sqrt(fParamP.GetParameter()[4]*fParamP.GetParameter()[4]+
422 fParamN.GetParameter()[4]*fParamN.GetParameter()[4]);
423 Double_t sigmaED0= TMath::Max(TMath::Sqrt(fRr)-fgkParams.fPSigmaRminDE,0.0)*fgkParams.fPSigmaCoefDE*p12*p12;
424 sigmaED0*= sigmaED0;
425 sigmaED0*= sigmaED0;
426 sigmaED0 = TMath::Sqrt(sigmaED0+fgkParams.fPSigmaOffsetDE*fgkParams.fPSigmaOffsetDE);
427 return (sigmaED0<fgkParams.fPSigmaMaxDE) ? sigmaED0: fgkParams.fPSigmaMaxDE;
428}
429
430
431Double_t AliESDv0::GetEffectiveSigmaAP0(){
432 //
433 // effective Sigma parameterization of point angle resolution
434 //
435 Double_t p12 = TMath::Sqrt(fParamP.GetParameter()[4]*fParamP.GetParameter()[4]+
436 fParamN.GetParameter()[4]*fParamN.GetParameter()[4]);
437 Double_t sigmaAPE= fgkParams.fPSigmaBase0APE;
438 sigmaAPE+= fgkParams.fPSigmaR0APE/(fgkParams.fPSigmaR1APE+fRr);
439 sigmaAPE*= (fgkParams.fPSigmaP0APE+fgkParams.fPSigmaP1APE*p12);
440 sigmaAPE = TMath::Min(sigmaAPE,fgkParams.fPSigmaMaxAPE);
441 return sigmaAPE;
442}
443
444
445Double_t AliESDv0::GetMinimaxSigmaAP0(){
446 //
447 // calculate mini-max effective sigma of point angle resolution
448 //
449 //compv0->fTree->SetAlias("SigmaAP2","max(min((SigmaAP0+SigmaAPE0)*0.5,1.5*SigmaAPE0),0.5*SigmaAPE0+0.003)");
450 Double_t effectiveSigma = GetEffectiveSigmaAP0();
451 Double_t sigmaMMAP = 0.5*(GetSigmaAP0()+effectiveSigma);
452 sigmaMMAP = TMath::Min(sigmaMMAP, fgkParams.fPMaxFractionAP0*effectiveSigma);
453 sigmaMMAP = TMath::Max(sigmaMMAP, fgkParams.fPMinFractionAP0*effectiveSigma+fgkParams.fPMinAP0);
454 return sigmaMMAP;
455}
456Double_t AliESDv0::GetMinimaxSigmaD0(){
457 //
458 // calculate mini-max sigma of dca resolution
459 //
460 //compv0->fTree->SetAlias("SigmaD2","max(min((SigmaD0+SigmaDE0)*0.5,1.5*SigmaDE0),0.5*SigmaDE0)");
461 Double_t effectiveSigma = GetEffectiveSigmaD0();
462 Double_t sigmaMMD0 = 0.5*(GetSigmaD0()+effectiveSigma);
463 sigmaMMD0 = TMath::Min(sigmaMMD0, fgkParams.fPMaxFractionD0*effectiveSigma);
464 sigmaMMD0 = TMath::Max(sigmaMMD0, fgkParams.fPMinFractionD0*effectiveSigma+fgkParams.fPMinD0);
465 return sigmaMMD0;
466}
467
468
469Double_t AliESDv0::GetLikelihoodAP(Int_t mode0, Int_t mode1){
470 //
471 // get likelihood for point angle
472 //
473 Double_t sigmaAP = 0.007; //default sigma
474 switch (mode0){
475 case 0:
476 sigmaAP = GetSigmaAP0(); // mode 0 - covariance matrix estimates used
477 break;
478 case 1:
479 sigmaAP = GetEffectiveSigmaAP0(); // mode 1 - effective sigma used
480 break;
481 case 2:
482 sigmaAP = GetMinimaxSigmaAP0(); // mode 2 - minimax sigma
483 break;
484 }
485 Double_t apNorm = TMath::Min(TMath::ACos(fPointAngle)/sigmaAP,50.);
486 //normalized point angle, restricted - because of overflow problems in Exp
487 Double_t likelihood = 0;
488 switch(mode1){
489 case 0:
490 likelihood = TMath::Exp(-0.5*apNorm*apNorm);
491 // one component
492 break;
493 case 1:
494 likelihood = (TMath::Exp(-0.5*apNorm*apNorm)+0.5* TMath::Exp(-0.25*apNorm*apNorm))/1.5;
495 // two components
496 break;
497 case 2:
498 likelihood = (TMath::Exp(-0.5*apNorm*apNorm)+0.5* TMath::Exp(-0.25*apNorm*apNorm)+0.25*TMath::Exp(-0.125*apNorm*apNorm))/1.75;
499 // three components
500 break;
501 }
502 return likelihood;
503}
504
505Double_t AliESDv0::GetLikelihoodD(Int_t mode0, Int_t mode1){
506 //
507 // get likelihood for DCA
508 //
509 Double_t sigmaD = 0.03; //default sigma
510 switch (mode0){
511 case 0:
512 sigmaD = GetSigmaD0(); // mode 0 - covariance matrix estimates used
513 break;
514 case 1:
515 sigmaD = GetEffectiveSigmaD0(); // mode 1 - effective sigma used
516 break;
517 case 2:
518 sigmaD = GetMinimaxSigmaD0(); // mode 2 - minimax sigma
519 break;
520 }
521 Double_t dNorm = TMath::Min(fDist2/sigmaD,50.);
522 //normalized point angle, restricted - because of overflow problems in Exp
523 Double_t likelihood = 0;
524 switch(mode1){
525 case 0:
526 likelihood = TMath::Exp(-2.*dNorm);
527 // one component
528 break;
529 case 1:
530 likelihood = (TMath::Exp(-2.*dNorm)+0.5* TMath::Exp(-dNorm))/1.5;
531 // two components
532 break;
533 case 2:
534 likelihood = (TMath::Exp(-2.*dNorm)+0.5* TMath::Exp(-dNorm)+0.25*TMath::Exp(-0.5*dNorm))/1.75;
535 // three components
536 break;
537 }
538 return likelihood;
539
540}
541
542Double_t AliESDv0::GetLikelihoodC(Int_t mode0, Int_t /*mode1*/){
543 //
544 // get likelihood for Causality
545 // !!! Causality variables defined in AliITStrackerMI !!!
546 // when more information was available
547 //
548 Double_t likelihood = 0.5;
549 Double_t minCausal = TMath::Min(fCausality[0],fCausality[1]);
550 Double_t maxCausal = TMath::Max(fCausality[0],fCausality[1]);
551 // minCausal = TMath::Max(minCausal,0.5*maxCausal);
552 //compv0->fTree->SetAlias("LCausal","(1.05-(2*(0.8-exp(-max(RC.fV0rec.fCausality[0],RC.fV0rec.fCausality[1])))+2*(0.8-exp(-min(RC.fV0rec.fCausality[0],RC.fV0rec.fCausality[1]))))/2)**4");
553
554 switch(mode0){
555 case 0:
556 //normalization
557 likelihood = TMath::Power((1.05-2*(0.8-TMath::Exp(-maxCausal))),4.);
558 break;
559 case 1:
560 likelihood = TMath::Power(1.05-(2*(0.8-TMath::Exp(-maxCausal))+(2*(0.8-TMath::Exp(-minCausal))))*0.5,4.);
561 break;
562 }
563 return likelihood;
564
565}
566
567void AliESDv0::SetCausality(Float_t pb0, Float_t pb1, Float_t pa0, Float_t pa1)
568{
569 //
570 // set probabilities
571 //
572 fCausality[0] = pb0; // probability - track 0 exist before vertex
573 fCausality[1] = pb1; // probability - track 1 exist before vertex
574 fCausality[2] = pa0; // probability - track 0 exist close after vertex
575 fCausality[3] = pa1; // probability - track 1 exist close after vertex
576}
577void AliESDv0::SetClusters(Int_t *clp, Int_t *clm)
578{
579 //
580 // Set its clusters indexes
581 //
582 for (Int_t i=0;i<6;i++) fClusters[0][i] = clp[i];
583 for (Int_t i=0;i<6;i++) fClusters[1][i] = clm[i];
584}
585
586
587void AliESDv0::SetP(const AliExternalTrackParam & paramp) {
588 //
589 // set track +
590 //
591 fParamP = paramp;
592}
593
594void AliESDv0::SetM(const AliExternalTrackParam & paramm){
595 //
596 //set track -
597 //
598 fParamN = paramm;
599}
600
601void AliESDv0::SetRp(const Double_t *rp){
602 //
603 // set pid +
604 //
605 for (Int_t i=0;i<5;i++) fRP[i]=rp[i];
606}
607
608void AliESDv0::SetRm(const Double_t *rm){
609 //
610 // set pid -
611 //
612 for (Int_t i=0;i<5;i++) fRM[i]=rm[i];
613}
614
615
616void AliESDv0::UpdatePID(Double_t pidp[5], Double_t pidm[5])
617{
618 //
619 // set PID hypothesy
620 //
621 // norm PID to 1
622 Float_t sump =0;
623 Float_t summ =0;
624 for (Int_t i=0;i<5;i++){
625 fRP[i]=pidp[i];
626 sump+=fRP[i];
627 fRM[i]=pidm[i];
628 summ+=fRM[i];
629 }
630 for (Int_t i=0;i<5;i++){
631 fRP[i]/=sump;
632 fRM[i]/=summ;
633 }
634}
635
636Float_t AliESDv0::GetProb(UInt_t p1, UInt_t p2){
637 //
638 //
639 //
640 //
641 return TMath::Max(fRP[p1]+fRM[p2], fRP[p2]+fRM[p1]);
642}
643
644Float_t AliESDv0::GetEffMass(UInt_t p1, UInt_t p2){
645 //
646 // calculate effective mass
647 //
648 const Float_t kpmass[5] = {5.10000000000000037e-04,1.05660000000000004e-01,1.39570000000000000e-01,
649 4.93599999999999983e-01, 9.38270000000000048e-01};
650 if (p1>4) return -1;
651 if (p2>4) return -1;
652 Float_t mass1 = kpmass[p1];
653 Float_t mass2 = kpmass[p2];
654 Double_t *m1 = fPP;
655 Double_t *m2 = fPM;
656 //
657 //if (fRP[p1]+fRM[p2]<fRP[p2]+fRM[p1]){
658 // m1 = fPM;
659 // m2 = fPP;
660 //}
661 //
662 Float_t e1 = TMath::Sqrt(mass1*mass1+
663 m1[0]*m1[0]+
664 m1[1]*m1[1]+
665 m1[2]*m1[2]);
666 Float_t e2 = TMath::Sqrt(mass2*mass2+
667 m2[0]*m2[0]+
668 m2[1]*m2[1]+
669 m2[2]*m2[2]);
670 Float_t mass =
671 (m2[0]+m1[0])*(m2[0]+m1[0])+
672 (m2[1]+m1[1])*(m2[1]+m1[1])+
673 (m2[2]+m1[2])*(m2[2]+m1[2]);
674
675 mass = TMath::Sqrt((e1+e2)*(e1+e2)-mass);
676 return mass;
677}
678