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