3 #include <TDirectory.h>
7 #include <Geant4/G4SystemOfUnits.hh>
9 #include <boost/tuple/tuple_comparison.hpp>
25 cout <<
"\n================ Begin Construct Mag Field =====================" << endl;
26 cout <<
"\n-----------------------------------------------------------"
27 <<
"\n Magnetic field Module - Verbosity:" <<
Verbosity()
28 <<
"\n-----------------------------------------------------------";
31 TFile *rootinput = TFile::Open(filename.c_str());
34 cout <<
"\n could not open " << filename <<
" exiting now" << endl;
38 "Reading the field grid from "
39 << filename <<
" ... " << endl;
43 TNtuple *field_map = (TNtuple *) gDirectory->Get(
"map");
44 Float_t ROOT_Z, ROOT_R, ROOT_PHI;
45 Float_t ROOT_BZ, ROOT_BR, ROOT_BPHI;
46 field_map->SetBranchAddress(
"z", &ROOT_Z);
47 field_map->SetBranchAddress(
"r", &ROOT_R);
48 field_map->SetBranchAddress(
"phi", &ROOT_PHI);
49 field_map->SetBranchAddress(
"bz", &ROOT_BZ);
50 field_map->SetBranchAddress(
"br", &ROOT_BR);
51 field_map->SetBranchAddress(
"bphi", &ROOT_BPHI);
55 nz = field_map->GetEntries(
"z>-1e6");
56 nr = field_map->GetEntries(
"r>-1e6");
57 nphi = field_map->GetEntries(
"phi>-1e6");
58 static const int NENTRIES = field_map->GetEntries();
61 cout <<
" ---> The field grid contained " << NENTRIES <<
" entries" << endl;
64 cout <<
"\n NENTRIES should be the same as the following values:"
65 <<
"\n [ Number of values r,phi,z: "
66 << nr <<
" " << nphi <<
" " << nz <<
" ]! " << endl;
69 if (nz != nr || nz != nphi || nr != nphi)
71 cout <<
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"
72 <<
"\n The file you entered is not a \"table\" of values"
73 <<
"\n Something very likely went oh so wrong"
74 <<
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
78 std::set<float> z_set, r_set, phi_set;
87 cout <<
" --> Sorting Entries..." << endl;
89 std::map<trio, trio> sorted_map;
90 for (
int i = 0; i < field_map->GetEntries(); i++)
92 field_map->GetEntry(i);
93 trio coord_key(ROOT_Z, ROOT_R, ROOT_PHI);
94 trio field_val(ROOT_BZ, ROOT_BR, ROOT_BPHI);
95 sorted_map[coord_key] = field_val;
97 z_set.insert(ROOT_Z *
cm);
98 r_set.insert(ROOT_R * cm);
99 phi_set.insert(ROOT_PHI *
deg);
105 map<trio, trio>::iterator
it = sorted_map.begin();
107 float last_z = it->first.get<0>();
108 for (it = sorted_map.begin(); it != sorted_map.end(); ++
it)
110 if (it->first.get<0>() != last_z)
112 last_z = it->first.get<0>();
120 cout <<
" --> Putting entries into containers... " << endl;
124 minz_ = *(z_set.begin());
125 std::set<float>::iterator ziter = z_set.end();
132 nphi = phi_set.size();
142 BFieldR_.resize(nz, vector<vector<float> >(nr, vector<float>(nphi, 0)));
143 BFieldPHI_.resize(nz, vector<vector<float> >(nr, vector<float>(nphi, 0)));
144 BFieldZ_.resize(nz, vector<vector<float> >(nr, vector<float>(nphi, 0)));
147 unsigned int ir = 0, iphi = 0, iz = 0;
148 map<trio, trio>::iterator iter = sorted_map.begin();
149 for (; iter != sorted_map.end(); ++iter)
152 float z = iter->first.get<0>() *
cm;
153 float r = iter->first.get<1>() *
cm;
154 float phi = iter->first.get<2>() *
deg;
155 float Bz = iter->second.get<0>() *
gauss;
156 float Br = iter->second.get<1>() *
gauss;
157 float Bphi = iter->second.get<2>() *
gauss;
180 if (iz > 0 && z <
z_map_[iz - 1])
182 cout <<
"!!!!!!!!! Your map isn't ordered.... z: " << z <<
" zprev: " <<
z_map_[iz - 1] << endl;
192 if (fabs(z) < 10 && ir < 10 &&
Verbosity() > 3)
202 <<
BFieldZ_[iz][ir][iphi] <<
")" << endl;
209 cout <<
"\n ---> ... read file successfully "
210 <<
"\n ---> Z Boundaries ~ zlow, zhigh: "
213 cout <<
"\n================= End Construct Mag Field ======================\n"
220 cout <<
"\nPHField3DCylindrical::GetFieldValue" << endl;
224 double r = sqrt(x * x + y * y);
228 phi = atan2(y, 0.00000000001);
234 if (phi < 0) phi += 2 *
M_PI;
240 double cylpoint[4] = {
z,
r,
phi, 0};
246 Bfield[0] = cos(phi) * BFieldCyl[1] - sin(phi) * BFieldCyl[2];
249 Bfield[1] = sin(phi) * BFieldCyl[1] + cos(phi) * BFieldCyl[2];
252 Bfield[2] = BFieldCyl[0];
261 cout <<
"!!!!!!!!!! Field point not in defined region (outside of z bounds)" << endl;
266 cout <<
"END PHField3DCylindrical::GetFieldValue\n"
267 <<
" ---> {Bx, By, Bz} : "
268 <<
"< " << Bfield[0] <<
", " << Bfield[1] <<
", " << Bfield[2] <<
" >" << endl;
276 float z = CylPoint[0];
277 float r = CylPoint[1];
278 float phi = CylPoint[2];
285 cout <<
"GetFieldCyl@ <z,r,phi>: {" << z <<
"," << r <<
"," << phi <<
"}" << endl;
290 cout <<
"!!!! Point not in defined region (|z| too large)" << endl;
297 cout <<
"!!!! Point not in defined region (radius too small in specific z-plane). Use min radius" << endl;
303 cout <<
"!!!! Point not in defined region (radius too large in specific z-plane)" << endl;
307 vector<float>::const_iterator ziter = upper_bound(
z_map_.begin(),
z_map_.end(),
z);
308 int z_index0 = distance(
z_map_.begin(), ziter) - 1;
309 int z_index1 = z_index0 + 1;
311 assert(z_index0 >= 0);
312 assert(z_index1 >= 0);
313 assert(z_index0 < (
int)
z_map_.size());
314 assert(z_index1 < (
int)
z_map_.size());
316 vector<float>::const_iterator riter = upper_bound(
r_map_.begin(),
r_map_.end(),
r);
317 int r_index0 = distance(
r_map_.begin(), riter) - 1;
318 if (r_index0 >= (
int)
r_map_.size())
321 cout <<
"!!!! Point not in defined region (radius too large in specific z-plane)" << endl;
325 int r_index1 = r_index0 + 1;
326 if (r_index1 >= (
int)
r_map_.size())
329 cout <<
"!!!! Point not in defined region (radius too large in specific z-plane)" << endl;
333 assert(r_index0 >= 0);
334 assert(r_index1 >= 0);
336 vector<float>::const_iterator phiiter = upper_bound(
phi_map_.begin(),
phi_map_.end(),
phi);
337 int phi_index0 = distance(
phi_map_.begin(), phiiter) - 1;
338 int phi_index1 = phi_index0 + 1;
339 if (phi_index1 >= (
int)
phi_map_.size())
342 assert(phi_index0 >= 0);
343 assert(phi_index0 < (
int)
phi_map_.size());
344 assert(phi_index1 >= 0);
346 double Br000 =
BFieldR_[z_index0][r_index0][phi_index0];
347 double Br001 =
BFieldR_[z_index0][r_index0][phi_index1];
348 double Br010 =
BFieldR_[z_index0][r_index1][phi_index0];
349 double Br011 =
BFieldR_[z_index0][r_index1][phi_index1];
350 double Br100 =
BFieldR_[z_index1][r_index0][phi_index0];
351 double Br101 =
BFieldR_[z_index1][r_index0][phi_index1];
352 double Br110 =
BFieldR_[z_index1][r_index1][phi_index0];
353 double Br111 =
BFieldR_[z_index1][r_index1][phi_index1];
355 double Bphi000 =
BFieldPHI_[z_index0][r_index0][phi_index0];
356 double Bphi001 =
BFieldPHI_[z_index0][r_index0][phi_index1];
357 double Bphi010 =
BFieldPHI_[z_index0][r_index1][phi_index0];
358 double Bphi011 =
BFieldPHI_[z_index0][r_index1][phi_index1];
359 double Bphi100 =
BFieldPHI_[z_index1][r_index0][phi_index0];
360 double Bphi101 =
BFieldPHI_[z_index1][r_index0][phi_index1];
361 double Bphi110 =
BFieldPHI_[z_index1][r_index1][phi_index0];
362 double Bphi111 =
BFieldPHI_[z_index1][r_index1][phi_index1];
364 double Bz000 =
BFieldZ_[z_index0][r_index0][phi_index0];
365 double Bz001 =
BFieldZ_[z_index0][r_index0][phi_index1];
366 double Bz100 =
BFieldZ_[z_index1][r_index0][phi_index0];
367 double Bz101 =
BFieldZ_[z_index1][r_index0][phi_index1];
368 double Bz010 =
BFieldZ_[z_index0][r_index1][phi_index0];
369 double Bz110 =
BFieldZ_[z_index1][r_index1][phi_index0];
370 double Bz011 =
BFieldZ_[z_index0][r_index1][phi_index1];
371 double Bz111 =
BFieldZ_[z_index1][r_index1][phi_index1];
373 double zweight = z -
z_map_[z_index0];
374 double zspacing = z_map_[z_index1] - z_map_[z_index0];
377 double rweight = r -
r_map_[r_index0];
378 double rspacing = r_map_[r_index1] - r_map_[r_index0];
381 double phiweight = phi -
phi_map_[phi_index0];
382 double phispacing = phi_map_[phi_index1] - phi_map_[phi_index0];
384 phispacing += 2 *
M_PI;
385 phiweight /= phispacing;
389 (1 - zweight) * ((1 - rweight) * ((1 - phiweight) * Bz000 + phiweight * Bz001) +
390 rweight * ((1 - phiweight) * Bz010 + phiweight * Bz011)) +
391 zweight * ((1 - rweight) * ((1 - phiweight) * Bz100 + phiweight * Bz101) +
392 rweight * ((1 - phiweight) * Bz110 + phiweight * Bz111));
396 (1 - zweight) * ((1 - rweight) * ((1 - phiweight) * Br000 + phiweight * Br001) +
397 rweight * ((1 - phiweight) * Br010 + phiweight * Br011)) +
398 zweight * ((1 - rweight) * ((1 - phiweight) * Br100 + phiweight * Br101) +
399 rweight * ((1 - phiweight) * Br110 + phiweight * Br111));
403 (1 - zweight) * ((1 - rweight) * ((1 - phiweight) * Bphi000 + phiweight * Bphi001) +
404 rweight * ((1 - phiweight) * Bphi010 + phiweight * Bphi011)) +
405 zweight * ((1 - rweight) * ((1 - phiweight) * Bphi100 + phiweight * Bphi101) +
406 rweight * ((1 - phiweight) * Bphi110 + phiweight * Bphi111));
421 cout <<
"End GFCyl Call: <bz,br,bphi> : {"
422 << BfieldCyl[0] /
gauss <<
"," << BfieldCyl[1] /
gauss <<
"," << BfieldCyl[2] /
gauss <<
"}"
442 unsigned middle = start + ((end -
start) / 2);
444 if (vec[middle] == key)
449 else if (vec[middle] > key)
451 return bin_search(vec, start, middle - 1, key, index);
453 return bin_search(vec, middle + 1, end, key, index);
460 << it->first.get<0>() *
cm <<
","
461 << it->first.get<1>() *
cm <<
","
462 << it->first.get<2>() *
deg <<
">"
465 << it->second.get<0>() *
gauss <<
","
466 << it->second.get<1>() *
gauss <<
","
467 << it->second.get<2>() *
gauss <<
">\n";