diff --git a/Detector/CMakeLists.txt b/Detector/CMakeLists.txt index e36d83ada704e5db05379fe1d617b525fcf5f8a8..770c4c2ba6a84c1be3917be1623b5de9485b634e 100644 --- a/Detector/CMakeLists.txt +++ b/Detector/CMakeLists.txt @@ -8,3 +8,4 @@ add_subdirectory(DetSegmentation) add_subdirectory(DetGeomSvc) add_subdirectory(DetIdentifier) add_subdirectory(MagneticFieldMap) +add_subdirectory(DetGeomMetaWriter) \ No newline at end of file diff --git a/Detector/DetCRD/compact/CRD_common_v01/Ecal_Crystal_Barrel_v02_02.xml b/Detector/DetCRD/compact/CRD_common_v01/Ecal_Crystal_Barrel_v02_02.xml index 30c16069d9a53f30d44f5db8a9a48b49cb0360ca..c7890968ccc1afa528f6787fcbc7c3bdc0dff758 100755 --- a/Detector/DetCRD/compact/CRD_common_v01/Ecal_Crystal_Barrel_v02_02.xml +++ b/Detector/DetCRD/compact/CRD_common_v01/Ecal_Crystal_Barrel_v02_02.xml @@ -50,7 +50,7 @@ grid_size_x="1*cm" grid_size_y="1*cm" grid_size_z="1*cm"/--> - <id>system:5,module:5,stave:4,dlayer:5,slayer:6,bar:15</id> + <id>system:5,module:5,stave:4,dlayer:5,slayer:6,bar:6</id> </readout> </readouts> diff --git a/Detector/DetCRD/compact/CRD_common_v01/Ecal_Crystal_Endcap_v02_02.xml b/Detector/DetCRD/compact/CRD_common_v01/Ecal_Crystal_Endcap_v02_02.xml index f84ab8c5b44055e6afd065b56af9d19fb7403100..faa9f1f11b33c20012e5f5d98a76ae13d5f0df6e 100755 --- a/Detector/DetCRD/compact/CRD_common_v01/Ecal_Crystal_Endcap_v02_02.xml +++ b/Detector/DetCRD/compact/CRD_common_v01/Ecal_Crystal_Endcap_v02_02.xml @@ -75,7 +75,7 @@ <readouts> <readout name="EcalEndcapsCollection"> <segmentation type="NoSegmentation"/> - <id>system:5,module:1,part:7,stave:7,type:4,dlayer:4,slayer:1,bar:7</id> + <id>system:5,module:1,part:7,stave:7,type:5,dlayer:4,slayer:1,bar:6</id> </readout> </readouts> diff --git a/Detector/DetCRD/compact/TDR_o1_v01/TDR_o1_v01-onlyCalo.xml b/Detector/DetCRD/compact/TDR_o1_v01/TDR_o1_v01-onlyCalo.xml new file mode 100644 index 0000000000000000000000000000000000000000..b7ddf149bf57708be51b2676ab6ec405c386b2f0 --- /dev/null +++ b/Detector/DetCRD/compact/TDR_o1_v01/TDR_o1_v01-onlyCalo.xml @@ -0,0 +1,53 @@ +<?xml version="1.0" encoding="UTF-8"?> +<lccdd xmlns:compact="http://www.lcsim.org/schemas/compact/1.0" + xmlns:xs="http://www.w3.org/2001/XMLSchema" + xs:noNamespaceSchemaLocation="http://www.lcsim.org/schemas/compact/1.0/compact.xsd"> + <info name="TDR_o1_v01" + title="CepC reference detctor for TDR" + author="" + url="http://cepc.ihep.ac.cn" + status="developing" + version="v01"> + <comment>CepC reference detector simulation models used for TDR </comment> + </info> + + <includes> + <gdmlFile ref="${DD4hepINSTALL}/DDDetectors/compact/elements.xml"/> + <gdmlFile ref="../CRD_common_v02/materials.xml"/> + </includes> + + <define> + <constant name="world_size" value="10*m"/> + <constant name="world_x" value="world_size"/> + <constant name="world_y" value="world_size"/> + <constant name="world_z" value="world_size"/> + + <include ref="${DD4hepINSTALL}/DDDetectors/compact/detector_types.xml"/> + </define> + + <include ref="./TDR_Dimensions_v01_01.xml"/> + + <include ref="../CRD_common_v01/Ecal_Crystal_Barrel_v02_02.xml"/> + <include ref="../CRD_common_v01/Ecal_Crystal_Endcap_v02_02.xml"/> + <include ref="../CRD_common_v01/SHcalGlass_Barrel_v05.xml"/> + <include ref="../CRD_common_v01/SHcalGlass_Endcaps_v01.xml"/> + + <fields> + <field name="InnerSolenoid" type="solenoid" + inner_field="Field_nominal_value" + outer_field="0" + zmax="SolenoidCoil_half_length" + inner_radius="SolenoidCoil_center_radius" + outer_radius="Solenoid_outer_radius"> + </field> + <!-- remove anti magnetic field in order to extrapolate to muon detector more easily--> + <!--field name="OuterSolenoid" type="solenoid" + inner_field="0" + outer_field="Field_outer_nominal_value" + zmax="SolenoidCoil_half_length" + inner_radius="Solenoid_outer_radius" + outer_radius="Yoke_barrel_inner_radius"> + </field--> + </fields> + +</lccdd> diff --git a/Detector/DetCRD/compact/TDR_o1_v01/TDR_o1_v01-onlyTracker.xml b/Detector/DetCRD/compact/TDR_o1_v01/TDR_o1_v01-onlyTracker.xml index 3bace0072e3c190c44b8d093d61d0ad93854cbf0..2d0180583e2c4a00b1d4ab0c41ea2c13ee203c4a 100644 --- a/Detector/DetCRD/compact/TDR_o1_v01/TDR_o1_v01-onlyTracker.xml +++ b/Detector/DetCRD/compact/TDR_o1_v01/TDR_o1_v01-onlyTracker.xml @@ -17,9 +17,10 @@ </includes> <define> - <constant name="world_x" value="1.83*m"/> - <constant name="world_y" value="1.83*m"/> - <constant name="world_z" value="8.0*m"/> + <constant name="world_size" value="10*m"/> + <constant name="world_x" value="world_size"/> + <constant name="world_y" value="world_size"/> + <constant name="world_z" value="world_size"/> <include ref="${DD4hepINSTALL}/DDDetectors/compact/detector_types.xml"/> </define> @@ -35,6 +36,13 @@ <include ref="../CRD_common_v01/OTKBarrel_v02.xml"/> <include ref="../CRD_common_v01/OTKEndcap_v02.xml"/> + <!--muon detector--> + <include ref="../CRD_common_v01/Muon_Barrel_v01_04.xml"/> + <include ref="../CRD_common_v01/Muon_Endcap_v01_02.xml"/> + <include ref="../CRD_common_v01/ParaffinEndcap_v01_01.xml"/> + <!--include ref="../CRD_common_v01/ConcreteWall_v01_01.xml"/--> + + <fields> <field name="InnerSolenoid" type="solenoid" inner_field="Field_nominal_value" @@ -43,13 +51,14 @@ inner_radius="SolenoidCoil_center_radius" outer_radius="Solenoid_outer_radius"> </field> - <field name="OuterSolenoid" type="solenoid" + <!-- remove anti magnetic field in order to extrapolate to muon detector more easily--> + <!--field name="OuterSolenoid" type="solenoid" inner_field="0" outer_field="Field_outer_nominal_value" zmax="SolenoidCoil_half_length" inner_radius="Solenoid_outer_radius" outer_radius="Yoke_barrel_inner_radius"> - </field> + </field--> </fields> </lccdd> diff --git a/Detector/DetCRD/scripts/TDR_o1_v01/calodigi.py b/Detector/DetCRD/scripts/TDR_o1_v01/calodigi.py index cdd1ed17abe460ff57573af004dafe02585b25e7..0c0483733fb0584e778d837ddba82504e5091404 100644 --- a/Detector/DetCRD/scripts/TDR_o1_v01/calodigi.py +++ b/Detector/DetCRD/scripts/TDR_o1_v01/calodigi.py @@ -15,7 +15,7 @@ rndmengine.Seeds = seed rndmgensvc = RndmGenSvc("RndmGenSvc") rndmgensvc.Engine = rndmengine.name() -geometry_option = "TDR_o1_v01/TDR_o1_v01.xml" +geometry_option = "TDR_o1_v01/TDR_o1_v01-onlyCalo.xml" if not os.getenv("DETCRDROOT"): print("Can't find the geometry. Please setup envvar DETCRDROOT." ) diff --git a/Detector/DetCRD/scripts/TDR_o1_v01/rec.py b/Detector/DetCRD/scripts/TDR_o1_v01/rec.py index 7f8a427bf3666f0986adbda659487c097b039711..593e9d5c285055830dc17a299af9722a1b184ce0 100644 --- a/Detector/DetCRD/scripts/TDR_o1_v01/rec.py +++ b/Detector/DetCRD/scripts/TDR_o1_v01/rec.py @@ -1,8 +1,13 @@ import os, sys from Gaudi.Configuration import * + +########## CEPCSWData ################# +cepcswdatatop ="/cvmfs/cepcsw.ihep.ac.cn/prototype/releases/data/latest" +####################################### + ############## GeomSvc ################# -geometry_option = "TDR_o1_v01/TDR_o1_v01.xml" +geometry_option = "TDR_o1_v01/TDR_o1_v01-onlyCalo.xml" if not os.getenv("DETCRDROOT"): print("Can't find the geometry. Please setup envvar DETCRDROOT." ) @@ -16,6 +21,8 @@ if not os.path.exists(geometry_path): from Configurables import DetGeomSvc geomsvc = DetGeomSvc("GeomSvc") geomsvc.compact = geometry_path +geomsvc.fastinit = True +geomsvc.metadata_path = os.path.join(cepcswdatatop, "CEPCSWData/offline-data/Detector/DetGeomMetaWriter/data/tdr25.5.0/GeometryMetaData.root") ####################################### ########### k4DataSvc #################### @@ -23,9 +30,6 @@ from Configurables import k4DataSvc podioevent = k4DataSvc("EventDataSvc", input="CaloDigi_TDR_o1_v01_00.root") ########################################## -########## CEPCSWData ################# -cepcswdatatop ="/cvmfs/cepcsw.ihep.ac.cn/prototype/releases/data/latest" -####################################### ########## CrystalEcalEnergyCorrectionSvc ######## from Configurables import CrystalEcalEnergyCorrectionSvc diff --git a/Detector/DetCRD/scripts/TDR_o1_v01/sim.py b/Detector/DetCRD/scripts/TDR_o1_v01/sim.py index 89159f9d7a9ab6f0521cc50f4d04ce444209f79a..51173d1dfbc516f3f85eafc8567c79fc964188eb 100644 --- a/Detector/DetCRD/scripts/TDR_o1_v01/sim.py +++ b/Detector/DetCRD/scripts/TDR_o1_v01/sim.py @@ -124,7 +124,7 @@ cal_sensdettool.CalNamesBirksConstants = [0.008415, 0.008415, 0.01, 0.01] # BGO from Configurables import MarlinEvtSeeder evtseeder = MarlinEvtSeeder("EventSeeder") - + # output from Configurables import PodioOutput out = PodioOutput("outputalg") diff --git a/Detector/DetCRD/scripts/TDR_o1_v01/tracking.py b/Detector/DetCRD/scripts/TDR_o1_v01/tracking.py index b4aae77a91ff4c01f8e03155af1be994d8158a5d..bcec4fc945c56f6b942561ffd5d1e6b0d9c43bb9 100644 --- a/Detector/DetCRD/scripts/TDR_o1_v01/tracking.py +++ b/Detector/DetCRD/scripts/TDR_o1_v01/tracking.py @@ -15,7 +15,7 @@ rndmengine.Seeds = seed rndmgensvc = RndmGenSvc("RndmGenSvc") rndmgensvc.Engine = rndmengine.name() -geometry_option = "TDR_o1_v01/TDR_o1_v01.xml" +geometry_option = "TDR_o1_v01/TDR_o1_v01-onlyTracker.xml" if not os.getenv("DETCRDROOT"): print("Can't find the geometry. Please setup envvar DETCRDROOT." ) diff --git a/Detector/DetGeomMetaWriter/CMakeLists.txt b/Detector/DetGeomMetaWriter/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..d995bc1fb9a4567e7c93dbd4ba80c6646fadaa5b --- /dev/null +++ b/Detector/DetGeomMetaWriter/CMakeLists.txt @@ -0,0 +1,9 @@ +gaudi_add_module(DetGeomMetaWriter + SOURCES src/DetGeomMetaWriter.h + src/DetGeomMetaWriter.cpp + LINK DetInterface + ${DD4hep_COMPONENT_LIBRARIES} + Gaudi::GaudiKernel + EDM4HEP::edm4hep EDM4HEP::edm4hepDict + k4FWCore::k4FWCore +) \ No newline at end of file diff --git a/Detector/DetGeomMetaWriter/script/dumpmeta.py b/Detector/DetGeomMetaWriter/script/dumpmeta.py new file mode 100644 index 0000000000000000000000000000000000000000..2927c8400ebb7fc9f8d47f57a90b8861148b91c9 --- /dev/null +++ b/Detector/DetGeomMetaWriter/script/dumpmeta.py @@ -0,0 +1,44 @@ +#!/usr/bin/env python +import os, sys +from Gaudi.Configuration import * + +from Configurables import k4DataSvc +dsvc = k4DataSvc("EventDataSvc") + +# option for standalone tracker study +#geometry_option = "TDR_o1_v01/TDR_o1_v01-oldVersion.xml" +#geometry_option = "TDR_o1_v01/TDR_o1_v01-patchOTK.xml" +geometry_option = "TDR_o1_v01/TDR_o1_v01.xml" + +if not os.getenv("DETCRDROOT"): + print("Can't find the geometry. Please setup envvar DETCRDROOT." ) + sys.exit(-1) + +geometry_path = os.path.join(os.getenv("DETCRDROOT"), "compact", geometry_option) +if not os.path.exists(geometry_path): + print("Can't find the compact geometry file: %s"%geometry_path) + sys.exit(-1) + +from Configurables import DetGeomSvc +geosvc = DetGeomSvc("GeomSvc") +geosvc.compact = geometry_path + +from Configurables import GeomMetaWriter +metaWriter = GeomMetaWriter("GeomMetaWriter") + +# output +from Configurables import PodioOutput +out = PodioOutput("outputalg") +out.filename = "GeometryMetaData.root" +out.outputCommands = ["keep *"] + +# ApplicationMgr +from Configurables import ApplicationMgr +mgr = ApplicationMgr( + TopAlg = [metaWriter, out], + EvtSel = 'NONE', + EvtMax = 0, + ExtSvc = [dsvc, geosvc], + HistogramPersistency = 'ROOT', + OutputLevel = ERROR +) diff --git a/Detector/DetGeomMetaWriter/src/DetGeomMetaWriter.cpp b/Detector/DetGeomMetaWriter/src/DetGeomMetaWriter.cpp new file mode 100644 index 0000000000000000000000000000000000000000..6536fd0464d9d982f03a2dc6ecbe5eac752a7d88 --- /dev/null +++ b/Detector/DetGeomMetaWriter/src/DetGeomMetaWriter.cpp @@ -0,0 +1,137 @@ +#include "DetGeomMetaWriter.h" +#include "DDRec/CellIDPositionConverter.h" +#include "DDSegmentation/BitFieldCoder.h" +#include "DetInterface/IGeomSvc.h" +#include "GaudiKernel/MsgStream.h" +#include <string> +#include <vector> + +DECLARE_COMPONENT(GeomMetaWriter) + +GeomMetaWriter::GeomMetaWriter(const std::string &aName, ISvcLocator *aSvcLoc) + : GaudiAlgorithm(aName, aSvcLoc) {} + +GeomMetaWriter::~GeomMetaWriter() {} + +void GeomMetaWriter::addDetIDNameMapping() { + auto dets = m_geomSvc->getDD4HepGeo().children(); + for (auto det : dets) { + int detId = det.second.id(); + if (detId != -1) { + _detIDVec.push_back(detId); + _detNameVec.push_back(det.first); + } + } + m_detIdMetaDataHandle.put(_detIDVec); + m_detNameMetaDataHandle.put(_detNameVec); +} + +void GeomMetaWriter::addReadoutDecoderMapping() { + + auto readouts = m_geomSvc->lcdd()->readouts(); + for (auto readout : readouts) { + auto readoutName = readout.first; + auto readoutIdSpec = m_geomSvc->lcdd()->readout(readoutName).idSpec(); + auto decoderStr = readoutIdSpec.decoder()->fieldDescription(); + + _readoutNameVec.push_back(readoutName); + _decoderVec.push_back(decoderStr); + } + + m_readoutNameMetaDataHandle.put(_readoutNameVec); + m_decoderMetaDataHandle.put(_decoderVec); +} + +void GeomMetaWriter::addEcalBarLengthMapping() { + // cellid -> bar length + // FIXME: HardCoded + // ECAL Barrel: system:5,module:5,stave:4,dlayer:5,slayer:6,bar:15 + // system = 20 + // module = 0/1 (odd=1/even=0) + // stave = 0 + // dlayer = 1-9 + // slayer = 0/1 + // bar = 0 + dd4hep::rec::CellID _cellId; + auto _barrel_decoder = m_geomSvc->getDecoder("EcalBarrelCollection"); + for (int module = 0; module <= 1; module++) { + for (int dlayer = 1; dlayer <= 9; dlayer++) { + for (int slayer = 0; slayer <= 1; slayer++) { + _barrel_decoder->set(_cellId, "system", 20); + _barrel_decoder->set(_cellId, "stave", 0); + _barrel_decoder->set(_cellId, "module", module); + _barrel_decoder->set(_cellId, "dlayer", dlayer); + _barrel_decoder->set(_cellId, "slayer", slayer); + _barrel_decoder->set(_cellId, "bar", 0); + debug() << "_cellId: " << _cellId << endmsg; + auto _length = m_geomSvc->getEcalBarLength(_cellId); + auto low32 = dd4hep::DDSegmentation::BitFieldCoder::lowWord(_cellId); + _ecalCellIdVec.push_back(static_cast<int>(low32)); + _ecalBarLengthVec.push_back(_length); + } + } + } + + // ECAL Endcap: + // system:5,module:1,part:7,stave:7,type:4,dlayer:4,slayer:1,bar:7 + // system = 29 + // module = 0/1 + // part = 0-10 + // stave = 10 + // type = 0-16 + // dlayer = 0-8 + // slayer = 0-1 + // bar = 0 + + auto _endcap_decoder = m_geomSvc->getDecoder("EcalEndcapsCollection"); + for (int module = 0; module <= 1; module++) { + for (int part = 0; part <= 10; part++) { + for (int stave = 0; stave <= 10; stave++) { + for (int type = 0; type <= 20; type++) { // FIXME: find accurate upper limit + for (int dlayer = 0; dlayer <= 9; dlayer++) { + for (int slayer = 0; slayer <= 1; slayer++) { + _endcap_decoder->set(_cellId, "system", 29); + _endcap_decoder->set(_cellId, "module", module); + _endcap_decoder->set(_cellId, "part", part); + _endcap_decoder->set(_cellId, "stave", stave); + _endcap_decoder->set(_cellId, "type", type); + _endcap_decoder->set(_cellId, "dlayer", dlayer); + _endcap_decoder->set(_cellId, "slayer", slayer); + _endcap_decoder->set(_cellId, "bar", 0); + try { + auto _length = m_geomSvc->getEcalBarLength(_cellId); + debug() << "_cellId: " << _cellId << endmsg; + debug() << "OK! " << "module: " << module << "part: " << part + << " stave: " << stave << " dlayer: " << dlayer + << " slayer: " << slayer << endmsg; + auto low32 = + dd4hep::DDSegmentation::BitFieldCoder::lowWord(_cellId); + _ecalCellIdVec.push_back(static_cast<int>(low32)); + _ecalBarLengthVec.push_back(_length); + } catch (const std::exception &e) { + continue; + } + } + } + } + } + } + } + + m_ecalCellIDMetaDataHandle.put(_ecalCellIdVec); + m_ecalBarLengthMetaDataHandle.put(_ecalBarLengthVec); +} + +StatusCode GeomMetaWriter::initialize() { + if (GaudiAlgorithm::initialize().isFailure()) { + return StatusCode::FAILURE; + } + addDetIDNameMapping(); + addReadoutDecoderMapping(); + addEcalBarLengthMapping(); + return StatusCode::SUCCESS; +} + +StatusCode GeomMetaWriter::execute() { return StatusCode::SUCCESS; } + +StatusCode GeomMetaWriter::finalize() { return GaudiAlgorithm::finalize(); } diff --git a/Detector/DetGeomMetaWriter/src/DetGeomMetaWriter.h b/Detector/DetGeomMetaWriter/src/DetGeomMetaWriter.h new file mode 100644 index 0000000000000000000000000000000000000000..3974e6ac9e7de3b50aa8801892cf767f62bd4c18 --- /dev/null +++ b/Detector/DetGeomMetaWriter/src/DetGeomMetaWriter.h @@ -0,0 +1,68 @@ +#ifndef GeomMetaWriter_h +#define GeomMetaWriter_h + +// GAUDI +#include "DetInterface/IGeomSvc.h" +#include "GaudiAlg/GaudiAlgorithm.h" +#include "GaudiKernel/SmartIF.h" + +// key4hep +#include <k4FWCore/DataHandle.h> +#include <k4FWCore/MetaDataHandle.h> + +#include <string> + +/** @class GeomMetaWriter + * Lightweight producer for edm data to test cellID + */ +class GeomMetaWriter : public GaudiAlgorithm { +public: + explicit GeomMetaWriter(const std::string &, ISvcLocator *); + virtual ~GeomMetaWriter(); + /** Initialize. + * @return status code + */ + virtual StatusCode initialize() final; + /** Execute. + * @return status code + */ + virtual StatusCode execute() final; + /** Finalize. + * @return status code + */ + virtual StatusCode finalize() final; + +private: + void addDetIDNameMapping(); + void addReadoutDecoderMapping(); + void addEcalBarLengthMapping(); + + MetaDataHandle<std::vector<int>> m_detIdMetaDataHandle{ + "DetIDVector", Gaudi::DataHandle::Writer}; + MetaDataHandle<std::vector<std::string>> + m_detNameMetaDataHandle{"DetNameVector", Gaudi::DataHandle::Writer}; + + MetaDataHandle<std::vector<std::string>> + m_readoutNameMetaDataHandle{"ReadoutNameVector", Gaudi::DataHandle::Writer}; + MetaDataHandle<std::vector<std::string>> + m_decoderMetaDataHandle{"CellIDDecoderStringVector", Gaudi::DataHandle::Writer}; + + MetaDataHandle<std::vector<int>> + m_ecalCellIDMetaDataHandle{"EcalCellIDVector", Gaudi::DataHandle::Writer}; + MetaDataHandle<std::vector<double>> + m_ecalBarLengthMetaDataHandle{"EcalBarLengthVector", Gaudi::DataHandle::Writer}; + + + std::vector<int> _detIDVec; + std::vector<std::string> _detNameVec; + + std::vector<std::string> _readoutNameVec; + std::vector<std::string> _decoderVec; + + std::vector<int> _ecalCellIdVec; + std::vector<double> _ecalBarLengthVec; + + SmartIF<IGeomSvc> m_geomSvc = service<IGeomSvc>("GeomSvc"); + dd4hep::VolumeManager m_volumeManager; +}; +#endif diff --git a/Detector/DetGeomSvc/CMakeLists.txt b/Detector/DetGeomSvc/CMakeLists.txt index 075bee7b5792b3d22f97d9bcf0b6aa8793cdc883..fb86121e97e23de0a538e9036af65f459f8f614b 100644 --- a/Detector/DetGeomSvc/CMakeLists.txt +++ b/Detector/DetGeomSvc/CMakeLists.txt @@ -12,6 +12,7 @@ gaudi_add_module(DetGeomSvc Gaudi::GaudiKernel ${GEAR_LIBRARIES} ${ROOT_LIBRARIES} + k4FWCore::k4FWCore ) target_include_directories(DetGeomSvc diff --git a/Detector/DetGeomSvc/src/DetGeomSvc.cpp b/Detector/DetGeomSvc/src/DetGeomSvc.cpp index 79847faf3c62360fc73fd408b9dcd65f16238efc..c8a2a0033bd1945d1bb97b6d6b04eea3105dda5d 100644 --- a/Detector/DetGeomSvc/src/DetGeomSvc.cpp +++ b/Detector/DetGeomSvc/src/DetGeomSvc.cpp @@ -1,141 +1,231 @@ #include "DetGeomSvc.h" -#include "gearimpl/GearParametersImpl.h" -#include "TMath.h" -#include "TMaterial.h" -#include "CLHEP/Units/SystemOfUnits.h" - -#include "G4Box.hh" -#include "G4Tubs.hh" -#include "G4LogicalVolume.hh" -#include "G4PVPlacement.hh" -#include "G4NistManager.hh" - #include "DD4hep/Detector.h" #include "DD4hep/Plugins.h" -#include "DDG4/Geant4Converter.h" -#include "DDG4/Geant4Mapping.h" -#include "DDRec/DetectorData.h" - -#include <iomanip> -#include <iostream> +#include "GaudiKernel/MsgStream.h" +#include "podio/FrameCategories.h" +#include <DD4hep/Segmentations.h> +#include <DDRec/CellIDPositionConverter.h> +#include <DDRec/DetectorData.h> +#include <string> DECLARE_COMPONENT(DetGeomSvc) -DetGeomSvc::DetGeomSvc(const std::string& name, ISvcLocator* svc) -: base_class(name, svc), m_dd4hep_geo(nullptr){ +DetGeomSvc::DetGeomSvc(const std::string &name, ISvcLocator *svc) + : base_class(name, svc), m_dd4hep_geo(nullptr) {} -} +DetGeomSvc::~DetGeomSvc() {} -DetGeomSvc::~DetGeomSvc() { +void DetGeomSvc::initDetIdToNames() { + if (m_enable_fastinit.value()) { + // fastinit enabled, no dd4hep geo initlized, read from metadata + auto _detIdVec = m_metadata->getValue<std::vector<int>>("DetIDVector"); + auto _detNameVec = + m_metadata->getValue<std::vector<std::string>>("DetNameVector"); + // Convert vectors to map (detID -> detName) + for (size_t i = 0; i < _detIdVec.size() && i < _detNameVec.size(); ++i) { + m_detIdToNames[_detIdVec[i]] = _detNameVec[i]; + } + } else { + auto world = getDD4HepGeo(); + auto subs = world.children(); + for (auto sub : subs) { + int detId = sub.second.id(); + if (detId != -1) + m_detIdToNames[detId] = sub.first; + info() << sub.second.id() << " " << sub.first << endmsg; + } + } +} +void DetGeomSvc::initReadoutNameToDecoder() { + if (m_enable_fastinit.value()) { + // fastinit enabled, no dd4hep geo initlized, read from metadata + auto _readoutNameVec = + m_metadata->getValue<std::vector<std::string>>("ReadoutNameVector"); + auto _decoderVec = m_metadata->getValue<std::vector<std::string>>( + "CellIDDecoderStringVector"); + // Convert vectors to map (detID -> detName) + for (size_t i = 0; i < _readoutNameVec.size() && i < _decoderVec.size(); + ++i) { + m_readoutNameToDecoder[_readoutNameVec[i]] = _decoderVec[i]; + } + } } -StatusCode -DetGeomSvc::initialize() { +StatusCode DetGeomSvc::initialize() { StatusCode sc = Service::initialize(); m_dd4hep_geo = &(dd4hep::Detector::getInstance()); // if DEBUG, set - //dd4hep::PrintLevel level = msgLevel(MSG::INFO) ? dd4hep::printLevel() : dd4hep::setPrintLevel(dd4hep::DEBUG); + // dd4hep::PrintLevel level = msgLevel(MSG::INFO) ? dd4hep::printLevel() : + // dd4hep::setPrintLevel(dd4hep::DEBUG); // if failed to load the compact, a runtime error will be thrown. + + // load the compact file m_dd4hep_geo->fromCompact(m_dd4hep_xmls.value()); // recover to old level, if not, too many DD4hep print - //dd4hep::setPrintLevel(level); - auto world = getDD4HepGeo(); - auto subs = world.children(); - for (auto sub : subs) { - int detId = sub.second.id(); - if (detId!=-1) m_detIdToNames[detId] = sub.first; - info() << sub.second.id() << " " << sub.first << endmsg; + // dd4hep::setPrintLevel(level); + if (m_enable_fastinit.value()) { + podio::ROOTFrameReader m_reader; + m_reader.openFile(m_metadata_path.value()); + auto frame = podio::Frame(m_reader.readEntry(podio::Category::Metadata, 0)); + m_metadata = + std::make_unique<podio::GenericParameters>(frame.getParameters()); } - return sc; } -StatusCode -DetGeomSvc::finalize() { +StatusCode DetGeomSvc::finalize() { StatusCode sc; // m_surface_manager has added as extension of Detector, so not delete? - //if (m_surface_manager != nullptr) delete m_surface_manager; - dd4hep::Detector::destroyInstance(); + // if (m_surface_manager != nullptr) delete m_surface_manager; + // if(m_dd4hep_geo) { + // dd4hep::Detector::destroyInstance(); + // } + // m_volumeManager.destroy(); return sc; } -dd4hep::DetElement -DetGeomSvc::getDD4HepGeo() { - if (lcdd()) { - return lcdd()->world(); - } - return dd4hep::DetElement(); -} - -dd4hep::Detector* -DetGeomSvc::lcdd() { - return m_dd4hep_geo; +dd4hep::DetElement DetGeomSvc::getDD4HepGeo() { + if (lcdd()) { + return lcdd()->world(); + } + return dd4hep::DetElement(); } -IGeomSvc::Decoder* -DetGeomSvc::getDecoder(const std::string& readout_name) { +dd4hep::Detector *DetGeomSvc::lcdd() { return m_dd4hep_geo; } - IGeomSvc::Decoder* decoder = nullptr; +IGeomSvc::Decoder *DetGeomSvc::getDecoder(const std::string &readout_name) { + IGeomSvc::Decoder *decoder = nullptr; + if (m_enable_fastinit.value()) { + if (m_readoutNameToDecoder.empty()) { + initReadoutNameToDecoder(); + } + auto it = m_readoutNameToDecoder.find(readout_name); + if (it != m_readoutNameToDecoder.end()) { + decoder = new IGeomSvc::Decoder(it->second); + } + } else { if (!lcdd()) { - error() << "Failed to get lcdd()" << endmsg; - return decoder; + error() << "Failed to get lcdd()" << endmsg; + return decoder; } auto readouts = m_dd4hep_geo->readouts(); if (readouts.find(readout_name) == readouts.end()) { - error() << "Failed to find readout name '" << readout_name << "'" - << " in DD4hep::readouts. " - << endmsg; - return decoder; + error() << "Failed to find readout name '" << readout_name << "'" + << " in DD4hep::readouts. " << endmsg; + return decoder; } - + dd4hep::Readout readout = lcdd()->readout(readout_name); - auto m_idspec = readout.idSpec(); + auto m_idspec = readout.idSpec(); decoder = m_idspec.decoder(); + } - if (!decoder) { - error() << "Failed to get the decoder with readout '" - << readout_name << "'" << endmsg; - } - - return decoder; + if (!decoder) { + error() << "Failed to get the decoder with readout '" << readout_name << "'" + << endmsg; + } + return decoder; } -const dd4hep::rec::SurfaceMap* -DetGeomSvc::getSurfaceMap(const std::string& det_name) { +const dd4hep::rec::SurfaceMap * +DetGeomSvc::getSurfaceMap(const std::string &det_name) { + if (!m_dd4hep_geo) { + return nullptr; + } if (m_surface_manager == nullptr) { - dd4hep::rec::SurfaceManager* surfaceMgr = nullptr; + dd4hep::rec::SurfaceManager *surfaceMgr = nullptr; // first check whether exist try { surfaceMgr = m_dd4hep_geo->extension<dd4hep::rec::SurfaceManager>(); - } - catch (std::runtime_error& e) { + } catch (std::runtime_error &e) { info() << e.what() << " " << surfaceMgr << endmsg; surfaceMgr = nullptr; } if (surfaceMgr) { m_surface_manager = surfaceMgr; - } - else { - m_dd4hep_geo->addExtension<dd4hep::rec::SurfaceManager>(m_surface_manager = new dd4hep::rec::SurfaceManager(*m_dd4hep_geo)); + } else { + m_dd4hep_geo->addExtension<dd4hep::rec::SurfaceManager>( + m_surface_manager = new dd4hep::rec::SurfaceManager(*m_dd4hep_geo)); } } return m_surface_manager->map(det_name); } -std::string -DetGeomSvc::getDetName(const int det_id) { +std::string DetGeomSvc::getDetName(const int det_id) { + if (m_detIdToNames.empty()) { + initDetIdToNames(); + } auto it = m_detIdToNames.find(det_id); - if (it!=m_detIdToNames.end()) + if (it != m_detIdToNames.end()) { return it->second; - else + } else { return lcdd()->world().name(); + } } + +double DetGeomSvc::getEcalBarLength(unsigned long cellId) { + if (m_enable_fastinit.value()) { + if (m_ecalCellIdToBarLengthMap.empty()) { + // Lazy init + _ecal_barrel_decoder = getDecoder("EcalBarrelCollection"); + _ecal_endcap_decoder = getDecoder("EcalEndcapsCollection"); + auto cellIDs = m_metadata->getValue<std::vector<int>>("EcalCellIDVector"); + auto barLengths = + m_metadata->getValue<std::vector<double>>("EcalBarLengthVector"); + for (size_t i = 0; i < cellIDs.size(); ++i) { + m_ecalCellIdToBarLengthMap[cellIDs[i]] = barLengths[i]; + } + } + dd4hep::rec::CellID _cellId = 0; + auto system = _ecal_barrel_decoder->get(cellId, "system"); + + // Configure decoder based on system type + auto decoder = (system == 20) ? _ecal_barrel_decoder : _ecal_endcap_decoder; + decoder->set(_cellId, "system", system); + decoder->set(_cellId, "bar", 0); + decoder->set(_cellId, "dlayer", decoder->get(cellId, "dlayer")); + decoder->set(_cellId, "slayer", decoder->get(cellId, "slayer")); + + if (system == 20) { + decoder->set(_cellId, "stave", 0); + decoder->set(_cellId, "module", decoder->get(cellId, "module") % 2); + } else if (system == 29) { + decoder->set(_cellId, "module", decoder->get(cellId, "module")); + decoder->set(_cellId, "part", decoder->get(cellId, "part")); + decoder->set(_cellId, "stave", decoder->get(cellId, "stave")); + decoder->set(_cellId, "type", decoder->get(cellId, "type")); + } + + // Lookup and return result + auto it = m_ecalCellIdToBarLengthMap.find(_cellId); + if (it != m_ecalCellIdToBarLengthMap.end()) { + return it->second; + } + + // Debug output for missing cell ID + fatal() << "GetEcalBarLength: Cell ID not found: " << _cellId + << " (system: " << system << ")" << endmsg; + throw GaudiException("GetEcalBarLength: Cell ID not found", name(), StatusCode::FAILURE); + } else { + if (!m_volumeManager) { + m_volumeManager = new dd4hep::VolumeManager( + dd4hep::VolumeManager::getVolumeManager(*m_dd4hep_geo)); + } + // use dd4hep volume + dd4hep::PlacedVolume ipv = m_volumeManager->lookupVolumePlacement(cellId); + dd4hep::Volume ivol = ipv.volume(); + std::vector<double> iVolParam = ivol.solid().dimensions(); + auto maxElement = std::max_element(iVolParam.begin(), iVolParam.end()); + iVolParam.clear(); + return *maxElement * 20; + } +} \ No newline at end of file diff --git a/Detector/DetGeomSvc/src/DetGeomSvc.h b/Detector/DetGeomSvc/src/DetGeomSvc.h index 6c8d00114f9884d3a4319762977bc691ecd401fa..8e497c1be2a5a8ab9624c21e7d1320bbb79bb9d6 100644 --- a/Detector/DetGeomSvc/src/DetGeomSvc.h +++ b/Detector/DetGeomSvc/src/DetGeomSvc.h @@ -5,12 +5,8 @@ #include "DetInterface/IGeomSvc.h" // Gaudi -#include "GaudiKernel/IIncidentListener.h" -#include "GaudiKernel/IIncidentSvc.h" -#include "GaudiKernel/Incident.h" -#include "GaudiKernel/MsgStream.h" +#include "Gaudi/Property.h" #include "GaudiKernel/Service.h" -#include "GaudiKernel/ServiceHandle.h" // DD4Hep #include "DD4hep/Detector.h" @@ -18,6 +14,8 @@ #include <gear/GEAR.h> #include <gearimpl/ZPlanarParametersImpl.h> #include <gearimpl/GearParametersImpl.h> +#include <string> +#include <k4FWCore/MetaDataHandle.h> class TGeoNode; @@ -38,16 +36,29 @@ class DetGeomSvc: public extends<Service, IGeomSvc> { Decoder* getDecoder(const std::string& readout_name) override; const dd4hep::rec::SurfaceMap* getSurfaceMap(const std::string& det_name) override; std::string getDetName(const int det_id) override; + double getEcalBarLength(unsigned long cellId) override; + void initDetIdToNames(); + void initReadoutNameToDecoder(); private: // DD4hep XML compact file path Gaudi::Property<std::string> m_dd4hep_xmls{this, "compact"}; + // init without loading dd4hep geometry + Gaudi::Property<bool> m_enable_fastinit{this, "fastinit", false}; - // dd4hep::Detector* m_dd4hep_geo; dd4hep::rec::SurfaceManager* m_surface_manager = nullptr; + dd4hep::VolumeManager* m_volumeManager = nullptr; - std::map<int, std::string> m_detIdToNames; + std::unordered_map<int, std::string> m_detIdToNames; + std::unordered_map<std::string, std::string> m_readoutNameToDecoder; + std::unordered_map<int, double> m_ecalCellIdToBarLengthMap; + + Gaudi::Property<std::string> m_metadata_path{this, "metadata_path"}; + std::unique_ptr<podio::GenericParameters> m_metadata; + + IGeomSvc::Decoder* _ecal_barrel_decoder; + IGeomSvc::Decoder* _ecal_endcap_decoder; }; #endif // GeomSvc_h diff --git a/Detector/DetInterface/include/DetInterface/IGeomSvc.h b/Detector/DetInterface/include/DetInterface/IGeomSvc.h index b314a2e103eb81dff8e4966fbe4e0dfc8525f589..4da2beb6bbca7f21e5d5a752da00c1048b7f3bfe 100644 --- a/Detector/DetInterface/include/DetInterface/IGeomSvc.h +++ b/Detector/DetInterface/include/DetInterface/IGeomSvc.h @@ -45,6 +45,7 @@ public: virtual const dd4hep::rec::SurfaceMap* getSurfaceMap(const std::string& det_name) = 0; virtual std::string getDetName(const int det_id) = 0; + virtual double getEcalBarLength(unsigned long cellId) = 0; virtual ~IGeomSvc() {} }; diff --git a/Reconstruction/RecPFACyber/include/CyberPFAlg.h b/Reconstruction/RecPFACyber/include/CyberPFAlg.h index 6f1d195ee7fd38f79ca2e882d4d640249e881638..1a13e1d583814d170ab8d62fe2a70d6894c42deb 100644 --- a/Reconstruction/RecPFACyber/include/CyberPFAlg.h +++ b/Reconstruction/RecPFACyber/include/CyberPFAlg.h @@ -93,8 +93,8 @@ protected: SmartIF<ICrystalEcalSvc> m_energycorsvc; std::map<std::string, dd4hep::DDSegmentation::BitFieldCoder*> map_readout_decoder; dd4hep::Detector* m_dd4hep; - dd4hep::rec::CellIDPositionConverter* m_cellIDConverter; - dd4hep::VolumeManager m_volumeManager; + //dd4hep::rec::CellIDPositionConverter* m_cellIDConverter; + //dd4hep::VolumeManager m_volumeManager; std::map<std::tuple<int, int, int, int, int>, int> barNumberMapEndcapMap; //DataCollection: moved into execute() to ensure everything can be cleand after one event. //CyberDataCol m_DataCol; diff --git a/Reconstruction/RecPFACyber/include/Tools/CaloHitsCreator.h b/Reconstruction/RecPFACyber/include/Tools/CaloHitsCreator.h index b4ffdeecb6af7c702d6c52c418ec47b4d69db0c4..05745e20634d49b045a9d659cb7d6018eed106e1 100644 --- a/Reconstruction/RecPFACyber/include/Tools/CaloHitsCreator.h +++ b/Reconstruction/RecPFACyber/include/Tools/CaloHitsCreator.h @@ -4,6 +4,7 @@ #include "k4FWCore/DataHandle.h" #include "CyberDataCol.h" #include "Tools/Algorithm.h" +#include "DetInterface/IGeomSvc.h" #include <DDRec/DetectorData.h> #include <DDRec/CellIDPositionConverter.h> #include <DD4hep/Segmentations.h> @@ -21,7 +22,7 @@ namespace Cyber{ std::vector<DataHandle<edm4hep::CalorimeterHitCollection>*>& r_CaloHitCols, std::map<std::string, dd4hep::DDSegmentation::BitFieldCoder*>& map_decoder, std::map<std::string, DataHandle<edm4hep::MCRecoCaloParticleAssociationCollection>*>& map_CaloParticleAssoCol, - const dd4hep::VolumeManager& m_volumeManager, + SmartIF<IGeomSvc>& m_geosvc, std::map<std::tuple<int, int, int, int, int>, int>& barNumberMapEndcapMap ); //StatusCode CreateMCParticleCaloHitsAsso( std::vector<DataHandle<edm4hep::CalorimeterHitCollection>*>& r_CaloHitCols, diff --git a/Reconstruction/RecPFACyber/src/CyberPFAlg.cpp b/Reconstruction/RecPFACyber/src/CyberPFAlg.cpp index 42f5421ca7a296badc5b9922f4b1f5ff0886babc..ddd0a1c2dc4ebcaa03abbcbe66487dc0a582935c 100644 --- a/Reconstruction/RecPFACyber/src/CyberPFAlg.cpp +++ b/Reconstruction/RecPFACyber/src/CyberPFAlg.cpp @@ -187,8 +187,8 @@ StatusCode CyberPFAlg::initialize() m_dd4hep = m_geosvc->lcdd(); if ( !m_dd4hep ) throw "CyberPFAlg :Failed to get dd4hep::Detector ..."; - m_cellIDConverter = new dd4hep::rec::CellIDPositionConverter(*m_dd4hep); - m_volumeManager = m_dd4hep->volumeManager(); + //m_cellIDConverter = new dd4hep::rec::CellIDPositionConverter(*m_dd4hep); + //m_volumeManager = m_dd4hep->volumeManager(); dd4hep::rec::ECALSystemInfoData* EcalEndcapData = m_geosvc->getDD4HepGeo().child("EcalEndcap").extension<dd4hep::rec::ECALSystemInfoData>(); @@ -742,7 +742,7 @@ StatusCode CyberPFAlg::execute() m_pMCParticleCreator->CreateMCParticle( m_DataCol, *r_MCParticleCol ); if(m_useMCPTrk) m_pTrackCreator->CreateTracksFromMCParticle(m_DataCol, *r_MCParticleCol); else m_pTrackCreator->CreateTracks( m_DataCol, r_TrackCols, r_MCPTrkAssoCol ); - m_pCaloHitsCreator->CreateCaloHits( m_DataCol, r_CaloHitCols, map_readout_decoder, map_CaloMCPAssoCols, m_volumeManager, barNumberMapEndcapMap); + m_pCaloHitsCreator->CreateCaloHits( m_DataCol, r_CaloHitCols, map_readout_decoder, map_CaloMCPAssoCols, m_geosvc, barNumberMapEndcapMap); //Perform PFA algorithm m_algorithmManager.RunAlgorithm( m_DataCol ); @@ -1734,7 +1734,7 @@ StatusCode CyberPFAlg::finalize() //r_HCalHitCols.clear(); r_CaloHitCols.clear(); //m_energycorsvc->finalize(); - delete m_cellIDConverter, m_geosvc; + //delete m_cellIDConverter, m_geosvc; info() << "Processed " << _nEvt << " events " << endmsg; return GaudiAlgorithm::finalize(); } diff --git a/Reconstruction/RecPFACyber/src/Tools/CaloHitsCreator.cpp b/Reconstruction/RecPFACyber/src/Tools/CaloHitsCreator.cpp index cad19fbedcd9f28d9e30da1dd1c1a5721cedd8e6..728d2b3b7b025a851e041d2441b9fdda63a48a7c 100644 --- a/Reconstruction/RecPFACyber/src/Tools/CaloHitsCreator.cpp +++ b/Reconstruction/RecPFACyber/src/Tools/CaloHitsCreator.cpp @@ -12,7 +12,7 @@ namespace Cyber{ std::vector<DataHandle<edm4hep::CalorimeterHitCollection>*>& r_CaloHitCols, std::map<std::string, dd4hep::DDSegmentation::BitFieldCoder*>& map_decoder, std::map<std::string, DataHandle<edm4hep::MCRecoCaloParticleAssociationCollection>*>& map_CaloParticleAssoCol, - const dd4hep::VolumeManager& m_volumeManager, + SmartIF<IGeomSvc>& m_geosvc, std::map<std::tuple<int, int, int, int, int>, int>& barNumberMapEndcapMap ) { if(r_CaloHitCols.size()==0 || settings.map_stringVecPars.at("CaloHitCollections").size()==0) StatusCode::SUCCESS; @@ -79,13 +79,14 @@ namespace Cyber{ m_bar.setQ(hit.getEnergy()/2., hit.getEnergy()/2.); m_bar.setT(hit.getTime(), hit.getTime()); - unsigned long long tmp_id = hit.getCellID(); - dd4hep::PlacedVolume ipv = m_volumeManager.lookupVolumePlacement(tmp_id); - dd4hep::Volume ivol = ipv.volume(); - std::vector< double > iVolParam = ivol.solid().dimensions(); - auto maxElement = std::max_element(iVolParam.begin(), iVolParam.end()); - iVolParam.clear(); - m_bar.setBarLength(*maxElement * 20); + //unsigned long long tmp_id = hit.getCellID(); + //dd4hep::PlacedVolume ipv = m_volumeManager.lookupVolumePlacement(tmp_id); + //dd4hep::Volume ivol = ipv.volume(); + //std::vector< double > iVolParam = ivol.solid().dimensions(); + //auto maxElement = std::max_element(iVolParam.begin(), iVolParam.end()); + //iVolParam.clear(); + //m_bar.setBarLength(*maxElement * 20); + m_bar.setBarLength(m_geosvc->getEcalBarLength(hit.getCellID())); for(int ilink=0; ilink<const_MCPCaloAssoCol->size(); ilink++){ if( hit == const_MCPCaloAssoCol->at(ilink).getRec() ) m_bar.addLinkedMCP( std::make_pair(const_MCPCaloAssoCol->at(ilink).getSim(), const_MCPCaloAssoCol->at(ilink).getWeight()) ); @@ -135,13 +136,14 @@ namespace Cyber{ m_bar.setQ(hit.getEnergy()/2., hit.getEnergy()/2.); m_bar.setT(hit.getTime(), hit.getTime()); - unsigned long long tmp_id = hit.getCellID(); - dd4hep::PlacedVolume ipv = m_volumeManager.lookupVolumePlacement(tmp_id); - dd4hep::Volume ivol = ipv.volume(); - std::vector< double > iVolParam = ivol.solid().dimensions(); - auto maxElement = std::max_element(iVolParam.begin(), iVolParam.end()); - iVolParam.clear(); - m_bar.setBarLength(*maxElement * 20); + //unsigned long long tmp_id = hit.getCellID(); + //dd4hep::PlacedVolume ipv = m_volumeManager.lookupVolumePlacement(tmp_id); + //dd4hep::Volume ivol = ipv.volume(); + //std::vector< double > iVolParam = ivol.solid().dimensions(); + //auto maxElement = std::max_element(iVolParam.begin(), iVolParam.end()); + //iVolParam.clear(); + //m_bar.setBarLength(*maxElement * 20); + m_bar.setBarLength(m_geosvc->getEcalBarLength(hit.getCellID())); for(int ilink=0; ilink<const_MCPCaloAssoCol->size(); ilink++){ if( hit == const_MCPCaloAssoCol->at(ilink).getRec() ) m_bar.addLinkedMCP( std::make_pair(const_MCPCaloAssoCol->at(ilink).getSim(), const_MCPCaloAssoCol->at(ilink).getWeight()) );