10#include <tracking/trackExtrapolateG4e/TrackExtrapolateG4e.h>
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>
47#include <CLHEP/Matrix/Vector.h>
48#include <CLHEP/Units/PhysicalConstants.h>
49#include <CLHEP/Units/SystemOfUnits.h>
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>
61#include <G4StepPoint.hh>
63#include <G4UImanager.hh>
64#include <G4VPhysicalVolume.hh>
65#include <G4VTouchable.hh>
66#include <G4ErrorSymMatrix.hh>
76#define TWOPI (2.0*M_PI)
77#define PI_8 (0.125*M_PI)
152 std::vector<Const::ChargedStable>& hypotheses)
166 m_MinPt = std::max(0.0, minPt) * CLHEP::GeV;
167 m_MinKE = std::max(0.0, minKE) * CLHEP::GeV;
182 B2FATAL(
"Coil geometry data are not available.");
192 G4ErrorPropagatorData::GetErrorPropagatorData()->SetTarget(
m_TargetExt);
195 B2FATAL(
"Beam pipe geometry data are not available.");
196 double beampipeRadius =
m_BeamPipeGeo->getParameter(
"Lv2OutBe.R2") * CLHEP::cm;
198 beampipeRadius =
m_BeamPipeGeo->getParameter(
"Lv2OutBe.R2") * CLHEP::cm;
205 double maxKLMTrackClusterDistance,
double maxECLTrackClusterDistance,
206 double minPt,
double minKE,
bool addHitsToRecoTrack,
207 std::vector<Const::ChargedStable>& hypotheses)
245 m_MinPt = std::max(0.0, minPt) * CLHEP::GeV;
246 m_MinKE = std::max(0.0, minKE) * CLHEP::GeV;
261 B2FATAL(
"Coil geometry data are not available.");
271 G4ErrorPropagatorData::GetErrorPropagatorData()->SetTarget(
m_TargetExt);
274 B2FATAL(
"Beam pipe geometry data are not available.");
275 double beampipeRadius =
m_BeamPipeGeo->getParameter(
"Lv2OutBe.R2") * CLHEP::cm;
277 beampipeRadius =
m_BeamPipeGeo->getParameter(
"Lv2OutBe.R2") * CLHEP::cm;
292 G4ErrorPropagatorData::GetErrorPropagatorData()->SetTarget(
m_TargetMuid);
311 width =
module->getPhiStripWidth();
313 width =
module->getZStripWidth();
320 double phi = M_PI_4 * (sector - 1);
328 bklmGeometry->
getActiveMiddleRadius(klmLayer.getSection(), klmLayer.getSector(), klmLayer.getLayer());
348 B2DEBUG(20, (byMuid ?
"muid" :
"ext"));
351 B2FATAL(
"KLM channel status data are not available.");
353 B2FATAL(
"KLM strip efficiency data are not available.");
355 B2FATAL(
"KLM likelihood parameters are not available.");
360 delete muidBuilder.second;
365 for (
int pdg : muidPdgCodes)
374 if (G4StateManager::GetStateManager()->GetCurrentState() == G4State_Idle) {
375 G4StateManager::GetStateManager()->SetNewState(G4State_GeomClosed);
378 G4ThreeVector directionAtIP, positionG4e, momentumG4e;
379 G4ErrorTrajErr covG4e(5);
384 G4ErrorPropagatorData::GetErrorPropagatorData()->SetTarget(
m_TargetMuid);
385 std::vector<std::pair<ECLCluster*, G4ThreeVector> > eclClusterInfo(
m_eclClusters.getEntries());
388 eclClusterInfo[c].second = G4ThreeVector(
m_eclClusters[c]->getClusterPosition().X(),
392 std::vector<std::pair<KLMCluster*, G4ThreeVector> > klmClusterInfo(
m_klmClusters.getEntries());
395 klmClusterInfo[c].second = G4ThreeVector(
m_klmClusters[c]->getClusterPosition().X(),
400 std::vector<std::map<const Track*, double> > bklmHitUsed(
m_klmHit2ds.getEntries());
403 int pdgCode = hypothesis.getPDGCode();
406 G4ErrorFreeTrajState g4eState(
"g4e_mu+", G4ThreeVector(), G4ThreeVector());
408 swim(extState, g4eState, &eclClusterInfo, &klmClusterInfo, &bklmHitUsed);
412 G4ErrorPropagatorData::GetErrorPropagatorData()->SetTarget(
m_TargetExt);
415 int pdgCode = hypothesis.getPDGCode();
417 G4ErrorFreeTrajState g4eState(
"g4e_mu+", G4ThreeVector(), G4ThreeVector());
419 swim(extState, g4eState);
440 delete muidBuilder.second;
458 const G4ThreeVector& position,
459 const G4ThreeVector& momentum,
460 const G4ErrorSymMatrix& covariance)
462 bool isCosmic =
false;
467 extMgr->
Initialize(
"Ext",
"default", 0.0, 0.25,
false, 0, std::vector<std::string>());
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");
480 if (G4StateManager::GetStateManager()->GetCurrentState() == G4State_Idle) {
481 G4StateManager::GetStateManager()->SetNewState(G4State_GeomClosed);
484 G4ErrorPropagatorData::GetErrorPropagatorData()->SetTarget(
m_TargetExt);
489 G4ThreeVector positionG4e = position * CLHEP::cm;
490 G4ThreeVector momentumG4e = momentum * CLHEP::GeV;
492 momentumG4e = -momentumG4e;
493 G4ErrorSymMatrix covarianceG4e(5, 0);
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,
498 momentumG4e.unit(), 0.0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
false
500 swim(extState, 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)
511 if (g4eState.GetMomentum().perp() <=
m_MinPt)
513 if (
m_TargetMuid->GetDistanceFromPoint(g4eState.GetPosition()) < 0.0)
515 G4ParticleDefinition* particle = G4ParticleTable::GetParticleTable()->FindParticle(extState.
pdgCode);
516 double mass = particle->GetPDGMass();
517 double minPSq = (mass +
m_MinKE) * (mass +
m_MinKE) - mass * mass;
519 std::vector<double> eclClusterDistance;
520 std::vector<ExtHit> eclHit1, eclHit2, eclHit3;
521 if (eclClusterInfo !=
nullptr) {
522 eclClusterDistance.resize(eclClusterInfo->size(), 1.0E10);
523 ExtHit tempExtHit(.0, extState.
pdgCode, Const::EDetector::ECL, 0, EXT_FIRST,
525 G4ThreeVector(), G4ThreeVector(), G4ErrorSymMatrix(6));
526 eclHit1.resize(eclClusterInfo->size(), tempExtHit);
527 eclHit2.resize(eclClusterInfo->size(), tempExtHit);
528 eclHit3.resize(eclClusterInfo->size(), tempExtHit);
531 std::vector<TrackClusterSeparation> klmHit;
532 if (klmClusterInfo !=
nullptr) {
533 klmHit.resize(klmClusterInfo->size());
537 if (extState.
track !=
nullptr)
539 G4ErrorMode propagationMode = (extState.
isCosmic ? G4ErrorMode_PropBackwards : G4ErrorMode_PropForwards);
540 m_ExtMgr->InitTrackPropagation(propagationMode);
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();
552 G4ThreeVector mom = track->GetMomentum();
555 if (step->GetStepLength() > 0.0) {
556 double dt = step->GetDeltaTime();
557 double dl = step->GetStepLength() / track->GetMaterial()->GetRadlen();
558 if (preStatus == fGeomBoundary) {
559 if (
m_TargetExt->GetDistanceFromPoint(pos) < 0.0) {
561 createExtHit(EXT_ENTER, extState, g4eState, preStepPoint, preTouch);
572 if (postStatus == fGeomBoundary) {
573 if (
m_TargetExt->GetDistanceFromPoint(pos) < 0.0) {
575 createExtHit(EXT_EXIT, extState, g4eState, postStepPoint, preTouch);
579 if (
createMuidHit(extState, g4eState, klmMuidLikelihood, bklmHitUsed)) {
581 m_ExtMgr->GetPropagator()->SetStepN(0);
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();
591 if (distance < eclClusterDistance[c]) {
592 eclClusterDistance[c] = distance;
593 G4ErrorSymMatrix covariance(6, 0);
595 eclHit3[c].update(EXT_ECLNEAR, extState.
tof, pos / CLHEP::cm, mom / CLHEP::GeV, covariance);
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);
606 G4ErrorSymMatrix covariance(6, 0);
608 eclHit1[c].update(EXT_ECLCROSS, tof, midPos / CLHEP::cm, mom / CLHEP::GeV, covariance);
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);
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);
626 eclHit2[c].update(EXT_ECLDL, length, midPos / CLHEP::cm, mom / CLHEP::GeV, covariance);
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));
649 if (errCode || (mom.mag2() < minPSq)) {
662 m_ExtMgr->EventTermination(propagationMode);
666 if (eclClusterInfo !=
nullptr) {
667 for (
unsigned int c = 0; c < eclClusterInfo->size(); ++c) {
668 if (eclHit1[c].getStatus() != EXT_FIRST) {
670 (*eclClusterInfo)[c].first->addRelationTo(h);
671 if (extState.
track !=
nullptr) {
675 if (eclHit2[c].getStatus() != EXT_FIRST) {
677 (*eclClusterInfo)[c].first->addRelationTo(h);
678 if (extState.
track !=
nullptr) {
682 if (eclHit3[c].getStatus() != EXT_FIRST) {
684 (*eclClusterInfo)[c].first->addRelationTo(h);
685 if (extState.
track !=
nullptr) {
692 if (klmClusterInfo !=
nullptr) {
696 unsigned int closestCluster = 0;
697 for (
unsigned int c = 0; c < klmClusterInfo->size(); ++c) {
698 if (klmHit[c].getDistance() > 1.0E9) {
703 (*klmClusterInfo)[c].first->addRelationTo(h);
704 if (extState.
track !=
nullptr) {
708 if (klmHit[c].getDistance() < minDistance) {
710 minDistance = klmHit[c].getDistance();
716 if (extState.
track !=
nullptr) {
717 extState.
track->
addRelationTo((*klmClusterInfo)[closestCluster].first, 1. / (minDistance / CLHEP::cm));
729 if (trackClusterSeparations.size() == 0) {
737 const auto closestSeparationIterator = std::min_element(trackClusterSeparations.begin(), trackClusterSeparations.end(),
738 [](
const auto & a,
const auto & b) {
739 return a.getDistance() < b.getDistance();
743 const auto& closestSeparation = *closestSeparationIterator;
744 klmCluster.setClusterTrackSeparation(closestSeparation.getDistance() / CLHEP::cm);
745 klmCluster.setClusterTrackSeparationAngle(closestSeparation.getTrackClusterSeparationAngle());
746 klmCluster.setClusterTrackRotationAngle(closestSeparation.getTrackRotationAngle());
755 if (g4eState.GetMomentum().perp() <=
m_MinPt)
757 if (
m_TargetExt->GetDistanceFromPoint(g4eState.GetPosition()) < 0.0)
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);
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();
775 G4ThreeVector mom = track->GetMomentum();
779 if (preStatus == fUndefined) {
781 createExtHit(EXT_FIRST, extState, g4eState, preStepPoint, preTouch);
785 if (step->GetStepLength() > 0.0) {
786 double dt = step->GetDeltaTime();
787 double dl = step->GetStepLength() / track->GetMaterial()->GetRadlen();
788 if (preStatus == fGeomBoundary) {
790 createExtHit(EXT_ENTER, extState, g4eState, preStepPoint, preTouch);
801 if (postStatus == fGeomBoundary) {
803 createExtHit(EXT_EXIT, extState, g4eState, postStepPoint, preTouch);
808 if (errCode || (mom.mag2() < minPSq)) {
810 createExtHit(EXT_STOP, extState, g4eState, postStepPoint, preTouch);
815 if (
m_TargetExt->GetDistanceFromPoint(pos) < 0.0) {
817 createExtHit(EXT_ESCAPE, extState, g4eState, postStepPoint, preTouch);
827 m_ExtMgr->EventTermination(propagationMode);
834 G4PhysicalVolumeStore* pvStore = G4PhysicalVolumeStore::GetInstance();
835 if (pvStore->size() == 0) {
836 B2FATAL(
"No geometry defined. Please create the geometry first.");
840 m_EnterExit =
new std::map<G4VPhysicalVolume*, enum VolTypes>;
842 for (std::vector<G4VPhysicalVolume*>::iterator iVol = pvStore->begin();
843 iVol != pvStore->end(); ++iVol) {
844 const G4String name = (*iVol)->GetName();
847 if (name.find(
"TOPModule") != std::string::npos) {
851 else if (name.find(
"_TOPPrism_") != std::string::npos ||
852 name.find(
"_TOPBarSegment") != std::string::npos ||
853 name.find(
"_TOPMirrorSegment") != std::string::npos) {
857 else if (name.find(
"TOPBarSegment1Glue") != std::string::npos ||
858 name.find(
"TOPBarSegment2Glue") != std::string::npos ||
859 name.find(
"TOPMirrorSegmentGlue") != std::string::npos) {
862 }
else if (name ==
"ARICH.AerogelSupportPlate") {
864 }
else if (name ==
"ARICH.AerogelImgPlate") {
866 }
else if (name.find(
"ARICH.HAPDWindow") != std::string::npos) {
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) {
877 else if (name.compare(0, 5,
"BKLM.") == 0) {
878 if (name.find(
"GasPhysical") != std::string::npos) {
880 }
else if (name.find(
"ScintActiveType") != std::string::npos) {
882 }
else if ((name.find(
"ScintType") != std::string::npos) ||
883 (name.find(
"ElectrodePhysical") != std::string::npos)) {
888 else if (name.compare(0, 14,
"StripSensitive") == 0) {
900 detID = Const::EDetector::invalidDetector;
903 G4VPhysicalVolume* pv = touch->GetVolume(0);
904 std::map<G4VPhysicalVolume*, enum VolTypes>::iterator it =
m_EnterExit->find(pv);
908 switch (it->second) {
910 detID = Const::EDetector::CDC;
911 copyID = pv->GetCopyNo();
914 detID = Const::EDetector::TOP;
915 copyID = -(pv->GetCopyNo());
918 detID = Const::EDetector::TOP;
919 if (touch->GetHistoryDepth() >= 1)
920 copyID = touch->GetVolume(1)->GetCopyNo();
923 detID = Const::EDetector::TOP;
924 if (touch->GetHistoryDepth() >= 2)
925 copyID = touch->GetVolume(2)->GetCopyNo();
928 detID = Const::EDetector::ARICH;
932 detID = Const::EDetector::ARICH;
936 detID = Const::EDetector::ARICH;
937 if (touch->GetHistoryDepth() >= 2)
938 copyID = touch->GetVolume(2)->GetCopyNo();
941 detID = Const::EDetector::ECL;
945 detID = Const::EDetector::BKLM;
946 if (touch->GetHistoryDepth() == DEPTH_RPC) {
948 int layer = touch->GetCopyNumber(4);
949 int sector = touch->GetCopyNumber(6);
950 int section = touch->GetCopyNumber(7);
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);
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());
984 ExtState extState = {&b2track, pdgCode,
false, 0.0, 0.0,
985 G4ThreeVector(0, 0, 1), 0.0, 0, 0, 0, -1, -1, -1, -1, 0, 0,
false
988 if (recoTrack ==
nullptr) {
989 B2WARNING(
"Track without associated RecoTrack: skipping extrapolation for this track.");
996 B2WARNING(
"RecoTrack fit failed for cardinal representation: skipping extrapolation for this track.");
1000 int charge = int(trackRep->getPDGCharge());
1006 TVector3 firstPosition, firstMomentum, lastPosition, lastMomentum;
1007 TMatrixDSym firstCov(6), lastCov(6);
1010 trackRep->getPosMomCov(firstState, firstPosition, firstMomentum, firstCov);
1012 trackRep->getPosMomCov(lastState, lastPosition, lastMomentum, lastCov);
1014 extState.
tof = lastState.getTime();
1015 if (lastPosition.Mag2() < firstPosition.Mag2()) {
1016 firstPosition = lastPosition;
1017 firstMomentum = -lastMomentum;
1019 trackRep->getPosMomCov(firstState, lastPosition, lastMomentum, lastCov);
1020 lastMomentum *= -1.0;
1022 extState.
tof = firstState.getTime();
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()));
1032 extState.
directionAtIP.set(firstMomentum.Unit().X(), firstMomentum.Unit().Y(), firstMomentum.Unit().Z());
1034 double radius = (firstMomentum.Perp() * CLHEP::GeV / CLHEP::eV) / (CLHEP::c_light * charge *
m_MagneticField);
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;
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;
1056 G4ThreeVector posG4e(lastPosition.X() * CLHEP::cm, lastPosition.Y() * CLHEP::cm,
1057 lastPosition.Z() * CLHEP::cm);
1058 G4ThreeVector momG4e(lastMomentum.X() * CLHEP::GeV, lastMomentum.Y() * CLHEP::GeV,
1059 lastMomentum.Z() * CLHEP::GeV);
1060 G4ErrorSymMatrix covG4e(5, 0);
1062 g4eState.SetData(
"g4e_" + particle->GetParticleName(), posG4e, momG4e);
1063 g4eState.SetParameters(posG4e, 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() <<
")");
1087 G4ErrorFreeTrajParam param = g4eState.GetParameters();
1088 double p = 1.0 / (param.GetInvP() * CLHEP::GeV);
1090 double lambda = param.GetLambda();
1091 double sinLambda = std::sin(lambda);
1092 double cosLambda = std::cos(lambda);
1093 double phi = param.GetPhi();
1094 double sinPhi = std::sin(phi);
1095 double cosPhi = std::cos(phi);
1099 G4ErrorMatrix jacobian(6, 5, 0);
1101 jacobian(4, 1) = -pSq * cosLambda * cosPhi;
1102 jacobian(5, 1) = -pSq * cosLambda * sinPhi;
1103 jacobian(6, 1) = -pSq * sinLambda;
1105 jacobian(4, 2) = -p * sinLambda * cosPhi;
1106 jacobian(5, 2) = -p * sinLambda * sinPhi;
1107 jacobian(6, 2) = p * cosLambda;
1109 jacobian(4, 3) = -p * cosLambda * sinPhi;
1110 jacobian(5, 3) = p * cosLambda * cosPhi;
1112 jacobian(1, 4) = -sinPhi;
1113 jacobian(2, 4) = cosPhi;
1115 jacobian(1, 5) = -sinLambda * cosPhi;
1116 jacobian(2, 5) = -sinLambda * sinPhi;
1117 jacobian(3, 5) = cosLambda;
1119 G4ErrorTrajErr g4eCov = g4eState.GetError();
1120 covariance.assign(g4eCov.similarity(jacobian));
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];
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);
1153 G4ErrorMatrix jacobian(5, 6, 0);
1155 jacobian(1, 4) = -pInvSq * cosLambda * cosPhi;
1156 jacobian(1, 5) = -pInvSq * cosLambda * sinPhi;
1157 jacobian(1, 6) = -pInvSq * sinLambda;
1159 jacobian(2, 4) = -pInv * sinLambda * cosPhi;
1160 jacobian(2, 5) = -pInv * sinLambda * sinPhi;
1161 jacobian(2, 6) = pInv * cosLambda;
1163 jacobian(3, 4) = -pPerpInv * sinPhi;
1164 jacobian(3, 5) = pPerpInv * cosPhi;
1166 jacobian(4, 1) = -sinPhi;
1167 jacobian(4, 2) = cosPhi;
1169 jacobian(5, 1) = -sinLambda * cosPhi;
1170 jacobian(5, 2) = -sinLambda * sinPhi;
1171 jacobian(5, 3) = cosLambda;
1173 covG4e = temp.similarity(jacobian);
1178 G4ErrorTrajErr& covG4e)
1190 G4ErrorSymMatrix temp(covariance);
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);
1202 G4ErrorMatrix jacobian(5, 6, 0);
1204 jacobian(1, 4) = -pInvSq * cosLambda * cosPhi;
1205 jacobian(1, 5) = -pInvSq * cosLambda * sinPhi;
1206 jacobian(1, 6) = -pInvSq * sinLambda;
1208 jacobian(2, 4) = -pInv * sinLambda * cosPhi;
1209 jacobian(2, 5) = -pInv * sinLambda * sinPhi;
1210 jacobian(2, 6) = pInv * cosLambda;
1212 jacobian(3, 4) = -pPerpInv * sinPhi;
1213 jacobian(3, 5) = pPerpInv * cosPhi;
1215 jacobian(4, 1) = -sinPhi;
1216 jacobian(4, 2) = cosPhi;
1218 jacobian(5, 1) = -sinLambda * cosPhi;
1219 jacobian(5, 2) = -sinLambda * sinPhi;
1220 jacobian(5, 3) = cosLambda;
1221 covG4e = temp.similarity(jacobian);
1227 const G4ErrorFreeTrajState& g4eState,
1228 const G4StepPoint* stepPoint,
const G4TouchableHandle& touch)
1233 G4ThreeVector pos(stepPoint->GetPosition() / CLHEP::cm);
1234 G4ThreeVector mom(stepPoint->GetMomentum() / CLHEP::GeV);
1237 G4ErrorSymMatrix covariance(6, 0);
1241 pos, mom, covariance);
1243 if (extState.
track !=
nullptr)
1251 std::vector<std::map<const Track*, double> >* bklmHitUsed)
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();
1270 (*bklmHitUsed)[intersection.
hit].insert(std::pair<const Track*, double>(extState.
track, intersection.
chi2));
1272 float layerBarrelEfficiency = 1.;
1276 intersection.
sector + 1, intersection.
layer + 1, plane, 1);
1290 intersection.
chi2 = -1.0;
1295 g4eState.GetG4Track()->GetVolume());
1301 int sector = intersection.
sector + 1;
1302 int layer = intersection.
layer + 1;
1305 const CLHEP::Hep3Vector localPosition = m->globalToLocal(intersection.
position);
1306 int zStrip = m->getZStripNumber(localPosition);
1307 int phiStrip = m->getPhiStripNumber(localPosition);
1308 if (zStrip >= 0 && phiStrip >= 0) {
1309 uint16_t channel1, channel2;
1311 section, sector, layer,
1314 section, sector, layer,
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));
1333 float layerBarrelEfficiency = 1.;
1337 intersection.
sector + 1, intersection.
layer + 1, plane, 1);
1358 float layerEndcapEfficiency = 1.;
1362 intersection.
sector + 1, intersection.
layer + 1, plane, 1);
1376 intersection.
chi2 = -1.0;
1380 int result, strip1, strip2;
1382 intersection.
position, &strip1, &strip2);
1384 uint16_t channel1, channel2;
1392 B2ERROR(
"Incomplete KLM channel status data.");
1398 float layerEndcapEfficiency = 1.;
1402 intersection.
sector + 1, intersection.
layer + 1, plane, 1);
1417 if (intersection.
chi2 >= 0.0) {
1418 ROOT::Math::XYZVector tpos(intersection.
position.x(),
1427 tpos, tposAtHitPlane,
1428 extState.
tof, intersection.
time,
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;
1441 g4eState.SetError(covG4e);
1442 extState.
chi2 += intersection.
chi2;
1458 double phi = intersection.
position.phi();
1461 if (phi > TWOPI - PI_8)
1463 int sector = (int)((phi + PI_8) / M_PI_4);
1476 intersection.
layer = layer;
1477 intersection.
sector = sector;
1493 double oldZ = std::fabs(oldPosition.z() -
m_OffsetZ);
1498 for (
int layer = extState.
firstEndcapLayer; layer <= outermostLayer; ++layer) {
1507 intersection.
layer = layer;
1509 isForward ? 2 : 1, intersection.
position) - 1;
1519 G4ThreeVector extPos0(intersection.
position);
1520 double diffBestMagSq = 1.0E60;
1522 int matchingLayer = intersection.
layer + 1;
1524 for (
int h = 0; h <
m_klmHit2ds.getEntries(); ++h) {
1528 if (hit->getLayer() != matchingLayer)
1530 if (hit->isOutOfTime())
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;
1538 if (std::fabs(dn) < 2.0) {
1541 if (diff.mag2() < diffBestMagSq) {
1542 diffBestMagSq = diff.mag2();
1548 if (std::fabs(dn) > 50.0)
1550 int sector = hit->getSector() - 1;
1551 int dSector = abs(intersection.
sector - sector);
1560 dn = diff * nHit + dn2;
1561 if (std::fabs(dn) > 1.0)
1564 G4ThreeVector extDir(intersection.
momentum.unit());
1565 double extDirA = extDir * nHit;
1566 if (std::fabs(extDirA) < 1.0E-6)
1568 G4ThreeVector projection = extDir * (dn2 / extDirA);
1569 if (projection.mag() > 15.0)
1571 diff += projection - nHit * dn;
1572 if (diff.mag2() < diffBestMagSq) {
1573 diffBestMagSq = diff.mag2();
1575 extPos0 = intersection.
position - projection;
1582 intersection.
isForward = (hit->getSection() == 1);
1583 intersection.
sector = hit->getSector() - 1;
1584 intersection.
time = hit->getTime();
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;
1591 nStrips = hit->getZStripMax() - hit->getZStripMin() + 1;
1593 factor = std::pow((0.9 + 0.4 * dn * dn), 1.5) * 0.55;
1596 int nStrips = hit->getPhiStripMax() - hit->getPhiStripMin() + 1;
1597 localVariance[0] *= (nStrips * nStrips);
1598 nStrips = hit->getZStripMax() - hit->getZStripMin() + 1;
1599 localVariance[1] *= (nStrips * nStrips);
1601 G4ThreeVector hitPos(hit->getPositionX(), hit->getPositionY(),
1602 hit->getPositionZ());
1604 if (intersection.
chi2 >= 0.0) {
1605 intersection.
hit = bestHit;
1606 hit->isOnTrack(
true);
1607 if (track !=
nullptr) {
1608 track->addRelationTo(hit, intersection.
chi2);
1616 return intersection.
chi2 >= 0.0;
1622 double diffBestMagSq = 1.0E60;
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) {
1631 if (hit->getLayer() != matchingLayer)
1633 if (hit->getSection() != matchingEndcap)
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;
1643 if (std::fabs(dn) > 2.0)
1646 if (diff.mag2() < diffBestMagSq) {
1647 diffBestMagSq = diff.mag2();
1654 intersection.
hit = bestHit;
1655 intersection.
isForward = (hit->getSection() == 2);
1656 intersection.
sector = hit->getSector() - 1;
1657 intersection.
time = hit->getTime();
1659 int nStrips = hit->getXStripMax() - hit->getXStripMin() + 1;
1660 localVariance[0] *= (nStrips * nStrips);
1661 nStrips = hit->getYStripMax() - hit->getYStripMin() + 1;
1662 localVariance[1] *= (nStrips * nStrips);
1663 G4ThreeVector hitPos(hit->getPositionX(), hit->getPositionY(), hit->getPositionZ());
1665 if (intersection.
chi2 >= 0.0) {
1668 if (track !=
nullptr) {
1669 track->addRelationTo(hit, intersection.
chi2);
1679 return intersection.
chi2 >= 0.0;
1684 const G4ThreeVector& hitPos,
const G4ThreeVector& extPos0)
1704 G4ThreeVector extPos(extPos0);
1705 G4ThreeVector extMom(intersection.
momentum);
1706 G4ThreeVector extDir(extMom.unit());
1707 G4ThreeVector diffPos(hitPos - extPos);
1708 G4ErrorSymMatrix extCov(intersection.
covariance);
1710 G4ErrorMatrix extPar(6, 1);
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();
1723 nC = G4ThreeVector(0.0, 0.0, 1.0);
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);
1731 double extDirA = extDir * nA;
1732 if (std::fabs(extDirA) < 1.0E-6)
1734 double extDirBA = extDir * nB / extDirA;
1735 double extDirCA = extDir * nC / extDirA;
1738 G4ThreeVector move = extDir * ((diffPos * nA) / extDirA);
1743 G4ErrorMatrix jacobian(2, 6);
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;
1751 G4ErrorMatrix residual(2, 1);
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];
1755 G4ErrorSymMatrix hitCov(2, 0);
1756 hitCov[0][0] = localVariance[0];
1757 hitCov[1][1] = localVariance[1];
1760 hitCov[0][0] *= 10.0;
1761 hitCov[1][1] *= 10.0;
1765 G4ErrorSymMatrix correction(extCov.similarity(jacobian) + hitCov);
1772 correction.invert(fail);
1778 intersection.
chi2 = (correction.similarityT(residual))[0][0];
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]);
1787 correction = hitCov - extCov.similarity(jacobian);
1788 correction.invert(fail);
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];
1801 intersection.
position = extPos + extDir * (((intersection.
position - extPos) * nA) / extDirA);
1827 if (outcome != MuidElementNumbers::c_NotReached) {
1829 int charge = klmMuidLikelihood->
getCharge();
1831 std::map<int, double> mapPdgPDF;
1832 for (
int pdg : signedPdgVector) {
1835 B2FATAL(
"Something went wrong: PDF for PDG code " << pdg <<
" not found!");
1836 double pdf = (search->second)->getPDF(klmMuidLikelihood);
1838 mapPdgPDF.insert(std::pair<int, double>(std::abs(pdg), pdf));
1840 if (denom < 1.0E-20)
1843 for (
auto const& [pdg, pdf] : mapPdgPDF) {
1844 klmMuidLikelihood->
setPDFValue(pdf, std::abs(pdg));
1846 klmMuidLikelihood->
setLogL(std::log(pdf), std::abs(pdg));
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).
@ c_ForwardSection
Forward.
@ c_BackwardSection
Backward.
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.
static const ChargedStable muon
muon particle
EDetector
Enum for identifying the detector components (detector and subdetector).
static const double doubleNaN
quiet_NaN
static const ChargedStable electron
electron particle
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.
@ c_ForwardSection
Forward.
@ c_BackwardSection
Backward.
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.
static const GeometryData & Instance(enum DataSource dataSource=c_Database, const GearDir *gearDir=nullptr)
Instantiation.
Store one Ext hit as a ROOT object.
@ c_IndexLevelLayer
Layer.
ChannelStatus
Channel status.
@ c_Dead
Dead channel (no signal).
@ c_Unknown
Unknown status (no data).
static const KLMElementNumbers & Instance()
Instantiation.
Store one muon-identification hit in the KLM as a ROOT object.
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.
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.
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.
bool wasFitSuccessful(const genfit::AbsTrackRep *representation=nullptr) const
Returns true if the last fit with the given representation was successful.
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.
genfit::AbsTrackRep * getCardinalRepresentation() const
Get a pointer to the cardinal track representation. You are not allowed to modify or delete it!
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...
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.
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...
unsigned int getNumberOfTotalHits() const
Return the number of cdc + svd + pxd + bklm + eklm hits.
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.
static ExtManager * GetManager()
Get pointer to the instance of this singleton class (create if needed)
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.
Class that bundles various TrackFitResults.
static const double T
[tesla]
Provides BKLM geometry parameters for simulation, reconstruction etc (from Gearbox or DataBase)
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.
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.
@ VOLTYPE_ARICH2
ARICH Img plate.
@ VOLTYPE_TOP2
TOP quartz.
@ VOLTYPE_ARICH3
ARICH HAPD window.
@ VOLTYPE_BKLM2
BKLM scintillator.
@ 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