GEMC  1.8
Geant4 Monte-Carlo Framework
MagneticField.cc
Go to the documentation of this file.
1 
8 // %%%%%%%%%%
9 // G4 headers
10 // %%%%%%%%%%
11 #include "G4UniformMagField.hh"
12 #include "G4PropagatorInField.hh"
13 #include "G4TransportationManager.hh"
14 #include "G4Mag_UsualEqRhs.hh"
15 #include "G4MagIntegratorStepper.hh"
16 #include "G4ChordFinder.hh"
17 #include "G4ClassicalRK4.hh"
18 #include "G4HelixSimpleRunge.hh"
19 
20 // %%%%%%%%%%
21 // Qt headers
22 // %%%%%%%%%%
23 #include <QtSql>
24 
25 // %%%%%%%%%%%%%
26 // gemc headers
27 // %%%%%%%%%%%%%
28 #include "detector.h"
29 #include "MagneticField.h"
30 #include "usage.h"
31 #include "string_utilities.h"
32 
33 map<string, MagneticField> get_magnetic_Fields(gemc_opts Opt)
34 {
35  string hd_msg = Opt.args["LOG_MSG"].args + " Magnetic Field >> " ;
36  double MGN_VERBOSITY = Opt.args["MGN_VERBOSITY"].arg;
37  string database = Opt.args["DATABASE"].args;
38  string dbtable = "magnetic_fields";
39  string dbhost = Opt.args["DBHOST"].args;
40 
41  string dbUser = Opt.args["DBUSER"].args;
42  string dbPswd = Opt.args["DBPSWD"].args;
43 
44  vector<opts> FIELD_SCALES = Opt.get_args("SCALE_FIELD");
45  map<string, double> SCALES;
46  for(unsigned int f=0; f<FIELD_SCALES.size(); f++)
47  {
48  string SCALEF = FIELD_SCALES[f].args;
49  string fieldname;
50  int ppos = SCALEF.find(",");
51  fieldname.assign(SCALEF, 0, ppos);
52  SCALEF.assign(SCALEF, ppos+1, SCALEF.size());
53  SCALES[fieldname] = atof(SCALEF.c_str());
54  }
55 
56  map<string, MagneticField> FieldMap;
57 
58  QSqlDatabase db = QSqlDatabase::addDatabase("QMYSQL");
59  db.setHostName(dbhost.c_str());
60  db.setDatabaseName(database.c_str());
61  db.setUserName( dbUser.c_str() );
62  db.setPassword( dbPswd.c_str() );
63 
64  bool ok = db.open();
65 
66  if(!ok)
67  {
68  cout << hd_msg << " Database not connected. Exiting." << endl;
69  exit(-1);
70  }
71 
72  MagneticField magneticField;
73  QSqlQuery q;
74  string dbexecute = "select name, type, magnitude, swim_method, description from " + dbtable ;
75  q.exec(dbexecute.c_str());
76 
77  if(MGN_VERBOSITY>2)
78  cout << hd_msg << " Available Magnetic Fields: " << endl << endl;
79 
80  while (q.next())
81  {
82  magneticField.name = TrimSpaces(gemc_tostring(q.value(0).toString()));
83  magneticField.type = gemc_tostring(q.value(1).toString());
84  magneticField.magnitude = gemc_tostring(q.value(2).toString());
85  magneticField.swim_method = gemc_tostring(q.value(3).toString());
86  magneticField.description = gemc_tostring(q.value(4).toString());
87 
88  // Sets MFM pointer to NULL
89  magneticField.init_MFM();
90  magneticField.gemcOpt = Opt;
91 
92  // Adjust Field Scale Factor if any
93  magneticField.scale_factor = 1;
94  if(SCALES.find(magneticField.name) != SCALES.end())
95  {
96  magneticField.scale_factor = SCALES[magneticField.name];
97  cout << hd_msg << " SCALE FACTOR: Magnetic Field " << magneticField.name << " scaled by: " << magneticField.scale_factor << endl;
98 
99  }
100 
101  FieldMap[magneticField.name] = magneticField;
102 
103  if(MGN_VERBOSITY>2)
104  {
105  cout << " ";
106  cout.width(15);
107  cout << magneticField.name << " | ";
108  cout << magneticField.description << endl;
109  cout << " ";
110  cout.width(15);
111  cout << magneticField.name << " | type: | ";
112  cout << magneticField.type << endl;
113  cout << " ";
114  cout.width(15);
115  cout << magneticField.name << " | Magnitude: | ";
116  cout << magneticField.magnitude << endl;
117  cout << " ";
118  cout.width(15);
119  cout << magneticField.name << " | Swim Method: | ";
120  cout << magneticField.swim_method << endl;
121  cout << " -------------------------------------------------------------- " << endl;
122  }
123  }
124  db.close();
125 
126  cout << endl;
127  db = QSqlDatabase();
128  db.removeDatabase("qt_sql_default_connection");
129 
130  return FieldMap;
131 }
132 
133 
134 void MagneticField::create_MFM()
135 {
136  string hd_msg = gemcOpt.args["LOG_MSG"].args + " Magnetic Field: >> ";
137  double MGN_VERBOSITY = gemcOpt.args["MGN_VERBOSITY"].arg;
138  string catch_v = gemcOpt.args["CATCH"].args;
139 
140  stringstream vars;
141  string var, format, symmetry, MapFile;
142  string var1, var2, var3, dim;
143  string dimSpc, dimFld;
144 
145  string integration_method;
146 
147  vars << type;
148  vars >> var >> format >> symmetry >> MapFile;
149 
150  mappedfield = NULL;
151 
152  // %%%%%%%%%%%%%%%%%%%%%%
153  // Uniform Magnetic Field
154  // %%%%%%%%%%%%%%%%%%%%%%
155  if(var == "uniform")
156  {
157  vars.clear();
158  vars << magnitude;
159  double x,y,z;
160  vars >> var; x = scale_factor*get_number(var);
161  vars >> var; y = scale_factor*get_number(var);
162  vars >> var; z = scale_factor*get_number(var);
163  if(MGN_VERBOSITY>3)
164  {
165  cout << hd_msg << " <" << name << "> is a uniform magnetic field type, scaled by " << scale_factor << endl;
166  cout << hd_msg << " <" << name << "> dimensions:" ;
167  cout << " (" << x/gauss << ", " << y/gauss << ", " << z/gauss << ") gauss." << endl;
168  }
169  G4UniformMagField* magField = new G4UniformMagField(G4ThreeVector((x/gauss)*gauss, (y/gauss)*gauss, (z/gauss)*gauss));
170 
171  G4Mag_UsualEqRhs* iEquation = new G4Mag_UsualEqRhs(magField);
172  G4MagIntegratorStepper* iStepper = new G4ClassicalRK4(iEquation);
173  G4ChordFinder* iChordFinder = new G4ChordFinder(magField, 1.0e-2*mm, iStepper);
174 
175  MFM = new G4FieldManager(magField, iChordFinder);
176 
177  return;
178  }
179 
180  // %%%%%%%%%%%%%%%%%%%%%%%
181  // Mapped Field
182  // 1D-dipole field
183  // Dependent on 2 transverse, longitudinal
184  // cartesian coordinates
185  // %%%%%%%%%%%%%%%%%%%%%%%
186  // For now it's "z-z" We'll see if in the future we need a more general
187  // code.
188  if(var == "mapped" && symmetry == "Dipole-y")
189  {
190  vars.clear();
191  vars << magnitude;
192  int TNPOINTS, LNPOINTS;
193  double tlimits[2], llimits[2];
194  double mapOrigin[3];
195  string units[5];
196 
197  vars >> var;
198  LNPOINTS = (int) get_number(var);
199  vars >> var1 >> var2 >> dim;
200  llimits[0] = get_number(var1 + "*" + dim);
201  llimits[1] = get_number(var2 + "*" + dim);
202  units[0].assign(dim);
203 
204  vars >> var;
205  TNPOINTS = (int) get_number(var);
206  vars >> var1 >> var2 >> dim;
207  tlimits[0] = get_number(var1 + "*" + dim);
208  tlimits[1] = get_number(var2 + "*" + dim);
209  units[1].assign(dim);
210 
211 
212  // Origin and Units
213  vars >> var1 >> var2 >> var3 >> dimSpc >> dimFld;
214  mapOrigin[0] = get_number(var1 + "*" + dimSpc);
215  mapOrigin[1] = get_number(var2 + "*" + dimSpc);
216  mapOrigin[2] = get_number(var3 + "*" + dimSpc);
217  units[3].assign( dimSpc );
218  units[4].assign( dimFld);
219 
220  vars.clear();
221  vars << swim_method;
222  vars >> integration_method;
223  vars.clear();
224  vars << swim_method;
225  vars >> integration_method;
226  if(MGN_VERBOSITY>3)
227  {
228  cout << hd_msg << " <" << name << "> is a " << symmetry << " mapped field in cartesian coordinates" << endl;;
229  cout << hd_msg << " <" << name << "> has (" << TNPOINTS << ", " << LNPOINTS << ") points." << endl;;
230  cout << hd_msg << " Tranverse Boundaries: (" << tlimits[0]/cm << ", " << tlimits[1]/cm << ") cm." << endl;
231  cout << hd_msg << " Longitudinal Boundaries: (" << llimits[0]/cm << ", " << llimits[1]/cm << ") cm." << endl;
232  cout << hd_msg << " Map Displacement: (" << mapOrigin[0]/cm << ", "
233  << mapOrigin[1]/cm << ", "
234  << mapOrigin[2]/cm << ") cm." << endl;
235  cout << hd_msg << " Integration Method: " << integration_method << endl;
236  cout << hd_msg << " Map Field Units: " << units[4] << endl;
237  cout << hd_msg << " Map File: " << MapFile << endl;
238  }
239  mappedfield = new MappedField(gemcOpt, TNPOINTS, LNPOINTS, tlimits, llimits, MapFile, mapOrigin, units, scale_factor, symmetry);
240  }
241 
242  // %%%%%%%%%%%%%%%%%%%%%%%
243  // Mapped Field
244  // phi-symmetric
245  // cylindrical coordinates
246  // %%%%%%%%%%%%%%%%%%%%%%%
247 
248  // by default the symmetry is around z. But it could be around x or y axis
249  if(var == "mapped" && (symmetry == "cylindrical" || symmetry == "cylindrical-x" || symmetry == "cylindrical-y"))
250  {
251  vars.clear();
252  vars << magnitude;
253  int TNPOINTS, LNPOINTS;
254  double tlimits[2], llimits[2];
255  double mapOrigin[3];
256  string units[5];
257 
258  vars >> var;
259  TNPOINTS = (int) get_number(var);
260  vars >> var1 >> var2 >> dim;
261  tlimits[0] = get_number(var1 + "*" + dim);
262  tlimits[1] = get_number(var2 + "*" + dim);
263  units[0].assign(dim);
264 
265 
266  vars >> var;
267  LNPOINTS = (int) get_number(var);
268  vars >> var1 >> var2 >> dim;
269  llimits[0] = get_number(var1 + "*" + dim);
270  llimits[1] = get_number(var2 + "*" + dim);
271  units[1].assign(dim);
272 
273  // Origin and Units
274  vars >> var1 >> var2 >> var3 >> dimSpc >> dimFld;
275  mapOrigin[0] = get_number(var1 + "*" + dimSpc);
276  mapOrigin[1] = get_number(var2 + "*" + dimSpc);
277  mapOrigin[2] = get_number(var3 + "*" + dimSpc);
278  units[3].assign( dimSpc );
279  units[4].assign( dimFld);
280 
281  vars.clear();
282  vars << swim_method;
283  vars >> integration_method;
284  vars.clear();
285  vars << swim_method;
286  vars >> integration_method;
287  if(MGN_VERBOSITY>3)
288  {
289  cout << hd_msg << " <" << name << "> is a phi-symmetric mapped field in cylindrical coordinates" << endl;;
290  cout << hd_msg << " <" << name << "> has (" << TNPOINTS << ", " << LNPOINTS << ") points." << endl;;
291  cout << hd_msg << " Tranverse Boundaries: (" << tlimits[0]/cm << ", " << tlimits[1]/cm << ") cm." << endl;
292  cout << hd_msg << " Longitudinal Boundaries: (" << llimits[0]/cm << ", " << llimits[1]/cm << ") cm." << endl;
293  cout << hd_msg << " Map Displacement: (" << mapOrigin[0]/cm << ", "
294  << mapOrigin[1]/cm << ", "
295  << mapOrigin[2]/cm << ") cm." << endl;
296  cout << hd_msg << " Integration Method: " << integration_method << endl;
297  cout << hd_msg << " Map Field Units: " << units[4] << endl;
298  cout << hd_msg << " Map File: " << MapFile << endl;
299  }
300  // last argument is dummy
301  mappedfield = new MappedField(gemcOpt, TNPOINTS, LNPOINTS, tlimits, llimits, MapFile, mapOrigin, units, scale_factor, symmetry, 0);
302  }
303 
304 
305 
306  // %%%%%%%%%%%%%%%%%%%%%%%
307  // Mapped Field
308  // phi-segmented
309  // cylindrical coordinates
310  // %%%%%%%%%%%%%%%%%%%%%%%
311  if(var == "mapped" && symmetry == "phi-segmented")
312  {
313  vars.clear();
314  vars << magnitude;
315  int TNPOINTS, PNPOINTS, LNPOINTS;
316  double plimits[2], tlimits[2], llimits[2];
317  double mapOrigin[3];
318  string units[5];
319 
320  vars >> var ;
321  PNPOINTS = (int) get_number(var);
322  vars >> var1 >> var2 >> dim;
323  plimits[0] = get_number(var1 + "*" + dim);
324  plimits[1] = get_number(var2 + "*" + dim);
325  units[0].assign(dim);
326 
327  vars >> var ;
328  TNPOINTS = (int) get_number(var);
329  vars >> var1 >> var2 >> dim;
330  tlimits[0] = get_number(var1 + "*" + dim);
331  tlimits[1] = get_number(var2 + "*" + dim);
332  units[1].assign(dim);
333 
334  vars >> var ;
335  LNPOINTS = (int) get_number(var);
336  vars >> var1 >> var2 >> dim;
337  llimits[0] = get_number(var1 + "*" + dim);
338  llimits[1] = get_number(var2 + "*" + dim);
339  units[2].assign(dim);
340 
341  // Origin and Units
342  vars >> var1 >> var2 >> var3 >> dimSpc >> dimFld;
343  mapOrigin[0] = get_number(var1 + "*" + dimSpc);
344  mapOrigin[1] = get_number(var2 + "*" + dimSpc);
345  mapOrigin[2] = get_number(var3 + "*" + dimSpc);
346  units[3].assign( dimSpc );
347  units[4].assign( dimFld);
348 
349  vars.clear();
350  vars << swim_method;
351  vars >> integration_method;
352  vars.clear();
353  vars << swim_method;
354  vars >> integration_method;
355  double segm_phi_span = 2*(plimits[1] - plimits[0]); // factor of two: the map is assumed to cover half the segment
356  int nsegments = (int) (360.0/(segm_phi_span/deg));
357  if( fabs( segm_phi_span*nsegments/deg - 360 ) > 1.0e-2 )
358  {
359  cout << hd_msg << " <" << name << "> segmentation is wrong: " << nsegments << " segments, each span "
360  << segm_phi_span/deg << ": it doesn't add to 360, but " << segm_phi_span*nsegments/deg << " Exiting." << endl;
361  exit(0);
362 
363  }
364  if(MGN_VERBOSITY>3)
365  {
366  cout << hd_msg << " <" << name << "> is a phi-segmented mapped field in cylindrical coordinates" << endl;;
367  cout << hd_msg << " <" << name << "> has (" << PNPOINTS << ", " << TNPOINTS << ", " << LNPOINTS << ") points" << endl;;
368  cout << hd_msg << " Azimuthal Boundaries: (" << plimits[0]/deg << ", " << plimits[1]/deg << ") deg" << endl;
369  cout << hd_msg << " Tranverse Boundaries: (" << tlimits[0]/cm << ", " << tlimits[1]/cm << ") cm" << endl;
370  cout << hd_msg << " Longitudinal Boundaries: (" << llimits[0]/cm << ", " << llimits[1]/cm << ") cm" << endl;
371  cout << hd_msg << " Phi segment span, number of segments: " << segm_phi_span/deg << " deg, " << nsegments << endl;
372  cout << hd_msg << " Map Displacement: (" << mapOrigin[0]/cm << ", "
373  << mapOrigin[1]/cm << ", "
374  << mapOrigin[2]/cm << ") cm" << endl;
375  cout << hd_msg << " Integration Method: " << integration_method << endl;
376  cout << hd_msg << " Map Field Units: " << units[4] << endl;
377  cout << hd_msg << " Map File: " << MapFile << endl;
378  }
379  mappedfield = new MappedField(gemcOpt, TNPOINTS, PNPOINTS, LNPOINTS, tlimits, plimits, llimits, MapFile, mapOrigin, units, scale_factor);
380  mappedfield->segm_phi_span = segm_phi_span;
381  mappedfield->nsegments = nsegments;
382  }
383 
384 
385 
386 
387  // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
388  // Mapped y-symmetric Field in Cartesian coordinates
389  // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
390  if(var == "mapped" && symmetry == "y-symmetric")
391  {
392  vars.clear();
393  vars << magnitude;
394  int XNPOINTS, YNPOINTS, ZNPOINTS;
395  double xlimits[2], ylimits[2], zlimits[2];
396  double mapOrigin[3];
397  string units[5];
398 
399  vars >> var ;
400  XNPOINTS = (int) get_number(var);
401  vars >> var1 >> var2 >> dim;
402  xlimits[0] = get_number(var1 + "*" + dim);
403  xlimits[1] = get_number(var2 + "*" + dim);
404  units[0].assign(dim);
405 
406  vars >> var ;
407  YNPOINTS = (int) get_number(var);
408  vars >> var1 >> var2 >> dim;
409  ylimits[0] = get_number(var1 + "*" + dim);
410  ylimits[1] = get_number(var2 + "*" + dim);
411  units[1].assign(dim);
412 
413  vars >> var ;
414  ZNPOINTS = (int) get_number(var);
415  vars >> var1 >> var2 >> dim;
416  zlimits[0] = get_number(var1 + "*" + dim);
417  zlimits[1] = get_number(var2 + "*" + dim);
418  units[2].assign(dim);
419 
420  // Origin and Units
421  vars >> var1 >> var2 >> var3 >> dimSpc >> dimFld;
422  mapOrigin[0] = get_number(var1 + "*" + dimSpc);
423  mapOrigin[1] = get_number(var2 + "*" + dimSpc);
424  mapOrigin[2] = get_number(var3 + "*" + dimSpc);
425  units[3].assign( dimSpc );
426  units[4].assign( dimFld);
427 
428  vars.clear();
429  vars << swim_method;
430  vars >> integration_method;
431  vars.clear();
432  vars << swim_method;
433  vars >> integration_method;
434 
435  if(MGN_VERBOSITY>3)
436  {
437  cout << hd_msg << " <" << name << "> is a cartesian y-symmetric mapped field in cartesian coordinates" << endl;;
438  cout << hd_msg << " <" << name << "> has (" << XNPOINTS << ", " << YNPOINTS << ", " << ZNPOINTS << ") points" << endl;;
439  cout << hd_msg << " X Boundaries: (" << xlimits[0]/cm << ", " << xlimits[1]/cm << ") cm" << endl;
440  cout << hd_msg << " Y Boundaries: (" << ylimits[0]/cm << ", " << ylimits[1]/cm << ") cm" << endl;
441  cout << hd_msg << " Z Boundaries: (" << zlimits[0]/cm << ", " << zlimits[1]/cm << ") cm" << endl;
442  cout << hd_msg << " Map Displacement: ("
443  << mapOrigin[0]/cm << ", " << mapOrigin[1]/cm << ", " << mapOrigin[2]/cm << ") cm" << endl;
444  cout << hd_msg << " Integration Method: " << integration_method << endl;
445  cout << hd_msg << " Map Field Units: `" << units[4] << "`" << endl;
446  cout << hd_msg << " Map File: `" << MapFile << "`" << endl;
447  }
448  mappedfield = new MappedField(gemcOpt, YNPOINTS, XNPOINTS, ZNPOINTS, ylimits, xlimits, zlimits,
449  MapFile, mapOrigin, units, scale_factor, true);
450  mappedfield->nsegments = 2;
451  }
452 
453 
454  if(integration_method == "RungeKutta")
455  {
456  // specialized equations for mapped magnetic field
457 
458  G4Mag_UsualEqRhs* iEquation = new G4Mag_UsualEqRhs(mappedfield);
459  G4MagIntegratorStepper* iStepper = new G4HelixSimpleRunge(iEquation);
460  G4ChordFinder* iChordFinder = new G4ChordFinder(mappedfield, 1.0e-2*mm, iStepper);
461 
462  MFM = new G4FieldManager(mappedfield, iChordFinder);
463  MFM->SetDeltaOneStep(1*mm);
464  MFM->SetDeltaIntersection(1*mm);
465  }
466 }
467 
468 MappedField::MappedField(){;}
469 MappedField::~MappedField(){;}
470 
471 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
472 // Constructor for Dipole field in cartesian coordinate
473 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
474 MappedField::MappedField(gemc_opts Opt, int TNPOINTS, int LNPOINTS, double tlimits[2], double llimits[2],
475  string filename, double origin[3], string units[5], double scale_factor, string symmetry)
476 {
477  gemcOpt = Opt;
478  Symmetry = symmetry;
479  string hd_msg = gemcOpt.args["LOG_MSG"].args + " Magnetic Field Constructor: >> ";
480  MGN_VERBOSITY = (int) gemcOpt.args["MGN_VERBOSITY"].arg ;
481  string field_dir = gemcOpt.args["FIELD_DIR"].args ;
482  if( getenv("JLAB_ROOT") == NULL){
483  cout << hd_msg << " ERROR: Please set JLAB_ROOT to point to the directory containing noarch/data \n";
484  exit(EXIT_FAILURE);
485  }
486  string field_file = (string) getenv("JLAB_ROOT") + "/noarch/data/" + filename;
487 
488  if(field_dir != "env")
489  field_file = field_dir + "/" + filename;
490 
491  mapOrigin[0] = origin[0];
492  mapOrigin[1] = origin[1];
493  mapOrigin[2] = origin[2];
494 
495  table_start[0] = tlimits[0];
496  table_start[1] = llimits[0];
497  cell_size[0] = (tlimits[1] - tlimits[0])/( TNPOINTS - 1);
498  cell_size[1] = (llimits[1] - llimits[0])/( LNPOINTS - 1);
499 
500  if(MGN_VERBOSITY>3)
501  {
502  cout << hd_msg << " The Transverse cell size is: " << cell_size[0]/cm << " cm" << endl
503  << hd_msg << " and the Longitudinal cell size is: " << cell_size[1]/cm << " cm" << endl;
504 
505  }
506 
507  ifstream IN(field_file.c_str());
508  if(!IN)
509  {
510  cout << hd_msg << " File " << field_file << " could not be opened. Exiting." << endl;
511  exit(0);
512  }
513  cout << hd_msg << " Reading map file: " << field_file << "... ";
514 
515  // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
516  // Setting up storage space for tables
517  // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
518  BDipole.resize( LNPOINTS );
519  for(int il=0; il<LNPOINTS; il++)
520  {
521  BDipole[il].resize(TNPOINTS);
522  }
523 
524  // %%%%%%%%%%%%%
525  // Filling table
526  // %%%%%%%%%%%%%
527  double TC, LC; // coordinates as read from the map
528  double B; // magnetic field as read from the map
529  unsigned int IT, IL; // indexes of the vector map
530 
531  for(int il=0; il<LNPOINTS; il++)
532  {
533  for(int it=0; it<TNPOINTS; it++)
534  {
535  IN >> LC >> TC >> B;
536  if(units[0] == "cm") LC *= cm;
537  else if(units[0] == "m") LC *= m;
538  else if(units[0] == "mm") LC *= mm;
539  else cout << hd_msg << " Dimension Unit " << units[1] << " not implemented yet." << endl;
540  if(units[1] == "cm") TC *= cm;
541  else if(units[1] == "m") TC *= m;
542  else if(units[1] == "mm") TC *= mm;
543  else cout << hd_msg << " Dimension Unit " << units[0] << " not implemented yet." << endl;
544 
545  // checking map consistency for transverse coordinate
546  if( (tlimits[0] + it*cell_size[0] - TC)/TC > 0.001)
547  {
548  cout << hd_msg << it << " transverse index wrong. Map point should be " << tlimits[0] + it*cell_size[0]
549  << " but it's " << TC << " instead." << endl;
550  }
551 
552  IT = (unsigned int) floor( ( TC/mm - table_start[0]/mm + cell_size[0]/mm/2 ) / ( cell_size[0]/mm ) ) ;
553  IL = (unsigned int) floor( ( LC/mm - table_start[1]/mm + cell_size[1]/mm/2 ) / ( cell_size[1]/mm ) ) ;
554 
555  B *= scale_factor;
556 
557  if(units[4] == "gauss")
558  {
559  BDipole[IL][IT] = B*gauss;
560  }
561  else if(units[4] == "kilogauss")
562  {
563  BDipole[IL][IT] = B*kilogauss;
564  }
565  else if(units[4] == "T")
566  {
567  BDipole[IL][IT] = B*tesla;
568  }
569  else
570  {
571  cout << hd_msg << " Field Unit " << units[4] << " not implemented yet." << endl;
572  }
573  }
574  // checking map consistency for longitudinal coordinate
575  if( (llimits[0] + il*cell_size[1] - LC)/LC > 0.001)
576  {
577  cout << hd_msg << il << " longitudinal index wrong. Map point should be " << llimits[0] + il*cell_size[1]
578  << " but it's " << LC << " instead." << endl;
579  }
580  }
581 
582  IN.close();
583  cout << " done!" << endl;
584 
585 }
586 
587 
588 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
589 // Constructor for phi-symmetric field in cylindrical coordinates
590 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
591 MappedField::MappedField(gemc_opts Opt, int TNPOINTS, int LNPOINTS, double tlimits[2], double llimits[2],
592  string filename, double origin[3], string units[5], double scale_factor, string symmetry, int dummy)
593 {
594  gemcOpt = Opt;
595  string hd_msg = gemcOpt.args["LOG_MSG"].args + " Magnetic Field Constructor: >> ";
596  MGN_VERBOSITY = (int) gemcOpt.args["MGN_VERBOSITY"].arg ;
597  string field_dir = gemcOpt.args["FIELD_DIR"].args ;
598  if( getenv("JLAB_ROOT") == NULL){
599  cout << hd_msg << " ERROR: Please set JLAB_ROOT to point to the directory containing noarch/data \n";
600  exit(EXIT_FAILURE);
601  }
602  string field_file = (string) getenv("JLAB_ROOT") + "/noarch/data/" + filename;
603  Symmetry = symmetry;
604 
605  if(field_dir != "env")
606  field_file = field_dir + "/" + filename;
607 
608 
609  mapOrigin[0] = origin[0];
610  mapOrigin[1] = origin[1];
611  mapOrigin[2] = origin[2];
612 
613  table_start[0] = tlimits[0];
614  table_start[1] = llimits[0];
615  cell_size[0] = (tlimits[1] - tlimits[0])/( TNPOINTS - 1);
616  cell_size[1] = (llimits[1] - llimits[0])/( LNPOINTS - 1);
617 
618  if(MGN_VERBOSITY>3)
619  {
620  cout << hd_msg << " The transverse cell size is: " << cell_size[0]/cm << " cm" << endl
621  << hd_msg << " and the longitudinal cell size is: " << cell_size[1]/cm << " cm" << endl;
622 
623  }
624 
625  ifstream IN(field_file.c_str());
626  if(!IN)
627  {
628  cout << hd_msg << " File " << field_file << " could not be opened. Exiting." << endl;
629  exit(0);
630  }
631  cout << hd_msg << " Reading map file: " << field_file << "... ";
632 
633  // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
634  // Setting up storage space for tables
635  // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
636  B2DCylT.resize( TNPOINTS );
637  B2DCylL.resize( TNPOINTS );
638  for(int it=0; it<TNPOINTS; it++)
639  {
640  B2DCylT[it].resize(LNPOINTS);
641  B2DCylL[it].resize(LNPOINTS);
642  }
643 
644  // %%%%%%%%%%%%%
645  // Filling table
646  // %%%%%%%%%%%%%
647  double TC, LC; // coordinates as read from the map
648  double BT, BL; // transverse and longitudinal magnetic field as read from the map
649  unsigned int IT, IL; // indexes of the vector map
650 
651  for(int it=0; it<TNPOINTS; it++)
652  {
653  for(int il=0; il<LNPOINTS; il++)
654  {
655  IN >> TC >> LC >> BT >> BL;
656  if(units[0] == "cm") TC *= cm;
657  else if(units[0] == "m") TC *= m;
658  else if(units[0] == "mm") TC *= mm;
659  else cout << hd_msg << " Dimension Unit " << units[0] << " not implemented yet." << endl;
660  if(units[1] == "cm") LC *= cm;
661  else if(units[1] == "m") LC *= m;
662  else if(units[1] == "mm") LC *= mm;
663  else cout << hd_msg << " Dimension Unit " << units[1] << " not implemented yet." << endl;
664 
665  // checking map consistency for longitudinal coordinate
666  if( (llimits[0] + il*cell_size[1] - LC)/LC > 0.001)
667  {
668  cout << hd_msg << il << " longitudinal index wrong. Map point should be " << llimits[0] + il*cell_size[1]
669  << " but it's " << LC << " instead." << endl;
670  }
671 
672  IT = (unsigned int) floor( ( TC/mm - table_start[0]/mm + cell_size[0]/mm/2 ) / ( cell_size[0]/mm ) ) ;
673  IL = (unsigned int) floor( ( LC/mm - table_start[1]/mm + cell_size[1]/mm/2 ) / ( cell_size[1]/mm ) ) ;
674 
675  BT *= scale_factor;
676  BL *= scale_factor;
677 
678  if(units[4] == "gauss")
679  {
680  B2DCylT[IT][IL] = BT*gauss;
681  B2DCylL[IT][IL] = BL*gauss;
682  }
683  else if(units[4] == "kilogauss")
684  {
685  B2DCylT[IT][IL] = BT*kilogauss;
686  B2DCylL[IT][IL] = BL*kilogauss;
687  }
688  else if(units[4] == "T")
689  {
690  B2DCylT[IT][IL] = BT*tesla;
691  B2DCylL[IT][IL] = BL*tesla;
692  }
693  else
694  {
695  cout << hd_msg << " Field Unit " << units[4] << " not implemented yet." << endl;
696  }
697  }
698  // checking map consistency for transverse coordinate
699  if( (tlimits[0] + it*cell_size[0] - TC)/TC > 0.001)
700  {
701  cout << hd_msg << it << " transverse index wrong. Map point should be " << tlimits[0] + it*cell_size[0]
702  << " but it's " << TC << " instead." << endl;
703  }
704  }
705 
706  IN.close();
707  cout << " done!" << endl;
708 
709 }
710 
711 
712 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
713 // Constructor for phi-segmented field in cylindrical coordinates. Field is in cartesian coordinates
714 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
715 MappedField::MappedField(gemc_opts Opt, int rNPOINTS, int pNPOINTS, int zNPOINTS, double tlimits[2], double plimits[2], double llimits[2],
716  string filename, double origin[3], string units[5], double scale_factor)
717 {
718  gemcOpt = Opt;
719  string hd_msg = gemcOpt.args["LOG_MSG"].args + " Magnetic Field Constructor: >> ";
720  MGN_VERBOSITY = (int) gemcOpt.args["MGN_VERBOSITY"].arg ;
721  string field_dir = gemcOpt.args["FIELD_DIR"].args ;
722  if( getenv("JLAB_ROOT") == NULL){
723  cout << hd_msg << " ERROR: Please set JLAB_ROOT to point to the directory containing noarch/data \n";
724  exit(EXIT_FAILURE);
725  }
726  string field_file = (string) getenv("JLAB_ROOT") + "/noarch/data/" + filename;
727 
728  if(field_dir != "env")
729  field_file = field_dir + "/" + filename;
730 
731  mapOrigin[0] = origin[0];
732  mapOrigin[1] = origin[1];
733  mapOrigin[2] = origin[2];
734 
735  table_start[0] = plimits[0];
736  table_start[1] = tlimits[0];
737  table_start[2] = llimits[0];
738  cell_size[0] = (plimits[1] - plimits[0])/( pNPOINTS - 1);
739  cell_size[1] = (tlimits[1] - tlimits[0])/( rNPOINTS - 1);
740  cell_size[2] = (llimits[1] - llimits[0])/( zNPOINTS - 1);
741 
742  if(MGN_VERBOSITY>3)
743  {
744  cout << hd_msg << " the phi cell size is: " << cell_size[0]/deg << " degrees" << endl
745  << hd_msg << " The radius cell size is: " << cell_size[1]/cm << " cm" << endl
746  << hd_msg << " and the z cell size is: " << cell_size[2]/cm << " cm" << endl;
747  }
748 
749  ifstream IN(field_file.c_str());
750  if(!IN)
751  {
752  cout << hd_msg << " File " << field_file << " could not be opened. Exiting." << endl;
753  exit(0);
754  }
755  cout << hd_msg << " Reading map file: " << field_file << "... ";
756 
757 
758  // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
759  // Setting up storage space for tables
760  // Field[phi][r][z]
761  // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
762 
763  B3DCylX.resize( pNPOINTS );
764  B3DCylY.resize( pNPOINTS );
765  B3DCylZ.resize( pNPOINTS );
766  for(int ip=0; ip<pNPOINTS; ip++)
767  {
768  B3DCylX[ip].resize( rNPOINTS );
769  B3DCylY[ip].resize( rNPOINTS );
770  B3DCylZ[ip].resize( rNPOINTS );
771  for(int ir=0; ir<rNPOINTS; ir++)
772  {
773  B3DCylX[ip][ir].resize( zNPOINTS );
774  B3DCylY[ip][ir].resize( zNPOINTS );
775  B3DCylZ[ip][ir].resize( zNPOINTS );
776  }
777  }
778 
779  // %%%%%%%%%%%%%
780  // Filling table
781  // %%%%%%%%%%%%%
782  double pC, tC, lC; // coordinates as read from the map
783  double Bx, By, Bz; // magnetic field components as read from the map
784  unsigned int It, Ip, Il; // indexes of the vector map
785  for(int ip=0; ip<pNPOINTS; ip++)
786  {
787  for(int it=0; it<rNPOINTS; it++)
788  {
789  for(int il=0; il<zNPOINTS; il++)
790  {
791  IN >> pC >> tC >> lC >> Bx >> By >> Bz;
792 
793  if(units[0] == "deg") pC = pC*deg;
794  else cout << hd_msg << " Dimension Unit " << units[0] << " not implemented yet." << endl;
795  if(units[1] == "cm") tC *= cm;
796  else cout << hd_msg << " Dimension Unit " << units[1] << " not implemented yet." << endl;
797  if(units[2] == "cm") lC *= cm;
798  else cout << hd_msg << " Dimension Unit " << units[2] << " not implemented yet." << endl;
799 
800  // checking map consistency for longitudinal coordinate
801  if( (llimits[0] + il*cell_size[2] - lC)/lC > 0.001)
802  {
803  cout << hd_msg << il << " longitudinal index wrong. Map point should be " << llimits[0] + il*cell_size[2]
804  << " but it's " << lC << " instead." << endl;
805  }
806 
807  Ip = (unsigned int) floor( ( pC/deg - table_start[0]/deg + cell_size[0]/deg/2 ) / ( cell_size[0]/deg ) ) ;
808  It = (unsigned int) floor( ( tC/mm - table_start[1]/mm + cell_size[1]/mm/2 ) / ( cell_size[1]/mm ) ) ;
809  Il = (unsigned int) floor( ( lC/mm - table_start[2]/mm + cell_size[2]/mm/2 ) / ( cell_size[2]/mm ) ) ;
810 
811  Bx *= scale_factor;
812  By *= scale_factor;
813  Bz *= scale_factor;
814 
815  if(units[4] == "gauss")
816  {
817  B3DCylX[Ip][It][Il] = Bx*gauss;
818  B3DCylY[Ip][It][Il] = By*gauss;
819  B3DCylZ[Ip][It][Il] = Bz*gauss;
820  }
821  else if(units[4] == "kilogauss")
822  {
823  B3DCylX[Ip][It][Il] = Bx*kilogauss;
824  B3DCylY[Ip][It][Il] = By*kilogauss;
825  B3DCylZ[Ip][It][Il] = Bz*kilogauss;
826  }
827  else if(units[4] == "T")
828  {
829  B3DCylX[Ip][It][Il] = Bx*tesla;
830  B3DCylY[Ip][It][Il] = By*tesla;
831  B3DCylZ[Ip][It][Il] = Bz*tesla;
832  }
833 
834  else
835  {
836  cout << hd_msg << " Field Unit " << units[4] << " not implemented yet." << endl;
837  }
838 
839  if(MGN_VERBOSITY>10 && MGN_VERBOSITY < 50)
840  {
841  cout << " phi=" << pC/deg << ", ip=" << Ip << " "
842  << " r=" << tC/cm << ", it=" << It << " "
843  << " z=" << lC/cm << ", iz=" << Il << " "
844  << " Bx=" << B3DCylX[Ip][It][Il]/kilogauss << " "
845  << " By=" << B3DCylY[Ip][It][Il]/kilogauss << " "
846  << " Bz=" << B3DCylZ[Ip][It][Il]/kilogauss << " kilogauss. Map Values: "
847  << " rBx=" << Bx << " "
848  << " rBy=" << By << " "
849  << " rBz=" << Bz << endl;
850 
851  }
852 
853  }
854 
855  // checking map consistency for transverse coordinate
856  if( (tlimits[0] + it*cell_size[1] - tC)/tC > 0.001)
857  {
858  cout << hd_msg << it << " transverse index wrong. Map point should be " << tlimits[0] + it*cell_size[0]
859  << " but it's " << tC << " instead." << endl;
860  }
861 
862  }
863 
864  // checking map consistency for azimuthal coordinate
865  if( (plimits[0] + ip*cell_size[0] - pC)/pC > 0.001)
866  {
867  cout << hd_msg << ip << " azimuthal index wrong. Map point should be " << plimits[0] + ip*cell_size[1]
868  << " but it's " << pC << " instead." << endl;
869  }
870 
871  }
872 
873  IN.close();
874  cout << " done!" << endl;
875 
876 }
877 
878 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
879 // Constructor for y-symmetric field in cartesian coordinates. Field is in cartesian coordinates
880 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
881 MappedField::MappedField(gemc_opts Opt, int yNPOINTS, int xNPOINTS, int zNPOINTS,
882  double ylimits[2], double xlimits[2], double zlimits[2],
883  string filename, double origin[3], string units[5], double scale_factor, bool ySymmetry)
884 {
885  gemcOpt = Opt;
886  string hd_msg = gemcOpt.args["LOG_MSG"].args + " Magnetic Field Constructor: >> ";
887  MGN_VERBOSITY = (int) gemcOpt.args["MGN_VERBOSITY"].arg ;
888  string field_dir = gemcOpt.args["FIELD_DIR"].args ;
889  if( getenv("JLAB_ROOT") == NULL){
890  cout << hd_msg << " ERROR: Please set JLAB_ROOT to point to the directory containing noarch/data \n";
891  exit(EXIT_FAILURE);
892  }
893  string field_file = (string) getenv("JLAB_ROOT") + "/noarch/data/" + filename;
894 
895  if(field_dir != "env")
896  field_file = field_dir + "/" + filename;
897 
898  double table_end[3];
899 
900  if( ySymmetry ) nsegments = 2;
901 
902  mapOrigin[0] = origin[0];
903  mapOrigin[1] = origin[1];
904  mapOrigin[2] = origin[2];
905 
906  table_start[0] = xlimits[0];
907  table_start[1] = ylimits[0];
908  table_start[2] = zlimits[0];
909 
910  table_end[0] = xlimits[1];
911  table_end[1] = ylimits[1];
912  table_end[2] = zlimits[1];
913 
914  if( ySymmetry )
915  {
916  table_start[1] = fabs( table_start[1] );
917  table_end[1] = fabs( table_end[1] );
918  }
919 
920  cell_size[0] = (table_end[0] - table_start[0])/( xNPOINTS - 1);
921  cell_size[1] = (table_end[1] - table_start[1])/( yNPOINTS - 1);
922  cell_size[2] = (table_end[2] - table_start[2])/( zNPOINTS - 1);
923 
924 
925  if(MGN_VERBOSITY>3)
926  {
927  cout << hd_msg << " the X cell size is: " << cell_size[0]/cm << " cm" << endl
928  << hd_msg << " The Y cell size is: " << cell_size[1]/cm << " cm" << endl
929  << hd_msg << " and the Z cell size is: " << cell_size[2]/cm << " cm" << endl;
930  }
931 
932  ifstream IN(field_file.c_str());
933  if(!IN)
934  {
935  cout << hd_msg << " File " << field_file << " could not be opened. Exiting." << endl;
936  exit(0);
937  }
938  cout << hd_msg << " Reading map file: " << field_file << "... ";
939 
940 
941  // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
942  // Setting up storage space for tables
943  // Field[x][y][z]
944  // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
945 
946  B3DCartX.resize( xNPOINTS );
947  B3DCartY.resize( xNPOINTS );
948  B3DCartZ.resize( xNPOINTS );
949  for(int ix=0; ix<xNPOINTS; ix++)
950  {
951  B3DCartX[ix].resize( yNPOINTS );
952  B3DCartY[ix].resize( yNPOINTS );
953  B3DCartZ[ix].resize( yNPOINTS );
954  for(int iy=0; iy<yNPOINTS; iy++)
955  {
956  B3DCartX[ix][iy].resize( zNPOINTS );
957  B3DCartY[ix][iy].resize( zNPOINTS );
958  B3DCartZ[ix][iy].resize( zNPOINTS );
959  }
960  }
961 
962  // %%%%%%%%%%%%%
963  // Filling table
964  // %%%%%%%%%%%%%
965  double xC, yC, zC; // coordinates as read from the map
966  double Bx, By, Bz; // magnetic field components as read from the map
967  unsigned int Ix, Iy, Iz; // indexes of the vector map
968  for(int ix=0; ix<xNPOINTS; ix++)
969  {
970  for(int iz=0; iz<zNPOINTS; iz++)
971  {
972  for(int iy=0; iy<yNPOINTS; iy++)
973  {
974  IN >> xC >> yC >> zC >> Bx >> By >> Bz;
975  if( ySymmetry ) yC = fabs( yC );
976 
977  if(units[0] == "cm") xC *= cm;
978  else if(units[0] == "m") xC *= m;
979  else if(units[0] == "mm") xC *= mm;
980  else cout << hd_msg << " Dimension Unit " << units[0] << " not implemented yet." << endl;
981  if(units[1] == "cm") yC *= cm;
982  else if(units[1] == "m") yC *= m;
983  else if(units[1] == "mm") yC *= mm;
984  else cout << hd_msg << " Dimension Unit " << units[1] << " not implemented yet." << endl;
985  if(units[2] == "cm") zC *= cm;
986  else if(units[2] == "m") zC *= m;
987  else if(units[2] == "mm") zC *= mm;
988  else cout << hd_msg << " Dimension Unit " << units[2] << " not implemented yet." << endl;
989 
990  // checking map consistency for Y coordinate
991  if( (table_start[1] + iy*cell_size[1] - yC)/yC > 0.001)
992  {
993  cout << hd_msg << iy << " Y index wrong. Map point should be " << table_start[1] + iy*cell_size[1]
994  << " but it's " << yC << " instead." << endl;
995  }
996 
997  Ix = (unsigned int) floor( ( xC/mm - table_start[0]/mm + cell_size[0]/mm/2 ) / ( cell_size[0]/mm ) ) ;
998  Iy = (unsigned int) floor( ( yC/mm - table_start[1]/mm + cell_size[1]/mm/2 ) / ( cell_size[1]/mm ) ) ;
999  Iz = (unsigned int) floor( ( zC/mm - table_start[2]/mm + cell_size[2]/mm/2 ) / ( cell_size[2]/mm ) ) ;
1000  Bx *= scale_factor;
1001  By *= scale_factor;
1002  Bz *= scale_factor;
1003 
1004  if(units[4] == "gauss")
1005  {
1006  B3DCartX[Ix][Iy][Iz] = Bx*gauss;
1007  B3DCartY[Ix][Iy][Iz] = By*gauss;
1008  B3DCartZ[Ix][Iy][Iz] = Bz*gauss;
1009  } else if(units[4] == "kilogauss")
1010  {
1011  B3DCartX[Ix][Iy][Iz] = Bx*kilogauss;
1012  B3DCartY[Ix][Iy][Iz] = By*kilogauss;
1013  B3DCartZ[Ix][Iy][Iz] = Bz*kilogauss;
1014  } else if(units[4] == "T")
1015  {
1016  B3DCartX[Ix][Iy][Iz] = Bx*tesla;
1017  B3DCartY[Ix][Iy][Iz] = By*tesla;
1018  B3DCartZ[Ix][Iy][Iz] = Bz*tesla;
1019  }
1020  else
1021  {
1022  cout << hd_msg << " Field Unit " << units[4] << " not implemented yet." << endl;
1023  }
1024 
1025  if(MGN_VERBOSITY>10 && MGN_VERBOSITY < 50)
1026  {
1027  cout << " x=" << xC/cm << ", ix=" << Ix << " "
1028  << " y=" << yC/cm << ", iy=" << Iy << " "
1029  << " z=" << zC/cm << ", iz=" << Iz << " "
1030  << " Bx=" << B3DCartX[Ix][Iy][Iz]/kilogauss << " "
1031  << " By=" << B3DCartY[Ix][Iy][Iz]/kilogauss << " "
1032  << " Bz=" << B3DCartZ[Ix][Iy][Iz]/kilogauss << " kilogauss. Map Values: "
1033  << " rBx=" << Bx << " "
1034  << " rBy=" << By << " "
1035  << " rBz=" << Bz << endl;
1036  }
1037  }
1038 
1039  // checking map consistency for Z coordinate
1040  if( (table_start[2] + iz*cell_size[2] - zC)/zC > 0.001)
1041  {
1042  cout << hd_msg << iz << " Z index wrong. Map point should be " << table_start[2] + iz*cell_size[2]
1043  << " but it's " << zC << " instead." << endl;
1044  }
1045 
1046  }
1047 
1048  // checking map consistency for X coordinate
1049  if( (table_start[0] + ix*cell_size[0] - xC)/xC > 0.001)
1050  {
1051  cout << hd_msg << ix << " X index wrong. Map point should be " << table_start[0] + ix*cell_size[0]
1052  << " but it's " << xC << " instead." << endl;
1053  }
1054  }
1055 
1056  IN.close();
1057  cout << " done!" << endl;
1058 
1059 }
1060 
1061 
1062 
1063 
1064 // %%%%%%%%%%%%%
1065 // GetFieldValue
1066 // %%%%%%%%%%%%%
1067 void MappedField::GetFieldValue(const double point[3], double *Bfield) const
1068 {
1069 
1070  double Point[3]; // global coordinates, in mm, after shifting the origin
1071  vector<double> Field[3];
1072  static int FIRST_ONLY;
1073 
1074  Point[0] = point[0] - mapOrigin[0]/mm;
1075  Point[1] = point[1] - mapOrigin[1]/mm;
1076  Point[2] = point[2] - mapOrigin[2]/mm;
1077 
1078  Bfield[0] = Bfield[1] = Bfield[2] = 0*gauss;
1079 
1080  // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1081  // Dipole field in cartesian coordinates
1082  // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1083  if(BDipole.size())
1084  {
1085  double psfield[3] = {0,0,0};
1086  unsigned int IT, IL;
1087 
1088  double LC; // longitudinal coordinate of the track in the global coordinate system
1089  double TC; // transverse and phy polar 2D coordinates: the map is phi-symmetric. phi is angle in respect to x
1090  if(Symmetry == "Dipole-x")
1091  {
1092  TC = Point[1];
1093  LC = Point[2];
1094  }
1095  if(Symmetry == "Dipole-y")
1096  {
1097  TC = Point[0];
1098  LC = Point[2];
1099  }
1100  if(Symmetry == "Dipole-z")
1101  {
1102  TC = Point[0];
1103  LC = Point[1];
1104  }
1105 
1106 
1107  IT = (unsigned int) floor( ( TC - table_start[0]/mm ) / (cell_size[0]/mm) ) ;
1108  IL = (unsigned int) floor( ( LC - table_start[1]/mm ) / (cell_size[1]/mm) ) ;
1109 
1110  if( fabs( table_start[0]/mm + IT*cell_size[0]/mm - TC) > fabs( table_start[0]/mm + (IT+1)*cell_size[0]/mm - TC) ) IT++;
1111  if( fabs( table_start[1]/mm + IL*cell_size[1]/mm - LC) > fabs( table_start[1]/mm + (IL+1)*cell_size[1]/mm - LC) ) IL++;
1112  // vector sizes are checked on both T and L components
1113  // (even if only one is enough)
1114  if(IL < BDipole.size())
1115  if(IT < BDipole[IL].size())
1116  {
1117  if(Symmetry == "Dipole-x") psfield[0] = BDipole[IL][IT];
1118  if(Symmetry == "Dipole-y") psfield[1] = BDipole[IL][IT];
1119  if(Symmetry == "Dipole-z") psfield[2] = BDipole[IL][IT];
1120  if(MGN_VERBOSITY>5 && FIRST_ONLY != 99)
1121  {
1122  cout << hd_msg << " Dipole Field: Cart. coordinates (cm), table indexes, magnetic field values (gauss):" << endl;
1123  cout << " x=" << point[0]/cm << " " ;
1124  cout << "y=" << point[1]/cm << " ";
1125  cout << "z=" << point[2]/cm << " ";
1126  cout << "IT=" << IT << " ";
1127  cout << "IL=" << IL << " ";
1128  cout << "Bx=" << psfield[0]/gauss << " ";
1129  cout << "By=" << psfield[1]/gauss << " ";
1130  cout << "Bz=" << psfield[2]/gauss << endl;
1131  }
1132  }
1133  for(int i=0; i<3; i++) Field[i].push_back(psfield[i]);
1134  }
1135 
1136 
1137  // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1138  // phi symmetric field in cylindrical coordinates
1139  // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1140  if(B2DCylT.size() && B2DCylL.size())
1141  {
1142  double psfield[3] = {0,0,0};
1143  unsigned int IT, IL;
1144 
1145  double LC; // longitudinal coordinate of the track in the global coordinate system
1146  double TC, phi; // transverse and phy polar 2D coordinates: the map is phi-symmetric. phi is angle in respect to x
1147  TC = sqrt(Point[0]*Point[0] + Point[1]*Point[1]);
1148  G4ThreeVector vpos(Point[0], Point[1], Point[2]);
1149  phi = vpos.phi();
1150  LC = Point[2];
1151 
1152  IT = (unsigned int) floor( ( TC - table_start[0]/mm ) / (cell_size[0]/mm) ) ;
1153  IL = (unsigned int) floor( ( LC - table_start[1]/mm ) / (cell_size[1]/mm) ) ;
1154 
1155  if( fabs( table_start[0]/mm + IT*cell_size[0]/mm - TC) > fabs( table_start[0]/mm + (IT+1)*cell_size[0]/mm - TC) ) IT++;
1156  if( fabs( table_start[1]/mm + IL*cell_size[1]/mm - LC) > fabs( table_start[1]/mm + (IL+1)*cell_size[1]/mm - LC) ) IL++;
1157 
1158 
1159  // vector sizes are checked on both T and L components
1160  // (even if only one is enough)
1161  if(IT < B2DCylT.size() && IT < B2DCylL.size())
1162  if(IL < B2DCylT[IT].size() && IL < B2DCylL[IT].size() )
1163  {
1164 
1165 
1166  psfield[0] = B2DCylT[IT][IL] * cos(phi);
1167  psfield[1] = B2DCylT[IT][IL] * sin(phi);
1168  psfield[2] = B2DCylL[IT][IL];
1169 
1170 
1171  if(Symmetry == "cylindrical-x")
1172  {
1173  psfield[2] = B2DCylT[IT][IL] * cos(phi);
1174  psfield[1] = B2DCylT[IT][IL] * sin(phi);
1175  psfield[0] = B2DCylL[IT][IL];
1176  }
1177 
1178  if(Symmetry == "cylindrical-y")
1179  {
1180  psfield[0] = B2DCylT[IT][IL] * cos(phi);
1181  psfield[2] = B2DCylT[IT][IL] * sin(phi);
1182  psfield[1] = B2DCylL[IT][IL];
1183  }
1184 
1185  if(MGN_VERBOSITY>5 && FIRST_ONLY != 99)
1186  {
1187  cout << hd_msg << " Phi-Simmetric Field: Cart. and Cyl. coordinates (cm), table indexes, magnetic field values (gauss):" << endl;
1188  cout << " x=" << point[0]/cm << " " ;
1189  cout << "y=" << point[1]/cm << " ";
1190  cout << "z=" << point[2]/cm << " ";
1191  cout << "r=" << TC/cm << " ";
1192  cout << "z=" << LC/cm << " ";
1193  cout << "phi=" << phi/deg << " ";
1194  cout << "IT=" << IT << " ";
1195  cout << "IL=" << IL << " ";
1196  cout << "Bx=" << psfield[0]/gauss << " ";
1197  cout << "By=" << psfield[1]/gauss << " ";
1198  cout << "Bz=" << psfield[2]/gauss << endl;
1199  }
1200  }
1201  for(int i=0; i<3; i++) Field[i].push_back(psfield[i]);
1202  }
1203 
1204 
1205  // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1206  // phi-segmented field in cylindrical coordinates
1207  // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1208  if(B3DCylX.size() && B3DCylY.size() && B3DCylZ.size())
1209  {
1210  double mfield[3] = {0,0,0};
1211  double rotpsfield[3] = {0,0,0};
1212 
1213  double tC, pC, lC;
1214  double pLC;
1215  unsigned int Ip, It, Il;
1216 
1217  pC = atan2( Point[1], Point[0] )*rad;
1218  if( pC < 0 ) pC += 360*deg;
1219 
1220  tC = sqrt(Point[0]*Point[0] + Point[1]*Point[1]);
1221  lC = Point[2];
1222 
1223  // Rotating the point to within the map limit
1224  int segment = 0;
1225  pLC = pC;
1226  while (pLC/deg > 30)
1227  {
1228  pLC -= 60*deg;
1229  segment++;
1230  }
1231 
1232  double dphi = pC - pLC;
1233  int sign = (pLC >= 0 ? 1 : -1);
1234  double apLC = fabs(pLC);
1235 
1236  Ip = (unsigned int) floor( ( apLC/deg - table_start[0]/deg ) / (cell_size[0]/deg ) ) ;
1237  It = (unsigned int) floor( ( tC/mm - table_start[1]/mm ) / (cell_size[1]/mm ) ) ;
1238  Il = (unsigned int) floor( ( lC/mm - table_start[2]/mm ) / (cell_size[2]/mm ) ) ;
1239 
1240  if( fabs( table_start[0]/mm + Ip*cell_size[0]/deg - apLC) > fabs( table_start[0]/mm + (Ip+1)*cell_size[0]/mm - apLC) ) Ip++;
1241  if( fabs( table_start[1]/mm + It*cell_size[1]/mm - tC) > fabs( table_start[1]/mm + (It+1)*cell_size[1]/mm - tC) ) It++;
1242  if( fabs( table_start[2]/mm + Il*cell_size[2]/mm - lC) > fabs( table_start[2]/mm + (Il+1)*cell_size[2]/mm - lC) ) Il++;
1243 
1244  // Getting Field at the rotated point
1245  // vector sizes are checked on all components
1246  // (even if only one is enough)
1247  if( Ip < B3DCylX.size() && Ip < B3DCylY.size() && Ip < B3DCylZ.size())
1248  if( It < B3DCylX[Ip].size() && It < B3DCylY[Ip].size() && It < B3DCylZ[Ip].size())
1249  if( Il < B3DCylX[Ip][It].size() && Il < B3DCylY[Ip][It].size() && Il < B3DCylZ[Ip][It].size())
1250  {
1251  // Field at local point
1252  mfield[0] = B3DCylX[Ip][It][Il];
1253  mfield[1] = B3DCylY[Ip][It][Il];
1254  mfield[2] = B3DCylZ[Ip][It][Il];
1255 
1256 
1257  // Rotating the field back to original point
1258  rotpsfield[0] = sign*mfield[0] * cos(dphi/rad) - mfield[1] * sin(dphi/rad);
1259  rotpsfield[1] = +sign*mfield[0] * sin(dphi/rad) + mfield[1] * cos(dphi/rad);
1260  rotpsfield[2] = sign*mfield[2];
1261 
1262  if(MGN_VERBOSITY>6 && FIRST_ONLY != 99)
1263  {
1264  cout << hd_msg << " Phi-Segmented Field: Cart. and Cyl. coord. (cm), indexes, local and rotated field values (gauss):" << endl;
1265  cout << " x=" << point[0]/cm << " " ;
1266  cout << "y=" << point[1]/cm << " ";
1267  cout << "z=" << point[2]/cm << " ";
1268  cout << "phi=" << pC/deg << " ";
1269  cout << "r=" << tC/cm << " ";
1270  cout << "z=" << lC/cm << " ";
1271  cout << "dphi=" << dphi/deg << " ";
1272  cout << "dphir=" << dphi/rad << " ";
1273  cout << "lphi=" << (sign == 1 ? "+ " : "- ") << pLC/deg << " ";
1274  cout << "segment=" << segment << " ";
1275  cout << "Ip=" << Ip << " ";
1276  cout << "It=" << It << " ";
1277  cout << "Il=" << Il << " ";
1278  cout << "lBx=" << mfield[0]/gauss << " ";
1279  cout << "lBy=" << mfield[1]/gauss << " ";
1280  cout << "lBz=" << mfield[2]/gauss << " " ;
1281  cout << "rBx=" << rotpsfield[0]/gauss << " ";
1282  cout << "rBy=" << rotpsfield[1]/gauss << " ";
1283  cout << "rBz=" << rotpsfield[2]/gauss << endl;
1284  }
1285  }
1286  for(int i=0; i<3; i++) Field[i].push_back(-rotpsfield[i]);
1287  }
1288 
1289  // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1290  // y-symmtric field in Cartesian coordinates
1291  // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1292  if(B3DCartX.size() && B3DCartY.size() && B3DCartZ.size())
1293  {
1294  double mfield[3] = {0,0,0};
1295  double xC, yC, zC;
1296  unsigned int Ix, Iy, Iz;
1297 
1298  xC = Point[0];
1299  yC = Point[1];
1300  zC = Point[2];
1301 
1302  // Use positive y-values for y-symmtric field
1303  if( nsegments == 2 ) yC = fabs( yC );
1304 
1305  Ix = (unsigned int) floor( ( xC/mm - table_start[0]/mm ) / (cell_size[0]/mm ) ) ;
1306  Iy = (unsigned int) floor( ( yC/mm - table_start[1]/mm ) / (cell_size[1]/mm ) ) ;
1307  Iz = (unsigned int) floor( ( zC/mm - table_start[2]/mm ) / (cell_size[2]/mm ) ) ;
1308 
1309  if( fabs( table_start[0]/mm + Ix*cell_size[0]/mm - xC) > fabs( table_start[0]/mm + (Ix+1)*cell_size[0]/mm - xC) ) Ix++;
1310  if( fabs( table_start[1]/mm + Iy*cell_size[1]/mm - yC) > fabs( table_start[1]/mm + (Iy+1)*cell_size[1]/mm - yC) ) Iy++;
1311  if( fabs( table_start[2]/mm + Iz*cell_size[2]/mm - zC) > fabs( table_start[2]/mm + (Iz+1)*cell_size[2]/mm - zC) ) Iz++;
1312 
1313  // Getting Field at the point in space
1314  // vector sizes are checked on all components
1315  // (even if only one is enough)
1316  if( Ix < B3DCartX.size() && Ix < B3DCartY.size() && Ix < B3DCartZ.size())
1317  if( Iy < B3DCartX[Ix].size() && Iy < B3DCartY[Ix].size() && Iy < B3DCartZ[Ix].size())
1318  if( Iz < B3DCartX[Ix][Iy].size() && Iz < B3DCartY[Ix][Iy].size() && Iz < B3DCartZ[Ix][Iy].size())
1319  {
1320  // Field at local point
1321  mfield[0] = B3DCartX[Ix][Iy][Iz];
1322  mfield[1] = B3DCartY[Ix][Iy][Iz];
1323  mfield[2] = B3DCartZ[Ix][Iy][Iz];
1324 
1325  if(MGN_VERBOSITY>6 && FIRST_ONLY != 99)
1326  {
1327  cout << hd_msg << " Segmented Field: Cart. coord. (cm), indexes and field values (gauss):" << endl;
1328  cout << " x=" << point[0]/cm << " " ;
1329  cout << "y=" << point[1]/cm << " ";
1330  cout << "z=" << point[2]/cm << " ";
1331  cout << "Ix=" << Ix << " ";
1332  cout << "Iy=" << Iy << " ";
1333  cout << "Iz=" << Iz << " ";
1334  cout << "lBx=" << mfield[0]/gauss << " ";
1335  cout << "lBy=" << mfield[1]/gauss << " ";
1336  cout << "lBz=" << mfield[2]/gauss << endl ;
1337  }
1338  }
1339  for(int i=0; i<3; i++) Field[i].push_back(mfield[i]);
1340  }
1341 
1342 
1343  // %%%%%%%%%%%%%%%%%%
1344  // Summing the Fields
1345  // %%%%%%%%%%%%%%%%%%
1346  for(unsigned int i=0; i<Field[0].size(); i++)
1347  for(int j=0; j<3; j++)
1348  Bfield[j] += Field[j][i];
1349 
1350  if(MGN_VERBOSITY>5 && FIRST_ONLY != 99)
1351  {
1352  cout << hd_msg << " Total Field: coordinates (cm), magnetic field values (gauss):" << endl ;
1353  cout << " x=" << point[0]/cm << " " ;
1354  cout << "y=" << point[1]/cm << " ";
1355  cout << "z=" << point[2]/cm << " ";
1356  cout << "Bx=" << Bfield[0]/gauss << " ";
1357  cout << "By=" << Bfield[1]/gauss << " ";
1358  cout << "Bz=" << Bfield[2]/gauss << endl << endl;
1359  }
1360 
1361  if(MGN_VERBOSITY == 99)
1362  FIRST_ONLY = 99;
1363 
1364 }
1365 
1366 
1367 
1368 
1369 
1370 
1371 
1372 
1373 
1374 
1375 
1376 
1377 
MagneticField::scale_factor
double scale_factor
Scale factor.
Definition: MagneticField.h:154
MagneticField::create_MFM
void create_MFM()
Creates the G4 Magnetic Field Manager.
Definition: MagneticField.cc:134
MappedField::segm_phi_span
double segm_phi_span
azimuthal span of the map for phi-segmented fields
Definition: MagneticField.h:109
get_magnetic_Fields
map< string, MagneticField > get_magnetic_Fields(gemc_opts Opt)
Fills magnetic field maps from Database.
Definition: MagneticField.cc:33
gemc_tostring
string gemc_tostring(QString input)
Definition: string_utilities.h:42
MagneticField::init_MFM
void init_MFM()
Initialize Field Manager Pointer to NULL.
Definition: MagneticField.h:163
MappedField::nsegments
int nsegments
number of segments
Definition: MagneticField.h:110
MagneticField::magnitude
string magnitude
Magnetic Field magnitude infos.
Definition: MagneticField.h:151
MagneticField.h
string_utilities.h
MagneticField::swim_method
string swim_method
Magnetic Field Swim Method.
Definition: MagneticField.h:152
MagneticField::gemcOpt
gemc_opts gemcOpt
gemc options map
Definition: MagneticField.h:148
gemc_opts::args
map< string, opts > args
Options map.
Definition: usage.h:68
MagneticField::name
string name
Magnetic Field identifier.
Definition: MagneticField.h:149
MappedField::GetFieldValue
void GetFieldValue(const double x[3], double *Bfield) const
Definition: MagneticField.cc:1067
TrimSpaces
string TrimSpaces(string in)
Removes leading and trailing spaces.
Definition: string_utilities.cc:25
MagneticField::mappedfield
MappedField * mappedfield
Mapped Magnetic Field.
Definition: MagneticField.h:156
get_number
double get_number(string)
Returns dimension from string, i.e. 100*cm.
Definition: string_utilities.cc:65
gemc_opts
Definition: usage.h:59
gemc_opts::get_args
vector< opts > get_args(string)
get a vector of arguments matching a string
Definition: usage.cc:983
MagneticField::type
string type
Type of Magnetic Field.
Definition: MagneticField.h:150
MagneticField::description
string description
Field Description.
Definition: MagneticField.h:153
Generated on Wed Jul 6 2016 11:12:55 for GEMC by   doxygen 1.8.11