LArVoxelReadout.h

Go to the documentation of this file.

#ifndef LArG4_LArVoxelReadout_h
#define LArG4_LArVoxelReadout_h

#include <stdint.h>
#include <vector>
#include <algorithm> // std::max()

#include "Geant4/G4VSensitiveDetector.hh"
#include "Geant4/G4PVPlacement.hh"
#include "Geant4/globals.hh"

#include "lardataobj/Simulation/SimChannel.h"
#include "larcore/Geometry/Geometry.h"
#include "larcorealg/Geometry/TPCGeo.h"
#include "larsim/Simulation/LArG4Parameters.h"
#include "larsim/LArG4/IonizationAndScintillation.h"
#include "lardata/DetectorInfoServices/LArPropertiesService.h"
#include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
#include "lardata/DetectorInfoServices/DetectorClocksService.h"

// Forward declarations
class G4HCofThisEvent;
class G4TouchableHistory;
class G4Step;
namespace CLHEP { class HEPRandomEngine; }

namespace larg4 {

  struct TPCID_t {
    unsigned short int Cryostat, TPC;
    bool operator< (const TPCID_t& than) const
      {
        return (Cryostat < than.Cryostat)
          || ((Cryostat == than.Cryostat) && (TPC < than.TPC));
      } // operator< ()
  }; // TPCID_t
  
  template <class IDTYPE>
  class G4PVPlacementWithID: public G4PVPlacement {
      public:
    typedef IDTYPE ID_t;
    
    ID_t ID; 
    
    G4PVPlacementWithID(const G4Transform3D& Transform3D, const G4String &pName,
      G4LogicalVolume* pLogical, G4VPhysicalVolume* pMother,
      G4bool pMany, G4int pCopyNo, G4bool pSurfChk = false,
      ID_t id = ID_t()
      ):
      G4PVPlacement
        (Transform3D, pName, pLogical, pMother, pMany, pCopyNo, pSurfChk),
      ID(id)
      {}
    
    virtual ~G4PVPlacementWithID() {}
  }; // G4PVPlacementWithID<>
  
  typedef G4PVPlacementWithID<TPCID_t> G4PVPlacementInTPC;
  
  
  class LArVoxelReadout : public G4VSensitiveDetector
  {
  public:
    typedef std::map<unsigned int, sim::SimChannel> ChannelMap_t;
    
    struct Setup_t {
      
      CLHEP::HepRandomEngine* propGen = nullptr;
      
      double offPlaneMargin = 0.0;
    }; // struct Setup_t
    
    
    LArVoxelReadout(std::string const& name);
    
    LArVoxelReadout
      (std::string const& name, unsigned int cryostat, unsigned int tpc);

    // Destructor
    virtual ~LArVoxelReadout();
    
    void Setup(Setup_t const& setupData);
    
    void SetSingleTPC(unsigned int cryostat, unsigned int tpc);
    
    void SetDiscoverTPC();
    
    // Required for classes that inherit from G4VSensitiveDetector.
    //
    // Called at start and end of each event.
    virtual void Initialize(G4HCofThisEvent*);
    virtual void EndOfEvent(G4HCofThisEvent*);

    // Called to clear any accumulated information.
    virtual void clear();

    // The key method of this class.  It's called by Geant4 for each
    // step within the read-out geometry.  It accumulates the energy
    // in the G4Step in the LArVoxelList.
    virtual G4bool ProcessHits( G4Step*, G4TouchableHistory* );

    // Empty methods; they have to be defined, but they're rarely
    // used in Geant4 applications.
    virtual void DrawAll();
    virtual void PrintAll();

    // Independent method; clears the vector of SimChannels as well as the
    // channel number to SimChannel map.  Has to be separate from the 
    // clear method above because that run is run for every G4 event, ie
    // each MCTruth in the art::Event, while we want to only run this at 
    // the end of the G4 processing for each art::Event.
    void ClearSimChannels();

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

    const ChannelMap_t& GetSimChannelMap() const;
    ChannelMap_t& GetSimChannelMap();
    
    const ChannelMap_t& GetSimChannelMap
      (unsigned short cryo, unsigned short tpc) const;
    ChannelMap_t& GetSimChannelMap(unsigned short cryo, unsigned short tpc);


  private:

    void SetOffPlaneChargeRecoveryMargin(double margin)
      { fOffPlaneMargin = std::max(margin, 0.0); }

    void SetRandomEngines(CLHEP::HepRandomEngine* pPropGen);
    
    
    geo::Point_t RecoverOffPlaneDeposit
      (geo::Point_t const& pos, geo::PlaneGeo const& plane) const;
    
    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
    {   
        for (auto c: v) if (c == tpc) return true;      
        return false;
    }

    // Used in electron-cluster calculations
    // External parameters for the electron-cluster calculation.
    // obtained from LArG4Parameters, LArProperties, and DetectorProperties services
    double                                    fDriftVelocity[3];
    double                                    fLongitudinalDiffusion;
    double                                    fTransverseDiffusion;
    double                                    fElectronLifetime;
    double                                    fElectronClusterSize;
    int                                       fMinNumberOfElCluster;
    // for c2: unused private data members
    //double                                    fSampleRate;
    //int                                       fTriggerOffset;
    bool                                      fDontDriftThem;
    std::vector<unsigned short int>           fSkipWireSignalInTPCs;
    double                                    fOffPlaneMargin = 0.0;

    std::vector<std::vector<ChannelMap_t>>    fChannelMaps; 
    art::ServiceHandle<geo::Geometry>         fGeoHandle;  
    art::ServiceHandle<sim::LArG4Parameters>  fLgpHandle;  
    unsigned int                              fTPC;        
    unsigned int                              fCstat;      
    bool                                      bSingleTPC;  
    
    CLHEP::HepRandomEngine*                   fPropGen = nullptr;  
    
    ::detinfo::ElecClock                         fClock;      

    //these are the things for doing the separated EDeps
    void ProcessStep(G4Step*);
    
    G4ThreeVector                             fStepStart;
    G4ThreeVector                             fStepEnd;
    size_t fNSteps;
  };

}

#endif // LArG4_LArVoxelReadout_h