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>
25 #include "AliVParticle.h"
26 #include "AliAODRecoDecay.h"
28 ClassImp(AliAODRecoDecay)
30 //--------------------------------------------------------------------------
31 AliAODRecoDecay::AliAODRecoDecay() :
34 fOwnSecondaryVtx(0x0),
36 fNProngs(0), fNDCA(0), fNPID(0),
37 fPx(0x0), fPy(0x0), fPz(0x0),
41 fEventNumber(-1),fRunNumber(-1)
44 // Default Constructor
47 //--------------------------------------------------------------------------
48 AliAODRecoDecay::AliAODRecoDecay(AliAODVertex *vtx2,Int_t nprongs,
50 Double_t *px,Double_t *py,Double_t *pz,
54 fOwnSecondaryVtx(0x0),
56 fNProngs(nprongs), fNDCA(0), fNPID(0),
57 fPx(0x0), fPy(0x0), fPz(0x0),
61 fEventNumber(-1),fRunNumber(-1)
64 // Constructor with AliAODVertex for decay vertex
67 fPx = new Double_t[GetNProngs()];
68 fPy = new Double_t[GetNProngs()];
69 fPz = new Double_t[GetNProngs()];
70 fd0 = new Double_t[GetNProngs()];
71 for(Int_t i=0; i<GetNProngs(); i++) {
78 //--------------------------------------------------------------------------
79 AliAODRecoDecay::AliAODRecoDecay(AliAODVertex *vtx2,Int_t nprongs,
84 fOwnSecondaryVtx(0x0),
86 fNProngs(nprongs), fNDCA(0), fNPID(0),
87 fPx(0x0), fPy(0x0), fPz(0x0),
91 fEventNumber(-1),fRunNumber(-1)
94 // Constructor with AliAODVertex for decay vertex and without prongs momenta
97 fd0 = new Double_t[GetNProngs()];
98 for(Int_t i=0; i<GetNProngs(); i++) fd0[i] = d0[i];
100 //--------------------------------------------------------------------------
101 AliAODRecoDecay::AliAODRecoDecay(const AliAODRecoDecay &source) :
102 AliVParticle(source),
103 fSecondaryVtx(source.fSecondaryVtx),
104 fOwnSecondaryVtx(source.fOwnSecondaryVtx),
105 fCharge(source.fCharge),
106 fNProngs(source.fNProngs), fNDCA(source.fNDCA), fNPID(source.fNPID),
107 fPx(0x0), fPy(0x0), fPz(0x0),
111 fEventNumber(source.fEventNumber),fRunNumber(source.fRunNumber)
116 if(source.GetNProngs()>0) {
117 fd0 = new Double_t[GetNProngs()];
118 memcpy(fd0,source.fd0,GetNProngs()*sizeof(Double_t));
120 fPx = new Double_t[GetNProngs()];
121 fPy = new Double_t[GetNProngs()];
122 fPz = new Double_t[GetNProngs()];
123 memcpy(fPx,source.fPx,GetNProngs()*sizeof(Double_t));
124 memcpy(fPy,source.fPy,GetNProngs()*sizeof(Double_t));
125 memcpy(fPz,source.fPz,GetNProngs()*sizeof(Double_t));
128 fPID = new Double_t[5*GetNProngs()];
129 memcpy(fPID,source.fPID,GetNProngs()*sizeof(Double_t));
132 fDCA = new Double_t[GetNProngs()*(GetNProngs()-1)/2];
133 memcpy(fDCA,source.fDCA,(GetNProngs()*(GetNProngs()-1)/2)*sizeof(Double_t));
137 //--------------------------------------------------------------------------
138 AliAODRecoDecay &AliAODRecoDecay::operator=(const AliAODRecoDecay &source)
141 // assignment operator
143 if(&source == this) return *this;
144 fSecondaryVtx = source.fSecondaryVtx;
145 fOwnSecondaryVtx = source.fOwnSecondaryVtx;
146 fCharge = source.fCharge;
147 fNProngs = source.fNProngs;
148 fNDCA = source.fNDCA;
149 fNPID = source.fNPID;
150 fEventNumber = source.fEventNumber;
151 fRunNumber = source.fRunNumber;
152 if(source.GetNProngs()>0) {
153 fd0 = new Double_t[GetNProngs()];
154 memcpy(fd0,source.fd0,GetNProngs()*sizeof(Double_t));
156 fPx = new Double_t[GetNProngs()];
157 fPy = new Double_t[GetNProngs()];
158 fPz = new Double_t[GetNProngs()];
159 memcpy(fPx,source.fPx,GetNProngs()*sizeof(Double_t));
160 memcpy(fPy,source.fPy,GetNProngs()*sizeof(Double_t));
161 memcpy(fPz,source.fPz,GetNProngs()*sizeof(Double_t));
164 fPID = new Double_t[5*GetNProngs()];
165 memcpy(fPID,source.fPID,GetNProngs()*sizeof(Double_t));
168 fDCA = new Double_t[GetNProngs()*(GetNProngs()-1)/2];
169 memcpy(fDCA,source.fDCA,(GetNProngs()*(GetNProngs()-1)/2)*sizeof(Double32_t));
174 //--------------------------------------------------------------------------
175 AliAODRecoDecay::~AliAODRecoDecay() {
177 // Default Destructor
179 if(fPx) { delete [] fPx; fPx=NULL; }
180 if(fPy) { delete [] fPy; fPy=NULL; }
181 if(fPz) { delete [] fPz; fPz=NULL; }
182 if(fd0) { delete [] fd0; fd0=NULL; }
183 if(fPID) { delete [] fPID; fPID=NULL; }
184 if(fDCA) { delete [] fDCA; fDCA=NULL; }
185 if(fOwnSecondaryVtx) { delete fOwnSecondaryVtx; fOwnSecondaryVtx=NULL; }
187 //----------------------------------------------------------------------------
188 Double_t AliAODRecoDecay::Alpha() const
191 // Armenteros-Podolanski alpha for 2-prong decays
193 if(GetNProngs()!=2) {
194 printf("Can be called only for 2-prong decays");
195 return (Double_t)-99999.;
197 return 1.-2./(1.+QlProng(0)/QlProng(1));
199 //----------------------------------------------------------------------------
200 Double_t AliAODRecoDecay::DecayLength(Double_t point[3]) const
203 // Decay length assuming it is produced at "point" [cm]
205 return TMath::Sqrt((point[0]-GetSecVtxX())
206 *(point[0]-GetSecVtxX())
207 +(point[1]-GetSecVtxY())
208 *(point[1]-GetSecVtxY())
209 +(point[2]-GetSecVtxZ())
210 *(point[2]-GetSecVtxZ()));
212 //----------------------------------------------------------------------------
213 Double_t AliAODRecoDecay::DecayLengthXY(Double_t point[3]) const
216 // Decay length in XY assuming it is produced at "point" [cm]
218 return TMath::Sqrt((point[0]-GetSecVtxX())
219 *(point[0]-GetSecVtxX())
220 +(point[1]-GetSecVtxY())
221 *(point[1]-GetSecVtxY()));
223 //----------------------------------------------------------------------------
224 Double_t AliAODRecoDecay::CosPointingAngle(Double_t point[3]) const
227 // Cosine of pointing angle in space assuming it is produced at "point"
229 TVector3 mom(Px(),Py(),Pz());
230 TVector3 fline(GetSecVtxX()-point[0],
231 GetSecVtxY()-point[1],
232 GetSecVtxZ()-point[2]);
234 Double_t pta = mom.Angle(fline);
236 return TMath::Cos(pta);
238 //----------------------------------------------------------------------------
239 Double_t AliAODRecoDecay::CosPointingAngleXY(Double_t point[3]) const
242 // Cosine of pointing angle in transverse plane assuming it is produced
245 TVector3 momXY(Px(),Py(),0.);
246 TVector3 flineXY(GetSecVtxX()-point[0],
247 GetSecVtxY()-point[1],
250 Double_t ptaXY = momXY.Angle(flineXY);
252 return TMath::Cos(ptaXY);
254 //----------------------------------------------------------------------------
255 Double_t AliAODRecoDecay::CosThetaStar(Int_t ip,UInt_t pdgvtx,UInt_t pdgprong0,UInt_t pdgprong1) const
258 // Only for 2-prong decays:
259 // Cosine of decay angle (theta*) in the rest frame of the mother particle
260 // for prong ip (0 or 1) with mass hypotheses pdgvtx for mother particle,
261 // pdgprong0 for prong 0 and pdgprong1 for prong1
263 if(GetNProngs()!=2) {
264 printf("Can be called only for 2-prong decays");
265 return (Double_t)-99999.;
267 Double_t massvtx = TDatabasePDG::Instance()->GetParticle(pdgvtx)->Mass();
269 massp[0] = TDatabasePDG::Instance()->GetParticle(pdgprong0)->Mass();
270 massp[1] = TDatabasePDG::Instance()->GetParticle(pdgprong1)->Mass();
272 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);
274 Double_t beta = P()/E(pdgvtx);
275 Double_t gamma = E(pdgvtx)/massvtx;
277 Double_t cts = (QlProng(ip)/gamma-beta*TMath::Sqrt(pStar*pStar+massp[ip]*massp[ip]))/pStar;
281 //---------------------------------------------------------------------------
282 Double_t AliAODRecoDecay::Ct(UInt_t pdg,Double_t point[3]) const
285 // Decay time * c assuming it is produced at "point" [cm]
287 Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
288 return DecayLength(point)*mass/P();
290 //---------------------------------------------------------------------------
291 Double_t AliAODRecoDecay::E(UInt_t pdg) const
296 Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
297 return TMath::Sqrt(mass*mass+P()*P());
299 //---------------------------------------------------------------------------
300 Double_t AliAODRecoDecay::EProng(Int_t ip,UInt_t pdg) const
303 // Energy of ip-th prong
305 Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
306 return TMath::Sqrt(mass*mass+PProng(ip)*PProng(ip));
308 //---------------------------------------------------------------------------
309 /*Int_t AliAODRecoDecay::GetIndexProng(Int_t ip) const
314 if(!GetNProngs()) return 999999;
315 UShort_t *indices = GetSecondaryVtx()->GetIndices();
318 //----------------------------------------------------------------------------
319 Double_t AliAODRecoDecay::ImpParXY(Double_t point[3]) const
322 // Impact parameter in the bending plane of the particle
325 Double_t k = -(GetSecVtxX()-point[0])*Px()-(GetSecVtxY()-point[1])*Py();
327 Double_t dx = GetSecVtxX()-point[0]+k*Px();
328 Double_t dy = GetSecVtxY()-point[1]+k*Py();
329 Double_t absImpPar = TMath::Sqrt(dx*dx+dy*dy);
330 TVector3 mom(Px(),Py(),Pz());
331 TVector3 fline(GetSecVtxX()-point[0],
332 GetSecVtxY()-point[1],
333 GetSecVtxZ()-point[2]);
334 TVector3 cross = mom.Cross(fline);
335 return (cross.Z()>0. ? absImpPar : -absImpPar);
337 //----------------------------------------------------------------------------
338 Double_t AliAODRecoDecay::InvMass(Int_t npdg,UInt_t *pdg) const
341 // Invariant mass for prongs mass hypotheses in pdg array
343 if(GetNProngs()!=npdg) {
344 printf("npdg != GetNProngs()");
345 return (Double_t)-99999.;
347 Double_t energysum = 0.;
349 for(Int_t i=0; i<GetNProngs(); i++) {
350 energysum += EProng(i,pdg[i]);
353 Double_t mass = TMath::Sqrt(energysum*energysum-P()*P());
357 //----------------------------------------------------------------------------
358 Double_t AliAODRecoDecay::InvMass2Prongs(Int_t ip1,Int_t ip2,
359 UInt_t pdg1,UInt_t pdg2) const
362 // 2-prong(ip1,ip2) invariant mass for prongs mass hypotheses in pdg1,2
364 Double_t energysum = EProng(ip1,pdg1) + EProng(ip2,pdg2);
365 Double_t psum2 = (PxProng(ip1)+PxProng(ip2))*(PxProng(ip1)+PxProng(ip2))
366 +(PyProng(ip1)+PyProng(ip2))*(PyProng(ip1)+PyProng(ip2))
367 +(PzProng(ip1)+PzProng(ip2))*(PzProng(ip1)+PzProng(ip2));
368 Double_t mass = TMath::Sqrt(energysum*energysum-psum2);
372 //---------------------------------------------------------------------------
373 void AliAODRecoDecay::Print(Option_t* /*option*/) const
376 // Print some information
378 printf("AliAODRecoDecay with %d prongs\n",GetNProngs());
379 printf("Secondary Vertex: (%f, %f, %f)\n",GetSecVtxX(),GetSecVtxY(),GetSecVtxZ());
383 //----------------------------------------------------------------------------
384 Double_t AliAODRecoDecay::ProngsRelAngle(Int_t ip1,Int_t ip2) const
387 // Relative angle between two prongs
389 TVector3 momA(PxProng(ip1),PyProng(ip1),PzProng(ip1));
390 TVector3 momB(PxProng(ip2),PyProng(ip2),PzProng(ip2));
392 Double_t angle = momA.Angle(momB);
396 //----------------------------------------------------------------------------
397 Double_t AliAODRecoDecay::QlProng(Int_t ip) const
400 // Longitudinal momentum of prong w.r.t. to total momentum
402 TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
403 TVector3 momTot(Px(),Py(),Pz());
405 return mom.Dot(momTot)/momTot.Mag();
407 //----------------------------------------------------------------------------
408 Double_t AliAODRecoDecay::QtProng(Int_t ip) const
411 // Transverse momentum of prong w.r.t. to total momentum
413 TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
414 TVector3 momTot(Px(),Py(),Pz());
416 return mom.Perp(momTot);
418 //----------------------------------------------------------------------------
419 Double_t AliAODRecoDecay::QlProngFlightLine(Int_t ip,Double_t point[3]) const
422 // Longitudinal momentum of prong w.r.t. to flight line between "point"
425 TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
426 TVector3 fline(GetSecVtxX()-point[0],
427 GetSecVtxY()-point[1],
428 GetSecVtxZ()-point[2]);
430 return mom.Dot(fline)/fline.Mag();
432 //----------------------------------------------------------------------------
433 Double_t AliAODRecoDecay::QtProngFlightLine(Int_t ip,Double_t point[3]) const
436 // Transverse momentum of prong w.r.t. to flight line between "point" and
439 TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
440 TVector3 fline(GetSecVtxX()-point[0],
441 GetSecVtxY()-point[1],
442 GetSecVtxZ()-point[2]);
444 return mom.Perp(fline);
446 //--------------------------------------------------------------------------