Belle II Software  release-05-01-25
TrackExtrapolateG4e.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010-2011 Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Leo Piilonen *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 /* Own header. */
12 #include <tracking/trackExtrapolateG4e/TrackExtrapolateG4e.h>
13 
14 /* Belle 2 headers. */
15 #include <ecl/geometry/ECLGeometryPar.h>
16 #include <framework/datastore/StoreArray.h>
17 #include <framework/datastore/StoreObjPtr.h>
18 #include <framework/gearbox/GearDir.h>
19 #include <framework/geometry/BFieldManager.h>
20 #include <framework/logging/Logger.h>
21 #include <genfit/Exception.h>
22 #include <klm/dataobjects/bklm/BKLMStatus.h>
23 #include <klm/bklm/geometry/GeometryPar.h>
24 #include <klm/bklm/geometry/Module.h>
25 #include <klm/dataobjects/KLMMuidLikelihood.h>
26 #include <klm/dataobjects/KLMMuidHit.h>
27 #include <klm/dataobjects/eklm/EKLMElementNumbers.h>
28 #include <klm/muid/MuidBuilder.h>
29 #include <klm/muid/MuidElementNumbers.h>
30 #include <mdst/dataobjects/ECLCluster.h>
31 #include <mdst/dataobjects/KLMCluster.h>
32 #include <mdst/dataobjects/Track.h>
33 #include <simulation/kernel/ExtCylSurfaceTarget.h>
34 #include <simulation/kernel/ExtManager.h>
35 
36 /* CLHEP headers. */
37 #include <CLHEP/Matrix/Vector.h>
38 #include <CLHEP/Units/PhysicalConstants.h>
39 #include <CLHEP/Units/SystemOfUnits.h>
40 
41 /* Geant4 headers. */
42 #include <G4ErrorFreeTrajState.hh>
43 #include <G4ErrorMatrix.hh>
44 #include <G4ErrorPropagatorData.hh>
45 #include <G4ErrorSymMatrix.hh>
46 #include <G4ParticleTable.hh>
47 #include <G4PhysicalVolumeStore.hh>
48 #include <G4Point3D.hh>
49 #include <G4StateManager.hh>
50 #include <G4Step.hh>
51 #include <G4StepPoint.hh>
52 #include <G4ThreeVector.hh>
53 #include <G4TouchableHandle.hh>
54 #include <G4Track.hh>
55 #include <G4UImanager.hh>
56 #include <G4VPhysicalVolume.hh>
57 #include <G4VTouchable.hh>
58 
59 /* C++ headers. */
60 #include <algorithm>
61 #include <cmath>
62 #include <iostream>
63 
64 #define TWOPI (2.0*M_PI)
65 #define PI_8 (0.125*M_PI)
66 #define DEPTH_RPC 9
67 #define DEPTH_SCINT 11
68 
69 using namespace std;
70 using namespace Belle2;
71 
72 TrackExtrapolateG4e* TrackExtrapolateG4e::m_Singleton = nullptr;
73 
74 TrackExtrapolateG4e* TrackExtrapolateG4e::getInstance()
75 {
76  if (m_Singleton == nullptr) m_Singleton = new TrackExtrapolateG4e;
77  return m_Singleton;
78 }
79 
80 TrackExtrapolateG4e::TrackExtrapolateG4e() :
81  m_ExtInitialized(false), // initialized later
82  m_MuidInitialized(false), // initialized later
83  m_MeanDt(0.0), // initialized later
84  m_MaxDt(0.0), // initialized later
85  m_MagneticField(0.0), // initialized later
86  m_MaxDistSqInVariances(0.0), // initialized later
87  m_MaxKLMTrackClusterDistance(0.0), // initialized later
88  m_MaxECLTrackClusterDistance(0.0), // initialized later
89  m_MinPt(0.0), // initialized later
90  m_MinKE(0.0), // initialized later
91  m_ExtMgr(nullptr), // initialized later
92  m_HypothesesExt(nullptr), // initialized later
93  m_HypothesesMuid(nullptr), // initialized later
94  m_DefaultHypotheses(nullptr), // initialized later
95  m_EnterExit(nullptr), // initialized later
96  m_BKLMVolumes(nullptr), // initialized later
97  m_TargetExt(nullptr), // initialized later
98  m_TargetMuid(nullptr), // initialized later
99  m_MinRadiusSq(0.0), // initialized later
100  m_OffsetZ(0.0), // initialized later
101  m_BarrelNSector(0), // initialized later
102  m_BarrelMaxR(0.0), // initialized later
103  m_BarrelMinR(0.0), // initialized later
104  m_BarrelHalfLength(0.0), // initialized later
105  m_OutermostActiveBarrelLayer(0), // initialized later
106  m_BarrelScintVariance(0.0), // initialized later
107  m_EndcapMaxR(0.0), // initialized later
108  m_EndcapMinR(0.0), // initialized later
109  m_EndcapMiddleZ(0.0), // initialized later
110  m_EndcapHalfLength(0.0), // initialized later
111  m_OutermostActiveForwardEndcapLayer(0), // initialized later
112  m_OutermostActiveBackwardEndcapLayer(0), // initialized later
113  m_EndcapScintVariance(0.0), // initialized later
114  m_eklmTransformData(nullptr) // initialized later
115 {
116  for (int j = 0; j < BKLMElementNumbers::getMaximalLayerNumber() + 1; ++j) {
117  m_BarrelPhiStripVariance[j] = 0.0;
118  m_BarrelZStripVariance[j] = 0.0;
119  m_BarrelPhiStripVariance[j] = 0.0;
120  m_EndcapModuleMiddleZ[j] = 0.0;
121  }
122  for (int s = 0; s < BKLMElementNumbers::getMaximalSectorNumber() + 1; ++s) {
123  for (int j = 0; j < BKLMElementNumbers::getMaximalLayerNumber() + 1; ++j) {
124  m_BarrelModuleMiddleRadius[0][s][j] = 0.0;
125  m_BarrelModuleMiddleRadius[1][s][j] = 0.0;
126  }
127  m_BarrelSectorPerp[s] = G4ThreeVector(0.0, 0.0, 0.0);
128  m_BarrelSectorPhi[s] = G4ThreeVector(0.0, 0.0, 0.0);
129  }
130  m_klmElementNumbers = &(KLMElementNumbers::Instance());
131 }
132 
133 TrackExtrapolateG4e::~TrackExtrapolateG4e()
134 {
135 }
136 
137 // Initialize for EXT
138 void TrackExtrapolateG4e::initialize(double minPt, double minKE,
139  std::vector<Const::ChargedStable>& hypotheses)
140 {
141  m_ExtInitialized = true;
142 
143  // Define required objects, register the new ones and relations
144  m_recoTracks.isRequired();
145  m_tracks.isRequired();
146  m_extHits.registerInDataStore();
147  m_tracks.registerRelationTo(m_extHits);
148 
149  // Save the magnetic field z component (gauss) at the origin
150  m_MagneticField = BFieldManager::getField(B2Vector3D(0, 0, 0)).Z() / Unit::T * CLHEP::tesla / CLHEP::gauss;
151 
152  // Convert user cutoff values to geant4 units
153  m_MinPt = max(0.0, minPt) * CLHEP::GeV;
154  m_MinKE = max(0.0, minKE) * CLHEP::GeV;
155 
156  // Save pointer to the list of particle hypotheses for EXT extrapolation
157  m_HypothesesExt = &hypotheses;
158 
159  // Define the list of volumes that will have their entry and/or
160  // exit points stored during the extrapolation.
161  registerVolumes();
162 
163  // Store the address of the ExtManager (used later)
164  m_ExtMgr = Simulation::ExtManager::GetManager();
165 
166  // Set up the EXT-specific geometry (might have already been done by MUID)
167  if (m_TargetExt == nullptr) {
168  GearDir coilContent = GearDir("/Detector/DetectorComponent[@name=\"COIL\"]/Content/");
169  double offsetZ = coilContent.getLength("OffsetZ") * CLHEP::cm;
170  double rMaxCoil = coilContent.getLength("Cryostat/Rmin") * CLHEP::cm;
171  double halfLength = coilContent.getLength("Cryostat/HalfLength") * CLHEP::cm;
172  m_TargetExt = new Simulation::ExtCylSurfaceTarget(rMaxCoil, offsetZ - halfLength, offsetZ + halfLength);
173  G4ErrorPropagatorData::GetErrorPropagatorData()->SetTarget(m_TargetExt);
174  GearDir beampipeContent = GearDir("/Detector/DetectorComponent[@name=\"BeamPipe\"]/Content/");
175  double beampipeRadius = beampipeContent.getLength("Lv2OutBe/R2", 1.20) * CLHEP::cm; // mm
176  m_MinRadiusSq = beampipeRadius * beampipeRadius; // mm^2
177  }
178 }
179 
180 // Initialize for MUID
181 void TrackExtrapolateG4e::initialize(double meanDt, double maxDt, double maxKLMTrackHitDistance,
182  double maxKLMTrackClusterDistance, double maxECLTrackClusterDistance,
183  double minPt, double minKE, bool addHitsToRecoTrack,
184  std::vector<Const::ChargedStable>& hypotheses)
185 {
186  m_MuidInitialized = true;
187  m_addHitsToRecoTrack = addHitsToRecoTrack;
188 
189  // Define required objects, register the new ones and relations
190  m_eclClusters.isRequired();
191  m_bklmHit2ds.isRequired();
192  m_eklmHit2ds.isRequired();
193  m_klmClusters.isRequired();
194  m_recoTracks.isRequired();
195  m_tracks.isRequired();
196  m_extHits.registerInDataStore();
197  m_klmMuidLikelihoods.registerInDataStore();
198  m_klmMuidHits.registerInDataStore();
199  m_trackClusterSeparations.registerInDataStore();
200  m_tracks.registerRelationTo(m_extHits);
201  m_tracks.registerRelationTo(m_klmMuidLikelihoods);
202  m_tracks.registerRelationTo(m_klmMuidHits);
203  m_tracks.registerRelationTo(m_bklmHit2ds);
204  m_tracks.registerRelationTo(m_eklmHit2ds);
205  m_tracks.registerRelationTo(m_trackClusterSeparations);
206  m_tracks.registerRelationTo(m_klmClusters);
207  m_klmClusters.registerRelationTo(m_trackClusterSeparations);
208  m_eclClusters.registerRelationTo(m_extHits);
209  RecoTrack::registerRequiredRelations(m_recoTracks);
210 
211  // Save the in-time cut's central value and width for valid hits
212  m_MeanDt = meanDt;
213  m_MaxDt = maxDt;
214 
215  // Save the magnetic field z component (gauss) at the origin
216  m_MagneticField = BFieldManager::getField(B2Vector3D(0, 0, 0)).Z() / Unit::T * CLHEP::tesla / CLHEP::gauss;
217 
218  // Convert from sigma to variance for hit-position uncertainty
219  m_MaxDistSqInVariances = maxKLMTrackHitDistance * maxKLMTrackHitDistance;
220 
221  // Convert user cutoff values to geant4 units
222  m_MaxKLMTrackClusterDistance = max(0.0, maxKLMTrackClusterDistance) * CLHEP::cm;
223  m_MaxECLTrackClusterDistance = max(0.0, maxECLTrackClusterDistance) * CLHEP::cm;
224  m_MinPt = max(0.0, minPt) * CLHEP::GeV;
225  m_MinKE = max(0.0, minKE) * CLHEP::GeV;
226 
227  // Save pointer to the list of particle hypotheses for EXT extrapolation
228  m_HypothesesMuid = &hypotheses;
229 
230  // Define the list of volumes that will have their entry and/or
231  // exit points stored during the extrapolation.
232  registerVolumes();
233 
234  // Store the address of the ExtManager (used later)
235  m_ExtMgr = Simulation::ExtManager::GetManager();
236 
237  // Set up the EXT-specific geometry (might have already been done by EXT)
238  if (m_TargetExt == nullptr) {
239  GearDir coilContent = GearDir("/Detector/DetectorComponent[@name=\"COIL\"]/Content/");
240  double offsetZ = coilContent.getLength("OffsetZ") * CLHEP::cm;
241  double rMaxCoil = coilContent.getLength("Cryostat/Rmin") * CLHEP::cm;
242  double halfLength = coilContent.getLength("Cryostat/HalfLength") * CLHEP::cm;
243  m_TargetExt = new Simulation::ExtCylSurfaceTarget(rMaxCoil, offsetZ - halfLength, offsetZ + halfLength);
244  GearDir beampipeContent = GearDir("/Detector/DetectorComponent[@name=\"BeamPipe\"]/Content/");
245  double beampipeRadius = beampipeContent.getLength("Lv2OutBe/R2", 1.20) * CLHEP::cm; // mm
246  m_MinRadiusSq = beampipeRadius * beampipeRadius; // mm^2
247  }
248 
249  // Set up the MUID-specific geometry
250  bklm::GeometryPar* bklmGeometry = bklm::GeometryPar::instance();
251  const EKLM::GeometryData& eklmGeometry = EKLM::GeometryData::Instance();
252  m_BarrelHalfLength = bklmGeometry->getHalfLength() * CLHEP::cm; // in G4 units (mm)
253  m_EndcapHalfLength = 0.5 * eklmGeometry.getSectionPosition()->getLength(); // in G4 units (mm)
254  m_OffsetZ = bklmGeometry->getOffsetZ() * CLHEP::cm; // in G4 units (mm)
255  double minZ = m_OffsetZ - (m_BarrelHalfLength + 2.0 * m_EndcapHalfLength);
256  double maxZ = m_OffsetZ + (m_BarrelHalfLength + 2.0 * m_EndcapHalfLength);
257  m_BarrelNSector = bklmGeometry->getNSector();
258  m_BarrelMaxR = bklmGeometry->getOuterRadius() * CLHEP::cm / std::cos(M_PI / m_BarrelNSector); // in G4 units (mm)
259  m_TargetMuid = new Simulation::ExtCylSurfaceTarget(m_BarrelMaxR, minZ, maxZ);
260  G4ErrorPropagatorData::GetErrorPropagatorData()->SetTarget(m_TargetMuid);
261 
262  m_BarrelHalfLength /= CLHEP::cm; // now in G4e units (cm)
263  m_EndcapHalfLength /= CLHEP::cm; // now in G4e units (cm)
264  m_OffsetZ /= CLHEP::cm; // now in G4e units (cm)
265  m_BarrelMinR = bklmGeometry->getGap1InnerRadius(); // in G4e units (cm)
266  m_BarrelMaxR /= CLHEP::cm; // now in G4e units (cm)
267  m_EndcapMinR = eklmGeometry.getSectionPosition()->getInnerR() / CLHEP::cm; // in G4e units (cm)
268  m_EndcapMaxR = eklmGeometry.getSectionPosition()->getOuterR() / CLHEP::cm; // in G4e units (cm)
269  m_EndcapMiddleZ = m_BarrelHalfLength + m_EndcapHalfLength; // in G4e units (cm)
270 
271  // Measurement uncertainties and acceptance windows
272  double width = eklmGeometry.getStripGeometry()->getWidth() / CLHEP::cm; // in G4e units (cm)
273  m_EndcapScintVariance = width * width / 12.0;
274  width = bklmGeometry->getScintHalfWidth() * 2.0; // in G4e units (cm)
275  m_BarrelScintVariance = width * width / 12.0;
276  int nBarrelLayers = bklmGeometry->getNLayer();
277  for (int layer = 1; layer <= nBarrelLayers; ++layer) {
278  const bklm::Module* module =
279  bklmGeometry->findModule(BKLMElementNumbers::c_ForwardSection, 1, layer);
280  width = module->getPhiStripWidth(); // in G4e units (cm)
281  m_BarrelPhiStripVariance[layer - 1] = width * width / 12.0;
282  width = module->getZStripWidth(); // in G4e units (cm)
283  m_BarrelZStripVariance[layer - 1] = width * width / 12.0;
284  }
285 
286  // KLM geometry (for associating KLM hit with extrapolated crossing point)
287  m_OutermostActiveBarrelLayer = nBarrelLayers - 1; // zero-based counting
288  for (int section = 0; section <= BKLMElementNumbers::getMaximalSectionNumber(); ++section) {
289  for (int sector = 1; sector <= m_BarrelNSector; ++sector) {
290  for (int layer = 1; layer <= nBarrelLayers; ++layer) {
291  m_BarrelModuleMiddleRadius[section][sector - 1][layer - 1] =
292  bklmGeometry->getActiveMiddleRadius(section, sector, layer); // in G4e units (cm)
293  }
294  }
295  }
296  double dz(eklmGeometry.getLayerShiftZ() / CLHEP::cm); // in G4e units (cm)
297  double z0((eklmGeometry.getSectionPosition()->getZ()
298  + eklmGeometry.getLayerShiftZ()
299  - 0.5 * eklmGeometry.getSectionPosition()->getLength()
300  - 0.5 * eklmGeometry.getLayerPosition()->getLength()) / CLHEP::cm); // in G4e units (cm)
301  int nEndcapLayers = eklmGeometry.getNLayers();
302  m_OutermostActiveForwardEndcapLayer = eklmGeometry.getNDetectorLayers(2) - 1; // zero-based counting
303  m_OutermostActiveBackwardEndcapLayer = eklmGeometry.getNDetectorLayers(1) - 1; // zero-based counting
304 
305  for (int layer = 1; layer <= nEndcapLayers; ++layer) {
306  m_EndcapModuleMiddleZ[layer - 1] = z0 + dz * (layer - 1); // in G4e units (cm)
307  }
308  for (int sector = 1; sector <= 8; ++sector) {
309  double phi = M_PI_4 * (sector - 1);
310  m_BarrelSectorPerp[sector - 1].set(std::cos(phi), std::sin(phi), 0.0);
311  m_BarrelSectorPhi[sector - 1].set(-std::sin(phi), std::cos(phi), 0.0);
312  }
313 
314  m_eklmTransformData = &(EKLM::TransformDataGlobalAligned::Instance());
315 }
316 
317 void TrackExtrapolateG4e::beginRun(bool byMuid)
318 {
319  B2DEBUG(20, (byMuid ? "muid" : "ext"));
320  if (byMuid) {
321  if (!m_klmChannelStatus.isValid())
322  B2FATAL("KLM channel status data are not available.");
323  if (!m_klmStripEfficiency.isValid())
324  B2FATAL("KLM strip efficiency data are not available.");
325  if (!m_klmLikelihoodParameters.isValid())
326  B2FATAL("KLM likelihood parameters are not available.");
327  std::vector<int> muidPdgCodes = MuidElementNumbers::getPDGVector();
328  if (!m_MuidBuilderMap.empty()) {
329  if (m_klmLikelihoodParameters.hasChanged()) { /* Clear m_MuidBuilderMap if KLMLikelihoodParameters payload changed. */
330  for (auto const& [pdg, muidBuilder] : m_MuidBuilderMap)
331  delete muidBuilder;
332  m_MuidBuilderMap.clear();
333  } else /* Return if m_MuidBuilderMap is already initialized. */
334  return;
335  }
336  for (int pdg : muidPdgCodes)
337  m_MuidBuilderMap.insert(std::pair<int, MuidBuilder*>(pdg, new MuidBuilder(pdg)));
338  }
339 }
340 
341 void TrackExtrapolateG4e::event(bool byMuid)
342 {
343 
344  // Put geant4 in proper state (in case this module is in a separate process)
345  if (G4StateManager::GetStateManager()->GetCurrentState() == G4State_Idle) {
346  G4StateManager::GetStateManager()->SetNewState(G4State_GeomClosed);
347  }
348 
349  G4ThreeVector directionAtIP, positionG4e, momentumG4e;
350  G4ErrorTrajErr covG4e(5); // initialized to zeroes
351 
352  // Loop over the reconstructed tracks
353  // Do extrapolation for each hypothesis of each reconstructed track.
354  if (byMuid) { // event() called by Muid module
355  G4ErrorPropagatorData::GetErrorPropagatorData()->SetTarget(m_TargetMuid);
356  std::vector<std::pair<ECLCluster*, G4ThreeVector> > eclClusterInfo(m_eclClusters.getEntries());
357  for (int c = 0; c < m_eclClusters.getEntries(); ++c) {
358  eclClusterInfo[c].first = m_eclClusters[c];
359  eclClusterInfo[c].second = G4ThreeVector(m_eclClusters[c]->getClusterPosition().X(),
360  m_eclClusters[c]->getClusterPosition().Y(),
361  m_eclClusters[c]->getClusterPosition().Z()) * CLHEP::cm;
362  }
363  std::vector<std::pair<KLMCluster*, G4ThreeVector> > klmClusterInfo(m_klmClusters.getEntries());
364  for (int c = 0; c < m_klmClusters.getEntries(); ++c) {
365  klmClusterInfo[c].first = m_klmClusters[c];
366  klmClusterInfo[c].second = G4ThreeVector(m_klmClusters[c]->getClusterPosition().X(),
367  m_klmClusters[c]->getClusterPosition().Y(),
368  m_klmClusters[c]->getClusterPosition().Z()) * CLHEP::cm;
369  }
370  // Keep track of (re-)use of BKLMHit2ds
371  std::vector<std::map<const Track*, double> > bklmHitUsed(m_bklmHit2ds.getEntries());
372  for (auto& b2track : m_tracks) {
373  for (const auto& hypothesis : *m_HypothesesMuid) {
374  int pdgCode = hypothesis.getPDGCode();
375  if (hypothesis == Const::electron || hypothesis == Const::muon) pdgCode = -pdgCode;
376  G4ErrorFreeTrajState g4eState("g4e_mu+", G4ThreeVector(), G4ThreeVector()); // will be updated
377  ExtState extState = getStartPoint(b2track, pdgCode, g4eState);
378  swim(extState, g4eState, &eclClusterInfo, &klmClusterInfo, &bklmHitUsed);
379  } // Muid hypothesis loop
380  } // Muid track loop
381  } else { // event() called by Ext module
382  G4ErrorPropagatorData::GetErrorPropagatorData()->SetTarget(m_TargetExt);
383  for (auto& b2track : m_tracks) {
384  for (const auto& hypothesis : *m_HypothesesExt) {
385  int pdgCode = hypothesis.getPDGCode();
386  if (hypothesis == Const::electron || hypothesis == Const::muon) pdgCode = -pdgCode;
387  G4ErrorFreeTrajState g4eState("g4e_mu+", G4ThreeVector(), G4ThreeVector()); // will be updated
388  ExtState extState = getStartPoint(b2track, pdgCode, g4eState);
389  swim(extState, g4eState);
390  } // Ext hypothesis loop
391  } // Ext track loop
392  } // byMuid
393 
394 }
395 
396 void TrackExtrapolateG4e::endRun(bool)
397 {
398 }
399 
400 void TrackExtrapolateG4e::terminate(bool byMuid)
401 {
402  if (m_DefaultHypotheses != nullptr) { delete m_DefaultHypotheses; }
403  if (byMuid) {
404  delete m_TargetMuid;
405  for (auto const& [pdg, muidBuilder] : m_MuidBuilderMap)
406  delete muidBuilder;
407  }
408  if (m_TargetExt != nullptr) {
409  delete m_TargetExt;
410  m_TargetExt = nullptr;
411  }
412  if (m_EnterExit != nullptr) {
413  delete m_EnterExit;
414  delete m_BKLMVolumes;
415  m_ExtMgr->RunTermination();
416  m_EnterExit = nullptr;
417  m_BKLMVolumes = nullptr;
418  }
419 }
420 
421 void TrackExtrapolateG4e::extrapolate(int pdgCode, // signed for charge
422  double tof, // in ns (from IP to position)
423  // DIVOT bool isCosmic, // true for back-extrapolation
424  const G4ThreeVector& position, // in cm (genfit2 units)
425  const G4ThreeVector& momentum, // in GeV/c (genfit2 units)
426  const G4ErrorSymMatrix& covariance) // (6x6) using cm, GeV/c (genfit2 units)
427 {
428  bool isCosmic = false; // DIVOT
429  if ((!m_ExtInitialized) && (!m_MuidInitialized)) {
430  // No EXT nor MUID module in analysis path ==> mimic ext::initialize() with reasonable defaults.
431  // The default values are taken from the EXT module's parameter definitions.
432  Simulation::ExtManager* extMgr = Simulation::ExtManager::GetManager();
433  extMgr->Initialize("Ext", "default", 0.0, 0.25, false, 0, vector<string>());
434  // Redefine geant4e step length, magnetic field step limitation (fraction of local curvature radius),
435  // and kinetic energy loss limitation (maximum fractional energy loss) by communicating with
436  // the geant4 UI. (Commands were defined in ExtMessenger when physics list was set up.)
437  // *NOTE* If module muid runs after this, its G4UImanager commands will override these.
438  G4UImanager::GetUIpointer()->ApplyCommand("/geant4e/limits/stepLength 250 mm");
439  G4UImanager::GetUIpointer()->ApplyCommand("/geant4e/limits/magField 0.001");
440  G4UImanager::GetUIpointer()->ApplyCommand("/geant4e/limits/energyLoss 0.05");
441  m_DefaultHypotheses = new std::vector<Const::ChargedStable>; // not used
442  initialize(0.1, 0.002, *m_DefaultHypotheses);
443  }
444 
445  // Put geant4 in proper state (in case this module is in a separate process)
446  if (G4StateManager::GetStateManager()->GetCurrentState() == G4State_Idle) {
447  G4StateManager::GetStateManager()->SetNewState(G4State_GeomClosed);
448  }
449 
450  G4ErrorPropagatorData::GetErrorPropagatorData()->SetTarget(m_TargetExt);
451 
452  // Do extrapolation for selected hypothesis (pion, electron, muon, kaon, proton,
453  // deuteron) for the selected track until calorimeter exit.
454 
455  G4ThreeVector positionG4e = position * CLHEP::cm; // convert from genfit2 units (cm) to geant4 units (mm)
456  G4ThreeVector momentumG4e = momentum * CLHEP::GeV; // convert from genfit2 units (GeV/c) to geant4 units (MeV/c)
457  // cppcheck-suppress knownConditionTrueFalse
458  if (isCosmic) momentumG4e = -momentumG4e;
459  G4ErrorSymMatrix covarianceG4e(5, 0); // in Geant4e units (GeV/c, cm)
460  fromPhasespaceToG4e(momentum, covariance, covarianceG4e);
461  G4String nameG4e("g4e_" + G4ParticleTable::GetParticleTable()->FindParticle(pdgCode)->GetParticleName());
462  G4ErrorFreeTrajState g4eState(nameG4e, positionG4e, momentumG4e, covarianceG4e);
463  ExtState extState = { nullptr, pdgCode, isCosmic, tof, 0.0, // for EXT and MUID
464  momentumG4e.unit(), 0.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, false // for MUID only
465  };
466  swim(extState, g4eState);
467 }
468 
469 void TrackExtrapolateG4e::identifyMuon(int pdgCode, // signed for charge
470  double tof, // in ns (from IP to position)
471  bool isCosmic, // true for back-extrapolation
472  const G4ThreeVector& position, // in cm (genfit2 units)
473  const G4ThreeVector& momentum, // in GeV/c (genfit2 units)
474  const G4ErrorSymMatrix& covariance) // (6x6) using cm, GeV/c (genfit2 units)
475 {
476 
477  if ((!m_ExtInitialized) && (!m_MuidInitialized)) {
478  // No EXT nor MUID module in analysis path ==> mimic ext::initialize() with reasonable defaults.
479  // The default values are taken from the MUID module's parameter definitions.
480  Simulation::ExtManager* extMgr = Simulation::ExtManager::GetManager();
481  extMgr->Initialize("Muid", "default", 0.0, 0.25, false, 0, vector<string>());
482  // Redefine geant4e step length, magnetic field step limitation (fraction of local curvature radius),
483  // and kinetic energy loss limitation (maximum fractional energy loss) by communicating with
484  // the geant4 UI. (Commands were defined in ExtMessenger when physics list was set up.)
485  // *NOTE* If module muid runs after this, its G4UImanager commands will override these.
486  G4UImanager::GetUIpointer()->ApplyCommand("/geant4e/limits/stepLength 250 mm");
487  G4UImanager::GetUIpointer()->ApplyCommand("/geant4e/limits/magField 0.001");
488  G4UImanager::GetUIpointer()->ApplyCommand("/geant4e/limits/energyLoss 0.05");
489  m_DefaultHypotheses = new std::vector<Const::ChargedStable>; // not used
490  initialize(0.0, 30.0, 3.5, 150.0, 100.0, 0.1, 0.002, false, *m_DefaultHypotheses);
491  }
492 
493  // Put geant4 in proper state (in case this module is in a separate process)
494  if (G4StateManager::GetStateManager()->GetCurrentState() == G4State_Idle) {
495  G4StateManager::GetStateManager()->SetNewState(G4State_GeomClosed);
496  }
497 
498  G4ErrorPropagatorData::GetErrorPropagatorData()->SetTarget(m_TargetMuid);
499 
500  // Do extrapolation for selected hypothesis (pion, electron, muon, kaon, proton,
501  // deuteron) for the selected track until calorimeter exit.
502 
503  G4ThreeVector positionG4e = position * CLHEP::cm; // from genfit2 units (cm) to geant4 units (mm)
504  G4ThreeVector momentumG4e = momentum * CLHEP::GeV; // from genfit2 units (GeV/c) to geant4 units (MeV/c)
505  if (isCosmic) momentumG4e = -momentumG4e;
506  G4ErrorSymMatrix covarianceG4e(5, 0); // in Geant4e units (GeV/c, cm)
507  fromPhasespaceToG4e(momentum, covariance, covarianceG4e);
508  G4String nameG4e("g4e_" + G4ParticleTable::GetParticleTable()->FindParticle(pdgCode)->GetParticleName());
509  G4ErrorFreeTrajState g4eState(nameG4e, positionG4e, momentumG4e, covarianceG4e);
510  ExtState extState = { nullptr, pdgCode, isCosmic, tof, 0.0, // for EXT and MUID
511  momentumG4e.unit(), 0.0, 0, 0, 0, -1, -1, -1, -1, 0, 0, false // for MUID only
512  };
513  swim(extState, g4eState, nullptr, nullptr, nullptr);
514 }
515 
516 // Swim one track for MUID until it stops or leaves the KLM-bounding cylinder
517 void TrackExtrapolateG4e::swim(ExtState& extState, G4ErrorFreeTrajState& g4eState,
518  const std::vector<std::pair<ECLCluster*, G4ThreeVector> >* eclClusterInfo,
519  const std::vector<std::pair<KLMCluster*, G4ThreeVector> >* klmClusterInfo,
520  std::vector<std::map<const Track*, double> >* bklmHitUsed)
521 {
522  if (extState.pdgCode == 0) return;
523  if (g4eState.GetMomentum().perp() <= m_MinPt) return;
524  if (m_TargetMuid->GetDistanceFromPoint(g4eState.GetPosition()) < 0.0) return;
525  G4ParticleDefinition* particle = G4ParticleTable::GetParticleTable()->FindParticle(extState.pdgCode);
526  double mass = particle->GetPDGMass();
527  double minPSq = (mass + m_MinKE) * (mass + m_MinKE) - mass * mass;
528  // Create structures for ECLCluster-Track matching
529  std::vector<double> eclClusterDistance;
530  std::vector<ExtHit> eclHit1, eclHit2, eclHit3;
531  if (eclClusterInfo != nullptr) {
532  eclClusterDistance.resize(eclClusterInfo->size(), 1.0E10); // "positive infinity"
533  ExtHit tempExtHit(extState.pdgCode, Const::EDetector::ECL, 0, EXT_FIRST,
534  extState.isCosmic, 0.0,
535  G4ThreeVector(), G4ThreeVector(), G4ErrorSymMatrix(6));
536  eclHit1.resize(eclClusterInfo->size(), tempExtHit);
537  eclHit2.resize(eclClusterInfo->size(), tempExtHit);
538  eclHit3.resize(eclClusterInfo->size(), tempExtHit);
539  }
540  // Create structures for KLMCluster-Track matching
541  std::vector<TrackClusterSeparation> klmHit;
542  if (klmClusterInfo != nullptr) {
543  klmHit.resize(klmClusterInfo->size()); // initialize each to huge distance
544  }
545  KLMMuidLikelihood* klmMuidLikelihood = m_klmMuidLikelihoods.appendNew(); // rest of this object will be filled later
546  klmMuidLikelihood->setPDGCode(extState.pdgCode);
547  if (extState.track != nullptr) { extState.track->addRelationTo(klmMuidLikelihood); }
548  G4ErrorMode propagationMode = (extState.isCosmic ? G4ErrorMode_PropBackwards : G4ErrorMode_PropForwards);
549  m_ExtMgr->InitTrackPropagation(propagationMode);
550  while (true) {
551  const G4int errCode = m_ExtMgr->PropagateOneStep(&g4eState, propagationMode);
552  G4Track* track = g4eState.GetG4Track();
553  const G4Step* step = track->GetStep();
554  const G4StepPoint* preStepPoint = step->GetPreStepPoint();
555  const G4StepPoint* postStepPoint = step->GetPostStepPoint();
556  G4TouchableHandle preTouch = preStepPoint->GetTouchableHandle();
557  G4VPhysicalVolume* pVol = preTouch->GetVolume();
558  const G4int preStatus = preStepPoint->GetStepStatus();
559  const G4int postStatus = postStepPoint->GetStepStatus();
560  G4ThreeVector pos = track->GetPosition(); // this is at postStepPoint
561  G4ThreeVector mom = track->GetMomentum(); // ditto
562  // Ignore the zero-length step by PropagateOneStep() at each boundary
563  if (extState.isCosmic) mom = -mom;
564  if (step->GetStepLength() > 0.0) {
565  double dt = step->GetDeltaTime();
566  double dl = step->GetStepLength() / track->GetMaterial()->GetRadlen();
567  if (preStatus == fGeomBoundary) { // first step in this volume?
568  if (m_TargetExt->GetDistanceFromPoint(pos) < 0.0) { // only hits outside ECL during muid
569  if (m_EnterExit->find(pVol) != m_EnterExit->end()) {
570  createExtHit(EXT_ENTER, extState, g4eState, preStepPoint, preTouch);
571  }
572  }
573  }
574  if (extState.isCosmic) {
575  extState.tof -= dt;
576  extState.length -= dl;
577  } else {
578  extState.tof += dt;
579  extState.length += dl;
580  }
581  if (postStatus == fGeomBoundary) { // last step in this volume?
582  if (m_TargetExt->GetDistanceFromPoint(pos) < 0.0) { // only hits outside ECL during muid
583  if (m_EnterExit->find(pVol) != m_EnterExit->end()) {
584  createExtHit(EXT_EXIT, extState, g4eState, postStepPoint, preTouch);
585  }
586  }
587  }
588  if (createMuidHit(extState, g4eState, klmMuidLikelihood, bklmHitUsed)) {
589  // Force geant4e to update its G4Track from the Kalman-updated state
590  m_ExtMgr->GetPropagator()->SetStepN(0);
591  }
592  if (eclClusterInfo != nullptr) {
593  for (unsigned int c = 0; c < eclClusterInfo->size(); ++c) {
594  G4ThreeVector eclPos((*eclClusterInfo)[c].second);
595  G4ThreeVector prePos(preStepPoint->GetPosition());
596  G4ThreeVector diff(prePos - eclPos);
597  double distance = diff.mag();
598  if (distance < m_MaxECLTrackClusterDistance) {
599  // fallback ECLNEAR in case no ECLCROSS is found
600  if (distance < eclClusterDistance[c]) {
601  eclClusterDistance[c] = distance;
602  G4ErrorSymMatrix covariance(6, 0);
603  fromG4eToPhasespace(g4eState, covariance);
604  eclHit3[c].update(EXT_ECLNEAR, extState.tof, pos / CLHEP::cm, mom / CLHEP::GeV, covariance);
605  }
606  // find position of crossing of the track with the ECLCluster's sphere
607  if (eclHit1[c].getStatus() == EXT_FIRST) {
608  if (pos.mag2() >= eclPos.mag2()) {
609  double r = eclPos.mag();
610  double preD = prePos.mag() - r;
611  double postD = pos.mag() - r;
612  double f = postD / (postD - preD);
613  G4ThreeVector midPos = pos + (prePos - pos) * f;
614  double tof = extState.tof + dt * f * (extState.isCosmic ? +1 : -1); // in ns, at end of step
615  G4ErrorSymMatrix covariance(6, 0);
616  fromG4eToPhasespace(g4eState, covariance);
617  eclHit1[c].update(EXT_ECLCROSS, tof, midPos / CLHEP::cm, mom / CLHEP::GeV, covariance);
618  }
619  }
620  }
621  // find closest distance to the radial line to the ECLCluster
622  if (eclHit2[c].getStatus() == EXT_FIRST) {
623  G4ThreeVector delta(pos - prePos);
624  G4ThreeVector perp(eclPos.cross(delta));
625  double perpMag2 = perp.mag2();
626  if (perpMag2 > 1.0E-10) {
627  double dist = std::fabs(diff * perp) / std::sqrt(perpMag2);
628  if (dist < m_MaxECLTrackClusterDistance) {
629  double f = eclPos * (prePos.cross(perp)) / perpMag2;
630  if ((f > -0.5) && (f <= 1.0)) {
631  G4ThreeVector midPos(prePos + f * delta);
632  double length = extState.length + dl * (1.0 - f) * (extState.isCosmic ? +1 : -1);
633  G4ErrorSymMatrix covariance(6, 0);
634  fromG4eToPhasespace(g4eState, covariance);
635  eclHit2[c].update(EXT_ECLDL, length, midPos / CLHEP::cm, mom / CLHEP::GeV, covariance);
636  }
637  }
638  }
639  }
640  }
641  }
642  if (klmClusterInfo != nullptr) {
643  for (unsigned int c = 0; c < klmClusterInfo->size(); ++c) {
644  G4ThreeVector klmPos = (*klmClusterInfo)[c].second;
645  G4ThreeVector separation = klmPos - pos;
646  double distance = separation.mag();
647  if (distance < klmHit[c].getDistance()) {
648  klmHit[c].setDistance(distance);
649  klmHit[c].setTrackClusterAngle(mom.angle(separation));
650  klmHit[c].setTrackClusterSeparationAngle(mom.angle(klmPos));
651  klmHit[c].setTrackRotationAngle(extState.directionAtIP.angle(mom));
652  klmHit[c].setTrackClusterInitialSeparationAngle(extState.directionAtIP.angle(klmPos));
653  }
654  }
655  }
656  }
657  // Post-step momentum too low?
658  if (errCode || (mom.mag2() < minPSq)) {
659  break;
660  }
661  // Detect escapes from the imaginary target cylinder.
662  if (m_TargetMuid->GetDistanceFromPoint(pos) < 0.0) {
663  break;
664  }
665  // Stop extrapolating as soon as the track curls inward too much
666  if (pos.perp2() < m_MinRadiusSq) {
667  break;
668  }
669  } // track-extrapolation "infinite" loop
670 
671  m_ExtMgr->EventTermination(propagationMode);
672 
673  finishTrack(extState, klmMuidLikelihood, (g4eState.GetPosition().z() > m_OffsetZ));
674 
675  if (eclClusterInfo != nullptr) {
676  for (unsigned int c = 0; c < eclClusterInfo->size(); ++c) {
677  if (eclHit1[c].getStatus() != EXT_FIRST) {
678  ExtHit* h = m_extHits.appendNew(eclHit1[c]);
679  (*eclClusterInfo)[c].first->addRelationTo(h);
680  extState.track->addRelationTo(h);
681  }
682  if (eclHit2[c].getStatus() != EXT_FIRST) {
683  ExtHit* h = m_extHits.appendNew(eclHit2[c]);
684  (*eclClusterInfo)[c].first->addRelationTo(h);
685  extState.track->addRelationTo(h);
686  }
687  if (eclHit3[c].getStatus() != EXT_FIRST) {
688  ExtHit* h = m_extHits.appendNew(eclHit3[c]);
689  (*eclClusterInfo)[c].first->addRelationTo(h);
690  extState.track->addRelationTo(h);
691  }
692  }
693  }
694 
695  if (klmClusterInfo != nullptr) {
696  // here we set a relation only to the closest KLMCluster
697  // and we don't set any relation if the distance is too large
698  double minDistance = m_MaxKLMTrackClusterDistance;
699  unsigned int closestCluster = 0;
700  for (unsigned int c = 0; c < klmClusterInfo->size(); ++c) {
701  if (klmHit[c].getDistance() > 1.0E9) {
702  continue;
703  }
704  TrackClusterSeparation* h = m_trackClusterSeparations.appendNew(klmHit[c]);
705  (*klmClusterInfo)[c].first->addRelationTo(h); // relation KLMCluster to TrackSep
706  extState.track->addRelationTo(h); // relation Track to TrackSep
707  if (klmHit[c].getDistance() < minDistance) {
708  closestCluster = c;
709  minDistance = klmHit[c].getDistance();
710  }
711  }
712  if (minDistance < m_MaxKLMTrackClusterDistance) {
713  // set the relation Track to KLMCluster, using the distance as weight
714  extState.track->addRelationTo((*klmClusterInfo)[closestCluster].first, 1. / minDistance);
715  }
716  }
717 
718 }
719 
720 // Swim one track for EXT until it stops or leaves the ECL-bounding cylinder
721 void TrackExtrapolateG4e::swim(ExtState& extState, G4ErrorFreeTrajState& g4eState)
722 {
723  if (extState.pdgCode == 0) return;
724  if (g4eState.GetMomentum().perp() <= m_MinPt) return;
725  if (m_TargetExt->GetDistanceFromPoint(g4eState.GetPosition()) < 0.0) return;
726  G4ParticleDefinition* particle = G4ParticleTable::GetParticleTable()->FindParticle(extState.pdgCode);
727  double mass = particle->GetPDGMass();
728  double minPSq = (mass + m_MinKE) * (mass + m_MinKE) - mass * mass;
729  G4ErrorMode propagationMode = (extState.isCosmic ? G4ErrorMode_PropBackwards : G4ErrorMode_PropForwards);
730  m_ExtMgr->InitTrackPropagation(propagationMode);
731  while (true) {
732  const G4int errCode = m_ExtMgr->PropagateOneStep(&g4eState, propagationMode);
733  G4Track* track = g4eState.GetG4Track();
734  const G4Step* step = track->GetStep();
735  const G4StepPoint* preStepPoint = step->GetPreStepPoint();
736  const G4StepPoint* postStepPoint = step->GetPostStepPoint();
737  G4TouchableHandle preTouch = preStepPoint->GetTouchableHandle();
738  G4VPhysicalVolume* pVol = preTouch->GetVolume();
739  const G4int preStatus = preStepPoint->GetStepStatus();
740  const G4int postStatus = postStepPoint->GetStepStatus();
741  G4ThreeVector pos = track->GetPosition(); // this is at postStepPoint
742  G4ThreeVector mom = track->GetMomentum(); // ditto
743  // First step on this track?
744  if (extState.isCosmic) mom = -mom;
745  if (preStatus == fUndefined) {
746  if (m_EnterExit->find(pVol) != m_EnterExit->end()) {
747  createExtHit(EXT_FIRST, extState, g4eState, preStepPoint, preTouch);
748  }
749  }
750  // Ignore the zero-length step by PropagateOneStep() at each boundary
751  if (step->GetStepLength() > 0.0) {
752  double dt = step->GetDeltaTime();
753  double dl = step->GetStepLength() / track->GetMaterial()->GetRadlen();
754  if (preStatus == fGeomBoundary) { // first step in this volume?
755  if (m_EnterExit->find(pVol) != m_EnterExit->end()) {
756  createExtHit(EXT_ENTER, extState, g4eState, preStepPoint, preTouch);
757  }
758  }
759  if (extState.isCosmic) {
760  extState.tof -= dt;
761  extState.length -= dl;
762  } else {
763  extState.tof += dt;
764  extState.length += dl;
765  }
766  // Last step in this volume?
767  if (postStatus == fGeomBoundary) {
768  if (m_EnterExit->find(pVol) != m_EnterExit->end()) {
769  createExtHit(EXT_EXIT, extState, g4eState, postStepPoint, preTouch);
770  }
771  }
772  }
773  // Post-step momentum too low?
774  if (errCode || (mom.mag2() < minPSq)) {
775  if (m_EnterExit->find(pVol) != m_EnterExit->end()) {
776  createExtHit(EXT_STOP, extState, g4eState, postStepPoint, preTouch);
777  }
778  break;
779  }
780  // Detect escapes from the imaginary target cylinder.
781  if (m_TargetExt->GetDistanceFromPoint(pos) < 0.0) {
782  if (m_EnterExit->find(pVol) != m_EnterExit->end()) {
783  createExtHit(EXT_ESCAPE, extState, g4eState, postStepPoint, preTouch);
784  }
785  break;
786  }
787  // Stop extrapolating as soon as the track curls inward too much
788  if (pos.perp2() < m_MinRadiusSq) {
789  break;
790  }
791  } // track-extrapolation "infinite" loop
792 
793  m_ExtMgr->EventTermination(propagationMode);
794 
795 }
796 
797 // Register the list of volumes for which entry/exit point is to be saved during extrapolation
798 void TrackExtrapolateG4e::registerVolumes()
799 {
800 
801  G4PhysicalVolumeStore* pvStore = G4PhysicalVolumeStore::GetInstance();
802  if (pvStore->size() == 0) {
803  B2FATAL("No geometry defined. Please create the geometry first.");
804  }
805 
806  if (m_EnterExit != nullptr) { return; } // Only do this once
807 
808  m_EnterExit = new map<G4VPhysicalVolume*, enum VolTypes>;
809  m_BKLMVolumes = new vector<G4VPhysicalVolume*>;
810  for (vector<G4VPhysicalVolume*>::iterator iVol = pvStore->begin();
811  iVol != pvStore->end(); ++iVol) {
812  const G4String name = (*iVol)->GetName();
813 
814  // CDC volumes have "CDC" in the name
815  if (name.find("CDC") != string::npos) {
816  // DO NOT STORE CDC HITS
817  // (*m_EnterExit)[*iVol] = VOLTYPE_CDC;
818  }
819  // TOP doesn't have one envelope; it has 16 "TOPModule"s
820  else if (name.find("TOPModule") != string::npos) {
821  (*m_EnterExit)[*iVol] = VOLTYPE_TOP1;
822  }
823  // TOP quartz bar (=sensitive)
824  else if (name.find("_TOPPrism_") != string::npos ||
825  name.find("_TOPBarSegment") != string::npos ||
826  name.find("_TOPMirrorSegment") != string::npos) {
827  (*m_EnterExit)[*iVol] = VOLTYPE_TOP2;
828  }
829  // TOP quartz glue (not sensitive?)
830  else if (name.find("TOPBarSegment1Glue") != string::npos ||
831  name.find("TOPBarSegment2Glue") != string::npos ||
832  name.find("TOPMirrorSegmentGlue") != string::npos) {
833  (*m_EnterExit)[*iVol] = VOLTYPE_TOP3;
834  // ARICH volumes
835  } else if (name == "ARICH.AerogelSupportPlate") {
836  (*m_EnterExit)[*iVol] = VOLTYPE_ARICH1;
837  } else if (name == "ARICH.AerogelImgPlate") {
838  (*m_EnterExit)[*iVol] = VOLTYPE_ARICH2;
839  } else if (name.find("ARICH.HAPDWindow") != string::npos) {
840  (*m_EnterExit)[*iVol] = VOLTYPE_ARICH3;
841  }
842 
843  // ECL crystal
844  else if (name.find("lv_barrel_crystal_") != string::npos ||
845  name.find("lv_forward_crystal_") != string::npos ||
846  name.find("lv_backward_crystal_") != string::npos) {
847  (*m_EnterExit)[*iVol] = VOLTYPE_ECL;
848  }
849  // Barrel KLM: BKLM.Layer**GasPhysical for RPCs or BKLM.Layer**ChimneyGasPhysical for RPCs
850  // BKLM.ScintActiveType*Physical for scintillator strips
851  else if (name.compare(0, 5, "BKLM.") == 0) {
852  if (name.find("GasPhysical") != string::npos) {
853  (*m_EnterExit)[*iVol] = VOLTYPE_BKLM1;
854  } else if (name.find("ScintActiveType") != string::npos) {
855  (*m_EnterExit)[*iVol] = VOLTYPE_BKLM2;
856  } else if ((name.find("ScintType") != string::npos) ||
857  (name.find("ElectrodePhysical") != string::npos)) {
858  m_BKLMVolumes->push_back(*iVol);
859  }
860  }
861  // Endcap KLM: StripSensitive_*
862  else if (name.compare(0, 14, "StripSensitive") == 0) {
863  (*m_EnterExit)[*iVol] = VOLTYPE_EKLM;
864  }
865  }
866 
867 }
868 
869 // Convert the physical volume to integer(-like) identifiers
870 void TrackExtrapolateG4e::getVolumeID(const G4TouchableHandle& touch, Const::EDetector& detID, int& copyID)
871 {
872 
873  // default values
874  detID = Const::EDetector::invalidDetector;
875  copyID = 0;
876 
877  G4VPhysicalVolume* pv = touch->GetVolume(0);
878  std::map<G4VPhysicalVolume*, enum VolTypes>::iterator it = m_EnterExit->find(pv);
879  if (it == m_EnterExit->end()) { return; }
880 
881  switch (it->second) {
882  case VOLTYPE_CDC:
883  detID = Const::EDetector::CDC;
884  copyID = pv->GetCopyNo();
885  return;
886  case VOLTYPE_TOP1:
887  detID = Const::EDetector::TOP;
888  copyID = -(pv->GetCopyNo()); // negative to distinguish module and quartz hits
889  return;
890  case VOLTYPE_TOP2:
891  detID = Const::EDetector::TOP;
892  if (touch->GetHistoryDepth() >= 1) copyID = touch->GetVolume(1)->GetCopyNo();
893  return;
894  case VOLTYPE_TOP3:
895  detID = Const::EDetector::TOP;
896  if (touch->GetHistoryDepth() >= 2) copyID = touch->GetVolume(2)->GetCopyNo();
897  return;
898  case VOLTYPE_ARICH1:
899  detID = Const::EDetector::ARICH;
900  copyID = 12345;
901  return;
902  case VOLTYPE_ARICH2:
903  detID = Const::EDetector::ARICH;
904  copyID = 6789;
905  return;
906  case VOLTYPE_ARICH3:
907  detID = Const::EDetector::ARICH;
908  if (touch->GetHistoryDepth() >= 2) copyID = touch->GetVolume(2)->GetCopyNo();
909  return;
910  case VOLTYPE_ECL:
911  detID = Const::EDetector::ECL;
912  copyID = ECL::ECLGeometryPar::Instance()->ECLVolumeToCellID(touch());
913  return;
914  case VOLTYPE_BKLM1: // BKLM RPCs
915  detID = Const::EDetector::BKLM;
916  if (touch->GetHistoryDepth() == DEPTH_RPC) {
917  // int plane = touch->GetCopyNumber(0);
918  int layer = touch->GetCopyNumber(4);
919  int sector = touch->GetCopyNumber(6);
920  int section = touch->GetCopyNumber(7);
921  copyID = BKLMElementNumbers::moduleNumber(section, sector, layer);
922  }
923  return;
924  case VOLTYPE_BKLM2: // BKLM scints
925  detID = Const::EDetector::BKLM;
926  if (touch->GetHistoryDepth() == DEPTH_SCINT) {
927  int strip = touch->GetCopyNumber(1);
928  int plane = (touch->GetCopyNumber(2) == BKLM_INNER) ?
929  BKLMElementNumbers::c_PhiPlane :
930  BKLMElementNumbers::c_ZPlane;
931  int layer = touch->GetCopyNumber(6);
932  int sector = touch->GetCopyNumber(8);
933  int section = touch->GetCopyNumber(9);
934  copyID = BKLMElementNumbers::channelNumber(
935  section, sector, layer, plane, strip);
936  BKLMStatus::setMaximalStrip(copyID, strip);
937  }
938  return;
939  case VOLTYPE_EKLM:
940  detID = Const::EDetector::EKLM;
941  copyID = EKLMElementNumbers::Instance().stripNumber(
942  touch->GetVolume(7)->GetCopyNo(),
943  touch->GetVolume(6)->GetCopyNo(),
944  touch->GetVolume(5)->GetCopyNo(),
945  touch->GetVolume(4)->GetCopyNo(),
946  touch->GetVolume(1)->GetCopyNo());
947  return;
948  }
949 
950 }
951 
952 ExtState TrackExtrapolateG4e::getStartPoint(const Track& b2track, int pdgCode, G4ErrorFreeTrajState& g4eState)
953 {
954  ExtState extState = {&b2track, pdgCode, false, 0.0, 0.0, // for EXT and MUID
955  G4ThreeVector(0, 0, 1), 0.0, 0, 0, 0, -1, -1, -1, -1, 0, 0, false // for MUID only
956  };
957  RecoTrack* recoTrack = b2track.getRelatedTo<RecoTrack>();
958  if (recoTrack == nullptr) {
959  B2WARNING("Track without associated RecoTrack: skipping extrapolation for this track.");
960  extState.pdgCode = 0; // prevent start of extrapolation in swim()
961  return extState;
962  }
963  const genfit::AbsTrackRep* trackRep = recoTrack->getCardinalRepresentation();
964  int charge = int(trackRep->getPDGCharge());
965  if (charge != 0) {
966  extState.pdgCode *= charge;
967  } else {
968  charge = 1; // should never happen but persist if it does
969  }
970  TVector3 firstPosition, firstMomentum, lastPosition, lastMomentum; // initialized to zeroes
971  TMatrixDSym firstCov(6), lastCov(6); // initialized to zeroes
972  try {
973  const genfit::MeasuredStateOnPlane& firstState = recoTrack->getMeasuredStateOnPlaneFromFirstHit(trackRep);
974  trackRep->getPosMomCov(firstState, firstPosition, firstMomentum, firstCov);
975  const genfit::MeasuredStateOnPlane& lastState = recoTrack->getMeasuredStateOnPlaneFromLastHit(trackRep);
976  trackRep->getPosMomCov(lastState, lastPosition, lastMomentum, lastCov);
977  // in genfit units (cm, GeV/c)
978  extState.tof = lastState.getTime(); // DIVOT: must be revised when IP profile (reconstructed beam spot) become available!
979  if (lastPosition.Mag2() < firstPosition.Mag2()) {
980  firstPosition = lastPosition;
981  firstMomentum = -lastMomentum;
982  firstCov = lastCov;
983  trackRep->getPosMomCov(firstState, lastPosition, lastMomentum, lastCov);
984  lastMomentum *= -1.0; // extrapolate backwards instead of forwards
985  extState.isCosmic = true;
986  extState.tof = firstState.getTime(); // DIVOT: must be revised when IP profile (reconstructed beam spot) become available!
987  }
988 
989  G4ParticleDefinition* particle = G4ParticleTable::GetParticleTable()->FindParticle(extState.pdgCode);
990  if (extState.pdgCode != trackRep->getPDG()) {
991  double pSq = lastMomentum.Mag2();
992  double mass = particle->GetPDGMass() / CLHEP::GeV;
993  extState.tof *= std::sqrt((pSq + mass * mass) / (pSq + lastState.getMass() * lastState.getMass()));
994  }
995 
996  extState.directionAtIP.set(firstMomentum.Unit().X(), firstMomentum.Unit().Y(), firstMomentum.Unit().Z());
997  if (m_MagneticField != 0.0) { // in gauss
998  double radius = (firstMomentum.Perp() * CLHEP::GeV / CLHEP::eV) / (CLHEP::c_light * charge * m_MagneticField); // in cm
999  double centerPhi = extState.directionAtIP.phi() - M_PI_2;
1000  double centerX = firstPosition.X() + radius * std::cos(centerPhi);
1001  double centerY = firstPosition.Y() + radius * std::sin(centerPhi);
1002  double pocaPhi = atan2(charge * centerY, charge * centerX) + M_PI;
1003  double ipPerp = extState.directionAtIP.perp();
1004  if (ipPerp > 0.0) {
1005  extState.directionAtIP.setX(+std::sin(pocaPhi) * ipPerp);
1006  extState.directionAtIP.setY(-std::cos(pocaPhi) * ipPerp);
1007  }
1008  } else {
1009  // No field: replace flaky covariance matrix with a diagonal one measured in 1.5T field
1010  // for a 10 GeV/c track ... and replace momentum magnitude with fixed 10 GeV/c
1011  lastCov *= 0.0;
1012  lastCov[0][0] = 5.0E-5;
1013  lastCov[1][1] = 1.0E-7;
1014  lastCov[2][2] = 5.0E-4;
1015  lastCov[3][3] = 3.5E-3;
1016  lastCov[4][4] = 3.5E-3;
1017  lastMomentum = lastMomentum.Unit() * 10.0;
1018  }
1019 
1020  G4ThreeVector posG4e(lastPosition.X() * CLHEP::cm, lastPosition.Y() * CLHEP::cm,
1021  lastPosition.Z() * CLHEP::cm); // in Geant4 units (mm)
1022  G4ThreeVector momG4e(lastMomentum.X() * CLHEP::GeV, lastMomentum.Y() * CLHEP::GeV,
1023  lastMomentum.Z() * CLHEP::GeV); // in Geant4 units (MeV/c)
1024  G4ErrorSymMatrix covG4e(5, 0); // in Geant4e units (GeV/c, cm)
1025  fromPhasespaceToG4e(lastMomentum, lastCov, covG4e); // in Geant4e units (GeV/c, cm)
1026  g4eState.SetData("g4e_" + particle->GetParticleName(), posG4e, momG4e);
1027  g4eState.SetParameters(posG4e, momG4e); // compute private-state parameters from momG4e
1028  g4eState.SetError(covG4e);
1029  } catch (const genfit::Exception&) {
1030  B2WARNING("genfit::MeasuredStateOnPlane() exception: skipping extrapolation for this track. initial momentum = ("
1031  << firstMomentum.X() << "," << firstMomentum.Y() << "," << firstMomentum.Z() << ")");
1032  extState.pdgCode = 0; // prevent start of extrapolation in swim()
1033  }
1034 
1035  return extState;
1036 }
1037 
1038 
1039 void TrackExtrapolateG4e::fromG4eToPhasespace(const G4ErrorFreeTrajState& g4eState, G4ErrorSymMatrix& covariance)
1040 {
1041 
1042  // Convert Geant4e covariance matrix with parameters 1/p, lambda, phi, yT, zT (in GeV/c, radians, cm)
1043  // to phase-space covariance matrix with parameters x, y, z, px, py, pz (in GeV/c, cm)
1044  // (1/p) = 1/sqrt( px^2 + py^2 + pz^2 )
1045  // phi = atan( py / px )
1046  // lambda = asin( pz / sqrt( px^2 + py^2 + pz^2 )
1047  // xT = x * cos(lambda) * cos(phi) + y * cos(lambda) * sin(phi) + z * sin(lambda)
1048  // yT = -x * sin(phi) + y * cos(phi)
1049  // zT = -x * sin(lambda) * cos(phi) - y * sin(lambda) * sin(phi) + z * cos(lambda)
1050 
1051  G4ErrorFreeTrajParam param = g4eState.GetParameters();
1052  double p = 1.0 / (param.GetInvP() * CLHEP::GeV); // in GeV/c
1053  double pSq = p * p;
1054  double lambda = param.GetLambda(); // in radians
1055  double sinLambda = std::sin(lambda);
1056  double cosLambda = std::cos(lambda);
1057  double phi = param.GetPhi(); // in radians
1058  double sinPhi = std::sin(phi);
1059  double cosPhi = std::cos(phi);
1060 
1061  // Transformation Jacobian 6x5 from Geant4e 5x5 to phase-space 6x6
1062 
1063  G4ErrorMatrix jacobian(6, 5, 0); // All entries are initialized to 0
1064 
1065  jacobian(4, 1) = -pSq * cosLambda * cosPhi; // @(px)/@(1/p)
1066  jacobian(5, 1) = -pSq * cosLambda * sinPhi; // @(py)/@(1/p)
1067  jacobian(6, 1) = -pSq * sinLambda; // @(pz)/@(1/p)
1068 
1069  jacobian(4, 2) = -p * sinLambda * cosPhi; // @(px)/@(lambda)
1070  jacobian(5, 2) = -p * sinLambda * sinPhi; // @(py)/@(lambda)
1071  jacobian(6, 2) = p * cosLambda; // @(pz)/@(lambda)
1072 
1073  jacobian(4, 3) = -p * cosLambda * sinPhi; // @(px)/@(phi)
1074  jacobian(5, 3) = p * cosLambda * cosPhi; // @(py)/@(phi)
1075 
1076  jacobian(1, 4) = -sinPhi; // @(x)/@(yT)
1077  jacobian(2, 4) = cosPhi; // @(y)/@(yT)
1078 
1079  jacobian(1, 5) = -sinLambda * cosPhi; // @(x)/@(zT)
1080  jacobian(2, 5) = -sinLambda * sinPhi; // @(y)/@(zT)
1081  jacobian(3, 5) = cosLambda; // @(z)/@(zT)
1082 
1083  G4ErrorTrajErr g4eCov = g4eState.GetError();
1084  covariance.assign(g4eCov.similarity(jacobian));
1085 
1086 }
1087 
1088 void TrackExtrapolateG4e::fromPhasespaceToG4e(const TVector3& momentum, const TMatrixDSym& covariance, G4ErrorTrajErr& covG4e)
1089 {
1090 
1091  // Convert phase-space covariance matrix with parameters x, y, z, px, py, pz (in genfit2 units: cm, GeV/c)
1092  // to Geant4e covariance matrix with parameters 1/p, lambda, phi, yT, zT (in geant4e units: GeV/c, radians, cm)
1093  // xT = x * cos(lambda) * cos(phi) + y * cos(lambda) * sin(phi) + z * sin(lambda)
1094  // yT = -x * sin(phi) + y * cos(phi)
1095  // zT = -x * sin(lambda) * cos(phi) - y * sin(lambda) * sin(phi) + z * cos(lambda)
1096  // (1/p) = 1/sqrt( px^2 + py^2 + pz^2 )
1097  // phi = atan( py / px )
1098  // lambda = asin( pz / sqrt( px^2 + py^2 + pz^2 )
1099 
1100  G4ErrorSymMatrix temp(6, 0);
1101  for (int k = 0; k < 6; ++k) {
1102  for (int j = k; j < 6; ++j) {
1103  temp[j][k] = covariance[j][k];
1104  }
1105  }
1106 
1107  double pInvSq = 1.0 / momentum.Mag2();
1108  double pInv = std::sqrt(pInvSq);
1109  double pPerpInv = 1.0 / momentum.Perp();
1110  double sinLambda = momentum.CosTheta();
1111  double cosLambda = std::sqrt(1.0 - sinLambda * sinLambda);
1112  double phi = momentum.Phi();
1113  double cosPhi = std::cos(phi);
1114  double sinPhi = std::sin(phi);
1115 
1116  // Transformation Jacobian 5x6 from phase-space 6x6 to Geant4e 5x5
1117  G4ErrorMatrix jacobian(5, 6, 0); // All entries are initialized to 0
1118 
1119  jacobian(1, 4) = -pInvSq * cosLambda * cosPhi; // @(1/p)/@(px)
1120  jacobian(1, 5) = -pInvSq * cosLambda * sinPhi; // @(1/p)/@(py)
1121  jacobian(1, 6) = -pInvSq * sinLambda; // @(1/p)/@(pz)
1122 
1123  jacobian(2, 4) = -pInv * sinLambda * cosPhi; // @(lambda)/@(px)
1124  jacobian(2, 5) = -pInv * sinLambda * sinPhi; // @(lambda)/@(py)
1125  jacobian(2, 6) = pInv * cosLambda; // @(lambda)/@(pz)
1126 
1127  jacobian(3, 4) = -pPerpInv * sinPhi; // @(phi)/@(px)
1128  jacobian(3, 5) = pPerpInv * cosPhi; // @(phi)/@(py)
1129 
1130  jacobian(4, 1) = -sinPhi; // @(yT)/@(x)
1131  jacobian(4, 2) = cosPhi; // @(yT)/@(y)
1132 
1133  jacobian(5, 1) = -sinLambda * cosPhi; // @(zT)/@(x)
1134  jacobian(5, 2) = -sinLambda * sinPhi; // @(zT)/@(y)
1135  jacobian(5, 3) = cosLambda; // @(zT)/@(z)
1136 
1137  covG4e = temp.similarity(jacobian);
1138 
1139 }
1140 
1141 void TrackExtrapolateG4e::fromPhasespaceToG4e(const G4ThreeVector& momentum, const G4ErrorSymMatrix& covariance,
1142  G4ErrorTrajErr& covG4e)
1143 {
1144 
1145  // Convert phase-space covariance matrix with parameters x, y, z, px, py, pz (in genfit2 units: cm, GeV/c)
1146  // to Geant4e covariance matrix with parameters 1/p, lambda, phi, yT, zT (in geant4e units: GeV/c, radians, cm)
1147  // xT = x * cos(lambda) * cos(phi) + y * cos(lambda) * sin(phi) + z * sin(lambda)
1148  // yT = -x * sin(phi) + y * cos(phi)
1149  // zT = -x * sin(lambda) * cos(phi) - y * sin(lambda) * sin(phi) + z * cos(lambda)
1150  // (1/p) = 1/sqrt( px^2 + py^2 + pz^2 )
1151  // phi = atan( py / px )
1152  // lambda = asin( pz / sqrt( px^2 + py^2 + pz^2 )
1153 
1154  G4ErrorSymMatrix temp(covariance);
1155 
1156  double pInvSq = 1.0 / momentum.mag2();
1157  double pInv = std::sqrt(pInvSq);
1158  double pPerpInv = 1.0 / momentum.perp();
1159  double sinLambda = momentum.cosTheta();
1160  double cosLambda = std::sqrt(1.0 - sinLambda * sinLambda);
1161  double phi = momentum.phi();
1162  double cosPhi = std::cos(phi);
1163  double sinPhi = std::sin(phi);
1164 
1165  // Transformation Jacobian 5x6 from phase-space 6x6 to Geant4e 5x5
1166  G4ErrorMatrix jacobian(5, 6, 0);
1167 
1168  jacobian(1, 4) = -pInvSq * cosLambda * cosPhi; // @(1/p)/@(px)
1169  jacobian(1, 5) = -pInvSq * cosLambda * sinPhi; // @(1/p)/@(py)
1170  jacobian(1, 6) = -pInvSq * sinLambda; // @(1/p)/@(pz)
1171 
1172  jacobian(2, 4) = -pInv * sinLambda * cosPhi; // @(lambda)/@(px)
1173  jacobian(2, 5) = -pInv * sinLambda * sinPhi; // @(lambda)/@(py)
1174  jacobian(2, 6) = pInv * cosLambda; // @(lambda)/@(pz)
1175 
1176  jacobian(3, 4) = -pPerpInv * sinPhi; // @(phi)/@(px)
1177  jacobian(3, 5) = pPerpInv * cosPhi; // @(phi)/@(py)
1178 
1179  jacobian(4, 1) = -sinPhi; // @(yT)/@(x)
1180  jacobian(4, 2) = cosPhi; // @(yT)/@(y)
1181 
1182  jacobian(5, 1) = -sinLambda * cosPhi; // @(zT)/@(x)
1183  jacobian(5, 2) = -sinLambda * sinPhi; // @(zT)/@(y)
1184  jacobian(5, 3) = cosLambda; // @(zT)/@(z)
1185  covG4e = temp.similarity(jacobian);
1186 
1187 }
1188 
1189 // write another volume-entry or volume-exit point on extrapolated track
1190 void TrackExtrapolateG4e::createExtHit(ExtHitStatus status, const ExtState& extState,
1191  const G4ErrorFreeTrajState& g4eState,
1192  const G4StepPoint* stepPoint, const G4TouchableHandle& touch)
1193 {
1194 
1195  Const::EDetector detID(Const::EDetector::invalidDetector);
1196  int copyID(0);
1197  getVolumeID(touch, detID, copyID);
1198  G4ThreeVector pos(stepPoint->GetPosition() / CLHEP::cm);
1199  G4ThreeVector mom(stepPoint->GetMomentum() / CLHEP::GeV);
1200  if (extState.isCosmic) mom = -mom;
1201  G4ErrorSymMatrix covariance(6, 0);
1202  fromG4eToPhasespace(g4eState, covariance);
1203  ExtHit* extHit = m_extHits.appendNew(extState.pdgCode, detID, copyID, status,
1204  extState.isCosmic, extState.tof,
1205  pos, mom, covariance);
1206  // If called standalone, there will be no associated track
1207  if (extState.track != nullptr) { extState.track->addRelationTo(extHit); }
1208 
1209 }
1210 
1211 // Write another volume-entry point on track.
1212 // The track state will be modified here by the Kalman fitter.
1213 
1214 bool TrackExtrapolateG4e::createMuidHit(ExtState& extState, G4ErrorFreeTrajState& g4eState, KLMMuidLikelihood* klmMuidLikelihood,
1215  std::vector<std::map<const Track*, double> >* bklmHitUsed)
1216 {
1217 
1218  Intersection intersection;
1219  intersection.hit = -1;
1220  intersection.chi2 = -1.0;
1221  intersection.position = g4eState.GetPosition() / CLHEP::cm;
1222  intersection.momentum = g4eState.GetMomentum() / CLHEP::GeV;
1223  G4ThreeVector prePos = g4eState.GetG4Track()->GetStep()->GetPreStepPoint()->GetPosition() / CLHEP::cm;
1224  G4ThreeVector oldPosition(prePos.x(), prePos.y(), prePos.z());
1225  double r = intersection.position.perp();
1226  double z = std::fabs(intersection.position.z() - m_OffsetZ);
1227 
1228  // Is the track in the barrel?
1229  if ((r > m_BarrelMinR) && (r < m_BarrelMaxR) && (z < m_BarrelHalfLength)) {
1230  // Did the track cross the inner midplane of a detector module?
1231  if (findBarrelIntersection(extState, oldPosition, intersection)) {
1232  fromG4eToPhasespace(g4eState, intersection.covariance);
1233  if (findMatchingBarrelHit(intersection, extState.track)) {
1234  (*bklmHitUsed)[intersection.hit].insert(std::pair<const Track*, double>(extState.track, intersection.chi2));
1235  extState.extLayerPattern |= (0x00000001 << intersection.layer);
1236  //efficiency implementation
1237  float phiBarrelEfficiency = m_klmStripEfficiency->getBarrelEfficiency((intersection.isForward ? 1 : 0), intersection.sector + 1,
1238  intersection.layer + 1, 1, 1);
1239  float zBarrelEfficiency = m_klmStripEfficiency->getBarrelEfficiency((intersection.isForward ? 1 : 0), intersection.sector + 1,
1240  intersection.layer + 1, 0, 1);
1241  klmMuidLikelihood->setExtBKLMEfficiencyValue(intersection.layer, phiBarrelEfficiency * zBarrelEfficiency);
1242 
1243  if (extState.lastBarrelExtLayer < intersection.layer) {
1244  extState.lastBarrelExtLayer = intersection.layer;
1245  }
1246  extState.hitLayerPattern |= (0x00000001 << intersection.layer);
1247  if (extState.lastBarrelHitLayer < intersection.layer) {
1248  extState.lastBarrelHitLayer = intersection.layer;
1249  }
1250  // If the updated point is outside the barrel, discard it and the Kalman-fitter adjustment
1251  r = intersection.position.perp();
1252  z = std::fabs(intersection.position.z() - m_OffsetZ);
1253  if ((r <= m_BarrelMinR) || (r >= m_BarrelMaxR) || (z >= m_BarrelHalfLength)) {
1254  intersection.chi2 = -1.0;
1255  }
1256  } else {
1257  // Record a no-hit track crossing if this step is strictly within a barrel sensitive volume
1258  vector<G4VPhysicalVolume*>::iterator j = find(m_BKLMVolumes->begin(), m_BKLMVolumes->end(), g4eState.GetG4Track()->GetVolume());
1259  if (j != m_BKLMVolumes->end()) {
1260  bool isDead = true; // by default, the nearest orthogonal strips are dead
1261  int section = intersection.isForward ?
1262  BKLMElementNumbers::c_ForwardSection :
1263  BKLMElementNumbers::c_BackwardSection;
1264  int sector = intersection.sector + 1; // from 0-based to 1-based enumeration
1265  int layer = intersection.layer + 1; // from 0-based to 1-based enumeration
1266  const bklm::Module* m = bklm::GeometryPar::instance()->findModule(section, sector, layer); // uses 1-based enumeration
1267  if (m) {
1268  const CLHEP::Hep3Vector localPosition = m->globalToLocal(intersection.position); // uses and returns position in cm
1269  int zStrip = m->getZStripNumber(localPosition);
1270  int phiStrip = m->getPhiStripNumber(localPosition);
1271  if (zStrip >= 0 && phiStrip >= 0) {
1272  uint16_t channel1, channel2;
1273  channel1 = m_klmElementNumbers->channelNumberBKLM(
1274  section, sector, layer,
1275  BKLMElementNumbers::c_ZPlane, zStrip);
1276  channel2 = m_klmElementNumbers->channelNumberBKLM(
1277  section, sector, layer,
1278  BKLMElementNumbers::c_PhiPlane, phiStrip);
1279  enum KLMChannelStatus::ChannelStatus status1, status2;
1280  status1 = m_klmChannelStatus->getChannelStatus(channel1);
1281  status2 = m_klmChannelStatus->getChannelStatus(channel2);
1282  if (status1 == KLMChannelStatus::c_Unknown ||
1283  status2 == KLMChannelStatus::c_Unknown)
1284  B2ERROR("No KLM channel status data."
1285  << LogVar("Section", section)
1286  << LogVar("Sector", sector)
1287  << LogVar("Layer", layer)
1288  << LogVar("Z strip", zStrip)
1289  << LogVar("Phi strip", phiStrip));
1290  isDead = (status1 == KLMChannelStatus::c_Dead ||
1291  status2 == KLMChannelStatus::c_Dead);
1292  }
1293  }
1294  if (!isDead) {
1295  extState.extLayerPattern |= (0x00000001 << intersection.layer); // valid extrapolation-crossing of the layer but no matching hit
1296  float phiBarrelEfficiency =
1297  m_klmStripEfficiency->getBarrelEfficiency(
1298  (intersection.isForward ? 1 : 0), intersection.sector + 1,
1299  intersection.layer + 1, BKLMElementNumbers::c_PhiPlane, 1);
1300  float zBarrelEfficiency =
1301  m_klmStripEfficiency->getBarrelEfficiency(
1302  (intersection.isForward ? 1 : 0), intersection.sector + 1,
1303  intersection.layer + 1, BKLMElementNumbers::c_ZPlane, 1);
1304  klmMuidLikelihood->setExtBKLMEfficiencyValue(intersection.layer, phiBarrelEfficiency * zBarrelEfficiency);
1305  } else {
1306  klmMuidLikelihood->setExtBKLMEfficiencyValue(intersection.layer, 0);
1307  }
1308  if (extState.lastBarrelExtLayer < intersection.layer) {
1309  extState.lastBarrelExtLayer = intersection.layer;
1310  }
1311  }
1312  }
1313  }
1314  }
1315 
1316  // Is the track in the endcap?
1317  if ((r > m_EndcapMinR) && (std::fabs(z - m_EndcapMiddleZ) < m_EndcapHalfLength)) {
1318  // Did the track cross the inner midplane of a detector module?
1319  if (findEndcapIntersection(extState, oldPosition, intersection)) {
1320  fromG4eToPhasespace(g4eState, intersection.covariance);
1321  if (findMatchingEndcapHit(intersection, extState.track)) {
1322  extState.extLayerPattern |= (0x00008000 << intersection.layer);
1323  //efficiency implementation
1324  float layerEndcapEfficiency = 1.;
1325  for (int plane = 1; plane <= EKLMElementNumbers::getMaximalPlaneNumber(); plane++) {
1326  layerEndcapEfficiency *= m_klmStripEfficiency->getEndcapEfficiency((intersection.isForward ? 1 : 0) + 1, intersection.sector + 1,
1327  intersection.layer + 1, plane, 1);
1328  }
1329  klmMuidLikelihood->setExtEKLMEfficiencyValue(intersection.layer, layerEndcapEfficiency);
1330 
1331  if (extState.lastEndcapExtLayer < intersection.layer) {
1332  extState.lastEndcapExtLayer = intersection.layer;
1333  }
1334  extState.hitLayerPattern |= (0x00008000 << intersection.layer);
1335  if (extState.lastEndcapHitLayer < intersection.layer) {
1336  extState.lastEndcapHitLayer = intersection.layer;
1337  }
1338  // If the updated point is outside the endcap, discard it and the Kalman-fitter adjustment
1339  r = intersection.position.perp();
1340  z = std::fabs(intersection.position.z() - m_OffsetZ);
1341  if ((r <= m_EndcapMinR) || (r >= m_EndcapMaxR) || (std::fabs(z - m_EndcapMiddleZ) >= m_EndcapHalfLength)) {
1342  intersection.chi2 = -1.0;
1343  }
1344  } else {
1345  bool isDead = true;
1346  int result, strip1, strip2;
1347  result = m_eklmTransformData->getStripsByIntersection(
1348  intersection.position, &strip1, &strip2);
1349  if (result == 0) {
1350  uint16_t channel1, channel2;
1351  channel1 = m_klmElementNumbers->channelNumberEKLM(strip1);
1352  channel2 = m_klmElementNumbers->channelNumberEKLM(strip2);
1353  enum KLMChannelStatus::ChannelStatus status1, status2;
1354  status1 = m_klmChannelStatus->getChannelStatus(channel1);
1355  status2 = m_klmChannelStatus->getChannelStatus(channel2);
1356  if (status1 == KLMChannelStatus::c_Unknown ||
1357  status2 == KLMChannelStatus::c_Unknown)
1358  B2ERROR("Incomplete KLM channel status data.");
1359  isDead = (status1 == KLMChannelStatus::c_Dead ||
1360  status2 == KLMChannelStatus::c_Dead);
1361  }
1362  if (!isDead) {
1363  extState.extLayerPattern |= (0x00008000 << intersection.layer); // valid extrapolation-crossing of the layer but no matching hit
1364  float layerEndcapEfficiency = 1.;
1365  for (int plane = 1; plane <= EKLMElementNumbers::getMaximalPlaneNumber(); plane++) {
1366  layerEndcapEfficiency *= m_klmStripEfficiency->getEndcapEfficiency((intersection.isForward ? 1 : 0) + 1, intersection.sector + 1,
1367  intersection.layer + 1, plane, 1);
1368  }
1369  klmMuidLikelihood->setExtEKLMEfficiencyValue(intersection.layer, layerEndcapEfficiency);
1370  } else {
1371  klmMuidLikelihood->setExtEKLMEfficiencyValue(intersection.layer, 0);
1372  }
1373  if (extState.lastEndcapExtLayer < intersection.layer) {
1374  extState.lastEndcapExtLayer = intersection.layer;
1375  }
1376  }
1377  }
1378  }
1379 
1380  // Create a new MuidHit and RelationEntry between it and the track.
1381  // Adjust geant4e's position, momentum and covariance based on matching hit and tell caller to update the geant4e state.
1382  if (intersection.chi2 >= 0.0) {
1383  TVector3 tpos(intersection.position.x(), intersection.position.y(), intersection.position.z());
1384  TVector3 tposAtHitPlane(intersection.positionAtHitPlane.x(),
1385  intersection.positionAtHitPlane.y(),
1386  intersection.positionAtHitPlane.z());
1387  KLMMuidHit* klmMuidHit = m_klmMuidHits.appendNew(extState.pdgCode, intersection.inBarrel, intersection.isForward,
1388  intersection.sector,
1389  intersection.layer, tpos,
1390  tposAtHitPlane, extState.tof, intersection.time, intersection.chi2);
1391  if (extState.track != nullptr) { extState.track->addRelationTo(klmMuidHit); }
1392  G4Point3D newPos(intersection.position.x() * CLHEP::cm,
1393  intersection.position.y() * CLHEP::cm,
1394  intersection.position.z() * CLHEP::cm);
1395  g4eState.SetPosition(newPos);
1396  G4Vector3D newMom(intersection.momentum.x() * CLHEP::GeV,
1397  intersection.momentum.y() * CLHEP::GeV,
1398  intersection.momentum.z() * CLHEP::GeV);
1399  g4eState.SetMomentum(newMom);
1400  G4ErrorTrajErr covG4e;
1401  fromPhasespaceToG4e(intersection.momentum, intersection.covariance, covG4e);
1402  g4eState.SetError(covG4e);
1403  extState.chi2 += intersection.chi2;
1404  extState.nPoint += 2; // two (orthogonal) independent hits per detector layer
1405  return true;
1406  }
1407 
1408  // Tell caller that the geant4e state was not modified.
1409  return false;
1410 
1411 }
1412 
1413 bool TrackExtrapolateG4e::findBarrelIntersection(ExtState& extState, const G4ThreeVector& oldPosition, Intersection& intersection)
1414 {
1415 
1416  // Be generous: allow outward-moving intersection to be in the dead space between
1417  // largest sensitive-volume Z and m_BarrelHalfLength, not necessarily in a geant4 sensitive volume
1418 
1419  if (std::fabs(intersection.position.z() - m_OffsetZ) > m_BarrelHalfLength) return false;
1420 
1421  double phi = intersection.position.phi();
1422  if (phi < 0.0) { phi += TWOPI; }
1423  if (phi > TWOPI - PI_8) { phi -= TWOPI; }
1424  int sector = (int)((phi + PI_8) / M_PI_4);
1425  int section = intersection.position.z() > m_OffsetZ ?
1426  BKLMElementNumbers::c_ForwardSection :
1427  BKLMElementNumbers::c_BackwardSection;
1428 
1429  double oldR = oldPosition * m_BarrelSectorPerp[sector];
1430  double newR = intersection.position * m_BarrelSectorPerp[sector];
1431 
1432  for (int layer = extState.firstBarrelLayer; layer <= m_OutermostActiveBarrelLayer; ++layer) {
1433  if (newR < m_BarrelModuleMiddleRadius[section][sector][layer]) break;
1434  if (oldR <= m_BarrelModuleMiddleRadius[section][sector][layer]) {
1435  extState.firstBarrelLayer = layer + 1; // ratchet outward for next call's loop starting value
1436  if (extState.firstBarrelLayer > m_OutermostActiveBarrelLayer) extState.escaped = true;
1437  intersection.inBarrel = true;
1438  intersection.isForward = intersection.position.z() > m_OffsetZ;
1439  intersection.layer = layer;
1440  intersection.sector = sector;
1441  return true;
1442  }
1443  }
1444 
1445  return false;
1446 
1447 }
1448 
1449 bool TrackExtrapolateG4e::findEndcapIntersection(ExtState& extState, const G4ThreeVector& oldPosition, Intersection& intersection)
1450 {
1451 
1452  // Be generous: allow intersection to be in the dead space between m_EndcapMinR and innermost
1453  // sensitive-volume radius or between outermost sensitive-volume radius and m_EndcapMaxR,
1454  // not necessarily in a geant4 sensitive volume
1455 
1456  if (oldPosition.perp() > m_EndcapMaxR) return false;
1457  if (intersection.position.perp() < m_EndcapMinR) return false;
1458 
1459  double oldZ = std::fabs(oldPosition.z() - m_OffsetZ);
1460  double newZ = std::fabs(intersection.position.z() - m_OffsetZ);
1461  bool isForward = intersection.position.z() > m_OffsetZ;
1462  int outermostLayer = isForward ? m_OutermostActiveForwardEndcapLayer
1463  : m_OutermostActiveBackwardEndcapLayer;
1464 
1465  for (int layer = extState.firstEndcapLayer; layer <= outermostLayer; ++layer) {
1466  if (newZ < m_EndcapModuleMiddleZ[layer]) break;
1467  if (oldZ <= m_EndcapModuleMiddleZ[layer]) {
1468  extState.firstEndcapLayer = layer + 1; // ratchet outward for next call's loop starting value
1469  if (extState.firstEndcapLayer > outermostLayer) extState.escaped = true;
1470  intersection.inBarrel = false;
1471  intersection.isForward = isForward;
1472  intersection.layer = layer;
1473  intersection.sector = m_eklmTransformData->getSectorByPosition(
1474  isForward ? 2 : 1, intersection.position) - 1;
1475  return true;
1476  }
1477  }
1478 
1479  return false;
1480 
1481 }
1482 
1483 bool TrackExtrapolateG4e::findMatchingBarrelHit(Intersection& intersection, const Track* track)
1484 
1485 {
1486  G4ThreeVector extPos0(intersection.position);
1487  double diffBestMagSq = 1.0E60;
1488  int bestHit = -1;
1489  int matchingLayer = intersection.layer + 1;
1490  G4ThreeVector n(m_BarrelSectorPerp[intersection.sector]);
1491  for (int h = 0; h < m_bklmHit2ds.getEntries(); ++h) {
1492  BKLMHit2d* hit = m_bklmHit2ds[h];
1493  if (hit->getLayer() != matchingLayer) continue;
1494  if (hit->isOutOfTime()) continue;
1495  if (std::fabs(hit->getTime() - m_MeanDt) > m_MaxDt) continue;
1496  G4ThreeVector diff(hit->getGlobalPositionX() - intersection.position.x(),
1497  hit->getGlobalPositionY() - intersection.position.y(),
1498  hit->getGlobalPositionZ() - intersection.position.z());
1499  double dn = diff * n; // in cm
1500  if (std::fabs(dn) < 2.0) {
1501  // Hit and extrapolated point are in the same sector
1502  diff -= n * dn;
1503  if (diff.mag2() < diffBestMagSq) {
1504  diffBestMagSq = diff.mag2();
1505  bestHit = h;
1506  extPos0 = intersection.position;
1507  }
1508  } else {
1509  // Accept a nearby hit in adjacent sector
1510  if (std::fabs(dn) > 50.0) continue;
1511  int sector = hit->getSector() - 1;
1512  int dSector = abs(intersection.sector - sector);
1513  if ((dSector != +1) && (dSector != m_BarrelNSector - 1)) continue;
1514  // Use the normal vector of the adjacent (hit's) sector
1515  G4ThreeVector nHit(m_BarrelSectorPerp[sector]);
1516  int section = intersection.isForward ?
1517  BKLMElementNumbers::c_ForwardSection :
1518  BKLMElementNumbers::c_BackwardSection;
1519  double dn2 = intersection.position * nHit - m_BarrelModuleMiddleRadius[section][sector][intersection.layer];
1520  dn = diff * nHit + dn2;
1521  if (std::fabs(dn) > 1.0) continue;
1522  // Project extrapolated track to the hit's plane in the adjacent sector
1523  G4ThreeVector extDir(intersection.momentum.unit());
1524  double extDirA = extDir * nHit;
1525  if (std::fabs(extDirA) < 1.0E-6) continue;
1526  G4ThreeVector projection = extDir * (dn2 / extDirA);
1527  if (projection.mag() > 15.0) continue;
1528  diff += projection - nHit * dn;
1529  if (diff.mag2() < diffBestMagSq) {
1530  diffBestMagSq = diff.mag2();
1531  bestHit = h;
1532  extPos0 = intersection.position - projection;
1533  }
1534  }
1535  }
1536 
1537  if (bestHit >= 0) {
1538  BKLMHit2d* hit = m_bklmHit2ds[bestHit];
1539  intersection.isForward = (hit->getSection() == 1);
1540  intersection.sector = hit->getSector() - 1;
1541  intersection.time = hit->getTime();
1542  double localVariance[2] = {m_BarrelScintVariance, m_BarrelScintVariance};
1543  if (hit->inRPC()) {
1544  int nStrips = hit->getPhiStripMax() - hit->getPhiStripMin() + 1;
1545  double dn = nStrips - 1.5;
1546  double factor = std::pow((0.9 + 0.4 * dn * dn), 1.5) * 0.60; // measured-in-Belle resolution
1547  localVariance[0] = m_BarrelPhiStripVariance[intersection.layer] * factor;
1548  nStrips = hit->getZStripMax() - hit->getZStripMin() + 1;
1549  dn = nStrips - 1.5;
1550  factor = std::pow((0.9 + 0.4 * dn * dn), 1.5) * 0.55; // measured-in-Belle resolution
1551  localVariance[1] = m_BarrelZStripVariance[intersection.layer] * factor;
1552  }
1553  G4ThreeVector hitPos(hit->getGlobalPositionX(), hit->getGlobalPositionY(), hit->getGlobalPositionZ());
1554  adjustIntersection(intersection, localVariance, hitPos, extPos0);
1555  if (intersection.chi2 >= 0.0) {
1556  intersection.hit = bestHit;
1557  hit->isOnTrack(true);
1558  if (track != nullptr) {
1559  track->addRelationTo(hit);
1560  RecoTrack* recoTrack = track->getRelatedTo<RecoTrack>();
1561  if (m_addHitsToRecoTrack) {
1562  recoTrack->addBKLMHit(hit, recoTrack->getNumberOfTotalHits() + 1);
1563  }
1564  }
1565  }
1566  }
1567  return intersection.chi2 >= 0.0;
1568 
1569 }
1570 
1571 bool TrackExtrapolateG4e::findMatchingEndcapHit(Intersection& intersection, const Track* track)
1572 {
1573  double diffBestMagSq = 1.0E60;
1574  int bestHit = -1;
1575  int matchingLayer = intersection.layer + 1;
1576  int matchingEndcap = (intersection.isForward ? 2 : 1);
1577  G4ThreeVector n(0.0, 0.0, (intersection.isForward ? 1.0 : -1.0));
1578  for (int h = 0; h < m_eklmHit2ds.getEntries(); ++h) {
1579  EKLMHit2d* hit = m_eklmHit2ds[h];
1580  if (hit->getLayer() != matchingLayer) continue;
1581  if (hit->getSection() != matchingEndcap) continue;
1582  // DIVOT no such function for EKLM!
1583  // if (hit->isOutOfTime()) continue;
1584  if (std::fabs(hit->getTime() - m_MeanDt) > m_MaxDt) continue;
1585  G4ThreeVector diff(hit->getPositionX() - intersection.position.x(),
1586  hit->getPositionY() - intersection.position.y(),
1587  hit->getPositionZ() - intersection.position.z());
1588  double dn = diff * n; // in cm
1589  if (std::fabs(dn) > 2.0) continue;
1590  diff -= n * dn;
1591  if (diff.mag2() < diffBestMagSq) {
1592  diffBestMagSq = diff.mag2();
1593  bestHit = h;
1594  }
1595  }
1596 
1597  if (bestHit >= 0) {
1598  EKLMHit2d* hit = m_eklmHit2ds[bestHit];
1599  intersection.hit = bestHit;
1600  intersection.isForward = (hit->getSection() == 2);
1601  intersection.sector = hit->getSector() - 1;
1602  intersection.time = hit->getTime();
1603  double localVariance[2] = {m_EndcapScintVariance, m_EndcapScintVariance};
1604  G4ThreeVector hitPos(hit->getPositionX(), hit->getPositionY(), hit->getPositionZ());
1605  adjustIntersection(intersection, localVariance, hitPos, intersection.position);
1606  if (intersection.chi2 >= 0.0) {
1607  // DIVOT no such function for EKLM!
1608  // hit->isOnTrack(true);
1609  if (track != nullptr) {
1610  RelationVector<EKLMAlignmentHit> eklmAlignmentHits = hit->getRelationsFrom<EKLMAlignmentHit>();
1611  track->addRelationTo(hit);
1612  RecoTrack* recoTrack = track->getRelatedTo<RecoTrack>();
1613  if (m_addHitsToRecoTrack) {
1614  for (unsigned int i = 0; i < eklmAlignmentHits.size(); ++i) {
1615  recoTrack->addEKLMHit(eklmAlignmentHits[i], recoTrack->getNumberOfTotalHits() + 1);
1616  }
1617  }
1618  }
1619  }
1620  }
1621  return intersection.chi2 >= 0.0;
1622 
1623 }
1624 
1625 void TrackExtrapolateG4e::adjustIntersection(Intersection& intersection, const double localVariance[2],
1626  const G4ThreeVector& hitPos, const G4ThreeVector& extPos0)
1627 {
1628 
1629 // Use the gain matrix formalism to get the corrected track parameters.
1630 // R. Fruhwirth, Application of Kalman Filtering, NIM A262 (1987) 444-450
1631 // Equations (7)
1632 // x_k^{k-1} = extPar[] 6 elements before filtering
1633 // C_k^{k-1} = extCov[] 6x6 symmetric before filtering
1634 // r_k^{k-1} = residual[] 2 elements before filtering
1635 // h_k = 2x6 projects cartesian coordinates to measurement-plane coordinates
1636 // H_k = @h_k/@x = jacobian[] 2x6 Jacobian of projection to measurement plane
1637 // R_k^{k-1} = correction[] 2x2 before Invert()
1638 // G_k = R^(-1) = correction[] 2x2 after Invert()
1639 // K_k = gain[] 6x2 Kalman gain matrix
1640 // x_k = extPar[] 6 elements after filtering
1641 // C_k = extCov[] 6x6 symmetric after filtering
1642 // r_k = residual[] 2 elements after filtering
1643 // Use the relation K*H*C = (C*H^T*R^-1)*H*C = C*(H^T*R^-1*H)*C^T
1644 
1645 // In most cases, extPos0 is the same as intersection.position. They differ only when
1646 // the nearest BKLM hit is in the sector adjacent to that of intersection.position.
1647 
1648  G4ThreeVector extPos(extPos0);
1649  G4ThreeVector extMom(intersection.momentum);
1650  G4ThreeVector extDir(extMom.unit());
1651  G4ThreeVector diffPos(hitPos - extPos);
1652  G4ErrorSymMatrix extCov(intersection.covariance);
1653 
1654 // Track parameters (x,y,z,px,py,pz) before correction
1655 
1656  G4ErrorMatrix extPar(6, 1); // initialized to all zeroes
1657  extPar[0][0] = extPos.x();
1658  extPar[1][0] = extPos.y();
1659  extPar[2][0] = extPos.z();
1660  extPar[3][0] = extMom.x();
1661  extPar[4][0] = extMom.y();
1662  extPar[5][0] = extMom.z();
1663 
1664  G4ThreeVector nA; // unit vector normal to the readout plane
1665  G4ThreeVector nB; // unit vector along phi- or x-readout direction (for barrel or endcap)
1666  G4ThreeVector nC; // unit vector along z- or y-readout direction (for barrel or endcap)
1667  if (intersection.inBarrel) {
1668  nA = m_BarrelSectorPerp[intersection.sector];
1669  nB = m_BarrelSectorPhi[intersection.sector];
1670  nC = G4ThreeVector(0.0, 0.0, 1.0);
1671  } else {
1672  double out = (intersection.isForward ? 1.0 : -1.0);
1673  nA = G4ThreeVector(0.0, 0.0, out);
1674  nB = G4ThreeVector(out, 0.0, 0.0);
1675  nC = G4ThreeVector(0.0, out, 0.0);
1676  }
1677 
1678 // Don't adjust the extrapolation if the track is nearly tangent to the readout plane.
1679 
1680  double extDirA = extDir * nA;
1681  if (std::fabs(extDirA) < 1.0E-6) return;
1682  double extDirBA = extDir * nB / extDirA;
1683  double extDirCA = extDir * nC / extDirA;
1684 
1685 // Move the extrapolated coordinate (at most a tiny amount!) to the plane of the hit.
1686 // If the moved point is outside the KLM, don't do Kalman filtering.
1687 
1688  G4ThreeVector move = extDir * ((diffPos * nA) / extDirA);
1689  extPos += move;
1690  diffPos -= move;
1691  intersection.positionAtHitPlane = extPos;
1692 
1693 // Projection jacobian onto the nB-nC measurement plane
1694 
1695  G4ErrorMatrix jacobian(2, 6); // initialized to all zeroes
1696  jacobian[0][0] = nB.x() - nA.x() * extDirBA;
1697  jacobian[0][1] = nB.y() - nA.y() * extDirBA;
1698  jacobian[0][2] = nB.z() - nA.z() * extDirBA;
1699  jacobian[1][0] = nC.x() - nA.x() * extDirCA;
1700  jacobian[1][1] = nC.y() - nA.y() * extDirCA;
1701  jacobian[1][2] = nC.z() - nA.z() * extDirCA;
1702 
1703 // Residuals of EXT track and KLM hit on the nB-nC measurement plane
1704 
1705  G4ErrorMatrix residual(2, 1); // initialized to all zeroes
1706  residual[0][0] = diffPos.x() * jacobian[0][0] + diffPos.y() * jacobian[0][1] + diffPos.z() * jacobian[0][2];
1707  residual[1][0] = diffPos.x() * jacobian[1][0] + diffPos.y() * jacobian[1][1] + diffPos.z() * jacobian[1][2];
1708 
1709 // Measurement errors in the detector plane
1710 
1711  G4ErrorSymMatrix hitCov(2, 0); // initialized to all zeroes
1712  hitCov[0][0] = localVariance[0];
1713  hitCov[1][1] = localVariance[1];
1714  // No magnetic field: increase the hit uncertainty
1715  if (m_MagneticField == 0.0) {
1716  hitCov[0][0] *= 10.0;
1717  hitCov[1][1] *= 10.0;
1718  }
1719 
1720 // Now get the correction matrix: combined covariance of EXT and KLM hit.
1721 // 1st dimension = nB, 2nd dimension = nC.
1722 
1723  G4ErrorSymMatrix correction(extCov.similarity(jacobian) + hitCov);
1724 
1725 // Ignore the best hit if it is too far from the extrapolated-track intersection in the hit's plane
1726 
1727  if (residual[0][0] * residual[0][0] > correction[0][0] * m_MaxDistSqInVariances) return;
1728  if (residual[1][0] * residual[1][0] > correction[1][1] * m_MaxDistSqInVariances) return;
1729 
1730  int fail = 0;
1731  correction.invert(fail);
1732  if (fail != 0) return;
1733 
1734 // Matrix inversion succeeeded and is reasonable.
1735 // Evaluate chi-squared increment assuming that the Kalman filter
1736 // won't be able to adjust the extrapolated track's position (fall-back).
1737 
1738  intersection.chi2 = (correction.similarityT(residual))[0][0];
1739 
1740 // Do the Kalman filtering
1741 
1742  G4ErrorMatrix gain((extCov * jacobian.T()) * correction);
1743  G4ErrorSymMatrix HRH(correction.similarityT(jacobian));
1744 
1745  extCov -= HRH.similarity(extCov);
1746  extPar += gain * residual;
1747  extPos.set(extPar[0][0], extPar[1][0], extPar[2][0]);
1748  extMom.set(extPar[3][0], extPar[4][0], extPar[5][0]);
1749 
1750 // Calculate the chi-squared increment using the Kalman-filtered state
1751 
1752  correction = hitCov - extCov.similarity(jacobian);
1753  correction.invert(fail);
1754  if (fail != 0) return;
1755 
1756  diffPos = hitPos - extPos;
1757  residual[0][0] = diffPos.x() * jacobian[0][0] + diffPos.y() * jacobian[0][1] + diffPos.z() * jacobian[0][2];
1758  residual[1][0] = diffPos.x() * jacobian[1][0] + diffPos.y() * jacobian[1][1] + diffPos.z() * jacobian[1][2];
1759  intersection.chi2 = (correction.similarityT(residual))[0][0];
1760 
1761 // Update the position, momentum and covariance of the point
1762 // Project the corrected extrapolation to the plane of the original
1763 // extrapolation's intersection.position. (Note: intersection.position is the same as
1764 // extPos0 in all cases except when nearest BKLM hit is in adjacent
1765 // sector, for which extPos0 is a projected position to the hit's plane.)
1766 // Also, leave the momentum magnitude unchanged.
1767 
1768  intersection.position = extPos + extDir * (((intersection.position - extPos) * nA) / extDirA);
1769  intersection.momentum = intersection.momentum.mag() * extMom.unit();
1770  intersection.covariance = extCov;
1771 
1772 }
1773 
1774 void TrackExtrapolateG4e::finishTrack(const ExtState& extState, KLMMuidLikelihood* klmMuidLikelihood, bool isForward)
1775 {
1776  /* Done with this track: compute KLM likelihoods and fill the relative dataobject. */
1777  int lastExtLayer = extState.lastBarrelExtLayer + extState.lastEndcapExtLayer + 1;
1778  unsigned int outcome = MuidElementNumbers::calculateExtrapolationOutcome(isForward, extState.escaped, extState.lastBarrelExtLayer,
1779  extState.lastEndcapExtLayer);
1780  klmMuidLikelihood->setOutcome(outcome);
1781  klmMuidLikelihood->setIsForward(isForward);
1782  klmMuidLikelihood->setBarrelExtLayer(extState.lastBarrelExtLayer);
1783  klmMuidLikelihood->setEndcapExtLayer(extState.lastEndcapExtLayer);
1784  klmMuidLikelihood->setBarrelHitLayer(extState.lastBarrelHitLayer);
1785  klmMuidLikelihood->setEndcapHitLayer(extState.lastEndcapHitLayer);
1786  klmMuidLikelihood->setExtLayer(lastExtLayer);
1787  klmMuidLikelihood->setHitLayer(((extState.lastEndcapHitLayer == -1) ?
1788  extState.lastBarrelHitLayer :
1789  extState.lastBarrelExtLayer + extState.lastEndcapHitLayer + 1));
1790  klmMuidLikelihood->setChiSquared(extState.chi2);
1791  klmMuidLikelihood->setDegreesOfFreedom(extState.nPoint);
1792  klmMuidLikelihood->setExtLayerPattern(extState.extLayerPattern);
1793  klmMuidLikelihood->setHitLayerPattern(extState.hitLayerPattern);
1794  /* Do KLM likelihood calculation. */
1795  if (outcome != MuidElementNumbers::c_NotReached) { /* Extrapolation reached KLM sensitive volume. */
1796  double denom = 0.0;
1797  int charge = klmMuidLikelihood->getCharge();
1798  std::vector<int> signedPdgVector = MuidElementNumbers::getPDGVector(charge);
1799  std::map<int, double> mapPdgPDF;
1800  for (int pdg : signedPdgVector) {
1801  auto search = m_MuidBuilderMap.find(pdg);
1802  if (search == m_MuidBuilderMap.end())
1803  B2FATAL("Something went wrong: PDF for PDG code " << pdg << " not found!");
1804  double pdf = (search->second)->getPDF(klmMuidLikelihood);
1805  denom += pdf;
1806  mapPdgPDF.insert(std::pair<int, double>(std::abs(pdg), pdf));
1807  }
1808  if (denom < 1.0E-20)
1809  klmMuidLikelihood->setJunkPDFValue(true); /* Anomaly: should be very rare. */
1810  else {
1811  for (auto const& [pdg, pdf] : mapPdgPDF) {
1812  klmMuidLikelihood->setPDFValue(pdf, std::abs(pdg));
1813  if (pdf > 0.0)
1814  klmMuidLikelihood->setLogL(std::log(pdf), std::abs(pdg));
1815  }
1816  }
1817  }
1818 }
Belle2::TrackExtrapolateG4e::m_BarrelScintVariance
double m_BarrelScintVariance
BKLM scintillator strip position variance (cm^2)
Definition: TrackExtrapolateG4e.h:389
Belle2::Intersection::layer
int layer
layer number (0..14 for barrel, 0..13 for endcap) of this point
Definition: TrackExtrapolateG4e.h:139
Belle2::RelationVector::size
size_t size() const
Get number of relations.
Definition: RelationVector.h:98
genfit::Exception
Exception class for error handling in GENFIT (provides storage for diagnostic information)
Definition: Exception.h:48
Belle2::ExtState::firstBarrelLayer
int firstBarrelLayer
MUID: outermost barrel layer encountered by the extrapolated track in the prior steps.
Definition: TrackExtrapolateG4e.h:109
Belle2::KLMMuidLikelihood::getCharge
int getCharge() const
Get the charge of the particle hypothesis used during the extrapolation.
Definition: KLMMuidLikelihood.cc:49
Belle2::TrackExtrapolateG4e::m_BarrelSectorPhi
G4ThreeVector m_BarrelSectorPhi[BKLMElementNumbers::getMaximalSectorNumber()+1]
azimuthal unit vector of each barrel sector
Definition: TrackExtrapolateG4e.h:399
Belle2::VOLTYPE_EKLM
@ VOLTYPE_EKLM
EKLM.
Definition: TrackExtrapolateG4e.h:87
Belle2::TrackExtrapolateG4e::m_MeanDt
double m_MeanDt
Mean hit - trigger time (ns)
Definition: TrackExtrapolateG4e.h:314
Belle2::TrackExtrapolateG4e::getStartPoint
ExtState getStartPoint(const Track &, int, G4ErrorFreeTrajState &)
Get the start point for a new reconstructed track with specific PDG hypothesis.
Definition: TrackExtrapolateG4e.cc:952
Belle2::RecoTrack::registerRequiredRelations
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:42
Belle2::BKLMStatus::setMaximalStrip
static void setMaximalStrip(int &module, int strip)
Set maximal strip number.
Definition: BKLMStatus.h:95
Belle2::MuidBuilder
Build the Muid likelihoods starting from the hit pattern and the transverse scattering in the KLM.
Definition: MuidBuilder.h:39
Belle2::EKLMGeometry::getNLayers
int getNLayers() const
Get number of layers.
Definition: EKLMGeometry.h:1695
Belle2::TrackExtrapolateG4e::createExtHit
void createExtHit(ExtHitStatus, const ExtState &, const G4ErrorFreeTrajState &, const G4StepPoint *, const G4TouchableHandle &)
Create another EXT extrapolation hit for a track candidate.
Definition: TrackExtrapolateG4e.cc:1190
Belle2::EKLMGeometry::ElementPosition::getZ
double getZ() const
Get Z coordinate.
Definition: EKLMGeometry.h:207
Belle2::RecoTrack::getCardinalRepresentation
genfit::AbsTrackRep * getCardinalRepresentation() const
Get a pointer to the cardinal track representation. You are not allowed to modify or delete it!
Definition: RecoTrack.h:547
Belle2::EKLM::TransformData::getSectorByPosition
int getSectorByPosition(int section, const HepGeom::Point3D< double > &position) const
Get sector by position.
Definition: TransformData.cc:418
Belle2::TrackExtrapolateG4e::m_EndcapMinR
double m_EndcapMinR
minimum radius (cm) of the endcaps
Definition: TrackExtrapolateG4e.h:405
Belle2::ExtState::extLayerPattern
int extLayerPattern
MUID: accumulated bit pattern of layers crossed by the extrapolated track.
Definition: TrackExtrapolateG4e.h:121
Belle2::BFieldManager::getField
static void getField(const double *pos, double *field)
return the magnetic field at a given position.
Definition: BFieldManager.h:110
genfit::MeasuredStateOnPlane
#StateOnPlane with additional covariance matrix.
Definition: MeasuredStateOnPlane.h:39
Belle2::Simulation::ExtManager::InitTrackPropagation
void InitTrackPropagation(G4ErrorMode)
Initialize for propagation of a track and set state to G4ErrorState_Propagating.
Definition: ExtManager.cc:109
Belle2::ExtState::lastEndcapExtLayer
int lastEndcapExtLayer
MUID: outermost endcap layer crossed by the extrapolated track.
Definition: TrackExtrapolateG4e.h:117
Belle2::TrackExtrapolateG4e::m_TargetExt
Simulation::ExtCylSurfaceTarget * m_TargetExt
virtual "target" cylinder for EXT (boundary beyond which extrapolation ends)
Definition: TrackExtrapolateG4e.h:356
Belle2::KLMMuidLikelihood::setExtLayer
void setExtLayer(int layer)
Set the outermost EKLM layer crossed in the extrapolation.
Definition: KLMMuidLikelihood.h:519
Belle2::Intersection::chi2
double chi2
chi-squared value of transverse deviation between extrapolated and measured hit positions
Definition: TrackExtrapolateG4e.h:151
genfit::AbsTrackRep::getPDG
int getPDG() const
Get the pdg code.
Definition: AbsTrackRep.h:272
Belle2::Const::electron
static const ChargedStable electron
electron particle
Definition: Const.h:533
Belle2::KLMChannelStatus::ChannelStatus
ChannelStatus
Channel status.
Definition: KLMChannelStatus.h:44
Belle2::TrackExtrapolateG4e::m_ExtMgr
Simulation::ExtManager * m_ExtMgr
Pointer to the ExtManager singleton.
Definition: TrackExtrapolateG4e.h:338
Belle2::RecoTrack::getMeasuredStateOnPlaneFromFirstHit
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:580
Belle2::VOLTYPE_TOP3
@ VOLTYPE_TOP3
TOP glue.
Definition: TrackExtrapolateG4e.h:73
Belle2::KLMMuidLikelihood::setBarrelExtLayer
void setBarrelExtLayer(int layer)
Set the outermost BKLM layer crossed in the extrapolation.
Definition: KLMMuidLikelihood.h:501
Belle2::TrackExtrapolateG4e::m_EndcapMaxR
double m_EndcapMaxR
maximum radius (cm) of the endcaps
Definition: TrackExtrapolateG4e.h:402
Belle2::TrackExtrapolateG4e::m_HypothesesMuid
const std::vector< Const::ChargedStable > * m_HypothesesMuid
ChargedStable hypotheses for MUID.
Definition: TrackExtrapolateG4e.h:344
Belle2::EKLMGeometry::getLayerPosition
const ElementPosition * getLayerPosition() const
Get position data for layers.
Definition: EKLMGeometry.h:1801
Belle2::VOLTYPE_BKLM1
@ VOLTYPE_BKLM1
BKLM RPC.
Definition: TrackExtrapolateG4e.h:83
Belle2::BKLMElementNumbers::c_ZPlane
@ c_ZPlane
Z.
Definition: BKLMElementNumbers.h:81
Belle2::KLMMuidLikelihood::setHitLayer
void setHitLayer(int layer)
Set the outermost KLM layer actually crossed by the track.
Definition: KLMMuidLikelihood.h:546
Belle2::TrackExtrapolateG4e::m_MaxDt
double m_MaxDt
Coincidence window half-width for in-time KLM hits (ns)
Definition: TrackExtrapolateG4e.h:317
Belle2::ExtState::hitLayerPattern
int hitLayerPattern
MUID: accumulated bit pattern of layers with matching hits.
Definition: TrackExtrapolateG4e.h:123
Belle2::TrackExtrapolateG4e::m_EndcapScintVariance
double m_EndcapScintVariance
EKLM scintillator strip position variance (cm^2)
Definition: TrackExtrapolateG4e.h:420
Belle2::KLMMuidLikelihood::setLogL
void setLogL(double logL, int pdg)
Set the log-likelihood.
Definition: KLMMuidLikelihood.h:400
Belle2::TrackExtrapolateG4e::m_bklmHit2ds
StoreArray< BKLMHit2d > m_bklmHit2ds
BKLM 2d hits.
Definition: TrackExtrapolateG4e.h:453
Belle2::KLMMuidLikelihood::setHitLayerPattern
void setHitLayerPattern(unsigned int pattern)
Set the pattern of the layers actually crossed by the track.
Definition: KLMMuidLikelihood.h:564
Belle2::KLMMuidLikelihood::setEndcapExtLayer
void setEndcapExtLayer(int layer)
Set the outermost EKLM layer crossed in the extrapolation.
Definition: KLMMuidLikelihood.h:510
Belle2::RelationsInterface::addRelationTo
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).
Definition: RelationsObject.h:144
Belle2::TrackExtrapolateG4e::m_klmMuidHits
StoreArray< KLMMuidHit > m_klmMuidHits
KLM muid hits.
Definition: TrackExtrapolateG4e.h:462
genfit::AbsTrackRep
Abstract base class for a track representation.
Definition: AbsTrackRep.h:66
Belle2::TrackExtrapolateG4e::findEndcapIntersection
bool findEndcapIntersection(ExtState &, const G4ThreeVector &, Intersection &)
Find the intersection point of the track with the crossed EKLM plane.
Definition: TrackExtrapolateG4e.cc:1449
Belle2::BKLMElementNumbers::c_ForwardSection
@ c_ForwardSection
Forward.
Definition: BKLMElementNumbers.h:47
Belle2::TrackExtrapolateG4e::findBarrelIntersection
bool findBarrelIntersection(ExtState &, const G4ThreeVector &, Intersection &)
Find the intersection point of the track with the crossed BKLM plane.
Definition: TrackExtrapolateG4e.cc:1413
Belle2::MuidElementNumbers::calculateExtrapolationOutcome
static unsigned int calculateExtrapolationOutcome(bool isForward, bool escaped, int lastBarrelLayer, int lastEndcapLayer)
Calculate the track extrapolation outcome.
Definition: MuidElementNumbers.cc:79
Belle2::RelationsInterface::getRelatedTo
TO * getRelatedTo(const std::string &name="", const std::string &namedRelation="") const
Get the object to which this object has a relation.
Definition: RelationsObject.h:250
Belle2::KLMElementNumbers::Instance
static const KLMElementNumbers & Instance()
Instantiation.
Definition: KLMElementNumbers.cc:31
Belle2::Simulation::ExtManager::EventTermination
void EventTermination(G4ErrorMode)
Terminate an event and set state to G4ErrorState_Init.
Definition: ExtManager.cc:139
Belle2::KLMMuidLikelihood::setDegreesOfFreedom
void setDegreesOfFreedom(int dof)
Set the number of degrees of freedom (= 2 times the number of KLM hits) for the chi-square computatio...
Definition: KLMMuidLikelihood.h:472
Belle2::Intersection::sector
int sector
sector number (0..7 for barrel, 0..3 for endcap) of this point
Definition: TrackExtrapolateG4e.h:137
Belle2::Simulation::ExtManager
It is the main interface for the user to define the setup and start the propagation.
Definition: ExtManager.h:53
Belle2::VOLTYPE_ARICH1
@ VOLTYPE_ARICH1
ARICH aerogel.
Definition: TrackExtrapolateG4e.h:75
Belle2::ExtState::lastBarrelHitLayer
int lastBarrelHitLayer
MUID: outermost barrel layer with a matching hit.
Definition: TrackExtrapolateG4e.h:115
Belle2::TrackExtrapolateG4e::m_tracks
StoreArray< Track > m_tracks
Tracks.
Definition: TrackExtrapolateG4e.h:471
Belle2::TrackExtrapolateG4e::findMatchingBarrelHit
bool findMatchingBarrelHit(Intersection &, const Track *)
Find the matching BKLM 2D hit nearest the intersection point of the track with the crossed BKLM plane...
Definition: TrackExtrapolateG4e.cc:1483
Belle2::Simulation::ExtManager::PropagateOneStep
G4int PropagateOneStep(G4ErrorTrajState *currentTS, G4ErrorMode mode=G4ErrorMode_PropForwards)
Propagate a track by one step.
Definition: ExtManager.cc:125
Belle2::TrackExtrapolateG4e::m_MinRadiusSq
double m_MinRadiusSq
Minimum squared radius (cm) outside of which extrapolation will continue.
Definition: TrackExtrapolateG4e.h:362
Belle2::TrackExtrapolateG4e::m_OffsetZ
double m_OffsetZ
offset (cm) along z axis of KLM midpoint from IP
Definition: TrackExtrapolateG4e.h:365
Belle2::KLMMuidLikelihood::setOutcome
void setOutcome(unsigned int outcome)
Set the outcome of this extrapolation.
Definition: KLMMuidLikelihood.h:483
Belle2::TrackExtrapolateG4e::m_MinKE
double m_MinKE
Minimum kinetic energy in MeV for extrapolation to continue.
Definition: TrackExtrapolateG4e.h:335
Belle2::ExtState::pdgCode
int pdgCode
Particle hypothesis that is being extrapolated.
Definition: TrackExtrapolateG4e.h:95
Belle2::KLMElementNumbers::channelNumberBKLM
uint16_t channelNumberBKLM(int section, int sector, int layer, int plane, int strip) const
Get channel number for BKLM.
Definition: KLMElementNumbers.cc:50
Belle2::VOLTYPE_TOP1
@ VOLTYPE_TOP1
TOP container.
Definition: TrackExtrapolateG4e.h:69
Belle2::Const::EDetector
EDetector
Enum for identifying the detector components (detector and subdetector).
Definition: Const.h:44
Belle2::EKLMAlignmentHit
This dataobject is used only for EKLM alignment.
Definition: EKLMAlignmentHit.h:39
Belle2::TrackExtrapolateG4e::m_OutermostActiveBarrelLayer
int m_OutermostActiveBarrelLayer
outermost barrel layer that is active for muon identification (user-defined)
Definition: TrackExtrapolateG4e.h:380
Belle2::TrackExtrapolateG4e
geant4e-based track extrapolation.
Definition: TrackExtrapolateG4e.h:168
Belle2::TrackExtrapolateG4e::m_TargetMuid
Simulation::ExtCylSurfaceTarget * m_TargetMuid
virtual "target" cylinder for MUID (boundary beyond which extrapolation ends)
Definition: TrackExtrapolateG4e.h:359
Belle2::bklm::Module
Define the geometry of a BKLM module Each sector [octant] contains Modules.
Definition: Module.h:86
Belle2::EKLMGeometry::ElementPosition::getInnerR
double getInnerR() const
Get inner radius.
Definition: EKLMGeometry.h:122
Belle2::KLMMuidLikelihood::setIsForward
void setIsForward(bool isForward)
Set if this extrapolation is in forward or backward B/EKLM.
Definition: KLMMuidLikelihood.h:492
Belle2::BKLMElementNumbers::getMaximalSectorNumber
static constexpr int getMaximalSectorNumber()
Get maximal sector number (1-based).
Definition: BKLMElementNumbers.h:243
Belle2::bklm::GeometryPar
Provides BKLM geometry parameters for simulation, reconstruction etc (from Gearbox or DataBase)
Definition: GeometryPar.h:48
Belle2::KLMMuidLikelihood::setJunkPDFValue
void setJunkPDFValue(bool flag)
Set the junk flag (1 if junk, 0 if not).
Definition: KLMMuidLikelihood.h:390
Belle2::TrackExtrapolateG4e::~TrackExtrapolateG4e
~TrackExtrapolateG4e()
destructor
Definition: TrackExtrapolateG4e.cc:133
Belle2::ExtHitStatus
ExtHitStatus
Define state of extrapolation for each recorded hit.
Definition: ExtHit.h:35
Belle2::ExtState
Data structure to define extrapolation state.
Definition: TrackExtrapolateG4e.h:91
Belle2::TrackExtrapolateG4e::registerVolumes
void registerVolumes()
Register the list of geant4 physical volumes whose entry/exit points will be saved during extrapolati...
Definition: TrackExtrapolateG4e.cc:798
Belle2::TrackExtrapolateG4e::m_MuidInitialized
bool m_MuidInitialized
Flag to indicate that MUID initialize() has been called.
Definition: TrackExtrapolateG4e.h:311
Belle2::TrackExtrapolateG4e::beginRun
void beginRun(bool flag)
Perform beginning-of-run actions.
Definition: TrackExtrapolateG4e.cc:317
Belle2::KLMMuidLikelihood::setExtEKLMEfficiencyValue
void setExtEKLMEfficiencyValue(int layer, float efficiency)
Set the efficiency of a given EKLM layer.
Definition: KLMMuidLikelihood.h:584
Belle2::Simulation::ExtManager::RunTermination
void RunTermination()
Terminate a run and set state to G4ErrorState_Init.
Definition: ExtManager.cc:153
Belle2::KLMChannelStatus::c_Dead
@ c_Dead
Dead channel (no signal).
Definition: KLMChannelStatus.h:53
Belle2::ExtState::lastEndcapHitLayer
int lastEndcapHitLayer
MUID: outermost endcap layer with a matching hit.
Definition: TrackExtrapolateG4e.h:119
Belle2::Simulation::ExtCylSurfaceTarget::GetDistanceFromPoint
virtual G4double GetDistanceFromPoint(const G4ThreeVector &point, const G4ThreeVector &direc) const
Get the distance from a point to the cylinder along direc.
Definition: ExtCylSurfaceTarget.cc:38
Belle2::TrackExtrapolateG4e::terminate
void terminate(bool flag)
Terminates this singleton.
Definition: TrackExtrapolateG4e.cc:400
Belle2::KLMMuidLikelihood::setExtLayerPattern
void setExtLayerPattern(unsigned int pattern)
Set the pattern of the layers crossed in the extrapolation.
Definition: KLMMuidLikelihood.h:555
Belle2::TrackExtrapolateG4e::event
void event(bool flag)
Performs track extrapolation for all tracks in one event.
Definition: TrackExtrapolateG4e.cc:341
Belle2::EKLM::TransformDataGlobalAligned::Instance
static const TransformDataGlobalAligned & Instance()
Instantiation.
Definition: TransformDataGlobalAligned.cc:17
Belle2::TrackExtrapolateG4e::m_BarrelPhiStripVariance
double m_BarrelPhiStripVariance[BKLMElementNumbers::getMaximalLayerNumber()+1]
BKLM RPC phi-measuring strip position variance (cm^2) by layer.
Definition: TrackExtrapolateG4e.h:383
genfit::AbsTrackRep::getPDGCharge
double getPDGCharge() const
Get the charge of the particle of the pdg code.
Definition: AbsTrackRep.cc:93
Belle2::Intersection::time
double time
time (ns) of matching BKLMHit2d
Definition: TrackExtrapolateG4e.h:149
Belle2::TrackExtrapolateG4e::m_ExtInitialized
bool m_ExtInitialized
Flag to indicate that EXT initialize() has been called.
Definition: TrackExtrapolateG4e.h:308
Belle2::ECL::ECLGeometryPar::Instance
static ECLGeometryPar * Instance()
Static method to get a reference to the ECLGeometryPar instance.
Definition: ECLGeometryPar.cc:151
Belle2::EKLMGeometry::getLayerShiftZ
double getLayerShiftZ() const
Get Z distance between two layers.
Definition: EKLMGeometry.h:1809
Belle2::bklm::GeometryPar::getHalfLength
double getHalfLength(void) const
Get the half-length along z of the BKLM.
Definition: GeometryPar.h:198
Belle2::RecoTrack::addEKLMHit
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:297
Belle2::KLMMuidLikelihood::setBarrelHitLayer
void setBarrelHitLayer(int layer)
Set the outermost BKLM layer actually crossed by the track.
Definition: KLMMuidLikelihood.h:528
Belle2::bklm::GeometryPar::getNLayer
int getNLayer(void) const
Get the number of modules in one sector.
Definition: GeometryPar.h:210
Belle2::TrackExtrapolateG4e::m_MagneticField
double m_MagneticField
Magnetic field z component (gauss) at origin.
Definition: TrackExtrapolateG4e.h:320
Belle2::RecoTrack
This is the Reconstruction Event-Data Model Track.
Definition: RecoTrack.h:78
Belle2::TrackExtrapolateG4e::m_HypothesesExt
const std::vector< Const::ChargedStable > * m_HypothesesExt
ChargedStable hypotheses for EXT.
Definition: TrackExtrapolateG4e.h:341
Belle2::TrackExtrapolateG4e::m_BarrelModuleMiddleRadius
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
Definition: TrackExtrapolateG4e.h:393
Belle2::Unit::T
static const double T
[tesla]
Definition: Unit.h:130
Belle2::ExtHit
Store one Ext hit as a ROOT object.
Definition: ExtHit.h:40
Belle2::TrackExtrapolateG4e::m_klmClusters
StoreArray< KLMCluster > m_klmClusters
KLM clusters.
Definition: TrackExtrapolateG4e.h:459
Belle2::BKLMElementNumbers::c_BackwardSection
@ c_BackwardSection
Backward.
Definition: BKLMElementNumbers.h:44
Belle2::TrackExtrapolateG4e::m_klmElementNumbers
const KLMElementNumbers * m_klmElementNumbers
KLM element numbers.
Definition: TrackExtrapolateG4e.h:432
Belle2::TrackExtrapolateG4e::m_klmChannelStatus
DBObjPtr< KLMChannelStatus > m_klmChannelStatus
Conditions-database object for KLM channel status.
Definition: TrackExtrapolateG4e.h:438
Belle2::Intersection::momentum
G4ThreeVector momentum
extrapolated-track momentum (GeV/c) at this intersection
Definition: TrackExtrapolateG4e.h:143
Belle2::TrackExtrapolateG4e::m_trackClusterSeparations
StoreArray< TrackClusterSeparation > m_trackClusterSeparations
Track cluster sepration.
Definition: TrackExtrapolateG4e.h:474
Belle2::BKLMElementNumbers::getMaximalLayerNumber
static constexpr int getMaximalLayerNumber()
Get maximal layer number (1-based).
Definition: BKLMElementNumbers.h:251
Belle2::RelationVector
Class for type safe access to objects that are referred to in relations.
Definition: DataStore.h:38
Belle2::Simulation::ExtManager::GetManager
static ExtManager * GetManager()
Get pointer to the instance of this singleton class (create if needed)
Definition: ExtManager.cc:75
Belle2::B2Vector3D
B2Vector3< double > B2Vector3D
typedef for common usage with double
Definition: B2Vector3.h:507
Belle2::bklm::GeometryPar::instance
static GeometryPar * instance(void)
Static method to get a reference to the singleton GeometryPar instance.
Definition: GeometryPar.cc:30
Belle2::TrackExtrapolateG4e::m_MaxECLTrackClusterDistance
double m_MaxECLTrackClusterDistance
user-defined maximum distance (mm) between ECLCluster and associated track (for EID)
Definition: TrackExtrapolateG4e.h:329
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
genfit::AbsTrackRep::getPosMomCov
virtual void getPosMomCov(const MeasuredStateOnPlane &state, TVector3 &pos, TVector3 &mom, TMatrixDSym &cov) const =0
Translates MeasuredStateOnPlane into 3D position, momentum and 6x6 covariance.
Belle2::TrackExtrapolateG4e::m_addHitsToRecoTrack
bool m_addHitsToRecoTrack
Parameter to add the found hits also to the reco tracks or not. Is turned off by default.
Definition: TrackExtrapolateG4e.h:426
Belle2::RecoTrack::addBKLMHit
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:283
Belle2::MuidElementNumbers::getPDGVector
static std::vector< int > getPDGVector()
Get a vector with all the hypothesis PDG codes used for Muid.
Definition: MuidElementNumbers.cc:167
Belle2::TrackExtrapolateG4e::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).
Definition: TrackExtrapolateG4e.cc:421
Belle2::TrackExtrapolateG4e::createMuidHit
bool createMuidHit(ExtState &, G4ErrorFreeTrajState &, KLMMuidLikelihood *, std::vector< std::map< const Track *, double > > *)
Create another MUID extrapolation hit for a track candidate.
Definition: TrackExtrapolateG4e.cc:1214
Belle2::TrackExtrapolateG4e::m_klmStripEfficiency
DBObjPtr< KLMStripEfficiency > m_klmStripEfficiency
Conditions-database object for KLM strip efficiency.
Definition: TrackExtrapolateG4e.h:441
Belle2::TrackExtrapolateG4e::m_recoTracks
StoreArray< RecoTrack > m_recoTracks
Reco tracks.
Definition: TrackExtrapolateG4e.h:468
Belle2::GearDir
GearDir is the basic class used for accessing the parameter store.
Definition: GearDir.h:41
Belle2::TrackExtrapolateG4e::m_BarrelMaxR
double m_BarrelMaxR
maximum radius (cm) of the barrel
Definition: TrackExtrapolateG4e.h:371
Belle2::ECL::ECLGeometryPar::ECLVolumeToCellID
int ECLVolumeToCellID(const G4VTouchable *)
Get Cell Id (LEP: new way)
Definition: ECLGeometryPar.cc:461
LogVar
Class to store variables with their name which were sent to the logging service.
Definition: LogVariableStream.h:24
Belle2::TrackExtrapolateG4e::getVolumeID
void getVolumeID(const G4TouchableHandle &, Const::EDetector &, int &)
Get the physical volume information for a geant4 physical volume.
Definition: TrackExtrapolateG4e.cc:870
Belle2::EKLMGeometry::getStripGeometry
const StripGeometry * getStripGeometry() const
Get strip geometry data.
Definition: EKLMGeometry.h:1873
Belle2::TrackExtrapolateG4e::m_MinPt
double m_MinPt
Minimum transverse momentum in MeV/c for extrapolation to be started.
Definition: TrackExtrapolateG4e.h:332
Belle2::EKLMElementNumbers::Instance
static const EKLMElementNumbers & Instance()
Instantiation.
Definition: EKLMElementNumbers.cc:20
Belle2::ExtState::isCosmic
bool isCosmic
True for back-propagation of a cosmic ray.
Definition: TrackExtrapolateG4e.h:97
Belle2::TrackExtrapolateG4e::m_BarrelNSector
int m_BarrelNSector
Number of barrel sectors.
Definition: TrackExtrapolateG4e.h:368
Belle2::KLMMuidLikelihood::setExtBKLMEfficiencyValue
void setExtBKLMEfficiencyValue(int layer, float efficiency)
Set the efficiency of a given BKLM layer.
Definition: KLMMuidLikelihood.h:574
Belle2::TrackExtrapolateG4e::m_OutermostActiveBackwardEndcapLayer
int m_OutermostActiveBackwardEndcapLayer
outermost backward-endcap layer that is active for muon identification (user-defined)
Definition: TrackExtrapolateG4e.h:417
Belle2::EKLMGeometry::getSectionPosition
const ElementPosition * getSectionPosition() const
Get position data for sections.
Definition: EKLMGeometry.h:1793
Belle2::EKLMGeometry::getNDetectorLayers
int getNDetectorLayers(int section) const
Get number of detector layers.
Definition: EKLMGeometry.cc:295
Belle2::EKLM::GeometryData
EKLM geometry data.
Definition: GeometryData.h:40
Belle2::EKLMElementNumbers::stripNumber
int stripNumber(int section, int layer, int sector, int plane, int strip) const
Get strip number.
Definition: EKLMElementNumbers.cc:212
Belle2::VOLTYPE_TOP2
@ VOLTYPE_TOP2
TOP quartz.
Definition: TrackExtrapolateG4e.h:71
Belle2::ExtState::escaped
bool escaped
MUID: flag to indicate that the extrapolated track escaped from the KLM.
Definition: TrackExtrapolateG4e.h:125
Belle2::TrackExtrapolateG4e::m_MuidBuilderMap
std::map< int, MuidBuilder * > m_MuidBuilderMap
PDF for the charged final state particle hypotheses.
Definition: TrackExtrapolateG4e.h:429
Belle2::TrackExtrapolateG4e::m_eklmTransformData
const EKLM::TransformDataGlobalAligned * m_eklmTransformData
EKLM transformation data.
Definition: TrackExtrapolateG4e.h:435
Belle2::TrackExtrapolateG4e::identifyMuon
void identifyMuon(int pdgCode, double tof, bool isCosmic, const G4ThreeVector &position, const G4ThreeVector &momentum, const G4ErrorSymMatrix &covariance)
Performs muon identification for a single track (specified in genfit2 units).
Definition: TrackExtrapolateG4e.cc:469
Belle2::KLMMuidLikelihood::setEndcapHitLayer
void setEndcapHitLayer(int layer)
Set the outermost EKLM layer actually crossed by the track.
Definition: KLMMuidLikelihood.h:537
Belle2::ExtState::firstEndcapLayer
int firstEndcapLayer
MUID: outermost barrel layer encountered by the extrapolated track in the prior steps.
Definition: TrackExtrapolateG4e.h:111
Belle2::Simulation::ExtManager::GetPropagator
G4ErrorPropagator * GetPropagator() const
Get the propagator.
Definition: ExtManager.h:79
Belle2::TrackExtrapolateG4e::m_EnterExit
std::map< G4VPhysicalVolume *, enum VolTypes > * m_EnterExit
Pointers to geant4 physical volumes whose entry/exit points will be saved.
Definition: TrackExtrapolateG4e.h:350
Belle2::EKLMElementNumbers::getMaximalPlaneNumber
static constexpr int getMaximalPlaneNumber()
Get maximal plane number.
Definition: EKLMElementNumbers.h:313
Belle2::KLMMuidLikelihood::setChiSquared
void setChiSquared(double chiSquared)
Set the chi-squared of the extrapolation.
Definition: KLMMuidLikelihood.h:463
Belle2::Intersection::position
G4ThreeVector position
extrapolated-track global position (cm) of this intersection
Definition: TrackExtrapolateG4e.h:141
Belle2::TrackExtrapolateG4e::m_EndcapModuleMiddleZ
double m_EndcapModuleMiddleZ[BKLMElementNumbers::getMaximalLayerNumber()+1]
hit-plane z (cm) of each IP layer relative to KLM midpoint
Definition: TrackExtrapolateG4e.h:423
Belle2::TrackExtrapolateG4e::swim
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.
Definition: TrackExtrapolateG4e.cc:517
Belle2::Intersection::covariance
G4ErrorSymMatrix covariance
extrapolated-track phase-space covariance matrix at this intersection
Definition: TrackExtrapolateG4e.h:145
Belle2::ExtState::length
double length
Length from start of extrapolation (rad lengths), updated during extrapolation.
Definition: TrackExtrapolateG4e.h:101
Belle2::EKLMGeometry::StripGeometry::getWidth
double getWidth() const
Get width.
Definition: EKLMGeometry.h:1199
Belle2::bklm::GeometryPar::getGap1InnerRadius
double getGap1InnerRadius(void) const
Get the radius of the inner tangent circle of gap 0 (innermost)
Definition: GeometryPar.h:228
Belle2::TrackExtrapolateG4e::findMatchingEndcapHit
bool findMatchingEndcapHit(Intersection &, const Track *)
Find the matching EKLM 2D hit nearest the intersection point of the track with the crossed EKLM plane...
Definition: TrackExtrapolateG4e.cc:1571
Belle2::gearbox::Interface::getLength
double getLength(const std::string &path="") const noexcept(false)
Get the parameter path as a double converted to the standard length unit.
Definition: Interface.h:261
Belle2::TrackExtrapolateG4e::m_eclClusters
StoreArray< ECLCluster > m_eclClusters
ECL clusters.
Definition: TrackExtrapolateG4e.h:447
Belle2::TrackExtrapolateG4e::m_klmLikelihoodParameters
DBObjPtr< KLMLikelihoodParameters > m_klmLikelihoodParameters
Conditions-database object for KLM likelihood parameters.
Definition: TrackExtrapolateG4e.h:444
Belle2::Simulation::ExtCylSurfaceTarget
Defines a closed cylinder for the geant4e "target", the surface that encloses the volume within which...
Definition: ExtCylSurfaceTarget.h:38
Belle2::TrackExtrapolateG4e::fromG4eToPhasespace
void fromG4eToPhasespace(const G4ErrorFreeTrajState &, G4ErrorSymMatrix &)
Convert the geant4e 5x5 covariance to phasespace 6x6 covariance.
Definition: TrackExtrapolateG4e.cc:1039
Belle2::TrackExtrapolateG4e::m_eklmHit2ds
StoreArray< EKLMHit2d > m_eklmHit2ds
EKLM 2d hits.
Definition: TrackExtrapolateG4e.h:456
Belle2::KLMElementNumbers::channelNumberEKLM
uint16_t channelNumberEKLM(int section, int sector, int layer, int plane, int strip) const
Get channel number for EKLM.
Definition: KLMElementNumbers.cc:64
Belle2::EKLM::TransformData::getStripsByIntersection
int getStripsByIntersection(const HepGeom::Point3D< double > &intersection, int *strip1, int *strip2) const
Find strips by intersection.
Definition: TransformData.cc:442
Belle2::VOLTYPE_ECL
@ VOLTYPE_ECL
ECL.
Definition: TrackExtrapolateG4e.h:81
Belle2::VOLTYPE_CDC
@ VOLTYPE_CDC
CDC.
Definition: TrackExtrapolateG4e.h:67
Belle2::KLMChannelStatus::c_Unknown
@ c_Unknown
Unknown status (no data).
Definition: KLMChannelStatus.h:47
Belle2::ExtState::chi2
double chi2
MUID: accumulated chi-squared of all in-plane transverse deviations between extrapolation and matchin...
Definition: TrackExtrapolateG4e.h:105
Belle2::ExtState::tof
double tof
Time of flight from IP (ns), updated during extrapolation.
Definition: TrackExtrapolateG4e.h:99
Belle2::TrackExtrapolateG4e::m_MaxDistSqInVariances
double m_MaxDistSqInVariances
user-defined maximum squared-distance (#variances) for matching hit to extrapolation
Definition: TrackExtrapolateG4e.h:323
Belle2::RecoTrack::getNumberOfTotalHits
unsigned int getNumberOfTotalHits() const
Return the number of cdc + svd + pxd + bklm + eklm hits.
Definition: RecoTrack.h:432
Belle2::BKLMElementNumbers::getMaximalSectionNumber
static constexpr int getMaximalSectionNumber()
Get maximal section number (0-based).
Definition: BKLMElementNumbers.h:235
Belle2::EKLMGeometry::ElementPosition::getOuterR
double getOuterR() const
Get outer radius.
Definition: EKLMGeometry.h:139
Belle2::bklm::GeometryPar::getScintHalfWidth
double getScintHalfWidth(void) const
Get the height of the entire volume of a scintillator strip (including TiO2 coating)
Definition: GeometryPar.h:135
Belle2::ExtState::track
const Track * track
Pointer to the reconstructed track.
Definition: TrackExtrapolateG4e.h:93
Belle2::bklm::GeometryPar::getNSector
int getNSector(void) const
Get the number of sectors of the BKLM.
Definition: GeometryPar.h:192
Belle2::KLMMuidHit
Store one muon-identification hit in the KLM as a ROOT object.
Definition: KLMMuidHit.h:34
Belle2::VOLTYPE_ARICH2
@ VOLTYPE_ARICH2
ARICH Img plate.
Definition: TrackExtrapolateG4e.h:77
Belle2::bklm::GeometryPar::getOuterRadius
double getOuterRadius(void) const
Get the radius of the inscribed circle of the outer polygon.
Definition: GeometryPar.h:204
Belle2::Intersection
intersection of muid-extrapolated track with a KLM layer
Definition: TrackExtrapolateG4e.h:129
Belle2::Const::muon
static const ChargedStable muon
muon particle
Definition: Const.h:534
Belle2::ExtState::nPoint
int nPoint
MUID: accumulated number of points with matching 2D hits.
Definition: TrackExtrapolateG4e.h:107
Belle2::ExtState::lastBarrelExtLayer
int lastBarrelExtLayer
MUID: outermost barrel layer crossed by the extrapolated track.
Definition: TrackExtrapolateG4e.h:113
Belle2::ExtState::directionAtIP
G4ThreeVector directionAtIP
MUID: initial direction of track, used for KLID.
Definition: TrackExtrapolateG4e.h:103
Belle2::Track
Class that bundles various TrackFitResults.
Definition: Track.h:35
Belle2::Intersection::inBarrel
bool inBarrel
flag to indicate if this point is in the barrel (true) or endcap (false)
Definition: TrackExtrapolateG4e.h:133
Belle2::VOLTYPE_BKLM2
@ VOLTYPE_BKLM2
BKLM scintillator.
Definition: TrackExtrapolateG4e.h:85
Belle2::KLMMuidLikelihood
Class to store the likelihoods from KLM with additional informations related to the extrapolation.
Definition: KLMMuidLikelihood.h:38
Belle2::TrackExtrapolateG4e::m_OutermostActiveForwardEndcapLayer
int m_OutermostActiveForwardEndcapLayer
outermost forward-endcap layer that is active for muon identification (user-defined)
Definition: TrackExtrapolateG4e.h:414
Belle2::TrackExtrapolateG4e::m_EndcapHalfLength
double m_EndcapHalfLength
half-length (cm) of either endcap
Definition: TrackExtrapolateG4e.h:411
Belle2::TrackExtrapolateG4e::m_BKLMVolumes
std::vector< G4VPhysicalVolume * > * m_BKLMVolumes
Pointers to BKLM geant4 sensitive (physical) volumes.
Definition: TrackExtrapolateG4e.h:353
Belle2::TrackExtrapolateG4e::m_klmMuidLikelihoods
StoreArray< KLMMuidLikelihood > m_klmMuidLikelihoods
KLM muid likelihoods.
Definition: TrackExtrapolateG4e.h:465
Belle2::EKLMGeometry::ElementPosition::getLength
double getLength() const
Get length.
Definition: EKLMGeometry.h:156
Belle2::TrackClusterSeparation
Store one Track-KLMCluster separation as a ROOT object.
Definition: TrackClusterSeparation.h:30
Belle2::Simulation::ExtManager::Initialize
void Initialize(const char[], const std::string &, double, double, bool, int, const std::vector< std::string > &)
Initialize Geant4 and Geant4e.
Definition: ExtManager.cc:179
Belle2::TrackExtrapolateG4e::m_extHits
StoreArray< ExtHit > m_extHits
Ext hits.
Definition: TrackExtrapolateG4e.h:450
Belle2::KLMMuidLikelihood::setPDFValue
void setPDFValue(double pdfValue, int pdg)
Set the normalized PDF.
Definition: KLMMuidLikelihood.h:327
Belle2::bklm::GeometryPar::findModule
const Module * findModule(int section, int sector, int layer) const
Get the pointer to the definition of a module.
Definition: GeometryPar.cc:716
Belle2::TrackExtrapolateG4e::endRun
void endRun(bool flag)
Perform end-of-run actions.
Definition: TrackExtrapolateG4e.cc:396
Belle2::BKLMElementNumbers::channelNumber
static uint16_t channelNumber(int section, int sector, int layer, int plane, int strip)
Get channel number.
Definition: BKLMElementNumbers.cc:27
Belle2::TrackExtrapolateG4e::m_BarrelZStripVariance
double m_BarrelZStripVariance[BKLMElementNumbers::getMaximalLayerNumber()+1]
BKLM RPC z-measuring strip position variance (cm^2) by layer.
Definition: TrackExtrapolateG4e.h:386
Belle2::TrackExtrapolateG4e::adjustIntersection
void adjustIntersection(Intersection &, const double *, const G4ThreeVector &, const G4ThreeVector &)
Nudge the track using the matching hit.
Definition: TrackExtrapolateG4e.cc:1625
Belle2::bklm::GeometryPar::getOffsetZ
double getOffsetZ(void) const
Get the global shift along a of the entire BKLM.
Definition: GeometryPar.h:174
Belle2::BKLMHit2d
Store one BKLM strip hit as a ROOT object.
Definition: BKLMHit2d.h:42
Belle2::TrackExtrapolateG4e::m_BarrelSectorPerp
G4ThreeVector m_BarrelSectorPerp[BKLMElementNumbers::getMaximalSectorNumber()+1]
normal unit vector of each barrel sector
Definition: TrackExtrapolateG4e.h:396
Belle2::TrackExtrapolateG4e::m_BarrelHalfLength
double m_BarrelHalfLength
half-length (cm) of the barrel
Definition: TrackExtrapolateG4e.h:377
Belle2::BKLMElementNumbers::c_PhiPlane
@ c_PhiPlane
Phi.
Definition: BKLMElementNumbers.h:84
Belle2::EKLM::GeometryData::Instance
static const GeometryData & Instance(enum DataSource dataSource=c_Database, const GearDir *gearDir=nullptr)
Instantiation.
Definition: GeometryData.cc:35
Belle2::Intersection::positionAtHitPlane
G4ThreeVector positionAtHitPlane
extrapolated-track position (cm) projected to the 2D hit's midplane
Definition: TrackExtrapolateG4e.h:147
Belle2::bklm::GeometryPar::getActiveMiddleRadius
double getActiveMiddleRadius(int section, int sector, int layer) const
Get the radial midpoint of the detector module's active volume of specified layer.
Definition: GeometryPar.cc:607
Belle2::TrackExtrapolateG4e::m_DefaultHypotheses
std::vector< Const::ChargedStable > * m_DefaultHypotheses
Default ChargedStable hypotheses (needed as call argument but not used)
Definition: TrackExtrapolateG4e.h:347
Belle2::EKLMHit2d
Class for 2d hits handling.
Definition: EKLMHit2d.h:39
Belle2::Intersection::isForward
bool isForward
flag to indicate if this point is in the forward (true) or backward (false) end
Definition: TrackExtrapolateG4e.h:135
Belle2::BKLMElementNumbers::moduleNumber
static uint16_t moduleNumber(int section, int sector, int layer, bool fatalError=true)
Get module number.
Definition: BKLMElementNumbers.cc:71
Belle2::VOLTYPE_ARICH3
@ VOLTYPE_ARICH3
ARICH HAPD window.
Definition: TrackExtrapolateG4e.h:79
Belle2::TrackExtrapolateG4e::m_EndcapMiddleZ
double m_EndcapMiddleZ
midpoint along z (cm) of the forward endcap from the KLM midpoint
Definition: TrackExtrapolateG4e.h:408
Belle2::TrackExtrapolateG4e::m_BarrelMinR
double m_BarrelMinR
minimum radius (cm) of the barrel
Definition: TrackExtrapolateG4e.h:374
Belle2::RecoTrack::getMeasuredStateOnPlaneFromLastHit
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:594
Belle2::Intersection::hit
int hit
index in {B,E}KLMHit2ds of matching hit
Definition: TrackExtrapolateG4e.h:131
Belle2::TrackExtrapolateG4e::initialize
void initialize(double minPt, double minKE, std::vector< Const::ChargedStable > &hypotheses)
Initialize for track extrapolation by the EXT module.
Definition: TrackExtrapolateG4e.cc:138
Belle2::TrackExtrapolateG4e::fromPhasespaceToG4e
void fromPhasespaceToG4e(const G4ThreeVector &, const G4ErrorSymMatrix &, G4ErrorTrajErr &)
Convert the phasespace covariance to geant4e covariance.
Definition: TrackExtrapolateG4e.cc:1141
Belle2::TrackExtrapolateG4e::finishTrack
void finishTrack(const ExtState &, KLMMuidLikelihood *, bool)
Complete muon identification after end of track extrapolation.
Definition: TrackExtrapolateG4e.cc:1774
Belle2::KLMMuidLikelihood::setPDGCode
void setPDGCode(int pdg)
Set the PDG code of the particle hypothesis used during the extrapolation.
Definition: KLMMuidLikelihood.h:317
Belle2::TrackExtrapolateG4e::m_MaxKLMTrackClusterDistance
double m_MaxKLMTrackClusterDistance
user-defined maximum distance (mm) between KLMCluster and associated track (for KLID)
Definition: TrackExtrapolateG4e.h:326