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