MpdFemtoPhysicalHelix.cxx 3.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147
  1. /**
  2. * \class MpdFemtoPhysicalHelix
  3. * \author Grigory Nigmatkulov, May 07 2018
  4. *
  5. * Parametrization of a physical helix (modification of StPhysicalHelix).
  6. *
  7. */
  8. // C++ headers
  9. #include <math.h>
  10. // PicoDst headers
  11. #include "MpdFemtoHelix.h"
  12. #include "MpdFemtoPhysicalHelix.h"
  13. //#ifdef _VANILLA_ROOT_
  14. #include "PhysicalConstants.h"
  15. #include "SystemOfUnits.h"
  16. //#else
  17. //#include "StarClassLibrary/PhysicalConstants.h"
  18. //#include "StarClassLibrary/SystemOfUnits.h"
  19. //#endif
  20. ClassImp(MpdFemtoPhysicalHelix)
  21. //_________________
  22. MpdFemtoPhysicalHelix::MpdFemtoPhysicalHelix() {
  23. // Default constructor
  24. /* no-op */
  25. }
  26. //_________________
  27. MpdFemtoPhysicalHelix::~MpdFemtoPhysicalHelix() {
  28. // Destructor
  29. /* no-op */
  30. }
  31. //_________________
  32. MpdFemtoPhysicalHelix::MpdFemtoPhysicalHelix(const TVector3& p,
  33. const TVector3& o,
  34. double B, double q) {
  35. mH = (q * B <= 0) ? 1 : -1;
  36. if (p.y() == 0 && p.x() == 0) {
  37. setPhase((M_PI / 4)*(1 - 2. * mH));
  38. } else {
  39. setPhase(atan2(p.y(), p.x()) - mH * M_PI / 2);
  40. }
  41. setDipAngle(atan2(p.z(), p.Perp()));
  42. mOrigin = o;
  43. #ifndef ST_NO_NAMESPACES
  44. {
  45. using namespace units;
  46. #endif
  47. setCurvature(::fabs((c_light * nanosecond / meter * q * B / tesla) /
  48. (p.Mag() / GeV * mCosDipAngle) / meter));
  49. #ifndef ST_NO_NAMESPACES
  50. }
  51. #endif
  52. }
  53. //_________________
  54. MpdFemtoPhysicalHelix::MpdFemtoPhysicalHelix(double c, double d,
  55. double phase, const TVector3& o,
  56. int h) : MpdFemtoHelix(c, d, phase, o, h) {
  57. /* no-op */
  58. }
  59. //_________________
  60. TVector3 MpdFemtoPhysicalHelix::momentum(double B) const {
  61. if (mSingularity) {
  62. return (TVector3(0, 0, 0));
  63. } else {
  64. #ifndef ST_NO_NAMESPACES
  65. {
  66. using namespace units;
  67. #endif
  68. double pt = GeV * fabs(c_light * nanosecond / meter * B / tesla) / (fabs(mCurvature) * meter);
  69. return ( TVector3(pt * cos(mPhase + mH * M_PI / 2), // pos part pos field
  70. pt * sin(mPhase + mH * M_PI / 2),
  71. pt * tan(mDipAngle)));
  72. #ifndef ST_NO_NAMESPACES
  73. }
  74. #endif
  75. } //else
  76. }
  77. //_________________
  78. TVector3 MpdFemtoPhysicalHelix::momentumAt(double S, double B) const {
  79. // Obtain phase-shifted momentum from phase-shift of origin
  80. MpdFemtoPhysicalHelix tmp(*this);
  81. tmp.moveOrigin(S);
  82. return tmp.momentum(B);
  83. }
  84. //_________________
  85. int MpdFemtoPhysicalHelix::charge(double B) const {
  86. return (B > 0 ? -mH : mH);
  87. }
  88. //_________________
  89. double MpdFemtoPhysicalHelix::geometricSignedDistance(double x, double y) {
  90. // Geometric signed distance
  91. double thePath = this->pathLength(x, y);
  92. TVector3 DCA2dPosition = this->at(thePath);
  93. DCA2dPosition.SetZ(0);
  94. TVector3 position(x, y, 0);
  95. TVector3 DCAVec = (DCA2dPosition - position);
  96. TVector3 momVec;
  97. /// Deal with straight tracks
  98. if (this->mSingularity) {
  99. momVec = this->at(1) - this->at(0);
  100. momVec.SetZ(0);
  101. } else {
  102. momVec = this->momentumAt(thePath, 1. / tesla); // Don't care about Bmag. Helicity is what matters.
  103. momVec.SetZ(0);
  104. }
  105. double cross = DCAVec.x() * momVec.y() - DCAVec.y() * momVec.x();
  106. double theSign = (cross >= 0) ? 1. : -1.;
  107. return theSign * DCAVec.Perp();
  108. }
  109. //_________________
  110. double MpdFemtoPhysicalHelix::curvatureSignedDistance(double x, double y) {
  111. // Protect against mH = 0 or zero field
  112. if (this->mSingularity || abs(this->mH) <= 0) {
  113. return (this->geometricSignedDistance(x, y));
  114. } else {
  115. return (this->geometricSignedDistance(x, y)) / (this->mH);
  116. }
  117. }
  118. //_________________
  119. double MpdFemtoPhysicalHelix::geometricSignedDistance(const TVector3& pos) {
  120. double sdca2d = this->geometricSignedDistance(pos.x(), pos.y());
  121. double theSign = (sdca2d >= 0) ? 1. : -1.;
  122. return (this->distance(pos))*theSign;
  123. }
  124. //_________________
  125. double MpdFemtoPhysicalHelix::curvatureSignedDistance(const TVector3& pos) {
  126. double sdca2d = this->curvatureSignedDistance(pos.x(), pos.y());
  127. double theSign = (sdca2d >= 0) ? 1. : -1.;
  128. return (this->distance(pos))*theSign;
  129. }