LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
detsim::SimDriftElectrons Class Reference
Inheritance diagram for detsim::SimDriftElectrons:
art::EDProducer art::detail::Producer art::detail::LegacyModule art::Modifier art::ModuleBase art::ProductRegistryHelper

Classes

struct  ChannelBookKeeping
 

Public Types

using ModuleType = EDProducer
 
template<typename UserConfig , typename KeysToIgnore = void>
using Table = Modifier::Table< UserConfig, KeysToIgnore >
 

Public Member Functions

 SimDriftElectrons (fhicl::ParameterSet const &pset)
 
void produce (art::Event &evt) override
 
void beginJob () override
 
void doBeginJob (SharedResources const &resources)
 
void doEndJob ()
 
void doRespondToOpenInputFile (FileBlock const &fb)
 
void doRespondToCloseInputFile (FileBlock const &fb)
 
void doRespondToOpenOutputFiles (FileBlock const &fb)
 
void doRespondToCloseOutputFiles (FileBlock const &fb)
 
bool doBeginRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doEndRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)
 
void fillProductDescriptions ()
 
void registerProducts (ProductDescriptions &productsToRegister)
 
ModuleDescription const & moduleDescription () const
 
void setModuleDescription (ModuleDescription const &)
 
std::array< std::vector< ProductInfo >, NumBranchTypes > const & getConsumables () const
 
void sortConsumables (std::string const &current_process_name)
 
std::unique_ptr< Worker > makeWorker (WorkerParams const &wp)
 
template<typename T , BranchType BT>
ViewToken< T > consumesView (InputTag const &tag)
 
template<typename T , BranchType BT>
ViewToken< T > mayConsumeView (InputTag const &tag)
 

Protected Member Functions

ConsumesCollector & consumesCollector ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > consumes (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > consumesView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void consumesMany ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > mayConsume (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > mayConsumeView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void mayConsumeMany ()
 

Private Types

typedef std::unordered_map< raw::ChannelID_t, ChannelBookKeepingChannelMap_t
 

Private Attributes

art::InputTag fSimModuleLabel
 
CLHEP::RandGauss fRandGauss
 
double fElectronLifetime
 
double fElectronClusterSize
 
int fMinNumberOfElCluster
 
double fLongitudinalDiffusion
 
double fTransverseDiffusion
 
double fLifetimeCorr_const
 
double fLDiff_const
 
double fTDiff_const
 
double fRecipDriftVel [3]
 
bool fStoreDriftedElectronClusters
 
std::vector< ChannelMap_tfChannelMaps
 
size_t fNCryostats
 
std::vector< size_t > fNTPCs
 
std::vector< double > fLongDiff
 
std::vector< double > fTransDiff1
 
std::vector< double > fTransDiff2
 
std::vector< double > fnElDiff
 
std::vector< double > fnEnDiff
 
double fDriftClusterPos [3]
 
art::ServiceHandle< geo::Geometry const > fGeometry
 Handle to the Geometry service. More...
 

Detailed Description

Definition at line 95 of file SimDriftElectrons_module.cc.

Member Typedef Documentation

Definition at line 135 of file SimDriftElectrons_module.cc.

Definition at line 17 of file EDProducer.h.

template<typename UserConfig , typename KeysToIgnore = void>
using art::detail::Producer::Table = Modifier::Table<UserConfig, KeysToIgnore>
inherited

Definition at line 26 of file Producer.h.

Constructor & Destructor Documentation

detsim::SimDriftElectrons::SimDriftElectrons ( fhicl::ParameterSet const &  pset)
explicit

Definition at line 161 of file SimDriftElectrons_module.cc.

References art::detail::EngineCreator::createEngine(), fRandGauss, fSimModuleLabel, and fStoreDriftedElectronClusters.

162  : art::EDProducer{pset}
163  , fSimModuleLabel{pset.get<art::InputTag>("SimulationLabel")}
164  // create a default random engine; obtain the random seed from
165  // NuRandomService, unless overridden in configuration with key
166  // "Seed"
168  -> registerAndSeedEngine(createEngine(0), pset, "Seed")}
169  , fStoreDriftedElectronClusters{pset.get<bool>("StoreDriftedElectronClusters", false)}
170  {
171  produces<std::vector<sim::SimChannel>>();
172  if (fStoreDriftedElectronClusters) { produces<std::vector<sim::SimDriftedElectronCluster>>(); }
173  }
base_engine_t & createEngine(seed_t seed)

Member Function Documentation

void detsim::SimDriftElectrons::beginJob ( )
overridevirtual

Reimplemented from art::EDProducer.

Definition at line 176 of file SimDriftElectrons_module.cc.

References e, sim::LArG4Parameters::ElectronClusterSize(), fElectronClusterSize, fElectronLifetime, fGeometry, fLDiff_const, fLifetimeCorr_const, fLongitudinalDiffusion, fMinNumberOfElCluster, fNCryostats, fNTPCs, fRecipDriftVel, fTDiff_const, fTransverseDiffusion, sim::LArG4Parameters::LongitudinalDiffusion(), MF_LOG_DEBUG, sim::LArG4Parameters::MinNumberOfElCluster(), n, geo::GeometryCore::Ncryostats(), geo::GeometryCore::NTPC(), and sim::LArG4Parameters::TransverseDiffusion().

177  {
178  // Define the physical constants we'll use.
179 
180  auto const detProp =
182  fElectronLifetime = detProp.ElectronLifetime(); // Electron lifetime as returned by the
183  // DetectorProperties service assumed to be in us;
184  for (int i = 0; i < 3; ++i) {
185  double driftVelocity = detProp.DriftVelocity(detProp.Efield(i),
186  detProp.Temperature()) *
187  1.e-3; // Drift velocity as returned by the DetectorProperties service
188  // assumed to be in cm/us. Multiply by 1.e-3 to convert into
189  // LArSoft standard velocity units, cm/ns;
190 
191  fRecipDriftVel[i] = 1. / driftVelocity;
192  }
193 
194  // To-do: Move the parameters we fetch from "LArG4" to detector
195  // properties.
197  fElectronClusterSize = paramHandle->ElectronClusterSize();
199  fLongitudinalDiffusion = paramHandle->LongitudinalDiffusion(); // cm^2/ns units
200  fTransverseDiffusion = paramHandle->TransverseDiffusion(); // cm^2/ns units
201 
202  MF_LOG_DEBUG("SimDriftElectrons")
203  << " e lifetime (ns): " << fElectronLifetime
204  << "\n Temperature (K): " << detProp.Temperature()
205  << "\n Drift velocity (cm/ns): " << 1. / fRecipDriftVel[0] << " " << 1. / fRecipDriftVel[1]
206  << " " << 1. / fRecipDriftVel[2];
207 
208  // Opposite of lifetime. Convert from us to standard LArSoft time units, ns;
210  fLDiff_const = std::sqrt(2. * fLongitudinalDiffusion);
211  fTDiff_const = std::sqrt(2. * fTransverseDiffusion);
212 
213  // For this detector's geometry, save the number of cryostats and
214  // the number of TPCs within each cryostat.
216  fNTPCs.resize(fNCryostats);
217  for (size_t n = 0; n < fNCryostats; ++n)
219  }
unsigned int NTPC(CryostatID const &cryoid=cryostat_zero) const
Returns the total number of TPCs in the specified cryostat.
Definition: GeometryCore.h:686
unsigned int Ncryostats() const
Returns the number of cryostats in the detector.
Definition: GeometryCore.h:430
art::ServiceHandle< geo::Geometry const > fGeometry
Handle to the Geometry service.
double TransverseDiffusion() const
double ElectronClusterSize() const
int MinNumberOfElCluster() const
#define MF_LOG_DEBUG(id)
Char_t n[5]
double LongitudinalDiffusion() const
Float_t e
Definition: plot.C:35
The data type to uniquely identify a cryostat.
Definition: geo_types.h:192
template<typename T , BranchType BT>
ProductToken< T > art::ModuleBase::consumes ( InputTag const &  tag)
protectedinherited

Definition at line 61 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::consumes().

62  {
63  return collector_.consumes<T, BT>(tag);
64  }
ConsumesCollector collector_
Definition: ModuleBase.h:56
ProductToken< T > consumes(InputTag const &)
ConsumesCollector & art::ModuleBase::consumesCollector ( )
protectedinherited

Definition at line 57 of file ModuleBase.cc.

References art::ModuleBase::collector_.

58  {
59  return collector_;
60  }
ConsumesCollector collector_
Definition: ModuleBase.h:56
template<typename T , BranchType BT>
void art::ModuleBase::consumesMany ( )
protectedinherited

Definition at line 75 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::consumesMany().

76  {
77  collector_.consumesMany<T, BT>();
78  }
ConsumesCollector collector_
Definition: ModuleBase.h:56
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::ModuleBase::consumesView ( InputTag const &  )
protectedinherited
template<typename T , BranchType BT>
ViewToken<T> art::ModuleBase::consumesView ( InputTag const &  tag)
inherited

Definition at line 68 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::consumesView().

69  {
70  return collector_.consumesView<T, BT>(tag);
71  }
ConsumesCollector collector_
Definition: ModuleBase.h:56
ViewToken< Element > consumesView(InputTag const &)
void art::detail::Producer::doBeginJob ( SharedResources const &  resources)
inherited

Definition at line 22 of file Producer.cc.

References art::detail::Producer::beginJobWithFrame(), and art::detail::Producer::setupQueues().

23  {
24  setupQueues(resources);
25  ProcessingFrame const frame{ScheduleID{}};
26  beginJobWithFrame(frame);
27  }
virtual void setupQueues(SharedResources const &)=0
virtual void beginJobWithFrame(ProcessingFrame const &)=0
bool art::detail::Producer::doBeginRun ( RunPrincipal rp,
ModuleContext const &  mc 
)
inherited

Definition at line 65 of file Producer.cc.

References art::detail::Producer::beginRunWithFrame(), art::RangeSet::forRun(), art::RunPrincipal::makeRun(), r, art::RunPrincipal::runID(), and art::ModuleContext::scheduleID().

66  {
67  auto r = rp.makeRun(mc, RangeSet::forRun(rp.runID()));
68  ProcessingFrame const frame{mc.scheduleID()};
69  beginRunWithFrame(r, frame);
70  r.commitProducts();
71  return true;
72  }
TRandom r
Definition: spectrum.C:23
virtual void beginRunWithFrame(Run &, ProcessingFrame const &)=0
static RangeSet forRun(RunID)
Definition: RangeSet.cc:51
bool art::detail::Producer::doBeginSubRun ( SubRunPrincipal srp,
ModuleContext const &  mc 
)
inherited

Definition at line 85 of file Producer.cc.

References art::detail::Producer::beginSubRunWithFrame(), art::RangeSet::forSubRun(), art::SubRunPrincipal::makeSubRun(), art::ModuleContext::scheduleID(), and art::SubRunPrincipal::subRunID().

86  {
87  auto sr = srp.makeSubRun(mc, RangeSet::forSubRun(srp.subRunID()));
88  ProcessingFrame const frame{mc.scheduleID()};
89  beginSubRunWithFrame(sr, frame);
90  sr.commitProducts();
91  return true;
92  }
virtual void beginSubRunWithFrame(SubRun &, ProcessingFrame const &)=0
static RangeSet forSubRun(SubRunID)
Definition: RangeSet.cc:57
void art::detail::Producer::doEndJob ( )
inherited

Definition at line 30 of file Producer.cc.

References art::detail::Producer::endJobWithFrame().

31  {
32  ProcessingFrame const frame{ScheduleID{}};
33  endJobWithFrame(frame);
34  }
virtual void endJobWithFrame(ProcessingFrame const &)=0
bool art::detail::Producer::doEndRun ( RunPrincipal rp,
ModuleContext const &  mc 
)
inherited

Definition at line 75 of file Producer.cc.

References art::detail::Producer::endRunWithFrame(), art::RunPrincipal::makeRun(), r, art::ModuleContext::scheduleID(), and art::Principal::seenRanges().

76  {
77  auto r = rp.makeRun(mc, rp.seenRanges());
78  ProcessingFrame const frame{mc.scheduleID()};
79  endRunWithFrame(r, frame);
80  r.commitProducts();
81  return true;
82  }
TRandom r
Definition: spectrum.C:23
virtual void endRunWithFrame(Run &, ProcessingFrame const &)=0
bool art::detail::Producer::doEndSubRun ( SubRunPrincipal srp,
ModuleContext const &  mc 
)
inherited

Definition at line 95 of file Producer.cc.

References art::detail::Producer::endSubRunWithFrame(), art::SubRunPrincipal::makeSubRun(), art::ModuleContext::scheduleID(), and art::Principal::seenRanges().

96  {
97  auto sr = srp.makeSubRun(mc, srp.seenRanges());
98  ProcessingFrame const frame{mc.scheduleID()};
99  endSubRunWithFrame(sr, frame);
100  sr.commitProducts();
101  return true;
102  }
virtual void endSubRunWithFrame(SubRun &, ProcessingFrame const &)=0
bool art::detail::Producer::doEvent ( EventPrincipal ep,
ModuleContext const &  mc,
std::atomic< std::size_t > &  counts_run,
std::atomic< std::size_t > &  counts_passed,
std::atomic< std::size_t > &  counts_failed 
)
inherited

Definition at line 105 of file Producer.cc.

References art::detail::Producer::checkPutProducts_, e, art::EventPrincipal::makeEvent(), art::detail::Producer::produceWithFrame(), and art::ModuleContext::scheduleID().

110  {
111  auto e = ep.makeEvent(mc);
112  ++counts_run;
113  ProcessingFrame const frame{mc.scheduleID()};
114  produceWithFrame(e, frame);
115  e.commitProducts(checkPutProducts_, &expectedProducts<InEvent>());
116  ++counts_passed;
117  return true;
118  }
bool const checkPutProducts_
Definition: Producer.h:70
Float_t e
Definition: plot.C:35
virtual void produceWithFrame(Event &, ProcessingFrame const &)=0
void art::detail::Producer::doRespondToCloseInputFile ( FileBlock const &  fb)
inherited

Definition at line 44 of file Producer.cc.

References art::detail::Producer::respondToCloseInputFileWithFrame().

45  {
46  ProcessingFrame const frame{ScheduleID{}};
48  }
virtual void respondToCloseInputFileWithFrame(FileBlock const &, ProcessingFrame const &)=0
TFile fb("Li6.root")
void art::detail::Producer::doRespondToCloseOutputFiles ( FileBlock const &  fb)
inherited

Definition at line 58 of file Producer.cc.

References art::detail::Producer::respondToCloseOutputFilesWithFrame().

59  {
60  ProcessingFrame const frame{ScheduleID{}};
62  }
virtual void respondToCloseOutputFilesWithFrame(FileBlock const &, ProcessingFrame const &)=0
TFile fb("Li6.root")
void art::detail::Producer::doRespondToOpenInputFile ( FileBlock const &  fb)
inherited

Definition at line 37 of file Producer.cc.

References art::detail::Producer::respondToOpenInputFileWithFrame().

38  {
39  ProcessingFrame const frame{ScheduleID{}};
41  }
virtual void respondToOpenInputFileWithFrame(FileBlock const &, ProcessingFrame const &)=0
TFile fb("Li6.root")
void art::detail::Producer::doRespondToOpenOutputFiles ( FileBlock const &  fb)
inherited

Definition at line 51 of file Producer.cc.

References art::detail::Producer::respondToOpenOutputFilesWithFrame().

52  {
53  ProcessingFrame const frame{ScheduleID{}};
55  }
virtual void respondToOpenOutputFilesWithFrame(FileBlock const &, ProcessingFrame const &)=0
TFile fb("Li6.root")
void art::Modifier::fillProductDescriptions ( )
inherited

Definition at line 10 of file Modifier.cc.

References art::ProductRegistryHelper::fillDescriptions(), and art::ModuleBase::moduleDescription().

11  {
13  }
void fillDescriptions(ModuleDescription const &md)
ModuleDescription const & moduleDescription() const
Definition: ModuleBase.cc:13
std::array< std::vector< ProductInfo >, NumBranchTypes > const & art::ModuleBase::getConsumables ( ) const
inherited

Definition at line 43 of file ModuleBase.cc.

References art::ModuleBase::collector_, and art::ConsumesCollector::getConsumables().

44  {
45  return collector_.getConsumables();
46  }
ConsumesCollector collector_
Definition: ModuleBase.h:56
std::array< std::vector< ProductInfo >, NumBranchTypes > const & getConsumables() const
std::unique_ptr< Worker > art::ModuleBase::makeWorker ( WorkerParams const &  wp)
inherited

Definition at line 37 of file ModuleBase.cc.

References art::ModuleBase::doMakeWorker(), and art::NumBranchTypes.

38  {
39  return doMakeWorker(wp);
40  }
virtual std::unique_ptr< Worker > doMakeWorker(WorkerParams const &wp)=0
template<typename T , BranchType BT>
ProductToken< T > art::ModuleBase::mayConsume ( InputTag const &  tag)
protectedinherited

Definition at line 82 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::mayConsume().

83  {
84  return collector_.mayConsume<T, BT>(tag);
85  }
ProductToken< T > mayConsume(InputTag const &)
ConsumesCollector collector_
Definition: ModuleBase.h:56
template<typename T , BranchType BT>
void art::ModuleBase::mayConsumeMany ( )
protectedinherited

Definition at line 96 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::mayConsumeMany().

97  {
98  collector_.mayConsumeMany<T, BT>();
99  }
ConsumesCollector collector_
Definition: ModuleBase.h:56
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::ModuleBase::mayConsumeView ( InputTag const &  )
protectedinherited
template<typename T , BranchType BT>
ViewToken<T> art::ModuleBase::mayConsumeView ( InputTag const &  tag)
inherited

Definition at line 89 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::mayConsumeView().

90  {
91  return collector_.mayConsumeView<T, BT>(tag);
92  }
ConsumesCollector collector_
Definition: ModuleBase.h:56
ViewToken< Element > mayConsumeView(InputTag const &)
ModuleDescription const & art::ModuleBase::moduleDescription ( ) const
inherited

Definition at line 13 of file ModuleBase.cc.

References art::errors::LogicError.

Referenced by art::OutputModule::doRespondToOpenInputFile(), art::OutputModule::doWriteEvent(), art::Modifier::fillProductDescriptions(), art::OutputModule::makePlugins_(), art::OutputWorker::OutputWorker(), reco::shower::LArPandoraModularShowerCreation::produce(), art::Modifier::registerProducts(), and art::OutputModule::registerProducts().

14  {
15  if (md_.has_value()) {
16  return *md_;
17  }
18 
20  "There was an error while calling moduleDescription().\n"}
21  << "The moduleDescription() base-class member function cannot be called\n"
22  "during module construction. To determine which module is "
23  "responsible\n"
24  "for calling it, find the '<module type>:<module "
25  "label>@Construction'\n"
26  "tag in the message prefix above. Please contact artists@fnal.gov\n"
27  "for guidance.\n";
28  }
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
std::optional< ModuleDescription > md_
Definition: ModuleBase.h:55
void detsim::SimDriftElectrons::produce ( art::Event evt)
overridevirtual
Todo:
think about effects of drift between planes.
Todo:
think about effects of drift between planes
Todo:
check on what happens if we allow the tdc value to be
Todo:
beyond the end of the expected number of ticks

Implements art::EDProducer.

Definition at line 222 of file SimDriftElectrons_module.cc.

References util::abs(), sim::SimChannel::AddIonizationElectrons(), detsim::SimDriftElectrons::ChannelBookKeeping::channelIndex, geo::vect::coord(), geo::CryostatID::Cryostat, DEFINE_ART_MODULE, geo::TPCGeo::DetectDriftDirection(), e, energy, fChannelMaps, fDriftClusterPos, fElectronClusterSize, fGeometry, fLDiff_const, fLifetimeCorr_const, fLongDiff, fMinNumberOfElCluster, fNCryostats, fnElDiff, fnEnDiff, fRandGauss, fRecipDriftVel, fSimModuleLabel, fStoreDriftedElectronClusters, fTDiff_const, fTransDiff1, fTransDiff2, art::ProductRetriever::getByLabel(), geo::PlaneGeo::GetCenter(), geo::CryostatID::getInvalidID(), geo::TPCID::getInvalidID(), MF_LOG_DEBUG, geo::GeometryCore::NearestChannel(), geo::TPCGeo::Nplanes(), larsim::Utils::SCE::out_of_bounds(), geo::TPCGeo::Plane(), geo::TPCGeo::PlanePitch(), geo::GeometryCore::PositionToCryostatID(), geo::GeometryCore::PositionToTPCID(), detsim::SimDriftElectrons::ChannelBookKeeping::stepList, geo::vect::toPoint(), geo::TPCID::TPC, geo::GeometryCore::TPC(), and xx.

223  {
224  // Fetch the SimEnergyDeposit objects for this event.
225  typedef art::Handle<std::vector<sim::SimEnergyDeposit>> energyDepositHandle_t;
226  energyDepositHandle_t energyDepositHandle;
227  // If there aren't any energy deposits for this event, don't
228  // panic. It's possible someone is doing a study with events
229  // outside the TPC, or where there are only non-ionizing
230  // particles, or something like that.
231  if (!event.getByLabel(fSimModuleLabel, energyDepositHandle)) return;
232 
233  // Define the container for the SimChannel objects that will be
234  // transferred to the art::Event after the put statement below.
235  std::unique_ptr<std::vector<sim::SimChannel>> channels(new std::vector<sim::SimChannel>);
236  // Container for the SimDriftedElectronCluster objects
237  std::unique_ptr<std::vector<sim::SimDriftedElectronCluster>>
238  SimDriftedElectronClusterCollection(new std::vector<sim::SimDriftedElectronCluster>);
239 
240  // Clear the channel maps from the last event. Remember,
241  // fChannelMaps is an array[cryo][tpc] of maps.
242  size_t cryo = 0;
243  fChannelMaps.resize(fNCryostats);
244  for (auto& cryoData : fChannelMaps) { // each, a vector of maps
245  cryoData.clear();
246  }
247 
248  auto const clockData =
250  auto const& tpcClock = clockData.TPCClock();
251 
252  auto const detProp =
254  // We're going through the input vector by index, rather than by
255  // iterator, because we need the index number to compute the
256  // associations near the end of this method.
257  auto const& energyDeposits = *energyDepositHandle;
258  auto energyDepositsSize = energyDeposits.size();
259 
260  // For each energy deposit in this event
261  for (size_t edIndex = 0; edIndex < energyDepositsSize; ++edIndex) {
262  auto const& energyDeposit = energyDeposits[edIndex];
263 
264  // "xyz" is the position of the energy deposit in world
265  // coordinates. Note that the units of distance in
266  // sim::SimEnergyDeposit are supposed to be cm.
267  auto const mp = energyDeposit.MidPoint();
268  std::array<double, 3> xyz;
269  mp.GetCoordinates(data(xyz));
270 
271  // From the position in world coordinates, determine the
272  // cryostat and tpc. If somehow the step is outside a tpc
273  // (e.g., cosmic rays in rock) just move on to the next one.
274  unsigned int cryostat = 0;
275  cryostat = fGeometry->PositionToCryostatID(mp).Cryostat;
276  if (cryostat == geo::CryostatID::getInvalidID()) {
277  mf::LogWarning("SimDriftElectrons") << "step " // << energyDeposit << "\n"
278  << "cannot be found in a cryostat\n";
279  continue;
280  }
281 
282  geo::TPCID tpcid;
283  tpcid = fGeometry->PositionToTPCID(mp);
284  if (tpcid.TPC == geo::TPCID::getInvalidID()) {
285  mf::LogWarning("SimDriftElectrons") << "step " // << energyDeposit << "\n"
286  << "cannot be found in a TPC\n";
287  continue;
288  }
289 
290  const geo::TPCGeo& tpcGeo = fGeometry->TPC(tpcid);
291 
292  // The drift direction can be either in the positive
293  // or negative direction in any coordinate x, y or z.
294  // Charge drift in ...
295  // +x: tpcGeo.DetectDriftDirection()==1
296  // -x: tpcGeo.DetectDriftDirection()==-1
297  // +y: tpcGeo.DetectDriftDirection()==2
298  // -y tpcGeo.DetectDriftDirection()==-2
299  // +z: tpcGeo.DetectDriftDirection()==3
300  // -z: tpcGeo.DetectDriftDirection()==-3
301 
302  // Define charge drift direction: driftcoordinate (x, y or z) and
303  // driftsign (positive or negative). Also define coordinates perpendicular
304  // to drift direction.
305  int driftcoordinate = std::abs(tpcGeo.DetectDriftDirection()) - 1; // x:0, y:1, z:2
306 
307  int transversecoordinate1 = 0;
308  int transversecoordinate2 = 0;
309  if (driftcoordinate == 0) {
310  transversecoordinate1 = 1;
311  transversecoordinate2 = 2;
312  }
313  else if (driftcoordinate == 1) {
314  transversecoordinate1 = 0;
315  transversecoordinate2 = 2;
316  }
317  else if (driftcoordinate == 2) {
318  transversecoordinate1 = 0;
319  transversecoordinate2 = 1;
320  }
321 
322  if (transversecoordinate1 == transversecoordinate2)
323  continue; // this is the case when driftcoordinate != 0, 1 or 2
324 
325  int driftsign = 0; // 1: +x, +y or +z, -1: -x, -y or -z
326  if (tpcGeo.DetectDriftDirection() > 0)
327  driftsign = 1;
328  else
329  driftsign = -1;
330 
331  // Check for charge deposits behind charge readout planes
332  auto const& plane_location = tpcGeo.Plane(0).GetCenter();
333  auto const plane_location_coord = geo::vect::coord(plane_location, driftcoordinate);
334  if (driftsign == 1 && plane_location_coord < xyz[driftcoordinate]) continue;
335  if (driftsign == -1 && plane_location_coord > xyz[driftcoordinate]) continue;
336 
338  // Center of plane is also returned in cm units
339  double DriftDistance = std::abs(xyz[driftcoordinate] - plane_location_coord);
340 
341  // Space-charge effect (SCE): Get SCE {x,y,z} offsets for
342  // particular location in TPC
343  geo::Vector_t posOffsets{0.0, 0.0, 0.0};
344  double posOffsetxyz[3] = {0.0, 0.0, 0.0}; // need this array for the driftcoordinate and
345  // transversecoordinates
346  auto const* SCE = lar::providerFrom<spacecharge::SpaceChargeService>();
347  if (SCE->EnableSimSpatialSCE() == true) {
348  posOffsets = SCE->GetPosOffsets(mp);
349  if (larsim::Utils::SCE::out_of_bounds(posOffsets)) { continue; }
350  posOffsetxyz[0] = posOffsets.X();
351  posOffsetxyz[1] = posOffsets.Y();
352  posOffsetxyz[2] = posOffsets.Z();
353  }
354 
355  double avegagetransversePos1 = 0.;
356  double avegagetransversePos2 = 0.;
357 
358  DriftDistance += -1. * posOffsetxyz[driftcoordinate];
359  avegagetransversePos1 = xyz[transversecoordinate1] + posOffsetxyz[transversecoordinate1];
360  avegagetransversePos2 = xyz[transversecoordinate2] + posOffsetxyz[transversecoordinate2];
361 
362  // Space charge distortion could push the energy deposit beyond the wire
363  // plane (see issue #15131). Given that we don't have any subtlety in the
364  // simulation of this region, bringing the deposit exactly on the plane
365  // should be enough for the time being.
366  if (DriftDistance < 0.) DriftDistance = 0.;
367 
368  // Drift time in ns
369  double TDrift = DriftDistance * fRecipDriftVel[0];
370 
371  if (tpcGeo.Nplanes() == 2 &&
372  driftcoordinate == 0) { // special case for ArgoNeuT (Nplanes = 2 and drift direction =
373  // x): plane 0 is the second wire plane
374  TDrift = ((DriftDistance - tpcGeo.PlanePitch(0, 1)) * fRecipDriftVel[0] +
375  tpcGeo.PlanePitch(0, 1) * fRecipDriftVel[1]);
376  }
377 
378  const int nIonizedElectrons = energyDeposit.NumElectrons();
379  const double lifetimecorrection = TMath::Exp(TDrift / fLifetimeCorr_const);
380  const double energy = energyDeposit.Energy();
381 
382  // if we have no electrons (too small energy or too large recombination)
383  // we are done already here
384  if (nIonizedElectrons <= 0) {
385  MF_LOG_DEBUG("SimDriftElectrons")
386  << "step " // << energyDeposit << "\n"
387  << "No electrons drifted to readout, " << energy << " MeV lost.";
388  continue;
389  }
390 
391  // includes the effect of lifetime: lifetimecorrection = exp[-tdrift/tau]
392  const double nElectrons = nIonizedElectrons * lifetimecorrection;
393 
394  // Longitudinal & transverse diffusion sigma (cm)
395  double SqrtT = std::sqrt(TDrift);
396  double LDiffSig = SqrtT * fLDiff_const;
397  double TDiffSig = SqrtT * fTDiff_const;
398  double electronclsize = fElectronClusterSize;
399 
400  // Number of electron clusters.
401  int nClus = (int)std::ceil(nElectrons / electronclsize);
402  if (nClus < fMinNumberOfElCluster) {
403  electronclsize = nElectrons / fMinNumberOfElCluster;
404  if (electronclsize < 1.0) { electronclsize = 1.0; }
405  nClus = (int)std::ceil(nElectrons / electronclsize);
406  }
407 
408  // Empty and resize the electron-cluster vectors.
409  fLongDiff.clear();
410  fTransDiff1.clear();
411  fTransDiff2.clear();
412  fnElDiff.clear();
413  fnEnDiff.clear();
414  fLongDiff.resize(nClus);
415  fTransDiff1.resize(nClus);
416  fTransDiff2.resize(nClus);
417  fnElDiff.resize(nClus, electronclsize);
418  fnEnDiff.resize(nClus);
419 
420  // fix the number of electrons in the last cluster, that has a smaller size
421  fnElDiff.back() = nElectrons - (nClus - 1) * electronclsize;
422 
423  for (size_t xx = 0; xx < fnElDiff.size(); ++xx) {
424  if (nElectrons > 0)
425  fnEnDiff[xx] = energy / nElectrons * fnElDiff[xx];
426  else
427  fnEnDiff[xx] = 0.;
428  }
429 
430  // Smear drift times by longitudinal diffusion
431  if (LDiffSig > 0.0)
432  fRandGauss.fireArray(nClus, &fLongDiff[0], 0., LDiffSig);
433  else
434  fLongDiff.assign(nClus, 0.0);
435 
436  if (TDiffSig > 0.0) {
437  // Smear the coordinates in plane perpendicular to drift direction by the transverse diffusion
438  fRandGauss.fireArray(nClus, &fTransDiff1[0], avegagetransversePos1, TDiffSig);
439  fRandGauss.fireArray(nClus, &fTransDiff2[0], avegagetransversePos2, TDiffSig);
440  }
441  else {
442  fTransDiff1.assign(nClus, avegagetransversePos1);
443  fTransDiff2.assign(nClus, avegagetransversePos2);
444  }
445 
446  // make a collection of electrons for each plane
447  for (unsigned int p = 0; p < tpcGeo.Nplanes(); ++p) {
448  // FIXME (KJK): Shouldn't this be the location for each plane
449  // we're iterating through? Right now it uses whatever value
450  // of plane_location_coord was set above.
451  fDriftClusterPos[driftcoordinate] = plane_location_coord;
452 
453  // Drift nClus electron clusters to the induction plane
454  for (int k = 0; k < nClus; ++k) {
455 
456  // Correct drift time for longitudinal diffusion and plane
457  double TDiff = TDrift + fLongDiff[k] * fRecipDriftVel[0];
458 
459  // Take into account different Efields between planes
460  // Also take into account special case for ArgoNeuT (Nplanes = 2 and
461  // drift direction = x): plane 0 is the second wire plane
462  for (unsigned int ip = 0; ip < p; ++ip) {
463  TDiff +=
464  tpcGeo.PlanePitch(ip + 1, ip) *
465  fRecipDriftVel[(tpcGeo.Nplanes() == 2 && driftcoordinate == 0) ? ip + 2 : ip + 1];
466  }
467 
468  fDriftClusterPos[transversecoordinate1] = fTransDiff1[k];
469  fDriftClusterPos[transversecoordinate2] = fTransDiff2[k];
470 
472 
473  // grab the nearest channel to the fDriftClusterPos position
474  try {
476  geo::vect::toPoint(fDriftClusterPos), geo::PlaneID{tpcid, p});
477 
480  // Add potential decay/capture/etc delay effect, simTime.
481  auto const simTime = energyDeposit.Time();
482  unsigned int tdc = tpcClock.Ticks(clockData.G4ToElecTime(TDiff + simTime));
483 
484  // Find whether we already have this channel in our map.
485  ChannelMap_t& channelDataMap = fChannelMaps[cryostat];
486  auto search = channelDataMap.find(channel);
487 
488  // We will find (or create) the pointer to a
489  // sim::SimChannel.
490  size_t channelIndex = 0;
491 
492  // Have we created the sim::SimChannel corresponding to
493  // channel ID?
494  if (search == channelDataMap.end()) {
495  // We haven't. Initialize the bookkeeping information
496  // for this channel.
497  ChannelBookKeeping bookKeeping;
498 
499  // Add a new channel to the end of the list we'll
500  // write out after we've processed this event.
501  bookKeeping.channelIndex = channels->size();
502  channels->emplace_back(channel);
503  channelIndex = bookKeeping.channelIndex;
504 
505  // Initialize a vector with the index of the step that
506  // created this channel.
507  bookKeeping.stepList.push_back(edIndex);
508 
509  // Save the bookkeeping information for this channel.
510  channelDataMap[channel] = bookKeeping;
511  }
512  else {
513  // We've created this SimChannel for a previous energy
514  // deposit. Get its address.
515 
516  auto& bookKeeping = search->second;
517  channelIndex = bookKeeping.channelIndex;
518 
519  // Has this step contributed to this channel before?
520  auto& stepList = bookKeeping.stepList;
521  if (!std::binary_search(stepList.begin(), stepList.end(), edIndex)) {
522  // No, so add this step's index to the list.
523  stepList.push_back(edIndex);
524  }
525  }
526 
527  sim::SimChannel* channelPtr = &(channels->at(channelIndex));
528 
529  // Add the electron clusters and energy to the
530  // sim::SimChannel
531  channelPtr->AddIonizationElectrons(energyDeposit.TrackID(),
532  tdc,
533  fnElDiff[k],
534  data(xyz),
535  fnEnDiff[k],
536  energyDeposit.OrigTrackID());
538  SimDriftedElectronClusterCollection->emplace_back(
539  fnElDiff[k],
540  TDiff + simTime, // timing
541  geo::Point_t{mp.X(), mp.Y(), mp.Z()}, // mean position of the deposited energy
542  geo::Point_t{fDriftClusterPos[0],
543  fDriftClusterPos[1],
544  fDriftClusterPos[2]}, // final position of the drifted cluster
545  geo::Point_t{
546  LDiffSig, TDiffSig, TDiffSig}, // Longitudinal (X) and transverse (Y,Z) diffusion
547  fnEnDiff[k], // deposited energy that originated this cluster
548  energyDeposit.TrackID());
549  }
550  catch (cet::exception& e) {
551  mf::LogDebug("SimDriftElectrons")
552  << "unable to drift electrons from point (" << xyz[0] << "," << xyz[1] << ","
553  << xyz[2] << ") with exception " << e;
554  } // end try to determine channel
555  } // end loop over clusters
556  } // end loop over planes
557  } // for each sim::SimEnergyDeposit
558 
559  // Write the sim::SimChannel collection.
560  event.put(std::move(channels));
561  if (fStoreDriftedElectronClusters) event.put(std::move(SimDriftedElectronClusterCollection));
562  }
double PlanePitch(unsigned int p1=0, unsigned int p2=1) const
Returns the center of the TPC volume in world coordinates [cm].
Definition: TPCGeo.cxx:361
Double_t xx
Definition: macro.C:12
Energy deposited on a readout channel by simulated tracks.
Definition: SimChannel.h:136
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
Definition: geo_vectors.h:160
Point_t const & GetCenter() const
Returns the centre of the wire plane in world coordinates [cm].
Definition: PlaneGeo.h:435
bool out_of_bounds(geo::Vector_t const &offset)
auto coord(Vector &v, unsigned int n) noexcept
Returns an object to manage the coordinate n of a vector.
unsigned int Nplanes() const
Number of planes in this tpc.
Definition: TPCGeo.h:137
The data type to uniquely identify a Plane.
Definition: geo_types.h:463
Geometry information for a single TPC.
Definition: TPCGeo.h:36
::geo::Point_t toPoint(Point const &p)
Convert the specified point into a geo::Point_t.
constexpr auto abs(T v)
Returns the absolute value of the argument.
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:211
art::ServiceHandle< geo::Geometry const > fGeometry
Handle to the Geometry service.
TPCGeo const & TPC(TPCID const &tpcid=tpc_zero) const
Returns the specified TPC.
Definition: GeometryCore.h:722
static constexpr CryostatID_t getInvalidID()
Return the value of the invalid ID as a r-value.
Definition: geo_types.h:285
static constexpr TPCID_t getInvalidID()
Return the value of the invalid TPC ID as a r-value.
Definition: geo_types.h:458
double energy
Definition: plottest35.C:25
The data type to uniquely identify a TPC.
Definition: geo_types.h:381
std::vector< ChannelMap_t > fChannelMaps
CryostatID PositionToCryostatID(Point_t const &point) const
Returns the ID of the cryostat at specified location.
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
Definition: geo_vectors.h:180
short int DetectDriftDirection() const
Returns the expected drift direction based on geometry.
Definition: TPCGeo.cxx:149
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
#define MF_LOG_DEBUG(id)
raw::ChannelID_t NearestChannel(Point_t const &worldLoc, PlaneID const &planeid) const
Returns the ID of the channel nearest to the specified position.
PlaneGeo const & Plane(geo::View_t view) const
Return the plane in the tpc with View_t view.
Definition: TPCGeo.cxx:252
void AddIonizationElectrons(TrackID_t trackID, TDC_t tdc, double numberElectrons, double const *xyz, double energy, TrackID_t origTrackID=util::kBogusI)
Add ionization electrons and energy to this channel.
Definition: SimChannel.cxx:49
std::unordered_map< raw::ChannelID_t, ChannelBookKeeping > ChannelMap_t
unsigned int ChannelID_t
Type representing the ID of a readout channel.
Definition: RawTypes.h:28
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:399
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
Event finding and building.
void art::Modifier::registerProducts ( ProductDescriptions productsToRegister)
inherited

Definition at line 16 of file Modifier.cc.

References art::ModuleBase::moduleDescription(), and art::ProductRegistryHelper::registerProducts().

17  {
18  ProductRegistryHelper::registerProducts(productsToRegister,
20  }
void registerProducts(ProductDescriptions &productsToRegister, ModuleDescription const &md)
ModuleDescription const & moduleDescription() const
Definition: ModuleBase.cc:13
void art::ModuleBase::setModuleDescription ( ModuleDescription const &  md)
inherited

Definition at line 31 of file ModuleBase.cc.

References art::ModuleBase::md_.

32  {
33  md_ = md;
34  }
std::optional< ModuleDescription > md_
Definition: ModuleBase.h:55
void art::ModuleBase::sortConsumables ( std::string const &  current_process_name)
inherited

Definition at line 49 of file ModuleBase.cc.

References art::ModuleBase::collector_, and art::ConsumesCollector::sortConsumables().

50  {
51  // Now that we know we have seen all the consumes declarations,
52  // sort the results for fast lookup later.
53  collector_.sortConsumables(current_process_name);
54  }
ConsumesCollector collector_
Definition: ModuleBase.h:56
void sortConsumables(std::string const &current_process_name)

Member Data Documentation

std::vector<ChannelMap_t> detsim::SimDriftElectrons::fChannelMaps
private

Definition at line 138 of file SimDriftElectrons_module.cc.

Referenced by produce().

double detsim::SimDriftElectrons::fDriftClusterPos[3]
private

Definition at line 154 of file SimDriftElectrons_module.cc.

Referenced by produce().

double detsim::SimDriftElectrons::fElectronClusterSize
private

Definition at line 112 of file SimDriftElectrons_module.cc.

Referenced by beginJob(), and produce().

double detsim::SimDriftElectrons::fElectronLifetime
private

Definition at line 111 of file SimDriftElectrons_module.cc.

Referenced by beginJob().

art::ServiceHandle<geo::Geometry const> detsim::SimDriftElectrons::fGeometry
private

Handle to the Geometry service.

Definition at line 156 of file SimDriftElectrons_module.cc.

Referenced by beginJob(), and produce().

double detsim::SimDriftElectrons::fLDiff_const
private

Definition at line 118 of file SimDriftElectrons_module.cc.

Referenced by beginJob(), and produce().

double detsim::SimDriftElectrons::fLifetimeCorr_const
private

Definition at line 117 of file SimDriftElectrons_module.cc.

Referenced by beginJob(), and produce().

std::vector<double> detsim::SimDriftElectrons::fLongDiff
private

Definition at line 148 of file SimDriftElectrons_module.cc.

Referenced by produce().

double detsim::SimDriftElectrons::fLongitudinalDiffusion
private

Definition at line 114 of file SimDriftElectrons_module.cc.

Referenced by beginJob().

int detsim::SimDriftElectrons::fMinNumberOfElCluster
private

Definition at line 113 of file SimDriftElectrons_module.cc.

Referenced by beginJob(), and produce().

size_t detsim::SimDriftElectrons::fNCryostats
private

Definition at line 144 of file SimDriftElectrons_module.cc.

Referenced by beginJob(), and produce().

std::vector<double> detsim::SimDriftElectrons::fnElDiff
private

Definition at line 151 of file SimDriftElectrons_module.cc.

Referenced by produce().

std::vector<double> detsim::SimDriftElectrons::fnEnDiff
private

Definition at line 152 of file SimDriftElectrons_module.cc.

Referenced by produce().

std::vector<size_t> detsim::SimDriftElectrons::fNTPCs
private

Definition at line 145 of file SimDriftElectrons_module.cc.

Referenced by beginJob().

CLHEP::RandGauss detsim::SimDriftElectrons::fRandGauss
private

Definition at line 109 of file SimDriftElectrons_module.cc.

Referenced by produce(), and SimDriftElectrons().

double detsim::SimDriftElectrons::fRecipDriftVel[3]
private

Definition at line 120 of file SimDriftElectrons_module.cc.

Referenced by beginJob(), and produce().

art::InputTag detsim::SimDriftElectrons::fSimModuleLabel
private

Definition at line 107 of file SimDriftElectrons_module.cc.

Referenced by produce(), and SimDriftElectrons().

bool detsim::SimDriftElectrons::fStoreDriftedElectronClusters
private

Definition at line 122 of file SimDriftElectrons_module.cc.

Referenced by produce(), and SimDriftElectrons().

double detsim::SimDriftElectrons::fTDiff_const
private

Definition at line 119 of file SimDriftElectrons_module.cc.

Referenced by beginJob(), and produce().

std::vector<double> detsim::SimDriftElectrons::fTransDiff1
private

Definition at line 149 of file SimDriftElectrons_module.cc.

Referenced by produce().

std::vector<double> detsim::SimDriftElectrons::fTransDiff2
private

Definition at line 150 of file SimDriftElectrons_module.cc.

Referenced by produce().

double detsim::SimDriftElectrons::fTransverseDiffusion
private

Definition at line 115 of file SimDriftElectrons_module.cc.

Referenced by beginJob().


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