EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PndFieldMap.cxx
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PndFieldMap.cxx
1 // -------------------------------------------------------------------------
2 // ----- PndFieldMap source file -----
3 // ----- Created 12/01/04 by M. Al/Turany (FairField.cxx) -----
4 // -------------------------------------------------------------------------
5 
6 
7 #include <iomanip>
8 #include <iostream>
9 #include <fstream>
10 #include "stdlib.h"
11 
12 #include "TArrayF.h"
13 #include "TFile.h"
14 #include "TMath.h"
15 
16 #include "PndFieldMap.h"
17 #include "PndFieldMapData.h"
18 #include "PndFieldPar.h"
19 
20 using namespace std;
21 
22 // ------------- Default constructor ----------------------------------
23 PndFieldMap::PndFieldMap()
24  : FairField(),
25  fFileName(""),
26  fScale(1.0),
27  funit(10.0),
28  fPosX(0), fPosY(0), fPosZ(0),
29  fXmin(0), fXmax(0), fXstep(0),
30  fYmin(0), fYmax(0), fYstep(0),
31  fZmin(0), fZmax(0), fZstep(0),
32  fNx(0),fNy(0),fNz(0),
33  fBx(NULL), fBy(NULL), fBz(NULL)
34 {
35  SetName("");
36  fType = 1;
37 }
38 // ------------------------------------------------------------------------
39 
40 
41 
42 // ------------- Standard constructor ---------------------------------
43 PndFieldMap::PndFieldMap(const char* mapName, const char* fileType)
44  : FairField(mapName),
45  fFileName(TString("")),
46  fScale(1.0),
47  funit(10.0),
48  fPosX(0), fPosY(0), fPosZ(0),
49  fXmin(0), fXmax(0), fXstep(0),
50  fYmin(0), fYmax(0), fYstep(0),
51  fZmin(0), fZmax(0), fZstep(0),
52  fNx(0),fNy(0),fNz(0),
53  fBx(NULL), fBy(NULL), fBz(NULL)
54 {
55  SetName(mapName);
56  TString dir = getenv("VMCWORKDIR");
57  fFileName = dir + "/input/" + mapName;
58  if ( fileType[0] == 'R' ) fFileName += ".root";
59  else fFileName += ".dat";
60  fType = 1;
61 }
62 // ------------------------------------------------------------------------
63 
64 
65 
66 // ------------ Constructor from PndFieldPar --------------------------
67 PndFieldMap::PndFieldMap(PndFieldPar* fieldPar)
68  : FairField(),
69  fFileName(TString("")),
70  fScale(1.0),
71  funit(10.0),
72  fPosX(0), fPosY(0), fPosZ(0),
73  fXmin(0), fXmax(0), fXstep(0),
74  fYmin(0), fYmax(0), fYstep(0),
75  fZmin(0), fZmax(0), fZstep(0),
76  fNx(0),fNy(0),fNz(0),
77  fBx(NULL), fBy(NULL), fBz(NULL)
78 {
79  fType = 1;
80  if ( ! fieldPar ) {
81  gLogger->Error(MESSAGE_ORIGIN, "PndConstField::PndConstField: empty parameter container!");
82  SetName("");
83  fType = -1;
84  }
85  else {
86  TString Name=GetName();
87  fieldPar->MapName(Name);
88  fPosX = fieldPar->GetPositionX();
89  fPosY = fieldPar->GetPositionY();
90  fPosZ = fieldPar->GetPositionZ();
91  fScale = fieldPar->GetScale();
92  TString dir = getenv("VMCWORKDIR");
93  fFileName = dir + "/input/" + Name + ".root";
94  fType = fieldPar->GetType();
95  }
96 }
97 // ------------------------------------------------------------------------
98 
99 
100 
101 // ------------ Destructor --------------------------------------------
102 PndFieldMap::~PndFieldMap() {
103 
104  //printf("PndFieldMap::~PndFieldMap() \n");
105  if ( fBx ) delete fBx;
106  if ( fBy ) delete fBy;
107  if ( fBz ) delete fBz;
108 }
109 // ------------------------------------------------------------------------
110 
111 
112 
113 // ----------- Intialisation ------------------------------------------
114 void PndFieldMap::Init() {
115  if (fFileName.EndsWith(".root")) ReadRootFile(fFileName, fName);
116  else if (fFileName.EndsWith(".dat")) ReadAsciiFile(fFileName);
117  else {
118  gLogger->Error(MESSAGE_ORIGIN,"-E- PndFieldMap::Init: No proper file name defined! (%s) ");
119  Fatal("Init", "No proper file name");
120  }
121 }
122 // ------------------------------------------------------------------------
123 
124 
125 
126 // ----------- Get x component of the field ---------------------------
127 Double_t PndFieldMap::GetBx(Double_t x, Double_t y, Double_t z) {
128 
129  Int_t ix = 0;
130  Int_t iy = 0;
131  Int_t iz = 0;
132  Double_t dx = 0.;
133  Double_t dy = 0.;
134  Double_t dz = 0.;
135 
136  if ( IsInside(x, y, z, ix, iy, iz, dx, dy, dz) ) {
137 
138  // Get Bx field values at grid cell corners
139  fHa[0][0][0] = fBx->At(ix *fNy*fNz + iy *fNz + iz);
140  fHa[1][0][0] = fBx->At((ix+1)*fNy*fNz + iy *fNz + iz);
141  fHa[0][1][0] = fBx->At(ix *fNy*fNz + (iy+1)*fNz + iz);
142  fHa[1][1][0] = fBx->At((ix+1)*fNy*fNz + (iy+1)*fNz + iz);
143  fHa[0][0][1] = fBx->At(ix *fNy*fNz + iy *fNz + (iz+1));
144  fHa[1][0][1] = fBx->At((ix+1)*fNy*fNz + iy *fNz + (iz+1));
145  fHa[0][1][1] = fBx->At(ix *fNy*fNz + (iy+1)*fNz + (iz+1));
146  fHa[1][1][1] = fBx->At((ix+1)*fNy*fNz + (iy+1)*fNz + (iz+1));
147 
148  // Return interpolated field value
149  return Interpolate(dx, dy, dz);
150 
151  }
152 
153  return 0.;
154 }
155 // ------------------------------------------------------------------------
156 
157 
158 
159 // ----------- Get y component of the field ---------------------------
160 Double_t PndFieldMap::GetBy(Double_t x, Double_t y, Double_t z) {
161 
162  Int_t ix = 0;
163  Int_t iy = 0;
164  Int_t iz = 0;
165  Double_t dx = 0.;
166  Double_t dy = 0.;
167  Double_t dz = 0.;
168 
169  if ( IsInside(x, y, z, ix, iy, iz, dx, dy, dz) ) {
170 
171  // Get By field values at grid cell corners
172  fHa[0][0][0] = fBy->At(ix *fNy*fNz + iy *fNz + iz);
173  fHa[1][0][0] = fBy->At((ix+1)*fNy*fNz + iy *fNz + iz);
174  fHa[0][1][0] = fBy->At(ix *fNy*fNz + (iy+1)*fNz + iz);
175  fHa[1][1][0] = fBy->At((ix+1)*fNy*fNz + (iy+1)*fNz + iz);
176  fHa[0][0][1] = fBy->At(ix *fNy*fNz + iy *fNz + (iz+1));
177  fHa[1][0][1] = fBy->At((ix+1)*fNy*fNz + iy *fNz + (iz+1));
178  fHa[0][1][1] = fBy->At(ix *fNy*fNz + (iy+1)*fNz + (iz+1));
179  fHa[1][1][1] = fBy->At((ix+1)*fNy*fNz + (iy+1)*fNz + (iz+1));
180 
181  // Return interpolated field value
182  return Interpolate(dx, dy, dz);
183 
184  }
185 
186  return 0.;
187 }
188 // ------------------------------------------------------------------------
189 
190 
191 
192 // ----------- Get z component of the field ---------------------------
193 Double_t PndFieldMap::GetBz(Double_t x, Double_t y, Double_t z) {
194 
195  Int_t ix = 0;
196  Int_t iy = 0;
197  Int_t iz = 0;
198  Double_t dx = 0.;
199  Double_t dy = 0.;
200  Double_t dz = 0.;
201 
202  if ( IsInside(x, y, z, ix, iy, iz, dx, dy, dz) ) {
203 
204  // Get Bz field values at grid cell corners
205  fHa[0][0][0] = fBz->At(ix *fNy*fNz + iy *fNz + iz);
206  fHa[1][0][0] = fBz->At((ix+1)*fNy*fNz + iy *fNz + iz);
207  fHa[0][1][0] = fBz->At(ix *fNy*fNz + (iy+1)*fNz + iz);
208  fHa[1][1][0] = fBz->At((ix+1)*fNy*fNz + (iy+1)*fNz + iz);
209  fHa[0][0][1] = fBz->At(ix *fNy*fNz + iy *fNz + (iz+1));
210  fHa[1][0][1] = fBz->At((ix+1)*fNy*fNz + iy *fNz + (iz+1));
211  fHa[0][1][1] = fBz->At(ix *fNy*fNz + (iy+1)*fNz + (iz+1));
212  fHa[1][1][1] = fBz->At((ix+1)*fNy*fNz + (iy+1)*fNz + (iz+1));
213 
214  // Return interpolated field value
215  return Interpolate(dx, dy, dz);
216 
217  }
218 
219  return 0.;
220 }
221 // ------------------------------------------------------------------------
222 
223 
224 
225 // ----------- Check whether a point is inside the map ----------------
226 Bool_t PndFieldMap::IsInside(Double_t x, Double_t y, Double_t z,
227  Int_t& ix, Int_t& iy, Int_t& iz,
228  Double_t& dx, Double_t& dy, Double_t& dz) {
229 
230  // --- Transform into local coordinate system
231  Double_t xl = x - fPosX;
232  Double_t yl = y - fPosY;
233  Double_t zl = z - fPosZ;
234 
235  // --- Check for being outside the map range
236  if ( ! ( xl >= fXmin && xl < fXmax && yl >= fYmin && yl < fYmax &&
237  zl >= fZmin && zl < fZmax ) ) {
238  ix = iy = iz = 0;
239  dx = dy = dz = 0.;
240  return kFALSE;
241  }
242 
243  // --- Determine grid cell
244  ix = Int_t( xl / fXstep );
245  iy = Int_t( yl / fYstep );
246  iz = Int_t( zl / fZstep );
247 
248 
249  // Relative distance from grid point (in units of cell size)
250  dx = xl / fXstep - Double_t(ix);
251  dy = yl / fYstep - Double_t(iy);
252  dz = zl / fZstep - Double_t(iz);
253 
254  return kTRUE;
255 
256 }
257 // ------------------------------------------------------------------------
258 
259 
260 
261 // ---------- Write the map to an ASCII file --------------------------
262 void PndFieldMap::WriteAsciiFile(const char* fileName) {
263 
264  // Open file
265  gLogger->Info(MESSAGE_ORIGIN, "PndFieldMap: Writing field map to ASCII file %s ",fileName);
266  ofstream mapFile(fileName);
267  if ( ! mapFile.is_open() ) {
268  gLogger->Error(MESSAGE_ORIGIN, "PndFieldMap:ReadAsciiFile: Could not open file! ");
269  return;
270  }
271 
272  // Write field map grid parameters
273  mapFile.precision(4);
274  mapFile << showpoint;
275  if ( fType == 1 ) mapFile << "nosym" << endl;
276  if ( fType == 2 ) mapFile << "Solenoid" << endl;
277  if ( fType == 3 ) mapFile << "Dipole" << endl;
278  if ( fType == 4 ) mapFile << "Trans" << endl;
279  if ( funit == 10.0 ) mapFile << "T" << endl;
280  if ( funit == 0.001 ) mapFile << "G" << endl;
281  if ( funit == 1.0 ) mapFile << "kG" << endl;
282 
283  mapFile << fXmin << " " << fXmax << " " << fNx << endl;
284  mapFile << fYmin << " " << fYmax << " " << fNy << endl;
285  mapFile << fZmin << " " << fZmax << " " << fNz << endl;
286 
287  // Write field values
288  Double_t factor = funit * fScale; // Takes out scaling
289  cout << right;
290  Int_t nTot = fNx * fNy * fNz;
291  cout << "-I- PndFieldMap: " << fNx*fNy*fNz << " entries to write... "
292  << setw(3) << 0 << " % ";
293  Int_t index=0;
294  div_t modul;
295  Int_t iDiv = TMath::Nint(nTot/100.);
296  for(Int_t ix=0; ix<fNx; ix++) {
297  for(Int_t iy=0; iy<fNy; iy++) {
298  for(Int_t iz=0; iz<fNz; iz++) {
299  index =ix*fNy*fNz + iy*fNz + iz;
300  if (iDiv!=0)
301  {
302  modul = div(index,iDiv);
303  if ( modul.rem == 0 ) {
304  Double_t perc = TMath::Nint(100.*index/nTot);
305  cout << "\b\b\b\b\b\b" << setw(3) << perc << " % " << flush;
306  }
307  }
308  mapFile << fBx->At(index)/factor << " " << fBy->At(index)/factor
309  << " " << fBz->At(index)/factor << endl;
310  } // z-Loop
311  } // y-Loop
312  } // x-Loop
313  cout << " " << index+1 << " written" << endl;
314  mapFile.close();
315 
316 }
317 // ------------------------------------------------------------------------
318 
319 
320 
321 // ------- Write field map to a ROOT file -----------------------------
322 void PndFieldMap::WriteRootFile(const char* fileName,
323  const char* mapName) {
324 
325  PndFieldMapData* data = new PndFieldMapData(mapName, *this);
326  TFile* oldFile = gFile;
327  TFile* file = new TFile(fileName, "RECREATE");
328  data->Write();
329  file->Close();
330  if(oldFile) oldFile->cd();
331 
332 }
333 // ------------------------------------------------------------------------
334 
335 
336 
337 // ----- Set the position of the field centre in global coordinates -----
338 void PndFieldMap::SetPosition(Double_t x, Double_t y, Double_t z) {
339  fPosX = x;
340  fPosY = y;
341  fPosZ = z;
342 }
343 // ------------------------------------------------------------------------
344 
345 
346 
347 // --------- Screen output --------------------------------------------
348 void PndFieldMap::Print() {
349  TString type = "Map";
350  if ( fType == 2 ) type = "Soleniod Map ";
351  if ( fType == 3 ) type = "Dipole Map ";
352  if ( fType == 4 ) type = "Trans Map ";
353  cout << "======================================================" << endl;
354  cout.precision(4);
355  cout << showpoint;
356  cout << "---- " << fTitle << " : " << fName << endl;
357  cout << "----" << endl;
358  cout << "---- Field type : " << type << endl;
359  cout << "----" << endl;
360  cout << "---- Field map grid : " << endl;
361  cout << "---- x = " << setw(4) << fXmin << " to " << setw(4) << fXmax
362  << " cm, " << fNx << " grid points, dx = " << fXstep << " cm" << endl;
363  cout << "---- y = " << setw(4) << fYmin << " to " << setw(4) << fYmax
364  << " cm, " << fNy << " grid points, dy = " << fYstep << " cm" << endl;
365  cout << "---- z = " << setw(4) << fZmin << " to " << setw(4) << fZmax
366  << " cm, " << fNz << " grid points, dz = " << fZstep << " cm" << endl;
367  cout << endl;
368  cout << "---- Field centre position: ( " << setw(6) << fPosX << ", "
369  << setw(6) << fPosY << ", " << setw(6) << fPosZ << ") cm" << endl;
370  cout << "---- Field scaling factor: " << fScale << endl;
371  Double_t bx = GetBx(0.,0.,0.);
372  Double_t by = GetBy(0.,0.,0.);
373  Double_t bz = GetBz(0.,0.,0.);
374  cout << "----" << endl;
375  cout << "---- Field at origin is ( " << setw(6) << bx << ", " << setw(6)
376  << by << ", " << setw(6) << bz << ") kG" << endl;
377  cout << "======================================================" << endl;
378 }
379 // ------------------------------------------------------------------------
380 
381 
382 
383 // --------- Reset parameters and data (private) ----------------------
384 void PndFieldMap::Reset() {
385  fPosX = fPosY = fPosZ = 0.;
386  fXmin = fYmin = fZmin = 0.;
387  fXmax = fYmax = fZmax = 0.;
388  fXstep = fYstep = fZstep = 0.;
389  fNx = fNy = fNz = 0;
390  fScale = 1.;
391  funit = 10.0;
392  if ( fBx ) { delete fBx; fBx = NULL; }
393  if ( fBy ) { delete fBy; fBy = NULL; }
394  if ( fBz ) { delete fBz; fBz = NULL; }
395 }
396 // ------------------------------------------------------------------------
397 
398 
399 
400 // ----- Read field map from ASCII file (private) ---------------------
401 void PndFieldMap::ReadAsciiFile(const char* fileName) {
402 
403  Double_t bx=0., by=0., bz=0.;
404  Double_t xx, yy, zz;
405  // Open file
406  cout << "-I- PndFieldMap: Reading field map from ASCII file "
407  << fileName << endl;
408  ifstream mapFile(fileName);
409  if ( ! mapFile.is_open() ) {
410  cerr << "-E- PndFieldMap:ReadAsciiFile: Could not open file! " << endl;
411  Fatal("ReadAsciiFile","Could not open file");
412  }
413 
414  // Read map type
415  TString type;
416  mapFile >> type;
417 
418  Int_t iType = 0;
419  if ( type == "nosym" ) iType = 1;
420  if ( type == "Solenoid") iType = 2;
421  if ( type == "Dipole" ) iType = 3;
422  if ( type == "Trans" ) iType = 4;
423  if ( fType != iType ) {
424  cout << "-E- PndFieldMap::ReadAsciiFile: Incompatible map types!"
425  << endl;
426  cout << " Field map is of type " << fType
427  << " but map on file is of type " << iType << endl;
428  Fatal("ReadAsciiFile","Incompatible map types");
429  }
430  // Read Units
431  TString unit;
432  mapFile >> unit;
433  if ( unit == "G" ) funit = 0.001;
434  else if ( unit == "T" ) funit = 10.0;
435  else if ( unit == "kG" ) funit=1.0;
436  else {
437  cout << "-E- FieldMap::ReadAsciiFile: No units!"
438  << endl;
439  Fatal("ReadAsciiFile","No units defined");
440  }
441 
442 
443  // Read grid parameters
444 
445  mapFile >>fXmin >> fXmax >> fNx;
446  mapFile >>fYmin >> fYmax >> fNy;
447  mapFile >>fZmin >> fZmax >> fNz;
448  fXstep = ( fXmax - fXmin ) / Double_t( fNx - 1 );
449  fYstep = ( fYmax - fYmin ) / Double_t( fNy - 1 );
450  fZstep = ( fZmax - fZmin ) / Double_t( fNz - 1 );
451 
452  // Create field arrays
453  fBx = new TArrayF(fNx * fNy * fNz);
454  fBy = new TArrayF(fNx * fNy * fNz);
455  fBz = new TArrayF(fNx * fNy * fNz);
456 
457  // Read the field values
458  Double_t factor = fScale * funit; // Factor 1/1000 for G -> kG
459  cout << right;
460  Int_t nTot = fNx * fNy * fNz;
461  cout << "-I- PndFieldMap: " << nTot << " entries to read... "
462  << setw(3) << 0 << " % ";
463  Int_t index = 0;
464  div_t modul;
465  Int_t iDiv = TMath::Nint(nTot/100.);
466  for (Int_t ix=0; ix<fNx; ix++) {
467  for (Int_t iy = 0; iy<fNy; iy++) {
468  for (Int_t iz = 0; iz<fNz; iz++) {
469  if (! mapFile.good()) cerr << "-E- PndFieldMap::ReadAsciiFile: "
470  << "I/O Error at " << ix << " "
471  << iy << " " << iz << endl;
472  index = ix*fNy*fNz + iy*fNz + iz;
473  if (iDiv!=0)
474  {
475  modul = div(index,iDiv);
476  if ( modul.rem == 0 ) {
477  Double_t perc = TMath::Nint(100.*index/nTot);
478  cout << "\b\b\b\b\b\b" << setw(3) << perc << " % " << flush;
479  }
480  }
481  //mapFile >> xx>>yy>>zz>> bx >> by >> bz ;
482  mapFile >> bx >> by >> bz ;
483  //cout << " x= " <<xx <<" y= " << yy<<" z= " << zz<<" bx= " << bx <<" by= " <<by <<" bz= " << bz<< endl;
484  fBx->AddAt(factor*bx, index);
485  fBy->AddAt(factor*by, index);
486  fBz->AddAt(factor*bz, index);
487  if ( mapFile.eof() ) {
488  cerr << endl << "-E- PndFieldMap::ReadAsciiFile: EOF"
489  << " reached at " << ix << " " << iy << " " << iz << endl;
490  mapFile.close();
491  break;
492  }
493  } // z-Loop
494  } // y-Loop0)
495  } // x-Loop
496 
497  cout << " " << index+1 << " read" << endl;
498 
499  mapFile.close();
500 
501 }
502 // ------------------------------------------------------------------------
503 
504 
505 
506 // ------------- Read field map from ROOT file (private) ---------------
507 void PndFieldMap::ReadRootFile(const char* fileName,
508  const char* mapName) {
509 
510  // Store gFile pointer
511  TFile* oldFile = gFile;
512 
513  // Open root file
514  gLogger->Info(MESSAGE_ORIGIN, "PndFieldMap: Reading field map from ROOT file %s ",fileName);
515  TFile* file = new TFile(fileName, "READ");
516  if (file->IsZombie()) {
517  gLogger->Error(MESSAGE_ORIGIN, "-E- PndFieldMap::ReadRootfile: Cannot read from file! ");
518  Fatal("ReadRootFile","Cannot read from file");
519  }
520 
521  // Get the field data object
522  PndFieldMapData* data = NULL;
523  file->GetObject(mapName, data);
524  if ( ! data ) {
525  gLogger->Error(MESSAGE_ORIGIN,"PndFieldMap::ReadRootFile: data object %s not found in file! ", fileName);
526  exit(-1);
527  }
528 
529  // Get the field parameters
530  SetField(data);
531 
532  // Close the root file and delete the data object
533  file->Close();
534  delete data;
535  delete file;
536  if ( oldFile ) oldFile->cd();
537 
538 }
539 // ------------------------------------------------------------------------
540 
541 
542 
543 // ------------ Set field parameters and data (private) ----------------
544 void PndFieldMap::SetField(const PndFieldMapData* data) {
545 
546  // Check compatibility
547  if ( data->GetType() != fType ) {
548  gLogger->Error(MESSAGE_ORIGIN,"PndFieldMap::SetField: Incompatible map types Field map is of type %s \n but map on file is of type %s ",fType,data->GetType());
549  Fatal("SetField","Incompatible map types");
550  }
551 
552 
553  fXmin = data->GetXmin();
554  fYmin = data->GetYmin();
555  fZmin = data->GetZmin();
556  fXmax = data->GetXmax();
557  fYmax = data->GetYmax();
558  fZmax = data->GetZmax();
559  fNx = data->GetNx();
560  fNy = data->GetNy();
561  fNz = data->GetNz();
562  fXstep = ( fXmax - fXmin ) / Double_t( fNx - 1 );
563  fYstep = ( fYmax - fYmin ) / Double_t( fNy - 1 );
564  fZstep = ( fZmax - fZmin ) / Double_t( fNz - 1 );
565  if ( fBx ) delete fBx;
566  if ( fBy ) delete fBy;
567  if ( fBz ) delete fBz;
568  fBx = new TArrayF(*(data->GetBx()));
569  fBy = new TArrayF(*(data->GetBy()));
570  fBz = new TArrayF(*(data->GetBz()));
571 
572  // Scale and convert from G(or T) to kG
573  Double_t factor = fScale * funit;
574  Int_t index = 0;
575  for (Int_t ix=0; ix<fNx; ix++) {
576  for (Int_t iy=0; iy<fNy; iy++) {
577  for (Int_t iz=0; iz<fNz; iz++) {
578  index = ix*fNy*fNz + iy*fNz + iz;
579  if ( fBx ) (*fBx)[index] = (*fBx)[index] * factor;
580  if ( fBy ) (*fBy)[index] = (*fBy)[index] * factor;
581  if ( fBz ) (*fBz)[index] = (*fBz)[index] * factor;
582  }
583  }
584  }
585 
586 }
587 // ------------------------------------------------------------------------
588 
589 
590 
591 // ------------ Interpolation in a grid cell (private) -----------------
592 Double_t PndFieldMap::Interpolate(Double_t dx, Double_t dy, Double_t dz) {
593 
594  // Interpolate in x coordinate
595  fHb[0][0] = fHa[0][0][0] + ( fHa[1][0][0]-fHa[0][0][0] ) * dx;
596  fHb[1][0] = fHa[0][1][0] + ( fHa[1][1][0]-fHa[0][1][0] ) * dx;
597  fHb[0][1] = fHa[0][0][1] + ( fHa[1][0][1]-fHa[0][0][1] ) * dx;
598  fHb[1][1] = fHa[0][1][1] + ( fHa[1][1][1]-fHa[0][1][1] ) * dx;
599 
600  // Interpolate in y coordinate
601  fHc[0] = fHb[0][0] + ( fHb[1][0] - fHb[0][0] ) * dy;
602  fHc[1] = fHb[0][1] + ( fHb[1][1] - fHb[0][1] ) * dy;
603 
604  // Interpolate in z coordinate
605  return fHc[0] + ( fHc[1] - fHc[0] ) * dz;
606 
607 }
608 // ------------------------------------------------------------------------
609 
610 
611 
612 ClassImp(PndFieldMap)