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