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

Last change on this file since 1201 was 807, checked in by garnier, 17 years ago

update
File size: 26.0 KB
RevLine 
[807]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.