Belle II Software prerelease-10-00-00a
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
68using namespace Belle2;
69
70TrackExtrapolateG4e* TrackExtrapolateG4e::m_Singleton = nullptr;
71
78
79TrackExtrapolateG4e::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
136
137// Initialize for EXT
138void 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
191void 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
333void 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
357void 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
413void TrackExtrapolateG4e::endRun(bool)
414{
415}
416
417void 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
439void 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
488void 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(.0, extState.pdgCode, Const::EDetector::ECL, 0, EXT_FIRST,
508 extState.isCosmic,
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 // Add the following three quantities to the KLMCluster data members and store some relations
687 // To be safe, add them only if h is a valid object
688 if (h != nullptr) {
689 (*klmClusterInfo)[c].first->addRelationTo(h); // relation KLMCluster to TrackClusterSeparation
690 (*klmClusterInfo)[c].first->setClusterTrackSeparation(h->getDistance() / CLHEP::cm);
691 (*klmClusterInfo)[c].first->setClusterTrackSeparationAngle(h->getTrackClusterSeparationAngle());
692 (*klmClusterInfo)[c].first->setClusterTrackRotationAngle(h->getTrackRotationAngle());
693 if (extState.track != nullptr) {
694 extState.track->addRelationTo(h); // relation Track to TrackClusterSeparation
695 }
696 } else {
697 (*klmClusterInfo)[c].first->setClusterTrackSeparation(Const::doubleNaN);
698 (*klmClusterInfo)[c].first->setClusterTrackSeparationAngle(Const::doubleNaN);
699 (*klmClusterInfo)[c].first->setClusterTrackRotationAngle(Const::doubleNaN);
700 }
701 if (klmHit[c].getDistance() < minDistance) {
702 closestCluster = c;
703 minDistance = klmHit[c].getDistance();
704 }
705 }
706 if (minDistance < m_MaxKLMTrackClusterDistance) {
707 // set the relation Track to KLMCluster, using the distance as weight
708 // but for being consistent with the basf2 conventions, store the distance in cm
709 if (extState.track != nullptr) {
710 extState.track->addRelationTo((*klmClusterInfo)[closestCluster].first, 1. / (minDistance / CLHEP::cm));
711 }
712 }
713 }
714}
715
716// Swim one track for EXT until it stops or leaves the ECL-bounding cylinder
717void TrackExtrapolateG4e::swim(ExtState& extState, G4ErrorFreeTrajState& g4eState)
718{
719 if (extState.pdgCode == 0)
720 return;
721 if (g4eState.GetMomentum().perp() <= m_MinPt)
722 return;
723 if (m_TargetExt->GetDistanceFromPoint(g4eState.GetPosition()) < 0.0)
724 return;
725 G4ParticleDefinition* particle = G4ParticleTable::GetParticleTable()->FindParticle(extState.pdgCode);
726 double mass = particle->GetPDGMass();
727 double minPSq = (mass + m_MinKE) * (mass + m_MinKE) - mass * mass;
728 G4ErrorMode propagationMode = (extState.isCosmic ? G4ErrorMode_PropBackwards : G4ErrorMode_PropForwards);
729 m_ExtMgr->InitTrackPropagation(propagationMode);
730 while (true) {
731 const G4int errCode = m_ExtMgr->PropagateOneStep(&g4eState, propagationMode);
732 G4Track* track = g4eState.GetG4Track();
733 const G4Step* step = track->GetStep();
734 const G4StepPoint* preStepPoint = step->GetPreStepPoint();
735 const G4StepPoint* postStepPoint = step->GetPostStepPoint();
736 G4TouchableHandle preTouch = preStepPoint->GetTouchableHandle();
737 G4VPhysicalVolume* pVol = preTouch->GetVolume();
738 const G4int preStatus = preStepPoint->GetStepStatus();
739 const G4int postStatus = postStepPoint->GetStepStatus();
740 G4ThreeVector pos = track->GetPosition(); // this is at postStepPoint
741 G4ThreeVector mom = track->GetMomentum(); // ditto
742 // First step on this track?
743 if (extState.isCosmic)
744 mom = -mom;
745 if (preStatus == fUndefined) {
746 if (m_EnterExit->find(pVol) != m_EnterExit->end()) {
747 createExtHit(EXT_FIRST, extState, g4eState, preStepPoint, preTouch);
748 }
749 }
750 // Ignore the zero-length step by PropagateOneStep() at each boundary
751 if (step->GetStepLength() > 0.0) {
752 double dt = step->GetDeltaTime();
753 double dl = step->GetStepLength() / track->GetMaterial()->GetRadlen();
754 if (preStatus == fGeomBoundary) { // first step in this volume?
755 if (m_EnterExit->find(pVol) != m_EnterExit->end()) {
756 createExtHit(EXT_ENTER, extState, g4eState, preStepPoint, preTouch);
757 }
758 }
759 if (extState.isCosmic) {
760 extState.tof -= dt;
761 extState.length -= dl;
762 } else {
763 extState.tof += dt;
764 extState.length += dl;
765 }
766 // Last step in this volume?
767 if (postStatus == fGeomBoundary) {
768 if (m_EnterExit->find(pVol) != m_EnterExit->end()) {
769 createExtHit(EXT_EXIT, extState, g4eState, postStepPoint, preTouch);
770 }
771 }
772 }
773 // Post-step momentum too low?
774 if (errCode || (mom.mag2() < minPSq)) {
775 if (m_EnterExit->find(pVol) != m_EnterExit->end()) {
776 createExtHit(EXT_STOP, extState, g4eState, postStepPoint, preTouch);
777 }
778 break;
779 }
780 // Detect escapes from the imaginary target cylinder.
781 if (m_TargetExt->GetDistanceFromPoint(pos) < 0.0) {
782 if (m_EnterExit->find(pVol) != m_EnterExit->end()) {
783 createExtHit(EXT_ESCAPE, extState, g4eState, postStepPoint, preTouch);
784 }
785 break;
786 }
787 // Stop extrapolating as soon as the track curls inward too much
788 if (pos.perp2() < m_MinRadiusSq) {
789 break;
790 }
791 } // track-extrapolation "infinite" loop
792
793 m_ExtMgr->EventTermination(propagationMode);
794
795}
796
797// Register the list of volumes for which entry/exit point is to be saved during extrapolation
798void TrackExtrapolateG4e::registerVolumes()
799{
800 G4PhysicalVolumeStore* pvStore = G4PhysicalVolumeStore::GetInstance();
801 if (pvStore->size() == 0) {
802 B2FATAL("No geometry defined. Please create the geometry first.");
803 }
804 if (m_EnterExit != nullptr) // Only do this once
805 return;
806 m_EnterExit = new std::map<G4VPhysicalVolume*, enum VolTypes>;
807 m_BKLMVolumes = new std::vector<G4VPhysicalVolume*>;
808 for (std::vector<G4VPhysicalVolume*>::iterator iVol = pvStore->begin();
809 iVol != pvStore->end(); ++iVol) {
810 const G4String name = (*iVol)->GetName();
811 // Do not store ExtHits in CDC, so let's start from TOP and ARICH.
812 // TOP doesn't have one envelope; it has 16 "TOPModule"s
813 if (name.find("TOPModule") != std::string::npos) {
814 (*m_EnterExit)[*iVol] = VOLTYPE_TOP1;
815 }
816 // TOP quartz bar (=sensitive)
817 else if (name.find("_TOPPrism_") != std::string::npos ||
818 name.find("_TOPBarSegment") != std::string::npos ||
819 name.find("_TOPMirrorSegment") != std::string::npos) {
820 (*m_EnterExit)[*iVol] = VOLTYPE_TOP2;
821 }
822 // TOP quartz glue (not sensitive?)
823 else if (name.find("TOPBarSegment1Glue") != std::string::npos ||
824 name.find("TOPBarSegment2Glue") != std::string::npos ||
825 name.find("TOPMirrorSegmentGlue") != std::string::npos) {
826 (*m_EnterExit)[*iVol] = VOLTYPE_TOP3;
827 // ARICH volumes
828 } else if (name == "ARICH.AerogelSupportPlate") {
829 (*m_EnterExit)[*iVol] = VOLTYPE_ARICH1;
830 } else if (name == "ARICH.AerogelImgPlate") {
831 (*m_EnterExit)[*iVol] = VOLTYPE_ARICH2;
832 } else if (name.find("ARICH.HAPDWindow") != std::string::npos) {
833 (*m_EnterExit)[*iVol] = VOLTYPE_ARICH3;
834 }
835 // ECL crystal
836 else if (name.find("lv_barrel_crystal_") != std::string::npos ||
837 name.find("lv_forward_crystal_") != std::string::npos ||
838 name.find("lv_backward_crystal_") != std::string::npos) {
839 (*m_EnterExit)[*iVol] = VOLTYPE_ECL;
840 }
841 // Barrel KLM: BKLM.Layer**GasPhysical for RPCs or BKLM.Layer**ChimneyGasPhysical for RPCs
842 // BKLM.ScintActiveType*Physical for scintillator strips
843 else if (name.compare(0, 5, "BKLM.") == 0) {
844 if (name.find("GasPhysical") != std::string::npos) {
845 (*m_EnterExit)[*iVol] = VOLTYPE_BKLM1;
846 } else if (name.find("ScintActiveType") != std::string::npos) {
847 (*m_EnterExit)[*iVol] = VOLTYPE_BKLM2;
848 } else if ((name.find("ScintType") != std::string::npos) ||
849 (name.find("ElectrodePhysical") != std::string::npos)) {
850 m_BKLMVolumes->push_back(*iVol);
851 }
852 }
853 // Endcap KLM: StripSensitive_*
854 else if (name.compare(0, 14, "StripSensitive") == 0) {
855 (*m_EnterExit)[*iVol] = VOLTYPE_EKLM;
856 }
857 }
858
859}
860
861// Convert the physical volume to integer(-like) identifiers
862void TrackExtrapolateG4e::getVolumeID(const G4TouchableHandle& touch, Const::EDetector& detID, int& copyID)
863{
864
865 // default values
866 detID = Const::EDetector::invalidDetector;
867 copyID = 0;
868
869 G4VPhysicalVolume* pv = touch->GetVolume(0);
870 std::map<G4VPhysicalVolume*, enum VolTypes>::iterator it = m_EnterExit->find(pv);
871 if (it == m_EnterExit->end())
872 return;
873
874 switch (it->second) {
875 case VOLTYPE_CDC:
876 detID = Const::EDetector::CDC;
877 copyID = pv->GetCopyNo();
878 return;
879 case VOLTYPE_TOP1:
880 detID = Const::EDetector::TOP;
881 copyID = -(pv->GetCopyNo()); // negative to distinguish module and quartz hits
882 return;
883 case VOLTYPE_TOP2:
884 detID = Const::EDetector::TOP;
885 if (touch->GetHistoryDepth() >= 1)
886 copyID = touch->GetVolume(1)->GetCopyNo();
887 return;
888 case VOLTYPE_TOP3:
889 detID = Const::EDetector::TOP;
890 if (touch->GetHistoryDepth() >= 2)
891 copyID = touch->GetVolume(2)->GetCopyNo();
892 return;
893 case VOLTYPE_ARICH1:
894 detID = Const::EDetector::ARICH;
895 copyID = 12345;
896 return;
897 case VOLTYPE_ARICH2:
898 detID = Const::EDetector::ARICH;
899 copyID = 6789;
900 return;
901 case VOLTYPE_ARICH3:
902 detID = Const::EDetector::ARICH;
903 if (touch->GetHistoryDepth() >= 2)
904 copyID = touch->GetVolume(2)->GetCopyNo();
905 return;
906 case VOLTYPE_ECL:
907 detID = Const::EDetector::ECL;
908 copyID = ECL::ECLGeometryPar::Instance()->ECLVolumeToCellID(touch());
909 return;
910 case VOLTYPE_BKLM1: // BKLM RPCs
911 detID = Const::EDetector::BKLM;
912 if (touch->GetHistoryDepth() == DEPTH_RPC) {
913 // int plane = touch->GetCopyNumber(0);
914 int layer = touch->GetCopyNumber(4);
915 int sector = touch->GetCopyNumber(6);
916 int section = touch->GetCopyNumber(7);
917 copyID = BKLMElementNumbers::moduleNumber(section, sector, layer);
918 }
919 return;
920 case VOLTYPE_BKLM2: // BKLM scints
921 detID = Const::EDetector::BKLM;
922 if (touch->GetHistoryDepth() == DEPTH_SCINT) {
923 int strip = touch->GetCopyNumber(1);
924 int plane = (touch->GetCopyNumber(2) == BKLM_INNER) ?
925 BKLMElementNumbers::c_PhiPlane :
926 BKLMElementNumbers::c_ZPlane;
927 int layer = touch->GetCopyNumber(6);
928 int sector = touch->GetCopyNumber(8);
929 int section = touch->GetCopyNumber(9);
930 copyID = BKLMElementNumbers::channelNumber(
931 section, sector, layer, plane, strip);
932 BKLMStatus::setMaximalStrip(copyID, strip);
933 }
934 return;
935 case VOLTYPE_EKLM:
936 detID = Const::EDetector::EKLM;
937 copyID = EKLMElementNumbers::Instance().stripNumber(
938 touch->GetVolume(7)->GetCopyNo(),
939 touch->GetVolume(6)->GetCopyNo(),
940 touch->GetVolume(5)->GetCopyNo(),
941 touch->GetVolume(4)->GetCopyNo(),
942 touch->GetVolume(1)->GetCopyNo());
943 return;
944 }
945
946}
947
948ExtState TrackExtrapolateG4e::getStartPoint(const Track& b2track, int pdgCode, G4ErrorFreeTrajState& g4eState)
949{
950 ExtState extState = {&b2track, pdgCode, false, 0.0, 0.0, // for EXT and MUID
951 G4ThreeVector(0, 0, 1), 0.0, 0, 0, 0, -1, -1, -1, -1, 0, 0, false // for MUID only
952 };
953 RecoTrack* recoTrack = b2track.getRelatedTo<RecoTrack>();
954 if (recoTrack == nullptr) {
955 B2WARNING("Track without associated RecoTrack: skipping extrapolation for this track.");
956 extState.pdgCode = 0; // prevent start of extrapolation in swim()
957 return extState;
958 }
959 const genfit::AbsTrackRep* trackRep = recoTrack->getCardinalRepresentation();
960 // check for a valid track fit
961 if (!recoTrack->wasFitSuccessful(trackRep)) {
962 B2WARNING("RecoTrack fit failed for cardinal representation: skipping extrapolation for this track.");
963 extState.pdgCode = 0; // prevent start of extrapolation in swim()
964 return extState;
965 }
966 int charge = int(trackRep->getPDGCharge());
967 if (charge != 0) {
968 extState.pdgCode *= charge;
969 } else {
970 charge = 1; // should never happen but persist if it does
971 }
972 TVector3 firstPosition, firstMomentum, lastPosition, lastMomentum; // initialized to zeroes
973 TMatrixDSym firstCov(6), lastCov(6); // initialized to zeroes
974 try {
975 const genfit::MeasuredStateOnPlane& firstState = recoTrack->getMeasuredStateOnPlaneFromFirstHit(trackRep);
976 trackRep->getPosMomCov(firstState, firstPosition, firstMomentum, firstCov);
977 const genfit::MeasuredStateOnPlane& lastState = recoTrack->getMeasuredStateOnPlaneFromLastHit(trackRep);
978 trackRep->getPosMomCov(lastState, lastPosition, lastMomentum, lastCov);
979 // in genfit units (cm, GeV/c)
980 extState.tof = lastState.getTime(); // DIVOT: must be revised when IP profile (reconstructed beam spot) become available!
981 if (lastPosition.Mag2() < firstPosition.Mag2()) {
982 firstPosition = lastPosition;
983 firstMomentum = -lastMomentum;
984 firstCov = lastCov;
985 trackRep->getPosMomCov(firstState, lastPosition, lastMomentum, lastCov);
986 lastMomentum *= -1.0; // extrapolate backwards instead of forwards
987 extState.isCosmic = true;
988 extState.tof = firstState.getTime(); // DIVOT: must be revised when IP profile (reconstructed beam spot) become available!
989 }
990
991 G4ParticleDefinition* particle = G4ParticleTable::GetParticleTable()->FindParticle(extState.pdgCode);
992 if (extState.pdgCode != trackRep->getPDG()) {
993 double pSq = lastMomentum.Mag2();
994 double mass = particle->GetPDGMass() / CLHEP::GeV;
995 extState.tof *= std::sqrt((pSq + mass * mass) / (pSq + lastState.getMass() * lastState.getMass()));
996 }
997
998 extState.directionAtIP.set(firstMomentum.Unit().X(), firstMomentum.Unit().Y(), firstMomentum.Unit().Z());
999 if (m_MagneticField != 0.0) { // in gauss
1000 double radius = (firstMomentum.Perp() * CLHEP::GeV / CLHEP::eV) / (CLHEP::c_light * charge * m_MagneticField); // in cm
1001 double centerPhi = extState.directionAtIP.phi() - M_PI_2;
1002 double centerX = firstPosition.X() + radius * std::cos(centerPhi);
1003 double centerY = firstPosition.Y() + radius * std::sin(centerPhi);
1004 double pocaPhi = atan2(charge * centerY, charge * centerX) + M_PI;
1005 double ipPerp = extState.directionAtIP.perp();
1006 if (ipPerp > 0.0) {
1007 extState.directionAtIP.setX(+std::sin(pocaPhi) * ipPerp);
1008 extState.directionAtIP.setY(-std::cos(pocaPhi) * ipPerp);
1009 }
1010 } else {
1011 // No field: replace flaky covariance matrix with a diagonal one measured in 1.5T field
1012 // for a 10 GeV/c track ... and replace momentum magnitude with fixed 10 GeV/c
1013 lastCov *= 0.0;
1014 lastCov[0][0] = 5.0E-5;
1015 lastCov[1][1] = 1.0E-7;
1016 lastCov[2][2] = 5.0E-4;
1017 lastCov[3][3] = 3.5E-3;
1018 lastCov[4][4] = 3.5E-3;
1019 lastMomentum = lastMomentum.Unit() * 10.0;
1020 }
1021
1022 G4ThreeVector posG4e(lastPosition.X() * CLHEP::cm, lastPosition.Y() * CLHEP::cm,
1023 lastPosition.Z() * CLHEP::cm); // in Geant4 units (mm)
1024 G4ThreeVector momG4e(lastMomentum.X() * CLHEP::GeV, lastMomentum.Y() * CLHEP::GeV,
1025 lastMomentum.Z() * CLHEP::GeV); // in Geant4 units (MeV/c)
1026 G4ErrorSymMatrix covG4e(5, 0); // in Geant4e units (GeV/c, cm)
1027 fromPhasespaceToG4e(lastMomentum, lastCov, covG4e); // in Geant4e units (GeV/c, cm)
1028 g4eState.SetData("g4e_" + particle->GetParticleName(), posG4e, momG4e);
1029 g4eState.SetParameters(posG4e, momG4e); // compute private-state parameters from momG4e
1030 g4eState.SetError(covG4e);
1031 } catch (const genfit::Exception&) {
1032 B2WARNING("genfit::MeasuredStateOnPlane() exception: skipping extrapolation for this track. initial momentum = ("
1033 << firstMomentum.X() << "," << firstMomentum.Y() << "," << firstMomentum.Z() << ")");
1034 extState.pdgCode = 0; // prevent start of extrapolation in swim()
1035 }
1036
1037 return extState;
1038}
1039
1040
1041void TrackExtrapolateG4e::fromG4eToPhasespace(const G4ErrorFreeTrajState& g4eState, G4ErrorSymMatrix& covariance)
1042{
1043
1044 // Convert Geant4e covariance matrix with parameters 1/p, lambda, phi, yT, zT (in GeV/c, radians, cm)
1045 // to phase-space covariance matrix with parameters x, y, z, px, py, pz (in GeV/c, cm)
1046 // (1/p) = 1/sqrt( px^2 + py^2 + pz^2 )
1047 // phi = atan( py / px )
1048 // lambda = asin( pz / sqrt( px^2 + py^2 + pz^2 )
1049 // xT = x * cos(lambda) * cos(phi) + y * cos(lambda) * sin(phi) + z * sin(lambda)
1050 // yT = -x * sin(phi) + y * cos(phi)
1051 // zT = -x * sin(lambda) * cos(phi) - y * sin(lambda) * sin(phi) + z * cos(lambda)
1052
1053 G4ErrorFreeTrajParam param = g4eState.GetParameters();
1054 double p = 1.0 / (param.GetInvP() * CLHEP::GeV); // in GeV/c
1055 double pSq = p * p;
1056 double lambda = param.GetLambda(); // in radians
1057 double sinLambda = std::sin(lambda);
1058 double cosLambda = std::cos(lambda);
1059 double phi = param.GetPhi(); // in radians
1060 double sinPhi = std::sin(phi);
1061 double cosPhi = std::cos(phi);
1062
1063 // Transformation Jacobian 6x5 from Geant4e 5x5 to phase-space 6x6
1064
1065 G4ErrorMatrix jacobian(6, 5, 0); // All entries are initialized to 0
1066
1067 jacobian(4, 1) = -pSq * cosLambda * cosPhi; // @(px)/@(1/p)
1068 jacobian(5, 1) = -pSq * cosLambda * sinPhi; // @(py)/@(1/p)
1069 jacobian(6, 1) = -pSq * sinLambda; // @(pz)/@(1/p)
1070
1071 jacobian(4, 2) = -p * sinLambda * cosPhi; // @(px)/@(lambda)
1072 jacobian(5, 2) = -p * sinLambda * sinPhi; // @(py)/@(lambda)
1073 jacobian(6, 2) = p * cosLambda; // @(pz)/@(lambda)
1074
1075 jacobian(4, 3) = -p * cosLambda * sinPhi; // @(px)/@(phi)
1076 jacobian(5, 3) = p * cosLambda * cosPhi; // @(py)/@(phi)
1077
1078 jacobian(1, 4) = -sinPhi; // @(x)/@(yT)
1079 jacobian(2, 4) = cosPhi; // @(y)/@(yT)
1080
1081 jacobian(1, 5) = -sinLambda * cosPhi; // @(x)/@(zT)
1082 jacobian(2, 5) = -sinLambda * sinPhi; // @(y)/@(zT)
1083 jacobian(3, 5) = cosLambda; // @(z)/@(zT)
1084
1085 G4ErrorTrajErr g4eCov = g4eState.GetError();
1086 covariance.assign(g4eCov.similarity(jacobian));
1087
1088}
1089
1090void TrackExtrapolateG4e::fromPhasespaceToG4e(const TVector3& momentum, const TMatrixDSym& covariance, G4ErrorTrajErr& covG4e)
1091{
1092
1093 // Convert phase-space covariance matrix with parameters x, y, z, px, py, pz (in genfit2 units: cm, GeV/c)
1094 // to Geant4e covariance matrix with parameters 1/p, lambda, phi, yT, zT (in geant4e units: GeV/c, radians, cm)
1095 // xT = x * cos(lambda) * cos(phi) + y * cos(lambda) * sin(phi) + z * sin(lambda)
1096 // yT = -x * sin(phi) + y * cos(phi)
1097 // zT = -x * sin(lambda) * cos(phi) - y * sin(lambda) * sin(phi) + z * cos(lambda)
1098 // (1/p) = 1/sqrt( px^2 + py^2 + pz^2 )
1099 // phi = atan( py / px )
1100 // lambda = asin( pz / sqrt( px^2 + py^2 + pz^2 )
1101
1102 G4ErrorSymMatrix temp(6, 0);
1103 for (int k = 0; k < 6; ++k) {
1104 for (int j = k; j < 6; ++j) {
1105 temp[j][k] = covariance[j][k];
1106 }
1107 }
1108
1109 double pInvSq = 1.0 / momentum.Mag2();
1110 double pInv = std::sqrt(pInvSq);
1111 double pPerpInv = 1.0 / momentum.Perp();
1112 double sinLambda = momentum.CosTheta();
1113 double cosLambda = std::sqrt(1.0 - sinLambda * sinLambda);
1114 double phi = momentum.Phi();
1115 double cosPhi = std::cos(phi);
1116 double sinPhi = std::sin(phi);
1117
1118 // Transformation Jacobian 5x6 from phase-space 6x6 to Geant4e 5x5
1119 G4ErrorMatrix jacobian(5, 6, 0); // All entries are initialized to 0
1120
1121 jacobian(1, 4) = -pInvSq * cosLambda * cosPhi; // @(1/p)/@(px)
1122 jacobian(1, 5) = -pInvSq * cosLambda * sinPhi; // @(1/p)/@(py)
1123 jacobian(1, 6) = -pInvSq * sinLambda; // @(1/p)/@(pz)
1124
1125 jacobian(2, 4) = -pInv * sinLambda * cosPhi; // @(lambda)/@(px)
1126 jacobian(2, 5) = -pInv * sinLambda * sinPhi; // @(lambda)/@(py)
1127 jacobian(2, 6) = pInv * cosLambda; // @(lambda)/@(pz)
1128
1129 jacobian(3, 4) = -pPerpInv * sinPhi; // @(phi)/@(px)
1130 jacobian(3, 5) = pPerpInv * cosPhi; // @(phi)/@(py)
1131
1132 jacobian(4, 1) = -sinPhi; // @(yT)/@(x)
1133 jacobian(4, 2) = cosPhi; // @(yT)/@(y)
1134
1135 jacobian(5, 1) = -sinLambda * cosPhi; // @(zT)/@(x)
1136 jacobian(5, 2) = -sinLambda * sinPhi; // @(zT)/@(y)
1137 jacobian(5, 3) = cosLambda; // @(zT)/@(z)
1138
1139 covG4e = temp.similarity(jacobian);
1140
1141}
1142
1143void TrackExtrapolateG4e::fromPhasespaceToG4e(const G4ThreeVector& momentum, const G4ErrorSymMatrix& covariance,
1144 G4ErrorTrajErr& covG4e)
1145{
1146
1147 // Convert phase-space covariance matrix with parameters x, y, z, px, py, pz (in genfit2 units: cm, GeV/c)
1148 // to Geant4e covariance matrix with parameters 1/p, lambda, phi, yT, zT (in geant4e units: GeV/c, radians, cm)
1149 // xT = x * cos(lambda) * cos(phi) + y * cos(lambda) * sin(phi) + z * sin(lambda)
1150 // yT = -x * sin(phi) + y * cos(phi)
1151 // zT = -x * sin(lambda) * cos(phi) - y * sin(lambda) * sin(phi) + z * cos(lambda)
1152 // (1/p) = 1/sqrt( px^2 + py^2 + pz^2 )
1153 // phi = atan( py / px )
1154 // lambda = asin( pz / sqrt( px^2 + py^2 + pz^2 )
1155
1156 G4ErrorSymMatrix temp(covariance);
1157
1158 double pInvSq = 1.0 / momentum.mag2();
1159 double pInv = std::sqrt(pInvSq);
1160 double pPerpInv = 1.0 / momentum.perp();
1161 double sinLambda = momentum.cosTheta();
1162 double cosLambda = std::sqrt(1.0 - sinLambda * sinLambda);
1163 double phi = momentum.phi();
1164 double cosPhi = std::cos(phi);
1165 double sinPhi = std::sin(phi);
1166
1167 // Transformation Jacobian 5x6 from phase-space 6x6 to Geant4e 5x5
1168 G4ErrorMatrix jacobian(5, 6, 0);
1169
1170 jacobian(1, 4) = -pInvSq * cosLambda * cosPhi; // @(1/p)/@(px)
1171 jacobian(1, 5) = -pInvSq * cosLambda * sinPhi; // @(1/p)/@(py)
1172 jacobian(1, 6) = -pInvSq * sinLambda; // @(1/p)/@(pz)
1173
1174 jacobian(2, 4) = -pInv * sinLambda * cosPhi; // @(lambda)/@(px)
1175 jacobian(2, 5) = -pInv * sinLambda * sinPhi; // @(lambda)/@(py)
1176 jacobian(2, 6) = pInv * cosLambda; // @(lambda)/@(pz)
1177
1178 jacobian(3, 4) = -pPerpInv * sinPhi; // @(phi)/@(px)
1179 jacobian(3, 5) = pPerpInv * cosPhi; // @(phi)/@(py)
1180
1181 jacobian(4, 1) = -sinPhi; // @(yT)/@(x)
1182 jacobian(4, 2) = cosPhi; // @(yT)/@(y)
1183
1184 jacobian(5, 1) = -sinLambda * cosPhi; // @(zT)/@(x)
1185 jacobian(5, 2) = -sinLambda * sinPhi; // @(zT)/@(y)
1186 jacobian(5, 3) = cosLambda; // @(zT)/@(z)
1187 covG4e = temp.similarity(jacobian);
1188
1189}
1190
1191// write another volume-entry or volume-exit point on extrapolated track
1192void TrackExtrapolateG4e::createExtHit(ExtHitStatus status, const ExtState& extState,
1193 const G4ErrorFreeTrajState& g4eState,
1194 const G4StepPoint* stepPoint, const G4TouchableHandle& touch)
1195{
1196 Const::EDetector detID(Const::EDetector::invalidDetector);
1197 int copyID(0);
1198 getVolumeID(touch, detID, copyID);
1199 G4ThreeVector pos(stepPoint->GetPosition() / CLHEP::cm);
1200 G4ThreeVector mom(stepPoint->GetMomentum() / CLHEP::GeV);
1201 if (extState.isCosmic)
1202 mom = -mom;
1203 G4ErrorSymMatrix covariance(6, 0);
1204 fromG4eToPhasespace(g4eState, covariance);
1205 ExtHit* extHit = m_extHits.appendNew(extState.tof, extState.pdgCode, detID, copyID, status,
1206 extState.isCosmic,
1207 pos, mom, covariance);
1208 // If called standalone, there will be no associated track
1209 if (extState.track != nullptr)
1210 extState.track->addRelationTo(extHit);
1211}
1212
1213// Write another volume-entry point on track.
1214// The track state will be modified here by the Kalman fitter.
1215
1216bool TrackExtrapolateG4e::createMuidHit(ExtState& extState, G4ErrorFreeTrajState& g4eState, KLMMuidLikelihood* klmMuidLikelihood,
1217 std::vector<std::map<const Track*, double> >* bklmHitUsed)
1218{
1219
1220 Intersection intersection;
1221 intersection.hit = -1;
1222 intersection.chi2 = -1.0;
1223 intersection.position = g4eState.GetPosition() / CLHEP::cm;
1224 intersection.momentum = g4eState.GetMomentum() / CLHEP::GeV;
1225 G4ThreeVector prePos = g4eState.GetG4Track()->GetStep()->GetPreStepPoint()->GetPosition() / CLHEP::cm;
1226 G4ThreeVector oldPosition(prePos.x(), prePos.y(), prePos.z());
1227 double r = intersection.position.perp();
1228 double z = std::fabs(intersection.position.z() - m_OffsetZ);
1229
1230 // Is the track in the barrel?
1231 if ((r > m_BarrelMinR) && (r < m_BarrelMaxR) && (z < m_BarrelHalfLength)) {
1232 // Did the track cross the inner midplane of a detector module?
1233 if (findBarrelIntersection(extState, oldPosition, intersection)) {
1234 fromG4eToPhasespace(g4eState, intersection.covariance);
1235 if (findMatchingBarrelHit(intersection, extState.track)) {
1236 (*bklmHitUsed)[intersection.hit].insert(std::pair<const Track*, double>(extState.track, intersection.chi2));
1237 extState.extLayerPattern |= (0x00000001 << intersection.layer);
1238 float layerBarrelEfficiency = 1.;
1239 for (int plane = 0; plane <= BKLMElementNumbers::getMaximalPlaneNumber(); plane++) {
1240 layerBarrelEfficiency *= m_klmStripEfficiency->getBarrelEfficiency(
1241 intersection.isForward ? BKLMElementNumbers::c_ForwardSection : BKLMElementNumbers::c_BackwardSection,
1242 intersection.sector + 1, intersection.layer + 1, plane, 1);
1243 }
1244 klmMuidLikelihood->setExtBKLMEfficiencyValue(intersection.layer, layerBarrelEfficiency);
1245 if (extState.lastBarrelExtLayer < intersection.layer) {
1246 extState.lastBarrelExtLayer = intersection.layer;
1247 }
1248 extState.hitLayerPattern |= (0x00000001 << intersection.layer);
1249 if (extState.lastBarrelHitLayer < intersection.layer) {
1250 extState.lastBarrelHitLayer = intersection.layer;
1251 }
1252 // If the updated point is outside the barrel, discard it and the Kalman-fitter adjustment
1253 r = intersection.position.perp();
1254 z = std::fabs(intersection.position.z() - m_OffsetZ);
1255 if ((r <= m_BarrelMinR) || (r >= m_BarrelMaxR) || (z >= m_BarrelHalfLength)) {
1256 intersection.chi2 = -1.0;
1257 }
1258 } else {
1259 // Record a no-hit track crossing if this step is strictly within a barrel sensitive volume
1260 std::vector<G4VPhysicalVolume*>::iterator j = find(m_BKLMVolumes->begin(), m_BKLMVolumes->end(),
1261 g4eState.GetG4Track()->GetVolume());
1262 if (j != m_BKLMVolumes->end()) {
1263 bool isDead = true; // by default, the nearest orthogonal strips are dead
1264 int section = intersection.isForward ?
1265 BKLMElementNumbers::c_ForwardSection :
1266 BKLMElementNumbers::c_BackwardSection;
1267 int sector = intersection.sector + 1; // from 0-based to 1-based enumeration
1268 int layer = intersection.layer + 1; // from 0-based to 1-based enumeration
1269 const bklm::Module* m = bklm::GeometryPar::instance()->findModule(section, sector, layer); // uses 1-based enumeration
1270 if (m) {
1271 const CLHEP::Hep3Vector localPosition = m->globalToLocal(intersection.position); // uses and returns position in cm
1272 int zStrip = m->getZStripNumber(localPosition);
1273 int phiStrip = m->getPhiStripNumber(localPosition);
1274 if (zStrip >= 0 && phiStrip >= 0) {
1275 uint16_t channel1, channel2;
1276 channel1 = m_klmElementNumbers->channelNumberBKLM(
1277 section, sector, layer,
1278 BKLMElementNumbers::c_ZPlane, zStrip);
1279 channel2 = m_klmElementNumbers->channelNumberBKLM(
1280 section, sector, layer,
1281 BKLMElementNumbers::c_PhiPlane, phiStrip);
1282 enum KLMChannelStatus::ChannelStatus status1, status2;
1283 status1 = m_klmChannelStatus->getChannelStatus(channel1);
1284 status2 = m_klmChannelStatus->getChannelStatus(channel2);
1285 if (status1 == KLMChannelStatus::c_Unknown ||
1286 status2 == KLMChannelStatus::c_Unknown)
1287 B2ERROR("No KLM channel status data."
1288 << LogVar("Section", section)
1289 << LogVar("Sector", sector)
1290 << LogVar("Layer", layer)
1291 << LogVar("Z strip", zStrip)
1292 << LogVar("Phi strip", phiStrip));
1293 isDead = (status1 == KLMChannelStatus::c_Dead ||
1294 status2 == KLMChannelStatus::c_Dead);
1295 }
1296 }
1297 if (!isDead) {
1298 extState.extLayerPattern |= (0x00000001 << intersection.layer); // valid extrapolation-crossing of the layer but no matching hit
1299 float layerBarrelEfficiency = 1.;
1300 for (int plane = 0; plane <= BKLMElementNumbers::getMaximalPlaneNumber(); plane++) {
1301 layerBarrelEfficiency *= m_klmStripEfficiency->getBarrelEfficiency(
1302 intersection.isForward ? BKLMElementNumbers::c_ForwardSection : BKLMElementNumbers::c_BackwardSection,
1303 intersection.sector + 1, intersection.layer + 1, plane, 1);
1304 }
1305 klmMuidLikelihood->setExtBKLMEfficiencyValue(intersection.layer, layerBarrelEfficiency);
1306 } else {
1307 klmMuidLikelihood->setExtBKLMEfficiencyValue(intersection.layer, 0);
1308 }
1309 if (extState.lastBarrelExtLayer < intersection.layer) {
1310 extState.lastBarrelExtLayer = intersection.layer;
1311 }
1312 }
1313 }
1314 }
1315 }
1316
1317 // Is the track in the endcap?
1318 if ((r > m_EndcapMinR) && (std::fabs(z - m_EndcapMiddleZ) < m_EndcapHalfLength)) {
1319 // Did the track cross the inner midplane of a detector module?
1320 if (findEndcapIntersection(extState, oldPosition, intersection)) {
1321 fromG4eToPhasespace(g4eState, intersection.covariance);
1322 if (findMatchingEndcapHit(intersection, extState.track)) {
1323 extState.extLayerPattern |= (0x00008000 << intersection.layer);
1324 float layerEndcapEfficiency = 1.;
1325 for (int plane = 1; plane <= EKLMElementNumbers::getMaximalPlaneNumber(); plane++) {
1326 layerEndcapEfficiency *= m_klmStripEfficiency->getEndcapEfficiency(
1327 intersection.isForward ? EKLMElementNumbers::c_ForwardSection : EKLMElementNumbers::c_BackwardSection,
1328 intersection.sector + 1, intersection.layer + 1, plane, 1);
1329 }
1330 klmMuidLikelihood->setExtEKLMEfficiencyValue(intersection.layer, layerEndcapEfficiency);
1331 if (extState.lastEndcapExtLayer < intersection.layer) {
1332 extState.lastEndcapExtLayer = intersection.layer;
1333 }
1334 extState.hitLayerPattern |= (0x00008000 << intersection.layer);
1335 if (extState.lastEndcapHitLayer < intersection.layer) {
1336 extState.lastEndcapHitLayer = intersection.layer;
1337 }
1338 // If the updated point is outside the endcap, discard it and the Kalman-fitter adjustment
1339 r = intersection.position.perp();
1340 z = std::fabs(intersection.position.z() - m_OffsetZ);
1341 if ((r <= m_EndcapMinR) || (r >= m_EndcapMaxR) || (std::fabs(z - m_EndcapMiddleZ) >= m_EndcapHalfLength)) {
1342 intersection.chi2 = -1.0;
1343 }
1344 } else {
1345 bool isDead = true;
1346 int result, strip1, strip2;
1347 result = m_eklmTransformData->getStripsByIntersection(
1348 intersection.position, &strip1, &strip2);
1349 if (result == 0) {
1350 uint16_t channel1, channel2;
1351 channel1 = m_klmElementNumbers->channelNumberEKLM(strip1);
1352 channel2 = m_klmElementNumbers->channelNumberEKLM(strip2);
1353 enum KLMChannelStatus::ChannelStatus status1, status2;
1354 status1 = m_klmChannelStatus->getChannelStatus(channel1);
1355 status2 = m_klmChannelStatus->getChannelStatus(channel2);
1356 if (status1 == KLMChannelStatus::c_Unknown ||
1357 status2 == KLMChannelStatus::c_Unknown)
1358 B2ERROR("Incomplete KLM channel status data.");
1359 isDead = (status1 == KLMChannelStatus::c_Dead ||
1360 status2 == KLMChannelStatus::c_Dead);
1361 }
1362 if (!isDead) {
1363 extState.extLayerPattern |= (0x00008000 << intersection.layer); // valid extrapolation-crossing of the layer but no matching hit
1364 float layerEndcapEfficiency = 1.;
1365 for (int plane = 1; plane <= EKLMElementNumbers::getMaximalPlaneNumber(); plane++) {
1366 layerEndcapEfficiency *= m_klmStripEfficiency->getEndcapEfficiency(
1367 intersection.isForward ? EKLMElementNumbers::c_ForwardSection : EKLMElementNumbers::c_BackwardSection,
1368 intersection.sector + 1, intersection.layer + 1, plane, 1);
1369 }
1370 klmMuidLikelihood->setExtEKLMEfficiencyValue(intersection.layer, layerEndcapEfficiency);
1371 } else {
1372 klmMuidLikelihood->setExtEKLMEfficiencyValue(intersection.layer, 0);
1373 }
1374 if (extState.lastEndcapExtLayer < intersection.layer) {
1375 extState.lastEndcapExtLayer = intersection.layer;
1376 }
1377 }
1378 }
1379 }
1380
1381 // Create a new MuidHit and RelationEntry between it and the track.
1382 // Adjust geant4e's position, momentum and covariance based on matching hit and tell caller to update the geant4e state.
1383 if (intersection.chi2 >= 0.0) {
1384 ROOT::Math::XYZVector tpos(intersection.position.x(),
1385 intersection.position.y(),
1386 intersection.position.z());
1387 ROOT::Math::XYZVector tposAtHitPlane(intersection.positionAtHitPlane.x(),
1388 intersection.positionAtHitPlane.y(),
1389 intersection.positionAtHitPlane.z());
1390 KLMMuidHit* klmMuidHit = m_klmMuidHits.appendNew(extState.pdgCode,
1391 intersection.inBarrel, intersection.isForward,
1392 intersection.sector, intersection.layer,
1393 tpos, tposAtHitPlane,
1394 extState.tof, intersection.time,
1395 intersection.chi2);
1396 if (extState.track != nullptr) { extState.track->addRelationTo(klmMuidHit); }
1397 G4Point3D newPos(intersection.position.x() * CLHEP::cm,
1398 intersection.position.y() * CLHEP::cm,
1399 intersection.position.z() * CLHEP::cm);
1400 g4eState.SetPosition(newPos);
1401 G4Vector3D newMom(intersection.momentum.x() * CLHEP::GeV,
1402 intersection.momentum.y() * CLHEP::GeV,
1403 intersection.momentum.z() * CLHEP::GeV);
1404 g4eState.SetMomentum(newMom);
1405 G4ErrorTrajErr covG4e;
1406 fromPhasespaceToG4e(intersection.momentum, intersection.covariance, covG4e);
1407 g4eState.SetError(covG4e);
1408 extState.chi2 += intersection.chi2;
1409 extState.nPoint += 2; // two (orthogonal) independent hits per detector layer
1410 return true;
1411 }
1412
1413 // Tell caller that the geant4e state was not modified.
1414 return false;
1415
1416}
1417
1418bool TrackExtrapolateG4e::findBarrelIntersection(ExtState& extState, const G4ThreeVector& oldPosition, Intersection& intersection)
1419{
1420 // Be generous: allow outward-moving intersection to be in the dead space between
1421 // largest sensitive-volume Z and m_BarrelHalfLength, not necessarily in a geant4 sensitive volume
1422 if (std::fabs(intersection.position.z() - m_OffsetZ) > m_BarrelHalfLength)
1423 return false;
1424 double phi = intersection.position.phi();
1425 if (phi < 0.0)
1426 phi += TWOPI;
1427 if (phi > TWOPI - PI_8)
1428 phi -= TWOPI;
1429 int sector = (int)((phi + PI_8) / M_PI_4);
1430 int section = intersection.position.z() > m_OffsetZ ?
1431 BKLMElementNumbers::c_ForwardSection :
1432 BKLMElementNumbers::c_BackwardSection;
1433 double oldR = oldPosition * m_BarrelSectorPerp[sector];
1434 double newR = intersection.position * m_BarrelSectorPerp[sector];
1435 for (int layer = extState.firstBarrelLayer; layer <= m_OutermostActiveBarrelLayer; ++layer) {
1436 if (newR < m_BarrelModuleMiddleRadius[section][sector][layer]) break;
1437 if (oldR <= m_BarrelModuleMiddleRadius[section][sector][layer]) {
1438 extState.firstBarrelLayer = layer + 1; // ratchet outward for next call's loop starting value
1439 if (extState.firstBarrelLayer > m_OutermostActiveBarrelLayer) extState.escaped = true;
1440 intersection.inBarrel = true;
1441 intersection.isForward = intersection.position.z() > m_OffsetZ;
1442 intersection.layer = layer;
1443 intersection.sector = sector;
1444 return true;
1445 }
1446 }
1447 return false;
1448}
1449
1450bool TrackExtrapolateG4e::findEndcapIntersection(ExtState& extState, const G4ThreeVector& oldPosition, Intersection& intersection)
1451{
1452 // Be generous: allow intersection to be in the dead space between m_EndcapMinR and innermost
1453 // sensitive-volume radius or between outermost sensitive-volume radius and m_EndcapMaxR,
1454 // not necessarily in a geant4 sensitive volume
1455 if (oldPosition.perp() > m_EndcapMaxR)
1456 return false;
1457 if (intersection.position.perp() < m_EndcapMinR)
1458 return false;
1459 double oldZ = std::fabs(oldPosition.z() - m_OffsetZ);
1460 double newZ = std::fabs(intersection.position.z() - m_OffsetZ);
1461 bool isForward = intersection.position.z() > m_OffsetZ;
1462 int outermostLayer = isForward ? m_OutermostActiveForwardEndcapLayer
1463 : m_OutermostActiveBackwardEndcapLayer;
1464 for (int layer = extState.firstEndcapLayer; layer <= outermostLayer; ++layer) {
1465 if (newZ < m_EndcapModuleMiddleZ[layer])
1466 break;
1467 if (oldZ <= m_EndcapModuleMiddleZ[layer]) {
1468 extState.firstEndcapLayer = layer + 1; // ratchet outward for next call's loop starting value
1469 if (extState.firstEndcapLayer > outermostLayer)
1470 extState.escaped = true;
1471 intersection.inBarrel = false;
1472 intersection.isForward = isForward;
1473 intersection.layer = layer;
1474 intersection.sector = m_eklmTransformData->getSectorByPosition(
1475 isForward ? 2 : 1, intersection.position) - 1;
1476 return true;
1477 }
1478 }
1479 return false;
1480}
1481
1482bool TrackExtrapolateG4e::findMatchingBarrelHit(Intersection& intersection, const Track* track)
1483
1484{
1485 G4ThreeVector extPos0(intersection.position);
1486 double diffBestMagSq = 1.0E60;
1487 int bestHit = -1;
1488 int matchingLayer = intersection.layer + 1;
1489 G4ThreeVector n(m_BarrelSectorPerp[intersection.sector]);
1490 for (int h = 0; h < m_klmHit2ds.getEntries(); ++h) {
1491 KLMHit2d* hit = m_klmHit2ds[h];
1492 if (hit->getSubdetector() != KLMElementNumbers::c_BKLM)
1493 continue;
1494 if (hit->getLayer() != matchingLayer)
1495 continue;
1496 if (hit->isOutOfTime())
1497 continue;
1498 if (std::fabs(hit->getTime() - m_MeanDt) > m_MaxDt)
1499 continue;
1500 G4ThreeVector diff(hit->getPositionX() - intersection.position.x(),
1501 hit->getPositionY() - intersection.position.y(),
1502 hit->getPositionZ() - intersection.position.z());
1503 double dn = diff * n; // in cm
1504 if (std::fabs(dn) < 2.0) {
1505 // Hit and extrapolated point are in the same sector
1506 diff -= n * dn;
1507 if (diff.mag2() < diffBestMagSq) {
1508 diffBestMagSq = diff.mag2();
1509 bestHit = h;
1510 extPos0 = intersection.position;
1511 }
1512 } else {
1513 // Accept a nearby hit in adjacent sector
1514 if (std::fabs(dn) > 50.0)
1515 continue;
1516 int sector = hit->getSector() - 1;
1517 int dSector = abs(intersection.sector - sector);
1518 if ((dSector != +1) && (dSector != m_BarrelNSector - 1))
1519 continue;
1520 // Use the normal vector of the adjacent (hit's) sector
1521 G4ThreeVector nHit(m_BarrelSectorPerp[sector]);
1522 int section = intersection.isForward ?
1523 BKLMElementNumbers::c_ForwardSection :
1524 BKLMElementNumbers::c_BackwardSection;
1525 double dn2 = intersection.position * nHit - m_BarrelModuleMiddleRadius[section][sector][intersection.layer];
1526 dn = diff * nHit + dn2;
1527 if (std::fabs(dn) > 1.0)
1528 continue;
1529 // Project extrapolated track to the hit's plane in the adjacent sector
1530 G4ThreeVector extDir(intersection.momentum.unit());
1531 double extDirA = extDir * nHit;
1532 if (std::fabs(extDirA) < 1.0E-6)
1533 continue;
1534 G4ThreeVector projection = extDir * (dn2 / extDirA);
1535 if (projection.mag() > 15.0)
1536 continue;
1537 diff += projection - nHit * dn;
1538 if (diff.mag2() < diffBestMagSq) {
1539 diffBestMagSq = diff.mag2();
1540 bestHit = h;
1541 extPos0 = intersection.position - projection;
1542 }
1543 }
1544 }
1545
1546 if (bestHit >= 0) {
1547 KLMHit2d* hit = m_klmHit2ds[bestHit];
1548 intersection.isForward = (hit->getSection() == 1);
1549 intersection.sector = hit->getSector() - 1;
1550 intersection.time = hit->getTime();
1551 double localVariance[2] = {m_BarrelScintVariance, m_BarrelScintVariance};
1552 if (hit->inRPC()) {
1553 int nStrips = hit->getPhiStripMax() - hit->getPhiStripMin() + 1;
1554 double dn = nStrips - 1.5;
1555 double factor = std::pow((0.9 + 0.4 * dn * dn), 1.5) * 0.60; // measured-in-Belle resolution
1556 localVariance[0] = m_BarrelPhiStripVariance[intersection.layer] * factor;
1557 nStrips = hit->getZStripMax() - hit->getZStripMin() + 1;
1558 dn = nStrips - 1.5;
1559 factor = std::pow((0.9 + 0.4 * dn * dn), 1.5) * 0.55; // measured-in-Belle resolution
1560 localVariance[1] = m_BarrelZStripVariance[intersection.layer] * factor;
1561 } else {
1562 int nStrips = hit->getPhiStripMax() - hit->getPhiStripMin() + 1;
1563 localVariance[0] *= (nStrips * nStrips); // variance inflated for multi-strip hit
1564 nStrips = hit->getZStripMax() - hit->getZStripMin() + 1;
1565 localVariance[1] *= (nStrips * nStrips); // variance inflated for multi-strip hit
1566 }
1567 G4ThreeVector hitPos(hit->getPositionX(), hit->getPositionY(),
1568 hit->getPositionZ());
1569 adjustIntersection(intersection, localVariance, hitPos, extPos0);
1570 if (intersection.chi2 >= 0.0) {
1571 intersection.hit = bestHit;
1572 hit->isOnTrack(true);
1573 if (track != nullptr) {
1574 track->addRelationTo(hit, intersection.chi2);
1575 RecoTrack* recoTrack = track->getRelatedTo<RecoTrack>();
1576 if (m_addHitsToRecoTrack) {
1577 recoTrack->addBKLMHit(hit, recoTrack->getNumberOfTotalHits() + 1);
1578 }
1579 }
1580 }
1581 }
1582 return intersection.chi2 >= 0.0;
1583
1584}
1585
1586bool TrackExtrapolateG4e::findMatchingEndcapHit(Intersection& intersection, const Track* track)
1587{
1588 double diffBestMagSq = 1.0E60;
1589 int bestHit = -1;
1590 int matchingLayer = intersection.layer + 1;
1591 int matchingEndcap = (intersection.isForward ? 2 : 1);
1592 G4ThreeVector n(0.0, 0.0, (intersection.isForward ? 1.0 : -1.0));
1593 for (int h = 0; h < m_klmHit2ds.getEntries(); ++h) {
1594 KLMHit2d* hit = m_klmHit2ds[h];
1595 if (hit->getSubdetector() != KLMElementNumbers::c_EKLM)
1596 continue;
1597 if (hit->getLayer() != matchingLayer)
1598 continue;
1599 if (hit->getSection() != matchingEndcap)
1600 continue;
1601 // DIVOT no such function for EKLM!
1602 // if (hit->isOutOfTime()) continue;
1603 if (std::fabs(hit->getTime() - m_MeanDt) > m_MaxDt)
1604 continue;
1605 G4ThreeVector diff(hit->getPositionX() - intersection.position.x(),
1606 hit->getPositionY() - intersection.position.y(),
1607 hit->getPositionZ() - intersection.position.z());
1608 double dn = diff * n; // in cm
1609 if (std::fabs(dn) > 2.0)
1610 continue;
1611 diff -= n * dn;
1612 if (diff.mag2() < diffBestMagSq) {
1613 diffBestMagSq = diff.mag2();
1614 bestHit = h;
1615 }
1616 }
1617
1618 if (bestHit >= 0) {
1619 KLMHit2d* hit = m_klmHit2ds[bestHit];
1620 intersection.hit = bestHit;
1621 intersection.isForward = (hit->getSection() == 2);
1622 intersection.sector = hit->getSector() - 1;
1623 intersection.time = hit->getTime();
1624 double localVariance[2] = {m_EndcapScintVariance, m_EndcapScintVariance};
1625 int nStrips = hit->getXStripMax() - hit->getXStripMin() + 1;
1626 localVariance[0] *= (nStrips * nStrips); // variance inflated for multi-strip hit
1627 nStrips = hit->getYStripMax() - hit->getYStripMin() + 1;
1628 localVariance[1] *= (nStrips * nStrips); // variance inflated for multi-strip hit
1629 G4ThreeVector hitPos(hit->getPositionX(), hit->getPositionY(), hit->getPositionZ());
1630 adjustIntersection(intersection, localVariance, hitPos, intersection.position);
1631 if (intersection.chi2 >= 0.0) {
1632 // DIVOT no such function for EKLM!
1633 // hit->isOnTrack(true);
1634 if (track != nullptr) {
1635 track->addRelationTo(hit, intersection.chi2);
1636 RecoTrack* recoTrack = track->getRelatedTo<RecoTrack>();
1637 if (m_addHitsToRecoTrack) {
1638 for (const EKLMAlignmentHit& alignmentHit : hit->getRelationsFrom<EKLMAlignmentHit>()) {
1639 recoTrack->addEKLMHit(&alignmentHit, recoTrack->getNumberOfTotalHits() + 1);
1640 }
1641 }
1642 }
1643 }
1644 }
1645 return intersection.chi2 >= 0.0;
1646
1647}
1648
1649void TrackExtrapolateG4e::adjustIntersection(Intersection& intersection, const double localVariance[2],
1650 const G4ThreeVector& hitPos, const G4ThreeVector& extPos0)
1651{
1652 // Use the gain matrix formalism to get the corrected track parameters.
1653 // R. Fruhwirth, Application of Kalman Filtering, NIM A262 (1987) 444-450
1654 // Equations (7)
1655 // x_k^{k-1} = extPar[] 6 elements before filtering
1656 // C_k^{k-1} = extCov[] 6x6 symmetric before filtering
1657 // r_k^{k-1} = residual[] 2 elements before filtering
1658 // h_k = 2x6 projects cartesian coordinates to measurement-plane coordinates
1659 // H_k = @h_k/@x = jacobian[] 2x6 Jacobian of projection to measurement plane
1660 // R_k^{k-1} = correction[] 2x2 before Invert()
1661 // G_k = R^(-1) = correction[] 2x2 after Invert()
1662 // K_k = gain[] 6x2 Kalman gain matrix
1663 // x_k = extPar[] 6 elements after filtering
1664 // C_k = extCov[] 6x6 symmetric after filtering
1665 // r_k = residual[] 2 elements after filtering
1666 // Use the relation K*H*C = (C*H^T*R^-1)*H*C = C*(H^T*R^-1*H)*C^T
1667
1668 // In most cases, extPos0 is the same as intersection.position. They differ only when
1669 // the nearest BKLM hit is in the sector adjacent to that of intersection.position.
1670 G4ThreeVector extPos(extPos0);
1671 G4ThreeVector extMom(intersection.momentum);
1672 G4ThreeVector extDir(extMom.unit());
1673 G4ThreeVector diffPos(hitPos - extPos);
1674 G4ErrorSymMatrix extCov(intersection.covariance);
1675 // Track parameters (x,y,z,px,py,pz) before correction
1676 G4ErrorMatrix extPar(6, 1); // initialized to all zeroes
1677 extPar[0][0] = extPos.x();
1678 extPar[1][0] = extPos.y();
1679 extPar[2][0] = extPos.z();
1680 extPar[3][0] = extMom.x();
1681 extPar[4][0] = extMom.y();
1682 extPar[5][0] = extMom.z();
1683 G4ThreeVector nA; // unit vector normal to the readout plane
1684 G4ThreeVector nB; // unit vector along phi- or x-readout direction (for barrel or endcap)
1685 G4ThreeVector nC; // unit vector along z- or y-readout direction (for barrel or endcap)
1686 if (intersection.inBarrel) {
1687 nA = m_BarrelSectorPerp[intersection.sector];
1688 nB = m_BarrelSectorPhi[intersection.sector];
1689 nC = G4ThreeVector(0.0, 0.0, 1.0);
1690 } else {
1691 double out = (intersection.isForward ? 1.0 : -1.0);
1692 nA = G4ThreeVector(0.0, 0.0, out);
1693 nB = G4ThreeVector(out, 0.0, 0.0);
1694 nC = G4ThreeVector(0.0, out, 0.0);
1695 }
1696 // Don't adjust the extrapolation if the track is nearly tangent to the readout plane.
1697 double extDirA = extDir * nA;
1698 if (std::fabs(extDirA) < 1.0E-6)
1699 return;
1700 double extDirBA = extDir * nB / extDirA;
1701 double extDirCA = extDir * nC / extDirA;
1702 // Move the extrapolated coordinate (at most a tiny amount!) to the plane of the hit.
1703 // If the moved point is outside the KLM, don't do Kalman filtering.
1704 G4ThreeVector move = extDir * ((diffPos * nA) / extDirA);
1705 extPos += move;
1706 diffPos -= move;
1707 intersection.positionAtHitPlane = extPos;
1708 // Projection jacobian onto the nB-nC measurement plane
1709 G4ErrorMatrix jacobian(2, 6); // initialized to all zeroes
1710 jacobian[0][0] = nB.x() - nA.x() * extDirBA;
1711 jacobian[0][1] = nB.y() - nA.y() * extDirBA;
1712 jacobian[0][2] = nB.z() - nA.z() * extDirBA;
1713 jacobian[1][0] = nC.x() - nA.x() * extDirCA;
1714 jacobian[1][1] = nC.y() - nA.y() * extDirCA;
1715 jacobian[1][2] = nC.z() - nA.z() * extDirCA;
1716 // Residuals of EXT track and KLM hit on the nB-nC measurement plane
1717 G4ErrorMatrix residual(2, 1); // initialized to all zeroes
1718 residual[0][0] = diffPos.x() * jacobian[0][0] + diffPos.y() * jacobian[0][1] + diffPos.z() * jacobian[0][2];
1719 residual[1][0] = diffPos.x() * jacobian[1][0] + diffPos.y() * jacobian[1][1] + diffPos.z() * jacobian[1][2];
1720 // Measurement errors in the detector plane
1721 G4ErrorSymMatrix hitCov(2, 0); // initialized to all zeroes
1722 hitCov[0][0] = localVariance[0];
1723 hitCov[1][1] = localVariance[1];
1724 // No magnetic field: increase the hit uncertainty
1725 if (m_MagneticField == 0.0) {
1726 hitCov[0][0] *= 10.0;
1727 hitCov[1][1] *= 10.0;
1728 }
1729 // Now get the correction matrix: combined covariance of EXT and KLM hit.
1730 // 1st dimension = nB, 2nd dimension = nC.
1731 G4ErrorSymMatrix correction(extCov.similarity(jacobian) + hitCov);
1732 // Ignore the best hit if it is too far from the extrapolated-track intersection in the hit's plane
1733 if (residual[0][0] * residual[0][0] > correction[0][0] * m_MaxDistSqInVariances)
1734 return;
1735 if (residual[1][0] * residual[1][0] > correction[1][1] * m_MaxDistSqInVariances)
1736 return;
1737 int fail = 0;
1738 correction.invert(fail);
1739 if (fail != 0)
1740 return;
1741 // Matrix inversion succeeded and is reasonable.
1742 // Evaluate chi-squared increment assuming that the Kalman filter
1743 // won't be able to adjust the extrapolated track's position (fall-back).
1744 intersection.chi2 = (correction.similarityT(residual))[0][0];
1745 // Do the Kalman filtering
1746 G4ErrorMatrix gain((extCov * jacobian.T()) * correction);
1747 G4ErrorSymMatrix HRH(correction.similarityT(jacobian));
1748 extCov -= HRH.similarity(extCov);
1749 extPar += gain * residual;
1750 extPos.set(extPar[0][0], extPar[1][0], extPar[2][0]);
1751 extMom.set(extPar[3][0], extPar[4][0], extPar[5][0]);
1752 // Calculate the chi-squared increment using the Kalman-filtered state
1753 correction = hitCov - extCov.similarity(jacobian);
1754 correction.invert(fail);
1755 if (fail != 0)
1756 return;
1757 diffPos = hitPos - extPos;
1758 residual[0][0] = diffPos.x() * jacobian[0][0] + diffPos.y() * jacobian[0][1] + diffPos.z() * jacobian[0][2];
1759 residual[1][0] = diffPos.x() * jacobian[1][0] + diffPos.y() * jacobian[1][1] + diffPos.z() * jacobian[1][2];
1760 intersection.chi2 = (correction.similarityT(residual))[0][0];
1761 // Update the position, momentum and covariance of the point
1762 // Project the corrected extrapolation to the plane of the original
1763 // extrapolation's intersection.position. (Note: intersection.position is the same as
1764 // extPos0 in all cases except when nearest BKLM hit is in adjacent
1765 // sector, for which extPos0 is a projected position to the hit's plane.)
1766 // Also, leave the momentum magnitude unchanged.
1767 intersection.position = extPos + extDir * (((intersection.position - extPos) * nA) / extDirA);
1768 intersection.momentum = intersection.momentum.mag() * extMom.unit();
1769 intersection.covariance = extCov;
1770}
1771
1772void TrackExtrapolateG4e::finishTrack(const ExtState& extState, KLMMuidLikelihood* klmMuidLikelihood, bool isForward)
1773{
1774 /* Done with this track: compute KLM likelihoods and fill the relative dataobject. */
1775 int lastExtLayer = extState.lastBarrelExtLayer + extState.lastEndcapExtLayer + 1;
1776 unsigned int outcome = MuidElementNumbers::calculateExtrapolationOutcome(
1777 isForward, extState.escaped, extState.lastBarrelExtLayer, extState.lastEndcapExtLayer);
1778 klmMuidLikelihood->setOutcome(outcome);
1779 klmMuidLikelihood->setIsForward(isForward);
1780 klmMuidLikelihood->setBarrelExtLayer(extState.lastBarrelExtLayer);
1781 klmMuidLikelihood->setEndcapExtLayer(extState.lastEndcapExtLayer);
1782 klmMuidLikelihood->setBarrelHitLayer(extState.lastBarrelHitLayer);
1783 klmMuidLikelihood->setEndcapHitLayer(extState.lastEndcapHitLayer);
1784 klmMuidLikelihood->setExtLayer(lastExtLayer);
1785 klmMuidLikelihood->setHitLayer(((extState.lastEndcapHitLayer == -1) ?
1786 extState.lastBarrelHitLayer :
1787 extState.lastBarrelExtLayer + extState.lastEndcapHitLayer + 1));
1788 klmMuidLikelihood->setChiSquared(extState.chi2);
1789 klmMuidLikelihood->setDegreesOfFreedom(extState.nPoint);
1790 klmMuidLikelihood->setExtLayerPattern(extState.extLayerPattern);
1791 klmMuidLikelihood->setHitLayerPattern(extState.hitLayerPattern);
1792 /* Do KLM likelihood calculation. */
1793 if (outcome != MuidElementNumbers::c_NotReached) { /* Extrapolation reached KLM sensitive volume. */
1794 double denom = 0.0;
1795 int charge = klmMuidLikelihood->getCharge();
1796 std::vector<int> signedPdgVector = MuidElementNumbers::getPDGVector(charge);
1797 std::map<int, double> mapPdgPDF;
1798 for (int pdg : signedPdgVector) {
1799 auto search = m_MuidBuilderMap.find(pdg);
1800 if (search == m_MuidBuilderMap.end())
1801 B2FATAL("Something went wrong: PDF for PDG code " << pdg << " not found!");
1802 double pdf = (search->second)->getPDF(klmMuidLikelihood);
1803 denom += pdf;
1804 mapPdgPDF.insert(std::pair<int, double>(std::abs(pdg), pdf));
1805 }
1806 if (denom < 1.0E-20)
1807 klmMuidLikelihood->setJunkPDFValue(true); /* Anomaly: should be very rare. */
1808 else {
1809 for (auto const& [pdg, pdf] : mapPdgPDF) {
1810 klmMuidLikelihood->setPDFValue(pdf, std::abs(pdg));
1811 if (pdf > 0.0)
1812 klmMuidLikelihood->setLogL(std::log(pdf), std::abs(pdg));
1813 }
1814 }
1815 }
1816}
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:660
Belle2::Const::EDetector
EDetector
Enum for identifying the detector components (detector and subdetector).
Definition Const.h:42
Belle2::Const::doubleNaN
static const double doubleNaN
quiet_NaN
Definition Const.h:703
Belle2::Const::electron
static const ChargedStable electron
electron particle
Definition Const.h:659
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::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::getLayerPosition
const ElementPosition * getLayerPosition() const
Get position data for layers.
Definition EKLMGeometry.h:1799
Belle2::EKLMGeometry::getLayerShiftZ
double getLayerShiftZ() const
Get Z distance between two layers.
Definition EKLMGeometry.h:1807
Belle2::EKLMGeometry::getStripGeometry
const StripGeometry * getStripGeometry() const
Get strip geometry data.
Definition EKLMGeometry.h:1871
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::ExtHit
Store one Ext hit as a ROOT object.
Definition ExtHit.h:31
Belle2::KLMChannelIndex
KLM channel index.
Definition KLMChannelIndex.h:23
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::Instance
static const KLMElementNumbers & Instance()
Instantiation.
Definition KLMElementNumbers.cc:29
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 information 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::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::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::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::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::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::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:948
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:1216
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:1450
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:1772
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:1192
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:1586
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:1143
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:1482
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:862
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 separation.
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:1418
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:798
Belle2::TrackExtrapolateG4e::fromG4eToPhasespace
void fromG4eToPhasespace(const G4ErrorFreeTrajState &, G4ErrorSymMatrix &)
Convert the geant4e 5x5 covariance to phasespace 6x6 covariance.
Definition TrackExtrapolateG4e.cc:1041
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:1649
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
Belle2::BFieldManager::getField
static void getField(const double *pos, double *field)
return the magnetic field at a given position.
Definition BFieldManager.h:90
Belle2::ExtHitStatus
ExtHitStatus
Define state of extrapolation for each recorded hit.
Definition ExtHit.h:26
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