LArSoft  v07_13_02
Liquid Argon Software toolkit - http://larsoft.org/
trkf::InteractPlane Class Reference

#include "InteractPlane.h"

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

Public Member Functions

 InteractPlane (double tcut)
 Constructor. More...
 
virtual ~InteractPlane ()
 Destructor. More...
 
Interactorclone () const
 Clone method. More...
 
virtual bool noise (const KTrack &trk, double s, TrackError &noise_matrix) const
 Calculate noise matrix. More...
 
double getTcut () const
 

Detailed Description

Definition at line 23 of file InteractPlane.h.

Constructor & Destructor Documentation

trkf::InteractPlane::InteractPlane ( double  tcut)

Constructor.

Constructor.

Arguments:

tcut - Maximum delta ray energy.

Definition at line 26 of file InteractPlane.cxx.

Referenced by clone().

26  :
27  Interactor(tcut)
28  {}
Interactor(double tcut)
Constructor.
Definition: Interactor.cxx:21
trkf::InteractPlane::~InteractPlane ( )
virtual

Destructor.

Definition at line 31 of file InteractPlane.cxx.

32  {}

Member Function Documentation

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

Clone method.

Implements trkf::Interactor.

Definition at line 36 of file InteractPlane.h.

References InteractPlane(), noise(), and s.

36 {return new InteractPlane(*this);}
InteractPlane(double tcut)
Constructor.
double trkf::Interactor::getTcut ( ) const
inlineinherited

Definition at line 39 of file Interactor.h.

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

Referenced by noise().

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

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 70 of file InteractPlane.cxx.

References trkf::Surface::BACKWARD, e, trkf::KTrack::getDirection(), trkf::KTrack::getSurface(), trkf::Interactor::getTcut(), trkf::KTrack::getVector(), and trkf::KTrack::Mass().

Referenced by clone(), and trkf::InteractGeneral::noise().

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

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