diff --git a/Tracking/Acts/ActsInterop/python/ActsConfigFlags.py b/Tracking/Acts/ActsInterop/python/ActsConfigFlags.py index b440f1f22843c918e706e589825922040689e6e7..4d121ed8abd55912b4da138b2e81cdb2faf4ec03 100644 --- a/Tracking/Acts/ActsInterop/python/ActsConfigFlags.py +++ b/Tracking/Acts/ActsInterop/python/ActsConfigFlags.py @@ -12,5 +12,6 @@ def createActsConfigFlags(): actscf.addFlag('Acts.TrackingGeometry.buildAllAvailableSubDetectors', False) # Build SCT, TRT and Calo if they are available # Track Finding Flags - TO BE ADDED + actscf.addFlag('Acts.TrackFinding.Seeding.WriteNtuple', False) # Write root file with ttree return actscf diff --git a/Tracking/Acts/ActsTrkFinding/CMakeLists.txt b/Tracking/Acts/ActsTrkFinding/CMakeLists.txt index becf15a5dd92f98cdb5e662a80349b5504c4857d..855276d611d4094e242ee12f358a361d16c43051 100644 --- a/Tracking/Acts/ActsTrkFinding/CMakeLists.txt +++ b/Tracking/Acts/ActsTrkFinding/CMakeLists.txt @@ -19,7 +19,7 @@ atlas_add_component( ActsTrkFinding src/*.h src/*.cxx src/components/*.cxx - LINK_LIBRARIES ActsTrkFindingLib ActsInteropLib MagFieldConditions) + LINK_LIBRARIES ActsTrkFindingLib ActsInteropLib MagFieldConditions CxxUtils) # Install files from the package: atlas_install_python_modules( python/*.py POST_BUILD_CMD ${ATLAS_FLAKE8} ) diff --git a/Tracking/Acts/ActsTrkFinding/python/ActsSeedingToolConfig.py b/Tracking/Acts/ActsTrkFinding/python/ActsSeedingToolConfig.py index 4d1fcec64dc3688fd18bdc9d25a9836e5efbeac7..e2a0b524d3ad032d286455b3b88042c458a8328a 100644 --- a/Tracking/Acts/ActsTrkFinding/python/ActsSeedingToolConfig.py +++ b/Tracking/Acts/ActsTrkFinding/python/ActsSeedingToolConfig.py @@ -14,8 +14,11 @@ ActsSeedToolProperties = { 'deltaRMax' : 400 * UnitConstants.mm }, 'ITkPixel' : { - 'rMax' : 320 * UnitConstants.mm, - 'deltaRMax' : 120 * UnitConstants.mm + 'rMax' : 360 * UnitConstants.mm, + 'deltaRMax' : 100 * UnitConstants.mm, + 'zMin' : -3000 * UnitConstants.mm, + 'zMax' : 3000 * UnitConstants.mm, + 'cotThetaMax' : 27.3 }, 'ITkStrip' : { 'rMax' : 1100 * UnitConstants.mm, @@ -69,6 +72,12 @@ def ActsSeedingToolBaseCfg(ConfigFlags, options.setdefault('zBinEdges', []) options.setdefault('zBinNeighborsTop', []) options.setdefault('zBinNeighborsBottom', []) + options.setdefault('WriteNtuple', ConfigFlags.Acts.TrackFinding.Seeding.WriteNtuple) + options.setdefault('TreeFolderName', 'actsValNtuples') + + if ConfigFlags.Acts.TrackFinding.Seeding.WriteNtuple: + HistService = CompFactory.THistSvc( Output = ["actsValNtuples DATAFILE='ActsSeedMakerValidation.root' OPT='RECREATE'"] ) + acc.addService(HistService) ActsTrk__ActsSeedingTool = CompFactory.getComp("ActsTrk::ActsSeedingTool") acc.setPrivateTools(ActsTrk__ActsSeedingTool(**options)) diff --git a/Tracking/Acts/ActsTrkFinding/src/ActsSeedingFromAthenaAlgorithm.cxx b/Tracking/Acts/ActsTrkFinding/src/ActsSeedingFromAthenaAlgorithm.cxx index 43430a7b6d0db330c34effb9955b2f4300b40fa7..2a70b02dc3e4f68f647b3e6ea3cbb2f8192d2569 100644 --- a/Tracking/Acts/ActsTrkFinding/src/ActsSeedingFromAthenaAlgorithm.cxx +++ b/Tracking/Acts/ActsTrkFinding/src/ActsSeedingFromAthenaAlgorithm.cxx @@ -1,5 +1,5 @@ /* - Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration */ #include "src/ActsSeedingFromAthenaAlgorithm.h" @@ -21,6 +21,9 @@ #include +#include "SiSPSeededTrackFinderData/ITkSiSpacePointForSeed.h" +#include "SiSPSeededTrackFinderData/ITkSiSpacePointsProSeed.h" + namespace ActsTrk { ActsSeedingFromAthenaAlgorithm::ActsSeedingFromAthenaAlgorithm( const std::string &name, @@ -43,7 +46,7 @@ namespace ActsTrk { ATH_CHECK( m_seedKey.initialize() ); ATH_CHECK( m_actsSpacePointKey.initialize() ); ATH_CHECK( m_actsSpacePointDataKey.initialize() ); - + return StatusCode::SUCCESS; } @@ -93,8 +96,14 @@ namespace ActsTrk { std::unique_ptr actsSpData = std::make_unique(); + // Convert result to Athena objects + std::unique_ptr< DataVector > spProSeeds = + std::make_unique< DataVector >(); + std::unique_ptr< DataVector > spForSeeds = + std::make_unique< DataVector >(); + // Conversion - // Supported collections are Pixel and SCT clusters + // Supported collections are Pixel and SCT/Strips clusters // The following piece of code converts Athena Clusters to Acts Space Points // This is a temporary solution used only for Acts validation purposes // @@ -106,32 +115,35 @@ namespace ActsTrk { for ( const SpacePointCollection* sp_collection: *spContainer) { for ( const Trk::SpacePoint* sp: *sp_collection ) { - - const Acts::Vector3 globalPos = sp->globalPosition(); - - InDet::SiSpacePointForSeed point; + + const Acts::Vector3 globalPos = sp->globalPosition(); float r[3] = { static_cast(globalPos[0]), static_cast(globalPos[1]), static_cast(globalPos[2]) }; - point.set( sp, r ); - Acts::Vector2 variance(point.covr(), point.covz()); + + std::unique_ptr< ITk::SiSpacePointForSeed > sp_for_seed = + std::make_unique< ITk::SiSpacePointForSeed >( sp, r ); - boost::container::static_vector indexes({counter++}); - std::unique_ptr toAdd = - std::make_unique( globalPos, - variance, + const Acts::Vector2 variance( sp_for_seed->covr(), sp_for_seed->covz() ); + + boost::container::static_vector indexes({counter++}); + std::unique_ptr toAdd = + std::make_unique( globalPos, + variance, *actsSpData.get(), - indexes ); - + indexes); + actsSpContainer->push_back( std::move(toAdd) ); + spForSeeds->push_back( std::move(sp_for_seed) ); + } } ATH_MSG_DEBUG( " \\__ Total of " << counter << " Space Points" ); - // ================================================== // + // ================================================== // // ===================== OUTPUTS ==================== // // ================================================== // @@ -144,7 +156,7 @@ namespace ActsTrk { SG::WriteHandle< ActsTrk::SpacePointData > spacePointDataHandle = SG::makeHandle( m_actsSpacePointDataKey, ctx ); ATH_MSG_DEBUG( " \\__ Space Point Data `" << m_actsSpacePointDataKey.key() << "` created ..." ); - + // ================================================== // // ===================== COMPUTATION ================ // // ================================================== // @@ -156,6 +168,26 @@ namespace ActsTrk { magFieldContext, *seedPtrs.get() ) ); ATH_MSG_DEBUG(" \\__ Created " << seedPtrs->size() << " seeds"); + + // Convert to Athena - TMP + spProSeeds->reserve( seedPtrs->size() ); + for (const ActsTrk::Seed* seed : *seedPtrs.get()) { + std::size_t bottom_idx = seed->sp()[0]->measurementIndexes()[0]; + std::size_t medium_idx = seed->sp()[1]->measurementIndexes()[0]; + std::size_t top_idx = seed->sp()[2]->measurementIndexes()[0]; + + ITk::SiSpacePointForSeed *bottom_sp = spForSeeds->at(bottom_idx); + ITk::SiSpacePointForSeed *medium_sp = spForSeeds->at(medium_idx); + ITk::SiSpacePointForSeed *top_sp = spForSeeds->at(top_idx); + + std::unique_ptr< ITk::SiSpacePointsProSeed > toAdd = + std::make_unique( bottom_sp, + medium_sp, + top_sp, + seed->z()); + + spProSeeds->push_back( std::move(toAdd) ); + } // ================================================== // // ===================== STORE OUTPUT =============== // diff --git a/Tracking/Acts/ActsTrkFinding/src/ActsSeedingTool.cxx b/Tracking/Acts/ActsTrkFinding/src/ActsSeedingTool.cxx index c49287dd6c44d122aa6301dff627836a9e317eaf..a275400997d82173f888549ce8c662c0d4119229 100644 --- a/Tracking/Acts/ActsTrkFinding/src/ActsSeedingTool.cxx +++ b/Tracking/Acts/ActsTrkFinding/src/ActsSeedingTool.cxx @@ -1,5 +1,5 @@ /* - Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration */ #include "src/ActsSeedingTool.h" @@ -40,6 +40,10 @@ namespace ActsTrk { ATH_MSG_INFO( " " << m_minPt ); ATH_MSG_INFO( " " << m_impactMax ); ATH_MSG_INFO( " " << m_numPhiNeighbors ); + ATH_MSG_INFO( " " << m_writeNtuple ); + ATH_MSG_INFO( " " << m_treeFolder ); + + ATH_CHECK( m_evtKey.initialize() ); if (m_zBinEdges.size() - 1 != m_zBinNeighborsTop.size() and @@ -54,12 +58,47 @@ namespace ActsTrk { ATH_MSG_ERROR("Inconsistent config zBinNeighborsBottom"); return StatusCode::FAILURE; } + + // Validation + if (m_writeNtuple) { + if (m_treeFolder.empty()) { + ATH_MSG_ERROR("Requested seed validation, but tree folder is not defined!"); + return StatusCode::FAILURE; + } + + ATH_CHECK( service("THistSvc", m_thistSvc) ); + std::string tree_name = std::string("SeedTree_") + name(); + std::replace( tree_name.begin(), tree_name.end(), '.', '_' ); + + m_outputTree = new TTree( tree_name.c_str() , "ActsSeedMakerValTool"); + + m_outputTree->Branch("eventNumber", &m_event_number, "eventNumber/L"); + m_outputTree->Branch("d0", &m_d0); + m_outputTree->Branch("z0", &m_z0); + m_outputTree->Branch("pt", &m_pt); + m_outputTree->Branch("eta", &m_eta); + m_outputTree->Branch("x1", &m_x1); + m_outputTree->Branch("x2", &m_x2); + m_outputTree->Branch("x3", &m_x3); + m_outputTree->Branch("y1", &m_y1); + m_outputTree->Branch("y2", &m_y2); + m_outputTree->Branch("y3", &m_y3); + m_outputTree->Branch("z1", &m_z1); + m_outputTree->Branch("z2", &m_z2); + m_outputTree->Branch("z3", &m_z3); + m_outputTree->Branch("r1", &m_r1); + m_outputTree->Branch("r2", &m_r2); + m_outputTree->Branch("r3", &m_r3); + + std::string full_tree_name = "/" + m_treeFolder + "/" + tree_name; + ATH_CHECK( m_thistSvc->regTree( full_tree_name.c_str(), m_outputTree ) ); + } return StatusCode::SUCCESS; } StatusCode - ActsSeedingTool::createSeeds(const EventContext& /*ctx*/, + ActsSeedingTool::createSeeds(const EventContext& ctx, const ActsTrk::SpacePointContainer& spContainer, const InDet::BeamSpotData& beamSpotData, const Acts::MagneticFieldContext& magFieldContext, @@ -67,7 +106,8 @@ namespace ActsTrk { { // Create Seeds auto groupSeeds = - createSeeds( spContainer.begin(), + createSeeds( ctx, + spContainer.begin(), spContainer.end(), beamSpotData, magFieldContext ); @@ -90,7 +130,8 @@ namespace ActsTrk { typename external_spacepoint_t > std::vector< Acts::Seed< external_spacepoint_t > > - ActsSeedingTool::createSeeds( spacepoint_iterator_t spBegin, + ActsSeedingTool::createSeeds( const EventContext& ctx, + spacepoint_iterator_t spBegin, spacepoint_iterator_t spEnd, const InDet::BeamSpotData& beamSpotData, const Acts::MagneticFieldContext& magFieldContext ) const { @@ -144,7 +185,50 @@ namespace ActsTrk { finder.createSeedsForGroup(state, std::back_inserter(seeds), group.bottom(), group.middle(), group.top()); } - + + // Validation + if (m_writeNtuple) { + std::lock_guard lock(m_mutex); + + long EvNumber = -1; + SG::ReadHandle eventInfo = SG::makeHandle( m_evtKey, ctx ); + if(eventInfo.isValid()) + EvNumber = eventInfo->eventNumber(); + + for (const auto& seed : seeds) { + auto [pt, eta, d0, z0] = estimateTrackParamsFromSeed(seed, bField); + + auto* bottom = seed.sp()[0]; + auto* medium = seed.sp()[1]; + auto* top = seed.sp()[2]; + + m_event_number = EvNumber; + + m_pt = pt; + m_eta = eta; + m_d0 = d0; + m_z0 = z0; + + // need to add r to sp EDM + m_x1 = bottom->x(); + m_y1 = bottom->y(); + m_z1 = bottom->z(); + m_r1 = std::sqrt(bottom->x()*bottom->x() + bottom->y()*bottom->y()); + + m_x2 = medium->x(); + m_y2 = medium->y(); + m_z2 = medium->z(); + m_r2 = std::sqrt(medium->x()*medium->x() + medium->y()*medium->y()); + + m_x3 = top->x(); + m_y3 = top->y(); + m_z3 = top->z(); + m_r3 = std::sqrt(top->x()*top->x() + top->y()*top->y()); + + m_outputTree->Fill(); + } + } + return seeds; } @@ -198,4 +282,85 @@ namespace ActsTrk { return std::make_pair( gridCfg,finderCfg ); } + template< typename external_spacepoint_t > + std::array + ActsSeedingTool::estimateTrackParamsFromSeed(const Acts::Seed< external_spacepoint_t >& seed, + const Acts::Vector3& magfield) const + { + // Computation is taken from what is done in ATLAS right now. + // We may want to change in the future to the ue of the dedicated function in ACTS Core instead + // Only 4 parameters for now, will add others in the future + + const auto* bottom = seed.sp()[0]; + const auto* medium = seed.sp()[1]; + const auto* top = seed.sp()[2]; + + // Central + float r_medium = std::sqrt(medium->x()*medium->x() + medium->y()*medium->y()); + float z_medium = medium->z(); + + float ax = medium->x() / r_medium; + float ay = medium->y() / r_medium; + + // Bottom + float dx_bottom = medium->x() - bottom->x(); + float dy_bottom = medium->y() - bottom->y(); + float dz_bottom = medium->z() - bottom->z(); + + float x_bottom = dx_bottom * ax + dy_bottom * ay; + float y_bottom = dy_bottom * ax - dx_bottom * ay; + float dxy_bottom = x_bottom * x_bottom + y_bottom * y_bottom; + float r2_bottom = 1. / dxy_bottom; + float dr_bottom = std::sqrt( 1/dxy_bottom ); + + float u_bottom = x_bottom * r2_bottom; + float v_bottom = y_bottom * r2_bottom; + + + // Top + float dx_top = top->x() - medium->x(); + float dy_top = top->y() - medium->y(); + float dz_top = top->z() - medium->z(); + + float x_top = dx_top * ax + dy_top * ay; + float y_top = dy_top * ax - dx_top * ay; + float dxy_top = x_top * x_top + y_top * y_top; + float r2_top = 1/dxy_top; + float dr_top = std::sqrt( 1/dxy_top ); + + float u_top = x_top * r2_top; + float v_top = y_top * r2_top; + + // Put together + float dU = u_top - u_bottom; + float A = (v_top - v_bottom) * dU; + float B = v_bottom - A * u_bottom; + + float K = 2. / (300. * magfield[2]); + float onePlusAsquare = 1. + A * A; + float BSquare = B * B; + + float tz_bottom = dz_bottom * dr_bottom; + float tz_top = dz_top * dr_top; + + float meanOneOverTanThetaSquare= tz_bottom * tz_top; + float meanOneOverTanTheta = std::sqrt(meanOneOverTanThetaSquare); + if (meanOneOverTanTheta < 1e-8) + meanOneOverTanTheta = 1e-8; + + float theta = std::atan(1. / meanOneOverTanTheta); + + + float pt = std::sqrt(onePlusAsquare / BSquare) / (1000 * K); + float eta = -std::log(std::tan(0.5 * theta)); + float d0 = std::abs( (A - B * r_medium) * r_medium ); + + float dZdR = ( top->z() - medium->z() ) * dr_top; + float z0 = z_medium - r_medium * dZdR; + + + std::array params = {pt, eta, d0, z0}; + return params; + } + } // namespace ActsTrk diff --git a/Tracking/Acts/ActsTrkFinding/src/ActsSeedingTool.h b/Tracking/Acts/ActsTrkFinding/src/ActsSeedingTool.h index cfa374316ab24183f86f5e4db677623100beea5b..b899427e069a06780d0fb91d013315e61ced6aed 100644 --- a/Tracking/Acts/ActsTrkFinding/src/ActsSeedingTool.h +++ b/Tracking/Acts/ActsTrkFinding/src/ActsSeedingTool.h @@ -1,5 +1,5 @@ /* - Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration */ #ifndef ACTSTRACKFINDING_ACTSSEEDINGTOOL_H @@ -16,6 +16,10 @@ #include "StoreGate/ReadCondHandleKey.h" #include "ActsGeometry/ATLASMagneticFieldWrapper.h" +#include "xAODEventInfo/EventInfo.h" + +#include "CxxUtils/checker_macros.h" + // ACTS CORE #include "Acts/Definitions/Units.hpp" #include "Acts/Definitions/Common.hpp" @@ -30,6 +34,9 @@ #include "Acts/Seeding/Seedfinder.hpp" #include "Acts/Seeding/Seed.hpp" +#include "GaudiKernel/ServiceHandle.h" +#include "GaudiKernel/ITHistSvc.h" +#include "TTree.h" namespace ActsTrk { @@ -61,7 +68,8 @@ namespace ActsTrk { >::type > std::vector< Acts::Seed< external_spacepoint_t > > - createSeeds( spacepoint_iterator_t spBegin, + createSeeds( const EventContext& ctx, + spacepoint_iterator_t spBegin, spacepoint_iterator_t spEnd, const InDet::BeamSpotData& beamSpotData, const Acts::MagneticFieldContext& magFieldContext ) const; @@ -75,11 +83,21 @@ namespace ActsTrk { > prepareConfiguration( const Acts::Vector2& beamPos, const Acts::Vector3& bField ) const; - + + + // For Validation + template< typename external_spacepoint_t > + std::array + estimateTrackParamsFromSeed(const Acts::Seed< external_spacepoint_t >& seed, + const Acts::Vector3& magfield) const; + // ********************************************************************* // ********************************************************************* private: + // Additional keys + SG::ReadHandleKey< xAOD::EventInfo > m_evtKey {this, "EventInfoKey", "EventInfo"}; + private: // Properties Gaudi::Property< float > m_rMax {this,"rMax",200. * Acts::UnitConstants::mm,""}; Gaudi::Property< float > m_deltaRMin {this,"deltaRMin",1. * Acts::UnitConstants::mm,""}; @@ -98,6 +116,34 @@ namespace ActsTrk { Gaudi::Property< std::vector > m_zBinEdges {this, "zBinEdges", {} , ""}; Gaudi::Property< std::vector> > m_zBinNeighborsTop{this, "zBinNeighborsTop", {}, ""}; Gaudi::Property< std::vector> > m_zBinNeighborsBottom{this, "zBinNeighborsBottom", {}, ""}; + Gaudi::Property< bool > m_writeNtuple {this, "WriteNtuple", false, ""}; + + private: + // Validation + mutable std::mutex m_mutex; + + ITHistSvc* m_thistSvc = nullptr; + TTree* m_outputTree = nullptr; + + Gaudi::Property< std::string > m_treeFolder {this, "TreeFolderName", "", ""}; + + mutable long m_event_number ATLAS_THREAD_SAFE = 0; + mutable float m_d0 ATLAS_THREAD_SAFE = 0; + mutable float m_z0 ATLAS_THREAD_SAFE = 0; + mutable float m_pt ATLAS_THREAD_SAFE = 0; + mutable float m_eta ATLAS_THREAD_SAFE = 0; + mutable double m_x1 ATLAS_THREAD_SAFE = 0; + mutable double m_x2 ATLAS_THREAD_SAFE = 0; + mutable double m_x3 ATLAS_THREAD_SAFE = 0; + mutable double m_y1 ATLAS_THREAD_SAFE = 0; + mutable double m_y2 ATLAS_THREAD_SAFE = 0; + mutable double m_y3 ATLAS_THREAD_SAFE = 0; + mutable double m_z1 ATLAS_THREAD_SAFE = 0; + mutable double m_z2 ATLAS_THREAD_SAFE = 0; + mutable double m_z3 ATLAS_THREAD_SAFE = 0; + mutable double m_r1 ATLAS_THREAD_SAFE = 0; + mutable double m_r2 ATLAS_THREAD_SAFE = 0; + mutable double m_r3 ATLAS_THREAD_SAFE = 0; }; } // namespace