| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // G4AdjointCrossSurfChecker class implementation
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* G4AdjointCrossSurfChecker::instance = nullptr;
46
47 //////////////////////////////////////////////////////////////////////////////
48 //
49 G4AdjointCrossSurfChecker::~G4AdjointCrossSurfChecker() { delete instance; }
50
51 //////////////////////////////////////////////////////////////////////////////
52 //
53 G4AdjointCrossSurfChecker* G4AdjointCrossSurfChecker::GetInstance()
54 {
55 if (instance == nullptr) instance = new G4AdjointCrossSurfChecker();
56 return instance;
57 }
58
59 //////////////////////////////////////////////////////////////////////////////
60 //
61 G4bool G4AdjointCrossSurfChecker::CrossingASphere(const G4Step* aStep, G4double sphere_radius,
62 G4ThreeVector sphere_center, G4ThreeVector& crossing_pos, G4double& cos_th, G4bool& GoingIn)
63 {
64 G4ThreeVector pos1 = aStep->GetPreStepPoint()->GetPosition() - sphere_center;
65 G4ThreeVector pos2 = aStep->GetPostStepPoint()->GetPosition() - sphere_center;
66 G4double r1 = pos1.mag();
67 G4double r2 = pos2.mag();
68 G4bool did_cross = false;
69
70 if (r1 <= sphere_radius && r2 > sphere_radius) {
71 did_cross = true;
72 GoingIn = false;
73 }
74 else if (r2 <= sphere_radius && r1 > sphere_radius) {
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::GoingInOrOutOfaVolume(const G4Step* aStep,
99 const G4String& volume_name, G4double&, G4bool& GoingIn) // from external surface
100 {
101 G4bool step_at_boundary = (aStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary);
102 G4bool did_cross = false;
103 if (step_at_boundary) {
104 const G4VTouchable* postStepTouchable = aStep->GetPostStepPoint()->GetTouchable();
105 const G4VTouchable* preStepTouchable = aStep->GetPreStepPoint()->GetTouchable();
106 if ((preStepTouchable != nullptr) && (postStepTouchable != nullptr) &&
107 (postStepTouchable->GetVolume() != nullptr) && (preStepTouchable->GetVolume() != nullptr))
108 {
109 G4String post_vol_name = postStepTouchable->GetVolume()->GetName();
110 G4String pre_vol_name = preStepTouchable->GetVolume()->GetName();
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 the cosine of the direction
123 }
124
125 /////////////////////////////////////////////////////////////////////////////
126 //
127 G4bool G4AdjointCrossSurfChecker::GoingInOrOutOfaVolumeByExtSurface(const G4Step* aStep,
128 const G4String& volume_name, const G4String& mother_logical_vol_name, G4double&,
129 G4bool& GoingIn) // from external surf.
130 {
131 G4bool step_at_boundary = (aStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary);
132 G4bool did_cross = false;
133 if (step_at_boundary) {
134 const G4VTouchable* postStepTouchable = aStep->GetPostStepPoint()->GetTouchable();
135 const G4VTouchable* preStepTouchable = aStep->GetPreStepPoint()->GetTouchable();
136 const G4VPhysicalVolume* postVol =
137 (postStepTouchable != nullptr) ? postStepTouchable->GetVolume() : nullptr;
138 const G4VPhysicalVolume* preVol =
139 (preStepTouchable != nullptr) ? preStepTouchable->GetVolume() : nullptr;
140 if (preStepTouchable != nullptr && postStepTouchable != nullptr && postVol != nullptr &&
141 preVol != nullptr)
142 {
143 G4String post_vol_name = postVol->GetName();
144 G4String post_log_vol_name = postVol->GetLogicalVolume()->GetName();
145 G4String pre_vol_name = preVol->GetName();
146 G4String pre_log_vol_name = preVol->GetLogicalVolume()->GetName();
147 if (post_vol_name == volume_name && pre_log_vol_name == mother_logical_vol_name) {
148 GoingIn = true;
149 did_cross = true;
150 }
151 else if (pre_vol_name == volume_name && post_log_vol_name == mother_logical_vol_name) {
152 GoingIn = false;
153 did_cross = true;
154 }
155 }
156 }
157 return did_cross; // still need to compute the cosine of the direction
158 }
159
160 //////////////////////////////////////////////////////////////////////////////
161 //
162 G4bool G4AdjointCrossSurfChecker::CrossingAGivenRegisteredSurface(const G4Step* aStep,
163 const G4String& surface_name, G4ThreeVector& crossing_pos, G4double& cos_to_surface,
164 G4bool& GoingIn)
165 {
166 G4int ind = FindRegisteredSurface(surface_name);
167 G4bool did_cross = false;
168 if (ind >= 0) {
169 did_cross = CrossingAGivenRegisteredSurface(aStep, ind, crossing_pos, cos_to_surface, GoingIn);
170 }
171 return did_cross;
172 }
173
174 //////////////////////////////////////////////////////////////////////////////
175 //
176 G4bool G4AdjointCrossSurfChecker::CrossingAGivenRegisteredSurface(const G4Step* aStep, G4int ind,
177 G4ThreeVector& crossing_pos, G4double& cos_to_surface, G4bool& GoingIn)
178 {
179 G4String surf_type = ListOfSurfaceType[ind];
180 G4double radius = ListOfSphereRadius[ind];
181 G4ThreeVector center = ListOfSphereCenter[ind];
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, center, crossing_pos, cos_to_surface, GoingIn);
188 }
189 else if (surf_type == "ExternalSurfaceOfAVolume") {
190 did_cross = GoingInOrOutOfaVolumeByExtSurface(aStep, vol1, vol2, cos_to_surface, GoingIn);
191 crossing_pos = aStep->GetPostStepPoint()->GetPosition();
192 }
193 else if (surf_type == "BoundaryBetweenTwoVolumes") {
194 did_cross = CrossingAnInterfaceBetweenTwoVolumes(
195 aStep, vol1, vol2, crossing_pos, cos_to_surface, GoingIn);
196 }
197 return did_cross;
198 }
199
200 //////////////////////////////////////////////////////////////////////////////
201 //
202 G4bool G4AdjointCrossSurfChecker::CrossingOneOfTheRegisteredSurface(const G4Step* aStep,
203 G4String& surface_name, G4ThreeVector& crossing_pos, G4double& cos_to_surface, G4bool& GoingIn)
204 {
205 for (std::size_t i = 0; i < ListOfSurfaceName.size(); ++i) {
206 if (CrossingAGivenRegisteredSurface(aStep, G4int(i), crossing_pos, cos_to_surface, GoingIn)) {
207 surface_name = ListOfSurfaceName[i];
208 return true;
209 }
210 }
211 return false;
212 }
213
214 //////////////////////////////////////////////////////////////////////////////
215 //
216 G4bool G4AdjointCrossSurfChecker::CrossingAnInterfaceBetweenTwoVolumes(const G4Step* aStep,
217 const G4String& vol1_name, const G4String& vol2_name, G4ThreeVector&, G4double&, G4bool& GoingIn)
218 {
219 G4bool step_at_boundary = (aStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary);
220 G4bool did_cross = false;
221 if (step_at_boundary) {
222 const G4VTouchable* postStepTouchable = aStep->GetPostStepPoint()->GetTouchable();
223 const G4VTouchable* preStepTouchable = aStep->GetPreStepPoint()->GetTouchable();
224 if ((preStepTouchable != nullptr) && (postStepTouchable != nullptr)) {
225 G4String post_vol_name = postStepTouchable->GetVolume()->GetName();
226 if (post_vol_name.empty()) {
227 post_vol_name = postStepTouchable->GetVolume()->GetLogicalVolume()->GetName();
228 }
229 G4String pre_vol_name = preStepTouchable->GetVolume()->GetName();
230 if (pre_vol_name.empty()) {
231 pre_vol_name = preStepTouchable->GetVolume()->GetLogicalVolume()->GetName();
232 }
233 if (pre_vol_name == vol1_name && post_vol_name == vol2_name) {
234 GoingIn = true;
235 did_cross = true;
236 }
237 else if (pre_vol_name == vol2_name && post_vol_name == vol1_name) {
238 GoingIn = false;
239 did_cross = true;
240 }
241 }
242 }
243 return did_cross; // still need to compute the cosine of the direction
244 }
245
246 //////////////////////////////////////////////////////////////////////////////
247 //
248 G4bool G4AdjointCrossSurfChecker::AddaSphericalSurface(
249 const G4String& SurfaceName, G4double radius, G4ThreeVector pos, G4double& Area)
250 {
251 G4int ind = FindRegisteredSurface(SurfaceName);
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::AddaSphericalSurfaceWithCenterAtTheCenterOfAVolume(
276 const G4String& SurfaceName, G4double radius, const G4String& volume_name, G4ThreeVector& center,
277 G4double& area)
278 {
279 G4VPhysicalVolume* thePhysicalVolume = nullptr;
280 G4PhysicalVolumeStore* thePhysVolStore = G4PhysicalVolumeStore::GetInstance();
281 thePhysicalVolume = thePhysVolStore->GetVolume(volume_name);
282 if (thePhysicalVolume != nullptr) {
283 G4VPhysicalVolume* daughter = thePhysicalVolume;
284 G4LogicalVolume* mother = thePhysicalVolume->GetMotherLogical();
285 G4AffineTransform theTransformationFromPhysVolToWorld = G4AffineTransform();
286 while (mother != nullptr) {
287 theTransformationFromPhysVolToWorld *=
288 G4AffineTransform(daughter->GetFrameRotation(), daughter->GetObjectTranslation());
289 for (std::size_t i = 0; i < thePhysVolStore->size(); ++i) {
290 if ((*thePhysVolStore)[i]->GetLogicalVolume() == mother) {
291 daughter = (*thePhysVolStore)[i];
292 mother = daughter->GetMotherLogical();
293 break;
294 }
295 }
296 }
297 center = theTransformationFromPhysVolToWorld.NetTranslation();
298 G4cout << "Center of the spherical surface is at the position: " << center / cm << " cm"
299 << G4endl;
300 }
301 else {
302 return false;
303 }
304 return AddaSphericalSurface(SurfaceName, radius, center, area);
305 }
306
307 //////////////////////////////////////////////////////////////////////////////
308 //
309 G4bool G4AdjointCrossSurfChecker::AddanExtSurfaceOfAvolume(
310 const G4String& SurfaceName, const G4String& volume_name, G4double& Area)
311 {
312 G4int ind = FindRegisteredSurface(SurfaceName);
313
314 G4VPhysicalVolume* thePhysicalVolume = nullptr;
315 G4PhysicalVolumeStore* thePhysVolStore = G4PhysicalVolumeStore::GetInstance();
316 thePhysicalVolume = thePhysVolStore->GetVolume(volume_name);
317 if (thePhysicalVolume == nullptr) {
318 return false;
319 }
320 Area = thePhysicalVolume->GetLogicalVolume()->GetSolid()->GetSurfaceArea();
321 G4String mother_vol_name = "";
322 G4LogicalVolume* theMother = thePhysicalVolume->GetMotherLogical();
323
324 if (theMother != nullptr) mother_vol_name = theMother->GetName();
325 if (ind >= 0) {
326 ListOfSurfaceType[ind] = "ExternalSurfaceOfAVolume";
327 ListOfSphereRadius[ind] = 0.;
328 ListOfSphereCenter[ind] = G4ThreeVector(0., 0., 0.);
329 ListOfVol1Name[ind] = volume_name;
330 ListOfVol2Name[ind] = std::move(mother_vol_name);
331 AreaOfSurface[ind] = Area;
332 }
333 else {
334 ListOfSurfaceName.push_back(SurfaceName);
335 ListOfSurfaceType.emplace_back("ExternalSurfaceOfAVolume");
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_vol_name));
340 AreaOfSurface.push_back(Area);
341 }
342 return true;
343 }
344
345 //////////////////////////////////////////////////////////////////////////////
346 //
347 G4bool G4AdjointCrossSurfChecker::AddanInterfaceBetweenTwoVolumes(const G4String& SurfaceName,
348 const G4String& volume_name1, const G4String& volume_name2, G4double& Area)
349 {
350 G4int ind = FindRegisteredSurface(SurfaceName);
351 Area = -1.; // the way to compute the surface is not known yet
352 if (ind >= 0) {
353 ListOfSurfaceType[ind] = "BoundaryBetweenTwoVolumes";
354 ListOfSphereRadius[ind] = 0.;
355 ListOfSphereCenter[ind] = G4ThreeVector(0., 0., 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("BoundaryBetweenTwoVolumes");
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::ClearListOfSelectedSurface()
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::FindRegisteredSurface(const G4String& name)
387 {
388 for (std::size_t i = 0; i < ListOfSurfaceName.size(); ++i) {
389 if (name == ListOfSurfaceName[i]) return G4int(i);
390 }
391 return -1;
392 }
393