LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
DriftElectronstoPlane_module.cc
Go to the documentation of this file.
1 
55 // LArSoft includes
64 
65 // Framework includes
71 #include "fhiclcpp/ParameterSet.h"
74 
75 // External libraries
76 #include "CLHEP/Random/RandGauss.h"
77 #include "TMath.h"
78 
79 #include <cmath>
80 #include <memory>
81 #include <utility>
82 #include <vector>
83 
84 namespace detsim {
85 
86  // Base class for creation of raw signals on wires.
88  public:
89  explicit DriftElectronstoPlane(fhicl::ParameterSet const& pset);
90 
91  // Methods that that are available for a module derived from
92  // art::EDProducer.
93  void produce(art::Event& evt) override;
94  void beginJob() override;
95 
96  private:
97  // The label of the module that created the sim::SimEnergyDeposit
98  // objects (as of Oct-2017, this is probably "largeant").
100 
101  CLHEP::RandGauss fRandGauss;
102 
109  double fRecombA;
110  double fRecombk;
111  double fModBoxA;
112  double fModBoxB;
114 
117  double fLDiff_const;
118  double fTDiff_const;
119  double fRecipDriftVel[3];
120 
121  // Save the number of cryostats, and the number of TPCs within
122  // each cryostat.
123  size_t fNCryostats;
124  std::vector<size_t> fNTPCs;
125 
126  // Per-cluster information.
127  std::vector<double> fLongDiff;
128  std::vector<double> fTransDiff1;
129  std::vector<double> fTransDiff2;
130  std::vector<double> fnElDiff;
131  std::vector<double> fnEnDiff;
132 
133  double fDriftClusterPos[3];
134 
136 
137  //IS calculationg
139 
140  }; // class DriftElectronstoPlane
141 
142  //-------------------------------------------------
144  : art::EDProducer{pset}
145  , fSimModuleLabel{pset.get<art::InputTag>("SimulationLabel")}
146  // create a default random engine; obtain the random seed from
147  // NuRandomService, unless overridden in configuration with key
148  // "Seed"
150  -> registerAndSeedEngine(createEngine(0), pset, "Seed")}
151  , fStoreDriftedElectronClusters{pset.get<bool>("StoreDriftedElectronClusters", true)}
152  , fLongitudinalDiffusion{pset.get<double>("LongitudinalDiffusion", 6.2e-9)}
153  , fTransverseDiffusion{pset.get<double>("TransverseDiffusion", 16.3e-9)}
154  , fElectronClusterSize{pset.get<double>("ElectronClusterSize", 600.0)}
155  , fMinNumberOfElCluster{pset.get<int>("MinNumberOfElCluster", 0)}
156  , fGeVToElectrons{pset.get<double>("GeVToElectrons", 4.237e+07)}
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)}
161  , fUseModBoxRecomb{pset.get<bool>("UseModBoxRecomb", true)}
162  , fISAlg{pset}
163  {
164  if (fStoreDriftedElectronClusters) { produces<std::vector<sim::SimDriftedElectronCluster>>(); }
165  }
166 
167  //-------------------------------------------------
169  {
170  // Define the physical constants we'll use.
171 
172  auto const detProp =
175  detProp
176  .ElectronLifetime(); // Electron lifetime as returned by the DetectorProperties service assumed to be in us;
177  for (int i = 0; i < 3; ++i) {
178  double driftVelocity =
179  detProp.DriftVelocity(detProp.Efield(i),
180  detProp.Temperature()) *
181  1.e-3; // Drift velocity as returned by the DetectorProperties service assumed to be in cm/us. Multiply by 1.e-3 to convert into LArSoft standard velocity units, cm/ns;
182 
183  fRecipDriftVel[i] = 1. / driftVelocity;
184  }
185  MF_LOG_DEBUG("DriftElectronstoPlane")
186  << " e lifetime (ns): " << fElectronLifetime
187  << "\n Temperature (K): " << detProp.Temperature()
188  << "\n Drift velocity (cm/ns): " << 1. / fRecipDriftVel[0] << " " << 1. / fRecipDriftVel[1]
189  << " " << 1. / fRecipDriftVel[2];
190 
191  // Opposite of lifetime. Convert from us to standard LArSoft time units, ns;
193  fLDiff_const = std::sqrt(2. * fLongitudinalDiffusion);
194  fTDiff_const = std::sqrt(2. * fTransverseDiffusion);
195 
196  // For this detector's geometry, save the number of cryostats and
197  // the number of TPCs within each cryostat.
199  fNTPCs.resize(fNCryostats);
200  for (size_t n = 0; n < fNCryostats; ++n)
202  }
203 
204  //-------------------------------------------------
206  {
207  // Fetch the SimEnergyDeposit objects for this event.
208  typedef art::Handle<std::vector<sim::SimEnergyDeposit>> energyDepositHandle_t;
209  energyDepositHandle_t energyDepositHandle;
210  // If there aren't any energy deposits for this event, don't
211  // panic. It's possible someone is doing a study with events
212  // outside the TPC, or where there are only non-ionizing
213  // particles, or something like that.
214  if (!event.getByLabel(fSimModuleLabel, energyDepositHandle)) return;
215  // Container for the SimDriftedElectronCluster objects
216  std::unique_ptr<std::vector<sim::SimDriftedElectronCluster>>
217  SimDriftedElectronClusterCollection(new std::vector<sim::SimDriftedElectronCluster>);
218 
219  // We're going through the input vector by index, rather than by
220  // iterator, because we need the index number to compute the
221  // associations near the end of this method.
222  auto const& energyDeposits = *energyDepositHandle;
223  auto energyDepositsSize = energyDeposits.size();
224 
225  auto const detProp =
227  // For each energy deposit in this event
228  for (size_t edIndex = 0; edIndex < energyDepositsSize; ++edIndex) {
229  auto const& energyDeposit = energyDeposits[edIndex];
230 
231  // "xyz" is the position of the energy deposit in world
232  // coordinates. Note that the units of distance in
233  // sim::SimEnergyDeposit are supposed to be cm.
234  auto const mp = energyDeposit.MidPoint();
235  double const xyz[3] = {mp.X(), mp.Y(), mp.Z()};
236 
237  // From the position in world coordinates, determine the
238  // cryostat and tpc. If somehow the step is outside a tpc
239  // (e.g., cosmic rays in rock) just move on to the next one.
240  const geo::TPCGeo& tpcGeo = fGeometry->TPC();
241 
242  // The drift direction can be either in the positive
243  // or negative direction in any coordinate x, y or z.
244  // Charge drift in ...
245  // +x: tpcGeo.DetectDriftDirection()==1
246  // -x: tpcGeo.DetectDriftDirection()==-1
247  // +y: tpcGeo.DetectDriftDirection()==2
248  // -y tpcGeo.DetectDriftDirection()==-2
249  // +z: tpcGeo.DetectDriftDirection()==3
250  // -z: tpcGeo.DetectDriftDirection()==-3
251 
252  //Define charge drift direction: driftcoordinate (x, y or z) and driftsign (positive or negative). Also define coordinates perpendicular to drift direction.
253  int driftcoordinate = std::abs(tpcGeo.DetectDriftDirection()) - 1; //x:0, y:1, z:2
254 
255  int transversecoordinate1 = 0;
256  int transversecoordinate2 = 0;
257  if (driftcoordinate == 0) {
258  transversecoordinate1 = 1;
259  transversecoordinate2 = 2;
260  }
261  else if (driftcoordinate == 1) {
262  transversecoordinate1 = 0;
263  transversecoordinate2 = 2;
264  }
265  else if (driftcoordinate == 2) {
266  transversecoordinate1 = 0;
267  transversecoordinate2 = 1;
268  }
269 
270  if (transversecoordinate1 == transversecoordinate2)
271  continue; //this is the case when driftcoordinate != 0, 1 or 2
272 
273  int driftsign = 0; //1: +x, +y or +z, -1: -x, -y or -z
274  if (tpcGeo.DetectDriftDirection() > 0)
275  driftsign = 1;
276  else
277  driftsign = -1;
278 
279  //Check for charge deposits behind charge readout planes
280  auto const plane_center = tpcGeo.Plane(0).GetCenter();
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;
284 
286  // Center of plane is also returned in cm units
287  double DriftDistance = std::abs(xyz[driftcoordinate] - plane_center_coord);
288 
289  // Space-charge effect (SCE): Get SCE {x,y,z} offsets for
290  // particular location in TPC
291  geo::Vector_t posOffsets{0.0, 0.0, 0.0};
292  double posOffsetxyz[3] = {
293  0.0, 0.0, 0.0}; //need this array for the driftcoordinate and transversecoordinates
294  auto const* SCE = lar::providerFrom<spacecharge::SpaceChargeService>();
295  if (SCE->EnableSimSpatialSCE() == true) {
296  posOffsets = SCE->GetPosOffsets(mp);
297  if (larsim::Utils::SCE::out_of_bounds(posOffsets)) continue;
298  posOffsetxyz[0] = posOffsets.X();
299  posOffsetxyz[1] = posOffsets.Y();
300  posOffsetxyz[2] = posOffsets.Z();
301  }
302 
303  double avegagetransversePos1 = 0.;
304  double avegagetransversePos2 = 0.;
305 
306  DriftDistance += -1. * posOffsetxyz[driftcoordinate];
307  avegagetransversePos1 = xyz[transversecoordinate1] + posOffsetxyz[transversecoordinate1];
308  avegagetransversePos2 = xyz[transversecoordinate2] + posOffsetxyz[transversecoordinate2];
309 
310  // Space charge distortion could push the energy deposit beyond the wire
311  // plane (see issue #15131). Given that we don't have any subtlety in the
312  // simulation of this region, bringing the deposit exactly on the plane
313  // should be enough for the time being.
314  if (DriftDistance < 0.) DriftDistance = 0.;
315 
316  // Drift time in ns
317  double TDrift = DriftDistance * fRecipDriftVel[0];
318 
319  if (
320  tpcGeo.Nplanes() == 2 &&
321  driftcoordinate ==
322  0) { // special case for ArgoNeuT (Nplanes = 2 and drift direction = x): plane 0 is the second wire plane
323  TDrift = ((DriftDistance - tpcGeo.PlanePitch(0, 1)) * fRecipDriftVel[0] +
324  tpcGeo.PlanePitch(0, 1) * fRecipDriftVel[1]);
325  }
326  const int nIonizedElectrons =
328 
329  const double lifetimecorrection = TMath::Exp(TDrift / fLifetimeCorr_const);
330  const double energy = energyDeposit.Energy();
331 
332  // if we have no electrons (too small energy or too large recombination)
333  // we are done already here
334  if (nIonizedElectrons <= 0) {
335  MF_LOG_DEBUG("DriftElectronstoPlane")
336  << "step " // << energyDeposit << "\n"
337  << "No electrons drifted to readout, " << energy << " MeV lost.";
338  continue;
339  }
340 
341  // includes the effect of lifetime: lifetimecorrection = exp[-tdrift/tau]
342  const double nElectrons = nIonizedElectrons * lifetimecorrection;
343 
344  // Longitudinal & transverse diffusion sigma (cm)
345  double SqrtT = std::sqrt(TDrift);
346  double LDiffSig = SqrtT * fLDiff_const;
347  double TDiffSig = SqrtT * fTDiff_const;
348  double electronclsize = fElectronClusterSize;
349 
350  // Number of electron clusters.
351  int nClus = (int)std::ceil(nElectrons / electronclsize);
352  if (nClus < fMinNumberOfElCluster) {
353  electronclsize = nElectrons / fMinNumberOfElCluster;
354  if (electronclsize < 1.0) { electronclsize = 1.0; }
355  nClus = (int)std::ceil(nElectrons / electronclsize);
356  }
357 
358  // Empty and resize the electron-cluster vectors.
359  fLongDiff.clear();
360  fTransDiff1.clear();
361  fTransDiff2.clear();
362  fnElDiff.clear();
363  fnEnDiff.clear();
364  fLongDiff.resize(nClus);
365  fTransDiff1.resize(nClus);
366  fTransDiff2.resize(nClus);
367  fnElDiff.resize(nClus, electronclsize);
368  fnEnDiff.resize(nClus);
369 
370  // fix the number of electrons in the last cluster, that has a smaller size
371  fnElDiff.back() = nElectrons - (nClus - 1) * electronclsize;
372 
373  for (size_t xx = 0; xx < fnElDiff.size(); ++xx) {
374  if (nElectrons > 0)
375  fnEnDiff[xx] = energy / nElectrons * fnElDiff[xx];
376  else
377  fnEnDiff[xx] = 0.;
378  }
379 
380  // Smear drift times by longitudinal diffusion
381  if (LDiffSig > 0.0)
382  fRandGauss.fireArray(nClus, &fLongDiff[0], 0., LDiffSig);
383  else
384  fLongDiff.assign(nClus, 0.0);
385 
386  if (TDiffSig > 0.0) {
387  // Smear the coordinates in plane perpendicular to drift direction by the transverse diffusion
388  fRandGauss.fireArray(nClus, &fTransDiff1[0], avegagetransversePos1, TDiffSig);
389  fRandGauss.fireArray(nClus, &fTransDiff2[0], avegagetransversePos2, TDiffSig);
390  }
391  else {
392  fTransDiff1.assign(nClus, avegagetransversePos1);
393  fTransDiff2.assign(nClus, avegagetransversePos2);
394  }
395 
396  // make a collection of electrons for each plane
397  for (size_t p = 0; p < tpcGeo.Nplanes(); ++p) {
398  auto const plane_center = tpcGeo.Plane(p).GetCenter();
399  fDriftClusterPos[driftcoordinate] = geo::vect::coord(plane_center, driftcoordinate);
400 
401  // Drift nClus electron clusters to the induction plane
402  for (int k = 0; k < nClus; ++k) {
403 
404  // Correct drift time for longitudinal diffusion and plane
405  double TDiff = TDrift + fLongDiff[k] * fRecipDriftVel[0];
406 
407  // Take into account different Efields between planes
408  // Also take into account special case for ArgoNeuT (Nplanes = 2 and drift direction = x): plane 0 is the second wire plane
409  for (size_t ip = 0; ip < p; ++ip) {
410  auto const plane_center = tpcGeo.Plane(ip).GetCenter();
411  auto const next_plane_center = tpcGeo.Plane(ip + 1).GetCenter();
412  TDiff +=
413  (geo::vect::coord(next_plane_center, driftcoordinate) -
414  geo::vect::coord(plane_center, driftcoordinate)) *
415  fRecipDriftVel[(tpcGeo.Nplanes() == 2 && driftcoordinate == 0) ? ip + 2 : ip + 1];
416  }
417 
418  fDriftClusterPos[transversecoordinate1] = fTransDiff1[k];
419  fDriftClusterPos[transversecoordinate2] = fTransDiff2[k];
420  auto const simTime = energyDeposit.Time();
422  SimDriftedElectronClusterCollection->emplace_back(
423  fnElDiff[k],
424  TDiff + simTime, // timing
425  geo::Point_t{mp.X(), mp.Y(), mp.Z()}, // mean position of the deposited energy
427  fDriftClusterPos[1],
428  fDriftClusterPos[2]}, // final position of the drifted cluster
429  geo::Point_t{
430  LDiffSig, TDiffSig, TDiffSig}, // Longitudinal (X) and transverse (Y,Z) diffusion
431  fnEnDiff[k], //deposited energy that originated this cluster
432  energyDeposit.TrackID());
433 
434  } // end loop over clusters
435  } // end loop over planes
436  } // for each sim::SimEnergyDeposit
437 
438  if (fStoreDriftedElectronClusters) event.put(std::move(SimDriftedElectronClusterCollection));
439  }
440 } // namespace detsim
441 
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.
Definition: GeometryCore.h:686
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
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.
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
art::ServiceHandle< geo::Geometry const > fGeometry
Handle to the Geometry service.
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.cc:6
Geometry information for a single TPC.
Definition: TPCGeo.h:36
Detector simulation of raw signals on wires.
constexpr auto abs(T v)
Returns the absolute value of the argument.
unsigned int Ncryostats() const
Returns the number of cryostats in the detector.
Definition: GeometryCore.h:430
TPCGeo const & TPC(TPCID const &tpcid=tpc_zero) const
Returns the specified TPC.
Definition: GeometryCore.h:722
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
Utility function for testing if Space Charge offsets are out of bounds.
double energy
Definition: plottest35.C:25
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 ...
void produce(art::Event &evt) override
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
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
short int DetectDriftDirection() const
Returns the expected drift direction based on geometry.
Definition: TPCGeo.cxx:149
contains information for a single step in the detector simulation
DriftElectronstoPlane(fhicl::ParameterSet const &pset)
#define MF_LOG_DEBUG(id)
Definition: MVAAlg.h:12
Char_t n[5]
PlaneGeo const & Plane(geo::View_t view) const
Return the plane in the tpc with View_t view.
Definition: TPCGeo.cxx:252
TCEvent evt
Definition: DataStructs.cxx:8
Float_t e
Definition: plot.C:35
art framework interface to geometry description
Event finding and building.
The data type to uniquely identify a cryostat.
Definition: geo_types.h:192