]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/AliAODRecoDecay.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
bug fixed
[u/mrichter/AliRoot.git] / STEER / AliAODRecoDecay.cxx
CommitLineData
7de7497b 1/**************************************************************************
2 * Copyright(c) 1998-2006, 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/////////////////////////////////////////////////////////////
17//
18// Base class for AOD reconstructed decay
19//
20// Author: A.Dainese, andrea.dainese@lnl.infn.it
21/////////////////////////////////////////////////////////////
22
23#include <TDatabasePDG.h>
24#include <TVector3.h>
e045ffda 25#include "AliLog.h"
26#include "AliVTrack.h"
7de7497b 27#include "AliAODRecoDecay.h"
28
29ClassImp(AliAODRecoDecay)
30
31//--------------------------------------------------------------------------
32AliAODRecoDecay::AliAODRecoDecay() :
e045ffda 33 AliVTrack(),
7de7497b 34 fSecondaryVtx(0x0),
a11de4a0 35 fOwnSecondaryVtx(0x0),
7de7497b 36 fCharge(0),
37 fNProngs(0), fNDCA(0), fNPID(0),
38 fPx(0x0), fPy(0x0), fPz(0x0),
39 fd0(0x0),
40 fDCA(0x0),
41 fPID(0x0),
42 fEventNumber(-1),fRunNumber(-1)
43{
44 //
45 // Default Constructor
46 //
7de7497b 47}
48//--------------------------------------------------------------------------
49AliAODRecoDecay::AliAODRecoDecay(AliAODVertex *vtx2,Int_t nprongs,
50 Short_t charge,
51 Double_t *px,Double_t *py,Double_t *pz,
52 Double_t *d0) :
e045ffda 53 AliVTrack(),
7de7497b 54 fSecondaryVtx(vtx2),
a11de4a0 55 fOwnSecondaryVtx(0x0),
7de7497b 56 fCharge(charge),
57 fNProngs(nprongs), fNDCA(0), fNPID(0),
58 fPx(0x0), fPy(0x0), fPz(0x0),
59 fd0(0x0),
60 fDCA(0x0),
61 fPID(0x0),
62 fEventNumber(-1),fRunNumber(-1)
63{
64 //
65 // Constructor with AliAODVertex for decay vertex
66 //
67
68 fPx = new Double_t[GetNProngs()];
69 fPy = new Double_t[GetNProngs()];
70 fPz = new Double_t[GetNProngs()];
71 fd0 = new Double_t[GetNProngs()];
72 for(Int_t i=0; i<GetNProngs(); i++) {
73 fPx[i] = px[i];
74 fPy[i] = py[i];
75 fPz[i] = pz[i];
76 fd0[i] = d0[i];
77 }
78}
79//--------------------------------------------------------------------------
80AliAODRecoDecay::AliAODRecoDecay(AliAODVertex *vtx2,Int_t nprongs,
81 Short_t charge,
82 Double_t *d0) :
e045ffda 83 AliVTrack(),
7de7497b 84 fSecondaryVtx(vtx2),
a11de4a0 85 fOwnSecondaryVtx(0x0),
7de7497b 86 fCharge(charge),
87 fNProngs(nprongs), fNDCA(0), fNPID(0),
88 fPx(0x0), fPy(0x0), fPz(0x0),
89 fd0(0x0),
90 fDCA(0x0),
91 fPID(0x0),
92 fEventNumber(-1),fRunNumber(-1)
93{
94 //
95 // Constructor with AliAODVertex for decay vertex and without prongs momenta
96 //
97
98 fd0 = new Double_t[GetNProngs()];
99 for(Int_t i=0; i<GetNProngs(); i++) fd0[i] = d0[i];
100}
101//--------------------------------------------------------------------------
102AliAODRecoDecay::AliAODRecoDecay(const AliAODRecoDecay &source) :
e045ffda 103 AliVTrack(source),
7de7497b 104 fSecondaryVtx(source.fSecondaryVtx),
a11de4a0 105 fOwnSecondaryVtx(source.fOwnSecondaryVtx),
7de7497b 106 fCharge(source.fCharge),
107 fNProngs(source.fNProngs), fNDCA(source.fNDCA), fNPID(source.fNPID),
108 fPx(0x0), fPy(0x0), fPz(0x0),
109 fd0(0x0),
110 fDCA(0x0),
111 fPID(0x0),
112 fEventNumber(source.fEventNumber),fRunNumber(source.fRunNumber)
113{
114 //
115 // Copy constructor
116 //
117 if(source.GetNProngs()>0) {
58b0186f 118 fd0 = new Double32_t[GetNProngs()];
119 memcpy(fd0,source.fd0,GetNProngs()*sizeof(Double32_t));
7de7497b 120 if(source.fPx) {
58b0186f 121 fPx = new Double32_t[GetNProngs()];
122 fPy = new Double32_t[GetNProngs()];
123 fPz = new Double32_t[GetNProngs()];
124 memcpy(fPx,source.fPx,GetNProngs()*sizeof(Double32_t));
125 memcpy(fPy,source.fPy,GetNProngs()*sizeof(Double32_t));
126 memcpy(fPz,source.fPz,GetNProngs()*sizeof(Double32_t));
7de7497b 127 }
128 if(source.fPID) {
58b0186f 129 fPID = new Double32_t[fNPID];
130 memcpy(fPID,source.fPID,fNPID*sizeof(Double32_t));
7de7497b 131 }
132 if(source.fDCA) {
58b0186f 133 fDCA = new Double32_t[fNDCA];
134 memcpy(fDCA,source.fDCA,fNDCA*sizeof(Double32_t));
7de7497b 135 }
136 }
137}
138//--------------------------------------------------------------------------
139AliAODRecoDecay &AliAODRecoDecay::operator=(const AliAODRecoDecay &source)
140{
141 //
142 // assignment operator
143 //
144 if(&source == this) return *this;
145 fSecondaryVtx = source.fSecondaryVtx;
a11de4a0 146 fOwnSecondaryVtx = source.fOwnSecondaryVtx;
7de7497b 147 fCharge = source.fCharge;
148 fNProngs = source.fNProngs;
149 fNDCA = source.fNDCA;
150 fNPID = source.fNPID;
151 fEventNumber = source.fEventNumber;
152 fRunNumber = source.fRunNumber;
153 if(source.GetNProngs()>0) {
341952b8 154 if(fd0)delete [] fd0;
58b0186f 155 fd0 = new Double32_t[GetNProngs()];
156 memcpy(fd0,source.fd0,GetNProngs()*sizeof(Double32_t));
7de7497b 157 if(source.fPx) {
341952b8 158 if(fPx) delete [] fPx;
58b0186f 159 fPx = new Double32_t[GetNProngs()];
341952b8 160 if(fPy) delete [] fPy;
58b0186f 161 fPy = new Double32_t[GetNProngs()];
341952b8 162 if(fPz) delete [] fPz;
58b0186f 163 fPz = new Double32_t[GetNProngs()];
164 memcpy(fPx,source.fPx,GetNProngs()*sizeof(Double32_t));
165 memcpy(fPy,source.fPy,GetNProngs()*sizeof(Double32_t));
166 memcpy(fPz,source.fPz,GetNProngs()*sizeof(Double32_t));
7de7497b 167 }
168 if(source.fPID) {
341952b8 169 if(fPID) delete [] fPID;
58b0186f 170 fPID = new Double32_t[fNPID];
171 memcpy(fPID,source.fPID,fNPID*sizeof(Double32_t));
7de7497b 172 }
173 if(source.fDCA) {
341952b8 174 if(fDCA) delete [] fDCA;
58b0186f 175 fDCA = new Double32_t[fNDCA];
176 memcpy(fDCA,source.fDCA,fNDCA*sizeof(Double32_t));
7de7497b 177 }
178 }
179 return *this;
180}
181//--------------------------------------------------------------------------
182AliAODRecoDecay::~AliAODRecoDecay() {
183 //
184 // Default Destructor
185 //
a11de4a0 186 if(fPx) { delete [] fPx; fPx=NULL; }
187 if(fPy) { delete [] fPy; fPy=NULL; }
188 if(fPz) { delete [] fPz; fPz=NULL; }
189 if(fd0) { delete [] fd0; fd0=NULL; }
190 if(fPID) { delete [] fPID; fPID=NULL; }
191 if(fDCA) { delete [] fDCA; fDCA=NULL; }
192 if(fOwnSecondaryVtx) { delete fOwnSecondaryVtx; fOwnSecondaryVtx=NULL; }
7de7497b 193}
194//----------------------------------------------------------------------------
195Double_t AliAODRecoDecay::Alpha() const
196{
197 //
198 // Armenteros-Podolanski alpha for 2-prong decays
199 //
200 if(GetNProngs()!=2) {
201 printf("Can be called only for 2-prong decays");
202 return (Double_t)-99999.;
203 }
204 return 1.-2./(1.+QlProng(0)/QlProng(1));
205}
206//----------------------------------------------------------------------------
207Double_t AliAODRecoDecay::DecayLength(Double_t point[3]) const
208{
209 //
210 // Decay length assuming it is produced at "point" [cm]
211 //
212 return TMath::Sqrt((point[0]-GetSecVtxX())
213 *(point[0]-GetSecVtxX())
214 +(point[1]-GetSecVtxY())
215 *(point[1]-GetSecVtxY())
216 +(point[2]-GetSecVtxZ())
217 *(point[2]-GetSecVtxZ()));
218}
219//----------------------------------------------------------------------------
220Double_t AliAODRecoDecay::DecayLengthXY(Double_t point[3]) const
221{
222 //
223 // Decay length in XY assuming it is produced at "point" [cm]
224 //
225 return TMath::Sqrt((point[0]-GetSecVtxX())
226 *(point[0]-GetSecVtxX())
227 +(point[1]-GetSecVtxY())
228 *(point[1]-GetSecVtxY()));
229}
230//----------------------------------------------------------------------------
231Double_t AliAODRecoDecay::CosPointingAngle(Double_t point[3]) const
232{
233 //
234 // Cosine of pointing angle in space assuming it is produced at "point"
235 //
236 TVector3 mom(Px(),Py(),Pz());
237 TVector3 fline(GetSecVtxX()-point[0],
238 GetSecVtxY()-point[1],
239 GetSecVtxZ()-point[2]);
240
241 Double_t pta = mom.Angle(fline);
242
243 return TMath::Cos(pta);
244}
245//----------------------------------------------------------------------------
246Double_t AliAODRecoDecay::CosPointingAngleXY(Double_t point[3]) const
247{
248 //
249 // Cosine of pointing angle in transverse plane assuming it is produced
250 // at "point"
251 //
252 TVector3 momXY(Px(),Py(),0.);
253 TVector3 flineXY(GetSecVtxX()-point[0],
254 GetSecVtxY()-point[1],
255 0.);
256
257 Double_t ptaXY = momXY.Angle(flineXY);
258
259 return TMath::Cos(ptaXY);
260}
261//----------------------------------------------------------------------------
262Double_t AliAODRecoDecay::CosThetaStar(Int_t ip,UInt_t pdgvtx,UInt_t pdgprong0,UInt_t pdgprong1) const
263{
264 //
265 // Only for 2-prong decays:
266 // Cosine of decay angle (theta*) in the rest frame of the mother particle
267 // for prong ip (0 or 1) with mass hypotheses pdgvtx for mother particle,
268 // pdgprong0 for prong 0 and pdgprong1 for prong1
269 //
270 if(GetNProngs()!=2) {
271 printf("Can be called only for 2-prong decays");
272 return (Double_t)-99999.;
273 }
274 Double_t massvtx = TDatabasePDG::Instance()->GetParticle(pdgvtx)->Mass();
275 Double_t massp[2];
276 massp[0] = TDatabasePDG::Instance()->GetParticle(pdgprong0)->Mass();
277 massp[1] = TDatabasePDG::Instance()->GetParticle(pdgprong1)->Mass();
278
279 Double_t pStar = TMath::Sqrt(TMath::Power(massvtx*massvtx-massp[0]*massp[0]-massp[1]*massp[1],2.)-4.*massp[0]*massp[0]*massp[1]*massp[1])/(2.*massvtx);
280
281 Double_t beta = P()/E(pdgvtx);
282 Double_t gamma = E(pdgvtx)/massvtx;
283
284 Double_t cts = (QlProng(ip)/gamma-beta*TMath::Sqrt(pStar*pStar+massp[ip]*massp[ip]))/pStar;
285
286 return cts;
287}
288//---------------------------------------------------------------------------
289Double_t AliAODRecoDecay::Ct(UInt_t pdg,Double_t point[3]) const
290{
291 //
292 // Decay time * c assuming it is produced at "point" [cm]
293 //
294 Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
295 return DecayLength(point)*mass/P();
296}
297//---------------------------------------------------------------------------
298Double_t AliAODRecoDecay::E(UInt_t pdg) const
299{
300 //
301 // Energy
302 //
303 Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
304 return TMath::Sqrt(mass*mass+P()*P());
305}
306//---------------------------------------------------------------------------
307Double_t AliAODRecoDecay::EProng(Int_t ip,UInt_t pdg) const
308{
309 //
310 // Energy of ip-th prong
311 //
312 Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
313 return TMath::Sqrt(mass*mass+PProng(ip)*PProng(ip));
314}
e045ffda 315//--------------------------------------------------------------------------
316Bool_t AliAODRecoDecay::GetCovarianceXYZPxPyPz(Double_t cv[21]) const {
317 //
318 // This function returns the global covariance matrix of the track params
319 //
320 // Cov(x,x) ... : cv[0]
321 // Cov(y,x) ... : cv[1] cv[2]
322 // Cov(z,x) ... : cv[3] cv[4] cv[5]
323 // Cov(px,x)... : cv[6] cv[7] cv[8] cv[9]
324 // Cov(py,x)... : cv[10] cv[11] cv[12] cv[13] cv[14]
325 // Cov(pz,x)... : cv[15] cv[16] cv[17] cv[18] cv[19] cv[20]
7de7497b 326 //
e045ffda 327 // For XYZ we take the cov of the vertex, for PxPyPz we take the
328 // sum of the covs of PxPyPz from the daughters, for the moment
329 // we set the cov between position and momentum as the sum of
330 // the same cov from the daughters.
7de7497b 331 //
e045ffda 332
333 Int_t j;
334 for(j=0;j<21;j++) cv[j]=0.;
335
336 if(!GetNDaughters()) {
337 AliError("No daughters available");
338 return kFALSE;
339 }
340
341 Double_t v[6];
342 AliAODVertex *secv=GetSecondaryVtx();
343 if(!secv) {
344 AliError("Vertex covariance matrix not available");
345 return kFALSE;
346 }
347 if(!secv->GetCovMatrix(v)) {
348 AliError("Vertex covariance matrix not available");
349 return kFALSE;
350 }
351
352 Double_t p[21]; for(j=0;j<21;j++) p[j]=0.;
353 Bool_t error=kFALSE;
354 for(Int_t i=1; i<GetNDaughters(); i++) {
355 AliVTrack *daugh = (AliVTrack*)GetDaughter(i);
356 Double_t dcov[21];
357 if(!daugh->GetCovarianceXYZPxPyPz(dcov)) error=kTRUE;
358 for(j=0;j<21;j++) p[j] += dcov[j];
359 }
360 if(error) {
361 AliError("No covariance for at least one daughter")
362 return kFALSE;
363 }
364
365 for(j=0; j<21; j++) {
366 if(j<6) {
367 cv[j] = v[j];
368 } else {
369 cv[j] = p[j];
370 }
371 }
372
373 return kTRUE;
374}
7de7497b 375//----------------------------------------------------------------------------
e045ffda 376UChar_t AliAODRecoDecay::GetITSClusterMap() const {
377 //
378 // We take the logical AND of the daughters cluster maps
379 // (only if all daughters have the bit for given layer, we set the bit)
380 //
381 UChar_t map=0;
382
383 if(!GetNDaughters()) {
384 AliError("No daughters available");
385 return map;
386 }
387
388 for(Int_t l=0; l<12; l++) { // loop on ITS layers (from here we cannot know how many they are; let's put 12 to be conservative)
389 Int_t bit = 1;
390 for(Int_t i=0; i<GetNDaughters(); i++) {
391 AliVTrack *daugh = (AliVTrack*)GetDaughter(i);
392 if(!TESTBIT(daugh->GetITSClusterMap(),l)) bit=0;
393 }
394 if(bit) SETBIT(map,l);
395 }
396
397 return map;
398}
399//--------------------------------------------------------------------------
400ULong_t AliAODRecoDecay::GetStatus() const {
401 //
402 // Same as for ITSClusterMap
403 //
404 ULong_t status=0;
405
406 if(!GetNDaughters()) {
407 AliError("No daughters available");
408 return status;
409 }
410
411 AliVTrack *daugh0 = (AliVTrack*)GetDaughter(0);
412 status = status&(daugh0->GetStatus());
413
414 for(Int_t i=1; i<GetNDaughters(); i++) {
415 AliVTrack *daugh = (AliVTrack*)GetDaughter(i);
416 status = status&(daugh->GetStatus());
417 }
418
419 return status;
420}
421//--------------------------------------------------------------------------
7de7497b 422Double_t AliAODRecoDecay::ImpParXY(Double_t point[3]) const
423{
424 //
425 // Impact parameter in the bending plane of the particle
426 // w.r.t. to "point"
427 //
428 Double_t k = -(GetSecVtxX()-point[0])*Px()-(GetSecVtxY()-point[1])*Py();
429 k /= Pt()*Pt();
430 Double_t dx = GetSecVtxX()-point[0]+k*Px();
431 Double_t dy = GetSecVtxY()-point[1]+k*Py();
432 Double_t absImpPar = TMath::Sqrt(dx*dx+dy*dy);
433 TVector3 mom(Px(),Py(),Pz());
434 TVector3 fline(GetSecVtxX()-point[0],
435 GetSecVtxY()-point[1],
436 GetSecVtxZ()-point[2]);
437 TVector3 cross = mom.Cross(fline);
438 return (cross.Z()>0. ? absImpPar : -absImpPar);
439}
440//----------------------------------------------------------------------------
441Double_t AliAODRecoDecay::InvMass(Int_t npdg,UInt_t *pdg) const
442{
443 //
444 // Invariant mass for prongs mass hypotheses in pdg array
445 //
446 if(GetNProngs()!=npdg) {
447 printf("npdg != GetNProngs()");
448 return (Double_t)-99999.;
449 }
450 Double_t energysum = 0.;
451
452 for(Int_t i=0; i<GetNProngs(); i++) {
453 energysum += EProng(i,pdg[i]);
454 }
455
456 Double_t mass = TMath::Sqrt(energysum*energysum-P()*P());
457
458 return mass;
459}
460//----------------------------------------------------------------------------
461Double_t AliAODRecoDecay::InvMass2Prongs(Int_t ip1,Int_t ip2,
462 UInt_t pdg1,UInt_t pdg2) const
463{
464 //
465 // 2-prong(ip1,ip2) invariant mass for prongs mass hypotheses in pdg1,2
466 //
467 Double_t energysum = EProng(ip1,pdg1) + EProng(ip2,pdg2);
468 Double_t psum2 = (PxProng(ip1)+PxProng(ip2))*(PxProng(ip1)+PxProng(ip2))
469 +(PyProng(ip1)+PyProng(ip2))*(PyProng(ip1)+PyProng(ip2))
470 +(PzProng(ip1)+PzProng(ip2))*(PzProng(ip1)+PzProng(ip2));
471 Double_t mass = TMath::Sqrt(energysum*energysum-psum2);
472
473 return mass;
474}
475//---------------------------------------------------------------------------
476void AliAODRecoDecay::Print(Option_t* /*option*/) const
477{
478 //
479 // Print some information
480 //
481 printf("AliAODRecoDecay with %d prongs\n",GetNProngs());
482 printf("Secondary Vertex: (%f, %f, %f)\n",GetSecVtxX(),GetSecVtxY(),GetSecVtxZ());
483
484 return;
485}
486//----------------------------------------------------------------------------
487Double_t AliAODRecoDecay::ProngsRelAngle(Int_t ip1,Int_t ip2) const
488{
489 //
490 // Relative angle between two prongs
491 //
492 TVector3 momA(PxProng(ip1),PyProng(ip1),PzProng(ip1));
493 TVector3 momB(PxProng(ip2),PyProng(ip2),PzProng(ip2));
494
495 Double_t angle = momA.Angle(momB);
496
497 return angle;
498}
499//----------------------------------------------------------------------------
500Double_t AliAODRecoDecay::QlProng(Int_t ip) const
501{
502 //
503 // Longitudinal momentum of prong w.r.t. to total momentum
504 //
505 TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
506 TVector3 momTot(Px(),Py(),Pz());
507
508 return mom.Dot(momTot)/momTot.Mag();
509}
510//----------------------------------------------------------------------------
511Double_t AliAODRecoDecay::QtProng(Int_t ip) const
512{
513 //
514 // Transverse momentum of prong w.r.t. to total momentum
515 //
516 TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
517 TVector3 momTot(Px(),Py(),Pz());
518
519 return mom.Perp(momTot);
520}
521//----------------------------------------------------------------------------
522Double_t AliAODRecoDecay::QlProngFlightLine(Int_t ip,Double_t point[3]) const
523{
524 //
525 // Longitudinal momentum of prong w.r.t. to flight line between "point"
526 // and fSecondaryVtx
527 //
528 TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
529 TVector3 fline(GetSecVtxX()-point[0],
530 GetSecVtxY()-point[1],
531 GetSecVtxZ()-point[2]);
532
533 return mom.Dot(fline)/fline.Mag();
534}
535//----------------------------------------------------------------------------
536Double_t AliAODRecoDecay::QtProngFlightLine(Int_t ip,Double_t point[3]) const
537{
538 //
539 // Transverse momentum of prong w.r.t. to flight line between "point" and
540 // fSecondaryVtx
541 //
542 TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
543 TVector3 fline(GetSecVtxX()-point[0],
544 GetSecVtxY()-point[1],
545 GetSecVtxZ()-point[2]);
546
547 return mom.Perp(fline);
548}
549//--------------------------------------------------------------------------