226 energyDepositHandle_t energyDepositHandle;
235 std::unique_ptr<std::vector<sim::SimChannel>> channels(
new std::vector<sim::SimChannel>);
237 std::unique_ptr<std::vector<sim::SimDriftedElectronCluster>>
238 SimDriftedElectronClusterCollection(
new std::vector<sim::SimDriftedElectronCluster>);
248 auto const clockData =
250 auto const& tpcClock = clockData.TPCClock();
257 auto const& energyDeposits = *energyDepositHandle;
258 auto energyDepositsSize = energyDeposits.size();
261 for (
size_t edIndex = 0; edIndex < energyDepositsSize; ++edIndex) {
262 auto const& energyDeposit = energyDeposits[edIndex];
267 auto const mp = energyDeposit.MidPoint();
268 std::array<double, 3> xyz;
269 mp.GetCoordinates(data(xyz));
274 unsigned int cryostat = 0;
278 <<
"cannot be found in a cryostat\n";
286 <<
"cannot be found in a TPC\n";
307 int transversecoordinate1 = 0;
308 int transversecoordinate2 = 0;
309 if (driftcoordinate == 0) {
310 transversecoordinate1 = 1;
311 transversecoordinate2 = 2;
313 else if (driftcoordinate == 1) {
314 transversecoordinate1 = 0;
315 transversecoordinate2 = 2;
317 else if (driftcoordinate == 2) {
318 transversecoordinate1 = 0;
319 transversecoordinate2 = 1;
322 if (transversecoordinate1 == transversecoordinate2)
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;
339 double DriftDistance =
std::abs(xyz[driftcoordinate] - plane_location_coord);
344 double posOffsetxyz[3] = {0.0, 0.0, 0.0};
346 auto const* SCE = lar::providerFrom<spacecharge::SpaceChargeService>();
347 if (SCE->EnableSimSpatialSCE() ==
true) {
348 posOffsets = SCE->GetPosOffsets(mp);
350 posOffsetxyz[0] = posOffsets.X();
351 posOffsetxyz[1] = posOffsets.Y();
352 posOffsetxyz[2] = posOffsets.Z();
355 double avegagetransversePos1 = 0.;
356 double avegagetransversePos2 = 0.;
358 DriftDistance += -1. * posOffsetxyz[driftcoordinate];
359 avegagetransversePos1 = xyz[transversecoordinate1] + posOffsetxyz[transversecoordinate1];
360 avegagetransversePos2 = xyz[transversecoordinate2] + posOffsetxyz[transversecoordinate2];
366 if (DriftDistance < 0.) DriftDistance = 0.;
372 driftcoordinate == 0) {
374 TDrift = ((DriftDistance - tpcGeo.
PlanePitch(0, 1)) * fRecipDriftVel[0] +
378 const int nIonizedElectrons = energyDeposit.NumElectrons();
380 const double energy = energyDeposit.Energy();
384 if (nIonizedElectrons <= 0) {
387 <<
"No electrons drifted to readout, " << energy <<
" MeV lost.";
392 const double nElectrons = nIonizedElectrons * lifetimecorrection;
395 double SqrtT = std::sqrt(TDrift);
401 int nClus = (int)std::ceil(nElectrons / electronclsize);
404 if (electronclsize < 1.0) { electronclsize = 1.0; }
405 nClus = (int)std::ceil(nElectrons / electronclsize);
417 fnElDiff.resize(nClus, electronclsize);
421 fnElDiff.back() = nElectrons - (nClus - 1) * electronclsize;
436 if (TDiffSig > 0.0) {
447 for (
unsigned int p = 0; p < tpcGeo.
Nplanes(); ++p) {
454 for (
int k = 0; k < nClus; ++k) {
457 double TDiff = TDrift +
fLongDiff[k] * fRecipDriftVel[0];
462 for (
unsigned int ip = 0; ip < p; ++ip) {
465 fRecipDriftVel[(tpcGeo.
Nplanes() == 2 && driftcoordinate == 0) ? ip + 2 : ip + 1];
469 fDriftClusterPos[transversecoordinate2] =
fTransDiff2[k];
481 auto const simTime = energyDeposit.Time();
482 unsigned int tdc = tpcClock.Ticks(clockData.G4ToElecTime(TDiff + simTime));
486 auto search = channelDataMap.find(channel);
490 size_t channelIndex = 0;
494 if (search == channelDataMap.end()) {
497 ChannelBookKeeping bookKeeping;
501 bookKeeping.channelIndex = channels->size();
502 channels->emplace_back(channel);
503 channelIndex = bookKeeping.channelIndex;
507 bookKeeping.stepList.push_back(edIndex);
510 channelDataMap[channel] = bookKeeping;
516 auto& bookKeeping = search->second;
517 channelIndex = bookKeeping.channelIndex;
520 auto& stepList = bookKeeping.stepList;
521 if (!std::binary_search(stepList.begin(), stepList.end(), edIndex)) {
523 stepList.push_back(edIndex);
536 energyDeposit.OrigTrackID());
538 SimDriftedElectronClusterCollection->emplace_back(
544 fDriftClusterPos[2]},
546 LDiffSig, TDiffSig, TDiffSig},
548 energyDeposit.TrackID());
552 <<
"unable to drift electrons from point (" << xyz[0] <<
"," << xyz[1] <<
"," 553 << xyz[2] <<
") with exception " <<
e;
560 event.put(std::move(channels));
std::vector< double > fTransDiff2
double PlanePitch(unsigned int p1=0, unsigned int p2=1) const
Returns the center of the TPC volume in world coordinates [cm].
double fElectronClusterSize
Energy deposited on a readout channel by simulated tracks.
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
Point_t const & GetCenter() const
Returns the centre of the wire plane in world coordinates [cm].
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.
The data type to uniquely identify a Plane.
Geometry information for a single TPC.
::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.
CLHEP::RandGauss fRandGauss
art::ServiceHandle< geo::Geometry const > fGeometry
Handle to the Geometry service.
TPCGeo const & TPC(TPCID const &tpcid=tpc_zero) const
Returns the specified TPC.
static constexpr CryostatID_t getInvalidID()
Return the value of the invalid ID as a r-value.
static constexpr TPCID_t getInvalidID()
Return the value of the invalid TPC ID as a r-value.
std::vector< double > fnEnDiff
The data type to uniquely identify a TPC.
std::vector< ChannelMap_t > fChannelMaps
art::InputTag fSimModuleLabel
double fLifetimeCorr_const
std::vector< double > fLongDiff
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.
double fDriftClusterPos[3]
std::vector< double > fTransDiff1
short int DetectDriftDirection() const
Returns the expected drift direction based on geometry.
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
std::vector< double > fnElDiff
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.
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.
std::unordered_map< raw::ChannelID_t, ChannelBookKeeping > ChannelMap_t
unsigned int ChannelID_t
Type representing the ID of a readout channel.
TPCID_t TPC
Index of the TPC within its cryostat.
bool fStoreDriftedElectronClusters
int fMinNumberOfElCluster
TPCID PositionToTPCID(Point_t const &point) const
Returns the ID of the TPC at specified location.
cet::coded_exception< error, detail::translate > exception
Event finding and building.