PWGDQ/dielectron/AliDielectronPair.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2009, 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 //  Dielectron Pair class. Internally it makes use of AliKFParticle.     //
  19 //                                                                       //
  20 ///////////////////////////////////////////////////////////////////////////
  21
  22
  23 #include <TDatabasePDG.h>
  24 #include <AliVTrack.h>
  25 #include <AliVVertex.h>
  26 #include <AliPID.h>
  27 #include <AliExternalTrackParam.h>
  28
  29 #include "AliDielectronPair.h"
  30
  31 ClassImp(AliDielectronPair)
  32
  33 AliDielectronPair::AliDielectronPair() :
  34   fType(-1),
  35   fLabel(-1),
  36   fPdgCode(0),
  37   fPair(),
  38   fD1(),
  39   fD2(),
  40   fRefD1(),
  41   fRefD2(),
  42   fKFUsage(kTRUE)
  43 {
  44   //
  45   // Default Constructor
  46   //
  47
  48 }
  49
  50 //______________________________________________
  51 AliDielectronPair::AliDielectronPair(AliVTrack * const particle1, Int_t pid1,
  52                                      AliVTrack * const particle2, Int_t pid2, Char_t type) :
  53   fType(type),
  54   fLabel(-1),
  55   fPdgCode(0),
  56   fPair(),
  57   fD1(),
  58   fD2(),
  59   fRefD1(),
  60   fRefD2(),
  61   fKFUsage(kTRUE)
  62 {
  63   //
  64   // Constructor with tracks
  65   //
  66   SetTracks(particle1, pid1, particle2, pid2);
  67 }
  68
  69 //______________________________________________
  70 AliDielectronPair::AliDielectronPair(const AliKFParticle * const particle1,
  71                                      const AliKFParticle * const particle2,
  72                                      AliVTrack * const refParticle1,
  73                                      AliVTrack * const refParticle2, Char_t type) :
  74   fType(type),
  75   fLabel(-1),
  76   fPdgCode(0),
  77   fPair(),
  78   fD1(),
  79   fD2(),
  80   fRefD1(),
  81   fRefD2(),
  82   fKFUsage(kTRUE)
  83 {
  84   //
  85   // Constructor with tracks
  86   //
  87   SetTracks(particle1, particle2,refParticle1,refParticle2);
  88 }
  89
  90 //______________________________________________
  91 AliDielectronPair::~AliDielectronPair()
  92 {
  93   //
  94   // Default Destructor
  95   //
  96
  97 }
  98
  99 //______________________________________________
 100 void AliDielectronPair::SetTracks(AliVTrack * const particle1, Int_t pid1,
 101                                   AliVTrack * const particle2, Int_t pid2)
 102 {
 103   //
 104   // Sort particles by pt, first particle larget Pt
 105   // set AliKF daughters and pair
 106   // refParticle1 and 2 are the original tracks. In the case of track rotation
 107   // they are needed in the framework
 108   //
 109   fPair.Initialize();
 110   fD1.Initialize();
 111   fD2.Initialize();
 112
 113   AliKFParticle kf1(*particle1,pid1);
 114   AliKFParticle kf2(*particle2,pid2);
 115
 116   fPair.AddDaughter(kf1);
 117   fPair.AddDaughter(kf2);
 118
 119   if (particle1->Pt()>particle2->Pt()){
 120     fRefD1 = particle1;
 121     fRefD2 = particle2;
 122     fD1+=kf1;
 123     fD2+=kf2;
 124   } else {
 125     fRefD1 = particle2;
 126     fRefD2 = particle1;
 127     fD1+=kf2;
 128     fD2+=kf1;
 129   }
 130 }
 131
 132 //______________________________________________
 133 void AliDielectronPair::SetTracks(const AliKFParticle * const particle1,
 134                                   const AliKFParticle * const particle2,
 135                                   AliVTrack * const refParticle1,
 136                                   AliVTrack * const refParticle2)
 137 {
 138   //
 139   // Sort particles by pt, first particle larget Pt
 140   // set AliKF daughters and pair
 141   // refParticle1 and 2 are the original tracks. In the case of track rotation
 142   // they are needed in the framework
 143   //
 144   fPair.Initialize();
 145   fD1.Initialize();
 146   fD2.Initialize();
 147
 148   AliKFParticle kf1(*particle1);
 149   AliKFParticle kf2(*particle2);
 150
 151   fPair.AddDaughter(kf1);
 152   fPair.AddDaughter(kf2);
 153
 154   if (kf1.GetPt()>kf2.GetPt()){
 155     fRefD1 = refParticle1;
 156     fRefD2 = refParticle2;
 157     fD1+=kf1;
 158     fD2+=kf2;
 159   } else {
 160     fRefD1 = refParticle2;
 161     fRefD2 = refParticle1;
 162     fD1+=kf2;
 163     fD2+=kf1;
 164   }
 165 }
 166
 167 //______________________________________________
 168 void AliDielectronPair::GetThetaPhiCM(Double_t &thetaHE, Double_t &phiHE, Double_t &thetaCS, Double_t &phiCS) const
 169 {
 170   //
 171   // Calculate theta and phi in helicity and Collins-Soper coordinate frame
 172   //
 173   const Double_t kBeamEnergy   = 3500.;
 174   Double_t pxyz1[3]={fD1.GetPx(),fD1.GetPy(),fD1.GetPz()};
 175   Double_t pxyz2[3]={fD2.GetPx(),fD2.GetPy(),fD2.GetPz()};
 176   Double_t eleMass=AliPID::ParticleMass(AliPID::kElectron);
 177   Double_t proMass=AliPID::ParticleMass(AliPID::kProton);
 178
 179 //   AliVParticle *d1 = static_cast<AliVParticle*>(fRefD1.GetObject());
 180 //   AliVParticle *d2 = static_cast<AliVParticle*>(fRefD2.GetObject());
 181
 182 //   d1->PxPyPz(pxyz1);
 183 //   d2->PxPyPz(pxyz2);
 184
 185   TLorentzVector projMom(0.,0.,-kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
 186   TLorentzVector targMom(0.,0., kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
 187
 188   // first & second daughter 4-mom
 189   TLorentzVector p1Mom(pxyz1[0],pxyz1[1],pxyz1[2],
 190                        TMath::Sqrt(pxyz1[0]*pxyz1[0]+pxyz1[1]*pxyz1[1]+pxyz1[2]*pxyz1[2]+eleMass*eleMass));
 191   TLorentzVector p2Mom(pxyz2[0],pxyz2[1],pxyz2[2],
 192                        TMath::Sqrt(pxyz2[0]*pxyz2[0]+pxyz2[1]*pxyz2[1]+pxyz2[2]*pxyz2[2]+eleMass*eleMass));
 193   // J/Psi 4-momentum vector
 194   TLorentzVector motherMom=p1Mom+p2Mom;
 195
 196   // boost all the 4-mom vectors to the mother rest frame
 197   TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect();
 198   p1Mom.Boost(beta);
 199   p2Mom.Boost(beta);
 200   projMom.Boost(beta);
 201   targMom.Boost(beta);
 202
 203     // x,y,z axes
 204   TVector3 zAxisHE = (motherMom.Vect()).Unit();
 205   TVector3 zAxisCS = ((projMom.Vect()).Unit()-(targMom.Vect()).Unit()).Unit();
 206   TVector3 yAxis = ((projMom.Vect()).Cross(targMom.Vect())).Unit();
 207   TVector3 xAxisHE = (yAxis.Cross(zAxisHE)).Unit();
 208   TVector3 xAxisCS = (yAxis.Cross(zAxisCS)).Unit();
 209
 210   // fill theta and phi
 211   if(fD1.GetQ()>0){
 212     thetaHE = zAxisHE.Dot((p1Mom.Vect()).Unit());
 213     thetaCS = zAxisCS.Dot((p1Mom.Vect()).Unit());
 214     phiHE   = TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxisHE));
 215     phiCS   = TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxisCS));
 216   } else {
 217     thetaHE = zAxisHE.Dot((p2Mom.Vect()).Unit());
 218     thetaCS = zAxisCS.Dot((p2Mom.Vect()).Unit());
 219     phiHE   = TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxisHE));
 220     phiCS   = TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxisCS));
 221   }
 222 }
 223
 224 //______________________________________________
 225 Double_t AliDielectronPair::PsiPair(Double_t MagField) const
 226 {
 227   //Following idea to use opening of colinear pairs in magnetic field from e.g. PHENIX
 228   //to ID conversions. Adapted from AliTRDv0Info class
 229   Double_t x, y, z;
 230   x = fPair.GetX();
 231   y = fPair.GetY();
 232   z = fPair.GetZ();
 233
 234   Double_t m1[3] = {0,0,0};
 235   Double_t m2[3] = {0,0,0};
 236
 237   m1[0] = fD1.GetPx();
 238   m1[1] = fD1.GetPy();
 239   m1[2] = fD1.GetPz();
 240
 241   m2[0] = fD2.GetPx();
 242   m2[1] = fD2.GetPy();
 243   m2[2] = fD2.GetPz();
 244
 245   Double_t deltat = 1.;
 246   deltat = TMath::ATan(m2[2]/(TMath::Sqrt(m2[0]*m2[0] + m2[1]*m2[1])+1.e-13))-
 247         TMath::ATan(m1[2]/(TMath::Sqrt(m1[0]*m1[0] + m1[1]*m1[1])+1.e-13));//difference of angles of the two daughter tracks with z-axis
 248
 249   Double_t radiussum = TMath::Sqrt(x*x + y*y) + 50;//radius to which tracks shall be propagated
 250
 251   Double_t mom1Prop[3];
 252   Double_t mom2Prop[3];
 253
 254   AliExternalTrackParam *d1 = static_cast<AliExternalTrackParam*>(fRefD1.GetObject());
 255   AliExternalTrackParam *d2 = static_cast<AliExternalTrackParam*>(fRefD2.GetObject());
 256
 257   AliExternalTrackParam nt(*d1), pt(*d2);
 258
 259   Double_t fPsiPair = 4.;
 260   if(nt.PropagateTo(radiussum,MagField) == 0)//propagate tracks to the outside
 261         fPsiPair =  -5.;
 262   if(pt.PropagateTo(radiussum,MagField) == 0)
 263         fPsiPair = -5.;
 264   pt.GetPxPyPz(mom1Prop);//Get momentum vectors of tracks after propagation
 265   nt.GetPxPyPz(mom2Prop);
 266
 267
 268
 269   Double_t pEle =
 270         TMath::Sqrt(mom2Prop[0]*mom2Prop[0]+mom2Prop[1]*mom2Prop[1]+mom2Prop[2]*mom2Prop[2]);//absolute momentum val
 271   Double_t pPos =
 272         TMath::Sqrt(mom1Prop[0]*mom1Prop[0]+mom1Prop[1]*mom1Prop[1]+mom1Prop[2]*mom1Prop[2]);//absolute momentum val
 273
 274   Double_t scalarproduct =
 275         mom1Prop[0]*mom2Prop[0]+mom1Prop[1]*mom2Prop[1]+mom1Prop[2]*mom2Prop[2];//scalar product of propagated posit
 276
 277   Double_t chipair = TMath::ACos(scalarproduct/(pEle*pPos));//Angle between propagated daughter tracks
 278
 279   fPsiPair =  TMath::Abs(TMath::ASin(deltat/chipair));
 280
 281   return fPsiPair;
 282
 283 }
 284
 285 //______________________________________________
 286 Double_t AliDielectronPair::ThetaPhiCM(const AliVParticle* d1, const AliVParticle* d2,
 287                                        const Bool_t isHE, const Bool_t isTheta)
 288 {
 289   // The function calculates theta and phi in the mother rest frame with
 290   // respect to the helicity coordinate system and Collins-Soper coordinate system
 291   // TO DO: generalize for different decays (only J/Psi->e+e- now)
 292
 293   // Laboratory frame 4-vectors:
 294   // projectile beam & target beam 4-mom
 295   // TODO: need to retrieve the beam energy from somewhere
 296   const Double_t kBeamEnergy   = 3500.;
 297   Double_t px1=d1->Px();
 298   Double_t py1=d1->Py();
 299   Double_t pz1=d1->Pz();
 300   Double_t px2=d2->Px();
 301   Double_t py2=d2->Py();
 302   Double_t pz2=d2->Pz();
 303   Double_t eleMass=AliPID::ParticleMass(AliPID::kElectron);
 304   Double_t proMass=AliPID::ParticleMass(AliPID::kProton);
 305
 306   TLorentzVector projMom(0.,0.,-kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
 307   TLorentzVector targMom(0.,0., kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
 308
 309   // first & second daughter 4-mom
 310   TLorentzVector p1Mom(px1,py1,pz1,TMath::Sqrt(px1*px1+py1*py1+pz1*pz1+eleMass*eleMass));
 311   TLorentzVector p2Mom(px2,py2,pz2,TMath::Sqrt(px2*px2+py2*py2+pz2*pz2+eleMass*eleMass));
 312   // J/Psi 4-momentum vector
 313   TLorentzVector motherMom=p1Mom+p2Mom;
 314
 315   // boost all the 4-mom vectors to the mother rest frame
 316   TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect();
 317   p1Mom.Boost(beta);
 318   p2Mom.Boost(beta);
 319   projMom.Boost(beta);
 320   targMom.Boost(beta);
 321
 322   // x,y,z axes
 323   TVector3 zAxis;
 324   if(isHE) zAxis = (motherMom.Vect()).Unit();
 325   else zAxis = ((projMom.Vect()).Unit()-(targMom.Vect()).Unit()).Unit();
 326   TVector3 yAxis = ((projMom.Vect()).Cross(targMom.Vect())).Unit();
 327   TVector3 xAxis = (yAxis.Cross(zAxis)).Unit();
 328
 329   // return either theta or phi
 330   if(isTheta) {
 331     if(d1->Charge()>0)
 332       return zAxis.Dot((p1Mom.Vect()).Unit());
 333     else
 334       return zAxis.Dot((p2Mom.Vect()).Unit());
 335
 336   }
 337   else {
 338     if(d1->Charge()>0)
 339       return TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxis));
 340     else
 341       return TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxis));
 342   }
 343 }
 344
 345 //______________________________________________
 346 Double_t AliDielectronPair::ThetaPhiCM(const Bool_t isHE, const Bool_t isTheta) const {
 347   // The function calculates theta and phi in the mother rest frame with
 348   // respect to the helicity coordinate system and Collins-Soper coordinate system
 349   // TO DO: generalize for different decays (only J/Psi->e+e- now)
 350
 351   // Laboratory frame 4-vectors:
 352   // projectile beam & target beam 4-mom
 353   AliVParticle *d1 = static_cast<AliVParticle*>(fRefD1.GetObject());
 354   AliVParticle *d2 = static_cast<AliVParticle*>(fRefD2.GetObject());
 355
 356   const Double_t kBeamEnergy   = 3500.;
 357   Double_t px1=d1->Px();
 358   Double_t py1=d1->Py();
 359   Double_t pz1=d1->Pz();
 360   Double_t px2=d2->Px();
 361   Double_t py2=d2->Py();
 362   Double_t pz2=d2->Pz();
 363   Double_t eleMass=AliPID::ParticleMass(AliPID::kElectron);
 364   Double_t proMass=AliPID::ParticleMass(AliPID::kProton);
 365
 366   TLorentzVector projMom(0.,0.,-kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
 367   TLorentzVector targMom(0.,0., kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
 368
 369   // first & second daughter 4-mom
 370   // first & second daughter 4-mom
 371   TLorentzVector p1Mom(px1,py1,pz1,TMath::Sqrt(px1*px1+py1*py1+pz1*pz1+eleMass*eleMass));
 372   TLorentzVector p2Mom(px2,py2,pz2,TMath::Sqrt(px2*px2+py2*py2+pz2*pz2+eleMass*eleMass));
 373   // J/Psi 4-momentum vector
 374   TLorentzVector motherMom=p1Mom+p2Mom;
 375
 376   // boost all the 4-mom vectors to the mother rest frame
 377   TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect();
 378   p1Mom.Boost(beta);
 379   p2Mom.Boost(beta);
 380   projMom.Boost(beta);
 381   targMom.Boost(beta);
 382
 383   // x,y,z axes
 384   TVector3 zAxis;
 385   if(isHE) zAxis = (motherMom.Vect()).Unit();
 386   else zAxis = ((projMom.Vect()).Unit()-(targMom.Vect()).Unit()).Unit();
 387   TVector3 yAxis = ((projMom.Vect()).Cross(targMom.Vect())).Unit();
 388   TVector3 xAxis = (yAxis.Cross(zAxis)).Unit();
 389
 390   // return either theta or phi
 391   if(isTheta) {
 392     if(fD1.GetQ()>0)
 393       return zAxis.Dot((p1Mom.Vect()).Unit());
 394     else
 395       return zAxis.Dot((p2Mom.Vect()).Unit());
 396   }
 397   else {
 398     if(fD1.GetQ()>0)
 399       return TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxis));
 400     else
 401       return TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxis));
 402   }
 403 }
 404 //______________________________________________
 405 Double_t AliDielectronPair::GetCosPointingAngle(const AliVVertex *primVtx) const
 406 {
 407   //
 408   // Calculate the poiting angle of the pair to the primary vertex and take the cosine
 409   //
 410   if(!primVtx) return -1.;
 411
 412   Double_t deltaPos[3]; //vector between the reference point and the V0 vertex
 413   deltaPos[0] = fPair.GetX() - primVtx->GetX();
 414   deltaPos[1] = fPair.GetY() - primVtx->GetY();
 415   deltaPos[2] = fPair.GetZ() - primVtx->GetZ();
 416
 417   Double_t momV02    = fPair.GetPx()*fPair.GetPx() + fPair.GetPy()*fPair.GetPy() + fPair.GetPz()*fPair.GetPz();
 418   Double_t deltaPos2 = deltaPos[0]*deltaPos[0] + deltaPos[1]*deltaPos[1] + deltaPos[2]*deltaPos[2];
 419
 420   Double_t cosinePointingAngle = (deltaPos[0]*fPair.GetPx() + deltaPos[1]*fPair.GetPy() + deltaPos[2]*fPair.GetPz()) / TMath::Sqrt(momV02 * deltaPos2);
 421
 422   return TMath::Abs(cosinePointingAngle);
 423
 424 }
 425
 426 // //______________________________________________
 427 // Double_t AliDielectronPair::GetLXY(const AliVVertex * const vtx) const
 428 // {
 429 //   //
 430 //   // Calculate the decay length in XY taking into account the primary vertex position
 431 //   //
 432 //   if(!vtx) return 0;
 433 //   return ( (Xv()-vtx->GetX()) * Px() + (Yv()-vtx->GetY()) * Py() )/Pt()  ;
 434 // }
 435
 436 // //______________________________________________
 437 // Double_t AliDielectronPair::GetPseudoProperTime(const AliVVertex * const vtx) const
 438 // {
 439 //   //
 440 //   // Calculate the pseudo proper time
 441 //   //
 442 //   Double_t lxy=GetLXY(vtx);
 443 //   Double_t psProperDecayLength = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/Pt();
 444 //   return psProperDecayLength;
 445 // }