]> 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
A little fix to the previous commit.
[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
37ClassImp(AliESDv0)
38
39const AliESDV0Params AliESDv0::fgkParams;
40
41AliESDv0::AliESDv0() :
42 AliVParticle(),
43 fParamN(),
44 fParamP(),
45 fEffMass(TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass()),
46 fDcaV0Daughters(0),
47 fChi2V0(0.),
48 fRr(0),
49 fDistSigma(0),
50 fChi2Before(0),
51 fChi2After(0),
52 fPointAngleFi(0),
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)
62{
63 //--------------------------------------------------------------------
64 // Default constructor (K0s)
65 //--------------------------------------------------------------------
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.;
75 }
76
77 for (Int_t i=0;i<6;i++){fClusters[0][i]=0; fClusters[1][i]=0;}
78 fNormDCAPrim[0]=fNormDCAPrim[1]=0;
79 for (Int_t i=0;i<3;i++){fAngle[i]=0;}
80 for (Int_t i=0;i<4;i++){fCausality[i]=0;}
81}
82
83AliESDv0::AliESDv0(const AliESDv0& v0) :
84 AliVParticle(v0),
85 fParamN(v0.fParamN),
86 fParamP(v0.fParamP),
87 fEffMass(v0.fEffMass),
88 fDcaV0Daughters(v0.fDcaV0Daughters),
89 fChi2V0(v0.fChi2V0),
90 fRr(v0.fRr),
91 fDistSigma(v0.fDistSigma),
92 fChi2Before(v0.fChi2Before),
93 fChi2After(v0.fChi2After),
94 fPointAngleFi(v0.fPointAngleFi),
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)
104{
105 //--------------------------------------------------------------------
106 // The copy constructor
107 //--------------------------------------------------------------------
108
109 for (int i=0; i<3; i++) {
110 fPos[i] = v0.fPos[i];
111 fNmom[i] = v0.fNmom[i];
112 fPmom[i] = v0.fPmom[i];
113 }
114 for (int i=0; i<6; i++) {
115 fPosCov[i] = v0.fPosCov[i];
116 }
117
118 for (Int_t i=0; i<2; i++) {
119 fNormDCAPrim[i]=v0.fNormDCAPrim[i];
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++){
126 fAngle[i]=v0.fAngle[i];
127 }
128 for (Int_t i=0;i<4;i++){fCausality[i]=v0.fCausality[i];}
129}
130
131
132AliESDv0::AliESDv0(const AliExternalTrackParam &t1, Int_t i1,
133 const AliExternalTrackParam &t2, Int_t i2) :
134 AliVParticle(),
135 fParamN(t1),
136 fParamP(t2),
137 fEffMass(TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass()),
138 fDcaV0Daughters(0),
139 fChi2V0(0.),
140 fRr(0),
141 fDistSigma(0),
142 fChi2Before(0),
143 fChi2After(0),
144 fPointAngleFi(0),
145 fPointAngleTh(0),
146 fPointAngle(0),
147 fPdgCode(kK0Short),
148 fNidx(i1),
149 fPidx(i2),
150 fStatus(0),
151 fNBefore(0),
152 fNAfter(0),
153 fOnFlyStatus(kFALSE)
154{
155 //--------------------------------------------------------------------
156 // Main constructor (K0s)
157 //--------------------------------------------------------------------
158
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 Int_t index=fNidx;
168 fNidx=fPidx;
169 fPidx=index;
170 }
171
172 for (Int_t i=0; i<6; i++) {
173 fPosCov[i]= 0.;
174 }
175
176 //Trivial estimation of the vertex parameters
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];
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;
185
186
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];
192 Double_t sx2=sn*sn*t2.GetSigmaY2()+ss, sy2=cs*cs*t2.GetSigmaY2()+ss;
193
194 Double_t sz1=t1.GetSigmaZ2(), sz2=t2.GetSigmaZ2();
195 Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
196 Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
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
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;}
208}
209
210AliESDv0& AliESDv0::operator=(const AliESDv0 &v0)
211{
212 //--------------------------------------------------------------------
213 // The assignment operator
214 //--------------------------------------------------------------------
215
216 if(this==&v0)return *this;
217 AliVParticle::operator=(v0);
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
263 // this overwrites the virtual TOBject::Copy()
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
273AliESDv0::~AliESDv0(){
274 //--------------------------------------------------------------------
275 // Empty destructor
276 //--------------------------------------------------------------------
277}
278
279// Start with AliVParticle functions
280Double_t AliESDv0::E() const {
281 //--------------------------------------------------------------------
282 // This gives the energy assuming the ChangeMassHypothesis was called
283 //--------------------------------------------------------------------
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();
300 return TMath::Sqrt(mass*mass+P()*P());
301}
302
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
333 Double_t lQlNeg = momNeg.Dot(momTot)/momTot.Mag();
334 Double_t lQlPos = momPos.Dot(momTot)/momTot.Mag();
335
336 //return 1.-2./(1.+lQlNeg/lQlPos);
337 return (lQlPos - lQlNeg)/(lQlPos + lQlNeg);
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
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 //--------------------------------------------------------------------
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;
366
367 fPdgCode=code;
368
369 switch (code) {
370 case kLambda0:
371 nmass=piMass; pmass=prMass; mass=l0Mass; ps=0.101; break;
372 case kLambda0Bar:
373 pmass=piMass; nmass=prMass; mass=l0Mass; ps=0.101; break;
374 case kK0Short:
375 break;
376 default:
377 AliError("invalide PDG code ! Assuming K0s...");
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
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
437Float_t AliESDv0::GetD(Double_t x0, Double_t y0, Double_t z0) const {
438 //--------------------------------------------------------------------
439 // This function returns V0's impact parameter calculated in 3D
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}
452
453Float_t AliESDv0::GetV0CosineOfPointingAngle(Double_t refPointX, Double_t refPointY, Double_t refPointZ) const {
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}
474
475
476// **** The following functions need to be revised
477
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
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 //
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;
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 }
645
646 //Bo: Double_t dNorm = TMath::Min(fDist2/sigmaD,50.);
647 Double_t dNorm = TMath::Min(fDcaV0Daughters/sigmaD,50.);//Bo:
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
668Double_t AliESDv0::GetLikelihoodC(Int_t mode0, Int_t /*mode1*/) const {
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}
703void AliESDv0::SetClusters(const Int_t *clp, const Int_t *clm)
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}
711
712Double_t AliESDv0::GetEffMass(UInt_t p1, UInt_t p2) const{
713 //
714 // calculate effective mass
715 //
716 const Double_t kpmass[5] = {TDatabasePDG::Instance()->GetParticle(kElectron)->Mass(),
717 TDatabasePDG::Instance()->GetParticle(kMuonMinus)->Mass(),
718 TDatabasePDG::Instance()->GetParticle(kPiPlus)->Mass(),
719 TDatabasePDG::Instance()->GetParticle(kKPlus)->Mass(),
720 TDatabasePDG::Instance()->GetParticle(kProton)->Mass()};
721 /*
722 if (p1>4) return -1;
723 if (p2>4) return -1;
724 Double_t mass1 = kpmass[p1];
725 Double_t mass2 = kpmass[p2];
726 const Double_t *m1 = fPmom;
727 const Double_t *m2 = fNmom;
728 //
729 //if (fRP[p1]+fRM[p2]<fRP[p2]+fRM[p1]){
730 // m1 = fPM;
731 // m2 = fPP;
732 //}
733 //
734 Double_t e1 = TMath::Sqrt(mass1*mass1+
735 m1[0]*m1[0]+
736 m1[1]*m1[1]+
737 m1[2]*m1[2]);
738 Double_t e2 = TMath::Sqrt(mass2*mass2+
739 m2[0]*m2[0]+
740 m2[1]*m2[1]+
741 m2[2]*m2[2]);
742 Double_t mass =
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
747 mass = (e1+e2)*(e1+e2)-mass;
748 if (mass < 0.) mass = 0.;
749 return (TMath::Sqrt(mass));
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
758}