MpdFieldMapSym3.cxx 5.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209
  1. // -------------------------------------------------------------------------
  2. // ----- MpdFieldMapSym3 source file -----
  3. // ----- Created 12/01/04 by M. Al/Turany (CbmField.cxx) -----
  4. // ----- Redesign 13/02/06 by V. Friese -----
  5. // -------------------------------------------------------------------------
  6. #include "MpdFieldMapSym3.h"
  7. #include "MpdFieldPar.h"
  8. #include "TArrayF.h"
  9. // ------------- Default constructor ----------------------------------
  10. MpdFieldMapSym3::MpdFieldMapSym3()
  11. : MpdFieldMap(),
  12. fHemiX(0.),
  13. fHemiY(0.),
  14. fHemiZ(0.)
  15. {
  16. fType = 3;
  17. }
  18. // ------------------------------------------------------------------------
  19. // ------------- Standard constructor ---------------------------------
  20. MpdFieldMapSym3::MpdFieldMapSym3(const char* mapName,
  21. const char* fileType)
  22. : MpdFieldMap(mapName, fileType),
  23. fHemiX(0.),
  24. fHemiY(0.),
  25. fHemiZ(0.)
  26. {
  27. fType = 3;
  28. }
  29. // ------------------------------------------------------------------------
  30. // ------------- Constructor from MpdFieldPar -------------------------
  31. MpdFieldMapSym3::MpdFieldMapSym3(MpdFieldPar* fieldPar)
  32. : MpdFieldMap(fieldPar),
  33. fHemiX(0.),
  34. fHemiY(0.),
  35. fHemiZ(0.)
  36. {
  37. fType = 3;
  38. }
  39. // ------------------------------------------------------------------------
  40. // ------------ Destructor --------------------------------------------
  41. MpdFieldMapSym3::~MpdFieldMapSym3() { }
  42. // ------------------------------------------------------------------------
  43. // ----------- Get x component of the field ---------------------------
  44. Double_t MpdFieldMapSym3::GetBx(Double_t x, Double_t y, Double_t z) {
  45. Int_t ix = 0;
  46. Int_t iy = 0;
  47. Int_t iz = 0;
  48. Double_t dx = 0.;
  49. Double_t dy = 0.;
  50. Double_t dz = 0.;
  51. if ( IsInside(x, y, z, ix, iy, iz, dx, dy, dz) ) {
  52. // Get Bx field values at grid cell corners
  53. fHa[0][0][0] = fBx->At(ix *fNy*fNz + iy *fNz + iz);
  54. fHa[1][0][0] = fBx->At((ix+1)*fNy*fNz + iy *fNz + iz);
  55. fHa[0][1][0] = fBx->At(ix *fNy*fNz + (iy+1)*fNz + iz);
  56. fHa[1][1][0] = fBx->At((ix+1)*fNy*fNz + (iy+1)*fNz + iz);
  57. fHa[0][0][1] = fBx->At(ix *fNy*fNz + iy *fNz + (iz+1));
  58. fHa[1][0][1] = fBx->At((ix+1)*fNy*fNz + iy *fNz + (iz+1));
  59. fHa[0][1][1] = fBx->At(ix *fNy*fNz + (iy+1)*fNz + (iz+1));
  60. fHa[1][1][1] = fBx->At((ix+1)*fNy*fNz + (iy+1)*fNz + (iz+1));
  61. // Return interpolated field value
  62. return Interpolate(dx, dy, dz) * fHemiX * fHemiY;
  63. }
  64. return 0.;
  65. }
  66. // ------------------------------------------------------------------------
  67. // ----------- Get y component of the field ---------------------------
  68. Double_t MpdFieldMapSym3::GetBy(Double_t x, Double_t y, Double_t z) {
  69. Int_t ix = 0;
  70. Int_t iy = 0;
  71. Int_t iz = 0;
  72. Double_t dx = 0.;
  73. Double_t dy = 0.;
  74. Double_t dz = 0.;
  75. if ( IsInside(x, y, z, ix, iy, iz, dx, dy, dz) ) {
  76. // Get By field values at grid cell corners
  77. fHa[0][0][0] = fBy->At(ix *fNy*fNz + iy *fNz + iz);
  78. fHa[1][0][0] = fBy->At((ix+1)*fNy*fNz + iy *fNz + iz);
  79. fHa[0][1][0] = fBy->At(ix *fNy*fNz + (iy+1)*fNz + iz);
  80. fHa[1][1][0] = fBy->At((ix+1)*fNy*fNz + (iy+1)*fNz + iz);
  81. fHa[0][0][1] = fBy->At(ix *fNy*fNz + iy *fNz + (iz+1));
  82. fHa[1][0][1] = fBy->At((ix+1)*fNy*fNz + iy *fNz + (iz+1));
  83. fHa[0][1][1] = fBy->At(ix *fNy*fNz + (iy+1)*fNz + (iz+1));
  84. fHa[1][1][1] = fBy->At((ix+1)*fNy*fNz + (iy+1)*fNz + (iz+1));
  85. // Return interpolated field value
  86. return Interpolate(dx, dy, dz);
  87. }
  88. return 0.;
  89. }
  90. // ------------------------------------------------------------------------
  91. // ----------- Get z component of the field ---------------------------
  92. Double_t MpdFieldMapSym3::GetBz(Double_t x, Double_t y, Double_t z) {
  93. Int_t ix = 0;
  94. Int_t iy = 0;
  95. Int_t iz = 0;
  96. Double_t dx = 0.;
  97. Double_t dy = 0.;
  98. Double_t dz = 0.;
  99. if ( IsInside(x, y, z, ix, iy, iz, dx, dy, dz) ) {
  100. // Get Bz field values at grid cell corners
  101. fHa[0][0][0] = fBz->At(ix *fNy*fNz + iy *fNz + iz);
  102. fHa[1][0][0] = fBz->At((ix+1)*fNy*fNz + iy *fNz + iz);
  103. fHa[0][1][0] = fBz->At(ix *fNy*fNz + (iy+1)*fNz + iz);
  104. fHa[1][1][0] = fBz->At((ix+1)*fNy*fNz + (iy+1)*fNz + iz);
  105. fHa[0][0][1] = fBz->At(ix *fNy*fNz + iy *fNz + (iz+1));
  106. fHa[1][0][1] = fBz->At((ix+1)*fNy*fNz + iy *fNz + (iz+1));
  107. fHa[0][1][1] = fBz->At(ix *fNy*fNz + (iy+1)*fNz + (iz+1));
  108. fHa[1][1][1] = fBz->At((ix+1)*fNy*fNz + (iy+1)*fNz + (iz+1));
  109. // Return interpolated field value
  110. return Interpolate(dx, dy, dz) * fHemiY * fHemiZ;
  111. }
  112. return 0.;
  113. }
  114. // ------------------------------------------------------------------------
  115. // ----------- Check whether a point is inside the map ----------------
  116. Bool_t MpdFieldMapSym3::IsInside(Double_t x, Double_t y, Double_t z,
  117. Int_t& ix, Int_t& iy, Int_t& iz,
  118. Double_t& dx, Double_t& dy,
  119. Double_t& dz) {
  120. // --- Transform into local coordinate system
  121. Double_t xl = x - fPosX;
  122. Double_t yl = y - fPosY;
  123. Double_t zl = z - fPosZ;
  124. // --- Reflect w.r.t. symmetry axes
  125. fHemiX = fHemiY = fHemiZ = 1.;
  126. if ( xl < 0. ) {
  127. fHemiX = -1.;
  128. xl = -1. * xl;
  129. }
  130. if ( yl < 0. ) {
  131. fHemiY = -1.;
  132. yl = -1. * yl;
  133. }
  134. if ( zl < 0. ) {
  135. fHemiZ = -1.;
  136. zl = -1. * zl;
  137. }
  138. // --- Check for being outside the map range
  139. if ( ! ( xl >= fXmin && xl < fXmax && yl >= fYmin && yl < fYmax &&
  140. zl >= fZmin && zl < fZmax ) ) {
  141. ix = iy = iz = 0;
  142. dx = dy = dz = 0.;
  143. return kFALSE;
  144. }
  145. // --- Determine grid cell
  146. ix = Int_t( (xl-fXmin) / fXstep );
  147. iy = Int_t( (yl-fYmin) / fYstep );
  148. iz = Int_t( (zl-fZmin) / fZstep );
  149. // Relative distance from grid point (in units of cell size)
  150. dx = (xl-fXmin) / fXstep - Double_t(ix);
  151. dy = (yl-fYmin) / fYstep - Double_t(iy);
  152. dz = (zl-fZmin) / fZstep - Double_t(iz);
  153. return kTRUE;
  154. }
  155. // ------------------------------------------------------------------------
  156. ClassImp(MpdFieldMapSym3)