Inheritance diagram for detsim::SimDriftElectrons:

Classes
struct	ChannelBookKeeping

Public Member Functions
	SimDriftElectrons (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 Types
typedef std::map< raw::ChannelID_t, ChannelBookKeeping >	ChannelMap_t

Private Attributes
art::InputTag	fSimModuleLabel

CLHEP::RandGauss	fRandGauss

double	fElectronLifetime

double	fElectronClusterSize

int	fMinNumberOfElCluster

double	fLongitudinalDiffusion

double	fTransverseDiffusion

double	fLifetimeCorr_const

double	fLDiff_const

double	fTDiff_const

double	fRecipDriftVel [3]

bool	fStoreDriftedElectronClusters

std::vector< ChannelMap_t >	fChannelMaps

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

larg4::ISCalcSeparate	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 94 of file SimDriftElectrons_module.cc.

Member Typedef Documentation

typedef std::map<raw::ChannelID_t, ChannelBookKeeping> detsim::SimDriftElectrons::ChannelMap_t

private

Definition at line 134 of file SimDriftElectrons_module.cc.

Constructor & Destructor Documentation

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

explicit

Definition at line 163 of file SimDriftElectrons_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", false)}
   {
     produces<std::vector<sim::SimChannel>>();
     if (fStoreDriftedElectronClusters) { produces<std::vector<sim::SimDriftedElectronCluster>>(); }
   }

Member Function Documentation

void detsim::SimDriftElectrons::beginJob ( )

overridevirtual

Reimplemented from art::EDProducer.

Definition at line 178 of file SimDriftElectrons_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;
     }
 
     // To-do: Move the parameters we fetch from "LArG4" to detector
     // properties.
     art::ServiceHandle<sim::LArG4Parameters const> paramHandle;
     fElectronClusterSize = paramHandle->ElectronClusterSize();
     fMinNumberOfElCluster = paramHandle->MinNumberOfElCluster();
     fLongitudinalDiffusion = paramHandle->LongitudinalDiffusion(); // cm^2/ns units
     fTransverseDiffusion = paramHandle->TransverseDiffusion();     // cm^2/ns units
 
     MF_LOG_DEBUG("SimDriftElectrons")
       << " 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::SimDriftElectrons::produce ( art::Event & evt )

overridevirtual

Todo:: think about effects of drift between planes.

Todo:: think about effects of drift between planes

Todo:: check on what happens if we allow the tdc value to be

Todo:: beyond the end of the expected number of ticks

Implements art::EDProducer.

Definition at line 225 of file SimDriftElectrons_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;
 
     // Define the container for the SimChannel objects that will be
     // transferred to the art::Event after the put statement below.
     std::unique_ptr<std::vector<sim::SimChannel>> channels(new std::vector<sim::SimChannel>);
     // Container for the SimDriftedElectronCluster objects
     std::unique_ptr<std::vector<sim::SimDriftedElectronCluster>>
       SimDriftedElectronClusterCollection(new std::vector<sim::SimDriftedElectronCluster>);
 
     // Clear the channel maps from the last event. Remember,
     // fChannelMaps is an array[cryo][tpc] of maps.
     size_t cryo = 0;
     fChannelMaps.resize(fNCryostats);
     for (auto& cryoData : fChannelMaps) { // each, a vector of maps
       cryoData.clear();
     }
 
     auto const clockData =
       art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(event);
     auto const& tpcClock = clockData.TPCClock();
 
     auto const detProp =
       art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(event, clockData);
     // 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();
 
     // 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;
       try {
         fGeometry->PositionToCryostat(xyz, cryostat);
       }
       catch (cet::exception& e) {
         mf::LogWarning("SimDriftElectrons") << "step " // << energyDeposit << "\n"
                                             << "cannot be found in a cryostat\n"
                                             << e;
         continue;
       }
 
       unsigned int tpc = 0;
       try {
         fGeometry->PositionToTPC(xyz, tpc, cryostat);
       }
       catch (cet::exception& e) {
         mf::LogWarning("SimDriftElectrons") << "step " // << energyDeposit << "\n"
                                             << "cannot be found in a TPC\n"
                                             << e;
         continue;
       }
 
       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.CalcIonAndScint(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("SimDriftElectrons")
           << "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.PlanePitch(ip+1, ip) *
                 fRecipDriftVel[(tpcGeo.Nplanes() == 2 && driftcoordinate == 0) ? ip + 2 : ip + 1];
           }
 
           fDriftClusterPos[transversecoordinate1] = fTransDiff1[k];
           fDriftClusterPos[transversecoordinate2] = fTransDiff2[k];
 
           /// \todo think about effects of drift between planes
 
           // grab the nearest channel to the fDriftClusterPos position
           try {
             raw::ChannelID_t channel =
               fGeometry->NearestChannel(fDriftClusterPos, p, tpc, cryostat);
 
             /// \todo check on what happens if we allow the tdc value to be
             /// \todo beyond the end of the expected number of ticks
             // Add potential decay/capture/etc delay effect, simTime.
             auto const simTime = energyDeposit.Time();
             unsigned int tdc = tpcClock.Ticks(clockData.G4ToElecTime(TDiff + simTime));
 
             // Find whether we already have this channel in our map.
             ChannelMap_t& channelDataMap = fChannelMaps[cryostat];
             auto search = channelDataMap.find(channel);
 
             // We will find (or create) the pointer to a
             // sim::SimChannel.
             size_t channelIndex = 0;
 
             // Have we created the sim::SimChannel corresponding to
             // channel ID?
             if (search == channelDataMap.end()) {
               // We haven't. Initialize the bookkeeping information
               // for this channel.
               ChannelBookKeeping bookKeeping;
 
               // Add a new channel to the end of the list we'll
               // write out after we've processed this event.
               bookKeeping.channelIndex = channels->size();
               channels->emplace_back(channel);
               channelIndex = bookKeeping.channelIndex;
 
               // Initialize a vector with the index of the step that
               // created this channel.
               bookKeeping.stepList.push_back(edIndex);
 
               // Save the bookkeeping information for this channel.
               channelDataMap[channel] = bookKeeping;
             }
             else {
               // We've created this SimChannel for a previous energy
               // deposit. Get its address.
 
               auto& bookKeeping = search->second;
               channelIndex = bookKeeping.channelIndex;
 
               // Has this step contributed to this channel before?
               auto& stepList = bookKeeping.stepList;
               auto stepSearch = std::find(stepList.begin(), stepList.end(), edIndex);
               if (stepSearch == stepList.end()) {
                 // No, so add this step's index to the list.
                 stepList.push_back(edIndex);
               }
             }
 
             sim::SimChannel* channelPtr = &(channels->at(channelIndex));
 
             // Add the electron clusters and energy to the
             // sim::SimChannel
             channelPtr->AddIonizationElectrons(
               energyDeposit.TrackID(), tdc, fnElDiff[k], xyz, fnEnDiff[k]);
 
             if (fStoreDriftedElectronClusters)
               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());
           }
           catch (cet::exception& e) {
             mf::LogDebug("SimDriftElectrons")
               << "unable to drift electrons from point (" << xyz[0] << "," << xyz[1] << ","
               << xyz[2] << ") with exception " << e;
           } // end try to determine channel
         }   // end loop over clusters
       }     // end loop over planes
     }       // for each sim::SimEnergyDeposit
 
     // Write the sim::SimChannel collection.
     event.put(std::move(channels));
     if (fStoreDriftedElectronClusters) event.put(std::move(SimDriftedElectronClusterCollection));
   }