Class defining a sensitive detector of type HCal. More...

#include <HcalSD.h>

Public Member Functions
	HcalSD (const std::string &name, simcore::ConditionsInterface &ci, const framework::config::Parameters &params)
	Constructor.

virtual	~HcalSD ()=default
	Destructor.

bool	isSensDet (G4LogicalVolume *volume) const override
	Check if the input logical volume is a part of the hcal sensitive volumes.

ldmx::HcalID	decodeCopyNumber (const std::uint32_t copyNumber, const G4ThreeVector &localPosition, const G4Box *scint)
	Decode copy number of scintillator bar.

virtual G4bool	ProcessHits (G4Step aStep, G4TouchableHistory ROhist) override
	Create a hit out of the energy deposition deposited during a step.

virtual void	saveHits (framework::Event &event) override
	Add our hits to the event bus and then reset the container.

virtual void	onFinishedEvent () override
	Cleanup SD and prepare a new-event state.

Public Member Functions inherited from simcore::SensitiveDetector
	SensitiveDetector (const std::string &name, simcore::ConditionsInterface &ci, const framework::config::Parameters &parameters)
	Constructor.

	DECLARE_FACTORY_WITH_WAREHOUSE (SensitiveDetector, SensitiveDetector *, const std::string &, simcore::ConditionsInterface &, const framework::config::Parameters &)
	The SD Factory.

virtual	~SensitiveDetector ()=default
	Destructor.

virtual void	EndOfEvent (G4HCofThisEvent *) override
	This is Geant4's handle to tell us the event is ending.

Static Public Attributes
static const std::string	COLLECTION_NAME = "HcalSimHits"
	name of collection to be added to event bus

Private Attributes
std::vector< std::string >	gdml_identifiers_

double	birksc1_

double	birksc2_

std::map< ldmx::HcalID, ldmx::SimCalorimeterHit >	hits_
	map of hits to add to the event (will be squashed)

bool	enable_hit_contribs_
	enable hit contribs

bool	compress_hit_contribs_
	compress hit contribs

int	max_origin_track_id_
	maximum track ID to be considered an "origin"

Additional Inherited Members
Protected Member Functions inherited from simcore::SensitiveDetector
template<class T >
const T &	getCondition (const std::string &condition_name)
	Record the configuration of this detector into the run header.

bool	isGeantino (const G4Step *step) const
	Check if the passed step is a step of a geantino.

const TrackMap &	getTrackMap () const
	Get a handle to the current track map.

	enableLogging ("SensitiveDetector")
	Enable logging for this class.

Detailed Description

Class defining a sensitive detector of type HCal.

Definition at line 24 of file HcalSD.h.

Constructor & Destructor Documentation

◆ HcalSD()

simcore::HcalSD::HcalSD	(	const std::string &	name,
		simcore::ConditionsInterface &	ci,
		const framework::config::Parameters &	params )

Constructor.

Parameters

name	The name of the sensitive detector.
ci	Conditions interface handle
params	python configuration parameters

Definition at line 22 of file HcalSD.cxx.

    : SensitiveDetector(name, ci, p), birksc1_(1.29e-2), birksc2_(9.59e-6) {
  gdml_identifiers_ = {p.get<std::vector<std::string>>("gdml_identifiers")};
  enable_hit_contribs_ = p.get<bool>("enable_hit_contribs");
  compress_hit_contribs_ = p.get<bool>("compress_hit_contribs");
  max_origin_track_id_ = p.get<int>("max_origin_track_id");
}

References compress_hit_contribs_, enable_hit_contribs_, framework::config::Parameters::get(), and max_origin_track_id_.

Member Function Documentation

◆ decodeCopyNumber()

ldmx::HcalID simcore::HcalSD::decodeCopyNumber	(	const std::uint32_t	copyNumber,
		const G4ThreeVector &	localPosition,
		const G4Box *	scint )

Decode copy number of scintillator bar.

Parameters

copyNumber	The copy number of the scintillator volume.
localPosition	The position of the hit (step mid-point).
scint	The G4Box of the scintillator volume.

Definition at line 31 of file HcalSD.cxx.

                                                          {
  const unsigned int version{copyNumber / 0x01000000};
  if (version != 0) {
    using Index = ldmx::PackedIndex<256, 256, 256>;
    return ldmx::HcalID{Index(copyNumber).field2(), Index(copyNumber).field1(),
                        Index(copyNumber).field0()};
  }
  const auto& geometry = getCondition<ldmx::HcalGeometry>(
      ldmx::HcalGeometry::CONDITIONS_OBJECT_NAME);
  unsigned int strip_id = 0;
  const unsigned int section = copyNumber / 1000;
  const unsigned int layer = copyNumber % 1000;
 
  // 5cm wide bars are HARD-CODED
  if (section == ldmx::HcalID::BACK) {
    if (geometry.backLayerIsHorizontal(layer)) {
      strip_id = int((localPosition.y() + scint->GetYHalfLength()) / 50.0);
    } else {
      strip_id = int((localPosition.x() + scint->GetXHalfLength()) / 50.0);
    }
  } else {
    strip_id = int((localPosition.z() + scint->GetZHalfLength()) / 50.0);
  }
  return ldmx::HcalID{section, layer, strip_id};
}

References ldmx::HcalGeometry::CONDITIONS_OBJECT_NAME, and simcore::SensitiveDetector::getCondition().

Referenced by ProcessHits().

◆ isSensDet()

bool simcore::HcalSD::isSensDet ( G4LogicalVolume * volume ) const

inlineoverridevirtual

Check if the input logical volume is a part of the hcal sensitive volumes.

Note: This will match if a) the volume has the auxiliary tag "Region" set to contain to "CalorimeterRegion" and b) the volume name contains one of the identifiers in the gdml_identifiers parameter

Implements simcore::SensitiveDetector.

Definition at line 50 of file HcalSD.h.

                                                         {
    auto region = volume->GetRegion();
    if (region and region->GetName().contains("CalorimeterRegion")) {
      const auto name{volume->GetName()};
      return std::find_if(std::begin(gdml_identifiers_),
                          std::end(gdml_identifiers_),
                          [&name](const auto& identifier) {
                            return name.contains(identifier);
                          }) != std::end(gdml_identifiers_);
    }
    return false;
  }

◆ onFinishedEvent()

virtual void simcore::HcalSD::onFinishedEvent ( )

inlineoverridevirtual

Cleanup SD and prepare a new-event state.

Implements simcore::SensitiveDetector.

Definition at line 89 of file HcalSD.h.

89{ hits_.clear(); }

simcore::HcalSD::hits_

std::map< ldmx::HcalID, ldmx::SimCalorimeterHit > hits_

map of hits to add to the event (will be squashed)

Definition HcalSD.h:105

References hits_.

◆ ProcessHits()

G4bool simcore::HcalSD::ProcessHits	(	G4Step *	aStep,
		G4TouchableHistory *	ROhist )

overridevirtual

Create a hit out of the energy deposition deposited during a step.

Parameters

[in]	step	The current step.
[in]	history	The readout history.

Implements simcore::SensitiveDetector.

Definition at line 59 of file HcalSD.cxx.

                                                                    {
  // Get the edep from the step.
  G4double edep = aStep->GetTotalEnergyDeposit();
 
  // Skip steps with no energy dep which come from non-Geantino particles.
  if (edep == 0.0 and not isGeantino(aStep)) {
    ldmx_log(trace) << "CalorimeterSD skipping step with zero edep.";
    return false;
  }
 
  //---------------------------------------------------------------------------------------------------
  //                Birks' Law
  //                ===========
  //
  //                      In the case of Scintillator as active medium, we can
  //              describe the quenching effects with the Birks' law,
  //              using the expression and the coefficients taken from
  //              the paper NIM 80 (1970) 239-244 for the organic
  //              scintillator NE-102:
  //                                     S*dE/dr
  //                    dL/dr = -----------------------------------
  //                               1 + C1*(dE/dr)
  //              with:
  //                    S=1
  //                    C1 = 1.29 x 10^-2  g*cm^-2*MeV^-1
  //                    C2 = 9.59 x 10^-6  g^2*cm^-4*MeV^-2
  //              These are the same values used by ATLAS TileCal
  //              and CMS HCAL (and also the default in Geant3).
  //              You can try different values for these parameters,
  //              to have an idea on the uncertainties due to them,
  //              by uncommenting one of the lines below.
  //              To get the "dE/dr" that appears in the formula,
  //              which has the dimensions
  //                    [ dE/dr ] = MeV * cm^2 / g
  //              we have to divide the energy deposit in MeV by the
  //              product of the step length (in cm) and the density
  //              of the scintillator:
 
  G4double birks_factor(1.0);
  G4double step_length = aStep->GetStepLength() / CLHEP::cm;
  // Do not apply Birks for gamma deposits!
  // Check, cut if necessary.
  if (step_length > 1.0e-6) {
    G4double rho = aStep->GetPreStepPoint()->GetMaterial()->GetDensity() /
                   (CLHEP::g / CLHEP::cm3);
    G4double dedx = edep / (rho * step_length);  //[MeV*cm^2/g]
    birks_factor = 1.0 / (1.0 + birksc1_ * dedx + birksc2_ * dedx * dedx);
    if (aStep->GetTrack()->GetDefinition() == G4Gamma::GammaDefinition())
      birks_factor = 1.0;
    if (aStep->GetTrack()->GetDefinition() == G4Neutron::NeutronDefinition())
      birks_factor = 1.0;
  }
 
  // update edep to include birksFactor
  edep *= birks_factor;
 
  // Get the scintillator solid box
  G4Box* scint = nullptr;
 
  if (aStep) {
    const auto* pre_step_point = aStep->GetPreStepPoint();
    if (pre_step_point) {
      const auto& touchable_handle = pre_step_point->GetTouchableHandle();
      if (touchable_handle) {
        const auto* volume = touchable_handle->GetVolume();
 
        if (volume) {
          const auto* logical_volume = volume->GetLogicalVolume();
          if (logical_volume) {
            auto* solid = logical_volume->GetSolid();
            if (solid) {
              scint = static_cast<G4Box*>(solid);
            }
          }
        }
      }
    }
  }
 
  // Set the step mid-point as the hit position.
  G4StepPoint* pre_point = aStep->GetPreStepPoint();
  G4StepPoint* post_point = aStep->GetPostStepPoint();
  // A Geant4 "touchable" is a way to uniquely identify a particular volume,
  // short for touchable detector element. See the detector definition and
  // response section of the Geant4 application developers manual for details.
  //
  // The TouchableHandle is just a reference counted pointer to a
  // G4TouchableHistory object, which is a concrete implementation of a
  // G4Touchable interface.
  //
  auto touchable_history{pre_point->GetTouchableHandle()->GetHistory()};
  // Affine transform for converting between local and global coordinates
  auto top_transform{touchable_history->GetTopTransform()};
  G4ThreeVector position =
      0.5 * (pre_point->GetPosition() + post_point->GetPosition());
  G4ThreeVector local_position = top_transform.TransformPoint(position);
 
  // Create the ID for the hit. Note 2 here corresponds to the "depth" of the
  // geometry tree. If this changes in the GDML, this would have to be updated
  // here. Currently, 0 corresponds to the world volume, 1 corresponds to the
  // Hcal, and 2 to the bars/absorbers
  int copy_num = touchable_history->GetVolume(2)->GetCopyNo();
  ldmx::HcalID id = decodeCopyNumber(copy_num, local_position, scint);
 
  if (hits_.find(id) == hits_.end()) {
    // hit in empty cell/bar
    auto& new_hit = hits_[id];
    new_hit.setID(id.raw());
    new_hit.setPosition(position[0], position[1], position[2]);
  }
 
  auto& hit = hits_[id];
 
  // hit variables
  const G4Track* track = aStep->GetTrack();
  auto time = track->GetGlobalTime();
  auto track_id = track->GetTrackID();
  auto pdg = track->GetParticleDefinition()->GetPDGEncoding();
 
  if (enable_hit_contribs_) {
    int contrib_i = hit.findContribIndex(track_id, pdg);
    if (compress_hit_contribs_ and contrib_i != -1) {
      hit.updateContrib(contrib_i, edep, time);
    } else {
      auto map{getTrackMap()};
      auto incident{map.findIncident(track_id)};
      // default "origin" is just the same as incident
      // "origin" checks if a hit "originates" from one of the earliest
      // track IDs (i.e. probably one of the primaries)
      int origin{incident};
      for (int i{1}; i < max_origin_track_id_; ++i) {
        if (map.isDescendant(track_id, i, 100)) {
          origin = i;
          break;
        }
      }
      hit.addContrib(incident, track_id, pdg, edep, time, origin);
    }
  } else {
    // no hit contribs and hit already exists
    hit.setEdep(hit.getEdep() + edep);
    if (time < hit.getTime() or hit.getTime() == 0) {
      hit.setTime(time);
    }
  }
  // Pre/post step details for scintillator response simulation
  // Note: These are set for each step, so the last step's details will be used
  // if multiple steps occur in the same bar
 
  // Convert back to mm
  hit.setPathLength(step_length * CLHEP::cm / CLHEP::mm);
  hit.setVelocity(track->GetVelocity());
  const auto& geometry = getCondition<ldmx::HcalGeometry>(
      ldmx::HcalGeometry::CONDITIONS_OBJECT_NAME);
  // Convert pre/post step position from global coordinates to coordinates
  // within the scintillator bar
  const auto local_pre_step_point{
      top_transform.TransformPoint(pre_point->GetPosition())};
  const auto local_post_step_point{
      top_transform.TransformPoint(post_point->GetPosition())};
 
  // And rotate them to a local coordinate system for the bar that always has
  // the same x/y/z definitions (see HcalGeometry for details)
  auto local_pre_position_rotated{geometry.rotateGlobalToLocalBarPosition(
      {local_pre_step_point[0], local_pre_step_point[1],
       local_pre_step_point[2]},
      id)};
 
  auto local_post_position_rotated{geometry.rotateGlobalToLocalBarPosition(
      {local_post_step_point[0], local_post_step_point[1],
       local_post_step_point[2]},
      id)};
  hit.setPreStepPosition(local_pre_position_rotated[0],
                         local_pre_position_rotated[1],
                         local_pre_position_rotated[2]);
  hit.setPostStepPosition(local_post_position_rotated[0],
                          local_post_position_rotated[1],
                          local_post_position_rotated[2]);
  hit.setPreStepTime(pre_point->GetGlobalTime());
  hit.setPostStepTime(post_point->GetGlobalTime());
 
  ldmx_log(trace) << hit;
 
  return true;
}

References compress_hit_contribs_, ldmx::HcalGeometry::CONDITIONS_OBJECT_NAME, decodeCopyNumber(), enable_hit_contribs_, simcore::SensitiveDetector::getCondition(), simcore::SensitiveDetector::getTrackMap(), hits_, simcore::SensitiveDetector::isGeantino(), and max_origin_track_id_.

◆ saveHits()

void simcore::HcalSD::saveHits ( framework::Event & event )

overridevirtual

Add our hits to the event bus and then reset the container.

Implements simcore::SensitiveDetector.

Definition at line 245 of file HcalSD.cxx.

                                           {
  // squash hits into list
  std::vector<ldmx::SimCalorimeterHit> hits;
  hits.reserve(hits_.size());
  for (const auto& [id, hit] : hits_) hits.push_back(hit);
  event.add(COLLECTION_NAME, hits);
}

References COLLECTION_NAME, and hits_.

Member Data Documentation

◆ birksc1_

double simcore::HcalSD::birksc1_

private

Definition at line 99 of file HcalSD.h.

◆ birksc2_

double simcore::HcalSD::birksc2_

private

Definition at line 102 of file HcalSD.h.

◆ COLLECTION_NAME

const std::string simcore::HcalSD::COLLECTION_NAME = "HcalSimHits"

static

name of collection to be added to event bus

Definition at line 27 of file HcalSD.h.

Referenced by saveHits().

◆ compress_hit_contribs_

bool simcore::HcalSD::compress_hit_contribs_

private

compress hit contribs

Definition at line 109 of file HcalSD.h.

Referenced by HcalSD(), and ProcessHits().

◆ enable_hit_contribs_

bool simcore::HcalSD::enable_hit_contribs_

private

enable hit contribs

Definition at line 107 of file HcalSD.h.

Referenced by HcalSD(), and ProcessHits().

◆ gdml_identifiers_

std::vector<std::string> simcore::HcalSD::gdml_identifiers_

private

Definition at line 97 of file HcalSD.h.

◆ hits_

std::map<ldmx::HcalID, ldmx::SimCalorimeterHit> simcore::HcalSD::hits_

private

map of hits to add to the event (will be squashed)

Definition at line 105 of file HcalSD.h.

Referenced by onFinishedEvent(), ProcessHits(), and saveHits().

◆ max_origin_track_id_

int simcore::HcalSD::max_origin_track_id_

private

maximum track ID to be considered an "origin"

Definition at line 111 of file HcalSD.h.

Referenced by HcalSD(), and ProcessHits().

The documentation for this class was generated from the following files:

SimCore/include/SimCore/SDs/HcalSD.h
SimCore/src/SimCore/SDs/HcalSD.cxx

Public Member Functions

Static Public Attributes

Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ HcalSD()

Member Function Documentation

◆ decodeCopyNumber()

◆ isSensDet()

◆ onFinishedEvent()

◆ ProcessHits()

◆ saveHits()

Member Data Documentation

◆ birksc1_

◆ birksc2_

◆ COLLECTION_NAME

◆ compress_hit_contribs_

◆ enable_hit_contribs_

◆ gdml_identifiers_

◆ hits_

◆ max_origin_track_id_