LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
trkf::InteractPlane Class Reference

#include "InteractPlane.h"

Inheritance diagram for trkf::InteractPlane:
trkf::Interactor

Public Member Functions

 InteractPlane (detinfo::DetectorPropertiesData const &detProp, double tcut)
 
Interactorclone () const override
 Clone method. More...
 
bool noise (const KTrack &trk, double s, TrackError &noise_matrix) const override
 Calculate noise matrix. More...
 
double getTcut () const
 

Private Attributes

detinfo::DetectorPropertiesData const & fDetProp
 

Detailed Description

Definition at line 27 of file InteractPlane.h.

Constructor & Destructor Documentation

trkf::InteractPlane::InteractPlane ( detinfo::DetectorPropertiesData const &  detProp,
double  tcut 
)

Constructor.

Arguments:

tcut - Maximum delta ray energy.

Definition at line 28 of file InteractPlane.cxx.

29  : Interactor(tcut), fDetProp{detProp}
30  {}
detinfo::DetectorPropertiesData const & fDetProp
Definition: InteractPlane.h:37
Interactor(double tcut)
Definition: Interactor.cxx:21

Member Function Documentation

Interactor* trkf::InteractPlane::clone ( ) const
inlineoverridevirtual

Clone method.

Implements trkf::Interactor.

Definition at line 31 of file InteractPlane.h.

References noise().

31 { return new InteractPlane(*this); }
InteractPlane(detinfo::DetectorPropertiesData const &detProp, double tcut)
double trkf::Interactor::getTcut ( ) const
inlineinherited

Definition at line 32 of file Interactor.h.

References trkf::Interactor::clone(), trkf::Interactor::fTcut, and trkf::Interactor::noise().

Referenced by noise().

32 { return fTcut; }
double fTcut
Maximum delta ray energy for dE/dx.
Definition: Interactor.h:41
bool trkf::InteractPlane::noise ( const KTrack trk,
double  s,
TrackError noise_matrix 
) const
overridevirtual

Calculate noise matrix.

Calculate noise matrix.

Arguments:

trk - Original track. s - Path distance. noise_matrix - Resultant noise matrix.

Returns: True if success.

Currently calculate noise from multiple scattering only.

Note about multiple scattering calculation:

In the case of normal incident track (u' = v' = 0), the multiple scattering calculations used in this class reduce to the familiar small-angle formulas found in the particle data book for a thick scatterer (meaning multiple scattering modifies both position and slope errors). However, the distance used in the logarithm factor of the rms scattering angle is an estimate of the total track length, rather than the incremental track length.

For non-normal incident track, the error ellipse is elongated in the radial direciton of the uv plane by factor sqrt(1 + u'^2 + v'^2). This is equivalent to expansion in the u direction by factor sqrt(1 + u'^2), and expansion in the v direction by sqrt(1 + v'^2), with uv correlation u' v' / sqrt((1 + u'^2)(1 + v'^2)).

Correlation between position and slope in the same view is sqrt(3)/2 regardless of normal incidence.

Correlation between position and slope in the opposite view is (sqrt(3)/2) u' v' / sqrt((1 + u'^2)(1 + v'^2))

Implements trkf::Interactor.

Definition at line 68 of file InteractPlane.cxx.

References util::abs(), trkf::Surface::BACKWARD, detinfo::DetectorPropertiesData::Density(), e, detinfo::DetectorPropertiesData::Eloss(), detinfo::DetectorPropertiesData::ElossVar(), fDetProp, trkf::KTrack::getDirection(), trkf::KTrack::getSurface(), trkf::Interactor::getTcut(), trkf::KTrack::getVector(), and trkf::KTrack::Mass().

Referenced by trkf::InteractGeneral::noise().

69  {
70  // Get LAr service.
71 
72  auto const* larprop = lar::providerFrom<detinfo::LArPropertiesService>();
73 
74  // Make sure we are on a plane surface (throw exception if not).
75 
76  const SurfPlane* psurf = dynamic_cast<const SurfPlane*>(&*trk.getSurface());
77  if (psurf == nullptr)
78  throw cet::exception("InteractPlane") << "InteractPlane called for non-planar surface.\n";
79 
80  // Clear noise matrix.
81 
82  noise_matrix.clear();
83 
84  // Unpack track parameters.
85 
86  const TrackVector& vec = trk.getVector();
87  double dudw = vec[2];
88  double dvdw = vec[3];
89  double pinv = vec[4];
90  double mass = trk.Mass();
91 
92  // If distance is zero, or momentum is infinite, return zero noise.
93 
94  if (pinv == 0. || s == 0.) return true;
95 
96  // Make a crude estimate of the range of the track.
97 
98  double p = 1. / std::abs(pinv);
99  double p2 = p * p;
100  double e2 = p2 + mass * mass;
101  double e = std::sqrt(e2);
102  double t = e - mass;
103  double dedx = 0.001 * fDetProp.Eloss(p, mass, getTcut());
104  double range = t / dedx;
105  if (range > 100.) range = 100.;
106 
107  // Calculate the radiation length in cm.
108 
109  double x0 = larprop->RadiationLength() / fDetProp.Density();
110 
111  // Calculate projected rms scattering angle.
112  // Use the estimted range in the logarithm factor.
113  // Use the incremental propagation distance in the square root factor.
114 
115  double betainv = std::sqrt(1. + pinv * pinv * mass * mass);
116  double theta_fact = (0.0136 * pinv * betainv) * (1. + 0.038 * std::log(range / x0));
117  double theta02 = theta_fact * theta_fact * std::abs(s / x0);
118 
119  // Calculate some sommon factors needed for multiple scattering.
120 
121  double ufact2 = 1. + dudw * dudw;
122  double vfact2 = 1. + dvdw * dvdw;
123  double uvfact2 = 1. + dudw * dudw + dvdw * dvdw;
124  double uvfact = std::sqrt(uvfact2);
125  double uv = dudw * dvdw;
126  double dist2_3 = s * s / 3.;
127  double dist_2 = std::abs(s) / 2.;
128  if (trk.getDirection() == Surface::BACKWARD) dist_2 = -dist_2;
129 
130  // Calculate energy loss fluctuations.
131 
132  double evar = 1.e-6 * fDetProp.ElossVar(p, mass) * std::abs(s); // E variance (GeV^2).
133  double pinvvar = evar * e2 / (p2 * p2 * p2); // Inv. p variance (1/GeV^2)
134 
135  // Fill elements of noise matrix.
136 
137  // Position submatrix.
138 
139  noise_matrix(0, 0) = dist2_3 * theta02 * ufact2; // sigma^2(u,u)
140  noise_matrix(1, 0) = dist2_3 * theta02 * uv; // sigma^2(u,v)
141  noise_matrix(1, 1) = dist2_3 * theta02 * vfact2; // sigma^2(v,v)
142 
143  // Slope submatrix.
144 
145  noise_matrix(2, 2) = theta02 * uvfact2 * ufact2; // sigma^2(u', u')
146  noise_matrix(3, 2) = theta02 * uvfact2 * uv; // sigma^2(v', u')
147  noise_matrix(3, 3) = theta02 * uvfact2 * vfact2; // sigma^2(v', v')
148 
149  // Same-view position-slope correlations.
150 
151  noise_matrix(2, 0) = dist_2 * theta02 * uvfact * ufact2; // sigma^2(u', u)
152  noise_matrix(3, 1) = dist_2 * theta02 * uvfact * vfact2; // sigma^2(v', v)
153 
154  // Opposite-view position-slope correlations.
155 
156  noise_matrix(2, 1) = dist_2 * theta02 * uvfact * uv; // sigma^2(u', v)
157  noise_matrix(3, 0) = dist_2 * theta02 * uvfact * uv; // sigma^2(v', u)
158 
159  // Momentum correlations (zero).
160 
161  noise_matrix(4, 0) = 0.; // sigma^2(pinv, u)
162  noise_matrix(4, 1) = 0.; // sigma^2(pinv, v)
163  noise_matrix(4, 2) = 0.; // sigma^2(pinv, u')
164  noise_matrix(4, 3) = 0.; // sigma^2(pinv, v')
165 
166  // Energy loss fluctuations.
167 
168  noise_matrix(4, 4) = pinvvar; // sigma^2(pinv, pinv)
169 
170  // Done (success).
171 
172  return true;
173  }
detinfo::DetectorPropertiesData const & fDetProp
Definition: InteractPlane.h:37
double ElossVar(double mom, double mass) const
Energy loss fluctuation ( )
constexpr auto abs(T v)
Returns the absolute value of the argument.
double getTcut() const
Definition: Interactor.h:32
double Density(double temperature=0.) const
Returns argon density at a given temperature.
KVector< 5 >::type TrackVector
Track state vector, dimension 5.
double Eloss(double mom, double mass, double tcut) const
Restricted mean energy loss (dE/dx)
Float_t e
Definition: plot.C:35
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33

Member Data Documentation

detinfo::DetectorPropertiesData const& trkf::InteractPlane::fDetProp
private

Definition at line 37 of file InteractPlane.h.

Referenced by noise().


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