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