174 auto const det_prop =
177 size_t nplanes = geom->
Nplanes();
180 std::unique_ptr< std::vector<anab::Calorimetry> > outputCalo(
new std::vector<anab::Calorimetry>);
185 std::vector<art::Ptr<recob::Track> > tracklist;
192 art::FindManyP<recob::Hit, recob::TrackHitMeta> fmHits(trackListHandle, evt, hitLabel);
198 for (
unsigned trk_i = 0; trk_i < tracklist.size(); trk_i++) {
202 const std::vector<art::Ptr<recob::Hit>> &hits = fmHits.at(trk_i);
203 const std::vector<const recob::TrackHitMeta *> &thms = fmHits.data(trk_i);
207 const std::vector<art::Ptr<anab::T0>> &this_t0s = fmT0s.at(trk_i);
208 if (this_t0s.size()) T0 = this_t0s.at(0)->Time();
212 std::vector<std::vector<unsigned>> hit_indices =
OrganizeHits(hits, thms, track, nplanes);
214 for (
unsigned plane_i = 0; plane_i < nplanes; plane_i++) {
216 float kinetic_energy = 0.;
217 std::vector<float> dEdxs;
218 std::vector<float> dQdxs;
219 std::vector<float> resranges;
220 std::vector<float> deadwireresranges;
222 std::vector<float> pitches;
223 std::vector<geo::Point_t> xyzs;
224 std::vector<size_t> tp_indices;
228 plane.
Plane = plane_i;
233 std::vector<float> lengths;
234 for (
unsigned hit_i = 0; hit_i < hit_indices[plane_i].size(); hit_i++) {
235 unsigned hit_index = hit_indices[plane_i][hit_i];
241 double pitch =
GetPitch(track, hits[hit_index], thms[hit_index]);
244 double charge =
GetCharge(hits[hit_index]);
247 double EField =
GetEfield(det_prop, track, hits[hit_index], thms[hit_index]);
249 double dQdx = charge / pitch;
253 fCaloAlg.
dEdx_AMP(clock_data, det_prop, dQdx, hits[hit_index]->PeakTime(), hits[hit_index]->WireID().Plane, T0, EField) : \
254 fCaloAlg.
dEdx_AREA(clock_data, det_prop, dQdx, hits[hit_index]->PeakTime(), hits[hit_index]->WireID().Plane, T0, EField);
257 if (xyzs.size() == 0) {
258 lengths.push_back(0.);
261 lengths.push_back((location - xyzs.back()).
r());
265 dEdxs.push_back(dEdx);
266 dQdxs.push_back(dQdx);
267 pitches.push_back(pitch);
268 xyzs.push_back(location);
269 kinetic_energy += dEdx * pitch;
278 tp_indices.push_back(hits[hit_index].
key());
285 if (lengths.size() > 1) {
286 range = std::accumulate(lengths.begin(), lengths.end(), 0.);
290 bool is_downstream = \
294 resranges.resize(lengths.size());
296 resranges[lengths.size() - 1] = lengths.back() / 2.;
297 for (
int i_len = lengths.size() - 2; i_len >= 0; i_len --) {
298 resranges[i_len] = resranges[i_len+1] + lengths[i_len+1];
302 resranges[0] = lengths[1] / 2.;
303 for (
unsigned i_len = 1; i_len < lengths.size(); i_len ++) {
304 resranges[i_len] = resranges[i_len-1] + lengths[i_len];
312 if (lengths.size() > 1) {
337 util::CreateAssn(*
this, evt, *outputCalo, tracklist[trk_i], *outputCaloAssn);
fhicl::Atom< unsigned > Cryostat
const recob::TrackTrajectory & Trajectory() const
Access to the stored recob::TrackTrajectory.
The data type to uniquely identify a Plane.
bool isValid
Whether this ID points to a valid element.
CryostatID_t Cryostat
Index of cryostat.
fhicl::Atom< std::string > TrackModuleLabel
fhicl::Atom< std::string > T0ModuleLabel
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
double dEdx(float dqdx, float Efield)
unsigned int Nplanes(unsigned int tpc=0, unsigned int cstat=0) const
Returns the total number of wire planes in the specified TPC.
double dEdx_AMP(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, recob::Hit const &hit, double pitch, double T0=0) const
double GetEfield(const detinfo::DetectorPropertiesData &dprop, const recob::Track &track, const art::Ptr< recob::Hit > hit, const recob::TrackHitMeta *meta)
ProductID put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
bool CreateAssn(PRODUCER const &prod, art::Event &evt, std::vector< T > const &a, art::Ptr< U > const &b, art::Assns< U, T > &assn, std::string a_instance, size_t indx=UINT_MAX)
Creates a single one-to-one association.
PlaneID_t Plane
Index of the plane within its TPC.
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
std::vector< std::vector< unsigned > > OrganizeHits(const std::vector< art::Ptr< recob::Hit >> &hits, const std::vector< const recob::TrackHitMeta * > &thms, const recob::Track &track, unsigned nplanes)
Point_t const & End() const
Returns the position of the last valid point of the trajectory [cm].
double dEdx_AREA(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, recob::Hit const &hit, double pitch, double T0=0) const
double GetPitch(const recob::Track &track, const art::Ptr< recob::Hit > hit, const recob::TrackHitMeta *meta)
double GetCharge(const art::Ptr< recob::Hit > hit)
constexpr double kBogusD
obviously bogus double value
geo::Point_t GetLocation(const recob::Track &track, const art::Ptr< recob::Hit > hit, const recob::TrackHitMeta *meta)
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
TPCID_t TPC
Index of the TPC within its cryostat.
fhicl::Atom< std::string > AssocHitModuleLabel
Point_t const & Start() const
Returns the position of the first valid point of the trajectory [cm].
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a "fitted" track: