diff --git a/CMakeLists.txt b/CMakeLists.txt index 31ddc44..2d41e51 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -19,7 +19,7 @@ limitations under the License. CMAKE_MINIMUM_REQUIRED(VERSION 3.12) -project(k4ProjectTemplate) +project(k4GaudiPandora) # please keep this layout for version setting, used by the automatic tagging script set(PACKAGE_VERSION_MAJOR 1) @@ -33,6 +33,7 @@ find_package(ROOT COMPONENTS RIO Tree) find_package(EDM4HEP) find_package(k4FWCore) find_package(Gaudi) +find_package(DD4hep) #--------------------------------------------------------------- include(GNUInstallDirs) @@ -48,5 +49,6 @@ endif() include(CTest) add_subdirectory(k4ProjectTemplate) +add_subdirectory(k4GaudiPandora) include(cmake/CreateProjectConfig.cmake) diff --git a/k4GaudiPandora/CMakeLists.txt b/k4GaudiPandora/CMakeLists.txt new file mode 100644 index 0000000..cfb3f18 --- /dev/null +++ b/k4GaudiPandora/CMakeLists.txt @@ -0,0 +1,54 @@ +#[[ +Copyright (c) 2020-2024 Key4hep-Project. + +This file is part of Key4hep. +See https://key4hep.github.io/key4hep-doc/ for further info. + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +]] + +set(_plugin_sources ${CMAKE_CURRENT_LIST_DIR}/src/DDSimpleMuonDigi.cc + src/CalorimeterHitType.cc) + +gaudi_add_module(k4GaudiPandoraPlugins + SOURCES ${_plugin_sources} + LINK Gaudi::GaudiKernel + k4FWCore::k4FWCore + k4FWCore::k4Interface + EDM4HEP::edm4hep + DD4hep::DDRec + DD4hep::DDCore + ) +target_include_directories(k4GaudiPandoraPlugins PRIVATE include) +install(TARGETS k4GaudiPandoraPlugins + EXPORT ${CMAKE_PROJECT_NAME}Targets + RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" COMPONENT bin + LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib + PUBLIC_HEADER DESTINATION "${CMAKE_INSTALL_INCLUDEDIR}/@{CMAKE_PROJECT_NAME}" + COMPONENT dev) + +include(CTest) + +#--- The genConf directory has been renamed to genConfDir in Gaudi 35r1 +#--- See https://gitlab.cern.ch/gaudi/Gaudi/-/merge_requests/1158 +set(GAUDI_GENCONF_DIR "genConfDir") + +function(set_test_env _testname) + set_property(TEST ${_testname} APPEND PROPERTY ENVIRONMENT + LD_LIBRARY_PATH=${CMAKE_BINARY_DIR}:$:$:$:$:$:$ENV{LD_LIBRARY_PATH} + PYTHONPATH=${CMAKE_BINARY_DIR}/${CMAKE_PROJECT_NAME}/${GAUDI_GENCONF_DIR}:$/../python:$ENV{PYTHONPATH} + PATH=$/../bin:$ENV{PATH} + K4PROJECTTEMPLATE=${CMAKE_CURRENT_LIST_DIR}/ + ) +endfunction() + diff --git a/k4GaudiPandora/include/CalorimeterHitType.h b/k4GaudiPandora/include/CalorimeterHitType.h new file mode 100644 index 0000000..940d1c4 --- /dev/null +++ b/k4GaudiPandora/include/CalorimeterHitType.h @@ -0,0 +1,129 @@ +/* + * Copyright (c) 2020-2024 Key4hep-Project. + * + * This file is part of Key4hep. + * See https://key4hep.github.io/key4hep-doc/ for further info. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +#ifndef CalorimeterHitType_h +#define CalorimeterHitType_h 1 + +#include + +/** Helper class for decoding/encoding lcio::CalorimeterHit types for the ILD + * detector. The encoding is: caloType + 10 * caloID + 1000 * layout + 10000 * layerNum
+ * (see enums: CaloType, CaloID and Layout for possible values).
+ * Example usage:
+ * 
+ *     lcio::CalorimeterHit* cHit = .... ;
+ *
+ *     // set the type (e.g. in digitization )
+ *     cHit->setType( CHT( CHT::em ,  CHT::ecal , CHT::plug , 12 ) ) ;
+ *
+ *     ...
+ *
+ *     CHT cht = cHit->getType() ;
+ *
+ *     //   sum energies for electromagentic, hadronic and tailcatcher:
+ *     if( cht.is( CHT::em ) )
+ *          e_em +=  cHit->getEnergy() ;
+ *     else
+ *       if ( cht.is(CHT::had ) )
+ *          e_had += cHit->getEnergy() ;
+ *       else
+ *          e_muon += cHit->getEnergy() ;
+ *
+ *     // use only EcalPlug hits:
+ *     if( cht.is( CHT::ecal) && cht.is( CHT::plug) )
+ *
+ *     // get the layer number (e.g. for calibration or clustering)
+ *     unsigned l = cht.layer() ;
+ *     // or directly :
+ *     unsigned l = CHT(  cHit->getType() ).layer()  ;
+ *
+ *     // detailed print:
+ *     std::cout <<  CHT(  cHit->getType() ) << std::endl ;
+ *
+ *
+ * + * F.Gaede, DESY, 12/2008 + */ + +class CHT { +public: + /** calorimeter types */ + enum class CaloType { em = 0, had = 1, muon = 2 }; + + /** calo ids - specific to ILD */ + enum class CaloID { unknown = 0, ecal = 1, hcal = 2, yoke = 3, lcal = 4, lhcal = 5, bcal = 6 }; + + /** calo layout / subdetector */ + enum class Layout { any = 0, barrel = 1, endcap = 2, plug = 3, ring = 4 }; + + /** C'tor for initialization from CalorimeterHit::getType() */ + CHT(int type) : m_type(type) {} + + /** C'tor for encoding the calo type inforamtion */ + CHT(CaloType c, CaloID n, Layout l, unsigned lay) + : m_type(c * fCaloType + n * fCaloID + l * fLayout + lay * fLayer) {} + + /** calorimeter type: CHT::em , CHT::had, CHT::muon */ + CaloType caloType() const { return (CaloType)(m_type % fCaloID); } + + /** calo ID - see enum CaloID for allowed values */ + CaloID caloID() const { return (CaloID)((m_type % fLayout) / fCaloID); } + + /** calo layout - see enum layout for allowed values */ + Layout layout() const { return (Layout)((m_type % fLayer) / fLayout); } + + /** calo layer of hit */ + unsigned layer() const { return unsigned(m_type) / fLayer; } + + bool is(CaloType t) const { return caloType() == t; } + + bool is(CaloID n) const { return caloID() == n; } + + bool is(Layout l) const { return layout() == l; } + + /** automatic conversion to int */ + operator int() const { return m_type; } + + /** explicit conversion to int */ + int toInt() const { return m_type; } + +protected: + int m_type; + + static const int fCaloType = 1; + static const int fCaloID = 10; + static const int fLayout = 1000; + static const int fLayer = 10000; +}; + +/** detailed string for calo type */ +std::ostream& operator<<(std::ostream& os, const CHT& cht); + +/** Return Layout based on the collection name, e.g. if name contains tolower("endcap") CHT::endcap is returned. In case no known layout + is found, CHT::any is returned.*/ +CHT::Layout layoutFromString(const std::string& name); + +/** Return caloID based on the collection name, e.g. if name contains tolower("HCal") CHT::hcal is returned. In case no known type + is found, CHT::unknown is returned.*/ +CHT::CaloID caloIDFromString(const std::string& name); + +/** Return caloType from string, e.g. if name contains tolower("Had") CHT::had is returned. In case no known type + is found, CHT::em is returned.*/ +CHT::CaloType caloTypeFromString(const std::string& name); + +#endif \ No newline at end of file diff --git a/k4GaudiPandora/include/DDSimpleMuonDigi.h b/k4GaudiPandora/include/DDSimpleMuonDigi.h new file mode 100644 index 0000000..8a55f9d --- /dev/null +++ b/k4GaudiPandora/include/DDSimpleMuonDigi.h @@ -0,0 +1,89 @@ +/* + * Copyright (c) 2020-2024 Key4hep-Project. + * + * This file is part of Key4hep. + * See https://key4hep.github.io/key4hep-doc/ for further info. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +#ifndef DDSimpleMuonDigi_H +#define DDSimpleMuonDigi_H + +#include "Gaudi/Property.h" +#include "edm4hep/CaloHitContributionCollection.h" +#include "edm4hep/CaloHitSimCaloHitLinkCollection.h" +#include "edm4hep/CalorimeterHitCollection.h" +#include "edm4hep/EventHeaderCollection.h" +#include "edm4hep/SimCalorimeterHitCollection.h" + +#include "CalorimeterHitType.h" +#include "DDRec/SurfaceManager.h" +#include "k4FWCore/Transformer.h" +#include "k4Interface/IGeoSvc.h" +#include "k4Interface/IUniqueIDGenSvc.h" + +#include +#include +#include + +struct DDSimpleMuonDigi final + : k4FWCore::MultiTransformer< + std::tuple( + const edm4hep::SimCalorimeterHitCollection&, const edm4hep::EventHeaderCollection&)> { + DDSimpleMuonDigi(const std::string& name, ISvcLocator* svcLoc); + + StatusCode initialize() override; + //StatusCode finalize() override; + + std::tuple operator()( + const edm4hep::SimCalorimeterHitCollection& simCaloHits, + const edm4hep::EventHeaderCollection& headers) const override; + +private: + Gaudi::Property m_subDetName{this, "SubDetectorName", "VXD", "Name of the subdetector"}; + Gaudi::Property> m_layersToKeepBarrelVec{ + this, "KeepBarrelLayersVec", {0}, "Vector of Barrel layers to be kept. Layers start at 1!"}; + Gaudi::Property> m_layersToKeepEndCapVec{ + this, "KeepEndcapLayersVec", {0}, "Vector of Endcap layers to be kept. Layers start at 1!"}; + //Gaudi::Property> useLayersBarrelVec{this, "useBarrelLayerVector", false, "whether to use the endcap layer vector"}; + //Gaudi::Property> useLayersEndcapVec{this, "useEndCapLayerVector", false, "whether to use the EndCap layer vector"}; + Gaudi::Property m_muonCollections{this, "muonCollections", "muonCollections", + "The input collection of Muons"}; + Gaudi::Property outputRelCollection{this, "outputRelCollection", "outputRelCollection", + "The output collection of relations"}; + Gaudi::Property outputMuonCollection{this, "outputMuonCollection", "outputMuonCollection", + "The output collection of muons"}; + Gaudi::Property m_encodingStringVariable{ + this, "EncodingStringParameterName", "GlobalTrackerReadoutID", + "The name of the DD4hep constant that contains the Encoding string for tracking detectors"}; + Gaudi::Property m_cellIDLayerString{this, "CellIDLayerString", "Layer", + "Name of the part of the cellID that holds the layer"}; + Gaudi::Property m_thresholdMuon{this, "MuonThreshold", {0.025}, "Threshold for muon"}; + Gaudi::Property m_timeThresholdMuon{this, "timethresholdMuon", {0.025}, "time threshold for muons"}; + Gaudi::Property m_calibrCoeffMuon{ + this, "calibrationCoeffmuon", {120000.0}, "Callibration coefficient of muons"}; + Gaudi::Property m_maxHitEnergyMuon{this, "maxMuonHitEnergy", {2.0}, "Threshold for maximum muon hit energy"}; + Gaudi::Property m_detectorNameBarrel{this, "detectornameB", "YokeBarrel", "Name of the subdetector"}; + Gaudi::Property m_detectorNameEndcap{this, "detectornameE", "YokeEndcap", + "Name of the second subdetector"}; + + std::string m_collName; + std::vector m_useLayersBarrelVec{}, m_useLayersEndcapVec{}; + SmartIF m_geoSvc; + SmartIF m_uidSvc; + + bool useLayer(CHT::Layout caloLayout, unsigned int layer) const; + float computeHitTime(const edm4hep::SimCalorimeterHit& h) const; +}; +DECLARE_COMPONENT(DDSimpleMuonDigi) +#endif \ No newline at end of file diff --git a/k4GaudiPandora/options/runDDSimpleMuonDigi.py b/k4GaudiPandora/options/runDDSimpleMuonDigi.py new file mode 100644 index 0000000..4b4c3ef --- /dev/null +++ b/k4GaudiPandora/options/runDDSimpleMuonDigi.py @@ -0,0 +1,67 @@ +# +# Copyright (c) 2020-2024 Key4hep-Project. +# +# This file is part of Key4hep. +# See https://key4hep.github.io/key4hep-doc/ for further info. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# +from Gaudi.Configuration import INFO +from k4FWCore import ApplicationMgr, IOSvc +from Configurables import EventDataSvc +from Configurables import DDSimpleMuonDigi + +from Configurables import GeoSvc +from Configurables import UniqueIDGenSvc +from Configurables import RootHistSvc +from Configurables import Gaudi__Histograming__Sink__Root as RootHistoSink +import os + +id_service = UniqueIDGenSvc("UniqueIDGenSvc") + +geoservice = GeoSvc("GeoSvc") +geoservice.detectors = [os.environ["K4GEO"]+"/FCCee/CLD/compact/CLD_o2_v06/CLD_o2_v06.xml"] +geoservice.OutputLevel = INFO +geoservice.EnableGeant4Geo = False + +digi = DDSimpleMuonDigi() + +digi.SubDetectorName = "VXD" +digi.KeepBarrelLayersVec = [] +digi.KeepEndcapLayersVec = [] +digi.muonCollections = "ECalBarrelCollection" # "ECalBarrelCollection","ECalEndcapCollection","HCalBarrelCollection","HCalEndcapCollection","HCalRingCollection","LumiCalCollection","YokeBarrelCollection","YokeEndcapCollection" +digi.outputRelCollection = "RelationMuonHit" +digi.outputMuonCollection = "CalorimeterHit" +digi.EncodingStringParameterName = "GlobalTrackerReadoutID" +digi.CellIDLayerString = "layer" +digi.MuonThreshold = 0.025 +digi.timethresholdMuon = 0.025 +digi.calibrationCoeffmuon = 120000.0 +digi.maxMuonHitEnergy = 2.0 +digi.detectornameB = "YokeBarrel" +digi.detectornameE = "YokeEndcap" + +iosvc = IOSvc() +iosvc.input = "../simulation/sim_partgun_1000.root" +iosvc.output = "../outputfiles/output_Gaudi.root" + +hps = RootHistSvc("HistogramPersistencySvc") +root_hist_svc = RootHistoSink("RootHistoSink") +root_hist_svc.FileName = "../outputfiles/ddmuondigi_hist.root" + +ApplicationMgr(TopAlg=[digi], + EvtSel="NONE", + EvtMax=-1, + ExtSvc=[EventDataSvc("EventDataSvc"), root_hist_svc], + OutputLevel=INFO, + ) diff --git a/k4GaudiPandora/src/CalorimeterHitType.cc b/k4GaudiPandora/src/CalorimeterHitType.cc new file mode 100644 index 0000000..610da06 --- /dev/null +++ b/k4GaudiPandora/src/CalorimeterHitType.cc @@ -0,0 +1,138 @@ +/* + * Copyright (c) 2020-2024 Key4hep-Project. + * + * This file is part of Key4hep. + * See https://key4hep.github.io/key4hep-doc/ for further info. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +#include "CalorimeterHitType.h" + +#include + +/** detailed string for calo type */ +std::ostream& operator<<(std::ostream& os, const CHT& cht) { + os << " calo hit type: "; + + switch (cht.caloType()) { + case CHT::em: + os << " em, "; + break; + case CHT::had: + os << " had, "; + break; + case CHT::muon: + os << " muon, "; + break; + default: + os << " - ,"; + } + switch (cht.caloID()) { + case CHT::ecal: + os << "ecal, "; + break; + case CHT::hcal: + os << "hcal, "; + break; + case CHT::yoke: + os << "yoke, "; + break; + case CHT::lcal: + os << "lcal, "; + break; + case CHT::lhcal: + os << "lhcal, "; + break; + case CHT::bcal: + os << "bcal, "; + break; + default: + os << " - ,"; + } + switch (cht.layout()) { + case CHT::any: + os << "any, "; + break; + case CHT::ring: + os << "ring, "; + break; + case CHT::endcap: + os << "endcap, "; + break; + case CHT::barrel: + os << "barrel, "; + break; + case CHT::plug: + os << "plug, "; + break; + default: + os << " - ,"; + } + os << " layer: " << cht.layer(); + + return os; +} + +/** Helper functions that should go to Marlinutil/CalorimeterHitTypes.hh */ + +CHT::Layout layoutFromString(const std::string& name) { + std::string str(name); + std::transform(str.begin(), str.end(), str.begin(), ::tolower); + + if (str.find("ring") != std::string::npos) + return CHT::ring; + if (str.find("plug") != std::string::npos) + return CHT::plug; + if (str.find("endcap") != std::string::npos) + return CHT::endcap; + if (str.find("barrel") != std::string::npos) + return CHT::barrel; + + std::cout << " not found :" << str << " in " << name << std::endl; + return CHT::any; +} + +CHT::CaloID caloIDFromString(const std::string& name) { + std::string str(name); + std::transform(str.begin(), str.end(), str.begin(), ::tolower); + + if (str.find("ecal") != std::string::npos) + return CHT::ecal; + if (str.find("hcal") != std::string::npos) + return CHT::hcal; + if (str.find("yoke") != std::string::npos) + return CHT::yoke; + if (str.find("lcal") != std::string::npos) + return CHT::lcal; + if (str.find("lhcal") != std::string::npos) + return CHT::lhcal; + if (str.find("bcal") != std::string::npos) + return CHT::bcal; + + return CHT::unknown; +} + +CHT::CaloType caloTypeFromString(const std::string& name) { + std::string str(name); + std::transform(str.begin(), str.end(), str.begin(), ::tolower); + + if (str.find("em") != std::string::npos) + return CHT::em; + if (str.find("had") != std::string::npos) + return CHT::had; + if (str.find("muon") != std::string::npos) + return CHT::muon; + + // jl: this should probably also have a separate "unknown" or "any" value? + return CHT::em; +} diff --git a/k4GaudiPandora/src/DDSimpleMuonDigi.cc b/k4GaudiPandora/src/DDSimpleMuonDigi.cc new file mode 100644 index 0000000..e80d75e --- /dev/null +++ b/k4GaudiPandora/src/DDSimpleMuonDigi.cc @@ -0,0 +1,183 @@ +/* + * Copyright (c) 2020-2024 Key4hep-Project. + * + * This file is part of Key4hep. + * See https://key4hep.github.io/key4hep-doc/ for further info. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +#include "DDSimpleMuonDigi.h" +#include +#include // abs +#include "DD4hep/DD4hepUnits.h" +#include "DD4hep/Detector.h" +#include "DDRec/DetectorData.h" +#include "GaudiKernel/MsgStream.h" +#include "edm4hep/CalorimeterHit.h" +#include "edm4hep/Constants.h" + +DECLARE_COMPONENT(DDSimpleMuonDigi) + +DDSimpleMuonDigi::DDSimpleMuonDigi(const std::string& aName, ISvcLocator* aSvcLoc) + : MultiTransformer(aName, aSvcLoc, + { + KeyValues("MUONCollection", {"ECalBarrelCollection"}), + KeyValues("HeaderName", {"EventHeader"}), + }, + {KeyValues("MUONOutputCollections", {"CalorimeterHit"}), + KeyValues("RelationOutputCollection", {"RelationMuonHit"})}) { + m_uidSvc = service("UniqueIDGenSvc", true); + if (!m_uidSvc) { + error() << "Unable to get UniqueIDGenSvc" << endmsg; + } + + m_geoSvc = serviceLocator()->service("GeoSvc"); // important to initialize m_geoSvc +} + +StatusCode DDSimpleMuonDigi::initialize() { + //Get the number of Layers in the Endcap and Barrel + int layersEndcap = 0, layersBarrel = 0; + try { + const auto mainDetector = m_geoSvc->getDetector(); + dd4hep::DetElement theDetector = mainDetector->detector(m_detectorNameBarrel); + const dd4hep::rec::LayeredCalorimeterData* yokeBarrelParameters = + theDetector.extension(); + + layersBarrel = yokeBarrelParameters->layers.size(); + if (!yokeBarrelParameters) { + error() << "oops - yokeBarrelParameters is a null pointer" << endmsg; + } + + } catch (std::exception& e) { + debug() << " oops - no Yoke Barrel available: " << e.what() << std::endl; + } + try { + const auto mainDetector = m_geoSvc->getDetector(); + dd4hep::DetElement theDetector = mainDetector->detector(m_detectorNameEndcap); + const dd4hep::rec::LayeredCalorimeterData* yokeEndcapParameters = + theDetector.extension(); + layersEndcap = yokeEndcapParameters->layers.size(); + } catch (std::exception& e) { + debug() << " oops - no Yoke Endcap available: " << e.what() << std::endl; + } + + //If the vectors are empty, we are keeping everything + if (m_layersToKeepBarrelVec.size() > 0) { + //layers start at 0 + m_useLayersBarrelVec = std::vector(layersBarrel, false); + for (int k : m_layersToKeepBarrelVec) { + m_useLayersBarrelVec[k - 1] = true; + } + } + if (m_layersToKeepEndCapVec.size() > 0) { + //layers start at 0 + m_useLayersEndcapVec = std::vector(layersEndcap, false); + for (int k : m_layersToKeepEndCapVec) { + m_useLayersEndcapVec[k - 1] = true; + } + } + return StatusCode::SUCCESS; +} +std::tuple DDSimpleMuonDigi::operator()( + const edm4hep::SimCalorimeterHitCollection& SimCaloHits, const edm4hep::EventHeaderCollection& headers) const { + debug() << " process event : " << headers[0].getEventNumber() << " - run " << headers[0].getRunNumber() + << endmsg; // headers[0].getRunNumber(),headers[0].getEventNumber() + + auto muoncol = edm4hep::CalorimeterHitCollection(); + auto muonRelcol = edm4hep::CaloHitSimCaloHitLinkCollection(); + + std::string initString; + + std::string colName = m_muonCollections; + CHT::Layout caloLayout = layoutFromString(colName); + + //auto col = headers[0].getCollection(m_muonCollections[i].c_str()); + initString = m_geoSvc->constantAsString(m_encodingStringVariable.value()); + dd4hep::DDSegmentation::BitFieldCoder bitFieldCoder(initString); // check! + // check if decoder contains "layer" + + for (const auto& hit : SimCaloHits) { + const int cellID = hit.getCellID(); + float energy = hit.getEnergy(); + //Get the layer number + unsigned int layer = bitFieldCoder.get(cellID, "layer"); + //Check if we want to use this later, else go to the next hit + if (!useLayer(caloLayout, layer)) + continue; + //Do the digitalization + float calibr_coeff = 1.; + calibr_coeff = m_calibrCoeffMuon; + float hitEnergy = calibr_coeff * energy; + if (hitEnergy > m_maxHitEnergyMuon) { + hitEnergy = m_maxHitEnergyMuon; + } + if (hitEnergy > m_thresholdMuon) { + edm4hep::MutableCalorimeterHit calHit = muoncol.create(); + calHit.setCellID(cellID); + calHit.setEnergy(hitEnergy); + calHit.setPosition(hit.getPosition()); + calHit.setType(CHT(CHT::muon, CHT::yoke, caloLayout, layer)); + calHit.setTime(computeHitTime(hit)); + auto muonRel = muonRelcol.create(); + muonRel.setFrom(calHit); + muonRel.setTo(hit); + } + } + + return std::make_tuple(std::move(muoncol), std::move(muonRelcol)); +} + +//StatusCode DDSimpleMuonDigi::finalize() { return StatusCode::SUCCESS; } + +//If the vectors are empty, we are keeping everything +bool DDSimpleMuonDigi::useLayer(CHT::Layout caloLayout, unsigned int layer) const { + switch (caloLayout) { + case CHT::barrel: + if (layer > m_useLayersBarrelVec.size() || m_useLayersBarrelVec.size() == 0) + return true; + return m_useLayersBarrelVec[layer]; //break not needed, because of return + case CHT::endcap: + if (layer > m_useLayersEndcapVec.size() || m_useLayersEndcapVec.size() == 0) + return true; + return m_useLayersEndcapVec[layer]; //break not needed, because of return + //For all other cases, always keep the hit + default: + return true; + } +} //useLayer + +float DDSimpleMuonDigi::computeHitTime(const edm4hep::SimCalorimeterHit& h) const { + // Sort sim hit MC contribution by time. + // Accumulate the energy from earliest time till the energy + // threshold is reached. The hit time is then estimated at this position in the array + using entry_type = std::pair; + std::vector timeToEnergyMapping{}; + auto singleHit = h.getContributions(); + + const unsigned int nContribs = singleHit.size(); + timeToEnergyMapping.reserve(nContribs); + + for (unsigned int c = 0; c < nContribs; ++c) { + timeToEnergyMapping.push_back({singleHit[c].getTime(), singleHit[c].getEnergy()}); + } + std::sort(timeToEnergyMapping.begin(), timeToEnergyMapping.end(), + [this](entry_type& lhs, entry_type& rhs) { return lhs.first < rhs.first; }); + float energySum = 0.f; + for (auto& entry : timeToEnergyMapping) { + energySum += entry.second * m_calibrCoeffMuon; + if (energySum > m_timeThresholdMuon) { + return entry.first; + } + } + return 0.f; +}