Inheritance diagram for detsim::DriftElectronstoPlane:

Public Member Functions
	DriftElectronstoPlane (fhicl::ParameterSet const &pset)

void	produce (art::Event &evt) override

void	beginJob () override

Public Member Functions inherited from art::EDProducer
	EDProducer (fhicl::ParameterSet const &pset)

template<typename Config >
	EDProducer (Table< Config > const &config)

std::string	workerType () const

Public Member Functions inherited from art::detail::Producer
virtual	~Producer () noexcept

	Producer (fhicl::ParameterSet const &)

	Producer (Producer const &)=delete

	Producer (Producer &&)=delete

Producer &	operator= (Producer const &)=delete

Producer &	operator= (Producer &&)=delete

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)

Public Member Functions inherited from art::Modifier
	~Modifier () noexcept

	Modifier ()

	Modifier (Modifier const &)=delete

	Modifier (Modifier &&)=delete

Modifier &	operator= (Modifier const &)=delete

Modifier &	operator= (Modifier &&)=delete

Public Member Functions inherited from art::ModuleBase
virtual	~ModuleBase () noexcept

	ModuleBase ()

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)

template<typename T , BranchType BT>
ViewToken< T >	consumesView (InputTag const &tag)

template<typename T , BranchType BT>
ViewToken< T >	mayConsumeView (InputTag const &tag)

Private Attributes
art::InputTag	fSimModuleLabel

CLHEP::RandGauss	fRandGauss

bool	fStoreDriftedElectronClusters

double	fLongitudinalDiffusion

double	fTransverseDiffusion

double	fElectronClusterSize

int	fMinNumberOfElCluster

double	fGeVToElectrons

double	fRecombA

double	fRecombk

double	fModBoxA

double	fModBoxB

bool	fUseModBoxRecomb

double	fElectronLifetime

double	fLifetimeCorr_const

double	fLDiff_const

double	fTDiff_const

double	fRecipDriftVel [3]

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...

ISCalculationSeparate	fISAlg

Additional Inherited Members
Public Types inherited from art::EDProducer
using	ModuleType = EDProducer

using	WorkerType = WorkerT< EDProducer >

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

Public Types inherited from art::Modifier
template<typename UserConfig , typename UserKeysToIgnore = void>
using	Table = ProducerTable< UserConfig, detail::ModuleConfig, UserKeysToIgnore >

Static Public Member Functions inherited from art::EDProducer
static void	commitEvent (EventPrincipal &ep, Event &e)

Protected Member Functions inherited from art::ModuleBase
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 ()

Detailed Description

Definition at line 81 of file DriftElectronstoPlane_module.cc.

Constructor & Destructor Documentation

detsim::DriftElectronstoPlane::DriftElectronstoPlane ( fhicl::ParameterSet const & pset )

explicit

Definition at line 137 of file DriftElectronstoPlane_module.cc.

     : art::EDProducer{pset}
     , fSimModuleLabel{pset.get<art::InputTag>("SimulationLabel")}
     // create a default random engine; obtain the random seed from
     // NuRandomService, unless overridden in configuration with key
     // "Seed"
     , fRandGauss{art::ServiceHandle<rndm::NuRandomService>{}->createEngine(*this, pset, "Seed")}
     , fStoreDriftedElectronClusters{pset.get<bool>("StoreDriftedElectronClusters", true)}
     , fLongitudinalDiffusion{pset.get<double>("LongitudinalDiffusion", 6.2e-9)}
     , fTransverseDiffusion{pset.get<double>("TransverseDiffusion", 16.3e-9)}
     , fElectronClusterSize{pset.get<double>("ElectronClusterSize", 600.0)}
     , fMinNumberOfElCluster{pset.get<int>("MinNumberOfElCluster", 0)}
     , fGeVToElectrons{pset.get<double>("GeVToElectrons", 4.237e+07)}
     , fRecombA{pset.get<double>("RecombA", 0.800)}
     , fRecombk{pset.get<double>("Recombk", 0.0486)}
     , fModBoxA{pset.get<double>("ModBoxA", 0.930)}
     , fModBoxB{pset.get<double>("ModBoxB", 0.212)}
     , fUseModBoxRecomb{pset.get<bool>("UseModBoxRecomb", true)}
     , fISAlg{pset}
   {
     if (fStoreDriftedElectronClusters) { produces<std::vector<sim::SimDriftedElectronCluster>>(); }
   }

Member Function Documentation

void detsim::DriftElectronstoPlane::beginJob ( )

overridevirtual

Reimplemented from art::EDProducer.

Definition at line 162 of file DriftElectronstoPlane_module.cc.

   {
     // Define the physical constants we'll use.
 
     auto const detProp =
       art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataForJob();
     fElectronLifetime =
       detProp
         .ElectronLifetime(); // Electron lifetime as returned by the DetectorProperties service assumed to be in us;
     for (int i = 0; i < 3; ++i) {
       double driftVelocity =
         detProp.DriftVelocity(detProp.Efield(i),
                               detProp.Temperature()) *
         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;
 
       fRecipDriftVel[i] = 1. / driftVelocity;
     }
     MF_LOG_DEBUG("DriftElectronstoPlane")
       << " e lifetime (ns): " << fElectronLifetime
       << "\n Temperature (K): " << detProp.Temperature()
       << "\n Drift velocity (cm/ns): " << 1. / fRecipDriftVel[0] << " " << 1. / fRecipDriftVel[1]
       << " " << 1. / fRecipDriftVel[2];
 
     // Opposite of lifetime. Convert from us to standard LArSoft time units, ns;
     fLifetimeCorr_const = -1000. * fElectronLifetime;
     fLDiff_const = std::sqrt(2. * fLongitudinalDiffusion);
     fTDiff_const = std::sqrt(2. * fTransverseDiffusion);
 
     // For this detector's geometry, save the number of cryostats and
     // the number of TPCs within each cryostat.
     fNCryostats = fGeometry->Ncryostats();
     fNTPCs.resize(fNCryostats);
     for (size_t n = 0; n < fNCryostats; ++n)
       fNTPCs[n] = fGeometry->NTPC(n);
   }

void detsim::DriftElectronstoPlane::produce ( art::Event & evt )

overridevirtual

Todo:: think about effects of drift between planes.

Todo:: think about effects of drift between planes

Implements art::EDProducer.

Definition at line 200 of file DriftElectronstoPlane_module.cc.

   {
     // Fetch the SimEnergyDeposit objects for this event.
     typedef art::Handle<std::vector<sim::SimEnergyDeposit>> energyDepositHandle_t;
     energyDepositHandle_t energyDepositHandle;
     // If there aren't any energy deposits for this event, don't
     // panic. It's possible someone is doing a study with events
     // outside the TPC, or where there are only non-ionizing
     // particles, or something like that.
     if (!event.getByLabel(fSimModuleLabel, energyDepositHandle)) return;
     // Container for the SimDriftedElectronCluster objects
     std::unique_ptr<std::vector<sim::SimDriftedElectronCluster>>
       SimDriftedElectronClusterCollection(new std::vector<sim::SimDriftedElectronCluster>);
 
     // We're going through the input vector by index, rather than by
     // iterator, because we need the index number to compute the
     // associations near the end of this method.
     auto const& energyDeposits = *energyDepositHandle;
     auto energyDepositsSize = energyDeposits.size();
 
     auto const detProp =
       art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(event);
     // For each energy deposit in this event
     for (size_t edIndex = 0; edIndex < energyDepositsSize; ++edIndex) {
       auto const& energyDeposit = energyDeposits[edIndex];
 
       // "xyz" is the position of the energy deposit in world
       // coordinates. Note that the units of distance in
       // sim::SimEnergyDeposit are supposed to be cm.
       auto const mp = energyDeposit.MidPoint();
       double const xyz[3] = {mp.X(), mp.Y(), mp.Z()};
 
       // From the position in world coordinates, determine the
       // cryostat and tpc. If somehow the step is outside a tpc
       // (e.g., cosmic rays in rock) just move on to the next one.
       unsigned int cryostat = 0;
       unsigned int tpc = 0;
       const geo::TPCGeo& tpcGeo = fGeometry->TPC(tpc, cryostat);
 
       // The drift direction can be either in the positive
       // or negative direction in any coordinate x, y or z.
       // Charge drift in ...
       // +x: tpcGeo.DetectDriftDirection()==1
       // -x: tpcGeo.DetectDriftDirection()==-1
       // +y: tpcGeo.DetectDriftDirection()==2
       // -y tpcGeo.DetectDriftDirection()==-2
       // +z: tpcGeo.DetectDriftDirection()==3
       // -z: tpcGeo.DetectDriftDirection()==-3
 
       //Define charge drift direction: driftcoordinate (x, y or z) and driftsign (positive or negative). Also define coordinates perpendicular to drift direction.
       int driftcoordinate = std::abs(tpcGeo.DetectDriftDirection()) - 1; //x:0, y:1, z:2
 
       int transversecoordinate1 = 0;
       int transversecoordinate2 = 0;
       if (driftcoordinate == 0) {
         transversecoordinate1 = 1;
         transversecoordinate2 = 2;
       }
       else if (driftcoordinate == 1) {
         transversecoordinate1 = 0;
         transversecoordinate2 = 2;
       }
       else if (driftcoordinate == 2) {
         transversecoordinate1 = 0;
         transversecoordinate2 = 1;
       }
 
       if (transversecoordinate1 == transversecoordinate2)
         continue; //this is the case when driftcoordinate != 0, 1 or 2
 
       int driftsign = 0; //1: +x, +y or +z, -1: -x, -y or -z
       if (tpcGeo.DetectDriftDirection() > 0)
         driftsign = 1;
       else
         driftsign = -1;
 
       //Check for charge deposits behind charge readout planes
       if (driftsign == 1 && tpcGeo.PlaneLocation(0)[driftcoordinate] < xyz[driftcoordinate])
         continue;
       if (driftsign == -1 && tpcGeo.PlaneLocation(0)[driftcoordinate] > xyz[driftcoordinate])
         continue;
 
       /// \todo think about effects of drift between planes.
       // Center of plane is also returned in cm units
       double DriftDistance =
         std::abs(xyz[driftcoordinate] - tpcGeo.PlaneLocation(0)[driftcoordinate]);
 
       // Space-charge effect (SCE): Get SCE {x,y,z} offsets for
       // particular location in TPC
       geo::Vector_t posOffsets{0.0, 0.0, 0.0};
       double posOffsetxyz[3] = {
         0.0, 0.0, 0.0}; //need this array for the driftcoordinate and transversecoordinates
       auto const* SCE = lar::providerFrom<spacecharge::SpaceChargeService>();
       if (SCE->EnableSimSpatialSCE() == true) {
         posOffsets = SCE->GetPosOffsets(mp);
         if (larsim::Utils::SCE::out_of_bounds(posOffsets)) continue;
         posOffsetxyz[0] = posOffsets.X();
         posOffsetxyz[1] = posOffsets.Y();
         posOffsetxyz[2] = posOffsets.Z();
       }
 
       double avegagetransversePos1 = 0.;
       double avegagetransversePos2 = 0.;
 
       DriftDistance += -1. * posOffsetxyz[driftcoordinate];
       avegagetransversePos1 = xyz[transversecoordinate1] + posOffsetxyz[transversecoordinate1];
       avegagetransversePos2 = xyz[transversecoordinate2] + posOffsetxyz[transversecoordinate2];
 
       // Space charge distortion could push the energy deposit beyond the wire
       // plane (see issue #15131). Given that we don't have any subtlety in the
       // simulation of this region, bringing the deposit exactly on the plane
       // should be enough for the time being.
       if (DriftDistance < 0.) DriftDistance = 0.;
 
       // Drift time in ns
       double TDrift = DriftDistance * fRecipDriftVel[0];
 
       if (
         tpcGeo.Nplanes() == 2 &&
         driftcoordinate ==
           0) { // special case for ArgoNeuT (Nplanes = 2 and drift direction = x): plane 0 is the second wire plane
         TDrift = ((DriftDistance - tpcGeo.PlanePitch(0, 1)) * fRecipDriftVel[0] +
                   tpcGeo.PlanePitch(0, 1) * fRecipDriftVel[1]);
       }
       const int nIonizedElectrons =
         fISAlg.CalculateIonizationAndScintillation(detProp, energyDeposit).numElectrons;
 
       const double lifetimecorrection = TMath::Exp(TDrift / fLifetimeCorr_const);
       const double energy = energyDeposit.Energy();
 
       // if we have no electrons (too small energy or too large recombination)
       // we are done already here
       if (nIonizedElectrons <= 0) {
         MF_LOG_DEBUG("DriftElectronstoPlane")
           << "step " // << energyDeposit << "\n"
           << "No electrons drifted to readout, " << energy << " MeV lost.";
         continue;
       }
 
       // includes the effect of lifetime: lifetimecorrection = exp[-tdrift/tau]
       const double nElectrons = nIonizedElectrons * lifetimecorrection;
 
       // Longitudinal & transverse diffusion sigma (cm)
       double SqrtT = std::sqrt(TDrift);
       double LDiffSig = SqrtT * fLDiff_const;
       double TDiffSig = SqrtT * fTDiff_const;
       double electronclsize = fElectronClusterSize;
 
       // Number of electron clusters.
       int nClus = (int)std::ceil(nElectrons / electronclsize);
       if (nClus < fMinNumberOfElCluster) {
         electronclsize = nElectrons / fMinNumberOfElCluster;
         if (electronclsize < 1.0) { electronclsize = 1.0; }
         nClus = (int)std::ceil(nElectrons / electronclsize);
       }
 
       // Empty and resize the electron-cluster vectors.
       fLongDiff.clear();
       fTransDiff1.clear();
       fTransDiff2.clear();
       fnElDiff.clear();
       fnEnDiff.clear();
       fLongDiff.resize(nClus);
       fTransDiff1.resize(nClus);
       fTransDiff2.resize(nClus);
       fnElDiff.resize(nClus, electronclsize);
       fnEnDiff.resize(nClus);
 
       // fix the number of electrons in the last cluster, that has a smaller size
       fnElDiff.back() = nElectrons - (nClus - 1) * electronclsize;
 
       for (size_t xx = 0; xx < fnElDiff.size(); ++xx) {
         if (nElectrons > 0)
           fnEnDiff[xx] = energy / nElectrons * fnElDiff[xx];
         else
           fnEnDiff[xx] = 0.;
       }
 
       // Smear drift times by longitudinal diffusion
       if (LDiffSig > 0.0)
         fRandGauss.fireArray(nClus, &fLongDiff[0], 0., LDiffSig);
       else
         fLongDiff.assign(nClus, 0.0);
 
       if (TDiffSig > 0.0) {
         // Smear the coordinates in plane perpendicular to drift direction by the transverse diffusion
         fRandGauss.fireArray(nClus, &fTransDiff1[0], avegagetransversePos1, TDiffSig);
         fRandGauss.fireArray(nClus, &fTransDiff2[0], avegagetransversePos2, TDiffSig);
       }
       else {
         fTransDiff1.assign(nClus, avegagetransversePos1);
         fTransDiff2.assign(nClus, avegagetransversePos2);
       }
 
       // make a collection of electrons for each plane
       for (size_t p = 0; p < tpcGeo.Nplanes(); ++p) {
 
         fDriftClusterPos[driftcoordinate] = tpcGeo.PlaneLocation(p)[driftcoordinate];
 
         // Drift nClus electron clusters to the induction plane
         for (int k = 0; k < nClus; ++k) {
 
           // Correct drift time for longitudinal diffusion and plane
           double TDiff = TDrift + fLongDiff[k] * fRecipDriftVel[0];
 
           // Take into account different Efields between planes
           // Also take into account special case for ArgoNeuT (Nplanes = 2 and drift direction = x): plane 0 is the second wire plane
           for (size_t ip = 0; ip < p; ++ip) {
             TDiff +=
               (tpcGeo.PlaneLocation(ip + 1)[driftcoordinate] -
                tpcGeo.PlaneLocation(ip)[driftcoordinate]) *
               fRecipDriftVel[(tpcGeo.Nplanes() == 2 && driftcoordinate == 0) ? ip + 2 : ip + 1];
           }
 
           fDriftClusterPos[transversecoordinate1] = fTransDiff1[k];
           fDriftClusterPos[transversecoordinate2] = fTransDiff2[k];
           auto const simTime = energyDeposit.Time();
           /// \todo think about effects of drift between planes
           SimDriftedElectronClusterCollection->emplace_back(
             fnElDiff[k],
             TDiff + simTime,                      // timing
             geo::Point_t{mp.X(), mp.Y(), mp.Z()}, // mean position of the deposited energy
             geo::Point_t{fDriftClusterPos[0],
                          fDriftClusterPos[1],
                          fDriftClusterPos[2]}, // final position of the drifted cluster
             geo::Point_t{
               LDiffSig, TDiffSig, TDiffSig}, // Longitudinal (X) and transverse (Y,Z) diffusion
             fnEnDiff[k],                     //deposited energy that originated this cluster
             energyDeposit.TrackID());
 
         } // end loop over clusters
       }   // end loop over planes
     }     // for each sim::SimEnergyDeposit
 
     if (fStoreDriftedElectronClusters) event.put(std::move(SimDriftedElectronClusterCollection));
   }