Rename AliVirtualParticle to AliVParticle (Markus)
[u/mrichter/AliRoot.git] / STEER / AliAODRecoDecay.cxx
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 "AliVParticle.h"
26 #include "AliAODRecoDecay.h"
27
28 ClassImp(AliAODRecoDecay)
29
30 //--------------------------------------------------------------------------
31 AliAODRecoDecay::AliAODRecoDecay() :
32   AliVParticle(),
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 //--------------------------------------------------------------------------
48 AliAODRecoDecay::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   AliVParticle(),
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 //--------------------------------------------------------------------------
78 AliAODRecoDecay::AliAODRecoDecay(AliAODVertex *vtx2,Int_t nprongs,
79                                  Short_t charge,
80                                  Double_t *d0) :
81   AliVParticle(),
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 //--------------------------------------------------------------------------
99 AliAODRecoDecay::AliAODRecoDecay(const AliAODRecoDecay &source) :
100   AliVParticle(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 //--------------------------------------------------------------------------
135 AliAODRecoDecay &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 //--------------------------------------------------------------------------
171 AliAODRecoDecay::~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 //----------------------------------------------------------------------------
184 Double_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 //----------------------------------------------------------------------------
196 Double_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 //----------------------------------------------------------------------------
209 Double_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 //----------------------------------------------------------------------------
220 Double_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 //----------------------------------------------------------------------------
235 Double_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 //----------------------------------------------------------------------------
251 Double_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 //---------------------------------------------------------------------------
278 Double_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 //---------------------------------------------------------------------------
287 Double_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 //---------------------------------------------------------------------------
296 Double_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 //----------------------------------------------------------------------------
315 Double_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 //----------------------------------------------------------------------------
334 Double_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 //----------------------------------------------------------------------------
354 Double_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 //---------------------------------------------------------------------------
369 void 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 //----------------------------------------------------------------------------
380 Double_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 //----------------------------------------------------------------------------
393 Double_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 //----------------------------------------------------------------------------
404 Double_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 //----------------------------------------------------------------------------
415 Double_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 //----------------------------------------------------------------------------
429 Double_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 //--------------------------------------------------------------------------