| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // G4AdjointCrossSurfChecker class implementat
27 //
28 // Author: L. Desorgher, SpaceIT GmbH
29 // Contract: ESA contract 21435/08/NL/AT
30 // Customer: ESA/ESTEC
31 // -------------------------------------------
32
33 #include "G4AdjointCrossSurfChecker.hh"
34
35 #include "G4AffineTransform.hh"
36 #include "G4PhysicalConstants.hh"
37 #include "G4PhysicalVolumeStore.hh"
38 #include "G4Step.hh"
39 #include "G4StepPoint.hh"
40 #include "G4SystemOfUnits.hh"
41 #include "G4VSolid.hh"
42
43 //////////////////////////////////////////////
44 //
45 G4ThreadLocal G4AdjointCrossSurfChecker* G4Adj
46
47 //////////////////////////////////////////////
48 //
49 G4AdjointCrossSurfChecker::~G4AdjointCrossSurf
50
51 //////////////////////////////////////////////
52 //
53 G4AdjointCrossSurfChecker* G4AdjointCrossSurfC
54 {
55 if (instance == nullptr) instance = new G4Ad
56 return instance;
57 }
58
59 //////////////////////////////////////////////
60 //
61 G4bool G4AdjointCrossSurfChecker::CrossingASph
62 G4ThreeVector sphere_center, G4ThreeVector&
63 {
64 G4ThreeVector pos1 = aStep->GetPreStepPoint(
65 G4ThreeVector pos2 = aStep->GetPostStepPoint
66 G4double r1 = pos1.mag();
67 G4double r2 = pos2.mag();
68 G4bool did_cross = false;
69
70 if (r1 <= sphere_radius && r2 > sphere_radiu
71 did_cross = true;
72 GoingIn = false;
73 }
74 else if (r2 <= sphere_radius && r1 > sphere_
75 did_cross = true;
76 GoingIn = true;
77 }
78
79 if (did_cross) {
80 G4ThreeVector dr = pos2 - pos1;
81 G4double r12 = r1 * r1;
82 G4double rdr = dr.mag();
83 G4double a, b, c, d;
84 a = rdr * rdr;
85 b = 2. * pos1.dot(dr);
86 c = r12 - sphere_radius * sphere_radius;
87 d = std::sqrt(b * b - 4. * a * c);
88 G4double l = (-b + d) / 2. / a;
89 if (l > 1.) l = (-b - d) / 2. / a;
90 crossing_pos = pos1 + l * dr;
91 cos_th = std::abs(dr.cosTheta(crossing_pos
92 }
93 return did_cross;
94 }
95
96 //////////////////////////////////////////////
97 //
98 G4bool G4AdjointCrossSurfChecker::GoingInOrOut
99 const G4String& volume_name, G4double&, G4bo
100 {
101 G4bool step_at_boundary = (aStep->GetPostSte
102 G4bool did_cross = false;
103 if (step_at_boundary) {
104 const G4VTouchable* postStepTouchable = aS
105 const G4VTouchable* preStepTouchable = aSt
106 if ((preStepTouchable != nullptr) && (post
107 (postStepTouchable->GetVolume() != nul
108 {
109 G4String post_vol_name = postStepTouchab
110 G4String pre_vol_name = preStepTouchable
111
112 if (post_vol_name == volume_name) {
113 GoingIn = true;
114 did_cross = true;
115 }
116 else if (pre_vol_name == volume_name) {
117 GoingIn = false;
118 did_cross = true;
119 }
120 }
121 }
122 return did_cross; // still need to compute
123 }
124
125 //////////////////////////////////////////////
126 //
127 G4bool G4AdjointCrossSurfChecker::GoingInOrOut
128 const G4String& volume_name, const G4String&
129 G4bool& GoingIn) // from external surf.
130 {
131 G4bool step_at_boundary = (aStep->GetPostSte
132 G4bool did_cross = false;
133 if (step_at_boundary) {
134 const G4VTouchable* postStepTouchable = aS
135 const G4VTouchable* preStepTouchable = aSt
136 const G4VPhysicalVolume* postVol =
137 (postStepTouchable != nullptr) ? postSte
138 const G4VPhysicalVolume* preVol =
139 (preStepTouchable != nullptr) ? preStepT
140 if (preStepTouchable != nullptr && postSte
141 preVol != nullptr)
142 {
143 G4String post_vol_name = postVol->GetNam
144 G4String post_log_vol_name = postVol->Ge
145 G4String pre_vol_name = preVol->GetName(
146 G4String pre_log_vol_name = preVol->GetL
147 if (post_vol_name == volume_name && pre_
148 GoingIn = true;
149 did_cross = true;
150 }
151 else if (pre_vol_name == volume_name &&
152 GoingIn = false;
153 did_cross = true;
154 }
155 }
156 }
157 return did_cross; // still need to compute
158 }
159
160 //////////////////////////////////////////////
161 //
162 G4bool G4AdjointCrossSurfChecker::CrossingAGiv
163 const G4String& surface_name, G4ThreeVector&
164 G4bool& GoingIn)
165 {
166 G4int ind = FindRegisteredSurface(surface_na
167 G4bool did_cross = false;
168 if (ind >= 0) {
169 did_cross = CrossingAGivenRegisteredSurfac
170 }
171 return did_cross;
172 }
173
174 //////////////////////////////////////////////
175 //
176 G4bool G4AdjointCrossSurfChecker::CrossingAGiv
177 G4ThreeVector& crossing_pos, G4double& cos_t
178 {
179 G4String surf_type = ListOfSurfaceType[ind];
180 G4double radius = ListOfSphereRadius[ind];
181 G4ThreeVector center = ListOfSphereCenter[in
182 G4String vol1 = ListOfVol1Name[ind];
183 G4String vol2 = ListOfVol2Name[ind];
184
185 G4bool did_cross = false;
186 if (surf_type == "Sphere") {
187 did_cross = CrossingASphere(aStep, radius,
188 }
189 else if (surf_type == "ExternalSurfaceOfAVol
190 did_cross = GoingInOrOutOfaVolumeByExtSurf
191 crossing_pos = aStep->GetPostStepPoint()->
192 }
193 else if (surf_type == "BoundaryBetweenTwoVol
194 did_cross = CrossingAnInterfaceBetweenTwoV
195 aStep, vol1, vol2, crossing_pos, cos_to_
196 }
197 return did_cross;
198 }
199
200 //////////////////////////////////////////////
201 //
202 G4bool G4AdjointCrossSurfChecker::CrossingOneO
203 G4String& surface_name, G4ThreeVector& cross
204 {
205 for (std::size_t i = 0; i < ListOfSurfaceNam
206 if (CrossingAGivenRegisteredSurface(aStep,
207 surface_name = ListOfSurfaceName[i];
208 return true;
209 }
210 }
211 return false;
212 }
213
214 //////////////////////////////////////////////
215 //
216 G4bool G4AdjointCrossSurfChecker::CrossingAnIn
217 const G4String& vol1_name, const G4String& v
218 {
219 G4bool step_at_boundary = (aStep->GetPostSte
220 G4bool did_cross = false;
221 if (step_at_boundary) {
222 const G4VTouchable* postStepTouchable = aS
223 const G4VTouchable* preStepTouchable = aSt
224 if ((preStepTouchable != nullptr) && (post
225 G4String post_vol_name = postStepTouchab
226 if (post_vol_name.empty()) {
227 post_vol_name = postStepTouchable->Get
228 }
229 G4String pre_vol_name = preStepTouchable
230 if (pre_vol_name.empty()) {
231 pre_vol_name = preStepTouchable->GetVo
232 }
233 if (pre_vol_name == vol1_name && post_vo
234 GoingIn = true;
235 did_cross = true;
236 }
237 else if (pre_vol_name == vol2_name && po
238 GoingIn = false;
239 did_cross = true;
240 }
241 }
242 }
243 return did_cross; // still need to compute
244 }
245
246 //////////////////////////////////////////////
247 //
248 G4bool G4AdjointCrossSurfChecker::AddaSpherica
249 const G4String& SurfaceName, G4double radius
250 {
251 G4int ind = FindRegisteredSurface(SurfaceNam
252 Area = 4. * pi * radius * radius;
253 if (ind >= 0) {
254 ListOfSurfaceType[ind] = "Sphere";
255 ListOfSphereRadius[ind] = radius;
256 ListOfSphereCenter[ind] = pos;
257 ListOfVol1Name[ind] = "";
258 ListOfVol2Name[ind] = "";
259 AreaOfSurface[ind] = Area;
260 }
261 else {
262 ListOfSurfaceName.push_back(SurfaceName);
263 ListOfSurfaceType.emplace_back("Sphere");
264 ListOfSphereRadius.push_back(radius);
265 ListOfSphereCenter.push_back(pos);
266 ListOfVol1Name.emplace_back("");
267 ListOfVol2Name.emplace_back("");
268 AreaOfSurface.push_back(Area);
269 }
270 return true;
271 }
272
273 //////////////////////////////////////////////
274 //
275 G4bool G4AdjointCrossSurfChecker::AddaSpherica
276 const G4String& SurfaceName, G4double radius
277 G4double& area)
278 {
279 G4VPhysicalVolume* thePhysicalVolume = nullp
280 G4PhysicalVolumeStore* thePhysVolStore = G4P
281 thePhysicalVolume = thePhysVolStore->GetVolu
282 if (thePhysicalVolume != nullptr) {
283 G4VPhysicalVolume* daughter = thePhysicalV
284 G4LogicalVolume* mother = thePhysicalVolum
285 G4AffineTransform theTransformationFromPhy
286 while (mother != nullptr) {
287 theTransformationFromPhysVolToWorld *=
288 G4AffineTransform(daughter->GetFrameRo
289 for (std::size_t i = 0; i < thePhysVolSt
290 if ((*thePhysVolStore)[i]->GetLogicalV
291 daughter = (*thePhysVolStore)[i];
292 mother = daughter->GetMotherLogical(
293 break;
294 }
295 }
296 }
297 center = theTransformationFromPhysVolToWor
298 G4cout << "Center of the spherical surface
299 << G4endl;
300 }
301 else {
302 return false;
303 }
304 return AddaSphericalSurface(SurfaceName, rad
305 }
306
307 //////////////////////////////////////////////
308 //
309 G4bool G4AdjointCrossSurfChecker::AddanExtSurf
310 const G4String& SurfaceName, const G4String&
311 {
312 G4int ind = FindRegisteredSurface(SurfaceNam
313
314 G4VPhysicalVolume* thePhysicalVolume = nullp
315 G4PhysicalVolumeStore* thePhysVolStore = G4P
316 thePhysicalVolume = thePhysVolStore->GetVolu
317 if (thePhysicalVolume == nullptr) {
318 return false;
319 }
320 Area = thePhysicalVolume->GetLogicalVolume()
321 G4String mother_vol_name = "";
322 G4LogicalVolume* theMother = thePhysicalVolu
323
324 if (theMother != nullptr) mother_vol_name =
325 if (ind >= 0) {
326 ListOfSurfaceType[ind] = "ExternalSurfaceO
327 ListOfSphereRadius[ind] = 0.;
328 ListOfSphereCenter[ind] = G4ThreeVector(0.
329 ListOfVol1Name[ind] = volume_name;
330 ListOfVol2Name[ind] = std::move(mother_vol
331 AreaOfSurface[ind] = Area;
332 }
333 else {
334 ListOfSurfaceName.push_back(SurfaceName);
335 ListOfSurfaceType.emplace_back("ExternalSu
336 ListOfSphereRadius.push_back(0.);
337 ListOfSphereCenter.emplace_back(0., 0., 0.
338 ListOfVol1Name.push_back(volume_name);
339 ListOfVol2Name.push_back(std::move(mother_
340 AreaOfSurface.push_back(Area);
341 }
342 return true;
343 }
344
345 //////////////////////////////////////////////
346 //
347 G4bool G4AdjointCrossSurfChecker::AddanInterfa
348 const G4String& volume_name1, const G4String
349 {
350 G4int ind = FindRegisteredSurface(SurfaceNam
351 Area = -1.; // the way to compute the surfa
352 if (ind >= 0) {
353 ListOfSurfaceType[ind] = "BoundaryBetweenT
354 ListOfSphereRadius[ind] = 0.;
355 ListOfSphereCenter[ind] = G4ThreeVector(0.
356 ListOfVol1Name[ind] = volume_name1;
357 ListOfVol2Name[ind] = volume_name2;
358 AreaOfSurface[ind] = Area;
359 }
360 else {
361 ListOfSurfaceName.push_back(SurfaceName);
362 ListOfSurfaceType.emplace_back("BoundaryBe
363 ListOfSphereRadius.push_back(0.);
364 ListOfSphereCenter.emplace_back(0., 0., 0.
365 ListOfVol1Name.push_back(volume_name1);
366 ListOfVol2Name.push_back(volume_name2);
367 AreaOfSurface.push_back(Area);
368 }
369 return true;
370 }
371
372 //////////////////////////////////////////////
373 //
374 void G4AdjointCrossSurfChecker::ClearListOfSel
375 {
376 ListOfSurfaceName.clear();
377 ListOfSurfaceType.clear();
378 ListOfSphereRadius.clear();
379 ListOfSphereCenter.clear();
380 ListOfVol1Name.clear();
381 ListOfVol2Name.clear();
382 }
383
384 //////////////////////////////////////////////
385 //
386 G4int G4AdjointCrossSurfChecker::FindRegistere
387 {
388 for (std::size_t i = 0; i < ListOfSurfaceNam
389 if (name == ListOfSurfaceName[i]) return G
390 }
391 return -1;
392 }
393