Transports energy depositions from GEANT4 to TPC channels. More...

#include "LArVoxelReadout.h"

Inheritance diagram for larg4::LArVoxelReadout:

Classes
struct	Setup_t
	Collection of what it takes to set a `LArVoxelReadout` up. More...

Public Types
typedef std::map< unsigned int, sim::SimChannel >	ChannelMap_t
	Type of map channel -> sim::SimChannel. More...

Public Member Functions
	LArVoxelReadout (std::string const &name)
	Constructor. Can detect which TPC to cover by the name. More...

	LArVoxelReadout (std::string const &name, unsigned int cryostat, unsigned int tpc)
	Constructor. Sets which TPC to work on. More...

virtual	~LArVoxelReadout ()

void	Setup (Setup_t const &setupData)
	Reads all the configuration elements from `setupData` More...

void	SetSingleTPC (unsigned int cryostat, unsigned int tpc)
	Associates this readout to one specific TPC. More...

void	SetDiscoverTPC ()
	Sets this readout to discover the TPC of each processed hit. More...

virtual void	Initialize (G4HCofThisEvent *)

virtual void	EndOfEvent (G4HCofThisEvent *)

virtual void	clear ()

virtual G4bool	ProcessHits (G4Step , G4TouchableHistory )

virtual void	DrawAll ()

virtual void	PrintAll ()

void	ClearSimChannels ()

std::vector< sim::SimChannel >	GetSimChannels () const
	Creates a list with the accumulated information for the single TPC. More...

std::vector< sim::SimChannel >	GetSimChannels (unsigned short cryo, unsigned short tpc) const
	Creates a list with the accumulated information for specified TPC. More...


const ChannelMap_t &	GetSimChannelMap () const
	Returns the accumulated channel -> SimChannel map for the single TPC. More...

ChannelMap_t &	GetSimChannelMap ()
	Returns the accumulated channel -> SimChannel map for the single TPC. More...


const ChannelMap_t &	GetSimChannelMap (unsigned short cryo, unsigned short tpc) const
	Returns the accumulated channel -> SimChannel map for the specified TPC. More...

ChannelMap_t &	GetSimChannelMap (unsigned short cryo, unsigned short tpc)
	Returns the accumulated channel -> SimChannel map for the specified TPC. More...

Private Member Functions
void	SetOffPlaneChargeRecoveryMargin (double margin)
	Sets the margin for recovery of charge drifted off-plane. More...

void	SetRandomEngines (CLHEP::HepRandomEngine *pPropGen)
	Sets the random generators to be used. More...

geo::Point_t	RecoverOffPlaneDeposit (geo::Point_t const &pos, geo::PlaneGeo const &plane) const
	Returns the point on the specified plane closest to position. More...

void	DriftIonizationElectrons (G4ThreeVector stepMidPoint, const double simTime, int trackID, unsigned short int cryostat, unsigned short int tpc)

bool	Has (std::vector< unsigned short int > v, unsigned short int tpc) const

void	ProcessStep (G4Step *)

Private Attributes
double	fDriftVelocity [3]

double	fLongitudinalDiffusion

double	fTransverseDiffusion

double	fElectronLifetime

double	fElectronClusterSize

int	fMinNumberOfElCluster

bool	fDontDriftThem

std::vector< unsigned short int >	fSkipWireSignalInTPCs

double	fOffPlaneMargin = 0.0
	Charge deposited within this many [cm] from the plane is lead onto it. More...

std::vector< std::vector< ChannelMap_t > >	fChannelMaps
	Maps of cryostat, tpc to channel data. More...

art::ServiceHandle< geo::Geometry >	fGeoHandle
	Handle to the Geometry service. More...

art::ServiceHandle< sim::LArG4Parameters >	fLgpHandle
	Handle to the LArG4 parameters service. More...

unsigned int	fTPC
	which TPC this LArVoxelReadout corresponds to More...

unsigned int	fCstat
	and in which cryostat (if bSingleTPC is true) More...

bool	bSingleTPC
	true if this readout is associated with a single TPC More...

CLHEP::HepRandomEngine *	fPropGen = nullptr
	random engine for charge propagation More...

::detinfo::ElecClock	fClock
	TPC electronics clock. More...

G4ThreeVector	fStepStart

G4ThreeVector	fStepEnd

size_t	fNSteps

Detailed Description

Transports energy depositions from GEANT4 to TPC channels.

This class acts on single energy depositions from GEANT4, simulating the transportation of the ensuing ionisation electrons to the readout channels:

the number of ionisation electrons is read from the current larg4::IonizationAndScintillation instance
space charge displacement is optionally applied
lifetime correction is applied
charge is split in small electron clusters
each cluster is subject to longitudinal and transverse diffusion
each cluster is assigned to one TPC channel for each wire plane
optionally, charge is forced to stay on the planes; otherwise charge drifting outside the plane is lost

For each energy deposition, entries on the appropriate sim::SimChannel are added, with the information of the position where the energy deposit happened (in global coordinates, centimeters), the ID of the Geant4 track which produced the deposition, and the quantized time of arrival to the channel (in global TDC tick units). At most one entry is added for each electron cluster, but entries from the same energy deposit can be compacted if falling on the same TDC tick.

The main entry point of this class is the method ProcessHits().

Options

A few optional behaviours are supported:

lead off-plane charge to the planes: regulated by RecoverOffPlaneDeposit(), if charge which reaches a wire plane is actually off it by less than the chosen margin, it's accounted for by that plane; by default the margin is 0 and all the charge off the plane is lost (with a warning)

Definition at line 157 of file LArVoxelReadout.h.

Member Typedef Documentation

typedef std::map<unsigned int, sim::SimChannel> larg4::LArVoxelReadout::ChannelMap_t

Type of map channel -> sim::SimChannel.

Definition at line 161 of file LArVoxelReadout.h.

Constructor & Destructor Documentation

larg4::LArVoxelReadout::LArVoxelReadout ( std::string const & name )

Constructor. Can detect which TPC to cover by the name.

Definition at line 47 of file LArVoxelReadout.cxx.

References ClearSimChannels(), fClock, SetDiscoverTPC(), SetSingleTPC(), and detinfo::DetectorClocks::TPCClock().

     : G4VSensitiveDetector(name)
   {
     // Initialize values for the electron-cluster calculation.
     ClearSimChannels();
 
     const detinfo::DetectorClocks* ts = lar::providerFrom<detinfo::DetectorClocksService>();
     fClock = ts->TPCClock();
 
     // the standard name contains cryostat and TPC;
     // if we don't find it, we will detect the TPC at each Geant hit
     unsigned int cryostat, tpc;
     if (std::sscanf(name.c_str(),"%*19s%1u%*4s%u",&cryostat, &tpc) == 2)
       SetSingleTPC(cryostat, tpc);
     else
       SetDiscoverTPC();
     
   } // LArVoxelReadout::LArVoxelReadout()

larg4::LArVoxelReadout::LArVoxelReadout	(	std::string const &	name,
		unsigned int	cryostat,
		unsigned int	tpc
	)

Constructor. Sets which TPC to work on.

Definition at line 70 of file LArVoxelReadout.cxx.

References SetSingleTPC().

71 : LArVoxelReadout(name)

72 { SetSingleTPC(cryostat, tpc); }

larg4::LArVoxelReadout::LArVoxelReadout

LArVoxelReadout(std::string const &name)

Constructor. Can detect which TPC to cover by the name.

Definition: LArVoxelReadout.cxx:47

larg4::LArVoxelReadout::SetSingleTPC

void SetSingleTPC(unsigned int cryostat, unsigned int tpc)

Associates this readout to one specific TPC.

Definition: LArVoxelReadout.cxx:85

larg4::LArVoxelReadout::~LArVoxelReadout ( )

virtual

Definition at line 102 of file LArVoxelReadout.cxx.

102 {}

Member Function Documentation

void larg4::LArVoxelReadout::clear ( void )

virtual

Definition at line 142 of file LArVoxelReadout.cxx.

143 {

144 }

void larg4::LArVoxelReadout::ClearSimChannels ( )

Definition at line 147 of file LArVoxelReadout.cxx.

References fChannelMaps, fGeoHandle, geo::GeometryCore::Ncryostats(), and geo::GeometryCore::NTPC().

Referenced by LArVoxelReadout().

                                          {
     fChannelMaps.resize(fGeoHandle->Ncryostats());
     size_t cryo = 0;
     for (auto& cryoData: fChannelMaps) { // each, a vector of maps
       cryoData.resize(fGeoHandle->NTPC(cryo++));
       for (auto& channelsMap: cryoData) channelsMap.clear(); // each, a map
     } // for cryostats
   } // LArVoxelReadout::ClearSimChannels()

void larg4::LArVoxelReadout::DrawAll ( )

virtual

Definition at line 582 of file LArVoxelReadout.cxx.

582 {}

void larg4::LArVoxelReadout::DriftIonizationElectrons	(	G4ThreeVector	stepMidPoint,
		const double	simTime,
		int	trackID,
		unsigned short int	cryostat,
		unsigned short int	tpc
	)

private

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

Definition at line 331 of file LArVoxelReadout.cxx.

References sim::SimChannel::AddIonizationElectrons(), geo::TPCGeo::DriftDirection(), e, energy, larg4::IonizationAndScintillation::EnergyDeposit(), fChannelMaps, fClock, fDriftVelocity, fElectronClusterSize, fElectronLifetime, fGeoHandle, fLongitudinalDiffusion, fMinNumberOfElCluster, fOffPlaneMargin, fPropGen, fTransverseDiffusion, larg4::IonizationAndScintillation::Instance(), geo::kNegX, geo::kPosX, LOG_DEBUG, geo::GeometryCore::NearestChannel(), geo::TPCGeo::Nplanes(), larg4::IonizationAndScintillation::NumberIonizationElectrons(), geo::TPCGeo::Plane(), geo::TPCGeo::Plane0Pitch(), geo::TPCGeo::PlaneLocation(), geo::TPCGeo::PlanePitch(), RecoverOffPlaneDeposit(), detinfo::ElecClock::Ticks(), geo::GeometryCore::TPC(), and xx.

Referenced by ProcessHits().

   {
     auto const * ts = lar::providerFrom<detinfo::DetectorClocksService>();
 
     // this must be always true, unless caller has been sloppy
     assert(fPropGen); // No propagation random generator provided?!
     
     CLHEP::RandGauss PropRand(*fPropGen);
 
     // This routine gets called frequently, once per every particle
     // traveling through every voxel. Use whatever tricks we can to
     // increase its execution speed.
 
     static double LifetimeCorr_const = -1000. * fElectronLifetime;
     static double LDiff_const        = std::sqrt(2.*fLongitudinalDiffusion);
     static double TDiff_const        = std::sqrt(2.*fTransverseDiffusion);
     static double RecipDriftVel[3]   = {1./fDriftVelocity[0], 
                                         1./fDriftVelocity[1], 
                                         1./fDriftVelocity[2]};
 
     struct Deposit_t {
       double energy = 0.;
       double electrons = 0.;
       
       void add(double more_energy, double more_electrons)
         { energy += more_energy; electrons += more_electrons; }
     }; // Deposit_t
     
     // Map of electrons to store - catalogued by map[channel][tdc]
     std::map<raw::ChannelID_t, std::map<unsigned int, Deposit_t>> DepositsToStore;
 
     double xyz1[3] = {0.};
 
     double const xyz[3] = {stepMidPoint.x() / CLHEP::cm,
                            stepMidPoint.y() / CLHEP::cm,
                            stepMidPoint.z() / CLHEP::cm};
 
     // Already know which TPC we're in because we have been told
 
     try{
       const geo::TPCGeo &tpcg = fGeoHandle->TPC(tpc, cryostat);
 
       // X drift distance - the drift direction can be either in
       // the positive or negative direction, so use std::abs
 
       double XDrift = std::abs(stepMidPoint.x()/CLHEP::cm - tpcg.PlaneLocation(0)[0]);
       //std::cout<<tpcg.DriftDirection()<<std::endl;
       if (tpcg.DriftDirection() == geo::kNegX)
         XDrift = stepMidPoint.x()/CLHEP::cm - tpcg.PlaneLocation(0)[0];
       else if (tpcg.DriftDirection() == geo::kPosX)
         XDrift = tpcg.PlaneLocation(0)[0] - stepMidPoint.x()/CLHEP::cm;
       
       if(XDrift < 0.) return;
 
       // Get SCE {x,y,z} offsets for particular location in TPC      
       geo::Vector_t posOffsets;
       auto const* SCE = lar::providerFrom<spacecharge::SpaceChargeService>();
       if (SCE->EnableSimSpatialSCE() == true)
       {
         posOffsets = SCE->GetPosOffsets({ xyz[0], xyz[1], xyz[2] });
       }
       posOffsets.SetX(-posOffsets.X());
 
       // Drift time (nano-sec)
       double TDrift;
       XDrift += posOffsets.X();
       
       // 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 (XDrift < 0.) XDrift = 0.;
       
       TDrift = XDrift * RecipDriftVel[0];
       if (tpcg.Nplanes() == 2){// special case for ArgoNeuT (plane 0 is the second wire plane)
         TDrift = ((XDrift - tpcg.PlanePitch(0,1)) * RecipDriftVel[0] 
                   + tpcg.PlanePitch(0,1) * RecipDriftVel[1]);
       }
           
       const double lifetimecorrection = TMath::Exp(TDrift / LifetimeCorr_const);
       const int    nIonizedElectrons  = larg4::IonizationAndScintillation::Instance()->NumberIonizationElectrons();
       const double energy             = larg4::IonizationAndScintillation::Instance()->EnergyDeposit();
       
       // if we have no electrons (too small energy or too large recombination)
       // we are done already here
       if (nIonizedElectrons <= 0) {
         LOG_DEBUG("LArVoxelReadout")
           << "No electrons drifted to readout, " << energy << " MeV lost.";
         return;
       }
       // includes the effect of lifetime
       const double nElectrons   = nIonizedElectrons * lifetimecorrection;
 
       // Longitudinal & transverse diffusion sigma (cm)
       double SqrtT    = std::sqrt(TDrift);
       double LDiffSig = SqrtT * LDiff_const;
       double TDiffSig = SqrtT * TDiff_const;
       double electronclsize = fElectronClusterSize;
       
       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);
       }
       
       // Compute arrays of values as quickly as possible.
       std::vector< double > XDiff(nClus);
       std::vector< double > YDiff(nClus);
       std::vector< double > ZDiff(nClus);
       std::vector< double > nElDiff(nClus, electronclsize);
       std::vector< double > nEnDiff(nClus);
 
       // fix the number of electrons in the last cluster, that has smaller size
       nElDiff.back() = nElectrons - (nClus-1)*electronclsize;
       
       for(size_t xx = 0; xx < nElDiff.size(); ++xx){
         if(nElectrons > 0) nEnDiff[xx] = energy/nElectrons*nElDiff[xx];
         else               nEnDiff[xx] = 0.;
       }
       
       double const avegageYtransversePos
         = (stepMidPoint.y()/CLHEP::cm) + posOffsets.Y();
       double const avegageZtransversePos
         = (stepMidPoint.z()/CLHEP::cm) + posOffsets.Z();
       
       // Smear drift times by x position and drift time
       if (LDiffSig > 0.0)
         PropRand.fireArray( nClus, &XDiff[0], 0., LDiffSig);
       else
         XDiff.assign(nClus, 0.0);
       
       if (TDiffSig > 0.0) {
         // Smear the Y,Z position by the transverse diffusion
         PropRand.fireArray( nClus, &YDiff[0], avegageYtransversePos, TDiffSig);
         PropRand.fireArray( nClus, &ZDiff[0], avegageZtransversePos, TDiffSig);
       }
       else {
         YDiff.assign(nClus, avegageYtransversePos);
         ZDiff.assign(nClus, avegageZtransversePos);
       }
 
       // make a collection of electrons for each plane
       for(size_t p = 0; p < tpcg.Nplanes(); ++p){
         
         geo::PlaneGeo const& plane = tpcg.Plane(p);
 
         double Plane0Pitch = tpcg.Plane0Pitch(p);
         
         // "-" sign is because Plane0Pitch output is positive. Andrzej
         xyz1[0] = tpcg.PlaneLocation(0)[0] - Plane0Pitch;
         
         // 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 + XDiff[k] * RecipDriftVel[0];
           // Take into account different Efields between planes
           // Also take into account special case for ArgoNeuT where Nplanes = 2.
           for (size_t ip = 0; ip<p; ++ip){
             TDiff += tpcg.PlanePitch(ip,ip+1) * RecipDriftVel[tpcg.Nplanes()==3?ip+1:ip+2];
           }
           xyz1[1] = YDiff[k];
           xyz1[2] = ZDiff[k];
           
           
           // grab the nearest channel to the xyz1 position
           try{
             if (fOffPlaneMargin != 0) {
               // get the effective position where to consider the charge landed;
               // 
               // Some optimisations are possible; in particular, this method
               // could be extended to inform us if the point was too far.
               // Currently, if that is the case the code will proceed, find the
               // point is off plane, emit a warning and skip the deposition.
               //
               auto const landingPos
                 = RecoverOffPlaneDeposit({ xyz1[0], xyz1[1], xyz1[2] }, plane);
               xyz1[0] = landingPos.X();
               xyz1[1] = landingPos.Y();
               xyz1[2] = landingPos.Z();
               
             } // if charge lands off plane
             uint32_t channel = fGeoHandle->NearestChannel(xyz1, p, tpc, cryostat);
             
             // Add potential decay/capture/etc delay effect, simTime.
             unsigned int tdc = fClock.Ticks(ts->G4ToElecTime(TDiff + simTime));
 
             // Add electrons produced by each cluster to the map
             DepositsToStore[channel][tdc].add(nEnDiff[k], nElDiff[k]);
           }
           catch(cet::exception &e){
             mf::LogWarning("LArVoxelReadout") << "unable to drift electrons from point ("
                                               << xyz[0] << "," << xyz[1] << "," << xyz[2]
                                               << ") with exception " << e;
           }
         } // end loop over clusters
       } // end loop over planes
       
       // Now store them in SimChannels
       ChannelMap_t& ChannelDataMap = fChannelMaps[cryostat][tpc];
       
       // browse deposited data on each channel: (channel; deposit data in time)
       for(auto const& deposit_per_channel: DepositsToStore){
         
         raw::ChannelID_t channel = deposit_per_channel.first;
         
         // find whether we already have this channel
         auto iChannelData = ChannelDataMap.find(channel);
         
         // channelData is the SimChannel these deposits are going to be added to
         // If there is such a channel already, use it (first beanch).
         // If it's a new channel, the inner assignment creates a new SimChannel
         // in the map, and we save its reference in channelData
         sim::SimChannel& channelData
           = (iChannelData == ChannelDataMap.end())
           ? (ChannelDataMap[channel] = sim::SimChannel(channel))
           : iChannelData->second;
         
         // go through all deposits, one for each TDC: (TDC, deposit data)
         for(auto const& deposit_per_tdc: deposit_per_channel.second) {
           channelData.AddIonizationElectrons(trackID,
                                              deposit_per_tdc.first,
                                              deposit_per_tdc.second.electrons,
                                              xyz,
                                              deposit_per_tdc.second.energy);
           
         } // for deposit on TDCs
       } // for deposit on channels
 
     } // end try intended to catch points where TPC can't be found
     catch(cet::exception &e){
       mf::LogWarning("LArVoxelReadout") << "step cannot be found in a TPC\n"
                                         << e;
     }
 
     return;
   }

void larg4::LArVoxelReadout::EndOfEvent ( G4HCofThisEvent * )

virtual

Definition at line 135 of file LArVoxelReadout.cxx.

References fNSteps.

   {
     mf::LogWarning("LArVoxelReadout") << "Total number of steps was " << fNSteps << std::endl;
 
   } // LArVoxelReadout::EndOfEvent()

const LArVoxelReadout::ChannelMap_t & larg4::LArVoxelReadout::GetSimChannelMap ( ) const

Returns the accumulated channel -> SimChannel map for the single TPC.

Definition at line 157 of file LArVoxelReadout.cxx.

References bSingleTPC, fCstat, and fTPC.

Referenced by GetSimChannelMap(), and larg4::LArG4::produce().

   {
     if (bSingleTPC) return GetSimChannelMap(fCstat, fTPC);
     throw cet::exception("LArVoxelReadout") << "TPC not specified";
   } // LArVoxelReadout::GetSimChannelMap() const

LArVoxelReadout::ChannelMap_t & larg4::LArVoxelReadout::GetSimChannelMap ( )

Returns the accumulated channel -> SimChannel map for the single TPC.

Definition at line 163 of file LArVoxelReadout.cxx.

References bSingleTPC, fCstat, fTPC, and GetSimChannelMap().

                                                                  {
     if (bSingleTPC) return GetSimChannelMap(fCstat, fTPC);
     throw cet::exception("LArVoxelReadout") << "TPC not specified";
   } // LArVoxelReadout::GetSimChannelMap()

const LArVoxelReadout::ChannelMap_t & larg4::LArVoxelReadout::GetSimChannelMap	(	unsigned short	cryo,
		unsigned short	tpc
	)		const

Returns the accumulated channel -> SimChannel map for the specified TPC.

Definition at line 170 of file LArVoxelReadout.cxx.

References fChannelMaps, and GetSimChannelMap().

171 { return fChannelMaps.at(cryo).at(tpc); }

larg4::LArVoxelReadout::fChannelMaps

std::vector< std::vector< ChannelMap_t > > fChannelMaps

Maps of cryostat, tpc to channel data.

Definition: LArVoxelReadout.h:327

LArVoxelReadout::ChannelMap_t & larg4::LArVoxelReadout::GetSimChannelMap	(	unsigned short	cryo,
		unsigned short	tpc
	)

Returns the accumulated channel -> SimChannel map for the specified TPC.

Definition at line 174 of file LArVoxelReadout.cxx.

References fChannelMaps.

175 { return fChannelMaps.at(cryo).at(tpc); }

larg4::LArVoxelReadout::fChannelMaps

std::vector< std::vector< ChannelMap_t > > fChannelMaps

Maps of cryostat, tpc to channel data.

Definition: LArVoxelReadout.h:327

std::vector< sim::SimChannel > larg4::LArVoxelReadout::GetSimChannels ( ) const

Creates a list with the accumulated information for the single TPC.

Definition at line 178 of file LArVoxelReadout.cxx.

References bSingleTPC, fCstat, and fTPC.

                                                                  {
     if (bSingleTPC) return GetSimChannels(fCstat, fTPC);
     throw cet::exception("LArVoxelReadout") << "TPC not specified";
   } // LArVoxelReadout::GetSimChannels()

std::vector< sim::SimChannel > larg4::LArVoxelReadout::GetSimChannels	(	unsigned short	cryo,
		unsigned short	tpc
	)		const

Creates a list with the accumulated information for specified TPC.

Definition at line 184 of file LArVoxelReadout.cxx.

References fChannelMaps.

   { 
     std::vector<sim::SimChannel> channels;
     const ChannelMap_t& chmap = fChannelMaps.at(cryo).at(tpc);
     channels.reserve(chmap.size());
     for(const auto& chpair: chmap) channels.push_back(chpair.second);
     return channels;
   } // LArVoxelReadout::GetSimChannels(short, short)

bool larg4::LArVoxelReadout::Has	(	std::vector< unsigned short int >	v,
		unsigned short int	tpc
	)		const

inlineprivate

Definition at line 304 of file LArVoxelReadout.h.

Referenced by ProcessHits().

     {   
         for (auto c: v) if (c == tpc) return true;      
         return false;
     }

void larg4::LArVoxelReadout::Initialize ( G4HCofThisEvent * )

virtual

Definition at line 106 of file LArVoxelReadout.cxx.

References sim::LArG4Parameters::DisableWireplanes(), sim::LArG4Parameters::ElectronClusterSize(), fDontDriftThem, fDriftVelocity, fElectronClusterSize, fElectronLifetime, fLgpHandle, fLongitudinalDiffusion, fMinNumberOfElCluster, fNSteps, fSkipWireSignalInTPCs, fTransverseDiffusion, LOG_DEBUG, sim::LArG4Parameters::LongitudinalDiffusion(), sim::LArG4Parameters::MinNumberOfElCluster(), sim::LArG4Parameters::NoElectronPropagation(), sim::LArG4Parameters::SkipWireSignalInTPCs(), and sim::LArG4Parameters::TransverseDiffusion().

   {
     // for c2: larp is unused
     //auto const * larp = lar::providerFrom<detinfo::LArPropertiesService>();
     auto const * detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();
     fElectronLifetime      = detprop->ElectronLifetime();
     for (int i = 0; i<3; ++i)
       fDriftVelocity[i]    = detprop->DriftVelocity(detprop->Efield(i),
                                                     detprop->Temperature())/1000.;
     
     fElectronClusterSize   = fLgpHandle->ElectronClusterSize();
     fMinNumberOfElCluster  = fLgpHandle->MinNumberOfElCluster();
     fLongitudinalDiffusion = fLgpHandle->LongitudinalDiffusion();
     fTransverseDiffusion   = fLgpHandle->TransverseDiffusion();
     fDontDriftThem         = fLgpHandle->DisableWireplanes();
     fSkipWireSignalInTPCs  = fLgpHandle->SkipWireSignalInTPCs();
 
     LOG_DEBUG("LArVoxelReadout")  << " e lifetime: "        << fElectronLifetime
                                   << "\n Temperature: "     << detprop->Temperature()
                                   << "\n Drift velocity: "  << fDriftVelocity[0]
                                   <<" "<<fDriftVelocity[1]<<" "<<fDriftVelocity[2];
 
     fDontDriftThem = (fDontDriftThem || fLgpHandle->NoElectronPropagation() );
 
     fNSteps=0;
   }

void larg4::LArVoxelReadout::PrintAll ( )

virtual

Definition at line 583 of file LArVoxelReadout.cxx.

583 {}

G4bool larg4::LArVoxelReadout::ProcessHits	(	G4Step *	step,
		G4TouchableHistory *	pHistory
	)

virtual

Definition at line 196 of file LArVoxelReadout.cxx.

References bSingleTPC, DriftIonizationElectrons(), fCstat, fDontDriftThem, fNSteps, fSkipWireSignalInTPCs, fTPC, larg4::ParticleListAction::GetCurrentTrackID(), Has(), larg4::G4PVPlacementWithID< IDTYPE >::ID, larg4::IonizationAndScintillation::Instance(), LOG_DEBUG, and larg4::IonizationAndScintillation::Reset().

   {
     // All work done for the "parallel world" "box of voxels" in
     // LArVoxelReadoutGeometry makes this a fairly simple routine.
     // First, the usual check for non-zero energy:
 
     // Only process the hit if the step is inside the active volume and
     // it has deposited energy.  The hit being inside the active volume
     // is virtually sure to happen because the LArVoxelReadoutGeometry
     // that this class makes use of only has voxels for inside the TPC.
 
     // The step can be no bigger than the size of the voxel,
     // because of the geometry set up in LArVoxelGeometry and the
     // transportation set up in PhysicsList.  Find the mid-point
     // of the step.
 
     if ( step->GetTotalEnergyDeposit() > 0 ){
       
       // Make sure we have the IonizationAndScintillation singleton
       // reset to this step
       larg4::IonizationAndScintillation::Instance()->Reset(step);
       fNSteps++;
       if( !fDontDriftThem ){
 
         G4ThreeVector midPoint = 0.5*( step->GetPreStepPoint()->GetPosition()
                                        + step->GetPostStepPoint()->GetPosition() );
         double g4time = step->GetPreStepPoint()->GetGlobalTime();
         
         // Find the Geant4 track ID for the particle responsible for depositing the
         // energy.  if we are only storing primary EM shower particles, and this energy
         // is from a secondary etc EM shower particle, the ID returned is the primary
         const int trackID = ParticleListAction::GetCurrentTrackID();
         
         // Find out which TPC we are in.
         // If this readout object covers just one, we already know it.
         // Otherwise, we have to ask Geant where we are.
         unsigned short int cryostat = 0, tpc = 0;
         if (bSingleTPC) {
           cryostat = fCstat;
           tpc = fTPC;
         }
         else {
           // detect the TPC we are in
           const G4VTouchable* pTouchable = step->GetPreStepPoint()->GetTouchable();
           if (!pTouchable) {
             throw cet::exception
               ("LArG4") << "Untouchable step in LArVoxelReadout::ProcessHits()";
           }
           
           // one of the ancestors of the touched volume is supposed to be
           // actually a G4PVPlacementInTPC that knows which TPC it covers;
           // currently, it's the 4th in the ladder:
           // [0] voxel [1] voxel tower [2] voxel plane [3] the full box;
           G4int depth = 0;
           while (depth < pTouchable->GetHistoryDepth()) {
             const G4PVPlacementInTPC* pPVinTPC = 
               dynamic_cast<const G4PVPlacementInTPC*>
               (pTouchable->GetVolume(depth++));
             if (!pPVinTPC) continue;
             cryostat = pPVinTPC->ID.Cryostat;
             tpc = pPVinTPC->ID.TPC;
             if (Has(fSkipWireSignalInTPCs, tpc))
             {
                 return true;
             }
             break;
           } // while
           if (depth < pTouchable->GetHistoryDepth()) {
             // this is a fundamental error where the step does not happen in
             // any TPC; this should not happen in the readout geometry!
             throw cet::exception
               ("LArG4") << "No TPC ID found in LArVoxelReadout::ProcessHits()";
           } // if
           LOG_DEBUG("LArVoxelReadoutHit") << " hit in C=" << cryostat << " T=" << tpc;
         } // if more than one TPC
         
         // Note that if there is no particle ID for this energy deposit, the
         // trackID will be sim::NoParticleId.
         
         DriftIonizationElectrons(midPoint, g4time, trackID, cryostat, tpc);
       } // end we are drifting
     } // end there is non-zero energy deposition
 
     return true;
   }

void larg4::LArVoxelReadout::ProcessStep ( G4Step * )

private

geo::Point_t larg4::LArVoxelReadout::RecoverOffPlaneDeposit	(	geo::Point_t const &	pos,
		geo::PlaneGeo const &	plane
	)		const

private

Returns the point on the specified plane closest to position.

Parameters

pos	the position to be tested (global coordinates, centimeters)
plane	the plane to test the position against

Returns: a position on plane, unless pos is too far from it

This method considers the distance of the position pos from the active part of the plane (see geo::Plane::DeltaFromActivePlane()). If the position is less than a configurable margin far from the plane, the closest point on the plane to that position is returned. Otherwise, the position itself is returned.

Ionization charge may be drifted so that when it arrives to the plane, it actually does not hit the area covered by wires. This can happen for many reasons:

space charge distortion led the point outside the fiducial volume (this may be prevented by specific code)
diffusion pushes the charge outside the instrumented region
the geometry of the wire planes is such that planes have different coverage and what one plane can cover, the next can't

The "recovery" consists in forcing the charge to the instrumented area covered by the plane wires. The distance of the drifted charge from each plane border is computed and compared to the margin. If that distance is smaller than the margin, it is neglected and the charge is assigned a new position on that border.

This method provides the position that should be used for the charge deposition.

This is a simplistic approach to the simulation of border effects, assuming that in fact the electric field, which is continuous and pointing to the collection wires, will drive the charge to the wires even when they are "off track". No correction is applied for the additional time that such deviation would take.

Definition at line 292 of file LArVoxelReadout.cxx.

References geo::PlaneGeo::ComposePoint(), geo::PlaneGeo::DeltaFromActivePlane(), geo::PlaneGeo::DistanceFromPlane(), fOffPlaneMargin, and geo::PlaneGeo::PointWidthDepthProjection().

Referenced by DriftIonizationElectrons(), and SetRandomEngines().

   {
     //
     // translate the landing position on the two frame coordinates
     // ("width" and "depth")
     //
     auto const landingPos = plane.PointWidthDepthProjection(pos);
     
     //
     // compute the distance of the landing position on the two frame
     // coordinates ("width" and "depth");
     // keep the point within 10 micrometers (0.001 cm) from the border
     //
     auto const offPlane = plane.DeltaFromActivePlane(landingPos, 0.001);
     
     //
     // if both the distances are below the margin, move the point to
     // the border
     //
     
     // nothing to recover: landing is inside
     if ((offPlane.X() == 0.0) && (offPlane.Y() == 0.0)) return pos;
     
     // won't recover: too far
     if ((std::abs(offPlane.X()) > fOffPlaneMargin)
       || (std::abs(offPlane.Y()) > fOffPlaneMargin))
       return pos;
     
     // we didn't fully decompose because it might be unnecessary;
     // now we need the full thing
     auto const distance = plane.DistanceFromPlane(pos);
     
     return plane.ComposePoint<geo::Point_t>(distance, landingPos + offPlane);
     
   } // LArVoxelReadout::RecoverOffPlaneDeposit()

void larg4::LArVoxelReadout::SetDiscoverTPC ( )

Sets this readout to discover the TPC of each processed hit.

Definition at line 93 of file LArVoxelReadout.cxx.

References bSingleTPC, fCstat, fTPC, and LOG_DEBUG.

Referenced by LArVoxelReadout().

                                        {
     bSingleTPC = false;
     fCstat = 0;
     fTPC = 0;
     LOG_DEBUG("LArVoxelReadout") << GetName() << " autodetects TPC";
   } // LArVoxelReadout::SetDiscoverTPC()

void larg4::LArVoxelReadout::SetOffPlaneChargeRecoveryMargin ( double margin )

inlineprivate

Sets the margin for recovery of charge drifted off-plane.

Parameters

margin the extent of the margin on each frame coordinate [cm]

This method sets the margin for the recovery of off-plane ionization charge. See RecoverOffPlaneDeposit() for a description of that feature.

This method is used by LArVoxelReadout::Setup().

Definition at line 251 of file LArVoxelReadout.h.

References max.

Referenced by Setup().

252 { fOffPlaneMargin = std::max(margin, 0.0); }

max

Int_t max

Definition: plot.C:27

larg4::LArVoxelReadout::fOffPlaneMargin

double fOffPlaneMargin

Charge deposited within this many [cm] from the plane is lead onto it.

Definition: LArVoxelReadout.h:325

void larg4::LArVoxelReadout::SetRandomEngines ( CLHEP::HepRandomEngine * pPropGen )

private

Sets the random generators to be used.

Definition at line 284 of file LArVoxelReadout.cxx.

References fPropGen, and RecoverOffPlaneDeposit().

Referenced by Setup().

                                                                        {
     assert(pPropGen); // random engine must be present
     fPropGen = pPropGen;
   } // LArVoxelReadout::SetRandomEngines()

void larg4::LArVoxelReadout::SetSingleTPC	(	unsigned int	cryostat,
		unsigned int	tpc
	)

Associates this readout to one specific TPC.

Definition at line 85 of file LArVoxelReadout.cxx.

References bSingleTPC, fCstat, fTPC, and LOG_DEBUG.

Referenced by larg4::LArVoxelReadoutGeometry::Construct(), and LArVoxelReadout().

                                                                             {
     bSingleTPC = true;
     fCstat = cryostat;
     fTPC = tpc;
     LOG_DEBUG("LArVoxelReadout")
       << GetName() << "covers C=" << fCstat << " T=" << fTPC;
   } // LArVoxelReadout::SetSingleTPC()

void larg4::LArVoxelReadout::Setup ( Setup_t const & setupData )

Reads all the configuration elements from setupData

Definition at line 76 of file LArVoxelReadout.cxx.

References larg4::LArVoxelReadout::Setup_t::offPlaneMargin, larg4::LArVoxelReadout::Setup_t::propGen, SetOffPlaneChargeRecoveryMargin(), and SetRandomEngines().

Referenced by larg4::LArVoxelReadoutGeometry::Construct().

                                                       {
     
     SetOffPlaneChargeRecoveryMargin(setupData.offPlaneMargin);
     SetRandomEngines(setupData.propGen);
     
   } // LArVoxelReadout::Setup()