76 #include "CLHEP/Random/RandGauss.h" 150 -> registerAndSeedEngine(
createEngine(0), pset,
"Seed")}
157 ,
fRecombA{pset.get<
double>(
"RecombA", 0.800)}
158 ,
fRecombk{pset.get<
double>(
"Recombk", 0.0486)}
159 ,
fModBoxA{pset.get<
double>(
"ModBoxA", 0.930)}
160 ,
fModBoxB{pset.get<
double>(
"ModBoxB", 0.212)}
177 for (
int i = 0; i < 3; ++i) {
178 double driftVelocity =
179 detProp.DriftVelocity(detProp.Efield(i),
180 detProp.Temperature()) *
187 <<
"\n Temperature (K): " << detProp.Temperature()
209 energyDepositHandle_t energyDepositHandle;
216 std::unique_ptr<std::vector<sim::SimDriftedElectronCluster>>
217 SimDriftedElectronClusterCollection(
new std::vector<sim::SimDriftedElectronCluster>);
222 auto const& energyDeposits = *energyDepositHandle;
223 auto energyDepositsSize = energyDeposits.size();
228 for (
size_t edIndex = 0; edIndex < energyDepositsSize; ++edIndex) {
229 auto const& energyDeposit = energyDeposits[edIndex];
234 auto const mp = energyDeposit.MidPoint();
235 double const xyz[3] = {mp.X(), mp.Y(), mp.Z()};
255 int transversecoordinate1 = 0;
256 int transversecoordinate2 = 0;
257 if (driftcoordinate == 0) {
258 transversecoordinate1 = 1;
259 transversecoordinate2 = 2;
261 else if (driftcoordinate == 1) {
262 transversecoordinate1 = 0;
263 transversecoordinate2 = 2;
265 else if (driftcoordinate == 2) {
266 transversecoordinate1 = 0;
267 transversecoordinate2 = 1;
270 if (transversecoordinate1 == transversecoordinate2)
281 auto const plane_center_coord =
geo::vect::coord(plane_center, driftcoordinate);
282 if (driftsign == 1 && plane_center_coord < xyz[driftcoordinate])
continue;
283 if (driftsign == -1 && plane_center_coord > xyz[driftcoordinate])
continue;
287 double DriftDistance =
std::abs(xyz[driftcoordinate] - plane_center_coord);
292 double posOffsetxyz[3] = {
294 auto const* SCE = lar::providerFrom<spacecharge::SpaceChargeService>();
295 if (SCE->EnableSimSpatialSCE() ==
true) {
296 posOffsets = SCE->GetPosOffsets(mp);
298 posOffsetxyz[0] = posOffsets.X();
299 posOffsetxyz[1] = posOffsets.Y();
300 posOffsetxyz[2] = posOffsets.Z();
303 double avegagetransversePos1 = 0.;
304 double avegagetransversePos2 = 0.;
306 DriftDistance += -1. * posOffsetxyz[driftcoordinate];
307 avegagetransversePos1 = xyz[transversecoordinate1] + posOffsetxyz[transversecoordinate1];
308 avegagetransversePos2 = xyz[transversecoordinate2] + posOffsetxyz[transversecoordinate2];
314 if (DriftDistance < 0.) DriftDistance = 0.;
323 TDrift = ((DriftDistance - tpcGeo.
PlanePitch(0, 1)) * fRecipDriftVel[0] +
326 const int nIonizedElectrons =
330 const double energy = energyDeposit.Energy();
334 if (nIonizedElectrons <= 0) {
337 <<
"No electrons drifted to readout, " << energy <<
" MeV lost.";
342 const double nElectrons = nIonizedElectrons * lifetimecorrection;
345 double SqrtT = std::sqrt(TDrift);
351 int nClus = (int)std::ceil(nElectrons / electronclsize);
354 if (electronclsize < 1.0) { electronclsize = 1.0; }
355 nClus = (int)std::ceil(nElectrons / electronclsize);
367 fnElDiff.resize(nClus, electronclsize);
371 fnElDiff.back() = nElectrons - (nClus - 1) * electronclsize;
386 if (TDiffSig > 0.0) {
397 for (
size_t p = 0; p < tpcGeo.
Nplanes(); ++p) {
402 for (
int k = 0; k < nClus; ++k) {
405 double TDiff = TDrift +
fLongDiff[k] * fRecipDriftVel[0];
409 for (
size_t ip = 0; ip < p; ++ip) {
415 fRecipDriftVel[(tpcGeo.
Nplanes() == 2 && driftcoordinate == 0) ? ip + 2 : ip + 1];
420 auto const simTime = energyDeposit.Time();
422 SimDriftedElectronClusterCollection->emplace_back(
428 fDriftClusterPos[2]},
430 LDiffSig, TDiffSig, TDiffSig},
432 energyDeposit.TrackID());
base_engine_t & createEngine(seed_t seed)
unsigned int NTPC(CryostatID const &cryoid=cryostat_zero) const
Returns the total number of TPCs in the specified cryostat.
double PlanePitch(unsigned int p1=0, unsigned int p2=1) const
Returns the center of the TPC volume in world coordinates [cm].
Utilities related to art service access.
contains objects relating to SimDriftedElectronCluster
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].
std::vector< double > fnElDiff
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.
art::ServiceHandle< geo::Geometry const > fGeometry
Handle to the Geometry service.
EDProducer(fhicl::ParameterSet const &pset)
Geometry information for a single TPC.
Detector simulation of raw signals on wires.
constexpr auto abs(T v)
Returns the absolute value of the argument.
double fLifetimeCorr_const
std::vector< double > fTransDiff1
unsigned int Ncryostats() const
Returns the number of cryostats in the detector.
std::vector< size_t > fNTPCs
double fElectronClusterSize
TPCGeo const & TPC(TPCID const &tpcid=tpc_zero) const
Returns the specified TPC.
bool fStoreDriftedElectronClusters
double fLongitudinalDiffusion
double fDriftClusterPos[3]
#define DEFINE_ART_MODULE(klass)
CLHEP::RandGauss fRandGauss
Utility function for testing if Space Charge offsets are out of bounds.
Data CalculateIonizationAndScintillation(detinfo::DetectorPropertiesData const &detProp, sim::SimEnergyDeposit const &edep) const
An art service to assist in the distribution of guaranteed unique seeds to all engines within an art ...
std::vector< double > fnEnDiff
double fTransverseDiffusion
void produce(art::Event &evt) override
ISCalculationSeparate fISAlg
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
art::InputTag fSimModuleLabel
short int DetectDriftDirection() const
Returns the expected drift direction based on geometry.
contains information for a single step in the detector simulation
DriftElectronstoPlane(fhicl::ParameterSet const &pset)
std::vector< double > fTransDiff2
PlaneGeo const & Plane(geo::View_t view) const
Return the plane in the tpc with View_t view.
int fMinNumberOfElCluster
art framework interface to geometry description
Event finding and building.
The data type to uniquely identify a cryostat.
std::vector< double > fLongDiff