Physics/Scattering/Scattering.hpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2018 CERN for the benefit of the Acts project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
0008
0009 #pragma once
0010
0011 #include <cmath>
0012
0013 #include "Acts/Utilities/Definitions.hpp"
0014 #include "Acts/Utilities/Helpers.hpp"
0015
0016 #include "Fatras/Kernel/detail/RandomNumberDistributions.hpp"
0017
0018 namespace Fatras {
0019
0020 /// @class Scatterering
0021 ///
0022 /// This is the (multiple) scattering plugin to the
0023 /// Physics list. It needs a scattering formula in order
0024 /// to provide the the scattering angle in 3D space.
0025 ///
0026 /// There's two options to apply the scattering
0027 /// - a parametric action that relates phi and theta (default: off)
0028 /// - an actuall out of direction scattering applying two random numbers
0029 template <typename formula_t> struct Scattering {
0030
0031   /// The flag to include scattering or not
0032   bool scattering = true;
0033
0034   /// Include the log term
0035   bool parametric = false;
0036   double projectionFactor = 1. / std::sqrt(2.);
0037
0038   /// The scattering formula
0039   formula_t angle;
0040
0041   /// This is the scattering call operator
0042   template <typename generator_t, typename detector_t, typename particle_t>
0043   std::vector<particle_t> operator()(generator_t &gen, const detector_t &det,
0044                                      particle_t &in) const {
0045
0046     // Do nothing if the flag is set to false
0047     if (not scattering) {
0048       return {};
0049     }
0050
0051     // 3D scattering angle
0052     double angle3D = angle(gen, det, in);
0053
0054     // parametric scattering
0055     if (parametric) {
0056
0057       // the initial values
0058       double theta = Acts::VectorHelpers::theta(in.momentum());
0059       double phi = Acts::VectorHelpers::phi(in.momentum());
0060       double sinTheta = (sin(theta) * sin(theta) > 10e-10) ? sin(theta) : 1.;
0061
0062       // sample them in an independent way
0063       double deltaTheta = projectionFactor * angle3D;
0064       double numDetlaPhi = 0.; //?? @THIS IS WRONG HERE !
0065       double deltaPhi = projectionFactor * numDetlaPhi / sinTheta;
0066
0067       // @todo: use bound parameter
0068       // (i) phi
0069       phi += deltaPhi;
0070       if (phi >= M_PI)
0071         phi -= M_PI;
0072       else if (phi < -M_PI)
0073         phi += M_PI;
0074       // (ii) theta
0075       theta += deltaTheta;
0076       if (theta > M_PI)
0077         theta -= M_PI;
0078       else if (theta < 0.)
0079         theta += M_PI;
0080
0081       double sphi = std::sin(phi);
0082       double cphi = std::cos(phi);
0083       double stheta = std::sin(theta);
0084       double ctheta = std::cos(theta);
0085
0086       // assign the new values
0087       in.scatter(in.p() * Acts::Vector3D(cphi * stheta, sphi * stheta, ctheta));
0088     } else {
0089
0090       /// uniform distribution
0091       UniformDist uniformDist = UniformDist(0., 1.);
0092
0093       // Create a random uniform distribution between in the intervall [0,1]
0094       double psi = 2. * M_PI * uniformDist(gen);
0095
0096       // more complex but "more true"
0097       Acts::Vector3D pDirection(in.momentum().normalized());
0098       double x = -pDirection.y();
0099       double y = pDirection.x();
0100       double z = 0.;
0101
0102       // if it runs along the z axis - no good ==> take the x axis
0103       if (pDirection.z() * pDirection.z() > 0.999999) {
0104         x = 1.;
0105         y = 0.;
0106       }
0107       // deflector direction
0108       Acts::Vector3D deflector(x, y, z);
0109       // rotate the new direction for scattering using theta and  psi
0110       Acts::RotationMatrix3D rotation;
0111       rotation = Acts::AngleAxis3D(psi, pDirection) *
0112                  Acts::AngleAxis3D(angle3D, deflector);
0113       // rotate and set a new direction to the cache
0114       in.scatter(in.p() * rotation * pDirection.normalized());
0115     }
0116     // scattering always returns an empty list
0117     // - it is a non-distructive process
0118     return {};
0119   }
0120 };
0121
0122 } // namespace Fatras