MpdFemtoSmearPair.cxx 2.5 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586
  1. //
  2. // Class that returns a MpdFemtoPair with smeared momentum
  3. //
  4. // MpdFemtoMaker headers
  5. #include "MpdFemtoSmearPair.h"
  6. // ROOT headers
  7. #include "TMath.h"
  8. ClassImp(MpdFemtoSmearPair)
  9. //_________________
  10. MpdFemtoSmearPair::MpdFemtoSmearPair() :
  11. mSmearedPair(),
  12. mParticle1(),
  13. mParticle2(),
  14. mFracPtRes(0),
  15. mPhi_a(0),
  16. mPhi_b(0),
  17. mPhi_alpha(0),
  18. mTheta_a(0),
  19. mTheta_b(0),
  20. mTheta_alpha(0),
  21. mHbtRandom(nullptr) {
  22. // Default constructor
  23. setup();
  24. if (mHbtRandom) {
  25. delete mHbtRandom;
  26. mHbtRandom = nullptr;
  27. }
  28. mHbtRandom = new TRandom3(0);
  29. }
  30. //_________________
  31. MpdFemtoSmearPair::MpdFemtoSmearPair(const MpdFemtoPair* unSmearedPair) {
  32. setup();
  33. setUnsmearedPair(unSmearedPair);
  34. }
  35. //_________________
  36. void MpdFemtoSmearPair::setup() {
  37. mSmearedPair.setTrack1(&mParticle1);
  38. mSmearedPair.setTrack2(&mParticle2);
  39. if (!mHbtRandom) {
  40. mHbtRandom = new TRandom3(0);
  41. }
  42. }
  43. //_________________
  44. void MpdFemtoSmearPair::setUnsmearedPair(const MpdFemtoPair* unSmearedPair) {
  45. mParticle1.resetFourMomentum(smearedMomentum(unSmearedPair->track1()->fourMomentum()));
  46. mParticle2.resetFourMomentum(smearedMomentum(unSmearedPair->track2()->fourMomentum()));
  47. }
  48. /**
  49. * philosophy behind this smearing is as follows:
  50. * what we *measure* is pT (via curvature), and the angles phi and theta
  51. * momentum is related to this via
  52. * px = pT*cos(phi) --> Dpx = DpT*cos(phi) - pT*sin(phi)*Dphi = px*DpT/pT - py*Dphi
  53. * py = pT*sin(phi) --> Dpy = DpT*sin(phi) + pT*cos(phi)*Dphi = py*DpT/pT + px*Dphi
  54. * pz = pT/tan(theta) --> Dpz = DpT/tan(theta) - pT*Dtheta/(sin(theta))^2
  55. * = pT*(DpT/pT)/tan(theta) - pT*Dtheta/(sin(theta))^2
  56. * = pz*DpT/pT - pT*Dtheta/(sin(theta))^2
  57. */
  58. TLorentzVector MpdFemtoSmearPair::smearedMomentum(TLorentzVector fourmom) {
  59. double pT = fourmom.Perp();
  60. double mass2 = fourmom.M2();
  61. double px = fourmom.X();
  62. double py = fourmom.Y();
  63. double pz = fourmom.Z();
  64. double sin2theta = TMath::Sin(fourmom.Theta());
  65. sin2theta = sin2theta * sin2theta;
  66. double DpT_div_pT = mHbtRandom->Gaus(0., mFracPtRes);
  67. double Dphi = mHbtRandom->Gaus(0., (mPhi_a + mPhi_b * TMath::Power(pT, mPhi_alpha)));
  68. double Dtheta = mHbtRandom->Gaus(0., (mTheta_a + mTheta_b * TMath::Power(pT, mTheta_alpha)));
  69. fourmom.SetX(px * (1.0 + DpT_div_pT) - py * Dphi);
  70. fourmom.SetY(py * (1.0 + DpT_div_pT) + px * Dphi);
  71. fourmom.SetZ(pz * (1.0 + DpT_div_pT) - pT * Dtheta / sin2theta);
  72. fourmom.SetE(TMath::Sqrt(mass2 + fourmom.X() * fourmom.X() +
  73. fourmom.Y() * fourmom.Y() + fourmom.Z() * fourmom.Z()));
  74. return fourmom;
  75. }