288 energyDepositHandle_t energyDepositHandle;
298 std::unique_ptr< std::vector<sim::SimChannel> > channels(
new std::vector<sim::SimChannel>);
300 std::unique_ptr< std::vector<sim::SimDriftedElectronCluster> > SimDriftedElectronClusterCollection(
new std::vector<sim::SimDriftedElectronCluster>);
307 cryoData.resize(
fNTPCs[cryo++]);
308 for (
auto& channelsMap: cryoData) channelsMap.clear();
314 auto const& energyDeposits = *energyDepositHandle;
315 auto energyDepositsSize = energyDeposits.size();
318 for (
size_t edIndex = 0; edIndex < energyDepositsSize; ++edIndex )
320 auto const& energyDeposit = energyDeposits[edIndex];
325 auto const mp = energyDeposit.MidPoint();
326 double const xyz[3] = { mp.X(), mp.Y(), mp.Z() };
331 unsigned int cryostat = 0;
337 <<
"cannot be found in a cryostat\n" 342 unsigned int tpc = 0;
348 <<
"cannot be found in a TPC\n" 369 int transversecoordinate1 = 0;
370 int transversecoordinate2 = 0;
371 if(driftcoordinate == 0)
373 transversecoordinate1 = 1;
374 transversecoordinate2 = 2;
376 else if(driftcoordinate == 1)
378 transversecoordinate1 = 0;
379 transversecoordinate2 = 2;
381 else if(driftcoordinate == 2)
383 transversecoordinate1 = 0;
384 transversecoordinate2 = 1;
387 if(transversecoordinate1 == transversecoordinate2)
continue;
394 if(driftsign == 1 && tpcGeo.
PlaneLocation(0)[driftcoordinate] < xyz[driftcoordinate] )
396 if(driftsign == -1 && tpcGeo.
PlaneLocation(0)[driftcoordinate] > xyz[driftcoordinate] )
403 double DriftDistance = std::abs(xyz[driftcoordinate] - tpcGeo.
PlaneLocation(0)[driftcoordinate]);
408 double posOffsetxyz[3] = {0.0,0.0,0.0};
409 auto const* SCE = lar::providerFrom<spacecharge::SpaceChargeService>();
410 if (SCE->EnableSimSpatialSCE() ==
true)
412 posOffsets = SCE->GetPosOffsets(mp);
413 posOffsetxyz[0] = posOffsets.X();
414 posOffsetxyz[1] = posOffsets.Y();
415 posOffsetxyz[2] = posOffsets.Z();
418 double avegagetransversePos1 = 0.;
419 double avegagetransversePos2 = 0.;
421 DriftDistance += -1.*posOffsetxyz[driftcoordinate];
422 avegagetransversePos1 = xyz[transversecoordinate1] + posOffsetxyz[transversecoordinate1];
423 avegagetransversePos2 = xyz[transversecoordinate2] + posOffsetxyz[transversecoordinate2];
430 if (DriftDistance < 0.) DriftDistance = 0.;
435 if (tpcGeo.
Nplanes() == 2 && driftcoordinate == 0){
436 TDrift = ((DriftDistance - tpcGeo.
PlanePitch(0,1)) * fRecipDriftVel[0]
437 + tpcGeo.
PlanePitch(0,1) * fRecipDriftVel[1]);
446 const double energy = energyDeposit.Energy();
450 if (nIonizedElectrons <= 0) {
453 <<
"No electrons drifted to readout, " << energy <<
" MeV lost.";
458 const double nElectrons = nIonizedElectrons * lifetimecorrection;
462 double SqrtT = std::sqrt(TDrift);
468 int nClus = (int) std::ceil(nElectrons / electronclsize);
472 if (electronclsize < 1.0)
474 electronclsize = 1.0;
476 nClus = (int) std::ceil(nElectrons / electronclsize);
488 fnElDiff.resize(nClus, electronclsize);
492 fnElDiff.back() = nElectrons - (nClus-1)*electronclsize;
507 if (TDiffSig > 0.0) {
520 for(
size_t p = 0; p < tpcGeo.
Nplanes(); ++p){
525 for(
int k = 0; k < nClus; ++k){
533 double TDiff = TDrift +
fLongDiff[k] * fRecipDriftVel[0];
537 for (
size_t ip = 0; ip<p; ++ip){
538 TDiff += (tpcGeo.
PlaneLocation(ip+1)[driftcoordinate] - tpcGeo.
PlaneLocation(ip)[driftcoordinate]) * fRecipDriftVel[(tpcGeo.
Nplanes() == 2 && driftcoordinate == 0)?ip+2:ip+1];
575 auto const simTime = energyDeposit.Time();
579 ChannelMap_t& channelDataMap = fChannelMaps[cryostat][tpc];
580 auto search = channelDataMap.find(channel);
585 size_t channelIndex=0;
589 if (search == channelDataMap.end())
595 ChannelBookKeeping_t bookKeeping;
599 bookKeeping.channelIndex = channels->size();
600 channels->emplace_back( channel );
601 channelIndex = bookKeeping.channelIndex;
610 bookKeeping.stepList.push_back( edIndex );
613 channelDataMap[channel] = bookKeeping;
621 auto& bookKeeping = search->second;
622 channelIndex = bookKeeping.channelIndex;
626 auto& stepList = bookKeeping.stepList;
627 auto stepSearch = std::find(stepList.begin(), stepList.end(), edIndex );
628 if ( stepSearch == stepList.end() ) {
630 stepList.push_back( edIndex );
657 {mp.X(),mp.Y(),mp.Z()},
659 {LDiffSig,TDiffSig,TDiffSig},
661 energyDeposit.TrackID()) );
667 mf::LogWarning(
"SimDriftElectrons") <<
"unable to drift electrons from point (" 668 << xyz[0] <<
"," << xyz[1] <<
"," << xyz[2]
669 <<
") with exception " <<
e;
711 event.put(std::move(channels));
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
double PlanePitch(unsigned int p1=0, unsigned int p2=1) const
Returns the center of the TPC volume in world coordinates [cm].
std::vector< double > fnEnDiff
double fElectronClusterSize
Energy deposited on a readout channel by simulated tracks.
CryostatGeo const & PositionToCryostat(geo::Point_t const &point) const
Returns the cryostat at specified location.
void CalculateIonizationAndScintillation(sim::SimEnergyDeposit const &edep)
unsigned int Nplanes() const
Number of planes in this tpc.
Geometry information for a single TPC.
std::vector< std::vector< ChannelMap_t > > fChannelMaps
art::ServiceHandle< geo::Geometry > fGeometry
Handle to the Geometry service.
int Ticks() const
Current clock tick (that is, the number of tick Time() falls in).
std::vector< size_t > fNTPCs
geo::TPCGeo const & PositionToTPC(geo::Point_t const &point) const
Returns the TPC at specified location.
std::vector< double > fTransDiff1
std::vector< double > fnElDiff
const detinfo::DetectorClocks * fTimeService
larg4::ISCalcSeparate fISAlg
virtual double G4ToElecTime(double g4_time) const =0
Given a simulation time [ns], converts it into electronics time [µs].
art::InputTag fSimModuleLabel
double fLifetimeCorr_const
void AddIonizationElectrons(TrackID_t trackID, TDC_t tdc, double numberElectrons, double const *xyz, double energy)
Add ionization electrons and energy to this channel.
std::map< raw::ChannelID_t, ChannelBookKeeping_t > ChannelMap_t
double fDriftClusterPos[3]
double NumberIonizationElectrons() const
short int DetectDriftDirection() const
Returns the expected drift direction based on geometry.
::detinfo::ElecClock fClock
TPC electronics clock.
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
std::unique_ptr< CLHEP::RandGauss > fRandGauss
TPCGeo const & TPC(unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified TPC.
unsigned int ChannelID_t
Type representing the ID of a readout channel.
bool fStoreDriftedElectronClusters
raw::ChannelID_t NearestChannel(geo::Point_t const &worldLoc, geo::PlaneID const &planeid) const
Returns the ID of the channel nearest to the specified position.
std::vector< double > fTransDiff2
int fMinNumberOfElCluster
const double * PlaneLocation(unsigned int p) const
Returns the coordinates of the center of the specified plane [cm].
std::vector< double > fLongDiff
cet::coded_exception< error, detail::translate > exception
Event finding and building.