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