MaterialScan implementation to shoot one ray along a defined direction and record the Material as a function of travel depth. More...

#include <MaterialScan.h>

Inheritance diagram for MaterialScanRay:

Public Member Functions
	MaterialScanRay (TFile *rootFile, const G4ThreeVector &origin, const G4ThreeVector &dir, double opening, int count, double sampleDepth, double maxDepth, bool splitByMaterials, const std::vector< std::string > &ignoredMaterials)
	Construct a new instance and set all parameters.

int	getNRays () const override
	Return the number of rays.

void	UserSteppingAction (const G4Step *step) override
	Record the material budget for each step of the particles.

bool	createNext (G4ThreeVector &origin, G4ThreeVector &direction) override
	Implement shooting along the ray.

std::string	getName () const
	Return the name of the scan.

Protected Member Functions
bool	checkStep (const G4Step *step)
	check for stuck tracks by looking at the step length

Protected Attributes
TFile *	m_rootFile
	Pointer to the root file for the histograms.

std::string	m_name
	Name of the scan, will be prefixed to all histogram names.

std::string	m_axisLabel
	Labels for the coordinate axes.

Private Member Functions
TH1D *	getHistogram (const std::string &name)
	get histogram for a given name, create if needed.

void	fillValue (const std::string &name, double value, double steplength)
	Fill the recorded material budget into the corresponding histogram.

Private Attributes
std::set< std::string >	m_ignoredMaterials
	Materials ignored when scanning.

std::map< std::string, std::unique_ptr< TH1D > >	m_regions
	Map holding pointers to all created histograms.

G4ThreeVector	m_origin
	Origin of the scan.

G4ThreeVector	m_dir
	Direction of the ray.

double	m_opening
	Opening angle in radian.

double	m_sampleDepth
	The ray length after which to sample.

double	m_maxDepth {0}
	Maximum depth for each ray after which it will be stopped.

double	m_curDepth {0}
	Current depth of the current ray.

double	m_scanDepth {0}
	The first ray does not record any material but just checks for the maximum useful depth to not get a plot which contains the PXD at the front and then continues for 500 more cm without any content.

int	m_count
	Amount of rays to shoot.

int	m_curRay { -1}
	Current Ray number: 0 = scan for maximum depth, 1..N = record materials.

bool	m_splitByMaterials
	If true Split by materials instead of regions.

int	m_zeroSteps {0}
	Count the number of steps with (almost) zero length.

Static Private Attributes
static constexpr double	c_zeroTolerance = 1e-6
	maximum Step length to be considered zero

static constexpr int	c_maxZeroStepsNudge = 10
	maximum number of consecutive zero steps before nudging the track along

static constexpr int	c_maxZeroStepsKill = 20
	maximum number of consecutive zero steps before killing the track

Detailed Description

MaterialScan implementation to shoot one ray along a defined direction and record the Material as a function of travel depth.

In contrast to the other implementations this produces a 1D histogram per region (plus a cumulative one) which contains the X_0 per bin so that the integral of the histogram ( $\sum_i binwidth_i * bincontent_i$ ) is equal to the total number of X_0 seen.

Definition at line 239 of file MaterialScan.h.

Constructor & Destructor Documentation

◆ MaterialScanRay()

MaterialScanRay	(	TFile *	rootFile,
		const G4ThreeVector &	origin,
		const G4ThreeVector &	dir,
		double	opening,
		int	count,
		double	sampleDepth,
		double	maxDepth,
		bool	splitByMaterials,
		const std::vector< std::string > &	ignoredMaterials
	)

inline

Construct a new instance and set all parameters.

Definition at line 242 of file MaterialScan.h.

                                                                 :
      MaterialScanBase(rootFile, "Ray", "ray length [cm]; material budget [X_0/cm]"),
      m_ignoredMaterials(ignoredMaterials.begin(), ignoredMaterials.end()),
      m_origin(origin), m_dir(dir), m_opening(opening), m_sampleDepth(sampleDepth),
      m_maxDepth(maxDepth), m_count(std::max(1, count)), m_splitByMaterials(splitByMaterials)
    {
      if (m_opening <= 0) m_count = 1;
    }

Member Function Documentation

◆ checkStep()

bool checkStep ( const G4Step * step )

protectedinherited

check for stuck tracks by looking at the step length

Definition at line 44 of file MaterialScan.cc.

{
  double stlen = step->GetStepLength();
  G4StepPoint* preStepPoint = step->GetPreStepPoint();
  G4Region* region = preStepPoint->GetPhysicalVolume()->GetLogicalVolume()->GetRegion();
  if (stlen < c_zeroTolerance) {
    ++m_zeroSteps;
  } else {
    m_zeroSteps = 0;
  }
  if (m_zeroSteps > c_maxZeroStepsNudge) {
    if (m_zeroSteps > c_maxZeroStepsKill) {
      B2ERROR("Track is stuck at " << preStepPoint->GetPosition() << " in volume '"
              << preStepPoint->GetPhysicalVolume()->GetLogicalVolume()->GetName()
              << " (" << region->GetName() << "): "
              << m_zeroSteps << " consecutive steps with length less then "
              << c_zeroTolerance << " mm, killing it");
      step->GetTrack()->SetTrackStatus(fStopAndKill);
    } else {
      B2WARNING("Track is stuck at " << preStepPoint->GetPosition() << " in volume '"
                << preStepPoint->GetPhysicalVolume()->GetLogicalVolume()->GetName()
                << " (" << region->GetName() << "): "
                << m_zeroSteps << " consecutive steps with length less then "
                << c_zeroTolerance << " mm, nudging it along");
      G4ThreeVector pos = step->GetTrack()->GetPosition();
      G4ThreeVector dir = step->GetTrack()->GetMomentumDirection();
      step->GetTrack()->SetPosition(pos + c_zeroTolerance * dir);
    }
    return false;
  }
  return true;
}

◆ createNext()

bool createNext	(	G4ThreeVector &	origin,
		G4ThreeVector &	direction
	)

overridevirtual

Implement shooting along the ray.

Implements MaterialScanBase.

Definition at line 215 of file MaterialScan.cc.

{
  if (m_curRay > m_count) return false;
 
  origin = m_origin;
  direction = m_dir;
  ++m_curRay;
  if (m_curRay == 0) {
    // First ray is just to check what the maximum flight length is
    return true;
  } else if (m_curRay == 1) {
    // this is the second ray, save the maximum depth for histograms
    m_maxDepth = m_maxDepth > 0 ? std::min(m_scanDepth, m_maxDepth) : m_scanDepth;
    B2INFO("Set maximum depth to " << m_maxDepth);
    // create summary histogram right now, not on demand, to make sure it is in the file
    if (m_splitByMaterials) {
      getHistogram("All_Materials_x0");
      getHistogram("All_Materials_lambda");
    } else {
      getHistogram("All_Regions_x0");
      getHistogram("All_Regions_lambda");
    }
  }
  m_curDepth = 0;
  if (m_opening > 0) {
    // add small angle to the ray where we want the opening angle to be covered
    // uniform by area so choose an angle uniform in cos(angle). We take half
    // the opening angle because it can go in both sides
    double theta = acos(gRandom->Uniform(cos(m_opening / 2.), 1));
    double phi = gRandom->Uniform(0, 2 * M_PI);
    auto orthogonal = direction.orthogonal();
    orthogonal.rotate(phi, direction);
    direction.rotate(theta, orthogonal);
  }
  B2DEBUG(10, "Create Ray " << m_curRay << " from " << origin << " direction " <<
          direction << " [theta=" << (direction.theta() / M_PI * 180) << ", phi=" << (direction.phi() / M_PI * 180) << "]");
  return true;
}

◆ fillValue()

void fillValue	(	const std::string &	name,
		double	value,
		double	steplength
	)

private

Fill the recorded material budget into the corresponding histogram.

Parameters

name	Name of the histogram
value	Value to store
steplength	The Steplength which produced the value (for correct subsampling)

Definition at line 266 of file MaterialScan.cc.

{
  TH1D* h = getHistogram(name);
  double x = m_curDepth;
  //The steps are not aligned to the bin sizes so make sure we sub-sample the
  //value and put it in the histogram in smaller steps. We choose 10 samples
  //per bin here which should be good enough.
  int samples = std::ceil(steplength / m_sampleDepth * 10);
  const double dx = steplength / samples;
  const double dv = value / samples / m_count / h->GetBinWidth(1);
  for (int i = 0; i < samples; ++i) {
    h->Fill(x * Unit::mm, dv);
    x += dx;
  }
}

◆ getHistogram()

TH1D * getHistogram ( const std::string & name )

private

get histogram for a given name, create if needed.

Parameters

name	Name of the histogram

Definition at line 254 of file MaterialScan.cc.

{
  std::unique_ptr<TH1D>& hist = m_regions[name];
  if (!hist) {
    //Create new histogram
    m_rootFile->cd(m_name.c_str());
    int bins = std::ceil(m_maxDepth / m_sampleDepth);
    hist.reset(new TH1D(name.c_str(), (name + ";" + m_axisLabel).c_str(), bins, 0, m_maxDepth * Unit::mm));
  }
  return hist.get();
}

◆ getName()

std::string getName ( ) const

inlineinherited

Return the name of the scan.

Definition at line 43 of file MaterialScan.h.

43{ return m_name; }

◆ getNRays()

int getNRays ( ) const

inlineoverridevirtual

Return the number of rays.

We have one extra ray to scan the maximum depth

Implements MaterialScanBase.

Definition at line 254 of file MaterialScan.h.

254{ return m_count + 1; }

◆ UserSteppingAction()

void UserSteppingAction ( const G4Step * step )

override

Record the material budget for each step of the particles.

Definition at line 282 of file MaterialScan.cc.

{
  //Get information about radiation and interaction length
  G4StepPoint* preStepPoint = step->GetPreStepPoint();
  G4Region* region = preStepPoint->GetPhysicalVolume()->GetLogicalVolume()->GetRegion();
  G4Material* material = preStepPoint->GetMaterial();
  const bool ignored = m_ignoredMaterials.count(material->GetName()) > 0;
  double stlen = step->GetStepLength();
  double x0 = stlen / (material->GetRadlen());
  double lambda = stlen / (material->GetNuclearInterLength());
  B2DEBUG(20, "Step in at " << preStepPoint->GetPosition() << " in volume '"
          << preStepPoint->GetPhysicalVolume()->GetLogicalVolume()->GetName()
          << " (" << region->GetName() << ")"
          << " with length=" << stlen << " mm");
  checkStep(step);
 
  //check if the depth is limited
  if (m_maxDepth > 0) {
    if (m_curDepth + stlen > m_maxDepth) {
      //Depth reached, kill track and add remaining part of the material budget
      G4Track* track = step->GetTrack();
      track->SetTrackStatus(fStopAndKill);
      double remaining = m_maxDepth - m_curDepth;
      x0 *= remaining / stlen;
      lambda *= remaining / stlen;
      stlen = remaining;
    }
  }
 
  if (m_curRay > 0 && !ignored) {
    //Fill x0 and lambda in a histogram for each region
    string x0_total = "All_Regions_x0";
    string lambda_total = "All_Regions_lambda";
    string x0_name = region->GetName() + "_x0";
    string lambda_name = region->GetName() + "_lambda";
    //or for each Material
    if (m_splitByMaterials) {
      x0_total = "All_Materials_x0";
      lambda_total = "All_Materials_lambda";
      x0_name = material->GetName() + "_x0";
      lambda_name = material->GetName() + "_lambda";
    }
    fillValue(x0_total, x0, stlen);
    fillValue(lambda_total, lambda, stlen);
    fillValue(x0_name, x0, stlen);
    fillValue(lambda_name, lambda, stlen);
  }
  m_curDepth += stlen;
  if (m_curRay == 0 && !ignored) {
    // update max depth only for materials we do not ignore ...
    m_scanDepth = m_curDepth;
  }
}

Member Data Documentation

◆ c_maxZeroStepsKill

constexpr int c_maxZeroStepsKill = 20

staticconstexprprivateinherited

maximum number of consecutive zero steps before killing the track

Definition at line 71 of file MaterialScan.h.

◆ c_maxZeroStepsNudge

constexpr int c_maxZeroStepsNudge = 10

staticconstexprprivateinherited

maximum number of consecutive zero steps before nudging the track along

Definition at line 69 of file MaterialScan.h.

◆ c_zeroTolerance

constexpr double c_zeroTolerance = 1e-6

staticconstexprprivateinherited

maximum Step length to be considered zero

Definition at line 67 of file MaterialScan.h.

◆ m_axisLabel

std::string m_axisLabel

protectedinherited

Labels for the coordinate axes.

Definition at line 64 of file MaterialScan.h.

◆ m_count

int m_count

private

Amount of rays to shoot.

Definition at line 294 of file MaterialScan.h.

◆ m_curDepth

double m_curDepth {0}

private

Current depth of the current ray.

Definition at line 286 of file MaterialScan.h.

◆ m_curRay

int m_curRay { -1}

private

Current Ray number: 0 = scan for maximum depth, 1..N = record materials.

Definition at line 296 of file MaterialScan.h.

◆ m_dir

G4ThreeVector m_dir

private

Direction of the ray.

Definition at line 278 of file MaterialScan.h.

◆ m_ignoredMaterials

std::set<std::string> m_ignoredMaterials

private

Materials ignored when scanning.

Definition at line 272 of file MaterialScan.h.

◆ m_maxDepth

double m_maxDepth {0}

private

Maximum depth for each ray after which it will be stopped.

0=no limit.

Definition at line 284 of file MaterialScan.h.

◆ m_name

std::string m_name

protectedinherited

Name of the scan, will be prefixed to all histogram names.

Definition at line 62 of file MaterialScan.h.

◆ m_opening

double m_opening

private

Opening angle in radian.

Definition at line 280 of file MaterialScan.h.

◆ m_origin

G4ThreeVector m_origin

private

Origin of the scan.

Definition at line 276 of file MaterialScan.h.

◆ m_regions

std::map<std::string, std::unique_ptr<TH1D> > m_regions

private

Map holding pointers to all created histograms.

Definition at line 274 of file MaterialScan.h.

◆ m_rootFile

TFile* m_rootFile

protectedinherited

Pointer to the root file for the histograms.

Definition at line 60 of file MaterialScan.h.

◆ m_sampleDepth

double m_sampleDepth

private

The ray length after which to sample.

Basically the bin width of the histogram

Definition at line 282 of file MaterialScan.h.

◆ m_scanDepth

double m_scanDepth {0}

private

The first ray does not record any material but just checks for the maximum useful depth to not get a plot which contains the PXD at the front and then continues for 500 more cm without any content.

This variable stores the maximum depth seen by this ray in non-ignored matierals

Definition at line 292 of file MaterialScan.h.

◆ m_splitByMaterials

bool m_splitByMaterials

private

If true Split by materials instead of regions.

Definition at line 298 of file MaterialScan.h.

◆ m_zeroSteps

int m_zeroSteps {0}

privateinherited

Count the number of steps with (almost) zero length.

Definition at line 73 of file MaterialScan.h.

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

simulation/modules/include/MaterialScan.h
simulation/modules/src/MaterialScan.cc

Public Member Functions

Protected Member Functions

Protected Attributes

Private Member Functions

Private Attributes

Static Private Attributes

Detailed Description

Constructor & Destructor Documentation

◆ MaterialScanRay()

Member Function Documentation

◆ checkStep()

◆ createNext()

◆ fillValue()

◆ getHistogram()

◆ getName()

◆ getNRays()

◆ UserSteppingAction()

Member Data Documentation

◆ c_maxZeroStepsKill

◆ c_maxZeroStepsNudge

◆ c_zeroTolerance

◆ m_axisLabel

◆ m_count

◆ m_curDepth

◆ m_curRay

◆ m_dir

◆ m_ignoredMaterials

◆ m_maxDepth

◆ m_name

◆ m_opening

◆ m_origin

◆ m_regions

◆ m_rootFile

◆ m_sampleDepth

◆ m_scanDepth

◆ m_splitByMaterials

◆ m_zeroSteps