]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/ESD/AliESDv0.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Updates in event mixing code for low-pt code
[u/mrichter/AliRoot.git] / STEER / ESD / AliESDv0.cxx
CommitLineData
e23730c7 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
4f679a16 16/* $Id$ */
17
e23730c7 18//-------------------------------------------------------------------------
19// Implementation of the ESD V0 vertex class
4f679a16 20// This class is part of the Event Data Summary
21// set of classes and contains information about
22// V0 kind vertexes generated by a neutral particle
e23730c7 23// Origin: Iouri Belikov, IReS, Strasbourg, Jouri.Belikov@cern.ch
d6a49f20 24// Modified by: Marian Ivanov, CERN, Marian.Ivanov@cern.ch
25// and Boris Hippolyte,IPHC, hippolyt@in2p3.fr
e23730c7 26//-------------------------------------------------------------------------
4f679a16 27
e23730c7 28#include <TMath.h>
90e48c0c 29#include <TDatabasePDG.h>
90e48c0c 30#include <TParticlePDG.h>
3e27a419 31#include <TVector3.h>
e23730c7 32
5f7789fc 33#include "AliLog.h"
e23730c7 34#include "AliESDv0.h"
97135b34 35#include "AliESDV0Params.h"
25f906db 36#include "AliKFParticle.h"
37#include "AliKFVertex.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
ed8f3bb5 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);
7935dc9e 168
169 Int_t index=fNidx;
170 fNidx=fPidx;
805508b1 171 fPidx=index;
ed8f3bb5 172 }
173
c7bafca9 174 for (Int_t i=0; i<6; i++) {
175 fPosCov[i]= 0.;
c7bafca9 176 }
177
178 //Trivial estimation of the vertex parameters
b75d63a7 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];
37919f0b 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;
c7bafca9 187
188
b75d63a7 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];
37919f0b 194 Double_t sx2=sn*sn*t2.GetSigmaY2()+ss, sy2=cs*cs*t2.GetSigmaY2()+ss;
c7bafca9 195
196 Double_t sz1=t1.GetSigmaZ2(), sz2=t2.GetSigmaZ2();
37919f0b 197 Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
198 Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
c7bafca9 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
b75d63a7 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;}
c028b974 210}
c7bafca9 211
3e27a419 212AliESDv0& AliESDv0::operator=(const AliESDv0 &v0)
213{
732a24fe 214 //--------------------------------------------------------------------
3e27a419 215 // The assignment operator
732a24fe 216 //--------------------------------------------------------------------
217
218 if(this==&v0)return *this;
646c9704 219 AliVParticle::operator=(v0);
732a24fe 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
3e27a419 265 // this overwrites the virtual TOBject::Copy()
732a24fe 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
c028b974 275AliESDv0::~AliESDv0(){
276 //--------------------------------------------------------------------
277 // Empty destructor
278 //--------------------------------------------------------------------
c7bafca9 279}
280
646c9704 281// Start with AliVParticle functions
282Double_t AliESDv0::E() const {
283 //--------------------------------------------------------------------
284 // This gives the energy assuming the ChangeMassHypothesis was called
285 //--------------------------------------------------------------------
3e27a419 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();
646c9704 302 return TMath::Sqrt(mass*mass+P()*P());
303}
c028b974 304
3e27a419 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
97135b34 335 Double_t lQlNeg = momNeg.Dot(momTot)/momTot.Mag();
ad07c38b 336 Double_t lQlPos = momPos.Dot(momTot)/momTot.Mag();
3e27a419 337
bc88fc12 338 //return 1.-2./(1.+lQlNeg/lQlPos);
339 return (lQlPos - lQlNeg)/(lQlPos + lQlNeg);
3e27a419 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
e23730c7 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 //--------------------------------------------------------------------
b75d63a7 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;
e23730c7 368
369 fPdgCode=code;
370
371 switch (code) {
372 case kLambda0:
b75d63a7 373 nmass=piMass; pmass=prMass; mass=l0Mass; ps=0.101; break;
e23730c7 374 case kLambda0Bar:
b75d63a7 375 pmass=piMass; nmass=prMass; mass=l0Mass; ps=0.101; break;
e23730c7 376 case kK0Short:
377 break;
378 default:
5f7789fc 379 AliError("invalide PDG code ! Assuming K0s...");
e23730c7 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
a0fd52b4 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
b75d63a7 439Float_t AliESDv0::GetD(Double_t x0, Double_t y0, Double_t z0) const {
e23730c7 440 //--------------------------------------------------------------------
a0fd52b4 441 // This function returns V0's impact parameter calculated in 3D
e23730c7 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}
c028b974 454
97135b34 455Float_t AliESDv0::GetV0CosineOfPointingAngle(Double_t refPointX, Double_t refPointY, Double_t refPointZ) const {
c028b974 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}
d6a49f20 476
477
478// **** The following functions need to be revised
479
074f017b 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
d6a49f20 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 //
b75d63a7 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;
d6a49f20 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 }
b75d63a7 647
648 //Bo: Double_t dNorm = TMath::Min(fDist2/sigmaD,50.);
649 Double_t dNorm = TMath::Min(fDcaV0Daughters/sigmaD,50.);//Bo:
d6a49f20 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
97135b34 670Double_t AliESDv0::GetLikelihoodC(Int_t mode0, Int_t /*mode1*/) const {
d6a49f20 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}
97135b34 705void AliESDv0::SetClusters(const Int_t *clp, const Int_t *clm)
d6a49f20 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}
9973a4bf 713
97135b34 714Double_t AliESDv0::GetEffMass(UInt_t p1, UInt_t p2) const{
9973a4bf 715 //
716 // calculate effective mass
717 //
89cb3a25 718 const Double_t kpmass[5] = {TDatabasePDG::Instance()->GetParticle(kElectron)->Mass(),
586d906c 719 TDatabasePDG::Instance()->GetParticle(kMuonMinus)->Mass(),
720 TDatabasePDG::Instance()->GetParticle(kPiPlus)->Mass(),
721 TDatabasePDG::Instance()->GetParticle(kKPlus)->Mass(),
722 TDatabasePDG::Instance()->GetParticle(kProton)->Mass()};
2278bbc8 723 /*
9973a4bf 724 if (p1>4) return -1;
725 if (p2>4) return -1;
89cb3a25 726 Double_t mass1 = kpmass[p1];
727 Double_t mass2 = kpmass[p2];
97135b34 728 const Double_t *m1 = fPmom;
729 const Double_t *m2 = fNmom;
9973a4bf 730 //
731 //if (fRP[p1]+fRM[p2]<fRP[p2]+fRM[p1]){
732 // m1 = fPM;
733 // m2 = fPP;
734 //}
735 //
89cb3a25 736 Double_t e1 = TMath::Sqrt(mass1*mass1+
9973a4bf 737 m1[0]*m1[0]+
738 m1[1]*m1[1]+
739 m1[2]*m1[2]);
89cb3a25 740 Double_t e2 = TMath::Sqrt(mass2*mass2+
9973a4bf 741 m2[0]*m2[0]+
742 m2[1]*m2[1]+
743 m2[2]*m2[2]);
89cb3a25 744 Double_t mass =
9973a4bf 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
89cb3a25 749 mass = (e1+e2)*(e1+e2)-mass;
750 if (mass < 0.) mass = 0.;
751 return (TMath::Sqrt(mass));
2278bbc8 752 */
753 if(p1>4 || p2>4) return -1;
25f906db 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];
2278bbc8 765 Double_t cmass = TMath::Sqrt(TMath::Max(kpmass[p1]*kpmass[p1]+kpmass[p2]*kpmass[p2]
25f906db 766 +2.*(TMath::Sqrt(e12*e22)-pmom[0]*nmom[0]-pmom[1]*nmom[1]-pmom[2]*nmom[2]),0.));
2278bbc8 767 return cmass;
768
9973a4bf 769}
25f906db 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) );
339fbe23 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();
25f906db 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) );
339fbe23 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();
25f906db 831 return 0;
832}