git://git.uio.no / u/mrichter/AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fix mismatched new/delete
[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>
ff7c57dd 25#include "AliVParticle.h"
7de7497b 26#include "AliAODRecoDecay.h"
27
28ClassImp(AliAODRecoDecay)
29
30//--------------------------------------------------------------------------
31AliAODRecoDecay::AliAODRecoDecay() :
ff7c57dd 32 AliVParticle(),
7de7497b 33 fSecondaryVtx(0x0),
a11de4a0 34 fOwnSecondaryVtx(0x0),
7de7497b 35 fCharge(0),
36 fNProngs(0), fNDCA(0), fNPID(0),
37 fPx(0x0), fPy(0x0), fPz(0x0),
38 fd0(0x0),
39 fDCA(0x0),
40 fPID(0x0),
41 fEventNumber(-1),fRunNumber(-1)
42{
43 //
44 // Default Constructor
45 //
7de7497b 46}
47//--------------------------------------------------------------------------
48AliAODRecoDecay::AliAODRecoDecay(AliAODVertex *vtx2,Int_t nprongs,
49 Short_t charge,
50 Double_t *px,Double_t *py,Double_t *pz,
51 Double_t *d0) :
ff7c57dd 52 AliVParticle(),
7de7497b 53 fSecondaryVtx(vtx2),
a11de4a0 54 fOwnSecondaryVtx(0x0),
7de7497b 55 fCharge(charge),
56 fNProngs(nprongs), fNDCA(0), fNPID(0),
57 fPx(0x0), fPy(0x0), fPz(0x0),
58 fd0(0x0),
59 fDCA(0x0),
60 fPID(0x0),
61 fEventNumber(-1),fRunNumber(-1)
62{
63 //
64 // Constructor with AliAODVertex for decay vertex
65 //
66
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++) {
72 fPx[i] = px[i];
73 fPy[i] = py[i];
74 fPz[i] = pz[i];
75 fd0[i] = d0[i];
76 }
77}
78//--------------------------------------------------------------------------
79AliAODRecoDecay::AliAODRecoDecay(AliAODVertex *vtx2,Int_t nprongs,
80 Short_t charge,
81 Double_t *d0) :
ff7c57dd 82 AliVParticle(),
7de7497b 83 fSecondaryVtx(vtx2),
a11de4a0 84 fOwnSecondaryVtx(0x0),
7de7497b 85 fCharge(charge),
86 fNProngs(nprongs), fNDCA(0), fNPID(0),
87 fPx(0x0), fPy(0x0), fPz(0x0),
88 fd0(0x0),
89 fDCA(0x0),
90 fPID(0x0),
91 fEventNumber(-1),fRunNumber(-1)
92{
93 //
94 // Constructor with AliAODVertex for decay vertex and without prongs momenta
95 //
96
97 fd0 = new Double_t[GetNProngs()];
98 for(Int_t i=0; i<GetNProngs(); i++) fd0[i] = d0[i];
99}
100//--------------------------------------------------------------------------
101AliAODRecoDecay::AliAODRecoDecay(const AliAODRecoDecay &source) :
ff7c57dd 102 AliVParticle(source),
7de7497b 103 fSecondaryVtx(source.fSecondaryVtx),
a11de4a0 104 fOwnSecondaryVtx(source.fOwnSecondaryVtx),
7de7497b 105 fCharge(source.fCharge),
106 fNProngs(source.fNProngs), fNDCA(source.fNDCA), fNPID(source.fNPID),
107 fPx(0x0), fPy(0x0), fPz(0x0),
108 fd0(0x0),
109 fDCA(0x0),
110 fPID(0x0),
111 fEventNumber(source.fEventNumber),fRunNumber(source.fRunNumber)
112{
113 //
114 // Copy constructor
115 //
116 if(source.GetNProngs()>0) {
117 fd0 = new Double_t[GetNProngs()];
118 memcpy(fd0,source.fd0,GetNProngs()*sizeof(Double_t));
119 if(source.fPx) {
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));
126 }
127 if(source.fPID) {
128 fPID = new Double_t[5*GetNProngs()];
129 memcpy(fPID,source.fPID,GetNProngs()*sizeof(Double_t));
130 }
131 if(source.fDCA) {
5cc73331 132 fDCA = new Double_t[GetNProngs()*(GetNProngs()-1)/2];
133 memcpy(fDCA,source.fDCA,(GetNProngs()*(GetNProngs()-1)/2)*sizeof(Double_t));
7de7497b 134 }
135 }
136}
137//--------------------------------------------------------------------------
138AliAODRecoDecay &AliAODRecoDecay::operator=(const AliAODRecoDecay &source)
139{
140 //
141 // assignment operator
142 //
143 if(&source == this) return *this;
144 fSecondaryVtx = source.fSecondaryVtx;
a11de4a0 145 fOwnSecondaryVtx = source.fOwnSecondaryVtx;
7de7497b 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) {
341952b8 153 if(fd0)delete [] fd0;
7de7497b 154 fd0 = new Double_t[GetNProngs()];
155 memcpy(fd0,source.fd0,GetNProngs()*sizeof(Double_t));
156 if(source.fPx) {
341952b8 157 if(fPx) delete [] fPx;
7de7497b 158 fPx = new Double_t[GetNProngs()];
341952b8 159 if(fPy) delete [] fPy;
7de7497b 160 fPy = new Double_t[GetNProngs()];
341952b8 161 if(fPz) delete [] fPz;
7de7497b 162 fPz = new Double_t[GetNProngs()];
163 memcpy(fPx,source.fPx,GetNProngs()*sizeof(Double_t));
164 memcpy(fPy,source.fPy,GetNProngs()*sizeof(Double_t));
165 memcpy(fPz,source.fPz,GetNProngs()*sizeof(Double_t));
166 }
167 if(source.fPID) {
341952b8 168 if(fPID) delete [] fPID;
7de7497b 169 fPID = new Double_t[5*GetNProngs()];
170 memcpy(fPID,source.fPID,GetNProngs()*sizeof(Double_t));
171 }
172 if(source.fDCA) {
341952b8 173 if(fDCA) delete [] fDCA;
5cc73331 174 fDCA = new Double_t[GetNProngs()*(GetNProngs()-1)/2];
175 memcpy(fDCA,source.fDCA,(GetNProngs()*(GetNProngs()-1)/2)*sizeof(Double32_t));
7de7497b 176 }
177 }
178 return *this;
179}
180//--------------------------------------------------------------------------
181AliAODRecoDecay::~AliAODRecoDecay() {
182 //
183 // Default Destructor
184 //
a11de4a0 185 if(fPx) { delete [] fPx; fPx=NULL; }
186 if(fPy) { delete [] fPy; fPy=NULL; }
187 if(fPz) { delete [] fPz; fPz=NULL; }
188 if(fd0) { delete [] fd0; fd0=NULL; }
189 if(fPID) { delete [] fPID; fPID=NULL; }
190 if(fDCA) { delete [] fDCA; fDCA=NULL; }
191 if(fOwnSecondaryVtx) { delete fOwnSecondaryVtx; fOwnSecondaryVtx=NULL; }
7de7497b 192}
193//----------------------------------------------------------------------------
194Double_t AliAODRecoDecay::Alpha() const
195{
196 //
197 // Armenteros-Podolanski alpha for 2-prong decays
198 //
199 if(GetNProngs()!=2) {
200 printf("Can be called only for 2-prong decays");
201 return (Double_t)-99999.;
202 }
203 return 1.-2./(1.+QlProng(0)/QlProng(1));
204}
205//----------------------------------------------------------------------------
206Double_t AliAODRecoDecay::DecayLength(Double_t point[3]) const
207{
208 //
209 // Decay length assuming it is produced at "point" [cm]
210 //
211 return TMath::Sqrt((point[0]-GetSecVtxX())
212 *(point[0]-GetSecVtxX())
213 +(point[1]-GetSecVtxY())
214 *(point[1]-GetSecVtxY())
215 +(point[2]-GetSecVtxZ())
216 *(point[2]-GetSecVtxZ()));
217}
218//----------------------------------------------------------------------------
219Double_t AliAODRecoDecay::DecayLengthXY(Double_t point[3]) const
220{
221 //
222 // Decay length in XY assuming it is produced at "point" [cm]
223 //
224 return TMath::Sqrt((point[0]-GetSecVtxX())
225 *(point[0]-GetSecVtxX())
226 +(point[1]-GetSecVtxY())
227 *(point[1]-GetSecVtxY()));
228}
229//----------------------------------------------------------------------------
230Double_t AliAODRecoDecay::CosPointingAngle(Double_t point[3]) const
231{
232 //
233 // Cosine of pointing angle in space assuming it is produced at "point"
234 //
235 TVector3 mom(Px(),Py(),Pz());
236 TVector3 fline(GetSecVtxX()-point[0],
237 GetSecVtxY()-point[1],
238 GetSecVtxZ()-point[2]);
239
240 Double_t pta = mom.Angle(fline);
241
242 return TMath::Cos(pta);
243}
244//----------------------------------------------------------------------------
245Double_t AliAODRecoDecay::CosPointingAngleXY(Double_t point[3]) const
246{
247 //
248 // Cosine of pointing angle in transverse plane assuming it is produced
249 // at "point"
250 //
251 TVector3 momXY(Px(),Py(),0.);
252 TVector3 flineXY(GetSecVtxX()-point[0],
253 GetSecVtxY()-point[1],
254 0.);
255
256 Double_t ptaXY = momXY.Angle(flineXY);
257
258 return TMath::Cos(ptaXY);
259}
260//----------------------------------------------------------------------------
261Double_t AliAODRecoDecay::CosThetaStar(Int_t ip,UInt_t pdgvtx,UInt_t pdgprong0,UInt_t pdgprong1) const
262{
263 //
264 // Only for 2-prong decays:
265 // Cosine of decay angle (theta*) in the rest frame of the mother particle
266 // for prong ip (0 or 1) with mass hypotheses pdgvtx for mother particle,
267 // pdgprong0 for prong 0 and pdgprong1 for prong1
268 //
269 if(GetNProngs()!=2) {
270 printf("Can be called only for 2-prong decays");
271 return (Double_t)-99999.;
272 }
273 Double_t massvtx = TDatabasePDG::Instance()->GetParticle(pdgvtx)->Mass();
274 Double_t massp[2];
275 massp[0] = TDatabasePDG::Instance()->GetParticle(pdgprong0)->Mass();
276 massp[1] = TDatabasePDG::Instance()->GetParticle(pdgprong1)->Mass();
277
278 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);
279
280 Double_t beta = P()/E(pdgvtx);
281 Double_t gamma = E(pdgvtx)/massvtx;
282
283 Double_t cts = (QlProng(ip)/gamma-beta*TMath::Sqrt(pStar*pStar+massp[ip]*massp[ip]))/pStar;
284
285 return cts;
286}
287//---------------------------------------------------------------------------
288Double_t AliAODRecoDecay::Ct(UInt_t pdg,Double_t point[3]) const
289{
290 //
291 // Decay time * c assuming it is produced at "point" [cm]
292 //
293 Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
294 return DecayLength(point)*mass/P();
295}
296//---------------------------------------------------------------------------
297Double_t AliAODRecoDecay::E(UInt_t pdg) const
298{
299 //
300 // Energy
301 //
302 Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
303 return TMath::Sqrt(mass*mass+P()*P());
304}
305//---------------------------------------------------------------------------
306Double_t AliAODRecoDecay::EProng(Int_t ip,UInt_t pdg) const
307{
308 //
309 // Energy of ip-th prong
310 //
311 Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
312 return TMath::Sqrt(mass*mass+PProng(ip)*PProng(ip));
313}
314//---------------------------------------------------------------------------
315/*Int_t AliAODRecoDecay::GetIndexProng(Int_t ip) const
316{
317 //
318 // Index of prong ip
319 //
320 if(!GetNProngs()) return 999999;
5cc73331 321UShort_t *indices = GetSecondaryVtx()->GetIndices();
7de7497b 322 return indices[ip];
323}*/
324//----------------------------------------------------------------------------
325Double_t AliAODRecoDecay::ImpParXY(Double_t point[3]) const
326{
327 //
328 // Impact parameter in the bending plane of the particle
329 // w.r.t. to "point"
330 //
331 Double_t k = -(GetSecVtxX()-point[0])*Px()-(GetSecVtxY()-point[1])*Py();
332 k /= Pt()*Pt();
333 Double_t dx = GetSecVtxX()-point[0]+k*Px();
334 Double_t dy = GetSecVtxY()-point[1]+k*Py();
335 Double_t absImpPar = TMath::Sqrt(dx*dx+dy*dy);
336 TVector3 mom(Px(),Py(),Pz());
337 TVector3 fline(GetSecVtxX()-point[0],
338 GetSecVtxY()-point[1],
339 GetSecVtxZ()-point[2]);
340 TVector3 cross = mom.Cross(fline);
341 return (cross.Z()>0. ? absImpPar : -absImpPar);
342}
343//----------------------------------------------------------------------------
344Double_t AliAODRecoDecay::InvMass(Int_t npdg,UInt_t *pdg) const
345{
346 //
347 // Invariant mass for prongs mass hypotheses in pdg array
348 //
349 if(GetNProngs()!=npdg) {
350 printf("npdg != GetNProngs()");
351 return (Double_t)-99999.;
352 }
353 Double_t energysum = 0.;
354
355 for(Int_t i=0; i<GetNProngs(); i++) {
356 energysum += EProng(i,pdg[i]);
357 }
358
359 Double_t mass = TMath::Sqrt(energysum*energysum-P()*P());
360
361 return mass;
362}
363//----------------------------------------------------------------------------
364Double_t AliAODRecoDecay::InvMass2Prongs(Int_t ip1,Int_t ip2,
365 UInt_t pdg1,UInt_t pdg2) const
366{
367 //
368 // 2-prong(ip1,ip2) invariant mass for prongs mass hypotheses in pdg1,2
369 //
370 Double_t energysum = EProng(ip1,pdg1) + EProng(ip2,pdg2);
371 Double_t psum2 = (PxProng(ip1)+PxProng(ip2))*(PxProng(ip1)+PxProng(ip2))
372 +(PyProng(ip1)+PyProng(ip2))*(PyProng(ip1)+PyProng(ip2))
373 +(PzProng(ip1)+PzProng(ip2))*(PzProng(ip1)+PzProng(ip2));
374 Double_t mass = TMath::Sqrt(energysum*energysum-psum2);
375
376 return mass;
377}
378//---------------------------------------------------------------------------
379void AliAODRecoDecay::Print(Option_t* /*option*/) const
380{
381 //
382 // Print some information
383 //
384 printf("AliAODRecoDecay with %d prongs\n",GetNProngs());
385 printf("Secondary Vertex: (%f, %f, %f)\n",GetSecVtxX(),GetSecVtxY(),GetSecVtxZ());
386
387 return;
388}
389//----------------------------------------------------------------------------
390Double_t AliAODRecoDecay::ProngsRelAngle(Int_t ip1,Int_t ip2) const
391{
392 //
393 // Relative angle between two prongs
394 //
395 TVector3 momA(PxProng(ip1),PyProng(ip1),PzProng(ip1));
396 TVector3 momB(PxProng(ip2),PyProng(ip2),PzProng(ip2));
397
398 Double_t angle = momA.Angle(momB);
399
400 return angle;
401}
402//----------------------------------------------------------------------------
403Double_t AliAODRecoDecay::QlProng(Int_t ip) const
404{
405 //
406 // Longitudinal momentum of prong w.r.t. to total momentum
407 //
408 TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
409 TVector3 momTot(Px(),Py(),Pz());
410
411 return mom.Dot(momTot)/momTot.Mag();
412}
413//----------------------------------------------------------------------------
414Double_t AliAODRecoDecay::QtProng(Int_t ip) const
415{
416 //
417 // Transverse momentum of prong w.r.t. to total momentum
418 //
419 TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
420 TVector3 momTot(Px(),Py(),Pz());
421
422 return mom.Perp(momTot);
423}
424//----------------------------------------------------------------------------
425Double_t AliAODRecoDecay::QlProngFlightLine(Int_t ip,Double_t point[3]) const
426{
427 //
428 // Longitudinal momentum of prong w.r.t. to flight line between "point"
429 // and fSecondaryVtx
430 //
431 TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
432 TVector3 fline(GetSecVtxX()-point[0],
433 GetSecVtxY()-point[1],
434 GetSecVtxZ()-point[2]);
435
436 return mom.Dot(fline)/fline.Mag();
437}
438//----------------------------------------------------------------------------
439Double_t AliAODRecoDecay::QtProngFlightLine(Int_t ip,Double_t point[3]) const
440{
441 //
442 // Transverse momentum of prong w.r.t. to flight line between "point" and
443 // fSecondaryVtx
444 //
445 TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
446 TVector3 fline(GetSecVtxX()-point[0],
447 GetSecVtxY()-point[1],
448 GetSecVtxZ()-point[2]);
449
450 return mom.Perp(fline);
451}
452//--------------------------------------------------------------------------