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 "AliVParticle.h"
  26 #include "AliAODRecoDecay.h"
  27
  28 ClassImp(AliAODRecoDecay)
  29
  30 //--------------------------------------------------------------------------
  31 AliAODRecoDecay::AliAODRecoDecay() :
  32   AliVParticle(),
  33   fSecondaryVtx(0x0),
  34   fOwnSecondaryVtx(0x0),
  35   fCharge(0),
  36   fNProngs(0), fNDCA(0), fNPID(0),
  37   fPx(0x0), fPy(0x0), fPz(0x0),
  38   fd0(0x0),
  39   fDCA(0x0),
  40   fPID(0x0),
  41   fEventNumber(-1),fRunNumber(-1)
  42 {
  43   //
  44   // Default Constructor
  45   //
  46 }
  47 //--------------------------------------------------------------------------
  48 AliAODRecoDecay::AliAODRecoDecay(AliAODVertex *vtx2,Int_t nprongs,
  49                                  Short_t charge,
  50                                  Double_t *px,Double_t *py,Double_t *pz,
  51                                  Double_t *d0) :
  52   AliVParticle(),
  53   fSecondaryVtx(vtx2),
  54   fOwnSecondaryVtx(0x0),
  55   fCharge(charge),
  56   fNProngs(nprongs), fNDCA(0), fNPID(0),
  57   fPx(0x0), fPy(0x0), fPz(0x0),
  58   fd0(0x0),
  59   fDCA(0x0),
  60   fPID(0x0),
  61   fEventNumber(-1),fRunNumber(-1)
  62 {
  63   //
  64   // Constructor with AliAODVertex for decay vertex
  65   //
  66
  67   fPx = new Double_t[GetNProngs()];
  68   fPy = new Double_t[GetNProngs()];
  69   fPz = new Double_t[GetNProngs()];
  70   fd0 = new Double_t[GetNProngs()];
  71   for(Int_t i=0; i<GetNProngs(); i++) {
  72     fPx[i] = px[i];
  73     fPy[i] = py[i];
  74     fPz[i] = pz[i];
  75     fd0[i] = d0[i];
  76   }
  77 }
  78 //--------------------------------------------------------------------------
  79 AliAODRecoDecay::AliAODRecoDecay(AliAODVertex *vtx2,Int_t nprongs,
  80                                  Short_t charge,
  81                                  Double_t *d0) :
  82   AliVParticle(),
  83   fSecondaryVtx(vtx2),
  84   fOwnSecondaryVtx(0x0),
  85   fCharge(charge),
  86   fNProngs(nprongs), fNDCA(0), fNPID(0),
  87   fPx(0x0), fPy(0x0), fPz(0x0),
  88   fd0(0x0),
  89   fDCA(0x0),
  90   fPID(0x0),
  91   fEventNumber(-1),fRunNumber(-1)
  92 {
  93   //
  94   // Constructor with AliAODVertex for decay vertex and without prongs momenta
  95   //
  96
  97   fd0 = new Double_t[GetNProngs()];
  98   for(Int_t i=0; i<GetNProngs(); i++) fd0[i] = d0[i];
  99 }
 100 //--------------------------------------------------------------------------
 101 AliAODRecoDecay::AliAODRecoDecay(const AliAODRecoDecay &source) :
 102   AliVParticle(source),
 103   fSecondaryVtx(source.fSecondaryVtx),
 104   fOwnSecondaryVtx(source.fOwnSecondaryVtx),
 105   fCharge(source.fCharge),
 106   fNProngs(source.fNProngs), fNDCA(source.fNDCA), fNPID(source.fNPID),
 107   fPx(0x0), fPy(0x0), fPz(0x0),
 108   fd0(0x0),
 109   fDCA(0x0),
 110   fPID(0x0),
 111   fEventNumber(source.fEventNumber),fRunNumber(source.fRunNumber)
 112 {
 113   //
 114   // Copy constructor
 115   //
 116   if(source.GetNProngs()>0) {
 117     fd0 = new Double_t[GetNProngs()];
 118     memcpy(fd0,source.fd0,GetNProngs()*sizeof(Double_t));
 119     if(source.fPx) {
 120       fPx = new Double_t[GetNProngs()];
 121       fPy = new Double_t[GetNProngs()];
 122       fPz = new Double_t[GetNProngs()];
 123       memcpy(fPx,source.fPx,GetNProngs()*sizeof(Double_t));
 124       memcpy(fPy,source.fPy,GetNProngs()*sizeof(Double_t));
 125       memcpy(fPz,source.fPz,GetNProngs()*sizeof(Double_t));
 126     }
 127     if(source.fPID) {
 128       fPID = new Double_t[5*GetNProngs()];
 129       memcpy(fPID,source.fPID,GetNProngs()*sizeof(Double_t));
 130     }
 131     if(source.fDCA) {
 132       fDCA = new Double_t[GetNProngs()*(GetNProngs()-1)/2];
 133       memcpy(fDCA,source.fDCA,(GetNProngs()*(GetNProngs()-1)/2)*sizeof(Double_t));
 134     }
 135   }
 136 }
 137 //--------------------------------------------------------------------------
 138 AliAODRecoDecay &AliAODRecoDecay::operator=(const AliAODRecoDecay &source)
 139 {
 140   //
 141   // assignment operator
 142   //
 143   if(&source == this) return *this;
 144   fSecondaryVtx = source.fSecondaryVtx;
 145   fOwnSecondaryVtx = source.fOwnSecondaryVtx;
 146   fCharge = source.fCharge;
 147   fNProngs = source.fNProngs;
 148   fNDCA = source.fNDCA;
 149   fNPID = source.fNPID;
 150   fEventNumber = source.fEventNumber;
 151   fRunNumber = source.fRunNumber;
 152   if(source.GetNProngs()>0) {
 153     fd0 = new Double_t[GetNProngs()];
 154     memcpy(fd0,source.fd0,GetNProngs()*sizeof(Double_t));
 155     if(source.fPx) {
 156       fPx = new Double_t[GetNProngs()];
 157       fPy = new Double_t[GetNProngs()];
 158       fPz = new Double_t[GetNProngs()];
 159       memcpy(fPx,source.fPx,GetNProngs()*sizeof(Double_t));
 160       memcpy(fPy,source.fPy,GetNProngs()*sizeof(Double_t));
 161       memcpy(fPz,source.fPz,GetNProngs()*sizeof(Double_t));
 162     }
 163     if(source.fPID) {
 164       fPID = new Double_t[5*GetNProngs()];
 165       memcpy(fPID,source.fPID,GetNProngs()*sizeof(Double_t));
 166     }
 167     if(source.fDCA) {
 168       fDCA = new Double_t[GetNProngs()*(GetNProngs()-1)/2];
 169       memcpy(fDCA,source.fDCA,(GetNProngs()*(GetNProngs()-1)/2)*sizeof(Double32_t));
 170     }
 171   }
 172   return *this;
 173 }
 174 //--------------------------------------------------------------------------
 175 AliAODRecoDecay::~AliAODRecoDecay() {
 176   //
 177   // Default Destructor
 178   //
 179   if(fPx) { delete [] fPx; fPx=NULL; }
 180   if(fPy) { delete [] fPy; fPy=NULL; }
 181   if(fPz) { delete [] fPz; fPz=NULL; }
 182   if(fd0) { delete [] fd0; fd0=NULL; }
 183   if(fPID) { delete [] fPID; fPID=NULL; }
 184   if(fDCA) { delete [] fDCA; fDCA=NULL; }
 185   if(fOwnSecondaryVtx) { delete fOwnSecondaryVtx; fOwnSecondaryVtx=NULL; }
 186 }
 187 //----------------------------------------------------------------------------
 188 Double_t AliAODRecoDecay::Alpha() const
 189 {
 190   //
 191   // Armenteros-Podolanski alpha for 2-prong decays
 192   //
 193   if(GetNProngs()!=2) {
 194     printf("Can be called only for 2-prong decays");
 195     return (Double_t)-99999.;
 196   }
 197   return 1.-2./(1.+QlProng(0)/QlProng(1));
 198 }
 199 //----------------------------------------------------------------------------
 200 Double_t AliAODRecoDecay::DecayLength(Double_t point[3]) const
 201 {
 202   //
 203   // Decay length assuming it is produced at "point" [cm]
 204   //
 205   return TMath::Sqrt((point[0]-GetSecVtxX())
 206                     *(point[0]-GetSecVtxX())
 207                     +(point[1]-GetSecVtxY())
 208                     *(point[1]-GetSecVtxY())
 209                     +(point[2]-GetSecVtxZ())
 210                     *(point[2]-GetSecVtxZ()));
 211 }
 212 //----------------------------------------------------------------------------
 213 Double_t AliAODRecoDecay::DecayLengthXY(Double_t point[3]) const
 214 {
 215   //
 216   // Decay length in XY assuming it is produced at "point" [cm]
 217   //
 218   return TMath::Sqrt((point[0]-GetSecVtxX())
 219                     *(point[0]-GetSecVtxX())
 220                     +(point[1]-GetSecVtxY())
 221                     *(point[1]-GetSecVtxY()));
 222 }
 223 //----------------------------------------------------------------------------
 224 Double_t AliAODRecoDecay::CosPointingAngle(Double_t point[3]) const
 225 {
 226   //
 227   // Cosine of pointing angle in space assuming it is produced at "point"
 228   //
 229   TVector3 mom(Px(),Py(),Pz());
 230   TVector3 fline(GetSecVtxX()-point[0],
 231                  GetSecVtxY()-point[1],
 232                  GetSecVtxZ()-point[2]);
 233
 234   Double_t pta = mom.Angle(fline);
 235
 236   return TMath::Cos(pta);
 237 }
 238 //----------------------------------------------------------------------------
 239 Double_t AliAODRecoDecay::CosPointingAngleXY(Double_t point[3]) const
 240 {
 241   //
 242   // Cosine of pointing angle in transverse plane assuming it is produced
 243   // at "point"
 244   //
 245   TVector3 momXY(Px(),Py(),0.);
 246   TVector3 flineXY(GetSecVtxX()-point[0],
 247                    GetSecVtxY()-point[1],
 248                    0.);
 249
 250   Double_t ptaXY = momXY.Angle(flineXY);
 251
 252   return TMath::Cos(ptaXY);
 253 }
 254 //----------------------------------------------------------------------------
 255 Double_t AliAODRecoDecay::CosThetaStar(Int_t ip,UInt_t pdgvtx,UInt_t pdgprong0,UInt_t pdgprong1) const
 256 {
 257   //
 258   // Only for 2-prong decays:
 259   // Cosine of decay angle (theta*) in the rest frame of the mother particle
 260   // for prong ip (0 or 1) with mass hypotheses pdgvtx for mother particle,
 261   // pdgprong0 for prong 0 and pdgprong1 for prong1
 262   //
 263   if(GetNProngs()!=2) {
 264     printf("Can be called only for 2-prong decays");
 265     return (Double_t)-99999.;
 266   }
 267   Double_t massvtx = TDatabasePDG::Instance()->GetParticle(pdgvtx)->Mass();
 268   Double_t massp[2];
 269   massp[0] = TDatabasePDG::Instance()->GetParticle(pdgprong0)->Mass();
 270   massp[1] = TDatabasePDG::Instance()->GetParticle(pdgprong1)->Mass();
 271
 272   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);
 273
 274   Double_t beta = P()/E(pdgvtx);
 275   Double_t gamma = E(pdgvtx)/massvtx;
 276
 277   Double_t cts = (QlProng(ip)/gamma-beta*TMath::Sqrt(pStar*pStar+massp[ip]*massp[ip]))/pStar;
 278
 279   return cts;
 280 }
 281 //---------------------------------------------------------------------------
 282 Double_t AliAODRecoDecay::Ct(UInt_t pdg,Double_t point[3]) const
 283 {
 284   //
 285   // Decay time * c assuming it is produced at "point" [cm]
 286   //
 287   Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
 288   return DecayLength(point)*mass/P();
 289 }
 290 //---------------------------------------------------------------------------
 291 Double_t AliAODRecoDecay::E(UInt_t pdg) const
 292 {
 293   //
 294   // Energy
 295   //
 296   Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
 297   return TMath::Sqrt(mass*mass+P()*P());
 298 }
 299 //---------------------------------------------------------------------------
 300 Double_t AliAODRecoDecay::EProng(Int_t ip,UInt_t pdg) const
 301 {
 302   //
 303   // Energy of ip-th prong
 304   //
 305   Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
 306   return TMath::Sqrt(mass*mass+PProng(ip)*PProng(ip));
 307 }
 308 //---------------------------------------------------------------------------
 309 /*Int_t AliAODRecoDecay::GetIndexProng(Int_t ip) const
 310 {
 311   //
 312   // Index of prong ip
 313   //
 314   if(!GetNProngs()) return 999999;
 315 UShort_t *indices = GetSecondaryVtx()->GetIndices();
 316   return indices[ip];
 317 }*/
 318 //----------------------------------------------------------------------------
 319 Double_t AliAODRecoDecay::ImpParXY(Double_t point[3]) const
 320 {
 321   //
 322   // Impact parameter in the bending plane of the particle
 323   // w.r.t. to "point"
 324   //
 325   Double_t k = -(GetSecVtxX()-point[0])*Px()-(GetSecVtxY()-point[1])*Py();
 326   k /= Pt()*Pt();
 327   Double_t dx = GetSecVtxX()-point[0]+k*Px();
 328   Double_t dy = GetSecVtxY()-point[1]+k*Py();
 329   Double_t absImpPar = TMath::Sqrt(dx*dx+dy*dy);
 330   TVector3 mom(Px(),Py(),Pz());
 331   TVector3 fline(GetSecVtxX()-point[0],
 332                  GetSecVtxY()-point[1],
 333                  GetSecVtxZ()-point[2]);
 334   TVector3 cross = mom.Cross(fline);
 335   return (cross.Z()>0. ? absImpPar : -absImpPar);
 336 }
 337 //----------------------------------------------------------------------------
 338 Double_t AliAODRecoDecay::InvMass(Int_t npdg,UInt_t *pdg) const
 339 {
 340   //
 341   // Invariant mass for prongs mass hypotheses in pdg array
 342   //
 343   if(GetNProngs()!=npdg) {
 344     printf("npdg != GetNProngs()");
 345     return (Double_t)-99999.;
 346   }
 347   Double_t energysum = 0.;
 348
 349   for(Int_t i=0; i<GetNProngs(); i++) {
 350     energysum += EProng(i,pdg[i]);
 351   }
 352
 353   Double_t mass = TMath::Sqrt(energysum*energysum-P()*P());
 354
 355   return mass;
 356 }
 357 //----------------------------------------------------------------------------
 358 Double_t AliAODRecoDecay::InvMass2Prongs(Int_t ip1,Int_t ip2,
 359                                       UInt_t pdg1,UInt_t pdg2) const
 360 {
 361   //
 362   // 2-prong(ip1,ip2) invariant mass for prongs mass hypotheses in pdg1,2
 363   //
 364   Double_t energysum = EProng(ip1,pdg1) + EProng(ip2,pdg2);
 365   Double_t psum2 = (PxProng(ip1)+PxProng(ip2))*(PxProng(ip1)+PxProng(ip2))
 366                   +(PyProng(ip1)+PyProng(ip2))*(PyProng(ip1)+PyProng(ip2))
 367                   +(PzProng(ip1)+PzProng(ip2))*(PzProng(ip1)+PzProng(ip2));
 368   Double_t mass = TMath::Sqrt(energysum*energysum-psum2);
 369
 370   return mass;
 371 }
 372 //---------------------------------------------------------------------------
 373 void AliAODRecoDecay::Print(Option_t* /*option*/) const
 374 {
 375   //
 376   // Print some information
 377   //
 378   printf("AliAODRecoDecay with %d prongs\n",GetNProngs());
 379   printf("Secondary Vertex: (%f, %f, %f)\n",GetSecVtxX(),GetSecVtxY(),GetSecVtxZ());
 380
 381   return;
 382 }
 383 //----------------------------------------------------------------------------
 384 Double_t AliAODRecoDecay::ProngsRelAngle(Int_t ip1,Int_t ip2) const
 385 {
 386   //
 387   // Relative angle between two prongs
 388   //
 389   TVector3 momA(PxProng(ip1),PyProng(ip1),PzProng(ip1));
 390   TVector3 momB(PxProng(ip2),PyProng(ip2),PzProng(ip2));
 391
 392   Double_t angle = momA.Angle(momB);
 393
 394   return angle;
 395 }
 396 //----------------------------------------------------------------------------
 397 Double_t AliAODRecoDecay::QlProng(Int_t ip) const
 398 {
 399   //
 400   // Longitudinal momentum of prong w.r.t. to total momentum
 401   //
 402   TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
 403   TVector3 momTot(Px(),Py(),Pz());
 404
 405   return mom.Dot(momTot)/momTot.Mag();
 406 }
 407 //----------------------------------------------------------------------------
 408 Double_t AliAODRecoDecay::QtProng(Int_t ip) const
 409 {
 410   //
 411   // Transverse momentum of prong w.r.t. to total momentum
 412   //
 413   TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
 414   TVector3 momTot(Px(),Py(),Pz());
 415
 416   return mom.Perp(momTot);
 417 }
 418 //----------------------------------------------------------------------------
 419 Double_t AliAODRecoDecay::QlProngFlightLine(Int_t ip,Double_t point[3]) const
 420 {
 421   //
 422   // Longitudinal momentum of prong w.r.t. to flight line between "point"
 423   // and fSecondaryVtx
 424   //
 425   TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
 426   TVector3 fline(GetSecVtxX()-point[0],
 427                  GetSecVtxY()-point[1],
 428                  GetSecVtxZ()-point[2]);
 429
 430   return mom.Dot(fline)/fline.Mag();
 431 }
 432 //----------------------------------------------------------------------------
 433 Double_t AliAODRecoDecay::QtProngFlightLine(Int_t ip,Double_t point[3]) const
 434 {
 435   //
 436   // Transverse momentum of prong w.r.t. to flight line between "point" and
 437   // fSecondaryVtx
 438   //
 439   TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
 440   TVector3 fline(GetSecVtxX()-point[0],
 441                  GetSecVtxY()-point[1],
 442                  GetSecVtxZ()-point[2]);
 443
 444   return mom.Perp(fline);
 445 }
 446 //--------------------------------------------------------------------------