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 //--------------------------------------------------------------------------