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