30 #include "larevt/CalibrationDBI/Interface/ChannelStatusProvider.h" 31 #include "larevt/CalibrationDBI/Interface/ChannelStatusService.h" 50 #include "cetlib/pow.h" 56 constexpr
unsigned int int_max_as_unsigned_int{std::numeric_limits<int>::max()};
57 typedef std::pair<unsigned, std::vector<unsigned>> OrganizedHits;
76 Comment(
"Module label for track producer.")};
79 Comment(
"Module label for T0 time producer."),
83 Name(
"AssocHitModuleLabel"),
84 Comment(
"Module label for association between tracks and hits. If not set, defaults to " 90 Comment(
"Method used to extract charge from a hit. Options: 0==Amplitude(), 1==Integral(), " 91 "2==SummedADC(), 3==SummedIntegral(). See the ChargeMethod enum.")};
94 Name(
"FieldDistortion"),
95 Comment(
"True if field distortion (i.e. from space charge) is included in the input.")};
98 Name(
"FieldDistortionEfield"),
99 Comment(
"True if field distortion (i.e. from space charge) is included in the input.")};
102 Name(
"TrackIsFieldDistortionCorrected"),
103 Comment(
"Whether the space-points on the input tracks have their points corrected for the " 104 "field distortions. " 105 "I.e. whether the track trajectory points represent charge as seen by wires or the " 106 "3D particle trajectory.")};
109 Comment(
"Which cryostat number the input tracks occupy.")};
112 Name(
"FieldDistortionCorrectionXSign"),
113 Comment(
"Sign of the field distortion correction to be applied in the X direction. " 114 "Positive by default."),
119 Comment(
"Configuration for the calo::CalorimetryAlg")};
123 Comment(
"List of INormalizeCharge tool configurations to use.")};
140 const std::vector<const recob::TrackHitMeta*>& thms,
145 const std::vector<const recob::TrackHitMeta*>& thms,
150 const std::vector<const recob::TrackHitMeta*>& thms,
166 const std::vector<recob::Hit const*>& sharedHits);
184 produces<std::vector<anab::Calorimetry>>();
185 produces<art::Assns<recob::Track, anab::Calorimetry>>();
186 std::vector<fhicl::ParameterSet> norm_tool_configs =
189 fNormTools.push_back(art::make_tool<INormalizeCharge>(p));
198 auto const det_prop =
201 size_t nplanes = geom->
Nplanes();
204 std::unique_ptr<std::vector<anab::Calorimetry>> outputCalo(
new std::vector<anab::Calorimetry>);
205 std::unique_ptr<art::Assns<recob::Track, anab::Calorimetry>> outputCaloAssn(
210 std::vector<art::Ptr<recob::Track>> tracklist;
228 for (
unsigned trk_i = 0; trk_i < tracklist.size(); trk_i++) {
232 const std::vector<art::Ptr<recob::Hit>>&
hits = fmHits.at(trk_i);
233 const std::vector<const recob::TrackHitMeta*>& thms = fmHits.
data(trk_i);
237 const std::vector<art::Ptr<anab::T0>>& this_t0s = fmT0s.at(trk_i);
238 if (this_t0s.size()) T0 = this_t0s.at(0)->Time();
242 std::vector<std::vector<OrganizedHits>> hit_indices =
OrganizeHits(hits, thms, track, nplanes);
244 for (
unsigned plane_i = 0; plane_i < nplanes; plane_i++) {
246 float kinetic_energy = 0.;
247 std::vector<float> dEdxs;
248 std::vector<float> dQdxs;
249 std::vector<float> resranges;
250 std::vector<float> deadwireresranges;
252 std::vector<float> pitches;
253 std::vector<geo::Point_t> xyzs;
254 std::vector<size_t> tp_indices;
258 plane.
Plane = plane_i;
263 std::vector<float> lengths;
264 for (
unsigned hit_i = 0; hit_i < hit_indices[plane_i].size(); hit_i++) {
265 unsigned hit_index = hit_indices[plane_i][hit_i].first;
267 std::vector<recob::Hit const*> sharedHits = {};
268 sharedHits.reserve(hit_indices[plane_i][hit_i].
second.size());
269 for (
const unsigned shared_hit_index : hit_indices[plane_i][hit_i].
second)
270 sharedHits.push_back(hits[shared_hit_index].
get());
276 double pitch =
GetPitch(track, hits[hit_index], thms[hit_index]);
279 double charge =
GetCharge(hits[hit_index], sharedHits);
282 double EField =
GetEfield(det_prop, track, hits[hit_index], thms[hit_index]);
286 double phi = acos(
abs(direction.x())) * 180 / M_PI;
288 double dQdx = charge / pitch;
303 hits[hit_index]->PeakTime(),
304 hits[hit_index]->WireID().Plane,
311 hits[hit_index]->PeakTime(),
312 hits[hit_index]->WireID().Plane,
318 if (xyzs.size() == 0) { lengths.push_back(0.); }
320 lengths.push_back((location - xyzs.back()).
r());
324 dEdxs.push_back(dEdx);
325 dQdxs.push_back(dQdx);
326 pitches.push_back(pitch);
327 xyzs.push_back(location);
328 kinetic_energy += dEdx * pitch;
337 tp_indices.push_back(hits[hit_index].key());
342 if (lengths.size() > 1) {
343 range = std::accumulate(lengths.begin(), lengths.end(), 0.);
352 resranges.resize(lengths.size());
354 resranges[lengths.size() - 1] = lengths.back() / 2.;
355 for (
int i_len = lengths.size() - 2; i_len >= 0; i_len--) {
356 resranges[i_len] = resranges[i_len + 1] + lengths[i_len + 1];
360 resranges[0] = lengths[1] / 2.;
361 for (
unsigned i_len = 1; i_len < lengths.size(); i_len++) {
362 resranges[i_len] = resranges[i_len - 1] + lengths[i_len];
370 if (lengths.size() > 1) {
383 outputCalo->push_back(
393 evt.
put(std::move(outputCalo));
394 evt.
put(std::move(outputCaloAssn));
401 const std::vector<const recob::TrackHitMeta*>& thms,
418 const std::vector<const recob::TrackHitMeta*>& thms,
422 std::vector<std::vector<OrganizedHits>> ret(nplanes);
423 for (
unsigned i = 0; i <
hits.size(); i++) {
424 if (
HitIsValid(thms[i], track)) { ret[
hits[i]->WireID().Plane].push_back({i, {}}); }
432 const std::vector<const recob::TrackHitMeta*>& thms,
441 struct HitIdentifier {
456 inline bool operator==(
const HitIdentifier& rhs)
const 458 return startTick == rhs.startTick && endTick == rhs.endTick && wire == rhs.wire;
462 inline bool operator>(
const HitIdentifier& rhs)
const {
return integral > rhs.integral; }
465 std::vector<std::vector<OrganizedHits>> ret(nplanes);
466 std::vector<std::vector<HitIdentifier>> hit_idents(nplanes);
467 for (
unsigned i = 0; i <
hits.size(); i++) {
469 HitIdentifier this_ident(*
hits[i]);
472 bool found_snippet =
false;
473 auto const plane =
hits[i]->WireID().Plane;
474 for (
unsigned j = 0; j < ret[plane].size(); j++) {
475 if (this_ident == hit_idents[plane][j]) {
476 found_snippet =
true;
477 if (this_ident > hit_idents[plane][j]) {
478 ret[plane][j].second.push_back(ret[plane][j].first);
479 ret[plane][j].first = i;
480 hit_idents[plane][j] = this_ident;
483 ret[plane][j].second.push_back(i);
488 if (!found_snippet) {
489 ret[plane].push_back({i, {}});
490 hit_idents[plane].push_back(this_ident);
499 if (thm->
Index() == int_max_as_unsigned_int)
return false;
515 auto const*
sce = lar::providerFrom<spacecharge::SpaceChargeService>();
523 ret.SetY(ret.Y() + offset.Y());
524 ret.SetZ(ret.Z() + offset.Z());
542 auto const*
sce = lar::providerFrom<spacecharge::SpaceChargeService>();
554 ret.SetY(ret.Y() + offset.Y());
555 ret.SetZ(ret.Z() + offset.Z());
566 auto const*
sce = lar::providerFrom<spacecharge::SpaceChargeService>();
589 dir = (loc_pdx - loc_mdx) / (loc_mdx - loc_pdx).r();
597 double cosgamma =
std::abs(std::sin(angleToVert) * dir.Y() + std::cos(angleToVert) * dir.Z());
599 if (cosgamma) { pitch = geom->
WirePitch(hit->
View()) / cosgamma; }
610 pitch = (locw_pdx_traj - loc).R();
616 const std::vector<recob::Hit const*>& sharedHits)
623 return std::accumulate(
627 [](
double sum,
recob::Hit const* sharedHit) {
return sum + sharedHit->Integral(); });
642 ret = nt->Normalize(ret, e, h, location, direction, t0);
653 auto const*
sce = lar::providerFrom<spacecharge::SpaceChargeService>();
655 double EField = dprop.
Efield();
662 EFieldOffsets = EField * EFieldOffsets;
664 EField = EFieldOffsets.r();
code to link reconstructed objects back to the MC truth information
std::vector< std::vector< OrganizedHits > > OrganizeHitsIndividual(const std::vector< art::Ptr< recob::Hit >> &hits, const std::vector< const recob::TrackHitMeta * > &thms, const recob::Track &track, unsigned nplanes)
GnocchiCalorimetry(Parameters const &config)
Functions to help with numbers.
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
fhicl::Atom< unsigned > Cryostat
void produce(art::Event &evt) override
constexpr bool operator>(Interval< Q, Cat > const a, Quantity< Args... > const b) noexcept
double GetCharge(const art::Ptr< recob::Hit > hit, const std::vector< recob::Hit const * > &sharedHits)
const recob::TrackTrajectory & Trajectory() const
Access to the stored recob::TrackTrajectory.
Declaration of signal hit object.
fhicl::DelegatedParameter NormTools
fhicl::Atom< float > FieldDistortionCorrectionXSign
Point_t const & LocationAtPoint(size_t i) const
Access to track position at different points.
EDProducer(fhicl::ParameterSet const &pset)
The data type to uniquely identify a Plane.
bool isValid
Whether this ID points to a valid element.
bool HasValidPoint(size_t i) const
Various functions related to the presence and the number of (valid) points.
constexpr auto abs(T v)
Returns the absolute value of the argument.
CryostatID_t Cryostat
Index of cryostat.
float Integral() const
Integral under the calibrated signal waveform of the hit, in tick x ADC units.
geo::View_t View() const
View for the plane of the hit.
WireID_t Wire
Index of the wire within its plane.
geo::WireID const & WireID() const
Initial tdc tick for hit.
data_const_reference data(size_type i) const
double Normalize(double dQdx, const art::Event &e, const recob::Hit &h, const geo::Point_t &location, const geo::Vector_t &direction, double t0)
geo::Point_t GetLocationAtWires(const recob::Track &track, const art::Ptr< recob::Hit > hit, const recob::TrackHitMeta *meta)
float PeakAmplitude() const
The estimated amplitude of the hit at its peak, in ADC units.
double Efield(unsigned int planegap=0) const
kV/cm
fhicl::Atom< std::string > TrackModuleLabel
TPCGeo const & TPC(TPCID const &tpcid=tpc_zero) const
Returns the specified TPC.
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
fhicl::Atom< std::string > T0ModuleLabel
double WireAngleToVertical(View_t view, TPCID const &tpcid) const
Returns the angle of the wires in the specified view from vertical.
double dEdx_AMP(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, recob::Hit const &hit, double pitch, double T0=0) const
#define DEFINE_ART_MODULE(klass)
std::vector< std::vector< OrganizedHits > > OrganizeHitsSnippets(const std::vector< art::Ptr< recob::Hit >> &hits, const std::vector< const recob::TrackHitMeta * > &thms, const recob::Track &track, unsigned nplanes)
std::vector< std::unique_ptr< INormalizeCharge > > fNormTools
double GetEfield(const detinfo::DetectorPropertiesData &dprop, const recob::Track &track, const art::Ptr< recob::Hit > hit, const recob::TrackHitMeta *meta)
fhicl::Atom< bool > FieldDistortion
Provides recob::Track data product.
geo::Point_t WireToTrajectoryPosition(const geo::Point_t &loc, const geo::TPCID &tpc)
constexpr PlaneID const & asPlaneID() const
Conversion to PlaneID (for convenience of notation).
The data type to uniquely identify a TPC.
PlaneID_t Plane
Index of the plane within its TPC.
raw::TDCtick_t StartTick() const
Initial tdc tick for hit.
float dEdx(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const TCSlice &slc, TP3D &tp3d)
fhicl::Atom< bool > TrackIsFieldDistortionCorrected
Detector simulation of raw signals on wires.
Encapsulate the geometry of a wire.
raw::TDCtick_t EndTick() const
Final tdc tick for hit.
fhicl::Table< calo::CalorimetryAlg::Config > CalorimetryAlgConfig
bool CreateAssn(art::Event &evt, std::vector< T > const &a, art::Ptr< U > const &b, art::Assns< U, T > &assn, std::string a_instance, size_t index=UINT_MAX)
Creates a single one-to-one association.
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Utility object to perform functions of association.
Encapsulate the construction of a single detector plane.
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)
geo::Point_t TrajectoryToWirePosition(const geo::Point_t &loc, const geo::TPCID &tpc)
float SummedADC() const
The sum of calibrated ADC counts of the hit (0. by default)
fhicl::Atom< bool > FieldDistortionEfield
unsigned int Nplanes(TPCID const &tpcid=tpc_zero) const
Returns the total number of planes in the specified TPC.
Vector_t DirectionAtPoint(size_t i) const
Access to track direction at different points.
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)
2D representation of charge deposited in the TDC/wire plane
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
TPCID_t TPC
Index of the TPC within its cryostat.
Collection of Physical constants used in LArSoft.
second_as<> second
Type of time stored in seconds, in double precision.
fhicl::Atom< std::string > AssocHitModuleLabel
Vector_t DriftDir() const
Returns the direction of the drift (vector pointing toward the planes).
Point_t const & Start() const
Returns the position of the first valid point of the trajectory [cm].
bool HitIsValid(const recob::TrackHitMeta *thm, const recob::Track &track)
Utility functions to extract information from recob::Track
bool operator==(infinite_endcount_iterator< T > const &, count_iterator< T > const &)
std::vector< std::vector< OrganizedHits > > OrganizeHits(const std::vector< art::Ptr< recob::Hit >> &hits, const std::vector< const recob::TrackHitMeta * > &thms, const recob::Track &track, unsigned nplanes)
Length_t WirePitch(PlaneID const &planeid=plane_zero) const
Returns the distance between two consecutive wires.
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a "fitted" track:
This is an interface for an art Tool which scales charge by some factor given information about its a...