Belle II Software development
TrackExtrapolateG4e Class Reference

geant4e-based track extrapolation. More...

#include <TrackExtrapolateG4e.h>

Public Member Functions

 ~TrackExtrapolateG4e ()
 destructor
 
void initialize (double minPt, double minKE, std::vector< Const::ChargedStable > &hypotheses)
 Initialize for track extrapolation by the EXT module.
 
void initialize (double meanDt, double maxDt, double maxSeparation, double maxKLMTrackClusterDistance, double maxECLTrackClusterDistance, double minPt, double minKE, bool addHitsToRecoTrack, std::vector< Const::ChargedStable > &hypotheses)
 Initialize for track extrapolation by the MUID module.
 
void beginRun (bool flag)
 Perform beginning-of-run actions.
 
void event (bool flag)
 Performs track extrapolation for all tracks in one event.
 
void endRun (bool flag)
 Perform end-of-run actions.
 
void terminate (bool flag)
 Terminates this singleton.
 
void extrapolate (int pdgCode, double tof, const G4ThreeVector &position, const G4ThreeVector &momentum, const G4ErrorSymMatrix &covariance)
 Performs track extrapolation for a single track (specified in genfit2 units).
 

Static Public Member Functions

static TrackExtrapolateG4egetInstance ()
 Get the singleton's address.
 

Private Member Functions

 TrackExtrapolateG4e ()
 constructor is hidden; user calls TrackExtrapolateG4e::getInstance() instead
 
 TrackExtrapolateG4e (TrackExtrapolateG4e &)
 copy constructor is hidden; user calls TrackExtrapolateG4e::getInstance() instead
 
void swim (ExtState &, G4ErrorFreeTrajState &, const std::vector< std::pair< ECLCluster *, G4ThreeVector > > *, const std::vector< std::pair< KLMCluster *, G4ThreeVector > > *, std::vector< std::map< const Track *, double > > *)
 Swim a single track (MUID) until it stops or leaves the target cylinder.
 
void swim (ExtState &, G4ErrorFreeTrajState &)
 Swim a single track (EXT) until it stops or leaves the target cylinder.
 
void findClosestTrackToKLMClusters ()
 Find the closest Track to each KLMCluster and fill the corresponding fields of KLMCluster objects.
 
void registerVolumes ()
 Register the list of geant4 physical volumes whose entry/exit points will be saved during extrapolation.
 
void getVolumeID (const G4TouchableHandle &, Const::EDetector &, int &)
 Get the physical volume information for a geant4 physical volume.
 
void fromG4eToPhasespace (const G4ErrorFreeTrajState &, G4ErrorSymMatrix &)
 Convert the geant4e 5x5 covariance to phasespace 6x6 covariance.
 
void fromPhasespaceToG4e (const G4ThreeVector &, const G4ErrorSymMatrix &, G4ErrorTrajErr &)
 Convert the phasespace covariance to geant4e covariance.
 
void fromPhasespaceToG4e (const TVector3 &, const TMatrixDSym &, G4ErrorTrajErr &)
 Convert the phasespace covariance to geant4e covariance.
 
ExtState getStartPoint (const Track &, int, G4ErrorFreeTrajState &)
 Get the start point for a new reconstructed track with specific PDG hypothesis.
 
void createExtHit (const ExtHitStatus, const ExtState &, const G4ErrorFreeTrajState &, const G4StepPoint *, const G4TouchableHandle &)
 Create another EXT extrapolation hit for a track candidate.
 
void createECLHit (const ExtState &, const G4ErrorFreeTrajState &, const G4StepPoint *, const G4StepPoint *, const G4TouchableHandle &, const std::pair< ECLCluster *, G4ThreeVector > &, double, double)
 Create another EXT ECL-crystal-crossing hit for a track candidate.
 
bool createMuidHit (ExtState &, G4ErrorFreeTrajState &, KLMMuidLikelihood *, std::vector< std::map< const Track *, double > > *)
 Create another MUID extrapolation hit for a track candidate.
 
bool findBarrelIntersection (ExtState &, const G4ThreeVector &, Intersection &)
 Find the intersection point of the track with the crossed BKLM plane.
 
bool findEndcapIntersection (ExtState &, const G4ThreeVector &, Intersection &)
 Find the intersection point of the track with the crossed EKLM plane.
 
bool findMatchingBarrelHit (Intersection &, const Track *)
 Find the matching BKLM 2D hit nearest the intersection point of the track with the crossed BKLM plane.
 
bool findMatchingEndcapHit (Intersection &, const Track *)
 Find the matching EKLM 2D hit nearest the intersection point of the track with the crossed EKLM plane.
 
void adjustIntersection (Intersection &, const double *, const G4ThreeVector &, const G4ThreeVector &)
 Nudge the track using the matching hit.
 
void finishTrack (const ExtState &, KLMMuidLikelihood *, bool)
 Complete muon identification after end of track extrapolation.
 

Private Attributes

bool m_ExtInitialized
 Flag to indicate that EXT initialize() has been called.
 
bool m_MuidInitialized
 Flag to indicate that MUID initialize() has been called.
 
double m_MeanDt
 Mean hit - trigger time (ns)
 
double m_MaxDt
 Coincidence window half-width for in-time KLM hits (ns)
 
double m_MagneticField
 Magnetic field z component (gauss) at origin.
 
double m_MaxDistSqInVariances
 user-defined maximum squared-distance (in number of variances) for matching hit to extrapolation
 
double m_MaxKLMTrackClusterDistance
 user-defined maximum distance (mm) between KLMCluster and associated track (for KLID)
 
double m_MaxECLTrackClusterDistance
 user-defined maximum distance (mm) between ECLCluster and associated track (for EID)
 
double m_MinPt
 Minimum transverse momentum in MeV/c for extrapolation to be started.
 
double m_MinKE
 Minimum kinetic energy in MeV for extrapolation to continue.
 
Simulation::ExtManagerm_ExtMgr
 Pointer to the ExtManager singleton.
 
const std::vector< Const::ChargedStable > * m_HypothesesExt
 ChargedStable hypotheses for EXT.
 
const std::vector< Const::ChargedStable > * m_HypothesesMuid
 ChargedStable hypotheses for MUID.
 
std::vector< Const::ChargedStable > * m_DefaultHypotheses
 Default ChargedStable hypotheses (needed as call argument but not used)
 
std::map< G4VPhysicalVolume *, enum VolTypes > * m_EnterExit
 Pointers to geant4 physical volumes whose entry/exit points will be saved.
 
std::vector< G4VPhysicalVolume * > * m_BKLMVolumes
 Pointers to BKLM geant4 sensitive (physical) volumes.
 
Simulation::ExtCylSurfaceTargetm_TargetExt
 virtual "target" cylinder for EXT (boundary beyond which extrapolation ends)
 
Simulation::ExtCylSurfaceTargetm_TargetMuid
 virtual "target" cylinder for MUID (boundary beyond which extrapolation ends)
 
DBObjPtr< COILGeometryParm_COILGeometryPar
 Conditions-database object for COIL geometry.
 
DBObjPtr< BeamPipeGeom_BeamPipeGeo
 Conditions-database object for beam pipe geometry.
 
double m_MinRadiusSq
 Minimum squared radius (cm) outside of which extrapolation will continue.
 
double m_OffsetZ
 offset (cm) along z axis of KLM midpoint from IP
 
int m_BarrelNSector
 Number of barrel sectors.
 
double m_BarrelMaxR
 maximum radius (cm) of the barrel
 
double m_BarrelMinR
 minimum radius (cm) of the barrel
 
double m_BarrelHalfLength
 half-length (cm) of the barrel
 
int m_OutermostActiveBarrelLayer
 outermost barrel layer that is active for muon identification (user-defined)
 
double m_BarrelPhiStripVariance [BKLMElementNumbers::getMaximalLayerNumber()+1]
 BKLM RPC phi-measuring strip position variance (cm^2) by layer.
 
double m_BarrelZStripVariance [BKLMElementNumbers::getMaximalLayerNumber()+1]
 BKLM RPC z-measuring strip position variance (cm^2) by layer.
 
double m_BarrelScintVariance
 BKLM scintillator strip position variance (cm^2)
 
double m_BarrelModuleMiddleRadius [2][BKLMElementNumbers::getMaximalSectorNumber()+1][BKLMElementNumbers::getMaximalLayerNumber()+1]
 hit-plane radius (cm) at closest distance to IP of each barrel end | sector | layer
 
G4ThreeVector m_BarrelSectorPerp [BKLMElementNumbers::getMaximalSectorNumber()+1]
 normal unit vector of each barrel sector
 
G4ThreeVector m_BarrelSectorPhi [BKLMElementNumbers::getMaximalSectorNumber()+1]
 azimuthal unit vector of each barrel sector
 
double m_EndcapMaxR
 maximum radius (cm) of the endcaps
 
double m_EndcapMinR
 minimum radius (cm) of the endcaps
 
double m_EndcapMiddleZ
 midpoint along z (cm) of the forward endcap from the KLM midpoint
 
double m_EndcapHalfLength
 half-length (cm) of either endcap
 
int m_OutermostActiveForwardEndcapLayer
 outermost forward-endcap layer that is active for muon identification (user-defined)
 
int m_OutermostActiveBackwardEndcapLayer
 outermost backward-endcap layer that is active for muon identification (user-defined)
 
double m_EndcapScintVariance
 EKLM scintillator strip position variance (cm^2)
 
double m_EndcapModuleMiddleZ [BKLMElementNumbers::getMaximalLayerNumber()+1]
 hit-plane z (cm) of each IP layer relative to KLM midpoint
 
bool m_addHitsToRecoTrack = false
 Parameter to add the found hits also to the reco tracks or not. Is turned off by default.
 
std::map< int, MuidBuilder * > m_MuidBuilderMap
 PDF for the charged final state particle hypotheses.
 
const EKLMElementNumbersm_eklmElementNumbers
 EKLM element numbers.
 
const KLMElementNumbersm_klmElementNumbers
 KLM element numbers.
 
const EKLM::TransformDataGlobalAlignedm_eklmTransformData
 EKLM transformation data.
 
DBObjPtr< KLMChannelStatusm_klmChannelStatus
 Conditions-database object for KLM channel status.
 
DBObjPtr< KLMStripEfficiencym_klmStripEfficiency
 Conditions-database object for KLM strip efficiency.
 
DBObjPtr< KLMLikelihoodParametersm_klmLikelihoodParameters
 Conditions-database object for KLM likelihood parameters.
 
StoreArray< ECLClusterm_eclClusters
 ECL clusters.
 
StoreArray< ExtHitm_extHits
 Ext hits.
 
StoreArray< KLMHit2dm_klmHit2ds
 KLM 2d hits.
 
StoreArray< KLMClusterm_klmClusters
 KLM clusters.
 
StoreArray< KLMMuidHitm_klmMuidHits
 KLM muid hits.
 
StoreArray< KLMMuidLikelihoodm_klmMuidLikelihoods
 KLM muid likelihoods.
 
StoreArray< RecoTrackm_recoTracks
 Reco tracks.
 
StoreArray< Trackm_tracks
 Tracks.
 
StoreArray< TrackClusterSeparationm_trackClusterSeparations
 Track cluster separation.
 

Static Private Attributes

static TrackExtrapolateG4em_Singleton = nullptr
 Stores pointer to the singleton class.
 

Detailed Description

geant4e-based track extrapolation.

This class extrapolates tracks outward from the outer perimeter of the CDC using geant4e.

This class requires a valid geometry in memory (gGeoManager). Therefore, a geometry building module should have been executed before this module is called.

This class has the same functions as a module - and these are called from the ExtModule's so-named functions - but also has an entry that can be called to extrapolate a single user-defined track.

Definition at line 173 of file TrackExtrapolateG4e.h.

Constructor & Destructor Documentation

◆ ~TrackExtrapolateG4e()

destructor

Definition at line 146 of file TrackExtrapolateG4e.cc.

147{
148}

◆ TrackExtrapolateG4e()

TrackExtrapolateG4e ( )
private

constructor is hidden; user calls TrackExtrapolateG4e::getInstance() instead

Definition at line 92 of file TrackExtrapolateG4e.cc.

92 :
93 m_ExtInitialized(false), // initialized later
94 m_MuidInitialized(false), // initialized later
95 m_MeanDt(0.0), // initialized later
96 m_MaxDt(0.0), // initialized later
97 m_MagneticField(0.0), // initialized later
98 m_MaxDistSqInVariances(0.0), // initialized later
99 m_MaxKLMTrackClusterDistance(0.0), // initialized later
100 m_MaxECLTrackClusterDistance(0.0), // initialized later
101 m_MinPt(0.0), // initialized later
102 m_MinKE(0.0), // initialized later
103 m_ExtMgr(nullptr), // initialized later
104 m_HypothesesExt(nullptr), // initialized later
105 m_HypothesesMuid(nullptr), // initialized later
106 m_DefaultHypotheses(nullptr), // initialized later
107 m_EnterExit(nullptr), // initialized later
108 m_BKLMVolumes(nullptr), // initialized later
109 m_TargetExt(nullptr), // initialized later
110 m_TargetMuid(nullptr), // initialized later
111 m_MinRadiusSq(0.0), // initialized later
112 m_OffsetZ(0.0), // initialized later
113 m_BarrelNSector(0), // initialized later
114 m_BarrelMaxR(0.0), // initialized later
115 m_BarrelMinR(0.0), // initialized later
116 m_BarrelHalfLength(0.0), // initialized later
117 m_OutermostActiveBarrelLayer(0), // initialized later
118 m_BarrelScintVariance(0.0), // initialized later
119 m_EndcapMaxR(0.0), // initialized later
120 m_EndcapMinR(0.0), // initialized later
121 m_EndcapMiddleZ(0.0), // initialized later
122 m_EndcapHalfLength(0.0), // initialized later
123 m_OutermostActiveForwardEndcapLayer(0), // initialized later
124 m_OutermostActiveBackwardEndcapLayer(0), // initialized later
125 m_EndcapScintVariance(0.0), // initialized later
126 m_eklmTransformData(nullptr) // initialized later
127{
128 for (int j = 0; j < BKLMElementNumbers::getMaximalLayerNumber() + 1; ++j) {
130 m_BarrelZStripVariance[j] = 0.0;
132 m_EndcapModuleMiddleZ[j] = 0.0;
133 }
134 for (int s = 0; s < BKLMElementNumbers::getMaximalSectorNumber() + 1; ++s) {
135 for (int j = 0; j < BKLMElementNumbers::getMaximalLayerNumber() + 1; ++j) {
136 m_BarrelModuleMiddleRadius[0][s][j] = 0.0;
137 m_BarrelModuleMiddleRadius[1][s][j] = 0.0;
138 }
139 m_BarrelSectorPerp[s] = G4ThreeVector(0.0, 0.0, 0.0);
140 m_BarrelSectorPhi[s] = G4ThreeVector(0.0, 0.0, 0.0);
141 }
144}
static constexpr int getMaximalLayerNumber()
Get maximal layer number (1-based).
static constexpr int getMaximalSectorNumber()
Get maximal sector number (1-based).
static const EKLMElementNumbers & Instance()
Instantiation.
static const KLMElementNumbers & Instance()
Instantiation.
double m_BarrelScintVariance
BKLM scintillator strip position variance (cm^2)
double m_EndcapHalfLength
half-length (cm) of either endcap
double m_MaxKLMTrackClusterDistance
user-defined maximum distance (mm) between KLMCluster and associated track (for KLID)
Simulation::ExtCylSurfaceTarget * m_TargetMuid
virtual "target" cylinder for MUID (boundary beyond which extrapolation ends)
int m_OutermostActiveForwardEndcapLayer
outermost forward-endcap layer that is active for muon identification (user-defined)
double m_BarrelMaxR
maximum radius (cm) of the barrel
G4ThreeVector m_BarrelSectorPhi[BKLMElementNumbers::getMaximalSectorNumber()+1]
azimuthal unit vector of each barrel sector
int m_OutermostActiveBackwardEndcapLayer
outermost backward-endcap layer that is active for muon identification (user-defined)
std::vector< Const::ChargedStable > * m_DefaultHypotheses
Default ChargedStable hypotheses (needed as call argument but not used)
const std::vector< Const::ChargedStable > * m_HypothesesExt
ChargedStable hypotheses for EXT.
double m_MagneticField
Magnetic field z component (gauss) at origin.
double m_BarrelHalfLength
half-length (cm) of the barrel
double m_MaxDt
Coincidence window half-width for in-time KLM hits (ns)
double m_BarrelMinR
minimum radius (cm) of the barrel
double m_BarrelModuleMiddleRadius[2][BKLMElementNumbers::getMaximalSectorNumber()+1][BKLMElementNumbers::getMaximalLayerNumber()+1]
hit-plane radius (cm) at closest distance to IP of each barrel end | sector | layer
double m_BarrelPhiStripVariance[BKLMElementNumbers::getMaximalLayerNumber()+1]
BKLM RPC phi-measuring strip position variance (cm^2) by layer.
double m_MinRadiusSq
Minimum squared radius (cm) outside of which extrapolation will continue.
double m_EndcapMiddleZ
midpoint along z (cm) of the forward endcap from the KLM midpoint
const EKLMElementNumbers * m_eklmElementNumbers
EKLM element numbers.
const std::vector< Const::ChargedStable > * m_HypothesesMuid
ChargedStable hypotheses for MUID.
double m_OffsetZ
offset (cm) along z axis of KLM midpoint from IP
const KLMElementNumbers * m_klmElementNumbers
KLM element numbers.
std::vector< G4VPhysicalVolume * > * m_BKLMVolumes
Pointers to BKLM geant4 sensitive (physical) volumes.
int m_OutermostActiveBarrelLayer
outermost barrel layer that is active for muon identification (user-defined)
G4ThreeVector m_BarrelSectorPerp[BKLMElementNumbers::getMaximalSectorNumber()+1]
normal unit vector of each barrel sector
double m_EndcapModuleMiddleZ[BKLMElementNumbers::getMaximalLayerNumber()+1]
hit-plane z (cm) of each IP layer relative to KLM midpoint
double m_MaxECLTrackClusterDistance
user-defined maximum distance (mm) between ECLCluster and associated track (for EID)
double m_EndcapScintVariance
EKLM scintillator strip position variance (cm^2)
double m_MinPt
Minimum transverse momentum in MeV/c for extrapolation to be started.
int m_BarrelNSector
Number of barrel sectors.
std::map< G4VPhysicalVolume *, enum VolTypes > * m_EnterExit
Pointers to geant4 physical volumes whose entry/exit points will be saved.
bool m_MuidInitialized
Flag to indicate that MUID initialize() has been called.
double m_EndcapMinR
minimum radius (cm) of the endcaps
double m_MeanDt
Mean hit - trigger time (ns)
double m_EndcapMaxR
maximum radius (cm) of the endcaps
Simulation::ExtCylSurfaceTarget * m_TargetExt
virtual "target" cylinder for EXT (boundary beyond which extrapolation ends)
bool m_ExtInitialized
Flag to indicate that EXT initialize() has been called.
double m_MinKE
Minimum kinetic energy in MeV for extrapolation to continue.
Simulation::ExtManager * m_ExtMgr
Pointer to the ExtManager singleton.
double m_MaxDistSqInVariances
user-defined maximum squared-distance (in number of variances) for matching hit to extrapolation
double m_BarrelZStripVariance[BKLMElementNumbers::getMaximalLayerNumber()+1]
BKLM RPC z-measuring strip position variance (cm^2) by layer.
const EKLM::TransformDataGlobalAligned * m_eklmTransformData
EKLM transformation data.

Member Function Documentation

◆ adjustIntersection()

void adjustIntersection ( Intersection & ,
const double * ,
const G4ThreeVector & ,
const G4ThreeVector &  )
private

Nudge the track using the matching hit.

Definition at line 1683 of file TrackExtrapolateG4e.cc.

1685{
1686 // Use the gain matrix formalism to get the corrected track parameters.
1687 // R. Fruhwirth, Application of Kalman Filtering, NIM A262 (1987) 444-450
1688 // Equations (7)
1689 // x_k^{k-1} = extPar[] 6 elements before filtering
1690 // C_k^{k-1} = extCov[] 6x6 symmetric before filtering
1691 // r_k^{k-1} = residual[] 2 elements before filtering
1692 // h_k = 2x6 projects cartesian coordinates to measurement-plane coordinates
1693 // H_k = @h_k/@x = jacobian[] 2x6 Jacobian of projection to measurement plane
1694 // R_k^{k-1} = correction[] 2x2 before Invert()
1695 // G_k = R^(-1) = correction[] 2x2 after Invert()
1696 // K_k = gain[] 6x2 Kalman gain matrix
1697 // x_k = extPar[] 6 elements after filtering
1698 // C_k = extCov[] 6x6 symmetric after filtering
1699 // r_k = residual[] 2 elements after filtering
1700 // Use the relation K*H*C = (C*H^T*R^-1)*H*C = C*(H^T*R^-1*H)*C^T
1701
1702 // In most cases, extPos0 is the same as intersection.position. They differ only when
1703 // the nearest BKLM hit is in the sector adjacent to that of intersection.position.
1704 G4ThreeVector extPos(extPos0);
1705 G4ThreeVector extMom(intersection.momentum);
1706 G4ThreeVector extDir(extMom.unit());
1707 G4ThreeVector diffPos(hitPos - extPos);
1708 G4ErrorSymMatrix extCov(intersection.covariance);
1709 // Track parameters (x,y,z,px,py,pz) before correction
1710 G4ErrorMatrix extPar(6, 1); // initialized to all zeroes
1711 extPar[0][0] = extPos.x();
1712 extPar[1][0] = extPos.y();
1713 extPar[2][0] = extPos.z();
1714 extPar[3][0] = extMom.x();
1715 extPar[4][0] = extMom.y();
1716 extPar[5][0] = extMom.z();
1717 G4ThreeVector nA; // unit vector normal to the readout plane
1718 G4ThreeVector nB; // unit vector along phi- or x-readout direction (for barrel or endcap)
1719 G4ThreeVector nC; // unit vector along z- or y-readout direction (for barrel or endcap)
1720 if (intersection.inBarrel) {
1721 nA = m_BarrelSectorPerp[intersection.sector];
1722 nB = m_BarrelSectorPhi[intersection.sector];
1723 nC = G4ThreeVector(0.0, 0.0, 1.0);
1724 } else {
1725 double out = (intersection.isForward ? 1.0 : -1.0);
1726 nA = G4ThreeVector(0.0, 0.0, out);
1727 nB = G4ThreeVector(out, 0.0, 0.0);
1728 nC = G4ThreeVector(0.0, out, 0.0);
1729 }
1730 // Don't adjust the extrapolation if the track is nearly tangent to the readout plane.
1731 double extDirA = extDir * nA;
1732 if (std::fabs(extDirA) < 1.0E-6)
1733 return;
1734 double extDirBA = extDir * nB / extDirA;
1735 double extDirCA = extDir * nC / extDirA;
1736 // Move the extrapolated coordinate (at most a tiny amount!) to the plane of the hit.
1737 // If the moved point is outside the KLM, don't do Kalman filtering.
1738 G4ThreeVector move = extDir * ((diffPos * nA) / extDirA);
1739 extPos += move;
1740 diffPos -= move;
1741 intersection.positionAtHitPlane = extPos;
1742 // Projection jacobian onto the nB-nC measurement plane
1743 G4ErrorMatrix jacobian(2, 6); // initialized to all zeroes
1744 jacobian[0][0] = nB.x() - nA.x() * extDirBA;
1745 jacobian[0][1] = nB.y() - nA.y() * extDirBA;
1746 jacobian[0][2] = nB.z() - nA.z() * extDirBA;
1747 jacobian[1][0] = nC.x() - nA.x() * extDirCA;
1748 jacobian[1][1] = nC.y() - nA.y() * extDirCA;
1749 jacobian[1][2] = nC.z() - nA.z() * extDirCA;
1750 // Residuals of EXT track and KLM hit on the nB-nC measurement plane
1751 G4ErrorMatrix residual(2, 1); // initialized to all zeroes
1752 residual[0][0] = diffPos.x() * jacobian[0][0] + diffPos.y() * jacobian[0][1] + diffPos.z() * jacobian[0][2];
1753 residual[1][0] = diffPos.x() * jacobian[1][0] + diffPos.y() * jacobian[1][1] + diffPos.z() * jacobian[1][2];
1754 // Measurement errors in the detector plane
1755 G4ErrorSymMatrix hitCov(2, 0); // initialized to all zeroes
1756 hitCov[0][0] = localVariance[0];
1757 hitCov[1][1] = localVariance[1];
1758 // No magnetic field: increase the hit uncertainty
1759 if (m_MagneticField == 0.0) {
1760 hitCov[0][0] *= 10.0;
1761 hitCov[1][1] *= 10.0;
1762 }
1763 // Now get the correction matrix: combined covariance of EXT and KLM hit.
1764 // 1st dimension = nB, 2nd dimension = nC.
1765 G4ErrorSymMatrix correction(extCov.similarity(jacobian) + hitCov);
1766 // Ignore the best hit if it is too far from the extrapolated-track intersection in the hit's plane
1767 if (residual[0][0] * residual[0][0] > correction[0][0] * m_MaxDistSqInVariances)
1768 return;
1769 if (residual[1][0] * residual[1][0] > correction[1][1] * m_MaxDistSqInVariances)
1770 return;
1771 int fail = 0;
1772 correction.invert(fail);
1773 if (fail != 0)
1774 return;
1775 // Matrix inversion succeeded and is reasonable.
1776 // Evaluate chi-squared increment assuming that the Kalman filter
1777 // won't be able to adjust the extrapolated track's position (fall-back).
1778 intersection.chi2 = (correction.similarityT(residual))[0][0];
1779 // Do the Kalman filtering
1780 G4ErrorMatrix gain((extCov * jacobian.T()) * correction);
1781 G4ErrorSymMatrix HRH(correction.similarityT(jacobian));
1782 extCov -= HRH.similarity(extCov);
1783 extPar += gain * residual;
1784 extPos.set(extPar[0][0], extPar[1][0], extPar[2][0]);
1785 extMom.set(extPar[3][0], extPar[4][0], extPar[5][0]);
1786 // Calculate the chi-squared increment using the Kalman-filtered state
1787 correction = hitCov - extCov.similarity(jacobian);
1788 correction.invert(fail);
1789 if (fail != 0)
1790 return;
1791 diffPos = hitPos - extPos;
1792 residual[0][0] = diffPos.x() * jacobian[0][0] + diffPos.y() * jacobian[0][1] + diffPos.z() * jacobian[0][2];
1793 residual[1][0] = diffPos.x() * jacobian[1][0] + diffPos.y() * jacobian[1][1] + diffPos.z() * jacobian[1][2];
1794 intersection.chi2 = (correction.similarityT(residual))[0][0];
1795 // Update the position, momentum and covariance of the point
1796 // Project the corrected extrapolation to the plane of the original
1797 // extrapolation's intersection.position. (Note: intersection.position is the same as
1798 // extPos0 in all cases except when nearest BKLM hit is in adjacent
1799 // sector, for which extPos0 is a projected position to the hit's plane.)
1800 // Also, leave the momentum magnitude unchanged.
1801 intersection.position = extPos + extDir * (((intersection.position - extPos) * nA) / extDirA);
1802 intersection.momentum = intersection.momentum.mag() * extMom.unit();
1803 intersection.covariance = extCov;
1804}

◆ beginRun()

void beginRun ( bool flag)

Perform beginning-of-run actions.

Parameters
flagTrue if called by Muid module, false if called by Ext module.

Definition at line 346 of file TrackExtrapolateG4e.cc.

347{
348 B2DEBUG(20, (byMuid ? "muid" : "ext"));
349 if (byMuid) {
350 if (!m_klmChannelStatus.isValid())
351 B2FATAL("KLM channel status data are not available.");
352 if (!m_klmStripEfficiency.isValid())
353 B2FATAL("KLM strip efficiency data are not available.");
354 if (!m_klmLikelihoodParameters.isValid())
355 B2FATAL("KLM likelihood parameters are not available.");
356 std::vector<int> muidPdgCodes = MuidElementNumbers::getPDGVector();
357 if (!m_MuidBuilderMap.empty()) {
358 if (m_klmLikelihoodParameters.hasChanged()) { /* Clear m_MuidBuilderMap if KLMLikelihoodParameters payload changed. */
359 for (auto const& muidBuilder : m_MuidBuilderMap)
360 delete muidBuilder.second;
361 m_MuidBuilderMap.clear();
362 } else /* Return if m_MuidBuilderMap is already initialized. */
363 return;
364 }
365 for (int pdg : muidPdgCodes)
366 m_MuidBuilderMap.insert(std::pair<int, MuidBuilder*>(pdg, new MuidBuilder(pdg)));
367 }
368}
static std::vector< int > getPDGVector()
Get a vector with all the hypothesis PDG codes used for Muid.
DBObjPtr< KLMLikelihoodParameters > m_klmLikelihoodParameters
Conditions-database object for KLM likelihood parameters.
std::map< int, MuidBuilder * > m_MuidBuilderMap
PDF for the charged final state particle hypotheses.
DBObjPtr< KLMChannelStatus > m_klmChannelStatus
Conditions-database object for KLM channel status.
DBObjPtr< KLMStripEfficiency > m_klmStripEfficiency
Conditions-database object for KLM strip efficiency.

◆ createExtHit()

void createExtHit ( const ExtHitStatus status,
const ExtState & extState,
const G4ErrorFreeTrajState & g4eState,
const G4StepPoint * stepPoint,
const G4TouchableHandle & touch )
private

Create another EXT extrapolation hit for a track candidate.

Definition at line 1226 of file TrackExtrapolateG4e.cc.

1229{
1230 Const::EDetector detID(Const::EDetector::invalidDetector);
1231 int copyID(0);
1232 getVolumeID(touch, detID, copyID);
1233 G4ThreeVector pos(stepPoint->GetPosition() / CLHEP::cm);
1234 G4ThreeVector mom(stepPoint->GetMomentum() / CLHEP::GeV);
1235 if (extState.isCosmic)
1236 mom = -mom;
1237 G4ErrorSymMatrix covariance(6, 0);
1238 fromG4eToPhasespace(g4eState, covariance);
1239 ExtHit* extHit = m_extHits.appendNew(extState.tof, extState.pdgCode, detID, copyID, status,
1240 extState.isCosmic,
1241 pos, mom, covariance);
1242 // If called standalone, there will be no associated track
1243 if (extState.track != nullptr)
1244 extState.track->addRelationTo(extHit);
1245}
EDetector
Enum for identifying the detector components (detector and subdetector).
Definition Const.h:42
void addRelationTo(const RelationsInterface< BASE > *object, float weight=1.0, const std::string &namedRelation="") const
Add a relation from this object to another object (with caching).
void getVolumeID(const G4TouchableHandle &, Const::EDetector &, int &)
Get the physical volume information for a geant4 physical volume.
StoreArray< ExtHit > m_extHits
Ext hits.
void fromG4eToPhasespace(const G4ErrorFreeTrajState &, G4ErrorSymMatrix &)
Convert the geant4e 5x5 covariance to phasespace 6x6 covariance.
bool isCosmic
True for back-propagation of a cosmic ray.
int pdgCode
Particle hypothesis that is being extrapolated.
double tof
Time of flight from IP (ns), updated during extrapolation.
const Track * track
Pointer to the reconstructed track.

◆ createMuidHit()

bool createMuidHit ( ExtState & extState,
G4ErrorFreeTrajState & g4eState,
KLMMuidLikelihood * klmMuidLikelihood,
std::vector< std::map< const Track *, double > > * bklmHitUsed )
private

Create another MUID extrapolation hit for a track candidate.

Definition at line 1250 of file TrackExtrapolateG4e.cc.

1252{
1253
1254 Intersection intersection;
1255 intersection.hit = -1;
1256 intersection.chi2 = -1.0;
1257 intersection.position = g4eState.GetPosition() / CLHEP::cm;
1258 intersection.momentum = g4eState.GetMomentum() / CLHEP::GeV;
1259 G4ThreeVector prePos = g4eState.GetG4Track()->GetStep()->GetPreStepPoint()->GetPosition() / CLHEP::cm;
1260 G4ThreeVector oldPosition(prePos.x(), prePos.y(), prePos.z());
1261 double r = intersection.position.perp();
1262 double z = std::fabs(intersection.position.z() - m_OffsetZ);
1263
1264 // Is the track in the barrel?
1265 if ((r > m_BarrelMinR) && (r < m_BarrelMaxR) && (z < m_BarrelHalfLength)) {
1266 // Did the track cross the inner midplane of a detector module?
1267 if (findBarrelIntersection(extState, oldPosition, intersection)) {
1268 fromG4eToPhasespace(g4eState, intersection.covariance);
1269 if (findMatchingBarrelHit(intersection, extState.track)) {
1270 (*bklmHitUsed)[intersection.hit].insert(std::pair<const Track*, double>(extState.track, intersection.chi2));
1271 extState.extLayerPattern |= (0x00000001 << intersection.layer);
1272 float layerBarrelEfficiency = 1.;
1273 for (int plane = 0; plane <= BKLMElementNumbers::getMaximalPlaneNumber(); plane++) {
1274 layerBarrelEfficiency *= m_klmStripEfficiency->getBarrelEfficiency(
1276 intersection.sector + 1, intersection.layer + 1, plane, 1);
1277 }
1278 klmMuidLikelihood->setExtBKLMEfficiencyValue(intersection.layer, layerBarrelEfficiency);
1279 if (extState.lastBarrelExtLayer < intersection.layer) {
1280 extState.lastBarrelExtLayer = intersection.layer;
1281 }
1282 extState.hitLayerPattern |= (0x00000001 << intersection.layer);
1283 if (extState.lastBarrelHitLayer < intersection.layer) {
1284 extState.lastBarrelHitLayer = intersection.layer;
1285 }
1286 // If the updated point is outside the barrel, discard it and the Kalman-fitter adjustment
1287 r = intersection.position.perp();
1288 z = std::fabs(intersection.position.z() - m_OffsetZ);
1289 if ((r <= m_BarrelMinR) || (r >= m_BarrelMaxR) || (z >= m_BarrelHalfLength)) {
1290 intersection.chi2 = -1.0;
1291 }
1292 } else {
1293 // Record a no-hit track crossing if this step is strictly within a barrel sensitive volume
1294 std::vector<G4VPhysicalVolume*>::iterator j = find(m_BKLMVolumes->begin(), m_BKLMVolumes->end(),
1295 g4eState.GetG4Track()->GetVolume());
1296 if (j != m_BKLMVolumes->end()) {
1297 bool isDead = true; // by default, the nearest orthogonal strips are dead
1298 int section = intersection.isForward ?
1301 int sector = intersection.sector + 1; // from 0-based to 1-based enumeration
1302 int layer = intersection.layer + 1; // from 0-based to 1-based enumeration
1303 const bklm::Module* m = bklm::GeometryPar::instance()->findModule(section, sector, layer); // uses 1-based enumeration
1304 if (m) {
1305 const CLHEP::Hep3Vector localPosition = m->globalToLocal(intersection.position); // uses and returns position in cm
1306 int zStrip = m->getZStripNumber(localPosition);
1307 int phiStrip = m->getPhiStripNumber(localPosition);
1308 if (zStrip >= 0 && phiStrip >= 0) {
1309 uint16_t channel1, channel2;
1310 channel1 = m_klmElementNumbers->channelNumberBKLM(
1311 section, sector, layer,
1313 channel2 = m_klmElementNumbers->channelNumberBKLM(
1314 section, sector, layer,
1316 enum KLMChannelStatus::ChannelStatus status1, status2;
1317 status1 = m_klmChannelStatus->getChannelStatus(channel1);
1318 status2 = m_klmChannelStatus->getChannelStatus(channel2);
1319 if (status1 == KLMChannelStatus::c_Unknown ||
1320 status2 == KLMChannelStatus::c_Unknown)
1321 B2ERROR("No KLM channel status data."
1322 << LogVar("Section", section)
1323 << LogVar("Sector", sector)
1324 << LogVar("Layer", layer)
1325 << LogVar("Z strip", zStrip)
1326 << LogVar("Phi strip", phiStrip));
1327 isDead = (status1 == KLMChannelStatus::c_Dead ||
1328 status2 == KLMChannelStatus::c_Dead);
1329 }
1330 }
1331 if (!isDead) {
1332 extState.extLayerPattern |= (0x00000001 << intersection.layer); // valid extrapolation-crossing of the layer but no matching hit
1333 float layerBarrelEfficiency = 1.;
1334 for (int plane = 0; plane <= BKLMElementNumbers::getMaximalPlaneNumber(); plane++) {
1335 layerBarrelEfficiency *= m_klmStripEfficiency->getBarrelEfficiency(
1337 intersection.sector + 1, intersection.layer + 1, plane, 1);
1338 }
1339 klmMuidLikelihood->setExtBKLMEfficiencyValue(intersection.layer, layerBarrelEfficiency);
1340 } else {
1341 klmMuidLikelihood->setExtBKLMEfficiencyValue(intersection.layer, 0);
1342 }
1343 if (extState.lastBarrelExtLayer < intersection.layer) {
1344 extState.lastBarrelExtLayer = intersection.layer;
1345 }
1346 }
1347 }
1348 }
1349 }
1350
1351 // Is the track in the endcap?
1352 if ((r > m_EndcapMinR) && (std::fabs(z - m_EndcapMiddleZ) < m_EndcapHalfLength)) {
1353 // Did the track cross the inner midplane of a detector module?
1354 if (findEndcapIntersection(extState, oldPosition, intersection)) {
1355 fromG4eToPhasespace(g4eState, intersection.covariance);
1356 if (findMatchingEndcapHit(intersection, extState.track)) {
1357 extState.extLayerPattern |= (0x00008000 << intersection.layer);
1358 float layerEndcapEfficiency = 1.;
1359 for (int plane = 1; plane <= EKLMElementNumbers::getMaximalPlaneNumber(); plane++) {
1360 layerEndcapEfficiency *= m_klmStripEfficiency->getEndcapEfficiency(
1362 intersection.sector + 1, intersection.layer + 1, plane, 1);
1363 }
1364 klmMuidLikelihood->setExtEKLMEfficiencyValue(intersection.layer, layerEndcapEfficiency);
1365 if (extState.lastEndcapExtLayer < intersection.layer) {
1366 extState.lastEndcapExtLayer = intersection.layer;
1367 }
1368 extState.hitLayerPattern |= (0x00008000 << intersection.layer);
1369 if (extState.lastEndcapHitLayer < intersection.layer) {
1370 extState.lastEndcapHitLayer = intersection.layer;
1371 }
1372 // If the updated point is outside the endcap, discard it and the Kalman-fitter adjustment
1373 r = intersection.position.perp();
1374 z = std::fabs(intersection.position.z() - m_OffsetZ);
1375 if ((r <= m_EndcapMinR) || (r >= m_EndcapMaxR) || (std::fabs(z - m_EndcapMiddleZ) >= m_EndcapHalfLength)) {
1376 intersection.chi2 = -1.0;
1377 }
1378 } else {
1379 bool isDead = true;
1380 int result, strip1, strip2;
1381 result = m_eklmTransformData->getStripsByIntersection(
1382 intersection.position, &strip1, &strip2);
1383 if (result == 0) {
1384 uint16_t channel1, channel2;
1385 channel1 = m_klmElementNumbers->channelNumberEKLM(strip1);
1386 channel2 = m_klmElementNumbers->channelNumberEKLM(strip2);
1387 enum KLMChannelStatus::ChannelStatus status1, status2;
1388 status1 = m_klmChannelStatus->getChannelStatus(channel1);
1389 status2 = m_klmChannelStatus->getChannelStatus(channel2);
1390 if (status1 == KLMChannelStatus::c_Unknown ||
1391 status2 == KLMChannelStatus::c_Unknown)
1392 B2ERROR("Incomplete KLM channel status data.");
1393 isDead = (status1 == KLMChannelStatus::c_Dead ||
1394 status2 == KLMChannelStatus::c_Dead);
1395 }
1396 if (!isDead) {
1397 extState.extLayerPattern |= (0x00008000 << intersection.layer); // valid extrapolation-crossing of the layer but no matching hit
1398 float layerEndcapEfficiency = 1.;
1399 for (int plane = 1; plane <= EKLMElementNumbers::getMaximalPlaneNumber(); plane++) {
1400 layerEndcapEfficiency *= m_klmStripEfficiency->getEndcapEfficiency(
1402 intersection.sector + 1, intersection.layer + 1, plane, 1);
1403 }
1404 klmMuidLikelihood->setExtEKLMEfficiencyValue(intersection.layer, layerEndcapEfficiency);
1405 } else {
1406 klmMuidLikelihood->setExtEKLMEfficiencyValue(intersection.layer, 0);
1407 }
1408 if (extState.lastEndcapExtLayer < intersection.layer) {
1409 extState.lastEndcapExtLayer = intersection.layer;
1410 }
1411 }
1412 }
1413 }
1414
1415 // Create a new MuidHit and RelationEntry between it and the track.
1416 // Adjust geant4e's position, momentum and covariance based on matching hit and tell caller to update the geant4e state.
1417 if (intersection.chi2 >= 0.0) {
1418 ROOT::Math::XYZVector tpos(intersection.position.x(),
1419 intersection.position.y(),
1420 intersection.position.z());
1421 ROOT::Math::XYZVector tposAtHitPlane(intersection.positionAtHitPlane.x(),
1422 intersection.positionAtHitPlane.y(),
1423 intersection.positionAtHitPlane.z());
1424 KLMMuidHit* klmMuidHit = m_klmMuidHits.appendNew(extState.pdgCode,
1425 intersection.inBarrel, intersection.isForward,
1426 intersection.sector, intersection.layer,
1427 tpos, tposAtHitPlane,
1428 extState.tof, intersection.time,
1429 intersection.chi2);
1430 if (extState.track != nullptr) { extState.track->addRelationTo(klmMuidHit); }
1431 G4Point3D newPos(intersection.position.x() * CLHEP::cm,
1432 intersection.position.y() * CLHEP::cm,
1433 intersection.position.z() * CLHEP::cm);
1434 g4eState.SetPosition(newPos);
1435 G4Vector3D newMom(intersection.momentum.x() * CLHEP::GeV,
1436 intersection.momentum.y() * CLHEP::GeV,
1437 intersection.momentum.z() * CLHEP::GeV);
1438 g4eState.SetMomentum(newMom);
1439 G4ErrorTrajErr covG4e;
1440 fromPhasespaceToG4e(intersection.momentum, intersection.covariance, covG4e);
1441 g4eState.SetError(covG4e);
1442 extState.chi2 += intersection.chi2;
1443 extState.nPoint += 2; // two (orthogonal) independent hits per detector layer
1444 return true;
1445 }
1446
1447 // Tell caller that the geant4e state was not modified.
1448 return false;
1449
1450}
static constexpr int getMaximalPlaneNumber()
Get maximal plane number (0-based).
static constexpr int getMaximalPlaneNumber()
Get maximal plane number.
@ c_Dead
Dead channel (no signal).
@ c_Unknown
Unknown status (no data).
void setExtEKLMEfficiencyValue(int layer, float efficiency)
Set the efficiency of a given EKLM layer.
void setExtBKLMEfficiencyValue(int layer, float efficiency)
Set the efficiency of a given BKLM layer.
bool findEndcapIntersection(ExtState &, const G4ThreeVector &, Intersection &)
Find the intersection point of the track with the crossed EKLM plane.
StoreArray< KLMMuidHit > m_klmMuidHits
KLM muid hits.
bool findMatchingEndcapHit(Intersection &, const Track *)
Find the matching EKLM 2D hit nearest the intersection point of the track with the crossed EKLM plane...
void fromPhasespaceToG4e(const G4ThreeVector &, const G4ErrorSymMatrix &, G4ErrorTrajErr &)
Convert the phasespace covariance to geant4e covariance.
bool findMatchingBarrelHit(Intersection &, const Track *)
Find the matching BKLM 2D hit nearest the intersection point of the track with the crossed BKLM plane...
bool findBarrelIntersection(ExtState &, const G4ThreeVector &, Intersection &)
Find the intersection point of the track with the crossed BKLM plane.
const Module * findModule(int section, int sector, int layer) const
Get the pointer to the definition of a module.
static GeometryPar * instance(void)
Static method to get a reference to the singleton GeometryPar instance.
int lastEndcapExtLayer
MUID: outermost endcap layer crossed by the extrapolated track.
double chi2
MUID: accumulated chi-squared of all in-plane transverse deviations between extrapolation and matchin...
int nPoint
MUID: accumulated number of points with matching 2D hits.
int lastEndcapHitLayer
MUID: outermost endcap layer with a matching hit.
int extLayerPattern
MUID: accumulated bit pattern of layers crossed by the extrapolated track.
int lastBarrelHitLayer
MUID: outermost barrel layer with a matching hit.
int hitLayerPattern
MUID: accumulated bit pattern of layers with matching hits.
int lastBarrelExtLayer
MUID: outermost barrel layer crossed by the extrapolated track.
int sector
sector number (0..7 for barrel, 0..3 for endcap) of this point
G4ThreeVector momentum
extrapolated-track momentum (GeV/c) at this intersection
double chi2
chi-squared value of transverse deviation between extrapolated and measured hit positions
bool inBarrel
flag to indicate if this point is in the barrel (true) or endcap (false)
int hit
index in {B,E}KLMHit2ds of matching hit
G4ThreeVector positionAtHitPlane
extrapolated-track position (cm) projected to the 2D hit's midplane
G4ThreeVector position
extrapolated-track global position (cm) of this intersection
G4ErrorSymMatrix covariance
extrapolated-track phase-space covariance matrix at this intersection
bool isForward
flag to indicate if this point is in the forward (true) or backward (false) end
double time
time (ns) of matching BKLMHit2d
int layer
layer number (0..14 for barrel, 0..13 for endcap) of this point

◆ endRun()

void endRun ( bool flag)

Perform end-of-run actions.

Parameters
flagTrue if called by Muid module, false if called by Ext module.

Definition at line 429 of file TrackExtrapolateG4e.cc.

430{
431}

◆ event()

void event ( bool flag)

Performs track extrapolation for all tracks in one event.

Parameters
flagTrue if called by Muid module, false if called by Ext module.

Definition at line 370 of file TrackExtrapolateG4e.cc.

371{
372
373 // Put geant4 in proper state (in case this module is in a separate process)
374 if (G4StateManager::GetStateManager()->GetCurrentState() == G4State_Idle) {
375 G4StateManager::GetStateManager()->SetNewState(G4State_GeomClosed);
376 }
377
378 G4ThreeVector directionAtIP, positionG4e, momentumG4e;
379 G4ErrorTrajErr covG4e(5); // initialized to zeroes
380
381 // Loop over the reconstructed tracks
382 // Do extrapolation for each hypothesis of each reconstructed track.
383 if (byMuid) { // event() called by Muid module
384 G4ErrorPropagatorData::GetErrorPropagatorData()->SetTarget(m_TargetMuid);
385 std::vector<std::pair<ECLCluster*, G4ThreeVector> > eclClusterInfo(m_eclClusters.getEntries());
386 for (int c = 0; c < m_eclClusters.getEntries(); ++c) {
387 eclClusterInfo[c].first = m_eclClusters[c];
388 eclClusterInfo[c].second = G4ThreeVector(m_eclClusters[c]->getClusterPosition().X(),
389 m_eclClusters[c]->getClusterPosition().Y(),
390 m_eclClusters[c]->getClusterPosition().Z()) * CLHEP::cm;
391 }
392 std::vector<std::pair<KLMCluster*, G4ThreeVector> > klmClusterInfo(m_klmClusters.getEntries());
393 for (int c = 0; c < m_klmClusters.getEntries(); ++c) {
394 klmClusterInfo[c].first = m_klmClusters[c];
395 klmClusterInfo[c].second = G4ThreeVector(m_klmClusters[c]->getClusterPosition().X(),
396 m_klmClusters[c]->getClusterPosition().Y(),
397 m_klmClusters[c]->getClusterPosition().Z()) * CLHEP::cm;
398 }
399 // Keep track of (re-)use of BKLMHit2ds
400 std::vector<std::map<const Track*, double> > bklmHitUsed(m_klmHit2ds.getEntries());
401 for (auto& b2track : m_tracks) {
402 for (const auto& hypothesis : *m_HypothesesMuid) {
403 int pdgCode = hypothesis.getPDGCode();
404 if (hypothesis == Const::electron || hypothesis == Const::muon)
405 pdgCode = -pdgCode;
406 G4ErrorFreeTrajState g4eState("g4e_mu+", G4ThreeVector(), G4ThreeVector()); // will be updated
407 ExtState extState = getStartPoint(b2track, pdgCode, g4eState);
408 swim(extState, g4eState, &eclClusterInfo, &klmClusterInfo, &bklmHitUsed);
409 } // Muid hypothesis loop
410 } // Muid track loop
411 } else { // event() called by Ext module
412 G4ErrorPropagatorData::GetErrorPropagatorData()->SetTarget(m_TargetExt);
413 for (auto& b2track : m_tracks) {
414 for (const auto& hypothesis : *m_HypothesesExt) {
415 int pdgCode = hypothesis.getPDGCode();
416 if (hypothesis == Const::electron || hypothesis == Const::muon) pdgCode = -pdgCode;
417 G4ErrorFreeTrajState g4eState("g4e_mu+", G4ThreeVector(), G4ThreeVector()); // will be updated
418 ExtState extState = getStartPoint(b2track, pdgCode, g4eState);
419 swim(extState, g4eState);
420 } // Ext hypothesis loop
421 } // Ext track loop
422 } // byMuid
423
424 if (byMuid) {
426 }
427}
static const ChargedStable muon
muon particle
Definition Const.h:660
static const ChargedStable electron
electron particle
Definition Const.h:659
ExtState getStartPoint(const Track &, int, G4ErrorFreeTrajState &)
Get the start point for a new reconstructed track with specific PDG hypothesis.
StoreArray< KLMCluster > m_klmClusters
KLM clusters.
StoreArray< Track > m_tracks
Tracks.
StoreArray< ECLCluster > m_eclClusters
ECL clusters.
StoreArray< KLMHit2d > m_klmHit2ds
KLM 2d hits.
void findClosestTrackToKLMClusters()
Find the closest Track to each KLMCluster and fill the corresponding fields of KLMCluster objects.
void swim(ExtState &, G4ErrorFreeTrajState &, const std::vector< std::pair< ECLCluster *, G4ThreeVector > > *, const std::vector< std::pair< KLMCluster *, G4ThreeVector > > *, std::vector< std::map< const Track *, double > > *)
Swim a single track (MUID) until it stops or leaves the target cylinder.

◆ extrapolate()

void extrapolate ( int pdgCode,
double tof,
const G4ThreeVector & position,
const G4ThreeVector & momentum,
const G4ErrorSymMatrix & covariance )

Performs track extrapolation for a single track (specified in genfit2 units).

Parameters
pdgCodeSigned PDG identifier of the particle hypothesis to be used for the extrapolation.
tofStarting time, i.e., time of flight from the IP, at the starting point (ns).
positionStarting point of the extrapolation (cm).
momentumMomentum of the track at the starting point (GeV/c).
covariancePhase-space covariance matrix (6x6) at the starting point (cm, GeV/c).

Definition at line 455 of file TrackExtrapolateG4e.cc.

461{
462 bool isCosmic = false; // DIVOT
463 if ((!m_ExtInitialized) && (!m_MuidInitialized)) {
464 // No EXT nor MUID module in analysis path ==> mimic ext::initialize() with reasonable defaults.
465 // The default values are taken from the EXT module's parameter definitions.
466 Simulation::ExtManager* extMgr = Simulation::ExtManager::GetManager();
467 extMgr->Initialize("Ext", "default", 0.0, 0.25, false, 0, std::vector<std::string>());
468 // Redefine geant4e step length, magnetic field step limitation (fraction of local curvature radius),
469 // and kinetic energy loss limitation (maximum fractional energy loss) by communicating with
470 // the geant4 UI. (Commands were defined in ExtMessenger when physics list was set up.)
471 // *NOTE* If module muid runs after this, its G4UImanager commands will override these.
472 G4UImanager::GetUIpointer()->ApplyCommand("/geant4e/limits/stepLength 250 mm");
473 G4UImanager::GetUIpointer()->ApplyCommand("/geant4e/limits/magField 0.001");
474 G4UImanager::GetUIpointer()->ApplyCommand("/geant4e/limits/energyLoss 0.05");
475 m_DefaultHypotheses = new std::vector<Const::ChargedStable>; // not used
476 initialize(0.1, 0.002, *m_DefaultHypotheses);
477 }
478
479 // Put geant4 in proper state (in case this module is in a separate process)
480 if (G4StateManager::GetStateManager()->GetCurrentState() == G4State_Idle) {
481 G4StateManager::GetStateManager()->SetNewState(G4State_GeomClosed);
482 }
483
484 G4ErrorPropagatorData::GetErrorPropagatorData()->SetTarget(m_TargetExt);
485
486 // Do extrapolation for selected hypothesis (pion, electron, muon, kaon, proton,
487 // deuteron) for the selected track until calorimeter exit.
488
489 G4ThreeVector positionG4e = position * CLHEP::cm; // convert from genfit2 units (cm) to geant4 units (mm)
490 G4ThreeVector momentumG4e = momentum * CLHEP::GeV; // convert from genfit2 units (GeV/c) to geant4 units (MeV/c)
491 if (isCosmic)
492 momentumG4e = -momentumG4e;
493 G4ErrorSymMatrix covarianceG4e(5, 0); // in Geant4e units (GeV/c, cm)
494 fromPhasespaceToG4e(momentum, covariance, covarianceG4e);
495 G4String nameG4e("g4e_" + G4ParticleTable::GetParticleTable()->FindParticle(pdgCode)->GetParticleName());
496 G4ErrorFreeTrajState g4eState(nameG4e, positionG4e, momentumG4e, covarianceG4e);
497 ExtState extState = { nullptr, pdgCode, isCosmic, tof, 0.0, // for EXT and MUID
498 momentumG4e.unit(), 0.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, false // for MUID only
499 };
500 swim(extState, g4eState);
501}
static ExtManager * GetManager()
Get pointer to the instance of this singleton class (create if needed)
Definition ExtManager.cc:71
void initialize(double minPt, double minKE, std::vector< Const::ChargedStable > &hypotheses)
Initialize for track extrapolation by the EXT module.

◆ findBarrelIntersection()

bool findBarrelIntersection ( ExtState & extState,
const G4ThreeVector & oldPosition,
Intersection & intersection )
private

Find the intersection point of the track with the crossed BKLM plane.

Definition at line 1452 of file TrackExtrapolateG4e.cc.

1453{
1454 // Be generous: allow outward-moving intersection to be in the dead space between
1455 // largest sensitive-volume Z and m_BarrelHalfLength, not necessarily in a geant4 sensitive volume
1456 if (std::fabs(intersection.position.z() - m_OffsetZ) > m_BarrelHalfLength)
1457 return false;
1458 double phi = intersection.position.phi();
1459 if (phi < 0.0)
1460 phi += TWOPI;
1461 if (phi > TWOPI - PI_8)
1462 phi -= TWOPI;
1463 int sector = (int)((phi + PI_8) / M_PI_4);
1464 int section = intersection.position.z() > m_OffsetZ ?
1467 double oldR = oldPosition * m_BarrelSectorPerp[sector];
1468 double newR = intersection.position * m_BarrelSectorPerp[sector];
1469 for (int layer = extState.firstBarrelLayer; layer <= m_OutermostActiveBarrelLayer; ++layer) {
1470 if (newR < m_BarrelModuleMiddleRadius[section][sector][layer]) break;
1471 if (oldR <= m_BarrelModuleMiddleRadius[section][sector][layer]) {
1472 extState.firstBarrelLayer = layer + 1; // ratchet outward for next call's loop starting value
1473 if (extState.firstBarrelLayer > m_OutermostActiveBarrelLayer) extState.escaped = true;
1474 intersection.inBarrel = true;
1475 intersection.isForward = intersection.position.z() > m_OffsetZ;
1476 intersection.layer = layer;
1477 intersection.sector = sector;
1478 return true;
1479 }
1480 }
1481 return false;
1482}
bool escaped
MUID: flag to indicate that the extrapolated track escaped from the KLM.
int firstBarrelLayer
MUID: outermost barrel layer encountered by the extrapolated track in the prior steps.

◆ findClosestTrackToKLMClusters()

void findClosestTrackToKLMClusters ( )
private

Find the closest Track to each KLMCluster and fill the corresponding fields of KLMCluster objects.

Definition at line 723 of file TrackExtrapolateG4e.cc.

724{
725 for (auto& klmCluster : m_klmClusters) {
726 const auto& trackClusterSeparations = klmCluster.getRelationsWith<TrackClusterSeparation>();
727
728 // If there are no TrackClusterSeparation objects related to this cluster, set NaN
729 if (trackClusterSeparations.size() == 0) {
730 klmCluster.setClusterTrackSeparation(Const::doubleNaN);
731 klmCluster.setClusterTrackSeparationAngle(Const::doubleNaN);
732 klmCluster.setClusterTrackRotationAngle(Const::doubleNaN);
733 continue;
734 }
735
736 // Look for the closest track by comparing the TrackClusterSeparation objects
737 const auto closestSeparationIterator = std::min_element(trackClusterSeparations.begin(), trackClusterSeparations.end(),
738 [](const auto & a, const auto & b) {
739 return a.getDistance() < b.getDistance();
740 });
741
742 // We found the closest track, let's set the cluster properties accordingly
743 const auto& closestSeparation = *closestSeparationIterator;
744 klmCluster.setClusterTrackSeparation(closestSeparation.getDistance() / CLHEP::cm);
745 klmCluster.setClusterTrackSeparationAngle(closestSeparation.getTrackClusterSeparationAngle());
746 klmCluster.setClusterTrackRotationAngle(closestSeparation.getTrackRotationAngle());
747 }
748}
static const double doubleNaN
quiet_NaN
Definition Const.h:703

◆ findEndcapIntersection()

bool findEndcapIntersection ( ExtState & extState,
const G4ThreeVector & oldPosition,
Intersection & intersection )
private

Find the intersection point of the track with the crossed EKLM plane.

Definition at line 1484 of file TrackExtrapolateG4e.cc.

1485{
1486 // Be generous: allow intersection to be in the dead space between m_EndcapMinR and innermost
1487 // sensitive-volume radius or between outermost sensitive-volume radius and m_EndcapMaxR,
1488 // not necessarily in a geant4 sensitive volume
1489 if (oldPosition.perp() > m_EndcapMaxR)
1490 return false;
1491 if (intersection.position.perp() < m_EndcapMinR)
1492 return false;
1493 double oldZ = std::fabs(oldPosition.z() - m_OffsetZ);
1494 double newZ = std::fabs(intersection.position.z() - m_OffsetZ);
1495 bool isForward = intersection.position.z() > m_OffsetZ;
1496 int outermostLayer = isForward ? m_OutermostActiveForwardEndcapLayer
1498 for (int layer = extState.firstEndcapLayer; layer <= outermostLayer; ++layer) {
1499 if (newZ < m_EndcapModuleMiddleZ[layer])
1500 break;
1501 if (oldZ <= m_EndcapModuleMiddleZ[layer]) {
1502 extState.firstEndcapLayer = layer + 1; // ratchet outward for next call's loop starting value
1503 if (extState.firstEndcapLayer > outermostLayer)
1504 extState.escaped = true;
1505 intersection.inBarrel = false;
1506 intersection.isForward = isForward;
1507 intersection.layer = layer;
1508 intersection.sector = m_eklmTransformData->getSectorByPosition(
1509 isForward ? 2 : 1, intersection.position) - 1;
1510 return true;
1511 }
1512 }
1513 return false;
1514}
bool isForward(int cellId)
Check whether the crystal is in forward ECL.
int firstEndcapLayer
MUID: outermost barrel layer encountered by the extrapolated track in the prior steps.

◆ findMatchingBarrelHit()

bool findMatchingBarrelHit ( Intersection & intersection,
const Track * track )
private

Find the matching BKLM 2D hit nearest the intersection point of the track with the crossed BKLM plane.

Definition at line 1516 of file TrackExtrapolateG4e.cc.

1518{
1519 G4ThreeVector extPos0(intersection.position);
1520 double diffBestMagSq = 1.0E60;
1521 int bestHit = -1;
1522 int matchingLayer = intersection.layer + 1;
1523 G4ThreeVector n(m_BarrelSectorPerp[intersection.sector]);
1524 for (int h = 0; h < m_klmHit2ds.getEntries(); ++h) {
1525 KLMHit2d* hit = m_klmHit2ds[h];
1526 if (hit->getSubdetector() != KLMElementNumbers::c_BKLM)
1527 continue;
1528 if (hit->getLayer() != matchingLayer)
1529 continue;
1530 if (hit->isOutOfTime())
1531 continue;
1532 if (std::fabs(hit->getTime() - m_MeanDt) > m_MaxDt)
1533 continue;
1534 G4ThreeVector diff(hit->getPositionX() - intersection.position.x(),
1535 hit->getPositionY() - intersection.position.y(),
1536 hit->getPositionZ() - intersection.position.z());
1537 double dn = diff * n; // in cm
1538 if (std::fabs(dn) < 2.0) {
1539 // Hit and extrapolated point are in the same sector
1540 diff -= n * dn;
1541 if (diff.mag2() < diffBestMagSq) {
1542 diffBestMagSq = diff.mag2();
1543 bestHit = h;
1544 extPos0 = intersection.position;
1545 }
1546 } else {
1547 // Accept a nearby hit in adjacent sector
1548 if (std::fabs(dn) > 50.0)
1549 continue;
1550 int sector = hit->getSector() - 1;
1551 int dSector = abs(intersection.sector - sector);
1552 if ((dSector != +1) && (dSector != m_BarrelNSector - 1))
1553 continue;
1554 // Use the normal vector of the adjacent (hit's) sector
1555 G4ThreeVector nHit(m_BarrelSectorPerp[sector]);
1556 int section = intersection.isForward ?
1559 double dn2 = intersection.position * nHit - m_BarrelModuleMiddleRadius[section][sector][intersection.layer];
1560 dn = diff * nHit + dn2;
1561 if (std::fabs(dn) > 1.0)
1562 continue;
1563 // Project extrapolated track to the hit's plane in the adjacent sector
1564 G4ThreeVector extDir(intersection.momentum.unit());
1565 double extDirA = extDir * nHit;
1566 if (std::fabs(extDirA) < 1.0E-6)
1567 continue;
1568 G4ThreeVector projection = extDir * (dn2 / extDirA);
1569 if (projection.mag() > 15.0)
1570 continue;
1571 diff += projection - nHit * dn;
1572 if (diff.mag2() < diffBestMagSq) {
1573 diffBestMagSq = diff.mag2();
1574 bestHit = h;
1575 extPos0 = intersection.position - projection;
1576 }
1577 }
1578 }
1579
1580 if (bestHit >= 0) {
1581 KLMHit2d* hit = m_klmHit2ds[bestHit];
1582 intersection.isForward = (hit->getSection() == 1);
1583 intersection.sector = hit->getSector() - 1;
1584 intersection.time = hit->getTime();
1585 double localVariance[2] = {m_BarrelScintVariance, m_BarrelScintVariance};
1586 if (hit->inRPC()) {
1587 int nStrips = hit->getPhiStripMax() - hit->getPhiStripMin() + 1;
1588 double dn = nStrips - 1.5;
1589 double factor = std::pow((0.9 + 0.4 * dn * dn), 1.5) * 0.60; // measured-in-Belle resolution
1590 localVariance[0] = m_BarrelPhiStripVariance[intersection.layer] * factor;
1591 nStrips = hit->getZStripMax() - hit->getZStripMin() + 1;
1592 dn = nStrips - 1.5;
1593 factor = std::pow((0.9 + 0.4 * dn * dn), 1.5) * 0.55; // measured-in-Belle resolution
1594 localVariance[1] = m_BarrelZStripVariance[intersection.layer] * factor;
1595 } else {
1596 int nStrips = hit->getPhiStripMax() - hit->getPhiStripMin() + 1;
1597 localVariance[0] *= (nStrips * nStrips); // variance inflated for multi-strip hit
1598 nStrips = hit->getZStripMax() - hit->getZStripMin() + 1;
1599 localVariance[1] *= (nStrips * nStrips); // variance inflated for multi-strip hit
1600 }
1601 G4ThreeVector hitPos(hit->getPositionX(), hit->getPositionY(),
1602 hit->getPositionZ());
1603 adjustIntersection(intersection, localVariance, hitPos, extPos0);
1604 if (intersection.chi2 >= 0.0) {
1605 intersection.hit = bestHit;
1606 hit->isOnTrack(true);
1607 if (track != nullptr) {
1608 track->addRelationTo(hit, intersection.chi2);
1609 RecoTrack* recoTrack = track->getRelatedTo<RecoTrack>();
1611 recoTrack->addBKLMHit(hit, recoTrack->getNumberOfTotalHits() + 1);
1612 }
1613 }
1614 }
1615 }
1616 return intersection.chi2 >= 0.0;
1617
1618}
bool addBKLMHit(const UsedBKLMHit *bklmHit, const unsigned int sortingParameter, OriginTrackFinder foundByTrackFinder=OriginTrackFinder::c_undefinedTrackFinder)
Adds a bklm hit with the given information to the reco track.
Definition RecoTrack.h:286
unsigned int getNumberOfTotalHits() const
Return the number of cdc + svd + pxd + bklm + eklm hits.
Definition RecoTrack.h:436
bool m_addHitsToRecoTrack
Parameter to add the found hits also to the reco tracks or not. Is turned off by default.
void adjustIntersection(Intersection &, const double *, const G4ThreeVector &, const G4ThreeVector &)
Nudge the track using the matching hit.

◆ findMatchingEndcapHit()

bool findMatchingEndcapHit ( Intersection & intersection,
const Track * track )
private

Find the matching EKLM 2D hit nearest the intersection point of the track with the crossed EKLM plane.

Definition at line 1620 of file TrackExtrapolateG4e.cc.

1621{
1622 double diffBestMagSq = 1.0E60;
1623 int bestHit = -1;
1624 int matchingLayer = intersection.layer + 1;
1625 int matchingEndcap = (intersection.isForward ? 2 : 1);
1626 G4ThreeVector n(0.0, 0.0, (intersection.isForward ? 1.0 : -1.0));
1627 for (int h = 0; h < m_klmHit2ds.getEntries(); ++h) {
1628 KLMHit2d* hit = m_klmHit2ds[h];
1629 if (hit->getSubdetector() != KLMElementNumbers::c_EKLM)
1630 continue;
1631 if (hit->getLayer() != matchingLayer)
1632 continue;
1633 if (hit->getSection() != matchingEndcap)
1634 continue;
1635 // DIVOT no such function for EKLM!
1636 // if (hit->isOutOfTime()) continue;
1637 if (std::fabs(hit->getTime() - m_MeanDt) > m_MaxDt)
1638 continue;
1639 G4ThreeVector diff(hit->getPositionX() - intersection.position.x(),
1640 hit->getPositionY() - intersection.position.y(),
1641 hit->getPositionZ() - intersection.position.z());
1642 double dn = diff * n; // in cm
1643 if (std::fabs(dn) > 2.0)
1644 continue;
1645 diff -= n * dn;
1646 if (diff.mag2() < diffBestMagSq) {
1647 diffBestMagSq = diff.mag2();
1648 bestHit = h;
1649 }
1650 }
1651
1652 if (bestHit >= 0) {
1653 KLMHit2d* hit = m_klmHit2ds[bestHit];
1654 intersection.hit = bestHit;
1655 intersection.isForward = (hit->getSection() == 2);
1656 intersection.sector = hit->getSector() - 1;
1657 intersection.time = hit->getTime();
1658 double localVariance[2] = {m_EndcapScintVariance, m_EndcapScintVariance};
1659 int nStrips = hit->getXStripMax() - hit->getXStripMin() + 1;
1660 localVariance[0] *= (nStrips * nStrips); // variance inflated for multi-strip hit
1661 nStrips = hit->getYStripMax() - hit->getYStripMin() + 1;
1662 localVariance[1] *= (nStrips * nStrips); // variance inflated for multi-strip hit
1663 G4ThreeVector hitPos(hit->getPositionX(), hit->getPositionY(), hit->getPositionZ());
1664 adjustIntersection(intersection, localVariance, hitPos, intersection.position);
1665 if (intersection.chi2 >= 0.0) {
1666 // DIVOT no such function for EKLM!
1667 // hit->isOnTrack(true);
1668 if (track != nullptr) {
1669 track->addRelationTo(hit, intersection.chi2);
1670 RecoTrack* recoTrack = track->getRelatedTo<RecoTrack>();
1672 for (const EKLMAlignmentHit& alignmentHit : hit->getRelationsFrom<EKLMAlignmentHit>()) {
1673 recoTrack->addEKLMHit(&alignmentHit, recoTrack->getNumberOfTotalHits() + 1);
1674 }
1675 }
1676 }
1677 }
1678 }
1679 return intersection.chi2 >= 0.0;
1680
1681}
bool addEKLMHit(const UsedEKLMHit *eklmHit, const unsigned int sortingParameter, OriginTrackFinder foundByTrackFinder=OriginTrackFinder::c_undefinedTrackFinder)
Adds an eklm hit with the given information to the reco track.
Definition RecoTrack.h:300

◆ finishTrack()

void finishTrack ( const ExtState & extState,
KLMMuidLikelihood * klmMuidLikelihood,
bool isForward )
private

Complete muon identification after end of track extrapolation.

Definition at line 1806 of file TrackExtrapolateG4e.cc.

1807{
1808 /* Done with this track: compute KLM likelihoods and fill the relative dataobject. */
1809 int lastExtLayer = extState.lastBarrelExtLayer + extState.lastEndcapExtLayer + 1;
1811 isForward, extState.escaped, extState.lastBarrelExtLayer, extState.lastEndcapExtLayer);
1812 klmMuidLikelihood->setOutcome(outcome);
1813 klmMuidLikelihood->setIsForward(isForward);
1814 klmMuidLikelihood->setBarrelExtLayer(extState.lastBarrelExtLayer);
1815 klmMuidLikelihood->setEndcapExtLayer(extState.lastEndcapExtLayer);
1816 klmMuidLikelihood->setBarrelHitLayer(extState.lastBarrelHitLayer);
1817 klmMuidLikelihood->setEndcapHitLayer(extState.lastEndcapHitLayer);
1818 klmMuidLikelihood->setExtLayer(lastExtLayer);
1819 klmMuidLikelihood->setHitLayer(((extState.lastEndcapHitLayer == -1) ?
1820 extState.lastBarrelHitLayer :
1821 extState.lastBarrelExtLayer + extState.lastEndcapHitLayer + 1));
1822 klmMuidLikelihood->setChiSquared(extState.chi2);
1823 klmMuidLikelihood->setDegreesOfFreedom(extState.nPoint);
1824 klmMuidLikelihood->setExtLayerPattern(extState.extLayerPattern);
1825 klmMuidLikelihood->setHitLayerPattern(extState.hitLayerPattern);
1826 /* Do KLM likelihood calculation. */
1827 if (outcome != MuidElementNumbers::c_NotReached) { /* Extrapolation reached KLM sensitive volume. */
1828 double denom = 0.0;
1829 int charge = klmMuidLikelihood->getCharge();
1830 std::vector<int> signedPdgVector = MuidElementNumbers::getPDGVector(charge);
1831 std::map<int, double> mapPdgPDF;
1832 for (int pdg : signedPdgVector) {
1833 auto search = m_MuidBuilderMap.find(pdg);
1834 if (search == m_MuidBuilderMap.end())
1835 B2FATAL("Something went wrong: PDF for PDG code " << pdg << " not found!");
1836 double pdf = (search->second)->getPDF(klmMuidLikelihood);
1837 denom += pdf;
1838 mapPdgPDF.insert(std::pair<int, double>(std::abs(pdg), pdf));
1839 }
1840 if (denom < 1.0E-20)
1841 klmMuidLikelihood->setJunkPDFValue(true); /* Anomaly: should be very rare. */
1842 else {
1843 for (auto const& [pdg, pdf] : mapPdgPDF) {
1844 klmMuidLikelihood->setPDFValue(pdf, std::abs(pdg));
1845 if (pdf > 0.0)
1846 klmMuidLikelihood->setLogL(std::log(pdf), std::abs(pdg));
1847 }
1848 }
1849 }
1850}
void setExtLayer(int layer)
Set the outermost EKLM layer crossed in the extrapolation.
void setExtLayerPattern(unsigned int pattern)
Set the pattern of the layers crossed in the extrapolation.
void setLogL(double logL, int pdg)
Set the log-likelihood.
void setChiSquared(double chiSquared)
Set the chi-squared of the extrapolation.
void setHitLayerPattern(unsigned int pattern)
Set the pattern of the layers actually crossed by the track.
void setBarrelExtLayer(int layer)
Set the outermost BKLM layer crossed in the extrapolation.
void setIsForward(bool isForward)
Set if this extrapolation is in forward or backward B/EKLM.
void setJunkPDFValue(bool flag)
Set the junk flag (1 if junk, 0 if not).
void setDegreesOfFreedom(int dof)
Set the number of degrees of freedom (= 2 times the number of KLM hits) for the chi-square computatio...
int getCharge() const
Get the charge of the particle hypothesis used during the extrapolation.
void setPDFValue(double pdfValue, int pdg)
Set the normalized PDF.
void setEndcapExtLayer(int layer)
Set the outermost EKLM layer crossed in the extrapolation.
void setBarrelHitLayer(int layer)
Set the outermost BKLM layer actually crossed by the track.
void setOutcome(unsigned int outcome)
Set the outcome of this extrapolation.
void setHitLayer(int layer)
Set the outermost KLM layer actually crossed by the track.
void setEndcapHitLayer(int layer)
Set the outermost EKLM layer actually crossed by the track.
static unsigned int calculateExtrapolationOutcome(bool isForward, bool escaped, int lastBarrelLayer, int lastEndcapLayer)
Calculate the track extrapolation outcome.
double charge(int pdgCode)
Returns electric charge of a particle with given pdg code.
Definition EvtPDLUtil.cc:44

◆ fromG4eToPhasespace()

void fromG4eToPhasespace ( const G4ErrorFreeTrajState & g4eState,
G4ErrorSymMatrix & covariance )
private

Convert the geant4e 5x5 covariance to phasespace 6x6 covariance.

Definition at line 1075 of file TrackExtrapolateG4e.cc.

1076{
1077
1078 // Convert Geant4e covariance matrix with parameters 1/p, lambda, phi, yT, zT (in GeV/c, radians, cm)
1079 // to phase-space covariance matrix with parameters x, y, z, px, py, pz (in GeV/c, cm)
1080 // (1/p) = 1/sqrt( px^2 + py^2 + pz^2 )
1081 // phi = atan( py / px )
1082 // lambda = asin( pz / sqrt( px^2 + py^2 + pz^2 )
1083 // xT = x * cos(lambda) * cos(phi) + y * cos(lambda) * sin(phi) + z * sin(lambda)
1084 // yT = -x * sin(phi) + y * cos(phi)
1085 // zT = -x * sin(lambda) * cos(phi) - y * sin(lambda) * sin(phi) + z * cos(lambda)
1086
1087 G4ErrorFreeTrajParam param = g4eState.GetParameters();
1088 double p = 1.0 / (param.GetInvP() * CLHEP::GeV); // in GeV/c
1089 double pSq = p * p;
1090 double lambda = param.GetLambda(); // in radians
1091 double sinLambda = std::sin(lambda);
1092 double cosLambda = std::cos(lambda);
1093 double phi = param.GetPhi(); // in radians
1094 double sinPhi = std::sin(phi);
1095 double cosPhi = std::cos(phi);
1096
1097 // Transformation Jacobian 6x5 from Geant4e 5x5 to phase-space 6x6
1098
1099 G4ErrorMatrix jacobian(6, 5, 0); // All entries are initialized to 0
1100
1101 jacobian(4, 1) = -pSq * cosLambda * cosPhi; // @(px)/@(1/p)
1102 jacobian(5, 1) = -pSq * cosLambda * sinPhi; // @(py)/@(1/p)
1103 jacobian(6, 1) = -pSq * sinLambda; // @(pz)/@(1/p)
1104
1105 jacobian(4, 2) = -p * sinLambda * cosPhi; // @(px)/@(lambda)
1106 jacobian(5, 2) = -p * sinLambda * sinPhi; // @(py)/@(lambda)
1107 jacobian(6, 2) = p * cosLambda; // @(pz)/@(lambda)
1108
1109 jacobian(4, 3) = -p * cosLambda * sinPhi; // @(px)/@(phi)
1110 jacobian(5, 3) = p * cosLambda * cosPhi; // @(py)/@(phi)
1111
1112 jacobian(1, 4) = -sinPhi; // @(x)/@(yT)
1113 jacobian(2, 4) = cosPhi; // @(y)/@(yT)
1114
1115 jacobian(1, 5) = -sinLambda * cosPhi; // @(x)/@(zT)
1116 jacobian(2, 5) = -sinLambda * sinPhi; // @(y)/@(zT)
1117 jacobian(3, 5) = cosLambda; // @(z)/@(zT)
1118
1119 G4ErrorTrajErr g4eCov = g4eState.GetError();
1120 covariance.assign(g4eCov.similarity(jacobian));
1121
1122}

◆ fromPhasespaceToG4e() [1/2]

void fromPhasespaceToG4e ( const G4ThreeVector & momentum,
const G4ErrorSymMatrix & covariance,
G4ErrorTrajErr & covG4e )
private

Convert the phasespace covariance to geant4e covariance.

Definition at line 1177 of file TrackExtrapolateG4e.cc.

1179{
1180
1181 // Convert phase-space covariance matrix with parameters x, y, z, px, py, pz (in genfit2 units: cm, GeV/c)
1182 // to Geant4e covariance matrix with parameters 1/p, lambda, phi, yT, zT (in geant4e units: GeV/c, radians, cm)
1183 // xT = x * cos(lambda) * cos(phi) + y * cos(lambda) * sin(phi) + z * sin(lambda)
1184 // yT = -x * sin(phi) + y * cos(phi)
1185 // zT = -x * sin(lambda) * cos(phi) - y * sin(lambda) * sin(phi) + z * cos(lambda)
1186 // (1/p) = 1/sqrt( px^2 + py^2 + pz^2 )
1187 // phi = atan( py / px )
1188 // lambda = asin( pz / sqrt( px^2 + py^2 + pz^2 )
1189
1190 G4ErrorSymMatrix temp(covariance);
1191
1192 double pInvSq = 1.0 / momentum.mag2();
1193 double pInv = std::sqrt(pInvSq);
1194 double pPerpInv = 1.0 / momentum.perp();
1195 double sinLambda = momentum.cosTheta();
1196 double cosLambda = std::sqrt(1.0 - sinLambda * sinLambda);
1197 double phi = momentum.phi();
1198 double cosPhi = std::cos(phi);
1199 double sinPhi = std::sin(phi);
1200
1201 // Transformation Jacobian 5x6 from phase-space 6x6 to Geant4e 5x5
1202 G4ErrorMatrix jacobian(5, 6, 0);
1203
1204 jacobian(1, 4) = -pInvSq * cosLambda * cosPhi; // @(1/p)/@(px)
1205 jacobian(1, 5) = -pInvSq * cosLambda * sinPhi; // @(1/p)/@(py)
1206 jacobian(1, 6) = -pInvSq * sinLambda; // @(1/p)/@(pz)
1207
1208 jacobian(2, 4) = -pInv * sinLambda * cosPhi; // @(lambda)/@(px)
1209 jacobian(2, 5) = -pInv * sinLambda * sinPhi; // @(lambda)/@(py)
1210 jacobian(2, 6) = pInv * cosLambda; // @(lambda)/@(pz)
1211
1212 jacobian(3, 4) = -pPerpInv * sinPhi; // @(phi)/@(px)
1213 jacobian(3, 5) = pPerpInv * cosPhi; // @(phi)/@(py)
1214
1215 jacobian(4, 1) = -sinPhi; // @(yT)/@(x)
1216 jacobian(4, 2) = cosPhi; // @(yT)/@(y)
1217
1218 jacobian(5, 1) = -sinLambda * cosPhi; // @(zT)/@(x)
1219 jacobian(5, 2) = -sinLambda * sinPhi; // @(zT)/@(y)
1220 jacobian(5, 3) = cosLambda; // @(zT)/@(z)
1221 covG4e = temp.similarity(jacobian);
1222
1223}

◆ fromPhasespaceToG4e() [2/2]

void fromPhasespaceToG4e ( const TVector3 & momentum,
const TMatrixDSym & covariance,
G4ErrorTrajErr & covG4e )
private

Convert the phasespace covariance to geant4e covariance.

Definition at line 1124 of file TrackExtrapolateG4e.cc.

1125{
1126
1127 // Convert phase-space covariance matrix with parameters x, y, z, px, py, pz (in genfit2 units: cm, GeV/c)
1128 // to Geant4e covariance matrix with parameters 1/p, lambda, phi, yT, zT (in geant4e units: GeV/c, radians, cm)
1129 // xT = x * cos(lambda) * cos(phi) + y * cos(lambda) * sin(phi) + z * sin(lambda)
1130 // yT = -x * sin(phi) + y * cos(phi)
1131 // zT = -x * sin(lambda) * cos(phi) - y * sin(lambda) * sin(phi) + z * cos(lambda)
1132 // (1/p) = 1/sqrt( px^2 + py^2 + pz^2 )
1133 // phi = atan( py / px )
1134 // lambda = asin( pz / sqrt( px^2 + py^2 + pz^2 )
1135
1136 G4ErrorSymMatrix temp(6, 0);
1137 for (int k = 0; k < 6; ++k) {
1138 for (int j = k; j < 6; ++j) {
1139 temp[j][k] = covariance[j][k];
1140 }
1141 }
1142
1143 double pInvSq = 1.0 / momentum.Mag2();
1144 double pInv = std::sqrt(pInvSq);
1145 double pPerpInv = 1.0 / momentum.Perp();
1146 double sinLambda = momentum.CosTheta();
1147 double cosLambda = std::sqrt(1.0 - sinLambda * sinLambda);
1148 double phi = momentum.Phi();
1149 double cosPhi = std::cos(phi);
1150 double sinPhi = std::sin(phi);
1151
1152 // Transformation Jacobian 5x6 from phase-space 6x6 to Geant4e 5x5
1153 G4ErrorMatrix jacobian(5, 6, 0); // All entries are initialized to 0
1154
1155 jacobian(1, 4) = -pInvSq * cosLambda * cosPhi; // @(1/p)/@(px)
1156 jacobian(1, 5) = -pInvSq * cosLambda * sinPhi; // @(1/p)/@(py)
1157 jacobian(1, 6) = -pInvSq * sinLambda; // @(1/p)/@(pz)
1158
1159 jacobian(2, 4) = -pInv * sinLambda * cosPhi; // @(lambda)/@(px)
1160 jacobian(2, 5) = -pInv * sinLambda * sinPhi; // @(lambda)/@(py)
1161 jacobian(2, 6) = pInv * cosLambda; // @(lambda)/@(pz)
1162
1163 jacobian(3, 4) = -pPerpInv * sinPhi; // @(phi)/@(px)
1164 jacobian(3, 5) = pPerpInv * cosPhi; // @(phi)/@(py)
1165
1166 jacobian(4, 1) = -sinPhi; // @(yT)/@(x)
1167 jacobian(4, 2) = cosPhi; // @(yT)/@(y)
1168
1169 jacobian(5, 1) = -sinLambda * cosPhi; // @(zT)/@(x)
1170 jacobian(5, 2) = -sinLambda * sinPhi; // @(zT)/@(y)
1171 jacobian(5, 3) = cosLambda; // @(zT)/@(z)
1172
1173 covG4e = temp.similarity(jacobian);
1174
1175}

◆ getInstance()

TrackExtrapolateG4e * getInstance ( )
static

Get the singleton's address.

Definition at line 85 of file TrackExtrapolateG4e.cc.

86{
87 if (m_Singleton == nullptr)
89 return m_Singleton;
90}
static TrackExtrapolateG4e * m_Singleton
Stores pointer to the singleton class.
TrackExtrapolateG4e()
constructor is hidden; user calls TrackExtrapolateG4e::getInstance() instead

◆ getStartPoint()

ExtState getStartPoint ( const Track & b2track,
int pdgCode,
G4ErrorFreeTrajState & g4eState )
private

Get the start point for a new reconstructed track with specific PDG hypothesis.

Definition at line 982 of file TrackExtrapolateG4e.cc.

983{
984 ExtState extState = {&b2track, pdgCode, false, 0.0, 0.0, // for EXT and MUID
985 G4ThreeVector(0, 0, 1), 0.0, 0, 0, 0, -1, -1, -1, -1, 0, 0, false // for MUID only
986 };
987 RecoTrack* recoTrack = b2track.getRelatedTo<RecoTrack>();
988 if (recoTrack == nullptr) {
989 B2WARNING("Track without associated RecoTrack: skipping extrapolation for this track.");
990 extState.pdgCode = 0; // prevent start of extrapolation in swim()
991 return extState;
992 }
993 const genfit::AbsTrackRep* trackRep = recoTrack->getCardinalRepresentation();
994 // check for a valid track fit
995 if (!recoTrack->wasFitSuccessful(trackRep)) {
996 B2WARNING("RecoTrack fit failed for cardinal representation: skipping extrapolation for this track.");
997 extState.pdgCode = 0; // prevent start of extrapolation in swim()
998 return extState;
999 }
1000 int charge = int(trackRep->getPDGCharge());
1001 if (charge != 0) {
1002 extState.pdgCode *= charge;
1003 } else {
1004 charge = 1; // should never happen but persist if it does
1005 }
1006 TVector3 firstPosition, firstMomentum, lastPosition, lastMomentum; // initialized to zeroes
1007 TMatrixDSym firstCov(6), lastCov(6); // initialized to zeroes
1008 try {
1009 const genfit::MeasuredStateOnPlane& firstState = recoTrack->getMeasuredStateOnPlaneFromFirstHit(trackRep);
1010 trackRep->getPosMomCov(firstState, firstPosition, firstMomentum, firstCov);
1011 const genfit::MeasuredStateOnPlane& lastState = recoTrack->getMeasuredStateOnPlaneFromLastHit(trackRep);
1012 trackRep->getPosMomCov(lastState, lastPosition, lastMomentum, lastCov);
1013 // in genfit units (cm, GeV/c)
1014 extState.tof = lastState.getTime(); // DIVOT: must be revised when IP profile (reconstructed beam spot) become available!
1015 if (lastPosition.Mag2() < firstPosition.Mag2()) {
1016 firstPosition = lastPosition;
1017 firstMomentum = -lastMomentum;
1018 firstCov = lastCov;
1019 trackRep->getPosMomCov(firstState, lastPosition, lastMomentum, lastCov);
1020 lastMomentum *= -1.0; // extrapolate backwards instead of forwards
1021 extState.isCosmic = true;
1022 extState.tof = firstState.getTime(); // DIVOT: must be revised when IP profile (reconstructed beam spot) become available!
1023 }
1024
1025 G4ParticleDefinition* particle = G4ParticleTable::GetParticleTable()->FindParticle(extState.pdgCode);
1026 if (extState.pdgCode != trackRep->getPDG()) {
1027 double pSq = lastMomentum.Mag2();
1028 double mass = particle->GetPDGMass() / CLHEP::GeV;
1029 extState.tof *= std::sqrt((pSq + mass * mass) / (pSq + lastState.getMass() * lastState.getMass()));
1030 }
1031
1032 extState.directionAtIP.set(firstMomentum.Unit().X(), firstMomentum.Unit().Y(), firstMomentum.Unit().Z());
1033 if (m_MagneticField != 0.0) { // in gauss
1034 double radius = (firstMomentum.Perp() * CLHEP::GeV / CLHEP::eV) / (CLHEP::c_light * charge * m_MagneticField); // in cm
1035 double centerPhi = extState.directionAtIP.phi() - M_PI_2;
1036 double centerX = firstPosition.X() + radius * std::cos(centerPhi);
1037 double centerY = firstPosition.Y() + radius * std::sin(centerPhi);
1038 double pocaPhi = atan2(charge * centerY, charge * centerX) + M_PI;
1039 double ipPerp = extState.directionAtIP.perp();
1040 if (ipPerp > 0.0) {
1041 extState.directionAtIP.setX(+std::sin(pocaPhi) * ipPerp);
1042 extState.directionAtIP.setY(-std::cos(pocaPhi) * ipPerp);
1043 }
1044 } else {
1045 // No field: replace flaky covariance matrix with a diagonal one measured in 1.5T field
1046 // for a 10 GeV/c track ... and replace momentum magnitude with fixed 10 GeV/c
1047 lastCov *= 0.0;
1048 lastCov[0][0] = 5.0E-5;
1049 lastCov[1][1] = 1.0E-7;
1050 lastCov[2][2] = 5.0E-4;
1051 lastCov[3][3] = 3.5E-3;
1052 lastCov[4][4] = 3.5E-3;
1053 lastMomentum = lastMomentum.Unit() * 10.0;
1054 }
1055
1056 G4ThreeVector posG4e(lastPosition.X() * CLHEP::cm, lastPosition.Y() * CLHEP::cm,
1057 lastPosition.Z() * CLHEP::cm); // in Geant4 units (mm)
1058 G4ThreeVector momG4e(lastMomentum.X() * CLHEP::GeV, lastMomentum.Y() * CLHEP::GeV,
1059 lastMomentum.Z() * CLHEP::GeV); // in Geant4 units (MeV/c)
1060 G4ErrorSymMatrix covG4e(5, 0); // in Geant4e units (GeV/c, cm)
1061 fromPhasespaceToG4e(lastMomentum, lastCov, covG4e); // in Geant4e units (GeV/c, cm)
1062 g4eState.SetData("g4e_" + particle->GetParticleName(), posG4e, momG4e);
1063 g4eState.SetParameters(posG4e, momG4e); // compute private-state parameters from momG4e
1064 g4eState.SetError(covG4e);
1065 } catch (const genfit::Exception&) {
1066 B2WARNING("genfit::MeasuredStateOnPlane() exception: skipping extrapolation for this track. initial momentum = ("
1067 << firstMomentum.X() << "," << firstMomentum.Y() << "," << firstMomentum.Z() << ")");
1068 extState.pdgCode = 0; // prevent start of extrapolation in swim()
1069 }
1070
1071 return extState;
1072}
bool wasFitSuccessful(const genfit::AbsTrackRep *representation=nullptr) const
Returns true if the last fit with the given representation was successful.
Definition RecoTrack.cc:336
genfit::AbsTrackRep * getCardinalRepresentation() const
Get a pointer to the cardinal track representation. You are not allowed to modify or delete it!
Definition RecoTrack.h:631
const genfit::MeasuredStateOnPlane & getMeasuredStateOnPlaneFromLastHit(const genfit::AbsTrackRep *representation=nullptr) const
Return genfit's MeasuredStateOnPlane for the last hit in a fit useful for extrapolation of measuremen...
Definition RecoTrack.cc:619
const genfit::MeasuredStateOnPlane & getMeasuredStateOnPlaneFromFirstHit(const genfit::AbsTrackRep *representation=nullptr) const
Return genfit's MeasuredStateOnPlane for the first hit in a fit useful for extrapolation of measureme...
Definition RecoTrack.cc:605
TO * getRelatedTo(const std::string &name="", const std::string &namedRelation="") const
Get the object to which this object has a relation.
G4ThreeVector directionAtIP
MUID: initial direction of track, used for KLID.

◆ getVolumeID()

void getVolumeID ( const G4TouchableHandle & touch,
Const::EDetector & detID,
int & copyID )
private

Get the physical volume information for a geant4 physical volume.

Definition at line 896 of file TrackExtrapolateG4e.cc.

897{
898
899 // default values
900 detID = Const::EDetector::invalidDetector;
901 copyID = 0;
902
903 G4VPhysicalVolume* pv = touch->GetVolume(0);
904 std::map<G4VPhysicalVolume*, enum VolTypes>::iterator it = m_EnterExit->find(pv);
905 if (it == m_EnterExit->end())
906 return;
907
908 switch (it->second) {
909 case VOLTYPE_CDC:
910 detID = Const::EDetector::CDC;
911 copyID = pv->GetCopyNo();
912 return;
913 case VOLTYPE_TOP1:
914 detID = Const::EDetector::TOP;
915 copyID = -(pv->GetCopyNo()); // negative to distinguish module and quartz hits
916 return;
917 case VOLTYPE_TOP2:
918 detID = Const::EDetector::TOP;
919 if (touch->GetHistoryDepth() >= 1)
920 copyID = touch->GetVolume(1)->GetCopyNo();
921 return;
922 case VOLTYPE_TOP3:
923 detID = Const::EDetector::TOP;
924 if (touch->GetHistoryDepth() >= 2)
925 copyID = touch->GetVolume(2)->GetCopyNo();
926 return;
927 case VOLTYPE_ARICH1:
928 detID = Const::EDetector::ARICH;
929 copyID = 12345;
930 return;
931 case VOLTYPE_ARICH2:
932 detID = Const::EDetector::ARICH;
933 copyID = 6789;
934 return;
935 case VOLTYPE_ARICH3:
936 detID = Const::EDetector::ARICH;
937 if (touch->GetHistoryDepth() >= 2)
938 copyID = touch->GetVolume(2)->GetCopyNo();
939 return;
940 case VOLTYPE_ECL:
941 detID = Const::EDetector::ECL;
943 return;
944 case VOLTYPE_BKLM1: // BKLM RPCs
945 detID = Const::EDetector::BKLM;
946 if (touch->GetHistoryDepth() == DEPTH_RPC) {
947 // int plane = touch->GetCopyNumber(0);
948 int layer = touch->GetCopyNumber(4);
949 int sector = touch->GetCopyNumber(6);
950 int section = touch->GetCopyNumber(7);
951 copyID = BKLMElementNumbers::moduleNumber(section, sector, layer);
952 }
953 return;
954 case VOLTYPE_BKLM2: // BKLM scints
955 detID = Const::EDetector::BKLM;
956 if (touch->GetHistoryDepth() == DEPTH_SCINT) {
957 int strip = touch->GetCopyNumber(1);
958 int plane = (touch->GetCopyNumber(2) == BKLM_INNER) ?
961 int layer = touch->GetCopyNumber(6);
962 int sector = touch->GetCopyNumber(8);
963 int section = touch->GetCopyNumber(9);
965 section, sector, layer, plane, strip);
966 BKLMStatus::setMaximalStrip(copyID, strip);
967 }
968 return;
969 case VOLTYPE_EKLM:
970 detID = Const::EDetector::EKLM;
972 touch->GetVolume(7)->GetCopyNo(),
973 touch->GetVolume(6)->GetCopyNo(),
974 touch->GetVolume(5)->GetCopyNo(),
975 touch->GetVolume(4)->GetCopyNo(),
976 touch->GetVolume(1)->GetCopyNo());
977 return;
978 }
979
980}
static KLMChannelNumber channelNumber(int section, int sector, int layer, int plane, int strip)
Get channel number.
static KLMModuleNumber moduleNumber(int section, int sector, int layer, bool fatalError=true)
Get module number.
static void setMaximalStrip(int &module, int strip)
Set maximal strip number.
Definition BKLMStatus.h:49
static ECLGeometryPar * Instance()
Static method to get a reference to the ECLGeometryPar instance.
int ECLVolumeToCellID(const G4VTouchable *)
Get Cell Id (LEP: new way)
int stripNumber(int section, int layer, int sector, int plane, int strip) const
Get strip number.
@ VOLTYPE_ARICH2
ARICH Img plate.
@ VOLTYPE_TOP2
TOP quartz.
@ VOLTYPE_TOP3
TOP glue.
@ VOLTYPE_ARICH3
ARICH HAPD window.
@ VOLTYPE_BKLM2
BKLM scintillator.
@ VOLTYPE_EKLM
EKLM.
@ VOLTYPE_BKLM1
BKLM RPC.
@ VOLTYPE_ARICH1
ARICH aerogel.
@ VOLTYPE_TOP1
TOP container.

◆ initialize() [1/2]

void initialize ( double meanDt,
double maxDt,
double maxSeparation,
double maxKLMTrackClusterDistance,
double maxECLTrackClusterDistance,
double minPt,
double minKE,
bool addHitsToRecoTrack,
std::vector< Const::ChargedStable > & hypotheses )

Initialize for track extrapolation by the MUID module.

Parameters
meanDtMean value of the in-time window (ns).
maxDtHalf-width of the in-time window (ns).
maxSeparationMaximum separation between track crossing and matching hit in detector plane (number of sigmas).
maxKLMTrackClusterDistanceMaximum distance between associated track and KLMCluster (cm), criterion for matching relation Track->KLMCluster on MDST.
maxECLTrackClusterDistanceMaximum distance between associated track and ECLCluster (cm).
minPtMinimum transverse momentum to begin extrapolation (GeV/c).
minKEMinimum kinetic energy to continue extrapolation (GeV/c).
addHitsToRecoTrackParameter to add the found hits also to the reco tracks or not. Is turned off by default.
hypothesesVector of charged-particle hypotheses used in extrapolation of each track.

Definition at line 204 of file TrackExtrapolateG4e.cc.

208{
209 m_MuidInitialized = true;
210 m_addHitsToRecoTrack = addHitsToRecoTrack;
211
212 // Define required objects, register the new ones and relations
213 m_eclClusters.isRequired();
214 m_klmHit2ds.isRequired();
215 m_klmClusters.isRequired();
216 m_recoTracks.isRequired();
217 m_tracks.isRequired();
218 m_extHits.registerInDataStore();
219 m_klmMuidLikelihoods.registerInDataStore();
220 m_klmMuidHits.registerInDataStore();
221 m_trackClusterSeparations.registerInDataStore();
222 m_tracks.registerRelationTo(m_extHits);
223 m_tracks.registerRelationTo(m_klmMuidLikelihoods);
224 m_tracks.registerRelationTo(m_klmMuidHits);
225 m_tracks.registerRelationTo(m_klmHit2ds);
226 m_tracks.registerRelationTo(m_trackClusterSeparations);
227 m_tracks.registerRelationTo(m_klmClusters);
228 m_klmClusters.registerRelationTo(m_trackClusterSeparations);
229 m_eclClusters.registerRelationTo(m_extHits);
231
232 // Save the in-time cut's central value and width for valid hits
233 m_MeanDt = meanDt;
234 m_MaxDt = maxDt;
235
236 // Save the magnetic field z component (gauss) at the origin
237 m_MagneticField = BFieldManager::getField(0, 0, 0).Z() / Unit::T * CLHEP::tesla / CLHEP::gauss;
238
239 // Convert from sigma to variance for hit-position uncertainty
240 m_MaxDistSqInVariances = maxKLMTrackHitDistance * maxKLMTrackHitDistance;
241
242 // Convert user cutoff values to geant4 units
243 m_MaxKLMTrackClusterDistance = std::max(0.0, maxKLMTrackClusterDistance) * CLHEP::cm; // mm
244 m_MaxECLTrackClusterDistance = std::max(0.0, maxECLTrackClusterDistance) * CLHEP::cm; // mm
245 m_MinPt = std::max(0.0, minPt) * CLHEP::GeV;
246 m_MinKE = std::max(0.0, minKE) * CLHEP::GeV;
247
248 // Save pointer to the list of particle hypotheses for EXT extrapolation
249 m_HypothesesMuid = &hypotheses;
250
251 // Define the list of volumes that will have their entry and/or
252 // exit points stored during the extrapolation.
254
255 // Store the address of the ExtManager (used later)
257
258 // Set up the EXT-specific geometry (might have already been done by EXT)
259 if (m_TargetExt == nullptr) {
260 if (!m_COILGeometryPar.isValid())
261 B2FATAL("Coil geometry data are not available.");
262 double offsetZ = m_COILGeometryPar->getGlobalOffsetZ();
263 double rMinCoil = m_COILGeometryPar->getCryoRmin();
264 double halfLength = m_COILGeometryPar->getCryoLength();
265 m_COILGeometryPar.addCallback([this, &offsetZ, &rMinCoil, &halfLength]() {
266 offsetZ = m_COILGeometryPar->getGlobalOffsetZ();
267 rMinCoil = m_COILGeometryPar->getCryoRmin();
268 halfLength = m_COILGeometryPar->getCryoLength();
269 });
270 m_TargetExt = new Simulation::ExtCylSurfaceTarget(rMinCoil, offsetZ - halfLength, offsetZ + halfLength);
271 G4ErrorPropagatorData::GetErrorPropagatorData()->SetTarget(m_TargetExt);
272 }
273 if (!m_BeamPipeGeo.isValid())
274 B2FATAL("Beam pipe geometry data are not available.");
275 double beampipeRadius = m_BeamPipeGeo->getParameter("Lv2OutBe.R2") * CLHEP::cm; // mm
276 m_BeamPipeGeo.addCallback([this, &beampipeRadius]() {
277 beampipeRadius = m_BeamPipeGeo->getParameter("Lv2OutBe.R2") * CLHEP::cm;
278 });
279 m_MinRadiusSq = beampipeRadius * beampipeRadius; // mm^2
280
281 // Set up the MUID-specific geometry
282 bklm::GeometryPar* bklmGeometry = bklm::GeometryPar::instance();
283 const EKLM::GeometryData& eklmGeometry = EKLM::GeometryData::Instance();
284 m_BarrelHalfLength = bklmGeometry->getHalfLength() * CLHEP::cm; // in G4 units (mm)
285 m_EndcapHalfLength = 0.5 * eklmGeometry.getSectionPosition()->getLength(); // in G4 units (mm)
286 m_OffsetZ = bklmGeometry->getOffsetZ() * CLHEP::cm; // in G4 units (mm)
287 double minZ = m_OffsetZ - (m_BarrelHalfLength + 2.0 * m_EndcapHalfLength);
288 double maxZ = m_OffsetZ + (m_BarrelHalfLength + 2.0 * m_EndcapHalfLength);
289 m_BarrelNSector = bklmGeometry->getNSector();
290 m_BarrelMaxR = bklmGeometry->getOuterRadius() * CLHEP::cm / std::cos(M_PI / m_BarrelNSector); // in G4 units (mm)
291 m_TargetMuid = new Simulation::ExtCylSurfaceTarget(m_BarrelMaxR, minZ, maxZ);
292 G4ErrorPropagatorData::GetErrorPropagatorData()->SetTarget(m_TargetMuid);
293
294 m_BarrelHalfLength /= CLHEP::cm; // now in G4e units (cm)
295 m_EndcapHalfLength /= CLHEP::cm; // now in G4e units (cm)
296 m_OffsetZ /= CLHEP::cm; // now in G4e units (cm)
297 m_BarrelMinR = bklmGeometry->getGap1InnerRadius(); // in G4e units (cm)
298 m_BarrelMaxR /= CLHEP::cm; // now in G4e units (cm)
299 m_EndcapMinR = eklmGeometry.getSectionPosition()->getInnerR() / CLHEP::cm; // in G4e units (cm)
300 m_EndcapMaxR = eklmGeometry.getSectionPosition()->getOuterR() / CLHEP::cm; // in G4e units (cm)
302
303 // Measurement uncertainties and acceptance windows
304 double width = eklmGeometry.getStripGeometry()->getWidth() / CLHEP::cm; // in G4e units (cm)
305 m_EndcapScintVariance = width * width / 12.0;
306 width = bklmGeometry->getScintHalfWidth() * 2.0; // in G4e units (cm)
307 m_BarrelScintVariance = width * width / 12.0;
308 for (int layer = 1; layer <= BKLMElementNumbers::getMaximalLayerNumber(); ++layer) {
309 const bklm::Module* module =
310 bklmGeometry->findModule(BKLMElementNumbers::c_ForwardSection, 1, layer);
311 width = module->getPhiStripWidth(); // in G4e units (cm)
312 m_BarrelPhiStripVariance[layer - 1] = width * width / 12.0;
313 width = module->getZStripWidth(); // in G4e units (cm)
314 m_BarrelZStripVariance[layer - 1] = width * width / 12.0;
315 }
316
317 // KLM geometry (for associating KLM hit with extrapolated crossing point)
319 for (int sector = 1; sector <= BKLMElementNumbers::getMaximalSectorNumber(); ++sector) {
320 double phi = M_PI_4 * (sector - 1);
321 m_BarrelSectorPerp[sector - 1].set(std::cos(phi), std::sin(phi), 0.0);
322 m_BarrelSectorPhi[sector - 1].set(-std::sin(phi), std::cos(phi), 0.0);
323 }
324 KLMChannelIndex klmLayers(KLMChannelIndex::c_IndexLevelLayer);
325 for (KLMChannelIndex& klmLayer : klmLayers) {
326 if (klmLayer.getSubdetector() == KLMElementNumbers::c_BKLM)
327 m_BarrelModuleMiddleRadius[klmLayer.getSection()][klmLayer.getSector() - 1][klmLayer.getLayer() - 1] =
328 bklmGeometry->getActiveMiddleRadius(klmLayer.getSection(), klmLayer.getSector(), klmLayer.getLayer()); // in G4e units (cm)
329 }
330 double dz(eklmGeometry.getLayerShiftZ() / CLHEP::cm); // in G4e units (cm)
331 double z0((eklmGeometry.getSectionPosition()->getZ()
332 + eklmGeometry.getLayerShiftZ()
333 - 0.5 * eklmGeometry.getSectionPosition()->getLength()
334 - 0.5 * eklmGeometry.getLayerPosition()->getLength()) / CLHEP::cm); // in G4e units (cm)
336 1; // zero-based counting
338 1; // zero-based counting
339 for (int layer = 1; layer <= EKLMElementNumbers::getMaximalLayerNumber(); ++layer) {
340 m_EndcapModuleMiddleZ[layer - 1] = z0 + dz * (layer - 1); // in G4e units (cm)
341 }
342
344}
static constexpr int getMaximalLayerNumber()
Get maximal layer number.
double getOuterR() const
Get outer radius.
double getInnerR() const
Get inner radius.
double getZ() const
Get Z coordinate.
double getLength() const
Get length.
double getWidth() const
Get width.
const ElementPosition * getLayerPosition() const
Get position data for layers.
double getLayerShiftZ() const
Get Z distance between two layers.
const StripGeometry * getStripGeometry() const
Get strip geometry data.
const ElementPosition * getSectionPosition() const
Get position data for sections.
static const GeometryData & Instance(enum DataSource dataSource=c_Database, const GearDir *gearDir=nullptr)
Instantiation.
static const TransformDataGlobalAligned & Instance()
Instantiation.
static void registerRequiredRelations(StoreArray< RecoTrack > &recoTracks, std::string const &pxdHitsStoreArrayName="", std::string const &svdHitsStoreArrayName="", std::string const &cdcHitsStoreArrayName="", std::string const &bklmHitsStoreArrayName="", std::string const &eklmHitsStoreArrayName="", std::string const &recoHitInformationStoreArrayName="")
Convenience method which registers all relations required to fully use a RecoTrack.
Definition RecoTrack.cc:53
StoreArray< KLMMuidLikelihood > m_klmMuidLikelihoods
KLM muid likelihoods.
StoreArray< TrackClusterSeparation > m_trackClusterSeparations
Track cluster separation.
StoreArray< RecoTrack > m_recoTracks
Reco tracks.
void registerVolumes()
Register the list of geant4 physical volumes whose entry/exit points will be saved during extrapolati...
DBObjPtr< COILGeometryPar > m_COILGeometryPar
Conditions-database object for COIL geometry.
DBObjPtr< BeamPipeGeo > m_BeamPipeGeo
Conditions-database object for beam pipe geometry.
static const double T
[tesla]
Definition Unit.h:120
double getScintHalfWidth(void) const
Get the height of the entire volume of a scintillator strip (including TiO2 coating)
double getGap1InnerRadius(void) const
Get the radius of the inner tangent circle of gap 0 (innermost)
double getOuterRadius(void) const
Get the radius of the inscribed circle of the outer polygon.
double getOffsetZ(void) const
Get the global shift along a of the entire BKLM.
int getNSector(void) const
Get the number of sectors of the BKLM.
double getHalfLength(void) const
Get the half-length along z of the BKLM.
double getActiveMiddleRadius(int section, int sector, int layer) const
Get the radial midpoint of the detector module's active volume of specified layer.
static void getField(const double *pos, double *field)
return the magnetic field at a given position.

◆ initialize() [2/2]

void initialize ( double minPt,
double minKE,
std::vector< Const::ChargedStable > & hypotheses )

Initialize for track extrapolation by the EXT module.

Parameters
minPtMinimum transverse momentum to begin extrapolation (GeV/c).
minKEMinimum kinetic energy to continue extrapolation (GeV/c).
hypothesesVector of charged-particle hypotheses used in extrapolation of each track.

Definition at line 151 of file TrackExtrapolateG4e.cc.

153{
154 m_ExtInitialized = true;
155
156 // Define required objects, register the new ones and relations
157 m_recoTracks.isRequired();
158 m_tracks.isRequired();
159 m_extHits.registerInDataStore();
160 m_tracks.registerRelationTo(m_extHits);
161
162 // Save the magnetic field z component (gauss) at the origin
163 m_MagneticField = BFieldManager::getField(0, 0, 0).Z() / Unit::T * CLHEP::tesla / CLHEP::gauss;
164
165 // Convert user cutoff values to geant4 units
166 m_MinPt = std::max(0.0, minPt) * CLHEP::GeV;
167 m_MinKE = std::max(0.0, minKE) * CLHEP::GeV;
168
169 // Save pointer to the list of particle hypotheses for EXT extrapolation
170 m_HypothesesExt = &hypotheses;
171
172 // Define the list of volumes that will have their entry and/or
173 // exit points stored during the extrapolation.
175
176 // Store the address of the ExtManager (used later)
178
179 // Set up the EXT-specific geometry (might have already been done by MUID)
180 if (m_TargetExt == nullptr) {
181 if (!m_COILGeometryPar.isValid())
182 B2FATAL("Coil geometry data are not available.");
183 double offsetZ = m_COILGeometryPar->getGlobalOffsetZ();
184 double rMinCoil = m_COILGeometryPar->getCryoRmin();
185 double halfLength = m_COILGeometryPar->getCryoLength();
186 m_COILGeometryPar.addCallback([this, &offsetZ, &rMinCoil, &halfLength]() {
187 offsetZ = m_COILGeometryPar->getGlobalOffsetZ();
188 rMinCoil = m_COILGeometryPar->getCryoRmin();
189 halfLength = m_COILGeometryPar->getCryoLength();
190 });
191 m_TargetExt = new Simulation::ExtCylSurfaceTarget(rMinCoil, offsetZ - halfLength, offsetZ + halfLength);
192 G4ErrorPropagatorData::GetErrorPropagatorData()->SetTarget(m_TargetExt);
193 }
194 if (!m_BeamPipeGeo.isValid())
195 B2FATAL("Beam pipe geometry data are not available.");
196 double beampipeRadius = m_BeamPipeGeo->getParameter("Lv2OutBe.R2") * CLHEP::cm; // mm
197 m_BeamPipeGeo.addCallback([this, &beampipeRadius]() {
198 beampipeRadius = m_BeamPipeGeo->getParameter("Lv2OutBe.R2") * CLHEP::cm;
199 });
200 m_MinRadiusSq = beampipeRadius * beampipeRadius; // mm^2
201}

◆ registerVolumes()

void registerVolumes ( )
private

Register the list of geant4 physical volumes whose entry/exit points will be saved during extrapolation.

Definition at line 832 of file TrackExtrapolateG4e.cc.

833{
834 G4PhysicalVolumeStore* pvStore = G4PhysicalVolumeStore::GetInstance();
835 if (pvStore->size() == 0) {
836 B2FATAL("No geometry defined. Please create the geometry first.");
837 }
838 if (m_EnterExit != nullptr) // Only do this once
839 return;
840 m_EnterExit = new std::map<G4VPhysicalVolume*, enum VolTypes>;
841 m_BKLMVolumes = new std::vector<G4VPhysicalVolume*>;
842 for (std::vector<G4VPhysicalVolume*>::iterator iVol = pvStore->begin();
843 iVol != pvStore->end(); ++iVol) {
844 const G4String name = (*iVol)->GetName();
845 // Do not store ExtHits in CDC, so let's start from TOP and ARICH.
846 // TOP doesn't have one envelope; it has 16 "TOPModule"s
847 if (name.find("TOPModule") != std::string::npos) {
848 (*m_EnterExit)[*iVol] = VOLTYPE_TOP1;
849 }
850 // TOP quartz bar (=sensitive)
851 else if (name.find("_TOPPrism_") != std::string::npos ||
852 name.find("_TOPBarSegment") != std::string::npos ||
853 name.find("_TOPMirrorSegment") != std::string::npos) {
854 (*m_EnterExit)[*iVol] = VOLTYPE_TOP2;
855 }
856 // TOP quartz glue (not sensitive?)
857 else if (name.find("TOPBarSegment1Glue") != std::string::npos ||
858 name.find("TOPBarSegment2Glue") != std::string::npos ||
859 name.find("TOPMirrorSegmentGlue") != std::string::npos) {
860 (*m_EnterExit)[*iVol] = VOLTYPE_TOP3;
861 // ARICH volumes
862 } else if (name == "ARICH.AerogelSupportPlate") {
863 (*m_EnterExit)[*iVol] = VOLTYPE_ARICH1;
864 } else if (name == "ARICH.AerogelImgPlate") {
865 (*m_EnterExit)[*iVol] = VOLTYPE_ARICH2;
866 } else if (name.find("ARICH.HAPDWindow") != std::string::npos) {
867 (*m_EnterExit)[*iVol] = VOLTYPE_ARICH3;
868 }
869 // ECL crystal
870 else if (name.find("lv_barrel_crystal_") != std::string::npos ||
871 name.find("lv_forward_crystal_") != std::string::npos ||
872 name.find("lv_backward_crystal_") != std::string::npos) {
873 (*m_EnterExit)[*iVol] = VOLTYPE_ECL;
874 }
875 // Barrel KLM: BKLM.Layer**GasPhysical for RPCs or BKLM.Layer**ChimneyGasPhysical for RPCs
876 // BKLM.ScintActiveType*Physical for scintillator strips
877 else if (name.compare(0, 5, "BKLM.") == 0) {
878 if (name.find("GasPhysical") != std::string::npos) {
879 (*m_EnterExit)[*iVol] = VOLTYPE_BKLM1;
880 } else if (name.find("ScintActiveType") != std::string::npos) {
881 (*m_EnterExit)[*iVol] = VOLTYPE_BKLM2;
882 } else if ((name.find("ScintType") != std::string::npos) ||
883 (name.find("ElectrodePhysical") != std::string::npos)) {
884 m_BKLMVolumes->push_back(*iVol);
885 }
886 }
887 // Endcap KLM: StripSensitive_*
888 else if (name.compare(0, 14, "StripSensitive") == 0) {
889 (*m_EnterExit)[*iVol] = VOLTYPE_EKLM;
890 }
891 }
892
893}

◆ swim() [1/2]

void swim ( ExtState & extState,
G4ErrorFreeTrajState & g4eState )
private

Swim a single track (EXT) until it stops or leaves the target cylinder.

Definition at line 751 of file TrackExtrapolateG4e.cc.

752{
753 if (extState.pdgCode == 0)
754 return;
755 if (g4eState.GetMomentum().perp() <= m_MinPt)
756 return;
757 if (m_TargetExt->GetDistanceFromPoint(g4eState.GetPosition()) < 0.0)
758 return;
759 G4ParticleDefinition* particle = G4ParticleTable::GetParticleTable()->FindParticle(extState.pdgCode);
760 double mass = particle->GetPDGMass();
761 double minPSq = (mass + m_MinKE) * (mass + m_MinKE) - mass * mass;
762 G4ErrorMode propagationMode = (extState.isCosmic ? G4ErrorMode_PropBackwards : G4ErrorMode_PropForwards);
763 m_ExtMgr->InitTrackPropagation(propagationMode);
764 while (true) {
765 const G4int errCode = m_ExtMgr->PropagateOneStep(&g4eState, propagationMode);
766 G4Track* track = g4eState.GetG4Track();
767 const G4Step* step = track->GetStep();
768 const G4StepPoint* preStepPoint = step->GetPreStepPoint();
769 const G4StepPoint* postStepPoint = step->GetPostStepPoint();
770 G4TouchableHandle preTouch = preStepPoint->GetTouchableHandle();
771 G4VPhysicalVolume* pVol = preTouch->GetVolume();
772 const G4int preStatus = preStepPoint->GetStepStatus();
773 const G4int postStatus = postStepPoint->GetStepStatus();
774 G4ThreeVector pos = track->GetPosition(); // this is at postStepPoint
775 G4ThreeVector mom = track->GetMomentum(); // ditto
776 // First step on this track?
777 if (extState.isCosmic)
778 mom = -mom;
779 if (preStatus == fUndefined) {
780 if (m_EnterExit->find(pVol) != m_EnterExit->end()) {
781 createExtHit(EXT_FIRST, extState, g4eState, preStepPoint, preTouch);
782 }
783 }
784 // Ignore the zero-length step by PropagateOneStep() at each boundary
785 if (step->GetStepLength() > 0.0) {
786 double dt = step->GetDeltaTime();
787 double dl = step->GetStepLength() / track->GetMaterial()->GetRadlen();
788 if (preStatus == fGeomBoundary) { // first step in this volume?
789 if (m_EnterExit->find(pVol) != m_EnterExit->end()) {
790 createExtHit(EXT_ENTER, extState, g4eState, preStepPoint, preTouch);
791 }
792 }
793 if (extState.isCosmic) {
794 extState.tof -= dt;
795 extState.length -= dl;
796 } else {
797 extState.tof += dt;
798 extState.length += dl;
799 }
800 // Last step in this volume?
801 if (postStatus == fGeomBoundary) {
802 if (m_EnterExit->find(pVol) != m_EnterExit->end()) {
803 createExtHit(EXT_EXIT, extState, g4eState, postStepPoint, preTouch);
804 }
805 }
806 }
807 // Post-step momentum too low?
808 if (errCode || (mom.mag2() < minPSq)) {
809 if (m_EnterExit->find(pVol) != m_EnterExit->end()) {
810 createExtHit(EXT_STOP, extState, g4eState, postStepPoint, preTouch);
811 }
812 break;
813 }
814 // Detect escapes from the imaginary target cylinder.
815 if (m_TargetExt->GetDistanceFromPoint(pos) < 0.0) {
816 if (m_EnterExit->find(pVol) != m_EnterExit->end()) {
817 createExtHit(EXT_ESCAPE, extState, g4eState, postStepPoint, preTouch);
818 }
819 break;
820 }
821 // Stop extrapolating as soon as the track curls inward too much
822 if (pos.perp2() < m_MinRadiusSq) {
823 break;
824 }
825 } // track-extrapolation "infinite" loop
826
827 m_ExtMgr->EventTermination(propagationMode);
828
829}
void createExtHit(const ExtHitStatus, const ExtState &, const G4ErrorFreeTrajState &, const G4StepPoint *, const G4TouchableHandle &)
Create another EXT extrapolation hit for a track candidate.
double length
Length from start of extrapolation (rad lengths), updated during extrapolation.

◆ swim() [2/2]

void swim ( ExtState & extState,
G4ErrorFreeTrajState & g4eState,
const std::vector< std::pair< ECLCluster *, G4ThreeVector > > * eclClusterInfo,
const std::vector< std::pair< KLMCluster *, G4ThreeVector > > * klmClusterInfo,
std::vector< std::map< const Track *, double > > * bklmHitUsed )
private

Swim a single track (MUID) until it stops or leaves the target cylinder.

Definition at line 504 of file TrackExtrapolateG4e.cc.

508{
509 if (extState.pdgCode == 0)
510 return;
511 if (g4eState.GetMomentum().perp() <= m_MinPt)
512 return;
513 if (m_TargetMuid->GetDistanceFromPoint(g4eState.GetPosition()) < 0.0)
514 return;
515 G4ParticleDefinition* particle = G4ParticleTable::GetParticleTable()->FindParticle(extState.pdgCode);
516 double mass = particle->GetPDGMass();
517 double minPSq = (mass + m_MinKE) * (mass + m_MinKE) - mass * mass;
518 // Create structures for ECLCluster-Track matching
519 std::vector<double> eclClusterDistance;
520 std::vector<ExtHit> eclHit1, eclHit2, eclHit3;
521 if (eclClusterInfo != nullptr) {
522 eclClusterDistance.resize(eclClusterInfo->size(), 1.0E10); // "positive infinity"
523 ExtHit tempExtHit(.0, extState.pdgCode, Const::EDetector::ECL, 0, EXT_FIRST,
524 extState.isCosmic,
525 G4ThreeVector(), G4ThreeVector(), G4ErrorSymMatrix(6));
526 eclHit1.resize(eclClusterInfo->size(), tempExtHit);
527 eclHit2.resize(eclClusterInfo->size(), tempExtHit);
528 eclHit3.resize(eclClusterInfo->size(), tempExtHit);
529 }
530 // Create structures for KLMCluster-Track matching
531 std::vector<TrackClusterSeparation> klmHit;
532 if (klmClusterInfo != nullptr) {
533 klmHit.resize(klmClusterInfo->size()); // initialize each to huge distance
534 }
535 KLMMuidLikelihood* klmMuidLikelihood = m_klmMuidLikelihoods.appendNew(); // rest of this object will be filled later
536 klmMuidLikelihood->setPDGCode(extState.pdgCode);
537 if (extState.track != nullptr)
538 extState.track->addRelationTo(klmMuidLikelihood);
539 G4ErrorMode propagationMode = (extState.isCosmic ? G4ErrorMode_PropBackwards : G4ErrorMode_PropForwards);
540 m_ExtMgr->InitTrackPropagation(propagationMode);
541 while (true) {
542 const G4int errCode = m_ExtMgr->PropagateOneStep(&g4eState, propagationMode);
543 G4Track* track = g4eState.GetG4Track();
544 const G4Step* step = track->GetStep();
545 const G4StepPoint* preStepPoint = step->GetPreStepPoint();
546 const G4StepPoint* postStepPoint = step->GetPostStepPoint();
547 G4TouchableHandle preTouch = preStepPoint->GetTouchableHandle();
548 G4VPhysicalVolume* pVol = preTouch->GetVolume();
549 const G4int preStatus = preStepPoint->GetStepStatus();
550 const G4int postStatus = postStepPoint->GetStepStatus();
551 G4ThreeVector pos = track->GetPosition(); // this is at postStepPoint
552 G4ThreeVector mom = track->GetMomentum(); // ditto
553 // Ignore the zero-length step by PropagateOneStep() at each boundary
554 if (extState.isCosmic) mom = -mom;
555 if (step->GetStepLength() > 0.0) {
556 double dt = step->GetDeltaTime();
557 double dl = step->GetStepLength() / track->GetMaterial()->GetRadlen();
558 if (preStatus == fGeomBoundary) { // first step in this volume?
559 if (m_TargetExt->GetDistanceFromPoint(pos) < 0.0) { // only hits outside ECL during muid
560 if (m_EnterExit->find(pVol) != m_EnterExit->end()) {
561 createExtHit(EXT_ENTER, extState, g4eState, preStepPoint, preTouch);
562 }
563 }
564 }
565 if (extState.isCosmic) {
566 extState.tof -= dt;
567 extState.length -= dl;
568 } else {
569 extState.tof += dt;
570 extState.length += dl;
571 }
572 if (postStatus == fGeomBoundary) { // last step in this volume?
573 if (m_TargetExt->GetDistanceFromPoint(pos) < 0.0) { // only hits outside ECL during muid
574 if (m_EnterExit->find(pVol) != m_EnterExit->end()) {
575 createExtHit(EXT_EXIT, extState, g4eState, postStepPoint, preTouch);
576 }
577 }
578 }
579 if (createMuidHit(extState, g4eState, klmMuidLikelihood, bklmHitUsed)) {
580 // Force geant4e to update its G4Track from the Kalman-updated state
581 m_ExtMgr->GetPropagator()->SetStepN(0);
582 }
583 if (eclClusterInfo != nullptr) {
584 for (unsigned int c = 0; c < eclClusterInfo->size(); ++c) {
585 G4ThreeVector eclPos((*eclClusterInfo)[c].second);
586 G4ThreeVector prePos(preStepPoint->GetPosition());
587 G4ThreeVector diff(prePos - eclPos);
588 double distance = diff.mag();
589 if (distance < m_MaxECLTrackClusterDistance) {
590 // fallback ECLNEAR in case no ECLCROSS is found
591 if (distance < eclClusterDistance[c]) {
592 eclClusterDistance[c] = distance;
593 G4ErrorSymMatrix covariance(6, 0);
594 fromG4eToPhasespace(g4eState, covariance);
595 eclHit3[c].update(EXT_ECLNEAR, extState.tof, pos / CLHEP::cm, mom / CLHEP::GeV, covariance);
596 }
597 // find position of crossing of the track with the ECLCluster's sphere
598 if (eclHit1[c].getStatus() == EXT_FIRST) {
599 if (pos.mag2() >= eclPos.mag2()) {
600 double r = eclPos.mag();
601 double preD = prePos.mag() - r;
602 double postD = pos.mag() - r;
603 double f = postD / (postD - preD);
604 G4ThreeVector midPos = pos + (prePos - pos) * f;
605 double tof = extState.tof + dt * f * (extState.isCosmic ? +1 : -1); // in ns, at end of step
606 G4ErrorSymMatrix covariance(6, 0);
607 fromG4eToPhasespace(g4eState, covariance);
608 eclHit1[c].update(EXT_ECLCROSS, tof, midPos / CLHEP::cm, mom / CLHEP::GeV, covariance);
609 }
610 }
611 }
612 // find closest distance to the radial line to the ECLCluster
613 if (eclHit2[c].getStatus() == EXT_FIRST) {
614 G4ThreeVector delta(pos - prePos);
615 G4ThreeVector perp(eclPos.cross(delta));
616 double perpMag2 = perp.mag2();
617 if (perpMag2 > 1.0E-10) {
618 double dist = std::fabs(diff * perp) / std::sqrt(perpMag2);
619 if (dist < m_MaxECLTrackClusterDistance) {
620 double f = eclPos * (prePos.cross(perp)) / perpMag2;
621 if ((f > -0.5) && (f <= 1.0)) {
622 G4ThreeVector midPos(prePos + f * delta);
623 double length = extState.length + dl * (1.0 - f) * (extState.isCosmic ? +1 : -1);
624 G4ErrorSymMatrix covariance(6, 0);
625 fromG4eToPhasespace(g4eState, covariance);
626 eclHit2[c].update(EXT_ECLDL, length, midPos / CLHEP::cm, mom / CLHEP::GeV, covariance);
627 }
628 }
629 }
630 }
631 }
632 }
633 if (klmClusterInfo != nullptr) {
634 for (unsigned int c = 0; c < klmClusterInfo->size(); ++c) {
635 G4ThreeVector klmPos = (*klmClusterInfo)[c].second;
636 G4ThreeVector separation = klmPos - pos;
637 double distance = separation.mag();
638 if (distance < klmHit[c].getDistance()) {
639 klmHit[c].setDistance(distance);
640 klmHit[c].setTrackClusterAngle(mom.angle(separation));
641 klmHit[c].setTrackClusterSeparationAngle(mom.angle(klmPos));
642 klmHit[c].setTrackRotationAngle(extState.directionAtIP.angle(mom));
643 klmHit[c].setTrackClusterInitialSeparationAngle(extState.directionAtIP.angle(klmPos));
644 }
645 }
646 }
647 }
648 // Post-step momentum too low?
649 if (errCode || (mom.mag2() < minPSq)) {
650 break;
651 }
652 // Detect escapes from the imaginary target cylinder.
653 if (m_TargetMuid->GetDistanceFromPoint(pos) < 0.0) {
654 break;
655 }
656 // Stop extrapolating as soon as the track curls inward too much
657 if (pos.perp2() < m_MinRadiusSq) {
658 break;
659 }
660 } // track-extrapolation "infinite" loop
661
662 m_ExtMgr->EventTermination(propagationMode);
663
664 finishTrack(extState, klmMuidLikelihood, (g4eState.GetPosition().z() > m_OffsetZ));
665
666 if (eclClusterInfo != nullptr) {
667 for (unsigned int c = 0; c < eclClusterInfo->size(); ++c) {
668 if (eclHit1[c].getStatus() != EXT_FIRST) {
669 ExtHit* h = m_extHits.appendNew(eclHit1[c]);
670 (*eclClusterInfo)[c].first->addRelationTo(h);
671 if (extState.track != nullptr) {
672 extState.track->addRelationTo(h);
673 }
674 }
675 if (eclHit2[c].getStatus() != EXT_FIRST) {
676 ExtHit* h = m_extHits.appendNew(eclHit2[c]);
677 (*eclClusterInfo)[c].first->addRelationTo(h);
678 if (extState.track != nullptr) {
679 extState.track->addRelationTo(h);
680 }
681 }
682 if (eclHit3[c].getStatus() != EXT_FIRST) {
683 ExtHit* h = m_extHits.appendNew(eclHit3[c]);
684 (*eclClusterInfo)[c].first->addRelationTo(h);
685 if (extState.track != nullptr) {
686 extState.track->addRelationTo(h);
687 }
688 }
689 }
690 }
691
692 if (klmClusterInfo != nullptr) {
693 // here we set a relation only to the closest KLMCluster
694 // and we don't set any relation if the distance is too large
695 double minDistance = m_MaxKLMTrackClusterDistance;
696 unsigned int closestCluster = 0;
697 for (unsigned int c = 0; c < klmClusterInfo->size(); ++c) {
698 if (klmHit[c].getDistance() > 1.0E9) {
699 continue;
700 }
701 TrackClusterSeparation* h = m_trackClusterSeparations.appendNew(klmHit[c]);
702 if (h != nullptr) {
703 (*klmClusterInfo)[c].first->addRelationTo(h); // relation KLMCluster to TrackClusterSeparation
704 if (extState.track != nullptr) {
705 extState.track->addRelationTo(h); // relation Track to TrackClusterSeparation
706 }
707 }
708 if (klmHit[c].getDistance() < minDistance) {
709 closestCluster = c;
710 minDistance = klmHit[c].getDistance();
711 }
712 }
713 if (minDistance < m_MaxKLMTrackClusterDistance) {
714 // set the relation Track to KLMCluster, using the distance as weight
715 // but for being consistent with the basf2 conventions, store the distance in cm
716 if (extState.track != nullptr) {
717 extState.track->addRelationTo((*klmClusterInfo)[closestCluster].first, 1. / (minDistance / CLHEP::cm));
718 }
719 }
720 }
721}
void setPDGCode(int pdg)
Set the PDG code of the particle hypothesis used during the extrapolation.
bool createMuidHit(ExtState &, G4ErrorFreeTrajState &, KLMMuidLikelihood *, std::vector< std::map< const Track *, double > > *)
Create another MUID extrapolation hit for a track candidate.
void finishTrack(const ExtState &, KLMMuidLikelihood *, bool)
Complete muon identification after end of track extrapolation.

◆ terminate()

void terminate ( bool flag)

Terminates this singleton.

Parameters
flagTrue if called by Muid module, false if called by Ext module.

Definition at line 433 of file TrackExtrapolateG4e.cc.

434{
435 if (m_DefaultHypotheses != nullptr)
436 delete m_DefaultHypotheses;
437 if (byMuid) {
438 delete m_TargetMuid;
439 for (auto const& muidBuilder : m_MuidBuilderMap)
440 delete muidBuilder.second;
441 }
442 if (m_TargetExt != nullptr) {
443 delete m_TargetExt;
444 m_TargetExt = nullptr;
445 }
446 if (m_EnterExit != nullptr) {
447 delete m_EnterExit;
448 delete m_BKLMVolumes;
449 m_ExtMgr->RunTermination();
450 m_EnterExit = nullptr;
451 m_BKLMVolumes = nullptr;
452 }
453}

Member Data Documentation

◆ m_addHitsToRecoTrack

bool m_addHitsToRecoTrack = false
private

Parameter to add the found hits also to the reco tracks or not. Is turned off by default.

Definition at line 429 of file TrackExtrapolateG4e.h.

◆ m_BarrelHalfLength

double m_BarrelHalfLength
private

half-length (cm) of the barrel

Definition at line 380 of file TrackExtrapolateG4e.h.

◆ m_BarrelMaxR

double m_BarrelMaxR
private

maximum radius (cm) of the barrel

Definition at line 374 of file TrackExtrapolateG4e.h.

◆ m_BarrelMinR

double m_BarrelMinR
private

minimum radius (cm) of the barrel

Definition at line 377 of file TrackExtrapolateG4e.h.

◆ m_BarrelModuleMiddleRadius

double m_BarrelModuleMiddleRadius[2][BKLMElementNumbers::getMaximalSectorNumber()+1][BKLMElementNumbers::getMaximalLayerNumber()+1]
private

hit-plane radius (cm) at closest distance to IP of each barrel end | sector | layer

Definition at line 395 of file TrackExtrapolateG4e.h.

◆ m_BarrelNSector

int m_BarrelNSector
private

Number of barrel sectors.

Definition at line 371 of file TrackExtrapolateG4e.h.

◆ m_BarrelPhiStripVariance

double m_BarrelPhiStripVariance[BKLMElementNumbers::getMaximalLayerNumber()+1]
private

BKLM RPC phi-measuring strip position variance (cm^2) by layer.

Definition at line 386 of file TrackExtrapolateG4e.h.

◆ m_BarrelScintVariance

double m_BarrelScintVariance
private

BKLM scintillator strip position variance (cm^2)

Definition at line 392 of file TrackExtrapolateG4e.h.

◆ m_BarrelSectorPerp

G4ThreeVector m_BarrelSectorPerp[BKLMElementNumbers::getMaximalSectorNumber()+1]
private

normal unit vector of each barrel sector

Definition at line 399 of file TrackExtrapolateG4e.h.

◆ m_BarrelSectorPhi

G4ThreeVector m_BarrelSectorPhi[BKLMElementNumbers::getMaximalSectorNumber()+1]
private

azimuthal unit vector of each barrel sector

Definition at line 402 of file TrackExtrapolateG4e.h.

◆ m_BarrelZStripVariance

double m_BarrelZStripVariance[BKLMElementNumbers::getMaximalLayerNumber()+1]
private

BKLM RPC z-measuring strip position variance (cm^2) by layer.

Definition at line 389 of file TrackExtrapolateG4e.h.

◆ m_BeamPipeGeo

DBObjPtr<BeamPipeGeo> m_BeamPipeGeo
private

Conditions-database object for beam pipe geometry.

Definition at line 362 of file TrackExtrapolateG4e.h.

◆ m_BKLMVolumes

std::vector<G4VPhysicalVolume*>* m_BKLMVolumes
private

Pointers to BKLM geant4 sensitive (physical) volumes.

Definition at line 350 of file TrackExtrapolateG4e.h.

◆ m_COILGeometryPar

DBObjPtr<COILGeometryPar> m_COILGeometryPar
private

Conditions-database object for COIL geometry.

Definition at line 359 of file TrackExtrapolateG4e.h.

◆ m_DefaultHypotheses

std::vector<Const::ChargedStable>* m_DefaultHypotheses
private

Default ChargedStable hypotheses (needed as call argument but not used)

Definition at line 344 of file TrackExtrapolateG4e.h.

◆ m_eclClusters

StoreArray<ECLCluster> m_eclClusters
private

ECL clusters.

Definition at line 453 of file TrackExtrapolateG4e.h.

◆ m_eklmElementNumbers

const EKLMElementNumbers* m_eklmElementNumbers
private

EKLM element numbers.

Definition at line 435 of file TrackExtrapolateG4e.h.

◆ m_eklmTransformData

const EKLM::TransformDataGlobalAligned* m_eklmTransformData
private

EKLM transformation data.

Definition at line 441 of file TrackExtrapolateG4e.h.

◆ m_EndcapHalfLength

double m_EndcapHalfLength
private

half-length (cm) of either endcap

Definition at line 414 of file TrackExtrapolateG4e.h.

◆ m_EndcapMaxR

double m_EndcapMaxR
private

maximum radius (cm) of the endcaps

Definition at line 405 of file TrackExtrapolateG4e.h.

◆ m_EndcapMiddleZ

double m_EndcapMiddleZ
private

midpoint along z (cm) of the forward endcap from the KLM midpoint

Definition at line 411 of file TrackExtrapolateG4e.h.

◆ m_EndcapMinR

double m_EndcapMinR
private

minimum radius (cm) of the endcaps

Definition at line 408 of file TrackExtrapolateG4e.h.

◆ m_EndcapModuleMiddleZ

double m_EndcapModuleMiddleZ[BKLMElementNumbers::getMaximalLayerNumber()+1]
private

hit-plane z (cm) of each IP layer relative to KLM midpoint

Definition at line 426 of file TrackExtrapolateG4e.h.

◆ m_EndcapScintVariance

double m_EndcapScintVariance
private

EKLM scintillator strip position variance (cm^2)

Definition at line 423 of file TrackExtrapolateG4e.h.

◆ m_EnterExit

std::map<G4VPhysicalVolume*, enum VolTypes>* m_EnterExit
private

Pointers to geant4 physical volumes whose entry/exit points will be saved.

Definition at line 347 of file TrackExtrapolateG4e.h.

◆ m_extHits

StoreArray<ExtHit> m_extHits
private

Ext hits.

Definition at line 456 of file TrackExtrapolateG4e.h.

◆ m_ExtInitialized

bool m_ExtInitialized
private

Flag to indicate that EXT initialize() has been called.

Definition at line 305 of file TrackExtrapolateG4e.h.

◆ m_ExtMgr

Simulation::ExtManager* m_ExtMgr
private

Pointer to the ExtManager singleton.

Definition at line 335 of file TrackExtrapolateG4e.h.

◆ m_HypothesesExt

const std::vector<Const::ChargedStable>* m_HypothesesExt
private

ChargedStable hypotheses for EXT.

Definition at line 338 of file TrackExtrapolateG4e.h.

◆ m_HypothesesMuid

const std::vector<Const::ChargedStable>* m_HypothesesMuid
private

ChargedStable hypotheses for MUID.

Definition at line 341 of file TrackExtrapolateG4e.h.

◆ m_klmChannelStatus

DBObjPtr<KLMChannelStatus> m_klmChannelStatus
private

Conditions-database object for KLM channel status.

Definition at line 444 of file TrackExtrapolateG4e.h.

◆ m_klmClusters

StoreArray<KLMCluster> m_klmClusters
private

KLM clusters.

Definition at line 462 of file TrackExtrapolateG4e.h.

◆ m_klmElementNumbers

const KLMElementNumbers* m_klmElementNumbers
private

KLM element numbers.

Definition at line 438 of file TrackExtrapolateG4e.h.

◆ m_klmHit2ds

StoreArray<KLMHit2d> m_klmHit2ds
private

KLM 2d hits.

Definition at line 459 of file TrackExtrapolateG4e.h.

◆ m_klmLikelihoodParameters

DBObjPtr<KLMLikelihoodParameters> m_klmLikelihoodParameters
private

Conditions-database object for KLM likelihood parameters.

Definition at line 450 of file TrackExtrapolateG4e.h.

◆ m_klmMuidHits

StoreArray<KLMMuidHit> m_klmMuidHits
private

KLM muid hits.

Definition at line 465 of file TrackExtrapolateG4e.h.

◆ m_klmMuidLikelihoods

StoreArray<KLMMuidLikelihood> m_klmMuidLikelihoods
private

KLM muid likelihoods.

Definition at line 468 of file TrackExtrapolateG4e.h.

◆ m_klmStripEfficiency

DBObjPtr<KLMStripEfficiency> m_klmStripEfficiency
private

Conditions-database object for KLM strip efficiency.

Definition at line 447 of file TrackExtrapolateG4e.h.

◆ m_MagneticField

double m_MagneticField
private

Magnetic field z component (gauss) at origin.

Definition at line 317 of file TrackExtrapolateG4e.h.

◆ m_MaxDistSqInVariances

double m_MaxDistSqInVariances
private

user-defined maximum squared-distance (in number of variances) for matching hit to extrapolation

Definition at line 320 of file TrackExtrapolateG4e.h.

◆ m_MaxDt

double m_MaxDt
private

Coincidence window half-width for in-time KLM hits (ns)

Definition at line 314 of file TrackExtrapolateG4e.h.

◆ m_MaxECLTrackClusterDistance

double m_MaxECLTrackClusterDistance
private

user-defined maximum distance (mm) between ECLCluster and associated track (for EID)

Definition at line 326 of file TrackExtrapolateG4e.h.

◆ m_MaxKLMTrackClusterDistance

double m_MaxKLMTrackClusterDistance
private

user-defined maximum distance (mm) between KLMCluster and associated track (for KLID)

Definition at line 323 of file TrackExtrapolateG4e.h.

◆ m_MeanDt

double m_MeanDt
private

Mean hit - trigger time (ns)

Definition at line 311 of file TrackExtrapolateG4e.h.

◆ m_MinKE

double m_MinKE
private

Minimum kinetic energy in MeV for extrapolation to continue.

Definition at line 332 of file TrackExtrapolateG4e.h.

◆ m_MinPt

double m_MinPt
private

Minimum transverse momentum in MeV/c for extrapolation to be started.

Definition at line 329 of file TrackExtrapolateG4e.h.

◆ m_MinRadiusSq

double m_MinRadiusSq
private

Minimum squared radius (cm) outside of which extrapolation will continue.

Definition at line 365 of file TrackExtrapolateG4e.h.

◆ m_MuidBuilderMap

std::map<int, MuidBuilder*> m_MuidBuilderMap
private

PDF for the charged final state particle hypotheses.

Definition at line 432 of file TrackExtrapolateG4e.h.

◆ m_MuidInitialized

bool m_MuidInitialized
private

Flag to indicate that MUID initialize() has been called.

Definition at line 308 of file TrackExtrapolateG4e.h.

◆ m_OffsetZ

double m_OffsetZ
private

offset (cm) along z axis of KLM midpoint from IP

Definition at line 368 of file TrackExtrapolateG4e.h.

◆ m_OutermostActiveBackwardEndcapLayer

int m_OutermostActiveBackwardEndcapLayer
private

outermost backward-endcap layer that is active for muon identification (user-defined)

Definition at line 420 of file TrackExtrapolateG4e.h.

◆ m_OutermostActiveBarrelLayer

int m_OutermostActiveBarrelLayer
private

outermost barrel layer that is active for muon identification (user-defined)

Definition at line 383 of file TrackExtrapolateG4e.h.

◆ m_OutermostActiveForwardEndcapLayer

int m_OutermostActiveForwardEndcapLayer
private

outermost forward-endcap layer that is active for muon identification (user-defined)

Definition at line 417 of file TrackExtrapolateG4e.h.

◆ m_recoTracks

StoreArray<RecoTrack> m_recoTracks
private

Reco tracks.

Definition at line 471 of file TrackExtrapolateG4e.h.

◆ m_Singleton

TrackExtrapolateG4e * m_Singleton = nullptr
staticprivate

Stores pointer to the singleton class.

Definition at line 302 of file TrackExtrapolateG4e.h.

◆ m_TargetExt

Simulation::ExtCylSurfaceTarget* m_TargetExt
private

virtual "target" cylinder for EXT (boundary beyond which extrapolation ends)

Definition at line 353 of file TrackExtrapolateG4e.h.

◆ m_TargetMuid

Simulation::ExtCylSurfaceTarget* m_TargetMuid
private

virtual "target" cylinder for MUID (boundary beyond which extrapolation ends)

Definition at line 356 of file TrackExtrapolateG4e.h.

◆ m_trackClusterSeparations

StoreArray<TrackClusterSeparation> m_trackClusterSeparations
private

Track cluster separation.

Definition at line 477 of file TrackExtrapolateG4e.h.

◆ m_tracks

StoreArray<Track> m_tracks
private

Tracks.

Definition at line 474 of file TrackExtrapolateG4e.h.


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