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