30 , fG4ModuleLabel(config.G4ModuleLabel())
31 , fSimChannelModuleLabel(config.SimChannelModuleLabel())
32 , fHitLabel(config.DefaultHitModuleLabel())
33 , fMinHitEnergyFraction(config.MinHitEnergyFraction())
34 , fOverrideRealData(config.OverrideRealData())
35 , fHitTimeRMS(config.HitTimeRMS())
63 std::vector<const sim::IDE*> ideps;
67 for (
auto mapitr = tdcidemap.begin(); mapitr != tdcidemap.end(); mapitr++) {
69 const std::vector<sim::IDE>& idevec =
72 for (
size_t iv = 0; iv < idevec.size(); ++iv) {
76 if (
abs(idevec[iv].trackID) == id) ideps.push_back(&(idevec[iv]));
87 std::vector<const sim::IDE*> ide_Ps;
92 for (
const auto& item :
sc->TDCIDEMap()) {
95 for (
const sim::IDE& ide : item.second) {
96 if (
abs(ide.
trackID) == id) ide_Ps.push_back(&ide);
111 return (a->Channel() < channel);
113 return ((ilb !=
fSimChannels.end()) && ((*ilb)->Channel() == channel)) ?
122 throw cet::exception(
"BackTracker") <<
"No sim::SimChannel corresponding " 123 <<
"to channel: " << channel <<
"\n";
130 const double hit_start_time,
131 const double hit_end_time)
const 134 if (!schannel)
return {};
136 std::vector<sim::TrackIDE> trackIDEs;
141 int start_tdc = clockData.
TPCTick2TDC(hit_start_time);
143 if (start_tdc < 0) start_tdc = 0;
144 if (end_tdc < 0) end_tdc = 0;
149 for (
size_t e = 0;
e < simides.size(); ++
e)
153 if (totalE < 1.
e-5) totalE = 1.;
157 for (
size_t e = 0;
e < simides.size(); ++
e) {
164 info.
energy = simides[
e].energy;
167 trackIDEs.push_back(info);
187 std::vector<int> retVec;
188 for (
auto const trackIDE : this->
HitToTrackIDEs(clockData, hit)) {
189 retVec.push_back(trackIDE.trackID);
200 std::vector<sim::TrackIDE> eveIDEs;
201 std::vector<sim::TrackIDE> trackIDEs = this->
HitToTrackIDEs(clockData, hit);
202 std::map<int, std::pair<double, double>> eveToEMap;
204 for (
const auto& trackIDE : trackIDEs) {
206 std::make_pair(trackIDE.energy, trackIDE.numElectrons));
207 if (check.second ==
false) {
208 check.first->second.first += trackIDE.energy;
209 check.first->second.second += trackIDE.numElectrons;
213 totalE += trackIDE.energy;
215 eveIDEs.reserve(eveToEMap.size());
216 for (
const auto& eveToE : eveToEMap) {
219 eveTrackIDE_tmp.
trackID = eveToE.first;
220 eveTrackIDE_tmp.
energy = eveToE.second.first;
224 eveIDEs.push_back(eveTrackIDE_tmp);
240 std::vector<art::Ptr<recob::Hit>> hitList;
241 std::vector<sim::TrackIDE> trackIDE;
242 for (
auto itr = hitsIn.begin(); itr != hitsIn.end(); ++itr) {
249 for (
auto itr_trakIDE = trackIDE.begin(); itr_trakIDE != trackIDE.end(); ++itr_trakIDE) {
251 hitList.push_back(hit);
261 std::vector<int>
const& tkIds,
269 std::vector<std::pair<int, art::Ptr<recob::Hit>>> hitList;
270 std::vector<sim::TrackIDE> tids;
271 for (
auto itr = hitsIn.begin(); itr != hitsIn.end(); ++itr) {
278 for (
auto itid = tids.begin(); itid != tids.end(); ++itid) {
279 for (
auto itkid = tkIds.begin(); itkid != tkIds.end(); ++itkid) {
280 if (itid->trackID == *itkid) {
282 hitList.push_back(std::make_pair(*itkid, hit));
289 std::vector<std::vector<art::Ptr<recob::Hit>>> truHits;
291 std::vector<art::Ptr<recob::Hit>> tmpHits;
292 for (
auto itkid = tkIds.begin(); itkid != tkIds.end(); ++itkid) {
294 for (
auto itr = hitList.begin(); itr != hitList.end(); ++itr) {
295 if (*itkid == (*itr).first) tmpHits.push_back((*itr).second);
297 truHits.push_back(tmpHits);
311 if (start_tdc < 0) start_tdc = 0;
312 if (end_tdc < 0) end_tdc = 0;
322 std::vector<const sim::IDE*> retVec;
325 if (start_tdc < 0) start_tdc = 0;
326 if (end_tdc < 0) end_tdc = 0;
328 if (start_tdc > end_tdc) {
throw; }
330 const std::vector<std::pair<unsigned short, std::vector<sim::IDE>>>& tdcIDEMap =
334 std::vector<const std::pair<unsigned short, std::vector<sim::IDE>>*> tdcIDEMap_SortedPointers;
335 for (
auto& pair : tdcIDEMap) {
336 tdcIDEMap_SortedPointers.push_back(&pair);
345 auto pairSort = [](
auto& a,
auto& b) {
return a->first < b->first; };
347 tdcIDEMap_SortedPointers.begin(), tdcIDEMap_SortedPointers.end(), pairSort)) {
348 std::sort(tdcIDEMap_SortedPointers.begin(), tdcIDEMap_SortedPointers.end(), pairSort);
351 std::vector<sim::IDE> dummyVec;
353 std::pair<double, std::vector<sim::IDE>> start_tdcPair =
354 std::make_pair(start_tdc, dummyVec);
356 std::pair<double, std::vector<sim::IDE>> end_tdcPair = std::make_pair(end_tdc, dummyVec);
357 auto start_tdcPair_P = &start_tdcPair;
358 auto end_tdcPair_P = &end_tdcPair;
359 auto mapFirst = std::lower_bound(
360 tdcIDEMap_SortedPointers.begin(), tdcIDEMap_SortedPointers.end(), start_tdcPair_P, pairSort);
363 auto mapLast = std::upper_bound(
364 tdcIDEMap_SortedPointers.begin(), tdcIDEMap_SortedPointers.end(), end_tdcPair_P, pairSort);
365 for (
auto& mapitr = mapFirst; mapitr != mapLast; ++mapitr) {
366 for (
auto& ide : (*mapitr)->second) {
367 retVec.push_back(&ide);
376 std::vector<double> xyz(3, 0.0);
378 for (
auto const& ide : ides) {
379 double weight = ide.numElectrons;
381 xyz[0] += (weight * ide.x);
382 xyz[1] += (weight * ide.y);
383 xyz[2] += (weight * ide.z);
386 throw cet::exception(
"BackTracker") <<
"No sim::IDEs providing non-zero number of electrons" 387 <<
" can't determine originating location from truth\n";
397 std::vector<sim::IDE> ides;
398 for (
auto ide_P : ide_Ps) {
399 ides.push_back(*ide_P);
408 std::vector<const sim::IDE*> ide_Ps = this->
HitToSimIDEs_Ps(clockData, hit);
414 std::set<int>
const& trackIds,
420 for (
const auto& tIDE : hitTrackIDEs) {
421 if (trackIds.find(tIDE.trackID) != trackIds.end()) {
427 if (hits.size() > 0) {
return double(
double(desired) /
double(hits.size())); }
433 std::set<int>
const& trackIds,
436 double totalCharge = 0., desired = 0.;
438 totalCharge +=
hit->Integral();
440 for (
const auto& trackIDE : trackIDEs) {
441 if (trackIds.find(trackIDE.trackID) != trackIds.end()) {
442 desired +=
hit->Integral();
447 if (totalCharge > 0.0) {
return (desired / totalCharge); }
453 std::set<int>
const& trackIds,
459 int desired = 0, total = 0;
463 for (
const auto& trackIDE : hitTrackIDEs) {
464 if (trackIds.find(trackIDE.trackID) != trackIds.end() &&
472 for (
const auto&
hit : allHits) {
477 for (
const auto& hitIDE : hitTrackIDEs) {
478 if (trackIds.find(hitIDE.trackID) != trackIds.end() &&
485 if (total >= 0) {
return double(
double(desired) /
double(total)); }
492 std::set<int>
const& trackIds,
497 double desired = 0., total = 0.;
500 for (
const auto& hitIDE : hitTrackIDEs) {
501 if (trackIds.find(hitIDE.trackID) != trackIds.end() &&
503 desired +=
hit->Integral();
509 for (
const auto&
hit : allHits) {
510 if (
hit->View() != view && view !=
geo::k3D) {
continue; }
512 for (
const auto& hitIDE : hitTrackIDEs) {
513 if (trackIds.find(hitIDE.trackID) != trackIds.end() &&
515 total +=
hit->Integral();
521 if (total > 0.) {
return desired / total; }
533 std::vector<sim::TrackIDE> trackIDEs =
535 for (
const auto& ide : trackIDEs) {
536 tids.insert(ide.trackID);
546 std::set<int> eveIds;
549 for (
const auto& ide : ides) {
550 eveIds.insert(ide.trackID);
562 std::vector<std::vector<std::vector<int>>> numHits(
fGeom->
Ncryostats());
563 std::vector<std::vector<std::vector<double>>> hitWeight(
fGeom->
Ncryostats());
564 std::vector<std::vector<std::vector<std::vector<double>>>> hitPos(
fGeom->
Ncryostats());
566 for (
size_t c = 0; c < numHits.size(); ++c) {
571 for (
size_t t = 0; t < numHits[c].size(); ++t) {
579 for (
auto const& hit_ptr :
hits) {
584 std::vector<double> hitOrigin = this->
HitToXYZ(clockData, hit_ptr);
585 geo::Point_t const worldLoc{hitOrigin[0], hitOrigin[1], hitOrigin[2]};
588 auto const [cstat, tpc] = std::make_tuple(tpcid.Cryostat, tpcid.TPC);
593 hitPos[cstat][tpc][hit.
WireID().
Plane] = hitOrigin;
600 std::vector<double> xyz(3);
601 for (
size_t c = 0; c < numHits.size(); ++c) {
602 for (
size_t t = 0; t < numHits[c].size(); ++t) {
603 for (
size_t p = 0; p < numHits[c][t].size(); ++p) {
605 if (numHits[c][t][p] == 1) {
607 xyz[0] += hitPos[c][t][p][0];
608 xyz[1] += hitPos[c][t][p][1];
609 xyz[2] += hitPos[c][t][p][2];
619 <<
"No hits to determine originating location from truth\n";
TrackID_t trackID
Geant4 supplied track ID.
const art::InputTag fHitLabel
unsigned int NTPC(CryostatID const &cryoid=cryostat_zero) const
Returns the total number of TPCs in the specified cryostat.
const geo::GeometryCore * fGeom
std::vector< art::Ptr< sim::SimChannel > > fSimChannels
std::vector< std::vector< art::Ptr< recob::Hit > > > TrackIdsToHits_Ps(detinfo::DetectorClocksData const &clockData, std::vector< int > const &tkIds, std::vector< art::Ptr< recob::Hit >> const &hitsIn) const
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
art::Ptr< sim::SimChannel > FindSimChannelPtr(raw::ChannelID_t channel) const
Returns the cached sim::SimChannel on the specified channel.
std::vector< sim::TrackIDE > HitToEveTrackIDEs(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
BackTracker(const fhiclConfig &config, const cheat::ParticleInventory *partInv, const geo::GeometryCore *geom)
const art::InputTag fSimChannelModuleLabel
float numElectrons
number of electrons from the particle detected on the wires
constexpr auto abs(T v)
Returns the absolute value of the argument.
const double fMinHitEnergyFraction
CryostatID_t Cryostat
Index of cryostat.
float Integral() const
Integral under the calibrated signal waveform of the hit, in tick x ADC units.
geo::WireID const & WireID() const
Initial tdc tick for hit.
float energy
energy from the particle with this trackID [MeV]
unsigned int Ncryostats() const
Returns the number of cryostats in the detector.
pure virtual base interface for detector clocks
std::vector< sim::TrackIDE > HitToTrackIDEs(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
3-dimensional objects, potentially hits, clusters, prongs, etc.
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Access the description of detector geometry.
std::vector< double > HitToXYZ(detinfo::DetectorClocksData const &clockData, const recob::Hit &hit) const
Ionization at a point of the TPC sensitive volume.
static const int NoParticleId
art::Ptr< sim::SimChannel > FindSimChannel(raw::ChannelID_t channel) const
Returns the cached sim::SimChannel on the specified channel.
std::vector< art::Ptr< recob::Hit > > TrackIdToHits_Ps(detinfo::DetectorClocksData const &clockData, int tkId, std::vector< art::Ptr< recob::Hit >> const &hitsIn) const
std::set< int > GetSetOfTrackIds() const
double HitCollectionPurity(detinfo::DetectorClocksData const &clockData, std::set< int > const &trackIds, std::vector< art::Ptr< recob::Hit >> const &hits) const
float energyFrac
fraction of hit energy from the particle with this trackID
const cheat::ParticleInventory * fPartInv
The data type to uniquely identify a TPC.
PlaneID_t Plane
Index of the plane within its TPC.
Description of geometry of one entire detector.
Definition of data types for geometry description.
float PeakTimeMinusRMS(float sigmas=+1.) const
Returns a time sigmas RMS away from the peak time.
Detector simulation of raw signals on wires.
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
std::vector< sim::TrackIDE > ChannelToTrackIDEs(detinfo::DetectorClocksData const &clockData, raw::ChannelID_t channel, const double hit_start_time, const double hit_end_time) const
std::set< int > GetSetOfEveIds() const
double HitChargeCollectionPurity(detinfo::DetectorClocksData const &clockData, std::set< int > const &trackIds, std::vector< art::Ptr< recob::Hit >> const &hits) const
decltype(auto) get(T &&obj)
ADL-aware version of std::to_string.
code to link reconstructed objects back to the MC truth information
int trackID
Geant4 supplied trackID.
const bool fOverrideRealData
Contains all timing reference information for the detector.
std::vector< sim::IDE > TrackIDsAndEnergies(TDC_t startTDC, TDC_t endTDC) const
Return all the recorded energy deposition within a time interval.
std::vector< int > HitToTrackIds(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
View_t View(PlaneID const &pid) const
Returns the view (wire orientation) on the channels of specified TPC plane.
double HitCollectionEfficiency(detinfo::DetectorClocksData const &clockData, std::set< int > const &trackIds, std::vector< art::Ptr< recob::Hit >> const &hits, std::vector< art::Ptr< recob::Hit >> const &allhits, geo::View_t const &view) const
std::vector< const sim::IDE * > TrackIdToSimIDEs_Ps(int const &id) const
const art::InputTag fG4ModuleLabel
unsigned int Nplanes(TPCID const &tpcid=tpc_zero) const
Returns the total number of planes in the specified TPC.
2D representation of charge deposited in the TDC/wire plane
float PeakTimePlusRMS(float sigmas=+1.) const
Returns a time sigmas RMS away from the peak time.
unsigned int ChannelID_t
Type representing the ID of a readout channel.
TPCID_t TPC
Index of the TPC within its cryostat.
std::vector< sim::IDE > HitToAvgSimIDEs(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
int TrackIdToEveTrackId(const int &tid) const
Ionization energy from a Geant4 track.
raw::ChannelID_t Channel() const
ID of the readout channel the hit was extracted from.
std::vector< double > SimIDEsToXYZ(std::vector< sim::IDE > const &ides) const
double HitChargeCollectionEfficiency(detinfo::DetectorClocksData const &clockData, std::set< int > const &trackIds, std::vector< art::Ptr< recob::Hit >> const &hits, std::vector< art::Ptr< recob::Hit >> const &allhits, geo::View_t const &view) const
double TPCTick2TDC(double const tick) const
Tools and modules for checking out the basics of the Monte Carlo.
TPCID PositionToTPCID(Point_t const &point) const
Returns the ID of the TPC at specified location.
cet::coded_exception< error, detail::translate > exception
The data type to uniquely identify a cryostat.
std::vector< const sim::IDE * > HitToSimIDEs_Ps(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
std::vector< double > SpacePointHitsToWeightedXYZ(detinfo::DetectorClocksData const &clockData, std::vector< art::Ptr< recob::Hit >> const &hits) const