Belle II Software  release-08-01-10
Geant4MaterialInterface Class Reference

AbsMaterialInterface implementation for use with Geant4's navigator. More...

#include <Geant4MaterialInterface.h>

Inheritance diagram for Geant4MaterialInterface:
Collaboration diagram for Geant4MaterialInterface:

Public Member Functions

bool initTrack (double posX, double posY, double posZ, double dirX, double dirY, double dirZ) override
 Initialize the navigator at given position and with given direction. More...
 
genfit::Material getMaterialParameters () override
 Get material parameters in current material.
 
double findNextBoundary (const genfit::RKTrackRep *rep, const genfit::M1x7 &state7, double sMax, bool varField=true) override
 Make a step (following the curvature) until step length sMax or the next boundary is reached. More...
 
virtual void setDebugLvl (unsigned int lvl=1)
 

Protected Attributes

unsigned int debugLvl_
 

Private Attributes

std::unique_ptr< class G4SafeNavigatornav_
 holds a object of G4SafeNavigator, which is located in Geant4MaterialInterface.cc
 
const class G4VPhysicalVolume * currentVolume_
 the volume the extraoplation is currenly located in
 
bool m_takingFullStep = false
 stores whether to call SetGeometricallyLimitedStep() because the full step length was taken.
 

Detailed Description

AbsMaterialInterface implementation for use with Geant4's navigator.

This allows to look up the material properties from the Geant4 geometry also used for simulation purposes.

Definition at line 29 of file Geant4MaterialInterface.h.

Member Function Documentation

◆ findNextBoundary()

double findNextBoundary ( const genfit::RKTrackRep rep,
const genfit::M1x7 state7,
double  sMax,
bool  varField = true 
)
overridevirtual

Make a step (following the curvature) until step length sMax or the next boundary is reached.

After making a step to a boundary, the position has to be beyond the boundary, i.e. the current material has to be that beyond the boundary. The actual step made is returned.

Implements AbsMaterialInterface.

Definition at line 339 of file Geant4MaterialInterface.cc.

343 {
344  // This assumes that sMax is small enough to take only a single RK
345  // step. This restriction comes about because RKPropagate only
346  // takes one step.
347  const double delta(1.E-2); // cm, distance limit beneath which straight-line steps are taken.
348  const double epsilon(1.E-1); // cm, allowed upper bound on arch deviation from straight line
349 
350  genfit::M1x3 SA;
351  genfit::M1x7 state7, oldState7;
352  oldState7 = stateOrig;
353 
354  int stepSign(sMax < 0 ? -1 : 1);
355 
356  G4ThreeVector pointOld(stateOrig[0] * CLHEP::cm,
357  stateOrig[1] * CLHEP::cm,
358  stateOrig[2] * CLHEP::cm);
359  G4ThreeVector dirOld(stepSign * stateOrig[3],
360  stepSign * stateOrig[4],
361  stepSign * stateOrig[5]);
362 
363  double s = 0; // trajectory length to boundary
364 
365  const unsigned maxIt = 300;
366  unsigned it = 0;
367 
368  // Initialize the geometry to the current location (set by caller).
369  double safety;
370 
371  m_takingFullStep = false;
372  double slDist = nav_->CheckNextStep(pointOld, dirOld, fabs(sMax) * CLHEP::cm, safety);
373  if (slDist == kInfinity) {
374  slDist = fabs(sMax);
375  m_takingFullStep = true;
376  } else {
377  slDist /= CLHEP::cm;
378  }
379  safety /= CLHEP::cm;
380 
381  // No boundary in sight?
382  if (safety > fabs(sMax)) {
383  B2DEBUG(20, " next boundary is farther away than sMax ");
384  return stepSign * safety; // sMax
385  }
386 
387  // Are we at the boundary?
388  if (slDist < delta) {
389  B2DEBUG(20, " very close to the boundary -> return @ it " << it
390  << " stepSign*slDist = "
391  << stepSign << "*" << slDist);
392  return stepSign * slDist;
393  }
394  double step = slDist;
395 
396  while (1) {
397  if (++it > maxIt) {
398  genfit::Exception exc("Geant4MaterialInterface::findNextBoundary ==> maximum number of iterations exceeded", __LINE__, __FILE__);
399  exc.setFatal();
400  throw exc;
401  }
402 
403  if (step < delta)
404  return stepSign * (s + step);
405 
406  // We have to find whether there's any boundary on our path.
407 
408  // Follow curved arch, then see if we may have missed a boundary.
409  // Always propagate complete way from original start to avoid
410  // inconsistent extrapolations. This is always a single RK step.
411  state7 = stateOrig;
412  rep->RKPropagate(state7, nullptr, SA, stepSign * (s + step), varField);
413 
414  G4ThreeVector pos(state7[0] * CLHEP::cm, state7[1] * CLHEP::cm, state7[2] * CLHEP::cm);
415  G4ThreeVector dir(stepSign * state7[3], stepSign * state7[4], stepSign * state7[5]);
416 
417  // Straight line distance between extrapolation finish and
418  // the end of the previously determined safe segment.
419  double dist2 = (pow(state7[0] - oldState7[0], 2)
420  + pow(state7[1] - oldState7[1], 2)
421  + pow(state7[2] - oldState7[2], 2));
422 
423  // If we moved less than safety, the volume cannot possibly have
424  // changed, so we skip further checks.
425  if (dist2 > safety * safety) {
426 
427  // Do we need to try again with a shorter step? There are two
428  // possible reasons:
429  //
430  // 1. Too much curvature?
431 
432  // Maximal lateral deviation² of the curved path from the
433  // straight line connecting beginning and end.
434  double maxDeviation2 = 0.25 * (step * step - dist2);
435  if (maxDeviation2 > epsilon * epsilon) {
436  // Need to take a shorter step to reliably estimate material,
437  // but only if we didn't move by safety.
438 
439  // Take a shorter step, but never shorter than safety.
440  step = safety + 0.5 * (step - safety);
441 
442  continue;
443  }
444 
445  // 2. Volume changed?
446  //
447  // Where are we after the step?
448 
449  if (m_takingFullStep) {
450  m_takingFullStep = false;
451  nav_->SetGeometricallyLimitedStep();
452  }
453 
454  std::unique_ptr<G4TouchableHistory> hist(nav_->CreateTouchableHistory());
455  G4VPhysicalVolume* newVolume = nav_->LocateGlobalPointAndSetup(pos, &dir);
456 
457  if (newVolume != currentVolume_) {
458 
459  // Volume changed during the extrapolation.
460 
461  // Extrapolation may not take the exact step length we asked
462  // for, so it can happen that a requested step < safety takes
463  // us across the boundary. This is then the best estimate we
464  // can get of the distance to the boundary with the stepper.
465  if (step <= safety)
466  return stepSign * (s + step);
467 
468  // Move back to last good point, but looking in the actual
469  // direction of the step.
470  G4ThreeVector dirCloser(pos - pointOld);
471  dirCloser.setMag(1.);
472  nav_->ResetHierarchyAndLocate(pointOld, dirCloser, *hist);
473 
474  // Look along the secant of the actual trajectory instead of
475  // the original direction. There should be a crossing within
476  // distance step.
477  double secantDist = nav_->CheckNextStep(pointOld, dirCloser,
478  step * CLHEP::cm, safety) / CLHEP::cm;
479  safety /= CLHEP::cm;
480  if (secantDist >= step) {
481  // Cannot be. Just take a shorter step, and hope that this
482  // works.
483  slDist = secantDist;
484  step = std::max(0.9 * step, safety);
485  } else {
486  slDist = step = std::max(secantDist, safety);
487  }
488 
489  // Are we at the boundary?
490  if (slDist < delta) {
491  B2DEBUG(20, " very close to the boundary -> return @ it " << it
492  << " stepSign*(s + slDist) = "
493  << stepSign << "*(" << s + slDist << ")");
494  return stepSign * (s + slDist);
495  }
496 
497  continue;
498  }
499  }
500 
501  // We're in the new place, the step was safe, advance.
502  s += step;
503 
504  oldState7 = state7;
505  pointOld = pos;
506  if (m_takingFullStep) {
507  m_takingFullStep = false;
508  // inform the navigator that the full geometrical step was taken. This is required for
509  // some Geant4 volume enter/exit optimizations to work.
510  nav_->SetGeometricallyLimitedStep();
511  }
512  nav_->LocateGlobalPointAndSetup(pos, &dir);
513  slDist = nav_->CheckNextStep(pos, dir,
514  (fabs(sMax) - s) * CLHEP::cm, safety);
515  if (slDist == kInfinity) {
516  m_takingFullStep = true;
517  slDist = fabs(sMax) - s;
518  } else {
519  slDist /= CLHEP::cm;
520  }
521  safety /= CLHEP::cm;
522  step = slDist;
523 
524  // No boundary in sight?
525  if (s + safety > fabs(sMax)) {
526  B2DEBUG(20, " next boundary is farther away than sMax");
527  return stepSign * (s + safety); // sMax
528  }
529 
530  // Are we at the boundary?
531  if (slDist < delta) {
532  B2DEBUG(20, " very close to the boundary -> return @ it " << it
533  << " stepSign*(s + slDist) = "
534  << stepSign << "*(" << s + slDist << ")");
535  return stepSign * (s + slDist);
536  }
537  }
538 }
bool m_takingFullStep
stores whether to call SetGeometricallyLimitedStep() because the full step length was taken.
const class G4VPhysicalVolume * currentVolume_
the volume the extraoplation is currenly located in
std::unique_ptr< class G4SafeNavigator > nav_
holds a object of G4SafeNavigator, which is located in Geant4MaterialInterface.cc
Exception class for error handling in GENFIT (provides storage for diagnostic information)
Definition: Exception.h:48
virtual double RKPropagate(M1x7 &state7, M7x7 *jacobian, M1x3 &SA, double S, bool varField=true, bool calcOnlyLastRowOfJ=false) const
The actual Runge Kutta propagation.
Definition: RKTrackRep.cc:1274

◆ initTrack()

bool initTrack ( double  posX,
double  posY,
double  posZ,
double  dirX,
double  dirY,
double  dirZ 
)
overridevirtual

Initialize the navigator at given position and with given direction.

Returns true if the volume changed.

Implements AbsMaterialInterface.

Definition at line 289 of file Geant4MaterialInterface.cc.


The documentation for this class was generated from the following files: