LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
TrackCalorimetryAlg.cxx
Go to the documentation of this file.
1 
10 #include "TrackCalorimetryAlg.h"
11 #include "fhiclcpp/ParameterSet.h"
12 
14 #include "lardata/ArtDataHelper/TrackUtils.h" // lar::util::TrackPitchInView()
21 
23  : caloAlg(p.get<fhicl::ParameterSet>("CalorimetryAlg"))
24 {
25  fNHitsToDetermineStart = p.get<unsigned int>("NHitsToDetermineStart", 3);
26 }
27 
29  detinfo::DetectorClocksData const& clock_data,
30  detinfo::DetectorPropertiesData const& det_prop,
31  std::vector<recob::Track> const& trackVector,
32  std::vector<recob::Hit> const& hitVector,
33  std::vector<std::vector<size_t>> const& hit_indices_per_track,
34  std::vector<anab::Calorimetry>& caloVector,
35  std::vector<size_t>& assnTrackCaloVector,
36  Providers_t providers)
37 {
38  auto const& geom = *providers.get<geo::GeometryCore>();
39 
40  //loop over the track list
41  for (size_t i_track = 0; i_track < trackVector.size(); i_track++) {
42 
43  recob::Track const& track = trackVector[i_track];
44  std::vector<float> path_length_fraction_vec(CreatePathLengthFractionVector(track));
45 
46  //sort hits into each plane
47  std::vector<std::vector<size_t>> hit_indices_per_plane(geom.Nplanes());
48  for (auto const& i_hit : hit_indices_per_track[i_track])
49  hit_indices_per_plane[hitVector[i_hit].WireID().Plane].push_back(i_hit);
50 
51  //loop over the planes
52  for (unsigned int i_plane = 0; i_plane < geom.Nplanes(); ++i_plane) {
53  geo::PlaneID const planeID{0, 0, i_plane};
54 
55  //project down the track into wire/tick space for this plane
56  std::vector<std::pair<geo::WireID, float>> traj_points_in_plane(
57  track.NumberTrajectoryPoints());
58  for (size_t i_trjpt = 0; i_trjpt < track.NumberTrajectoryPoints(); i_trjpt++) {
59  double x_pos = track.LocationAtPoint(i_trjpt).X();
60  float tick = det_prop.ConvertXToTicks(x_pos, (int)i_plane, 0, 0);
61  traj_points_in_plane[i_trjpt] =
62  std::make_pair(geom.NearestWireID(track.LocationAtPoint(i_trjpt), planeID), tick);
63  }
64 
65  HitPropertiesMultiset_t HitPropertiesMultiset;
66  //now loop through hits
67  for (auto const& i_hit : hit_indices_per_plane[i_plane])
68  AnalyzeHit(clock_data,
69  det_prop,
70  hitVector[i_hit],
71  track,
72  traj_points_in_plane,
73  path_length_fraction_vec,
74  HitPropertiesMultiset,
75  geom);
76 
77  //PrintHitPropertiesMultiset(HitPropertiesMultiset);
79  HitPropertiesMultiset, track, i_track, caloVector, assnTrackCaloVector, planeID);
80 
81  } //end loop over planes
82 
83  } //end loop over tracks
84 
85 } //end ExtractCalorimetry
86 
88 public:
89  dist_projected(recob::Hit const& h, geo::GeometryCore const& g) : hit(h), geom(g) {}
90  bool operator()(std::pair<geo::WireID, float> const& i, std::pair<geo::WireID, float> const& j)
91  {
92  float dw_i =
93  ((int)(i.first.Wire) - (int)(hit.WireID().Wire)) * geom.WirePitch(i.first.asPlaneID());
94  float dw_j =
95  ((int)(j.first.Wire) - (int)(hit.WireID().Wire)) * geom.WirePitch(j.first.asPlaneID());
96  float dt_i = i.second - hit.PeakTime();
97  float dt_j = j.second - hit.PeakTime();
98  return std::hypot(dw_i, dt_i) < std::hypot(dw_j, dt_j);
99  }
100 
101 private:
102  recob::Hit const& hit;
104 };
105 
107  recob::Track const& track)
108 {
109 
110  std::vector<float> trk_path_length_frac_vec(track.NumberTrajectoryPoints());
111 
112  float cumulative_path_length = 0;
113  const float total_path_length = track.Length();
114  for (size_t i_trj = 1; i_trj < track.NumberTrajectoryPoints(); i_trj++) {
115  cumulative_path_length += (track.LocationAtPoint(i_trj) - track.LocationAtPoint(i_trj - 1)).R();
116  trk_path_length_frac_vec[i_trj] = cumulative_path_length / total_path_length;
117  }
118 
119  return trk_path_length_frac_vec;
120 }
121 
123  detinfo::DetectorClocksData const& clock_data,
124  detinfo::DetectorPropertiesData const& det_prop,
125  recob::Hit const& hit,
126  recob::Track const& track,
127  std::vector<std::pair<geo::WireID, float>> const& traj_points_in_plane,
128  std::vector<float> const& path_length_fraction_vec,
129  HitPropertiesMultiset_t& HitPropertiesMultiset,
130  geo::GeometryCore const& geom)
131 {
132 
133  size_t traj_iter = std::distance(traj_points_in_plane.begin(),
134  std::min_element(traj_points_in_plane.begin(),
135  traj_points_in_plane.end(),
136  dist_projected(hit, geom)));
137  float pitch = lar::util::TrackPitchInView(track, geom.View(hit.WireID().Plane), traj_iter);
138 
139  HitPropertiesMultiset.emplace(hit.Integral(),
140  hit.Integral() / pitch,
141  caloAlg.dEdx_AREA(clock_data, det_prop, hit, pitch),
142  pitch,
143  track.LocationAtPoint<TVector3>(traj_iter),
144  path_length_fraction_vec[traj_iter]);
145 }
146 
148 {
149 
150  if (hpm.size() <= fNHitsToDetermineStart) return false;
151 
152  float charge_start = 0, charge_end = 0;
153  unsigned int counter = 0;
154  for (HitPropertiesMultiset_t::iterator it_hpm = hpm.begin(); it_hpm != hpm.end(); it_hpm++) {
155  charge_start += (*it_hpm).charge;
156  counter++;
157  if (counter == fNHitsToDetermineStart) break;
158  }
159 
160  counter = 0;
161  for (HitPropertiesMultiset_t::reverse_iterator it_hpm = hpm.rbegin(); it_hpm != hpm.rend();
162  it_hpm++) {
163  charge_end += (*it_hpm).charge;
164  counter++;
165  if (counter == fNHitsToDetermineStart) break;
166  }
167 
168  return (charge_start > charge_end);
169 }
170 
172  recob::Track const& track,
173  size_t const& i_track,
174  std::vector<anab::Calorimetry>& caloVector,
175  std::vector<size_t>& assnTrackCaloVector,
176  geo::PlaneID const& planeID)
177 {
178  size_t n_hits = hpm.size();
179  std::vector<float> dEdxVector, dQdxVector, resRangeVector, deadWireVector, pitchVector;
180  std::vector<TVector3> XYZVector;
181 
182  dEdxVector.reserve(n_hits);
183  dQdxVector.reserve(n_hits);
184  resRangeVector.reserve(n_hits);
185  deadWireVector.reserve(n_hits);
186  pitchVector.reserve(n_hits);
187  XYZVector.reserve(n_hits);
188 
189  float kinetic_energy = 0, track_length = track.Length();
190  if (IsInvertedTrack(hpm)) {
191 
192  for (HitPropertiesMultiset_t::iterator it_hpm = hpm.begin(); it_hpm != hpm.end(); it_hpm++) {
193  dEdxVector.push_back((*it_hpm).dEdx);
194  dQdxVector.push_back((*it_hpm).dQdx);
195  resRangeVector.push_back((*it_hpm).path_fraction * track_length);
196  pitchVector.push_back((*it_hpm).pitch);
197  XYZVector.push_back((*it_hpm).xyz);
198  kinetic_energy += dEdxVector.back() * pitchVector.back();
199  }
200  }
201  else {
202 
203  for (HitPropertiesMultiset_t::reverse_iterator it_hpm = hpm.rbegin(); it_hpm != hpm.rend();
204  it_hpm++) {
205  dEdxVector.push_back((*it_hpm).dEdx);
206  dQdxVector.push_back((*it_hpm).dQdx);
207  resRangeVector.push_back((1 - (*it_hpm).path_fraction) * track_length);
208  pitchVector.push_back((*it_hpm).pitch);
209  XYZVector.push_back((*it_hpm).xyz);
210  kinetic_energy += dEdxVector.back() * pitchVector.back();
211  }
212  }
213 
214  caloVector.emplace_back(kinetic_energy,
215  dEdxVector,
216  dQdxVector,
217  resRangeVector,
218  deadWireVector,
219  track_length,
220  pitchVector,
222  planeID);
223  assnTrackCaloVector.emplace_back(i_track);
224 }
225 
227 {
228 
229  for (auto const& hit : hpm)
230  hit.Print();
231 
232  std::cout << "Inverted? " << IsInvertedTrack(hpm) << std::endl;
233 }
intermediate_table::iterator iterator
geo::GeometryCore const & geom
bool IsInvertedTrack(HitPropertiesMultiset_t const &)
Declaration of signal hit object.
Point_t const & LocationAtPoint(size_t i) const
Access to track position at different points.
Definition: Track.h:160
void MakeCalorimetryObject(HitPropertiesMultiset_t const &hpm, recob::Track const &track, size_t const &i_track, std::vector< anab::Calorimetry > &caloVector, std::vector< size_t > &assnTrackCaloVector, geo::PlaneID const &planeID)
void ExtractCalorimetry(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, std::vector< recob::Track > const &, std::vector< recob::Hit > const &, std::vector< std::vector< size_t >> const &, std::vector< anab::Calorimetry > &, std::vector< size_t > &, Providers_t providers)
The data type to uniquely identify a Plane.
Definition: geo_types.h:463
size_t NumberTrajectoryPoints() const
Various functions related to the presence and the number of (valid) points.
Definition: Track.h:136
bool operator()(std::pair< geo::WireID, float > const &i, std::pair< geo::WireID, float > const &j)
Provider const * get() const
Returns the provider with the specified type.
Definition: ProviderPack.h:182
float Integral() const
Integral under the calibrated signal waveform of the hit, in tick x ADC units.
Definition: Hit.h:244
geo::WireID const & WireID() const
Initial tdc tick for hit.
Definition: Hit.h:280
std::multiset< HitProperties, HitPropertySorter > HitPropertiesMultiset_t
dist_projected(recob::Hit const &h, geo::GeometryCore const &g)
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
std::vector< float > CreatePathLengthFractionVector(recob::Track const &track)
auto counter(T begin, T end)
Returns an object to iterate values from begin to end in a range-for loop.
Definition: counter.h:295
double Length(size_t p=0) const
Access to various track properties.
Definition: Track.h:207
parameter set interface
void PrintHitPropertiesMultiset(HitPropertiesMultiset_t const &hpm)
T get(std::string const &key) const
Definition: ParameterSet.h:314
tick_as<> tick
Tick number, represented by std::ptrdiff_t.
Definition: electronics.h:73
double ConvertXToTicks(double X, int p, int t, int c) const
Provides recob::Track data product.
std::vector< Point_t > convertCollToPoint(std::vector< Point > const &coll)
Definition: TrackingTypes.h:68
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:481
Description of geometry of one entire detector.
Definition: GeometryCore.h:119
Definition of data types for geometry description.
Detector simulation of raw signals on wires.
void AnalyzeHit(detinfo::DetectorClocksData const &, detinfo::DetectorPropertiesData const &, recob::Hit const &, recob::Track const &, std::vector< std::pair< geo::WireID, float >> const &, std::vector< float > const &, HitPropertiesMultiset_t &, geo::GeometryCore const &)
decltype(auto) get(T &&obj)
ADL-aware version of std::to_string.
Definition: StdUtils.h:120
TrackCalorimetryAlg(fhicl::ParameterSet const &p)
Contains all timing reference information for the detector.
Container for a list of pointers to providers.
Definition: ProviderPack.h:111
View_t View(PlaneID const &pid) const
Returns the view (wire orientation) on the channels of specified TPC plane.
double TrackPitchInView(recob::Track const &track, geo::View_t view, size_t trajectory_point=0U)
Returns the projected length of track on a wire pitch step [cm].
Definition: TrackUtils.cxx:75
double dEdx_AREA(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, recob::Hit const &hit, double pitch, double T0=0) const
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:46
recob::Hit const & hit
Float_t track
Definition: plot.C:35
Utility functions to extract information from recob::Track
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a "fitted" track:
Definition: Track.h:49