source: trunk/examples/novice/gemc/src/MagneticField.cc @ 1282

Last change on this file since 1282 was 807, checked in by garnier, 16 years ago

update
File size: 26.0 KB
Line 
1
2
3
4
5
6
7
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
32map<string, MagneticField> get_magnetic_Fields(gemc_opts Opt)
33{
34 string hd_msg        = Opt.args["LOG_MSG"].args + " Magnetic Field >> " ;
35 double MGN_VERBOSITY = Opt.args["MGN_VERBOSITY"].arg;
36 string database      = Opt.args["DATABASE"].args;
37 string dbtable       = "magnetic_fields";
38
39 map<string, MagneticField> FieldMap;
40
41 QSqlDatabase db = QSqlDatabase::addDatabase("QMYSQL");
42 db.setHostName("clasdb.jlab.org");
43 db.setDatabaseName(database.c_str());
44 db.setUserName("clasuser");
45 bool ok = db.open();
46
47 if(!ok)
48 {
49    cout << hd_msg << " Database not connected. Exiting." << endl;
50    exit(-1);
51 }
52
53 MagneticField magneticField;
54 QSqlQuery q;
55 string dbexecute = "select name, type, magnitude, swim_method, description from " + dbtable ;
56 q.exec(dbexecute.c_str());
57
58 if(MGN_VERBOSITY>2)
59    cout << hd_msg << " Available Magnetic Fields: " << endl  << endl;
60
61 while (q.next())
62 {
63    magneticField.name        = TrimSpaces(q.value(0).toString().toStdString());
64    magneticField.type        = q.value(1).toString().toStdString();
65    magneticField.magnitude   = q.value(2).toString().toStdString();
66    magneticField.swim_method = q.value(3).toString().toStdString();
67    magneticField.description = q.value(4).toString().toStdString();
68
69    // Sets MFM pointer to NULL
70    magneticField.init_MFM();
71    magneticField.gemcOpt = Opt;
72
73    FieldMap[magneticField.name] = magneticField;
74
75    if(MGN_VERBOSITY>2)
76    {
77       cout << "      ";
78       cout.width(15);
79       cout << magneticField.name << "   |   ";
80       cout << magneticField.description << endl;
81       cout << "      ";
82       cout.width(15);
83       cout << magneticField.name << "   |  type:        |  ";
84       cout << magneticField.type << endl;
85       cout << "      ";
86       cout.width(15);
87       cout << magneticField.name << "   |  Magnitude:   |  ";
88       cout << magneticField.magnitude << endl;
89       cout << "      ";
90       cout.width(15);
91       cout << magneticField.name << "   |  Swim Method: |  ";
92       cout << magneticField.swim_method << endl;
93       cout << "          --------------------------------------------------------------  " << endl;
94    }
95 }
96 db.close();
97
98 cout << endl;
99 db = QSqlDatabase();
100 db.removeDatabase("qt_sql_default_connection");
101
102 return FieldMap;
103}
104
105
106void MagneticField::create_MFM()
107{
108 string hd_msg        = gemcOpt.args["LOG_MSG"].args + " Magnetic Field: >> ";
109 double MGN_VERBOSITY = gemcOpt.args["MGN_VERBOSITY"].arg ;
110 string catch_v       = gemcOpt.args["CATCH"].args;
111
112 stringstream vars;
113 string var, format, symmetry, MapFile;
114 string var1, var2, dim;
115
116 string integration_method;       
117
118 vars << type;
119 vars >> var >> format >> symmetry >> MapFile;
120
121 mappedfield = NULL;
122
123 // %%%%%%%%%%%%%%%%%%%%%%
124 // Uniform Magnetic Field
125 // %%%%%%%%%%%%%%%%%%%%%%
126 if(var == "uniform")
127 {
128    vars.clear();
129    vars << magnitude;
130    double x,y,z;
131    vars >> var; x = get_number(var);
132    vars >> var; y = get_number(var);
133    vars >> var; z = get_number(var);
134    if(MGN_VERBOSITY>3)
135    {
136      cout << hd_msg << " <" << name << "> is a uniform magnetic field type." << endl;
137      cout << hd_msg << " <" << name << "> dimensions:" ;
138      cout << "      (" << x/gauss << ", " << y/gauss << ", " << z/gauss << ") gauss." << endl;
139    }
140    G4UniformMagField* magField = new G4UniformMagField(G4ThreeVector((x/gauss)*gauss, (y/gauss)*gauss, (z/gauss)*gauss));
141
142    G4Mag_UsualEqRhs*       iEquation    = new G4Mag_UsualEqRhs(magField);
143    G4MagIntegratorStepper* iStepper     = new G4ClassicalRK4(iEquation);
144    G4ChordFinder*          iChordFinder = new G4ChordFinder(magField, 1.0e-2*mm, iStepper);
145
146    MFM = new G4FieldManager(magField, iChordFinder);
147
148    return;
149 }
150
151 // %%%%%%%%%%%%%%%%%%%%%%%
152 // Mapped Field
153 // phi-symmetric
154 // cylindrical coordinates
155 // %%%%%%%%%%%%%%%%%%%%%%%
156 if(var == "mapped" && symmetry == "cylindrical")
157 {
158    vars.clear();
159    vars << magnitude;
160    int TNPOINTS, LNPOINTS;         
161    double tlimits[2], llimits[2];   
162    double mapOrigin[3];             
163    string units[5];                 
164
165    vars >> var;
166    TNPOINTS  = (int) get_number(var);
167    vars >> var1 >> var2 >> dim;
168    tlimits[0]   =  get_number(var1 + "*" + dim);
169    tlimits[1]   =  get_number(var2 + "*" + dim);
170    units[0].assign(dim);
171
172
173    vars >> var;
174    LNPOINTS  = (int) get_number(var);
175    vars >> var1 >> var2 >> dim;
176    llimits[0]   =  get_number(var1 + "*" + dim);
177    llimits[1]   =  get_number(var2 + "*" + dim);
178    units[1].assign(dim);
179
180    // Origin
181    vars >> var ; mapOrigin[0] =       get_number(var);
182    vars >> var ; mapOrigin[1] =       get_number(var);
183    vars >> var ; mapOrigin[2] =       get_number(var);
184    vars >> var ; units[3].assign(var);
185
186    // Field Units
187    vars >> var ; units[4].assign(var);
188
189    vars.clear();
190    vars << swim_method;
191    vars >> integration_method;
192    vars.clear();
193    vars << swim_method;
194    vars >> integration_method;
195    if(MGN_VERBOSITY>3)
196    {
197      cout << hd_msg << " <" << name << "> is a phi-symmetric mapped field in cylindrical coordinates" << endl;;
198      cout << hd_msg << " <" << name << "> has (" << TNPOINTS << ", " << LNPOINTS << ") points." << endl;;
199      cout << hd_msg << " Tranverse Boundaries:     (" << tlimits[0]/cm << ", " << tlimits[1]/cm << ") cm." << endl;
200      cout << hd_msg << " Longitudinal Boundaries:  (" << llimits[0]/cm << ", " << llimits[1]/cm << ") cm." << endl;
201      cout << hd_msg << " Map Displacement:  (" << mapOrigin[0]/cm << ", "
202                                                << mapOrigin[1]/cm << ", "
203                                                << mapOrigin[2]/cm << ") cm." << endl;
204      cout << hd_msg << " Integration Method: " << integration_method << endl;
205      cout << hd_msg << " Map Field Units: " << units[4] << endl;
206      cout << hd_msg << " Map File: " << MapFile << endl;
207    }
208    mappedfield = new MappedField(gemcOpt, TNPOINTS, LNPOINTS, tlimits, llimits, MapFile, mapOrigin, units);
209
210 }
211
212 // %%%%%%%%%%%%%%%%%%%%%%%
213 // Mapped Field
214 // phi-segmented
215 // cylindrical coordinates
216 // %%%%%%%%%%%%%%%%%%%%%%%
217 if(var == "mapped" && symmetry == "phi-segmented")
218 {
219    vars.clear();
220    vars << magnitude;
221    int TNPOINTS, PNPOINTS, LNPOINTS;             
222    double plimits[2], tlimits[2], llimits[2];   
223    double mapOrigin[3];                         
224    string units[5];                             
225
226    vars >> var ;
227    PNPOINTS  = (int) get_number(var);
228    vars >> var1 >> var2 >> dim;
229    plimits[0]   =  get_number(var1 + "*" + dim);
230    plimits[1]   =  get_number(var2 + "*" + dim);
231    units[0].assign(dim);
232
233    vars >> var ;
234    TNPOINTS  = (int) get_number(var);
235    vars >> var1 >> var2 >> dim;
236    tlimits[0]   =  get_number(var1 + "*" + dim);
237    tlimits[1]   =  get_number(var2 + "*" + dim);
238    units[1].assign(dim);
239
240    vars >> var ;
241    LNPOINTS  = (int) get_number(var);
242    vars >> var1 >> var2 >> dim;
243    llimits[0]   =  get_number(var1 + "*" + dim);
244    llimits[1]   =  get_number(var2 + "*" + dim);
245    units[2].assign(dim);
246
247    // Origin
248    vars >> var ; mapOrigin[0] =       get_number(var);
249    vars >> var ; mapOrigin[1] =       get_number(var);
250    vars >> var ; mapOrigin[2] =       get_number(var);
251    vars >> var ; units[3].assign(var);
252
253    // Field Units
254    vars >> var ; units[4].assign(var);
255
256    vars.clear();
257    vars << swim_method;
258    vars >> integration_method;
259    vars.clear();
260    vars << swim_method;
261    vars >> integration_method;
262    double segm_phi_span = 2*(plimits[1] - plimits[0]);   // factor of two: the map is assumed to cover half the segment
263    int nsegments = (int) (360.0/(segm_phi_span/deg)); 
264    if( fabs( segm_phi_span*nsegments/deg - 360 ) > 1.0e-2 )
265    {
266      cout << hd_msg << " <" << name << "> segmentation is wrong:  " << nsegments << " segments, each span "
267           << segm_phi_span/deg << ": it doesn't add to 360, but " << segm_phi_span*nsegments/deg << " Exiting." << endl;
268      exit(0);
269
270    }
271    if(MGN_VERBOSITY>3)
272    {
273      cout << hd_msg << " <" << name << "> is a phi-segmented mapped field in cylindrical coordinates" << endl;;
274      cout << hd_msg << " <" << name << "> has (" << PNPOINTS << ", " << TNPOINTS << ", " << LNPOINTS << ") points" << endl;;
275      cout << hd_msg << " Azimuthal Boundaries:     (" << plimits[0]/deg << ", " << plimits[1]/deg << ") deg" << endl;
276      cout << hd_msg << " Tranverse Boundaries:     (" << tlimits[0]/cm <<  ", " << tlimits[1]/cm  << ") cm"  << endl;
277      cout << hd_msg << " Longitudinal Boundaries:  (" << llimits[0]/cm <<  ", " << llimits[1]/cm  << ") cm"  << endl;
278      cout << hd_msg << " Phi segment span, number of segments:  " << segm_phi_span/deg <<  " deg, " << nsegments << endl;
279      cout << hd_msg << " Map Displacement:  (" << mapOrigin[0]/cm << ", "
280                                                << mapOrigin[1]/cm << ", "
281                                                << mapOrigin[2]/cm << ") cm" << endl;
282      cout << hd_msg << " Integration Method: " << integration_method << endl;
283      cout << hd_msg << " Map Field Units: " << units[4] << endl;
284      cout << hd_msg << " Map File: " << MapFile << endl;
285    }
286    mappedfield = new MappedField(gemcOpt, TNPOINTS, PNPOINTS, LNPOINTS, tlimits, plimits, llimits, MapFile, mapOrigin, units);
287    mappedfield->segm_phi_span = segm_phi_span;
288    mappedfield->nsegments     = nsegments;
289 }
290
291
292 if(integration_method == "RungeKutta")
293 {
294    // specialized equations for mapped magnetic field
295
296    G4Mag_UsualEqRhs*       iEquation    = new G4Mag_UsualEqRhs(mappedfield);
297    G4MagIntegratorStepper* iStepper     = new G4HelixSimpleRunge(iEquation);
298    G4ChordFinder*          iChordFinder = new G4ChordFinder(mappedfield, 1.0e-2*mm, iStepper);
299
300    MFM = new G4FieldManager(mappedfield, iChordFinder);
301    MFM->SetDeltaOneStep(1*mm);
302    MFM->SetDeltaIntersection(1*mm);
303 }
304}
305
306MappedField::MappedField(){;}
307MappedField::~MappedField(){;}
308
309// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
310// Constructor for phi-symmetric field in cylindrical coordinates
311// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
312MappedField::MappedField(gemc_opts Opt, int TNPOINTS, int LNPOINTS, double tlimits[2], double llimits[2],
313                         string filename, double origin[3], string units[5])
314{
315 gemcOpt = Opt;
316 string hd_msg  = gemcOpt.args["LOG_MSG"].args + " Magnetic Field Constructor: >> ";
317 MGN_VERBOSITY           = (int) gemcOpt.args["MGN_VERBOSITY"].arg ;
318
319 mapOrigin[0] = origin[0];
320 mapOrigin[1] = origin[1];
321 mapOrigin[2] = origin[2];
322
323 table_start[0] = tlimits[0];
324 table_start[1] = llimits[0];
325 cell_size[0]   = (tlimits[1] - tlimits[0])/( TNPOINTS - 1);
326 cell_size[1]   = (llimits[1] - llimits[0])/( LNPOINTS - 1);
327
328 if(MGN_VERBOSITY>3)
329 {
330    cout << hd_msg << " The transverse cell size is: "       << cell_size[0]/cm << " cm" << endl
331         << hd_msg << " and the longitudinal cell size is: " << cell_size[1]/cm << " cm" << endl;
332
333 }
334
335 ifstream IN(filename.c_str());
336 if(!IN)
337 {
338    cout << hd_msg << " File " << filename << " could not be opened. Exiting." << endl;
339    exit(0);
340 }
341 cout << hd_msg << " Reading map file: " << filename << "... ";
342
343 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
344 // Setting up storage space for tables
345 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
346 B2DCylT.resize( TNPOINTS );
347 B2DCylL.resize( TNPOINTS );
348 for(int it=0; it<TNPOINTS; it++)
349 {
350    B2DCylT[it].resize(LNPOINTS);
351    B2DCylL[it].resize(LNPOINTS);
352 }
353
354 // %%%%%%%%%%%%%
355 // Filling table
356 // %%%%%%%%%%%%%
357 double TC, LC;         // coordinates as read from the map
358 double BT, BL;         // transverse and longitudinal magnetic field as read from the map
359 unsigned int IT, IL;   // indexes of the vector map
360
361 for(int it=0; it<TNPOINTS; it++)
362 {
363    for(int il=0; il<LNPOINTS; il++)
364    {
365       IN >> TC >> LC >> BT >>  BL;
366       if(units[0] == "cm") TC *= cm;
367       else                 cout << hd_msg << " Dimension Unit " << units[0] << " not implemented yet." << endl;
368       if(units[1] == "cm") LC *= cm;
369       else                 cout << hd_msg << " Dimension Unit " << units[1] << " not implemented yet." << endl;
370
371       // checking map consistency for longitudinal coordinate
372       if( (llimits[0] + il*cell_size[1] - LC)/LC > 0.001)
373       {
374          cout << hd_msg << il << " longitudinal index wrong. Map point should be " <<  llimits[0] + il*cell_size[1]
375               << " but it's  " << LC << " instead." << endl;
376       }
377
378       IT = (unsigned int) floor( ( TC/mm - table_start[0]/mm + cell_size[0]/mm/2 ) / ( cell_size[0]/mm ) ) ;
379       IL = (unsigned int) floor( ( LC/mm - table_start[1]/mm + cell_size[1]/mm/2 ) / ( cell_size[1]/mm ) ) ;
380
381       if(units[4] == "gauss")
382       {
383          B2DCylT[IT][IL] = BT*gauss;
384          B2DCylL[IT][IL] = BL*gauss;
385       }
386       if(units[4] == "kilogauss")
387       {
388          B2DCylT[IT][IL] = BT*kilogauss;
389          B2DCylL[IT][IL] = BL*kilogauss;
390       }
391       else
392       {
393          cout << hd_msg << " Field Unit " << units[4] << " not implemented yet." << endl;
394       }
395    }
396    // checking map consistency for transverse coordinate
397    if( (tlimits[0]  + it*cell_size[0] - TC)/TC > 0.001)
398    {
399       cout << hd_msg << it << " transverse index wrong. Map point should be " <<  tlimits[0] + it*cell_size[0]
400            << " but it's  " << TC << " instead." << endl;
401    }
402 }
403
404 IN.close();
405 cout << " done!" << endl;
406
407}
408
409
410// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
411// Constructor for phi-segmented field in cylindrical coordinates. Field is in cartesian coordinates
412// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
413MappedField::MappedField(gemc_opts Opt, int rNPOINTS, int pNPOINTS, int zNPOINTS, double tlimits[2], double plimits[2], double llimits[2],
414                         string filename, double origin[3], string units[5])
415{
416 gemcOpt        = Opt;
417 string hd_msg  = gemcOpt.args["LOG_MSG"].args + " Magnetic Field Constructor: >> ";
418 MGN_VERBOSITY  = (int) gemcOpt.args["MGN_VERBOSITY"].arg ;
419
420 mapOrigin[0] = origin[0];
421 mapOrigin[1] = origin[1];
422 mapOrigin[2] = origin[2];
423
424 table_start[0] = plimits[0];
425 table_start[1] = tlimits[0];
426 table_start[2] = llimits[0];
427 cell_size[0]   = (plimits[1] - plimits[0])/( pNPOINTS - 1);
428 cell_size[1]   = (tlimits[1] - tlimits[0])/( rNPOINTS - 1);
429 cell_size[2]   = (llimits[1] - llimits[0])/( zNPOINTS - 1);
430
431 if(MGN_VERBOSITY>3)
432 {
433    cout << hd_msg << " the phi cell size is: "    << cell_size[0]/deg << " degrees" << endl
434         << hd_msg << " The radius cell size is: " << cell_size[1]/cm  << " cm" << endl
435         << hd_msg << " and the z cell size is: "  << cell_size[2]/cm  << " cm" << endl;
436 }
437
438 ifstream IN(filename.c_str());
439 if(!IN)
440 {
441    cout << hd_msg << " File " << filename << " could not be opened. Exiting." << endl;
442    exit(0);
443 }
444 cout << hd_msg << " Reading map file: " << filename << "... ";
445
446
447 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
448 // Setting up storage space for tables
449 // Field[phi][r][z]
450 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
451
452 B3DCylX.resize( pNPOINTS );
453 B3DCylY.resize( pNPOINTS );
454 B3DCylZ.resize( pNPOINTS );
455 for(int ip=0; ip<pNPOINTS; ip++)
456 {
457    B3DCylX[ip].resize( rNPOINTS );
458    B3DCylY[ip].resize( rNPOINTS );
459    B3DCylZ[ip].resize( rNPOINTS );
460    for(int ir=0; ir<rNPOINTS; ir++)
461    {
462       B3DCylX[ip][ir].resize( zNPOINTS );
463       B3DCylY[ip][ir].resize( zNPOINTS );
464       B3DCylZ[ip][ir].resize( zNPOINTS );
465    }
466 }
467
468 // %%%%%%%%%%%%%
469 // Filling table
470 // %%%%%%%%%%%%%
471 double pC, tC, lC;         // coordinates as read from the map
472 double Bx, By, Bz;         // magnetic field components as read from the map
473 unsigned int It, Ip, Il;   // indexes of the vector map
474 for(int ip=0; ip<pNPOINTS; ip++)
475 {
476    for(int it=0; it<rNPOINTS; it++)
477    {
478       for(int il=0; il<zNPOINTS; il++)
479       {
480          IN >> pC >> tC >> lC >> Bx >> By >> Bz;
481
482          if(units[0] == "deg") pC = pC*deg;
483          else                  cout << hd_msg << " Dimension Unit " << units[0] << " not implemented yet." << endl;
484          if(units[1] == "cm")  tC *= cm;
485          else                  cout << hd_msg << " Dimension Unit " << units[1] << " not implemented yet." << endl;
486          if(units[2] == "cm")  lC *= cm;
487          else                  cout << hd_msg << " Dimension Unit " << units[2] << " not implemented yet." << endl;
488
489          // checking map consistency for longitudinal coordinate
490          if( (llimits[0] + il*cell_size[2] - lC)/lC > 0.001)
491          {
492             cout << hd_msg << il << " longitudinal index wrong. Map point should be " <<  llimits[0] + il*cell_size[2]
493                  << " but it's  " << lC << " instead." << endl;
494          }
495
496          Ip = (unsigned int) floor( ( pC/deg - table_start[0]/deg + cell_size[0]/deg/2 )  / ( cell_size[0]/deg ) ) ;
497          It = (unsigned int) floor( ( tC/mm  - table_start[1]/mm  + cell_size[1]/mm/2  )  / ( cell_size[1]/mm  ) ) ;
498          Il = (unsigned int) floor( ( lC/mm  - table_start[2]/mm  + cell_size[2]/mm/2  )  / ( cell_size[2]/mm  ) ) ;
499
500          if(units[4] == "gauss")
501          {
502             B3DCylX[Ip][It][Il] = Bx*gauss;
503             B3DCylY[Ip][It][Il] = By*gauss;
504             B3DCylZ[Ip][It][Il] = Bz*gauss;
505          }
506          if(units[4] == "kilogauss")
507          {
508             B3DCylX[Ip][It][Il] = Bx*kilogauss;
509             B3DCylY[Ip][It][Il] = By*kilogauss;
510             B3DCylZ[Ip][It][Il] = Bz*kilogauss;
511          }
512          else
513          {
514             cout << hd_msg << " Field Unit " << units[4] << " not implemented yet." << endl;
515          }
516
517          if(MGN_VERBOSITY>10)
518          {
519             cout << " phi=" << pC/deg << ", ip=" << Ip << "   "
520                  << " r="   << tC/cm  << ", it=" << It << "   "
521                  << " z="   << lC/cm  << ", iz=" << Il << "   "
522                  << " Bx=" << B3DCylX[Ip][It][Il]/kilogauss  << " "
523                  << " By=" << B3DCylY[Ip][It][Il]/kilogauss  << " "
524                  << " Bz=" << B3DCylZ[Ip][It][Il]/kilogauss  << " kilogauss.   Map Values: "
525                  << " rBx=" << Bx  << " "
526                  << " rBy=" << By  << " "
527                  << " rBz=" << Bz  << endl;
528
529          }
530
531       }
532
533       // checking map consistency for transverse coordinate
534       if( (tlimits[0]  + it*cell_size[1] - tC)/tC > 0.001)
535       {
536          cout << hd_msg << it << " transverse index wrong. Map point should be " <<  tlimits[0] + it*cell_size[0]
537               << " but it's  " << tC << " instead." << endl;
538       }
539
540    }
541
542    // checking map consistency for azimuthal coordinate
543    if( (plimits[0] + ip*cell_size[0] - pC)/pC > 0.001)
544    {
545       cout << hd_msg << ip << " azimuthal index wrong. Map point should be " <<  plimits[0] + ip*cell_size[1]
546            << " but it's  " << pC << " instead." << endl;
547    }
548
549 }
550
551 IN.close();
552 cout << " done!" << endl;
553
554}
555
556
557
558
559// %%%%%%%%%%%%%
560// GetFieldValue
561// %%%%%%%%%%%%%
562void MappedField::GetFieldValue(const double point[3], double *Bfield) const
563{
564
565 double Point[3]; // global coordinates, in mm, after shifting the origin
566 vector<double> Field[3];
567
568 Point[0] = point[0] - mapOrigin[0]/mm;
569 Point[1] = point[1] - mapOrigin[1]/mm;
570 Point[2] = point[2] - mapOrigin[2]/mm;
571
572 Bfield[0] = Bfield[1] = Bfield[2] = 0*gauss;
573
574 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
575 // phi symmetric field in cylindrical coordinates
576 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
577 if(B2DCylT.size() && B2DCylL.size())
578 {
579    double psfield[3] = {0,0,0};
580    unsigned int IT, IL;
581
582    double LC;            // longitudinal coordinate of the track in the global coordinate system
583    double TC, phi;       // transverse and phy polar 2D coordinates: the map is phi-symmetric. phi is angle in respect to x
584    TC              = sqrt(Point[0]*Point[0] + Point[1]*Point[1]);
585    if( TC!=0 ) phi = acos(Point[0]/TC);
586    else        phi = 0;
587    LC    = Point[2];
588
589    IT = (unsigned int) floor( ( TC - table_start[0]/mm ) / (cell_size[0]/mm) ) ;
590    IL = (unsigned int) floor( ( LC - table_start[1]/mm ) / (cell_size[1]/mm) ) ;
591
592    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++;
593    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++;
594
595
596    // vector sizes are checked on both T and L components
597    // (even if only one is enough)
598    if(IT < B2DCylT.size() && IT < B2DCylL.size())
599       if(IL < B2DCylT[IT].size() && IL < B2DCylL[IT].size()  )
600       {
601          psfield[0] = B2DCylT[IT][IL] * cos(phi);
602          psfield[1] = B2DCylT[IT][IL] * sin(phi);
603          psfield[2] = B2DCylL[IT][IL];
604          if(MGN_VERBOSITY>5)
605           {
606             cout << hd_msg << " Phi-Simmetric Field: Cart. and Cyl. coordinates (cm), table indexes, magnetic field values (gauss):" << endl;
607             cout << " x="    << point[0]/cm << "   " ;
608             cout << "y="     << point[1]/cm << "   ";
609             cout << "z="     << point[2]/cm << "   ";
610             cout << "r="     << TC/cm << "   ";
611             cout << "z="     << LC/cm << "  ";
612             cout << "phi="   << phi*180/3.141592 << "   ";
613             cout << "IT="    << IT << "   ";
614             cout << "IL="    << IL << "  ";
615             cout << "Bx="    << psfield[0]/gauss << "  ";
616             cout << "By="    << psfield[1]/gauss << "  ";
617             cout << "Bz="    << psfield[2]/gauss << endl;
618          }
619       }
620    for(int i=0; i<3; i++) Field[i].push_back(psfield[i]);
621 }
622
623
624 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
625 // phi-segmented field in cylindrical coordinates
626 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
627 if(B3DCylX.size() && B3DCylY.size() && B3DCylZ.size())
628 {
629    double mfield[3]      = {0,0,0};
630    double rotpsfield[3]  = {0,0,0};
631
632    double tC, pC, lC;         
633    double pLC;               
634    double Bx, By, Bz;         
635    unsigned int Ip, It, Il;   
636
637    pC = atan2( Point[1], Point[0] )*rad;
638    if( pC < 0 ) pC += 360*deg;
639
640    tC = sqrt(Point[0]*Point[0] + Point[1]*Point[1]); 
641    lC = Point[2];                                     
642
643    // Rotating the point to within the map limit
644    int segment = 0;
645    pLC = pC;
646    while (pLC/deg > 30)
647    {
648       pLC -= 60*deg;
649       segment++;
650    }
651
652    double dphi = pC - pLC;
653    int    sign = (pLC >= 0 ? 1 : -1);
654    double apLC = fabs(pLC);
655
656    Ip = (unsigned int) floor( ( apLC/deg - table_start[0]/deg ) / (cell_size[0]/deg ) ) ;
657    It = (unsigned int) floor( ( tC/mm    - table_start[1]/mm  ) / (cell_size[1]/mm  ) ) ;
658    Il = (unsigned int) floor( ( lC/mm    - table_start[2]/mm  ) / (cell_size[2]/mm  ) ) ;
659
660    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++;
661    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++;
662    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++;
663
664    // Getting Field at the rotated point
665    // vector sizes are checked on all components
666    // (even if only one is enough)
667    if(      Ip < B3DCylX.size()         && Ip < B3DCylY.size()         && Ip < B3DCylZ.size())
668       if(   It < B3DCylX[Ip].size()     && It < B3DCylY[Ip].size()     && It < B3DCylZ[Ip].size())
669          if(Il < B3DCylX[Ip][It].size() && Il < B3DCylY[Ip][It].size() && Il < B3DCylZ[Ip][It].size())
670          {
671             // Field at local point
672             mfield[0] = B3DCylX[Ip][It][Il];
673             mfield[1] = B3DCylY[Ip][It][Il];
674             mfield[2] = B3DCylZ[Ip][It][Il];
675
676
677             // Rotating the field back to original point
678             rotpsfield[0] =  sign*mfield[0] * cos(dphi/rad) - mfield[1] * sin(dphi/rad);
679             rotpsfield[1] = +sign*mfield[0] * sin(dphi/rad) + mfield[1] * cos(dphi/rad);
680             rotpsfield[2] =  sign*mfield[2];
681
682             if(MGN_VERBOSITY>6)
683              {
684                cout << hd_msg << " Phi-Segmented Field: Cart. and Cyl. coord. (cm), indexes, local and rotated field values (gauss):" << endl;
685                cout << " x="       << point[0]/cm << "   " ;
686                cout << "y="        << point[1]/cm << "   ";
687                cout << "z="        << point[2]/cm << "   ";
688                cout << "phi="      << pC/deg << "   ";
689                cout << "r="        << tC/cm << "   ";
690                cout << "z="        << lC/cm << "  ";
691                cout << "dphi="     << dphi/deg << "   ";
692                cout << "dphir="     << dphi/rad << "   ";
693                cout << "lphi="     << (sign == 1 ? "+ " : "- ") << pLC/deg << "   ";
694                cout << "segment="  << segment << "   ";
695                cout << "Ip="       << Ip << "  ";
696                cout << "It="       << It << "   ";
697                cout << "Il="       << Il << "  ";
698                cout << "lBx="      << mfield[0]/gauss << "  ";
699                cout << "lBy="      << mfield[1]/gauss << "  ";
700                cout << "lBz="      << mfield[2]/gauss << " " ;
701                cout << "rBx="      << rotpsfield[0]/gauss << "  ";
702                cout << "rBy="      << rotpsfield[1]/gauss << "  ";
703                cout << "rBz="      << rotpsfield[2]/gauss << endl;
704             }
705          }
706    for(int i=0; i<3; i++) Field[i].push_back(-rotpsfield[i]);
707 }
708
709
710 // %%%%%%%%%%%%%%%%%%
711 // Summing the Fields
712 // %%%%%%%%%%%%%%%%%%
713 for(int i=0; i<Field[0].size(); i++)
714     for(int j=0; j<3; j++)
715       Bfield[j] +=  Field[j][i];
716
717 if(MGN_VERBOSITY>5)
718 {
719    cout << hd_msg << " Total Field: coordinates (cm), magnetic field values (gauss):" << endl ;
720    cout << " x="    << point[0]/cm << "   " ;
721    cout << "y="     << point[1]/cm << "   ";
722    cout << "z="     << point[2]/cm << "   ";
723    cout << "Bx="    << Bfield[0]/gauss << "  ";
724    cout << "By="    << Bfield[1]/gauss << "  ";
725    cout << "Bz="    << Bfield[2]/gauss << endl << endl;
726 }
727
728}
729
730
731
732
733
734
735
736
737
738
739
740
741
Note: See TracBrowser for help on using the repository browser.