| 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 // class G4PVPlacement Implementation 26 // class G4PVPlacement Implementation
27 // 27 //
28 // 24.07.95 P.Kent, First non-stub version. <<
29 // ------------------------------------------- 28 // ----------------------------------------------------------------------
30 29
31 #include "G4PVPlacement.hh" 30 #include "G4PVPlacement.hh"
32 #include "G4AffineTransform.hh" 31 #include "G4AffineTransform.hh"
33 #include "G4UnitsTable.hh" 32 #include "G4UnitsTable.hh"
34 #include "G4LogicalVolume.hh" 33 #include "G4LogicalVolume.hh"
35 #include "G4VSolid.hh" 34 #include "G4VSolid.hh"
36 35
37 // ------------------------------------------- 36 // ----------------------------------------------------------------------
38 // Constructor 37 // Constructor
39 // 38 //
40 G4PVPlacement::G4PVPlacement( G4RotationMatrix 39 G4PVPlacement::G4PVPlacement( G4RotationMatrix* pRot,
41 const G4ThreeVector& t 40 const G4ThreeVector& tlate,
42 const G4String& pName, 41 const G4String& pName,
43 G4LogicalVolume* 42 G4LogicalVolume* pLogical,
44 G4VPhysicalVolum 43 G4VPhysicalVolume* pMother,
45 G4bool pMany, 44 G4bool pMany,
46 G4int pCopyNo, 45 G4int pCopyNo,
47 G4bool pSurfChk 46 G4bool pSurfChk )
48 : G4VPhysicalVolume(pRot, tlate, pName, pLog 47 : G4VPhysicalVolume(pRot, tlate, pName, pLogical, pMother),
49 fmany(pMany), fcopyNo(pCopyNo) 48 fmany(pMany), fcopyNo(pCopyNo)
50 { 49 {
51 if (pMother != nullptr) << 50 if (pMother)
52 { 51 {
53 G4LogicalVolume* motherLogical = pMother-> 52 G4LogicalVolume* motherLogical = pMother->GetLogicalVolume();
54 if (pLogical == motherLogical) 53 if (pLogical == motherLogical)
55 { 54 {
56 G4Exception("G4PVPlacement::G4PVPlacemen 55 G4Exception("G4PVPlacement::G4PVPlacement()", "GeomVol0002",
57 FatalException, "Cannot plac 56 FatalException, "Cannot place a volume inside itself!");
58 } 57 }
59 SetMotherLogical(motherLogical); 58 SetMotherLogical(motherLogical);
60 motherLogical->AddDaughter(this); 59 motherLogical->AddDaughter(this);
61 if (pSurfChk) { CheckOverlaps(); } 60 if (pSurfChk) { CheckOverlaps(); }
62 } 61 }
63 } 62 }
64 63
65 // ------------------------------------------- 64 // ----------------------------------------------------------------------
66 // Constructor 65 // Constructor
67 // 66 //
68 G4PVPlacement::G4PVPlacement( const G4Transfor 67 G4PVPlacement::G4PVPlacement( const G4Transform3D& Transform3D,
69 const G4String& 68 const G4String& pName,
70 G4LogicalV 69 G4LogicalVolume* pLogical,
71 G4VPhysica 70 G4VPhysicalVolume* pMother,
72 G4bool pMa 71 G4bool pMany,
73 G4int pCop 72 G4int pCopyNo,
74 G4bool pSu 73 G4bool pSurfChk )
75 : G4VPhysicalVolume(NewPtrRotMatrix(Transfor 74 : G4VPhysicalVolume(NewPtrRotMatrix(Transform3D.getRotation().inverse()),
76 Transform3D.getTranslati 75 Transform3D.getTranslation(), pName, pLogical, pMother),
77 fmany(pMany), fcopyNo(pCopyNo) 76 fmany(pMany), fcopyNo(pCopyNo)
78 { 77 {
79 fallocatedRotM = (GetRotation() != nullptr); << 78 fallocatedRotM = (GetRotation() != 0);
80 if (pMother != nullptr) << 79 if (pMother)
81 { 80 {
82 G4LogicalVolume* motherLogical = pMother-> 81 G4LogicalVolume* motherLogical = pMother->GetLogicalVolume();
83 if (pLogical == motherLogical) 82 if (pLogical == motherLogical)
84 G4Exception("G4PVPlacement::G4PVPlacemen 83 G4Exception("G4PVPlacement::G4PVPlacement()", "GeomVol0002",
85 FatalException, "Cannot plac 84 FatalException, "Cannot place a volume inside itself!");
86 SetMotherLogical(motherLogical); 85 SetMotherLogical(motherLogical);
87 motherLogical->AddDaughter(this); 86 motherLogical->AddDaughter(this);
88 if (pSurfChk) { CheckOverlaps(); } 87 if (pSurfChk) { CheckOverlaps(); }
89 } 88 }
90 } 89 }
91 90
92 // ------------------------------------------- 91 // ----------------------------------------------------------------------
93 // Constructor 92 // Constructor
94 // 93 //
95 // The logical volume of the mother is utilise 94 // The logical volume of the mother is utilised (not the physical)
96 // 95 //
97 G4PVPlacement::G4PVPlacement( G4RotationMatrix 96 G4PVPlacement::G4PVPlacement( G4RotationMatrix* pRot,
98 const G4ThreeVector& t 97 const G4ThreeVector& tlate,
99 G4LogicalVolume* 98 G4LogicalVolume* pCurrentLogical,
100 const G4String& pName, 99 const G4String& pName,
101 G4LogicalVolume* 100 G4LogicalVolume* pMotherLogical,
102 G4bool pMany, 101 G4bool pMany,
103 G4int pCopyNo, 102 G4int pCopyNo,
104 G4bool pSurfChk 103 G4bool pSurfChk )
105 : G4VPhysicalVolume(pRot, tlate, pName, pCur 104 : G4VPhysicalVolume(pRot, tlate, pName, pCurrentLogical, nullptr),
106 fmany(pMany), fcopyNo(pCopyNo) 105 fmany(pMany), fcopyNo(pCopyNo)
107 { 106 {
108 if (pCurrentLogical == pMotherLogical) 107 if (pCurrentLogical == pMotherLogical)
109 { 108 {
110 G4Exception("G4PVPlacement::G4PVPlacement( 109 G4Exception("G4PVPlacement::G4PVPlacement()", "GeomVol0002",
111 FatalException, "Cannot place 110 FatalException, "Cannot place a volume inside itself!");
112 } 111 }
113 SetMotherLogical(pMotherLogical); 112 SetMotherLogical(pMotherLogical);
114 if (pMotherLogical != nullptr) { pMotherLogi << 113 if (pMotherLogical) { pMotherLogical->AddDaughter(this); }
115 if ((pSurfChk) && ((pMotherLogical) != nullp << 114 if ((pSurfChk) && (pMotherLogical)) { CheckOverlaps(); }
116 } 115 }
117 116
>> 117
118 // ------------------------------------------- 118 // ----------------------------------------------------------------------
119 // Constructor 119 // Constructor
120 // 120 //
121 G4PVPlacement::G4PVPlacement( const G4Transfor 121 G4PVPlacement::G4PVPlacement( const G4Transform3D& Transform3D,
122 G4LogicalV 122 G4LogicalVolume* pCurrentLogical,
123 const G4String& 123 const G4String& pName,
124 G4LogicalV 124 G4LogicalVolume* pMotherLogical,
125 G4bool pMa 125 G4bool pMany,
126 G4int pCop 126 G4int pCopyNo,
127 G4bool pSu 127 G4bool pSurfChk )
128 : G4VPhysicalVolume(nullptr, Transform3D.get 128 : G4VPhysicalVolume(nullptr, Transform3D.getTranslation(),
129 pName, pCurrentLogical, 129 pName, pCurrentLogical, nullptr),
130 fmany(pMany), fcopyNo(pCopyNo) 130 fmany(pMany), fcopyNo(pCopyNo)
131 { 131 {
132 if (pCurrentLogical == pMotherLogical) 132 if (pCurrentLogical == pMotherLogical)
133 { 133 {
134 G4Exception("G4PVPlacement::G4PVPlacement( 134 G4Exception("G4PVPlacement::G4PVPlacement()", "GeomVol0002",
135 FatalException, "Cannot place 135 FatalException, "Cannot place a volume inside itself!");
136 } 136 }
137 SetRotation( NewPtrRotMatrix(Transform3D.get 137 SetRotation( NewPtrRotMatrix(Transform3D.getRotation().inverse()) );
138 fallocatedRotM = (GetRotation() != nullptr); 138 fallocatedRotM = (GetRotation() != nullptr);
139 SetMotherLogical(pMotherLogical); 139 SetMotherLogical(pMotherLogical);
140 if (pMotherLogical != nullptr) { pMotherLogi << 140 if (pMotherLogical) { pMotherLogical->AddDaughter(this); }
141 if ((pSurfChk) && ((pMotherLogical) != nullp << 141 if ((pSurfChk) && (pMotherLogical)) { CheckOverlaps(); }
142 } 142 }
143 143
144 // ------------------------------------------- 144 // ----------------------------------------------------------------------
145 // Fake default constructor - sets only member 145 // Fake default constructor - sets only member data and allocates memory
146 // for usage restri 146 // for usage restricted to object persistency.
147 // 147 //
148 G4PVPlacement::G4PVPlacement( __void__& a ) 148 G4PVPlacement::G4PVPlacement( __void__& a )
149 : G4VPhysicalVolume(a) 149 : G4VPhysicalVolume(a)
150 { 150 {
151 } 151 }
152 152
153 // ------------------------------------------- 153 // ----------------------------------------------------------------------
154 // Destructor 154 // Destructor
155 // 155 //
156 G4PVPlacement::~G4PVPlacement() 156 G4PVPlacement::~G4PVPlacement()
157 { 157 {
158 if( fallocatedRotM ){ delete this->GetRotati 158 if( fallocatedRotM ){ delete this->GetRotation() ; }
159 } 159 }
160 160
161 // ------------------------------------------- 161 // ----------------------------------------------------------------------
162 // IsMany 162 // IsMany
163 // 163 //
164 G4bool G4PVPlacement::IsMany() const 164 G4bool G4PVPlacement::IsMany() const
165 { 165 {
166 return fmany; 166 return fmany;
167 } 167 }
168 168
169 // ------------------------------------------- 169 // ----------------------------------------------------------------------
170 // SetCopyNo 170 // SetCopyNo
171 // 171 //
172 void G4PVPlacement::SetCopyNo(G4int newCopyNo) 172 void G4PVPlacement::SetCopyNo(G4int newCopyNo)
173 { 173 {
174 fcopyNo = newCopyNo; 174 fcopyNo = newCopyNo;
175 } 175 }
176 176
177 // ------------------------------------------- 177 // ----------------------------------------------------------------------
178 // IsReplicated 178 // IsReplicated
179 // 179 //
180 G4bool G4PVPlacement::IsReplicated() const 180 G4bool G4PVPlacement::IsReplicated() const
181 { 181 {
182 return false; 182 return false;
183 } 183 }
184 184
185 // ------------------------------------------- 185 // ----------------------------------------------------------------------
186 // IsParameterised 186 // IsParameterised
187 // 187 //
188 G4bool G4PVPlacement::IsParameterised() const 188 G4bool G4PVPlacement::IsParameterised() const
189 { 189 {
190 return false; 190 return false;
191 } 191 }
192 192
193 // ------------------------------------------- 193 // ----------------------------------------------------------------------
194 // GetParameterisation 194 // GetParameterisation
195 // 195 //
196 G4VPVParameterisation* G4PVPlacement::GetParam 196 G4VPVParameterisation* G4PVPlacement::GetParameterisation() const
197 { 197 {
198 return nullptr; 198 return nullptr;
199 } 199 }
200 200
201 // ------------------------------------------- 201 // ----------------------------------------------------------------------
202 // GetReplicationData 202 // GetReplicationData
203 // 203 //
204 void G4PVPlacement:: 204 void G4PVPlacement::
205 GetReplicationData( EAxis&, G4int&, G4double&, 205 GetReplicationData( EAxis&, G4int&, G4double&, G4double&, G4bool& ) const
206 { 206 {
207 // No-operations 207 // No-operations
208 } 208 }
209 209
210 // ------------------------------------------- 210 // ----------------------------------------------------------------------
211 // IsRegularRepeatedStructure 211 // IsRegularRepeatedStructure
212 // 212 //
213 // This is for specialised repeated volumes (r 213 // This is for specialised repeated volumes (replicas, parameterised vol.)
214 // 214 //
215 G4bool G4PVPlacement::IsRegularStructure() con 215 G4bool G4PVPlacement::IsRegularStructure() const
216 { 216 {
217 return false; 217 return false;
218 } 218 }
219 219
220 // ------------------------------------------- 220 // ----------------------------------------------------------------------
221 // IsRegularRepeatedStructure 221 // IsRegularRepeatedStructure
222 // 222 //
223 // This is for specialised repeated volumes (r 223 // This is for specialised repeated volumes (replicas, parameterised vol.)
224 // 224 //
225 G4int G4PVPlacement::GetRegularStructureId() c 225 G4int G4PVPlacement::GetRegularStructureId() const
226 { 226 {
227 return 0; 227 return 0;
228 } 228 }
229 229
230 // ------------------------------------------- 230 // ----------------------------------------------------------------------
231 // VolumeType 231 // VolumeType
232 // 232 //
233 // Information to help identify sub-navigator 233 // Information to help identify sub-navigator which will be used
234 // 234 //
235 EVolume G4PVPlacement::VolumeType() const 235 EVolume G4PVPlacement::VolumeType() const
236 { 236 {
237 return kNormal; 237 return kNormal;
238 } 238 }
239 239
240 // ------------------------------------------- 240 // ----------------------------------------------------------------------
241 // CheckOverlaps 241 // CheckOverlaps
242 // 242 //
243 G4bool G4PVPlacement::CheckOverlaps(G4int res, 243 G4bool G4PVPlacement::CheckOverlaps(G4int res, G4double tol,
244 G4bool ver 244 G4bool verbose, G4int maxErr)
245 { 245 {
246 if (res <= 0) { return false; } 246 if (res <= 0) { return false; }
247 247
248 G4VSolid* solid = GetLogicalVolume()->GetSol 248 G4VSolid* solid = GetLogicalVolume()->GetSolid();
249 G4LogicalVolume* motherLog = GetMotherLogica 249 G4LogicalVolume* motherLog = GetMotherLogical();
250 if (motherLog == nullptr) { return false; } 250 if (motherLog == nullptr) { return false; }
251 251
252 G4int trials = 0; 252 G4int trials = 0;
253 G4bool retval = false; 253 G4bool retval = false;
254 254
255 if (verbose) 255 if (verbose)
256 { 256 {
257 G4cout << "Checking overlaps for volume " 257 G4cout << "Checking overlaps for volume "
258 << GetName() << ':' << GetCopyNo() 258 << GetName() << ':' << GetCopyNo()
259 << " (" << solid->GetEntityType() < 259 << " (" << solid->GetEntityType() << ") ... ";
260 } 260 }
261 261
262 // Check that random points are gererated co 262 // Check that random points are gererated correctly
263 // 263 //
264 G4ThreeVector ptmp = solid->GetPointOnSurfac 264 G4ThreeVector ptmp = solid->GetPointOnSurface();
265 if (solid->Inside(ptmp) != kSurface) 265 if (solid->Inside(ptmp) != kSurface)
266 { 266 {
267 G4String position[3] = { "outside", "surfa 267 G4String position[3] = { "outside", "surface", "inside" };
268 std::ostringstream message; 268 std::ostringstream message;
269 message << "Sample point is not on the sur 269 message << "Sample point is not on the surface !" << G4endl
270 << " The issue is detecte 270 << " The issue is detected for volume "
271 << GetName() << ':' << GetCopyNo() 271 << GetName() << ':' << GetCopyNo()
272 << " (" << solid->GetEntityType() 272 << " (" << solid->GetEntityType() << ")" << G4endl
273 << " generated point " << 273 << " generated point " << ptmp
274 << " is " << position[solid->Insid 274 << " is " << position[solid->Inside(ptmp)];
275 G4Exception("G4PVPlacement::CheckOverlaps( 275 G4Exception("G4PVPlacement::CheckOverlaps()",
276 "GeomVol1002", JustWarning, me 276 "GeomVol1002", JustWarning, message);
277 return false; 277 return false;
278 } 278 }
279 279
280 // Generate random points on the surface of 280 // Generate random points on the surface of the solid,
281 // transform them into the mother volume coo 281 // transform them into the mother volume coordinate system
282 // and find the bonding box 282 // and find the bonding box
283 // 283 //
284 std::vector<G4ThreeVector> points(res); 284 std::vector<G4ThreeVector> points(res);
285 G4double xmin = kInfinity, ymin = kInfinit 285 G4double xmin = kInfinity, ymin = kInfinity, zmin = kInfinity;
286 G4double xmax = -kInfinity, ymax = -kInfinit 286 G4double xmax = -kInfinity, ymax = -kInfinity, zmax = -kInfinity;
287 G4AffineTransform Tm(GetRotation(), GetTrans 287 G4AffineTransform Tm(GetRotation(), GetTranslation());
288 for (G4int i = 0; i < res; ++i) 288 for (G4int i = 0; i < res; ++i)
289 { 289 {
290 points[i] = Tm.TransformPoint(solid->GetPo 290 points[i] = Tm.TransformPoint(solid->GetPointOnSurface());
291 xmin = std::min(xmin, points[i].x()); 291 xmin = std::min(xmin, points[i].x());
292 ymin = std::min(ymin, points[i].y()); 292 ymin = std::min(ymin, points[i].y());
293 zmin = std::min(zmin, points[i].z()); 293 zmin = std::min(zmin, points[i].z());
294 xmax = std::max(xmax, points[i].x()); 294 xmax = std::max(xmax, points[i].x());
295 ymax = std::max(ymax, points[i].y()); 295 ymax = std::max(ymax, points[i].y());
296 zmax = std::max(zmax, points[i].z()); 296 zmax = std::max(zmax, points[i].z());
297 } 297 }
298 G4ThreeVector scenter(0.5*(xmax+xmin), 0.5*( 298 G4ThreeVector scenter(0.5*(xmax+xmin), 0.5*(ymax+ymin), 0.5*(zmax+zmin));
299 G4double sradius = 0.5*G4ThreeVector(xmax-xm 299 G4double sradius = 0.5*G4ThreeVector(xmax-xmin, ymax-ymin, zmax-zmin).mag();
300 300
301 // Check overlap with the mother volume 301 // Check overlap with the mother volume
302 // 302 //
303 G4int overlapCount = 0; 303 G4int overlapCount = 0;
304 G4double overlapSize = -kInfinity; 304 G4double overlapSize = -kInfinity;
305 G4ThreeVector overlapPoint; 305 G4ThreeVector overlapPoint;
306 G4VSolid* motherSolid = motherLog->GetSolid( 306 G4VSolid* motherSolid = motherLog->GetSolid();
307 for (G4int i = 0; i < res; ++i) 307 for (G4int i = 0; i < res; ++i)
308 { 308 {
309 G4ThreeVector mp = points[i]; 309 G4ThreeVector mp = points[i];
310 if (motherSolid->Inside(mp) != kOutside) c 310 if (motherSolid->Inside(mp) != kOutside) continue;
311 G4double distin = motherSolid->DistanceToI 311 G4double distin = motherSolid->DistanceToIn(mp);
312 if (distin < tol) continue; // too small o 312 if (distin < tol) continue; // too small overlap
313 ++overlapCount; 313 ++overlapCount;
314 if (distin <= overlapSize) continue; 314 if (distin <= overlapSize) continue;
315 overlapSize = distin; 315 overlapSize = distin;
316 overlapPoint = mp; 316 overlapPoint = mp;
317 } 317 }
318 318
319 // Print information on overlap 319 // Print information on overlap
320 // 320 //
321 if (overlapCount > 0) 321 if (overlapCount > 0)
322 { 322 {
323 ++trials; 323 ++trials;
324 retval = true; 324 retval = true;
325 std::ostringstream message; 325 std::ostringstream message;
326 message << "Overlap with mother volume !" 326 message << "Overlap with mother volume !" << G4endl
327 << " Overlap is detected 327 << " Overlap is detected for volume "
328 << GetName() << ':' << GetCopyNo() 328 << GetName() << ':' << GetCopyNo()
329 << " (" << solid->GetEntityType() 329 << " (" << solid->GetEntityType() << ")"
330 << " with its mother volume " << m 330 << " with its mother volume " << motherLog->GetName()
331 << " (" << motherSolid->GetEntityT 331 << " (" << motherSolid->GetEntityType() << ")" << G4endl
332 << " protrusion at mother 332 << " protrusion at mother local point " << overlapPoint
333 << " by " << G4BestUnit(overlapSiz 333 << " by " << G4BestUnit(overlapSize, "Length")
334 << " (max of " << overlapCount << 334 << " (max of " << overlapCount << " cases)";
335 if (trials >= maxErr) 335 if (trials >= maxErr)
336 { 336 {
337 message << G4endl 337 message << G4endl
338 << "NOTE: Reached maximum fixed 338 << "NOTE: Reached maximum fixed number -" << maxErr
339 << "- of overlaps reports for th 339 << "- of overlaps reports for this volume !";
340 } 340 }
341 G4Exception("G4PVPlacement::CheckOverlaps( 341 G4Exception("G4PVPlacement::CheckOverlaps()",
342 "GeomVol1002", JustWarning, me 342 "GeomVol1002", JustWarning, message);
343 if (trials >= maxErr) { return true; } 343 if (trials >= maxErr) { return true; }
344 } 344 }
345 345
346 // Checking overlaps with each 'sister' volu 346 // Checking overlaps with each 'sister' volumes
347 // 347 //
348 G4VSolid* previous = nullptr; 348 G4VSolid* previous = nullptr;
349 G4ThreeVector pmin_local(0.,0.,0.), pmax_loc 349 G4ThreeVector pmin_local(0.,0.,0.), pmax_local(0.,0.,0.);
350 350
351 for (std::size_t k = 0; k < motherLog->GetNo 351 for (std::size_t k = 0; k < motherLog->GetNoDaughters(); ++k)
352 { 352 {
353 G4VPhysicalVolume* daughter = motherLog->G 353 G4VPhysicalVolume* daughter = motherLog->GetDaughter((G4int)k);
354 if (daughter == this) continue; 354 if (daughter == this) continue;
355 G4bool check_encapsulation = true; 355 G4bool check_encapsulation = true;
356 356
357 G4AffineTransform Td(daughter->GetRotation 357 G4AffineTransform Td(daughter->GetRotation(), daughter->GetTranslation());
358 G4VSolid* daughterSolid = daughter->GetLog 358 G4VSolid* daughterSolid = daughter->GetLogicalVolume()->GetSolid();
359 if (previous != daughterSolid) 359 if (previous != daughterSolid)
360 { 360 {
361 daughterSolid->BoundingLimits(pmin_local 361 daughterSolid->BoundingLimits(pmin_local, pmax_local);
362 previous = daughterSolid; 362 previous = daughterSolid;
363 } 363 }
364 overlapCount = 0; 364 overlapCount = 0;
365 overlapSize = -kInfinity; 365 overlapSize = -kInfinity;
366 if (!Td.IsRotated()) { // no rotation, onl 366 if (!Td.IsRotated()) { // no rotation, only translation
367 G4ThreeVector offset = Td.NetTranslation 367 G4ThreeVector offset = Td.NetTranslation();
368 G4ThreeVector pmin(pmin_local + offset); 368 G4ThreeVector pmin(pmin_local + offset);
369 G4ThreeVector pmax(pmax_local + offset); 369 G4ThreeVector pmax(pmax_local + offset);
370 if (pmin.x() >= xmax) continue; 370 if (pmin.x() >= xmax) continue;
371 if (pmin.y() >= ymax) continue; 371 if (pmin.y() >= ymax) continue;
372 if (pmin.z() >= zmax) continue; 372 if (pmin.z() >= zmax) continue;
373 if (pmax.x() <= xmin) continue; 373 if (pmax.x() <= xmin) continue;
374 if (pmax.y() <= ymin) continue; 374 if (pmax.y() <= ymin) continue;
375 if (pmax.z() <= zmin) continue; 375 if (pmax.z() <= zmin) continue;
376 for (G4int i = 0; i < res; ++i) 376 for (G4int i = 0; i < res; ++i)
377 { 377 {
378 G4ThreeVector p = points[i]; 378 G4ThreeVector p = points[i];
379 if (p.x() <= pmin.x()) continue; 379 if (p.x() <= pmin.x()) continue;
380 if (p.x() >= pmax.x()) continue; 380 if (p.x() >= pmax.x()) continue;
381 if (p.y() <= pmin.y()) continue; 381 if (p.y() <= pmin.y()) continue;
382 if (p.y() >= pmax.y()) continue; 382 if (p.y() >= pmax.y()) continue;
383 if (p.z() <= pmin.z()) continue; 383 if (p.z() <= pmin.z()) continue;
384 if (p.z() >= pmax.z()) continue; 384 if (p.z() >= pmax.z()) continue;
385 G4ThreeVector md = p - offset; 385 G4ThreeVector md = p - offset;
386 if (daughterSolid->Inside(md) == kInsi 386 if (daughterSolid->Inside(md) == kInside)
387 { 387 {
388 check_encapsulation = false; 388 check_encapsulation = false;
389 G4double distout = daughterSolid->Di 389 G4double distout = daughterSolid->DistanceToOut(md);
390 if (distout < tol) continue; // too 390 if (distout < tol) continue; // too small overlap
391 ++overlapCount; 391 ++overlapCount;
392 if (distout <= overlapSize) continue 392 if (distout <= overlapSize) continue;
393 overlapSize = distout; 393 overlapSize = distout;
394 overlapPoint = md; 394 overlapPoint = md;
395 } 395 }
396 } 396 }
397 } 397 }
398 else // transformation with rotation 398 else // transformation with rotation
399 { 399 {
400 G4ThreeVector pmin(pmin_local), pmax(pma 400 G4ThreeVector pmin(pmin_local), pmax(pmax_local);
401 G4ThreeVector dcenter = Td.TransformPoin 401 G4ThreeVector dcenter = Td.TransformPoint(0.5*(pmin + pmax));
402 G4double dradius = 0.5*((pmax - pmin).ma 402 G4double dradius = 0.5*((pmax - pmin).mag());
403 if ((scenter - dcenter).mag2() >= (sradi 403 if ((scenter - dcenter).mag2() >= (sradius + dradius)*(sradius + dradius)) continue;
404 if (dcenter.x() - dradius >= xmax) conti 404 if (dcenter.x() - dradius >= xmax) continue;
405 if (dcenter.y() - dradius >= ymax) conti 405 if (dcenter.y() - dradius >= ymax) continue;
406 if (dcenter.z() - dradius >= zmax) conti 406 if (dcenter.z() - dradius >= zmax) continue;
407 if (dcenter.x() + dradius <= xmin) conti 407 if (dcenter.x() + dradius <= xmin) continue;
408 if (dcenter.y() + dradius <= ymin) conti 408 if (dcenter.y() + dradius <= ymin) continue;
409 if (dcenter.z() + dradius <= zmin) conti 409 if (dcenter.z() + dradius <= zmin) continue;
410 410
411 G4ThreeVector pbox[8] = { 411 G4ThreeVector pbox[8] = {
412 G4ThreeVector(pmin.x(), pmin.y(), pmin 412 G4ThreeVector(pmin.x(), pmin.y(), pmin.z()),
413 G4ThreeVector(pmax.x(), pmin.y(), pmin 413 G4ThreeVector(pmax.x(), pmin.y(), pmin.z()),
414 G4ThreeVector(pmin.x(), pmax.y(), pmin 414 G4ThreeVector(pmin.x(), pmax.y(), pmin.z()),
415 G4ThreeVector(pmax.x(), pmax.y(), pmin 415 G4ThreeVector(pmax.x(), pmax.y(), pmin.z()),
416 G4ThreeVector(pmin.x(), pmin.y(), pmax 416 G4ThreeVector(pmin.x(), pmin.y(), pmax.z()),
417 G4ThreeVector(pmax.x(), pmin.y(), pmax 417 G4ThreeVector(pmax.x(), pmin.y(), pmax.z()),
418 G4ThreeVector(pmin.x(), pmax.y(), pmax 418 G4ThreeVector(pmin.x(), pmax.y(), pmax.z()),
419 G4ThreeVector(pmax.x(), pmax.y(), pmax 419 G4ThreeVector(pmax.x(), pmax.y(), pmax.z())
420 }; 420 };
421 G4double dxmin = kInfinity, dymin = kI 421 G4double dxmin = kInfinity, dymin = kInfinity, dzmin = kInfinity;
422 G4double dxmax = -kInfinity, dymax = -kI 422 G4double dxmax = -kInfinity, dymax = -kInfinity, dzmax = -kInfinity;
423 for (const auto & i : pbox) << 423 for (G4int i = 0; i < 8; ++i)
424 { 424 {
425 G4ThreeVector p = Td.TransformPoint(i) << 425 G4ThreeVector p = Td.TransformPoint(pbox[i]);
426 dxmin = std::min(dxmin, p.x()); 426 dxmin = std::min(dxmin, p.x());
427 dymin = std::min(dymin, p.y()); 427 dymin = std::min(dymin, p.y());
428 dzmin = std::min(dzmin, p.z()); 428 dzmin = std::min(dzmin, p.z());
429 dxmax = std::max(dxmax, p.x()); 429 dxmax = std::max(dxmax, p.x());
430 dymax = std::max(dymax, p.y()); 430 dymax = std::max(dymax, p.y());
431 dzmax = std::max(dzmax, p.z()); 431 dzmax = std::max(dzmax, p.z());
432 } 432 }
433 if (dxmin >= xmax) continue; 433 if (dxmin >= xmax) continue;
434 if (dymin >= ymax) continue; 434 if (dymin >= ymax) continue;
435 if (dzmin >= zmax) continue; 435 if (dzmin >= zmax) continue;
436 if (dxmax <= xmin) continue; 436 if (dxmax <= xmin) continue;
437 if (dymax <= ymin) continue; 437 if (dymax <= ymin) continue;
438 if (dzmax <= zmin) continue; 438 if (dzmax <= zmin) continue;
439 for (G4int i = 0; i < res; ++i) 439 for (G4int i = 0; i < res; ++i)
440 { 440 {
441 G4ThreeVector p = points[i]; 441 G4ThreeVector p = points[i];
442 if (p.x() >= dxmax) continue; 442 if (p.x() >= dxmax) continue;
443 if (p.x() <= dxmin) continue; 443 if (p.x() <= dxmin) continue;
444 if (p.y() >= dymax) continue; 444 if (p.y() >= dymax) continue;
445 if (p.y() <= dymin) continue; 445 if (p.y() <= dymin) continue;
446 if (p.z() >= dzmax) continue; 446 if (p.z() >= dzmax) continue;
447 if (p.z() <= dzmin) continue; 447 if (p.z() <= dzmin) continue;
448 G4ThreeVector md = Td.InverseTransform 448 G4ThreeVector md = Td.InverseTransformPoint(p);
449 if (daughterSolid->Inside(md) == kInsi 449 if (daughterSolid->Inside(md) == kInside)
450 { 450 {
451 check_encapsulation = false; 451 check_encapsulation = false;
452 G4double distout = daughterSolid->Di 452 G4double distout = daughterSolid->DistanceToOut(md);
453 if (distout < tol) continue; // too 453 if (distout < tol) continue; // too small overlap
454 ++overlapCount; 454 ++overlapCount;
455 if (distout <= overlapSize) continue 455 if (distout <= overlapSize) continue;
456 overlapSize = distout; 456 overlapSize = distout;
457 overlapPoint = md; 457 overlapPoint = md;
458 } 458 }
459 } 459 }
460 } 460 }
461 461
462 // Print information on overlap 462 // Print information on overlap
463 // 463 //
464 if (overlapCount > 0) 464 if (overlapCount > 0)
465 { 465 {
466 ++trials; 466 ++trials;
467 retval = true; 467 retval = true;
468 std::ostringstream message; 468 std::ostringstream message;
469 message << "Overlap with volume already 469 message << "Overlap with volume already placed !" << G4endl
470 << " Overlap is detecte 470 << " Overlap is detected for volume "
471 << GetName() << ':' << GetCopyNo 471 << GetName() << ':' << GetCopyNo()
472 << " (" << solid->GetEntityType( 472 << " (" << solid->GetEntityType() << ") with "
473 << daughter->GetName() << ':' << 473 << daughter->GetName() << ':' << daughter->GetCopyNo()
474 << " (" << daughterSolid->GetEnt 474 << " (" << daughterSolid->GetEntityType() << ")" << G4endl
475 << " overlap at local p 475 << " overlap at local point " << overlapPoint
476 << " by " << G4BestUnit(overlapS 476 << " by " << G4BestUnit(overlapSize, "Length")
477 << " (max of " << overlapCount < 477 << " (max of " << overlapCount << " cases)";
478 if (trials >= maxErr) 478 if (trials >= maxErr)
479 { 479 {
480 message << G4endl 480 message << G4endl
481 << "NOTE: Reached maximum fixe 481 << "NOTE: Reached maximum fixed number -" << maxErr
482 << "- of overlaps reports for 482 << "- of overlaps reports for this volume !";
483 } 483 }
484 G4Exception("G4PVPlacement::CheckOverlap 484 G4Exception("G4PVPlacement::CheckOverlaps()",
485 "GeomVol1002", JustWarning, 485 "GeomVol1002", JustWarning, message);
486 if (trials >= maxErr) { return true; } 486 if (trials >= maxErr) { return true; }
487 } 487 }
488 else if (check_encapsulation) 488 else if (check_encapsulation)
489 { 489 {
490 // Now checking that 'sister' volume is 490 // Now checking that 'sister' volume is not totally included
491 // and overlapping. Generate a single po 491 // and overlapping. Generate a single point inside of
492 // the 'sister' volume and verify that t 492 // the 'sister' volume and verify that the point in NOT inside
493 // the current volume 493 // the current volume
494 // 494 //
495 G4ThreeVector pSurface = daughterSolid-> 495 G4ThreeVector pSurface = daughterSolid->GetPointOnSurface();
496 G4ThreeVector normal = daughterSolid->Su 496 G4ThreeVector normal = daughterSolid->SurfaceNormal(pSurface);
497 G4ThreeVector pInside = pSurface - norma 497 G4ThreeVector pInside = pSurface - normal*1.e-4; // move point to inside
498 G4ThreeVector dPoint = (daughterSolid->I 498 G4ThreeVector dPoint = (daughterSolid->Inside(pInside) == kInside) ?
499 pInside : pSurface; 499 pInside : pSurface;
500 500
501 // Transform the generated point to the 501 // Transform the generated point to the mother's coordinate system
502 // and then to current volume's coordina 502 // and then to current volume's coordinate system
503 // 503 //
504 G4ThreeVector mp2 = Td.TransformPoint(dP 504 G4ThreeVector mp2 = Td.TransformPoint(dPoint);
505 G4ThreeVector msi = Tm.InverseTransformP 505 G4ThreeVector msi = Tm.InverseTransformPoint(mp2);
506 506
507 if (solid->Inside(msi) == kInside) 507 if (solid->Inside(msi) == kInside)
508 { 508 {
509 ++trials; 509 ++trials;
510 retval = true; 510 retval = true;
511 std::ostringstream message; 511 std::ostringstream message;
512 message << "Overlap with volume alread 512 message << "Overlap with volume already placed !" << G4endl
513 << " Overlap is detec 513 << " Overlap is detected for volume "
514 << GetName() << ':' << GetCopy 514 << GetName() << ':' << GetCopyNo()
515 << " (" << solid->GetEntityTyp 515 << " (" << solid->GetEntityType() << ")" << G4endl
516 << " apparently fully 516 << " apparently fully encapsulating volume "
517 << daughter->GetName() << ':' 517 << daughter->GetName() << ':' << daughter->GetCopyNo()
518 << " (" << daughterSolid->GetE 518 << " (" << daughterSolid->GetEntityType() << ")"
519 << " at the same level!"; 519 << " at the same level!";
520 if (trials >= maxErr) 520 if (trials >= maxErr)
521 { 521 {
522 message << G4endl 522 message << G4endl
523 << "NOTE: Reached maximum fi 523 << "NOTE: Reached maximum fixed number -" << maxErr
524 << "- of overlaps reports fo 524 << "- of overlaps reports for this volume !";
525 } 525 }
526 G4Exception("G4PVPlacement::CheckOverl 526 G4Exception("G4PVPlacement::CheckOverlaps()",
527 "GeomVol1002", JustWarning 527 "GeomVol1002", JustWarning, message);
528 if (trials >= maxErr) { return true; 528 if (trials >= maxErr) { return true; }
529 } 529 }
530 } 530 }
531 } 531 }
532 532
533 if (verbose && trials == 0) { G4cout << "OK! 533 if (verbose && trials == 0) { G4cout << "OK! " << G4endl; }
534 return retval; 534 return retval;
535 } 535 }
536 536
537 // ------------------------------------------- 537 // ----------------------------------------------------------------------
538 // NewPtrRotMatrix 538 // NewPtrRotMatrix
539 // 539 //
540 // Auxiliary function for 2nd & 4th constructo 540 // Auxiliary function for 2nd & 4th constructors (those with G4Transform3D)
541 // Creates a new rotation matrix on the heap ( 541 // Creates a new rotation matrix on the heap (using "new") and copies its
542 // argument into it. 542 // argument into it.
543 // 543 //
544 // NOTE: Ownership of the returned pointer is 544 // NOTE: Ownership of the returned pointer is left to the caller !
545 // No entity is currently responsible to 545 // No entity is currently responsible to delete this memory.
546 // 546 //
547 G4RotationMatrix* 547 G4RotationMatrix*
548 G4PVPlacement::NewPtrRotMatrix(const G4Rotatio 548 G4PVPlacement::NewPtrRotMatrix(const G4RotationMatrix &RotMat)
549 { 549 {
550 G4RotationMatrix* pRotMatrix; 550 G4RotationMatrix* pRotMatrix;
551 if ( RotMat.isIdentity() ) 551 if ( RotMat.isIdentity() )
552 { 552 {
553 pRotMatrix = nullptr; 553 pRotMatrix = nullptr;
554 } 554 }
555 else 555 else
556 { 556 {
557 pRotMatrix = new G4RotationMatrix(RotMat) 557 pRotMatrix = new G4RotationMatrix(RotMat);
558 } 558 }
559 return pRotMatrix; 559 return pRotMatrix;
560 } 560 }
561 561