From 40aedf46c5d50d763a2f618ffad5bd6a0fa8c095 Mon Sep 17 00:00:00 2001 From: Federico Meloni Date: Fri, 14 Nov 2025 23:32:49 +0100 Subject: [PATCH] added cut on outliers --- source/Refitting/include/RefitFinal.h | 1 + source/Refitting/src/RefitFinal.cc | 59 +++++++++------------------ source/Utils/include/FilterTracks.h | 9 ++-- source/Utils/src/FilterTracks.cc | 47 +++++++++++++++++---- 4 files changed, 66 insertions(+), 50 deletions(-) diff --git a/source/Refitting/include/RefitFinal.h b/source/Refitting/include/RefitFinal.h index 1a6bffd..a2ed3ee 100644 --- a/source/Refitting/include/RefitFinal.h +++ b/source/Refitting/include/RefitFinal.h @@ -91,6 +91,7 @@ class RefitFinal : public marlin::Processor { bool _extrapolateForward = true; int _minClustersOnTrackAfterFit = 0; + int _maxOutliersAllowed = 99; float _initialTrackError_d0{}; float _initialTrackError_phi0{}; diff --git a/source/Refitting/src/RefitFinal.cc b/source/Refitting/src/RefitFinal.cc index e4ec0f3..e40654b 100644 --- a/source/Refitting/src/RefitFinal.cc +++ b/source/Refitting/src/RefitFinal.cc @@ -8,6 +8,8 @@ #include #include +#include + #include #include #include @@ -85,6 +87,9 @@ RefitFinal::RefitFinal() : Processor("RefitFinal") { registerProcessorParameter("MinClustersOnTrackAfterFit", "Final minimum number of track clusters", _minClustersOnTrackAfterFit, int(4)); + registerProcessorParameter("MaxOutliersAllowed", "Maximum number of outliers allowed on the refitted track", + _maxOutliersAllowed, int(99)); + registerProcessorParameter("ReducedChi2Cut", "Cut on maximum allowed reduced chi2", _ReducedChi2Cut, double(-1.0)); } @@ -177,24 +182,6 @@ void RefitFinal::processEvent(LCEvent* evt) { ++hitInSubDet[_encoder->operator[](UTIL::LCTrackerCellID::subdet())]; } - //int init_status = FitInit2(track, marlin_trk.get()); - - //if (init_status != 0) { - // continue; - //} - - //streamlog_out(DEBUG4) << "Refit: Trackstate after initialisation\n" << marlin_trk->toString() << std::endl; - - //streamlog_out(DEBUG5) << "track initialised " << std::endl; - - //int fit_status = marlin_trk->fit(); - - //streamlog_out(DEBUG4) << "RefitHit: Trackstate after fit()\n" << marlin_trk->toString() << std::endl; - - //if (fit_status != 0) { - // continue; - //} - auto lcio_trk = std::unique_ptr(new IMPL::TrackImpl()); // setup initial dummy covariance matrix @@ -211,6 +198,13 @@ void RefitFinal::processEvent(LCEvent* evt) { covMatrix[9] = (_initialTrackError_z0); // sigma_z0^2 covMatrix[14] = (_initialTrackError_tanL); // sigma_tanl^2 + // get magnetic field at the origin + dd4hep::Detector& lcdd = dd4hep::Detector::getInstance(); + const double position[3] = {0, 0, 0}; // position to calculate magnetic field at (the origin in this case) + double magneticFieldVector[3] = {0, 0, 0}; // initialise object to hold magnetic field + lcdd.field().magneticField(position, magneticFieldVector); // get the magnetic field vector from DD4hep + _bField = magneticFieldVector[2] / dd4hep::tesla; + const bool fit_direction = _extrapolateForward ? MarlinTrk::IMarlinTrack::forward : MarlinTrk::IMarlinTrack::backward; int return_code = MarlinTrk::createFinalisedLCIOTrack(marlin_trk.get(), trkHits, lcio_trk.get(), fit_direction, covMatrix, _bField, _Max_Chi2_Incr); @@ -242,6 +236,14 @@ void RefitFinal::processEvent(LCEvent* evt) { marlin_trk->getOutliers(outliers); + if (int(outliers.size()) > _maxOutliersAllowed) { + streamlog_out(DEBUG3) << "More than " << _maxOutliersAllowed + << " outliers: Track " + "Discarded. Number of outliers = " + << outliers.size() << std::endl; + continue; + } + std::vector all_hits; all_hits.reserve(hits_in_fit.size() + outliers.size()); @@ -316,24 +318,3 @@ LCCollection* RefitFinal::GetCollection(LCEvent* evt, std::string colName) { return col; } - -int RefitFinal::FitInit2(Track* track, MarlinTrk::IMarlinTrack* marlinTrk) { - TrackStateImpl trackState; - - if (_refPoint == -1) { - trackState = TrackStateImpl(TrackState::AtOther, track->getD0(), track->getPhi(), track->getOmega(), track->getZ0(), - track->getTanLambda(), track->getCovMatrix(), track->getReferencePoint()); - } else { - const TrackState* trackAtHit = track->getTrackState(_refPoint); - if (not trackAtHit) { - streamlog_out(ERROR) << "Cannot find trackstate for " << _refPoint << std::endl; - return MarlinTrk::IMarlinTrack::error; - } - trackState = TrackStateImpl(*trackAtHit); - } - - const bool direction = _extrapolateForward ? MarlinTrk::IMarlinTrack::forward : MarlinTrk::IMarlinTrack::backward; - marlinTrk->initialise(trackState, _bField, direction); - - return MarlinTrk::IMarlinTrack::success; -} diff --git a/source/Utils/include/FilterTracks.h b/source/Utils/include/FilterTracks.h index 67a995c..636eb73 100644 --- a/source/Utils/include/FilterTracks.h +++ b/source/Utils/include/FilterTracks.h @@ -52,9 +52,11 @@ class FilterTracks : public marlin::Processor { private: //! Input track collection std::string _InputTrackCollection{}; + std::string _InputTrackRelationCollection{}; //! Output track collection std::string _OutputTrackCollection{}; + std::string _OutputTrackRelationCollection{}; bool _BarrelOnly = false; bool _HasCaloState = false; @@ -71,15 +73,16 @@ class FilterTracks : public marlin::Processor { //! Cut off for maximum number of holes on track int _MaxHoles = 0; - //! Cut off for maximum number of outliers on track - int _MaxOutliers = 20; - //! Cut off for momentum (GeV) float _MinPt = 1.0; // units GeV //! Cut off for spatial and temporal chi squared values float _Chi2Spatial = 0; + //! Cut off for impact parameters [mm] + float _MaxD0 = 999.; + float _MaxZ0 = 999.; + //! Default magnetic field value (Tesla) float _Bz = 5.0; }; diff --git a/source/Utils/src/FilterTracks.cc b/source/Utils/src/FilterTracks.cc index 3db287a..7aba9f9 100644 --- a/source/Utils/src/FilterTracks.cc +++ b/source/Utils/src/FilterTracks.cc @@ -9,6 +9,10 @@ #include #include #include +#include +#include + +#include FilterTracks aFilterTracks; @@ -18,6 +22,20 @@ FilterTracks::FilterTracks() : Processor("FilterTracks") { "filtered collection"; // register steering parameters: name, description, class-variable, default value + registerInputCollection(LCIO::TRACK, "InputTrackCollectionName", "Name of the input collection", + _InputTrackCollection, _InputTrackCollection); + + registerOutputCollection(LCIO::TRACK, "OutputTrackCollectionName", "Name of output collection", + _OutputTrackCollection, std::string("FilteredTracks")); + + registerInputCollection(LCIO::LCRELATION, "InputRelationCollectionName", + "Name of the input track to MCParticle relation collection", _InputTrackRelationCollection, + _InputTrackRelationCollection); + + registerOutputCollection(LCIO::LCRELATION, "OutputRelationCollectionName", + "Refit Track to MCParticle relation collection Name", _OutputTrackRelationCollection, + _OutputTrackRelationCollection); + registerProcessorParameter("BarrelOnly", "If true, just keep tracks with only barrel hits", _BarrelOnly, _BarrelOnly); registerProcessorParameter("HasCaloState", @@ -34,17 +52,13 @@ FilterTracks::FilterTracks() : Processor("FilterTracks") { registerProcessorParameter("MaxHoles", "Maximum number of holes on track", _MaxHoles, _MaxHoles); - registerProcessorParameter("MaxOutliers", "Max number of outliers", _MaxOutliers, _MaxOutliers); - registerProcessorParameter("MinPt", "Minimum transverse momentum", _MinPt, _MinPt); - registerProcessorParameter("Chi2Spatial", "Spatial chi squared", _Chi2Spatial, _Chi2Spatial); + registerProcessorParameter("MaxD0", "Max d0", _MaxD0, _MaxD0); - registerInputCollection(LCIO::TRACK, "InputTrackCollectionName", "Name of the input collection", - _InputTrackCollection, _InputTrackCollection); + registerProcessorParameter("MaxZ0", "Max z0", _MaxZ0, _MaxZ0); - registerOutputCollection(LCIO::TRACK, "OutputTrackCollectionName", "Name of output collection", - _OutputTrackCollection, std::string("FilteredTracks")); + registerProcessorParameter("Chi2Spatial", "Spatial chi squared", _Chi2Spatial, _Chi2Spatial); } void FilterTracks::init() { @@ -70,6 +84,9 @@ void FilterTracks::processEvent(LCEvent* evt) { LCCollectionVec* OutputTrackCollection = new LCCollectionVec(LCIO::TRACK); OutputTrackCollection->setSubset(true); + // track-mcpart relation output collections + UTIL::LCRelationNavigator myNav = UTIL::LCRelationNavigator(LCIO::TRACK, LCIO::MCPARTICLE); + // Get input collection LCCollection* InputTrackCollection = evt->getCollection(_InputTrackCollection); @@ -77,6 +94,10 @@ void FilterTracks::processEvent(LCEvent* evt) { throw EVENT::Exception("Invalid collection type: " + InputTrackCollection->getTypeName()); } + // Get input relations + LCCollection* InputTrackRelationCollection = evt->getCollection(_InputTrackRelationCollection); + std::shared_ptr relation = std::make_shared(InputTrackRelationCollection); + // Filter std::string encoderString = lcio::LCTrackerCellID::encoding_string(); UTIL::CellIDDecoder decoder(encoderString); @@ -122,14 +143,24 @@ void FilterTracks::processEvent(LCEvent* evt) { } } else { // track property cuts if (nhittotal > _NHitsTotal && nhitvertex > _NHitsVertex && nhitinner > _NHitsInner && nhitouter > _NHitsOuter && - pt > _MinPt && chi2spatial > _Chi2Spatial && nholes <= _MaxHoles) { + pt > _MinPt && chi2spatial > _Chi2Spatial && nholes <= _MaxHoles && + fabs(trk->getD0()) < _MaxD0 && fabs(trk->getZ0()) < _MaxZ0) { OutputTrackCollection->addElement(trk); + + auto mcParticleVec = relation->getRelatedToObjects(trk); + auto weightVec = relation->getRelatedToWeights(trk); + for (size_t irel = 0; irel < mcParticleVec.size(); ++irel) { + myNav.addRelation(trk, mcParticleVec[irel], weightVec[irel]); + } + } } } // Save output track collection evt->addCollection(OutputTrackCollection, _OutputTrackCollection); + LCCollection* outputTrackHitRel = myNav.createLCCollection(); + evt->addCollection(outputTrackHitRel, _OutputTrackRelationCollection); } void FilterTracks::end() {}