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 #include <AliESDEvent.h>
  29
  30 #include "AliDielectronPair.h"
  31
  32 ClassImp(AliDielectronPair)
  33
  34 Double_t AliDielectronPair::fBeamEnergy=-1.;
  35
  36 AliDielectronPair::AliDielectronPair() :
  37   fType(-1),
  38   fLabel(-1),
  39   fPdgCode(0),
  40   fPair(),
  41   fD1(),
  42   fD2(),
  43   fRefD1(),
  44   fRefD2(),
  45   fKFUsage(kTRUE)
  46 {
  47   //
  48   // Default Constructor
  49   //
  50
  51 }
  52
  53 //______________________________________________
  54 AliDielectronPair::AliDielectronPair(AliVTrack * const particle1, Int_t pid1,
  55                                      AliVTrack * const particle2, Int_t pid2, Char_t type) :
  56   fType(type),
  57   fLabel(-1),
  58   fPdgCode(0),
  59   fPair(),
  60   fD1(),
  61   fD2(),
  62   fRefD1(),
  63   fRefD2(),
  64   fKFUsage(kTRUE)
  65 {
  66   //
  67   // Constructor with tracks
  68   //
  69   SetTracks(particle1, pid1, particle2, pid2);
  70 }
  71
  72 //______________________________________________
  73 AliDielectronPair::AliDielectronPair(const AliKFParticle * const particle1,
  74                                      const AliKFParticle * const particle2,
  75                                      AliVTrack * const refParticle1,
  76                                      AliVTrack * const refParticle2, Char_t type) :
  77   fType(type),
  78   fLabel(-1),
  79   fPdgCode(0),
  80   fPair(),
  81   fD1(),
  82   fD2(),
  83   fRefD1(),
  84   fRefD2(),
  85   fKFUsage(kTRUE)
  86 {
  87   //
  88   // Constructor with tracks
  89   //
  90   SetTracks(particle1, particle2,refParticle1,refParticle2);
  91 }
  92
  93 //______________________________________________
  94 AliDielectronPair::~AliDielectronPair()
  95 {
  96   //
  97   // Default Destructor
  98   //
  99
 100 }
 101
 102 //______________________________________________
 103 void AliDielectronPair::SetTracks(AliVTrack * const particle1, Int_t pid1,
 104                                   AliVTrack * const particle2, Int_t pid2)
 105 {
 106   //
 107   // Sort particles by pt, first particle larget Pt
 108   // set AliKF daughters and pair
 109   // refParticle1 and 2 are the original tracks. In the case of track rotation
 110   // they are needed in the framework
 111   //
 112   fPair.Initialize();
 113   fD1.Initialize();
 114   fD2.Initialize();
 115
 116   AliKFParticle kf1(*particle1,pid1);
 117   AliKFParticle kf2(*particle2,pid2);
 118
 119   fPair.AddDaughter(kf1);
 120   fPair.AddDaughter(kf2);
 121
 122   if (particle1->Pt()>particle2->Pt()){
 123     fRefD1 = particle1;
 124     fRefD2 = particle2;
 125     fD1+=kf1;
 126     fD2+=kf2;
 127   } else {
 128     fRefD1 = particle2;
 129     fRefD2 = particle1;
 130     fD1+=kf2;
 131     fD2+=kf1;
 132   }
 133 }
 134 //______________________________________________
 135 void AliDielectronPair::SetGammaTracks(AliVTrack * const particle1, Int_t pid1,
 136                                        AliVTrack * const particle2, Int_t pid2)
 137 {
 138   //
 139   // Sort particles by pt, first particle larget Pt
 140   // set AliKF daughters and a GAMMA pair
 141   // refParticle1 and 2 are the original tracks. In the case of track rotation
 142   // they are needed in the framework
 143   //
 144   fD1.Initialize();
 145   fD2.Initialize();
 146
 147   AliKFParticle kf1(*particle1,pid1);
 148   AliKFParticle kf2(*particle2,pid2);
 149   fPair.ConstructGamma(kf1,kf2);
 150
 151   if (particle1->Pt()>particle2->Pt()){
 152     fRefD1 = particle1;
 153     fRefD2 = particle2;
 154     fD1+=kf1;
 155     fD2+=kf2;
 156   } else {
 157     fRefD1 = particle2;
 158     fRefD2 = particle1;
 159     fD1+=kf2;
 160     fD2+=kf1;
 161   }
 162 }
 163
 164 //______________________________________________
 165 void AliDielectronPair::SetTracks(const AliKFParticle * const particle1,
 166                                   const AliKFParticle * const particle2,
 167                                   AliVTrack * const refParticle1,
 168                                   AliVTrack * const refParticle2)
 169 {
 170   //
 171   // Sort particles by pt, first particle larget Pt
 172   // set AliKF daughters and pair
 173   // refParticle1 and 2 are the original tracks. In the case of track rotation
 174   // they are needed in the framework
 175   //
 176   fPair.Initialize();
 177   fD1.Initialize();
 178   fD2.Initialize();
 179
 180   AliKFParticle kf1(*particle1);
 181   AliKFParticle kf2(*particle2);
 182
 183   fPair.AddDaughter(kf1);
 184   fPair.AddDaughter(kf2);
 185
 186   if (kf1.GetPt()>kf2.GetPt()){
 187     fRefD1 = refParticle1;
 188     fRefD2 = refParticle2;
 189     fD1+=kf1;
 190     fD2+=kf2;
 191   } else {
 192     fRefD1 = refParticle2;
 193     fRefD2 = refParticle1;
 194     fD1+=kf2;
 195     fD2+=kf1;
 196   }
 197 }
 198
 199 //______________________________________________
 200 void AliDielectronPair::GetThetaPhiCM(Double_t &thetaHE, Double_t &phiHE, Double_t &thetaCS, Double_t &phiCS) const
 201 {
 202   //
 203   // Calculate theta and phi in helicity and Collins-Soper coordinate frame
 204   //
 205   Double_t pxyz1[3]={fD1.GetPx(),fD1.GetPy(),fD1.GetPz()};
 206   Double_t pxyz2[3]={fD2.GetPx(),fD2.GetPy(),fD2.GetPz()};
 207   Double_t eleMass=AliPID::ParticleMass(AliPID::kElectron);
 208   Double_t proMass=AliPID::ParticleMass(AliPID::kProton);
 209
 210 //   AliVParticle *d1 = static_cast<AliVParticle*>(fRefD1.GetObject());
 211 //   AliVParticle *d2 = static_cast<AliVParticle*>(fRefD2.GetObject());
 212
 213 //   d1->PxPyPz(pxyz1);
 214 //   d2->PxPyPz(pxyz2);
 215
 216   TLorentzVector projMom(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+proMass*proMass));
 217   TLorentzVector targMom(0.,0., fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+proMass*proMass));
 218
 219   // first & second daughter 4-mom
 220   TLorentzVector p1Mom(pxyz1[0],pxyz1[1],pxyz1[2],
 221                        TMath::Sqrt(pxyz1[0]*pxyz1[0]+pxyz1[1]*pxyz1[1]+pxyz1[2]*pxyz1[2]+eleMass*eleMass));
 222   TLorentzVector p2Mom(pxyz2[0],pxyz2[1],pxyz2[2],
 223                        TMath::Sqrt(pxyz2[0]*pxyz2[0]+pxyz2[1]*pxyz2[1]+pxyz2[2]*pxyz2[2]+eleMass*eleMass));
 224   // J/Psi 4-momentum vector
 225   TLorentzVector motherMom=p1Mom+p2Mom;
 226
 227   // boost all the 4-mom vectors to the mother rest frame
 228   TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect();
 229   p1Mom.Boost(beta);
 230   p2Mom.Boost(beta);
 231   projMom.Boost(beta);
 232   targMom.Boost(beta);
 233
 234     // x,y,z axes
 235   TVector3 zAxisHE = (motherMom.Vect()).Unit();
 236   TVector3 zAxisCS = ((projMom.Vect()).Unit()-(targMom.Vect()).Unit()).Unit();
 237   TVector3 yAxis = ((projMom.Vect()).Cross(targMom.Vect())).Unit();
 238   TVector3 xAxisHE = (yAxis.Cross(zAxisHE)).Unit();
 239   TVector3 xAxisCS = (yAxis.Cross(zAxisCS)).Unit();
 240
 241   // fill theta and phi
 242   if(fD1.GetQ()>0){
 243     thetaHE = zAxisHE.Dot((p1Mom.Vect()).Unit());
 244     thetaCS = zAxisCS.Dot((p1Mom.Vect()).Unit());
 245     phiHE   = TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxisHE));
 246     phiCS   = TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxisCS));
 247   } else {
 248     thetaHE = zAxisHE.Dot((p2Mom.Vect()).Unit());
 249     thetaCS = zAxisCS.Dot((p2Mom.Vect()).Unit());
 250     phiHE   = TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxisHE));
 251     phiCS   = TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxisCS));
 252   }
 253 }
 254
 255
 256 //______________________________________________
 257 Double_t AliDielectronPair::PsiPair(Double_t MagField) const
 258 {
 259   //Following idea to use opening of colinear pairs in magnetic field from e.g. PHENIX
 260   //to ID conversions. Adapted from AliTRDv0Info class
 261   Double_t x, y;//, z;
 262   x = fPair.GetX();
 263   y = fPair.GetY();
 264   //  z = fPair.GetZ();
 265
 266   Double_t m1[3] = {0,0,0};
 267   Double_t m2[3] = {0,0,0};
 268
 269   m1[0] = fD1.GetPx();
 270   m1[1] = fD1.GetPy();
 271   m1[2] = fD1.GetPz();
 272
 273   m2[0] = fD2.GetPx();
 274   m2[1] = fD2.GetPy();
 275   m2[2] = fD2.GetPz();
 276
 277   Double_t deltat = 1.;
 278   deltat = TMath::ATan(m2[2]/(TMath::Sqrt(m2[0]*m2[0] + m2[1]*m2[1])+1.e-13))-
 279         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
 280
 281   Double_t radiussum = TMath::Sqrt(x*x + y*y) + 50;//radius to which tracks shall be propagated
 282
 283   Double_t mom1Prop[3];
 284   Double_t mom2Prop[3];
 285
 286   AliExternalTrackParam *d1 = static_cast<AliExternalTrackParam*>(fRefD1.GetObject());
 287   AliExternalTrackParam *d2 = static_cast<AliExternalTrackParam*>(fRefD2.GetObject());
 288
 289   AliExternalTrackParam nt(*d1), pt(*d2);
 290
 291   Double_t fPsiPair = 4.;
 292   if(nt.PropagateTo(radiussum,MagField) == 0)//propagate tracks to the outside
 293         fPsiPair =  -5.;
 294   if(pt.PropagateTo(radiussum,MagField) == 0)
 295         fPsiPair = -5.;
 296   pt.GetPxPyPz(mom1Prop);//Get momentum vectors of tracks after propagation
 297   nt.GetPxPyPz(mom2Prop);
 298
 299
 300
 301   Double_t pEle =
 302         TMath::Sqrt(mom2Prop[0]*mom2Prop[0]+mom2Prop[1]*mom2Prop[1]+mom2Prop[2]*mom2Prop[2]);//absolute momentum val
 303   Double_t pPos =
 304         TMath::Sqrt(mom1Prop[0]*mom1Prop[0]+mom1Prop[1]*mom1Prop[1]+mom1Prop[2]*mom1Prop[2]);//absolute momentum val
 305
 306   Double_t scalarproduct =
 307         mom1Prop[0]*mom2Prop[0]+mom1Prop[1]*mom2Prop[1]+mom1Prop[2]*mom2Prop[2];//scalar product of propagated posit
 308
 309   Double_t chipair = TMath::ACos(scalarproduct/(pEle*pPos));//Angle between propagated daughter tracks
 310
 311   fPsiPair =  TMath::Abs(TMath::ASin(deltat/chipair));
 312
 313   return fPsiPair;
 314
 315 }
 316
 317 //______________________________________________
 318 Double_t AliDielectronPair::ThetaPhiCM(const AliVParticle* d1, const AliVParticle* d2,
 319                                        Bool_t isHE, Bool_t isTheta)
 320 {
 321   // The function calculates theta and phi in the mother rest frame with
 322   // respect to the helicity coordinate system and Collins-Soper coordinate system
 323   // TO DO: generalize for different decays (only J/Psi->e+e- now)
 324
 325   // Laboratory frame 4-vectors:
 326   // projectile beam & target beam 4-mom
 327   Double_t px1=d1->Px();
 328   Double_t py1=d1->Py();
 329   Double_t pz1=d1->Pz();
 330   Double_t px2=d2->Px();
 331   Double_t py2=d2->Py();
 332   Double_t pz2=d2->Pz();
 333   Double_t eleMass=AliPID::ParticleMass(AliPID::kElectron);
 334   Double_t proMass=AliPID::ParticleMass(AliPID::kProton);
 335 //   printf(" beam energy %f \n ", fBeamEnergy);
 336   TLorentzVector projMom(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+proMass*proMass));
 337   TLorentzVector targMom(0.,0., fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+proMass*proMass));
 338
 339   // first & second daughter 4-mom
 340   TLorentzVector p1Mom(px1,py1,pz1,TMath::Sqrt(px1*px1+py1*py1+pz1*pz1+eleMass*eleMass));
 341   TLorentzVector p2Mom(px2,py2,pz2,TMath::Sqrt(px2*px2+py2*py2+pz2*pz2+eleMass*eleMass));
 342   // J/Psi 4-momentum vector
 343   TLorentzVector motherMom=p1Mom+p2Mom;
 344
 345   // boost all the 4-mom vectors to the mother rest frame
 346   TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect();
 347   p1Mom.Boost(beta);
 348   p2Mom.Boost(beta);
 349   projMom.Boost(beta);
 350   targMom.Boost(beta);
 351
 352   // x,y,z axes
 353   TVector3 zAxis;
 354   if(isHE) zAxis = (motherMom.Vect()).Unit();
 355   else zAxis = ((projMom.Vect()).Unit()-(targMom.Vect()).Unit()).Unit();
 356   TVector3 yAxis = ((projMom.Vect()).Cross(targMom.Vect())).Unit();
 357   TVector3 xAxis = (yAxis.Cross(zAxis)).Unit();
 358
 359   // return either theta or phi
 360   if(isTheta) {
 361     if(d1->Charge()>0)
 362       return zAxis.Dot((p1Mom.Vect()).Unit());
 363     else
 364       return zAxis.Dot((p2Mom.Vect()).Unit());
 365
 366   }
 367   else {
 368     if(d1->Charge()>0)
 369       return TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxis));
 370     else
 371       return TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxis));
 372   }
 373 }
 374
 375 //______________________________________________
 376 Double_t AliDielectronPair::ThetaPhiCM(Bool_t isHE, Bool_t isTheta) const {
 377   // The function calculates theta and phi in the mother rest frame with
 378   // respect to the helicity coordinate system and Collins-Soper coordinate system
 379   // TO DO: generalize for different decays (only J/Psi->e+e- now)
 380
 381   // Laboratory frame 4-vectors:
 382   // projectile beam & target beam 4-mom
 383   AliVParticle *d1 = static_cast<AliVParticle*>(fRefD1.GetObject());
 384   AliVParticle *d2 = static_cast<AliVParticle*>(fRefD2.GetObject());
 385
 386   Double_t px1=d1->Px();
 387   Double_t py1=d1->Py();
 388   Double_t pz1=d1->Pz();
 389   Double_t px2=d2->Px();
 390   Double_t py2=d2->Py();
 391   Double_t pz2=d2->Pz();
 392   Double_t eleMass=AliPID::ParticleMass(AliPID::kElectron);
 393   Double_t proMass=AliPID::ParticleMass(AliPID::kProton);
 394
 395   TLorentzVector projMom(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+proMass*proMass));
 396   TLorentzVector targMom(0.,0., fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+proMass*proMass));
 397
 398   // first & second daughter 4-mom
 399   // first & second daughter 4-mom
 400   TLorentzVector p1Mom(px1,py1,pz1,TMath::Sqrt(px1*px1+py1*py1+pz1*pz1+eleMass*eleMass));
 401   TLorentzVector p2Mom(px2,py2,pz2,TMath::Sqrt(px2*px2+py2*py2+pz2*pz2+eleMass*eleMass));
 402   // J/Psi 4-momentum vector
 403   TLorentzVector motherMom=p1Mom+p2Mom;
 404
 405   // boost all the 4-mom vectors to the mother rest frame
 406   TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect();
 407   p1Mom.Boost(beta);
 408   p2Mom.Boost(beta);
 409   projMom.Boost(beta);
 410   targMom.Boost(beta);
 411
 412   // x,y,z axes
 413   TVector3 zAxis;
 414   if(isHE) zAxis = (motherMom.Vect()).Unit();
 415   else zAxis = ((projMom.Vect()).Unit()-(targMom.Vect()).Unit()).Unit();
 416   TVector3 yAxis = ((projMom.Vect()).Cross(targMom.Vect())).Unit();
 417   TVector3 xAxis = (yAxis.Cross(zAxis)).Unit();
 418
 419   // return either theta or phi
 420   if(isTheta) {
 421     if(fD1.GetQ()>0)
 422       return zAxis.Dot((p1Mom.Vect()).Unit());
 423     else
 424       return zAxis.Dot((p2Mom.Vect()).Unit());
 425   }
 426   else {
 427     if(fD1.GetQ()>0)
 428       return TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxis));
 429     else
 430       return TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxis));
 431   }
 432 }
 433 //______________________________________________
 434 Double_t AliDielectronPair::GetCosPointingAngle(const AliVVertex *primVtx) const
 435 {
 436   //
 437   // Calculate the poiting angle of the pair to the primary vertex and take the cosine
 438   //
 439   if(!primVtx) return -1.;
 440
 441   Double_t deltaPos[3]; //vector between the reference point and the V0 vertex
 442   deltaPos[0] = fPair.GetX() - primVtx->GetX();
 443   deltaPos[1] = fPair.GetY() - primVtx->GetY();
 444   deltaPos[2] = fPair.GetZ() - primVtx->GetZ();
 445
 446   Double_t momV02    = Px()*Px() + Py()*Py() + Pz()*Pz();
 447   Double_t deltaPos2 = deltaPos[0]*deltaPos[0] + deltaPos[1]*deltaPos[1] + deltaPos[2]*deltaPos[2];
 448
 449   Double_t cosinePointingAngle = (deltaPos[0]*Px() + deltaPos[1]*Py() + deltaPos[2]*Pz()) / TMath::Sqrt(momV02 * deltaPos2);
 450
 451   return TMath::Abs(cosinePointingAngle);
 452
 453 }
 454
 455 //______________________________________________
 456 Double_t AliDielectronPair::GetArmAlpha() const
 457 {
 458   //
 459   // Calculate the Armenteros-Podolanski Alpha
 460   //
 461   Int_t qD1 = fD1.GetQ();
 462
 463   TVector3 momNeg( (qD1<0?fD1.GetPx():fD2.GetPx()),
 464                    (qD1<0?fD1.GetPy():fD2.GetPy()),
 465                    (qD1<0?fD1.GetPz():fD2.GetPz()) );
 466   TVector3 momPos( (qD1<0?fD2.GetPx():fD1.GetPx()),
 467                    (qD1<0?fD2.GetPy():fD1.GetPy()),
 468                    (qD1<0?fD2.GetPz():fD1.GetPz()) );
 469   TVector3 momTot(Px(),Py(),Pz());
 470
 471   Double_t lQlNeg = momNeg.Dot(momTot)/momTot.Mag();
 472   Double_t lQlPos = momPos.Dot(momTot)/momTot.Mag();
 473
 474   return ((lQlPos - lQlNeg)/(lQlPos + lQlNeg));
 475 }
 476
 477 //______________________________________________
 478 Double_t AliDielectronPair::GetArmPt() const
 479 {
 480   //
 481   // Calculate the Armenteros-Podolanski Pt
 482   //
 483   Int_t qD1 = fD1.GetQ();
 484
 485   TVector3 momNeg( (qD1<0?fD1.GetPx():fD2.GetPx()),
 486                    (qD1<0?fD1.GetPy():fD2.GetPy()),
 487                    (qD1<0?fD1.GetPz():fD2.GetPz()) );
 488   TVector3 momTot(Px(),Py(),Pz());
 489
 490   return (momNeg.Perp(momTot));
 491 }
 492
 493 //______________________________________________
 494 void AliDielectronPair::GetDCA(const AliVVertex *primVtx, Double_t d0z0[2]) const
 495 {
 496   //
 497   // Calculate the dca of the mother with respect to the primary vertex
 498   //
 499   if(!primVtx) return;
 500
 501   d0z0[0] = TMath::Sqrt(TMath::Power(Xv()-primVtx->GetX(),2) +
 502                         TMath::Power(Yv()-primVtx->GetY(),2) );
 503
 504   d0z0[1] = Zv() - primVtx->GetZ();
 505   return;
 506
 507 }
 508
 509 // //______________________________________________
 510 // Double_t AliDielectronPair::GetLXY(const AliVVertex * const vtx) const
 511 // {
 512 //   //
 513 //   // Calculate the decay length in XY taking into account the primary vertex position
 514 //   //
 515 //   if(!vtx) return 0;
 516 //   return ( (Xv()-vtx->GetX()) * Px() + (Yv()-vtx->GetY()) * Py() )/Pt()  ;
 517 // }
 518
 519 // //______________________________________________
 520 // Double_t AliDielectronPair::GetPseudoProperTime(const AliVVertex * const vtx) const
 521 // {
 522 //   //
 523 //   // Calculate the pseudo proper time
 524 //   //
 525 //   Double_t lxy=GetLXY(vtx);
 526 //   Double_t psProperDecayLength = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/Pt();
 527 //   return psProperDecayLength;
 528 // }
 529
 530
 531 //______________________________________________
 532 Double_t AliDielectronPair::PhivPair(Double_t MagField) const
 533 {
 534   //Following idea to use opening of colinear pairs in magnetic field from e.g. PHENIX
 535   //to ID conversions. Angle between ee plane and magnetic field is calculated.
 536
 537   //Define local buffer variables for leg properties
 538   Double_t px1=-9999.,py1=-9999.,pz1=-9999.;
 539   Double_t px2=-9999.,py2=-9999.,pz2=-9999.;
 540
 541   if(MagField>0){
 542     if(fD1.GetQ()>0){
 543       px1 = fD1.GetPx();
 544       py1 = fD1.GetPy();
 545       pz1 = fD1.GetPz();
 546
 547       px2 = fD2.GetPx();
 548       py2 = fD2.GetPy();
 549       pz2 = fD2.GetPz();
 550     }else{
 551       px1 = fD2.GetPx();
 552       py1 = fD2.GetPy();
 553       pz1 = fD2.GetPz();
 554
 555       px2 = fD1.GetPx();
 556       py2 = fD1.GetPy();
 557       pz2 = fD1.GetPz();
 558     }
 559   }else{
 560     if(fD1.GetQ()>0){
 561       px1 = fD2.GetPx();
 562       py1 = fD2.GetPy();
 563       pz1 = fD2.GetPz();
 564
 565       px2 = fD1.GetPx();
 566       py2 = fD1.GetPy();
 567       pz2 = fD1.GetPz();
 568     }else{
 569       px1 = fD1.GetPx();
 570       py1 = fD1.GetPy();
 571       pz1 = fD1.GetPz();
 572
 573       px2 = fD2.GetPx();
 574       py2 = fD2.GetPy();
 575       pz2 = fD2.GetPz();
 576     }
 577   }
 578
 579   Double_t px = px1+px2;
 580   Double_t py = py1+py2;
 581   Double_t pz = pz1+pz2;
 582   Double_t dppair = TMath::Sqrt(px*px+py*py+pz*pz);
 583
 584   //unit vector of (pep+pem)
 585   Double_t pl = dppair;
 586   Double_t ux = px/pl;
 587   Double_t uy = py/pl;
 588   Double_t uz = pz/pl;
 589   Double_t ax = uy/TMath::Sqrt(ux*ux+uy*uy);
 590   Double_t ay = -ux/TMath::Sqrt(ux*ux+uy*uy);
 591
 592   //momentum of e+ and e- in (ax,ay,az) axis. Note that az=0 by
 593   //definition.
 594   //Double_t ptep = iep->Px()*ax + iep->Py()*ay;
 595   //Double_t ptem = iem->Px()*ax + iem->Py()*ay;
 596
 597   Double_t pxep = px1;
 598   Double_t pyep = py1;
 599   Double_t pzep = pz1;
 600   Double_t pxem = px2;
 601   Double_t pyem = py2;
 602   Double_t pzem = pz2;
 603
 604   //vector product of pep X pem
 605   Double_t vpx = pyep*pzem - pzep*pyem;
 606   Double_t vpy = pzep*pxem - pxep*pzem;
 607   Double_t vpz = pxep*pyem - pyep*pxem;
 608   Double_t vp = sqrt(vpx*vpx+vpy*vpy+vpz*vpz);
 609   //Double_t thev = acos(vpz/vp);
 610
 611   //unit vector of pep X pem
 612   Double_t vx = vpx/vp;
 613   Double_t vy = vpy/vp;
 614   Double_t vz = vpz/vp;
 615
 616   //The third axis defined by vector product (ux,uy,uz)X(vx,vy,vz)
 617   Double_t wx = uy*vz - uz*vy;
 618   Double_t wy = uz*vx - ux*vz;
 619   //Double_t wz = ux*vy - uy*vx;
 620   //Double_t wl = sqrt(wx*wx+wy*wy+wz*wz);
 621   // by construction, (wx,wy,wz) must be a unit vector.
 622   // measure angle between (wx,wy,wz) and (ax,ay,0). The angle between them
 623   // should be small if the pair is conversion
 624   //
 625   Double_t cosPhiV = wx*ax + wy*ay;
 626   Double_t phiv = TMath::ACos(cosPhiV);
 627
 628   return phiv;
 629
 630 }
 631
 632 //______________________________________________
 633 Double_t AliDielectronPair::GetPairPlaneAngle(Double_t v0rpH2, Int_t VariNum)const
 634 {
 635
 636   // Calculate the angle between electron pair plane and variables
 637   // kv0rpH2 is reaction plane angle using V0-A,C,AC,Random
 638
 639   Double_t px1=-9999.,py1=-9999.,pz1=-9999.;
 640   Double_t px2=-9999.,py2=-9999.,pz2=-9999.;
 641
 642   px1 = fD1.GetPx();
 643   py1 = fD1.GetPy();
 644   pz1 = fD1.GetPz();
 645
 646   px2 = fD2.GetPx();
 647   py2 = fD2.GetPy();
 648   pz2 = fD2.GetPz();
 649
 650   //p1+p2
 651   Double_t px = px1+px2;
 652   Double_t py = py1+py2;
 653   Double_t pz = pz1+pz2;
 654
 655   // normal vector of ee plane
 656   Double_t pnorx = py1*pz2 - pz1*py2;
 657   Double_t pnory = pz1*px2 - px1*pz2;
 658   Double_t pnorz = px1*py2 - py1*px2;
 659   Double_t pnor  = TMath::Sqrt( pnorx*pnorx + pnory*pnory + pnorz*pnorz );
 660
 661   //unit vector
 662   Double_t upnx = -9999.;
 663   Double_t upny = -9999.;
 664   Double_t upnz = -9999.;
 665
 666   if (pnor !=0)
 667     {
 668       upnx= pnorx/pnor;
 669       upny= pnory/pnor;
 670       upnz= pnorz/pnor;
 671     }
 672
 673
 674   Double_t ax = -9999.;
 675   Double_t ay = -9999.;
 676   Double_t az = -9999.;
 677
 678   //variable 1
 679   //seeing the angle between ee decay plane and reaction plane by using V0-A,C,AC,Random
 680           if(VariNum == 1){
 681                 ax = TMath::Sin(v0rpH2);
 682                 ay = -TMath::Cos(v0rpH2);
 683                 az = 0.0;
 684           }
 685
 686
 687         //variable 2
 688         //seeing the angle between ee decay plane and (p1+p2) rot ez
 689           else if (VariNum == 2 ){
 690                 ax = py;
 691                 ay = -px;
 692                 az = 0.0;
 693           }
 694
 695         //variable 3
 696         //seeing the angle between ee decay plane and (p1+p2) rot (p1+p2)x'z
 697           else if (VariNum == 3 ){
 698                 Double_t rotpx = px*TMath::Cos(v0rpH2)+py*TMath::Sin(v0rpH2);
 699                 //Double_t rotpy = 0.0;
 700                 // Double_t rotpz = pz;
 701
 702                 ax = py*pz;
 703                 ay = pz*rotpx-pz*px;
 704                 az = -rotpx*py;
 705           }
 706
 707         //variable 4
 708         //seeing the angle between ee decay plane and (p1+p2) rot ey'
 709           else if (VariNum == 4){
 710                 ax = 0.0;
 711                 ay = 0.0;
 712                 az = pz;
 713           }
 714
 715         Double_t denomHelper = ax*ax + ay*ay +az*az;
 716         Double_t uax = -9999.;
 717         Double_t uay = -9999.;
 718         Double_t uaz = -9999.;
 719
 720         if (denomHelper !=0) {
 721           uax = ax/TMath::Sqrt(denomHelper);
 722           uay = ay/TMath::Sqrt(denomHelper);
 723           uaz = az/TMath::Sqrt(denomHelper);
 724         }
 725
 726         //PM is the angle between Pair plane and a plane decided by using variable 1-4
 727
 728         Double_t cosPM = upnx*uax + upny*uay + upnz*uaz;
 729         Double_t PM = TMath::ACos(cosPM);
 730
 731         //keep interval [0,pi/2]
 732         if(PM > TMath::Pi()/2){
 733           PM -= TMath::Pi();
 734           PM *= -1.0;
 735
 736         }
 737         return PM;
 738 }
 739
 740
 741 //_______________________________________________
 742 Double_t AliDielectronPair::PairPlaneMagInnerProduct(Double_t ZDCrpH1) const
 743 {
 744
 745   // Calculate inner product of the strong magnetic field and electron pair plane
 746
 747   if(ZDCrpH1 == 0.) return -9999.;
 748
 749   Double_t px1=-9999.,py1=-9999.,pz1=-9999.;
 750   Double_t px2=-9999.,py2=-9999.,pz2=-9999.;
 751
 752   px1 = fD1.GetPx();
 753   py1 = fD1.GetPy();
 754   pz1 = fD1.GetPz();
 755
 756
 757   px2 = fD2.GetPx();
 758   py2 = fD2.GetPy();
 759   pz2 = fD2.GetPz();
 760
 761   // normal vector of ee plane
 762   Double_t pnorx = py2*pz1 - pz2*py1;
 763   Double_t pnory = pz2*px1 - px2*pz1;
 764   Double_t pnorz = px2*py1 - py2*px1;
 765   Double_t pnor  = TMath::Sqrt( pnorx*pnorx + pnory*pnory + pnorz*pnorz );
 766
 767   //unit vector
 768   Double_t upnx = -9999.;
 769   Double_t upny = -9999.;
 770   //Double_t upnz = -9999.;
 771
 772   if (pnor == 0) return -9999.;
 773   upnx= pnorx/pnor;
 774   upny= pnory/pnor;
 775   //upnz= pnorz/pnor;
 776
 777   //direction of strong magnetic field
 778   Double_t magx = TMath::Cos(ZDCrpH1+(TMath::Pi()/2));
 779   Double_t magy = TMath::Sin(ZDCrpH1+(TMath::Pi()/2));
 780
 781   //inner product of strong magnetic field and  ee plane
 782   Double_t upnmag = upnx*magx + upny*magy;
 783
 784   return upnmag;
 785 }
 786
 787
 788 //______________________________________________
 789 void AliDielectronPair::SetBeamEnergy(AliVEvent *ev, Double_t beamEbyHand)
 790 {
 791   //
 792   // set the beam energy (by hand in case of AODs)
 793   //
 794   if(ev->IsA()==AliESDEvent::Class())
 795     fBeamEnergy = ((AliESDEvent*)ev)->GetBeamEnergy();
 796   else
 797     fBeamEnergy = beamEbyHand;
 798 }
 799
 800