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