diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx index e1053fa6f84d9242e6010ba756981f9ccc05d9dd..41e74b75b8f3c770709b64e400b5600bd6c1c1d1 100644 --- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx +++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx @@ -106,6 +106,7 @@ StatusCode NtupleDumperAlg::initialize() ATH_CHECK(m_preshowerCalibratedContainer.initialize()); ATH_CHECK(m_ecalCalibratedContainer.initialize()); + ATH_CHECK(m_eventInfoKey.initialize()); ATH_CHECK(detStore()->retrieve(m_sctHelper, "FaserSCT_ID")); ATH_CHECK(detStore()->retrieve(m_vetoNuHelper, "VetoNuID")); @@ -160,8 +161,8 @@ StatusCode NtupleDumperAlg::initialize() //WAVEFORMS addWaveBranches("VetoNu",2,4); - addWaveBranches("VetoSt1",2,6); - addWaveBranches("VetoSt2",1,14); + addWaveBranches("VetoSt1",1,14); + addWaveBranches("VetoSt2",2,6); addWaveBranches("Timing",4,8); addWaveBranches("Preshower",2,12); addWaveBranches("Calo",4,0); @@ -197,6 +198,7 @@ StatusCode NtupleDumperAlg::initialize() m_tree->Branch("TrackSegment_pz", &m_trackseg_pz); //TrackCollection + m_tree->Branch("Track_PropagationError", &m_propagationError, "Track_PropagationError/I"); m_tree->Branch("longTracks", &m_longTracks, "longTracks/I"); m_tree->Branch("Track_Chi2", &m_Chi2); m_tree->Branch("Track_nDoF", &m_DoF); @@ -304,6 +306,24 @@ StatusCode NtupleDumperAlg::initialize() m_tree->Branch("truthPdg", &m_truthPdg, "truthPdg/I"); m_tree->Branch("CrossSection", &m_crossSection, "crossSection/D"); + + // first 10 truth particles + + m_tree->Branch("truth_P", &m_truth_P); + m_tree->Branch("truth_px", &m_truth_px); + m_tree->Branch("truth_py", &m_truth_py); + m_tree->Branch("truth_pz", &m_truth_pz); + m_tree->Branch("truth_m", &m_truth_m); + m_tree->Branch("truth_pdg", &m_truth_pdg); + + m_tree->Branch("truth_prod_x", &m_truth_prod_x); + m_tree->Branch("truth_prod_y", &m_truth_prod_y); + m_tree->Branch("truth_prod_z", &m_truth_prod_z); + + m_tree->Branch("truth_dec_x", &m_truth_dec_x); + m_tree->Branch("truth_dec_y", &m_truth_dec_y); + m_tree->Branch("truth_dec_z", &m_truth_dec_z); + // for mother + daughter particle truth infomation m_tree->Branch("truthM_P", &m_truthM_P); @@ -546,30 +566,64 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const SG::ReadHandle<xAOD::TruthParticleContainer> truthParticleContainer { m_truthParticleContainer, ctx }; if (truthParticleContainer.isValid() && truthParticleContainer->size() > 0) { - for (auto particle : *truthParticleContainer) - { - if ( particle->barcode() == 1 || particle->barcode() == 2 || particle->barcode() == 3 ) - { - - if ( particle->pdgId() == 32) // mother particle (A') - { - m_truthM_P.push_back(particle->p4().P()); - m_truthM_px.push_back(particle->p4().X()); - m_truthM_py.push_back(particle->p4().Y()); - m_truthM_pz.push_back(particle->p4().Z()); - - if ( particle->hasProdVtx()) { - m_truthM_x.push_back(particle->prodVtx()->x()); - m_truthM_y.push_back(particle->prodVtx()->y()); - m_truthM_z.push_back(particle->prodVtx()->z()); - } else { + int ipart(0); + for (auto particle : *truthParticleContainer) + { + + // loop over first 10 truth particles (for non A' samples) + + if (ipart++ < 10) { + + m_truth_P.push_back(particle->p4().P()); + m_truth_px.push_back(particle->p4().X()); + m_truth_py.push_back(particle->p4().Y()); + m_truth_pz.push_back(particle->p4().Z()); + m_truth_m.push_back(particle->m()); + m_truth_pdg.push_back(particle->pdgId()); + + if ( particle->hasProdVtx()) { + m_truth_prod_x.push_back(particle->prodVtx()->x()); + m_truth_prod_y.push_back(particle->prodVtx()->y()); + m_truth_prod_z.push_back(particle->prodVtx()->z()); + } else { + m_truth_prod_x.push_back(NaN); + m_truth_prod_y.push_back(NaN); + m_truth_prod_z.push_back(NaN); + } + + if ( particle->hasDecayVtx()) { + m_truth_dec_x.push_back(particle->decayVtx()->x()); + m_truth_dec_y.push_back(particle->decayVtx()->y()); + m_truth_dec_z.push_back(particle->decayVtx()->z()); + } else { + m_truth_dec_x.push_back(NaN); + m_truth_dec_y.push_back(NaN); + m_truth_dec_z.push_back(NaN); + } + } + + + if ( particle->barcode() == 1 || particle->barcode() == 2 || particle->barcode() == 3 ) { + + if ( particle->pdgId() == 32) { // mother particle (A') + + m_truthM_P.push_back(particle->p4().P()); + m_truthM_px.push_back(particle->p4().X()); + m_truthM_py.push_back(particle->p4().Y()); + m_truthM_pz.push_back(particle->p4().Z()); + + if ( particle->hasDecayVtx()) { // decay vertex for A' particle + m_truthM_x.push_back(particle->decayVtx()->x()); + m_truthM_y.push_back(particle->decayVtx()->y()); + m_truthM_z.push_back(particle->decayVtx()->z()); + } else { m_truthM_x.push_back(NaN); m_truthM_y.push_back(NaN); m_truthM_z.push_back(NaN); } + } - } if ( particle->pdgId() == 11) // daughter particle (positron) { m_truthd0_P.push_back(particle->p4().P()); @@ -604,9 +658,9 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const m_truthd1_z.push_back(NaN); } } - } + } } - } + } } // load in calibrated calo container @@ -683,6 +737,14 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const } } + SG::ReadDecorHandle<xAOD::EventInfo,uint32_t> eventInfo (m_eventInfoKey, ctx); + if (eventInfo->errorState(xAOD::EventInfo_v1::SCT) == xAOD::EventInfo::Error) { + m_propagationError = 1; + ATH_MSG_DEBUG("NtupleDumper: xAOD::EventInfo::SCT::Error"); + } else { + m_propagationError = 0; + } + // get geometry context FaserActsGeometryContext faserGeometryContext = m_trackingGeometryTool->getNominalGeometryContext(); auto gctx = faserGeometryContext.context(); @@ -1229,7 +1291,20 @@ NtupleDumperAlg::clearTree() const m_truthBarcode = 0; m_truthPdg = 0; - m_truthM_P.clear(); + m_truth_P.clear(); + m_truth_px.clear(); + m_truth_py.clear(); + m_truth_pz.clear(); + m_truth_m.clear(); + m_truth_pdg.clear(); + m_truth_prod_x.clear(); + m_truth_prod_y.clear(); + m_truth_prod_z.clear(); + m_truth_dec_x.clear(); + m_truth_dec_y.clear(); + m_truth_dec_z.clear(); + + m_truthM_P.clear(); m_truthM_px.clear(); m_truthM_py.clear(); m_truthM_pz.clear(); diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h index c970480b2586af42fe046144905473aa505cbd35..c5b21b97622c83a6cf3dea34af4ad928951c8279 100644 --- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h +++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h @@ -21,6 +21,8 @@ #include "FaserActsKalmanFilter/IFiducialParticleTool.h" #include "FaserActsKalmanFilter/ITrackTruthMatchingTool.h" #include "TrackerSimEvent/FaserSiHitCollection.h" +#include "xAODEventInfo/EventInfo.h" +#include "StoreGate/ReadDecorHandle.h" #include <vector> @@ -84,7 +86,8 @@ private: SG::ReadHandleKey<xAOD::FaserTriggerData> m_FaserTriggerData { this, "FaserTriggerDataKey", "FaserTriggerData", "ReadHandleKey for xAOD::FaserTriggerData"}; SG::ReadHandleKey<xAOD::WaveformClock> m_ClockWaveformContainer { this, "WaveformClockKey", "WaveformClock", "ReadHandleKey for ClockWaveforms Container"}; - ToolHandle<IFaserActsExtrapolationTool> m_extrapolationTool { this, "ExtrapolationTool", "FaserActsExtrapolationTool" }; + SG::ReadDecorHandleKey<xAOD::EventInfo> m_eventInfoKey{this, "EventInfoKey", "EventInfo"}; + ToolHandle<IFaserActsExtrapolationTool> m_extrapolationTool { this, "ExtrapolationTool", "FaserActsExtrapolationTool" }; ToolHandle<IFaserActsTrackingGeometryTool> m_trackingGeometryTool {this, "TrackingGeometryTool", "FaserActsTrackingGeometryTool"}; ToolHandle<ITrackTruthMatchingTool> m_trackTruthMatchingTool {this, "TrackTruthMatchingTool", "TrackTruthMatchingTool"}; ToolHandle<IFiducialParticleTool> m_fiducialParticleTool {this, "FiducialParticleTool", "FiducialParticleTool"}; @@ -191,6 +194,7 @@ private: mutable std::vector<double> m_trackseg_pz; mutable int m_longTracks; + mutable int m_propagationError; mutable std::vector<double> m_Chi2; mutable std::vector<double> m_DoF; mutable std::vector<double> m_xup; @@ -276,7 +280,7 @@ private: mutable std::vector<double> m_truthM_py; mutable std::vector<double> m_truthM_pz; - mutable std::vector<double> m_truthM_x; + mutable std::vector<double> m_truthM_x; // decay vertex of A' mutable std::vector<double> m_truthM_y; mutable std::vector<double> m_truthM_z; @@ -286,7 +290,7 @@ private: mutable std::vector<double> m_truthd0_py; mutable std::vector<double> m_truthd0_pz; - mutable std::vector<double> m_truthd0_x; + mutable std::vector<double> m_truthd0_x; // production vertex for daughter particles mutable std::vector<double> m_truthd0_y; mutable std::vector<double> m_truthd0_z; @@ -300,6 +304,27 @@ private: mutable std::vector<double> m_truthd1_y; mutable std::vector<double> m_truthd1_z; + // first 10 truth particles + + mutable std::vector<double> m_truth_P; + mutable std::vector<double> m_truth_px; + mutable std::vector<double> m_truth_py; + mutable std::vector<double> m_truth_pz; + mutable std::vector<double> m_truth_m; + + mutable std::vector<double> m_truth_dec_x; // components of decay vertex (mm) + mutable std::vector<double> m_truth_dec_y; + mutable std::vector<double> m_truth_dec_z; + + mutable std::vector<double> m_truth_prod_x; // components of production vertex (mm) + mutable std::vector<double> m_truth_prod_y; + mutable std::vector<double> m_truth_prod_z; + + mutable std::vector<int> m_truth_pdg; // pdg of first 10 truth particles + + + + mutable double m_truthLeptonMomentum; mutable int m_truthBarcode; mutable int m_truthPdg; diff --git a/Tracking/Acts/FaserActsKalmanFilter/python/CKF2Config.py b/Tracking/Acts/FaserActsKalmanFilter/python/CKF2Config.py index 0c890e1eb0542037393884a0d181abe9799c6664..5d6ed5c4f3170cfc57b058d6bc9432329c1da561 100644 --- a/Tracking/Acts/FaserActsKalmanFilter/python/CKF2Config.py +++ b/Tracking/Acts/FaserActsKalmanFilter/python/CKF2Config.py @@ -107,9 +107,7 @@ def CKF2Cfg(flags, **kwargs): ckf.PerformanceWriter = False ckf.nMax = 10 - # Use larger chi2 cut until alignment improves - # ckf.chi2Max = 25 - ckf.chi2Max = 100000 + ckf.chi2Max = 25 # Protect against running out of memory on busy events ckf.maxSteps = 5000 acc.addEventAlgo(ckf) diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx index 3647ac98abdd9e5176fb0902fb368e89dc741dcb..f990be9458e0ee2ac0a5bdef8b670703f1424c26 100644 --- a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx +++ b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx @@ -2,6 +2,7 @@ #include "StoreGate/ReadHandle.h" #include "StoreGate/ReadCondHandleKey.h" +#include "StoreGate/WriteDecorHandle.h" #include "TrackerSpacePoint/FaserSCT_SpacePointCollection.h" #include "TrackerSpacePoint/FaserSCT_SpacePoint.h" #include "TrackerIdentifier/FaserSCT_ID.h" @@ -28,9 +29,9 @@ #include "Acts/Propagator/Propagator.hpp" #include "Acts/TrackFitting/GainMatrixSmoother.hpp" #include "Acts/TrackFitting/GainMatrixUpdater.hpp" - - #include "Acts/EventData/Measurement.hpp" +#include "Acts/Propagator/PropagatorError.hpp" +#include "Acts/TrackFinding/CombinatorialKalmanFilterError.hpp" size_t CKF2::TrajectoryInfo::nClusters {0}; @@ -49,6 +50,7 @@ StatusCode CKF2::initialize() { ATH_CHECK(m_kalmanFitterTool1.retrieve()); ATH_CHECK(m_createTrkTrackTool.retrieve()); ATH_CHECK(m_trackCollection.initialize()); + ATH_CHECK(m_eventInfoKey.initialize()); if (m_performanceWriter && !m_noDiagnostics) { ATH_CHECK(m_performanceWriterTool.retrieve()); } @@ -78,6 +80,8 @@ StatusCode CKF2::execute() { const EventContext& ctx = Gaudi::Hive::currentContext(); m_numberOfEvents++; + SG::WriteDecorHandle<xAOD::EventInfo,uint32_t> eventInfo (m_eventInfoKey, ctx); + SG::WriteHandle trackContainer{m_trackCollection, ctx}; std::unique_ptr<TrackCollection> outputTracks = std::make_unique<TrackCollection>(); @@ -150,6 +154,14 @@ StatusCode CKF2::execute() { std::list<TrajectoryInfo> allTrajectories; for (auto &result : results) { if (not result.ok()) { + // TODO use status bits for different errors + // result.error() == Acts::CombinatorialKalmanFilterError::NoTrackFound + if (result.error() == Acts::PropagatorError::StepCountLimitReached || + result.error() == Acts::CombinatorialKalmanFilterError::PropagationReachesMaxSteps) { + if (!eventInfo->updateErrorState(xAOD::EventInfo::SCT, xAOD::EventInfo::Error)) { + ATH_MSG_WARNING (" cannot set error state for SCT"); + } + } continue; } CKFResult ckfResult = result.value(); diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.h b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.h index 2111afcf6b52d46b45e1e2539d105290b937689e..170f0353904bef82095146c50bf98d5721ef2138 100644 --- a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.h +++ b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.h @@ -21,6 +21,7 @@ #include "KalmanFitterTool.h" #include <boost/dynamic_bitset.hpp> #include "CreateTrkTrackTool.h" +#include "xAODEventInfo/EventInfo.h" using ConstTrackStateProxy = Acts::detail_lt::TrackStateProxy<IndexSourceLink, 6, true>; using ClusterSet = boost::dynamic_bitset<>; @@ -148,6 +149,7 @@ private: ToolHandle<CreateTrkTrackTool> m_createTrkTrackTool {this, "CreateTrkTrackTool", "CreateTrkTrackTool"}; Gaudi::Property<bool> m_isMC {this, "isMC", false}; SG::WriteHandleKey<TrackCollection> m_trackCollection { this, "OutputCollection", "CKFTrackCollection", "Output track collection name" }; + SG::WriteDecorHandleKey<xAOD::EventInfo> m_eventInfoKey{this, "EventInfoKey", "EventInfo"}; }; #endif // FASERACTSKALMANFILTER_CKF2_H