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