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 "AliLog.h"
  26 #include "AliVTrack.h"
  27 #include "AliAODRecoDecay.h"
  28
  29 ClassImp(AliAODRecoDecay)
  30
  31 //--------------------------------------------------------------------------
  32 AliAODRecoDecay::AliAODRecoDecay() :
  33   AliVTrack(),
  34   fSecondaryVtx(0x0),
  35   fOwnSecondaryVtx(0x0),
  36   fCharge(0),
  37   fNProngs(0), fNDCA(0), fNPID(0),
  38   fPx(0x0), fPy(0x0), fPz(0x0),
  39   fd0(0x0),
  40   fDCA(0x0),
  41   fPID(0x0),
  42   fEventNumber(-1),fRunNumber(-1)
  43 {
  44   //
  45   // Default Constructor
  46   //
  47 }
  48 //--------------------------------------------------------------------------
  49 AliAODRecoDecay::AliAODRecoDecay(AliAODVertex *vtx2,Int_t nprongs,
  50                                  Short_t charge,
  51                                  Double_t *px,Double_t *py,Double_t *pz,
  52                                  Double_t *d0) :
  53   AliVTrack(),
  54   fSecondaryVtx(vtx2),
  55   fOwnSecondaryVtx(0x0),
  56   fCharge(charge),
  57   fNProngs(nprongs), fNDCA(0), fNPID(0),
  58   fPx(0x0), fPy(0x0), fPz(0x0),
  59   fd0(0x0),
  60   fDCA(0x0),
  61   fPID(0x0),
  62   fEventNumber(-1),fRunNumber(-1)
  63 {
  64   //
  65   // Constructor with AliAODVertex for decay vertex
  66   //
  67
  68   fPx = new Double_t[GetNProngs()];
  69   fPy = new Double_t[GetNProngs()];
  70   fPz = new Double_t[GetNProngs()];
  71   fd0 = new Double_t[GetNProngs()];
  72   for(Int_t i=0; i<GetNProngs(); i++) {
  73     fPx[i] = px[i];
  74     fPy[i] = py[i];
  75     fPz[i] = pz[i];
  76     fd0[i] = d0[i];
  77   }
  78 }
  79 //--------------------------------------------------------------------------
  80 AliAODRecoDecay::AliAODRecoDecay(AliAODVertex *vtx2,Int_t nprongs,
  81                                  Short_t charge,
  82                                  Double_t *d0) :
  83   AliVTrack(),
  84   fSecondaryVtx(vtx2),
  85   fOwnSecondaryVtx(0x0),
  86   fCharge(charge),
  87   fNProngs(nprongs), fNDCA(0), fNPID(0),
  88   fPx(0x0), fPy(0x0), fPz(0x0),
  89   fd0(0x0),
  90   fDCA(0x0),
  91   fPID(0x0),
  92   fEventNumber(-1),fRunNumber(-1)
  93 {
  94   //
  95   // Constructor with AliAODVertex for decay vertex and without prongs momenta
  96   //
  97
  98   fd0 = new Double_t[GetNProngs()];
  99   for(Int_t i=0; i<GetNProngs(); i++) fd0[i] = d0[i];
 100 }
 101 //--------------------------------------------------------------------------
 102 AliAODRecoDecay::AliAODRecoDecay(const AliAODRecoDecay &source) :
 103   AliVTrack(source),
 104   fSecondaryVtx(source.fSecondaryVtx),
 105   fOwnSecondaryVtx(source.fOwnSecondaryVtx),
 106   fCharge(source.fCharge),
 107   fNProngs(source.fNProngs), fNDCA(source.fNDCA), fNPID(source.fNPID),
 108   fPx(0x0), fPy(0x0), fPz(0x0),
 109   fd0(0x0),
 110   fDCA(0x0),
 111   fPID(0x0),
 112   fEventNumber(source.fEventNumber),fRunNumber(source.fRunNumber)
 113 {
 114   //
 115   // Copy constructor
 116   //
 117   if(source.GetNProngs()>0) {
 118     fd0 = new Double32_t[GetNProngs()];
 119     memcpy(fd0,source.fd0,GetNProngs()*sizeof(Double32_t));
 120     if(source.fPx) {
 121       fPx = new Double32_t[GetNProngs()];
 122       fPy = new Double32_t[GetNProngs()];
 123       fPz = new Double32_t[GetNProngs()];
 124       memcpy(fPx,source.fPx,GetNProngs()*sizeof(Double32_t));
 125       memcpy(fPy,source.fPy,GetNProngs()*sizeof(Double32_t));
 126       memcpy(fPz,source.fPz,GetNProngs()*sizeof(Double32_t));
 127     }
 128     if(source.fPID) {
 129       fPID = new Double32_t[fNPID];
 130       memcpy(fPID,source.fPID,fNPID*sizeof(Double32_t));
 131     }
 132     if(source.fDCA) {
 133       fDCA = new Double32_t[fNDCA];
 134       memcpy(fDCA,source.fDCA,fNDCA*sizeof(Double32_t));
 135     }
 136   }
 137 }
 138 //--------------------------------------------------------------------------
 139 AliAODRecoDecay &AliAODRecoDecay::operator=(const AliAODRecoDecay &source)
 140 {
 141   //
 142   // assignment operator
 143   //
 144   if(&source == this) return *this;
 145   fSecondaryVtx = source.fSecondaryVtx;
 146   fOwnSecondaryVtx = source.fOwnSecondaryVtx;
 147   fCharge = source.fCharge;
 148   fNProngs = source.fNProngs;
 149   fNDCA = source.fNDCA;
 150   fNPID = source.fNPID;
 151   fEventNumber = source.fEventNumber;
 152   fRunNumber = source.fRunNumber;
 153   if(source.GetNProngs()>0) {
 154     if(fd0)delete [] fd0;
 155     fd0 = new Double32_t[GetNProngs()];
 156     memcpy(fd0,source.fd0,GetNProngs()*sizeof(Double32_t));
 157     if(source.fPx) {
 158       if(fPx) delete [] fPx;
 159       fPx = new Double32_t[GetNProngs()];
 160       if(fPy) delete [] fPy;
 161       fPy = new Double32_t[GetNProngs()];
 162       if(fPz) delete [] fPz;
 163       fPz = new Double32_t[GetNProngs()];
 164       memcpy(fPx,source.fPx,GetNProngs()*sizeof(Double32_t));
 165       memcpy(fPy,source.fPy,GetNProngs()*sizeof(Double32_t));
 166       memcpy(fPz,source.fPz,GetNProngs()*sizeof(Double32_t));
 167     }
 168     if(source.fPID) {
 169       if(fPID) delete [] fPID;
 170       fPID = new Double32_t[fNPID];
 171       memcpy(fPID,source.fPID,fNPID*sizeof(Double32_t));
 172     }
 173     if(source.fDCA) {
 174       if(fDCA) delete [] fDCA;
 175       fDCA = new Double32_t[fNDCA];
 176       memcpy(fDCA,source.fDCA,fNDCA*sizeof(Double32_t));
 177     }
 178   }
 179   return *this;
 180 }
 181 //--------------------------------------------------------------------------
 182 AliAODRecoDecay::~AliAODRecoDecay() {
 183   //
 184   // Default Destructor
 185   //
 186   if(fPx) { delete [] fPx; fPx=NULL; }
 187   if(fPy) { delete [] fPy; fPy=NULL; }
 188   if(fPz) { delete [] fPz; fPz=NULL; }
 189   if(fd0) { delete [] fd0; fd0=NULL; }
 190   if(fPID) { delete [] fPID; fPID=NULL; }
 191   if(fDCA) { delete [] fDCA; fDCA=NULL; }
 192   if(fOwnSecondaryVtx) { delete fOwnSecondaryVtx; fOwnSecondaryVtx=NULL; }
 193 }
 194 //----------------------------------------------------------------------------
 195 Double_t AliAODRecoDecay::Alpha() const
 196 {
 197   //
 198   // Armenteros-Podolanski alpha for 2-prong decays
 199   //
 200   if(GetNProngs()!=2) {
 201     printf("Can be called only for 2-prong decays");
 202     return (Double_t)-99999.;
 203   }
 204   return 1.-2./(1.+QlProng(0)/QlProng(1));
 205 }
 206 //----------------------------------------------------------------------------
 207 Double_t AliAODRecoDecay::DecayLength(Double_t point[3]) const
 208 {
 209   //
 210   // Decay length 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                     +(point[2]-GetSecVtxZ())
 217                     *(point[2]-GetSecVtxZ()));
 218 }
 219 //----------------------------------------------------------------------------
 220 Double_t AliAODRecoDecay::DecayLengthXY(Double_t point[3]) const
 221 {
 222   //
 223   // Decay length in XY assuming it is produced at "point" [cm]
 224   //
 225   return TMath::Sqrt((point[0]-GetSecVtxX())
 226                     *(point[0]-GetSecVtxX())
 227                     +(point[1]-GetSecVtxY())
 228                     *(point[1]-GetSecVtxY()));
 229 }
 230 //----------------------------------------------------------------------------
 231 Double_t AliAODRecoDecay::CosPointingAngle(Double_t point[3]) const
 232 {
 233   //
 234   // Cosine of pointing angle in space assuming it is produced at "point"
 235   //
 236   TVector3 mom(Px(),Py(),Pz());
 237   TVector3 fline(GetSecVtxX()-point[0],
 238                  GetSecVtxY()-point[1],
 239                  GetSecVtxZ()-point[2]);
 240
 241   Double_t pta = mom.Angle(fline);
 242
 243   return TMath::Cos(pta);
 244 }
 245 //----------------------------------------------------------------------------
 246 Double_t AliAODRecoDecay::CosPointingAngleXY(Double_t point[3]) const
 247 {
 248   //
 249   // Cosine of pointing angle in transverse plane assuming it is produced
 250   // at "point"
 251   //
 252   TVector3 momXY(Px(),Py(),0.);
 253   TVector3 flineXY(GetSecVtxX()-point[0],
 254                    GetSecVtxY()-point[1],
 255                    0.);
 256
 257   Double_t ptaXY = momXY.Angle(flineXY);
 258
 259   return TMath::Cos(ptaXY);
 260 }
 261 //----------------------------------------------------------------------------
 262 Double_t AliAODRecoDecay::CosThetaStar(Int_t ip,UInt_t pdgvtx,UInt_t pdgprong0,UInt_t pdgprong1) const
 263 {
 264   //
 265   // Only for 2-prong decays:
 266   // Cosine of decay angle (theta*) in the rest frame of the mother particle
 267   // for prong ip (0 or 1) with mass hypotheses pdgvtx for mother particle,
 268   // pdgprong0 for prong 0 and pdgprong1 for prong1
 269   //
 270   if(GetNProngs()!=2) {
 271     printf("Can be called only for 2-prong decays");
 272     return (Double_t)-99999.;
 273   }
 274   Double_t massvtx = TDatabasePDG::Instance()->GetParticle(pdgvtx)->Mass();
 275   Double_t massp[2];
 276   massp[0] = TDatabasePDG::Instance()->GetParticle(pdgprong0)->Mass();
 277   massp[1] = TDatabasePDG::Instance()->GetParticle(pdgprong1)->Mass();
 278
 279   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);
 280
 281   Double_t beta = P()/E(pdgvtx);
 282   Double_t gamma = E(pdgvtx)/massvtx;
 283
 284   Double_t cts = (QlProng(ip)/gamma-beta*TMath::Sqrt(pStar*pStar+massp[ip]*massp[ip]))/pStar;
 285
 286   return cts;
 287 }
 288 //---------------------------------------------------------------------------
 289 Double_t AliAODRecoDecay::Ct(UInt_t pdg,Double_t point[3]) const
 290 {
 291   //
 292   // Decay time * c assuming it is produced at "point" [cm]
 293   //
 294   Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
 295   return DecayLength(point)*mass/P();
 296 }
 297 //---------------------------------------------------------------------------
 298 Double_t AliAODRecoDecay::E(UInt_t pdg) const
 299 {
 300   //
 301   // Energy
 302   //
 303   Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
 304   return TMath::Sqrt(mass*mass+P()*P());
 305 }
 306 //---------------------------------------------------------------------------
 307 Double_t AliAODRecoDecay::EProng(Int_t ip,UInt_t pdg) const
 308 {
 309   //
 310   // Energy of ip-th prong
 311   //
 312   Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
 313   return TMath::Sqrt(mass*mass+PProng(ip)*PProng(ip));
 314 }
 315 //--------------------------------------------------------------------------
 316 Bool_t AliAODRecoDecay::GetCovarianceXYZPxPyPz(Double_t cv[21]) const {
 317   //
 318   // This function returns the global covariance matrix of the track params
 319   //
 320   // Cov(x,x) ... :   cv[0]
 321   // Cov(y,x) ... :   cv[1]  cv[2]
 322   // Cov(z,x) ... :   cv[3]  cv[4]  cv[5]
 323   // Cov(px,x)... :   cv[6]  cv[7]  cv[8]  cv[9]
 324   // Cov(py,x)... :   cv[10] cv[11] cv[12] cv[13] cv[14]
 325   // Cov(pz,x)... :   cv[15] cv[16] cv[17] cv[18] cv[19] cv[20]
 326   //
 327   // For XYZ we take the cov of the vertex, for PxPyPz we take the
 328   // sum of the covs of PxPyPz from the daughters, for the moment
 329   // we set the cov between position and momentum as the sum of
 330   // the same cov from the daughters.
 331   //
 332
 333   Int_t j;
 334   for(j=0;j<21;j++) cv[j]=0.;
 335
 336   if(!GetNDaughters()) {
 337     AliError("No daughters available");
 338     return kFALSE;
 339   }
 340
 341   Double_t v[6];
 342   AliAODVertex *secv=GetSecondaryVtx();
 343   if(!secv) {
 344     AliError("Vertex covariance matrix not available");
 345     return kFALSE;
 346   }
 347   if(!secv->GetCovMatrix(v)) {
 348     AliError("Vertex covariance matrix not available");
 349     return kFALSE;
 350   }
 351
 352   Double_t p[21]; for(j=0;j<21;j++) p[j]=0.;
 353   Bool_t error=kFALSE;
 354   for(Int_t i=1; i<GetNDaughters(); i++) {
 355     AliVTrack *daugh = (AliVTrack*)GetDaughter(i);
 356     Double_t dcov[21];
 357     if(!daugh->GetCovarianceXYZPxPyPz(dcov)) error=kTRUE;
 358     for(j=0;j<21;j++) p[j] += dcov[j];
 359   }
 360   if(error) {
 361     AliError("No covariance for at least one daughter")
 362     return kFALSE;
 363   }
 364
 365   for(j=0; j<21; j++) {
 366     if(j<6) {
 367       cv[j] = v[j];
 368     } else {
 369       cv[j] = p[j];
 370     }
 371   }
 372
 373   return kTRUE;
 374 }
 375 //----------------------------------------------------------------------------
 376 UChar_t  AliAODRecoDecay::GetITSClusterMap() const {
 377   //
 378   // We take the logical AND of the daughters cluster maps
 379   // (only if all daughters have the bit for given layer, we set the bit)
 380   //
 381   UChar_t map=0;
 382
 383   if(!GetNDaughters()) {
 384     AliError("No daughters available");
 385     return map;
 386   }
 387
 388   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)
 389     Int_t bit = 1;
 390     for(Int_t i=0; i<GetNDaughters(); i++) {
 391       AliVTrack *daugh = (AliVTrack*)GetDaughter(i);
 392       if(!TESTBIT(daugh->GetITSClusterMap(),l)) bit=0;
 393     }
 394     if(bit) SETBIT(map,l);
 395   }
 396
 397   return map;
 398 }
 399 //--------------------------------------------------------------------------
 400 ULong_t AliAODRecoDecay::GetStatus() const {
 401   //
 402   // Same as for ITSClusterMap
 403   //
 404   ULong_t status=0;
 405
 406   if(!GetNDaughters()) {
 407     AliError("No daughters available");
 408     return status;
 409   }
 410
 411   AliVTrack *daugh0 = (AliVTrack*)GetDaughter(0);
 412   status = status&(daugh0->GetStatus());
 413
 414   for(Int_t i=1; i<GetNDaughters(); i++) {
 415     AliVTrack *daugh = (AliVTrack*)GetDaughter(i);
 416     status = status&(daugh->GetStatus());
 417   }
 418
 419   return status;
 420 }
 421 //--------------------------------------------------------------------------
 422 Double_t AliAODRecoDecay::ImpParXY(Double_t point[3]) const
 423 {
 424   //
 425   // Impact parameter in the bending plane of the particle
 426   // w.r.t. to "point"
 427   //
 428   Double_t k = -(GetSecVtxX()-point[0])*Px()-(GetSecVtxY()-point[1])*Py();
 429   k /= Pt()*Pt();
 430   Double_t dx = GetSecVtxX()-point[0]+k*Px();
 431   Double_t dy = GetSecVtxY()-point[1]+k*Py();
 432   Double_t absImpPar = TMath::Sqrt(dx*dx+dy*dy);
 433   TVector3 mom(Px(),Py(),Pz());
 434   TVector3 fline(GetSecVtxX()-point[0],
 435                  GetSecVtxY()-point[1],
 436                  GetSecVtxZ()-point[2]);
 437   TVector3 cross = mom.Cross(fline);
 438   return (cross.Z()>0. ? absImpPar : -absImpPar);
 439 }
 440 //----------------------------------------------------------------------------
 441 Double_t AliAODRecoDecay::InvMass(Int_t npdg,UInt_t *pdg) const
 442 {
 443   //
 444   // Invariant mass for prongs mass hypotheses in pdg array
 445   //
 446   if(GetNProngs()!=npdg) {
 447     printf("npdg != GetNProngs()");
 448     return (Double_t)-99999.;
 449   }
 450   Double_t energysum = 0.;
 451
 452   for(Int_t i=0; i<GetNProngs(); i++) {
 453     energysum += EProng(i,pdg[i]);
 454   }
 455
 456   Double_t mass = TMath::Sqrt(energysum*energysum-P()*P());
 457
 458   return mass;
 459 }
 460 //----------------------------------------------------------------------------
 461 Double_t AliAODRecoDecay::InvMass2Prongs(Int_t ip1,Int_t ip2,
 462                                       UInt_t pdg1,UInt_t pdg2) const
 463 {
 464   //
 465   // 2-prong(ip1,ip2) invariant mass for prongs mass hypotheses in pdg1,2
 466   //
 467   Double_t energysum = EProng(ip1,pdg1) + EProng(ip2,pdg2);
 468   Double_t psum2 = (PxProng(ip1)+PxProng(ip2))*(PxProng(ip1)+PxProng(ip2))
 469                   +(PyProng(ip1)+PyProng(ip2))*(PyProng(ip1)+PyProng(ip2))
 470                   +(PzProng(ip1)+PzProng(ip2))*(PzProng(ip1)+PzProng(ip2));
 471   Double_t mass = TMath::Sqrt(energysum*energysum-psum2);
 472
 473   return mass;
 474 }
 475 //---------------------------------------------------------------------------
 476 void AliAODRecoDecay::Print(Option_t* /*option*/) const
 477 {
 478   //
 479   // Print some information
 480   //
 481   printf("AliAODRecoDecay with %d prongs\n",GetNProngs());
 482   printf("Secondary Vertex: (%f, %f, %f)\n",GetSecVtxX(),GetSecVtxY(),GetSecVtxZ());
 483
 484   return;
 485 }
 486 //----------------------------------------------------------------------------
 487 Double_t AliAODRecoDecay::ProngsRelAngle(Int_t ip1,Int_t ip2) const
 488 {
 489   //
 490   // Relative angle between two prongs
 491   //
 492   TVector3 momA(PxProng(ip1),PyProng(ip1),PzProng(ip1));
 493   TVector3 momB(PxProng(ip2),PyProng(ip2),PzProng(ip2));
 494
 495   Double_t angle = momA.Angle(momB);
 496
 497   return angle;
 498 }
 499 //----------------------------------------------------------------------------
 500 Double_t AliAODRecoDecay::QlProng(Int_t ip) const
 501 {
 502   //
 503   // Longitudinal momentum of prong w.r.t. to total momentum
 504   //
 505   TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
 506   TVector3 momTot(Px(),Py(),Pz());
 507
 508   return mom.Dot(momTot)/momTot.Mag();
 509 }
 510 //----------------------------------------------------------------------------
 511 Double_t AliAODRecoDecay::QtProng(Int_t ip) const
 512 {
 513   //
 514   // Transverse momentum of prong w.r.t. to total momentum
 515   //
 516   TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
 517   TVector3 momTot(Px(),Py(),Pz());
 518
 519   return mom.Perp(momTot);
 520 }
 521 //----------------------------------------------------------------------------
 522 Double_t AliAODRecoDecay::QlProngFlightLine(Int_t ip,Double_t point[3]) const
 523 {
 524   //
 525   // Longitudinal momentum of prong w.r.t. to flight line between "point"
 526   // and fSecondaryVtx
 527   //
 528   TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
 529   TVector3 fline(GetSecVtxX()-point[0],
 530                  GetSecVtxY()-point[1],
 531                  GetSecVtxZ()-point[2]);
 532
 533   return mom.Dot(fline)/fline.Mag();
 534 }
 535 //----------------------------------------------------------------------------
 536 Double_t AliAODRecoDecay::QtProngFlightLine(Int_t ip,Double_t point[3]) const
 537 {
 538   //
 539   // Transverse momentum of prong w.r.t. to flight line between "point" and
 540   // fSecondaryVtx
 541   //
 542   TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
 543   TVector3 fline(GetSecVtxX()-point[0],
 544                  GetSecVtxY()-point[1],
 545                  GetSecVtxZ()-point[2]);
 546
 547   return mom.Perp(fline);
 548 }
 549 //--------------------------------------------------------------------------