Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/navigation/src/G4VIntersectionLocator.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 G4VIntersectionLocator implementation
 27 //
 28 // 27.10.08 - John Apostolakis, Tatiana Nikitina.
 29 // ---------------------------------------------------------------------------
 30  
 31 #include <iomanip>
 32 #include <sstream>
 33 
 34 #include "globals.hh"
 35 #include "G4ios.hh"
 36 #include "G4AutoDelete.hh"
 37 #include "G4SystemOfUnits.hh"
 38 #include "G4VIntersectionLocator.hh"
 39 #include "G4TouchableHandle.hh"
 40 #include "G4GeometryTolerance.hh"
 41 
 42 ///////////////////////////////////////////////////////////////////////////
 43 //
 44 // Constructor
 45 //
 46 G4VIntersectionLocator::G4VIntersectionLocator(G4Navigator* theNavigator)
 47  : fiNavigator(theNavigator)
 48 {
 49   kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
 50 
 51   if( fiNavigator->GetExternalNavigation() == nullptr )
 52   {
 53     fHelpingNavigator = new G4Navigator();
 54   }
 55   else // Must clone the navigator, together with External Navigation
 56   {
 57     fHelpingNavigator = fiNavigator->Clone();
 58   }
 59 }      
 60 
 61 ///////////////////////////////////////////////////////////////////////////
 62 //
 63 // Destructor.
 64 //
 65 G4VIntersectionLocator::~G4VIntersectionLocator()
 66 {
 67   delete fHelpingNavigator;
 68   delete fpTouchable;
 69 }
 70 
 71 ///////////////////////////////////////////////////////////////////////////
 72 //
 73 // Dump status of propagator to cout (old method)
 74 //
 75 void
 76 G4VIntersectionLocator::printStatus( const G4FieldTrack& StartFT,
 77                                      const G4FieldTrack& CurrentFT, 
 78                                            G4double      requestStep, 
 79                                            G4double      safety,
 80                                            G4int         stepNo)
 81 {
 82   std::ostringstream os; 
 83   printStatus( StartFT,CurrentFT,requestStep,safety,stepNo,os,fVerboseLevel);
 84   G4cout << os.str();
 85 }
 86 
 87 ///////////////////////////////////////////////////////////////////////////
 88 //
 89 // Dumps status of propagator.
 90 //
 91 void
 92 G4VIntersectionLocator::printStatus( const G4FieldTrack& StartFT,
 93                                      const G4FieldTrack& CurrentFT, 
 94                                            G4double      requestStep, 
 95                                            G4double      safety,
 96                                            G4int         stepNo,
 97                                            std::ostream& os,
 98                                            G4int         verboseLevel)
 99 {
100   // const G4int verboseLevel= fVerboseLevel;
101   const G4ThreeVector StartPosition       = StartFT.GetPosition();
102   const G4ThreeVector StartUnitVelocity   = StartFT.GetMomentumDir();
103   const G4ThreeVector CurrentPosition     = CurrentFT.GetPosition();
104   const G4ThreeVector CurrentUnitVelocity = CurrentFT.GetMomentumDir();
105 
106   G4double step_len = CurrentFT.GetCurveLength() - StartFT.GetCurveLength();
107   G4long oldprc;  // cout/cerr precision settings
108 
109   if( ((stepNo == 0) && (verboseLevel <3)) || (verboseLevel >= 3) )
110   {
111     oldprc = os.precision(4);
112     os << std::setw( 6)  << " " 
113            << std::setw( 25) << " Current Position  and  Direction" << " "
114            << G4endl; 
115     os << std::setw( 5) << "Step#" 
116            << std::setw(10) << "  s  " << " "
117            << std::setw(10) << "X(mm)" << " "
118            << std::setw(10) << "Y(mm)" << " "  
119            << std::setw(10) << "Z(mm)" << " "
120            << std::setw( 7) << " N_x " << " "
121            << std::setw( 7) << " N_y " << " "
122            << std::setw( 7) << " N_z " << " " ;
123     os << std::setw( 7) << " Delta|N|" << " "
124            << std::setw( 9) << "StepLen" << " "  
125            << std::setw(12) << "StartSafety" << " "  
126            << std::setw( 9) << "PhsStep" << " ";  
127     os << G4endl;
128     os.precision(oldprc);
129   }
130   if((stepNo == 0) && (verboseLevel <=3))
131   {
132     // Recurse to print the start values
133     //
134     printStatus( StartFT, StartFT, -1.0, safety, -1, os, verboseLevel);
135   }
136   if( verboseLevel <= 3 )
137   {
138     if( stepNo >= 0)
139     {
140        os << std::setw( 4) << stepNo << " ";
141     }
142     else
143     {
144        os << std::setw( 5) << "Start" ;
145     }
146     oldprc = os.precision(8);
147     os << std::setw(10) << CurrentFT.GetCurveLength() << " "; 
148     os << std::setw(10) << CurrentPosition.x() << " "
149            << std::setw(10) << CurrentPosition.y() << " "
150            << std::setw(10) << CurrentPosition.z() << " ";
151     os.precision(4);
152     os << std::setw( 7) << CurrentUnitVelocity.x() << " "
153            << std::setw( 7) << CurrentUnitVelocity.y() << " "
154            << std::setw( 7) << CurrentUnitVelocity.z() << " ";
155     os.precision(3); 
156     os << std::setw( 7)
157            << CurrentFT.GetMomentum().mag()- StartFT.GetMomentum().mag()
158            << " "; 
159     os << std::setw( 9) << step_len << " "; 
160     os << std::setw(12) << safety << " ";
161     if( requestStep != -1.0 )
162     {
163       os << std::setw( 9) << requestStep << " ";
164     }
165     else
166     {
167       os << std::setw( 9) << "Init/NotKnown" << " "; 
168     }
169     os << G4endl;
170     os.precision(oldprc);
171   }
172   else // if( verboseLevel > 3 )
173   {
174     //  Multi-line output
175        
176     os << "Step taken was " << step_len  
177            << " out of PhysicalStep= " <<  requestStep << G4endl;
178     os << "Final safety is: " << safety << G4endl;
179     os << "Chord length = " << (CurrentPosition-StartPosition).mag()
180            << G4endl;
181     os << G4endl; 
182   }
183 }
184 
185 ///////////////////////////////////////////////////////////////////////////
186 //
187 // ReEstimateEndPoint.
188 //
189 G4FieldTrack G4VIntersectionLocator::
190 ReEstimateEndpoint( const G4FieldTrack& CurrentStateA,  
191                     const G4FieldTrack& EstimatedEndStateB,
192                           G4double      , // linearDistSq,  // NOT used
193                           G4double
194 #ifdef G4DEBUG_FIELD
195   curveDist
196 #endif
197                                    )
198 {  
199   G4FieldTrack newEndPoint( CurrentStateA );
200   auto integrDriver = GetChordFinderFor()->GetIntegrationDriver(); 
201 
202   G4FieldTrack retEndPoint( CurrentStateA );
203   G4bool goodAdvance;
204   G4int  itrial = 0;
205   const G4int no_trials = 20;
206 
207 
208   G4double endCurveLen= EstimatedEndStateB.GetCurveLength();
209 
210   do  // Loop checking, 07.10.2016, JA
211   {
212     G4double currentCurveLen = newEndPoint.GetCurveLength();
213     G4double advanceLength = endCurveLen - currentCurveLen ; 
214     if (std::abs(advanceLength)<kCarTolerance)
215     {
216       goodAdvance=true;
217     }
218     else
219     {
220       goodAdvance = integrDriver->AccurateAdvance(newEndPoint, advanceLength,
221                                                   GetEpsilonStepFor());
222     }
223   }
224   while( !goodAdvance && (++itrial < no_trials) );
225 
226   if( goodAdvance )
227   {
228     retEndPoint = newEndPoint; 
229   }
230   else
231   {
232     retEndPoint = EstimatedEndStateB; // Could not improve without major work !!
233   }
234 
235   //  All the work is done
236   //  below are some diagnostics only -- before the return!
237   // 
238   const G4String MethodName("G4VIntersectionLocator::ReEstimateEndpoint()");
239 
240 #ifdef G4VERBOSE
241   G4int  latest_good_trials = 0;
242   if( itrial > 1)
243   {
244     if( fVerboseLevel > 0 )
245     {
246       G4cout << MethodName << " called - goodAdv= " << goodAdvance
247              << " trials = " << itrial
248              << " previous good= " << latest_good_trials
249              << G4endl;
250     }
251     latest_good_trials = 0; 
252   }
253   else
254   {   
255     ++latest_good_trials; 
256   }
257 #endif
258 
259 #ifdef G4DEBUG_FIELD
260   G4double lengthDone = newEndPoint.GetCurveLength() 
261                       - CurrentStateA.GetCurveLength(); 
262   if( !goodAdvance )
263   {
264     if( fVerboseLevel >= 3 )
265     {
266       G4cout << MethodName << "> AccurateAdvance failed " ;
267       G4cout << " in " << itrial << " integration trials/steps. " << G4endl;
268       G4cout << " It went only " << lengthDone << " instead of " << curveDist
269              << " -- a difference of " << curveDist - lengthDone  << G4endl;
270       G4cout << " ReEstimateEndpoint> Reset endPoint to original value!"
271              << G4endl;
272     }
273   }
274   G4double linearDist = ( EstimatedEndStateB.GetPosition() 
275                       - CurrentStateA.GetPosition() ).mag(); 
276   static G4int noInaccuracyWarnings = 0; 
277   G4int maxNoWarnings = 10;
278   if (  (noInaccuracyWarnings < maxNoWarnings ) 
279        || (fVerboseLevel > 1) )
280     {
281       G4ThreeVector move = newEndPoint.GetPosition()
282                          - EstimatedEndStateB.GetPosition();
283       std::ostringstream message;
284       message.precision(12);
285       message << " Integration inaccuracy requires" 
286               << " an adjustment in the step's endpoint."  << G4endl
287               << "   Two mid-points are further apart than their"
288               << " curve length difference"                << G4endl 
289               << "   Dist = "       << linearDist
290               << " curve length = " << curveDist             << G4endl; 
291       message << " Correction applied is " << move.mag() << G4endl
292               << "  Old Estimated B position= "
293               << EstimatedEndStateB.GetPosition() << G4endl
294               << "  Recalculated    Position= "
295               << newEndPoint.GetPosition() << G4endl
296               << "   Change ( new - old )   = " << move;
297       G4Exception("G4VIntersectionLocator::ReEstimateEndpoint()",
298                   "GeomNav1002", JustWarning, message);
299     }
300 /*
301 #else
302   // Statistics on the RMS value of the corrections
303 
304   static G4ThreadLocal G4int noCorrections = 0;
305   ++noCorrections; 
306   if( goodAdvance )
307   {
308     static G4ThreadLocal G4double sumCorrectionsSq;
309     sumCorrectionsSq += (EstimatedEndStateB.GetPosition() - 
310                          newEndPoint.GetPosition()).mag2();
311   }
312 */
313 #endif
314 
315   return retEndPoint;
316 }
317 
318 
319 ///////////////////////////////////////////////////////////////////////////
320 //
321 // ReEstimateEndPoint.
322 //
323 //   The following values are returned in curveError 
324 //       0: Normal - no problem
325 //       1: Unexpected co-incidence - milder mixup
326 //       2: Real mixup - EndB is NOT past StartA  
327 //            ( ie. StartA's curve-lengh is bigger than EndB's)
328 
329 
330 G4bool G4VIntersectionLocator::
331 CheckAndReEstimateEndpoint( const G4FieldTrack& CurrentStartA,  
332                             const G4FieldTrack& EstimatedEndB,
333                                   G4FieldTrack& RevisedEndPoint,
334                                   G4int&        curveError)
335 {
336   G4double linDistSq, curveDist; 
337 
338   G4bool recalculated = false;
339   curveError= 0;
340 
341   linDistSq = ( EstimatedEndB.GetPosition() 
342               - CurrentStartA.GetPosition() ).mag2(); 
343   curveDist = EstimatedEndB.GetCurveLength()
344             - CurrentStartA.GetCurveLength();
345   if(  (curveDist>=0.0) 
346      && (curveDist*curveDist *(1.0+2.0*fiEpsilonStep ) < linDistSq ) )
347   {
348     G4FieldTrack newEndPointFT = EstimatedEndB;  // Unused
349 
350     if (curveDist>0.0) 
351     {
352        // Re-integrate to obtain a new B
353        RevisedEndPoint = ReEstimateEndpoint( CurrentStartA,
354                                              EstimatedEndB,
355                                              linDistSq,    
356                                              curveDist );
357        recalculated = true;
358     }
359     else
360     {
361        // Zero length -> no advance!
362        newEndPointFT = CurrentStartA;
363        recalculated = true;
364        curveError = 1;  // Unexpected co-incidence - milder mixup
365 
366        G4Exception("G4MultiLevelLocator::EstimateIntersectionPoint()",
367            "GeomNav1002", JustWarning, 
368            "A & B are at equal distance in 2nd half. A & B will coincide." );
369     }
370   }
371 
372   // Sanity check
373   //
374   if( curveDist < 0.0 )
375   {
376     curveError = 2;  //  Real mixup
377   }
378   return recalculated;
379 }
380 
381 ///////////////////////////////////////////////////////////////////////////
382 //
383 // Method for finding SurfaceNormal of Intersecting Solid 
384 //
385 G4ThreeVector G4VIntersectionLocator::
386 GetLocalSurfaceNormal(const G4ThreeVector& CurrentE_Point, G4bool& validNormal)
387 {
388   G4ThreeVector Normal(G4ThreeVector(0.0,0.0,0.0));
389   G4VPhysicalVolume* located;
390 
391   validNormal = false;
392   fHelpingNavigator->SetWorldVolume(GetNavigatorFor()->GetWorldVolume());
393   located = fHelpingNavigator->LocateGlobalPointAndSetup( CurrentE_Point );
394 
395   delete fpTouchable;
396   fpTouchable = fHelpingNavigator->CreateTouchableHistory();
397 
398   // To check if we can use GetGlobalExitNormal() 
399   //
400   G4ThreeVector localPosition = fpTouchable->GetHistory()
401                 ->GetTopTransform().TransformPoint(CurrentE_Point);
402 
403   // Issue: in the case of coincident surfaces, this version does not recognise 
404   //        which side you are located onto (can return vector with wrong sign.)
405   // TO-DO: use direction (of chord) to identify volume we will be "entering"
406 
407   if( located != nullptr)
408   { 
409     G4LogicalVolume* pLogical= located->GetLogicalVolume(); 
410     G4VSolid*        pSolid; 
411 
412     if( (pLogical != nullptr) && ( (pSolid=pLogical->GetSolid()) != nullptr ) )
413     {
414       if ( ( pSolid->Inside(localPosition)==kSurface )
415         || ( pSolid->DistanceToOut(localPosition) < 1000.0 * kCarTolerance ) )
416       {
417         Normal = pSolid->SurfaceNormal(localPosition);
418         validNormal = true;
419 
420 #ifdef G4DEBUG_FIELD
421         if( std::fabs(Normal.mag2() - 1.0 ) > CLHEP::perThousand) 
422         {
423           G4cerr << "PROBLEM in G4VIntersectionLocator::GetLocalSurfaceNormal."
424                  << G4endl;
425           G4cerr << "  Normal is not unit - mag=" << Normal.mag() << G4endl; 
426           G4cerr << "  at trial local point " << CurrentE_Point << G4endl;
427           G4cerr <<  "  Solid is " << *pSolid << G4endl;
428         }
429 #endif
430       }
431     }
432   }
433   return Normal;
434 }
435 
436 ///////////////////////////////////////////////////////////////////////////
437 //
438 // Adjustment of Found Intersection
439 //
440 G4bool G4VIntersectionLocator::
441 AdjustmentOfFoundIntersection( const G4ThreeVector& CurrentA_Point,
442                                const G4ThreeVector& CurrentE_Point, 
443                                const G4ThreeVector& CurrentF_Point,
444                                const G4ThreeVector& MomentumDir,
445                                const G4bool         IntersectAF,
446                                      G4ThreeVector& IntersectionPoint,  // I/O
447                                      G4double&      NewSafety,          // I/O 
448                                      G4double&      fPreviousSafety,    // I/O
449                                      G4ThreeVector& fPreviousSftOrigin )// I/O
450 {     
451   G4double dist,lambda;
452   G4ThreeVector Normal, NewPoint, Point_G;
453   G4bool goodAdjust = false, Intersects_FP = false, validNormal = false;
454 
455   // Get SurfaceNormal of Intersecting Solid
456   //
457   Normal = GetGlobalSurfaceNormal(CurrentE_Point,validNormal);
458   if(!validNormal) { return false; }
459 
460   // Intersection between Line and Plane
461   //
462   G4double n_d_m = Normal.dot(MomentumDir);
463   if ( std::abs(n_d_m)>kCarTolerance )
464   {
465 #ifdef G4VERBOSE
466     if ( fVerboseLevel>1 )
467     {
468       G4Exception("G4VIntersectionLocator::AdjustmentOfFoundIntersection()",
469                   "GeomNav0003", JustWarning,
470                   "No intersection. Parallels lines!");
471     }
472 #endif
473     lambda =- Normal.dot(CurrentF_Point-CurrentE_Point)/n_d_m;
474 
475     // New candidate for Intersection
476     //
477     NewPoint = CurrentF_Point+lambda*MomentumDir;
478 
479     // Distance from CurrentF to Calculated Intersection
480     //
481     dist = std::abs(lambda);
482 
483     if ( dist<kCarTolerance*0.001 )  { return false; }
484 
485     // Calculation of new intersection point on the path.
486     //
487     if ( IntersectAF )  //  First part intersects
488     {
489       G4double stepLengthFP; 
490       G4ThreeVector Point_P = CurrentA_Point;
491       GetNavigatorFor()->LocateGlobalPointWithinVolume(Point_P);
492       Intersects_FP = IntersectChord( Point_P, NewPoint, NewSafety,
493                                       fPreviousSafety, fPreviousSftOrigin,
494                                       stepLengthFP, Point_G );
495 
496     }
497     else   // Second part intersects
498     {      
499       G4double stepLengthFP; 
500       GetNavigatorFor()->LocateGlobalPointWithinVolume(CurrentF_Point );
501       Intersects_FP = IntersectChord( CurrentF_Point, NewPoint, NewSafety,
502                                       fPreviousSafety, fPreviousSftOrigin,
503                                       stepLengthFP, Point_G );
504     }
505     if ( Intersects_FP )
506     {
507       goodAdjust = true;
508       IntersectionPoint = Point_G;              
509     }
510   }
511 
512   return goodAdjust;
513 }
514 
515 ///////////////////////////////////////////////////////////////////////////
516 //
517 // GetSurfaceNormal.
518 //
519 G4ThreeVector G4VIntersectionLocator::
520 GetSurfaceNormal(const G4ThreeVector& CurrentInt_Point,
521                        G4bool& validNormal)
522 {
523   G4ThreeVector NormalAtEntry; // ( -10. , -10., -10. ); 
524 
525   G4ThreeVector NormalAtEntryLast, NormalAtEntryGlobal, diffNormals;
526   G4bool validNormalLast; 
527 
528   // Relies on a call to Navigator::ComputeStep in IntersectChord before
529   // this call
530   //
531   NormalAtEntryLast = GetLastSurfaceNormal( CurrentInt_Point, validNormalLast );
532     // May return valid=false in cases, including
533     //  - if the candidate volume was not found (eg exiting world), or
534     //  - a replica was involved -- determined the step size.
535     // (This list is not complete.) 
536 
537 #ifdef G4DEBUG_FIELD
538   if  ( validNormalLast
539    && ( std::fabs(NormalAtEntryLast.mag2() - 1.0) > perThousand ) )
540   {
541     std::ostringstream message; 
542     message << "PROBLEM: Normal is not unit - magnitude = "
543             << NormalAtEntryLast.mag() << G4endl; 
544     message << "   at trial intersection point " << CurrentInt_Point << G4endl;
545     message << "   Obtained from Get *Last* Surface Normal."; 
546     G4Exception("G4VIntersectionLocator::GetSurfaceNormal()",
547                 "GeomNav1002", JustWarning, message);
548   }
549 #endif
550 
551   if( validNormalLast ) 
552   {
553     NormalAtEntry = NormalAtEntryLast;  
554   }
555   validNormal = validNormalLast;
556 
557   return NormalAtEntry;
558 }
559 
560 ///////////////////////////////////////////////////////////////////////////
561 //
562 // GetGlobalSurfaceNormal.
563 //
564 G4ThreeVector G4VIntersectionLocator::
565 GetGlobalSurfaceNormal(const G4ThreeVector& CurrentE_Point,
566                              G4bool& validNormal)
567 {
568   G4ThreeVector localNormal = GetLocalSurfaceNormal(CurrentE_Point,validNormal);
569   G4AffineTransform localToGlobal =          //  Must use the same Navigator !!
570       fHelpingNavigator->GetLocalToGlobalTransform();
571   G4ThreeVector globalNormal = localToGlobal.TransformAxis( localNormal );
572 
573 #ifdef G4DEBUG_FIELD
574   if( validNormal && ( std::fabs(globalNormal.mag2() - 1.0) > perThousand ) ) 
575   {
576     std::ostringstream message; 
577     message << "**************************************************************"
578             << G4endl;
579     message << " Bad Normal in G4VIntersectionLocator::GetGlobalSurfaceNormal "
580             << G4endl;
581     message << "  * Constituents: " << G4endl;
582     message << "    Local  Normal= " << localNormal << G4endl;
583     message << "    Transform: " << G4endl
584             << "      Net Translation= " << localToGlobal.NetTranslation()
585             << G4endl
586             << "      Net Rotation   = " << localToGlobal.NetRotation()
587             << G4endl;
588     message << "  * Result: " << G4endl;
589     message << "     Global Normal= " << localNormal << G4endl;
590     message << "**************************************************************";
591     G4Exception("G4VIntersectionLocator::GetGlobalSurfaceNormal()",
592                 "GeomNav1002", JustWarning, message);
593   }
594 #endif
595 
596   return globalNormal;
597 }
598 
599 ///////////////////////////////////////////////////////////////////////////
600 //
601 // GetLastSurfaceNormal.
602 //
603 G4ThreeVector G4VIntersectionLocator::
604 GetLastSurfaceNormal( const G4ThreeVector& intersectPoint,
605                             G4bool& normalIsValid) const
606 {
607   G4ThreeVector normalVec;
608   G4bool validNorm;
609   normalVec = fiNavigator->GetGlobalExitNormal( intersectPoint, &validNorm ); 
610   normalIsValid = validNorm;
611 
612   return normalVec;
613 }
614 
615 ///////////////////////////////////////////////////////////////////////////
616 //
617 // ReportTrialStep.
618 //
619 void G4VIntersectionLocator::ReportTrialStep( G4int step_no, 
620                                         const G4ThreeVector& ChordAB_v,
621                                         const G4ThreeVector& ChordEF_v,
622                                         const G4ThreeVector& NewMomentumDir,
623                                         const G4ThreeVector& NormalAtEntry,
624                                               G4bool validNormal )
625 {
626   G4double ABchord_length  = ChordAB_v.mag(); 
627   G4double MomDir_dot_Norm = NewMomentumDir.dot( NormalAtEntry );
628   G4double MomDir_dot_ABchord;
629   MomDir_dot_ABchord = (1.0 / ABchord_length) * NewMomentumDir.dot( ChordAB_v );
630 
631   std::ostringstream  outStream; 
632   outStream << std::setw(6)  << " Step# "
633     << std::setw(17) << " |ChordEF|(mag)" << "  "
634     << std::setw(18) << " uMomentum.Normal" << "  "
635     << std::setw(18) << " uMomentum.ABdir " << "  " 
636     << std::setw(16) << " AB-dist         " << " " 
637     << " Chord Vector (EF) " 
638     << G4endl;
639   outStream.precision(7); 
640   outStream << " " << std::setw(5) << step_no           
641     << " " << std::setw(18) << ChordEF_v.mag() 
642     << " " << std::setw(18) << MomDir_dot_Norm    
643     << " " << std::setw(18) << MomDir_dot_ABchord 
644     << " " << std::setw(12) << ABchord_length     
645     << " " << ChordEF_v
646     << G4endl;
647   outStream << " MomentumDir= " << " " << NewMomentumDir 
648     << " Normal at Entry E= " << NormalAtEntry
649     << " AB chord =   " << ChordAB_v
650     << G4endl;
651   G4cout << outStream.str();
652 
653   if( ( std::fabs(NormalAtEntry.mag2() - 1.0) > perThousand ) ) 
654   {
655     std::ostringstream message; 
656     message << "Normal is not unit - mag= " << NormalAtEntry.mag() << G4endl
657             << "         ValidNormalAtE = " << validNormal;
658     G4Exception("G4VIntersectionLocator::ReportTrialStep()",
659                 "GeomNav1002", JustWarning, message);
660   }
661   return; 
662 }
663 
664 ///////////////////////////////////////////////////////////////////////////
665 //
666 // LocateGlobalPointWithinVolumeAndCheck.
667 //
668 // Locate point using navigator: updates state of Navigator
669 // By default, it assumes that the point is inside the current volume,
670 // and returns true.
671 // In check mode, checks whether the point is *inside* the volume.
672 //   If it is inside, it returns true
673 //   If not, issues a warning and returns false.
674 //
675 G4bool G4VIntersectionLocator::
676 LocateGlobalPointWithinVolumeAndCheck( const G4ThreeVector& position )
677 {
678   G4bool good = true;
679   G4Navigator* nav = GetNavigatorFor();
680   const G4String
681   MethodName("G4VIntersectionLocator::LocateGlobalPointWithinVolumeAndCheck()");
682 
683   if( fCheckMode )
684   {
685     G4bool navCheck= nav->IsCheckModeActive();  // Recover original value
686     nav->CheckMode(true);
687 
688     // Identify the current volume
689     
690     G4TouchableHandle startTH= nav->CreateTouchableHistoryHandle();
691     G4VPhysicalVolume* motherPhys = startTH->GetVolume();
692     G4VSolid*          motherSolid = startTH->GetSolid();
693     G4AffineTransform transform = nav->GetGlobalToLocalTransform();
694     G4int motherCopyNo = motherPhys->GetCopyNo();
695     
696     // Let's check that the point is inside the current solid
697     G4ThreeVector  localPosition = transform.TransformPoint(position);
698     EInside        inMother = motherSolid->Inside( localPosition );
699     if( inMother != kInside )
700     {
701       std::ostringstream message; 
702       message << "Position located "
703               << ( inMother == kSurface ? " on Surface " : " outside " )
704               << "expected volume" << G4endl
705               << "  Safety (from Outside) = "
706               << motherSolid->DistanceToIn(localPosition);
707       G4Exception("G4VIntersectionLocator::LocateGlobalPointWithinVolumeAndCheck()",
708                   "GeomNav1002", JustWarning, message);
709     }
710     
711     // 1. Simple next step - quick relocation and check result.
712     // nav->LocateGlobalPointWithinVolume( position );
713     
714     // 2. Full relocation - to cross-check answer !
715     G4VPhysicalVolume* nextPhysical= nav->LocateGlobalPointAndSetup(position);
716     if(    (nextPhysical != motherPhys)
717         || (nextPhysical->GetCopyNo() != motherCopyNo )
718        )
719     {
720       G4Exception("G4VIntersectionLocator::LocateGlobalPointWithinVolumeAndCheck()",
721                   "GeomNav1002", JustWarning,
722                   "Position located outside expected volume.");
723     }
724     nav->CheckMode(navCheck);  // Recover original value
725   }
726   else
727   {
728     nav->LocateGlobalPointWithinVolume( position );
729   }
730   return good;
731 }
732 
733 ///////////////////////////////////////////////////////////////////////////
734 //
735 // LocateGlobalPointWithinVolumeCheckAndReport.
736 //
737 void G4VIntersectionLocator::
738 LocateGlobalPointWithinVolumeCheckAndReport( const G4ThreeVector& position,
739                                              const G4String& CodeLocationInfo,
740                                              G4int /* CheckMode */)
741 {
742   // Save value of Check mode first
743   G4bool oldCheck = GetCheckMode();
744   
745   G4bool ok = LocateGlobalPointWithinVolumeAndCheck( position );
746   if( !ok )
747   {
748     std::ostringstream message; 
749     message << "Failed point location." << G4endl
750             << "   Code Location info: " << CodeLocationInfo;
751     G4Exception("G4VIntersectionLocator::LocateGlobalPointWithinVolumeCheckAndReport()",
752                 "GeomNav1002", JustWarning, message);
753   }
754   
755   SetCheckMode( oldCheck );
756 }
757 
758 ///////////////////////////////////////////////////////////////////////////
759 //
760 // ReportReversedPoints.
761 //
762 void G4VIntersectionLocator::
763 ReportReversedPoints( std::ostringstream& msg,
764                       const G4FieldTrack& StartPointVel, 
765                       const G4FieldTrack& EndPointVel,
766                             G4double NewSafety, G4double epsStep,
767                       const G4FieldTrack& A_PtVel,
768                       const G4FieldTrack& B_PtVel,
769                       const G4FieldTrack& SubStart_PtVel,
770                       const G4ThreeVector& E_Point,
771                       const G4FieldTrack& ApproxIntersecPointV,
772                             G4int  substep_no, G4int substep_no_p, G4int depth )
773 {
774    // Expect that 'msg' can hold the name of the calling method
775 
776    // FieldTrack 'points' A and B have been tangled
777    // Whereas A should be before B, it is found that curveLen(B) < curveLen(A)
778    G4int verboseLevel= 5; 
779    G4double curveDist = B_PtVel.GetCurveLength() - A_PtVel.GetCurveLength();
780    G4VIntersectionLocator::printStatus( A_PtVel,  B_PtVel,
781                            -1.0, NewSafety,  substep_no, msg, verboseLevel );
782    msg << "Error in advancing propagation." << G4endl
783        << "   The final curve point is NOT further along"
784        << "  than the original!" << G4endl
785        << "   Going *backwards* from len(A) = " << A_PtVel.GetCurveLength()
786        << "  to len(B) = " << B_PtVel.GetCurveLength() << G4endl
787        << "      Curve distance is " << curveDist / CLHEP::millimeter << " mm "
788        << G4endl
789        << "      Point A' (start) is " << A_PtVel  << G4endl
790        << "      Point B' (end)   is " << B_PtVel << G4endl;
791    msg << "      fEpsStep= " << epsStep << G4endl << G4endl;
792 
793    G4long oldprc = msg.precision(20);
794    msg << " In full precision, the position, momentum, E_kin, length, rest mass "
795        << " ... are: " << G4endl;
796    msg << " Point A[0] (Curve   start) is " << StartPointVel << G4endl
797        << " Point S    (Sub     start) is " << SubStart_PtVel
798        << " Point A'   (Current start) is " << A_PtVel << G4endl
799        << " Point E    (Trial Point)   is " << E_Point << G4endl
800        << " Point F    (Intersection)  is " << ApproxIntersecPointV << G4endl
801        << " Point B'   (Current end)   is " << B_PtVel << G4endl
802        << " Point B[0] (Curve   end)   is " << EndPointVel << G4endl
803        << G4endl
804        << " LocateIntersection parameters are : " << G4endl
805        << "      Substep no (total) = "  << substep_no << G4endl
806        << "      Substep no         = "  << substep_no_p << " at depth= " << depth;
807    msg.precision(oldprc);
808 }
809 
810 ///////////////////////////////////////////////////////////////////////////
811 //
812 // ReportProgress.
813 //
814 void G4VIntersectionLocator::ReportProgress( std::ostream& oss,
815                                     const G4FieldTrack& StartPointVel, 
816                                     const G4FieldTrack& EndPointVel,
817                                           G4int         substep_no, 
818                                     const G4FieldTrack& A_PtVel,
819                                     const G4FieldTrack& B_PtVel,
820                                           G4double      safetyLast,
821                                           G4int         depth )
822 
823 {
824   oss << "ReportProgress: Current status of intersection search: " << G4endl;
825   if( depth > 0 ) { oss << " Depth= " << depth; }
826   oss << " Substep no = " << substep_no << G4endl;
827   G4int  verboseLevel = 5; 
828   G4double safetyPrev = -1.0;  // Add as argument ?
829 
830   printStatus( StartPointVel, EndPointVel, -1.0, -1.0, -1, 
831               oss, verboseLevel);  
832   oss << " * Start and end-point of requested Step:" << G4endl;
833   oss << " ** State of point A: "; 
834   printStatus( A_PtVel, A_PtVel, -1.0, safetyPrev, substep_no-1,
835                oss, verboseLevel);  
836   oss << " ** State of point B: "; 
837   printStatus( A_PtVel, B_PtVel, -1.0, safetyLast, substep_no, 
838                oss, verboseLevel);
839 }
840 
841 ///////////////////////////////////////////////////////////////////////////
842 //
843 // ReportImmediateHit.
844 //
845 void
846 G4VIntersectionLocator::ReportImmediateHit( const char*          MethodName, 
847                                             const G4ThreeVector& StartPosition, 
848                                             const G4ThreeVector& TrialPoint, 
849                                                   G4double       tolerance,
850                                             unsigned long int    numCalls )
851 {
852   static G4ThreadLocal unsigned int occurredOnTop= 0;
853   static G4ThreadLocal G4ThreeVector* ptrLast = nullptr;
854   if( ptrLast == nullptr )
855   {
856      ptrLast= new G4ThreeVector( DBL_MAX, DBL_MAX, DBL_MAX );
857      G4AutoDelete::Register(ptrLast);
858   }
859   G4ThreeVector &lastStart= *ptrLast;
860 
861   if( (TrialPoint - StartPosition).mag2() < tolerance*tolerance) 
862   {
863      static G4ThreadLocal unsigned int numUnmoved = 0;
864      static G4ThreadLocal unsigned int numStill = 0;    // Still at same point
865 
866      G4cout << "Intersection F == start A in " << MethodName;
867      G4cout << "Start Point: " << StartPosition << G4endl;
868      G4cout << " Start-Trial: " << TrialPoint - StartPosition;
869      G4cout << " Start-last: " << StartPosition - lastStart;
870 
871      if( (StartPosition - lastStart).mag() < tolerance )
872      {
873         // We are at position of last 'Start' position - ie unmoved
874         ++numUnmoved;
875         ++numStill; 
876         G4cout << " { Unmoved: "  << " still#= " << numStill
877                << " total # = " << numUnmoved << " } - ";
878      }
879      else
880      {
881         numStill = 0; 
882      }
883      G4cout << " Occurred: " << ++occurredOnTop;  
884      G4cout <<  " out of total calls= " << numCalls;
885      G4cout << G4endl;
886      lastStart = StartPosition;
887   }
888 }  // End of ReportImmediateHit() 
889