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