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