MpdFieldMap.cxx 18 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562
  1. // -------------------------------------------------------------------------
  2. // MpdFieldMap source file -----
  3. // Created 23/07/13 by P. Batyuk (MPD) <batyuk@jinr.ru> -----
  4. // from PndFieldMap (PNDROOT) -----
  5. // -------------------------------------------------------------------------
  6. /// Last modified: 24.07.2013, P.B.
  7. /// Last modified: 30.07.2013, P.B.
  8. /// Last modified: 31.07.2013, P.B.
  9. #include <iomanip>
  10. #include <iostream>
  11. #include <fstream>
  12. #include "stdlib.h"
  13. #include "TArrayF.h"
  14. #include "TFile.h"
  15. #include "TMath.h"
  16. #include "MpdFieldMap.h"
  17. #include "MpdFieldMapData.h"
  18. #include "MpdFieldPar.h"
  19. using namespace std;
  20. // ------------- Default constructor ----------------------------------
  21. MpdFieldMap::MpdFieldMap() {
  22. fPosX = fPosY = fPosZ = 0.;
  23. fXmin = fYmin = fZmin = 0.;
  24. fXmax = fYmax = fZmax = 0.;
  25. fXstep = fYstep = fZstep = 0.;
  26. fNx = fNy = fNz = 0;
  27. fScale = 1.;
  28. funit = 10.0;
  29. fBx = fBy = fBz = NULL;
  30. fPosX = fPosY = fPosZ = 0.;
  31. fName = "";
  32. fFileName = "";
  33. fType = 1;
  34. }
  35. // ------------- Standard constructor ---------------------------------
  36. MpdFieldMap::MpdFieldMap(const char* mapName, const char* fileType)
  37. : FairField(mapName) {
  38. fPosX = fPosY = fPosZ = 0.;
  39. fXmin = fYmin = fZmin = 0.;
  40. fXmax = fYmax = fZmax = 0.;
  41. fXstep = fYstep = fZstep = 0.;
  42. fNx = fNy = fNz = 0;
  43. fScale = 1.;
  44. funit = 10.0;
  45. fBx = fBy = fBz = NULL;
  46. fName = mapName;
  47. TString dir = getenv("VMCWORKDIR");
  48. fFileName = dir + "/input/" + mapName;
  49. if ( fileType[0] == 'R' ) fFileName += ".root";
  50. else fFileName += ".dat";
  51. fType = 1;
  52. }
  53. // ------------ Constructor from MpdFieldPar --------------------------
  54. MpdFieldMap::MpdFieldMap(MpdFieldPar* fieldPar) {
  55. fType = 1;
  56. fPosX = fPosY = fPosZ = 0.;
  57. fXmin = fYmin = fZmin = 0.;
  58. fXmax = fYmax = fZmax = 0.;
  59. fXstep = fYstep = fZstep = 0.;
  60. fNx = fNy = fNz = 0;
  61. fScale = 1.;
  62. funit = 10.0;
  63. fBx = fBy = fBz = NULL;
  64. if ( ! fieldPar ) {
  65. cerr << "-W- MpdConstField::MpdConstField: empty parameter container!"
  66. << endl;
  67. fName = "";
  68. fFileName = "";
  69. fType = -1;
  70. }
  71. else {
  72. fieldPar->MapName(fName);
  73. fPosX = fieldPar->GetPositionX();
  74. fPosY = fieldPar->GetPositionY();
  75. fPosZ = fieldPar->GetPositionZ();
  76. fScale = fieldPar->GetScale();
  77. TString dir = getenv("VMCWORKDIR");
  78. fFileName = dir + "/input/" + fName + ".root";
  79. fType = fieldPar->GetType();
  80. }
  81. }
  82. // ------------ Destructor --------------------------------------------
  83. MpdFieldMap::~MpdFieldMap() {
  84. if ( fBx ) delete fBx;
  85. if ( fBy ) delete fBy;
  86. if ( fBz ) delete fBz;
  87. }
  88. // ----------- Intialization ------------------------------------------
  89. void MpdFieldMap::Init() {
  90. if (fFileName.EndsWith(".root")) ReadRootFile(fFileName, fName);
  91. else if (fFileName.EndsWith(".dat")) ReadAsciiFile(fFileName);
  92. else {
  93. cerr << "-E- MpdFieldMap::Init: No proper file name defined! ("
  94. << fFileName << ")" << endl;
  95. Fatal("Init", "No proper file name");
  96. }
  97. }
  98. // ----------- Check whether a point is inside the map ----------------
  99. Bool_t MpdFieldMap::IsInside(Double_t x, Double_t y, Double_t z) {
  100. /// Define geometry sizes of TPC to interpolate mag. field....
  101. /// Extrapolation used to define mag. field components in vicinity of the TPC-volume
  102. /// r1extrap <= R <= r2extrap
  103. /// r1 - 2cm <= R <= r2 + 2 cm
  104. Double_t delta_r(2.);
  105. Double_t r1(40.3), r2(120.3); /// Real internal and external radiuses of the TPC
  106. //Double_t fZmin(-170.), fZmax(170.);
  107. //Double_t fXmin_r2 = -r2, fXmax_r2 = r2;
  108. //Double_t fYmin_r2 = fXmin_r2, fYmax_r2 = fXmax_r2;
  109. //Double_t fXmin_r1 = -r1, fXmax_r1 = r1;
  110. //Double_t fYmin_r1 = fXmin_r1, fYmax_r1 = fXmax_r1;
  111. /// Checking only correct field map initialization
  112. /// In fact, the mag. field used at the point (0, 0, 0) is considered to be (0, 0, 5) kG
  113. if (x == 0 && y == 0 && z== 0) return kTRUE;
  114. /// Checking a point to be inside the TPC (B-field Map) (r1 - delta_r = 40.3 cm - 2 cm, r2 + delta_r = 120.3 cm + 2 cm)
  115. else if ( TMath::Power(x, 2) + TMath::Power(y, 2) >= TMath::Power(r1 - delta_r, 2) &&
  116. TMath::Power(x, 2) + TMath::Power(y, 2) <= TMath::Power(r2 + delta_r, 2)) return kTRUE;
  117. else return kFALSE;
  118. }
  119. // ---------- Write the map to an ASCII file --------------------------
  120. void MpdFieldMap::WriteAsciiFile(const char* fileName) {
  121. // Open file
  122. cout << "-I- MpdFieldMap: Writing field map to ASCII file "
  123. << fileName << endl;
  124. ofstream mapFile(fileName);
  125. if ( ! mapFile.is_open() ) {
  126. cerr << "-E- MpdFieldMap:ReadAsciiFile: Could not open file! " << endl;
  127. return;
  128. }
  129. // Write field map grid parameters
  130. mapFile.precision(4);
  131. mapFile << showpoint;
  132. if ( fType == 1 ) mapFile << "nosym" << endl;
  133. if ( funit == 10.0 ) mapFile << "T" << endl;
  134. else if ( funit == 0.001 ) mapFile << "G" << endl;
  135. else if ( funit == 1.0 ) mapFile << "kG" << endl;
  136. mapFile << fXmin << " " << fXmax << " " << fNx << endl;
  137. mapFile << fYmin << " " << fYmax << " " << fNy << endl;
  138. mapFile << fZmin << " " << fZmax << " " << fNz << endl;
  139. // Write field values
  140. Double_t factor = funit * fScale; // Takes out scaling
  141. cout << right;
  142. Int_t nTot = fNx * fNy * fNz;
  143. cout << "-I- MpdFieldMap: " << fNx*fNy*fNz << " entries to write... "
  144. << setw(3) << 0 << " % ";
  145. Int_t index=0;
  146. div_t modul;
  147. Int_t iDiv = TMath::Nint(nTot/100.) + 1;
  148. for(Int_t ix=0; ix<fNx; ix++) {
  149. for(Int_t iy=0; iy<fNy; iy++) {
  150. for(Int_t iz=0; iz<fNz; iz++) {
  151. index =ix*fNy*fNz + iy*fNz + iz;
  152. modul = div(index,iDiv);
  153. if ( modul.rem == 0 ) {
  154. Double_t perc = TMath::Nint(100.*index/nTot);
  155. cout << "\b\b\b\b\b\b" << setw(3) << perc << " % " << flush;
  156. }
  157. mapFile << fBx->At(index)/factor << " " << fBy->At(index)/factor
  158. << " " << fBz->At(index)/factor << endl;
  159. } // z-Loop
  160. } // y-Loop
  161. } // x-Loop
  162. cout << " " << index+1 << " written" << endl;
  163. mapFile.close();
  164. }
  165. // ------- Write field map to a ROOT file -----------------------------
  166. void MpdFieldMap::WriteRootFile(const char* fileName,
  167. const char* mapName) {
  168. MpdFieldMapData* data = new MpdFieldMapData(mapName, *this);
  169. TFile* oldFile = gFile;
  170. TFile* file = new TFile(fileName, "RECREATE");
  171. data->Write();
  172. file->Close();
  173. if(oldFile) oldFile->cd();
  174. }
  175. // ----- Set the position of the field centre in global coordinates -----
  176. void MpdFieldMap::SetPosition(Double_t x, Double_t y, Double_t z) {
  177. fPosX = x;
  178. fPosY = y;
  179. fPosZ = z;
  180. }
  181. // --------- Screen output --------------------------------------------
  182. void MpdFieldMap::Print() {
  183. TString type = "Map";
  184. if (fType == 1 ) type = "Nosym Map";
  185. else if ( fType == 2 ) type = "Solenoid Map ";
  186. else if ( fType == 3 ) type = "Dipole Map ";
  187. else if ( fType == 4 ) type = "Trans Map ";
  188. cout << "======================================================" << endl;
  189. cout.precision(4);
  190. cout << showpoint;
  191. cout << "---- " << fTitle << " : " << fName << endl;
  192. cout << "----" << endl;
  193. cout << "---- Field type : " << type << endl;
  194. cout << "----" << endl;
  195. cout << "---- Field map grid : " << endl;
  196. cout << "---- x = " << setw(4) << fXmin << " to " << setw(4) << fXmax
  197. << " cm, " << fNx << " grid points, dx = " << fXstep << " cm" << endl;
  198. cout << "---- y = " << setw(4) << fYmin << " to " << setw(4) << fYmax
  199. << " cm, " << fNy << " grid points, dy = " << fYstep << " cm" << endl;
  200. cout << "---- z = " << setw(4) << fZmin << " to " << setw(4) << fZmax
  201. << " cm, " << fNz << " grid points, dz = " << fZstep << " cm" << endl;
  202. cout << endl;
  203. cout << "---- Field centre position: ( " << setw(6) << fPosX << ", "
  204. << setw(6) << fPosY << ", " << setw(6) << fPosZ << ") cm" << endl;
  205. cout << "---- Field scaling factor: " << fScale << endl;
  206. Double_t bx = GetBx(0.,0.,0.);
  207. Double_t by = GetBy(0.,0.,0.);
  208. Double_t bz = GetBz(0.,0.,0.);
  209. cout << "----" << endl;
  210. cout << "---- Field at origin is ( " << setw(6) << bx << ", " << setw(6)
  211. << by << ", " << setw(6) << bz << ") kG" << endl;
  212. cout << "======================================================" << endl;
  213. }
  214. // --------- Reset parameters and data (private) ----------------------
  215. void MpdFieldMap::Reset() {
  216. fPosX = fPosY = fPosZ = 0.;
  217. fXmin = fYmin = fZmin = 0.;
  218. fXmax = fYmax = fZmax = 0.;
  219. fXstep = fYstep = fZstep = 0.;
  220. fNx = fNy = fNz = 0;
  221. fScale = 1.;
  222. funit = 10.0;
  223. if ( fBx ) { delete fBx; fBx = NULL; }
  224. if ( fBy ) { delete fBy; fBy = NULL; }
  225. if ( fBz ) { delete fBz; fBz = NULL; }
  226. }
  227. // ----- Read field map from ASCII file (private) ---------------------
  228. void MpdFieldMap::ReadAsciiFile(const char* fileName) {
  229. Double_t xx=0., yy=0., zz=0.;
  230. Double_t bx=0., by=0., bz=0.;
  231. // Open file
  232. cout << "-I- MpdFieldMap: Reading field map from ASCII file "
  233. << fileName << endl;
  234. ifstream mapFile(fileName);
  235. if ( ! mapFile.is_open() ) {
  236. cerr << "-E- MpdFieldMap:ReadAsciiFile: Could not open file! " << endl;
  237. Fatal("ReadAsciiFile","Could not open file");
  238. }
  239. // Read map type
  240. TString type;
  241. mapFile >> type;
  242. Int_t iType = 0;
  243. if ( type == "nosym" ) iType = 1;
  244. else if ( type == "Solenoid") iType = 2;
  245. else if ( type == "Dipole" ) iType = 3;
  246. else if ( type == "Trans" ) iType = 4;
  247. if ( fType != iType ) {
  248. cout << "-E- MpdFieldMap::ReadAsciiFile: Incompatible map types!"
  249. << endl;
  250. cout << " Field map is of type " << fType
  251. << " but map on file is of type " << iType << endl;
  252. Fatal("ReadAsciiFile","Incompatible map types");
  253. }
  254. // Read Units
  255. TString unit;
  256. mapFile >> unit;
  257. if ( unit == "G" ) funit = 0.001;
  258. else if ( unit == "T" ) funit = 10.0;
  259. else if ( unit == "kG" ) funit=1.0;
  260. else {
  261. cout << "-E- FieldMap::ReadAsciiFile: No units!"
  262. << endl;
  263. Fatal("ReadAsciiFile","No units defined");
  264. }
  265. // Read grid parameters
  266. mapFile >>fXmin >> fXmax >> fNx;
  267. mapFile >>fYmin >> fYmax >> fNy;
  268. mapFile >>fZmin >> fZmax >> fNz;
  269. fXstep = ( fXmax - fXmin ) / Double_t( fNx - 1 );
  270. fYstep = ( fYmax - fYmin ) / Double_t( fNy - 1 );
  271. fZstep = ( fZmax - fZmin ) / Double_t( fNz - 1 );
  272. // Create field arrays
  273. fBx = new TArrayF(fNx * fNy * fNz);
  274. fBy = new TArrayF(fNx * fNy * fNz);
  275. fBz = new TArrayF(fNx * fNy * fNz);
  276. // Read the field values
  277. Double_t factor = fScale * funit; // Factor 1/1000 for G -> kG
  278. cout << right;
  279. Int_t nTot = fNx * fNy * fNz;
  280. cout << "-I- MpdFieldMap: " << nTot << " entries to read... "
  281. << setw(3) << 0 << " % ";
  282. Int_t index = 0;
  283. div_t modul;
  284. Int_t iDiv = TMath::Nint(nTot/100.) + 1;
  285. for (Int_t ix=0; ix<fNx; ix++) {
  286. for (Int_t iy = 0; iy<fNy; iy++) {
  287. for (Int_t iz = 0; iz<fNz; iz++) {
  288. if (! mapFile.good()) cerr << "-E- MpdFieldMap::ReadAsciiFile: "
  289. << "I/O Error at " << ix << " "
  290. << iy << " " << iz << endl;
  291. index = ix*fNy*fNz + iy*fNz + iz;
  292. modul = div(index,iDiv);
  293. if ( modul.rem == 0 ) {
  294. Double_t perc = TMath::Nint(100.*index/nTot);
  295. cout << "\b\b\b\b\b\b" << setw(3) << perc << " % " << flush;
  296. }
  297. //mapFile >> bx >> by >> bz;
  298. mapFile >> xx >> yy >> zz >> bx >> by >> bz ;
  299. //cout << " x= " <<xx <<" y= " << yy<<" z= " << zz<<" bx= " << bx <<" by= " <<by <<" bz= " << bz<< endl;
  300. fBx->AddAt(factor*bx, index);
  301. fBy->AddAt(factor*by, index);
  302. fBz->AddAt(factor*bz, index);
  303. if ( mapFile.eof() ) {
  304. cerr << endl << "-E- MpdFieldMap::ReadAsciiFile: EOF"
  305. << " reached at " << ix << " " << iy << " " << iz << endl;
  306. mapFile.close();
  307. break;
  308. }
  309. } // z-Loop
  310. } // y-Loop0)
  311. } // x-Loop
  312. cout << " " << index+1 << " read" << endl;
  313. mapFile.close();
  314. }
  315. // ------------- Read field map from ROOT file (private) ---------------
  316. void MpdFieldMap::ReadRootFile(const char* fileName,
  317. const char* mapName) {
  318. cout<<"Field Map is read from ASCII file ..."<<endl;
  319. }
  320. // ------------ Set field parameters and data (private) ----------------
  321. void MpdFieldMap::SetField(const MpdFieldMapData* data) {
  322. // Check compatibility
  323. if ( data->GetType() != fType ) {
  324. cout << "-E- MpdFieldMap::SetField: Incompatible map types!"
  325. << endl;
  326. cout << " Field map is of type " << fType
  327. << " but map on file is of type " << data->GetType() << endl;
  328. Fatal("SetField","Incompatible map types");
  329. }
  330. fXmin = data->GetXmin();
  331. fYmin = data->GetYmin();
  332. fZmin = data->GetZmin();
  333. fXmax = data->GetXmax();
  334. fYmax = data->GetYmax();
  335. fZmax = data->GetZmax();
  336. fNx = data->GetNx();
  337. fNy = data->GetNy();
  338. fNz = data->GetNz();
  339. fXstep = ( fXmax - fXmin ) / Double_t( fNx - 1 );
  340. fYstep = ( fYmax - fYmin ) / Double_t( fNy - 1 );
  341. fZstep = ( fZmax - fZmin ) / Double_t( fNz - 1 );
  342. if ( fBx ) delete fBx;
  343. if ( fBy ) delete fBy;
  344. if ( fBz ) delete fBz;
  345. fBx = new TArrayF(*(data->GetBx()));
  346. fBy = new TArrayF(*(data->GetBy()));
  347. fBz = new TArrayF(*(data->GetBz()));
  348. // Scale and convert from G(or T) to kG
  349. Double_t factor = fScale * funit;
  350. Int_t index = 0;
  351. for (Int_t ix=0; ix<fNx; ix++) {
  352. for (Int_t iy=0; iy<fNy; iy++) {
  353. for (Int_t iz=0; iz<fNz; iz++) {
  354. index = ix*fNy*fNz + iy*fNz + iz;
  355. if ( fBx ) (*fBx)[index] = (*fBx)[index] * factor;
  356. if ( fBy ) (*fBy)[index] = (*fBy)[index] * factor;
  357. if ( fBz ) (*fBz)[index] = (*fBz)[index] * factor;
  358. }
  359. }
  360. }
  361. }
  362. // ------------------------------------------------------------------------
  363. Double_t MpdFieldMap::MPDCalcPhi(Double_t x, Double_t y) {
  364. /// We assume that Phi is located in the interval [0; 2pi), P.B.
  365. Double_t phi(0.);
  366. if (x > 0. && y>=0.) phi = TMath::ATan(y/x);
  367. else if (x > 0. && y < 0.) phi = TMath::ATan(y/x) + 2*TMath::Pi();
  368. else if (x < 0.) phi = TMath::ATan(y/x) + TMath::Pi();
  369. else if (x == 0. && y > 0.) phi = 0.5*TMath::Pi();
  370. else if (x == 0. && y < 0.) phi = 1.5*TMath::Pi();
  371. else if (x == 0. && y == 0.) phi = 0.;
  372. return phi;
  373. }
  374. //-----------------------------------------------------------------------------------
  375. Double_t MpdFieldMap::MPDCalcBr(Double_t x, Double_t y) {
  376. /// Mag. Field in T!!!
  377. Double_t br_coeff[25] = {3.55683796461738e-005, -1.45190262790289e-005, -4.22589216944936e-005, -3.7825589770834e-006,
  378. 5.31016888137445e-006, -0.000357993220605977, 0.000104590154606646, 0.000421949428873677,
  379. 1.14589001127912e-005, 7.04773281719065e-007, 0.000289792998812542, -0.000269810303185643,
  380. -0.000367922801279625, 8.60754473338966e-005, 0.000114899755175272, -0.000301080291666739,
  381. 0.000447886344136466, 0.000250622861695885, -0.000333504572726095, -0.000111012107349364,
  382. 7.82187500000153e-005, -0.00015617776565864, -8.1476324428263e-005, 0.000159232809881339,
  383. 2.93802936497351e-005};
  384. Double_t sum(0.);
  385. for (int i(0); i<5; i++) sum += br_coeff[i]*pow(y,i);
  386. for (int i(5); i<10; i++) sum += br_coeff[i]*x*pow(y,i-5);
  387. for (int i(10); i<15; i++) sum += br_coeff[i]*pow(x,2)*pow(y,i-10);
  388. for (int i(15); i<20; i++) sum += br_coeff[i]*pow(x,3)*pow(y,i-15);
  389. for (int i(20); i<25; i++) sum += br_coeff[i]*pow(x,4)*pow(y,i-20);
  390. return sum;
  391. }
  392. //-----------------------------------------------------------------------------------
  393. Double_t MpdFieldMap::MPDCalcBz(Double_t x, Double_t y) {
  394. /// Mag. Field in T!!!
  395. Double_t bz_coeff[25] = {0.500047838143456, 0.00044540880088384, -5.35189259780466e-005, -0.000247723183378179,
  396. 2.63504568969453e-005, 1.23804543665921e-005, 0.000177653223608631, 1.66011921364229e-005,
  397. 6.08797488443891e-005, -8.64585638589155e-005, -3.187574503429e-005, -0.000208086128594243,
  398. -9.1879573536513e-005, -1.4640341932548e-005, 0.000284560265815692, 7.60046874443087e-005,
  399. 0.000394401600524841, -0.000120094441661944, 0.000101870491589828, -0.000203451132764876,
  400. -1.8619791659889e-005, -0.000151059550349997, 8.41518325664925e-005, -3.90510111913933e-005,
  401. 2.77271406987634e-005};
  402. Double_t sum(0.);
  403. for (int i(0); i<5; i++) sum += bz_coeff[i]*pow(y,i);
  404. for (int i(5); i<10; i++) sum += bz_coeff[i]*x*pow(y,i-5);
  405. for (int i(10); i<15; i++) sum += bz_coeff[i]*pow(x,2)*pow(y,i-10);
  406. for (int i(15); i<20; i++) sum += bz_coeff[i]*pow(x,3)*pow(y,i-15);
  407. for (int i(20); i<25; i++) sum += bz_coeff[i]*pow(x,4)*pow(y,i-20);
  408. return sum;
  409. }
  410. /// MPDGetBx(y, z) return Mag. Field in kG!
  411. //-------------------------------------------------------------------------------------
  412. Double_t MpdFieldMap::GetBx(Double_t x, Double_t y, Double_t z) {
  413. if (IsInside(x, y, z) == kTRUE) {
  414. x /= 100; y /= 100; z /=100;
  415. Double_t phi = MPDCalcPhi(x, y);
  416. Double_t r = TMath::Sqrt(x*x + y*y);
  417. Double_t Br = MPDCalcBr(r,z);
  418. return Br*TMath::Cos(phi)*fScale*funit;
  419. }
  420. else return 0;
  421. }
  422. //--------------------------------------------------------------------------------------
  423. Double_t MpdFieldMap::GetBy(Double_t x, Double_t y, Double_t z) {
  424. if (IsInside(x, y, z) == kTRUE) {
  425. x /= 100; y /= 100; z /=100;
  426. Double_t phi = MPDCalcPhi(x, y);
  427. Double_t r = TMath::Sqrt(x*x + y*y);
  428. Double_t Br = MPDCalcBr(r,z);
  429. return Br*TMath::Sin(phi)*fScale*funit;
  430. }
  431. else return 0;
  432. }
  433. //--------------------------------------------------------------------------------------
  434. Double_t MpdFieldMap::GetBz(Double_t x, Double_t y, Double_t z) {
  435. if (IsInside(x, y, z) == kTRUE) {
  436. x /= 100; y /= 100; z /= 100;
  437. Double_t phi = MPDCalcPhi(x, y);
  438. Double_t r = TMath::Sqrt(x*x + y*y);
  439. return MPDCalcBz(r,z)*fScale*funit;
  440. }
  441. else return 5.;
  442. }
  443. //---------------------------------------------------------------------------------------
  444. void MpdFieldMap::GetBxyz(const Double_t point[3], Double_t* bField) {
  445. bField[0] = GetBx(point[0], point[1], point[2]);
  446. bField[1] = GetBy(point[0], point[1], point[2]);
  447. bField[2] = GetBz(point[0], point[1], point[2]);
  448. //cout<<"x = " <<point[0] <<" y = " <<point[1] <<" z = " <<point[2];
  449. //cout<<" Bx = "<<bField[0]<<" By = "<<bField[1]<<" Bz = "<<bField[2]<<endl;
  450. }
  451. //---------------------------------------------------------------------------------------
  452. void MpdFieldMap::GetFieldValue(const Double_t point[3], Double_t* bField) {
  453. bField[0] = GetBx(point[0], point[1], point[2]);
  454. bField[1] = GetBy(point[0], point[1], point[2]);
  455. bField[2] = GetBz(point[0], point[1], point[2]);
  456. }
  457. Double_t MpdFieldMap::Interpolate(Double_t dx, Double_t dy, Double_t dz)
  458. {
  459. cout<<"MpdFieldMap::Interpolate( is not implemented!!!"<<endl;
  460. return 0;
  461. }
  462. ClassImp(MpdFieldMap)