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