4 #include <TDirectory.h>
9 #include <Geant4/G4SystemOfUnits.hh>
28 std::cout <<
" ------------- PHField2D::PHField2D() ------------------" << std::endl;
31 TFile *rootinput = TFile::Open(filename.c_str());
34 std::cout <<
" could not open " << filename <<
" exiting now" << std::endl;
40 std::cout <<
" Field grid file: " << filename << std::endl;
44 Float_t ROOT_Z, ROOT_R;
45 Float_t ROOT_BZ, ROOT_BR;
47 TNtuple *field_map = (TNtuple *) gDirectory->Get(
"fieldmap");
54 field_map = (TNtuple *) gDirectory->Get(
"map");
57 std::cout <<
"PHField2D: could not locate ntuple of name map or fieldmap, exiting now" << std::endl;
62 field_map->SetBranchAddress(
"z", &ROOT_Z);
63 field_map->SetBranchAddress(
"r", &ROOT_R);
64 field_map->SetBranchAddress(
"bz", &ROOT_BZ);
65 field_map->SetBranchAddress(
"br", &ROOT_BR);
69 nz = field_map->GetEntries(
"z>-1e6");
70 nr = field_map->GetEntries(
"r>-1e6");
71 static const int NENTRIES = field_map->GetEntries();
76 std::cout <<
" The field grid contained " << NENTRIES <<
" entries" << std::endl;
80 std::cout <<
"\n NENTRIES should be the same as the following values:"
81 <<
"\n [ Number of values r,z: "
82 << nr <<
" " << nz <<
" ]! " << std::endl;
87 std::cout <<
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"
88 <<
"\n The file you entered is not a \"table\" of values"
89 <<
"\n Something very likely went oh so wrong"
90 <<
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
94 std::set<float> z_set, r_set;
103 std::cout <<
" --> Sorting Entries..." << std::endl;
105 std::map<trio, trio> sorted_map;
106 for (
int i = 0;
i < field_map->GetEntries();
i++)
108 field_map->GetEntry(
i);
109 trio coord_key(ROOT_Z, ROOT_R);
110 trio field_val(ROOT_BZ, ROOT_BR);
111 sorted_map[coord_key] = field_val;
113 z_set.insert(ROOT_Z *
cm);
114 r_set.insert(ROOT_R * cm);
120 std::map<trio, trio>::iterator
it = sorted_map.begin();
122 float last_z = std::get<0>(it->first);
123 for (it = sorted_map.begin(); it != sorted_map.end(); ++
it)
125 if (std::get<0>(it->first) != last_z)
127 last_z = std::get<0>(it->first);
135 std::cout <<
" --> Putting entries into containers... " << std::endl;
139 minz_ = *(z_set.begin());
140 std::set<float>::iterator ziter = z_set.end();
150 std::copy(z_set.begin(), z_set.end(),
z_map_.begin());
151 std::copy(r_set.begin(), r_set.end(),
r_map_.begin());
154 BFieldR_.resize(nz, std::vector<float>(nr, 0));
155 BFieldZ_.resize(nz, std::vector<float>(nr, 0));
158 unsigned int ir = 0, iz = 0;
159 std::map<trio, trio>::iterator iter = sorted_map.begin();
160 for (; iter != sorted_map.end(); ++iter)
163 float z = std::get<0>(iter->first) *
cm;
164 float r = std::get<1>(iter->first) *
cm;
189 if (iz > 0 && z <
z_map_[iz - 1])
191 std::cout <<
"!!!!!!!!! Your map isn't ordered.... z: " << z <<
" zprev: " <<
z_map_[iz - 1] << std::endl;
200 if (std::fabs(z) < 10 && ir < 10 &&
Verbosity() > 3)
208 <<
BFieldZ_[iz][ir] <<
")" << std::endl;
217 std::cout <<
" Mag field z boundaries (min,max): (" <<
minz_ /
cm <<
", " <<
maxz_ /
cm <<
") cm" << std::endl;
221 std::cout <<
" Mag field r max boundary: " <<
r_map_.back() /
cm <<
" cm" << std::endl;
226 std::cout <<
" -----------------------------------------------------------" << std::endl;
234 std::cout <<
"\nPHField2D::GetFieldValue" << std::endl;
239 double r = sqrt(x * x + y * y);
251 double cylpoint[4] = {
z,
r,
phi, 0};
257 Bfield[0] = cos(phi) * BFieldCyl[1] - sin(phi) * BFieldCyl[2];
260 Bfield[1] = sin(phi) * BFieldCyl[1] + cos(phi) * BFieldCyl[2];
263 Bfield[2] = BFieldCyl[0];
272 std::cout <<
"!!!!!!!!!! Field point not in defined region (outside of z bounds)" << std::endl;
278 std::cout <<
"END PHField2D::GetFieldValue\n"
279 <<
" ---> {Bx, By, Bz} : "
280 <<
"< " << Bfield[0] <<
", " << Bfield[1] <<
", " << Bfield[2] <<
" >" << std::endl;
288 float z = CylPoint[0];
289 float r = CylPoint[1];
297 std::cout <<
"GetFieldCyl@ <z,r>: {" << z <<
"," << r <<
"}" << std::endl;
304 std::cout <<
"!!!! Point not in defined region (radius too large in specific z-plane)" << std::endl;
316 if (!((r >
r_map_[r_index0]) && (r <
r_map_[r_index1])))
319 std::vector<float>::const_iterator riter = upper_bound(
r_map_.begin(),
r_map_.end(),
r);
321 if (r_index0 >=
r_map_.size())
325 std::cout <<
"!!!! Point not in defined region (radius too large in specific z-plane)" << std::endl;
330 r_index1 = r_index0 + 1;
331 if (r_index1 >=
r_map_.size())
335 std::cout <<
"!!!! Point not in defined region (radius too large in specific z-plane)" << std::endl;
348 if (!((z >
z_map_[z_index0]) && (z <
z_map_[z_index1])))
351 std::vector<float>::const_iterator ziter = upper_bound(
z_map_.begin(),
z_map_.end(),
z);
353 z_index1 = z_index0 + 1;
354 if (z_index1 >=
z_map_.size())
358 std::cout <<
"!!!! Point not in defined region (z too large in specific r-plane)" << std::endl;
368 double Br000 =
BFieldR_[z_index0][r_index0];
369 double Br010 =
BFieldR_[z_index0][r_index1];
370 double Br100 =
BFieldR_[z_index1][r_index0];
371 double Br110 =
BFieldR_[z_index1][r_index1];
373 double Bz000 =
BFieldZ_[z_index0][r_index0];
374 double Bz100 =
BFieldZ_[z_index1][r_index0];
375 double Bz010 =
BFieldZ_[z_index0][r_index1];
376 double Bz110 =
BFieldZ_[z_index1][r_index1];
378 double zweight = z -
z_map_[z_index0];
379 double zspacing = z_map_[z_index1] - z_map_[z_index0];
382 double rweight = r -
r_map_[r_index0];
383 double rspacing = r_map_[r_index1] - r_map_[r_index0];
388 (1 - zweight) * ((1 - rweight) * Bz000 +
390 zweight * ((1 - rweight) * Bz100 +
395 (1 - zweight) * ((1 - rweight) * Br000 +
397 zweight * ((1 - rweight) * Br100 +
405 std::cout <<
"End GFCyl Call: <bz,br,bphi> : {"
406 << BfieldCyl[0] /
gauss <<
"," << BfieldCyl[1] /
gauss <<
"," << BfieldCyl[2] /
gauss <<
"}"
416 std::cout <<
" Key: <"
417 << std::get<0>(it->first) /
cm <<
","
418 << std::get<1>(it->first) /
cm <<
">"
422 << std::get<1>(it->second) /
magfield_unit <<
">" << std::endl;