From ea909bf10e36cf785a2582d9a5982d0d3482f882 Mon Sep 17 00:00:00 2001 From: Chengdong Fu Date: Sun, 31 Mar 2024 22:04:10 +0800 Subject: [PATCH 01/10] fix bug --- Reconstruction/SiliconTracking/src/SiliconTrackingAlg.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.cpp b/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.cpp index 1ef104fa4..37110a924 100644 --- a/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.cpp +++ b/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.cpp @@ -1674,7 +1674,8 @@ void SiliconTrackingAlg::CreateTrack(TrackExtended * trackAR ) { for (int itrk=0; itrkgetTrackerHitExtendedVec(); + // fucd: TrackerHitExtendedVec& will change after merge split tracks, so must TrackerHitExtendedVec + TrackerHitExtendedVec hitVecOld = trackOld->getTrackerHitExtendedVec(); float phiNew = trackAR->getPhi(); float phiOld = trackOld->getPhi(); From 767b1bf1b22e03a5214d23d05c4f97814ac2dbf8 Mon Sep 17 00:00:00 2001 From: Chengdong Fu Date: Mon, 1 Apr 2024 08:09:42 +0800 Subject: [PATCH 02/10] import as assocaition maker between MCPartilce and Track --- .../RecAssociationMaker/CMakeLists.txt | 20 +++ .../src/TrackParticleRelationAlg.cpp | 150 ++++++++++++++++++ .../src/TrackParticleRelationAlg.h | 38 +++++ 3 files changed, 208 insertions(+) create mode 100644 Reconstruction/RecAssociationMaker/CMakeLists.txt create mode 100644 Reconstruction/RecAssociationMaker/src/TrackParticleRelationAlg.cpp create mode 100644 Reconstruction/RecAssociationMaker/src/TrackParticleRelationAlg.h diff --git a/Reconstruction/RecAssociationMaker/CMakeLists.txt b/Reconstruction/RecAssociationMaker/CMakeLists.txt new file mode 100644 index 000000000..d854356b8 --- /dev/null +++ b/Reconstruction/RecAssociationMaker/CMakeLists.txt @@ -0,0 +1,20 @@ +# Modules +gaudi_add_module(RecAssociationMaker + SOURCES src/TrackParticleRelationAlg.cpp + + LINK GearSvc + EventSeeder + TrackSystemSvcLib + DataHelperLib + KiTrackLib + Gaudi::GaudiKernel + k4FWCore::k4FWCore + ${GEAR_LIBRARIES} + ${GSL_LIBRARIES} + ${LCIO_LIBRARIES} +) +install(TARGETS RecAssociationMaker + EXPORT CEPCSWTargets + RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" COMPONENT bin + LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib + COMPONENT dev) diff --git a/Reconstruction/RecAssociationMaker/src/TrackParticleRelationAlg.cpp b/Reconstruction/RecAssociationMaker/src/TrackParticleRelationAlg.cpp new file mode 100644 index 000000000..1c764257d --- /dev/null +++ b/Reconstruction/RecAssociationMaker/src/TrackParticleRelationAlg.cpp @@ -0,0 +1,150 @@ +#include "TrackParticleRelationAlg.h" + +DECLARE_COMPONENT(TrackParticleRelationAlg) + +TrackParticleRelationAlg::TrackParticleRelationAlg(const std::string& name, ISvcLocator* svcLoc) +: GaudiAlgorithm(name, svcLoc){ + declareProperty("MCParticleCollection", m_inMCParticleColHdl, "Handle of the Input MCParticle collection"); +} + +StatusCode TrackParticleRelationAlg::initialize() { + for(auto name : m_inTrackCollectionNames) { + m_inTrackColHdls.push_back(new DataHandle (name, Gaudi::DataHandle::Reader, this)); + m_outColHdls.push_back(new DataHandle (name+"ParticleAssociation", Gaudi::DataHandle::Writer, this)); + } + + for(unsigned i=0; i (m_inAssociationCollectionNames[i], Gaudi::DataHandle::Reader, this)); + } + + if(m_inAssociationColHdls.size()==0) { + warning() << "no Association Collection input" << endmsg; + return StatusCode::FAILURE; + } + + m_nEvt = 0; + return GaudiAlgorithm::initialize(); +} + +StatusCode TrackParticleRelationAlg::execute() { + info() << "start Event " << m_nEvt << endmsg; + + // MCParticle + const edm4hep::MCParticleCollection* mcpCol = nullptr; + try { + mcpCol = m_inMCParticleColHdl.get(); + for (auto mcp : *mcpCol) { + debug() << "id=" << mcp.id() << " pdgid=" << mcp.getPDG() << " charge=" << mcp.getCharge() << " genstat=" << mcp.getGeneratorStatus() << " simstat=" << mcp.getSimulatorStatus() + << " p=" << mcp.getMomentum() << endmsg; + int nparents = mcp.parents_size(); + int ndaughters = mcp.daughters_size(); + for (int i=0; i mapHitParticle; + debug() << "reading Association" << endmsg; + for (auto hdl : m_inAssociationColHdls) { + const edm4hep::MCRecoTrackerAssociationCollection* assCol = nullptr; + try { + assCol = hdl->get(); + } + catch ( GaudiException &e ) { + debug() << "Collection " << hdl->fullKey() << " is unavailable in event " << m_nEvt << endmsg; + return StatusCode::FAILURE; + } + for (auto ass: *assCol) { + auto hit = ass.getRec(); + auto particle = ass.getSim().getMCParticle(); + mapHitParticle[hit] = particle; + } + } + + // Handle all input TrackCollection + for (unsigned icol=0; icolget(); + } + catch ( GaudiException &e ) { + debug() << "Collection " << hdl->fullKey() << " is unavailable in event " << m_nEvt << endmsg; + } + + auto outCol = m_outColHdls[icol]->createAndPut(); + if(!outCol) { + error() << "Collection " << m_outColHdls[icol]->fullKey() << " cannot be created" << endmsg; + return StatusCode::FAILURE; + } + + if(trkCol) { + std::map > mapParticleHits; + + for (auto track: *trkCol) { + std::map mapParticleNHits; + + // Count hits related to MCParticle + int nhits = track.trackerHits_size(); + for (int ihit=0; ihitsecond; + if(!particle.isAvailable()) continue; + debug() << "track " << track.id() << "'s hit " << hit.id() << " link to MCParticle " << particle.id() << endmsg; + auto itPN = mapParticleNHits.find(particle); + if (itPN!=mapParticleNHits.end()) itPN->second++; + else mapParticleNHits[particle] = 1; + + if (msgLevel(MSG::DEBUG)) mapParticleHits[particle].push_back(hit); + } + } + debug() << "Total " << mapParticleNHits.size() << " particles link to the track " << track.id() << endmsg; + + // Save to collection + for (std::map::iterator it=mapParticleNHits.begin(); it!=mapParticleNHits.end(); it++) { + auto ass = outCol->create(); + ass.setWeight(it->second); + ass.setRec(track); + ass.setSim(it->first); + } + } + + if (msgLevel(MSG::DEBUG)) { + for (std::map >::iterator it=mapParticleHits.begin(); it!=mapParticleHits.end(); it++) { + auto particle = it->first; + auto hits = it->second; + debug() << "=== MCPaticle ===" << particle << endmsg; + for (auto hit : hits) { + debug() << hit << endmsg; + } + } + } + } + } + + m_nEvt++; + return StatusCode::SUCCESS; +} + +StatusCode TrackParticleRelationAlg::finalize() { + + return StatusCode::SUCCESS; +} diff --git a/Reconstruction/RecAssociationMaker/src/TrackParticleRelationAlg.h b/Reconstruction/RecAssociationMaker/src/TrackParticleRelationAlg.h new file mode 100644 index 000000000..9f19212bd --- /dev/null +++ b/Reconstruction/RecAssociationMaker/src/TrackParticleRelationAlg.h @@ -0,0 +1,38 @@ +#ifndef TrackParticleRelationAlg_h +#define TrackParticleRelationAlg_h 1 + +#include "k4FWCore/DataHandle.h" +#include "GaudiAlg/GaudiAlgorithm.h" + +#include "edm4hep/MCParticleCollection.h" +#include "edm4hep/TrackCollection.h" +#include "edm4hep/MCRecoTrackerAssociationCollection.h" +#include "edm4hep/MCRecoTrackParticleAssociationCollection.h" + +class TrackParticleRelationAlg : public GaudiAlgorithm { + public: + // Constructor of this form must be provided + TrackParticleRelationAlg( const std::string& name, ISvcLocator* pSvcLocator ); + + // Three mandatory member functions of any algorithm + StatusCode initialize() override; + StatusCode execute() override; + StatusCode finalize() override; + + private: + // input MCParticle + DataHandle m_inMCParticleColHdl{"MCParticle", Gaudi::DataHandle::Reader, this}; + // input Tracks to make relation + std::vector* > m_inTrackColHdls; + Gaudi::Property > m_inTrackCollectionNames{this, "TrackList", {"SiTracks"}}; + // input TrackerAssociation to link TrackerHit and SimTrackerHit + std::vector* > m_inAssociationColHdls; + Gaudi::Property > m_inAssociationCollectionNames{this, "TrackerAssociationList", {"VXDTrackerHitAssociation", + "SITTrackerHitAssociation", "SETTrackerHitAssociation", "FTDTrackerHitAssociation"}}; + + // output TrackParticleAssociation + std::vector* > m_outColHdls; + + int m_nEvt; +}; +#endif From 0106792a112272b751409f8fefdecb52cee8ebd3 Mon Sep 17 00:00:00 2001 From: Chengdong Fu Date: Mon, 1 Apr 2024 08:16:22 +0800 Subject: [PATCH 03/10] import RecAssociationMaker --- Reconstruction/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/Reconstruction/CMakeLists.txt b/Reconstruction/CMakeLists.txt index 6dd63b02c..6bb53e061 100644 --- a/Reconstruction/CMakeLists.txt +++ b/Reconstruction/CMakeLists.txt @@ -4,3 +4,4 @@ add_subdirectory(PFA) add_subdirectory(SiliconTracking) add_subdirectory(Tracking) add_subdirectory(RecGenfitAlg) +add_subdirectory(RecAssociationMaker) From 1904ab2faacab2c400d242dc5ce58fb41e30421e Mon Sep 17 00:00:00 2001 From: Chengdong Fu Date: Mon, 1 Apr 2024 08:19:58 +0800 Subject: [PATCH 04/10] import fitter tool --- Reconstruction/Tracking/CMakeLists.txt | 3 + .../include/Tracking/ITrackFitterTool.h | 40 +++++++ .../Tracking/src/FitterTool/KalTestTool.cpp | 103 ++++++++++++++++++ .../Tracking/src/FitterTool/KalTestTool.h | 43 ++++++++ 4 files changed, 189 insertions(+) create mode 100644 Reconstruction/Tracking/include/Tracking/ITrackFitterTool.h create mode 100644 Reconstruction/Tracking/src/FitterTool/KalTestTool.cpp create mode 100644 Reconstruction/Tracking/src/FitterTool/KalTestTool.h diff --git a/Reconstruction/Tracking/CMakeLists.txt b/Reconstruction/Tracking/CMakeLists.txt index f3e77b472..9ea9ecd67 100644 --- a/Reconstruction/Tracking/CMakeLists.txt +++ b/Reconstruction/Tracking/CMakeLists.txt @@ -1,9 +1,12 @@ +gaudi_add_header_only_library(TrackingLib) + # Modules gaudi_add_module(Tracking SOURCES src/Clupatra/ClupatraAlg.cpp src/Clupatra/clupatra_new.cpp src/FullLDCTracking/FullLDCTrackingAlg.cpp src/TruthTracker/TruthTrackerAlg.cpp + src/FitterTool/KalTestTool.cpp LINK GearSvc EventSeeder TrackSystemSvcLib diff --git a/Reconstruction/Tracking/include/Tracking/ITrackFitterTool.h b/Reconstruction/Tracking/include/Tracking/ITrackFitterTool.h new file mode 100644 index 000000000..4de2b46e8 --- /dev/null +++ b/Reconstruction/Tracking/include/Tracking/ITrackFitterTool.h @@ -0,0 +1,40 @@ +#ifndef ITrackFitterTool_h +#define ITrackFitterTool_h + +/* + * Description: + * ITrackFitterTool is used to fit a track candidate to obtain track parameter + * + * The interface: + * * Fit: peform on tracker hits according prepared geometry + * + * Author: FU Chengdong + */ + +#include "GaudiKernel/IAlgTool.h" +#include "edm4hep/TrackState.h" +#include + + + +namespace edm4hep{ + class MutableTrack; + class TrackerHit; +} + +class ITrackFitterTool: virtual public IAlgTool { + public: + + DeclareInterfaceID(ITrackFitterTool, 0, 1); + virtual ~ITrackFitterTool() {} + + virtual int Fit(edm4hep::MutableTrack track, std::vector& trackHits, + const decltype(edm4hep::TrackState::covMatrix)& covMatrix, double maxChi2perHit, bool backward) = 0; + virtual int Fit(edm4hep::MutableTrack track, std::vector& trackHits, + edm4hep::TrackState trackState, double maxChi2perHit, bool backward) = 0; + virtual std::vector >& GetHitsInFit() = 0; + virtual std::vector >& GetOutliers() = 0; + virtual void Clear() = 0; +}; + +#endif diff --git a/Reconstruction/Tracking/src/FitterTool/KalTestTool.cpp b/Reconstruction/Tracking/src/FitterTool/KalTestTool.cpp new file mode 100644 index 000000000..773f1998a --- /dev/null +++ b/Reconstruction/Tracking/src/FitterTool/KalTestTool.cpp @@ -0,0 +1,103 @@ +#include "KalTestTool.h" + +#include "TrackSystemSvc/ITrackSystemSvc.h" +#include "TrackSystemSvc/MarlinTrkUtils.h" +#include "TrackSystemSvc/IMarlinTrack.h" +#include "DetInterface/IGeomSvc.h" + +#include "DD4hep/Detector.h" +#include "DD4hep/DD4hepUnits.h" + +#include "edm4hep/TrackerHit.h" +#include "edm4hep/TrackState.h" +#include "edm4hep/MutableTrack.h" + +DECLARE_COMPONENT(KalTestTool) + +StatusCode KalTestTool::initialize() { + StatusCode sc; + always() << m_fitterName << endmsg; + if (m_fitterName=="KalTest") { + auto _trackSystemSvc = service("TrackSystemSvc"); + if (!_trackSystemSvc) { + error() << "Failed to find TrackSystemSvc ..." << endmsg; + return StatusCode::FAILURE; + } + m_factoryMarlinTrk = _trackSystemSvc->getTrackSystem(this); + m_factoryMarlinTrk->setOption(MarlinTrk::IMarlinTrkSystem::CFG::useQMS, m_useQMS); + m_factoryMarlinTrk->setOption(MarlinTrk::IMarlinTrkSystem::CFG::usedEdx, m_usedEdx); + m_factoryMarlinTrk->setOption(MarlinTrk::IMarlinTrkSystem::CFG::useSmoothing, m_useSmoothing); + m_factoryMarlinTrk->init(); + } + else { + error() << "fitter " << m_fitterName << " has not been imported" << endmsg; + return StatusCode::FAILURE; + } + + auto _geomSvc = service("GeomSvc"); + if ( !_geomSvc ) { + error() << "Failed to find GeomSvc ..." << endmsg; + return StatusCode::FAILURE; + } + + if (m_magneticField.value()==0) { + const dd4hep::Direction& field = _geomSvc->lcdd()->field().magneticField(dd4hep::Position(0,0,0)); + double Bz = field.z()/dd4hep::tesla; + if (Bz==0) { + error() << "magnetic field = 0, KalmanFilter cannot run" << endmsg; + return StatusCode::FAILURE; + } + m_magneticField = Bz; + } + + return sc; +} + +StatusCode KalTestTool::finalize() { + StatusCode sc; + return sc; +} + +int KalTestTool::Fit(edm4hep::MutableTrack track, std::vector& trackHits, + const decltype(edm4hep::TrackState::covMatrix)& covMatrix, double maxChi2perHit, bool backward) { + if (m_hitsInFit.size()!=0 || m_outliers.size()!=0) { + error() << "Important! vector not clear, still store the data of last event!" << endmsg; + return 0; + } + + if (m_fitterName=="KalTest") { + debug() << "start..." << endmsg; + std::shared_ptr marlinTrack(m_factoryMarlinTrk->createTrack()); + debug() << "created MarlinKalTestTrack" << endmsg; + int status = MarlinTrk::createFinalisedLCIOTrack(marlinTrack.get(), trackHits, &track, backward, covMatrix, m_magneticField, maxChi2perHit); + + marlinTrack->getHitsInFit(m_hitsInFit); + marlinTrack->getOutliers(m_outliers); + return status; + } + + error() << "Don't support the Fitter " << m_fitterName << endmsg; + return 0; +} + +int KalTestTool::Fit(edm4hep::MutableTrack track, std::vector& trackHits, + edm4hep::TrackState trackState, double maxChi2perHit, bool backward) { + if (m_hitsInFit.size()!=0 || m_outliers.size()!=0) { + error() << "Important! vector not clear, still store the data of last event!" << endmsg; + return 0; + } + + if (m_fitterName=="KalTest") { + debug() << "start..." << endmsg; + std::shared_ptr marlinTrack(m_factoryMarlinTrk->createTrack()); + debug() << "created MarlinKalTestTrack" << endmsg; + int status = MarlinTrk::createFinalisedLCIOTrack(marlinTrack.get(), trackHits, &track, backward, &trackState, m_magneticField, maxChi2perHit); + + marlinTrack->getHitsInFit(m_hitsInFit); + marlinTrack->getOutliers(m_outliers); + return status; + } + + error() << "Don't support the Fitter " << m_fitterName << endmsg; + return 0; +} diff --git a/Reconstruction/Tracking/src/FitterTool/KalTestTool.h b/Reconstruction/Tracking/src/FitterTool/KalTestTool.h new file mode 100644 index 000000000..0cce58bc4 --- /dev/null +++ b/Reconstruction/Tracking/src/FitterTool/KalTestTool.h @@ -0,0 +1,43 @@ +#ifndef KalTestTool_h +#define KalTestTool_h + +#include "GaudiKernel/AlgTool.h" +#include "Tracking/ITrackFitterTool.h" +#include + +namespace MarlinTrk { + class IMarlinTrkSystem ; +} + +class KalTestTool : public extends { + public: + using extends::extends; + //KalTestTool(void* p) { m_pAlgUsing=p; }; + + virtual int Fit(edm4hep::MutableTrack track, std::vector& trackHits, + const decltype(edm4hep::TrackState::covMatrix)& covMatrix, double maxChi2perHit, bool backward = true) override; + virtual int Fit(edm4hep::MutableTrack track, std::vector& trackHits, + edm4hep::TrackState trackState, double maxChi2perHit, bool backward = true) override; + + StatusCode initialize() override; + StatusCode finalize() override; + + std::vector >& GetHitsInFit() override {return m_hitsInFit;}; + std::vector >& GetOutliers() override {return m_outliers;}; + void Clear() override {m_hitsInFit.clear(); m_outliers.clear();}; + + private: + Gaudi::Property m_fitterName{this, "Fitter", "KalTest"}; + Gaudi::Property m_useQMS{this, "MSOn", true}; + Gaudi::Property m_usedEdx{this, "EnergyLossOn", true}; + Gaudi::Property m_useSmoothing{this, "Smooth", true}; + Gaudi::Property m_magneticField{this, "MagneticField", 0}; + + MarlinTrk::IMarlinTrkSystem* m_factoryMarlinTrk = nullptr; + + void* m_pAlgUsing = nullptr; + std::vector > m_hitsInFit ; + std::vector > m_outliers ; +}; + +#endif From a8ca1ffae80c79de0089b594bd3e5921b527059c Mon Sep 17 00:00:00 2001 From: Chengdong Fu Date: Mon, 1 Apr 2024 14:15:56 +0800 Subject: [PATCH 05/10] change to ITrackFitterTool --- Reconstruction/SiliconTracking/CMakeLists.txt | 1 + .../src/ForwardTrackingAlg.cpp | 17 +- .../SiliconTracking/src/ForwardTrackingAlg.h | 4 + .../src/SiliconTrackingAlg.cpp | 237 ++++++++++++------ .../SiliconTracking/src/SiliconTrackingAlg.h | 16 +- .../SiliconTracking/src/TrackSubsetAlg.cpp | 29 ++- .../SiliconTracking/src/TrackSubsetAlg.h | 4 + .../Tracking/src/Clupatra/ClupatraAlg.cpp | 3 + .../Tracking/src/Clupatra/ClupatraAlg.h | 3 + .../FullLDCTracking/FullLDCTrackingAlg.cpp | 37 ++- .../src/FullLDCTracking/FullLDCTrackingAlg.h | 3 + Service/GearSvc/src/GearSvc.h | 3 +- 12 files changed, 256 insertions(+), 101 deletions(-) mode change 100644 => 100755 Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp diff --git a/Reconstruction/SiliconTracking/CMakeLists.txt b/Reconstruction/SiliconTracking/CMakeLists.txt index 9eb873504..8fd9be94f 100644 --- a/Reconstruction/SiliconTracking/CMakeLists.txt +++ b/Reconstruction/SiliconTracking/CMakeLists.txt @@ -7,6 +7,7 @@ gaudi_add_module(SiliconTracking src/TrackSubsetAlg.cpp LINK GearSvc EventSeeder + TrackingLib TrackSystemSvcLib DataHelperLib KiTrackLib diff --git a/Reconstruction/SiliconTracking/src/ForwardTrackingAlg.cpp b/Reconstruction/SiliconTracking/src/ForwardTrackingAlg.cpp index 3eb06e297..3184563ff 100644 --- a/Reconstruction/SiliconTracking/src/ForwardTrackingAlg.cpp +++ b/Reconstruction/SiliconTracking/src/ForwardTrackingAlg.cpp @@ -4,7 +4,7 @@ #include "DataHelper/Navigation.h" #include "edm4hep/TrackerHit.h" -#include "edm4hep/TrackerHit.h" +#include "edm4hep/TrackerHitPlane.h" #include "edm4hep/Track.h" #if __has_include("edm4hep/EDM4hepVersion.h") #include "edm4hep/EDM4hepVersion.h" @@ -74,12 +74,16 @@ StatusCode ForwardTrackingAlg::initialize(){ // can be printed. As this is mainly used for debugging it is not a steerable parameter. //if( _useCED )MarlinCED::init(this) ; //CED + // Set up the track fit tool + m_fitTool = ToolHandle(m_fitToolName.value()); + if(m_dumpTime){ NTuplePtr nt1(ntupleSvc(), "MyTuples/Time"+name()); if ( !nt1 ) { m_tuple = ntupleSvc()->book("MyTuples/Time"+name(),CLID_ColumnWiseTuple,"Tracking time"); if ( 0 != m_tuple ) { m_tuple->addItem ("timeTotal", m_timeTotal ).ignore(); + m_tuple->addItem ("timeKalman", m_timeKalman ).ignore(); } else { fatal() << "Cannot book MyTuples/Time"+name() <getSector() ].push_back( ftdHit ); } } - + + if(m_dumpTime&&m_tuple) m_timeKalman = 0; + if( !_map_sector_hits.empty() ){ /**********************************************************************************************/ /* Check if no sector is overflowing with hits */ @@ -661,7 +667,7 @@ StatusCode ForwardTrackingAlg::execute(){ debug() << "\t\t---Save Tracks---" << endmsg ; //auto trkCol = _outColHdl.createAndPut(); - + for (unsigned int i=0; i < tracks.size(); i++){ FTDTrack* myTrack = dynamic_cast< FTDTrack* >( tracks[i] ); @@ -942,7 +948,7 @@ bool ForwardTrackingAlg::setCriteria( unsigned round ){ } void ForwardTrackingAlg::finaliseTrack( edm4hep::MutableTrack* trackImpl ){ - + auto stopwatch = TStopwatch(); Fitter fitter( trackImpl , _trkSystem ); //trackImpl->trackStates().clear(); @@ -1024,6 +1030,9 @@ void ForwardTrackingAlg::finaliseTrack( edm4hep::MutableTrack* trackImpl ){ trackImpl->addToSubDetectorHitNumbers(hitNumbers[UTIL::ILDDetID::SET]); trackImpl->addToSubDetectorHitNumbers(hitNumbers[UTIL::ILDDetID::ETD]); #endif + + if(m_dumpTime&&m_tuple) m_timeKalman += stopwatch.RealTime()*1000; + return; } diff --git a/Reconstruction/SiliconTracking/src/ForwardTrackingAlg.h b/Reconstruction/SiliconTracking/src/ForwardTrackingAlg.h index bd7a5173b..fee28994b 100644 --- a/Reconstruction/SiliconTracking/src/ForwardTrackingAlg.h +++ b/Reconstruction/SiliconTracking/src/ForwardTrackingAlg.h @@ -14,6 +14,7 @@ #include "edm4hep/TrackerHit.h" #include "TrackSystemSvc/IMarlinTrkSystem.h" +#include "Tracking/ITrackFitterTool.h" //#include "gear/BField.h" #include "KiTrack/Segment.h" @@ -239,6 +240,7 @@ class ForwardTrackingAlg : public GaudiAlgorithm { NTuple::Tuple* m_tuple; NTuple::Item m_timeTotal; + NTuple::Item m_timeKalman; /** B field in z direction */ double _Bz; @@ -261,6 +263,7 @@ class ForwardTrackingAlg : public GaudiAlgorithm { Gaudi::Property > _critMinimaInit{this, "CriteriaMin", {} }; Gaudi::Property > _critMaximaInit{this, "CriteriaMax", {} }; Gaudi::Property m_dumpTime{this, "DumpTime", true}; + Gaudi::Property m_fitToolName{this, "FitterTool", "KalTestTool/KalTest110"}; std::map > _critMinima; std::map > _critMaxima; @@ -291,6 +294,7 @@ class ForwardTrackingAlg : public GaudiAlgorithm { unsigned _nTrackCandidatesPlus; MarlinTrk::IMarlinTrkSystem* _trkSystem; + ToolHandle m_fitTool; /** The quality of the output track collection */ int _output_track_col_quality ; diff --git a/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.cpp b/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.cpp index 37110a924..c10689121 100644 --- a/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.cpp +++ b/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.cpp @@ -113,13 +113,17 @@ StatusCode SiliconTrackingAlg::initialize() { _nRun = -1 ; _nEvt = 0 ; //printParameters() ; - if(m_dumpTime){ + + // Set up the track fit tool + m_fitTool = ToolHandle(m_fitToolName.value()); + + if (m_dumpTime) { NTuplePtr nt1(ntupleSvc(), "MyTuples/Time"+name()); if ( !nt1 ) { m_tuple = ntupleSvc()->book("MyTuples/Time"+name(),CLID_ColumnWiseTuple,"Tracking time"); if ( 0 != m_tuple ) { - m_tuple->addItem ("timeTotal", m_timeTotal ).ignore(); - m_tuple->addItem ("timeKalman", m_timeKalman ).ignore(); + m_tuple->addItem("timeTotal", m_timeTotal).ignore(); + m_tuple->addItem("timeKalman", m_timeKalman).ignore(); } else { fatal() << "Cannot book MyTuples/Time"+name() <book("MyTuples/"+name(),CLID_ColumnWiseTuple,"Tracking time"); + if ( 0 != m_tuple ) { + m_tupleDebug->addItem("nhit", m_nhit, 0, 100).ignore(); + m_tupleDebug->addIndexedItem("layer", m_nhit, m_layer).ignore(); + m_tupleDebug->addIndexedItem("chi2", m_nhit, m_chi2).ignore(); + m_tupleDebug->addIndexedItem("used", m_nhit, m_used).ignore(); + } + else { + fatal() << "Cannot book MyTuples/"+name() <("TrackSystemSvc"); @@ -192,7 +217,7 @@ StatusCode SiliconTrackingAlg::execute(){ // zero triplet counters _ntriplets = _ntriplets_good = _ntriplets_2MCP = _ntriplets_3MCP = _ntriplets_1MCP_Bad = _ntriplets_bad = 0; - + // Clearing the working containers from the previous event // FIXME: partly done at the end of the event, in CleanUp. Make it consistent. //_tracksWithNHitsContainer.clear(); @@ -205,14 +230,14 @@ StatusCode SiliconTrackingAlg::execute(){ //int runNo = header.getRunNumber(); //debug() << "Processing Run[" << runNo << "]::Event[" << evtNo << "]" << endmsg; - _trackImplVec.reserve(100); + //_trackImplVec.reserve(100); int successVTX = InitialiseVTX(); //int successFTD = 0; int successFTD = InitialiseFTD(); if (successVTX == 1) { - debug() << " phi theta layer nh o : m : i :: o*m*i " << endmsg; + debug() << " --phi-- --theta-- --layer-- nh o : m : i :: o*m*i " << endmsg; for (int iPhi=0; iPhi<_nDivisionsInPhi; ++iPhi) { for (int iTheta=0; iTheta<_nDivisionsInTheta;++iTheta) { @@ -225,7 +250,7 @@ StatusCode SiliconTrackingAlg::execute(){ } if (successFTD == 1) { - debug() << " phi side layer nh o : m : i :: o*m*i " << endmsg; + debug() << " --phi-- --side-- --layer-- nh o : m : i :: o*m*i " << endmsg; TrackingInFTD(); // Perform tracking in the FTD debug() << "End of Processing FTD sectors" << endmsg; } @@ -250,7 +275,17 @@ StatusCode SiliconTrackingAlg::execute(){ trackIter < tracksWithNHits.end(); trackIter++) { CreateTrack( *trackIter ); } - debug() << "End of creating "<< nHits << " hits tracks " << endmsg; + debug() << "End of creating "<< nHits << " hits tracks, total " << _trackImplVec.size() << " tracks now" << endmsg; + + if (msgLevel(MSG::DEBUG)) { + for (auto trackAR : _trackImplVec) { + debug() << "track " << trackAR << " has hits:" << endmsg; + TrackerHitExtendedVec& hitVec = trackAR->getTrackerHitExtendedVec(); + for (auto hit : hitVec) { + debug() << " " << hit->getTrackerHit().id().collectionID << "-" << hit->getTrackerHit().id().index << endmsg; + } + } + } } if (_attachFast == 0) { @@ -539,7 +574,7 @@ int SiliconTrackingAlg::InitialiseFTD() { int iCode = iSemiSphere + 2*layer + 2*_nlayersFTD*iPhi; _sectorsFTD[iCode].push_back( hitExt ); - debug() << " FTD Pixel Hit added : @ " << pos[0] << " " << pos[1] << " " << pos[2] << " drphi " << hitExt->getResolutionRPhi() << " dz " << hitExt->getResolutionZ() << " iPhi = " << iPhi << " iSemiSphere " << iSemiSphere << " iCode = " << iCode << " layer = " << layer << " cellID = " << hit.getCellID() << endmsg; + debug() << " FTD Pixel Hit " << hit.id().index << " added : @ " << pos[0] << " " << pos[1] << " " << pos[2] << " drphi " << hitExt->getResolutionRPhi() << " dz " << hitExt->getResolutionZ() << " iPhi = " << iPhi << " iSemiSphere " << iSemiSphere << " iCode = " << iCode << " layer = " << layer << " cellID = " << hit.getCellID() << endmsg; } } // Reading out FTD SpacePoint Collection @@ -603,7 +638,7 @@ int SiliconTrackingAlg::InitialiseFTD() { double Phi = atan2(pos[1],pos[0]); if (Phi < 0.) Phi = Phi + TWOPI; - + // get the layer number unsigned int layer = static_cast(getLayerID(hit)); unsigned int petalIndex = static_cast(getModuleID(hit)); @@ -636,7 +671,7 @@ int SiliconTrackingAlg::InitialiseFTD() { int iCode = iSemiSphere + 2*layer + 2*_nlayersFTD*iPhi; _sectorsFTD[iCode].push_back( hitExt ); - debug() << " FTD SpacePoint Hit added : @ " << pos[0] << " " << pos[1] << " " << pos[2] << " drphi " << hitExt->getResolutionRPhi() << " dz " << hitExt->getResolutionZ() << " iPhi = " << iPhi << " iSemiSphere " << iSemiSphere << " iCode = " << iCode << " layer = " << layer << " cellID = " << hit.getCellID() << endmsg; + debug() << " FTD SpacePoint Hit added : @ " << pos[0] << " " << pos[1] << " " << pos[2] << " drphi " << hitExt->getResolutionRPhi() << " dz " << hitExt->getResolutionZ() << " iPhi = " << iPhi << " iSemiSphere " << iSemiSphere << " iCode = " << iCode << " layer = " << layer << " cellID = " << hit.getCellID() << endmsg; } } @@ -713,7 +748,7 @@ int SiliconTrackingAlg::InitialiseVTX() { // SJA:FIXME: just use planar res for now hitExt->setResolutionRPhi(hit.getCovMatrix()[2]); hitExt->setResolutionZ(hit.getCovMatrix()[5]); - + // set type is now only used in one place where it is set to 0 to reject hits from a fit, set to INT_MAX to try and catch any missuse hitExt->setType(int(INT_MAX)); // det is no longer used set to INT_MAX to try and catch any missuse @@ -733,7 +768,7 @@ int SiliconTrackingAlg::InitialiseVTX() { double Phi = atan2(pos[1],pos[0]); if (Phi < 0.) Phi = Phi + TWOPI; - + // get the layer number int layer = getLayerID(hit); @@ -747,7 +782,7 @@ int SiliconTrackingAlg::InitialiseVTX() { int iCode = layer + _nLayers*iPhi + _nLayers*_nDivisionsInPhi*iTheta; _sectors[iCode].push_back( hitExt ); - debug() << " VXD Hit " << hit.id() << " added : @ " << pos[0] << " " << pos[1] << " " << pos[2] << " drphi " << hitExt->getResolutionRPhi() << " dz " << hitExt->getResolutionZ() << " iPhi = " << iPhi << " iTheta " << iTheta << " iCode = " << iCode << " layer = " << layer << endmsg; + debug() << " VXD Hit " << hit.id().index << " added : @ " << pos[0] << " " << pos[1] << " " << pos[2] << " drphi " << hitExt->getResolutionRPhi() << " dz " << hitExt->getResolutionZ() << " iPhi = " << iPhi << " iTheta " << iTheta << " iCode = " << iCode << " layer = " << layer << endmsg; } } @@ -906,7 +941,7 @@ int SiliconTrackingAlg::InitialiseVTX() { int iCode = layer + _nLayers*iPhi + _nLayers*_nDivisionsInPhi*iTheta; _sectors[iCode].push_back( hitExt ); - debug() << " SIT Hit " << trkhit.id() << " added : @ " << pos[0] << " " << pos[1] << " " << pos[2] << " drphi " << hitExt->getResolutionRPhi() << " dz " << hitExt->getResolutionZ() << " iPhi = " << iPhi << " iTheta " << iTheta << " iCode = " << iCode << " layer = " << layer << endmsg; + debug() << " SIT Hit " << trkhit.id().index << " added : @ " << pos[0] << " " << pos[1] << " " << pos[2] << " drphi " << hitExt->getResolutionRPhi() << " dz " << hitExt->getResolutionZ() << " iPhi = " << iPhi << " iTheta " << iTheta << " iCode = " << iCode << " layer = " << layer << endmsg; } } @@ -1063,13 +1098,13 @@ void SiliconTrackingAlg::ProcessOneSector(int iPhi, int iTheta) { if (nHitsInner > 0) { - debug() << " " - << std::setw(3) << iPhi << " " << std::setw(3) << ipMiddle << " " << std::setw(3) << ipInner << " " - << std::setw(3) << iTheta << " " << std::setw(3) << itMiddle << " " << std::setw(3) << itInner << " " - << std::setw(3) << nLR[0] << " " << std::setw(3) << nLR[1] << " " << std::setw(3) << nLR[2] << " " - << std::setw(3) << nHitsOuter << " : " << std::setw(3) << nHitsMiddle << " : " << std::setw(3) << nHitsInner << " :: " - << std::setw(3) << nHitsOuter*nHitsMiddle* nHitsInner << endmsg; - + debug() << "pattern " + << std::setw(3) << iPhi << " " << std::setw(3) << ipMiddle << " " << std::setw(3) << ipInner << " " + << std::setw(3) << iTheta << " " << std::setw(3) << itMiddle << " " << std::setw(3) << itInner << " " + << std::setw(3) << nLR[0] << " " << std::setw(3) << nLR[1] << " " << std::setw(3) << nLR[2] << " " + << std::setw(3) << nHitsOuter << " : " << std::setw(3) << nHitsMiddle << " : " << std::setw(3) << nHitsInner << " :: " + << std::setw(3) << nHitsOuter*nHitsMiddle*nHitsInner << endmsg; + // test all triplets for (int iOuter=0; iOutergetTrackerHit().id().index << " " << middleHit->getTrackerHit().id().index << " " << innerHit->getTrackerHit().id().index << " innerLayer=" << innerLayer << endmsg; for (int layer = innerLayer-1; layer>=0; layer--) { // loop over remaining layers float distMin = 1.0e+20; @@ -1571,7 +1609,7 @@ int SiliconTrackingAlg::BuildTrack(TrackerHitExtended * outerHit, float chi2RPhi; float chi2Z; - //debug() << "######## number of hits to fit with _fastfitter = " << NPT << endmsg; + //debug() << "######## number of hits to fit with _fastfitter = " << NPT << endmsg; _fastfitter->fastHelixFit(NPT, xh, yh, rh, ph, wrh, zh, wzh,iopt, par, epar, chi2RPhi, chi2Z); par[3] = par[3]*par[0]/fabs(par[0]); @@ -1623,7 +1661,8 @@ int SiliconTrackingAlg::BuildTrack(TrackerHitExtended * outerHit, } // endloop over remaining layers TrackerHitExtendedVec& hvec = trackAR->getTrackerHitExtendedVec(); int nTotalHits = int(hvec.size()); - debug() << "######## number of hits to return = " << nTotalHits << endmsg; + debug() << "######## number of hits to return = " << nTotalHits << endmsg; + return nTotalHits; } @@ -1654,15 +1693,20 @@ void SiliconTrackingAlg::CreateTrack(TrackExtended * trackAR ) { Method which creates Track out of TrackExtended objects. Checks for possible track splitting (separate track segments in VXD and FTD). */ - + debug() << "CreateTrack " << trackAR << endmsg; TrackerHitExtendedVec& hitVec = trackAR->getTrackerHitExtendedVec(); int nHits = int(hitVec.size()); for (int i=0; igetTrackExtendedVec(); - if (trackVec.size() != 0) + //edm4hep::TrackerHit trkHit = hitVec[i]->getTrackerHit(); + //debug() << "Hit " << trkHit.id().collectionID << "-" << trkHit.id().index << endmsg; + if (trackVec.size() != 0) { + //debug() << " used in TrackExtended " << endmsg; return ; + //continue; + } } // First check if the current track is piece of the split one @@ -1674,7 +1718,6 @@ void SiliconTrackingAlg::CreateTrack(TrackExtended * trackAR ) { for (int itrk=0; itrkgetTrackerHitExtendedVec(); float phiNew = trackAR->getPhi(); @@ -1714,6 +1757,7 @@ void SiliconTrackingAlg::CreateTrack(TrackExtended * trackAR ) { } for (int ih=0;ihgetTrackerHit(); + //debug() << "old track's hit " << trkHit.id().collectionID << " " << trkHit.id().index << endmsg; xh[ih+nHits] = trkHit.getPosition()[0]; yh[ih+nHits] = trkHit.getPosition()[1]; zh[ih+nHits] = float(trkHit.getPosition()[2]); @@ -1810,6 +1854,15 @@ void SiliconTrackingAlg::CreateTrack(TrackExtended * trackAR ) { // Attach hits belonging to the current track segment to // the track already created if (found == 1) { + if (iBad!=-1) { + if (iBad>=nHits) + debug() << "Bad hit " << iBad << " " << hitVecOld[iBad-nHits]->getTrackerHit().id().collectionID << "-" << hitVecOld[iBad-nHits]->getTrackerHit().id().index << endmsg; + else + debug() << "Bad hit " << iBad << " " << hitVec[iBad]->getTrackerHit().id().collectionID << "-" << hitVec[iBad]->getTrackerHit().id().index << endmsg; + } + else + debug() << "Bad hit " << iBad << endmsg; + trackOld->ClearTrackerHitExtendedVec(); for (int i=0;igetTrackerHit().id().collectionID << "-" << trkHit->getTrackerHit().id().index << endmsg; trkHit->clearTrackVec(); if (icur == iBad) { } @@ -1848,7 +1902,8 @@ void SiliconTrackingAlg::CreateTrack(TrackExtended * trackAR ) { trackOld->setChi2(chi2Min*float(ndf)); trackOld->setNDF(ndf); trackOld->setCovMatrix(eparmin); - // trackOld->setReferencePoint(refPointMin); + // TODO: by fucd, removed why in ILCSoft + //trackOld->setReferencePoint(refPointMin); } delete[] xh; @@ -2027,7 +2082,8 @@ void SiliconTrackingAlg::AttachRemainingVTXHitsFast() { } void SiliconTrackingAlg::AttachRemainingVTXHitsSlow() { - + debug() << "AttachRemainingVTXHitsSlow()" << endmsg; + TrackerHitExtendedVec nonAttachedHits; nonAttachedHits.clear(); @@ -2050,11 +2106,9 @@ void SiliconTrackingAlg::AttachRemainingVTXHitsSlow() { if(isize>maxTrackSize)maxTrackSize = isize; } if (maxTrackSize<=3) { - debug() << " Add non attached hit to list: id = " << hit->getTrackerHit().id() << endmsg; + debug() << " Add non attached hit to list: id = " << hit->getTrackerHit().id().index << endmsg; nonAttachedHits.push_back( hit ); } - - } } } @@ -2065,7 +2119,7 @@ void SiliconTrackingAlg::AttachRemainingVTXHitsSlow() { int nTrk = int(_trackImplVec.size()); for (int iHit=0; iHitgetZ0(); float omega = trackAR->getOmega(); float tanlambda = trackAR->getTanLambda(); + //debug() << "attach hit now, field = " << _bField << " d0 = " << d0 << " phi0 = " << phi0 << " z0 = " << z0 << " omega = " << omega << " tanlambda = " << tanlambda << endmsg; helix.Initialize_Canonical(phi0,d0,z0,omega,tanlambda,_bField); float distance[3]; float time = helix.getDistanceToPoint(pos,distance); if (time < 1.0e+10) { if (distance[2] < minDist) { - minDist = distance[2]; + minDist = distance[2]; trackToAttach = trackAR; } } @@ -2123,16 +2178,18 @@ void SiliconTrackingAlg::AttachRemainingVTXHitsSlow() { } if (minDist < _minDistCutAttach && trackToAttach != NULL) { int iopt = 2; - debug() << " Hit: id = " << hit->getTrackerHit().id() << " : try attachement"<< endmsg; + debug() << " Hit: id = " << hit->getTrackerHit().id().index << " : try attachement"<< endmsg; AttachHitToTrack(trackToAttach,hit,iopt); } else { - debug() << " Hit: id = " << hit->getTrackerHit().id() << " rejected due to distance cut of " <<_minDistCutAttach<< " min distance = " << minDist << endmsg; + debug() << " Hit: id = " << hit->getTrackerHit().id().index << " rejected due to distance cut of " <<_minDistCutAttach<< " min distance = " << minDist << endmsg; } } } } void SiliconTrackingAlg::AttachRemainingFTDHitsSlow() { + debug() << "AttachRemainingFTDHitsSlow()" << endmsg; + TrackerHitExtendedVec nonAttachedHits; nonAttachedHits.clear(); @@ -2146,6 +2203,7 @@ void SiliconTrackingAlg::AttachRemainingFTDHitsSlow() { TrackerHitExtended * hit = hitVec[iH]; TrackExtendedVec& trackVec = hit->getTrackExtendedVec(); if (trackVec.size()==0) + debug() << " Add non attached hit to list: id = " << hit->getTrackerHit().id().index << endmsg; nonAttachedHits.push_back( hit ); } } @@ -2173,7 +2231,7 @@ void SiliconTrackingAlg::AttachRemainingFTDHitsSlow() { // SJA:FIXME: check to see if allowing no hits in the same sensor vs no hits in the same layer works // if (hit->getTrackerHit()->getType() == hitVector[IHIT]->getTrackerHit()->getType()) { if (hit->getTrackerHit().getCellID() == hitVector[IHIT]->getTrackerHit().getCellID()) { - + debug() << " hit: id = " << hit->getTrackerHit().id().index << " condsidered delta together with hit " << hitVector[IHIT]->getTrackerHit().id().index << endmsg; consider = false; break; } @@ -2193,7 +2251,7 @@ void SiliconTrackingAlg::AttachRemainingFTDHitsSlow() { float time = helix.getDistanceToPoint(pos,distance); if (time < 1.0e+10) { if (distance[2] < minDist) { - minDist = distance[2]; + minDist = distance[2]; trackToAttach = trackAR; } } @@ -2202,8 +2260,12 @@ void SiliconTrackingAlg::AttachRemainingFTDHitsSlow() { } if (minDist < _minDistCutAttach && trackToAttach != NULL) { int iopt = 2; + debug() << " Hit: id = " << hit->getTrackerHit().id().index << " : try attachement"<< endmsg; AttachHitToTrack(trackToAttach,hit,iopt); - } + } + else { + debug() << " Hit: id = " << hit->getTrackerHit().id().index << " rejected due to distance cut of " <<_minDistCutAttach<< " min distance = " << minDist << endmsg; + } } } @@ -2321,6 +2383,7 @@ void SiliconTrackingAlg::TrackingInFTD() { if( iCodeOuter >= _sectorsFTD.size()){ error() << "iCodeOuter index out of range: iCodeOuter = " << iCodeOuter << " _sectorsFTD.size() = " << _sectorsFTD.size() << " exit(1) called from file " << __FILE__ << " line " << __LINE__<< endmsg; + info() << "iS=" << iS << " nLS=" << nLS[0] << " _nlayersFTD=" << _nlayersFTD << " ipOuter=" << ipOuter << endmsg; exit(1); } @@ -2371,12 +2434,12 @@ void SiliconTrackingAlg::TrackingInFTD() { // << "\nMiddle Hit Type "<< hitMiddle->getTrackerHit()->getType() << " z = " << hitMiddle->getTrackerHit().getPosition()[2] // << "\nInner Hit Type "<< hitInner->getTrackerHit()->getType() << " z = " << hitInner->getTrackerHit().getPosition()[2] << endmsg; - debug() << " " - << std::setw(3) << ipOuter << " " << std::setw(3) << ipMiddle << " " << std::setw(3) << ipInner << " " - << std::setw(3) << iS << " " - << std::setw(3) << nLS[0] << " " << std::setw(3) << nLS[1] << " " << std::setw(3) << nLS[2] << " " - << std::setw(3) << nOuter << " : " << std::setw(3) << nMiddle << " : " << std::setw(3) << nInner << " :: " - << std::setw(3) << nOuter*nMiddle* nInner << endmsg; + debug() << "pattern " + << std::setw(3) << ipOuter << " " << std::setw(3) << ipMiddle << " " << std::setw(3) << ipInner << " " + << ((iS==0)?" +z ":" -z ") + << std::setw(3) << nLS[0] << " " << std::setw(3) << nLS[1] << " " << std::setw(3) << nLS[2] << " " + << std::setw(3) << nOuter << " : " << std::setw(3) << nMiddle << " : " << std::setw(3) << nInner << " :: " + << std::setw(3) << nOuter*nMiddle*nInner << endmsg; TrackExtended * trackAR = TestTriplet(hitOuter,hitMiddle,hitInner,helix); @@ -2583,7 +2646,7 @@ void SiliconTrackingAlg::FinalRefit(edm4hep::TrackCollection* trk_col) { TrackerHitExtendedVec& hitVec = trackAR->getTrackerHitExtendedVec(); int nHits = int(hitVec.size()); - //debug() << " Track " << iTrk << " has hits: " << nHits << endmsg; + debug() << " Track " << iTrk << " has hits: " << nHits << endmsg; if( nHits >= _minimalHits) { // int * lh = new int[nHits]; std::vector lh; @@ -2701,10 +2764,12 @@ void SiliconTrackingAlg::FinalRefit(edm4hep::TrackCollection* trk_col) { int nFit = 0; for (int i=0; igetTrackerHit(); + debug() << "TrackerHit " << i << " id = " << trkHit.id().collectionID << "-" << trkHit.id().index << endmsg; // check if the hit has been rejected as being on the same layer and further from the helix lh==0 if (lh[i] == 1) { - edm4hep::TrackerHit trkHit = hitVec[i]->getTrackerHit(); - debug() << "TrackerHit " << i << " id = " << trkHit.id() << endmsg; + //edm4hep::TrackerHit trkHit = hitVec[i]->getTrackerHit(); + debug() << " accept" << endmsg; nFit++; if(trkHit.isAvailable()) { trkHits.push_back(trkHit); @@ -2723,6 +2788,10 @@ void SiliconTrackingAlg::FinalRefit(edm4hep::TrackCollection* trk_col) { debug() << "SiliconTrackingAlg::FinalRefit: Cannot fit less than 3 hits. Number of hits = " << trkHits.size() << endmsg; continue ; } + + // test w/o fit + //continue; + //TrackImpl* Track = new TrackImpl ; //auto track = trk_col->create(); //fucd @@ -2786,6 +2855,7 @@ void SiliconTrackingAlg::FinalRefit(edm4hep::TrackCollection* trk_col) { */ bool fit_backwards = IMarlinTrack::backward; + /* MarlinTrk::IMarlinTrack* marlinTrk = nullptr; try{ marlinTrk = _trksystem->createTrack(); @@ -2794,11 +2864,13 @@ void SiliconTrackingAlg::FinalRefit(edm4hep::TrackCollection* trk_col) { error() << "Cannot create MarlinTrack ! " << endmsg; return; } + */ int status = 0; debug() << "call createFinalisedLCIOTrack now" << endmsg; try { - status = MarlinTrk::createFinalisedLCIOTrack(marlinTrk, trkHits, &track, fit_backwards, covMatrix, _bField, _maxChi2PerHit); + //status = MarlinTrk::createFinalisedLCIOTrack(marlinTrk, trkHits, &track, fit_backwards, covMatrix, _bField, _maxChi2PerHit); + status = m_fitTool->Fit(track, trkHits, covMatrix, _maxChi2PerHit, fit_backwards); } catch (...) { // delete Track; // delete marlinTrk; @@ -2819,29 +2891,49 @@ void SiliconTrackingAlg::FinalRefit(edm4hep::TrackCollection* trk_col) { std::vector > hits_in_fit ; std::vector > outliers ; std::vector all_hits; - all_hits.reserve(300); - - marlinTrk->getHitsInFit(hits_in_fit); - + //all_hits.reserve(300); + + UTIL::BitField64 cellID_encoder( UTIL::ILDCellID0::encoder_string ) ; + + //marlinTrk->getHitsInFit(hits_in_fit); + hits_in_fit = m_fitTool->GetHitsInFit(); + for ( unsigned ihit = 0; ihit < hits_in_fit.size(); ++ihit) { - debug() << "Hit id =" << hits_in_fit[ihit].first.id() << endmsg; - edm4hep::TrackerHit trk = hits_in_fit[ihit].first; - all_hits.push_back(trk);//hits_in_fit[ihit].first); + edm4hep::TrackerHit hit = hits_in_fit[ihit].first; + //all_hits.push_back(hit);//hits_in_fit[ihit].first); + if (m_tupleDebug) { + cellID_encoder.setValue(hit.getCellID()); + m_layer[m_nhit] = cellID_encoder[lcio::ILDCellID0::layer]; + m_chi2[m_nhit] = hits_in_fit[ihit].second; + m_used[m_nhit] = true; + debug() << "Hit id =" << hit.id().index << " layer = " << m_layer[m_nhit] << " det = " << cellID_encoder[lcio::ILDCellID0::subdet] << endmsg; + m_nhit++; + } + all_hits.push_back(hit); } - UTIL::BitField64 cellID_encoder( UTIL::ILDCellID0::encoder_string ) ; - MarlinTrk::addHitNumbersToTrack(&track, all_hits, true, cellID_encoder); - marlinTrk->getOutliers(outliers); - + //marlinTrk->getOutliers(outliers); + outliers = m_fitTool->GetOutliers(); + for ( unsigned ihit = 0; ihit < outliers.size(); ++ihit) { all_hits.push_back(outliers[ihit].first); + if (m_tupleDebug) { + cellID_encoder.setValue(outliers[ihit].first.getCellID()); + m_layer[m_nhit] = cellID_encoder[lcio::ILDCellID0::layer]; + m_chi2[m_nhit] = outliers[ihit].second; + m_used[m_nhit] = false; + m_nhit++; + } } MarlinTrk::addHitNumbersToTrack(&track, all_hits, false, cellID_encoder); - - delete marlinTrk; + if (m_tupleDebug) m_tupleDebug->write(); + + //delete marlinTrk; + m_fitTool->Clear(); + #if EDM4HEP_BUILD_VERSION > EDM4HEP_VERSION(0, 9, 0) int nhits_in_vxd = track.getSubdetectorHitNumbers(0); int nhits_in_ftd = track.getSubdetectorHitNumbers(1); @@ -2852,8 +2944,7 @@ void SiliconTrackingAlg::FinalRefit(edm4hep::TrackCollection* trk_col) { int nhits_in_sit = track.getSubDetectorHitNumbers(2); #endif - //debug() << " Hit numbers for Track "<< track.id() << ": " - debug() << " Hit numbers for Track "<< iTrk <<": " + debug() << " Hit numbers for Track "<< iTrk <<" (id=" << track.id().index << "): " << " vxd hits = " << nhits_in_vxd << " ftd hits = " << nhits_in_ftd << " sit hits = " << nhits_in_sit @@ -2865,13 +2956,13 @@ void SiliconTrackingAlg::FinalRefit(edm4hep::TrackCollection* trk_col) { if( status != IMarlinTrack::success ) { //delete track; - debug() << "FinalRefit: Track fit failed with error code " << status << " track dropped. Number of hits = "<< trkHits.size() << endmsg; + debug() << "FinalRefit: Track fit failed with error code " << status << " track dropped. Number of hits = "<< trkHits.size() << endmsg; continue ; } if( track.getNdf() < 0) { //delete track; - debug() << "FinalRefit: Track fit returns " << track.getNdf() << " degress of freedom track dropped. Number of hits = "<< trkHits.size() << endmsg; + debug() << "FinalRefit: Track fit returns " << track.getNdf() << " degress of freedom track dropped. Number of hits = "<< trkHits.size() << endmsg; //delete track; continue ; } @@ -2885,7 +2976,7 @@ void SiliconTrackingAlg::FinalRefit(edm4hep::TrackCollection* trk_col) { if(trkStateIP.location !=1) continue; /* if (trkStateIP == 0) { - debug() << "SiliconTrackingAlg::FinalRefit: Track fit returns " << track->getNdf() << " degress of freedom track dropped. Number of hits = "<< trkHits.size() << endmsg; + debug() << "SiliconTrackingAlg::FinalRefit: Track fit returns " << track->getNdf() << " degress of freedom track dropped. Number of hits = "<< trkHits.size() << endmsg; throw EVENT::Exception( std::string("SiliconTrackingAlg::FinalRefit: trkStateIP pointer == NULL ") ) ; } */ @@ -2898,8 +2989,8 @@ void SiliconTrackingAlg::FinalRefit(edm4hep::TrackCollection* trk_col) { float cov[15]; - for (int i = 0 ; i<15 ; ++i) { - cov[i] = trkStateIP.covMatrix.operator[](i); + for (int icov = 0 ; icov < 15 ; ++icov) { + cov[icov] = trkStateIP.covMatrix.operator[](icov); } trackAR->setCovMatrix(cov); @@ -2925,8 +3016,7 @@ void SiliconTrackingAlg::FinalRefit(edm4hep::TrackCollection* trk_col) { } if(m_dumpTime&&m_tuple) m_timeKalman = stopwatch.RealTime()*1000; - debug() << " -> run " << _nRun << " event " << _nEvt << endmsg; - debug() << "Number of Si tracks = " << nSiSegments << endmsg; + info() << " -> run " << _nRun << " event " << _nEvt << " Number of Si tracks = " << nSiSegments << endmsg; debug() << "Total 4-momentum of Si Tracks : E = " << std::setprecision(7) << eTot << " Px = " << pxTot << " Py = " << pyTot @@ -3076,6 +3166,9 @@ StatusCode SiliconTrackingAlg::setupGearGeom(){ } } + + info() << "nvxd = " << _nLayersVTX << " nsit = " << _nLayersSIT << " nftd = " << _nlayersFTD << endmsg; + return StatusCode::SUCCESS; } diff --git a/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.h b/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.h index e8e915b56..c457d446b 100644 --- a/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.h +++ b/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.h @@ -23,6 +23,7 @@ #include "DataHelper/HelixClass.h" #include "TrackSystemSvc/IMarlinTrack.h" +#include "Tracking/ITrackFitterTool.h" #include #include @@ -211,6 +212,7 @@ class SiliconTrackingAlg : public GaudiAlgorithm { /** pointer to the IMarlinTrkSystem instance */ MarlinTrk::IMarlinTrkSystem* _trksystem ; + ToolHandle m_fitTool; //bool _runMarlinTrkDiagnostics; //std::string _MarlinTrkDiagnosticsName; typedef std::vector IntVec; @@ -238,7 +240,7 @@ class SiliconTrackingAlg : public GaudiAlgorithm { Gaudi::Property _cutOnPt{this, "CutOnPt", 0.05}; Gaudi::Property _minimalHits{this, "MinimalHits",3}; Gaudi::Property _nHitsChi2{this, "NHitsChi2", 5}; - Gaudi::Property _max_hits_per_sector{this, "MaxHitsPerSector", 100}; + Gaudi::Property _max_hits_per_sector{this, "MaxHitsPerSector", 200}; Gaudi::Property _attachFast{this, "FastAttachment", 0}; Gaudi::Property _useSIT{this, "UseSIT", true}; Gaudi::Property _initialTrackError_d0{this, "InitialTrackErrorD0",1e6}; @@ -253,8 +255,11 @@ class SiliconTrackingAlg : public GaudiAlgorithm { Gaudi::Property _ElossOn{this, "EnergyLossOn", true}; Gaudi::Property _SmoothOn{this, "SmoothOn", true}; Gaudi::Property _helix_max_r{this, "HelixMaxR", 2000.}; - Gaudi::Property m_dumpTime{this, "DumpTime", true}; + Gaudi::Property m_dumpTime{this, "DumpTime", true}; + Gaudi::Property m_debug{this, "Debug", false}; + //std::vector _colours; + Gaudi::Property m_fitToolName{this, "FitterTool", "KalTestTool/KalTest111"}; /** helper function to get collection using try catch block */ //LCCollection* GetCollection( LCEvent * evt, std::string colName ) ; @@ -302,10 +307,15 @@ class SiliconTrackingAlg : public GaudiAlgorithm { std::vector _sectors; std::vector _sectorsFTD; - NTuple::Tuple* m_tuple; + NTuple::Tuple* m_tuple = nullptr; NTuple::Item m_timeTotal; NTuple::Item m_timeKalman; + NTuple::Tuple* m_tupleDebug = nullptr; + NTuple::Item m_nhit; + NTuple::Array m_layer; + NTuple::Array m_chi2; + NTuple::Array m_used; /** * A helper class to allow good code readability by accessing tracks with N hits. * As the smalest valid track contains three hits, but the first index in a vector is 0, diff --git a/Reconstruction/SiliconTracking/src/TrackSubsetAlg.cpp b/Reconstruction/SiliconTracking/src/TrackSubsetAlg.cpp index d4741f2eb..c8b6d1ea9 100644 --- a/Reconstruction/SiliconTracking/src/TrackSubsetAlg.cpp +++ b/Reconstruction/SiliconTracking/src/TrackSubsetAlg.cpp @@ -41,12 +41,16 @@ StatusCode TrackSubsetAlg::initialize() { _nRun = 0 ; _nEvt = 0 ; + // Set up the track fit tool + m_fitTool = ToolHandle(m_fitToolName.value()); + if(m_dumpTime){ NTuplePtr nt1(ntupleSvc(), "MyTuples/Time"+name()); if ( !nt1 ) { m_tuple = ntupleSvc()->book("MyTuples/Time"+name(),CLID_ColumnWiseTuple,"Tracking time"); if ( 0 != m_tuple ) { m_tuple->addItem ("timeTotal", m_timeTotal ).ignore(); + m_tuple->addItem ("timeKalman", m_timeKalman ).ignore(); } else { fatal() << "Cannot book MyTuples/Time"+name() <createTrack(); + //MarlinTrk::IMarlinTrack* marlinTrk = _trkSystem->createTrack(); int error = 0; try { - error = MarlinTrk::createFinalisedLCIOTrack(marlinTrk, trackerHits, &trackImpl, fit_backwards, covMatrix, _bField, _maxChi2PerHit); - + //error = MarlinTrk::createFinalisedLCIOTrack(marlinTrk, trackerHits, &trackImpl, fit_backwards, covMatrix, _bField, _maxChi2PerHit); + error = m_fitTool->Fit(trackImpl, trackerHits, covMatrix, _maxChi2PerHit, fit_backwards); + } catch (...) { // delete Track; @@ -288,10 +293,11 @@ StatusCode TrackSubsetAlg::execute(){ std::vector > hits_in_fit ; std::vector > outliers ; std::vector all_hits; - all_hits.reserve(300); - - marlinTrk->getHitsInFit(hits_in_fit); + //all_hits.reserve(300); + //marlinTrk->getHitsInFit(hits_in_fit); + hits_in_fit = m_fitTool->GetHitsInFit(); + for ( unsigned ihit = 0; ihit < hits_in_fit.size(); ++ihit) { all_hits.push_back(hits_in_fit[ihit].first); } @@ -300,15 +306,17 @@ StatusCode TrackSubsetAlg::execute(){ MarlinTrk::addHitNumbersToTrack(&trackImpl, all_hits, true, cellID_encoder); - marlinTrk->getOutliers(outliers); - + //marlinTrk->getOutliers(outliers); + outliers = m_fitTool->GetOutliers(); + for ( unsigned ihit = 0; ihit < outliers.size(); ++ihit) { all_hits.push_back(outliers[ihit].first); } MarlinTrk::addHitNumbersToTrack(&trackImpl, all_hits, false, cellID_encoder); - delete marlinTrk; + //delete marlinTrk; + m_fitTool->Clear(); if( error != MarlinTrk::IMarlinTrack::success ) { //delete trackImpl; @@ -355,6 +363,7 @@ StatusCode TrackSubsetAlg::execute(){ if(m_dumpTime&&m_tuple){ m_timeTotal = stopwatch.RealTime()*1000; + m_timeKalman = stopwatch_kalman.RealTime()*1000; m_tuple->write(); } diff --git a/Reconstruction/SiliconTracking/src/TrackSubsetAlg.h b/Reconstruction/SiliconTracking/src/TrackSubsetAlg.h index 1023eba97..534524f9b 100644 --- a/Reconstruction/SiliconTracking/src/TrackSubsetAlg.h +++ b/Reconstruction/SiliconTracking/src/TrackSubsetAlg.h @@ -8,6 +8,7 @@ #include "edm4hep/TrackerHitCollection.h" //#include "edm4hep/Track.h" #include "TrackSystemSvc/IMarlinTrkSystem.h" +#include "Tracking/ITrackFitterTool.h" #include "Math/ProbFunc.h" @@ -56,6 +57,7 @@ class TrackSubsetAlg : public GaudiAlgorithm { protected: MarlinTrk::IMarlinTrkSystem* _trkSystem; + ToolHandle m_fitTool; /* Input collection */ std::vector* > _inTrackColHdls; std::vector* > _inTrackerHitColHdls; @@ -76,6 +78,7 @@ class TrackSubsetAlg : public GaudiAlgorithm { Gaudi::Property _maxChi2PerHit{this, "MaxChi2PerHit", 1e2}; Gaudi::Property _omega{this, "Omega", 0.75}; Gaudi::Property m_dumpTime{this, "DumpTime", true}; + Gaudi::Property m_fitToolName{this, "FitterTool", "KalTestTool/KalTest111"}; float _bField; @@ -84,6 +87,7 @@ class TrackSubsetAlg : public GaudiAlgorithm { NTuple::Tuple* m_tuple; NTuple::Item m_timeTotal; + NTuple::Item m_timeKalman; }; /** A functor to return whether two tracks are compatible: The criterion is if the share a TrackerHit or more */ diff --git a/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.cpp b/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.cpp index d218e041b..937adba6a 100644 --- a/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.cpp +++ b/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.cpp @@ -201,6 +201,9 @@ StatusCode ClupatraAlg::initialize() { // don't need, since Gaudi Algorithm will print all Property //printParameters() ; + // Set up the track fit tool + m_fitTool = ToolHandle(m_fitToolName.value()); + auto _trackSystemSvc = service("TrackSystemSvc"); if ( !_trackSystemSvc ) { error() << "Fail to find TrackSystemSvc ..." << endmsg; diff --git a/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.h b/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.h index 962b1c544..fb451483d 100644 --- a/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.h +++ b/Reconstruction/Tracking/src/Clupatra/ClupatraAlg.h @@ -22,6 +22,7 @@ #include "TrackSystemSvc/HelixTrack.h" #include "TrackSystemSvc/HelixFit.h" #include "TrackSystemSvc/IMarlinTrack.h" +#include "Tracking/ITrackFitterTool.h" namespace MarlinTrk{ class IMarlinTrkSystem ; @@ -139,6 +140,7 @@ class ClupatraAlg : public GaudiAlgorithm { Gaudi::Property _MSOn{this, "MultipleScatteringOn", false}; Gaudi::Property _ElossOn{this, "EnergyLossOn", true}; Gaudi::Property _SmoothOn{this, "SmoothOn", false}; + Gaudi::Property m_fitToolName{this, "FitterTool", "KalTestTool/KalTest010"}; DataHandle _TPCHitCollectionHandle{"TPCTrackerHits", Gaudi::DataHandle::Reader, this}; @@ -154,6 +156,7 @@ class ClupatraAlg : public GaudiAlgorithm { int _nEvt ; MarlinTrk::IMarlinTrkSystem* _trksystem ; + ToolHandle m_fitTool; const gear::TPCParameters* _gearTPC ; diff --git a/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp b/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp old mode 100644 new mode 100755 index ee5f1feb2..95265a4cc --- a/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp +++ b/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp @@ -6,7 +6,7 @@ #include #include -#include +#include #include #if __has_include("edm4hep/EDM4hepVersion.h") #include "edm4hep/EDM4hepVersion.h" @@ -94,7 +94,11 @@ std::string toString( int iTrk, edm4hep::Track tpcTrack, float bField=3.5 ) { float pzTPC = helixTPC.getMomentum()[2]; const float ptot = sqrt(pxTPC*pxTPC+pyTPC*pyTPC+pzTPC*pzTPC); +#if EDM4HEP_BUILD_VERSION > EDM4HEP_VERSION(0, 9, 0) + sprintf(strg,"%3i %5i %9.3f %9.3f %9.3f %7.2f %7.2f %7.2f %4i %4i %8.3f %8i",iTrk,tpcTrack.id().index, +#else sprintf(strg,"%3i %5i %9.3f %9.3f %9.3f %7.2f %7.2f %7.2f %4i %4i %8.3f %8i",iTrk,tpcTrack.id(), +#endif ptot, d0TPC,z0TPC,pxTPC,pyTPC,pzTPC,nHits,ndfTPC,Chi2TPC,nlinkedTracks); return std::string( strg ) ; @@ -154,6 +158,9 @@ StatusCode FullLDCTrackingAlg::initialize() { PIOVER2 = 0.5*PI; TWOPI = 2*PI; + // Set up the track fit tool + m_fitTool = ToolHandle(m_fitToolName.value()); + if(m_dumpTime){ NTuplePtr nt1(ntupleSvc(), "MyTuples/Time"+name()); if ( !nt1 ) { @@ -414,13 +421,14 @@ void FullLDCTrackingAlg::AddTrackColToEvt(TrackExtendedVec & trkVec, edm4hep::Tr debug() << "trkHits ready" << endmsg; bool fit_backwards = IMarlinTrack::backward; - MarlinTrk::IMarlinTrack* marlinTrk = _trksystem->createTrack(); - debug() << "marlinTrk created " << endmsg; + //MarlinTrk::IMarlinTrack* marlinTrk = _trksystem->createTrack(); + //debug() << "marlinTrk created " << endmsg; int error_code = 0; try { - error_code = MarlinTrk::createFinalisedLCIOTrack(marlinTrk, trkHits, &track, fit_backwards, &ts_initial, _bField, _maxChi2PerHit); + //error_code = MarlinTrk::createFinalisedLCIOTrack(marlinTrk, trkHits, &track, fit_backwards, &ts_initial, _bField, _maxChi2PerHit); + error_code = m_fitTool->Fit(track, trkHits, ts_initial, _maxChi2PerHit, fit_backwards); } catch (...) { // delete track; @@ -445,8 +453,9 @@ void FullLDCTrackingAlg::AddTrackColToEvt(TrackExtendedVec & trkVec, edm4hep::Tr std::vector all_hits; all_hits.reserve(hits_in_fit.size()); - marlinTrk->getHitsInFit(hits_in_fit); - + //marlinTrk->getHitsInFit(hits_in_fit); + hits_in_fit = m_fitTool->GetHitsInFit(); + for ( unsigned ihit = 0; ihit < hits_in_fit.size(); ++ihit) { all_hits.push_back(hits_in_fit[ihit].first); } @@ -455,17 +464,18 @@ void FullLDCTrackingAlg::AddTrackColToEvt(TrackExtendedVec & trkVec, edm4hep::Tr MarlinTrk::addHitNumbersToTrack(&track, all_hits, true, cellID_encoder); - marlinTrk->getOutliers(outliers); - + //marlinTrk->getOutliers(outliers); + outliers = m_fitTool->GetOutliers(); + for ( unsigned ihit = 0; ihit < outliers.size(); ++ihit) { all_hits.push_back(outliers[ihit].first); } MarlinTrk::addHitNumbersToTrack(&track, all_hits, false, cellID_encoder); - - delete marlinTrk; - + //delete marlinTrk; + m_fitTool->Clear(); + if( error_code != IMarlinTrack::success ) { debug() << "FullLDCTrackingAlg::AddTrackColToEvt: Track fit failed with error code " << error_code << " track dropped. Number of hits = "<< trkHits.size() << endmsg; //delete Track; @@ -534,6 +544,7 @@ void FullLDCTrackingAlg::AddTrackColToEvt(TrackExtendedVec & trkVec, edm4hep::Tr int nhits_in_tpc = track.getSubDetectorHitNumbers(3); int nhits_in_set = track.getSubDetectorHitNumbers(4); #endif + //int nhits_in_vxd = Track->subdetectorHitNumbers()[ 2 * lcio::ILDDetID::VXD - 2 ]; //int nhits_in_ftd = Track->subdetectorHitNumbers()[ 2 * lcio::ILDDetID::FTD - 2 ]; //int nhits_in_sit = Track->subdetectorHitNumbers()[ 2 * lcio::ILDDetID::SIT - 2 ]; @@ -1264,7 +1275,11 @@ void FullLDCTrackingAlg::prepareVectors() { const float pTot = sqrt(pxSi*pxSi+pySi*pySi+pzSi*pzSi); const int ndfSi = trackExt->getNDF(); float Chi2Si = trackExt->getChi2()/float(trackExt->getNDF()); +#if EDM4HEP_BUILD_VERSION > EDM4HEP_VERSION(0, 9, 0) + sprintf(strg,"%3i %5i %9.3f %9.3f %9.3f %7.2f %7.2f %7.2f %4i %4i %8.3f",iTrk, siTrack.id().index, pTot, d0Si,z0Si,pxSi,pySi,pzSi,nHits, ndfSi, Chi2Si); +#else sprintf(strg,"%3i %5i %9.3f %9.3f %9.3f %7.2f %7.2f %7.2f %4i %4i %8.3f",iTrk, siTrack.id(), pTot, d0Si,z0Si,pxSi,pySi,pzSi,nHits, ndfSi, Chi2Si); +#endif debug() << strg << endmsg; if(nHits>0){ diff --git a/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.h b/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.h index cfc804879..2e9b9586b 100755 --- a/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.h +++ b/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.h @@ -19,6 +19,7 @@ #include "TrackSystemSvc/MarlinTrkUtils.h" #include "TrackSystemSvc/IMarlinTrack.h" +#include "Tracking/ITrackFitterTool.h" #include "edm4hep/SimTrackerHitCollection.h" #include "edm4hep/TrackerHitCollection.h" @@ -319,6 +320,7 @@ class FullLDCTrackingAlg : public GaudiAlgorithm { /** pointer to the IMarlinTrkSystem instance */ MarlinTrk::IMarlinTrkSystem* _trksystem ; + ToolHandle m_fitTool; Gaudi::Property _MSOn{this, "MultipleScatteringOn", true}; Gaudi::Property _ElossOn{this, "EnergyLossOn", true}; @@ -397,6 +399,7 @@ class FullLDCTrackingAlg : public GaudiAlgorithm { Gaudi::Property _maxAllowedSiHitRejectionsForTrackCombination{this, "MaxAllowedSiHitRejectionsForTrackCombination", 2}; Gaudi::Property m_dumpTime{this, "DumpTime", true}; //float _dPCutForForcedMerging; + Gaudi::Property m_fitToolName{this, "FitterTool", "KalTestTool/KalTest111"}; bool _runMarlinTrkDiagnostics = false; std::string _MarlinTrkDiagnosticsName; diff --git a/Service/GearSvc/src/GearSvc.h b/Service/GearSvc/src/GearSvc.h index 3048526cb..312491eb4 100644 --- a/Service/GearSvc/src/GearSvc.h +++ b/Service/GearSvc/src/GearSvc.h @@ -10,7 +10,7 @@ class GearSvc : public extends { public: GearSvc(const std::string& name, ISvcLocator* svc); - ~GearSvc(); + virtual ~GearSvc(); gear::GearMgr* getGearMgr() override; @@ -28,6 +28,7 @@ class GearSvc : public extends TGeoNode* FindNode(TGeoNode* mother, char* name); Gaudi::Property m_gearFile{this, "GearXMLFile", ""}; + Gaudi::Property m_field{this, "MagneticField", 0}; gear::GearMgr* m_gearMgr; }; From 5541344dce88070a9cf861a43acf213703f3711a Mon Sep 17 00:00:00 2001 From: Chengdong Fu Date: Mon, 1 Apr 2024 14:22:36 +0800 Subject: [PATCH 06/10] remove very short path --- Simulation/DetSimSD/src/GenericTrackerSensitiveDetector.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Simulation/DetSimSD/src/GenericTrackerSensitiveDetector.cpp b/Simulation/DetSimSD/src/GenericTrackerSensitiveDetector.cpp index 5e9ba1a6a..c1b8aca9a 100644 --- a/Simulation/DetSimSD/src/GenericTrackerSensitiveDetector.cpp +++ b/Simulation/DetSimSD/src/GenericTrackerSensitiveDetector.cpp @@ -33,10 +33,12 @@ G4bool GenericTrackerSensitiveDetector::ProcessHits(G4Step* step, G4TouchableHis dd4hep::Position direction = postPos - prePos; dd4hep::Position position = mean_direction(prePos,postPos); double hit_len = direction.R(); + if (hit_len < 1E-9) return true; if (hit_len > 0) { double new_len = mean_length(h.preMom(),h.postMom())/hit_len; direction *= new_len/hit_len; } + dd4hep::sim::Geant4TrackerHit* hit = nullptr; hit = new dd4hep::sim::Geant4TrackerHit(h.track->GetTrackID(), h.track->GetDefinition()->GetPDGEncoding(), From 949853d2a419d00eec5e43afd05afe09c7886544 Mon Sep 17 00:00:00 2001 From: Chengdong Fu Date: Mon, 1 Apr 2024 14:26:13 +0800 Subject: [PATCH 07/10] add TStraightTrack case --- Utilities/KalTest/src/kaltracklib/TKalDetCradle.cxx | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/Utilities/KalTest/src/kaltracklib/TKalDetCradle.cxx b/Utilities/KalTest/src/kaltracklib/TKalDetCradle.cxx index 0ff61c707..010face67 100644 --- a/Utilities/KalTest/src/kaltracklib/TKalDetCradle.cxx +++ b/Utilities/KalTest/src/kaltracklib/TKalDetCradle.cxx @@ -100,8 +100,11 @@ void TKalDetCradle::Transport(const TKalTrackSite &from, // site from const TVMeasLayer& ml_to = to.GetHit().GetMeasLayer() ; TVector3 x0; this->Transport(from, ml_to, x0, sv, F, Q ) ; - - THelicalTrack hel(sv, x0, to.GetHit().GetBfield()) ; + + double bfield = to.GetHit().GetBfield(); + TVTrack* trk = 0; + if (bfield==0) trk = new TStraightTrack(sv, x0); + else trk = new THelicalTrack(sv, x0, bfield); // --------------------------------------------------------------------- // Move pivot from last expected hit to actural hit at site to @@ -113,14 +116,14 @@ void TKalDetCradle::Transport(const TKalTrackSite &from, // site from Int_t sdim = sv.GetNrows(); // number of track parameters TKalMatrix DF(sdim, sdim); // propagator matrix segment - hel.MoveTo(to.GetPivot(), fid, &DF); // move pivot to actual hit (to) + trk->MoveTo(to.GetPivot(), fid, &DF); // move pivot to actual hit (to) F = DF * F; // update F accordingly - hel.PutInto(sv); // save updated hel to sv + trk->PutInto(sv); // save updated hel to sv } else { to.SetPivot(x0); // if it is a 1-dim hit } - + delete trk; } // From 548a99f5cae1aae71238c12563d3b92753493dca Mon Sep 17 00:00:00 2001 From: Chengdong Fu Date: Mon, 1 Apr 2024 14:56:13 +0800 Subject: [PATCH 08/10] Error to Warning --- Service/GearSvc/src/GearSvc.cpp | 40 ++++++++++++++++--- .../KalDet/src/ild/common/MaterialDataBase.cc | 12 +++++- .../KalDet/src/ild/vxd/ILDVXDKalDetector.cc | 20 +++++++--- 3 files changed, 59 insertions(+), 13 deletions(-) diff --git a/Service/GearSvc/src/GearSvc.cpp b/Service/GearSvc/src/GearSvc.cpp index 695e45549..cf0d6ec9b 100644 --- a/Service/GearSvc/src/GearSvc.cpp +++ b/Service/GearSvc/src/GearSvc.cpp @@ -73,9 +73,15 @@ StatusCode GearSvc::initialize() info() << "Fill GEAR data from GeomSvc" << endmsg; m_gearMgr->setDetectorName("CRD_o1_v01"); - const dd4hep::Direction& field = geomSvc->lcdd()->field().magneticField(dd4hep::Position(0,0,0)); - gear::ConstantBField* b = new gear::ConstantBField(gear::Vector3D(field.x()/dd4hep::tesla, field.y()/dd4hep::tesla, field.z()/dd4hep::tesla)); - m_gearMgr->setBField(b); + if (m_field.value()==0) { + const dd4hep::Direction& field = geomSvc->lcdd()->field().magneticField(dd4hep::Position(0,0,0)); + gear::ConstantBField* b = new gear::ConstantBField(gear::Vector3D(field.x()/dd4hep::tesla, field.y()/dd4hep::tesla, field.z()/dd4hep::tesla)); + m_gearMgr->setBField(b); + } + else { + gear::ConstantBField* b = new gear::ConstantBField(gear::Vector3D(0, 0, m_field.value())); + m_gearMgr->setBField(b); + } dd4hep::DetElement world = geomSvc->getDD4HepGeo(); const std::map& subs = world.children(); @@ -211,8 +217,11 @@ StatusCode GearSvc::convertVXD(dd4hep::DetElement& vxd){ dd4hep::rec::MaterialManager matMgr( dd4hep::Detector::getInstance().world().volume() ) ; const dd4hep::rec::ZPlanarData::LayerLayout& l = vxdData->layers[0] ; - dd4hep::rec::Vector3D a( l.distanceSensitive + l.thicknessSensitive, l.phi0 , 0. , dd4hep::rec::Vector3D::cylindrical ) ; - dd4hep::rec::Vector3D b( l.distanceSupport + l.thicknessSupport, l.phi0 , 0. , dd4hep::rec::Vector3D::cylindrical ) ; + double offset = l.offsetSupport; + //dd4hep::rec::Vector3D a( l.distanceSensitive + l.thicknessSensitive, l.phi0 , 0. , dd4hep::rec::Vector3D::cylindrical ) ; + //dd4hep::rec::Vector3D b( l.distanceSupport + l.thicknessSupport, l.phi0 , 0. , dd4hep::rec::Vector3D::cylindrical ) ; + dd4hep::rec::Vector3D a( l.distanceSensitive + l.thicknessSensitive, l.offsetSupport, 2.*dd4hep::mm); + dd4hep::rec::Vector3D b( l.distanceSupport + l.thicknessSupport, l.offsetSupport, 2.*dd4hep::mm); const dd4hep::rec::MaterialVec& materials = matMgr.materialsBetween( a , b ) ; dd4hep::rec::MaterialData mat = ( materials.size() > 1 ? matMgr.createAveragedMaterial( materials ) : materials[0].first ) ; @@ -228,6 +237,25 @@ StatusCode GearSvc::convertVXD(dd4hep::DetElement& vxd){ mat.interactionLength()/dd4hep::mm); m_gearMgr->registerSimpleMaterial(VXDSupportMaterial); + if (vxdData->rOuterShell>vxdData->rInnerShell) { + dd4hep::rec::Vector3D a1( vxdData->rInnerShell, 0, 2.*dd4hep::mm); + dd4hep::rec::Vector3D b1( vxdData->rOuterShell, 0, 2.*dd4hep::mm); + const dd4hep::rec::MaterialVec& materials1 = matMgr.materialsBetween( a1 , b1 ) ; + dd4hep::rec::MaterialData mat1 = ( materials1.size() > 1 ? matMgr.createAveragedMaterial( materials1 ) : materials1[0].first ) ; + + std::cout << " ####### found materials between points : " << a1 << " and " << b1 << " : " ; + for( unsigned i=0,n=materials1.size();iregisterSimpleMaterial(VXDShellMaterial); + } + info() << vxdData->rInnerShell << " " << vxdData->rOuterShell << " " << vxdData->zHalfShell << " " << vxdData->gapShell << endmsg; for(int i=0,n=vxdData->layers.size(); ilayers[i]; @@ -702,7 +730,7 @@ StatusCode GearSvc::convertFTD(dd4hep::DetElement& ftd){ senRinner, senThickness, senLengthMin, senLengthMax, senWidth, 0); } m_gearMgr->setFTDParameters(ftdParam); - + info() << "nftd = " << nLayers << endmsg; return StatusCode::SUCCESS; } diff --git a/Utilities/KalDet/src/ild/common/MaterialDataBase.cc b/Utilities/KalDet/src/ild/common/MaterialDataBase.cc index f181b8b53..277f17aea 100644 --- a/Utilities/KalDet/src/ild/common/MaterialDataBase.cc +++ b/Utilities/KalDet/src/ild/common/MaterialDataBase.cc @@ -202,7 +202,6 @@ void MaterialDataBase::createMaterials(const gear::GearMgr& gearMgr, IGeomSvc* g TMaterial &ftdsupport = *new TMaterial(name.c_str(), "", A, Z, density, radlen, 0.) ; this->addMaterial(&ftdsupport, name); - // VXD Support Material if(geoSvc){ //obselete @@ -222,9 +221,18 @@ void MaterialDataBase::createMaterials(const gear::GearMgr& gearMgr, IGeomSvc* g name = vxd_sup_mat.getName() ; TMaterial &vxdsupport = *new TMaterial(name.c_str(), "", A, Z, density, radlen, 0.); this->addMaterial(&vxdsupport, name); + + const gear::SimpleMaterial& vxd_shell_mat = gearMgr.getSimpleMaterial("VXDShellMaterial"); + A = vxd_shell_mat.getA(); + Z = vxd_shell_mat.getZ(); + density = vxd_shell_mat.getDensity() * (1000.0/ 1000000.0); // kg/m^3 -> g/cm^3 + radlen = vxd_shell_mat.getRadLength() / 10.0 ; // mm -> cm + name = vxd_shell_mat.getName() ; + TMaterial &vxdshell = *new TMaterial(name.c_str(), "", A, Z, density, radlen, 0.); + this->addMaterial(&vxdshell, name); } catch( gear::UnknownParameterException& e){ - std::cout << "Error while read material from GeomSvc!" << std::endl; + std::cout << "Warning! cannot get material VXDSupportMaterial and VXDShellMaterial from GeomSvc! GearMgr=" << &gearMgr << std::endl; } } } diff --git a/Utilities/KalDet/src/ild/vxd/ILDVXDKalDetector.cc b/Utilities/KalDet/src/ild/vxd/ILDVXDKalDetector.cc index 3043a7f39..3cf666966 100644 --- a/Utilities/KalDet/src/ild/vxd/ILDVXDKalDetector.cc +++ b/Utilities/KalDet/src/ild/vxd/ILDVXDKalDetector.cc @@ -45,7 +45,6 @@ ILDVXDKalDetector::ILDVXDKalDetector( const gear::GearMgr& gearMgr, IGeomSvc* ge TMaterial & beryllium = *MaterialDataBase::Instance().getMaterial("beryllium"); // needed for cryostat - TMaterial & aluminium = *MaterialDataBase::Instance().getMaterial("aluminium"); _vxd_Cryostat.exists = false; @@ -215,6 +214,12 @@ ILDVXDKalDetector::ILDVXDKalDetector( const gear::GearMgr& gearMgr, IGeomSvc* ge } } + if (_vxd_Cryostat.shellInnerR>0&&_vxd_Cryostat.shellThickness>0) { + TMaterial & shell = *MaterialDataBase::Instance().getMaterial("VXDShellMaterial"); + Add( new ILDCylinderMeasLayer(air, shell , _vxd_Cryostat.shellInnerR, _vxd_Cryostat.shelllHalfZ, 0, 0, 0, _bZ, dummy,-1,"VXDShellInnerWall" ) ); + Add( new ILDCylinderMeasLayer(shell, air , _vxd_Cryostat.shellInnerR+_vxd_Cryostat.shellThickness, _vxd_Cryostat.shelllHalfZ, 0, 0, 0, _bZ, dummy,-1,"VXDShellOuterWall" ) ); + } + if (_vxd_Cryostat.exists) { // build Cryostat according to mokka driver vxd04.cc @@ -357,6 +362,11 @@ void ILDVXDKalDetector::setupGearGeom( const gear::GearMgr& gearMgr ){ _relative_position_of_measurement_surface = pVXDDetMain.getDoubleVal( "relative_position_of_measurement_surface" ); } catch (gear::UnknownParameterException& e) {} + + _vxd_Cryostat.shellInnerR = pVXDDetMain.getShellInnerRadius(); + _vxd_Cryostat.shellThickness = pVXDDetMain.getShellOuterRadius() - pVXDDetMain.getShellInnerRadius(); + _vxd_Cryostat.shelllHalfZ = pVXDDetMain.getShellHalfLength(); + try { const gear::GearParameters& pVXDInfra = gearMgr.getGearParameters("VXDInfra"); _vxd_Cryostat.alRadius = pVXDInfra.getDoubleVal( "CryostatAlRadius" ); @@ -365,16 +375,16 @@ void ILDVXDKalDetector::setupGearGeom( const gear::GearMgr& gearMgr ){ _vxd_Cryostat.alZEndCap = pVXDInfra.getDoubleVal( "CryostatAlZEndCap" ); _vxd_Cryostat.alHalfZ = pVXDInfra.getDoubleVal( "CryostatAlHalfZ" ); - _vxd_Cryostat.shellInnerR = pVXDDetMain.getShellInnerRadius(); - _vxd_Cryostat.shellThickness = pVXDDetMain.getShellOuterRadius() - _vxd_Cryostat.shellInnerR; - _vxd_Cryostat.shelllHalfZ = pVXDDetMain.getShellHalfLength(); + //_vxd_Cryostat.shellInnerR = pVXDDetMain.getShellInnerRadius(); + //_vxd_Cryostat.shellThickness = pVXDDetMain.getShellOuterRadius() - _vxd_Cryostat.shellInnerR; + //_vxd_Cryostat.shelllHalfZ = pVXDDetMain.getShellHalfLength(); _vxd_Cryostat.exists = true; //std::cout << "VXDInfra: " << _vxd_Cryostat.alRadius << " " << _vxd_Cryostat.alThickness << " " << _vxd_Cryostat.alInnerR << " " << _vxd_Cryostat.alZEndCap << " " // << _vxd_Cryostat.alHalfZ << " " << _vxd_Cryostat.shellInnerR << " " << _vxd_Cryostat.shellThickness << " " << _vxd_Cryostat.shelllHalfZ << std::endl; } catch (gear::UnknownParameterException& e) { - std::cout << e.what() << std::endl ; + //std::cout << "ILDVXDKalDetector: " << e.what() << ", will not be built" << std::endl ; _vxd_Cryostat.exists = false; } From fccb196f15b148192cac709ff2d8cb2a2c0ed6f1 Mon Sep 17 00:00:00 2001 From: Chengdong Fu Date: Mon, 1 Apr 2024 20:29:49 +0800 Subject: [PATCH 09/10] remove discard code --- .../src/SiliconTrackingAlg.cpp | 40 ------------------- 1 file changed, 40 deletions(-) diff --git a/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.cpp b/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.cpp index c10689121..df203553f 100644 --- a/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.cpp +++ b/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.cpp @@ -233,7 +233,6 @@ StatusCode SiliconTrackingAlg::execute(){ //_trackImplVec.reserve(100); int successVTX = InitialiseVTX(); - //int successFTD = 0; int successFTD = InitialiseFTD(); if (successVTX == 1) { @@ -299,10 +298,6 @@ StatusCode SiliconTrackingAlg::execute(){ debug() << "End of picking up remaining hits " << endmsg; - //edm4hep::TrackCollection* trkCol = nullptr; - //edm4hep::LCRelationCollection* relCol = nullptr; - //auto trkCol = _outColHdl.createAndPut(); - //auto relCol = _outRelColHdl.createAndPut(); /* LCCollectionVec * trkCol = new LCCollectionVec(LCIO::TRACK); // if we want to point back to the hits we need to set the flag @@ -2829,43 +2824,9 @@ void SiliconTrackingAlg::FinalRefit(edm4hep::TrackCollection* trk_col) { for (std::vector< std::pair >::iterator it=r2_values.begin(); it!=r2_values.end(); ++it) { trkHits.push_back(it->second); } - //std::cout << "fucd------------------3 " << _trksystem << std::endl; - //for (unsigned ihit_indx=0 ; ihit_indx < trkHits.size(); ++ihit_indx) { - // std::cout << "fucd trk hit " << *trkHits[ihit_indx] << " " << trkHits[ihit_indx]->getCovMatrix()[0] - // << " " << BitSet32(trkHits[ihit_indx]->getType())[ ILDTrkHitTypeBit::COMPOSITE_SPACEPOINT ] << endmsg; - //} - /* - auto _trackSystemSvc = service("TrackSystemSvc"); - if ( !_trackSystemSvc ) { - error() << "Failed to find TrackSystemSvc ..." << endmsg; - return; - } - _trksystem = _trackSystemSvc->getTrackSystem(); - if( _trksystem == 0 ){ - error() << "Cannot initialize MarlinTrkSystem of Type: KalTest" <setOption( IMarlinTrkSystem::CFG::useQMS, _MSOn ) ; - _trksystem->setOption( IMarlinTrkSystem::CFG::usedEdx, _ElossOn) ; - _trksystem->setOption( IMarlinTrkSystem::CFG::useSmoothing, _SmoothOn) ; - _trksystem->init() ; - */ bool fit_backwards = IMarlinTrack::backward; - /* - MarlinTrk::IMarlinTrack* marlinTrk = nullptr; - try{ - marlinTrk = _trksystem->createTrack(); - } - catch(...){ - error() << "Cannot create MarlinTrack ! " << endmsg; - return; - } - */ - int status = 0; debug() << "call createFinalisedLCIOTrack now" << endmsg; try { @@ -2891,7 +2852,6 @@ void SiliconTrackingAlg::FinalRefit(edm4hep::TrackCollection* trk_col) { std::vector > hits_in_fit ; std::vector > outliers ; std::vector all_hits; - //all_hits.reserve(300); UTIL::BitField64 cellID_encoder( UTIL::ILDCellID0::encoder_string ) ; From 13c93c857d731c2482d0138aa6f0b3d13eb8e74c Mon Sep 17 00:00:00 2001 From: Chengdong Fu Date: Mon, 1 Apr 2024 20:48:11 +0800 Subject: [PATCH 10/10] add switch for debug --- Service/TrackSystemSvc/src/MarlinTrkUtils.cc | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/Service/TrackSystemSvc/src/MarlinTrkUtils.cc b/Service/TrackSystemSvc/src/MarlinTrkUtils.cc index fcacb19cf..bf2ab05d9 100644 --- a/Service/TrackSystemSvc/src/MarlinTrkUtils.cc +++ b/Service/TrackSystemSvc/src/MarlinTrkUtils.cc @@ -608,7 +608,16 @@ namespace MarlinTrk { trkStateCalo.location = MarlinTrk::Location::AtCalorimeter; track->addToTrackStates(trkStateCalo); } else { +#ifdef DEBUG std::cout << " >>>>>>>>>>> MarlinTrk::finaliseLCIOTrack: could not get TrackState at Calo Face " << std::endl ; + if (last_constrained_hit.isAvailable()) { + auto pos = last_constrained_hit.getPosition(); + std::cout << " last_constrained_hit = " << pos.x << "," << pos.y << "," << pos.z << std::endl; + } + else { + std::cout << " last_constrained_hit not Available" << std::endl; + } +#endif //delete trkStateCalo; } } else {