Deriving from VTrack now. (A. Dainese)
[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 "AliLog.h"
26 #include "AliVTrack.h"
27 #include "AliAODRecoDecay.h"
28
29 ClassImp(AliAODRecoDecay)
30
31 //--------------------------------------------------------------------------
32 AliAODRecoDecay::AliAODRecoDecay() :
33   AliVTrack(),
34   fSecondaryVtx(0x0),
35   fOwnSecondaryVtx(0x0),
36   fCharge(0),
37   fNProngs(0), fNDCA(0), fNPID(0),
38   fPx(0x0), fPy(0x0), fPz(0x0),
39   fd0(0x0),
40   fDCA(0x0),
41   fPID(0x0), 
42   fEventNumber(-1),fRunNumber(-1)
43 {
44   //
45   // Default Constructor
46   //
47 }
48 //--------------------------------------------------------------------------
49 AliAODRecoDecay::AliAODRecoDecay(AliAODVertex *vtx2,Int_t nprongs,
50                                  Short_t charge,
51                                  Double_t *px,Double_t *py,Double_t *pz,
52                                  Double_t *d0) :
53   AliVTrack(),
54   fSecondaryVtx(vtx2),
55   fOwnSecondaryVtx(0x0),
56   fCharge(charge),
57   fNProngs(nprongs), fNDCA(0), fNPID(0),
58   fPx(0x0), fPy(0x0), fPz(0x0),
59   fd0(0x0),
60   fDCA(0x0),
61   fPID(0x0), 
62   fEventNumber(-1),fRunNumber(-1)
63 {
64   //
65   // Constructor with AliAODVertex for decay vertex
66   //
67
68   fPx = new Double_t[GetNProngs()];
69   fPy = new Double_t[GetNProngs()];
70   fPz = new Double_t[GetNProngs()];
71   fd0 = new Double_t[GetNProngs()];
72   for(Int_t i=0; i<GetNProngs(); i++) {
73     fPx[i] = px[i];
74     fPy[i] = py[i];
75     fPz[i] = pz[i];
76     fd0[i] = d0[i];
77   }
78 }
79 //--------------------------------------------------------------------------
80 AliAODRecoDecay::AliAODRecoDecay(AliAODVertex *vtx2,Int_t nprongs,
81                                  Short_t charge,
82                                  Double_t *d0) :
83   AliVTrack(),
84   fSecondaryVtx(vtx2),
85   fOwnSecondaryVtx(0x0),
86   fCharge(charge),
87   fNProngs(nprongs), fNDCA(0), fNPID(0),
88   fPx(0x0), fPy(0x0), fPz(0x0),
89   fd0(0x0),
90   fDCA(0x0),
91   fPID(0x0), 
92   fEventNumber(-1),fRunNumber(-1)
93 {
94   //
95   // Constructor with AliAODVertex for decay vertex and without prongs momenta
96   //
97
98   fd0 = new Double_t[GetNProngs()];
99   for(Int_t i=0; i<GetNProngs(); i++) fd0[i] = d0[i];
100 }
101 //--------------------------------------------------------------------------
102 AliAODRecoDecay::AliAODRecoDecay(const AliAODRecoDecay &source) :
103   AliVTrack(source),
104   fSecondaryVtx(source.fSecondaryVtx),
105   fOwnSecondaryVtx(source.fOwnSecondaryVtx),
106   fCharge(source.fCharge),
107   fNProngs(source.fNProngs), fNDCA(source.fNDCA), fNPID(source.fNPID),
108   fPx(0x0), fPy(0x0), fPz(0x0),
109   fd0(0x0), 
110   fDCA(0x0),
111   fPID(0x0), 
112   fEventNumber(source.fEventNumber),fRunNumber(source.fRunNumber)
113 {
114   //
115   // Copy constructor
116   //
117   if(source.GetNProngs()>0) {
118     fd0 = new Double32_t[GetNProngs()];
119     memcpy(fd0,source.fd0,GetNProngs()*sizeof(Double32_t));
120     if(source.fPx) {
121       fPx = new Double32_t[GetNProngs()];
122       fPy = new Double32_t[GetNProngs()];
123       fPz = new Double32_t[GetNProngs()];
124       memcpy(fPx,source.fPx,GetNProngs()*sizeof(Double32_t));
125       memcpy(fPy,source.fPy,GetNProngs()*sizeof(Double32_t));
126       memcpy(fPz,source.fPz,GetNProngs()*sizeof(Double32_t));
127     }
128     if(source.fPID) {
129       fPID = new Double32_t[fNPID];
130       memcpy(fPID,source.fPID,fNPID*sizeof(Double32_t));
131     }
132     if(source.fDCA) {
133       fDCA = new Double32_t[fNDCA];
134       memcpy(fDCA,source.fDCA,fNDCA*sizeof(Double32_t));
135     }
136   }
137 }
138 //--------------------------------------------------------------------------
139 AliAODRecoDecay &AliAODRecoDecay::operator=(const AliAODRecoDecay &source)
140 {
141   //
142   // assignment operator
143   //
144   if(&source == this) return *this;
145   fSecondaryVtx = source.fSecondaryVtx;
146   fOwnSecondaryVtx = source.fOwnSecondaryVtx;
147   fCharge = source.fCharge;
148   fNProngs = source.fNProngs;
149   fNDCA = source.fNDCA;
150   fNPID = source.fNPID;
151   fEventNumber = source.fEventNumber;
152   fRunNumber = source.fRunNumber;
153   if(source.GetNProngs()>0) {
154     if(fd0)delete [] fd0; 
155     fd0 = new Double32_t[GetNProngs()];
156     memcpy(fd0,source.fd0,GetNProngs()*sizeof(Double32_t));
157     if(source.fPx) {
158       if(fPx) delete [] fPx; 
159       fPx = new Double32_t[GetNProngs()];
160       if(fPy) delete [] fPy; 
161       fPy = new Double32_t[GetNProngs()];
162       if(fPz) delete [] fPz; 
163       fPz = new Double32_t[GetNProngs()];
164       memcpy(fPx,source.fPx,GetNProngs()*sizeof(Double32_t));
165       memcpy(fPy,source.fPy,GetNProngs()*sizeof(Double32_t));
166       memcpy(fPz,source.fPz,GetNProngs()*sizeof(Double32_t));
167     }
168     if(source.fPID) {
169       if(fPID) delete [] fPID; 
170       fPID = new Double32_t[fNPID];
171       memcpy(fPID,source.fPID,fNPID*sizeof(Double32_t));
172     }
173     if(source.fDCA) {
174       if(fDCA) delete [] fDCA; 
175       fDCA = new Double32_t[fNDCA];
176       memcpy(fDCA,source.fDCA,fNDCA*sizeof(Double32_t));
177     }
178   }
179   return *this;
180 }
181 //--------------------------------------------------------------------------
182 AliAODRecoDecay::~AliAODRecoDecay() {
183   //  
184   // Default Destructor
185   //
186   if(fPx) { delete [] fPx; fPx=NULL; } 
187   if(fPy) { delete [] fPy; fPy=NULL; }
188   if(fPz) { delete [] fPz; fPz=NULL; }
189   if(fd0) { delete [] fd0; fd0=NULL; }
190   if(fPID) { delete [] fPID; fPID=NULL; }
191   if(fDCA) { delete [] fDCA; fDCA=NULL; }
192   if(fOwnSecondaryVtx) { delete fOwnSecondaryVtx; fOwnSecondaryVtx=NULL; }
193 }
194 //----------------------------------------------------------------------------
195 Double_t AliAODRecoDecay::Alpha() const 
196 {
197   //
198   // Armenteros-Podolanski alpha for 2-prong decays
199   //
200   if(GetNProngs()!=2) {
201     printf("Can be called only for 2-prong decays");
202     return (Double_t)-99999.;
203   }
204   return 1.-2./(1.+QlProng(0)/QlProng(1));
205 }
206 //----------------------------------------------------------------------------
207 Double_t AliAODRecoDecay::DecayLength(Double_t point[3]) const 
208 {
209   //
210   // Decay length 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                     +(point[2]-GetSecVtxZ())
217                     *(point[2]-GetSecVtxZ()));  
218 }
219 //----------------------------------------------------------------------------
220 Double_t AliAODRecoDecay::DecayLengthXY(Double_t point[3]) const 
221 {
222   //
223   // Decay length in XY assuming it is produced at "point" [cm]
224   //
225   return TMath::Sqrt((point[0]-GetSecVtxX())
226                     *(point[0]-GetSecVtxX())
227                     +(point[1]-GetSecVtxY())
228                     *(point[1]-GetSecVtxY()));  
229 }
230 //----------------------------------------------------------------------------
231 Double_t AliAODRecoDecay::CosPointingAngle(Double_t point[3]) const 
232 {
233   //
234   // Cosine of pointing angle in space assuming it is produced at "point"
235   //
236   TVector3 mom(Px(),Py(),Pz());
237   TVector3 fline(GetSecVtxX()-point[0],
238                  GetSecVtxY()-point[1],
239                  GetSecVtxZ()-point[2]);
240
241   Double_t pta = mom.Angle(fline);
242
243   return TMath::Cos(pta); 
244 }
245 //----------------------------------------------------------------------------
246 Double_t AliAODRecoDecay::CosPointingAngleXY(Double_t point[3]) const 
247 {
248   //
249   // Cosine of pointing angle in transverse plane assuming it is produced 
250   // at "point"
251   //
252   TVector3 momXY(Px(),Py(),0.);
253   TVector3 flineXY(GetSecVtxX()-point[0],
254                    GetSecVtxY()-point[1],
255                    0.);
256
257   Double_t ptaXY = momXY.Angle(flineXY);
258
259   return TMath::Cos(ptaXY); 
260 }
261 //----------------------------------------------------------------------------
262 Double_t AliAODRecoDecay::CosThetaStar(Int_t ip,UInt_t pdgvtx,UInt_t pdgprong0,UInt_t pdgprong1) const 
263 {
264   //
265   // Only for 2-prong decays: 
266   // Cosine of decay angle (theta*) in the rest frame of the mother particle
267   // for prong ip (0 or 1) with mass hypotheses pdgvtx for mother particle,
268   // pdgprong0 for prong 0 and pdgprong1 for prong1
269   //
270   if(GetNProngs()!=2) {
271     printf("Can be called only for 2-prong decays");
272     return (Double_t)-99999.;
273   }
274   Double_t massvtx = TDatabasePDG::Instance()->GetParticle(pdgvtx)->Mass();
275   Double_t massp[2];
276   massp[0] = TDatabasePDG::Instance()->GetParticle(pdgprong0)->Mass();
277   massp[1] = TDatabasePDG::Instance()->GetParticle(pdgprong1)->Mass();
278
279   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);
280
281   Double_t beta = P()/E(pdgvtx);
282   Double_t gamma = E(pdgvtx)/massvtx;
283
284   Double_t cts = (QlProng(ip)/gamma-beta*TMath::Sqrt(pStar*pStar+massp[ip]*massp[ip]))/pStar;
285
286   return cts;
287 }
288 //---------------------------------------------------------------------------
289 Double_t AliAODRecoDecay::Ct(UInt_t pdg,Double_t point[3]) const
290 {
291   //
292   // Decay time * c assuming it is produced at "point" [cm]
293   //
294   Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
295   return DecayLength(point)*mass/P();
296 }
297 //---------------------------------------------------------------------------
298 Double_t AliAODRecoDecay::E(UInt_t pdg) const 
299 {
300   //
301   // Energy
302   //
303   Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
304   return TMath::Sqrt(mass*mass+P()*P());
305 }
306 //---------------------------------------------------------------------------
307 Double_t AliAODRecoDecay::EProng(Int_t ip,UInt_t pdg) const 
308 {
309   //
310   // Energy of ip-th prong 
311   //
312   Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
313   return TMath::Sqrt(mass*mass+PProng(ip)*PProng(ip));
314 }
315 //--------------------------------------------------------------------------
316 Bool_t AliAODRecoDecay::GetCovarianceXYZPxPyPz(Double_t cv[21]) const {
317   //
318   // This function returns the global covariance matrix of the track params
319   // 
320   // Cov(x,x) ... :   cv[0]
321   // Cov(y,x) ... :   cv[1]  cv[2]
322   // Cov(z,x) ... :   cv[3]  cv[4]  cv[5]
323   // Cov(px,x)... :   cv[6]  cv[7]  cv[8]  cv[9]
324   // Cov(py,x)... :   cv[10] cv[11] cv[12] cv[13] cv[14]
325   // Cov(pz,x)... :   cv[15] cv[16] cv[17] cv[18] cv[19] cv[20]
326   //
327   // For XYZ we take the cov of the vertex, for PxPyPz we take the 
328   // sum of the covs of PxPyPz from the daughters, for the moment 
329   // we set the cov between position and momentum as the sum of 
330   // the same cov from the daughters.
331   //
332
333   Int_t j;
334   for(j=0;j<21;j++) cv[j]=0.;
335
336   if(!GetNDaughters()) {
337     AliError("No daughters available");
338     return kFALSE;
339   }
340
341   Double_t v[6];
342   AliAODVertex *secv=GetSecondaryVtx();
343   if(!secv) {
344     AliError("Vertex covariance matrix not available");
345     return kFALSE;
346   }
347   if(!secv->GetCovMatrix(v)) {
348     AliError("Vertex covariance matrix not available");
349     return kFALSE;
350   }
351
352   Double_t p[21]; for(j=0;j<21;j++) p[j]=0.;
353   Bool_t error=kFALSE;
354   for(Int_t i=1; i<GetNDaughters(); i++) {
355     AliVTrack *daugh = (AliVTrack*)GetDaughter(i);
356     Double_t dcov[21];
357     if(!daugh->GetCovarianceXYZPxPyPz(dcov)) error=kTRUE;
358     for(j=0;j<21;j++) p[j] += dcov[j];
359   }
360   if(error) {
361     AliError("No covariance for at least one daughter")
362     return kFALSE;
363   }
364
365   for(j=0; j<21; j++) {
366     if(j<6) {
367       cv[j] = v[j];
368     } else {
369       cv[j] = p[j];
370     }
371   }
372
373   return kTRUE;
374 }
375 //----------------------------------------------------------------------------
376 UChar_t  AliAODRecoDecay::GetITSClusterMap() const {
377   //
378   // We take the logical AND of the daughters cluster maps 
379   // (only if all daughters have the bit for given layer, we set the bit)
380   //
381   UChar_t map=0;
382
383   if(!GetNDaughters()) {
384     AliError("No daughters available");
385     return map;
386   }
387
388   for(Int_t l=0; l<12; l++) { // loop on ITS layers (from here we cannot know how many they are; let's put 12 to be conservative)
389     Int_t bit = 1;
390     for(Int_t i=0; i<GetNDaughters(); i++) {
391       AliVTrack *daugh = (AliVTrack*)GetDaughter(i);
392       if(!TESTBIT(daugh->GetITSClusterMap(),l)) bit=0; 
393     }
394     if(bit) SETBIT(map,l);
395   }
396
397   return map;
398 }
399 //--------------------------------------------------------------------------
400 ULong_t AliAODRecoDecay::GetStatus() const {
401   // 
402   // Same as for ITSClusterMap
403   //
404   ULong_t status=0;
405
406   if(!GetNDaughters()) {
407     AliError("No daughters available");
408     return status;
409   }
410
411   AliVTrack *daugh0 = (AliVTrack*)GetDaughter(0);
412   status = status&(daugh0->GetStatus());
413
414   for(Int_t i=1; i<GetNDaughters(); i++) {
415     AliVTrack *daugh = (AliVTrack*)GetDaughter(i);
416     status = status&(daugh->GetStatus());
417   }
418
419   return status;
420 }
421 //--------------------------------------------------------------------------
422 Double_t AliAODRecoDecay::ImpParXY(Double_t point[3]) const 
423 {
424   //
425   // Impact parameter in the bending plane of the particle 
426   // w.r.t. to "point"
427   //
428   Double_t k = -(GetSecVtxX()-point[0])*Px()-(GetSecVtxY()-point[1])*Py();
429   k /= Pt()*Pt();
430   Double_t dx = GetSecVtxX()-point[0]+k*Px();
431   Double_t dy = GetSecVtxY()-point[1]+k*Py();
432   Double_t absImpPar = TMath::Sqrt(dx*dx+dy*dy);
433   TVector3 mom(Px(),Py(),Pz());
434   TVector3 fline(GetSecVtxX()-point[0],
435                  GetSecVtxY()-point[1],
436                  GetSecVtxZ()-point[2]);
437   TVector3 cross = mom.Cross(fline);
438   return (cross.Z()>0. ? absImpPar : -absImpPar);
439 }
440 //----------------------------------------------------------------------------
441 Double_t AliAODRecoDecay::InvMass(Int_t npdg,UInt_t *pdg) const 
442 {
443   //
444   // Invariant mass for prongs mass hypotheses in pdg array
445   //
446   if(GetNProngs()!=npdg) {
447     printf("npdg != GetNProngs()");
448     return (Double_t)-99999.;
449   }
450   Double_t energysum = 0.;
451
452   for(Int_t i=0; i<GetNProngs(); i++) {
453     energysum += EProng(i,pdg[i]);
454   }
455
456   Double_t mass = TMath::Sqrt(energysum*energysum-P()*P());
457
458   return mass;
459 }
460 //----------------------------------------------------------------------------
461 Double_t AliAODRecoDecay::InvMass2Prongs(Int_t ip1,Int_t ip2,
462                                       UInt_t pdg1,UInt_t pdg2) const 
463 {
464   //
465   // 2-prong(ip1,ip2) invariant mass for prongs mass hypotheses in pdg1,2
466   //
467   Double_t energysum = EProng(ip1,pdg1) + EProng(ip2,pdg2);
468   Double_t psum2 = (PxProng(ip1)+PxProng(ip2))*(PxProng(ip1)+PxProng(ip2))
469                   +(PyProng(ip1)+PyProng(ip2))*(PyProng(ip1)+PyProng(ip2))
470                   +(PzProng(ip1)+PzProng(ip2))*(PzProng(ip1)+PzProng(ip2));
471   Double_t mass = TMath::Sqrt(energysum*energysum-psum2);
472
473   return mass;
474 }
475 //---------------------------------------------------------------------------
476 void AliAODRecoDecay::Print(Option_t* /*option*/) const 
477 {
478   //
479   // Print some information
480   //
481   printf("AliAODRecoDecay with %d prongs\n",GetNProngs());
482   printf("Secondary Vertex: (%f, %f, %f)\n",GetSecVtxX(),GetSecVtxY(),GetSecVtxZ());
483
484   return;
485 }
486 //----------------------------------------------------------------------------
487 Double_t AliAODRecoDecay::ProngsRelAngle(Int_t ip1,Int_t ip2) const 
488 {
489   //
490   // Relative angle between two prongs
491   //
492   TVector3 momA(PxProng(ip1),PyProng(ip1),PzProng(ip1));
493   TVector3 momB(PxProng(ip2),PyProng(ip2),PzProng(ip2));
494
495   Double_t angle = momA.Angle(momB);
496
497   return angle; 
498 }
499 //----------------------------------------------------------------------------
500 Double_t AliAODRecoDecay::QlProng(Int_t ip) const 
501 {
502   //
503   // Longitudinal momentum of prong w.r.t. to total momentum
504   //
505   TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
506   TVector3 momTot(Px(),Py(),Pz());
507
508   return mom.Dot(momTot)/momTot.Mag();
509 }
510 //----------------------------------------------------------------------------
511 Double_t AliAODRecoDecay::QtProng(Int_t ip) const 
512 {
513   //
514   // Transverse momentum of prong w.r.t. to total momentum  
515   //
516   TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
517   TVector3 momTot(Px(),Py(),Pz());
518
519   return mom.Perp(momTot);
520 }
521 //----------------------------------------------------------------------------
522 Double_t AliAODRecoDecay::QlProngFlightLine(Int_t ip,Double_t point[3]) const 
523 {
524   //
525   // Longitudinal momentum of prong w.r.t. to flight line between "point"
526   // and fSecondaryVtx
527   //
528   TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
529   TVector3 fline(GetSecVtxX()-point[0],
530                  GetSecVtxY()-point[1],
531                  GetSecVtxZ()-point[2]);
532
533   return mom.Dot(fline)/fline.Mag();
534 }
535 //----------------------------------------------------------------------------
536 Double_t AliAODRecoDecay::QtProngFlightLine(Int_t ip,Double_t point[3]) const 
537 {
538   //
539   // Transverse momentum of prong w.r.t. to flight line between "point" and 
540   // fSecondaryVtx 
541   //
542   TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
543   TVector3 fline(GetSecVtxX()-point[0],
544                  GetSecVtxY()-point[1],
545                  GetSecVtxZ()-point[2]);
546
547   return mom.Perp(fline);
548 }
549 //--------------------------------------------------------------------------