LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
sim::MCRecoEdep Class Reference

#include "MCRecoEdep.h"

Public Member Functions

 MCRecoEdep (fhicl::ParameterSet const &pset)
 Default constructor with fhicl parameters. More...
 
void MakeMCEdep (const std::vector< sim::SimChannel > &schArray)
 
void MakeMCEdep (const std::vector< sim::SimEnergyDeposit > &sedArray)
 
void MakeMCEdep (const std::vector< sim::SimEnergyDepositLite > &sedArray)
 
bool ExistTrack (const unsigned int track_id) const
 
int TrackToEdepIndex (unsigned int track_id) const
 Converts a track ID to MCEdep array index. Returns -1 if no corresponding array found . More...
 
const std::vector< sim::MCEdep > & GetEdepArrayAt (size_t edep_index) const
 Returns a vector of MCEdep object at the given index. More...
 
const std::map< unsigned int, size_t > TrackIndexMap () const
 Returns a map of track id <-> MCEdep vector index. More...
 
void Clear ()
 

Protected Member Functions

std::vector< sim::MCEdep > & __GetEdepArray__ (unsigned int track_id)
 

Protected Attributes

bool _debug_mode
 
bool _save_mchit
 
std::map< unsigned int, size_t > _track_index
 
std::vector< std::vector< sim::MCEdep > > _mc_edeps
 

Detailed Description

Definition at line 94 of file MCRecoEdep.h.

Constructor & Destructor Documentation

sim::MCRecoEdep::MCRecoEdep ( fhicl::ParameterSet const &  pset)

Default constructor with fhicl parameters.

Definition at line 40 of file MCRecoEdep.cxx.

References fhicl::ParameterSet::get().

42  {
43  _debug_mode = pset.get<bool>("DebugMode");
44  _save_mchit = pset.get<bool>("SaveMCHit");
45  }

Member Function Documentation

std::vector< sim::MCEdep > & sim::MCRecoEdep::__GetEdepArray__ ( unsigned int  track_id)
protected

Definition at line 54 of file MCRecoEdep.cxx.

55  {
56  if (ExistTrack(track_id)) return _mc_edeps.at((*_track_index.find(track_id)).second);
57  _track_index.insert(std::pair<unsigned int, size_t>(track_id, _mc_edeps.size()));
58  _mc_edeps.push_back(std::vector<sim::MCEdep>());
59  return (*(_mc_edeps.rbegin()));
60  }
bool ExistTrack(const unsigned int track_id) const
Definition: MCRecoEdep.h:107
std::map< unsigned int, size_t > _track_index
Definition: MCRecoEdep.h:138
std::vector< std::vector< sim::MCEdep > > _mc_edeps
Definition: MCRecoEdep.h:139
second_as<> second
Type of time stored in seconds, in double precision.
Definition: spacetime.h:82
void sim::MCRecoEdep::Clear ( )
inline

Definition at line 125 of file MCRecoEdep.h.

Referenced by MCReco::produce().

126  {
127  _mc_edeps.clear();
128  _track_index.clear();
129  std::vector<std::vector<sim::MCEdep>>().swap(_mc_edeps);
130  std::map<unsigned int, size_t>().swap(_track_index);
131  }
std::map< unsigned int, size_t > _track_index
Definition: MCRecoEdep.h:138
std::vector< std::vector< sim::MCEdep > > _mc_edeps
Definition: MCRecoEdep.h:139
void swap(lar::deep_const_fwd_iterator_nested< CITER, INNERCONTEXTRACT > &a, lar::deep_const_fwd_iterator_nested< CITER, INNERCONTEXTRACT > &b)
bool sim::MCRecoEdep::ExistTrack ( const unsigned int  track_id) const
inline

Definition at line 107 of file MCRecoEdep.h.

108  {
109  return (_track_index.find(track_id) != _track_index.end());
110  }
std::map< unsigned int, size_t > _track_index
Definition: MCRecoEdep.h:138
const std::vector< sim::MCEdep > & sim::MCRecoEdep::GetEdepArrayAt ( size_t  edep_index) const

Returns a vector of MCEdep object at the given index.

Definition at line 47 of file MCRecoEdep.cxx.

Referenced by sim::MCTrackRecoAlg::Reconstruct(), and sim::MCShowerRecoAlg::Reconstruct().

48  {
49  if (edep_index >= _mc_edeps.size())
50  throw cet::exception(__FUNCTION__) << Form("Track ID %zu not found!", edep_index);
51  return _mc_edeps.at(edep_index);
52  }
std::vector< std::vector< sim::MCEdep > > _mc_edeps
Definition: MCRecoEdep.h:139
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
void sim::MCRecoEdep::MakeMCEdep ( const std::vector< sim::SimChannel > &  schArray)

Definition at line 62 of file MCRecoEdep.cxx.

References geo::GeometryCore::ChannelToWire(), sim::details::createPlaneIndexMap(), sim::MCEdep::deps, and edep.

Referenced by MCReco::MakeMCEdep().

63  {
64  _mc_edeps.clear();
65  _track_index.clear();
66 
68 
69  // Key map to identify a unique particle energy deposition point
70  std::map<std::pair<UniquePosition, unsigned int>, int> hit_index_m;
71 
72  auto pindex = details::createPlaneIndexMap();
73 
74  if (_debug_mode) std::cout << "Processing " << schArray.size() << " channels..." << std::endl;
75  // Loop over channels
76  for (size_t i = 0; i < schArray.size(); ++i) {
77 
78  // Get data to loop over
79  auto const& sch = schArray[i];
80  const auto& sch_map(sch.TDCIDEMap());
81  // Channel
82  UInt_t ch = sch.Channel();
83  // Loop over ticks
84  for (auto tdc_iter = sch_map.begin(); tdc_iter != sch_map.end(); ++tdc_iter) {
85  // Loop over IDEs
86  for (auto const& ide : (*tdc_iter).second) {
87 
88  int track_id = ide.trackID;
89  if (track_id < 0) track_id = track_id * (-1);
90  unsigned int real_track_id = track_id;
91 
92  UniquePosition pos(ide.x, ide.y, ide.z);
93 
94  int hit_index = -1;
95  auto key = std::make_pair(pos, real_track_id);
96  auto hit_index_track_iter = hit_index_m.find(key);
97  if (hit_index_track_iter == hit_index_m.end()) {
98  int new_hit_index = this->__GetEdepArray__(real_track_id).size();
99  hit_index_m[key] = new_hit_index;
100  }
101  else {
102  hit_index = (*hit_index_track_iter).second;
103  }
104  auto const pid = geom->ChannelToWire(ch)[0].planeID();
105  auto const channel_id = pindex[pid];
106  double charge = ide.numElectrons;
107  if (hit_index < 0) {
108  // This particle energy deposition is never recorded so far. Create a new Edep
109  //float charge = ide.numElectrons * detp->ElectronsToADC();
110  this->__GetEdepArray__(real_track_id)
111  .emplace_back(pos, pid, pindex.size(), ide.energy, charge, channel_id);
112  }
113  else {
114  // Append charge to the relevant edep (@ hit_index)
115  //float charge = ide.numElectrons * detp->ElectronsToADC();
116  MCEdep& edep = this->__GetEdepArray__(real_track_id).at(hit_index);
117  edep.deps[channel_id].charge += charge;
118  edep.deps[channel_id].energy += ide.energy;
119  }
120  } // end looping over ides in this tick
121  } // end looping over ticks in this channel
122  } // end looping over channels
123 
124  if (_debug_mode) {
125  std::cout << Form(" Collected %zu particles' energy depositions...", _mc_edeps.size())
126  << std::endl;
127  }
128  }
std::vector< WireID > ChannelToWire(raw::ChannelID_t const channel) const
Returns a list of wires connected to the specified TPC channel.
std::map< geo::PlaneID, size_t > createPlaneIndexMap()
Definition: MCRecoEdep.cxx:26
std::map< unsigned int, size_t > _track_index
Definition: MCRecoEdep.h:138
Double_t edep
Definition: macro.C:13
std::vector< std::vector< sim::MCEdep > > _mc_edeps
Definition: MCRecoEdep.h:139
std::vector< sim::MCEdep > & __GetEdepArray__(unsigned int track_id)
Definition: MCRecoEdep.cxx:54
void sim::MCRecoEdep::MakeMCEdep ( const std::vector< sim::SimEnergyDeposit > &  sedArray)

Definition at line 130 of file MCRecoEdep.cxx.

References geo::GeometryCore::ChannelToWire(), sim::details::createPlaneIndexMap(), sim::MCEdep::deps, e, edep, geo::GeometryCore::Iterate(), geo::GeometryCore::NearestChannel(), geo::GeometryCore::PositionToCryostatID(), and geo::GeometryCore::PositionToTPCID().

131  {
132  _mc_edeps.clear();
133  _track_index.clear();
134 
136 
137  // Key map to identify a unique particle energy deposition point
138  std::map<std::pair<UniquePosition, unsigned int>, int> hit_index_m;
139 
140  auto pindex = details::createPlaneIndexMap();
141 
142  if (_debug_mode)
143  std::cout << "Processing " << sedArray.size() << " energy deposits..." << std::endl;
144  // Loop over channels
145  for (size_t i = 0; i < sedArray.size(); ++i) {
146 
147  // Get data to loop over
148  auto const& sed = sedArray[i];
149 
150  // David Caratelli:
151  // much of the code below is taken from the module:
152  // https://cdcvs.fnal.gov/redmine/projects/larsim/repository/revisions/develop/entry/larsim/ElectronDrift/SimDriftElectrons_module.cc
153 
154  // given this SimEnergyDeposit, find the TPC channel information
155  // "xyz" is the position of the energy deposit in world
156  // coordinates. Note that the units of distance in
157  // sim::SimEnergyDeposit are supposed to be cm.
158  auto const mp = sed.MidPoint();
159  // From the position in world coordinates, determine the
160  // cryostat and tpc. If somehow the step is outside a tpc
161  // (e.g., cosmic rays in rock) just move on to the next one.
162  unsigned int cryostat = 0;
163  try {
164  geom->PositionToCryostatID(mp);
165  }
166  catch (cet::exception& e) {
167  mf::LogWarning("SimDriftElectrons") << "step " // << energyDeposit << "\n"
168  << "cannot be found in a cryostat\n"
169  << e;
170  continue;
171  }
172  unsigned int tpc = 0;
173  try {
174  geom->PositionToTPCID(mp);
175  }
176  catch (cet::exception& e) {
177  mf::LogWarning("SimDriftElectrons") << "step " // << energyDeposit << "\n"
178  << "cannot be found in a TPC\n"
179  << e;
180  continue;
181  }
182  geo::TPCID const tpcid{cryostat, tpc};
183 
184  //Define charge drift direction: driftcoordinate (x, y or z) and driftsign (positive or negative). Also define coordinates perpendicular to drift direction.
185  // unused int driftcoordinate = std::abs(tpcGeo.DetectDriftDirection())-1; //x:0, y:1, z:2
186 
187  // make a collection of electrons for each plane
188  for (auto const& planeid : geom->Iterate<geo::PlaneID>()) {
189 
190  // grab the nearest channel to the fDriftClusterPos position
191  // David Caratelli, comment begin:
192  // NOTE: the below code works only when the drift coordinate is indeed in x (i.e. 0th coordinate)
193  // see code linked above for a much more careful treatment of the coordinate system
194  // David Caratelli, comment end.
195  raw::ChannelID_t ch;
196  try {
197  ch = geom->NearestChannel(mp, planeid);
198  }
199  catch (cet::exception& e) {
200  mf::LogWarning("SimDriftElectrons") << "step " // << energyDeposit << "\n"
201  << "nearest wire not in TPC\n"
202  << e;
203  continue;
204  }
205 
206  int track_id = sed.TrackID();
207 
208  if (track_id < 0) track_id = track_id * (-1);
209  unsigned int real_track_id = track_id;
210 
211  UniquePosition pos(mp.X(), mp.Y(), mp.Z());
212 
213  int hit_index = -1;
214  auto key = std::make_pair(pos, real_track_id);
215  auto hit_index_track_iter = hit_index_m.find(key);
216  if (hit_index_track_iter == hit_index_m.end()) {
217  int new_hit_index = this->__GetEdepArray__(real_track_id).size();
218  hit_index_m[key] = new_hit_index;
219  }
220  else {
221  hit_index = (*hit_index_track_iter).second;
222  }
223  auto const pid = geom->ChannelToWire(ch)[0].planeID();
224  auto const channel_id = pindex[pid];
225  double charge = sed.NumElectrons();
226  if (hit_index < 0) {
227  // This particle energy deposition is never recorded so far. Create a new Edep
228  //float charge = ide.numElectrons * detp->ElectronsToADC();
229  this->__GetEdepArray__(real_track_id)
230  .emplace_back(pos, pid, pindex.size(), sed.Energy(), charge, channel_id);
231  }
232  else {
233  // Append charge to the relevant edep (@ hit_index)
234  //float charge = ide.numElectrons * detp->ElectronsToADC();
235  MCEdep& edep = this->__GetEdepArray__(real_track_id).at(hit_index);
236  edep.deps[channel_id].charge += charge;
237  edep.deps[channel_id].energy += sed.Energy();
238  }
239  } // end looping over planes
240  } // end looping over SimEnergyDeposits
241 
242  if (_debug_mode) {
243  std::cout << Form(" Collected %zu particles' energy depositions...", _mc_edeps.size())
244  << std::endl;
245  }
246  }
details::range_type< T > Iterate() const
Initializes the specified ID with the ID of the first cryostat.
Definition: GeometryCore.h:541
std::vector< WireID > ChannelToWire(raw::ChannelID_t const channel) const
Returns a list of wires connected to the specified TPC channel.
The data type to uniquely identify a Plane.
Definition: geo_types.h:463
std::map< geo::PlaneID, size_t > createPlaneIndexMap()
Definition: MCRecoEdep.cxx:26
std::map< unsigned int, size_t > _track_index
Definition: MCRecoEdep.h:138
The data type to uniquely identify a TPC.
Definition: geo_types.h:381
Double_t edep
Definition: macro.C:13
std::vector< std::vector< sim::MCEdep > > _mc_edeps
Definition: MCRecoEdep.h:139
CryostatID PositionToCryostatID(Point_t const &point) const
Returns the ID of the cryostat at specified location.
std::vector< sim::MCEdep > & __GetEdepArray__(unsigned int track_id)
Definition: MCRecoEdep.cxx:54
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
raw::ChannelID_t NearestChannel(Point_t const &worldLoc, PlaneID const &planeid) const
Returns the ID of the channel nearest to the specified position.
unsigned int ChannelID_t
Type representing the ID of a readout channel.
Definition: RawTypes.h:28
Float_t e
Definition: plot.C:35
TPCID PositionToTPCID(Point_t const &point) const
Returns the ID of the TPC at specified location.
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
void sim::MCRecoEdep::MakeMCEdep ( const std::vector< sim::SimEnergyDepositLite > &  sedArray)

Definition at line 248 of file MCRecoEdep.cxx.

249  {
250  // Create a substitute array of sim::SimEnergyDeposit to avoid duplicating code...
251  // Note that this makes use of the explicit conversion operator
252  // defined in SimEnergyDepositLite class. Information will be partial.
253  // Most notably for MakeMCEdep, charge (numElectrons) will be 0.
254  std::vector<sim::SimEnergyDeposit> new_sedArray(sedArray.begin(), sedArray.end());
255  MakeMCEdep(new_sedArray);
256  }
void MakeMCEdep(const std::vector< sim::SimChannel > &schArray)
Definition: MCRecoEdep.cxx:62
const std::map<unsigned int, size_t> sim::MCRecoEdep::TrackIndexMap ( ) const
inline

Returns a map of track id <-> MCEdep vector index.

Definition at line 123 of file MCRecoEdep.h.

123 { return _track_index; }
std::map< unsigned int, size_t > _track_index
Definition: MCRecoEdep.h:138
int sim::MCRecoEdep::TrackToEdepIndex ( unsigned int  track_id) const
inline

Converts a track ID to MCEdep array index. Returns -1 if no corresponding array found .

Definition at line 113 of file MCRecoEdep.h.

Referenced by sim::MCTrackRecoAlg::Reconstruct(), and sim::MCShowerRecoAlg::Reconstruct().

114  {
115  auto iter = _track_index.find(track_id);
116  return (iter == _track_index.end() ? -1 : (int)((*iter).second));
117  }
std::map< unsigned int, size_t > _track_index
Definition: MCRecoEdep.h:138

Member Data Documentation

bool sim::MCRecoEdep::_debug_mode
protected

Definition at line 136 of file MCRecoEdep.h.

std::vector<std::vector<sim::MCEdep> > sim::MCRecoEdep::_mc_edeps
protected

Definition at line 139 of file MCRecoEdep.h.

bool sim::MCRecoEdep::_save_mchit
protected

Definition at line 137 of file MCRecoEdep.h.

std::map<unsigned int, size_t> sim::MCRecoEdep::_track_index
protected

Definition at line 138 of file MCRecoEdep.h.


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