sTgcDigitMaker Class Reference

Handles the digitization of sTGC hits, converting simulated hits into sTGC digits. More...

#include <sTgcDigitMaker.h>

Inheritance diagram for sTgcDigitMaker:

Collaboration diagram for sTgcDigitMaker:

Classes
struct	DigiConditions
	Digitize a given hit, determining the time and charge spread on wires, pads and strips. More...
struct	GammaParameter
	Parameters of a gamma probability distribution function, required for estimating wire digit's time of arrival. More...
struct	Ionization
	Ionization object with distance, position on the hit segment and position on the wire. More...

Public Types
enum	digitMode { StripsOnly = 1, StripsAndPads = 2, AllChType = 3 }
	Constructor initializing digitization parameters. More...
using	sTgcDigitVec = std::vector< std::unique_ptr< sTgcDigit > >
using	TimedHit = TimedHitPtr< xAOD::MuonSimHit >
using	ReadoutChannelType = sTgcIdHelper::sTgcChannelTypes
using	sTgcDigitVec = std::vector< std::unique_ptr< sTgcDigit > >
	Digitize a given hit. More...

Public Member Functions
	sTgcDigitMaker (const Muon::IMuonIdHelperSvc *idHelperSvc, const int channelTypes, double meanGasGain, bool doPadChargeSharing, bool applyAsBuiltBLines)
virtual	~sTgcDigitMaker ()
StatusCode	initialize ()
	Initializes sTgcHitIdHelper and sTgcIdHelper, and call the functions to read files containing digitization parameters. More...
sTgcDigitVec	executeDigi (const DigiConditions &condContainers, const sTGCSimHit &hit) const
	sTgcDigitMaker (const Muon::IMuonIdHelperSvc *idHelperSvc, digitMode mode, double meanGasGain, bool doPadChargeSharing)
virtual	~sTgcDigitMaker ()
	Destructor. More...
StatusCode	initialize ()
	Initialize digitization parameters, including reading necessary data files. More...
sTgcDigitVec	executeDigi (const DigiConditions &condContainers, const TimedHit &hit) const
bool	msgLvl (const MSG::Level lvl) const
	Test the output level. More...
MsgStream &	msg () const
	The standard message stream. More...
MsgStream &	msg (const MSG::Level lvl) const
	The standard message stream. More...
void	setLevel (MSG::Level lvl)
	Change the current logging level. More...

Private Member Functions
StatusCode	readFileOfTimeArrival ()
	Read share/sTGC_Digitization_timeArrival.dat. More...
StatusCode	readFileOfTimeOffsetStrip ()
	Read share/sTGC_Digitization_timeOffsetStrip.dat. More...
Ionization	pointClosestApproach (const MuonGM::sTgcReadoutElement *readoutEle, const Identifier &id, int wireNumber, const Amg::Vector3D &preStepPos, const Amg::Vector3D &postStepPos) const
	Determine the points where the distance between two segments is smallest. More...
double	getTimeOffsetStrip (size_t neighbor_index) const
	Get digit time offset of a strip depending on its relative position to the strip at the centre of the cluster. More...
GammaParameter	getGammaParameter (double distance) const
	Find the gamma pdf parameters of a given distance. More...
double	getMostProbableArrivalTime (double distance) const
	Get the most probable time of arrival. More...
StatusCode	readFileOfTimeArrival ()
	Reads time arrival data file. More...
StatusCode	readFileOfTimeOffsetStrip ()
	Reads strip time offset data file. More...
Ionization	pointClosestApproach (const MuonGMR4::WireGroupDesign &wireDesign, int wireNumber, const Amg::Vector3D &locHitPos, const Amg::Vector3D &locHitDir, const double stepLength) const
	Computes the closest approach between a trajectory and a wire segment. More...
double	getTimeOffsetStrip (size_t neighbor_index) const
	Gets the time offset for a strip cluster. More...
GammaParameter	getGammaParameter (double distance) const
	Retrieves gamma distribution parameters based on distance. More...
double	getMostProbableArrivalTime (double distance) const
	Computes the most probable arrival time based on the distance of closest approach. More...
double	chargeIntegral (double N, double M) const
void	initMessaging () const
	Initialize our message level and MessageSvc. More...

Static Private Member Functions
static void	addDigit (sTgcDigitVec &digits, const Identifier &id, const uint16_t bctag, const double digittime, const double charge)
static double	getPadChargeFraction (double distance)
static void	addDigit (sTgcDigitVec &digits, const Identifier &id, uint16_t bctag, double digittime, double charge)
	Adds a digit to the appropriate cache. More...
static double	getPadChargeFraction (double distance)
	Computes charge fraction shared among pads. More...

Private Attributes
std::vector< GammaParameter >	m_gammaParameter
std::array< double, 5 >	m_mostProbableArrivalTime {make_array<double, 5>(0.)}
std::array< double, 6 >	m_timeOffsetStrip {make_array<double, 6>(0.)}
const Muon::IMuonIdHelperSvc *	m_idHelperSvc {nullptr}
int	m_channelTypes {3}
	define offsets and widths of time windows for signals from wiregangs and strips. More...
double	m_theta {10}
double	m_meanGasGain {5.e4}
bool	m_doPadSharing {false}
bool	m_doTimeOffsetStrip {false}
double	m_StripResolution {0.0949}
double	m_posResIncident {1.}
double	m_posResAngular {0.305/m_StripResolution}
double	m_chargeAngularFactor {4.0}
bool	m_applyAsBuiltBLines {false}
digitMode	m_digitMode = AllChType
	define offsets and widths of time windows for signals from wiregroups and strips. More...
std::string	m_nm
	Message source name. More...
boost::thread_specific_ptr< MsgStream >	m_msg_tls
	MsgStream instance (a std::cout like with print-out levels) More...
std::atomic< IMessageSvc * >	m_imsg { nullptr }
	MessageSvc pointer. More...
std::atomic< MSG::Level >	m_lvl { MSG::NIL }
	Current logging level. More...
std::atomic_flag m_initialized	ATLAS_THREAD_SAFE = ATOMIC_FLAG_INIT
	Messaging initialized (initMessaging) More...

Static Private Attributes
static constexpr std::array< double, 2 >	m_clusterProfile {0.573, 1.092}
static constexpr std::array< double, 2 >	m_clusterParams {0.573, 1.092}

Detailed Description

Handles the digitization of sTGC hits, converting simulated hits into sTGC digits.

for SCT_ChargeTrappingTool

All functionality of sTGC digitization is implemented to this class.

Definition at line 36 of file MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h.

Member Typedef Documentation

ReadoutChannelType

using sTgcDigitMaker::ReadoutChannelType = sTgcIdHelper::sTgcChannelTypes

Definition at line 72 of file MuonPhaseII/MuonDigitization/sTgcDigitizationR4/sTgcDigitizationR4/sTgcDigitMaker.h.

sTgcDigitVec [1/2]

using sTgcDigitMaker::sTgcDigitVec = std::vector<std::unique_ptr<sTgcDigit> >

Definition at line 64 of file MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h.

sTgcDigitVec [2/2]

using sTgcDigitMaker::sTgcDigitVec = std::vector<std::unique_ptr<sTgcDigit> >

Digitize a given hit.

Parameters

condContainers	Conditions required for digitization.
hit	The simulated hit to be digitized.

Definition at line 79 of file MuonPhaseII/MuonDigitization/sTgcDigitizationR4/sTgcDigitizationR4/sTgcDigitMaker.h.

TimedHit

using sTgcDigitMaker::TimedHit = TimedHitPtr<xAOD::MuonSimHit>

Definition at line 71 of file MuonPhaseII/MuonDigitization/sTgcDigitizationR4/sTgcDigitizationR4/sTgcDigitMaker.h.

Member Enumeration Documentation

digitMode

enum sTgcDigitMaker::digitMode

Constructor initializing digitization parameters.

Enumerator
StripsOnly
StripsAndPads
AllChType

Definition at line 40 of file MuonPhaseII/MuonDigitization/sTgcDigitizationR4/sTgcDigitizationR4/sTgcDigitMaker.h.

                 {
    StripsOnly = 1,
    StripsAndPads = 2,
    AllChType = 3,
  };

Constructor & Destructor Documentation

sTgcDigitMaker() [1/2]

sTgcDigitMaker::sTgcDigitMaker	(	const Muon::IMuonIdHelperSvc *	idHelperSvc,
		const int	channelTypes,
		double	meanGasGain,
		bool	doPadChargeSharing,
		bool	applyAsBuiltBLines
	)

Definition at line 33 of file MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx.

  : AthMessaging ("sTgcDigitMaker"),
  m_idHelperSvc{idHelperSvc},
  m_channelTypes{channelTypes},
  m_meanGasGain{meanGasGain},
  m_doPadSharing{doPadChargeSharing},
  m_applyAsBuiltBLines(applyAsBuiltBLines) {}

~sTgcDigitMaker() [1/2]

sTgcDigitMaker::~sTgcDigitMaker ( )

virtualdefault

sTgcDigitMaker() [2/2]

sTgcDigitMaker::sTgcDigitMaker	(	const Muon::IMuonIdHelperSvc *	idHelperSvc,
		digitMode	mode,
		double	meanGasGain,
		bool	doPadChargeSharing
	)

Definition at line 23 of file MuonPhaseII/MuonDigitization/sTgcDigitizationR4/src/sTgcDigitMaker.cxx.

  : AthMessaging("sTgcDigitMaker"),
    m_idHelperSvc{idHelperSvc},
    m_digitMode(mode),
    m_meanGasGain{meanGasGain},
    m_doPadSharing{doPadChargeSharing}{}

~sTgcDigitMaker() [2/2]

virtual sTgcDigitMaker::~sTgcDigitMaker ( )

virtual

Destructor.

Member Function Documentation

addDigit() [1/2]

void sTgcDigitMaker::addDigit ( sTgcDigitVec &  digits, const Identifier &  id, const uint16_t  bctag, const double  digittime, const double  charge  )

staticprivate

Definition at line 735 of file MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx.

                                                   {
  
  constexpr double tolerance = 0.1;
  if (std::find_if(digits.begin(),digits.end(), [&](std::unique_ptr<sTgcDigit>& known) {
      return known->identify() == id && std::abs(digittime - known->time()) < tolerance;
  }) == digits.end()) {
    digits.push_back(std::make_unique<sTgcDigit>(id, bctag, digittime, charge, 0, 0));
  }
}

addDigit() [2/2]

static void sTgcDigitMaker::addDigit ( sTgcDigitVec &  digits, const Identifier &  id, uint16_t  bctag, double  digittime, double  charge  )

staticprivate

Adds a digit to the appropriate cache.

chargeIntegral()

double sTgcDigitMaker::chargeIntegral ( double  N, double  M  ) const

private

Definition at line 763 of file MuonPhaseII/MuonDigitization/sTgcDigitizationR4/src/sTgcDigitMaker.cxx.

                                                              {
 
    double term1 = 0.25 * std::erf( M / (std::sqrt(2) * m_clusterParams[0]));
    double term2 = 0.25 * std::erf( N / (std::sqrt(2) * m_clusterParams[0]));
    double term3 = 0.25 * std::erf( M / (std::sqrt(2) * m_clusterParams[1]));
    double term4 = 0.25 * std::erf( N / (std::sqrt(2) * m_clusterParams[1]));
 
    return (term1 - term2 + term3 - term4);
}

executeDigi() [1/2]

sTgcDigitMaker::sTgcDigitVec sTgcDigitMaker::executeDigi	(	const DigiConditions &	condContainers,
		const sTGCSimHit &	hit
	)		const

Definition at line 69 of file MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx.

                                                                                      {
  // SimHits without energy loss are not recorded.
  double energyDeposit = hit.depositEnergy(); // Energy deposit in MeV
  if(energyDeposit==0.) return {};
 
  const sTgcIdHelper& idHelper{m_idHelperSvc->stgcIdHelper()};
  sTgcSimIdToOfflineId simToOffline{&idHelper};
  const Identifier offlineId = simToOffline.convert(hit.sTGCId());
  const Identifier layid = idHelper.channelID(offlineId, 
                                              idHelper.multilayer(offlineId),
                                              idHelper.gasGap(offlineId), 
                                              sTgcIdHelper::sTgcChannelTypes::Wire, 1);
 
  ATH_MSG_VERBOSE("sTgc hit:  time " << hit.globalTime() 
              << " position " << Amg::toString(hit.globalPosition(), 2) 
              << " mclink " << hit.particleLink() << " PDG ID " << hit.particleEncoding() );
 
  
  ATH_MSG_DEBUG("Retrieving detector element for: "<< m_idHelperSvc->toStringDetEl(layid) << " energyDeposit "<<energyDeposit );
 
  const MuonGM::sTgcReadoutElement* detEl = condContainers.detMgr->getsTgcReadoutElement(layid);
  if( !detEl ){ return {}; }
  
 
  // DO THE DIGITIZATTION HERE ////////
 
  // Required precision on length in mm
  constexpr double length_tolerance = 0.01;
 
  // Retrieve the wire surface
  int surfHash_wire = detEl->surfaceHash(layid);
  const Trk::PlaneSurface& SURF_WIRE = detEl->surface(surfHash_wire); // get the wire surface
 
  // Hit post-Step global position
  const Amg::Vector3D& GPOS{hit.globalPosition()};
  // Hit pre-Step global position
  const Amg::Vector3D& pre_pos{hit.globalPrePosition()};
  // Hit global direction
  const Amg::Vector3D& GLODIRE{hit.globalDirection()};
 
  // Hit position in the wire surface's coordinate frame
  Amg::Vector3D hitOnSurface_wire{SURF_WIRE.transform().inverse() * GPOS};
  Amg::Vector3D pre_pos_wire_surf{SURF_WIRE.transform().inverse() * pre_pos};
  Amg::Vector2D posOnSurf_wire{0.5 * (hitOnSurface_wire + pre_pos_wire_surf).block<2,1>(0,0)};
 
  /* Determine the closest wire and the distance of closest approach
   * Since most particles pass through the the wire plane between two wires,
   * the nearest wire should be one of these two wire. Otherwise, the particle's
   * trajectory is uncommon, and such rare case is not supported yet.
   *
   * Finding that nearest wire follows the following steps:
   * - Compute the distance to the wire at the center of the current wire pitch
   * - Compute the distance to the other adjacent wire
   */
 
  // Wire number of the current wire pitch
  int wire_number = detEl->getDesign(layid)->wireNumber(posOnSurf_wire);
 
  // If wire number is invalid, verify if hit is near the edge of the chamber.
  const int number_wires = detEl->numberOfWires(layid);
  if ((wire_number < 1) || (wire_number > number_wires)) {
    // Compute the wire number using either the pos-step position or
    // the pre-step position, whichever yields a valid wire number.
    wire_number = detEl->getDesign(layid)->wireNumber(hitOnSurface_wire.block<2,1>(0,0));
    if ((wire_number < 1) || (wire_number > number_wires)) {
      wire_number = detEl->getDesign(layid)->wireNumber(pre_pos_wire_surf.block<2,1>(0,0));
    }
  }
 
  // Compute the position of the ionization and its distance relative to a given sTGC wire.
  Ionization ionization = pointClosestApproach(detEl, layid, wire_number, pre_pos_wire_surf, hitOnSurface_wire);
  double dist_wire = ionization.distance;
  if (dist_wire > 0.) {
    // Determine the other adjacent wire, which is
    //  -1 if particle crosses the wire surface between wire_number-1 and wire_number
    //  +1 if particle crosses the wire surface between wire_number and wire_number+1
    int adjacent = 1;
    if (ionization.posOnSegment.x() < ionization.posOnWire.x()) {adjacent = -1;}
 
    // Find the position of the ionization with respect to the adjacent wire
    Ionization ion_adj = pointClosestApproach(detEl, layid, wire_number+adjacent, pre_pos_wire_surf, hitOnSurface_wire);
    // Keep the closest point
    double dist_wire_adj = ion_adj.distance;
    if ((dist_wire_adj > 0.) && (dist_wire_adj < dist_wire)) {
      dist_wire = dist_wire_adj;
      wire_number += adjacent;
      ionization = std::move(ion_adj);
    }
  } else {
    ATH_MSG_DEBUG("Failed to get the distance between the wire number = " << wire_number
                    << " and hit at " <<Amg::toString(hitOnSurface_wire, 2)
                    << ". Number of wires = " << number_wires
                    << ", "<<m_idHelperSvc->toStringGasGap(layid));
    return {};
  }
 
  // Update the position of ionization on the wire surface
  posOnSurf_wire = ionization.posOnWire.block<2,1>(0,0);
  // Position of the ionization in the global coordinate frame
  const Amg::Vector3D glob_ionization_pos = SURF_WIRE.transform() * ionization.posOnWire;
 
  ATH_MSG_VERBOSE("Ionization_info: distance: " << ionization.distance
    << " posOnTrack: " <<Amg::toString(ionization.posOnSegment, 3)
    << " posOnWire: " << Amg::toString(ionization.posOnWire, 3)
    << " hitGPos: " << Amg::toString(GPOS,3)
    << " hitPrePos: " <<Amg::toString(pre_pos, 3)
    << " EDep: " << hit.depositEnergy() << " EKin: " << hit.kineticEnergy()
    << " pdgId: " << hit.particleEncoding()
    <<" "<<m_idHelperSvc->toStringGasGap(layid));
 
  // Distance should be in the range [0, 0.9] mm, excepting
  // - particles pass through the wire plane near the edges
  // - secondary particles created inside the gas gap that go through the gas gap partially.
  //   Most of such particles are not muons and have low kinetic energy.
  // - particle with trajectory parallel to the sTGC wire plane
  const double wire_pitch = detEl->wirePitch();
  if ((dist_wire > 0.) && (std::abs(hit.particleEncoding()) == 13) && (dist_wire > (wire_pitch/2))) {
    ATH_MSG_DEBUG("Distance to the nearest wire (" << dist_wire << ") is greater than expected.");
    ATH_MSG_DEBUG("Hit globalPos: "<<Amg::toString(GPOS,2)
                 << " globalPrePos: " << Amg::toString(pre_pos, 2)
                 << " EDeposited: " << hit.depositEnergy() 
                 << " EKinetic: " << hit.kineticEnergy()
                 << " pdgID: " << hit.particleEncoding()
                 << " "<<m_idHelperSvc->toStringGasGap(layid));
  }
 
  // Do not digitize hits that are too far from the nearest wire
  if (dist_wire > wire_pitch) {
    return {};
  }
 
  // Get the gamma pdf parameters associated with the distance of closest approach.
  const GammaParameter gamParam = getGammaParameter(dist_wire);
  const double par_kappa = gamParam.kParameter;
  const double par_theta = gamParam.thetaParameter;
  const double most_prob_time = getMostProbableArrivalTime(dist_wire);
  // Compute the most probable value of the gamma pdf
  double gamma_mpv = (par_kappa - 1) * par_theta;
  // If the most probable value is less than zero, then set it to zero
  if (gamma_mpv < 0.) {gamma_mpv = 0.;}
  const double t0_par = most_prob_time - gamma_mpv;
 
  // Digit time follows a gamma distribution, so a value val is
  // chosen using a gamma random generator then is shifted by t0
  // to account for drift time.
  // Note: CLHEP::RandGamma takes the parameters k and lambda,
  // where lambda = 1 / theta.
  double digit_time = t0_par + CLHEP::RandGamma::shoot(condContainers.rndmEngine, par_kappa, 1/par_theta);
 
  // Sometimes, digit_time is negative because t0_par can be negative.
  // In such case, discard the negative value and shoot RandGamma for another value.
  // However, if that has already been done many times then set digit_time to zero
  // in order to avoid runaway loop.
  constexpr int shoot_limit = 4;
  int shoot_counter = 0;
  while (digit_time < 0.) {
    if (shoot_counter > shoot_limit) {
      digit_time = 0.;
      break;
    }
    digit_time = t0_par + CLHEP::RandGamma::shoot(condContainers.rndmEngine, par_kappa, 1/par_theta);
    ++shoot_counter;
  }
 
  ATH_MSG_DEBUG("sTgcDigitMaker distance = " << dist_wire
                << ", time = " << digit_time
                << ", k parameter = " << par_kappa
                << ", theta parameter = " << par_theta
                << ", most probable time = " << most_prob_time);
 
  bool isValid = false;
  if (condContainers.efficiencies) {
    const double efficiency = condContainers.efficiencies->getEfficiency(layid);
    // Lose Hits to match HV efficiency
    if (CLHEP::RandFlat::shoot(condContainers.rndmEngine,0.0,1.0) > efficiency) return {};
  }
 
  
  sTgcDigitVec digits{};
 
  //const double stripPropagationTime = 3.3*CLHEP::ns/CLHEP::m * detEl->distanceToReadout(posOnSurf_strip, elemId); // 8.5*ns/m was used until MC10.
  const double stripPropagationTime = 0.; // 8.5*ns/m was used until MC10.
 
  double sDigitTimeWire = digit_time;
  double sDigitTimePad = sDigitTimeWire;
  double sDigitTimeStrip = sDigitTimeWire + stripPropagationTime;
 
  uint16_t bctag = 0;
 
  //##################################################################################
  //######################################### strip readout ##########################
  //##################################################################################
  ATH_MSG_DEBUG("sTgcDigitMaker::strip response ");
  int channelType = sTgcIdHelper::sTgcChannelTypes::Strip;
 
  Identifier newId = idHelper.channelID(layid,idHelper.multilayer(layid), 
                                        idHelper.gasGap(layid), channelType, 1, isValid);
 
  // get strip surface
  int surfHash_strip =  detEl->surfaceHash(newId);
  const Trk::PlaneSurface& SURF_STRIP = detEl->surface(surfHash_strip); // get the strip surface
 
  const Amg::Vector3D hitOnSurface_strip = SURF_STRIP.transform().inverse()*glob_ionization_pos;
  Amg::Vector2D posOnSurf_strip {Amg::Vector2D::Zero()};
 
  if(m_applyAsBuiltBLines){
    //This block is used to apply As-Built and BLine corrections for dedicated studies.
    Amg::Vector3D posAfterAsBuilt {Amg::Vector3D::Zero()};
    detEl->spacePointPosition(newId, hitOnSurface_strip.x(), hitOnSurface_strip.y(), posAfterAsBuilt);
    posOnSurf_strip = posAfterAsBuilt.block<2,1>(0,0);
  } else {
    posOnSurf_strip = hitOnSurface_strip.block<2,1>(0,0);
  }
 
  bool insideBounds = SURF_STRIP.insideBounds(posOnSurf_strip);
  if(!insideBounds) {
    ATH_MSG_DEBUG("Outside of the strip surface boundary : " <<  m_idHelperSvc->toString(newId) << "; local position " <<posOnSurf_strip );
    return {};
  }
 
  //************************************ find the nearest readout element **************************************
 
  int stripNumber = detEl->stripNumber(posOnSurf_strip, newId);
  if( stripNumber == -1 ){
    // Verify if the energy deposit is at the boundary
    const double new_posX = (posOnSurf_strip.x() > 0.0)? posOnSurf_strip.x() - length_tolerance
                                                      : posOnSurf_strip.x() + length_tolerance;
    const Amg::Vector2D new_position(new_posX, posOnSurf_strip.y());
    stripNumber = detEl->stripNumber(new_position, newId);
    // Skip hit if still unable to obtain strip number
    if (stripNumber < 1) {
      ATH_MSG_WARNING("Failed to obtain strip number " << m_idHelperSvc->toString(newId) );
      ATH_MSG_WARNING("Position on strip surface = (" << posOnSurf_strip.x() << ", " << posOnSurf_strip.y() << ")");
      return {};
    }
  }
  newId = idHelper.channelID(layid, idHelper.multilayer(layid), 
                             idHelper.gasGap(layid), channelType, stripNumber, isValid);
  if(!isValid && stripNumber != -1) {
    ATH_MSG_ERROR("Failed to obtain identifier " << m_idHelperSvc->toString(newId) );
    return {};
  }
 
  const int NumberOfStrips = detEl->numberOfStrips(newId);
  const double stripHalfPitch = detEl->channelPitch(newId)*0.5; // 3.2/2 = 1.6 mm
 
  //************************************ conversion of energy to charge **************************************
 
  // Typical ionized charge in pC per keV deposited. The constant is determined from ionization
  // study with Garfield program. A note titled "Charge Energy Relation" which outlines
  // conversion can be found here:
  //   https://cernbox.cern.ch/index.php/apps/files/?dir=/__myshares/Documents (id:274113) // link is dead
  const double ionized_charge = (5.65E-6)*energyDeposit/CLHEP::keV;
 
 // To get avalanche gain, polya function is taken from Blum paper https://inspirehep.net/literature/807304
 // m_polyaFunction = new TF1("m_polyaFunction","(1.0/[1])*(TMath::Power([0]+1,[0]+1)/TMath::Gamma([0]+1))*TMath::Power(x/[1],[0])*TMath::Exp(-([0]+1)*x/[1])",0,3000000);
 
  // Mean value for total gain due to E field;
  // To calculate this gain from polya distibution, we replace in gamma PDF:
  //     alpha = 1+theta and
  //     beta = 1+theta/mean
  // With these substitutions, gamma PDF gives the same sampling values as those from polya PDF.
  const double gain =  CLHEP::RandGamma::shoot(condContainers.rndmEngine, 1. + m_theta, (1. + m_theta)/m_meanGasGain);
 
  // total charge after avalanche
  const double total_charge = gain*ionized_charge;
 
  //************************************ spread charge among readout element **************************************
  // Charge Spread including tan(theta) resolution term.
  const double tan_theta = GLODIRE.perp()/GLODIRE.z();
  // The angle dependance on strip resolution goes as tan^2(angle)
  const double angle_dependency = std::hypot(m_posResIncident, m_posResAngular * tan_theta);
 
  const double cluster_posX = posOnSurf_strip.x();
  double peak_position = CLHEP::RandGaussZiggurat::shoot(condContainers.rndmEngine, cluster_posX, m_StripResolution*angle_dependency);
 
  // Each readout plane reads about half the total charge produced on the wire,
  // including a tan(theta) term to describe the increase of charge with incident angle
  const double norm = 0.5 * total_charge;
  // Strip cluster charge profile described by a double Gaussian.
  // The independent parameter x is in the units of strip channel, one strip channel = 3.2 mm,
  //   so convert position from mm to strip channel if it is not already.
std::unique_ptr<TF1> clusterProfile = std::make_unique<TF1>("fgaus",
    "0.5 * [2] / (sqrt(2 * TMath::Pi()) * [0]) * exp(-0.5 * (x / [0])^2) + "
    "0.5 * [2] / (sqrt(2 * TMath::Pi()) * [1]) * exp(-0.5 * (x / [1])^2)", -300., 300.);
  clusterProfile->SetParameters(m_clusterProfile[0], // sigma of 1st Gaussian
                                m_clusterProfile[1], // sigma of 2nd Gaussian
                                norm);               // normalization factor
 
  // Lower limit on strip charge (arbitrary limit), in pC, which has the same units as the parameter ionized_charge. 
  constexpr double tolerance_charge = 0.0005;
  // Set a maximum number of neighbour strips to avoid very long loop. This is an arbitrary number.
  constexpr unsigned int max_neighbor = 10;
 
  // Spread charge on the strips that are on the upper half of the strip cluster
  for (unsigned int iStrip = 0; iStrip <= max_neighbor; ++iStrip) {
    int currentStrip = stripNumber + iStrip;
    if (currentStrip > NumberOfStrips) break;
 
    // Get the strip identifier and create the digit
    newId = idHelper.channelID(layid, idHelper.multilayer(layid), idHelper.gasGap(layid), 
                               channelType, currentStrip, isValid);
    if (isValid) {
      Amg::Vector2D locpos{Amg::Vector2D::Zero()};
      if (!detEl->stripPosition(newId, locpos)) {
          ATH_MSG_WARNING("Failed to obtain local position for identifier " << m_idHelperSvc->toString(newId) );
      }
 
      // Estimate the digit charge
      // Position with respect to the peak of the charge curve
      double x_relative = locpos.x() - peak_position;
      // In clusterProfile curve, position should be in the units of strip channel
      double charge = std::hypot(1, m_chargeAngularFactor * tan_theta) * clusterProfile->Integral(x_relative/(2*stripHalfPitch) - 0.5, x_relative/(2*stripHalfPitch) + 0.5);
      // If charge is too small, stop creating neighbor strip
      if (charge < tolerance_charge) break;
 
      // Estimate digit time
      double strip_time = sDigitTimeStrip;
      // Strip time response can be delayed due to the resistive layer.
      // A correction would be required if the actual VMM front-end doesn't re-align the strip timing.
      if (m_doTimeOffsetStrip) {
        // Determine how far the current strip is from the middle strip
        int indexFromMiddleStrip = iStrip;
        if ((iStrip > 0) && ((x_relative + (0.75 - iStrip * 2) * stripHalfPitch) < 0))
           indexFromMiddleStrip = iStrip - 1;
        // Add time delay due to resistive layer
        strip_time += getTimeOffsetStrip(indexFromMiddleStrip);
      }
 
      addDigit(digits, newId, bctag, strip_time, charge);
  
      ATH_MSG_VERBOSE("Created a strip digit: strip number = " << currentStrip << ", charge = " << charge
                      << ", time = " << strip_time << ", time offset = " << strip_time-sDigitTimeStrip
                      << ", neighbor index = " << iStrip
                      << ", strip position = "<<Amg::toString(locpos, 2));
    }
  }
 
  // The lower half of the strip cluster
  for (unsigned int iStrip = 1; iStrip <= max_neighbor; ++iStrip) {
    int currentStrip = stripNumber - iStrip;
    if (currentStrip < 1) break;
 
    newId = idHelper.channelID(layid,idHelper.multilayer(layid), 
                              idHelper.gasGap(layid), channelType, currentStrip, isValid);
    if (isValid) {
      Amg::Vector2D locpos{Amg::Vector2D::Zero()};
      if (!detEl->stripPosition(newId, locpos)) {
        ATH_MSG_WARNING("Failed to obtain local position for identifier " << m_idHelperSvc->toString(newId) );
      }
 
      // Estimate the digit charge
      double x_relative = locpos.x() - peak_position;
      double charge = std::hypot(1, m_chargeAngularFactor * tan_theta) * clusterProfile->Integral(x_relative/(2*stripHalfPitch) - 0.5, x_relative/(2*stripHalfPitch) + 0.5);
      if (charge < tolerance_charge) break;
 
      // Estimate digit time
      double strip_time = sDigitTimeStrip;
      // Time delay due to resistive layer
      if (m_doTimeOffsetStrip) {
        int indexFromMiddleStrip = ((x_relative + (iStrip * 2 - 0.75) * stripHalfPitch) > 0)? iStrip - 1 : iStrip;
        strip_time += getTimeOffsetStrip(indexFromMiddleStrip);
      }
 
      addDigit(digits, newId, bctag, strip_time, charge);
 
      ATH_MSG_VERBOSE("Created a strip digit: strip number = " << currentStrip << ", charge = " << charge
                      << ", time = " << strip_time << ", time offset = " << strip_time-sDigitTimeStrip
                      << ", neighbor index = " << iStrip
                      << ", strip position = " << Amg::toString(locpos, 2));
    }
  }
  // end of strip digitization
 
  if(m_channelTypes==1) {
    ATH_MSG_WARNING("Only digitize strip response !");
    return digits;
  }
 
  //##################################################################################
  //######################################### pad readout ##########################
  //##################################################################################
  ATH_MSG_DEBUG("sTgcDigitMaker::pad response ");
  channelType = sTgcIdHelper::sTgcChannelTypes::Pad;
 
  Identifier PAD_ID = idHelper.channelID(layid, idHelper.multilayer(layid), 
                                        idHelper.gasGap(layid), channelType, 1);// find the a pad id
  //************************************ find the nearest readout element **************************************
  int  surfHash_pad =  detEl->surfaceHash(PAD_ID);
  const Trk::PlaneSurface& SURF_PAD = detEl->surface(surfHash_pad); // get the pad surface
 
  const Amg::Vector3D hitOnSurface_pad = SURF_PAD.transform().inverse()*glob_ionization_pos;
  const Amg::Vector2D posOnSurf_pad(hitOnSurface_pad.x(), hitOnSurface_pad.y());
 
  
  insideBounds = SURF_PAD.insideBounds(posOnSurf_pad);
 
  if(insideBounds) {
    int padNumber = detEl->stripNumber(posOnSurf_pad, PAD_ID);
    if( padNumber == -1 ){
      // Verify if the energy deposit is at the boundary
      const double new_posX = (posOnSurf_pad.x()>0.0)? posOnSurf_pad.x()-length_tolerance
                                                    : posOnSurf_pad.x()+length_tolerance;
      const double new_posY = (posOnSurf_pad.y()>0.0)? posOnSurf_pad.y()-length_tolerance
                                                    : posOnSurf_pad.y()+length_tolerance;
      const Amg::Vector2D new_position(new_posX, new_posY);
      padNumber = detEl->stripNumber(new_position, PAD_ID);
      // Skip hit if still unable to obtain pad number
      if (padNumber < 1) {
        ATH_MSG_WARNING("Failed to obtain pad number " << m_idHelperSvc->toString(PAD_ID) );
        ATH_MSG_WARNING("Position on pad surface = (" << posOnSurf_pad.x() << ", " << posOnSurf_pad.y() << ")");
        return digits;
      }
    }
    newId = idHelper.channelID(layid, idHelper.multilayer(layid), 
                               idHelper.gasGap(layid), channelType, padNumber, isValid);
    if(isValid) {
      // Find centre position of pad
      Amg::Vector2D padPos{Amg::Vector2D::Zero()};
      if (!detEl->stripPosition(newId, padPos)) {
        ATH_MSG_ERROR("Failed to obtain local position for neighbor pad with identifier " << m_idHelperSvc->toString(newId) );
        return digits;
      }
      // Pad sharing needs to look at position on hit vs pad boundaries
      const Amg::Vector2D diff = posOnSurf_pad - padPos;
      double halfPadWidthX = 0.5*detEl->getPadDesign(newId)->channelWidth(posOnSurf_pad, true, true);
      double halfPadWidthY = 0.5*detEl->getPadDesign(newId)->channelWidth(posOnSurf_pad, false);
 
      // Charge sharing happens within 4mm window for pads
      // i.e. the charge is spread over a 2D gaussian of total radius about 4mm
      // This value depends on actual width of the distribution, but we don't have the
      // width of the charge distribution for pads. So lets consider 2.5*sigma, where
      // sigma is the width of the charge distribution for strips.
      double deltaX = halfPadWidthX - std::abs(diff.x());
      double deltaY= halfPadWidthY - std::abs(diff.y());
      bool isNeighX = deltaX < 2.5*m_clusterProfile[1];
      bool isNeighY = deltaY < 2.5*m_clusterProfile[1];
      // Pad width can be calculated to be very slightly larger than it should due to rounding errors
      // So if a hit falls on a given pad but is "outside" the width, just define it to be on the boundary of 2 pads.
      if (deltaX < 0.) deltaX = 0.1;
      if (deltaY < 0.) deltaY = 0.1;
 
      if (m_doPadSharing && (isNeighX || isNeighY)){
        std::pair<int, int> padEtaPhi(detEl->getPadDesign(newId)->channelNumber(padPos));
        // Phi == 1 at the right in the local geometry
        // In local coordinates, the pad to the right has phi' = phi - 1
        int newPhi = diff.x() > 0 ? padEtaPhi.second-1 : padEtaPhi.second+1;
        int newEta = diff.y() > 0 ? padEtaPhi.first+1 : padEtaPhi.first-1;
        bool validEta = newEta > 0 && newEta < detEl->getPadDesign(newId)->nPadH+1;
        bool validPhi = newPhi > 0 && newPhi < detEl->getPadDesign(newId)->nPadColumns+1;
    
        if (isNeighX && isNeighY && validEta && validPhi){
          // 4 pads total; makes life a bit harder. Corner of 4 valid pads
          Identifier neigh_ID_X = idHelper.padID(newId, idHelper.multilayer(newId), idHelper.gasGap(newId),
                                                  channelType, padEtaPhi.first, newPhi, isValid);
          Identifier neigh_ID_Y = idHelper.padID(newId, idHelper.multilayer(newId), idHelper.gasGap(newId),
                                                  channelType, newEta, padEtaPhi.second, isValid);
          Identifier neigh_ID_XY = idHelper.padID(newId, idHelper.multilayer(newId), idHelper.gasGap(newId),
                                                  channelType, newEta, newPhi, isValid);
          double xQfraction = getPadChargeFraction(deltaX);
          double yQfraction = getPadChargeFraction(deltaY);
 
          // Main pad gets 1 - Qfraction of total
          addDigit(digits, neigh_ID_X, bctag, sDigitTimePad, xQfraction*(1.-yQfraction)*0.5*total_charge);
          addDigit(digits, neigh_ID_Y, bctag, sDigitTimePad, yQfraction*(1.-xQfraction)*0.5*total_charge);
          addDigit(digits, neigh_ID_XY, bctag, sDigitTimePad, xQfraction*yQfraction*0.5*total_charge);
          addDigit(digits, newId, bctag, sDigitTimePad, (1.-xQfraction-yQfraction+xQfraction*yQfraction)*0.5*total_charge);
        } else if (isNeighX && validPhi){
          // There is only 1 neighbor, immediately to the left or right.
          Identifier neigh_ID = idHelper.padID(newId, idHelper.multilayer(newId), idHelper.gasGap(newId), 
                                               channelType, padEtaPhi.first, newPhi, isValid);
 
          // NeighborPad gets Qfraction of the total pad charge: 0.5*total_charge
          // Main pad gets 1 - Qfraction of total
          double xQfraction = getPadChargeFraction(deltaX);
          addDigit(digits, newId, bctag, sDigitTimePad, (1.-xQfraction)*0.5*total_charge);
          addDigit(digits, neigh_ID, bctag, sDigitTimePad, xQfraction*0.5*total_charge);
        }
        else if (isNeighY && validEta){
          // There is only 1 neighbor, immediately above or below
          Identifier neigh_ID = idHelper.padID(newId, idHelper.multilayer(newId), idHelper.gasGap(newId), 
                                               channelType, newEta, padEtaPhi.second, isValid);
 
          // NeighborPad gets Qfraction of the total pad charge: 0.5*total_charge
          // Main pad gets 1 - Qfraction of total
          double yQfraction = getPadChargeFraction(deltaX);
          addDigit(digits, newId, bctag, sDigitTimePad, (1.-yQfraction)*0.5*total_charge);
          addDigit(digits, neigh_ID, bctag, sDigitTimePad, yQfraction*0.5*total_charge);
        }
 
      }
      else{
        // No charge sharing: hit is nicely isolated within pad
        addDigit(digits, newId, bctag, sDigitTimePad, 0.5*total_charge);
      }
 
    }
    else if(padNumber != -1) {
      ATH_MSG_ERROR("Failed to obtain identifier " << m_idHelperSvc->toString(newId) );
    }
  }
  else {
    ATH_MSG_DEBUG("Outside of the pad surface boundary :" << m_idHelperSvc->toString(PAD_ID)<< " local position " <<posOnSurf_pad );
  }
 
  if(m_channelTypes==2) {
    ATH_MSG_WARNING("Only digitize strip/pad response !");
    return digits;
  }
 
 
  //##################################################################################
  //######################################### wire readout ##########################
  //##################################################################################
  ATH_MSG_DEBUG("sTgcDigitMaker::wire response ");
  channelType = sTgcIdHelper::sTgcChannelTypes::Wire;
 
    // Find the ID of the first wiregroup
  Identifier WIREGP_ID = idHelper.channelID(layid, idHelper.multilayer(layid),
                                            idHelper.gasGap(layid), channelType, 1);
 
   //************************************ find the nearest readout element **************************************
   insideBounds = SURF_WIRE.insideBounds(posOnSurf_wire);
 
    if(insideBounds) {
        // Determine the wire number
 
        int wiregroupNumber = detEl->stripNumber(posOnSurf_wire, WIREGP_ID);
        if( wiregroupNumber == -1 ){
          // Verify if the energy deposit is at the boundary
          const double new_posX = (posOnSurf_wire.x() > 0.0)? posOnSurf_wire.x() - length_tolerance
                                                           : posOnSurf_wire.x() + length_tolerance;
          const Amg::Vector2D new_position(new_posX, posOnSurf_wire.y());
          wiregroupNumber = detEl->stripNumber(new_position, WIREGP_ID);
          // Skip hit if still unable to obtain pad number
          if (wiregroupNumber < 1) {
            ATH_MSG_WARNING("Failed to obtain wire number " << m_idHelperSvc->toString(WIREGP_ID) );
            ATH_MSG_WARNING("Position on wire surface = (" << posOnSurf_wire.x() << ", " << posOnSurf_wire.y() << ")");
            return digits;
          }
        }
 
        // Find ID of the actual wiregroup
        newId = idHelper.channelID(layid, idHelper.multilayer(layid), idHelper.gasGap(layid), 
                                   channelType, wiregroupNumber, isValid);
 
        if(isValid) {
          int NumberOfWiregroups = detEl->numberOfStrips(newId);
          if(wiregroupNumber>=1&&wiregroupNumber<=NumberOfWiregroups) {
            addDigit(digits, newId, bctag, sDigitTimeWire, total_charge);
          }
        } // end of if(isValid)
        else if (wiregroupNumber != -1){
          ATH_MSG_ERROR("Failed to obtain wiregroup identifier " << m_idHelperSvc->toString(newId) );
        }
    }
    else {
      ATH_MSG_DEBUG("Outside of the wire surface boundary :" << m_idHelperSvc->toString(WIREGP_ID)<< " local position " <<posOnSurf_wire );
    }
    // end of wire digitization
 
  return digits;
}

executeDigi() [2/2]

sTgcDigitMaker::sTgcDigitVec sTgcDigitMaker::executeDigi	(	const DigiConditions &	condContainers,
		const TimedHit &	hit
	)		const

Definition at line 53 of file MuonPhaseII/MuonDigitization/sTgcDigitizationR4/src/sTgcDigitMaker.cxx.

                                                                                                                      {
  // Extract energy deposit from the hit
  double energyDeposit = hit->energyDeposit();
  if (energyDeposit == 0.) return {};  // Ignore hits with no energy deposit
 
  // Retrieve the detector element for the given hit
  const MuonGMR4::MuonDetectorManager* detMgr = condContainers.detMgr;
 
  // Convert hit to identifier and retrieve sTgc readout element
  const Identifier hitId = hit->identify();
  ATH_MSG_DEBUG("Retrieving detector element for: "<< m_idHelperSvc->toStringDetEl(hitId) << " energyDeposit "<< energyDeposit );
  const MuonGMR4::sTgcReadoutElement* readoutElement = detMgr->getsTgcReadoutElement(hitId);
 
  // Projecting the hit position on wire surface given the hit position
  // and direction w.r.t. the wire plane
  Amg::Vector3D locHitDir = xAOD::toEigen(hit->localDirection());
  Amg::Vector3D locHitPos = xAOD::toEigen(hit->localPosition());
  ATH_MSG_VERBOSE("sTgc hit:  time " << hit->globalTime() 
  << " position " << Amg::toString(locHitPos, 2) 
  << " direction" << Amg::toString(locHitDir, 2) 
  << " mclink " << hit->genParticleLink() << " PDG ID " << hit->pdgId() );
  ATH_MSG_VERBOSE("Projecting hit to Wire Surface" );
  const double scale = Amg::intersect<3>(locHitPos, locHitDir, Amg::Vector3D::UnitZ(), 0.).value_or(0);
  Amg::Vector3D hitOnWireSurf = locHitPos + scale * locHitDir;
  ATH_MSG_VERBOSE("Hit on Wire Surface: " << Amg::toString(hitOnWireSurf, 2));
 
  /* Determine the closest wire and the distance of closest approach
   * Since most particles pass through the the wire plane between two wires,
   * the nearest wire should be one of these two wire. Otherwise, the particle's
   * trajectory is uncommon, and such rare case is not supported yet.
   *
   * Finding that nearest wire follows the following steps:
   * - Compute the distance to the wire at the center of the current wire pitch
   * - Compute the distance to the other adjacent wire
   */
  const sTgcIdHelper& idHelper{m_idHelperSvc->stgcIdHelper()};
  const Identifier wireLayId = idHelper.channelID(hitId, 
                                                idHelper.multilayer(hitId),
                                                idHelper.gasGap(hitId),
                                                ReadoutChannelType::Wire,
                                                1);
 
  const MuonGMR4::WireGroupDesign& wireDesign{readoutElement->wireDesign(wireLayId)};
  Amg::Vector2D hitOnWire2D = hitOnWireSurf.block<2,1>(0,0);
  if(!wireDesign.insideTrapezoid(hitOnWire2D)) {
    return {};
  }
  std::pair<int, int> wireGrpWireNum = wireDesign.wireNumber(hitOnWire2D);
  int wireNumber = -1;
  wireNumber = wireDesign.numPitchesToGroup(wireGrpWireNum.first) + wireGrpWireNum.second;
  // If wire number is invalid, verify if hit is near the edge of the chamber.
  // Try to get the wire number from the x and y coordinates of the hit position rather than the hit on wire surface.
  const int numWires = wireDesign.nAllWires();
  if((wireNumber < 1) || (wireNumber > numWires)) {
    wireGrpWireNum = wireDesign.wireNumber(locHitPos.block<2,1>(0,0));
    wireNumber = wireDesign.numPitchesToGroup(wireGrpWireNum.first) + wireGrpWireNum.second;
    if((wireNumber < 1) || (wireNumber > numWires)) {
      ATH_MSG_WARNING("Unable to retrieve the wire number, skipping the hit: " << m_idHelperSvc->toString(hitId));
      return {};
    }
  }
  // Compute the position of the ionization and its distance to the closest wire.
  Ionization ionization = pointClosestApproach(wireDesign, wireNumber, locHitPos, locHitDir, hit->stepLength());
  double distToWire = ionization.distance;
  if(distToWire > 0.) {
    // Determine on which side of the wire does the particle cross
    //  -1 if particle crosses the wire surface between wireNumber-1 and wireNumber
    //  +1 if particle crosses the wire surface between wireNumber and wireNumber+1
    int adjacent = 1;
    if (ionization.posOnSegment.x() < ionization.posOnWire.x()) {adjacent = -1;}
  // Find the position of the ionization with respect to the adjacent wire
  Ionization ionizationAdj = pointClosestApproach(wireDesign, wireNumber+adjacent, locHitPos, locHitDir, hit->stepLength());
  // Determine if the adjacent wire is closer to the hit on wire surface
  double distToWireAdj = ionizationAdj.distance;
  if ((distToWireAdj > 0.) && (distToWireAdj < distToWire)) {
    distToWire = distToWireAdj;
    wireNumber += adjacent;
    ionization = std::move(ionizationAdj);
  }
  } else {
    ATH_MSG_DEBUG("Failed to get the distance between the wire number = " << wireNumber
                  << " and hit at " <<Amg::toString(hitOnWireSurf, 2)
                  << ". Number of wires = " << numWires
                  << ", "<<m_idHelperSvc->toStringGasGap(wireLayId));
    return {};
  }
 
  // Update the position of ionization on the wire surface
  hitOnWireSurf = ionization.posOnWire;
  ATH_MSG_VERBOSE("Ionization_info: distance: " << ionization.distance
    << " posOnTrack: " <<Amg::toString(ionization.posOnSegment, 3)
    << " posOnWire: " << Amg::toString(ionization.posOnWire, 3)
    << " hit local Pos: " << Amg::toString(locHitPos,2)
    << " hit local Dir: " <<Amg::toString(locHitDir, 2)
    << " EDep: " << hit->energyDeposit() << " EKin: " << hit->kineticEnergy()
    << " pdgId: " << hit->pdgId()
    <<  "gasGap: "<<m_idHelperSvc->toStringGasGap(wireLayId));
 
  // Distance should be in the range [0, 0.9] mm, excepting
  // - particles pass through the wire plane near the edges
  // - secondary particles created inside the gas gap that go through the gas gap partially.
  //   Most of such particles are not muons and have low kinetic energy.
  // - particle with trajectory parallel to the sTGC wire plane
  const double wirePitch = wireDesign.stripPitch();
  if ((distToWire > 0.) && (std::abs(hit->pdgId()) == 13) && (distToWire > (0.5 * wirePitch))) {
    ATH_MSG_DEBUG("Distance to the nearest wire (" << distToWire << ") is greater than expected.");
    ATH_MSG_DEBUG("Hit local Pos: "<<Amg::toString(locHitPos, 2)
                 << " hit local Dir: " <<Amg::toString(locHitDir, 2)
                 << " EDeposited: " << hit->energyDeposit() 
                 << " EKinetic: " << hit->kineticEnergy()
                 << " pdgID: " << hit->pdgId()
                 << " gasGap: "<<m_idHelperSvc->toStringGasGap(wireLayId));
  }
  // Do not digitize hits that are too far from the nearest wire
  if (distToWire > wirePitch) {
    return {};
  }
 
  // Get the gamma pdf parameters associated with the distance of closest approach.
  const GammaParameter gamParam = getGammaParameter(distToWire);
  const double par_kappa = gamParam.kParameter;
  const double par_theta = gamParam.thetaParameter;
  const double most_prob_time = getMostProbableArrivalTime(distToWire);
  // Compute the most probable value of the gamma pdf
  double gamma_mpv = (par_kappa - 1) * par_theta;
  // If the most probable value is less than zero, then set it to zero
  if (gamma_mpv < 0.) {gamma_mpv = 0.;}
  const double t0_par = most_prob_time - gamma_mpv;
 
  // Digit time follows a gamma distribution, so a value val is
  // chosen using a gamma random generator then is shifted by t0
  // to account for drift time.
  // Note: CLHEP::RandGamma takes the parameters k and lambda,
  // where lambda = 1 / theta.
  double digitTime = t0_par + CLHEP::RandGamma::shoot(condContainers.rndEngine, par_kappa, 1/par_theta);
 
  // Sometimes, digitTime is negative because t0_par can be negative.
  // In such case, discard the negative value and shoot RandGamma for another value.
  // However, if that has already been done many times then set digitTime to zero
  // in order to avoid runaway loop.
 
  constexpr int shoot_limit = 4;
  int shoot_counter = 0;
  while (digitTime < 0.) {
    if (shoot_counter > shoot_limit) {
      digitTime = 0.;
      break;
    }
    digitTime = t0_par + CLHEP::RandGamma::shoot(condContainers.rndEngine, par_kappa, 1/par_theta);
    ++shoot_counter;
  }
 
  ATH_MSG_DEBUG("sTgcDigitMaker distance = " << distToWire
    << ", time = " << digitTime
    << ", k parameter = " << par_kappa
    << ", theta parameter = " << par_theta
    << ", most probable time = " << most_prob_time);
 
  if (condContainers.efficiencies) {
    const double efficiency = condContainers.efficiencies->getEfficiency(wireLayId);
    // Lose Hits to match HV efficiency
    if (CLHEP::RandFlat::shoot(condContainers.rndEngine,0.0,1.0) > efficiency) return {};
  }
 
  sTgcDigitVec digits{};
 
  double sDigitTimeWire = digitTime;
  double sDigitTimePad = sDigitTimeWire;
  double sDigitTimeStrip = sDigitTimeWire;
 
  uint16_t bctag = 0;
 
  //##################################################################################
  //######################################### strip readout ##########################
  //##################################################################################
  ATH_MSG_DEBUG("sTgcDigitMaker::strip response ");
  bool isValid = false;
  int channelType = ReadoutChannelType::Strip;
  const Identifier stripLayId = idHelper.channelID(wireLayId, idHelper.multilayer(wireLayId), 
                                                    idHelper.gasGap(wireLayId), channelType, 1, isValid);
 
  const MuonGMR4::StripDesign& stripDesign{readoutElement->stripDesign(stripLayId)};
  // Rotate clockwise 90 in xy plane to get to strip plane orientation, and then extract x,y
  const Amg::Vector2D hitOnStripSurf = (Amg::getRotateZ3D(-90. * Gaudi::Units::deg) * hitOnWireSurf).block<2,1>(0,0);
 
  bool insideBounds = stripDesign.insideTrapezoid(hitOnStripSurf);
  if(!insideBounds) {
    ATH_MSG_DEBUG("Outside of the strip surface boundary : " <<  m_idHelperSvc->toString(stripLayId) << "; local position " <<Amg::toString(hitOnStripSurf, 2));
    return {};
  }
 
  //************************************ find the nearest readout element **************************************
  // Required precision on length in mm
  constexpr double tolerance_length = 0.01*Gaudi::Units::mm;
  int stripNumber = stripDesign.stripNumber(hitOnStripSurf);
  if( stripNumber < 0){
    // Verify if the energy deposit is at the boundary
    const double newPosX = (hitOnStripSurf.x() > 0.0)? hitOnStripSurf.x() - tolerance_length
                                                      : hitOnStripSurf.x() + tolerance_length;
    const Amg::Vector2D newPos(newPosX, hitOnStripSurf.y());
    stripNumber = stripDesign.stripNumber(newPos);
    // Skip hit if still unable to obtain strip number
    if (stripNumber < 0) {
      ATH_MSG_WARNING("Failed to obtain strip number " << m_idHelperSvc->toString(stripLayId) );
      ATH_MSG_WARNING("Position on strip surface = (" << hitOnStripSurf.x() << ", " << hitOnStripSurf.y() << ")");
      return {};
    }
  }
  isValid = false;
  const Identifier stripHitId = idHelper.channelID(stripLayId, idHelper.multilayer(stripLayId), 
                          idHelper.gasGap(stripLayId), channelType, stripNumber, isValid);
  if(!isValid && stripNumber != -1) {
  ATH_MSG_ERROR("Failed to obtain identifier " << m_idHelperSvc->toString(stripHitId) );
  return {};
  }
 
  const int numStrips = stripDesign.numStrips();
  const double stripHalfPitch = 0.5 * stripDesign.stripPitch(); 
 
  //************************************ conversion of energy to charge **************************************
 
  // Typical ionized charge in pC per keV deposited. The constant is determined from ionization
  // study with Garfield program. A note titled "Charge Energy Relation" which outlines
  // conversion can be found here:
  // https://gitlab.cern.ch/carleton-ATLAS/carleton-nsw/stgc-ionization-gain/-/blob/master/Charge_Energy_Relation.pdf?ref_type=heads 
  const double ionized_charge = (5.65E-6)*energyDeposit/CLHEP::keV;
 
  // To get avalanche gain, polya function is taken from Blum paper https://inspirehep.net/literature/807304
  // m_polyaFunction = new TF1("m_polyaFunction","(1.0/[1])*(TMath::Power([0]+1,[0]+1)/TMath::Gamma([0]+1))*TMath::Power(x/[1],[0])*TMath::Exp(-([0]+1)*x/[1])",0,3000000);
 
  // Mean value for total gain due to E field;
  // To calculate this gain from polya distibution, we replace in gamma PDF:
  //     alpha = 1+theta and
  //     beta = 1+theta/mean
  // With these substitutions, gamma PDF gives the same sampling values as those from polya PDF.
  const double gain =  CLHEP::RandGamma::shoot(condContainers.rndEngine, 1. + m_theta, (1. + m_theta)/m_meanGasGain);
 
  // total charge after avalanche
  const double total_charge = gain*ionized_charge;
 
  // Charge Spread including tan(theta) resolution term.
  const double tan_theta = locHitDir.perp()/locHitDir.z();
  // The angle dependance on strip resolution goes as tan^2(angle)
  const double angle_dependency = std::hypot(m_posResIncident, m_posResAngular * tan_theta);
 
  const double cluster_posX = hitOnStripSurf.x();
  double peak_position = CLHEP::RandGaussZiggurat::shoot(condContainers.rndEngine, cluster_posX, m_StripResolution*angle_dependency);
 
  // Each readout plane reads about half the total charge produced on the wire,
  // including a tan(theta) term to describe the increase of charge with incident angle
  const double norm = 0.5 * total_charge;
 
  // Lower limit on strip charge (arbitrary limit), in pC, which has the same units as the parameter ionized_charge. 
  constexpr double tolerance_charge = 0.0005;
  // Set a maximum number of neighbour strips to avoid very long loop. This is an arbitrary number.
  constexpr unsigned int max_neighbor = 10;
 
  // Spread charge on the strips that are on the upper half of the strip cluster
  for (unsigned int iStrip = 0; iStrip <= max_neighbor; ++iStrip) {
    int currentStrip = stripNumber + iStrip;
    if (currentStrip > numStrips) {
      ATH_MSG_VERBOSE("Breaking the upper half strip loop stripNumber: " << currentStrip << " > " << numStrips);
      break;
    }
 
    // Get the strip identifier and create the digit
    isValid = false;
    const Identifier currentStripId = idHelper.channelID(stripLayId, idHelper.multilayer(stripLayId), idHelper.gasGap(stripLayId), 
                              channelType, currentStrip, isValid);
    if (isValid) {
      Amg::Vector2D currentStripPos = readoutElement->localChannelPosition(currentStripId);
      if (currentStripPos == Amg::Vector2D::Zero()) {
          ATH_MSG_VERBOSE("Failed to obtain local position for identifier " << m_idHelperSvc->toString(currentStripId) );
      }
 
      // Estimate the digit charge
      // Position with respect to the peak of the charge curve
      double x_relative = currentStripPos.x() - peak_position;
      // In clusterProfile curve, position should be in the units of strip channel
      double charge = std::hypot(1, m_chargeAngularFactor * tan_theta) * norm * chargeIntegral(x_relative/(2*stripHalfPitch) - 0.5, x_relative/(2*stripHalfPitch) + 0.5);
      // If charge is too small, stop creating neighbor strip
      if (charge < tolerance_charge) {
        ATH_MSG_VERBOSE("Breaking the upper half strip loop stripCharge: " << charge << " < " << tolerance_charge);
        break;
      }
 
      // Estimate digit time
      double strip_time = sDigitTimeStrip;
      // Strip time response can be delayed due to the resistive layer.
      // A correction would be required if the actual VMM front-end doesn't re-align the strip timing.
      if (m_doTimeOffsetStrip) {
        // Determine how far the current strip is from the middle strip
        int indexFromMiddleStrip = iStrip;
        if ((iStrip > 0) && ((x_relative + (0.75 - iStrip * 2) * stripHalfPitch) < 0))
          indexFromMiddleStrip = iStrip - 1;
        // Add time delay due to resistive layer
        strip_time += getTimeOffsetStrip(indexFromMiddleStrip);
      }
 
      addDigit(digits, currentStripId, bctag, strip_time, charge);
 
      ATH_MSG_VERBOSE("Created a strip digit: strip number = " << currentStrip << ", charge = " << charge
                      << ", time = " << strip_time << ", time offset = " << strip_time-sDigitTimeStrip
                      << ", neighbor index = " << iStrip
                      << ", strip position = "<<Amg::toString(currentStripPos, 2));
    }
  }
 
  // The lower half of the strip cluster
  for (unsigned int iStrip = 1; iStrip <= max_neighbor; ++iStrip) {
    int currentStrip = stripNumber - iStrip;
    if (currentStrip < 1) {
      ATH_MSG_VERBOSE("Breaking the lower half strip loop stripNumber: " << currentStrip << " < 1 ");
      break;
    } 
 
    isValid = false;
    const Identifier currentStripId = idHelper.channelID(stripLayId, idHelper.multilayer(stripLayId), idHelper.gasGap(stripLayId), 
                              channelType, currentStrip, isValid);
    if (isValid) {
      Amg::Vector2D currentStripPos = readoutElement->localChannelPosition(currentStripId);
      if (currentStripPos == Amg::Vector2D::Zero()) {
          ATH_MSG_WARNING("Failed to obtain local position for identifier " << m_idHelperSvc->toString(currentStripId) );
      }
 
      // Estimate the digit charge
      double x_relative = currentStripPos.x() - peak_position;
      double charge = std::hypot(1, m_chargeAngularFactor * tan_theta) * norm * chargeIntegral((x_relative/(2*stripHalfPitch) - 0.5), (x_relative/(2*stripHalfPitch) + 0.5));
      if (charge < tolerance_charge) {
        ATH_MSG_VERBOSE("Breaking the lower half strip loop stripCharge: " << charge << " < " << tolerance_charge);
        break;
      }
 
      // Estimate digit time
      double strip_time = sDigitTimeStrip;
      // Time delay due to resistive layer
      if (m_doTimeOffsetStrip) {
        int indexFromMiddleStrip = ((x_relative + (iStrip * 2 - 0.75) * stripHalfPitch) > 0)? iStrip - 1 : iStrip;
        strip_time += getTimeOffsetStrip(indexFromMiddleStrip);
      }
 
      addDigit(digits, currentStripId, bctag, strip_time, charge);
 
      ATH_MSG_VERBOSE("Created a strip digit: strip number = " << currentStrip << ", charge = " << charge
                      << ", time = " << strip_time << ", time offset = " << strip_time-sDigitTimeStrip
                      << ", neighbor index = " << iStrip
                      << ", strip position = " << Amg::toString(currentStripPos, 2));
    }
  }
  // end of strip digitization
 
  if(m_digitMode == digitMode::StripsOnly) {
    ATH_MSG_WARNING("Only digitize strip response !");
    return digits;
  }
 
  //##################################################################################
  //######################################### pad readout ##########################
  //##################################################################################
  ATH_MSG_DEBUG("sTgcDigitMaker::pad response ");
  channelType = ReadoutChannelType::Pad;
 
  Identifier padLayId = idHelper.channelID(wireLayId, idHelper.multilayer(wireLayId), 
                                        idHelper.gasGap(wireLayId), channelType, 1);
  
  const MuonGMR4::PadDesign& padDesign{readoutElement->padDesign(padLayId)};
  //************************************ find the nearest readout element **************************************
  
  // Since the orientation of pad layer is the same as wirelayer, the x,y can be taken directly
  // from the hit position on wire surface
  const Amg::Vector2D hitOnPadSurf = hitOnWireSurf.block<2,1>(0,0);
 
  
  insideBounds = padDesign.insideTrapezoid(hitOnPadSurf);
 
  if(insideBounds) {
    std::pair<uint, uint> padEtaPhi = padDesign.channelNumber(hitOnPadSurf);
    unsigned int padEta = padEtaPhi.first;
    unsigned int padPhi = padEtaPhi.second;
    if( padEta < 1 || padPhi < 1){
      // Verify if the energy deposit is at the boundary
      const double newPosX = (hitOnPadSurf.x()>0.0)? hitOnPadSurf.x() - tolerance_length
                                                    : hitOnPadSurf.x() + tolerance_length;
      const double newPosY = (hitOnPadSurf.y()>0.0)? hitOnPadSurf.y() - tolerance_length
                                                    : hitOnPadSurf.y() + tolerance_length;
      const Amg::Vector2D newPosition(newPosX, newPosY);
      padEtaPhi = padDesign.channelNumber(newPosition);
      padEta = padEtaPhi.first;
      padPhi = padEtaPhi.second;
      // Skip hit if still unable to obtain pad number
      if( padEta < 1 || padPhi < 1){
        ATH_MSG_WARNING("Failed to obtain pad number for " << m_idHelperSvc->toString(padLayId) );
        ATH_MSG_WARNING("Position on pad surface = (" << hitOnPadSurf.x() << ", " << hitOnPadSurf.y() << ")");
        return digits;
      }
    }
    isValid = false;
    const Identifier padHitId = idHelper.padID(padLayId, idHelper.multilayer(padLayId), 
                               idHelper.gasGap(padLayId), channelType, padEta, padPhi, isValid);
    if(isValid) {
      // Find centre position of pad
      Amg::Vector2D padPos = readoutElement->localChannelPosition(padHitId);
      if (padPos == Amg::Vector2D::Zero()) {
        ATH_MSG_ERROR("Failed to obtain local position for neighbor pad with identifier " << m_idHelperSvc->toString(padHitId) );
        return digits;
      }
 
      // Pad sharing needs to look at position on hit vs pad boundaries
      const Amg::Vector2D diff = hitOnPadSurf - padPos;
      double halfPadHeight = 0.5 * readoutElement->padHeight(padHitId);
      // Extracting pad corners to get the width of the pad at hit position
      std::array<Amg::Vector2D, 4> padHitCorners = padDesign.padCorners(padEtaPhi);
      double padBottomBase = (padHitCorners[0] - padHitCorners[1]).norm();
      double padTopBase = (padHitCorners[2] - padHitCorners[3]).norm();
      // Linearly interpolating pad width (as done by channelWidth function in R3)
      // Interpolating pad width at the y-position of the hit (diff.y()) using pad top and bottom bases.
      // The expression below computes the width at a relative vertical position inside the pad.
      // Specifically: (1 + diff.y() / halfPadHeight) * 0.25 = (y_from_bottom / padHeight)
      double halfPadWidth = 0.5 * padBottomBase + 0.25 * (1 + diff.y()/halfPadHeight) * (padTopBase - padBottomBase);
 
      // Charge sharing happens within 4mm window for pads
      // i.e. the charge is spread over a 2D gaussian of total radius about 4mm
      // This value depends on actual width of the distribution, but we don't have the
      // width of the charge distribution for pads. So lets consider 2.5*sigma, where
      // sigma is the width of the charge distribution for strips.
      double deltaX = halfPadWidth - std::abs(diff.x());
      double deltaY = halfPadHeight - std::abs(diff.y());
      bool isNeighX = deltaX < 2.5*m_clusterParams[1];
      bool isNeighY = deltaY < 2.5*m_clusterParams[1];
      // Pad width can be calculated to be very slightly larger than it should due to rounding errors
      // So if a hit falls on a given pad but is "outside" the width, just define it to be on the boundary of 2 pads.
      if (deltaX < 0.) deltaX = 0.1;
      if (deltaY < 0.) deltaY = 0.1;
 
      if (m_doPadSharing && (isNeighX || isNeighY)){
        // Phi == 1 at the right in the local geometry
        // In local coordinates, the pad to the right has phi' = phi - 1
        unsigned int newPhi = diff.x() > 0 ? padEtaPhi.second-1 : padEtaPhi.second+1;
        unsigned int newEta = diff.y() > 0 ? padEtaPhi.first+1 : padEtaPhi.first-1;
        bool validEta = newEta > 0 && newEta < readoutElement->numPadEta(padHitId) + 1;
        bool validPhi = newPhi > 0 && newPhi < readoutElement->numPadPhi(padHitId) + 1;
    
        if (isNeighX && isNeighY && validEta && validPhi){
          // 4 pads total; makes life a bit harder. Corner of 4 valid pads
          isValid = false;
          const Identifier neigh_ID_X = idHelper.padID(padHitId, idHelper.multilayer(padHitId), idHelper.gasGap(padHitId),
                                                  channelType, padEta, newPhi, isValid);
          isValid = false;
          const Identifier neigh_ID_Y = idHelper.padID(padHitId, idHelper.multilayer(padHitId), idHelper.gasGap(padHitId),
                                                  channelType, newEta, padPhi, isValid);
          isValid = false;
          const Identifier neigh_ID_XY = idHelper.padID(padHitId, idHelper.multilayer(padHitId), idHelper.gasGap(padHitId),
                                                  channelType, newEta, newPhi, isValid);
          double xQfraction = getPadChargeFraction(deltaX);
          double yQfraction = getPadChargeFraction(deltaY);
 
          // Main pad gets 1 - Qfraction of total
          addDigit(digits, neigh_ID_X, bctag, sDigitTimePad, xQfraction*(1.-yQfraction)*0.5*total_charge);
          addDigit(digits, neigh_ID_Y, bctag, sDigitTimePad, yQfraction*(1.-xQfraction)*0.5*total_charge);
          addDigit(digits, neigh_ID_XY, bctag, sDigitTimePad, xQfraction*yQfraction*0.5*total_charge);
          addDigit(digits, padHitId, bctag, sDigitTimePad, (1.-xQfraction-yQfraction+xQfraction*yQfraction)*0.5*total_charge);
        } else if (isNeighX && validPhi){
          // There is only 1 neighbor, immediately to the left or right.
          isValid = false;
          const Identifier neigh_ID = idHelper.padID(padHitId, idHelper.multilayer(padHitId), idHelper.gasGap(padHitId), 
                                               channelType, padEta, newPhi, isValid);
 
          // NeighborPad gets Qfraction of the total pad charge: 0.5*total_charge
          // Main pad gets 1 - Qfraction of total
          double xQfraction = getPadChargeFraction(deltaX);
          addDigit(digits, padHitId, bctag, sDigitTimePad, (1.-xQfraction)*0.5*total_charge);
          addDigit(digits, neigh_ID, bctag, sDigitTimePad, xQfraction*0.5*total_charge);
        }
        else if (isNeighY && validEta){
          // There is only 1 neighbor, immediately above or below
          isValid = false;
          const Identifier neigh_ID = idHelper.padID(padHitId, idHelper.multilayer(padHitId), idHelper.gasGap(padHitId), 
                                               channelType, newEta, padPhi, isValid);
 
          // NeighborPad gets Qfraction of the total pad charge: 0.5*total_charge
          // Main pad gets 1 - Qfraction of total
          double yQfraction = getPadChargeFraction(deltaX);
          addDigit(digits, padHitId, bctag, sDigitTimePad, (1.-yQfraction)*0.5*total_charge);
          addDigit(digits, neigh_ID, bctag, sDigitTimePad, yQfraction*0.5*total_charge);
        }
 
      }
      else{
        // No charge sharing: hit is nicely isolated within pad
        addDigit(digits, padHitId, bctag, sDigitTimePad, 0.5*total_charge);
      }
 
    }
    else if(padEta != 0 || padPhi != 0) {
      ATH_MSG_ERROR("Failed to obtain identifier " << m_idHelperSvc->toString(padHitId) );
    }
  }
  else {
    ATH_MSG_DEBUG("Outside of the pad surface boundary :" << m_idHelperSvc->toString(padLayId)<< " local position " <<Amg::toString(hitOnPadSurf, 2));
  }
 
  if(m_digitMode == digitMode::StripsAndPads) {
    ATH_MSG_WARNING("Only digitize strip/pad response !");
    return digits;
  }
 
  //##################################################################################
  //######################################### wire readout ##########################
  //##################################################################################
  ATH_MSG_DEBUG("sTgcDigitMaker::wire response ");
  channelType = ReadoutChannelType::Wire;
 
  // wireLayId already defined at the start
  //************************************ find the nearest readout element **************************************
  insideBounds = wireDesign.insideTrapezoid(hitOnWireSurf.block<2,1>(0,0));
 
  if(insideBounds) {
      // Determine the wire number
      int wiregroupNumber = (wireDesign.wireNumber(hitOnWireSurf.block<2,1>(0,0))).first;
      if( wiregroupNumber < 1){
        // Verify if the energy deposit is at the boundary
        const double newPosX = (hitOnWireSurf.x() > 0.0)? hitOnWireSurf.x() - tolerance_length
                                                          : hitOnWireSurf.x() + tolerance_length;
        const Amg::Vector2D newPos(newPosX, hitOnWireSurf.y());
        wiregroupNumber = (wireDesign.wireNumber(newPos)).first;
        // Skip hit if still unable to obtain wire number
        if (wiregroupNumber < 1) {
          ATH_MSG_WARNING("Failed to obtain wire number " << m_idHelperSvc->toString(wireLayId) );
          ATH_MSG_WARNING("Position on wire surface = (" << hitOnWireSurf.x() << ", " << hitOnWireSurf.y() << ")");
          return digits;
        }
      }
 
      // Find ID of the actual wiregroup
      isValid = false;
      const Identifier wireGroupHitId = idHelper.channelID(wireLayId, idHelper.multilayer(wireLayId), idHelper.gasGap(wireLayId), 
                                  channelType, wiregroupNumber, isValid);
 
      if(isValid) {
        int numWireGroups = readoutElement->numChannels(wireLayId);
        if(wiregroupNumber >=1 && wiregroupNumber <= numWireGroups) {
          addDigit(digits, wireGroupHitId, bctag, sDigitTimeWire, total_charge);
        }
      } // end of if(isValid)
      else if (wiregroupNumber != -1){
        ATH_MSG_ERROR("Failed to obtain wiregroup identifier " << m_idHelperSvc->toString(wireGroupHitId));
      }
  }
  else {
    ATH_MSG_DEBUG("Outside of the wire surface boundary :" << m_idHelperSvc->toString(wireLayId)<< " local position " <<Amg::toString(hitOnWireSurf, 2));
  }
  // end of wire digitization
 
  return digits;
}

getGammaParameter() [1/2]

sTgcDigitMaker::GammaParameter sTgcDigitMaker::getGammaParameter ( double  distance ) const

private

Find the gamma pdf parameters of a given distance.

Definition at line 790 of file MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx.

                                                                                    {
 
  const double d = std::abs(distance); 
  // Find the parameters assuming the container is sorted in ascending order of 'lowEdge'
  int index{-1};
  for (const auto& par: m_gammaParameter) {
    if (d < par.lowEdge) {
      break;
    }
    ++index;
  }
  return m_gammaParameter.at(index);
}

getGammaParameter() [2/2]

GammaParameter sTgcDigitMaker::getGammaParameter ( double  distance ) const

private

Retrieves gamma distribution parameters based on distance.

getMostProbableArrivalTime() [1/2]

double sTgcDigitMaker::getMostProbableArrivalTime ( double  distance ) const

private

Get the most probable time of arrival.

Definition at line 805 of file MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx.

                                                                       {
 
  const double d{std::abs(distance)};
  double mpt{0};
  for (size_t t = 0 ; t < m_mostProbableArrivalTime.size(); ++t){
    mpt += m_mostProbableArrivalTime[t] * std::pow(d, t);
  }
  return mpt;
}

getMostProbableArrivalTime() [2/2]

double sTgcDigitMaker::getMostProbableArrivalTime ( double  distance ) const

private

Computes the most probable arrival time based on the distance of closest approach.

getPadChargeFraction() [1/2]

double sTgcDigitMaker::getPadChargeFraction ( double  distance )

staticprivate

Definition at line 860 of file MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx.

                                                           {
  // The charge fraction that is found past a distance x away from the
  // centre of a 2D gaussian distribution of width of cluster profile is
  // described by a modified error function.
 
  // The modified error function perfectly describes
  // the pad charge sharing distribution figure 16 of the sTGC
  // testbeam paper https://arxiv.org/pdf/1509.06329.pdf
 
  return 0.5 * (1.0 - std::erf( distance / (std::sqrt(2) * m_clusterProfile[1])));
}

getPadChargeFraction() [2/2]

static double sTgcDigitMaker::getPadChargeFraction ( double  distance )

staticprivate

Computes charge fraction shared among pads.

getTimeOffsetStrip() [1/2]

double sTgcDigitMaker::getTimeOffsetStrip ( size_t  neighbor_index ) const

private

Get digit time offset of a strip depending on its relative position to the strip at the centre of the cluster.

It returns 0 ns by default, as well as when it fails or container is empty.

Definition at line 855 of file MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx.

                                                                     {
    return m_timeOffsetStrip.at(std::min(neighbor_index, m_timeOffsetStrip.size() -1));
}

getTimeOffsetStrip() [2/2]

double sTgcDigitMaker::getTimeOffsetStrip ( size_t  neighbor_index ) const

private

Gets the time offset for a strip cluster.

Get digit time offset of a strip depending on its relative position to the strip at the centre of the cluster. It returns 0 ns by default, as well as when it fails or container is empty.

initialize() [1/2]

StatusCode sTgcDigitMaker::initialize

(

)

Initializes sTgcHitIdHelper and sTgcIdHelper, and call the functions to read files containing digitization parameters.

Definition at line 49 of file MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx.

                                      {
 
  // Read share/sTGC_Digitization_timeArrivale.dat, containing the digit time of arrival
  ATH_CHECK(readFileOfTimeArrival());
 
  // Read share/sTGC_Digitization_timeOffsetStrip.dat if the the strip time correction is enable
  if (m_doTimeOffsetStrip) {
    ATH_CHECK(readFileOfTimeOffsetStrip());
  }
  // check the digitization channel type
  if(m_channelTypes!=1 && m_channelTypes!=2 && m_channelTypes!=3) {
      ATH_MSG_ERROR("Invalid ChannelTypes : " << m_channelTypes << " (valid values are : 1 --> strips ; 2 --> strips / wires ; 3 --> strips / wires / pads)");
      return StatusCode::FAILURE;
  }
  return StatusCode::SUCCESS;
}

initialize() [2/2]

StatusCode sTgcDigitMaker::initialize

(

)

Initialize digitization parameters, including reading necessary data files.

initMessaging()

void AthMessaging::initMessaging ( ) const

privateinherited

Initialize our message level and MessageSvc.

This method should only be called once.

Definition at line 39 of file AthMessaging.cxx.

{
  m_imsg = Athena::getMessageSvc();
  m_lvl = m_imsg ?
    static_cast<MSG::Level>( m_imsg.load()->outputLevel(m_nm) ) :
    MSG::INFO;
}

msg() [1/2]

MsgStream & AthMessaging::msg ( ) const

inlineinherited

The standard message stream.

Returns a reference to the default message stream May not be invoked before sysInitialize() has been invoked.

Definition at line 164 of file AthMessaging.h.

{
  MsgStream* ms = m_msg_tls.get();
  if (!ms) {
    if (!m_initialized.test_and_set()) initMessaging();
    ms = new MsgStream(m_imsg,m_nm);
    m_msg_tls.reset( ms );
  }
 
  ms->setLevel (m_lvl);
  return *ms;
}

msg() [2/2]

MsgStream & AthMessaging::msg ( const MSG::Level  lvl ) const

inlineinherited

The standard message stream.

Returns a reference to the default message stream May not be invoked before sysInitialize() has been invoked.

Definition at line 179 of file AthMessaging.h.

{ return msg() << lvl; }

msgLvl()

bool AthMessaging::msgLvl ( const MSG::Level  lvl ) const

inlineinherited

Test the output level.

Parameters

lvl	The message level to test against

Returns

boolean Indicating if messages at given level will be printed

Return values

true	Messages at level "lvl" will be printed

Definition at line 151 of file AthMessaging.h.

{
  if (!m_initialized.test_and_set()) initMessaging();
  if (m_lvl <= lvl) {
    msg() << lvl;
    return true;
  } else {
    return false;
  }
}

pointClosestApproach() [1/2]

sTgcDigitMaker::Ionization sTgcDigitMaker::pointClosestApproach ( const MuonGM::sTgcReadoutElement *  readoutEle, const Identifier &  id, int  wireNumber, const Amg::Vector3D &  preStepPos, const Amg::Vector3D &  postStepPos  ) const

private

Determine the points where the distance between two segments is smallest.

Given two segments, e.g. a particle trajectory and a sTGC wire, solve for the two points, the point on the trajectory and the point on the wire, where the distance between the two segments is the smallest.

Positions returned are in the local coordinate frame of the wire plane. Returns an object with distance of -9.99 in case of error.

Definition at line 637 of file MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx.

                                                                                                      {
  // tolerance in mm
  constexpr double tolerance = 0.001;
  // minimum segment length
  constexpr double min_length = 0.1;
 
  // Position of the ionization
  Ionization ionization;
 
  // Wire number should be one or greater
  if (wireNumber < 1) return ionization;
 
  // Wire number too large
  if (wireNumber > detEl->numberOfWires(id)) return ionization;
 
  // Finding smallest distance and the points at the smallest distance.
  //  The smallest distance between two lines is perpendicular to both lines.
  //  We can construct two lines in the wire surface local coordinate frame:
  //  - one for the hit segment with equation h0 + t * v_h, where h0 is a point
  //    and v_h is the unit vector of the hit segment
  //  - another for the wire with similar equation w0 + s * v_w, where w0 is a
  //    point and v_w is the unit vector of the wire line
  //  Then it is possible to determine the factors s and t by requiring the
  //  dot product to be zero:
  //   1. (h0 - w0) \dot v_h = 0
  //   2. (h0 - w0) \dot v_w = 0
 
  // Wire pitch
  const double wire_pitch = detEl->wirePitch();
  // Wire local position on the wire plane, the y-coordinate is arbitrary and z-coordinate is zero
  const double wire_posX = detEl->positionFirstWire(id) + (wireNumber - 1) * wire_pitch;
  const Amg::Vector3D wire_position(wire_posX, postStepPos.y(), 0.);
  // The wires are parallel to Y in the wire plane's coordinate frame
  const Amg::Vector3D wire_direction{Amg::Vector3D::UnitY()};
 
  // particle trajectory
  Amg::Vector3D hit_direction{postStepPos - preStepPos};
  const double seg_length = hit_direction.mag();
  if (seg_length > tolerance) hit_direction /= seg_length;
 
  // Find the point on the track segment that is closest to the wire
  if (seg_length < min_length) {
    ionization.posOnSegment = postStepPos;
    ionization.posOnWire = wire_position;
    ionization.distance = std::hypot(postStepPos.x() - wire_posX, postStepPos.z());
    return ionization;
  }
 
  // Dot product between the wire and hit direction
  const double cos_theta = wire_direction.dot(hit_direction);
  // distance between the hit and wire
  const Amg::Vector3D dist_wire_hit = postStepPos - wire_position;
 
  // Verifier the special case where the two lines are parallel
  if (std::abs(cos_theta - 1.0) < tolerance) {
    ATH_MSG_DEBUG("The track segment is parallel to the wire, position of digit is undefined");
    ionization.posOnSegment = 0.5 * (postStepPos + preStepPos);
    ionization.posOnWire = Amg::Vector3D(wire_posX, ionization.posOnSegment.y(), 0.0);
    ionization.distance = std::hypot(ionization.posOnSegment.x() - wire_posX,
                                     ionization.posOnSegment.z());
    return ionization;
  }
 
  // Perpendicular component squared
  const double sin_theta_2 = 1.0 - cos_theta * cos_theta;
 
  //* Point on the hit segment
  const double factor_hit = (-dist_wire_hit.dot(hit_direction)
                           + dist_wire_hit.dot(wire_direction) * cos_theta) / sin_theta_2;
  Amg::Vector3D ionization_pos = postStepPos + factor_hit * hit_direction;
 
  // If the point is on the track segment, then compute the other point on the wire.
  // Otherwise, set the ionization at the pre-step position and compute where it
  // should be on the wire.
  Amg::Vector3D pos_on_wire{Amg::Vector3D::Zero()};
  if (factor_hit < seg_length) {
    //* Point on the wire line
    const double factor_wire = (dist_wire_hit.dot(wire_direction)
                              - dist_wire_hit.dot(hit_direction) * cos_theta) / sin_theta_2;
    pos_on_wire = wire_position + factor_wire * wire_direction;
  } else {
    ionization_pos = preStepPos;
    pos_on_wire = wire_position - (seg_length * cos_theta) * wire_direction;
  }
 
  // Save the distance and ionization position
  ionization.distance = (ionization_pos - pos_on_wire).mag();
  ionization.posOnSegment = ionization_pos;
  ionization.posOnWire = pos_on_wire;
 
  return ionization;
}

pointClosestApproach() [2/2]

sTgcDigitMaker::Ionization sTgcDigitMaker::pointClosestApproach ( const MuonGMR4::WireGroupDesign &  wireDesign, int  wireNumber, const Amg::Vector3D &  locHitPos, const Amg::Vector3D &  locHitDir, const double  stepLength  ) const

private

Computes the closest approach between a trajectory and a wire segment.

Given two segments, e.g. a particle trajectory and a sTGC wire, solve for the two points, the point on the trajectory and the point on the wire, where the distance between the two segments is the smallest.

Positions returned are in the local coordinate frame of the wire plane. Returns an object with distance of -9.99 in case of error.

Definition at line 613 of file MuonPhaseII/MuonDigitization/sTgcDigitizationR4/src/sTgcDigitMaker.cxx.

                                                                                               {
 
  constexpr double angular_tolerance = 1e-3;
  // Position of the ionization
  Ionization ionization;
  // Finding smallest distance and the points at the smallest distance.
  //  The smallest distance between two lines is perpendicular to both lines.
  // Previous logic was to find the perpendicular distance between two lines using projection geometry
  //  We can construct two lines in the wire surface local coordinate frame:
  //  - one for the hit segment with equation h0 + t * v_h, where h0 is a point
  //    and v_h is the unit vector of the hit segment
  //  - another for the wire with similar equation w0 + s * v_w, where w0 is a
  //    point and v_w is the unit vector of the wire line
  //  Then it is possible to determine the closest points on each line
  //  by requiring that the vector between them is perpendicular to both:
  //   1. (h0 + t*v_h - w0 - s*v_w) · v_h = 0
  //   2. (h0 + t*v_h - w0 - s*v_w) · v_w = 0
 
  // We have replaced this logic with using Amg::intersect<3> method
 
  // Geometry setup
  const double wirePitch = wireDesign.stripPitch();
  const double wirePosX = wireDesign.firstStripPos().x() + (wireNumber - 1) * wirePitch;
  const Amg::Vector3D wireDir(wireDesign.stripDir().x(), wireDesign.stripDir().y(), 0.0);
  const Amg::Vector3D wirePos(wirePosX, locHitPos.y(), 0.);
 
  // Use Amg::intersect to find closest point on hit segment to wire plane
  std::optional<double> scaleHit = Amg::intersect<3>(locHitPos, locHitDir, Amg::Vector3D::UnitZ(), 0);
  if (!scaleHit || std::abs(std::abs(wireDir.dot(locHitDir)) - 1.0) < angular_tolerance) {
    ATH_MSG_DEBUG("The track segment is parallel to the wire, position of digit is undefined");
    ionization.posOnSegment = locHitPos;
    ionization.posOnWire = wirePos;
    ionization.distance = std::hypot(locHitPos.x() - wirePosX, locHitPos.z());
    return ionization;
  }
  // Position on hit segment
  Amg::Vector3D ionizationPos = locHitPos + scaleHit.value() * locHitDir;
 
  if (scaleHit.value() > stepLength) {
    ionization.posOnSegment = locHitPos;
    const Amg::Vector3D closestPointToWirePlane = locHitPos + stepLength * locHitDir;
    ionization.posOnWire = wirePos;
    ionization.distance = std::abs(closestPointToWirePlane.z());
    return ionization;
  }
 
  // Project ionization position onto the wire line
  const double scaleWire = (ionizationPos - wirePos).dot(wireDir);
  const Amg::Vector3D closestPointOnWire = wirePos + scaleWire * wireDir;
 
  // Fill ionization result
  ionization.posOnSegment = ionizationPos;
  ionization.posOnWire = closestPointOnWire;
  ionization.distance = (ionizationPos - closestPointOnWire).mag();
 
  return ionization;
}

readFileOfTimeArrival() [1/2]

StatusCode sTgcDigitMaker::readFileOfTimeArrival ( )

private

Read share/sTGC_Digitization_timeArrival.dat.

Definition at line 750 of file MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx.

                                                 {
  // Verify the file sTGC_Digitization_timeArrival.dat exists
  const std::string file_name = "sTGC_Digitization_timeArrival.dat";
  std::string file_path = PathResolver::find_file(file_name, "DATAPATH");
  if(file_path.empty()) {
    ATH_MSG_FATAL("readFileOfTimeWindowOffset(): Could not find file " << file_name );
    return StatusCode::FAILURE;
  }
 
  // Open the sTGC_Digitization_timeArrival.dat file
  std::ifstream ifs{file_path, std::ios::in};
  if(ifs.bad()) {
    ATH_MSG_FATAL("sTgcDigitMaker: Failed to open time of arrival file " << file_name );
    return StatusCode::FAILURE;
  }
 
  // Read the sTGC_Digitization_timeWindowOffset.dat file
  std::string line;
  GammaParameter param{};
 
  while (std::getline(ifs, line)) {
    std::string key;
    std::istringstream iss(line);
    iss >> key;
    if (key.compare("bin") == 0) {
      iss >> param.lowEdge >> param.kParameter >> param.thetaParameter;
      m_gammaParameter.push_back(param);
    } else if (key.compare("mpv") == 0)  {
      double mpt;
      int idx{0};
      while (iss >> mpt) {m_mostProbableArrivalTime[idx++] = mpt;}
    }
  }
 
  // Close the file
  ifs.close();
  return StatusCode::SUCCESS;
}

readFileOfTimeArrival() [2/2]

StatusCode sTgcDigitMaker::readFileOfTimeArrival ( )

private

Reads time arrival data file.

readFileOfTimeOffsetStrip() [1/2]

StatusCode sTgcDigitMaker::readFileOfTimeOffsetStrip ( )

private

Read share/sTGC_Digitization_timeOffsetStrip.dat.

Definition at line 816 of file MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx.

                                                     {
  // Verify the file sTGC_Digitization_timeOffsetStrip.dat exists
  const std::string file_name = "sTGC_Digitization_timeOffsetStrip.dat";
  std::string file_path = PathResolver::find_file(file_name, "DATAPATH");
  if(file_path.empty()) {
    ATH_MSG_FATAL("readFileOfTimeWindowOffset(): Could not find file " << file_name );
    return StatusCode::FAILURE;
  }
 
  // Open the sTGC_Digitization_timeOffsetStrip.dat file
  std::ifstream ifs{file_path, std::ios::in};
  if(ifs.bad()) {
    ATH_MSG_FATAL("Failed to open time of arrival file " << file_name );
    return StatusCode::FAILURE;
  }
 
  // Initialize the container to store the time offset.
  // The number of parameters, 6, corresponds to the number of lines to be read
  // from sTGC_Digitization_timeOffsetStrip.dat.
  // Setting the default offset to 0 ns.
 
  // Read the input file
  std::string line;
  size_t index{0};
  double value{0.0};
  while (std::getline(ifs, line)) {
    std::string key;
    std::istringstream iss(line);
    iss >> key;
    if (key.compare("strip") == 0) {
      iss >> index >> value;
      if (index >= m_timeOffsetStrip.size()) continue;
      m_timeOffsetStrip.at(index) = value;
    }
  }
  return StatusCode::SUCCESS;
}

readFileOfTimeOffsetStrip() [2/2]

StatusCode sTgcDigitMaker::readFileOfTimeOffsetStrip ( )

private

Reads strip time offset data file.

setLevel()

void AthMessaging::setLevel ( MSG::Level  lvl )

inherited

Change the current logging level.

Use this rather than msg().setLevel() for proper operation with MT.

Definition at line 28 of file AthMessaging.cxx.

{
  m_lvl = lvl;
}

Member Data Documentation

ATLAS_THREAD_SAFE

std::atomic_flag m_initialized AthMessaging::ATLAS_THREAD_SAFE = ATOMIC_FLAG_INIT

mutableprivateinherited

Messaging initialized (initMessaging)

Definition at line 141 of file AthMessaging.h.

m_applyAsBuiltBLines

bool sTgcDigitMaker::m_applyAsBuiltBLines {false}

private

Definition at line 165 of file MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h.

m_channelTypes

int sTgcDigitMaker::m_channelTypes {3}

private

define offsets and widths of time windows for signals from wiregangs and strips.

The offsets are defined as relative time diffference with respect to the time after TOF and cable length corrections. Bunch crossing time is specified.

Definition at line 147 of file MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h.

m_chargeAngularFactor

double sTgcDigitMaker::m_chargeAngularFactor {4.0}

private

Definition at line 163 of file MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h.

m_clusterParams

constexpr std::array<double, 2> sTgcDigitMaker::m_clusterParams {0.573, 1.092}

staticconstexprprivate

Definition at line 175 of file MuonPhaseII/MuonDigitization/sTgcDigitizationR4/sTgcDigitizationR4/sTgcDigitMaker.h.

m_clusterProfile

constexpr std::array<double, 2> sTgcDigitMaker::m_clusterProfile {0.573, 1.092}

staticconstexprprivate

Definition at line 161 of file MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h.

m_digitMode

digitMode sTgcDigitMaker::m_digitMode = AllChType

private

define offsets and widths of time windows for signals from wiregroups and strips.

The offsets are defined as relative time diffference with respect to the time after TOF and cable length corrections. Bunch crossing time is specified.

Definition at line 185 of file MuonPhaseII/MuonDigitization/sTgcDigitizationR4/sTgcDigitizationR4/sTgcDigitMaker.h.

m_doPadSharing

bool sTgcDigitMaker::m_doPadSharing {false}

private

Definition at line 150 of file MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h.

m_doTimeOffsetStrip

bool sTgcDigitMaker::m_doTimeOffsetStrip {false}

private

Definition at line 154 of file MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h.

m_gammaParameter

std::vector< GammaParameter > sTgcDigitMaker::m_gammaParameter

private

Definition at line 131 of file MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h.

m_idHelperSvc

const Muon::IMuonIdHelperSvc * sTgcDigitMaker::m_idHelperSvc {nullptr}

private

Definition at line 138 of file MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h.

m_imsg

std::atomic<IMessageSvc*> AthMessaging::m_imsg { nullptr }

mutableprivateinherited

MessageSvc pointer.

Definition at line 135 of file AthMessaging.h.

m_lvl

std::atomic<MSG::Level> AthMessaging::m_lvl { MSG::NIL }

mutableprivateinherited

Current logging level.

Definition at line 138 of file AthMessaging.h.

m_meanGasGain

double sTgcDigitMaker::m_meanGasGain {5.e4}

private

Definition at line 149 of file MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h.

m_mostProbableArrivalTime

std::array< double, 5 > sTgcDigitMaker::m_mostProbableArrivalTime {make_array<double, 5>(0.)}

private

Definition at line 133 of file MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h.

m_msg_tls

boost::thread_specific_ptr<MsgStream> AthMessaging::m_msg_tls

mutableprivateinherited

MsgStream instance (a std::cout like with print-out levels)

Definition at line 132 of file AthMessaging.h.

m_nm

std::string AthMessaging::m_nm

privateinherited

Message source name.

Definition at line 129 of file AthMessaging.h.

m_posResAngular

double sTgcDigitMaker::m_posResAngular {0.305/m_StripResolution}

private

Definition at line 158 of file MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h.

m_posResIncident

double sTgcDigitMaker::m_posResIncident {1.}

private

Definition at line 157 of file MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h.

m_StripResolution

double sTgcDigitMaker::m_StripResolution {0.0949}

private

Definition at line 156 of file MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h.

m_theta

double sTgcDigitMaker::m_theta {10}

private

Definition at line 148 of file MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h.

m_timeOffsetStrip

std::array< double, 6 > sTgcDigitMaker::m_timeOffsetStrip {make_array<double, 6>(0.)}

private

Definition at line 136 of file MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h.

---

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

• MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h
• MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx