PWG3/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 "AliDielectronPair.h"
  24 #include "AliVTrack.h"
  25 #include "AliPID.h"
  26
  27 ClassImp(AliDielectronPair)
  28
  29 AliDielectronPair::AliDielectronPair() :
  30   fType(-1),
  31   fLabel(-1),
  32   fPair(),
  33   fD1(),
  34   fD2(),
  35   fRefD1(),
  36   fRefD2()
  37 {
  38   //
  39   // Default Constructor
  40   //
  41
  42 }
  43
  44 //______________________________________________
  45 AliDielectronPair::AliDielectronPair(AliVTrack * const particle1, Int_t pid1,
  46                                      AliVTrack * const particle2, Int_t pid2, Char_t type) :
  47   fType(type),
  48   fLabel(-1),
  49   fPair(),
  50   fD1(),
  51   fD2(),
  52   fRefD1(),
  53   fRefD2()
  54 {
  55   //
  56   // Constructor with tracks
  57   //
  58   SetTracks(particle1, pid1, particle2, pid2);
  59 }
  60
  61 //______________________________________________
  62 AliDielectronPair::~AliDielectronPair()
  63 {
  64   //
  65   // Default Destructor
  66   //
  67
  68 }
  69
  70 //______________________________________________
  71 void AliDielectronPair::SetTracks(AliVTrack * const particle1, Int_t pid1,
  72                                   AliVTrack * const particle2, Int_t pid2)
  73 {
  74   //
  75   // Sort particles by pt, first particle larget Pt
  76   // set AliKF daughters and pair
  77   //
  78   fPair.Initialize();
  79   fD1.Initialize();
  80   fD2.Initialize();
  81
  82   AliKFParticle kf1(*particle1,pid1);
  83   AliKFParticle kf2(*particle2,pid2);
  84
  85   fPair.AddDaughter(kf1);
  86   fPair.AddDaughter(kf2);
  87
  88   if (particle1->Pt()>particle2->Pt()){
  89     fRefD1 = particle1;
  90     fRefD2 = particle2;
  91     fD1+=kf1;
  92     fD2+=kf2;
  93   } else {
  94     fRefD1 = particle2;
  95     fRefD2 = particle1;
  96     fD1+=kf2;
  97     fD2+=kf1;
  98   }
  99 }
 100
 101 //______________________________________________
 102 void AliDielectronPair::GetThetaPhiCM(Double_t &thetaHE, Double_t &phiHE, Double_t &thetaCS, Double_t &phiCS) const
 103 {
 104   //
 105   // Calculate theta and phi in helicity and Collins-Soper coordinate frame
 106   //
 107   const Double_t kBeamEnergy   = 3500.;
 108   Double_t pxyz1[3]={0,0,0};
 109   Double_t pxyz2[3]={0,0,0};
 110   Double_t eleMass=AliPID::ParticleMass(AliPID::kElectron);
 111   Double_t proMass=AliPID::ParticleMass(AliPID::kProton);
 112
 113   AliVParticle *d1 = static_cast<AliVParticle*>(fRefD1.GetObject());
 114   AliVParticle *d2 = static_cast<AliVParticle*>(fRefD2.GetObject());
 115
 116   d1->PxPyPz(pxyz1);
 117   d2->PxPyPz(pxyz2);
 118
 119   TLorentzVector projMom(0.,0.,-kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
 120   TLorentzVector targMom(0.,0., kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
 121
 122   // first & second daughter 4-mom
 123   TLorentzVector p1Mom(pxyz1[0],pxyz1[1],pxyz1[2],
 124                        TMath::Sqrt(pxyz1[0]*pxyz1[0]+pxyz1[1]*pxyz1[1]+pxyz1[2]*pxyz1[2]+eleMass*eleMass));
 125   TLorentzVector p2Mom(pxyz2[0],pxyz2[1],pxyz2[2],
 126                        TMath::Sqrt(pxyz2[0]*pxyz2[0]+pxyz2[1]*pxyz2[1]+pxyz2[2]*pxyz2[2]+eleMass*eleMass));
 127   // J/Psi 4-momentum vector
 128   TLorentzVector motherMom=p1Mom+p2Mom;
 129
 130   // boost all the 4-mom vectors to the mother rest frame
 131   TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect();
 132   p1Mom.Boost(beta);
 133   p2Mom.Boost(beta);
 134   projMom.Boost(beta);
 135   targMom.Boost(beta);
 136
 137     // x,y,z axes
 138   TVector3 zAxisHE = (motherMom.Vect()).Unit();
 139   TVector3 zAxisCS = ((projMom.Vect()).Unit()-(targMom.Vect()).Unit()).Unit();
 140   TVector3 yAxis = ((projMom.Vect()).Cross(targMom.Vect())).Unit();
 141   TVector3 xAxisHE = (yAxis.Cross(zAxisHE)).Unit();
 142   TVector3 xAxisCS = (yAxis.Cross(zAxisCS)).Unit();
 143
 144   // fill theta and phi
 145   if(d1->Charge()>0){
 146     thetaHE = zAxisHE.Dot((p1Mom.Vect()).Unit());
 147     thetaCS = zAxisCS.Dot((p1Mom.Vect()).Unit());
 148     phiHE   = TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxisHE));
 149     phiCS   = TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxisCS));
 150   } else {
 151     thetaHE = zAxisHE.Dot((p2Mom.Vect()).Unit());
 152     thetaCS = zAxisCS.Dot((p2Mom.Vect()).Unit());
 153     phiHE   = TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxisHE));
 154     phiCS   = TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxisCS));
 155   }
 156 }
 157
 158 //______________________________________________
 159 Double_t AliDielectronPair::ThetaPhiCM(const AliVParticle* d1, const AliVParticle* d2,
 160                                        const Bool_t isHE, const Bool_t isTheta)
 161 {
 162   // The function calculates theta and phi in the mother rest frame with
 163   // respect to the helicity coordinate system and Collins-Soper coordinate system
 164   // TO DO: generalize for different decays (only J/Psi->e+e- now)
 165
 166   // Laboratory frame 4-vectors:
 167   // projectile beam & target beam 4-mom
 168   // TODO: need to retrieve the beam energy from somewhere
 169   const Double_t kBeamEnergy   = 3500.;
 170   Double_t px1=d1->Px();
 171   Double_t py1=d1->Py();
 172   Double_t pz1=d1->Pz();
 173   Double_t px2=d2->Px();
 174   Double_t py2=d2->Py();
 175   Double_t pz2=d2->Pz();
 176   Double_t eleMass=AliPID::ParticleMass(AliPID::kElectron);
 177   Double_t proMass=AliPID::ParticleMass(AliPID::kProton);
 178
 179   TLorentzVector projMom(0.,0.,-kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
 180   TLorentzVector targMom(0.,0., kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
 181
 182   // first & second daughter 4-mom
 183   TLorentzVector p1Mom(px1,py1,pz1,TMath::Sqrt(px1*px1+py1*py1+pz1*pz1+eleMass*eleMass));
 184   TLorentzVector p2Mom(px2,py2,pz2,TMath::Sqrt(px2*px2+py2*py2+pz2*pz2+eleMass*eleMass));
 185   // J/Psi 4-momentum vector
 186   TLorentzVector motherMom=p1Mom+p2Mom;
 187
 188   // boost all the 4-mom vectors to the mother rest frame
 189   TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect();
 190   p1Mom.Boost(beta);
 191   p2Mom.Boost(beta);
 192   projMom.Boost(beta);
 193   targMom.Boost(beta);
 194
 195   // x,y,z axes
 196   TVector3 zAxis;
 197   if(isHE) zAxis = (motherMom.Vect()).Unit();
 198   else zAxis = ((projMom.Vect()).Unit()-(targMom.Vect()).Unit()).Unit();
 199   TVector3 yAxis = ((projMom.Vect()).Cross(targMom.Vect())).Unit();
 200   TVector3 xAxis = (yAxis.Cross(zAxis)).Unit();
 201
 202   // return either theta or phi
 203   if(isTheta) {
 204     if(d1->Charge()>0)
 205       return zAxis.Dot((p1Mom.Vect()).Unit());
 206     else
 207       return zAxis.Dot((p2Mom.Vect()).Unit());
 208
 209   }
 210   else {
 211     if(d1->Charge()>0)
 212       return TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxis));
 213     else
 214       return TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxis));
 215   }
 216 }
 217
 218 //______________________________________________
 219 Double_t AliDielectronPair::ThetaPhiCM(const Bool_t isHE, const Bool_t isTheta) const {
 220   // The function calculates theta and phi in the mother rest frame with
 221   // respect to the helicity coordinate system and Collins-Soper coordinate system
 222   // TO DO: generalize for different decays (only J/Psi->e+e- now)
 223
 224   // Laboratory frame 4-vectors:
 225   // projectile beam & target beam 4-mom
 226   AliVParticle *d1 = dynamic_cast<AliVParticle*>(fRefD1.GetObject());
 227   AliVParticle *d2 = dynamic_cast<AliVParticle*>(fRefD2.GetObject());
 228
 229   const Double_t kBeamEnergy   = 3500.;
 230   Double_t px1=d1->Px();
 231   Double_t py1=d1->Py();
 232   Double_t pz1=d1->Pz();
 233   Double_t px2=d2->Px();
 234   Double_t py2=d2->Py();
 235   Double_t pz2=d2->Pz();
 236   Double_t eleMass=AliPID::ParticleMass(AliPID::kElectron);
 237   Double_t proMass=AliPID::ParticleMass(AliPID::kProton);
 238
 239   TLorentzVector projMom(0.,0.,-kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
 240   TLorentzVector targMom(0.,0., kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+proMass*proMass));
 241
 242   // first & second daughter 4-mom
 243   // first & second daughter 4-mom
 244   TLorentzVector p1Mom(px1,py1,pz1,TMath::Sqrt(px1*px1+py1*py1+pz1*pz1+eleMass*eleMass));
 245   TLorentzVector p2Mom(px2,py2,pz2,TMath::Sqrt(px2*px2+py2*py2+pz2*pz2+eleMass*eleMass));
 246   // J/Psi 4-momentum vector
 247   TLorentzVector motherMom=p1Mom+p2Mom;
 248
 249   // boost all the 4-mom vectors to the mother rest frame
 250   TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect();
 251   p1Mom.Boost(beta);
 252   p2Mom.Boost(beta);
 253   projMom.Boost(beta);
 254   targMom.Boost(beta);
 255
 256   // x,y,z axes
 257   TVector3 zAxis;
 258   if(isHE) zAxis = (motherMom.Vect()).Unit();
 259   else zAxis = ((projMom.Vect()).Unit()-(targMom.Vect()).Unit()).Unit();
 260   TVector3 yAxis = ((projMom.Vect()).Cross(targMom.Vect())).Unit();
 261   TVector3 xAxis = (yAxis.Cross(zAxis)).Unit();
 262
 263   // return either theta or phi
 264   if(isTheta) {
 265     if(fD1.GetQ()>0)
 266       return zAxis.Dot((p1Mom.Vect()).Unit());
 267     else
 268       return zAxis.Dot((p2Mom.Vect()).Unit());
 269   }
 270   else {
 271     if(fD1.GetQ()>0)
 272       return TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxis));
 273     else
 274       return TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxis));
 275   }
 276 }