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