]> 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
Protection against div. by 0 in the Set(xyz,p..) for tracks with momentum along X...
[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 <TMath.h>
90e48c0c 29#include <TDatabasePDG.h>
90e48c0c 30#include <TParticlePDG.h>
3e27a419 31#include <TVector3.h>
e23730c7 32
5f7789fc 33#include "AliLog.h"
e23730c7 34#include "AliESDv0.h"
97135b34 35#include "AliESDV0Params.h"
e23730c7 36
37ClassImp(AliESDv0)
38
562dd0b4 39const AliESDV0Params AliESDv0::fgkParams;
d6a49f20 40
90e48c0c 41AliESDv0::AliESDv0() :
646c9704 42 AliVParticle(),
d6a49f20 43 fParamN(),
b75d63a7 44 fParamP(),
562dd0b4 45 fEffMass(TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass()),
46 fDcaV0Daughters(0),
47 fChi2V0(0.),
48 fRr(0),
d6a49f20 49 fDistSigma(0),
50 fChi2Before(0),
d6a49f20 51 fChi2After(0),
d6a49f20 52 fPointAngleFi(0),
562dd0b4 53 fPointAngleTh(0),
54 fPointAngle(0),
55 fPdgCode(kK0Short),
56 fNidx(0),
57 fPidx(0),
58 fStatus(0),
59 fNBefore(0),
60 fNAfter(0),
61 fOnFlyStatus(kFALSE)
90e48c0c 62{
e23730c7 63 //--------------------------------------------------------------------
64 // Default constructor (K0s)
65 //--------------------------------------------------------------------
6605de26 66
67 for (Int_t i=0; i<3; i++) {
68 fPos[i] = 0.;
69 fNmom[i] = 0.;
70 fPmom[i] = 0.;
71 }
72
73 for (Int_t i=0; i<6; i++) {
74 fPosCov[i]= 0.;
6605de26 75 }
d6a49f20 76
d6a49f20 77 for (Int_t i=0;i<6;i++){fClusters[0][i]=0; fClusters[1][i]=0;}
78 fNormDCAPrim[0]=fNormDCAPrim[1]=0;
b75d63a7 79 for (Int_t i=0;i<3;i++){fAngle[i]=0;}
d6a49f20 80 for (Int_t i=0;i<4;i++){fCausality[i]=0;}
e23730c7 81}
82
d6a49f20 83AliESDv0::AliESDv0(const AliESDv0& v0) :
646c9704 84 AliVParticle(v0),
562dd0b4 85 fParamN(v0.fParamN),
86 fParamP(v0.fParamP),
d6a49f20 87 fEffMass(v0.fEffMass),
88 fDcaV0Daughters(v0.fDcaV0Daughters),
89 fChi2V0(v0.fChi2V0),
d6a49f20 90 fRr(v0.fRr),
d6a49f20 91 fDistSigma(v0.fDistSigma),
92 fChi2Before(v0.fChi2Before),
d6a49f20 93 fChi2After(v0.fChi2After),
d6a49f20 94 fPointAngleFi(v0.fPointAngleFi),
562dd0b4 95 fPointAngleTh(v0.fPointAngleTh),
96 fPointAngle(v0.fPointAngle),
97 fPdgCode(v0.fPdgCode),
98 fNidx(v0.fNidx),
99 fPidx(v0.fPidx),
100 fStatus(v0.fStatus),
101 fNBefore(v0.fNBefore),
102 fNAfter(v0.fNAfter),
103 fOnFlyStatus(v0.fOnFlyStatus)
c028b974 104{
d6a49f20 105 //--------------------------------------------------------------------
106 // The copy constructor
107 //--------------------------------------------------------------------
c028b974 108
109 for (int i=0; i<3; i++) {
d6a49f20 110 fPos[i] = v0.fPos[i];
111 fNmom[i] = v0.fNmom[i];
112 fPmom[i] = v0.fPmom[i];
c028b974 113 }
114 for (int i=0; i<6; i++) {
d6a49f20 115 fPosCov[i] = v0.fPosCov[i];
c028b974 116 }
c028b974 117
d6a49f20 118 for (Int_t i=0; i<2; i++) {
b75d63a7 119 fNormDCAPrim[i]=v0.fNormDCAPrim[i];
d6a49f20 120 }
121 for (Int_t i=0;i<6;i++){
122 fClusters[0][i]=v0.fClusters[0][i];
123 fClusters[1][i]=v0.fClusters[1][i];
124 }
125 for (Int_t i=0;i<3;i++){
d6a49f20 126 fAngle[i]=v0.fAngle[i];
c028b974 127 }
d6a49f20 128 for (Int_t i=0;i<4;i++){fCausality[i]=v0.fCausality[i];}
c028b974 129}
130
732a24fe 131
c7bafca9 132AliESDv0::AliESDv0(const AliExternalTrackParam &t1, Int_t i1,
133 const AliExternalTrackParam &t2, Int_t i2) :
646c9704 134 AliVParticle(),
562dd0b4 135 fParamN(t1),
136 fParamP(t2),
c7bafca9 137 fEffMass(TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass()),
c028b974 138 fDcaV0Daughters(0),
562dd0b4 139 fChi2V0(0.),
140 fRr(0),
141 fDistSigma(0),
142 fChi2Before(0),
143 fChi2After(0),
144 fPointAngleFi(0),
145 fPointAngleTh(0),
b75d63a7 146 fPointAngle(0),
562dd0b4 147 fPdgCode(kK0Short),
c7bafca9 148 fNidx(i1),
d6a49f20 149 fPidx(i2),
d6a49f20 150 fStatus(0),
d6a49f20 151 fNBefore(0),
d6a49f20 152 fNAfter(0),
562dd0b4 153 fOnFlyStatus(kFALSE)
c7bafca9 154{
155 //--------------------------------------------------------------------
156 // Main constructor (K0s)
157 //--------------------------------------------------------------------
158
ed8f3bb5 159 //Make sure the daughters are ordered (needed for the on-the-fly V0s)
160 Short_t cN=t1.Charge(), cP=t2.Charge();
161 if ((cN>0) && (cN != cP)) {
162 fParamN.~AliExternalTrackParam();
163 new (&fParamN) AliExternalTrackParam(t2);
164 fParamP.~AliExternalTrackParam();
165 new (&fParamP) AliExternalTrackParam(t1);
166 }
167
c7bafca9 168 for (Int_t i=0; i<6; i++) {
169 fPosCov[i]= 0.;
c7bafca9 170 }
171
172 //Trivial estimation of the vertex parameters
b75d63a7 173 Double_t alpha=t1.GetAlpha(), cs=TMath::Cos(alpha), sn=TMath::Sin(alpha);
174 Double_t tmp[3];
175 t1.GetPxPyPz(tmp);
176 Double_t px1=tmp[0], py1=tmp[1], pz1=tmp[2];
177 t1.GetXYZ(tmp);
178 Double_t x1=tmp[0], y1=tmp[1], z1=tmp[2];
37919f0b 179 const Double_t ss=0.0005*0.0005;//a kind of a residual misalignment precision
180 Double_t sx1=sn*sn*t1.GetSigmaY2()+ss, sy1=cs*cs*t1.GetSigmaY2()+ss;
c7bafca9 181
182
b75d63a7 183 alpha=t2.GetAlpha(); cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
184 t2.GetPxPyPz(tmp);
185 Double_t px2=tmp[0], py2=tmp[1], pz2=tmp[2];
186 t2.GetXYZ(tmp);
187 Double_t x2=tmp[0], y2=tmp[1], z2=tmp[2];
37919f0b 188 Double_t sx2=sn*sn*t2.GetSigmaY2()+ss, sy2=cs*cs*t2.GetSigmaY2()+ss;
c7bafca9 189
190 Double_t sz1=t1.GetSigmaZ2(), sz2=t2.GetSigmaZ2();
37919f0b 191 Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
192 Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
c7bafca9 193 Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
194 fPos[0]=wx1*x1 + wx2*x2; fPos[1]=wy1*y1 + wy2*y2; fPos[2]=wz1*z1 + wz2*z2;
195
196 //fPos[0]=0.5*(x1+x2); fPos[1]=0.5*(y1+y2); fPos[2]=0.5*(z1+z2);
197 fNmom[0]=px1; fNmom[1]=py1; fNmom[2]=pz1;
198 fPmom[0]=px2; fPmom[1]=py2; fPmom[2]=pz2;
199
b75d63a7 200 for (Int_t i=0;i<6;i++){fClusters[0][i]=0; fClusters[1][i]=0;}
201 fNormDCAPrim[0]=fNormDCAPrim[1]=0;
202 for (Int_t i=0;i<3;i++){fAngle[i]=0;}
203 for (Int_t i=0;i<4;i++){fCausality[i]=0;}
c028b974 204}
c7bafca9 205
3e27a419 206AliESDv0& AliESDv0::operator=(const AliESDv0 &v0)
207{
732a24fe 208 //--------------------------------------------------------------------
3e27a419 209 // The assignment operator
732a24fe 210 //--------------------------------------------------------------------
211
212 if(this==&v0)return *this;
646c9704 213 AliVParticle::operator=(v0);
732a24fe 214 fParamN = v0.fParamN;
215 fParamP = v0.fParamP;
216 fEffMass = v0.fEffMass;
217 fDcaV0Daughters = v0.fDcaV0Daughters;
218 fChi2V0 = v0.fChi2V0;
219 fRr = v0.fRr;
220 fDistSigma = v0.fDistSigma;
221 fChi2Before = v0.fChi2Before;
222 fChi2After = v0.fChi2After;
223 fPointAngleFi = v0.fPointAngleFi;
224 fPointAngleTh = v0.fPointAngleTh;
225 fPointAngle = v0.fPointAngle;
226 fPdgCode = v0.fPdgCode;
227 fNidx = v0.fNidx;
228 fPidx = v0.fPidx;
229 fStatus = v0.fStatus;
230 fNBefore = v0.fNBefore;
231 fNAfter = v0.fNAfter;
232 fOnFlyStatus = v0.fOnFlyStatus;
233
234 for (int i=0; i<3; i++) {
235 fPos[i] = v0.fPos[i];
236 fNmom[i] = v0.fNmom[i];
237 fPmom[i] = v0.fPmom[i];
238 }
239 for (int i=0; i<6; i++) {
240 fPosCov[i] = v0.fPosCov[i];
241 }
242 for (Int_t i=0; i<2; i++) {
243 fNormDCAPrim[i]=v0.fNormDCAPrim[i];
244 }
245 for (Int_t i=0;i<6;i++){
246 fClusters[0][i]=v0.fClusters[0][i];
247 fClusters[1][i]=v0.fClusters[1][i];
248 }
249 for (Int_t i=0;i<3;i++){
250 fAngle[i]=v0.fAngle[i];
251 }
252 for (Int_t i=0;i<4;i++){fCausality[i]=v0.fCausality[i];}
253
254 return *this;
255}
256
257void AliESDv0::Copy(TObject& obj) const {
258
3e27a419 259 // this overwrites the virtual TOBject::Copy()
732a24fe 260 // to allow run time copying without casting
261 // in AliESDEvent
262
263 if(this==&obj)return;
264 AliESDv0 *robj = dynamic_cast<AliESDv0*>(&obj);
265 if(!robj)return; // not an aliesesv0
266 *robj = *this;
267}
268
c028b974 269AliESDv0::~AliESDv0(){
270 //--------------------------------------------------------------------
271 // Empty destructor
272 //--------------------------------------------------------------------
c7bafca9 273}
274
646c9704 275// Start with AliVParticle functions
276Double_t AliESDv0::E() const {
277 //--------------------------------------------------------------------
278 // This gives the energy assuming the ChangeMassHypothesis was called
279 //--------------------------------------------------------------------
3e27a419 280 return E(fPdgCode);
281}
282
283Double_t AliESDv0::Y() const {
284 //--------------------------------------------------------------------
285 // This gives the energy assuming the ChangeMassHypothesis was called
286 //--------------------------------------------------------------------
287 return Y(fPdgCode);
288}
289
290// Then extend AliVParticle functions
291Double_t AliESDv0::E(Int_t pdg) const {
292 //--------------------------------------------------------------------
293 // This gives the energy with the particle hypothesis as argument
294 //--------------------------------------------------------------------
295 Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
646c9704 296 return TMath::Sqrt(mass*mass+P()*P());
297}
c028b974 298
3e27a419 299Double_t AliESDv0::Y(Int_t pdg) const {
300 //--------------------------------------------------------------------
301 // This gives the rapidity with the particle hypothesis as argument
302 //--------------------------------------------------------------------
303 return 0.5*TMath::Log((E(pdg)+Pz())/(E(pdg)-Pz()+1.e-13));
304}
305
306// Now the functions for analysis consistency
307Double_t AliESDv0::RapK0Short() const {
308 //--------------------------------------------------------------------
309 // This gives the pseudorapidity assuming a K0s particle
310 //--------------------------------------------------------------------
311 return Y(kK0Short);
312}
313
314Double_t AliESDv0::RapLambda() const {
315 //--------------------------------------------------------------------
316 // This gives the pseudorapidity assuming a (Anti) Lambda particle
317 //--------------------------------------------------------------------
318 return Y(kLambda0);
319}
320
321Double_t AliESDv0::AlphaV0() const {
322 //--------------------------------------------------------------------
323 // This gives the Armenteros-Podolanski alpha
324 //--------------------------------------------------------------------
325 TVector3 momNeg(fNmom[0],fNmom[1],fNmom[2]);
326 TVector3 momPos(fPmom[0],fPmom[1],fPmom[2]);
327 TVector3 momTot(Px(),Py(),Pz());
328
97135b34 329 Double_t lQlNeg = momNeg.Dot(momTot)/momTot.Mag();
ad07c38b 330 Double_t lQlPos = momPos.Dot(momTot)/momTot.Mag();
3e27a419 331
bc88fc12 332 //return 1.-2./(1.+lQlNeg/lQlPos);
333 return (lQlPos - lQlNeg)/(lQlPos + lQlNeg);
3e27a419 334}
335
336Double_t AliESDv0::PtArmV0() const {
337 //--------------------------------------------------------------------
338 // This gives the Armenteros-Podolanski ptarm
339 //--------------------------------------------------------------------
340 TVector3 momNeg(fNmom[0],fNmom[1],fNmom[2]);
341 TVector3 momTot(Px(),Py(),Pz());
342
343 return momNeg.Perp(momTot);
344}
345
346// Eventually the older functions
e23730c7 347Double_t AliESDv0::ChangeMassHypothesis(Int_t code) {
348 //--------------------------------------------------------------------
349 // This function changes the mass hypothesis for this V0
350 // and returns the "kinematical quality" of this hypothesis
351 //--------------------------------------------------------------------
b75d63a7 352 static
353 Double_t piMass=TDatabasePDG::Instance()->GetParticle(kPiPlus)->Mass();
354 static
355 Double_t prMass=TDatabasePDG::Instance()->GetParticle(kProton)->Mass();
356 static
357 Double_t k0Mass=TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass();
358 static
359 Double_t l0Mass=TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass();
360
361 Double_t nmass=piMass, pmass=piMass, mass=k0Mass, ps=0.206;
e23730c7 362
363 fPdgCode=code;
364
365 switch (code) {
366 case kLambda0:
b75d63a7 367 nmass=piMass; pmass=prMass; mass=l0Mass; ps=0.101; break;
e23730c7 368 case kLambda0Bar:
b75d63a7 369 pmass=piMass; nmass=prMass; mass=l0Mass; ps=0.101; break;
e23730c7 370 case kK0Short:
371 break;
372 default:
5f7789fc 373 AliError("invalide PDG code ! Assuming K0s...");
e23730c7 374 fPdgCode=kK0Short;
375 break;
376 }
377
378 Double_t pxn=fNmom[0], pyn=fNmom[1], pzn=fNmom[2];
379 Double_t pxp=fPmom[0], pyp=fPmom[1], pzp=fPmom[2];
380
381 Double_t en=TMath::Sqrt(nmass*nmass + pxn*pxn + pyn*pyn + pzn*pzn);
382 Double_t ep=TMath::Sqrt(pmass*pmass + pxp*pxp + pyp*pyp + pzp*pzp);
383 Double_t pxl=pxn+pxp, pyl=pyn+pyp, pzl=pzn+pzp;
384 Double_t pl=TMath::Sqrt(pxl*pxl + pyl*pyl + pzl*pzl);
385
386 fEffMass=TMath::Sqrt((en+ep)*(en+ep)-pl*pl);
387
388 Double_t beta=pl/(en+ep);
389 Double_t pln=(pxn*pxl + pyn*pyl + pzn*pzl)/pl;
390 Double_t plp=(pxp*pxl + pyp*pyl + pzp*pzl)/pl;
391
392 Double_t pt2=pxp*pxp + pyp*pyp + pzp*pzp - plp*plp;
393
394 Double_t a=(plp-pln)/(plp+pln);
395 a -= (pmass*pmass-nmass*nmass)/(mass*mass);
396 a = 0.25*beta*beta*mass*mass*a*a + pt2;
397
398 return (a - ps*ps);
399
400}
401
402void AliESDv0::GetPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const {
403 //--------------------------------------------------------------------
404 // This function returns V0's momentum (global)
405 //--------------------------------------------------------------------
406 px=fNmom[0]+fPmom[0];
407 py=fNmom[1]+fPmom[1];
408 pz=fNmom[2]+fPmom[2];
409}
410
411void AliESDv0::GetXYZ(Double_t &x, Double_t &y, Double_t &z) const {
412 //--------------------------------------------------------------------
413 // This function returns V0's position (global)
414 //--------------------------------------------------------------------
415 x=fPos[0];
416 y=fPos[1];
417 z=fPos[2];
418}
419
b75d63a7 420Float_t AliESDv0::GetD(Double_t x0, Double_t y0, Double_t z0) const {
e23730c7 421 //--------------------------------------------------------------------
422 // This function returns V0's impact parameter
423 //--------------------------------------------------------------------
424 Double_t x=fPos[0],y=fPos[1],z=fPos[2];
425 Double_t px=fNmom[0]+fPmom[0];
426 Double_t py=fNmom[1]+fPmom[1];
427 Double_t pz=fNmom[2]+fPmom[2];
428
429 Double_t dx=(y0-y)*pz - (z0-z)*py;
430 Double_t dy=(x0-x)*pz - (z0-z)*px;
431 Double_t dz=(x0-x)*py - (y0-y)*px;
432 Double_t d=TMath::Sqrt((dx*dx+dy*dy+dz*dz)/(px*px+py*py+pz*pz));
433 return d;
434}
c028b974 435
97135b34 436Float_t AliESDv0::GetV0CosineOfPointingAngle(Double_t refPointX, Double_t refPointY, Double_t refPointZ) const {
c028b974 437 // calculates the pointing angle of the V0 wrt a reference point
438
439 Double_t momV0[3]; //momentum of the V0
440 GetPxPyPz(momV0[0],momV0[1],momV0[2]);
441
442 Double_t deltaPos[3]; //vector between the reference point and the V0 vertex
443 deltaPos[0] = fPos[0] - refPointX;
444 deltaPos[1] = fPos[1] - refPointY;
445 deltaPos[2] = fPos[2] - refPointZ;
446
447 Double_t momV02 = momV0[0]*momV0[0] + momV0[1]*momV0[1] + momV0[2]*momV0[2];
448 Double_t deltaPos2 = deltaPos[0]*deltaPos[0] + deltaPos[1]*deltaPos[1] + deltaPos[2]*deltaPos[2];
449
450 Double_t cosinePointingAngle = (deltaPos[0]*momV0[0] +
451 deltaPos[1]*momV0[1] +
452 deltaPos[2]*momV0[2] ) /
453 TMath::Sqrt(momV02 * deltaPos2);
454
455 return cosinePointingAngle;
456}
d6a49f20 457
458
459// **** The following functions need to be revised
460
074f017b 461void AliESDv0::GetPosCov(Double_t cov[6]) const {
462
463 for (Int_t i=0; i<6; ++i) cov[i] = fPosCov[i];
464
465}
466
d6a49f20 467Double_t AliESDv0::GetSigmaY(){
468 //
469 // return sigmay in y at vertex position using covariance matrix
470 //
471 const Double_t * cp = fParamP.GetCovariance();
472 const Double_t * cm = fParamN.GetCovariance();
473 Double_t sigmay = cp[0]+cm[0]+ cp[5]*(fParamP.GetX()-fRr)*(fParamP.GetX()-fRr)+ cm[5]*(fParamN.GetX()-fRr)*(fParamN.GetX()-fRr);
474 return (sigmay>0) ? TMath::Sqrt(sigmay):100;
475}
476
477Double_t AliESDv0::GetSigmaZ(){
478 //
479 // return sigmay in y at vertex position using covariance matrix
480 //
481 const Double_t * cp = fParamP.GetCovariance();
482 const Double_t * cm = fParamN.GetCovariance();
483 Double_t sigmaz = cp[2]+cm[2]+ cp[9]*(fParamP.GetX()-fRr)*(fParamP.GetX()-fRr)+ cm[9]*(fParamN.GetX()-fRr)*(fParamN.GetX()-fRr);
484 return (sigmaz>0) ? TMath::Sqrt(sigmaz):100;
485}
486
487Double_t AliESDv0::GetSigmaD0(){
488 //
489 // Sigma parameterization using covariance matrix
490 //
491 // sigma of distance between two tracks in vertex position
492 // sigma of DCA is proportianal to sigmaD0
493 // factor 2 difference is explained by the fact that the DCA is calculated at the position
494 // where the tracks as closest together ( not exact position of the vertex)
495 //
496 const Double_t * cp = fParamP.GetCovariance();
497 const Double_t * cm = fParamN.GetCovariance();
498 Double_t sigmaD0 = cp[0]+cm[0]+cp[2]+cm[2]+fgkParams.fPSigmaOffsetD0*fgkParams.fPSigmaOffsetD0;
499 sigmaD0 += ((fParamP.GetX()-fRr)*(fParamP.GetX()-fRr))*(cp[5]+cp[9]);
500 sigmaD0 += ((fParamN.GetX()-fRr)*(fParamN.GetX()-fRr))*(cm[5]+cm[9]);
501 return (sigmaD0>0)? TMath::Sqrt(sigmaD0):100;
502}
503
504
505Double_t AliESDv0::GetSigmaAP0(){
506 //
507 //Sigma parameterization using covariance matrices
508 //
b75d63a7 509 Double_t prec = TMath::Sqrt((fNmom[0]+fPmom[0])*(fNmom[0]+fPmom[0])
510 +(fNmom[1]+fPmom[1])*(fNmom[1]+fPmom[1])
511 +(fNmom[2]+fPmom[2])*(fNmom[2]+fPmom[2]));
512 Double_t normp = TMath::Sqrt(fPmom[0]*fPmom[0]+fPmom[1]*fPmom[1]+fPmom[2]*fPmom[2])/prec; // fraction of the momenta
513 Double_t normm = TMath::Sqrt(fNmom[0]*fNmom[0]+fNmom[1]*fNmom[1]+fNmom[2]*fNmom[2])/prec;
d6a49f20 514 const Double_t * cp = fParamP.GetCovariance();
515 const Double_t * cm = fParamN.GetCovariance();
516 Double_t sigmaAP0 = fgkParams.fPSigmaOffsetAP0*fgkParams.fPSigmaOffsetAP0; // minimal part
517 sigmaAP0 += (cp[5]+cp[9])*(normp*normp)+(cm[5]+cm[9])*(normm*normm); // angular resolution part
518 Double_t sigmaAP1 = GetSigmaD0()/(TMath::Abs(fRr)+0.01); // vertex position part
519 sigmaAP0 += 0.5*sigmaAP1*sigmaAP1;
520 return (sigmaAP0>0)? TMath::Sqrt(sigmaAP0):100;
521}
522
523Double_t AliESDv0::GetEffectiveSigmaD0(){
524 //
525 // minimax - effective Sigma parameterization
526 // p12 effective curvature and v0 radius postion used as parameters
527 //
528 Double_t p12 = TMath::Sqrt(fParamP.GetParameter()[4]*fParamP.GetParameter()[4]+
529 fParamN.GetParameter()[4]*fParamN.GetParameter()[4]);
530 Double_t sigmaED0= TMath::Max(TMath::Sqrt(fRr)-fgkParams.fPSigmaRminDE,0.0)*fgkParams.fPSigmaCoefDE*p12*p12;
531 sigmaED0*= sigmaED0;
532 sigmaED0*= sigmaED0;
533 sigmaED0 = TMath::Sqrt(sigmaED0+fgkParams.fPSigmaOffsetDE*fgkParams.fPSigmaOffsetDE);
534 return (sigmaED0<fgkParams.fPSigmaMaxDE) ? sigmaED0: fgkParams.fPSigmaMaxDE;
535}
536
537
538Double_t AliESDv0::GetEffectiveSigmaAP0(){
539 //
540 // effective Sigma parameterization of point angle resolution
541 //
542 Double_t p12 = TMath::Sqrt(fParamP.GetParameter()[4]*fParamP.GetParameter()[4]+
543 fParamN.GetParameter()[4]*fParamN.GetParameter()[4]);
544 Double_t sigmaAPE= fgkParams.fPSigmaBase0APE;
545 sigmaAPE+= fgkParams.fPSigmaR0APE/(fgkParams.fPSigmaR1APE+fRr);
546 sigmaAPE*= (fgkParams.fPSigmaP0APE+fgkParams.fPSigmaP1APE*p12);
547 sigmaAPE = TMath::Min(sigmaAPE,fgkParams.fPSigmaMaxAPE);
548 return sigmaAPE;
549}
550
551
552Double_t AliESDv0::GetMinimaxSigmaAP0(){
553 //
554 // calculate mini-max effective sigma of point angle resolution
555 //
556 //compv0->fTree->SetAlias("SigmaAP2","max(min((SigmaAP0+SigmaAPE0)*0.5,1.5*SigmaAPE0),0.5*SigmaAPE0+0.003)");
557 Double_t effectiveSigma = GetEffectiveSigmaAP0();
558 Double_t sigmaMMAP = 0.5*(GetSigmaAP0()+effectiveSigma);
559 sigmaMMAP = TMath::Min(sigmaMMAP, fgkParams.fPMaxFractionAP0*effectiveSigma);
560 sigmaMMAP = TMath::Max(sigmaMMAP, fgkParams.fPMinFractionAP0*effectiveSigma+fgkParams.fPMinAP0);
561 return sigmaMMAP;
562}
563Double_t AliESDv0::GetMinimaxSigmaD0(){
564 //
565 // calculate mini-max sigma of dca resolution
566 //
567 //compv0->fTree->SetAlias("SigmaD2","max(min((SigmaD0+SigmaDE0)*0.5,1.5*SigmaDE0),0.5*SigmaDE0)");
568 Double_t effectiveSigma = GetEffectiveSigmaD0();
569 Double_t sigmaMMD0 = 0.5*(GetSigmaD0()+effectiveSigma);
570 sigmaMMD0 = TMath::Min(sigmaMMD0, fgkParams.fPMaxFractionD0*effectiveSigma);
571 sigmaMMD0 = TMath::Max(sigmaMMD0, fgkParams.fPMinFractionD0*effectiveSigma+fgkParams.fPMinD0);
572 return sigmaMMD0;
573}
574
575
576Double_t AliESDv0::GetLikelihoodAP(Int_t mode0, Int_t mode1){
577 //
578 // get likelihood for point angle
579 //
580 Double_t sigmaAP = 0.007; //default sigma
581 switch (mode0){
582 case 0:
583 sigmaAP = GetSigmaAP0(); // mode 0 - covariance matrix estimates used
584 break;
585 case 1:
586 sigmaAP = GetEffectiveSigmaAP0(); // mode 1 - effective sigma used
587 break;
588 case 2:
589 sigmaAP = GetMinimaxSigmaAP0(); // mode 2 - minimax sigma
590 break;
591 }
592 Double_t apNorm = TMath::Min(TMath::ACos(fPointAngle)/sigmaAP,50.);
593 //normalized point angle, restricted - because of overflow problems in Exp
594 Double_t likelihood = 0;
595 switch(mode1){
596 case 0:
597 likelihood = TMath::Exp(-0.5*apNorm*apNorm);
598 // one component
599 break;
600 case 1:
601 likelihood = (TMath::Exp(-0.5*apNorm*apNorm)+0.5* TMath::Exp(-0.25*apNorm*apNorm))/1.5;
602 // two components
603 break;
604 case 2:
605 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;
606 // three components
607 break;
608 }
609 return likelihood;
610}
611
612Double_t AliESDv0::GetLikelihoodD(Int_t mode0, Int_t mode1){
613 //
614 // get likelihood for DCA
615 //
616 Double_t sigmaD = 0.03; //default sigma
617 switch (mode0){
618 case 0:
619 sigmaD = GetSigmaD0(); // mode 0 - covariance matrix estimates used
620 break;
621 case 1:
622 sigmaD = GetEffectiveSigmaD0(); // mode 1 - effective sigma used
623 break;
624 case 2:
625 sigmaD = GetMinimaxSigmaD0(); // mode 2 - minimax sigma
626 break;
627 }
b75d63a7 628
629 //Bo: Double_t dNorm = TMath::Min(fDist2/sigmaD,50.);
630 Double_t dNorm = TMath::Min(fDcaV0Daughters/sigmaD,50.);//Bo:
d6a49f20 631 //normalized point angle, restricted - because of overflow problems in Exp
632 Double_t likelihood = 0;
633 switch(mode1){
634 case 0:
635 likelihood = TMath::Exp(-2.*dNorm);
636 // one component
637 break;
638 case 1:
639 likelihood = (TMath::Exp(-2.*dNorm)+0.5* TMath::Exp(-dNorm))/1.5;
640 // two components
641 break;
642 case 2:
643 likelihood = (TMath::Exp(-2.*dNorm)+0.5* TMath::Exp(-dNorm)+0.25*TMath::Exp(-0.5*dNorm))/1.75;
644 // three components
645 break;
646 }
647 return likelihood;
648
649}
650
97135b34 651Double_t AliESDv0::GetLikelihoodC(Int_t mode0, Int_t /*mode1*/) const {
d6a49f20 652 //
653 // get likelihood for Causality
654 // !!! Causality variables defined in AliITStrackerMI !!!
655 // when more information was available
656 //
657 Double_t likelihood = 0.5;
658 Double_t minCausal = TMath::Min(fCausality[0],fCausality[1]);
659 Double_t maxCausal = TMath::Max(fCausality[0],fCausality[1]);
660 // minCausal = TMath::Max(minCausal,0.5*maxCausal);
661 //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");
662
663 switch(mode0){
664 case 0:
665 //normalization
666 likelihood = TMath::Power((1.05-2*(0.8-TMath::Exp(-maxCausal))),4.);
667 break;
668 case 1:
669 likelihood = TMath::Power(1.05-(2*(0.8-TMath::Exp(-maxCausal))+(2*(0.8-TMath::Exp(-minCausal))))*0.5,4.);
670 break;
671 }
672 return likelihood;
673
674}
675
676void AliESDv0::SetCausality(Float_t pb0, Float_t pb1, Float_t pa0, Float_t pa1)
677{
678 //
679 // set probabilities
680 //
681 fCausality[0] = pb0; // probability - track 0 exist before vertex
682 fCausality[1] = pb1; // probability - track 1 exist before vertex
683 fCausality[2] = pa0; // probability - track 0 exist close after vertex
684 fCausality[3] = pa1; // probability - track 1 exist close after vertex
685}
97135b34 686void AliESDv0::SetClusters(const Int_t *clp, const Int_t *clm)
d6a49f20 687{
688 //
689 // Set its clusters indexes
690 //
691 for (Int_t i=0;i<6;i++) fClusters[0][i] = clp[i];
692 for (Int_t i=0;i<6;i++) fClusters[1][i] = clm[i];
693}
9973a4bf 694
97135b34 695Double_t AliESDv0::GetEffMass(UInt_t p1, UInt_t p2) const{
9973a4bf 696 //
697 // calculate effective mass
698 //
89cb3a25 699 const Double_t kpmass[5] = {TDatabasePDG::Instance()->GetParticle(kElectron)->Mass(),
586d906c 700 TDatabasePDG::Instance()->GetParticle(kMuonMinus)->Mass(),
701 TDatabasePDG::Instance()->GetParticle(kPiPlus)->Mass(),
702 TDatabasePDG::Instance()->GetParticle(kKPlus)->Mass(),
703 TDatabasePDG::Instance()->GetParticle(kProton)->Mass()};
2278bbc8 704 /*
9973a4bf 705 if (p1>4) return -1;
706 if (p2>4) return -1;
89cb3a25 707 Double_t mass1 = kpmass[p1];
708 Double_t mass2 = kpmass[p2];
97135b34 709 const Double_t *m1 = fPmom;
710 const Double_t *m2 = fNmom;
9973a4bf 711 //
712 //if (fRP[p1]+fRM[p2]<fRP[p2]+fRM[p1]){
713 // m1 = fPM;
714 // m2 = fPP;
715 //}
716 //
89cb3a25 717 Double_t e1 = TMath::Sqrt(mass1*mass1+
9973a4bf 718 m1[0]*m1[0]+
719 m1[1]*m1[1]+
720 m1[2]*m1[2]);
89cb3a25 721 Double_t e2 = TMath::Sqrt(mass2*mass2+
9973a4bf 722 m2[0]*m2[0]+
723 m2[1]*m2[1]+
724 m2[2]*m2[2]);
89cb3a25 725 Double_t mass =
9973a4bf 726 (m2[0]+m1[0])*(m2[0]+m1[0])+
727 (m2[1]+m1[1])*(m2[1]+m1[1])+
728 (m2[2]+m1[2])*(m2[2]+m1[2]);
729
89cb3a25 730 mass = (e1+e2)*(e1+e2)-mass;
731 if (mass < 0.) mass = 0.;
732 return (TMath::Sqrt(mass));
2278bbc8 733 */
734 if(p1>4 || p2>4) return -1;
735 Double_t e12 = kpmass[p1]*kpmass[p1]+fPmom[0]*fPmom[0]+fPmom[1]*fPmom[1]+fPmom[2]*fPmom[2];
736 Double_t e22 = kpmass[p2]*kpmass[p2]+fNmom[0]*fNmom[0]+fNmom[1]*fNmom[1]+fNmom[2]*fNmom[2];
737 Double_t cmass = TMath::Sqrt(TMath::Max(kpmass[p1]*kpmass[p1]+kpmass[p2]*kpmass[p2]
738 +2.*(TMath::Sqrt(e12*e22)-fPmom[0]*fNmom[0]-fPmom[1]*fNmom[1]-fPmom[2]*fNmom[2]),0.));
739 return cmass;
740
9973a4bf 741}