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