Geant4MaterialInterface Class Reference

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

#include <Geant4MaterialInterface.h>

Inheritance 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.
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.

Private Attributes
std::unique_ptr< class G4SafeNavigator >	nav_
	holds a object of G4SafeNavigator, which is located in Geant4MaterialInterface.cc
const class G4VPhysicalVolume *	currentVolume_
	the volume the extrapolation is currently 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.

Constructor & Destructor Documentation

Geant4MaterialInterface()

Geant4MaterialInterface

(

)

Definition at line 276 of file Geant4MaterialInterface.cc.

  : nav_(new G4SafeNavigator()), currentVolume_(0)
{
  G4VPhysicalVolume* world = geometry::GeometryManager::getInstance().getTopVolume();
  nav_->SetWorldVolume(world);
}

~Geant4MaterialInterface()

~Geant4MaterialInterface

(

)

Definition at line 283 of file Geant4MaterialInterface.cc.

{
}

Member Function Documentation

findNextBoundary()

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.

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.

Definition at line 339 of file Geant4MaterialInterface.cc.

{
  // This assumes that sMax is small enough to take only a single RK
  // step.  This restriction comes about because RKPropagate only
  // takes one step.
  const double delta(1.E-2); // cm, distance limit beneath which straight-line steps are taken.
  const double epsilon(1.E-1); // cm, allowed upper bound on arch deviation from straight line
 
  genfit::M1x3 SA;
  genfit::M1x7 state7, oldState7;
  oldState7 = stateOrig;
 
  int stepSign(sMax < 0 ? -1 : 1);
 
  G4ThreeVector pointOld(stateOrig[0] * CLHEP::cm,
                         stateOrig[1] * CLHEP::cm,
                         stateOrig[2] * CLHEP::cm);
  G4ThreeVector dirOld(stepSign * stateOrig[3],
                       stepSign * stateOrig[4],
                       stepSign * stateOrig[5]);
 
  double s = 0;  // trajectory length to boundary
 
  const unsigned maxIt = 300;
  unsigned it = 0;
 
  // Initialize the geometry to the current location (set by caller).
  double safety;
 
  m_takingFullStep = false;
  double slDist = nav_->CheckNextStep(pointOld, dirOld, fabs(sMax) * CLHEP::cm, safety);
  if (slDist == kInfinity) {
    slDist = fabs(sMax);
    m_takingFullStep = true;
  } else {
    slDist /= CLHEP::cm;
  }
  safety /= CLHEP::cm;
 
  // No boundary in sight?
  if (safety > fabs(sMax)) {
    B2DEBUG(20, "   next boundary is farther away than sMax ");
    return stepSign * safety; // sMax
  }
 
  // Are we at the boundary?
  if (slDist < delta) {
    B2DEBUG(20, "   very close to the boundary -> return @ it " << it
            << " stepSign*slDist = "
            << stepSign << "*" << slDist);
    return stepSign * slDist;
  }
  double step = slDist;
 
  while (1) {
    if (++it > maxIt) {
      genfit::Exception exc("Geant4MaterialInterface::findNextBoundary ==> maximum number of iterations exceeded", __LINE__, __FILE__);
      exc.setFatal();
      throw exc;
    }
 
    if (step < delta)
      return stepSign * (s + step);
 
    // We have to find whether there's any boundary on our path.
 
    // Follow curved arch, then see if we may have missed a boundary.
    // Always propagate complete way from original start to avoid
    // inconsistent extrapolations.  This is always a single RK step.
    state7 = stateOrig;
    rep->RKPropagate(state7, nullptr, SA, stepSign * (s + step), varField);
 
    G4ThreeVector pos(state7[0] * CLHEP::cm, state7[1] * CLHEP::cm, state7[2] * CLHEP::cm);
    G4ThreeVector dir(stepSign * state7[3], stepSign * state7[4], stepSign * state7[5]);
 
    // Straight line distance between extrapolation finish and
    // the end of the previously determined safe segment.
    double dist2 = (pow(state7[0] - oldState7[0], 2)
                    + pow(state7[1] - oldState7[1], 2)
                    + pow(state7[2] - oldState7[2], 2));
 
    // If we moved less than safety, the volume cannot possibly have
    // changed, so we skip further checks.
    if (dist2 > safety * safety) {
 
      // Do we need to try again with a shorter step?  There are two
      // possible reasons:
      //
      // 1. Too much curvature?
 
      // Maximal lateral deviation² of the curved path from the
      // straight line connecting beginning and end.
      double maxDeviation2 = 0.25 * (step * step - dist2);
      if (maxDeviation2 > epsilon * epsilon) {
        // Need to take a shorter step to reliably estimate material,
        // but only if we didn't move by safety.
 
        // Take a shorter step, but never shorter than safety.
        step = safety + 0.5 * (step - safety);
 
        continue;
      }
 
      // 2. Volume changed?
      //
      // Where are we after the step?
 
      if (m_takingFullStep) {
        m_takingFullStep = false;
        nav_->SetGeometricallyLimitedStep();
      }
 
      std::unique_ptr<G4TouchableHistory> hist(nav_->CreateTouchableHistory());
      G4VPhysicalVolume* newVolume = nav_->LocateGlobalPointAndSetup(pos, &dir);
 
      if (newVolume != currentVolume_) {
 
        // Volume changed during the extrapolation.
 
        // Extrapolation may not take the exact step length we asked
        // for, so it can happen that a requested step < safety takes
        // us across the boundary.  This is then the best estimate we
        // can get of the distance to the boundary with the stepper.
        if (step <= safety)
          return stepSign * (s + step);
 
        // Move back to last good point, but looking in the actual
        // direction of the step.
        G4ThreeVector dirCloser(pos - pointOld);
        dirCloser.setMag(1.);
        nav_->ResetHierarchyAndLocate(pointOld, dirCloser, *hist);
 
        // Look along the secant of the actual trajectory instead of
        // the original direction. There should be a crossing within
        // distance step.
        double secantDist = nav_->CheckNextStep(pointOld, dirCloser,
                                                step * CLHEP::cm, safety) / CLHEP::cm;
        safety /= CLHEP::cm;
        if (secantDist >= step) {
          // Cannot be.  Just take a shorter step, and hope that this
          // works.
          slDist = secantDist;
          step = std::max(0.9 * step, safety);
        } else {
          slDist = step = std::max(secantDist, safety);
        }
 
        // Are we at the boundary?
        if (slDist < delta) {
          B2DEBUG(20, "   very close to the boundary -> return @ it " << it
                  << " stepSign*(s + slDist) = "
                  << stepSign << "*(" << s + slDist << ")");
          return stepSign * (s + slDist);
        }
 
        continue;
      }
    }
 
    // We're in the new place, the step was safe, advance.
    s += step;
 
    oldState7 = state7;
    pointOld = pos;
    if (m_takingFullStep) {
      m_takingFullStep = false;
      // inform the navigator that the full geometrical step was taken. This is required for
      // some Geant4 volume enter/exit optimizations to work.
      nav_->SetGeometricallyLimitedStep();
    }
    nav_->LocateGlobalPointAndSetup(pos, &dir);
    slDist = nav_->CheckNextStep(pos, dir,
                                 (fabs(sMax) - s) * CLHEP::cm, safety);
    if (slDist == kInfinity) {
      m_takingFullStep = true;
      slDist = fabs(sMax) - s;
    } else {
      slDist /= CLHEP::cm;
    }
    safety /= CLHEP::cm;
    step = slDist;
 
    // No boundary in sight?
    if (s + safety > fabs(sMax)) {
      B2DEBUG(20, "   next boundary is farther away than sMax");
      return stepSign * (s + safety); // sMax
    }
 
    // Are we at the boundary?
    if (slDist < delta) {
      B2DEBUG(20, "   very close to the boundary -> return @ it " << it
              << " stepSign*(s + slDist) = "
              << stepSign << "*(" << s + slDist << ")");
      return stepSign * (s + slDist);
    }
  }
}

getMaterialParameters()

genfit::Material getMaterialParameters ( )

override

Get material parameters in current material.

Definition at line 309 of file Geant4MaterialInterface.cc.

{
  assert(currentVolume_);
 
  const G4Material* mat = currentVolume_->GetLogicalVolume()->GetMaterial();
 
  double density, Z, A, radiationLength, mEE;
  if (mat->GetNumberOfElements() == 1) {
    Z = mat->GetZ();
    A = mat->GetA();
  } else {
    // Calculate weight-averaged A, Z
    A = Z = 0;
    for (unsigned i = 0; i < mat->GetNumberOfElements(); ++i) {
      const G4Element* element = (*mat->GetElementVector())[i];
      Z += element->GetZ() * mat->GetFractionVector()[i];
      A += element->GetA() * mat->GetFractionVector()[i];
    }
  }
 
  density = mat->GetDensity() / CLHEP::g * CLHEP::cm3;
  // Z has correct units
  A *= CLHEP::mole / CLHEP::g;
  radiationLength = mat->GetRadlen() / CLHEP::cm;
  mEE = mat->GetIonisation()->GetMeanExcitationEnergy() / CLHEP::eV;
  return genfit::Material(density, Z, A, radiationLength, mEE);
}

initTrack()

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.

Returns true if the volume changed.

Definition at line 289 of file Geant4MaterialInterface.cc.

{
  G4ThreeVector pos(posX * CLHEP::cm, posY * CLHEP::cm, posZ * CLHEP::cm);
  G4ThreeVector dir(dirX, dirY, dirZ);
 
  if (m_takingFullStep) {
    m_takingFullStep = false;
    nav_->SetGeometricallyLimitedStep();
  }
 
  const G4VPhysicalVolume* newVolume = nav_->LocateGlobalPointAndSetup(pos, &dir);
  bool volChanged = newVolume != currentVolume_;
  currentVolume_ = newVolume;
 
  return volChanged;
}

Member Data Documentation

currentVolume_

const class G4VPhysicalVolume* currentVolume_

private

the volume the extrapolation is currently located in

Definition at line 63 of file Geant4MaterialInterface.h.

m_takingFullStep

bool m_takingFullStep = false

private

stores whether to call SetGeometricallyLimitedStep() because the full step length was taken.

Definition at line 68 of file Geant4MaterialInterface.h.

nav_

std::unique_ptr<class G4SafeNavigator> nav_

private

holds a object of G4SafeNavigator, which is located in Geant4MaterialInterface.cc

Definition at line 60 of file Geant4MaterialInterface.h.

---

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

• tracking/modules/genfitUtilities/include/Geant4MaterialInterface.h
• tracking/modules/genfitUtilities/src/Geant4MaterialInterface.cc