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