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