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,
 479                                  Int_t ndgCk,Int_t *pdgDg) const
 480 {
 481   //
 482   // Check if this candidate is matched to a MC signal
 483   // If no, return -1
 484   // If yes, return label (>=0) of the AliAODMCParticle
 485   //
 486   // if ndgCk>0, checks also daughters PDGs
 487   //
 488   Int_t ndg=GetNDaughters();
 489   if(!ndg) {
 490     AliError("No daughters available");
 491     return -1;
 492   }
 493   if(ndg>10) {
 494     AliError("Only decays with <10 daughters supported");
 495     return -1;
 496   }
 497   if(ndgCk>0 && ndgCk!=ndg) {
 498     AliError("Wrong number of daughter PDGs passed");
 499     return -1;
 500   }
 501
 502   Int_t dgLabels[10];
 503
 504   // loop on daughters and write the labels
 505   for(Int_t i=0; i<ndg; i++) {
 506     AliAODTrack *trk = (AliAODTrack*)GetDaughter(i);
 507     dgLabels[i] = trk->GetLabel();
 508   }
 509
 510   return MatchToMC(pdgabs,mcArray,dgLabels,ndg,ndgCk,pdgDg);
 511 }
 512 //----------------------------------------------------------------------------
 513 Int_t AliAODRecoDecay::MatchToMC(Int_t pdgabs,TClonesArray *mcArray,
 514                                  Int_t dgLabels[10],Int_t ndg,
 515                                  Int_t ndgCk,Int_t *pdgDg) const
 516 {
 517   //
 518   // Check if this candidate is matched to a MC signal
 519   // If no, return -1
 520   // If yes, return label (>=0) of the AliAODMCParticle
 521   //
 522
 523   Int_t labMom[10]={0,0,0,0,0,0,0,0,0,0};
 524   Int_t i,j,lab,labMother,pdgMother,pdgPart;
 525   AliAODMCParticle *part=0;
 526   AliAODMCParticle *mother=0;
 527   Double_t pxSumDgs=0.,pySumDgs=0.,pzSumDgs=0.;
 528   Bool_t pdgUsed[10]={kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE};
 529
 530   // loop on daughter labels
 531   for(i=0; i<ndg; i++) {
 532     labMom[i]=-1;
 533     lab = dgLabels[i];
 534     if(lab<0) {
 535       printf("daughter with negative label %d\n",lab);
 536       return -1;
 537     }
 538     part = (AliAODMCParticle*)mcArray->At(lab);
 539     if(!part) {
 540       printf("no MC particle\n");
 541       return -1;
 542     }
 543
 544     // check the PDG of the daughter, if requested
 545     if(ndgCk>0) {
 546       pdgPart=TMath::Abs(part->GetPdgCode());
 547       for(j=0; j<ndg; j++) {
 548         if(!pdgUsed[j] && pdgPart==pdgDg[j]) {
 549           pdgUsed[j]=kTRUE;
 550           break;
 551         }
 552       }
 553     }
 554
 555     // for the J/psi, check that the daughters are electrons
 556     if(pdgabs==443) {
 557       if(TMath::Abs(part->GetPdgCode())!=11) return -1;
 558     }
 559
 560     mother = part;
 561     while(mother->GetMother()>=0) {
 562       labMother=mother->GetMother();
 563       mother = (AliAODMCParticle*)mcArray->At(labMother);
 564       if(!mother) {
 565         printf("no MC mother particle\n");
 566         break;
 567       }
 568       pdgMother = TMath::Abs(mother->GetPdgCode());
 569       if(pdgMother==pdgabs) {
 570         labMom[i]=labMother;
 571         // keep sum of daughters' momenta, to check for mom conservation
 572         pxSumDgs += part->Px();
 573         pySumDgs += part->Py();
 574         pzSumDgs += part->Pz();
 575         break;
 576       } else if(pdgMother>pdgabs || pdgMother<10) {
 577         break;
 578       }
 579     }
 580     if(labMom[i]==-1) return -1; // mother PDG not ok for this daughter
 581   } // end loop on daughters
 582
 583   // check if the candidate is signal
 584   labMother=labMom[0];
 585   // all labels have to be the same and !=-1
 586   for(i=0; i<ndg; i++) {
 587     if(labMom[i]==-1)        return -1;
 588     if(labMom[i]!=labMother) return -1;
 589   }
 590
 591   // check that all daughter PDGs are matched
 592   if(ndgCk>0) {
 593     for(i=0; i<ndg; i++) {
 594       if(pdgUsed[i]==kFALSE) return -1;
 595     }
 596   }
 597
 598   /*
 599   // check that the mother decayed in <GetNDaughters()> prongs
 600   Int_t ndg2 = TMath::Abs(mother->GetDaughter(1)-mother->GetDaughter(0))+1;
 601   printf("  MC daughters %d\n",ndg2);
 602   //if(ndg!=GetNDaughters()) return -1;
 603   AliAODMCParticle* p1=(AliAODMCParticle*)(mcArray->At(mother->GetDaughter(1)));
 604   AliAODMCParticle* p0=(AliAODMCParticle*)(mcArray->At(mother->GetDaughter(0)));
 605   printf(">>>>>>>> pdg %d  %d %d %d   %d %d\n",mother->GetDaughter(0),mother->GetDaughter(1),dgLabels[0],dgLabels[1],p0->GetPdgCode(),p1->GetPdgCode());
 606   */
 607
 608   // the above works only for non-resonant decays,
 609   // it's better to check for mom conservation
 610   mother = (AliAODMCParticle*)mcArray->At(labMother);
 611   Double_t pxMother = mother->Px();
 612   Double_t pyMother = mother->Py();
 613   Double_t pzMother = mother->Pz();
 614   // within 0.1%
 615   if((TMath::Abs(pxMother-pxSumDgs)/(TMath::Abs(pxMother)+1.e-13)) > 0.00001 &&
 616      (TMath::Abs(pyMother-pySumDgs)/(TMath::Abs(pyMother)+1.e-13)) > 0.00001 &&
 617      (TMath::Abs(pzMother-pzSumDgs)/(TMath::Abs(pzMother)+1.e-13)) > 0.00001)
 618     return -1;
 619
 620   return labMother;
 621 }
 622 //---------------------------------------------------------------------------
 623 void AliAODRecoDecay::Print(Option_t* /*option*/) const
 624 {
 625   //
 626   // Print some information
 627   //
 628   printf("AliAODRecoDecay with %d prongs\n",GetNProngs());
 629   printf("Secondary Vertex: (%f, %f, %f)\n",GetSecVtxX(),GetSecVtxY(),GetSecVtxZ());
 630
 631   return;
 632 }
 633 //----------------------------------------------------------------------------
 634 Double_t AliAODRecoDecay::ProngsRelAngle(Int_t ip1,Int_t ip2) const
 635 {
 636   //
 637   // Relative angle between two prongs
 638   //
 639   TVector3 momA(PxProng(ip1),PyProng(ip1),PzProng(ip1));
 640   TVector3 momB(PxProng(ip2),PyProng(ip2),PzProng(ip2));
 641
 642   Double_t angle = momA.Angle(momB);
 643
 644   return angle;
 645 }
 646 //----------------------------------------------------------------------------
 647 Double_t AliAODRecoDecay::QlProng(Int_t ip) const
 648 {
 649   //
 650   // Longitudinal momentum of prong w.r.t. to total momentum
 651   //
 652   TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
 653   TVector3 momTot(Px(),Py(),Pz());
 654
 655   return mom.Dot(momTot)/momTot.Mag();
 656 }
 657 //----------------------------------------------------------------------------
 658 Double_t AliAODRecoDecay::QtProng(Int_t ip) const
 659 {
 660   //
 661   // Transverse momentum of prong w.r.t. to total momentum
 662   //
 663   TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
 664   TVector3 momTot(Px(),Py(),Pz());
 665
 666   return mom.Perp(momTot);
 667 }
 668 //----------------------------------------------------------------------------
 669 Double_t AliAODRecoDecay::QlProngFlightLine(Int_t ip,Double_t point[3]) const
 670 {
 671   //
 672   // Longitudinal momentum of prong w.r.t. to flight line between "point"
 673   // and fSecondaryVtx
 674   //
 675   TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
 676   TVector3 fline(GetSecVtxX()-point[0],
 677                  GetSecVtxY()-point[1],
 678                  GetSecVtxZ()-point[2]);
 679
 680   return mom.Dot(fline)/fline.Mag();
 681 }
 682 //----------------------------------------------------------------------------
 683 Double_t AliAODRecoDecay::QtProngFlightLine(Int_t ip,Double_t point[3]) const
 684 {
 685   //
 686   // Transverse momentum of prong w.r.t. to flight line between "point" and
 687   // fSecondaryVtx
 688   //
 689   TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
 690   TVector3 fline(GetSecVtxX()-point[0],
 691                  GetSecVtxY()-point[1],
 692                  GetSecVtxZ()-point[2]);
 693
 694   return mom.Perp(fline);
 695 }
 696 //--------------------------------------------------------------------------