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