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 "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::DecayLength(Double_t point[3]) const
 212 {
 213   //
 214   // Decay length assuming it is produced at "point" [cm]
 215   //
 216   return TMath::Sqrt((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 pta = mom.Angle(fline);
 246
 247   return TMath::Cos(pta);
 248 }
 249 //----------------------------------------------------------------------------
 250 Double_t AliAODRecoDecay::CosPointingAngleXY(Double_t point[3]) const
 251 {
 252   //
 253   // Cosine of pointing angle in transverse plane assuming it is produced
 254   // at "point"
 255   //
 256   TVector3 momXY(Px(),Py(),0.);
 257   TVector3 flineXY(GetSecVtxX()-point[0],
 258                    GetSecVtxY()-point[1],
 259                    0.);
 260
 261   Double_t ptaXY = momXY.Angle(flineXY);
 262
 263   return TMath::Cos(ptaXY);
 264 }
 265 //----------------------------------------------------------------------------
 266 Double_t AliAODRecoDecay::CosThetaStar(Int_t ip,UInt_t pdgvtx,UInt_t pdgprong0,UInt_t pdgprong1) const
 267 {
 268   //
 269   // Only for 2-prong decays:
 270   // Cosine of decay angle (theta*) in the rest frame of the mother particle
 271   // for prong ip (0 or 1) with mass hypotheses pdgvtx for mother particle,
 272   // pdgprong0 for prong 0 and pdgprong1 for prong1
 273   //
 274   if(GetNProngs()!=2) {
 275     printf("Can be called only for 2-prong decays");
 276     return (Double_t)-99999.;
 277   }
 278   Double_t massvtx = TDatabasePDG::Instance()->GetParticle(pdgvtx)->Mass();
 279   Double_t massp[2];
 280   massp[0] = TDatabasePDG::Instance()->GetParticle(pdgprong0)->Mass();
 281   massp[1] = TDatabasePDG::Instance()->GetParticle(pdgprong1)->Mass();
 282
 283   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);
 284
 285   Double_t beta = P()/E(pdgvtx);
 286   Double_t gamma = E(pdgvtx)/massvtx;
 287
 288   Double_t cts = (QlProng(ip)/gamma-beta*TMath::Sqrt(pStar*pStar+massp[ip]*massp[ip]))/pStar;
 289
 290   return cts;
 291 }
 292 //---------------------------------------------------------------------------
 293 Double_t AliAODRecoDecay::Ct(UInt_t pdg,Double_t point[3]) const
 294 {
 295   //
 296   // Decay time * c assuming it is produced at "point" [cm]
 297   //
 298   Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
 299   return DecayLength(point)*mass/P();
 300 }
 301 //---------------------------------------------------------------------------
 302 Double_t AliAODRecoDecay::E(UInt_t pdg) const
 303 {
 304   //
 305   // Energy
 306   //
 307   Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
 308   return TMath::Sqrt(mass*mass+P()*P());
 309 }
 310 //---------------------------------------------------------------------------
 311 Double_t AliAODRecoDecay::EProng(Int_t ip,UInt_t pdg) const
 312 {
 313   //
 314   // Energy of ip-th prong
 315   //
 316   Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
 317   return TMath::Sqrt(mass*mass+PProng(ip)*PProng(ip));
 318 }
 319 //--------------------------------------------------------------------------
 320 Bool_t AliAODRecoDecay::GetCovarianceXYZPxPyPz(Double_t cv[21]) const {
 321   //
 322   // This function returns the global covariance matrix of the track params
 323   //
 324   // Cov(x,x) ... :   cv[0]
 325   // Cov(y,x) ... :   cv[1]  cv[2]
 326   // Cov(z,x) ... :   cv[3]  cv[4]  cv[5]
 327   // Cov(px,x)... :   cv[6]  cv[7]  cv[8]  cv[9]
 328   // Cov(py,x)... :   cv[10] cv[11] cv[12] cv[13] cv[14]
 329   // Cov(pz,x)... :   cv[15] cv[16] cv[17] cv[18] cv[19] cv[20]
 330   //
 331   // For XYZ we take the cov of the vertex, for PxPyPz we take the
 332   // sum of the covs of PxPyPz from the daughters, for the moment
 333   // we set the cov between position and momentum as the sum of
 334   // the same cov from the daughters.
 335   //
 336
 337   Int_t j;
 338   for(j=0;j<21;j++) cv[j]=0.;
 339
 340   if(!GetNDaughters()) {
 341     AliError("No daughters available");
 342     return kFALSE;
 343   }
 344
 345   Double_t v[6];
 346   AliAODVertex *secv=GetSecondaryVtx();
 347   if(!secv) {
 348     AliError("Vertex covariance matrix not available");
 349     return kFALSE;
 350   }
 351   if(!secv->GetCovMatrix(v)) {
 352     AliError("Vertex covariance matrix not available");
 353     return kFALSE;
 354   }
 355
 356   Double_t p[21]; for(j=0;j<21;j++) p[j]=0.;
 357   Bool_t error=kFALSE;
 358   for(Int_t i=1; i<GetNDaughters(); i++) {
 359     AliVTrack *daugh = (AliVTrack*)GetDaughter(i);
 360     Double_t dcov[21];
 361     if(!daugh->GetCovarianceXYZPxPyPz(dcov)) error=kTRUE;
 362     for(j=0;j<21;j++) p[j] += dcov[j];
 363   }
 364   if(error) {
 365     AliError("No covariance for at least one daughter")
 366     return kFALSE;
 367   }
 368
 369   for(j=0; j<21; j++) {
 370     if(j<6) {
 371       cv[j] = v[j];
 372     } else {
 373       cv[j] = p[j];
 374     }
 375   }
 376
 377   return kTRUE;
 378 }
 379 //----------------------------------------------------------------------------
 380 UChar_t  AliAODRecoDecay::GetITSClusterMap() const {
 381   //
 382   // We take the logical AND of the daughters cluster maps
 383   // (only if all daughters have the bit for given layer, we set the bit)
 384   //
 385   UChar_t map=0;
 386
 387   if(!GetNDaughters()) {
 388     AliError("No daughters available");
 389     return map;
 390   }
 391
 392   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)
 393     Int_t bit = 1;
 394     for(Int_t i=0; i<GetNDaughters(); i++) {
 395       AliVTrack *daugh = (AliVTrack*)GetDaughter(i);
 396       if(!TESTBIT(daugh->GetITSClusterMap(),l)) bit=0;
 397     }
 398     if(bit) SETBIT(map,l);
 399   }
 400
 401   return map;
 402 }
 403 //--------------------------------------------------------------------------
 404 ULong_t AliAODRecoDecay::GetStatus() const {
 405   //
 406   // Same as for ITSClusterMap
 407   //
 408   ULong_t status=0;
 409
 410   if(!GetNDaughters()) {
 411     AliError("No daughters available");
 412     return status;
 413   }
 414
 415   AliVTrack *daugh0 = (AliVTrack*)GetDaughter(0);
 416   status = status&(daugh0->GetStatus());
 417
 418   for(Int_t i=1; i<GetNDaughters(); i++) {
 419     AliVTrack *daugh = (AliVTrack*)GetDaughter(i);
 420     status = status&(daugh->GetStatus());
 421   }
 422
 423   return status;
 424 }
 425 //--------------------------------------------------------------------------
 426 Bool_t AliAODRecoDecay::PropagateToDCA(const AliVVertex* vtx,Double_t b,Double_t maxd,Double_t dz[2],Double_t covar[3])
 427 {
 428   // compute impact parameters to the vertex vtx and their covariance matrix
 429   // b is the Bz, needed to propagate correctly the track to vertex
 430   // return kFALSE is something went wrong
 431
 432   AliWarning("The AliAODRecoDecay momentum is not updated to the DCA");
 433
 434   Bool_t retval=kTRUE;
 435   if(Charge()==0) {
 436     // convert to AliNeutralTrackParam
 437     AliNeutralTrackParam ntp(this);
 438     retval = ntp.PropagateToDCA(vtx,b,maxd,dz,covar);
 439   } else {
 440     // convert to AliExternalTrackParam
 441     AliExternalTrackParam etp(this);
 442     retval = etp.PropagateToDCA(vtx,b,maxd,dz,covar);
 443   }
 444   return retval;
 445 }
 446 //--------------------------------------------------------------------------
 447 Double_t AliAODRecoDecay::ImpParXY(Double_t point[3]) const
 448 {
 449   //
 450   // Impact parameter in the bending plane of the particle
 451   // w.r.t. to "point"
 452   //
 453   Double_t k = -(GetSecVtxX()-point[0])*Px()-(GetSecVtxY()-point[1])*Py();
 454   k /= Pt()*Pt();
 455   Double_t dx = GetSecVtxX()-point[0]+k*Px();
 456   Double_t dy = GetSecVtxY()-point[1]+k*Py();
 457   Double_t absImpPar = TMath::Sqrt(dx*dx+dy*dy);
 458   TVector3 mom(Px(),Py(),Pz());
 459   TVector3 fline(GetSecVtxX()-point[0],
 460                  GetSecVtxY()-point[1],
 461                  GetSecVtxZ()-point[2]);
 462   TVector3 cross = mom.Cross(fline);
 463   return (cross.Z()>0. ? absImpPar : -absImpPar);
 464 }
 465 //----------------------------------------------------------------------------
 466 Double_t AliAODRecoDecay::InvMass(Int_t npdg,UInt_t *pdg) const
 467 {
 468   //
 469   // Invariant mass for prongs mass hypotheses in pdg array
 470   //
 471   if(GetNProngs()!=npdg) {
 472     printf("npdg != GetNProngs()");
 473     return (Double_t)-99999.;
 474   }
 475   Double_t energysum = 0.;
 476
 477   for(Int_t i=0; i<GetNProngs(); i++) {
 478     energysum += EProng(i,pdg[i]);
 479   }
 480
 481   Double_t mass = TMath::Sqrt(energysum*energysum-P()*P());
 482
 483   return mass;
 484 }
 485 //----------------------------------------------------------------------------
 486 Double_t AliAODRecoDecay::InvMass2Prongs(Int_t ip1,Int_t ip2,
 487                                       UInt_t pdg1,UInt_t pdg2) const
 488 {
 489   //
 490   // 2-prong(ip1,ip2) invariant mass for prongs mass hypotheses in pdg1,2
 491   //
 492   Double_t energysum = EProng(ip1,pdg1) + EProng(ip2,pdg2);
 493   Double_t psum2 = (PxProng(ip1)+PxProng(ip2))*(PxProng(ip1)+PxProng(ip2))
 494                   +(PyProng(ip1)+PyProng(ip2))*(PyProng(ip1)+PyProng(ip2))
 495                   +(PzProng(ip1)+PzProng(ip2))*(PzProng(ip1)+PzProng(ip2));
 496   Double_t mass = TMath::Sqrt(energysum*energysum-psum2);
 497
 498   return mass;
 499 }
 500 //----------------------------------------------------------------------------
 501 Int_t AliAODRecoDecay::MatchToMC(Int_t pdgabs,TClonesArray *mcArray,
 502                                  Int_t ndgCk,Int_t *pdgDg) const
 503 {
 504   //
 505   // Check if this candidate is matched to a MC signal
 506   // If no, return -1
 507   // If yes, return label (>=0) of the AliAODMCParticle
 508   //
 509   // if ndgCk>0, checks also daughters PDGs
 510   //
 511   Int_t ndg=GetNDaughters();
 512   if(!ndg) {
 513     AliError("No daughters available");
 514     return -1;
 515   }
 516   if(ndg>10) {
 517     AliError("Only decays with <10 daughters supported");
 518     return -1;
 519   }
 520   if(ndgCk>0 && ndgCk!=ndg) {
 521     AliError("Wrong number of daughter PDGs passed");
 522     return -1;
 523   }
 524
 525   Int_t dgLabels[10] = {0};
 526
 527   // loop on daughters and write the labels
 528   for(Int_t i=0; i<ndg; i++) {
 529     AliAODTrack *trk = (AliAODTrack*)GetDaughter(i);
 530     dgLabels[i] = trk->GetLabel();
 531   }
 532
 533   return MatchToMC(pdgabs,mcArray,dgLabels,ndg,ndgCk,pdgDg);
 534 }
 535 //----------------------------------------------------------------------------
 536 Int_t AliAODRecoDecay::MatchToMC(Int_t pdgabs,TClonesArray *mcArray,
 537                                  Int_t dgLabels[10],Int_t ndg,
 538                                  Int_t ndgCk,Int_t *pdgDg) const
 539 {
 540   //
 541   // Check if this candidate is matched to a MC signal
 542   // If no, return -1
 543   // If yes, return label (>=0) of the AliAODMCParticle
 544   //
 545
 546   Int_t labMom[10]={0,0,0,0,0,0,0,0,0,0};
 547   Int_t i,j,lab,labMother,pdgMother,pdgPart;
 548   AliAODMCParticle *part=0;
 549   AliAODMCParticle *mother=0;
 550   Double_t pxSumDgs=0.,pySumDgs=0.,pzSumDgs=0.;
 551   Bool_t pdgUsed[10]={kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE};
 552
 553   // loop on daughter labels
 554   for(i=0; i<ndg; i++) {
 555     labMom[i]=-1;
 556     lab = dgLabels[i];
 557     if(lab<0) {
 558       printf("daughter with negative label %d\n",lab);
 559       return -1;
 560     }
 561     part = (AliAODMCParticle*)mcArray->At(lab);
 562     if(!part) {
 563       printf("no MC particle\n");
 564       return -1;
 565     }
 566
 567     // check the PDG of the daughter, if requested
 568     if(ndgCk>0) {
 569       pdgPart=TMath::Abs(part->GetPdgCode());
 570       for(j=0; j<ndg; j++) {
 571         if(!pdgUsed[j] && pdgPart==pdgDg[j]) {
 572           pdgUsed[j]=kTRUE;
 573           break;
 574         }
 575       }
 576     }
 577
 578     // for the J/psi, check that the daughters are electrons
 579     if(pdgabs==443) {
 580       if(TMath::Abs(part->GetPdgCode())!=11) return -1;
 581     }
 582
 583     mother = part;
 584     while(mother->GetMother()>=0) {
 585       labMother=mother->GetMother();
 586       mother = (AliAODMCParticle*)mcArray->At(labMother);
 587       if(!mother) {
 588         printf("no MC mother particle\n");
 589         break;
 590       }
 591       pdgMother = TMath::Abs(mother->GetPdgCode());
 592       if(pdgMother==pdgabs) {
 593         labMom[i]=labMother;
 594         // keep sum of daughters' momenta, to check for mom conservation
 595         pxSumDgs += part->Px();
 596         pySumDgs += part->Py();
 597         pzSumDgs += part->Pz();
 598         break;
 599       } else if(pdgMother>pdgabs || pdgMother<10) {
 600         break;
 601       }
 602     }
 603     if(labMom[i]==-1) return -1; // mother PDG not ok for this daughter
 604   } // end loop on daughters
 605
 606   // check if the candidate is signal
 607   labMother=labMom[0];
 608   // all labels have to be the same and !=-1
 609   for(i=0; i<ndg; i++) {
 610     if(labMom[i]==-1)        return -1;
 611     if(labMom[i]!=labMother) return -1;
 612   }
 613
 614   // check that all daughter PDGs are matched
 615   if(ndgCk>0) {
 616     for(i=0; i<ndg; i++) {
 617       if(pdgUsed[i]==kFALSE) return -1;
 618     }
 619   }
 620
 621   /*
 622   // check that the mother decayed in <GetNDaughters()> prongs
 623   Int_t ndg2 = TMath::Abs(mother->GetDaughter(1)-mother->GetDaughter(0))+1;
 624   printf("  MC daughters %d\n",ndg2);
 625   //if(ndg!=GetNDaughters()) return -1;
 626   AliAODMCParticle* p1=(AliAODMCParticle*)(mcArray->At(mother->GetDaughter(1)));
 627   AliAODMCParticle* p0=(AliAODMCParticle*)(mcArray->At(mother->GetDaughter(0)));
 628   printf(">>>>>>>> pdg %d  %d %d %d   %d %d\n",mother->GetDaughter(0),mother->GetDaughter(1),dgLabels[0],dgLabels[1],p0->GetPdgCode(),p1->GetPdgCode());
 629   */
 630
 631   // the above works only for non-resonant decays,
 632   // it's better to check for mom conservation
 633   mother = (AliAODMCParticle*)mcArray->At(labMother);
 634   Double_t pxMother = mother->Px();
 635   Double_t pyMother = mother->Py();
 636   Double_t pzMother = mother->Pz();
 637   // within 0.1%
 638   if((TMath::Abs(pxMother-pxSumDgs)/(TMath::Abs(pxMother)+1.e-13)) > 0.00001 &&
 639      (TMath::Abs(pyMother-pySumDgs)/(TMath::Abs(pyMother)+1.e-13)) > 0.00001 &&
 640      (TMath::Abs(pzMother-pzSumDgs)/(TMath::Abs(pzMother)+1.e-13)) > 0.00001)
 641     return -1;
 642
 643   return labMother;
 644 }
 645 //---------------------------------------------------------------------------
 646 void AliAODRecoDecay::Print(Option_t* /*option*/) const
 647 {
 648   //
 649   // Print some information
 650   //
 651   printf("AliAODRecoDecay with %d prongs\n",GetNProngs());
 652   printf("Secondary Vertex: (%f, %f, %f)\n",GetSecVtxX(),GetSecVtxY(),GetSecVtxZ());
 653
 654   return;
 655 }
 656 //----------------------------------------------------------------------------
 657 Double_t AliAODRecoDecay::ProngsRelAngle(Int_t ip1,Int_t ip2) const
 658 {
 659   //
 660   // Relative angle between two prongs
 661   //
 662   TVector3 momA(PxProng(ip1),PyProng(ip1),PzProng(ip1));
 663   TVector3 momB(PxProng(ip2),PyProng(ip2),PzProng(ip2));
 664
 665   Double_t angle = momA.Angle(momB);
 666
 667   return angle;
 668 }
 669 //----------------------------------------------------------------------------
 670 Double_t AliAODRecoDecay::QlProng(Int_t ip) const
 671 {
 672   //
 673   // Longitudinal momentum of prong w.r.t. to total momentum
 674   //
 675   TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
 676   TVector3 momTot(Px(),Py(),Pz());
 677
 678   return mom.Dot(momTot)/momTot.Mag();
 679 }
 680 //----------------------------------------------------------------------------
 681 Double_t AliAODRecoDecay::QtProng(Int_t ip) const
 682 {
 683   //
 684   // Transverse momentum of prong w.r.t. to total momentum
 685   //
 686   TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
 687   TVector3 momTot(Px(),Py(),Pz());
 688
 689   return mom.Perp(momTot);
 690 }
 691 //----------------------------------------------------------------------------
 692 Double_t AliAODRecoDecay::QlProngFlightLine(Int_t ip,Double_t point[3]) const
 693 {
 694   //
 695   // Longitudinal momentum of prong w.r.t. to flight line between "point"
 696   // and fSecondaryVtx
 697   //
 698   TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
 699   TVector3 fline(GetSecVtxX()-point[0],
 700                  GetSecVtxY()-point[1],
 701                  GetSecVtxZ()-point[2]);
 702
 703   return mom.Dot(fline)/fline.Mag();
 704 }
 705 //----------------------------------------------------------------------------
 706 Double_t AliAODRecoDecay::QtProngFlightLine(Int_t ip,Double_t point[3]) const
 707 {
 708   //
 709   // Transverse momentum of prong w.r.t. to flight line between "point" and
 710   // fSecondaryVtx
 711   //
 712   TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
 713   TVector3 fline(GetSecVtxX()-point[0],
 714                  GetSecVtxY()-point[1],
 715                  GetSecVtxZ()-point[2]);
 716
 717   return mom.Perp(fline);
 718 }
 719 //--------------------------------------------------------------------------