CbmRoot
CbmFieldMap.cxx
Go to the documentation of this file.
1 // -------------------------------------------------------------------------
2 // ----- CbmFieldMap source file -----
3 // ----- Created 12/01/04 by M. Al/Turany (CbmField.cxx) -----
4 // ----- Redesign 13/02/06 by V. Friese -----
5 // -------------------------------------------------------------------------
6 #include "CbmFieldMap.h"
7 
8 #include "CbmFieldMapCreator.h" // for CbmFieldMapCreator
9 #include "CbmFieldMapData.h" // for CbmFieldMapData
10 #include "CbmFieldPar.h" // for CbmFieldPar
11 
12 #include <FairField.h> // for FairField
13 #include <FairLogger.h> // for LOG, Logger
14 
15 #include <TArrayF.h> // for TArrayF
16 #include <TFile.h> // for TFile, gFile
17 #include <TMath.h> // for Nint
18 
19 #include <assert.h> // for assert
20 #include <iomanip> // for operator<<, setw
21 #include <iostream> // for operator<<, basic_ostream, endl, ost...
22 #include <stdlib.h> // for div, div_t, getenv, exit
23 
24 using std::cerr;
25 using std::cout;
26 using std::endl;
27 using std::flush;
28 using std::right;
29 using std::setw;
30 using std::showpoint;
31 
32 
33 // ------------- Default constructor ----------------------------------
34 CbmFieldMap::CbmFieldMap()
35  : FairField()
36  , fFileName("")
37  , fScale(1.)
38  , fPosX(0.)
39  , fPosY(0.)
40  , fPosZ(0.)
41  , fXmin(0.)
42  , fXmax(0.)
43  , fXstep(0.)
44  , fYmin(0.)
45  , fYmax(0.)
46  , fYstep(0.)
47  , fZmin(0.)
48  , fZmax(0.)
49  , fZstep(0.)
50  , fNx(0)
51  , fNy(0)
52  , fNz(0)
53  , fBx(nullptr)
54  , fBy(nullptr)
55  , fBz(nullptr)
56  , fBxOrigin(0.)
57  , fByOrigin(0.)
58  , fBzOrigin(0.) {
59  // Initilization of arrays is to my knowledge not
60  // possible in member initalization lists
61  for (Int_t i = 0; i < 2; i++) {
62  fHc[i] = 0;
63  for (Int_t j = 0; j < 2; j++) {
64  fHb[i][j] = 0;
65  for (Int_t k = 0; k < 2; k++) {
66  fHa[i][j][k] = 0;
67  }
68  }
69  }
70  // Assign values to data members of base classes
71  // There is no appropriate constructor of the base
72  // class.
73  fName = "";
74  fType = 1;
75 }
76 // ------------------------------------------------------------------------
77 
78 
79 // ------------- Standard constructor ---------------------------------
80 CbmFieldMap::CbmFieldMap(const char* mapName, const char* fileType)
81  : FairField()
82  , fFileName("")
83  , fScale(1.)
84  , fPosX(0.)
85  , fPosY(0.)
86  , fPosZ(0.)
87  , fXmin(0.)
88  , fXmax(0.)
89  , fXstep(0.)
90  , fYmin(0.)
91  , fYmax(0.)
92  , fYstep(0.)
93  , fZmin(0.)
94  , fZmax(0.)
95  , fZstep(0.)
96  , fNx(0)
97  , fNy(0)
98  , fNz(0)
99  , fBx(nullptr)
100  , fBy(nullptr)
101  , fBz(nullptr)
102  , fBxOrigin(0.)
103  , fByOrigin(0.)
104  , fBzOrigin(0.) {
105  // Initilization of arrays is to my knowledge not
106  // possible in member initalization lists
107  for (Int_t i = 0; i < 2; i++) {
108  fHc[i] = 0;
109  for (Int_t j = 0; j < 2; j++) {
110  fHb[i][j] = 0;
111  for (Int_t k = 0; k < 2; k++) {
112  fHa[i][j][k] = 0;
113  }
114  }
115  }
116  // Assign values to data members of base classes
117  // There is no appropriate constructor of the base
118  // class.
119  fName = mapName;
120  fTitle = "CbmFieldMap";
121  TString dir = getenv("VMCWORKDIR");
122  fFileName = dir + "/input/" + mapName;
123  if (fileType[0] == 'R') {
124  fFileName += ".root";
125  } else {
126  fFileName += ".dat";
127  }
128  fType = 1;
129  LOG(info) << "Filename is " << fFileName;
130 }
131 // ------------------------------------------------------------------------
132 
133 
134 // ------------ Constructor from CbmFieldPar --------------------------
135 CbmFieldMap::CbmFieldMap(CbmFieldPar* fieldPar)
136  : FairField()
137  , fFileName("")
138  , fScale(1.)
139  , fPosX(0.)
140  , fPosY(0.)
141  , fPosZ(0.)
142  , fXmin(0.)
143  , fXmax(0.)
144  , fXstep(0.)
145  , fYmin(0.)
146  , fYmax(0.)
147  , fYstep(0.)
148  , fZmin(0.)
149  , fZmax(0.)
150  , fZstep(0.)
151  , fNx(0)
152  , fNy(0)
153  , fNz(0)
154  , fBx(nullptr)
155  , fBy(nullptr)
156  , fBz(nullptr)
157  , fBxOrigin(0.)
158  , fByOrigin(0.)
159  , fBzOrigin(0.) {
160  // Initilization of arrays is to my knowledge not
161  // possible in member initalization lists
162  for (Int_t i = 0; i < 2; i++) {
163  fHc[i] = 0;
164  for (Int_t j = 0; j < 2; j++) {
165  fHb[i][j] = 0;
166  for (Int_t k = 0; k < 2; k++) {
167  fHa[i][j][k] = 0;
168  }
169  }
170  }
171  // Assign values to data members of base classes
172  // There is no appropriate constructor of the base
173  // class.
174  fName = "";
175  fType = 1;
176  if (!fieldPar) {
177  cerr << "-W- CbmFieldConst::CbmFieldMap: empty parameter container!"
178  << endl;
179  } else {
180  fieldPar->MapName(fName);
181  fPosX = fieldPar->GetPositionX();
182  fPosY = fieldPar->GetPositionY();
183  fPosZ = fieldPar->GetPositionZ();
184  fScale = fieldPar->GetScale();
185  TString dir = getenv("VMCWORKDIR");
186  fFileName = dir + "/input/" + fName + ".root";
187  fType = fieldPar->GetType();
188  }
189 }
190 // ------------------------------------------------------------------------
191 
192 
193 // ------------ Constructor from CbmFieldMapCreator ---------------------
194 CbmFieldMap::CbmFieldMap(CbmFieldMapCreator* creator)
195  : FairField()
196  , fFileName("")
197  , fScale(1.)
198  , fPosX(0.)
199  , fPosY(0.)
200  , fPosZ(0.)
201  , fXmin(0.)
202  , fXmax(0.)
203  , fXstep(0.)
204  , fYmin(0.)
205  , fYmax(0.)
206  , fYstep(0.)
207  , fZmin(0.)
208  , fZmax(0.)
209  , fZstep(0.)
210  , fNx(0)
211  , fNy(0)
212  , fNz(0)
213  , fBx(nullptr)
214  , fBy(nullptr)
215  , fBz(nullptr)
216  , fBxOrigin(0.)
217  , fByOrigin(0.)
218  , fBzOrigin(0.) {
219  // Initilization of arrays is to my knowledge not
220  // possible in member initalization lists
221  for (Int_t i = 0; i < 2; i++) {
222  fHc[i] = 0;
223  for (Int_t j = 0; j < 2; j++) {
224  fHb[i][j] = 0;
225  for (Int_t k = 0; k < 2; k++) {
226  fHa[i][j][k] = 0;
227  }
228  }
229  }
230  // Assign values to data members of base classes
231  // There is no appropriate constructor of the base
232  // class.
233  fName = "";
234  fType = 1;
235  if (!creator) {
236  cerr << "-W- CbmFieldMap: no creator given!" << endl;
237  } else {
238  fType = 1;
239  fName = creator->GetMapName();
240  fXmin = creator->GetXmin();
241  fXmax = creator->GetXmax();
242  fYmin = creator->GetYmin();
243  fYmax = creator->GetYmax();
244  fZmin = creator->GetZmin();
245  fZmax = creator->GetZmax();
246  fNx = creator->GetNx();
247  fNy = creator->GetNy();
248  fNz = creator->GetNz();
249  fXstep = (fXmax - fXmin) / Double_t(fNx - 1);
250  fYstep = (fYmax - fYmin) / Double_t(fNy - 1);
251  fZstep = (fZmax - fZmin) / Double_t(fNz - 1);
252  fBx = creator->GetBx();
253  fBy = creator->GetBy();
254  fBz = creator->GetBz();
255  }
256 }
257 
258 // ------------------------------------------------------------------------
259 
260 
261 // ------------ Destructor --------------------------------------------
262 CbmFieldMap::~CbmFieldMap() {
263  if (fBx) delete fBx;
264  if (fBy) delete fBy;
265  if (fBz) delete fBz;
266 }
267 // ------------------------------------------------------------------------
268 
269 
270 // ----------- Intialisation ------------------------------------------
271 void CbmFieldMap::Init() {
272  if (fFileName.EndsWith(".root"))
273  ReadRootFile(fFileName, fName);
274  else if (fFileName.EndsWith(".dat"))
275  ReadAsciiFile(fFileName);
276  else {
277  cerr << "-E- CbmFieldMap::Init: No proper file name defined! (" << fFileName
278  << ")" << endl;
279  Fatal("Init", "No proper file name");
280  }
281  // Fill values needed in the Print() function. This is needed to allow
282  // a constant Print() function.
283  fBxOrigin = GetBx(0., 0., 0.);
284  fByOrigin = GetBy(0., 0., 0.);
285  fBzOrigin = GetBz(0., 0., 0.);
286 
287  Print();
288 }
289 // ------------------------------------------------------------------------
290 
291 
292 // ----- Initialisation from arrays -------------------------------------
293 void CbmFieldMap::Init(Int_t nX,
294  Double_t xMin,
295  Double_t xMax,
296  Int_t nY,
297  Double_t yMin,
298  Double_t yMax,
299  Int_t nZ,
300  Double_t zMin,
301  Double_t zMax,
302  TArrayF* bx,
303  TArrayF* by,
304  TArrayF* bz) {
305 
306  Reset();
307  fNx = nX;
308  fXmin = xMin;
309  fXmax = xMax;
310  fNy = nY;
311  fYmin = yMin;
312  fYmax = yMax;
313  fNz = nZ;
314  fZmin = zMin;
315  fZmax = zMax;
316  Int_t nPoints = nX * nY * nZ;
317  assert(bx->GetSize() == by->GetSize());
318  assert(bz->GetSize() == bx->GetSize());
319  if (bx->GetSize() != nPoints) {
320  LOG(error) << "CbmFieldMap: array size " << bx->GetSize()
321  << " does not match number of grid points.";
322  return;
323  }
324  fBx = new TArrayF(nPoints);
325  fBy = new TArrayF(nPoints);
326  fBz = new TArrayF(nPoints);
327  Int_t index = 0;
328  for (Int_t ix = 0; ix < fNx; ix++) {
329  for (Int_t iy = 0; iy < fNy; iy++) {
330  for (Int_t iz = 0; iz < fNz; iz++) {
331  index = ix * fNy * fNz + iy * fNz + iz;
332  (*fBx)[index] = (*bx)[index];
333  (*fBy)[index] = (*by)[index];
334  (*fBz)[index] = (*bz)[index];
335  } //# iz
336  } //# iy
337  } //# ix
338 
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  fBxOrigin = GetBx(0., 0., 0.);
343  fByOrigin = GetBy(0., 0., 0.);
344  fBzOrigin = GetBz(0., 0., 0.);
345 
346  LOG(info) << "CbmFieldMap: Initialised from " << nPoints << " grid points";
347  Print();
348 }
349 // ------------------------------------------------------------------------
350 
351 
352 // ----------- Get x component of the field ---------------------------
353 Double_t CbmFieldMap::GetBx(Double_t x, Double_t y, Double_t z) {
354 
355  Int_t ix = 0;
356  Int_t iy = 0;
357  Int_t iz = 0;
358  Double_t dx = 0.;
359  Double_t dy = 0.;
360  Double_t dz = 0.;
361 
362  if (IsInside(x, y, z, ix, iy, iz, dx, dy, dz)) {
363 
364  // Get Bx field values at grid cell corners
365  fHa[0][0][0] = fBx->At(ix * fNy * fNz + iy * fNz + iz);
366  fHa[1][0][0] = fBx->At((ix + 1) * fNy * fNz + iy * fNz + iz);
367  fHa[0][1][0] = fBx->At(ix * fNy * fNz + (iy + 1) * fNz + iz);
368  fHa[1][1][0] = fBx->At((ix + 1) * fNy * fNz + (iy + 1) * fNz + iz);
369  fHa[0][0][1] = fBx->At(ix * fNy * fNz + iy * fNz + (iz + 1));
370  fHa[1][0][1] = fBx->At((ix + 1) * fNy * fNz + iy * fNz + (iz + 1));
371  fHa[0][1][1] = fBx->At(ix * fNy * fNz + (iy + 1) * fNz + (iz + 1));
372  fHa[1][1][1] = fBx->At((ix + 1) * fNy * fNz + (iy + 1) * fNz + (iz + 1));
373 
374  // Return interpolated field value
375  return Interpolate(dx, dy, dz);
376  }
377 
378  return 0.;
379 }
380 // ------------------------------------------------------------------------
381 
382 
383 // ----------- Get y component of the field ---------------------------
384 Double_t CbmFieldMap::GetBy(Double_t x, Double_t y, Double_t z) {
385 
386  Int_t ix = 0;
387  Int_t iy = 0;
388  Int_t iz = 0;
389  Double_t dx = 0.;
390  Double_t dy = 0.;
391  Double_t dz = 0.;
392 
393  if (IsInside(x, y, z, ix, iy, iz, dx, dy, dz)) {
394 
395  // Get By field values at grid cell corners
396  fHa[0][0][0] = fBy->At(ix * fNy * fNz + iy * fNz + iz);
397  fHa[1][0][0] = fBy->At((ix + 1) * fNy * fNz + iy * fNz + iz);
398  fHa[0][1][0] = fBy->At(ix * fNy * fNz + (iy + 1) * fNz + iz);
399  fHa[1][1][0] = fBy->At((ix + 1) * fNy * fNz + (iy + 1) * fNz + iz);
400  fHa[0][0][1] = fBy->At(ix * fNy * fNz + iy * fNz + (iz + 1));
401  fHa[1][0][1] = fBy->At((ix + 1) * fNy * fNz + iy * fNz + (iz + 1));
402  fHa[0][1][1] = fBy->At(ix * fNy * fNz + (iy + 1) * fNz + (iz + 1));
403  fHa[1][1][1] = fBy->At((ix + 1) * fNy * fNz + (iy + 1) * fNz + (iz + 1));
404 
405  // Return interpolated field value
406  return Interpolate(dx, dy, dz);
407  }
408 
409  return 0.;
410 }
411 // ------------------------------------------------------------------------
412 
413 
414 // ----------- Get z component of the field ---------------------------
415 Double_t CbmFieldMap::GetBz(Double_t x, Double_t y, Double_t z) {
416 
417  Int_t ix = 0;
418  Int_t iy = 0;
419  Int_t iz = 0;
420  Double_t dx = 0.;
421  Double_t dy = 0.;
422  Double_t dz = 0.;
423 
424  if (IsInside(x, y, z, ix, iy, iz, dx, dy, dz)) {
425 
426  // Get Bz field values at grid cell corners
427  fHa[0][0][0] = fBz->At(ix * fNy * fNz + iy * fNz + iz);
428  fHa[1][0][0] = fBz->At((ix + 1) * fNy * fNz + iy * fNz + iz);
429  fHa[0][1][0] = fBz->At(ix * fNy * fNz + (iy + 1) * fNz + iz);
430  fHa[1][1][0] = fBz->At((ix + 1) * fNy * fNz + (iy + 1) * fNz + iz);
431  fHa[0][0][1] = fBz->At(ix * fNy * fNz + iy * fNz + (iz + 1));
432  fHa[1][0][1] = fBz->At((ix + 1) * fNy * fNz + iy * fNz + (iz + 1));
433  fHa[0][1][1] = fBz->At(ix * fNy * fNz + (iy + 1) * fNz + (iz + 1));
434  fHa[1][1][1] = fBz->At((ix + 1) * fNy * fNz + (iy + 1) * fNz + (iz + 1));
435 
436  // Return interpolated field value
437  return Interpolate(dx, dy, dz);
438  }
439 
440  return 0.;
441 }
442 // ------------------------------------------------------------------------
443 
444 
445 // ----------- Check whether a point is inside the map ----------------
446 Bool_t CbmFieldMap::IsInside(Double_t x,
447  Double_t y,
448  Double_t z,
449  Int_t& ix,
450  Int_t& iy,
451  Int_t& iz,
452  Double_t& dx,
453  Double_t& dy,
454  Double_t& dz) {
455 
456  // --- Transform into local coordinate system
457  Double_t xl = x - fPosX;
458  Double_t yl = y - fPosY;
459  Double_t zl = z - fPosZ;
460 
461  // --- Check for being outside the map range
462  if (!(xl >= fXmin && xl < fXmax && yl >= fYmin && yl < fYmax && zl >= fZmin
463  && zl < fZmax)) {
464  ix = iy = iz = 0;
465  dx = dy = dz = 0.;
466  return kFALSE;
467  }
468 
469  // --- Determine grid cell
470  ix = Int_t((xl - fXmin) / fXstep);
471  iy = Int_t((yl - fYmin) / fYstep);
472  iz = Int_t((zl - fZmin) / fZstep);
473 
474 
475  // Relative distance from grid point (in units of cell size)
476  dx = (xl - fXmin) / fXstep - Double_t(ix);
477  dy = (yl - fYmin) / fYstep - Double_t(iy);
478  dz = (zl - fZmin) / fZstep - Double_t(iz);
479 
480  return kTRUE;
481 }
482 // ------------------------------------------------------------------------
483 
484 
485 // ---------- Write the map to an ASCII file --------------------------
486 void CbmFieldMap::WriteAsciiFile(const char* fileName) {
487 
488  // Open file
489  LOG(info) << "Writing field map to ASCII file " << fileName;
490  std::ofstream mapFile(fileName);
491  if (!mapFile.is_open()) {
492  LOG(error) << "Could not open file! ";
493  return;
494  }
495 
496  // Write field map grid parameters
497  mapFile.precision(6);
498  mapFile << showpoint;
499  if (fType == 1) mapFile << "nosym" << endl;
500  if (fType == 2) mapFile << "sym2" << endl;
501  if (fType == 3) mapFile << "sym3" << endl;
502  mapFile << fXmin << " " << fXmax << " " << fNx << endl;
503  mapFile << fYmin << " " << fYmax << " " << fNy << endl;
504  mapFile << fZmin << " " << fZmax << " " << fNz << endl;
505 
506  // Write field values
507  Double_t factor = 10. * fScale; // Takes out scaling and converts kG->T
508  cout << right;
509  Int_t nTot = fNx * fNy * fNz;
510  cout << "-I- CbmFieldMap: " << fNx * fNy * fNz << " entries to write... "
511  << setw(3) << 0 << " % ";
512  Int_t index = 0;
513  div_t modul;
514  Int_t iDiv = TMath::Nint(nTot / 100.);
515  for (Int_t ix = 0; ix < fNx; ix++) {
516  for (Int_t iy = 0; iy < fNy; iy++) {
517  for (Int_t iz = 0; iz < fNz; iz++) {
518  index = ix * fNy * fNz + iy * fNz + iz;
519  modul = div(index, iDiv);
520  if (modul.rem == 0) {
521  Double_t perc = TMath::Nint(100. * index / nTot);
522  cout << "\b\b\b\b\b\b" << setw(3) << perc << " % " << flush;
523  }
524  mapFile << fBx->At(index) / factor << " " << fBy->At(index) / factor
525  << " " << fBz->At(index) / factor << endl;
526  } // z-Loop
527  } // y-Loop
528  } // x-Loop
529  cout << " " << index + 1 << " written" << endl;
530  mapFile.close();
531 }
532 // ------------------------------------------------------------------------
533 
534 
535 // ------- Write field map to a ROOT file -----------------------------
536 void CbmFieldMap::WriteRootFile(const char* fileName, const char* mapName) {
537 
538  CbmFieldMapData* data = new CbmFieldMapData(mapName, *this);
539  TFile* oldFile = gFile;
540  TFile* file = new TFile(fileName, "RECREATE");
541  data->Write();
542  file->Close();
543  if (oldFile) oldFile->cd();
544 }
545 // ------------------------------------------------------------------------
546 
547 
548 // ----- Set the position of the field centre in global coordinates -----
549 void CbmFieldMap::SetPosition(Double_t x, Double_t y, Double_t z) {
550  fPosX = x;
551  fPosY = y;
552  fPosZ = z;
553 }
554 // ------------------------------------------------------------------------
555 
556 
557 // --------- Screen output --------------------------------------------
558 void CbmFieldMap::Print(Option_t*) const {
559  TString type = "Map";
560  if (fType == 2) type = "Map sym2";
561  if (fType == 3) type = "Map sym3";
562  cout << "======================================================" << endl;
563  cout.precision(4);
564  cout << showpoint;
565  cout << "---- " << fTitle << " : " << fName << endl;
566  cout << "----" << endl;
567  cout << "---- Field type : " << type << endl;
568  cout << "----" << endl;
569  cout << "---- Field map grid : " << endl;
570  cout << "---- x = " << setw(4) << fXmin << " to " << setw(4) << fXmax
571  << " cm, " << fNx << " grid points, dx = " << fXstep << " cm" << endl;
572  cout << "---- y = " << setw(4) << fYmin << " to " << setw(4) << fYmax
573  << " cm, " << fNy << " grid points, dy = " << fYstep << " cm" << endl;
574  cout << "---- z = " << setw(4) << fZmin << " to " << setw(4) << fZmax
575  << " cm, " << fNz << " grid points, dz = " << fZstep << " cm" << endl;
576  cout << endl;
577  cout << "---- Field centre position: ( " << setw(6) << fPosX << ", "
578  << setw(6) << fPosY << ", " << setw(6) << fPosZ << ") cm" << endl;
579  cout << "---- Field scaling factor: " << fScale << endl;
580  cout << "----" << endl;
581  cout << "---- Field at origin is ( " << setw(6) << fBxOrigin << ", "
582  << setw(6) << fByOrigin << ", " << setw(6) << fBzOrigin << ") kG"
583  << endl;
584 
585  cout << "======================================================" << endl;
586 }
587 // ------------------------------------------------------------------------
588 
589 
590 // --------- Reset parameters and data (private) ----------------------
591 void CbmFieldMap::Reset() {
592  fPosX = fPosY = fPosZ = 0.;
593  fXmin = fYmin = fZmin = 0.;
594  fXmax = fYmax = fZmax = 0.;
595  fXstep = fYstep = fZstep = 0.;
596  fNx = fNy = fNz = 0;
597  fScale = 1.;
598  if (fBx) {
599  delete fBx;
600  fBx = nullptr;
601  }
602  if (fBy) {
603  delete fBy;
604  fBy = nullptr;
605  }
606  if (fBz) {
607  delete fBz;
608  fBz = nullptr;
609  }
610 }
611 // ------------------------------------------------------------------------
612 
613 
614 // ----- Read field map from ASCII file (private) ---------------------
615 void CbmFieldMap::ReadAsciiFile(const char* fileName) {
616 
617  Double_t bx = 0., by = 0., bz = 0.;
618 
619  // Open file
620  LOG(info) << "CbmFieldMap: Reading field map from ASCII file " << fileName;
621  std::ifstream mapFile(fileName);
622  if (!mapFile.is_open()) {
623  LOG(error) << "CbmFieldMap:ReadAsciiFile: Could not open file!";
624  LOG(fatal) << "CbmFieldMap:ReadAsciiFile: Could not open file!";
625  }
626 
627  // Read map type
628  TString type;
629  mapFile >> type;
630  Int_t iType = 0;
631  if (type == "nosym") iType = 1;
632  if (type == "sym2") iType = 2;
633  if (type == "sym3") iType = 3;
634  if (fType != iType) {
635  cout << "-E- CbmFieldMap::ReadAsciiFile: Incompatible map types!" << endl;
636  cout << " Field map is of type " << fType
637  << " but map on file is of type " << iType << endl;
638  Fatal("ReadAsciiFile", "Incompatible map types");
639  }
640 
641  // Read grid parameters
642  mapFile >> fXmin >> fXmax >> fNx;
643  mapFile >> fYmin >> fYmax >> fNy;
644  mapFile >> fZmin >> fZmax >> fNz;
645  fXstep = (fXmax - fXmin) / Double_t(fNx - 1);
646  fYstep = (fYmax - fYmin) / Double_t(fNy - 1);
647  fZstep = (fZmax - fZmin) / Double_t(fNz - 1);
648 
649  // Create field arrays
650  fBx = new TArrayF(fNx * fNy * fNz);
651  fBy = new TArrayF(fNx * fNy * fNz);
652  fBz = new TArrayF(fNx * fNy * fNz);
653 
654  // Read the field values
655  Double_t factor = fScale * 10.; // Factor 10 for T -> kG
656  cout << right;
657  Int_t nTot = fNx * fNy * fNz;
658  cout << "-I- CbmFieldMap: " << nTot << " entries to read... " << setw(3) << 0
659  << " % ";
660  Int_t index = 0;
661  div_t modul;
662  Int_t iDiv = TMath::Nint(nTot / 100.);
663  for (Int_t ix = 0; ix < fNx; ix++) {
664  for (Int_t iy = 0; iy < fNy; iy++) {
665  for (Int_t iz = 0; iz < fNz; iz++) {
666  if (!mapFile.good())
667  cerr << "-E- CbmFieldMap::ReadAsciiFile: "
668  << "I/O Error at " << ix << " " << iy << " " << iz << endl;
669  index = ix * fNy * fNz + iy * fNz + iz;
670  modul = div(index, iDiv);
671  if (modul.rem == 0) {
672  Double_t perc = TMath::Nint(100. * index / nTot);
673  cout << "\b\b\b\b\b\b" << setw(3) << perc << " % " << flush;
674  }
675  mapFile >> bx >> by >> bz;
676  fBx->AddAt(factor * bx, index);
677  fBy->AddAt(factor * by, index);
678  fBz->AddAt(factor * bz, index);
679  if (mapFile.eof()) {
680  cerr << endl
681  << "-E- CbmFieldMap::ReadAsciiFile: EOF"
682  << " reached at " << ix << " " << iy << " " << iz << endl;
683  mapFile.close();
684  break;
685  }
686  } // z-Loop
687  } // y-Loop0)
688  } // x-Loop
689 
690  cout << " " << index + 1 << " read" << endl;
691 
692  mapFile.close();
693 }
694 // ------------------------------------------------------------------------
695 
696 
697 // ------------- Read field map from ROOT file (private) ---------------
698 void CbmFieldMap::ReadRootFile(const char* fileName, const char* mapName) {
699 
700  // Store gFile pointer
701  TFile* oldFile = gFile;
702 
703  // Open root file
704  LOG(info) << "CbmFieldMap: Reading field map from ROOT file " << fileName;
705  TFile* file = new TFile(fileName, "READ");
706  if (!(file->IsOpen())) {
707  LOG(error) << "CbmFieldMap:ReadRootFile: Could not open file!";
708  LOG(fatal) << "CbmFieldMap:ReadRootFile: Could not open file!";
709  }
710 
711  // Get the field data object
712  CbmFieldMapData* data = nullptr;
713  file->GetObject(mapName, data);
714  if (!data) {
715  cout << "-E- CbmFieldMap::ReadRootFile: data object " << fileName
716  << " not found in file! " << endl;
717  exit(-1);
718  }
719 
720  // Get the field parameters
721  SetField(data);
722 
723  // Close the root file and delete the data object
724  file->Close();
725  delete data;
726  if (oldFile) oldFile->cd();
727 }
728 // ------------------------------------------------------------------------
729 
730 
731 // ------------ Set field parameters and data (private) ----------------
732 void CbmFieldMap::SetField(const CbmFieldMapData* data) {
733 
734  // Check compatibility
735  if (data->GetType() != fType) {
736  if (!((data->GetType() == 3) && (fType == 5))) // E.Litvinenko
737  {
738  cout << "-E- CbmFieldMap::SetField: Incompatible map types!" << endl;
739  cout << " Field map is of type " << fType
740  << " but map on file is of type " << data->GetType() << endl;
741  Fatal("SetField", "Incompatible map types");
742  } else
743  cout << " CbmFieldMap::SetField: Warning: You are using PosDepScaled "
744  "map (original map type = 3)"
745  << endl;
746  }
747 
748 
749  fXmin = data->GetXmin();
750  fYmin = data->GetYmin();
751  fZmin = data->GetZmin();
752  fXmax = data->GetXmax();
753  fYmax = data->GetYmax();
754  fZmax = data->GetZmax();
755  fNx = data->GetNx();
756  fNy = data->GetNy();
757  fNz = data->GetNz();
758  fXstep = (fXmax - fXmin) / Double_t(fNx - 1);
759  fYstep = (fYmax - fYmin) / Double_t(fNy - 1);
760  fZstep = (fZmax - fZmin) / Double_t(fNz - 1);
761  if (fBx) delete fBx;
762  if (fBy) delete fBy;
763  if (fBz) delete fBz;
764  fBx = new TArrayF(*(data->GetBx()));
765  fBy = new TArrayF(*(data->GetBy()));
766  fBz = new TArrayF(*(data->GetBz()));
767 
768  // Scale and convert from T to kG
769  Double_t factor = fScale * 10.;
770  Int_t index = 0;
771  for (Int_t ix = 0; ix < fNx; ix++) {
772  for (Int_t iy = 0; iy < fNy; iy++) {
773  for (Int_t iz = 0; iz < fNz; iz++) {
774  index = ix * fNy * fNz + iy * fNz + iz;
775  if (fBx) (*fBx)[index] = (*fBx)[index] * factor;
776  if (fBy) (*fBy)[index] = (*fBy)[index] * factor;
777  if (fBz) (*fBz)[index] = (*fBz)[index] * factor;
778  }
779  }
780  }
781 }
782 // ------------------------------------------------------------------------
783 
784 
785 // ------------ Interpolation in a grid cell (private) -----------------
786 Double_t CbmFieldMap::Interpolate(Double_t dx, Double_t dy, Double_t dz) {
787 
788  // Interpolate in x coordinate
789  fHb[0][0] = fHa[0][0][0] + (fHa[1][0][0] - fHa[0][0][0]) * dx;
790  fHb[1][0] = fHa[0][1][0] + (fHa[1][1][0] - fHa[0][1][0]) * dx;
791  fHb[0][1] = fHa[0][0][1] + (fHa[1][0][1] - fHa[0][0][1]) * dx;
792  fHb[1][1] = fHa[0][1][1] + (fHa[1][1][1] - fHa[0][1][1]) * dx;
793 
794  // Interpolate in y coordinate
795  fHc[0] = fHb[0][0] + (fHb[1][0] - fHb[0][0]) * dy;
796  fHc[1] = fHb[0][1] + (fHb[1][1] - fHb[0][1]) * dy;
797 
798  // Interpolate in z coordinate
799  return fHc[0] + (fHc[1] - fHc[0]) * dz;
800 }
801 // ------------------------------------------------------------------------
802 
803 
804 ClassImp(CbmFieldMap)
CbmFieldMap::fFileName
TString fFileName
Definition: CbmFieldMap.h:197
CbmFieldMapData::GetBz
TArrayF * GetBz() const
Definition: CbmFieldMapData.h:64
CbmFieldMapData::GetYmin
Double_t GetYmin() const
Definition: CbmFieldMapData.h:51
CbmFieldMap::Interpolate
Double_t Interpolate(Double_t dx, Double_t dy, Double_t dz)
Definition: CbmFieldMap.cxx:786
CbmFieldMap::ReadRootFile
void ReadRootFile(const char *fileName, const char *mapName)
Definition: CbmFieldMap.cxx:698
CbmFieldPar::GetPositionY
Double_t GetPositionY() const
Definition: CbmFieldPar.h:69
CbmFieldMapData::GetYmax
Double_t GetYmax() const
Definition: CbmFieldMapData.h:54
CbmFieldMap::fYmin
Double_t fYmin
Definition: CbmFieldMap.h:210
CbmFieldMap::fScale
Double_t fScale
Definition: CbmFieldMap.h:201
CbmFieldMapData::GetXmax
Double_t GetXmax() const
Definition: CbmFieldMapData.h:53
CbmFieldMap::fBxOrigin
Double_t fBxOrigin
Interpolated field (1-dim)
Definition: CbmFieldMap.h:230
CbmFieldMapCreator::GetXmax
Double_t GetXmax() const
Definition: CbmFieldMapCreator.h:83
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
CbmFieldMapCreator::GetYmax
Double_t GetYmax() const
Definition: CbmFieldMapCreator.h:85
CbmFieldMap::GetBy
TArrayF * GetBy() const
Definition: CbmFieldMap.h:159
CbmFieldMapCreator
Definition: CbmFieldMapCreator.h:30
CbmFieldMap::fXmin
Double_t fXmin
Definition: CbmFieldMap.h:209
CbmFieldMapCreator::GetZmin
Double_t GetZmin() const
Definition: CbmFieldMapCreator.h:86
CbmFieldMap::fPosZ
Double_t fPosZ
Definition: CbmFieldMap.h:205
CbmFieldMapCreator::GetXmin
Double_t GetXmin() const
Definition: CbmFieldMapCreator.h:82
CbmFieldMap::Print
virtual void Print(Option_t *="") const
Definition: CbmFieldMap.cxx:558
CbmFieldMap::fNz
Int_t fNz
Definition: CbmFieldMap.h:215
CbmFieldMapData.h
CbmFieldMap::fNy
Int_t fNy
Definition: CbmFieldMap.h:215
CbmFieldMapData::GetZmin
Double_t GetZmin() const
Definition: CbmFieldMapData.h:52
CbmFieldPar::GetPositionX
Double_t GetPositionX() const
Definition: CbmFieldPar.h:68
CbmFieldMapData::GetNx
Int_t GetNx() const
Definition: CbmFieldMapData.h:56
CbmFieldMap::fXstep
Double_t fXstep
Definition: CbmFieldMap.h:209
CbmFieldMapData::GetBy
TArrayF * GetBy() const
Definition: CbmFieldMapData.h:63
CbmFieldMap::SetPosition
virtual void SetPosition(Double_t x, Double_t y, Double_t z)
Definition: CbmFieldMap.cxx:549
CbmFieldMap::fPosX
Double_t fPosX
Definition: CbmFieldMap.h:205
CbmFieldMap::fBzOrigin
Double_t fBzOrigin
y-component of the field at the origin
Definition: CbmFieldMap.h:232
CbmFieldMapCreator::GetMapName
TString GetMapName() const
Definition: CbmFieldMapCreator.h:78
CbmFieldMapCreator::GetBx
TArrayF * GetBx() const
Definition: CbmFieldMapCreator.h:88
CbmFieldMap::fZstep
Double_t fZstep
Definition: CbmFieldMap.h:211
CbmFieldPar::GetScale
Double_t GetScale() const
Definition: CbmFieldPar.h:71
CbmFieldMap::SetField
void SetField(const CbmFieldMapData *data)
Definition: CbmFieldMap.cxx:732
CbmFieldMap::fZmin
Double_t fZmin
Definition: CbmFieldMap.h:211
CbmFieldMap::fBz
TArrayF * fBz
Definition: CbmFieldMap.h:221
CbmFieldMapCreator::GetNz
Int_t GetNz() const
Definition: CbmFieldMapCreator.h:81
CbmFieldMap::fBy
TArrayF * fBy
Definition: CbmFieldMap.h:220
CbmFieldMap::GetBz
TArrayF * GetBz() const
Definition: CbmFieldMap.h:160
CbmFieldMap
Definition: CbmFieldMap.h:34
CbmFieldMap::fHc
Double_t fHc[2]
Interpolated field (2-dim)
Definition: CbmFieldMap.h:228
CbmFieldMapCreator::GetBy
TArrayF * GetBy() const
Definition: CbmFieldMapCreator.h:89
CbmFieldMapData::GetType
Int_t GetType() const
Definition: CbmFieldMapData.h:49
CbmFieldMap::fPosY
Double_t fPosY
Definition: CbmFieldMap.h:205
CbmFieldMapCreator::GetNy
Int_t GetNy() const
Definition: CbmFieldMapCreator.h:80
CbmFieldPar.h
CbmFieldMapCreator.h
CbmFieldMapCreator::GetBz
TArrayF * GetBz() const
Definition: CbmFieldMapCreator.h:90
ClassImp
ClassImp(CbmConverterManager) InitStatus CbmConverterManager
Definition: CbmConverterManager.cxx:12
CbmFieldMap::fBx
TArrayF * fBx
Definition: CbmFieldMap.h:219
CbmFieldMapData::GetXmin
Double_t GetXmin() const
Definition: CbmFieldMapData.h:50
CbmFieldMapCreator::GetNx
Int_t GetNx() const
Definition: CbmFieldMapCreator.h:79
CbmFieldPar::GetPositionZ
Double_t GetPositionZ() const
Definition: CbmFieldPar.h:70
CbmFieldMap::IsInside
virtual Bool_t IsInside(Double_t x, Double_t y, Double_t z, Int_t &ix, Int_t &iy, Int_t &iz, Double_t &dx, Double_t &dy, Double_t &dz)
Definition: CbmFieldMap.cxx:446
CbmFieldMap::WriteRootFile
void WriteRootFile(const char *fileName, const char *mapName)
Definition: CbmFieldMap.cxx:536
CbmFieldMapCreator::GetZmax
Double_t GetZmax() const
Definition: CbmFieldMapCreator.h:87
CbmFieldMap::ReadAsciiFile
void ReadAsciiFile(const char *fileName)
Definition: CbmFieldMap.cxx:615
CbmFieldPar::MapName
void MapName(TString &name)
Definition: CbmFieldPar.h:67
CbmFieldMap::fYstep
Double_t fYstep
Definition: CbmFieldMap.h:210
CbmFieldMap::fYmax
Double_t fYmax
Definition: CbmFieldMap.h:210
CbmFieldMapData::GetNy
Int_t GetNy() const
Definition: CbmFieldMapData.h:57
CbmFieldMap::fZmax
Double_t fZmax
Definition: CbmFieldMap.h:211
CbmFieldPar::GetType
Int_t GetType() const
Definition: CbmFieldPar.h:57
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmFieldMapData
Definition: CbmFieldMapData.h:29
CbmFieldMap::WriteAsciiFile
void WriteAsciiFile(const char *fileName)
Definition: CbmFieldMap.cxx:486
y
Double_t y
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmFieldMap::fHb
Double_t fHb[2][2]
Field at corners of a grid cell.
Definition: CbmFieldMap.h:227
CbmFieldMapData::GetZmax
Double_t GetZmax() const
Definition: CbmFieldMapData.h:55
CbmFieldMapCreator::GetYmin
Double_t GetYmin() const
Definition: CbmFieldMapCreator.h:84
CbmFieldMap::Init
virtual void Init()
Definition: CbmFieldMap.cxx:271
CbmFieldMap::fHa
Double_t fHa[2][2][2]
Definition: CbmFieldMap.h:226
CbmFieldMap::fNx
Int_t fNx
Definition: CbmFieldMap.h:215
CbmFieldMap.h
CbmFieldMap::CbmFieldMap
CbmFieldMap()
Definition: CbmFieldMap.cxx:34
CbmFieldPar
Definition: CbmFieldPar.h:31
CbmFieldMap::fByOrigin
Double_t fByOrigin
x-component of the field at the origin
Definition: CbmFieldMap.h:231
CbmFieldMapData::GetNz
Int_t GetNz() const
Definition: CbmFieldMapData.h:58
CbmFieldMap::Reset
void Reset()
Definition: CbmFieldMap.cxx:591
CbmFieldMapData::GetBx
TArrayF * GetBx() const
Definition: CbmFieldMapData.h:62
CbmFieldMap::fXmax
Double_t fXmax
Definition: CbmFieldMap.h:209
CbmFieldMap::~CbmFieldMap
virtual ~CbmFieldMap()
Definition: CbmFieldMap.cxx:262
CbmFieldMap::GetBx
TArrayF * GetBx() const
Definition: CbmFieldMap.h:158