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
  28 #include "AliDielectronPair.h"
  29
  30 ClassImp(AliDielectronPair)
  31
  32 AliDielectronPair::AliDielectronPair() :
  33   fType(-1),
  34   fLabel(-1),
  35   fPair(),
  36   fD1(),
  37   fD2(),
  38   fRefD1(),
  39   fRefD2()
  40 {
  41   //
  42   // Default Constructor
  43   //
  44
  45 }
  46
  47 //______________________________________________
  48 AliDielectronPair::AliDielectronPair(AliVTrack * const particle1, Int_t pid1,
  49                                      AliVTrack * const particle2, Int_t pid2, Char_t type) :
  50   fType(type),
  51   fLabel(-1),
  52   fPair(),
  53   fD1(),
  54   fD2(),
  55   fRefD1(),
  56   fRefD2()
  57 {
  58   //
  59   // Constructor with tracks
  60   //
  61   SetTracks(particle1, pid1, particle2, pid2);
  62 }
  63
  64 //______________________________________________
  65 AliDielectronPair::AliDielectronPair(const AliKFParticle * const particle1,
  66                                      const AliKFParticle * const particle2,
  67                                      AliVTrack * const refParticle1,
  68                                      AliVTrack * const refParticle2, Char_t type) :
  69   fType(type),
  70   fLabel(-1),
  71   fPair(),
  72   fD1(),
  73   fD2(),
  74   fRefD1(),
  75   fRefD2()
  76 {
  77   //
  78   // Constructor with tracks
  79   //
  80   SetTracks(particle1, particle2,refParticle1,refParticle2);
  81 }
  82
  83 //______________________________________________
  84 AliDielectronPair::~AliDielectronPair()
  85 {
  86   //
  87   // Default Destructor
  88   //
  89
  90 }
  91
  92 //______________________________________________
  93 void AliDielectronPair::SetTracks(AliVTrack * const particle1, Int_t pid1,
  94                                   AliVTrack * const particle2, Int_t pid2)
  95 {
  96   //
  97   // Sort particles by pt, first particle larget Pt
  98   // set AliKF daughters and pair
  99   // refParticle1 and 2 are the original tracks. In the case of track rotation
 100   // they are needed in the framework
 101   //
 102   fPair.Initialize();
 103   fD1.Initialize();
 104   fD2.Initialize();
 105
 106   AliKFParticle kf1(*particle1,pid1);
 107   AliKFParticle kf2(*particle2,pid2);
 108
 109   fPair.AddDaughter(kf1);
 110   fPair.AddDaughter(kf2);
 111
 112   if (particle1->Pt()>particle2->Pt()){
 113     fRefD1 = particle1;
 114     fRefD2 = particle2;
 115     fD1+=kf1;
 116     fD2+=kf2;
 117   } else {
 118     fRefD1 = particle2;
 119     fRefD2 = particle1;
 120     fD1+=kf2;
 121     fD2+=kf1;
 122   }
 123 }
 124
 125 //______________________________________________
 126 void AliDielectronPair::SetTracks(const AliKFParticle * const particle1,
 127                                   const AliKFParticle * const particle2,
 128                                   AliVTrack * const refParticle1,
 129                                   AliVTrack * const refParticle2)
 130 {
 131   //
 132   // Sort particles by pt, first particle larget Pt
 133   // set AliKF daughters and pair
 134   // refParticle1 and 2 are the original tracks. In the case of track rotation
 135   // they are needed in the framework
 136   //
 137   fPair.Initialize();
 138   fD1.Initialize();
 139   fD2.Initialize();
 140
 141   AliKFParticle kf1(*particle1);
 142   AliKFParticle kf2(*particle2);
 143
 144   fPair.AddDaughter(kf1);
 145   fPair.AddDaughter(kf2);
 146
 147   if (kf1.GetPt()>kf2.GetPt()){
 148     fRefD1 = refParticle1;
 149     fRefD2 = refParticle2;
 150     fD1+=kf1;
 151     fD2+=kf2;
 152   } else {
 153     fRefD1 = refParticle2;
 154     fRefD2 = refParticle1;
 155     fD1+=kf2;
 156     fD2+=kf1;
 157   }
 158 }
 159
 160 //______________________________________________
 161 void AliDielectronPair::GetThetaPhiCM(Double_t &thetaHE, Double_t &phiHE, Double_t &thetaCS, Double_t &phiCS) const
 162 {
 163   //
 164   // Calculate theta and phi in helicity and Collins-Soper coordinate frame
 165   //
 166   const Double_t kBeamEnergy   = 3500.;
 167   Double_t pxyz1[3]={fD1.GetPx(),fD1.GetPy(),fD1.GetPz()};
 168   Double_t pxyz2[3]={fD2.GetPx(),fD2.GetPy(),fD2.GetPz()};
 169   Double_t eleMass=AliPID::ParticleMass(AliPID::kElectron);
 170   Double_t proMass=AliPID::ParticleMass(AliPID::kProton);
 171
 172 //   AliVParticle *d1 = static_cast<AliVParticle*>(fRefD1.GetObject());
 173 //   AliVParticle *d2 = static_cast<AliVParticle*>(fRefD2.GetObject());
 174
 175 //   d1->PxPyPz(pxyz1);
 176 //   d2->PxPyPz(pxyz2);
 177
 178   TLorentzVector projMom(0.,0.,-kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
 179   TLorentzVector targMom(0.,0., kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
 180
 181   // first & second daughter 4-mom
 182   TLorentzVector p1Mom(pxyz1[0],pxyz1[1],pxyz1[2],
 183                        TMath::Sqrt(pxyz1[0]*pxyz1[0]+pxyz1[1]*pxyz1[1]+pxyz1[2]*pxyz1[2]+eleMass*eleMass));
 184   TLorentzVector p2Mom(pxyz2[0],pxyz2[1],pxyz2[2],
 185                        TMath::Sqrt(pxyz2[0]*pxyz2[0]+pxyz2[1]*pxyz2[1]+pxyz2[2]*pxyz2[2]+eleMass*eleMass));
 186   // J/Psi 4-momentum vector
 187   TLorentzVector motherMom=p1Mom+p2Mom;
 188
 189   // boost all the 4-mom vectors to the mother rest frame
 190   TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect();
 191   p1Mom.Boost(beta);
 192   p2Mom.Boost(beta);
 193   projMom.Boost(beta);
 194   targMom.Boost(beta);
 195
 196     // x,y,z axes
 197   TVector3 zAxisHE = (motherMom.Vect()).Unit();
 198   TVector3 zAxisCS = ((projMom.Vect()).Unit()-(targMom.Vect()).Unit()).Unit();
 199   TVector3 yAxis = ((projMom.Vect()).Cross(targMom.Vect())).Unit();
 200   TVector3 xAxisHE = (yAxis.Cross(zAxisHE)).Unit();
 201   TVector3 xAxisCS = (yAxis.Cross(zAxisCS)).Unit();
 202
 203   // fill theta and phi
 204   if(fD1.GetQ()>0){
 205     thetaHE = zAxisHE.Dot((p1Mom.Vect()).Unit());
 206     thetaCS = zAxisCS.Dot((p1Mom.Vect()).Unit());
 207     phiHE   = TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxisHE));
 208     phiCS   = TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxisCS));
 209   } else {
 210     thetaHE = zAxisHE.Dot((p2Mom.Vect()).Unit());
 211     thetaCS = zAxisCS.Dot((p2Mom.Vect()).Unit());
 212     phiHE   = TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxisHE));
 213     phiCS   = TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxisCS));
 214   }
 215 }
 216
 217 //______________________________________________
 218 Double_t AliDielectronPair::ThetaPhiCM(const AliVParticle* d1, const AliVParticle* d2,
 219                                        const Bool_t isHE, const Bool_t isTheta)
 220 {
 221   // The function calculates theta and phi in the mother rest frame with
 222   // respect to the helicity coordinate system and Collins-Soper coordinate system
 223   // TO DO: generalize for different decays (only J/Psi->e+e- now)
 224
 225   // Laboratory frame 4-vectors:
 226   // projectile beam & target beam 4-mom
 227   // TODO: need to retrieve the beam energy from somewhere
 228   const Double_t kBeamEnergy   = 3500.;
 229   Double_t px1=d1->Px();
 230   Double_t py1=d1->Py();
 231   Double_t pz1=d1->Pz();
 232   Double_t px2=d2->Px();
 233   Double_t py2=d2->Py();
 234   Double_t pz2=d2->Pz();
 235   Double_t eleMass=AliPID::ParticleMass(AliPID::kElectron);
 236   Double_t proMass=AliPID::ParticleMass(AliPID::kProton);
 237
 238   TLorentzVector projMom(0.,0.,-kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
 239   TLorentzVector targMom(0.,0., kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
 240
 241   // first & second daughter 4-mom
 242   TLorentzVector p1Mom(px1,py1,pz1,TMath::Sqrt(px1*px1+py1*py1+pz1*pz1+eleMass*eleMass));
 243   TLorentzVector p2Mom(px2,py2,pz2,TMath::Sqrt(px2*px2+py2*py2+pz2*pz2+eleMass*eleMass));
 244   // J/Psi 4-momentum vector
 245   TLorentzVector motherMom=p1Mom+p2Mom;
 246
 247   // boost all the 4-mom vectors to the mother rest frame
 248   TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect();
 249   p1Mom.Boost(beta);
 250   p2Mom.Boost(beta);
 251   projMom.Boost(beta);
 252   targMom.Boost(beta);
 253
 254   // x,y,z axes
 255   TVector3 zAxis;
 256   if(isHE) zAxis = (motherMom.Vect()).Unit();
 257   else zAxis = ((projMom.Vect()).Unit()-(targMom.Vect()).Unit()).Unit();
 258   TVector3 yAxis = ((projMom.Vect()).Cross(targMom.Vect())).Unit();
 259   TVector3 xAxis = (yAxis.Cross(zAxis)).Unit();
 260
 261   // return either theta or phi
 262   if(isTheta) {
 263     if(d1->Charge()>0)
 264       return zAxis.Dot((p1Mom.Vect()).Unit());
 265     else
 266       return zAxis.Dot((p2Mom.Vect()).Unit());
 267
 268   }
 269   else {
 270     if(d1->Charge()>0)
 271       return TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxis));
 272     else
 273       return TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxis));
 274   }
 275 }
 276
 277 //______________________________________________
 278 Double_t AliDielectronPair::ThetaPhiCM(const Bool_t isHE, const Bool_t isTheta) const {
 279   // The function calculates theta and phi in the mother rest frame with
 280   // respect to the helicity coordinate system and Collins-Soper coordinate system
 281   // TO DO: generalize for different decays (only J/Psi->e+e- now)
 282
 283   // Laboratory frame 4-vectors:
 284   // projectile beam & target beam 4-mom
 285   AliVParticle *d1 = static_cast<AliVParticle*>(fRefD1.GetObject());
 286   AliVParticle *d2 = static_cast<AliVParticle*>(fRefD2.GetObject());
 287
 288   const Double_t kBeamEnergy   = 3500.;
 289   Double_t px1=d1->Px();
 290   Double_t py1=d1->Py();
 291   Double_t pz1=d1->Pz();
 292   Double_t px2=d2->Px();
 293   Double_t py2=d2->Py();
 294   Double_t pz2=d2->Pz();
 295   Double_t eleMass=AliPID::ParticleMass(AliPID::kElectron);
 296   Double_t proMass=AliPID::ParticleMass(AliPID::kProton);
 297
 298   TLorentzVector projMom(0.,0.,-kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
 299   TLorentzVector targMom(0.,0., kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
 300
 301   // first & second daughter 4-mom
 302   // first & second daughter 4-mom
 303   TLorentzVector p1Mom(px1,py1,pz1,TMath::Sqrt(px1*px1+py1*py1+pz1*pz1+eleMass*eleMass));
 304   TLorentzVector p2Mom(px2,py2,pz2,TMath::Sqrt(px2*px2+py2*py2+pz2*pz2+eleMass*eleMass));
 305   // J/Psi 4-momentum vector
 306   TLorentzVector motherMom=p1Mom+p2Mom;
 307
 308   // boost all the 4-mom vectors to the mother rest frame
 309   TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect();
 310   p1Mom.Boost(beta);
 311   p2Mom.Boost(beta);
 312   projMom.Boost(beta);
 313   targMom.Boost(beta);
 314
 315   // x,y,z axes
 316   TVector3 zAxis;
 317   if(isHE) zAxis = (motherMom.Vect()).Unit();
 318   else zAxis = ((projMom.Vect()).Unit()-(targMom.Vect()).Unit()).Unit();
 319   TVector3 yAxis = ((projMom.Vect()).Cross(targMom.Vect())).Unit();
 320   TVector3 xAxis = (yAxis.Cross(zAxis)).Unit();
 321
 322   // return either theta or phi
 323   if(isTheta) {
 324     if(fD1.GetQ()>0)
 325       return zAxis.Dot((p1Mom.Vect()).Unit());
 326     else
 327       return zAxis.Dot((p2Mom.Vect()).Unit());
 328   }
 329   else {
 330     if(fD1.GetQ()>0)
 331       return TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxis));
 332     else
 333       return TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxis));
 334   }
 335 }
 336
 337 // //______________________________________________
 338 // Double_t AliDielectronPair::GetLXY(const AliVVertex * const vtx) const
 339 // {
 340 //   //
 341 //   // Calculate the decay length in XY taking into account the primary vertex position
 342 //   //
 343 //   if(!vtx) return 0;
 344 //   return ( (Xv()-vtx->GetX()) * Px() + (Yv()-vtx->GetY()) * Py() )/Pt()  ;
 345 // }
 346
 347 // //______________________________________________
 348 // Double_t AliDielectronPair::GetPseudoProperTime(const AliVVertex * const vtx) const
 349 // {
 350 //   //
 351 //   // Calculate the pseudo proper time
 352 //   //
 353 //   Double_t lxy=GetLXY(vtx);
 354 //   Double_t psProperDecayLength = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/Pt();
 355 //   return psProperDecayLength;
 356 // }