Some init values and deletes
[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   fOwnSecondaryVtx(0x0),
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   //
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   fOwnSecondaryVtx(0x0),
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 //--------------------------------------------------------------------------
79 AliAODRecoDecay::AliAODRecoDecay(AliAODVertex *vtx2,Int_t nprongs,
80                                  Short_t charge,
81                                  Double_t *d0) :
82   AliVParticle(),
83   fSecondaryVtx(vtx2),
84   fOwnSecondaryVtx(0x0),
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 //--------------------------------------------------------------------------
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),
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) {
132       fDCA = new Double_t[GetNProngs()*(GetNProngs()-1)/2];
133       memcpy(fDCA,source.fDCA,(GetNProngs()*(GetNProngs()-1)/2)*sizeof(Double_t));
134     }
135   }
136 }
137 //--------------------------------------------------------------------------
138 AliAODRecoDecay &AliAODRecoDecay::operator=(const AliAODRecoDecay &source)
139 {
140   //
141   // assignment operator
142   //
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     if(fd0)delete [] fd0; 
154     fd0 = new Double_t[GetNProngs()];
155     memcpy(fd0,source.fd0,GetNProngs()*sizeof(Double_t));
156     if(source.fPx) {
157       if(fPx) delete [] fPx; 
158       fPx = new Double_t[GetNProngs()];
159       if(fPy) delete [] fPy; 
160       fPy = new Double_t[GetNProngs()];
161       if(fPz) delete [] fPz; 
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) {
168       if(fPID) delete [] fPID; 
169       fPID = new Double_t[5*GetNProngs()];
170       memcpy(fPID,source.fPID,GetNProngs()*sizeof(Double_t));
171     }
172     if(source.fDCA) {
173       if(fDCA) delete [] fDCA; 
174       fDCA = new Double_t[GetNProngs()*(GetNProngs()-1)/2];
175       memcpy(fDCA,source.fDCA,(GetNProngs()*(GetNProngs()-1)/2)*sizeof(Double32_t));
176     }
177   }
178   return *this;
179 }
180 //--------------------------------------------------------------------------
181 AliAODRecoDecay::~AliAODRecoDecay() {
182   //  
183   // Default Destructor
184   //
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; }
192 }
193 //----------------------------------------------------------------------------
194 Double_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 //----------------------------------------------------------------------------
206 Double_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 //----------------------------------------------------------------------------
219 Double_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 //----------------------------------------------------------------------------
230 Double_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 //----------------------------------------------------------------------------
245 Double_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 //----------------------------------------------------------------------------
261 Double_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 //---------------------------------------------------------------------------
288 Double_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 //---------------------------------------------------------------------------
297 Double_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 //---------------------------------------------------------------------------
306 Double_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;
321 UShort_t *indices = GetSecondaryVtx()->GetIndices();
322   return indices[ip];
323 }*/
324 //----------------------------------------------------------------------------
325 Double_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 //----------------------------------------------------------------------------
344 Double_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 //----------------------------------------------------------------------------
364 Double_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 //---------------------------------------------------------------------------
379 void 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 //----------------------------------------------------------------------------
390 Double_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 //----------------------------------------------------------------------------
403 Double_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 //----------------------------------------------------------------------------
414 Double_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 //----------------------------------------------------------------------------
425 Double_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 //----------------------------------------------------------------------------
439 Double_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 //--------------------------------------------------------------------------