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     if(fd0)delete [] fd0;
 154     fd0 = new Double_t[GetNProngs()];
 155     memcpy(fd0,source.fd0,GetNProngs()*sizeof(Double_t));
 156     if(source.fPx) {
 157       if(fPx) delete [] fPx;
 158       fPx = new Double_t[GetNProngs()];
 159       if(fPy) delete [] fPy;
 160       fPy = new Double_t[GetNProngs()];
 161       if(fPz) delete [] fPz;
 162       fPz = new Double_t[GetNProngs()];
 163       memcpy(fPx,source.fPx,GetNProngs()*sizeof(Double_t));
 164       memcpy(fPy,source.fPy,GetNProngs()*sizeof(Double_t));
 165       memcpy(fPz,source.fPz,GetNProngs()*sizeof(Double_t));
 166     }
 167     if(source.fPID) {
 168       if(fPID) delete [] fPID;
 169       fPID = new Double_t[5*GetNProngs()];
 170       memcpy(fPID,source.fPID,GetNProngs()*sizeof(Double_t));
 171     }
 172     if(source.fDCA) {
 173       if(fDCA) delete [] fDCA;
 174       fDCA = new Double_t[GetNProngs()*(GetNProngs()-1)/2];
 175       memcpy(fDCA,source.fDCA,(GetNProngs()*(GetNProngs()-1)/2)*sizeof(Double32_t));
 176     }
 177   }
 178   return *this;
 179 }
 180 //--------------------------------------------------------------------------
 181 AliAODRecoDecay::~AliAODRecoDecay() {
 182   //
 183   // Default Destructor
 184   //
 185   if(fPx) { delete [] fPx; fPx=NULL; }
 186   if(fPy) { delete [] fPy; fPy=NULL; }
 187   if(fPz) { delete [] fPz; fPz=NULL; }
 188   if(fd0) { delete [] fd0; fd0=NULL; }
 189   if(fPID) { delete [] fPID; fPID=NULL; }
 190   if(fDCA) { delete [] fDCA; fDCA=NULL; }
 191   if(fOwnSecondaryVtx) { delete fOwnSecondaryVtx; fOwnSecondaryVtx=NULL; }
 192 }
 193 //----------------------------------------------------------------------------
 194 Double_t AliAODRecoDecay::Alpha() const
 195 {
 196   //
 197   // Armenteros-Podolanski alpha for 2-prong decays
 198   //
 199   if(GetNProngs()!=2) {
 200     printf("Can be called only for 2-prong decays");
 201     return (Double_t)-99999.;
 202   }
 203   return 1.-2./(1.+QlProng(0)/QlProng(1));
 204 }
 205 //----------------------------------------------------------------------------
 206 Double_t AliAODRecoDecay::DecayLength(Double_t point[3]) const
 207 {
 208   //
 209   // Decay length assuming it is produced at "point" [cm]
 210   //
 211   return TMath::Sqrt((point[0]-GetSecVtxX())
 212                     *(point[0]-GetSecVtxX())
 213                     +(point[1]-GetSecVtxY())
 214                     *(point[1]-GetSecVtxY())
 215                     +(point[2]-GetSecVtxZ())
 216                     *(point[2]-GetSecVtxZ()));
 217 }
 218 //----------------------------------------------------------------------------
 219 Double_t AliAODRecoDecay::DecayLengthXY(Double_t point[3]) const
 220 {
 221   //
 222   // Decay length in XY assuming it is produced at "point" [cm]
 223   //
 224   return TMath::Sqrt((point[0]-GetSecVtxX())
 225                     *(point[0]-GetSecVtxX())
 226                     +(point[1]-GetSecVtxY())
 227                     *(point[1]-GetSecVtxY()));
 228 }
 229 //----------------------------------------------------------------------------
 230 Double_t AliAODRecoDecay::CosPointingAngle(Double_t point[3]) const
 231 {
 232   //
 233   // Cosine of pointing angle in space assuming it is produced at "point"
 234   //
 235   TVector3 mom(Px(),Py(),Pz());
 236   TVector3 fline(GetSecVtxX()-point[0],
 237                  GetSecVtxY()-point[1],
 238                  GetSecVtxZ()-point[2]);
 239
 240   Double_t pta = mom.Angle(fline);
 241
 242   return TMath::Cos(pta);
 243 }
 244 //----------------------------------------------------------------------------
 245 Double_t AliAODRecoDecay::CosPointingAngleXY(Double_t point[3]) const
 246 {
 247   //
 248   // Cosine of pointing angle in transverse plane assuming it is produced
 249   // at "point"
 250   //
 251   TVector3 momXY(Px(),Py(),0.);
 252   TVector3 flineXY(GetSecVtxX()-point[0],
 253                    GetSecVtxY()-point[1],
 254                    0.);
 255
 256   Double_t ptaXY = momXY.Angle(flineXY);
 257
 258   return TMath::Cos(ptaXY);
 259 }
 260 //----------------------------------------------------------------------------
 261 Double_t AliAODRecoDecay::CosThetaStar(Int_t ip,UInt_t pdgvtx,UInt_t pdgprong0,UInt_t pdgprong1) const
 262 {
 263   //
 264   // Only for 2-prong decays:
 265   // Cosine of decay angle (theta*) in the rest frame of the mother particle
 266   // for prong ip (0 or 1) with mass hypotheses pdgvtx for mother particle,
 267   // pdgprong0 for prong 0 and pdgprong1 for prong1
 268   //
 269   if(GetNProngs()!=2) {
 270     printf("Can be called only for 2-prong decays");
 271     return (Double_t)-99999.;
 272   }
 273   Double_t massvtx = TDatabasePDG::Instance()->GetParticle(pdgvtx)->Mass();
 274   Double_t massp[2];
 275   massp[0] = TDatabasePDG::Instance()->GetParticle(pdgprong0)->Mass();
 276   massp[1] = TDatabasePDG::Instance()->GetParticle(pdgprong1)->Mass();
 277
 278   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);
 279
 280   Double_t beta = P()/E(pdgvtx);
 281   Double_t gamma = E(pdgvtx)/massvtx;
 282
 283   Double_t cts = (QlProng(ip)/gamma-beta*TMath::Sqrt(pStar*pStar+massp[ip]*massp[ip]))/pStar;
 284
 285   return cts;
 286 }
 287 //---------------------------------------------------------------------------
 288 Double_t AliAODRecoDecay::Ct(UInt_t pdg,Double_t point[3]) const
 289 {
 290   //
 291   // Decay time * c assuming it is produced at "point" [cm]
 292   //
 293   Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
 294   return DecayLength(point)*mass/P();
 295 }
 296 //---------------------------------------------------------------------------
 297 Double_t AliAODRecoDecay::E(UInt_t pdg) const
 298 {
 299   //
 300   // Energy
 301   //
 302   Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
 303   return TMath::Sqrt(mass*mass+P()*P());
 304 }
 305 //---------------------------------------------------------------------------
 306 Double_t AliAODRecoDecay::EProng(Int_t ip,UInt_t pdg) const
 307 {
 308   //
 309   // Energy of ip-th prong
 310   //
 311   Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
 312   return TMath::Sqrt(mass*mass+PProng(ip)*PProng(ip));
 313 }
 314 //---------------------------------------------------------------------------
 315 /*Int_t AliAODRecoDecay::GetIndexProng(Int_t ip) const
 316 {
 317   //
 318   // Index of prong ip
 319   //
 320   if(!GetNProngs()) return 999999;
 321 UShort_t *indices = GetSecondaryVtx()->GetIndices();
 322   return indices[ip];
 323 }*/
 324 //----------------------------------------------------------------------------
 325 Double_t AliAODRecoDecay::ImpParXY(Double_t point[3]) const
 326 {
 327   //
 328   // Impact parameter in the bending plane of the particle
 329   // w.r.t. to "point"
 330   //
 331   Double_t k = -(GetSecVtxX()-point[0])*Px()-(GetSecVtxY()-point[1])*Py();
 332   k /= Pt()*Pt();
 333   Double_t dx = GetSecVtxX()-point[0]+k*Px();
 334   Double_t dy = GetSecVtxY()-point[1]+k*Py();
 335   Double_t absImpPar = TMath::Sqrt(dx*dx+dy*dy);
 336   TVector3 mom(Px(),Py(),Pz());
 337   TVector3 fline(GetSecVtxX()-point[0],
 338                  GetSecVtxY()-point[1],
 339                  GetSecVtxZ()-point[2]);
 340   TVector3 cross = mom.Cross(fline);
 341   return (cross.Z()>0. ? absImpPar : -absImpPar);
 342 }
 343 //----------------------------------------------------------------------------
 344 Double_t AliAODRecoDecay::InvMass(Int_t npdg,UInt_t *pdg) const
 345 {
 346   //
 347   // Invariant mass for prongs mass hypotheses in pdg array
 348   //
 349   if(GetNProngs()!=npdg) {
 350     printf("npdg != GetNProngs()");
 351     return (Double_t)-99999.;
 352   }
 353   Double_t energysum = 0.;
 354
 355   for(Int_t i=0; i<GetNProngs(); i++) {
 356     energysum += EProng(i,pdg[i]);
 357   }
 358
 359   Double_t mass = TMath::Sqrt(energysum*energysum-P()*P());
 360
 361   return mass;
 362 }
 363 //----------------------------------------------------------------------------
 364 Double_t AliAODRecoDecay::InvMass2Prongs(Int_t ip1,Int_t ip2,
 365                                       UInt_t pdg1,UInt_t pdg2) const
 366 {
 367   //
 368   // 2-prong(ip1,ip2) invariant mass for prongs mass hypotheses in pdg1,2
 369   //
 370   Double_t energysum = EProng(ip1,pdg1) + EProng(ip2,pdg2);
 371   Double_t psum2 = (PxProng(ip1)+PxProng(ip2))*(PxProng(ip1)+PxProng(ip2))
 372                   +(PyProng(ip1)+PyProng(ip2))*(PyProng(ip1)+PyProng(ip2))
 373                   +(PzProng(ip1)+PzProng(ip2))*(PzProng(ip1)+PzProng(ip2));
 374   Double_t mass = TMath::Sqrt(energysum*energysum-psum2);
 375
 376   return mass;
 377 }
 378 //---------------------------------------------------------------------------
 379 void AliAODRecoDecay::Print(Option_t* /*option*/) const
 380 {
 381   //
 382   // Print some information
 383   //
 384   printf("AliAODRecoDecay with %d prongs\n",GetNProngs());
 385   printf("Secondary Vertex: (%f, %f, %f)\n",GetSecVtxX(),GetSecVtxY(),GetSecVtxZ());
 386
 387   return;
 388 }
 389 //----------------------------------------------------------------------------
 390 Double_t AliAODRecoDecay::ProngsRelAngle(Int_t ip1,Int_t ip2) const
 391 {
 392   //
 393   // Relative angle between two prongs
 394   //
 395   TVector3 momA(PxProng(ip1),PyProng(ip1),PzProng(ip1));
 396   TVector3 momB(PxProng(ip2),PyProng(ip2),PzProng(ip2));
 397
 398   Double_t angle = momA.Angle(momB);
 399
 400   return angle;
 401 }
 402 //----------------------------------------------------------------------------
 403 Double_t AliAODRecoDecay::QlProng(Int_t ip) const
 404 {
 405   //
 406   // Longitudinal momentum of prong w.r.t. to total momentum
 407   //
 408   TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
 409   TVector3 momTot(Px(),Py(),Pz());
 410
 411   return mom.Dot(momTot)/momTot.Mag();
 412 }
 413 //----------------------------------------------------------------------------
 414 Double_t AliAODRecoDecay::QtProng(Int_t ip) const
 415 {
 416   //
 417   // Transverse momentum of prong w.r.t. to total momentum
 418   //
 419   TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
 420   TVector3 momTot(Px(),Py(),Pz());
 421
 422   return mom.Perp(momTot);
 423 }
 424 //----------------------------------------------------------------------------
 425 Double_t AliAODRecoDecay::QlProngFlightLine(Int_t ip,Double_t point[3]) const
 426 {
 427   //
 428   // Longitudinal momentum of prong w.r.t. to flight line between "point"
 429   // and fSecondaryVtx
 430   //
 431   TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
 432   TVector3 fline(GetSecVtxX()-point[0],
 433                  GetSecVtxY()-point[1],
 434                  GetSecVtxZ()-point[2]);
 435
 436   return mom.Dot(fline)/fline.Mag();
 437 }
 438 //----------------------------------------------------------------------------
 439 Double_t AliAODRecoDecay::QtProngFlightLine(Int_t ip,Double_t point[3]) const
 440 {
 441   //
 442   // Transverse momentum of prong w.r.t. to flight line between "point" and
 443   // fSecondaryVtx
 444   //
 445   TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
 446   TVector3 fline(GetSecVtxX()-point[0],
 447                  GetSecVtxY()-point[1],
 448                  GetSecVtxZ()-point[2]);
 449
 450   return mom.Perp(fline);
 451 }
 452 //--------------------------------------------------------------------------