Public Member Functions | Private Attributes | List of all members
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...
 
- Public Member Functions inherited from trkf::Interactor
 Interactor (double tcut)
 
virtual ~Interactor ()
 
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:41
Interactor(double tcut)
Definition: Interactor.cxx:21

Member Function Documentation

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

Clone method.

Implements trkf::Interactor.

Definition at line 32 of file InteractPlane.h.

33  {
34  return new InteractPlane(*this);
35  }
InteractPlane(detinfo::DetectorPropertiesData const &detProp, double tcut)
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 69 of file InteractPlane.cxx.

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

Member Data Documentation

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

Definition at line 41 of file InteractPlane.h.


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