1 /**************************************************************************
2 * Copyright(c) 1998-2006, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 /////////////////////////////////////////////////////////////
18 // Base class for AOD reconstructed decay
20 // Author: A.Dainese, andrea.dainese@lnl.infn.it
21 /////////////////////////////////////////////////////////////
23 #include <TDatabasePDG.h>
26 #include "AliVTrack.h"
27 #include "AliAODRecoDecay.h"
29 ClassImp(AliAODRecoDecay)
31 //--------------------------------------------------------------------------
32 AliAODRecoDecay::AliAODRecoDecay() :
35 fOwnSecondaryVtx(0x0),
37 fNProngs(0), fNDCA(0), fNPID(0),
38 fPx(0x0), fPy(0x0), fPz(0x0),
42 fEventNumber(-1),fRunNumber(-1)
45 // Default Constructor
48 //--------------------------------------------------------------------------
49 AliAODRecoDecay::AliAODRecoDecay(AliAODVertex *vtx2,Int_t nprongs,
51 Double_t *px,Double_t *py,Double_t *pz,
55 fOwnSecondaryVtx(0x0),
57 fNProngs(nprongs), fNDCA(0), fNPID(0),
58 fPx(0x0), fPy(0x0), fPz(0x0),
62 fEventNumber(-1),fRunNumber(-1)
65 // Constructor with AliAODVertex for decay vertex
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++) {
79 //--------------------------------------------------------------------------
80 AliAODRecoDecay::AliAODRecoDecay(AliAODVertex *vtx2,Int_t nprongs,
85 fOwnSecondaryVtx(0x0),
87 fNProngs(nprongs), fNDCA(0), fNPID(0),
88 fPx(0x0), fPy(0x0), fPz(0x0),
92 fEventNumber(-1),fRunNumber(-1)
95 // Constructor with AliAODVertex for decay vertex and without prongs momenta
98 fd0 = new Double_t[GetNProngs()];
99 for(Int_t i=0; i<GetNProngs(); i++) fd0[i] = d0[i];
101 //--------------------------------------------------------------------------
102 AliAODRecoDecay::AliAODRecoDecay(const AliAODRecoDecay &source) :
104 fSecondaryVtx(source.fSecondaryVtx),
105 fOwnSecondaryVtx(source.fOwnSecondaryVtx),
106 fCharge(source.fCharge),
107 fNProngs(source.fNProngs), fNDCA(source.fNDCA), fNPID(source.fNPID),
108 fPx(0x0), fPy(0x0), fPz(0x0),
112 fEventNumber(source.fEventNumber),fRunNumber(source.fRunNumber)
117 if(source.GetNProngs()>0) {
118 fd0 = new Double32_t[GetNProngs()];
119 memcpy(fd0,source.fd0,GetNProngs()*sizeof(Double32_t));
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));
129 fPID = new Double32_t[fNPID];
130 memcpy(fPID,source.fPID,fNPID*sizeof(Double32_t));
133 fDCA = new Double32_t[fNDCA];
134 memcpy(fDCA,source.fDCA,fNDCA*sizeof(Double32_t));
138 //--------------------------------------------------------------------------
139 AliAODRecoDecay &AliAODRecoDecay::operator=(const AliAODRecoDecay &source)
142 // assignment operator
144 if(&source == this) return *this;
145 fSecondaryVtx = source.fSecondaryVtx;
146 fOwnSecondaryVtx = source.fOwnSecondaryVtx;
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) {
154 if(fd0)delete [] fd0;
155 fd0 = new Double32_t[GetNProngs()];
156 memcpy(fd0,source.fd0,GetNProngs()*sizeof(Double32_t));
158 if(fPx) delete [] fPx;
159 fPx = new Double32_t[GetNProngs()];
160 if(fPy) delete [] fPy;
161 fPy = new Double32_t[GetNProngs()];
162 if(fPz) delete [] fPz;
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));
169 if(fPID) delete [] fPID;
170 fPID = new Double32_t[fNPID];
171 memcpy(fPID,source.fPID,fNPID*sizeof(Double32_t));
174 if(fDCA) delete [] fDCA;
175 fDCA = new Double32_t[fNDCA];
176 memcpy(fDCA,source.fDCA,fNDCA*sizeof(Double32_t));
181 //--------------------------------------------------------------------------
182 AliAODRecoDecay::~AliAODRecoDecay() {
184 // Default Destructor
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; }
194 //----------------------------------------------------------------------------
195 Double_t AliAODRecoDecay::Alpha() const
198 // Armenteros-Podolanski alpha for 2-prong decays
200 if(GetNProngs()!=2) {
201 printf("Can be called only for 2-prong decays");
202 return (Double_t)-99999.;
204 return 1.-2./(1.+QlProng(0)/QlProng(1));
206 //----------------------------------------------------------------------------
207 Double_t AliAODRecoDecay::DecayLength(Double_t point[3]) const
210 // Decay length assuming it is produced at "point" [cm]
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()));
219 //----------------------------------------------------------------------------
220 Double_t AliAODRecoDecay::DecayLengthXY(Double_t point[3]) const
223 // Decay length in XY assuming it is produced at "point" [cm]
225 return TMath::Sqrt((point[0]-GetSecVtxX())
226 *(point[0]-GetSecVtxX())
227 +(point[1]-GetSecVtxY())
228 *(point[1]-GetSecVtxY()));
230 //----------------------------------------------------------------------------
231 Double_t AliAODRecoDecay::CosPointingAngle(Double_t point[3]) const
234 // Cosine of pointing angle in space assuming it is produced at "point"
236 TVector3 mom(Px(),Py(),Pz());
237 TVector3 fline(GetSecVtxX()-point[0],
238 GetSecVtxY()-point[1],
239 GetSecVtxZ()-point[2]);
241 Double_t pta = mom.Angle(fline);
243 return TMath::Cos(pta);
245 //----------------------------------------------------------------------------
246 Double_t AliAODRecoDecay::CosPointingAngleXY(Double_t point[3]) const
249 // Cosine of pointing angle in transverse plane assuming it is produced
252 TVector3 momXY(Px(),Py(),0.);
253 TVector3 flineXY(GetSecVtxX()-point[0],
254 GetSecVtxY()-point[1],
257 Double_t ptaXY = momXY.Angle(flineXY);
259 return TMath::Cos(ptaXY);
261 //----------------------------------------------------------------------------
262 Double_t AliAODRecoDecay::CosThetaStar(Int_t ip,UInt_t pdgvtx,UInt_t pdgprong0,UInt_t pdgprong1) const
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
270 if(GetNProngs()!=2) {
271 printf("Can be called only for 2-prong decays");
272 return (Double_t)-99999.;
274 Double_t massvtx = TDatabasePDG::Instance()->GetParticle(pdgvtx)->Mass();
276 massp[0] = TDatabasePDG::Instance()->GetParticle(pdgprong0)->Mass();
277 massp[1] = TDatabasePDG::Instance()->GetParticle(pdgprong1)->Mass();
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);
281 Double_t beta = P()/E(pdgvtx);
282 Double_t gamma = E(pdgvtx)/massvtx;
284 Double_t cts = (QlProng(ip)/gamma-beta*TMath::Sqrt(pStar*pStar+massp[ip]*massp[ip]))/pStar;
288 //---------------------------------------------------------------------------
289 Double_t AliAODRecoDecay::Ct(UInt_t pdg,Double_t point[3]) const
292 // Decay time * c assuming it is produced at "point" [cm]
294 Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
295 return DecayLength(point)*mass/P();
297 //---------------------------------------------------------------------------
298 Double_t AliAODRecoDecay::E(UInt_t pdg) const
303 Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
304 return TMath::Sqrt(mass*mass+P()*P());
306 //---------------------------------------------------------------------------
307 Double_t AliAODRecoDecay::EProng(Int_t ip,UInt_t pdg) const
310 // Energy of ip-th prong
312 Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
313 return TMath::Sqrt(mass*mass+PProng(ip)*PProng(ip));
315 //--------------------------------------------------------------------------
316 Bool_t AliAODRecoDecay::GetCovarianceXYZPxPyPz(Double_t cv[21]) const {
318 // This function returns the global covariance matrix of the track params
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]
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.
334 for(j=0;j<21;j++) cv[j]=0.;
336 if(!GetNDaughters()) {
337 AliError("No daughters available");
342 AliAODVertex *secv=GetSecondaryVtx();
344 AliError("Vertex covariance matrix not available");
347 if(!secv->GetCovMatrix(v)) {
348 AliError("Vertex covariance matrix not available");
352 Double_t p[21]; for(j=0;j<21;j++) p[j]=0.;
354 for(Int_t i=1; i<GetNDaughters(); i++) {
355 AliVTrack *daugh = (AliVTrack*)GetDaughter(i);
357 if(!daugh->GetCovarianceXYZPxPyPz(dcov)) error=kTRUE;
358 for(j=0;j<21;j++) p[j] += dcov[j];
361 AliError("No covariance for at least one daughter")
365 for(j=0; j<21; j++) {
375 //----------------------------------------------------------------------------
376 UChar_t AliAODRecoDecay::GetITSClusterMap() const {
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)
383 if(!GetNDaughters()) {
384 AliError("No daughters available");
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)
390 for(Int_t i=0; i<GetNDaughters(); i++) {
391 AliVTrack *daugh = (AliVTrack*)GetDaughter(i);
392 if(!TESTBIT(daugh->GetITSClusterMap(),l)) bit=0;
394 if(bit) SETBIT(map,l);
399 //--------------------------------------------------------------------------
400 ULong_t AliAODRecoDecay::GetStatus() const {
402 // Same as for ITSClusterMap
406 if(!GetNDaughters()) {
407 AliError("No daughters available");
411 AliVTrack *daugh0 = (AliVTrack*)GetDaughter(0);
412 status = status&(daugh0->GetStatus());
414 for(Int_t i=1; i<GetNDaughters(); i++) {
415 AliVTrack *daugh = (AliVTrack*)GetDaughter(i);
416 status = status&(daugh->GetStatus());
421 //--------------------------------------------------------------------------
422 Double_t AliAODRecoDecay::ImpParXY(Double_t point[3]) const
425 // Impact parameter in the bending plane of the particle
428 Double_t k = -(GetSecVtxX()-point[0])*Px()-(GetSecVtxY()-point[1])*Py();
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);
440 //----------------------------------------------------------------------------
441 Double_t AliAODRecoDecay::InvMass(Int_t npdg,UInt_t *pdg) const
444 // Invariant mass for prongs mass hypotheses in pdg array
446 if(GetNProngs()!=npdg) {
447 printf("npdg != GetNProngs()");
448 return (Double_t)-99999.;
450 Double_t energysum = 0.;
452 for(Int_t i=0; i<GetNProngs(); i++) {
453 energysum += EProng(i,pdg[i]);
456 Double_t mass = TMath::Sqrt(energysum*energysum-P()*P());
460 //----------------------------------------------------------------------------
461 Double_t AliAODRecoDecay::InvMass2Prongs(Int_t ip1,Int_t ip2,
462 UInt_t pdg1,UInt_t pdg2) const
465 // 2-prong(ip1,ip2) invariant mass for prongs mass hypotheses in pdg1,2
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);
475 //---------------------------------------------------------------------------
476 void AliAODRecoDecay::Print(Option_t* /*option*/) const
479 // Print some information
481 printf("AliAODRecoDecay with %d prongs\n",GetNProngs());
482 printf("Secondary Vertex: (%f, %f, %f)\n",GetSecVtxX(),GetSecVtxY(),GetSecVtxZ());
486 //----------------------------------------------------------------------------
487 Double_t AliAODRecoDecay::ProngsRelAngle(Int_t ip1,Int_t ip2) const
490 // Relative angle between two prongs
492 TVector3 momA(PxProng(ip1),PyProng(ip1),PzProng(ip1));
493 TVector3 momB(PxProng(ip2),PyProng(ip2),PzProng(ip2));
495 Double_t angle = momA.Angle(momB);
499 //----------------------------------------------------------------------------
500 Double_t AliAODRecoDecay::QlProng(Int_t ip) const
503 // Longitudinal momentum of prong w.r.t. to total momentum
505 TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
506 TVector3 momTot(Px(),Py(),Pz());
508 return mom.Dot(momTot)/momTot.Mag();
510 //----------------------------------------------------------------------------
511 Double_t AliAODRecoDecay::QtProng(Int_t ip) const
514 // Transverse momentum of prong w.r.t. to total momentum
516 TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
517 TVector3 momTot(Px(),Py(),Pz());
519 return mom.Perp(momTot);
521 //----------------------------------------------------------------------------
522 Double_t AliAODRecoDecay::QlProngFlightLine(Int_t ip,Double_t point[3]) const
525 // Longitudinal momentum of prong w.r.t. to flight line between "point"
528 TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
529 TVector3 fline(GetSecVtxX()-point[0],
530 GetSecVtxY()-point[1],
531 GetSecVtxZ()-point[2]);
533 return mom.Dot(fline)/fline.Mag();
535 //----------------------------------------------------------------------------
536 Double_t AliAODRecoDecay::QtProngFlightLine(Int_t ip,Double_t point[3]) const
539 // Transverse momentum of prong w.r.t. to flight line between "point" and
542 TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
543 TVector3 fline(GetSecVtxX()-point[0],
544 GetSecVtxY()-point[1],
545 GetSecVtxZ()-point[2]);
547 return mom.Perp(fline);
549 //--------------------------------------------------------------------------