DenseEnvironmentExtension.hpp

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

// This file is part of the Acts project.
//
// Copyright (C) 2018-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 <functional>

#include "Acts/Geometry/GeometryContext.hpp"
#include "Acts/MagneticField/MagneticFieldContext.hpp"
#include "Acts/Material/Interactions.hpp"
#include "Acts/Propagator/Propagator.hpp"
#include "Acts/Utilities/Helpers.hpp"

namespace Acts {

struct DenseEnvironmentExtension {
  double currentMomentum = 0.;
  double initialMomentum = 0.;
  Material material;
  std::array<double, 4> dLdl;
  std::array<double, 4> qop;
  std::array<double, 4> dPds;
  double dgdqopValue = 0.;
  double g = 0.;
  std::array<double, 4> tKi;
  std::array<double, 4> Lambdappi;
  std::array<double, 4> energy;

  DenseEnvironmentExtension() = default;

  template <typename propagator_state_t, typename stepper_t>
  int bid(const propagator_state_t& state, const stepper_t& stepper) const {
    // Check for valid particle properties
    if (stepper.charge(state.stepping) == 0. || state.options.mass == 0. ||
        stepper.momentum(state.stepping) < state.options.momentumCutOff) {
      return 0;
    }

    // Check existence of a volume with material
    if (!state.navigation.currentVolume ||
        !state.navigation.currentVolume->volumeMaterial()) {
      return 0;
    }
    return 2;
  }

  template <typename propagator_state_t, typename stepper_t>
  bool k(const propagator_state_t& state, const stepper_t& stepper,
         Vector3D& knew, const Vector3D& bField, std::array<double, 4>& kQoP,
         const int i = 0, const double h = 0.,
         const Vector3D& kprev = Vector3D()) {
    // i = 0 is used for setup and evaluation of k
    if (i == 0) {
      // Set up container for energy loss
      auto volumeMaterial = state.navigation.currentVolume->volumeMaterial();
      Vector3D position = stepper.position(state.stepping);
      material = (volumeMaterial->material(position));
      initialMomentum = stepper.momentum(state.stepping);
      currentMomentum = initialMomentum;
      qop[0] = stepper.charge(state.stepping) / initialMomentum;
      initializeEnergyLoss(state);
      // Evaluate k
      knew = qop[0] * stepper.direction(state.stepping).cross(bField);
      // Evaluate k for the time propagation
      Lambdappi[0] =
          -qop[0] * qop[0] * qop[0] * g * energy[0] /
          (stepper.charge(state.stepping) * stepper.charge(state.stepping));
      //~ tKi[0] = std::hypot(1, state.options.mass / initialMomentum);
      tKi[0] = std::hypot(1, state.options.mass * qop[0]);
      kQoP[0] = Lambdappi[0];
    } else {
      // Update parameters and check for momentum condition
      updateEnergyLoss(state.options.mass, h, state.stepping, stepper, i);
      if (currentMomentum < state.options.momentumCutOff) {
        return false;
      }
      // Evaluate k
      knew = qop[i] *
             (stepper.direction(state.stepping) + h * kprev).cross(bField);
      // Evaluate k_i for the time propagation
      double qopNew = qop[0] + h * Lambdappi[i - 1];
      Lambdappi[i] =
          -qopNew * qopNew * qopNew * g * energy[i] /
          (stepper.charge(state.stepping) * stepper.charge(state.stepping) *
           UnitConstants::C * UnitConstants::C);
      tKi[i] = std::hypot(1, state.options.mass * qopNew);
      kQoP[i] = Lambdappi[i];
    }
    return true;
  }

  template <typename propagator_state_t, typename stepper_t>
  bool finalize(propagator_state_t& state, const stepper_t& stepper,
                const double h) const {
    // Evaluate the new momentum
    double newMomentum =
        stepper.momentum(state.stepping) +
        (h / 6.) * (dPds[0] + 2. * (dPds[1] + dPds[2]) + dPds[3]);

    // Break propagation if momentum becomes below cut-off
    if (newMomentum < state.options.momentumCutOff) {
      return false;
    }

    // Add derivative dlambda/ds = Lambda''
    state.stepping.derivative(7) =
        -std::sqrt(state.options.mass * state.options.mass +
                   newMomentum * newMomentum) *
        g / (newMomentum * newMomentum * newMomentum);

    // Update momentum
    state.stepping.p = newMomentum;
    // Add derivative dt/ds = 1/(beta * c) = sqrt(m^2 * p^{-2} + c^{-2})
    state.stepping.derivative(3) =
        std::hypot(1, state.options.mass / newMomentum);
    // Update time
    state.stepping.t += (h / 6.) * (tKi[0] + 2. * (tKi[1] + tKi[2]) + tKi[3]);

    return true;
  }

  template <typename propagator_state_t, typename stepper_t>
  bool finalize(propagator_state_t& state, const stepper_t& stepper,
                const double h, FreeMatrix& D) const {
    return finalize(state, stepper, h) && transportMatrix(state, stepper, h, D);
  }

 private:
  template <typename propagator_state_t, typename stepper_t>
  bool transportMatrix(propagator_state_t& state, const stepper_t& stepper,
                       const double h, FreeMatrix& D) const {

    auto& sd = state.stepping.stepData;
    auto dir = stepper.direction(state.stepping);

    D = FreeMatrix::Identity();
    const double half_h = h * 0.5;

    // This sets the reference to the sub matrices
    // dFdx is already initialised as (3x3) zero
    auto dFdT = D.block<3, 3>(0, 4);
    auto dFdL = D.block<3, 1>(0, 7);
    // dGdx is already initialised as (3x3) identity
    auto dGdT = D.block<3, 3>(4, 4);
    auto dGdL = D.block<3, 1>(4, 7);

    ActsMatrixD<3, 3> dk1dT = ActsMatrixD<3, 3>::Zero();
    ActsMatrixD<3, 3> dk2dT = ActsMatrixD<3, 3>::Identity();
    ActsMatrixD<3, 3> dk3dT = ActsMatrixD<3, 3>::Identity();
    ActsMatrixD<3, 3> dk4dT = ActsMatrixD<3, 3>::Identity();

    ActsVectorD<3> dk1dL = ActsVectorD<3>::Zero();
    ActsVectorD<3> dk2dL = ActsVectorD<3>::Zero();
    ActsVectorD<3> dk3dL = ActsVectorD<3>::Zero();
    ActsVectorD<3> dk4dL = ActsVectorD<3>::Zero();

    std::array<double, 4> jdL;

    // Evaluation of the rightmost column without the last term.
    jdL[0] = dLdl[0];
    dk1dL = dir.cross(sd.B_first);

    jdL[1] = dLdl[1] * (1. + half_h * jdL[0]);
    dk2dL = (1. + half_h * jdL[0]) * (dir + half_h * sd.k1).cross(sd.B_middle) +
            qop[1] * half_h * dk1dL.cross(sd.B_middle);

    jdL[2] = dLdl[2] * (1. + half_h * jdL[1]);
    dk3dL = (1. + half_h * jdL[1]) * (dir + half_h * sd.k2).cross(sd.B_middle) +
            qop[2] * half_h * dk2dL.cross(sd.B_middle);

    jdL[3] = dLdl[3] * (1. + h * jdL[2]);
    dk4dL = (1. + h * jdL[2]) * (dir + h * sd.k3).cross(sd.B_last) +
            qop[3] * h * dk3dL.cross(sd.B_last);

    dk1dT(0, 1) = sd.B_first.z();
    dk1dT(0, 2) = -sd.B_first.y();
    dk1dT(1, 0) = -sd.B_first.z();
    dk1dT(1, 2) = sd.B_first.x();
    dk1dT(2, 0) = sd.B_first.y();
    dk1dT(2, 1) = -sd.B_first.x();
    dk1dT *= qop[0];

    dk2dT += half_h * dk1dT;
    dk2dT = qop[1] * VectorHelpers::cross(dk2dT, sd.B_middle);

    dk3dT += half_h * dk2dT;
    dk3dT = qop[2] * VectorHelpers::cross(dk3dT, sd.B_middle);

    dk4dT += h * dk3dT;
    dk4dT = qop[3] * VectorHelpers::cross(dk4dT, sd.B_last);

    dFdT.setIdentity();
    dFdT += h / 6. * (dk1dT + dk2dT + dk3dT);
    dFdT *= h;

    dFdL = h * h / 6. * (dk1dL + dk2dL + dk3dL);

    dGdT += h / 6. * (dk1dT + 2. * (dk2dT + dk3dT) + dk4dT);

    dGdL = h / 6. * (dk1dL + 2. * (dk2dL + dk3dL) + dk4dL);

    // Evaluation of the dLambda''/dlambda term
    D(7, 7) += (h / 6.) * (jdL[0] + 2. * (jdL[1] + jdL[2]) + jdL[3]);

    // The following comment lines refer to the application of the time being
    // treated as a position. Since t and qop are treated independently for now,
    // this just serves as entry point for building their relation
    //~ double dtpp1dl = -state.options.mass * state.options.mass * qop[0] *
    //~ qop[0] *
    //~ (3. * g + qop[0] * dgdqop(energy[0], state.options.mass,
    //~ state.options.absPdgCode,
    //~ state.options.meanEnergyLoss));

    double dtp1dl = qop[0] * state.options.mass * state.options.mass /
                    std::hypot(1, qop[0] * state.options.mass);
    double qopNew = qop[0] + half_h * Lambdappi[0];

    //~ double dtpp2dl = -state.options.mass * state.options.mass * qopNew *
    //~ qopNew *
    //~ (3. * g * (1. + half_h * jdL[0]) +
    //~ qopNew * dgdqop(energy[1], state.options.mass,
    //~ state.options.absPdgCode,
    //~ state.options.meanEnergyLoss));

    double dtp2dl = qopNew * state.options.mass * state.options.mass /
                    std::hypot(1, qopNew * state.options.mass);
    qopNew = qop[0] + half_h * Lambdappi[1];

    //~ double dtpp3dl = -state.options.mass * state.options.mass * qopNew *
    //~ qopNew *
    //~ (3. * g * (1. + half_h * jdL[1]) +
    //~ qopNew * dgdqop(energy[2], state.options.mass,
    //~ state.options.absPdgCode,
    //~ state.options.meanEnergyLoss));

    double dtp3dl = qopNew * state.options.mass * state.options.mass /
                    std::hypot(1, qopNew * state.options.mass);
    qopNew = qop[0] + half_h * Lambdappi[2];
    double dtp4dl = qopNew * state.options.mass * state.options.mass /
                    std::hypot(1, qopNew * state.options.mass);

    //~ D(3, 7) = h * state.options.mass * state.options.mass * qop[0] /
    //~ std::hypot(1., state.options.mass * qop[0])
    //~ + h * h / 6. * (dtpp1dl + dtpp2dl + dtpp3dl);

    D(3, 7) = (h / 6.) * (dtp1dl + 2. * (dtp2dl + dtp3dl) + dtp4dl);
    return true;
  }

  template <typename propagator_state_t>
  void initializeEnergyLoss(const propagator_state_t& state) {
    energy[0] = std::hypot(initialMomentum, state.options.mass);
    // use unit length as thickness to compute the energy loss per unit length
    Acts::MaterialProperties slab(material, 1);
    // Use the same energy loss throughout the step.
    if (state.options.meanEnergyLoss) {
      g = -computeEnergyLossMean(slab, state.options.absPdgCode,
                                 state.options.mass, qop[0]);
    } else {
      // TODO using the unit path length is not quite right since the most
      //      probably energy loss is not independent from the path length.
      g = -computeEnergyLossMode(slab, state.options.absPdgCode,
                                 state.options.mass, qop[0]);
    }
    // Change of the momentum per path length
    // dPds = dPdE * dEds
    dPds[0] = g * energy[0] / initialMomentum;
    if (state.stepping.covTransport) {
      // Calculate the change of the energy loss per path length and
      // inverse momentum
      if (state.options.includeGgradient) {
        if (state.options.meanEnergyLoss) {
          dgdqopValue = deriveEnergyLossMeanQOverP(
              slab, state.options.absPdgCode, state.options.mass, qop[0]);
        } else {
          // TODO path length dependence; see above
          dgdqopValue = deriveEnergyLossModeQOverP(
              slab, state.options.absPdgCode, state.options.mass, qop[0]);
        }
      }
      // Calculate term for later error propagation
      dLdl[0] = (-qop[0] * qop[0] * g * energy[0] *
                     (3. - (initialMomentum * initialMomentum) /
                               (energy[0] * energy[0])) -
                 qop[0] * qop[0] * qop[0] * energy[0] * dgdqopValue);
    }
  }

  template <typename stepper_state_t, typename stepper_t>
  void updateEnergyLoss(const double mass, const double h,
                        const stepper_state_t& state, const stepper_t& stepper,
                        const int i) {
    // Update parameters related to a changed momentum
    currentMomentum = initialMomentum + h * dPds[i - 1];
    energy[i] = std::sqrt(currentMomentum * currentMomentum + mass * mass);
    dPds[i] = g * energy[i] / currentMomentum;
    qop[i] = stepper.charge(state) / currentMomentum;
    // Calculate term for later error propagation
    if (state.covTransport) {
      dLdl[i] = (-qop[i] * qop[i] * g * energy[i] *
                     (3. - (currentMomentum * currentMomentum) /
                               (energy[i] * energy[i])) -
                 qop[i] * qop[i] * qop[i] * energy[i] * dgdqopValue);
    }
  }
};

template <typename action_list_t = ActionList<>,
          typename aborter_list_t = AbortList<>>
struct DenseStepperPropagatorOptions
    : public PropagatorOptions<action_list_t, aborter_list_t> {
  DenseStepperPropagatorOptions(
      const DenseStepperPropagatorOptions<action_list_t, aborter_list_t>&
          dspo) = default;

  DenseStepperPropagatorOptions(
      std::reference_wrapper<const GeometryContext> gctx,
      std::reference_wrapper<const MagneticFieldContext> mctx)
      : PropagatorOptions<action_list_t, aborter_list_t>(gctx, mctx) {}

  bool meanEnergyLoss = true;

  bool includeGgradient = true;

  double momentumCutOff = 0.;

  template <typename extended_aborter_list_t>
  DenseStepperPropagatorOptions<action_list_t, extended_aborter_list_t> extend(
      extended_aborter_list_t aborters) const {
    DenseStepperPropagatorOptions<action_list_t, extended_aborter_list_t>
        eoptions(this->geoContext, this->magFieldContext);
    // Copy the options over
    eoptions.direction = this->direction;
    eoptions.absPdgCode = this->absPdgCode;
    eoptions.mass = this->mass;
    eoptions.maxSteps = this->maxSteps;
    eoptions.maxStepSize = this->maxStepSize;
    eoptions.targetTolerance = this->targetTolerance;
    eoptions.pathLimit = this->pathLimit;
    eoptions.loopProtection = this->loopProtection;
    eoptions.loopFraction = this->loopFraction;
    // Output option
    eoptions.debug = this->debug;
    eoptions.debugString = this->debugString;
    eoptions.debugPfxWidth = this->debugPfxWidth;
    eoptions.debugMsgWidth = this->debugMsgWidth;
    // Stepper options
    eoptions.tolerance = this->tolerance;
    eoptions.stepSizeCutOff = this->stepSizeCutOff;
    // Action / abort list
    eoptions.actionList = this->actionList;
    eoptions.abortList = std::move(aborters);
    // Copy dense environment specific parameters
    eoptions.meanEnergyLoss = meanEnergyLoss;
    eoptions.includeGgradient = includeGgradient;
    eoptions.momentumCutOff = momentumCutOff;
    // And return the options
    return eoptions;
  }
};

}  // namespace Acts