| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // class G4ParameterisedNavigation Implementat
27 //
28 // Initial Author: P.Kent, 1996
29 // Revisions:
30 // J. Apostolakis 24 Nov 2005, Revised/fixed
31 // J. Apostolakis 4 Feb 2005, Reintroducting
32 // for materials
33 // G. Cosmo 11 Mar 2004, Added Check mo
34 // G. Cosmo 15 May 2002, Extended to 3-
35 // J. Apostolakis 5 Mar 1998, Enabled parame
36 // -------------------------------------------
37
38 // Note 1: Design/implementation note for exte
39 // We cannot make the solid, dimensions and tr
40 // parent because the voxelisation will not ha
41 // So the following can NOT be done:
42 // sampleSolid = curParam->ComputeSolid(num,
43 // sampleSolid->ComputeDimensions(curParam,
44 // curParam->ComputeTransformation(num, curP
45
46 #include "G4ParameterisedNavigation.hh"
47 #include "G4TouchableHistory.hh"
48 #include "G4VNestedParameterisation.hh"
49
50 #include "G4AuxiliaryNavServices.hh"
51
52 #include <cassert>
53
54 // *******************************************
55 // Constructor
56 // *******************************************
57 //
58 G4ParameterisedNavigation::G4ParameterisedNavi
59
60 // *******************************************
61 // Destructor
62 // *******************************************
63 //
64 G4ParameterisedNavigation::~G4ParameterisedNav
65
66 // *******************************************
67 // ComputeStep
68 // *******************************************
69 //
70 G4double G4ParameterisedNavigation::
71 ComputeStep(const G4ThreeV
72 const G4ThreeV
73 const G4double
74 G4double
75 G4Naviga
76 G4bool&
77 G4ThreeV
78 G4bool&
79 G4bool&
80 G4VPhysi
81 G4int& b
82 {
83 G4VPhysicalVolume *motherPhysical, *samplePh
84 G4VPVParameterisation *sampleParam;
85 G4LogicalVolume *motherLogical;
86 G4VSolid *motherSolid, *sampleSolid;
87 G4ThreeVector sampleDirection;
88 G4double ourStep=currentProposedStepLength,
89 G4double motherSafety, motherStep = DBL_MAX;
90 G4bool motherValidExitNormal = false;
91 G4ThreeVector motherExitNormal;
92
93 G4int sampleNo;
94
95 G4bool initialNode, noStep;
96 G4SmartVoxelNode *curVoxelNode;
97 G4long curNoVolumes, contentNo;
98 G4double voxelSafety;
99
100 // Replication data
101 //
102 EAxis axis;
103 G4int nReplicas;
104 G4double width, offset;
105 G4bool consuming;
106
107 motherPhysical = history.GetTopVolume();
108 motherLogical = motherPhysical->GetLogicalVo
109 motherSolid = motherLogical->GetSolid();
110
111 //
112 // Compute mother safety
113 //
114
115 motherSafety = motherSolid->DistanceToOut(lo
116 ourSafety = motherSafety; // Wo
117
118 #ifdef G4VERBOSE
119 if ( fCheck )
120 {
121 if( motherSafety < 0.0 )
122 {
123 motherSolid->DumpInfo();
124 std::ostringstream message;
125 message << "Negative Safety In Voxel Nav
126 << " Current solid " << m
127 << " gave negative safety: " <<
128 << " for the current (loc
129 G4Exception("G4ParameterisedNavigation::
130 "GeomNav0003", FatalExceptio
131 }
132 if( motherSolid->Inside(localPoint) == kOu
133 {
134 std::ostringstream message;
135 message << "Point is outside Current Vol
136 << " Point " << localPo
137 << " is outside current volume "
138 << G4endl;
139 G4double estDistToSolid = motherSolid->D
140 G4cout << " Estimated isotropic
141 << estDistToSolid;
142 if( estDistToSolid > 100.0 * motherSolid
143 {
144 motherSolid->DumpInfo();
145 G4Exception("G4ParameterisedNavigation
146 "GeomNav0003", FatalExcept
147 "Point is far outside Curr
148 }
149 else
150 {
151 G4Exception("G4ParameterisedNavigation
152 "GeomNav1002", JustWarning
153 "Point is a little outside
154 }
155 }
156
157 // Compute early:
158 // a) to check whether point is (wrongly)
159 // (signaled if step < 0 or
160 // b) to check value against answer of da
161 //
162 motherStep = motherSolid->DistanceToOut(lo
163 lo
164 tr
165 &mo
166 &mo
167
168 if( (motherStep >= kInfinity) || (motherSt
169 {
170 // Error - indication of being outside s
171 //
172 fLogger->ReportOutsideMother(localPoint,
173
174 ourStep = motherStep = 0.0;
175 exiting = true;
176 entering = false;
177
178 // If we are outside the solid does the
179 validExitNormal = motherValidExitNormal;
180 exitNormal = motherExitNormal;
181
182 *pBlockedPhysical = nullptr; // or mothe
183 blockedReplicaNo = 0; // or motherRepli
184
185 newSafety = 0.0;
186 return ourStep;
187 }
188 }
189 #endif
190
191 initialNode = true;
192 noStep = true;
193
194 // By definition, the parameterised volume i
195 // (and only) daughter of the mother volume
196 //
197 samplePhysical = motherLogical->GetDaughter(
198 samplePhysical->GetReplicationData(axis,nRep
199 fBList.Enlarge(nReplicas);
200 fBList.Reset();
201
202 // Exiting normal optimisation
203 //
204 if (exiting && (*pBlockedPhysical==samplePhy
205 {
206 if (localDirection.dot(exitNormal)>=kMinEx
207 {
208 // Block exited daughter replica; Must b
209 //
210 fBList.BlockVolume(blockedReplicaNo);
211 ourSafety = 0;
212 }
213 }
214 exiting = false;
215 entering = false;
216
217 sampleParam = samplePhysical->GetParameteris
218
219 // Loop over voxels & compute daughter safet
220
221 do
222 {
223 curVoxelNode = fVoxelNode;
224 curNoVolumes = curVoxelNode->GetNoContaine
225
226 for ( contentNo=curNoVolumes-1; contentNo>
227 {
228 sampleNo = curVoxelNode->GetVolume((G4in
229 if ( !fBList.IsBlocked(sampleNo) )
230 {
231 fBList.BlockVolume(sampleNo);
232
233 // Call virtual methods, and copy info
234 //
235 sampleSolid = IdentifyAndPlaceSolid( s
236 s
237
238 G4AffineTransform sampleTf(samplePhysi
239 samplePhysi
240 sampleTf.Invert();
241 const G4ThreeVector samplePoint = samp
242 const G4double sampleSafety = sampleSo
243 if ( sampleSafety<ourSafety )
244 {
245 ourSafety = sampleSafety;
246 }
247 if ( sampleSafety<=ourStep )
248 {
249 sampleDirection = sampleTf.Transform
250 G4double sampleStep =
251 sampleSolid->DistanceToIn(s
252 if ( sampleStep<=ourStep )
253 {
254 ourStep = sampleStep;
255 entering = true;
256 exiting = false;
257 *pBlockedPhysical = samplePhysical
258 blockedReplicaNo = sampleNo;
259 #ifdef G4VERBOSE
260 // Check to see that the resulti
261 // This check could eventually b
262 // candidate.
263
264 if ( ( fCheck ) && ( sampleStep
265 {
266 G4ThreeVector intersectionPoin
267 intersectionPoint = samplePoin
268 EInside insideIntPt = sampleSo
269 if( insideIntPt != kSurface )
270 {
271 G4long oldcoutPrec = G4cout.
272 std::ostringstream message;
273 message << "Navigator gets c
274 << G4endl
275 << " Inaccu
276 << " for solid " <<
277 << " Solid
278 << sampleStep << " y
279 if( insideIntPt == kInside )
280 {
281 message << "-kInside-";
282 }
283 else if( insideIntPt == kOut
284 {
285 message << "-kOutside-";
286 }
287 else
288 {
289 message << "-kSurface-";
290 }
291 message << " for this point
292 << " Point
293 << G4endl;
294 if ( insideIntPt != kInside
295 {
296 message << " Distan
297 << sampleSolid->Di
298 }
299 if ( insideIntPt != kOutside
300 {
301 message << " Distan
302 << sampleSolid->Di
303 }
304 G4Exception("G4Parameterised
305 "GeomNav1002", J
306 G4cout.precision(oldcoutPrec
307 }
308 }
309 #endif
310 }
311 }
312 }
313 }
314
315 if ( initialNode )
316 {
317 initialNode = false;
318 voxelSafety = ComputeVoxelSafety(localPo
319 if ( voxelSafety<ourSafety )
320 {
321 ourSafety = voxelSafety;
322 }
323 if ( currentProposedStepLength<ourSafety
324 {
325 // Guaranteed physics limited
326 //
327 noStep = false;
328 entering = false;
329 exiting = false;
330 *pBlockedPhysical = nullptr;
331 ourStep = kInfinity;
332 }
333 else
334 {
335 // Consider intersection with mother s
336 //
337 if ( motherSafety<=ourStep )
338 {
339 if ( !fCheck )
340 {
341 motherStep = motherSolid->Distance
342
343
344
345
346 }
347
348 if( ( motherStep < 0.0 ) || ( mother
349 {
350 #ifdef G4VERBOSE
351 fLogger->ReportOutsideMother(local
352 mothe
353 #endif
354 ourStep = motherStep = 0.0;
355 // Rely on the code below to set t
356 // exiting, entering, exitNormal
357 // pBlockedPhysical etc.
358 }
359 #ifdef G4VERBOSE
360 if( motherValidExitNormal && ( fChec
361 {
362 fLogger->CheckAndReportBadNormal(m
363 l
364 m
365 "
366 }
367 #endif
368 if ( motherStep<=ourStep )
369 {
370 ourStep = motherStep;
371 exiting = true;
372 entering = false;
373 if ( validExitNormal )
374 {
375 const G4RotationMatrix* rot = mo
376 if (rot != nullptr)
377 {
378 exitNormal *= rot->inverse();
379 }
380 }
381 }
382 else
383 {
384 validExitNormal = false;
385 }
386 }
387 }
388 newSafety = ourSafety;
389 }
390 if (noStep)
391 {
392 noStep = LocateNextVoxel(localPoint, loc
393 }
394 } while (noStep);
395
396 return ourStep;
397 }
398
399 // *******************************************
400 // ComputeSafety
401 // *******************************************
402 //
403 G4double
404 G4ParameterisedNavigation::ComputeSafety(const
405 const
406 const
407 {
408 G4VPhysicalVolume *motherPhysical, *samplePh
409 G4VPVParameterisation *sampleParam;
410 G4LogicalVolume *motherLogical;
411 G4VSolid *motherSolid, *sampleSolid;
412 G4double motherSafety, ourSafety;
413 G4int sampleNo, curVoxelNodeNo;
414
415 G4SmartVoxelNode *curVoxelNode;
416 G4long curNoVolumes, contentNo;
417 G4double voxelSafety;
418
419 // Replication data
420 //
421 EAxis axis;
422 G4int nReplicas;
423 G4double width, offset;
424 G4bool consuming;
425
426 motherPhysical = history.GetTopVolume();
427 motherLogical = motherPhysical->GetLogicalVo
428 motherSolid = motherLogical->GetSolid();
429
430 //
431 // Compute mother safety
432 //
433
434 motherSafety = motherSolid->DistanceToOut(lo
435 ourSafety = motherSafety;
436
437 //
438 // Compute daughter safeties
439 //
440
441 // By definition, parameterised volumes exis
442 // daughter of the mother volume
443 //
444 samplePhysical = motherLogical->GetDaughter(
445 samplePhysical->GetReplicationData(axis, nRe
446 width, of
447 sampleParam = samplePhysical->GetParameteris
448
449 // Look inside the current Voxel only at the
450 //
451 if ( axis==kUndefined ) // 3D case: cur
452 { // fro
453 curVoxelNode = fVoxelNode;
454 }
455 else // 1D case: cur
456 {
457 curVoxelNodeNo = G4int((localPoint(fVoxelA
458 -fVoxelHeader->GetM
459 curVoxelNode = fVoxelHeader->GetSlice(curV
460 fVoxelNodeNo = curVoxelNodeNo;
461 fVoxelNode = curVoxelNode;
462 }
463 curNoVolumes = curVoxelNode->GetNoContained(
464
465 for ( contentNo=curNoVolumes-1; contentNo>=0
466 {
467 sampleNo = curVoxelNode->GetVolume((G4int)
468
469 // Call virtual methods, and copy informat
470 //
471 sampleSolid= IdentifyAndPlaceSolid( sample
472
473 G4AffineTransform sampleTf(samplePhysical-
474 samplePhysical-
475 sampleTf.Invert();
476 const G4ThreeVector samplePoint = sampleTf
477 G4double sampleSafety = sampleSolid->Dista
478 if ( sampleSafety<ourSafety )
479 {
480 ourSafety = sampleSafety;
481 }
482 }
483
484 voxelSafety = ComputeVoxelSafety(localPoint,
485 if ( voxelSafety<ourSafety )
486 {
487 ourSafety=voxelSafety;
488 }
489
490 return ourSafety;
491 }
492
493 // *******************************************
494 // ComputeVoxelSafety
495 //
496 // Computes safety from specified point to col
497 // using already located point.
498 // *******************************************
499 //
500 G4double G4ParameterisedNavigation::
501 ComputeVoxelSafety(const G4ThreeVector& localP
502 const EAxis pAxis) const
503 {
504 // If no best axis is specified, adopt defau
505 // strategy as for placements
506 //
507 if ( pAxis==kUndefined )
508 {
509 return G4VoxelNavigation::ComputeVoxelSafe
510 }
511
512 G4double voxelSafety, plusVoxelSafety, minus
513 G4double curNodeOffset, minCurCommonDelta, m
514 G4long minCurNodeNoDelta, maxCurNodeNoDelta;
515
516 // Compute linear intersection distance to b
517 // to collected nodes at current level
518 //
519 curNodeOffset = fVoxelNodeNo*fVoxelSliceWidt
520 minCurCommonDelta = localPoint(fVoxelAxis)
521 - fVoxelHeader->GetMinExte
522 maxCurNodeNoDelta = fVoxelNode->GetMaxEquiva
523 minCurNodeNoDelta = fVoxelNodeNo-fVoxelNode-
524 maxCurCommonDelta = fVoxelSliceWidth-minCurC
525 plusVoxelSafety = minCurNodeNoDelta*fVoxel
526 minusVoxelSafety = maxCurNodeNoDelta*fVoxel
527 voxelSafety = std::min(plusVoxelSafety,minus
528
529 if ( voxelSafety<0 )
530 {
531 voxelSafety = 0;
532 }
533
534 return voxelSafety;
535 }
536
537 // *******************************************
538 // LocateNextVoxel
539 //
540 // Finds the next voxel from the current voxel
541 // in the specified direction.
542 //
543 // Returns false if all voxels considered
544 // true otherwise
545 // [current Step ends inside same voxel or lea
546 // *******************************************
547 //
548 G4bool G4ParameterisedNavigation::
549 LocateNextVoxel( const G4ThreeVector& localPoi
550 const G4ThreeVector& localDir
551 const G4double currentStep,
552 const EAxis pAxis)
553 {
554 // If no best axis is specified, adopt defau
555 // location strategy as for placements
556 //
557 if ( pAxis==kUndefined )
558 {
559 return G4VoxelNavigation::LocateNextVoxel(
560
561
562 }
563
564 G4bool isNewVoxel;
565 G4int newNodeNo;
566 G4double minVal, maxVal, curMinExtent, curCo
567
568 curMinExtent = fVoxelHeader->GetMinExtent();
569 curCoord = localPoint(fVoxelAxis)+currentSte
570 minVal = curMinExtent+fVoxelNode->GetMinEqui
571 isNewVoxel = false;
572
573 if ( minVal<=curCoord )
574 {
575 maxVal = curMinExtent
576 + (fVoxelNode->GetMaxEquivalentSlic
577 if ( maxVal<curCoord )
578 {
579 newNodeNo = fVoxelNode->GetMaxEquivalent
580 if ( newNodeNo<G4int(fVoxelHeader->GetNo
581 {
582 fVoxelNodeNo = newNodeNo;
583 fVoxelNode = fVoxelHeader->GetSlice(ne
584 isNewVoxel = true;
585 }
586 }
587 }
588 else
589 {
590 newNodeNo = fVoxelNode->GetMinEquivalentSl
591
592 // Must locate from newNodeNo no and down
593 // Repeat or earlier code...
594 //
595 if ( newNodeNo>=0 )
596 {
597 fVoxelNodeNo = newNodeNo;
598 fVoxelNode = fVoxelHeader->GetSlice(newN
599 isNewVoxel = true;
600 }
601 }
602 return isNewVoxel;
603 }
604
605 // *******************************************
606 // LevelLocate
607 // *******************************************
608 //
609 G4bool
610 G4ParameterisedNavigation::LevelLocate( G4Navi
611 const G4VPhy
612 const G4int
613 const G4Thre
614 const G4Thre
615 const G4bool
616 G4Thre
617 {
618 G4SmartVoxelHeader *motherVoxelHeader;
619 G4SmartVoxelNode *motherVoxelNode;
620 G4VPhysicalVolume *motherPhysical, *pPhysica
621 G4VPVParameterisation *pParam;
622 G4LogicalVolume *motherLogical;
623 G4VSolid *pSolid;
624 G4ThreeVector samplePoint;
625 G4int voxelNoDaughters, replicaNo;
626
627 motherPhysical = history.GetTopVolume();
628 motherLogical = motherPhysical->GetLogicalVo
629 motherVoxelHeader = motherLogical->GetVoxelH
630
631 // Find the voxel containing the point
632 //
633 motherVoxelNode = ParamVoxelLocate(motherVox
634
635 voxelNoDaughters = (G4int)motherVoxelNode->G
636 if ( voxelNoDaughters==0 ) { return false;
637
638 pPhysical = motherLogical->GetDaughter(0);
639 pParam = pPhysical->GetParameterisation();
640
641 // Save parent history in touchable history
642 // ... for use as parent t-h in ComputeMat
643 //
644 G4TouchableHistory parentTouchable( history
645
646 // Search replicated daughter volume
647 //
648 for ( auto sampleNo=voxelNoDaughters-1; samp
649 {
650 replicaNo = motherVoxelNode->GetVolume(sam
651 if ( (replicaNo!=blockedNum) || (pPhysical
652 {
653 // Obtain solid (as it can vary) and obt
654 //
655 pSolid = IdentifyAndPlaceSolid( replicaN
656
657 // Setup history
658 //
659 history.NewLevel(pPhysical, kParameteris
660 samplePoint = history.GetTopTransform().
661 if ( !G4AuxiliaryNavServices::CheckPoint
662 globalDirection, history.GetTopTra
663 {
664 history.BackLevel();
665 }
666 else
667 {
668 // Enter this daughter
669 //
670 localPoint = samplePoint;
671
672 // Set the correct copy number in phys
673 //
674 pPhysical->SetCopyNo(replicaNo);
675
676 // Set the correct solid and material
677 //
678 G4LogicalVolume *pLogical = pPhysical-
679 pLogical->SetSolid(pSolid);
680 pLogical->UpdateMaterial(pParam->Compu
681 pPhysical, &p
682 return true;
683 }
684 }
685 }
686 return false;
687 }
688
689 void G4ParameterisedNavigation::RelocateWithin
690
691 {
692 auto motherLogical = motherPhysical->GetLogi
693
694 /* this should only be called on parameteriz
695 assert(motherPhysical->GetRegularStructureId
696 assert(motherLogical->GetNoDaughters() == 1)
697
698 if ( auto pVoxelHeader = motherLogical->GetV
699 ParamVoxelLocate( pVoxelHeader, localPoint
700 }
701