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