MaterialInteractor.hpp

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file MaterialInteractor.hpp

// This file is part of the Acts project.
//
// Copyright (C) 2019 CERN for the benefit of the Acts project
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.

#pragma once

#include "Acts/Definitions/Units.hpp"
#include "Acts/Geometry/TrackingVolume.hpp"
#include "Acts/Material/MaterialInteraction.hpp"
#include "Acts/Material/MaterialSlab.hpp"
#include "Acts/Propagator/detail/PointwiseMaterialInteraction.hpp"
#include "Acts/Propagator/detail/VolumeMaterialInteraction.hpp"
#include "Acts/Surfaces/Surface.hpp"

#include <sstream>

namespace Acts {

struct MaterialInteractor {
  bool multipleScattering = true;
  bool energyLoss = true;
  bool recordInteractions = false;

  using result_type = RecordedMaterial;

  template <typename propagator_state_t, typename stepper_t,
            typename navigator_t>
  void operator()(propagator_state_t& state, const stepper_t& stepper,
                  const navigator_t& navigator, result_type& result,
                  const Logger& logger) const {
    // In case of Volume material update the result of the previous step
    if (recordInteractions and not result.materialInteractions.empty() and
        not result.materialInteractions.back().volume.empty() and
        result.materialInteractions.back().updatedVolumeStep == false) {
      updateResult(state, stepper, result);
    }

    // If we are on target, everything should have been done
    if (navigator.targetReached(state.navigation)) {
      return;
    }
    // Do nothing if nothing is what is requested.
    if (not(multipleScattering or energyLoss or recordInteractions)) {
      return;
    }
    // We only have material interactions if there is potential material
    const Surface* surface = navigator.currentSurface(state.navigation);
    const TrackingVolume* volume = navigator.currentVolume(state.navigation);

    if (not(surface and surface->surfaceMaterial()) and
        not(volume and volume->volumeMaterial())) {
      return;
    }

    if (surface and surface->surfaceMaterial()) {
      // Prepare relevant input particle properties
      detail::PointwiseMaterialInteraction d(surface, state, stepper);

      // Determine the effective traversed material and its properties
      // Material exists but it's not real, i.e. vacuum; there is nothing to do
      if (not d.evaluateMaterialSlab(state, navigator)) {
        return;
      }

      // Evaluate the material effects
      d.evaluatePointwiseMaterialInteraction(multipleScattering, energyLoss);

      if (energyLoss) {
        using namespace UnitLiterals;
        ACTS_VERBOSE(d.slab << " absPdg=" << d.absPdg
                            << " mass=" << d.mass / 1_MeV << "MeV"
                            << " momentum=" << d.momentum / 1_GeV << "GeV"
                            << " energyloss=" << d.Eloss / 1_MeV << "MeV");
      }

      // To integrate process noise, we need to transport
      // the covariance to the current position in space
      if (d.performCovarianceTransport) {
        stepper.transportCovarianceToCurvilinear(state.stepping);
      }
      // Change the noise updater depending on the navigation direction
      NoiseUpdateMode mode = (state.options.direction == Direction::Forward)
                                 ? addNoise
                                 : removeNoise;
      // Apply the material interactions
      d.updateState(state, stepper, mode);
      // Record the result
      recordResult(d, result);
    } else if (recordInteractions && volume and volume->volumeMaterial()) {
      // Prepare relevant input particle properties
      detail::VolumeMaterialInteraction d(volume, state, stepper);
      // Determine the effective traversed material and its properties
      // Material exists but it's not real, i.e. vacuum; there is nothing to do
      if (not d.evaluateMaterialSlab(state, navigator)) {
        return;
      }
      // Record the result
      recordResult(d, result);
    }
  }

 private:
  void recordResult(const detail::PointwiseMaterialInteraction& d,
                    result_type& result) const {
    result.materialInX0 += d.slab.thicknessInX0();
    result.materialInL0 += d.slab.thicknessInL0();
    // Record the interaction if requested
    if (recordInteractions) {
      MaterialInteraction mi;
      mi.position = d.pos;
      mi.time = d.time;
      mi.direction = d.dir;
      mi.deltaP = d.nextP - d.momentum;
      mi.sigmaPhi2 = d.variancePhi;
      mi.sigmaTheta2 = d.varianceTheta;
      mi.sigmaQoP2 = d.varianceQoverP;
      mi.surface = d.surface;
      mi.volume = InteractionVolume();
      mi.pathCorrection = d.pathCorrection;
      mi.materialSlab = d.slab;
      result.materialInteractions.push_back(std::move(mi));
    }
  }

  void recordResult(const detail::VolumeMaterialInteraction& d,
                    result_type& result) const {
    // Record the interaction
    MaterialInteraction mi;
    mi.position = d.pos;
    mi.time = d.time;
    mi.direction = d.dir;
    mi.surface = nullptr;
    mi.volume = d.volume;
    mi.pathCorrection = d.pathCorrection;
    mi.materialSlab = d.slab;
    result.materialInteractions.push_back(std::move(mi));
  }

  template <typename propagator_state_t, typename stepper_t>
  void updateResult(propagator_state_t& state, const stepper_t& stepper,
                    result_type& result) const {
    // Update the previous interaction
    Vector3 shift = stepper.position(state.stepping) -
                    result.materialInteractions.back().position;
    double momentum = stepper.direction(state.stepping).norm();
    result.materialInteractions.back().deltaP =
        momentum - result.materialInteractions.back().direction.norm();
    result.materialInteractions.back().materialSlab.scaleThickness(
        shift.norm());
    result.materialInteractions.back().updatedVolumeStep = true;
    result.materialInX0 +=
        result.materialInteractions.back().materialSlab.thicknessInX0();
    result.materialInL0 +=
        result.materialInteractions.back().materialSlab.thicknessInL0();
  }
};

}  // namespace Acts