From 9d5250567b50b64ad6f815fecd0fdc583bd1f66b Mon Sep 17 00:00:00 2001 From: Federico Meloni Date: Wed, 6 Aug 2025 14:27:42 +0200 Subject: [PATCH 01/14] update podio req --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 819860d..0911731 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -29,7 +29,7 @@ set(${PROJECT_NAME}_VERSION "${${PROJECT_NAME}_VERSION_MAJOR}.${${PROJECT_NAME}_ option(BUILD_TRACKING "Build conformal tracking" ON) -find_package(podio 1.2.99 REQUIRED) +find_package(podio 1.3 REQUIRED) find_package(EDM4HEP REQUIRED) find_package(DD4hep REQUIRED) find_package(Gaudi REQUIRED) From 63650ffb74645eb1f024656fcd562a9cacbf11f7 Mon Sep 17 00:00:00 2001 From: samf25 Date: Thu, 31 Jul 2025 13:30:00 -0400 Subject: [PATCH 02/14] Calo and OTRM --- k4Reco/CMakeLists.txt | 5 + k4Reco/CaloDigi/README.md | 0 k4Reco/CaloDigi/include/CalorimeterHitType.h | 129 +++++ .../CaloDigi/include/FilterDoubleLayerHits.h | 95 ++++ k4Reco/CaloDigi/include/RealisticCaloDigi.h | 118 +++++ .../include/RealisticCaloDigiScinPpd.h | 31 ++ .../include/RealisticCaloDigiSilicon.h | 27 + k4Reco/CaloDigi/include/RealisticCaloReco.h | 72 +++ .../include/RealisticCaloRecoScinPpd.h | 23 + .../include/RealisticCaloRecoSilicon.h | 23 + k4Reco/CaloDigi/src/CalorimeterHitType.cc | 138 +++++ k4Reco/CaloDigi/src/FilterDoubleLayerHits.cc | 293 +++++++++++ k4Reco/CaloDigi/src/RealisticCaloDigi.cc | 302 +++++++++++ .../CaloDigi/src/RealisticCaloDigiScinPpd.cc | 49 ++ .../CaloDigi/src/RealisticCaloDigiSilicon.cc | 44 ++ k4Reco/CaloDigi/src/RealisticCaloReco.cc | 92 ++++ .../CaloDigi/src/RealisticCaloRecoScinPpd.cc | 31 ++ .../CaloDigi/src/RealisticCaloRecoSilicon.cc | 17 + .../components/OverlayTimingRandomMix.cpp | 471 ++++++++++++++++++ .../components/OverlayTimingRandomMix.h | 152 ++++++ 20 files changed, 2112 insertions(+) create mode 100644 k4Reco/CaloDigi/README.md create mode 100644 k4Reco/CaloDigi/include/CalorimeterHitType.h create mode 100644 k4Reco/CaloDigi/include/FilterDoubleLayerHits.h create mode 100644 k4Reco/CaloDigi/include/RealisticCaloDigi.h create mode 100644 k4Reco/CaloDigi/include/RealisticCaloDigiScinPpd.h create mode 100644 k4Reco/CaloDigi/include/RealisticCaloDigiSilicon.h create mode 100644 k4Reco/CaloDigi/include/RealisticCaloReco.h create mode 100644 k4Reco/CaloDigi/include/RealisticCaloRecoScinPpd.h create mode 100644 k4Reco/CaloDigi/include/RealisticCaloRecoSilicon.h create mode 100644 k4Reco/CaloDigi/src/CalorimeterHitType.cc create mode 100644 k4Reco/CaloDigi/src/FilterDoubleLayerHits.cc create mode 100644 k4Reco/CaloDigi/src/RealisticCaloDigi.cc create mode 100644 k4Reco/CaloDigi/src/RealisticCaloDigiScinPpd.cc create mode 100644 k4Reco/CaloDigi/src/RealisticCaloDigiSilicon.cc create mode 100644 k4Reco/CaloDigi/src/RealisticCaloReco.cc create mode 100644 k4Reco/CaloDigi/src/RealisticCaloRecoScinPpd.cc create mode 100644 k4Reco/CaloDigi/src/RealisticCaloRecoSilicon.cc create mode 100644 k4Reco/Overlay/components/OverlayTimingRandomMix.cpp create mode 100644 k4Reco/Overlay/components/OverlayTimingRandomMix.h diff --git a/k4Reco/CMakeLists.txt b/k4Reco/CMakeLists.txt index e002025..b56433a 100644 --- a/k4Reco/CMakeLists.txt +++ b/k4Reco/CMakeLists.txt @@ -20,6 +20,11 @@ limitations under the License. set(sources DDPlanarDigi/components/DDPlanarDigi.cpp Overlay/components/OverlayTiming.cpp + Overlay/components/OverlayTimingRandomMix.cpp + CaloDigi/src/RealisticCaloDigiScinPpd.cc + CaloDigi/src/RealisticCaloDigiSilicon.cc + CaloDigi/src/RealisticCaloRecoScinPpd.cc + CaloDigi/src/RealisticCaloRecoSilicon.cc GaudiLumiCalClusterer/components/GaudiLumiCalClusterer.cpp GaudiLumiCalClusterer/src/LumiCalClusterer.cpp GaudiLumiCalClusterer/src/LumiCalClusterer_getCalHits.cpp diff --git a/k4Reco/CaloDigi/README.md b/k4Reco/CaloDigi/README.md new file mode 100644 index 0000000..e69de29 diff --git a/k4Reco/CaloDigi/include/CalorimeterHitType.h b/k4Reco/CaloDigi/include/CalorimeterHitType.h new file mode 100644 index 0000000..614a104 --- /dev/null +++ b/k4Reco/CaloDigi/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 CaloType { em = 0, had = 1, muon = 2 }; + + /** calo ids - specific to ILD */ + enum CaloID { unknown = 0, ecal = 1, hcal = 2, yoke = 3, lcal = 4, lhcal = 5, bcal = 6 }; + + /** calo layout / subdetector */ + enum 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/k4Reco/CaloDigi/include/FilterDoubleLayerHits.h b/k4Reco/CaloDigi/include/FilterDoubleLayerHits.h new file mode 100644 index 0000000..5a109c0 --- /dev/null +++ b/k4Reco/CaloDigi/include/FilterDoubleLayerHits.h @@ -0,0 +1,95 @@ +#ifndef FilterDoubleLayerHits_h +#define FilterDoubleLayerHits_h 1 + +#include +#include "k4Interface/IGeoSvc.h" +#include + +#include "DDRec/SurfaceManager.h" +#include +#include + +#include + +#include + + + +/** Utility processor that removes tracker hits in double layers if they don't have + * a corresponding close-by hit in the other sublayer. + * Pairs of considered layers are configurable and extracted from the cellID word. + * Works for all four lcio hit classes. + * + * @parameter InputCollection name of the hit collection with (Sim)TrackerHits/(Sim)CalorimeterHits + * @parameter OutputCollections ( ColName StartLayer EndLayer ) + * + * @author N. Bartosik, INFN Torino, S. Ferraro + * @date 17 June 2020 + * @version $Id: $ + */ + +struct FilterDoubleLayerHits : public k4FWCore::Transformer { + +protected: + + static const size_t NHITS_MAX = 10000000; + + struct SensorPosition{ + unsigned int layer; + unsigned int side; + unsigned int ladder; + unsigned int module; + + bool operator<(const SensorPosition& rhs) const { + return std::tie(layer, side, ladder, module) < std::tie(rhs.layer, rhs.side, rhs.ladder, rhs.module); + } + }; + + /// Double layer cut struct + struct DoubleLayerCut{ + unsigned int layer0 ; + unsigned int layer1 ; + double dPhi_max ; + double dTheta_max ; + }; + + + public: + FilterDoubleLayerHits(const std::string& name, ISvcLocator* svcLoc) ; + + /** Called at the begin of the job before anything is read. + * Use to initialize the processor, e.g. book histograms. + */ + virtual StatusCode initialize() ; + + /** Called for every event - the working horse. + */ + edm4hep::TrackerHitPlaneCollection operator()(const edm4hep::TrackerHitPlaneCollection& inputTrackerHitCollection) const override; + + /** Called after data processing for clean up. + */ + virtual StatusCode finalize(); + + + protected: + dd4hep::rec::Vector2D globalToLocal(long int cellID, const dd4hep::rec::Vector3D& posGlobal, dd4hep::rec::ISurface** surf) const; + + Gaudi::Property m_subDetName{this, "SubDetectorName", "Vertex", "Name of sub detector"}; + Gaudi::Property m_fillHistos{this, "FillHistograms", false, "Whether to fill diagnostic histograms"}; + Gaudi::Property m_dtMax{this, "DeltaTimeMax", -1.0, "Maximum time difference between hits in a doublet [ns]"}; + Gaudi::Property m_encodingStringVariable{this, "EncodingStringParameterName", "GlobalTrackerReadoutID", "The name of the DD4hep constant that contains the Encoding string for tracking detectors"}; + Gaudi::Property> m_dlCutConfigs{this, "DoubleLayerCuts" , {"0", "1", "0.5", "0.05"}, "Layer IDs and angular cuts [mrad] to be applied: "}; + + ////Double layer cuts configuration + std::vector m_dlCuts {}; + ////Surface map for getting local hit positions at sensor surface + const dd4hep::rec::SurfaceMap* m_map {nullptr}; + ////Monitoring histograms + std::map m_histos {}; + // GeoSvc + SmartIF m_geoSvc; + SmartIF m_histSvc; +} ; + +#endif \ No newline at end of file diff --git a/k4Reco/CaloDigi/include/RealisticCaloDigi.h b/k4Reco/CaloDigi/include/RealisticCaloDigi.h new file mode 100644 index 0000000..be2bc16 --- /dev/null +++ b/k4Reco/CaloDigi/include/RealisticCaloDigi.h @@ -0,0 +1,118 @@ +#ifndef REALISTICCALODIGI_H +#define REALISTICCALODIGI_H 1 + +#include +#include +#include +#include +#include +#include + +#include "k4Interface/IGeoSvc.h" +#include "k4Interface/IUniqueIDGenSvc.h" + +#include "TRandom2.h" + +#include +#include +#include +#include +#include + + +/** === RealisticCaloDigi Processor ===
+ Digitisation of calorimeter hits + e.g. timing, dead cells, miscalibrations + this is virtual class, technology-blind + technology-specific classes can inherit from this one + D. Jeans 02/2016, rewrite of parts of ILDCaloDigi, DDCaloDigi + R. Ete 11/2020, rewrite of charge integration and extension of timing treatment + */ + +struct RealisticCaloDigi : k4FWCore::MultiTransformer< + std::tuple( + const edm4hep::SimCalorimeterHitCollection&, + const edm4hep::EventHeaderCollection&)> { + public: + RealisticCaloDigi(const std::string& name, ISvcLocator* svcLoc); + /** Called at the begin of the job before anything is read. + * Use to initialize the processor, e.g. book histograms. + */ + StatusCode initialize(); + + /** Called for every run. + */ + std::tuple operator()( + const edm4hep::SimCalorimeterHitCollection& inputSim, + const edm4hep::EventHeaderCollection& headers) const; + + /** Called after data processing for clean up. + */ + StatusCode finalize(); + + + protected: + + // energy scales we know about + enum { MIP, GEVDEP, NPE }; + // integration result types + using integr_res = std::pair; + using integr_res_opt = std::optional; + using integr_function = std::function; + + virtual float EnergyDigi(float energy, float event_correl_miscalib) const; + virtual integr_res_opt Integrate( const edm4hep::SimCalorimeterHit * hit ) const; + + integr_res_opt StandardIntegration( const edm4hep::SimCalorimeterHit * hit ) const ; + integr_res_opt ROCIntegration( const edm4hep::SimCalorimeterHit * hit ) const ; + float SmearTime(float time) const; + + // virtual methods to be be overloaded in tech-specific derived classes + virtual int getMyUnit() const = 0 ; + virtual float digitiseDetectorEnergy(float energy) const = 0 ; + virtual float convertEnergy( float energy, int inScale ) const = 0; // convert energy from input to output scale + + // timing + Gaudi::Property m_time_apply{this, "timingCut", 0, "Use hit times"}; + Gaudi::Property m_time_correctForPropagation{this, "timingCorrectForPropagation", 0, "Correct hit times for propagation: radial distance/c"}; + Gaudi::Property m_time_windowMin{this, "timingWindowMin", -10.0f, "Time Window minimum time in ns"}; + Gaudi::Property m_time_windowMax{this, "timingWindowMax", 100.0f, "Time Window maximum time in ns"}; + Gaudi::Property m_integration_method{this, "integrationMethod", "Standard", "Energy integration and time calculation method. Options: Standard, ROC"}; + Gaudi::Property m_fast_shaper{this, "fastShaper", 0.f, "Fast shaper value. Unit in ns"}; + Gaudi::Property m_slow_shaper{this, "slowShaper", 0.f, "Slow shaper value. Unit in ns"}; + Gaudi::Property m_time_resol{this, "timingResolution", 0.f, "Time resolution to apply (gaussian smearing). Unit in ns"}; + // additional digi effects + Gaudi::Property m_calib_mip{this, "calibration_mip", 1.0e-4f, "Average G4 deposited energy by MIP for calibration"}; + Gaudi::Property m_misCalib_uncorrel{this, "miscalibration_uncorrel", 0.0f, "Uncorrelated random Gaussian miscalibration (as a fraction: 1.0 = 100%)"}; + Gaudi::Property m_misCalib_correl{this, "miscalibration_correl", 0.0f, "Correlated random Gaussian miscalibration (as a fraction: 1.0 = 100%)"}; + Gaudi::Property m_deadCell_fraction{this, "deadCell_fraction", 0.0f, "Random dead cell fraction (as a fraction: 0->1)"}; + // simple model of electronics properties + Gaudi::Property m_elec_noiseMip{this, "elec_noise_mip", 0.0f, "Typical electronics noise (in MIP units)"}; + Gaudi::Property m_elec_rangeMip{this, "elec_range_mip", 2500.0f, "Maximum of dynamic range of electronics (in MIPs)"}; + // code for layer info for cellID decoder + Gaudi::Property m_encodingStringVariable{this, "EncodingStringParameterName", "GlobalTrackerReadoutID", "The name of the DD4hep constant that contains the Encoding string for tracking detectors"}; + // energy threshold + Gaudi::Property m_threshold_value{this, "threshold", 0.5f, "Threshold for Hit"}; + Gaudi::Property m_threshold_unit{this, "thresholdUnit", std::string("MIP"), "Unit for threshold. Can be \"GeV\", \"MIP\" or \"px\". MIP and px need properly set calibration constants"}; + // id parameters + Gaudi::Property m_calo_type {this, "CaloType", "em", "Calorimeter Type: em, had, mu"}; + Gaudi::Property m_calo_id {this, "CaloID", "ecal", "Calorimeter ID: ecal, hcal, yoke, lcal, lhcal, bcal"}; + Gaudi::Property m_calo_layout {this, "CaloLayout", "barrel", "Calorimeter Layout: barrel, endcap, ring, plug"}; + + + int m_threshold_iunit{}; + inline static thread_local TRandom2 m_engine; + SmartIF m_geoSvc; + SmartIF m_uidSvc; + + + integr_function m_integr_function{}; + +} ; + +#endif + + + diff --git a/k4Reco/CaloDigi/include/RealisticCaloDigiScinPpd.h b/k4Reco/CaloDigi/include/RealisticCaloDigiScinPpd.h new file mode 100644 index 0000000..8a509ca --- /dev/null +++ b/k4Reco/CaloDigi/include/RealisticCaloDigiScinPpd.h @@ -0,0 +1,31 @@ +#ifndef DIGITIZER_DDCCALODIGISCINT_H +#define DIGITIZER_DDCCALODIGISCINT_H 1 + +#include "RealisticCaloDigi.h" + + +/** === RealisticCaloDigiScinPpd Processor ===
+ realistic digitisation of scint+PPD (SiPM, MPPC) calorimeter hits + D.Jeans 02/2016. +*/ + +struct RealisticCaloDigiScinPpd : public RealisticCaloDigi { + + public: + RealisticCaloDigiScinPpd(const std::string& name, ISvcLocator* svcLoc); + + protected: + int getMyUnit() const {return NPE;} + float digitiseDetectorEnergy(float energy) const ; // apply scin+PPD specific effects + float convertEnergy( float energy, int inputUnit ) const; // convert energy from input to output scale + + Gaudi::Property m_PPD_pe_per_mip{this, "ppd_mipPe", 10.0f, "# Photo-electrons per MIP (scintillator): used to Poisson smear #PEs if >0"}; + Gaudi::Property m_PPD_n_pixels{this, "ppd_npix", 10000, "Total number of MPPC/SiPM pixels for implementation of saturation effect"}; + Gaudi::Property m_misCalibNpix{this, "ppd_npix_uncert", 0.05f, "Fractional uncertainty of effective total number of MPPC/SiPM pixels"}; + Gaudi::Property m_pixSpread{this, "ppd_pix_spread", 0.05f, "Variation of PPD pixel signal (as a fraction: 0.01=1%)"}; +} ; + +#endif + + + diff --git a/k4Reco/CaloDigi/include/RealisticCaloDigiSilicon.h b/k4Reco/CaloDigi/include/RealisticCaloDigiSilicon.h new file mode 100644 index 0000000..9954e0a --- /dev/null +++ b/k4Reco/CaloDigi/include/RealisticCaloDigiSilicon.h @@ -0,0 +1,27 @@ +#ifndef DIGITIZER_REALISTICCALODIGISILICON_H +#define DIGITIZER_REALISTICCALODIGISILICON_H 1 + +#include "RealisticCaloDigi.h" + +/** === RealisticCaloDigiSilicon Processor ===
+ realistic digitisation of silicon calorimeter hits + D.Jeans 02/2016. +*/ + +struct RealisticCaloDigiSilicon : public RealisticCaloDigi { + + public: + RealisticCaloDigiSilicon(const std::string& name, ISvcLocator* svcLoc) ; + + protected: + int getMyUnit() const {return MIP;} + float convertEnergy( float energy, int inputUnit )const; // convert energy from input to output (MIP) scale + float digitiseDetectorEnergy(float energy)const; // apply silicon-specific realistic digitisation + + Gaudi::Property m_ehEnergy{this, "silicon_pairEnergy" , 3.6f, "energy required to create e-h pair in silicon (in eV)"}; +} ; + +#endif + + + diff --git a/k4Reco/CaloDigi/include/RealisticCaloReco.h b/k4Reco/CaloDigi/include/RealisticCaloReco.h new file mode 100644 index 0000000..99d83f1 --- /dev/null +++ b/k4Reco/CaloDigi/include/RealisticCaloReco.h @@ -0,0 +1,72 @@ +#ifndef DIGITIZER_REALISTICCALORECO_H +#define DIGITIZER_REALISTICCALORECO_H 1 + +#include + +#include +#include +#include + +#include "CalorimeterHitType.h" +#include "k4Interface/IGeoSvc.h" + +#include +#include + + +/** === RealisticCaloReco Processor ===
+ realistic reconstruction of calorimeter hits + e.g. apply sampling fraction correction + virtual class, technology indenpendent + D.Jeans 02/2016. + + 24 March 2016: removed gap corrections - to be put into separate processor + changed relations: now keep relation between reconstructed and simulated hits. +*/ + + +struct RealisticCaloReco : k4FWCore::MultiTransformer( + const edm4hep::CaloHitSimCaloHitLinkCollection&)> { + + public: + RealisticCaloReco(const std::string& name, ISvcLocator* svcLoc); + /** Called at the begin of the job before anything is read. + * Use to initialize the processor, e.g. book histograms. + */ + StatusCode initialize(); + + /** Called for every run. + */ + std::tuple operator()( + const edm4hep::CaloHitSimCaloHitLinkCollection& inputLinks) const; + + /** Called after data processing for clean up. + */ + StatusCode finalize(); + + + + protected: + + float getLayerCalib( int ilayer ) const; + virtual float reconstructEnergy(const edm4hep::CalorimeterHit* hit, int layer) const = 0; // to be overloaded, technology-specific + + // parameters + // Grouping of calo layers + Gaudi::Property> m_calLayers{this, "calibration_layergroups", {}, "Grouping of calo layers"}; + // Calibration coefficients for layers groups + Gaudi::Property> m_calibrCoeff{this, "calibration_factorsMipGev", {}, "Calibration coefficients (MIP->shower GeV) of layers groups"}; + // Cell ID layer string + Gaudi::Property m_encodingStringVariable{this, "EncodingStringParameterName", "GlobalTrackerReadoutID", "The name of the DD4hep constant that contains the Encoding string for tracking detectors"}; + + SmartIF m_geoSvc; + +} ; + +#endif + + + diff --git a/k4Reco/CaloDigi/include/RealisticCaloRecoScinPpd.h b/k4Reco/CaloDigi/include/RealisticCaloRecoScinPpd.h new file mode 100644 index 0000000..3c4707b --- /dev/null +++ b/k4Reco/CaloDigi/include/RealisticCaloRecoScinPpd.h @@ -0,0 +1,23 @@ +#ifndef REALISTICCALORECOSCINPPD_H +#define REALISTICCALORECOSCINPPD_H 1 + +#include "RealisticCaloReco.h" + +/** === RealisticCaloRecoSilicon Processor ===
+ realistic reconstruction of scint+PPD calorimeter hits + D.Jeans 02/2016. +*/ + +struct RealisticCaloRecoScinPpd final : RealisticCaloReco { + public: + RealisticCaloRecoScinPpd(const std::string& name, ISvcLocator* svcLoc); + + protected: + float reconstructEnergy(const edm4hep::CalorimeterHit* hit, int layer) const override; + + Gaudi::Property m_PPD_pe_per_mip{this, "ppd_mipPe", 10.0f, "# Photo-electrons per MIP (scintillator): used to Poisson smear #PEs if >0"}; + Gaudi::Property m_PPD_n_pixels{this, "ppd_npix", 10000, "Total number of MPPC/SiPM pixels for implementation of saturation effect"}; + +} ; + +#endif diff --git a/k4Reco/CaloDigi/include/RealisticCaloRecoSilicon.h b/k4Reco/CaloDigi/include/RealisticCaloRecoSilicon.h new file mode 100644 index 0000000..682bd30 --- /dev/null +++ b/k4Reco/CaloDigi/include/RealisticCaloRecoSilicon.h @@ -0,0 +1,23 @@ +#ifndef REALISTICCALORECOSILICON_H +#define REALISTICCALORECOSILICON_H 1 + +#include "RealisticCaloReco.h" + +/** === RealisticCaloRecoSilicon Processor ===
+ realistic reconstruction of silicon calorimeter hits + D.Jeans 02/2016. + + 24 March 2016: removed gap corrections - to be put into separate processor + +*/ + +struct RealisticCaloRecoSilicon final : RealisticCaloReco { + + public: + RealisticCaloRecoSilicon(const std::string& name, ISvcLocator* svcLoc); + + protected: + float reconstructEnergy(const edm4hep::CalorimeterHit* hit, int layer) const override; +} ; + +#endif diff --git a/k4Reco/CaloDigi/src/CalorimeterHitType.cc b/k4Reco/CaloDigi/src/CalorimeterHitType.cc new file mode 100644 index 0000000..610da06 --- /dev/null +++ b/k4Reco/CaloDigi/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/k4Reco/CaloDigi/src/FilterDoubleLayerHits.cc b/k4Reco/CaloDigi/src/FilterDoubleLayerHits.cc new file mode 100644 index 0000000..dd58cce --- /dev/null +++ b/k4Reco/CaloDigi/src/FilterDoubleLayerHits.cc @@ -0,0 +1,293 @@ +#include "FilterDoubleLayerHits.h" +#include +#include +#include + +#include +#include + +#include + +#include "DD4hep/DD4hepUnits.h" + +#include +#include + +DECLARE_COMPONENT(FilterDoubleLayerHits) + +FilterDoubleLayerHits::FilterDoubleLayerHits(const std::string& name, ISvcLocator* svcLoc) : Transformer(name, svcLoc, + KeyValues("InputCollection", {"VXDTrackerHitPlanes"}), + KeyValues("OutputCollection", {"VXDTrackerHitPlanes_DLFiltered"})) {} + + +dd4hep::rec::Vector2D FilterDoubleLayerHits::globalToLocal(long int cellID, const dd4hep::rec::Vector3D& posGlobal, dd4hep::rec::ISurface** surfptr=nullptr) const{ + dd4hep::rec::ISurface* surf; + // Using directly the provided surface object if available + if (surfptr && *surfptr) { + surf = *surfptr; + } else { + // Finding the surface corresponding to the cellID + dd4hep::rec::SurfaceMap::const_iterator surfIt = m_map->find( cellID ); + if( surfIt == m_map->end() ){ + throw GaudiException(" FilterDoubleLayerHits::processEvent(): no surface found for cellID: " + cellID, "Fail", StatusCode::FAILURE); + } + surf = surfIt->second; + // Saving the surface object outside the function to be reused for the same cellID + if (surfptr) *surfptr = surf; + } + // Converting global position to local in [cm] + dd4hep::rec::Vector2D posLocal = surf->globalToLocal( dd4hep::mm * posGlobal ); + + return dd4hep::rec::Vector2D( posLocal.u() / dd4hep::mm, posLocal.v() / dd4hep::mm ); +} + + +StatusCode FilterDoubleLayerHits::initialize() { + m_geoSvc = serviceLocator()->service("GeoSvc"); + if (!m_geoSvc) { + error() << "Unable to retrieve the GeoSvc" << endmsg; + return StatusCode::FAILURE; + } + // Get Histogram and Data Services + SmartIF histSvc; + histSvc = serviceLocator()->service("THistSvc"); + + debug() << " init called " << endmsg; + + // Extracting double-layer cut configurations + m_dlCuts.resize( m_dlCutConfigs.size() / 4 ) ; + + unsigned i=0,index=0 ; + while( i < m_dlCutConfigs.size() ){ + m_dlCuts[index].layer0 = std::atoi( m_dlCutConfigs[ i++ ].c_str() ) ; + m_dlCuts[index].layer1 = std::atoi( m_dlCutConfigs[ i++ ].c_str() ) ; + m_dlCuts[index].dPhi_max = std::atof( m_dlCutConfigs[ i++ ].c_str() ) / 1e3 ; // converting mrad -> rad + m_dlCuts[index].dTheta_max = std::atof( m_dlCutConfigs[ i++ ].c_str() ) / 1e3 ; // converting mrad -> rad + ++index ; + } + + // Get the surface map from the SurfaceManager + dd4hep::Detector& theDetector = *(m_geoSvc->getDetector()); + dd4hep::rec::SurfaceManager& surfMan = *theDetector.extension(); + dd4hep::DetElement det = theDetector.detector( m_subDetName ); + + m_map = surfMan.map( det.name() ); + + if( !m_map ) { + throw GaudiException( " Could not find surface map for detector: " + m_subDetName + " in SurfaceManager", "Fail", StatusCode::FAILURE); + } + + // Booking diagnostic histograms for each configured cut + char hname[100]; + for(size_t iCut=0, nCuts=m_dlCuts.size() ; iCut layers{cut.layer0, cut.layer1}; + for (auto layer : layers) { + sprintf(hname, "h2_posUV_rejected_layer_%d", layer); + m_histos[ std::string(hname) ] = new TH2F( hname , ";U [mm]; V [mm]", 500, -100, 100, 1000, -200, 200 ); + } + } + + for (const auto& [name, histo] : m_histos) { + (void)m_histSvc->regHist("/histos/" + name, histo); + } + + // Printing the configured cut + debug() << iCut << ". layers: " << cut.layer0 << " >> " << cut.layer1 << "; dPhi: " + << cut.dPhi_max << " rad; dTheta: " << cut.dTheta_max << " rad" << endmsg; + } + + return StatusCode::SUCCESS; + +} + + +edm4hep::TrackerHitPlaneCollection FilterDoubleLayerHits::operator()(const edm4hep::TrackerHitPlaneCollection& inputTrackerHitCollection) const{ + std::string initString; + initString = m_geoSvc->constantAsString(m_encodingStringVariable.value()); + dd4hep::DDSegmentation::BitFieldCoder bitFieldCoder(initString); // check! + + //---- create the output collection + edm4hep::TrackerHitPlaneCollection outCol; + outCol.setSubsetCollection(); + + // Set acceptance flags for all hits to FALSE + const size_t nHit = inputTrackerHitCollection.size(); + bool hitAccepted[NHITS_MAX]; ////Array of flags for hits to be accepted + memset(&hitAccepted, false, nHit); + + ////Map of vectors of hits grouped by position in the detector + std::map > hitsGrouped{}; + + // Splitting hits by sensor ids for faster association + for (size_t iHit = 0; iHit < nHit ; iHit++) { + + edm4hep::TrackerHitPlane h = inputTrackerHitCollection.at( iHit ); + + unsigned int layerID = bitFieldCoder.get(h.getCellID(), "layer"); + unsigned int sideID = bitFieldCoder.get(h.getCellID(), "side"); + unsigned int ladderID = bitFieldCoder.get(h.getCellID(), "module"); + unsigned int moduleID = bitFieldCoder.get(h.getCellID(), "sensor"); + + SensorPosition sensPos = {layerID, sideID, ladderID, moduleID}; + if (hitsGrouped.find(sensPos) == hitsGrouped.end()) { + hitsGrouped[sensPos] = std::vector(); + hitsGrouped[sensPos].reserve(nHit); + } + hitsGrouped[sensPos].push_back(iHit); + } + + //---- loop over hits + for (size_t iHit = 0; iHit < nHit ; iHit++) { + + // Skipping hits that are already accepted + if (hitAccepted[iHit]) continue; + + edm4hep::TrackerHitPlane h = inputTrackerHitCollection.at( iHit ); + + unsigned int layerID = bitFieldCoder.get(h.getCellID(), "layer"); + unsigned int sideID = bitFieldCoder.get(h.getCellID(), "side"); + unsigned int ladderID = bitFieldCoder.get(h.getCellID(), "module"); + unsigned int moduleID = bitFieldCoder.get(h.getCellID(), "sensor"); + debug() << " Checking 1st hit " << iHit << " / " << nHit << " at layer: " << layerID << " ladder: " << ladderID << " module: " << moduleID << endmsg ; + + const SensorPosition sensPos = {layerID, sideID, ladderID, moduleID}; + + // Checking if the hit is at the inner double layer to be filtered + const DoubleLayerCut* dlCut(0); + for (int iCut=0, nCuts=m_dlCuts.size(); iCutlayer1) continue; + + // Getting local and global hit positions + dd4hep::rec::Vector3D posGlobal( h.getPosition().x, h.getPosition().y, h.getPosition().z ); + dd4hep::rec::Vector2D posLocal = globalToLocal( h.getCellID(), posGlobal ); + + // Setting the values for closest hits + double dR_min(999.0); + double dU_closest(0.0); + double dTheta_closest(0.0); + double dPhi_closest(0.0); + double dt_closest(0.0); + + // Looking for the compliment hits in the 2nd sublayer + size_t nCompatibleHits(0); + SensorPosition sensPos2 = sensPos; + sensPos2.layer = dlCut->layer1; + dd4hep::rec::ISurface* surf=nullptr; + // Checking if there are any hits in the corresponding sensor at the other sublayer + if (hitsGrouped.find(sensPos2) == hitsGrouped.end()) continue; + for (size_t iHit2 : hitsGrouped.at(sensPos2)) { + edm4hep::TrackerHitPlane h2 = inputTrackerHitCollection.at( iHit2 ); + unsigned int layerID2 = bitFieldCoder.get(h2.getCellID(), "layer"); + + // Checking whether hit is in the time acceptance window + double dt = h2.getTime() - h.getTime(); + if (m_dtMax >= 0.0 && std::fabs(dt) > m_dtMax) continue; + + // Getting the local and global hit positions + dd4hep::rec::Vector3D posGlobal2( h2.getPosition().x, h2.getPosition().y, h2.getPosition().z ); + dd4hep::rec::Vector2D posLocal2 = globalToLocal( h2.getCellID(), posGlobal2, &surf ); + + // Checking whether hit is close enough to the 1st one + double dU = posLocal2.u() - posLocal.u(); + double dTheta = posGlobal2.theta() - posGlobal.theta(); + double dPhi = std::fabs(posGlobal2.phi() - posGlobal.phi()); + if (dPhi > dd4hep::pi) dPhi = dd4hep::twopi - dPhi; + double dR = sqrt(dPhi*dPhi + dTheta*dTheta); + debug() << " Checking 2nd hit at layer: " << layerID2 << "; dPhi: " << dPhi << "; dTheta: " << dTheta << endmsg; + + // Updating the minimal values + if (dR < dR_min) { + dR_min = dR; + dU_closest = dU; + dTheta_closest = dTheta; + dPhi_closest = dPhi; + dt_closest = dt; + } + + // Skipping if the hit is outside the cut window + if (std::fabs(dPhi) > dlCut->dPhi_max) continue; + if (std::fabs(dTheta) > dlCut->dTheta_max) continue; + + nCompatibleHits++; + hitAccepted[iHit2] = true; + debug() << " Accepted 2nd hit at layer: " << layerID2 << "; dPhi: " << dPhi << "; dTheta: " << dTheta << endmsg; + } + // Filling diagnostic histograms + if (m_fillHistos && dR_min < 998) { + char hname[100]; + sprintf(hname, "h_dU_layers_%d_%d", dlCut->layer0, dlCut->layer1); + m_histos.find(hname)->second->Fill(dU_closest); + sprintf(hname, "h2_dU_dPhi_layers_%d_%d", dlCut->layer0, dlCut->layer1); + m_histos.find(hname)->second->Fill(dU_closest, dPhi_closest*1e3); + sprintf(hname, "h_dTheta_layers_%d_%d", dlCut->layer0, dlCut->layer1); + m_histos.find(hname)->second->Fill(dTheta_closest*1e3); + sprintf(hname, "h_dPhi_layers_%d_%d", dlCut->layer0, dlCut->layer1); + m_histos.find(hname)->second->Fill(dPhi_closest*1e3); + sprintf(hname, "h_dt_layers_%d_%d", dlCut->layer0, dlCut->layer1); + m_histos.find(hname)->second->Fill(dt_closest); + } + + // Accepting the first hit if it has at least one compatible pair + if (nCompatibleHits > 0) { + hitAccepted[iHit] = true; + debug() << " Accepted 1st hit at layer: " << layerID << endmsg; + } + } + + // Adding accepted hits to the output collection + size_t nHitsAccepted(0); + for (size_t iHit = 0; iHit < nHit; iHit++) { + if (!hitAccepted[iHit]) { + // Filling the positions of rejected hits + if (m_fillHistos) { + edm4hep::TrackerHitPlane h = inputTrackerHitCollection.at( iHit ); + unsigned int layerID = bitFieldCoder.get(h.getCellID(), "layer"); + + // Getting local hit position + dd4hep::rec::Vector3D posGlobal( h.getPosition().x, h.getPosition().y, h.getPosition().z ); + dd4hep::rec::Vector2D posLocal = globalToLocal( h.getCellID(), posGlobal ); + + char hname[100]; + sprintf(hname, "h2_posUV_rejected_layer_%d", layerID); + m_histos.find(hname)->second->Fill(posLocal.u(), posLocal.v()); + } + continue; + } + outCol.push_back( inputTrackerHitCollection.at( iHit ) ); + nHitsAccepted++; + } + info() << " " << nHitsAccepted << " hits added to collection." << endmsg; + + return outCol; +} + +StatusCode FilterDoubleLayerHits::finalize(){ + return StatusCode::SUCCESS; +} \ No newline at end of file diff --git a/k4Reco/CaloDigi/src/RealisticCaloDigi.cc b/k4Reco/CaloDigi/src/RealisticCaloDigi.cc new file mode 100644 index 0000000..4ec1458 --- /dev/null +++ b/k4Reco/CaloDigi/src/RealisticCaloDigi.cc @@ -0,0 +1,302 @@ +// Calorimeter digitiser +#include "RealisticCaloDigi.h" + +#include +#include +#include +#include + +#include + +#include + +#include +#include +#include +#include +#include +#include + +#include "CLHEP/Units/PhysicalConstants.h" + + +using namespace std; +using namespace std::placeholders; + +struct MCC { + float energy {0.f}; + float time {0.f}; +}; + +RealisticCaloDigi::RealisticCaloDigi(const std::string& name, ISvcLocator* svcLoc) : MultiTransformer(name, svcLoc, + { KeyValues("inputHitCollections", {"SimCalorimeterHits"}), + KeyValues("inputHeaderCollections", {"EventHeader"}) }, + { KeyValues("outputHitCollections", {"CalorimeterHits"}), + KeyValues("outputRelationCollections", {"CaloHitLinks"}) }) {} + +StatusCode RealisticCaloDigi::initialize() { + m_geoSvc = serviceLocator()->service("GeoSvc"); + if (!m_geoSvc) { + error() << "Unable to retrieve the GeoSvc" << endmsg; + return StatusCode::FAILURE; + } + m_uidSvc = service("UniqueIDGenSvc", true); + if (!m_uidSvc) { + error() << "Unable to get UniqueIDGenSvc" << endmsg; + } + + + // unit in which threshold is specified + if (m_threshold_unit.value().compare("MIP") == 0){ + m_threshold_iunit=MIP; + } else if (m_threshold_unit.value().compare("GeV") == 0){ + m_threshold_iunit=GEVDEP; + } else if (m_threshold_unit.value().compare("px") == 0){ + m_threshold_iunit=NPE; + } else { + error() << "could not identify threshold unit. Please use \"GeV\", \"MIP\" or \"px\"! Aborting." << endmsg; + } + + // convert the threshold to the approriate units (i.e. MIP for silicon, NPE for scint) + m_threshold_value = convertEnergy( m_threshold_value, m_threshold_iunit ); + + // deal with timing calculations + std::map integrations = { + {"Standard", std::bind(&RealisticCaloDigi::StandardIntegration, this, _1)}, + {"ROC", std::bind(&RealisticCaloDigi::ROCIntegration, this, _1)} + }; + auto findIter = integrations.find( m_integration_method ) ; + if(integrations.end() == findIter) { + error() << "Could not guess timing calculation method!" << endmsg; + error() << "Available are: Standard, ROC. Provided: " << m_integration_method << endmsg; + error() << "Aborting..." << endmsg; + } + m_integr_function = findIter->second; + + // check if parameters are correctly set for the ROC integration + if("ROC" == m_integration_method) { + if(m_fast_shaper == 0.0f || m_slow_shaper == 0.0f) { + error() << "Fast/slow shaper parameter(s) not set. Required for ROC integration!" << endmsg; + error() << "Aborting..." << endmsg; + } + } + + return StatusCode::SUCCESS; +} + + + +std::tuple RealisticCaloDigi::operator()( + const edm4hep::SimCalorimeterHitCollection& inputSim, + const edm4hep::EventHeaderCollection& headers) const { + auto seed = m_uidSvc->getUniqueID(headers[0].getEventNumber(), headers[0].getRunNumber(), this->name()); + debug() << "Using seed " << seed << " for event " << headers[0].getEventNumber() << " and run " + << headers[0].getRunNumber() << endmsg; + m_engine.SetSeed(seed); + + // decide on this event's correlated miscalibration + float event_correl_miscalib = ( m_misCalib_correl>0 ) ? m_engine.Gaus(1.0, m_misCalib_correl) : 0; + + edm4hep::CalorimeterHitCollection newcol; + edm4hep::CaloHitSimCaloHitLinkCollection relcol; + + CHT::CaloType cht_type = caloTypeFromString(m_calo_type); + CHT::CaloID cht_id = caloIDFromString(m_calo_id); + CHT::Layout cht_lay = layoutFromString(m_calo_layout); + + std::string initString; + initString = m_geoSvc->constantAsString(m_encodingStringVariable.value()); + dd4hep::DDSegmentation::BitFieldCoder bitFieldCoder(initString); // check! + + debug() << "Number of elements = " << inputSim.size() << endmsg; + // loop over input hits + for (int j=0; j < inputSim.size(); ++j) { + edm4hep::SimCalorimeterHit simhit0 = inputSim.at( j ); + edm4hep::SimCalorimeterHit *simhit = &simhit0; + + // deal with energy integration and timing aspects + auto integrationResult = Integrate(simhit); + if( ! integrationResult.has_value() ) { + continue; + } + float time = integrationResult.value().first; + float energyDep = integrationResult.value().second; + // apply extra energy digitisation onto the energy + float energyDig = EnergyDigi(energyDep, event_correl_miscalib); + + if (energyDig > m_threshold_value) { // write out this hit + edm4hep::MutableCalorimeterHit newhit = newcol.create(); + newhit.setCellID( simhit->getCellID() ); + newhit.setTime( time ); + newhit.setPosition( simhit->getPosition() ); + newhit.setEnergy( energyDig ); + + int layer = bitFieldCoder.get(simhit->getCellID(), "layer"); + newhit.setType( CHT( cht_type, cht_id, cht_lay, layer ) ); + + debug() << "orig/new hit energy: " << simhit->getEnergy() << " " << newhit.getEnergy() << endmsg; + + edm4hep::MutableCaloHitSimCaloHitLink rel = relcol.create(); + rel.setTo(simhit0); + rel.setFrom(newhit); + rel.setWeight(1.0); + + } // theshold + } // input hits + + return std::make_tuple(std::move(newcol), std::move(relcol)); +} + +//------------------------------------------------------------------------------ + +StatusCode RealisticCaloDigi::finalize(){ + return StatusCode::SUCCESS; +} + +//------------------------------------------------------------------------------ + +RealisticCaloDigi::integr_res_opt RealisticCaloDigi::Integrate( const edm4hep::SimCalorimeterHit * hit ) const { + return m_integr_function(hit); +} + +//------------------------------------------------------------------------------ + +float RealisticCaloDigi::EnergyDigi(float energy, float event_correl_miscalib) const{ + // some extra digi effects + // controlled by _applyDigi = 0 (none), 1 (apply) + // input parameters: hit energy ( in any unit: effects are all relative ) + // returns energy ( in units determined by the overloaded digitiseDetectorEnergy ) + + float e_out(energy); + e_out = digitiseDetectorEnergy(energy); // this is an overloaded method, provides energy in technology-dependent units + + // the following make only relative changes to the energy + + // random miscalib, uncorrelated in cells + if (m_misCalib_uncorrel>0) { + float miscal(0); + miscal = m_engine.Gaus(1.0, m_misCalib_uncorrel); + e_out*=miscal; + } + + // random miscalib, correlated across cells in one event + if (m_misCalib_correl>0) e_out*=event_correl_miscalib; + + float oneMipInMyUnits = convertEnergy( 1.0, MIP ); + // limited electronics dynamic range + if ( m_elec_rangeMip > 0 ) e_out = std::min ( e_out, m_elec_rangeMip*oneMipInMyUnits ); + // add electronics noise + if ( m_elec_noiseMip > 0 ) { + e_out += m_engine.Gaus(0, m_elec_noiseMip*oneMipInMyUnits); + } + + // random cell kill + if (m_deadCell_fraction>0) { + if (m_engine.Uniform(0., 1.) < m_deadCell_fraction ) e_out=0; + } + return e_out; +} + +//------------------------------------------------------------------------------ + +RealisticCaloDigi::integr_res_opt RealisticCaloDigi::StandardIntegration( const edm4hep::SimCalorimeterHit * hit ) const { + // apply timing cuts on simhit contributions + // outputs a (time,energy) pair + float timeCorrection(0); + if ( m_time_correctForPropagation ) { // time of flight from IP to this point + float r = pow(hit->getPosition().x,2) + pow(hit->getPosition().y,2) + pow(hit->getPosition().z,2); + timeCorrection = sqrt(r)/CLHEP::c_light; // [speed of light in mm/ns] + } + // this is Oskar's simple (and probably the most correct) method for treatment of timing + // - collect energy in some predefined time window around collision time (possibly corrected for TOF) + // - assign time of earliest contribution to hit + float energySum = 0; + float earliestTime=std::numeric_limits::max(); + for(edm4hep::CaloHitContribution contribution : hit->getContributions()){ // loop over all contributions + float timei = contribution.getTime(); //absolute hit timing of current subhit + float energyi = contribution.getEnergy(); //energy of current subhit + float relativetime = timei - timeCorrection; // wrt time of flight + if (relativetime>m_time_windowMin && relativetime m_time_windowMin && earliestTime < m_time_windowMax){ //accept this hit + return integr_res{SmearTime(earliestTime), energySum}; + } + return std::nullopt; +} + +//------------------------------------------------------------------------------ + +RealisticCaloDigi::integr_res_opt RealisticCaloDigi::ROCIntegration( const edm4hep::SimCalorimeterHit * hit ) const { + const unsigned int ncontrib = hit->contributions_size() ; + // Sort MC contribution by time + std::vector mcconts{ncontrib}; + for(int i=0; igetContributions().size();i++){ + mcconts[i].energy = hit->getContributions(i).getEnergy(); + mcconts[i].time = hit->getContributions(i).getTime(); + } + std::sort(mcconts.begin(), mcconts.end(), [](auto lhs, auto rhs){ + return (lhs.time < rhs.time); + }); + // Accumulate energy until threshold is reached. + // The first MC contriubtion after the threshold has been reached sets the hit time + bool passThreshold = false; + float epar=0.f, hitTime=0.f; + unsigned int thresholdIndex=0; + // First determine the hit time (hitTime) and the initial hit index + // at which we need to start the integration (thresholdIndex) + for(unsigned int i=0; i m_threshold_value ) { + hitTime = timej; + passThreshold = true ; + break; + } + } + if(passThreshold) { + break; + } + } + // check hit time + const float thresholdTime = mcconts[thresholdIndex].time; + if( not (thresholdTime>m_time_windowMin && thresholdTime0.f ? time + m_engine.Gaus(0, m_time_resol) : time; +} + diff --git a/k4Reco/CaloDigi/src/RealisticCaloDigiScinPpd.cc b/k4Reco/CaloDigi/src/RealisticCaloDigiScinPpd.cc new file mode 100644 index 0000000..f93d467 --- /dev/null +++ b/k4Reco/CaloDigi/src/RealisticCaloDigiScinPpd.cc @@ -0,0 +1,49 @@ +// Calorimeter digitiser for the IDC ECAL and HCAL +// For other detectors/models SimpleCaloDigi should be used + +#include "RealisticCaloDigiScinPpd.h" + +#include +#include +#include +#include + +#include "CLHEP/Random/RandGauss.h" +#include "CLHEP/Random/RandBinomial.h" + +using namespace std; + +DECLARE_COMPONENT(RealisticCaloDigiScinPpd) + +RealisticCaloDigiScinPpd::RealisticCaloDigiScinPpd(const std::string& name, ISvcLocator* svcLoc) + : RealisticCaloDigi(name, svcLoc) {} + +float RealisticCaloDigiScinPpd::convertEnergy( float energy, int inUnit ) const { // convert energy from input to output scale (NPE) + if ( inUnit==NPE ) return energy; + else if ( inUnit==MIP ) return m_PPD_pe_per_mip*energy; + else if ( inUnit==GEVDEP ) return m_PPD_pe_per_mip*energy/m_calib_mip; + + throw std::runtime_error("RealisticCaloDigiScinPpd::convertEnergy - unknown unit " + std::to_string(inUnit)); +} + +float RealisticCaloDigiScinPpd::digitiseDetectorEnergy(float energy) const { + // input energy in deposited GeV + // output in npe + float npe = energy*m_PPD_pe_per_mip/m_calib_mip; // convert to pe scale + + if (m_PPD_n_pixels>0){ + // apply average sipm saturation behaviour + npe = m_PPD_n_pixels*(1.0 - exp( -npe/m_PPD_n_pixels ) ); + //apply binomial smearing + float p = npe/m_PPD_n_pixels; // fraction of hit pixels on SiPM + npe = m_engine.Binomial(m_PPD_n_pixels, p); //npe now quantised to integer pixels + + if (m_pixSpread>0) { + // variations in pixel capacitance + npe *= m_engine.Gaus(1, m_pixSpread/sqrt(npe) ); + } + } + + return npe; +} + diff --git a/k4Reco/CaloDigi/src/RealisticCaloDigiSilicon.cc b/k4Reco/CaloDigi/src/RealisticCaloDigiSilicon.cc new file mode 100644 index 0000000..77b73c5 --- /dev/null +++ b/k4Reco/CaloDigi/src/RealisticCaloDigiSilicon.cc @@ -0,0 +1,44 @@ +// Calorimeter digitiser for the IDC ECAL and HCAL +// For other detectors/models SimpleCaloDigi should be used +#include "RealisticCaloDigiSilicon.h" + +#include +#include +#include +#include + +#include "CLHEP/Random/RandPoisson.h" +#include "CLHEP/Random/RandGauss.h" + +using namespace std; + +DECLARE_COMPONENT(RealisticCaloDigiSilicon) + +RealisticCaloDigiSilicon::RealisticCaloDigiSilicon(const std::string& name, ISvcLocator* svcLoc) + : RealisticCaloDigi(name, svcLoc) {} + +float RealisticCaloDigiSilicon::convertEnergy( float energy, int inUnit ) const { // convert energy from input to output scale (MIP) + // converts input energy to MIP scale + if ( inUnit==MIP ) return energy; + else if ( inUnit==GEVDEP ) return energy/m_calib_mip; + + throw std::runtime_error("RealisticCaloDigiSilicon::convertEnergy - unknown unit " + std::to_string(inUnit)); +} + + +float RealisticCaloDigiSilicon::digitiseDetectorEnergy(float energy) const { + // applies extra digitisation to silicon hits + // input energy in deposited GeV + // output is MIP scale + float smeared_energy(energy); + if ( m_ehEnergy>0 ) { + // calculate #e-h pairs + float nehpairs = 1e9*energy/m_ehEnergy; // check units of energy! _ehEnergy is in eV, energy in GeV + // fluctuate it by Poisson (actually an overestimate: Fano factor actually makes it smaller, however even this overstimated effect is tiny for our purposes) + smeared_energy *= m_engine.Poisson( nehpairs )/nehpairs; + } + + return smeared_energy/m_calib_mip; // convert to MIP units +} + + diff --git a/k4Reco/CaloDigi/src/RealisticCaloReco.cc b/k4Reco/CaloDigi/src/RealisticCaloReco.cc new file mode 100644 index 0000000..213a75b --- /dev/null +++ b/k4Reco/CaloDigi/src/RealisticCaloReco.cc @@ -0,0 +1,92 @@ +#include "RealisticCaloReco.h" + +#include + +#include +#include +#include + +#include +#include +#include +#include +#include + +using namespace std; + +RealisticCaloReco::RealisticCaloReco(const std::string& name, ISvcLocator* svcLoc) : MultiTransformer(name, svcLoc, + { KeyValues("inputLinkCollections", {"CaloHitLinks"}) }, + { KeyValues("outputHitCollections", {"CalorimeterHitsRec"}), + KeyValues("outputRelationCollections", {"CaloHitLinksRec"}) }) {} + +StatusCode RealisticCaloReco::initialize() { + m_geoSvc = serviceLocator()->service("GeoSvc"); + if (!m_geoSvc) { + error() << "Unable to retrieve the GeoSvc" << endmsg; + return StatusCode::FAILURE; + } + + assert ( m_calibrCoeff.size()>0 ); + assert ( m_calibrCoeff.size() == m_calLayers.size() ); + + return StatusCode::SUCCESS; +} + +//----------------------------------------------------------------------------------------------- + +std::tuple RealisticCaloReco::operator()( + const edm4hep::CaloHitSimCaloHitLinkCollection& inputLinks) const { + // * Reading Collections of digitised calorimeter Hits * + std::string initString; + initString = m_geoSvc->constantAsString(m_encodingStringVariable.value()); + dd4hep::DDSegmentation::BitFieldCoder bitFieldCoder(initString); // check! + + edm4hep::CalorimeterHitCollection newcol; + edm4hep::CaloHitSimCaloHitLinkCollection relcol; + debug() << " number of elements = " << inputLinks.size() << endmsg; + + for (int j(0); j < inputLinks.size(); ++j) { + edm4hep::CaloHitSimCaloHitLink link = inputLinks.at( j ) ; + edm4hep::CalorimeterHit hit0 = link.getFrom(); + edm4hep::CalorimeterHit *hit = &hit0; + edm4hep::MutableCalorimeterHit calhit = newcol.create(); // make new hit + + int cellid = hit->getCellID(); + float energy = reconstructEnergy( hit, bitFieldCoder.get(cellid, "layer") ); // overloaded method, technology dependent + + calhit.setCellID(cellid); + calhit.setEnergy(energy); + calhit.setTime( hit->getTime() ); + calhit.setPosition( hit->getPosition() ); + calhit.setType( hit->getType() ); + + edm4hep::MutableCaloHitSimCaloHitLink newLink = relcol.create(); + newLink.setFrom(calhit); + newLink.setTo(link.getTo()); + newLink.setWeight(1.0); + } + + return std::make_tuple(std::move(newcol), std::move(relcol)); +} + +float RealisticCaloReco::getLayerCalib( int ilayer ) const{ + float calib_coeff = 0; + // retrieve calibration constants + // Fixed the following logic (DJeans, June 2016) + int min(0),max(0); + for (unsigned int k(0); k < m_calLayers.size(); ++k) { + if ( k > 0 ) min+=m_calLayers[k-1]; + max+=m_calLayers[k]; + if (ilayer >= min && ilayer < max) { + calib_coeff = m_calibrCoeff[k]; + break; + } + } + assert( calib_coeff>0 ); + return calib_coeff; +} + + +StatusCode RealisticCaloReco::finalize(){ return StatusCode::SUCCESS; } + diff --git a/k4Reco/CaloDigi/src/RealisticCaloRecoScinPpd.cc b/k4Reco/CaloDigi/src/RealisticCaloRecoScinPpd.cc new file mode 100644 index 0000000..e25e120 --- /dev/null +++ b/k4Reco/CaloDigi/src/RealisticCaloRecoScinPpd.cc @@ -0,0 +1,31 @@ +#include "RealisticCaloRecoScinPpd.h" +#include +#include + +DECLARE_COMPONENT(RealisticCaloRecoScinPpd) + +RealisticCaloRecoScinPpd::RealisticCaloRecoScinPpd(const std::string& name, ISvcLocator* svcLoc) : RealisticCaloReco(name, svcLoc) {} + +float RealisticCaloRecoScinPpd::reconstructEnergy(const edm4hep::CalorimeterHit* hit, int layer) const{ + // here the input energy should be in NPE + float energy = hit->getEnergy(); + + // first de-saturate PPD response + // this is the fraction of SiPM pixels fired above which a linear continuation of the saturation-reconstruction function is used. + // 0.95 of nPixel corresponds to a energy correction of factor ~3. + const float r = 0.95; + if (energy < r*m_PPD_n_pixels){ //current hit below linearisation threshold, reconstruct energy normally: + energy = -m_PPD_n_pixels * std::log ( 1. - ( energy / m_PPD_n_pixels ) ); + } else { //current hit is aove linearisation threshold, reconstruct using linear continuation function: + energy = 1/(1-r)*(energy-r*m_PPD_n_pixels)-m_PPD_n_pixels*std::log(1-r); + } + // then go back to MIP scale + energy/=m_PPD_pe_per_mip; + + // now correct for sampling fraction (calibration from MIP -> shower GeV) + energy *= getLayerCalib( layer ); + + return energy; +} + + diff --git a/k4Reco/CaloDigi/src/RealisticCaloRecoSilicon.cc b/k4Reco/CaloDigi/src/RealisticCaloRecoSilicon.cc new file mode 100644 index 0000000..09b56c9 --- /dev/null +++ b/k4Reco/CaloDigi/src/RealisticCaloRecoSilicon.cc @@ -0,0 +1,17 @@ +#include "RealisticCaloRecoSilicon.h" +#include +#include +#include + + +DECLARE_COMPONENT(RealisticCaloRecoSilicon) + +RealisticCaloRecoSilicon::RealisticCaloRecoSilicon(const std::string& name, ISvcLocator* svcLoc) : RealisticCaloReco(name, svcLoc) {} + +float RealisticCaloRecoSilicon::reconstructEnergy(const edm4hep::CalorimeterHit* hit, int layer) const{ + // here the input energy should be in MIPs + float energy = hit->getEnergy(); + // now correct for sampling fraction + energy *= getLayerCalib( layer ); + return energy; +} diff --git a/k4Reco/Overlay/components/OverlayTimingRandomMix.cpp b/k4Reco/Overlay/components/OverlayTimingRandomMix.cpp new file mode 100644 index 0000000..d0fbd92 --- /dev/null +++ b/k4Reco/Overlay/components/OverlayTimingRandomMix.cpp @@ -0,0 +1,471 @@ +/* + * 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 "OverlayTimingRandomMix.h" +#include + +#include "podio/FrameCategories.h" + +#include "edm4hep/CaloHitContributionCollection.h" +#include "edm4hep/Constants.h" +#include "edm4hep/EventHeaderCollection.h" +#include "edm4hep/MCParticleCollection.h" +#include "edm4hep/SimCalorimeterHitCollection.h" +#include "edm4hep/SimTrackerHitCollection.h" + +#include "k4FWCore/MetadataUtils.h" + +#include + +#include +#include +#include +#include +namespace fs = std::filesystem; + +std::vector filesInFolder(const std::string &folderPath) +{ + std::vector files; + try + { + for (const auto &entry : fs::directory_iterator(folderPath)) + { + if (fs::is_regular_file(entry.path()) && entry.path().extension() == ".root") + { + files.push_back(entry.path().string()); + } + } + } + catch (const std::exception &ex) + { + std::cerr << "Error: " << ex.what() << std::endl; + } + return files; +} + + +template inline float time_of_flight(const T& pos) { + // Returns the time of flight to the radius in ns + // Assumming positions in mm, then mm/m/s = 10^-3 s = 10^6 ns + return std::sqrt((pos[0] * pos[0]) + (pos[1] * pos[1]) + (pos[2] * pos[2])) / TMath::C() * 1e6; +} + +std::pair OverlayTimingRandomMix::define_time_windows(const std::string& collection_name) const { + try { + return {m_timeWindows.value().at(collection_name)[0], m_timeWindows.value().at(collection_name)[1]}; + } catch (const std::out_of_range& e) { + error() << "No time window defined for collection " << collection_name << endmsg; + throw e; + } +} + +StatusCode OverlayTimingRandomMix::initialize() { + m_uidSvc = service("UniqueIDGenSvc", true); + if (!m_uidSvc) { + error() << "Unable to get UniqueIDGenSvc" << endmsg; + } + + std::vector> inputFiles; + if (m_inputFileNames.value()[0][0].find("root") != std::string::npos) { + inputFiles = m_inputFileNames.value(); + } else { + for (auto group : m_inputFileNames.value()) { + for (auto directory : group) { + inputFiles.push_back(filesInFolder(directory)); + } + } + } + + m_bkgEvents = make_unique(inputFiles); + for (auto& valset : m_bkgEvents->m_totalNumberOfEvents) { + for (auto& val : valset) { + if (val == 0) { + std::string err = "No events found in the background files"; + for (auto& file : m_inputFileNames.value()) { + err += " " + file[0]; + } + error() << err << endmsg; + return StatusCode::FAILURE; + } + } + if (std::any_of(valset.begin(), valset.end(), + [this](const int& val) { return this->m_startWithBackgroundEvent >= val; })) { + throw GaudiException("StartWithBackgroundEvent is larger than the number of events in the background files", name(), + StatusCode::FAILURE); + } + } + + if (m_Noverlay.empty()) { + info() << "Using the default number of overlay events (1) for each group, since none was specified with " + "NumberBackground " + << endmsg; + m_Noverlay = std::vector(m_bkgEvents->size(), 1); + } + + if (m_Poisson.empty()) { + info() << "Using the default overlay mode (no Poission distribution) for each group, since none was specified with " + "Poisson_random_NOverlay" + << endmsg; + m_Poisson = std::vector(m_bkgEvents->size(), false); + } + + return StatusCode::SUCCESS; +} + +retType OverlayTimingRandomMix::operator()( + const edm4hep::EventHeaderCollection& headers, + const edm4hep::MCParticleCollection& particles, + const std::vector& simTrackerHits, + const std::vector& simCaloHits) const { + + const auto seed = m_uidSvc->getUniqueID(headers[0].getEventNumber(), headers[0].getRunNumber(), this->name()); + m_engine.seed(seed); + + // Output collections + auto oparticles = edm4hep::MCParticleCollection(); + auto osimTrackerHits = std::vector(); + auto osimCaloHits = std::vector(); + auto ocaloHitContribs = std::vector(); + for (size_t i = 0; i < simCaloHits.size(); ++i) { + ocaloHitContribs.emplace_back(edm4hep::CaloHitContributionCollection()); + } + + // Copy MCParticles for physics event into a new collection + for (const auto&& part : particles) { + oparticles->push_back(part.clone(false)); + } + // Fix relations to point to the new particles + for (size_t i = 0; i < particles.size(); ++i) { + for (const auto& parent : particles[i].getParents()) { + oparticles[i].addToParents(oparticles[parent.getObjectID().index]); + } + for (const auto& daughter : particles[i].getDaughters()) { + oparticles[i].addToDaughters(oparticles[daughter.getObjectID().index]); + } + } + + // Copy the SimTrackerHits and crop them + for (size_t i = 0; i < simTrackerHits.size(); ++i) { + const auto& coll = simTrackerHits[i]; + const auto name = inputLocations(SIMTRACKERHIT_INDEX_POSITION)[i]; + const auto [this_start, this_stop] = define_time_windows(name); + auto ocoll = edm4hep::SimTrackerHitCollection(); + for (const auto&& simTrackerHit : *coll) { + const float tof = time_of_flight(simTrackerHit.getPosition()); + if ((simTrackerHit.getTime() > this_start + tof) && (simTrackerHit.getTime() < this_stop + tof)) { + auto nhit = simTrackerHit.clone(false); + if (simTrackerHit.getParticle().getObjectID().index != -1) + nhit.setParticle(oparticles[simTrackerHit.getParticle().getObjectID().index]); + ocoll->push_back(nhit); + } + } + osimTrackerHits.emplace_back(std::move(ocoll)); + } + + // Copy the SimCalorimeterHits and crop them together with the contributions + std::map> cellIDsMap; + for (size_t i = 0; i < simCaloHits.size(); ++i) { + const auto& coll = simCaloHits[i]; + const auto name = inputLocations(SIMCALOHIT_INDEX_POSITION)[i]; + const auto [this_start, this_stop] = define_time_windows(name); + auto& calHitMap = cellIDsMap[i]; + auto& caloHitContribs = ocaloHitContribs[i]; + for (const auto&& simCaloHit : *coll) { + const float tof = time_of_flight(simCaloHit.getPosition()); + bool within_time_window = false; + std::vector thisContribs; + for (const auto& contrib : simCaloHit.getContributions()) { + if (!((contrib.getTime() > this_start + tof) && (contrib.getTime() < this_stop + tof))) + continue; + within_time_window = true; + // TODO: Make sure a contribution is not added twice + auto newContrib = contrib.clone(false); + newContrib.setParticle(oparticles[contrib.getParticle().getObjectID().index]); + thisContribs.push_back(caloHitContribs.size()); + caloHitContribs.push_back(std::move(newContrib)); + } + if (within_time_window) { + auto newhit = simCaloHit.clone(false); + for (const auto& contrib : thisContribs) { + newhit.addToContributions(caloHitContribs[contrib]); + } + calHitMap.emplace(simCaloHit.getCellID(), std::move(newhit)); + } + } + } + + // Iterate over each group of files and parameters + for (size_t groupIndex = 0; groupIndex < m_bkgEvents->size(); groupIndex++) { + if (m_randomBX) { + m_physBX = std::uniform_int_distribution(0, m_NBunchTrain - 1)(m_engine); + debug() << "Physics Event was placed in the " << m_physBX << " bunch crossing!" << endmsg; + } + + // define a permutation for the events to overlay -- the physics event is per definition at position 0 + std::vector permutation; + std::vector v_file_indices(m_bkgEvents->m_fileNames[groupIndex].size()); // vector of indices + + // Permutation has negative values and the last one is 0 + // if (!m_randomBX) then m_physBX (default = 1) + for (int i = -(m_physBX - 1); i < m_NBunchTrain - (m_physBX - 1); ++i) { + permutation.push_back(i); + } + std::iota(std::begin(v_file_indices), std::end(v_file_indices), 0); // Fill with 0, 1, ... + std::shuffle(permutation.begin(), permutation.end(), m_engine); + std::shuffle(v_file_indices.begin(), v_file_indices.end(), m_engine); + + // TODO: Check that there is anything to overlay + + debug() << "Starting overlay at event: " << m_bkgEvents->m_nextEntry[groupIndex] << " for the background group " + << groupIndex << endmsg; + + if (m_startWithBackgroundEvent >= 0) { + info() << "Skipping to event: " << m_startWithBackgroundEvent << endmsg; + for (auto& entry : m_bkgEvents->m_nextEntry[groupIndex]) { + entry = m_startWithBackgroundEvent; + } + } + + // Overlay the background events to each bunchcrossing in the bunch train + for (int bxInTrain = 0; bxInTrain < m_NBunchTrain; ++bxInTrain) { + const int BX_number_in_train = permutation.at(bxInTrain); + + int NOverlay_to_this_BX = 0; + + if (m_Poisson[groupIndex]) { + NOverlay_to_this_BX = std::poisson_distribution<>(m_Noverlay[groupIndex])(m_engine); + } else { + NOverlay_to_this_BX = m_Noverlay[groupIndex]; + } + + debug() << "Will overlay " << NOverlay_to_this_BX << " events to BX number " << BX_number_in_train + m_physBX + << endmsg; + + for (int k = 0; k < NOverlay_to_this_BX; ++k) { + info() << "Overlaying background event " << m_bkgEvents->m_nextEntry[groupIndex][k] << " from group " << groupIndex + << " to BX " << bxInTrain << endmsg; + if (m_bkgEvents->m_nextEntry[groupIndex][k] >= m_bkgEvents->m_totalNumberOfEvents[groupIndex][k] && + !m_allowReusingBackgroundFiles) { + throw GaudiException("No more events in background file", name(), StatusCode::FAILURE); + } + podio::Reader reader = m_bkgEvents->open(groupIndex, v_file_indices[k]); + debug() << "File: " << m_bkgEvents->m_fileNames[groupIndex][v_file_indices[k]] + <<"\nNumber of Events: "<< reader.getEvents() << endmsg; + const auto backgroundEvent = reader.readEvent(m_bkgEvents->m_nextEntry[groupIndex][k]); + m_bkgEvents->m_nextEntry[groupIndex][k]++; + m_bkgEvents->m_nextEntry[groupIndex][k] %= m_bkgEvents->m_totalNumberOfEvents[groupIndex][k]; + const auto availableCollections = backgroundEvent.getAvailableCollections(); + + // Either 0 or negative + const auto timeOffset = BX_number_in_train * m_deltaT; + + if (std::find(availableCollections.begin(), availableCollections.end(), m_MCParticleCollectionName) == + availableCollections.end()) { + warning() << "Collection " << m_MCParticleCollectionName << " not found in background event" << endmsg; + } + + // To fix the relations we will need to have a map from old to new particle index + std::map oldToNewMap; + std::map, std::vector>> parentDaughterMap; + + if (m_mergeMCParticles) { + const auto& bgParticles = backgroundEvent.get(m_MCParticleCollectionName); + int j = oparticles.size(); + for (size_t i = 0; i < bgParticles.size(); ++i) { + auto npart = bgParticles[i].clone(false); + + npart.setTime(bgParticles[i].getTime() + timeOffset); + npart.setOverlay(true); + oparticles->push_back(npart); + for (const auto& parent : bgParticles[i].getParents()) { + parentDaughterMap[j].first.push_back(parent.getObjectID().index); + } + for (const auto& daughter : bgParticles[i].getDaughters()) { + parentDaughterMap[j].second.push_back(daughter.getObjectID().index); + } + oldToNewMap[i] = j; + j++; + } + for (const auto& [index, parentsDaughters] : parentDaughterMap) { + const auto& [parents, daughters] = parentsDaughters; + for (const auto& parent : parents) { + if (parentDaughterMap.find(oldToNewMap[parent]) == parentDaughterMap.end()) { + // warning() << "Parent " << parent << " not found in background event" << endmsg; + continue; + } + oparticles[index].addToParents(oparticles[oldToNewMap[parent]]); + } + for (const auto& daughter : daughters) { + if (parentDaughterMap.find(oldToNewMap[daughter]) == parentDaughterMap.end()) { + // warning() << "Parent " << daughter << " not found in background event" << endmsg; + continue; + } + // info() << "Adding (daughter) " << daughter << " to " << index << endmsg; + oparticles[index].addToDaughters(oparticles[oldToNewMap[daughter]]); + } + } + } + + for (size_t i = 0; i < simTrackerHits.size(); ++i) { + const auto name = inputLocations(SIMTRACKERHIT_INDEX_POSITION)[i]; + debug() << "Processing collection " << name << endmsg; + if (std::find(availableCollections.begin(), availableCollections.end(), name) == availableCollections.end()) { + warning() << "Collection " << name << " not found in background event" << endmsg; + continue; + } + const auto [this_start, this_stop] = define_time_windows(name); + // There are only contributions to the readout if the hits are in the integration window + if (this_stop <= (BX_number_in_train - m_physBX) * m_deltaT) { + info() << "Skipping collection " << name << " as it is not in the integration window" << endmsg; + continue; + } + auto& ocoll = osimTrackerHits[i]; + for (const auto&& simTrackerHit : backgroundEvent.get(name)) { + const float tof = time_of_flight(simTrackerHit.getPosition()); + + if (!((simTrackerHit.getTime() + timeOffset > this_start + tof) && + (simTrackerHit.getTime() + timeOffset < this_stop + tof))) + continue; + auto nhit = simTrackerHit.clone(false); + nhit.setOverlay(true); + nhit.setTime(simTrackerHit.getTime() + timeOffset); + if (m_mergeMCParticles) { + nhit.setParticle(oparticles[oldToNewMap[simTrackerHit.getParticle().getObjectID().index]]); + } else { + edm4hep::MCParticle mcp = simTrackerHit.getParticle(); + if (mcp.isAvailable()) { + // Preserve Momentum + edm4hep::Vector3d mom = mcp.getMomentum(); + nhit.setMomentum({(float)mom.x, (float)mom.y, (float)mom.z}); + } + } + ocoll->push_back(nhit); + } + } + + for (size_t i = 0; i < simCaloHits.size(); ++i) { + const auto name = inputLocations(SIMCALOHIT_INDEX_POSITION)[i]; + debug() << "Processing collection " << name << endmsg; + if (std::find(availableCollections.begin(), availableCollections.end(), name) == availableCollections.end()) { + warning() << "Collection " << name << " not found in background event" << endmsg; + continue; + } + const auto [this_start, this_stop] = define_time_windows(name); + // There are only contributions to the readout if the hits are in the integration window + if (this_stop <= (BX_number_in_train - m_physBX) * m_deltaT) { + info() << "Skipping collection " << name << " as it is not in the integration window" << endmsg; + continue; + } + + auto& calHitMap = cellIDsMap[i]; + auto& calHitContribs = ocaloHitContribs[i]; + for (const auto&& simCaloHit : backgroundEvent.get(name)) { + if (calHitMap.find(simCaloHit.getCellID()) == calHitMap.end()) { + // There is no hit at this position. The new hit can be added, if it is not outside the window + auto calhit = edm4hep::MutableSimCalorimeterHit(); + bool add = false; + for (const auto& contrib : simCaloHit.getContributions()) { + if ((contrib.getTime() + timeOffset > this_start) && (contrib.getTime() + timeOffset < this_stop)) { + add = true; + // TODO: Make sure a contribution is not added twice + auto newContrib = contrib.clone(false); + if (m_mergeMCParticles) { + newContrib.setParticle(oparticles[oldToNewMap[contrib.getParticle().getObjectID().index]]); + } else { + newContrib.setParticle(edm4hep::MCParticle()); + } + newContrib.setTime(contrib.getTime() + timeOffset); + calhit.addToContributions(newContrib); + calHitContribs.push_back(newContrib); + } + } + if (add) { + calhit.setCellID(simCaloHit.getCellID()); + calhit.setEnergy(simCaloHit.getEnergy()); + calhit.setPosition(simCaloHit.getPosition()); + calHitMap[calhit.getCellID()] = calhit; + } + } else { + // there is already a hit at this position + auto& calhit = calHitMap[simCaloHit.getCellID()]; + for (const auto& contrib : simCaloHit.getContributions()) { + if ((contrib.getTime() + timeOffset > this_start) && (contrib.getTime() + timeOffset < this_stop)) { + // TODO: Make sure a contribution is not added twice + auto newContrib = contrib.clone(false); + if (m_mergeMCParticles) { + newContrib.setParticle(oparticles[oldToNewMap[contrib.getParticle().getObjectID().index]]); + } else { + newContrib.setParticle(edm4hep::MCParticle()); + } + newContrib.setTime(contrib.getTime() + timeOffset); + calhit.addToContributions(newContrib); + calHitContribs.push_back(newContrib); + } + } + } + } + } + } + } + } + // Move the SimCalorimeterHitCollections to the output vector + // So far they are stored in a map with the cellID as key + // but they don't belong to any collection yet + for (const auto& [index, calHitMap] : cellIDsMap) { + auto ocoll = edm4hep::SimCalorimeterHitCollection(); + for (const auto& [cellID, hit] : calHitMap) { + ocoll->push_back(std::move(hit)); + } + osimCaloHits.emplace_back(std::move(ocoll)); + } + + debug() << "\n\t\tCollection\t\t|\t\tPre BIB\t\t|\t\tPost BIB\t\t\n--------------------------------------------------------------------------\n"; + for (int trkCol = 0; trkCol < simTrackerHits.size(); trkCol++) { + debug() << "\tTrackerHits " << trkCol << "\t|\t\t" << simTrackerHits[trkCol]->size() << "\t\t|\t\t" << osimTrackerHits[trkCol]->size() << "\n"; + } + for (int calCol = 0; calCol < simCaloHits.size(); calCol++) { + debug() << "\tCaloHits " << calCol << "\t|\t\t" << simCaloHits[calCol]->size() << "\t\t|\t\t" << osimCaloHits[calCol]->size() << "\n"; + } + debug() << endmsg; + + return std::make_tuple(std::move(oparticles), std::move(osimTrackerHits), std::move(osimCaloHits), + std::move(ocaloHitContribs)); +} + +StatusCode OverlayTimingRandomMix::finalize() { + if (m_copyCellIDMetadata) { + for (const auto& [input, output] : + {std::make_pair(inputLocations("SimTrackerHits"), outputLocations("OutputSimTrackerHits")), + std::make_pair(inputLocations("SimCalorimeterHits"), outputLocations("OutputSimCalorimeterHits"))}) { + for (size_t i = 0; i < input.size(); ++i) { + const auto value = k4FWCore::getParameter( + podio::collMetadataParamName(input[i], edm4hep::labels::CellIDEncoding), this); + if (value.has_value()) { + k4FWCore::putParameter(podio::collMetadataParamName(output[i], edm4hep::labels::CellIDEncoding), + value.value(), this); + } else { + warning() << "No metadata found for " << input[i] << " when copying CellID metadata was requested" << endmsg; + } + } + } + } + + return Gaudi::Algorithm::finalize(); +} + +DECLARE_COMPONENT(OverlayTimingRandomMix) diff --git a/k4Reco/Overlay/components/OverlayTimingRandomMix.h b/k4Reco/Overlay/components/OverlayTimingRandomMix.h new file mode 100644 index 0000000..afab73b --- /dev/null +++ b/k4Reco/Overlay/components/OverlayTimingRandomMix.h @@ -0,0 +1,152 @@ +/* + * 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. + */ + +/** Background overlay algorithm + + This algorithm overlays background events on top of the signal events. The + background events are read from sets of input files, and the signal events + are the input in the main event loop. + + The MCParticleCollection in signal are background are overlaid into one + collection. The SimTrackerHit collections are cropped and overlayed if they + are in the time window. The SimCalorimeterHit collections are overlayed + based on the cellID. If a signal hit has the same cellID as a background + hit, they are combined into a single hit. Only hits that have + CaloHitContributions in the time range are considered. + +**/ + +#include "podio/Frame.h" +#include "podio/Reader.h" + +#include "edm4hep/CaloHitContributionCollection.h" +#include "edm4hep/EventHeaderCollection.h" +#include "edm4hep/MCParticleCollection.h" +#include "edm4hep/SimCalorimeterHitCollection.h" +#include "edm4hep/SimTrackerHitCollection.h" + +#include "k4FWCore/Transformer.h" +#include "k4Interface/IUniqueIDGenSvc.h" + +// Needed for some of the more complex properties +#include "Gaudi/Parsers/Factory.h" +#include "Gaudi/Property.h" + +#include +#include +#include +#include + +struct EventHolder { + std::vector> m_fileNames; + std::vector> m_totalNumberOfEvents; + + std::vector> m_nextEntry; + + EventHolder(const std::vector>& fileNames) : m_fileNames(fileNames) { + m_totalNumberOfEvents.resize(m_fileNames.size()); + m_nextEntry.resize(m_fileNames.size()); + + for (int group = 0; group < m_fileNames.size(); group++) { + m_nextEntry[group].resize(m_fileNames[group].size()); + for (auto& name : m_fileNames[group]) { + m_totalNumberOfEvents[group].push_back(1);//m_rootFileReaders[group].back().getEntries("events")); + } + } + } + EventHolder() = default; + + podio::Reader open(int groupIndex, int index) { + return podio::makeReader(m_fileNames[groupIndex][index]); + } + + // TODO: Cache functionality + // podio::Frame& read + + size_t size() const { return m_fileNames.size(); } +}; + +using retType = + std::tuple, + std::vector, std::vector>; + +struct OverlayTimingRandomMix : public k4FWCore::MultiTransformer&, + const std::vector&)> { + OverlayTimingRandomMix(const std::string& name, ISvcLocator* svcLoc) + : MultiTransformer( + name, svcLoc, + {KeyValues("EventHeader", {"EventHeader"}), KeyValues("MCParticles", {"DefaultMCParticles"}), + KeyValues("SimTrackerHits", {"DefaultSimTrackerHits"}), + KeyValues("SimCalorimeterHits", {"DefaultSimCalorimeterHits"})}, + {KeyValues("OutputMCParticles", {"NewMCParticles"}), KeyValues("OutputSimTrackerHits", {"NewSimTrackerHits"}), + KeyValues("OutputSimCalorimeterHits", {"NewSimCalorimeterHits"}), + KeyValues("OutputCaloHitContributions", {"OverlayCaloHitContributions"})}) {} + + template void overlayCollection(std::string collName, const podio::CollectionBase& inColl); + + virtual StatusCode initialize() final; + virtual StatusCode finalize() final; + + retType virtual operator()( + const edm4hep::EventHeaderCollection& headers, const edm4hep::MCParticleCollection& mcParticles, + const std::vector& simTrackerHits, + const std::vector& simCalorimeterHits) const final; + + std::pair define_time_windows(const std::string& Collection_name) const; + +private: + // These correspond to the index position in the argument list + constexpr static int SIMTRACKERHIT_INDEX_POSITION = 2; + constexpr static int SIMCALOHIT_INDEX_POSITION = 3; + + Gaudi::Property m_randomBX{this, "RandomBx", false, "Place the physics event at an random position in the train: overrides PhysicsBX"}; + mutable Gaudi::Property m_physBX{this, "PhysicsBX", 1, "Number of the Bunch crossing of the physics event"}; + Gaudi::Property m_NBunchTrain{this, "NBunchtrain", 1, "Number of bunches in a bunch train"}; + // Gaudi::Property m_startWithBackgroundFile{this, "StartBackgroundFileIndex", -1, + // "Which background file to startWith"}; + Gaudi::Property m_startWithBackgroundEvent{this, "StartBackgroundEventIndex", -1, "Which background event to startWith"}; + + Gaudi::Property>> m_inputFileNames{this, "BackgroundFileNames", {}, "Name of the edm4hep input file(s) with background."}; + //Gaudi::Property m_separateEventFiles{this, "SeparateEventFiles", false, "For setups where each event is in a different file (like Muon Collider)."} + + Gaudi::Property> m_Noverlay{this, "NumberBackground", {}, "Number of Background events to overlay - either fixed or Poisson mean"}; + + Gaudi::Property> m_Poisson{this, "Poisson_random_NOverlay", {}, "Draw random number of Events to overlay from Poisson distribution with mean value NumberBackground"}; + + Gaudi::Property m_MCParticleCollectionName{this, "BackgroundMCParticleCollectionName", "MCParticle", "The name of the MCParticle collection in the background files"}; + + Gaudi::Property m_deltaT{this, "Delta_t", float(0.5), "Time difference between BXs in the BXtrain"}; + + mutable std::unique_ptr m_bkgEvents{}; + + Gaudi::Property>> m_timeWindows{this, "TimeWindows", std::map>(), "Time windows for the different collections"}; + Gaudi::Property m_allowReusingBackgroundFiles{this, "AllowReusingBackgroundFiles", false, "If true the same background file can be used for the same event"}; + Gaudi::Property m_copyCellIDMetadata{this, "CopyCellIDMetadata", false, "If metadata is found in the signal file, copy it to the output file, replacing the old names with the new names"}; + + Gaudi::Property m_mergeMCParticles{this, "MergeMCParticles", true, "Merge the MC Particle collections"}; + + // Gaudi::Property m_maxCachedFrames{ + // this, "MaxCachedFrames", 0, "Maximum number of frames cached from background files"}; + +private: + inline static thread_local std::mt19937 m_engine; + SmartIF m_uidSvc; +}; From a6ae4e71f41852bad5736bc905369115207955a4 Mon Sep 17 00:00:00 2001 From: samf25 Date: Thu, 31 Jul 2025 13:32:57 -0400 Subject: [PATCH 03/14] Move Filter --- k4Reco/CMakeLists.txt | 1 + k4Reco/{CaloDigi => Tracking}/include/FilterDoubleLayerHits.h | 0 k4Reco/{CaloDigi => Tracking}/src/FilterDoubleLayerHits.cc | 0 3 files changed, 1 insertion(+) rename k4Reco/{CaloDigi => Tracking}/include/FilterDoubleLayerHits.h (100%) rename k4Reco/{CaloDigi => Tracking}/src/FilterDoubleLayerHits.cc (100%) diff --git a/k4Reco/CMakeLists.txt b/k4Reco/CMakeLists.txt index b56433a..7e679e9 100644 --- a/k4Reco/CMakeLists.txt +++ b/k4Reco/CMakeLists.txt @@ -102,6 +102,7 @@ if(BUILD_TRACKING) Tracking/src/ClonesAndSplitTracksFinder.cpp Tracking/src/RefitFinal.cpp Tracking/src/TruthTrackFinder.cpp + Tracking/src/FilterDoubleLayerHits.cc ) list(APPEND libs diff --git a/k4Reco/CaloDigi/include/FilterDoubleLayerHits.h b/k4Reco/Tracking/include/FilterDoubleLayerHits.h similarity index 100% rename from k4Reco/CaloDigi/include/FilterDoubleLayerHits.h rename to k4Reco/Tracking/include/FilterDoubleLayerHits.h diff --git a/k4Reco/CaloDigi/src/FilterDoubleLayerHits.cc b/k4Reco/Tracking/src/FilterDoubleLayerHits.cc similarity index 100% rename from k4Reco/CaloDigi/src/FilterDoubleLayerHits.cc rename to k4Reco/Tracking/src/FilterDoubleLayerHits.cc From 0a8788875079f5fea70907d6e62dc0f1a7134465 Mon Sep 17 00:00:00 2001 From: Federico Meloni Date: Wed, 6 Aug 2025 15:55:33 +0200 Subject: [PATCH 04/14] depend on k4geo --- CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index 0911731..7127736 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -36,6 +36,7 @@ find_package(Gaudi REQUIRED) find_package(k4FWCore 1.4 REQUIRED) find_package(ROOT REQUIRED COMPONENTS RIO Tree) find_package(k4SimGeant4 REQUIRED) +find_package(k4Geo REQUIRED) include(cmake/Key4hepConfig.cmake) From e47537e3efc53e24fa78d132ebb47688d6c23f40 Mon Sep 17 00:00:00 2001 From: Federico Meloni Date: Sun, 10 Aug 2025 22:56:01 +0200 Subject: [PATCH 05/14] fix cmake find --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 7127736..af6d33f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -36,7 +36,7 @@ find_package(Gaudi REQUIRED) find_package(k4FWCore 1.4 REQUIRED) find_package(ROOT REQUIRED COMPONENTS RIO Tree) find_package(k4SimGeant4 REQUIRED) -find_package(k4Geo REQUIRED) +find_package(k4geo REQUIRED) include(cmake/Key4hepConfig.cmake) From 8e691a7579b6a58150e449447cd0371ec837a37a Mon Sep 17 00:00:00 2001 From: samf25 Date: Mon, 11 Aug 2025 13:02:24 -0400 Subject: [PATCH 06/14] CMake changes to allow for compilation in spack image --- CMakeLists.txt | 1 + k4Reco/CMakeLists.txt | 4 ++++ 2 files changed, 5 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index af6d33f..63faae9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -37,6 +37,7 @@ find_package(k4FWCore 1.4 REQUIRED) find_package(ROOT REQUIRED COMPONENTS RIO Tree) find_package(k4SimGeant4 REQUIRED) find_package(k4geo REQUIRED) +find_package(CLHEP REQUIRED) include(cmake/Key4hepConfig.cmake) diff --git a/k4Reco/CMakeLists.txt b/k4Reco/CMakeLists.txt index 7e679e9..a34cb7e 100644 --- a/k4Reco/CMakeLists.txt +++ b/k4Reco/CMakeLists.txt @@ -37,6 +37,8 @@ set(sources GaudiLumiCalClusterer/src/LCCluster.cpp GaudiLumiCalClusterer/src/SuperTrueClusterWeights.cpp GaudiLumiCalClusterer/src/ProjectionInfo.cpp + CaloDigi/src/RealisticCaloDigi.cc + CaloDigi/src/RealisticCaloReco.cc ) set(libs @@ -50,6 +52,7 @@ set(libs ROOT::Core ROOT::MathCore ROOT::Tree + CLHEP::Random ) if(BUILD_TRACKING) @@ -128,6 +131,7 @@ target_include_directories(k4RecoPlugins PUBLIC $ $ $ + $ $) if(BUILD_TRACKING) From f4d65fa575a1b0c4e172a51f53f7c175812e94dc Mon Sep 17 00:00:00 2001 From: Thomas Madlener Date: Tue, 12 Aug 2025 14:15:06 +0200 Subject: [PATCH 07/14] Use the correct header files The others are no longer present since the switch to templated links --- k4Reco/CaloDigi/src/RealisticCaloDigi.cc | 2 +- k4Reco/CaloDigi/src/RealisticCaloReco.cc | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/k4Reco/CaloDigi/src/RealisticCaloDigi.cc b/k4Reco/CaloDigi/src/RealisticCaloDigi.cc index 4ec1458..5d4dcc4 100644 --- a/k4Reco/CaloDigi/src/RealisticCaloDigi.cc +++ b/k4Reco/CaloDigi/src/RealisticCaloDigi.cc @@ -2,7 +2,7 @@ #include "RealisticCaloDigi.h" #include -#include +#include #include #include diff --git a/k4Reco/CaloDigi/src/RealisticCaloReco.cc b/k4Reco/CaloDigi/src/RealisticCaloReco.cc index 213a75b..f6f88d1 100644 --- a/k4Reco/CaloDigi/src/RealisticCaloReco.cc +++ b/k4Reco/CaloDigi/src/RealisticCaloReco.cc @@ -4,7 +4,7 @@ #include #include -#include +#include #include #include From 3890430e63d52395252ea08ffe9aaa12c3e64733 Mon Sep 17 00:00:00 2001 From: Samuel Ferraro <113254671+samf25@users.noreply.github.com> Date: Mon, 18 Aug 2025 14:33:01 -0400 Subject: [PATCH 08/14] Encoding String was wrong --- k4Reco/CaloDigi/include/RealisticCaloDigi.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/k4Reco/CaloDigi/include/RealisticCaloDigi.h b/k4Reco/CaloDigi/include/RealisticCaloDigi.h index be2bc16..ab4fc72 100644 --- a/k4Reco/CaloDigi/include/RealisticCaloDigi.h +++ b/k4Reco/CaloDigi/include/RealisticCaloDigi.h @@ -92,7 +92,7 @@ struct RealisticCaloDigi : k4FWCore::MultiTransformer< Gaudi::Property m_elec_noiseMip{this, "elec_noise_mip", 0.0f, "Typical electronics noise (in MIP units)"}; Gaudi::Property m_elec_rangeMip{this, "elec_range_mip", 2500.0f, "Maximum of dynamic range of electronics (in MIPs)"}; // code for layer info for cellID decoder - Gaudi::Property m_encodingStringVariable{this, "EncodingStringParameterName", "GlobalTrackerReadoutID", "The name of the DD4hep constant that contains the Encoding string for tracking detectors"}; + Gaudi::Property m_encodingStringVariable{this, "EncodingStringParameterName", "GlobalCalorimeterReadoutID", "The name of the DD4hep constant that contains the Encoding string for tracking detectors"}; // energy threshold Gaudi::Property m_threshold_value{this, "threshold", 0.5f, "Threshold for Hit"}; Gaudi::Property m_threshold_unit{this, "thresholdUnit", std::string("MIP"), "Unit for threshold. Can be \"GeV\", \"MIP\" or \"px\". MIP and px need properly set calibration constants"}; From 36654bf7271b026959179c186ea74434dd013b4d Mon Sep 17 00:00:00 2001 From: Samuel Ferraro <113254671+samf25@users.noreply.github.com> Date: Mon, 18 Aug 2025 14:33:40 -0400 Subject: [PATCH 09/14] Encoding String was wrong here, too --- k4Reco/CaloDigi/include/RealisticCaloReco.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/k4Reco/CaloDigi/include/RealisticCaloReco.h b/k4Reco/CaloDigi/include/RealisticCaloReco.h index 99d83f1..f8fb77d 100644 --- a/k4Reco/CaloDigi/include/RealisticCaloReco.h +++ b/k4Reco/CaloDigi/include/RealisticCaloReco.h @@ -60,7 +60,7 @@ struct RealisticCaloReco : k4FWCore::MultiTransformer> m_calibrCoeff{this, "calibration_factorsMipGev", {}, "Calibration coefficients (MIP->shower GeV) of layers groups"}; // Cell ID layer string - Gaudi::Property m_encodingStringVariable{this, "EncodingStringParameterName", "GlobalTrackerReadoutID", "The name of the DD4hep constant that contains the Encoding string for tracking detectors"}; + Gaudi::Property m_encodingStringVariable{this, "EncodingStringParameterName", "GlobalCalorimeterReadoutID", "The name of the DD4hep constant that contains the Encoding string for tracking detectors"}; SmartIF m_geoSvc; From 92246c955ed2d9f4b06266153931f6f779a364f1 Mon Sep 17 00:00:00 2001 From: Samuel Ferraro Date: Tue, 25 Nov 2025 10:13:35 -0800 Subject: [PATCH 10/14] Update Refitter --- k4Reco/GaudiTrkUtils/include/GaudiDDKalTest.h | 2 +- k4Reco/GaudiTrkUtils/src/GaudiDDKalTest.cpp | 7 +- k4Reco/Tracking/include/RefitFinal.h | 89 ++++------ k4Reco/Tracking/src/RefitFinal.cpp | 160 ++++++++---------- 4 files changed, 104 insertions(+), 154 deletions(-) diff --git a/k4Reco/GaudiTrkUtils/include/GaudiDDKalTest.h b/k4Reco/GaudiTrkUtils/include/GaudiDDKalTest.h index b502d03..8137711 100644 --- a/k4Reco/GaudiTrkUtils/include/GaudiDDKalTest.h +++ b/k4Reco/GaudiTrkUtils/include/GaudiDDKalTest.h @@ -77,7 +77,7 @@ class GaudiDDKalTest { // virtual void setOption(unsigned CFGOption, bool val); /** initialise track fitter system */ - void init(); + void init(bool msOn = true, bool energyLossOn = true); // Copy the encoder void setEncoder(const dd4hep::DDSegmentation::BitFieldCoder& encoder) { m_encoder = encoder; } diff --git a/k4Reco/GaudiTrkUtils/src/GaudiDDKalTest.cpp b/k4Reco/GaudiTrkUtils/src/GaudiDDKalTest.cpp index 8dd91c4..9bc0f4a 100644 --- a/k4Reco/GaudiTrkUtils/src/GaudiDDKalTest.cpp +++ b/k4Reco/GaudiTrkUtils/src/GaudiDDKalTest.cpp @@ -69,10 +69,9 @@ GaudiDDKalTest::~GaudiDDKalTest() { // } // } -void GaudiDDKalTest::init() { - // TODO: Don't hardcode the options - this->includeMultipleScattering(true); - this->includeEnergyLoss(true); +void GaudiDDKalTest::init(bool msOn, bool energyLossOn) { + this->includeMultipleScattering(msOn); + this->includeEnergyLoss(energyLossOn); // this->sm // std::cout << " -------------------------------------------------------------------------------- " << endmsg; diff --git a/k4Reco/Tracking/include/RefitFinal.h b/k4Reco/Tracking/include/RefitFinal.h index d695b27..68eccd0 100644 --- a/k4Reco/Tracking/include/RefitFinal.h +++ b/k4Reco/Tracking/include/RefitFinal.h @@ -35,80 +35,53 @@ #include #include -struct RefitFinal final - : k4FWCore::MultiTransformer( - const edm4hep::TrackCollection&, const std::vector&)> { +struct RefitFinal final : k4FWCore::MultiTransformer< + std::tuple( + const edm4hep::TrackCollection&, const std::vector&)> { RefitFinal(const std::string& name, ISvcLocator* svcLoc); StatusCode initialize() override; - std::tuple - operator()(const edm4hep::TrackCollection&, - const std::vector&) const override; - - int FitInit2(const edm4hep::Track& track, GaudiDDKalTestTrack& _marlinTrk) const; - - // /* helper function to get collection using try catch block */ - // lcio::LCCollection* GetCollection(lcio::LCEvent* evt, std::string colName); - - // /** Input track collection name for refitting. - // */ - // std::string _input_track_col_name = "TruthTracks"; - - // /** output track collection name. - // */ - // std::string _output_track_col_name = "RefittedTracks"; - - // /** Input track relations name. - // */ - // std::string _input_track_rel_name = "SiTrackRelations"; - - // /** Output track relations name for refitting. - // */ - // std::string _output_track_rel_name = "RefittedRelation"; - - // /** pointer to the IMarlinTrkSystem instance - // */ - // MarlinTrk::IMarlinTrkSystem* _trksystem = nullptr; - - // int _n_run = -1; - - // bool _MSOn = true; - // bool _ElossOn = true; - // bool _SmoothOn = false; - // double _Max_Chi2_Incr = DBL_MAX; - // int _refPoint = -1; - - // bool _extrapolateForward = true; - // int _minClustersOnTrackAfterFit = 0; - - // std::shared_ptr _encoder{}; - - // registerProcessorParameter("Max_Chi2_Incr", "maximum allowable chi2 increment when moving from one site to - // another", - // _Max_Chi2_Incr, _Max_Chi2_Incr); + std::tuple operator()( + const edm4hep::TrackCollection&, + const std::vector& + ) const override; Gaudi::Property m_MSOn{this, "MultipleScatteringOn", true, "Use MultipleScattering in Fit"}; Gaudi::Property m_ElossOn{this, "EnergyLossOn", true, "Use Energy Loss in Fit"}; Gaudi::Property m_SmoothOn{this, "SmoothOn", false, "Smooth All Mesurement Sites in Fit"}; + + Gaudi::Property m_initialTrackError_d0{this, "InitialTrackErrorD0", 1.e6f, + "Value used for the initial d0 variance of the trackfit"}; + Gaudi::Property m_initialTrackError_phi0{this, "InitialTrackErrorPhi0", 1.e2f, + "Value used for the initial phi0 variance of the trackfit"}; + Gaudi::Property m_initialTrackError_omega{this, "InitialTrackErrorOmega", 1.e-4f, + "Value used for the initial omega variance of the trackfit"}; + Gaudi::Property m_initialTrackError_z0{this, "InitialTrackErrorZ0", 1.e6f, + "Value used for the initial z0 variance of the trackfit"}; + Gaudi::Property m_initialTrackError_tanL{this, "InitialTrackErrorTanL", 1.e2f, + "Value used for the initial tanL variance of the trackfit"}; Gaudi::Property m_Max_Chi2_Incr{this, "Max_Chi2_Incr", std::numeric_limits::max(), - "maximum allowable chi2 increment when moving from one site to another"}; - Gaudi::Property m_refPoint{ - this, "ReferencePoint", -1, - "Identifier of the reference point to use for the fit initialisation, -1 means at 0 0 0"}; - Gaudi::Property m_extrapolateForward{ - this, "extrapolateForward", true, - "if true extrapolation in the forward direction (in-out), otherwise backward (out-in)"}; + "maximum allowable chi2 increment when moving from one site to another"}; + Gaudi::Property m_refPoint{this, "ReferencePoint", -1, + "Identifier of the reference point to use for the fit initialisation, -1 means at 0 0 0"}; + Gaudi::Property m_extrapolateForward{this, "extrapolateForward", true, + "if true extrapolation in the forward direction (in-out), otherwise backward (out-in)"}; Gaudi::Property m_minClustersOnTrackAfterFit{this, "MinClustersOnTrackAfterFit", 4, - "Final minimum number of track clusters"}; + "Final minimum number of track clusters"}; + Gaudi::Property m_maxOutliersAllowed{this, "MaxOutliersAllowed", 99, + "Maximum number of outliers allowed on the refitted track"}; + Gaudi::Property m_ReducedChi2Cut{this, "ReducedChi2Cut", -1.0, + "Cut on maximum allowed reduced chi2"}; Gaudi::Property m_geoSvcName{this, "GeoSvcName", "GeoSvc", "The name of the GeoSvc instance"}; - Gaudi::Property m_encodingStringVariable{ - this, "EncodingStringParameterName", "GlobalTrackerReadoutID", + Gaudi::Property m_encodingStringVariable{this, "EncodingStringParameterName", "GlobalTrackerReadoutID", "The name of the DD4hep constant that contains the Encoding string for tracking detectors"}; + GaudiDDKalTest m_ddkaltest{this}; SmartIF m_geoSvc; + float m_bField = 5.0; dd4hep::DDSegmentation::BitFieldCoder m_encoder; }; diff --git a/k4Reco/Tracking/src/RefitFinal.cpp b/k4Reco/Tracking/src/RefitFinal.cpp index a63c0ab..2d0d6cf 100644 --- a/k4Reco/Tracking/src/RefitFinal.cpp +++ b/k4Reco/Tracking/src/RefitFinal.cpp @@ -53,11 +53,6 @@ StatusCode RefitFinal::initialize() { streamlog::logscope* scope = new streamlog::logscope(streamlog::out); setStreamlogOutputLevel(this, scope); - // usually a good idea to - // printParameters(); - - // _trksystem = MarlinTrk::Factory::createMarlinTrkSystem("DDKalTest", nullptr, ""); - m_geoSvc = serviceLocator()->service(m_geoSvcName); if (!m_geoSvc) { error() << "Unable to retrieve GeoSvc" << endmsg; @@ -66,35 +61,16 @@ StatusCode RefitFinal::initialize() { std::string cellIDEncodingString = m_geoSvc->constantAsString(m_encodingStringVariable.value()); m_encoder = dd4hep::DDSegmentation::BitFieldCoder(cellIDEncodingString); - m_ddkaltest.init(); + m_ddkaltest.init(m_MSOn, m_ElossOn); m_ddkaltest.setEncoder(m_encoder); - // /////////////////////////////// - - // _encoder = std::make_shared(lcio::LCTrackerCellID::encoding_string()); - - // if (not _trksystem) { - // throw EVENT::Exception("Cannot initialize MarlinTrkSystem of Type: DDKalTest"); - // } - - // _trksystem->setOption(MarlinTrk::IMarlinTrkSystem::CFG::useQMS, _MSOn); - // _trksystem->setOption(MarlinTrk::IMarlinTrkSystem::CFG::usedEdx, _ElossOn); - // _trksystem->setOption(MarlinTrk::IMarlinTrkSystem::CFG::useSmoothing, _SmoothOn); - // _trksystem->init(); - - // _n_run = 0; - return StatusCode::SUCCESS; } -std::tuple -RefitFinal::operator()(const edm4hep::TrackCollection& input_track_col, - const std::vector& input_rel_col) const { - // // set the correct configuration for the tracking system for this event - // MarlinTrk::TrkSysConfig mson(_trksystem, _MSOn); - // MarlinTrk::TrkSysConfig elosson(_trksystem, _ElossOn); - // MarlinTrk::TrkSysConfig smoothon(_trksystem, _SmoothOn); - +std::tuple RefitFinal::operator()( + const edm4hep::TrackCollection& input_track_col, + const std::vector& input_rel_col +) const { // establish the track collection that will be created edm4hep::TrackCollection trackVec; edm4hep::TrackMCParticleLinkCollection trackRelationCollection; @@ -103,11 +79,22 @@ RefitFinal::operator()(const edm4hep::TrackCollection& input_track_col, debug() << "No input relation collection, not creating one either" << endmsg; } + std::map trackIndexToMCParticle; + if (!input_rel_col.empty()) { + const auto& relationCol = *input_rel_col[0]; + for (const auto& rel : relationCol) { + edm4hep::Track trk = rel.getFrom(); + trackIndexToMCParticle[trk.getObjectID().index] = rel.getTo(); + } + } + const size_t nTracks = input_track_col.size(); debug() << " Number of Tracks " << nTracks << endmsg; // loop over the input tracks and refit + int counter = 0; + std::map hitInSubDet; for (size_t iTrack = 0; iTrack < nTracks; ++iTrack) { const auto& track = input_track_col.at(iTrack); const auto& trkHits = track.getTrackerHits(); @@ -118,51 +105,56 @@ RefitFinal::operator()(const edm4hep::TrackCollection& input_track_col, trkHitsPtr.push_back(&hit); } - auto marlin_trk = GaudiDDKalTestTrack(this, const_cast(&m_ddkaltest)); + auto gaudi_trk = GaudiDDKalTestTrack(this, const_cast(&m_ddkaltest)); debug() << "---- track n = " << iTrack << " n hits = " << trkHits.size() << endmsg; - - for (const auto& ptr : trkHitsPtr) { - marlin_trk.addHit(ptr); - } - - int init_status = FitInit2(track, marlin_trk); - - if (init_status != 0) { + if (trkHits.size() < 3) { + debug() << "Track " << iTrack << " has less than 3 hits, skipping" << endmsg; continue; } - // debug() << "Refit: Trackstate after initialisation\n" << marlin_trk.toString() << endmsg; + hitInSubDet.clear(); - debug() << "track initialised " << endmsg; + for (const auto& ptr : trkHitsPtr) { + gaudi_trk.addHit(ptr); + ++hitInSubDet[m_encoder.get(ptr->getCellID(), "system")]; + } - int fit_status = marlin_trk.fit(); - // debug() << "RefitHit: Trackstate after fit()\n" << marlin_trk.toString() << endmsg; + edm4hep::MutableTrack edm4hep_trk = trackVec.create(); - if (fit_status != 0) { - continue; - } + edm4hep::CovMatrix6f initialCov; + initialCov[0] = m_initialTrackError_d0; + initialCov[2] = m_initialTrackError_phi0; + initialCov[5] = m_initialTrackError_omega; + initialCov[9] = m_initialTrackError_z0; + initialCov[14] = m_initialTrackError_tanL; - edm4hep::MutableTrack lcio_trk; + const bool fit_direction = m_extrapolateForward; - GaudiTrkUtils trkUtils(static_cast(this), m_ddkaltest, m_geoSvc, - m_encodingStringVariable.value()); + GaudiTrkUtils trkUtils( + static_cast(this), + m_ddkaltest, m_geoSvc, + m_encodingStringVariable.value() + ); - int return_code = trkUtils.finaliseLCIOTrack(marlin_trk, lcio_trk, trkHitsPtr, true); + int return_code = trkUtils.createFinalisedLCIOTrack( + gaudi_trk, trkHitsPtr, edm4hep_trk, + fit_direction, initialCov, m_bField, + m_Max_Chi2_Incr.value() + ); if (return_code != 0) { debug() << "finaliseLCIOTrack failed" << endmsg; continue; } - // debug() << " *** created finalized LCIO track - return code " << return_code << std::endl << *lcio_trk << endmsg; + debug() << " *** created finalized LCIO track - return code " << return_code << endmsg; // fit finished - get hits in the fit - // remember the hits are ordered in the order in which they were fitted - const auto hits_in_fit = marlin_trk.getHitsInFit(); + const auto hits_in_fit = gaudi_trk.getHitsInFit(); if (int(hits_in_fit.size()) < m_minClustersOnTrackAfterFit) { debug() << "Less than " << m_minClustersOnTrackAfterFit @@ -172,7 +164,7 @@ RefitFinal::operator()(const edm4hep::TrackCollection& input_track_col, continue; } - const auto outliers = marlin_trk.getOutliers(); + const auto outliers = gaudi_trk.getOutliers(); std::vector all_hits; std::vector hits_in_fit_ptr; @@ -194,52 +186,38 @@ RefitFinal::operator()(const edm4hep::TrackCollection& input_track_col, trkUtils.addHitNumbersToTrack(subdetectorHitNumbers, all_hits, false, encoder2); trkUtils.addHitNumbersToTrack(subdetectorHitNumbers, hits_in_fit_ptr, true, encoder2); for (const auto num : subdetectorHitNumbers) { - lcio_trk.addToSubdetectorHitNumbers(num); + edm4hep_trk.addToSubdetectorHitNumbers(num); } - // debug() << "processEvent: Hit numbers for track " << lcio_trk.id() << ": " << endmsg; + // debug() << "processEvent: Hit numbers for track " << edm4hep_trk.id() << ": " << endmsg; int detID = 0; - for (size_t ip = 0; ip < lcio_trk.getSubdetectorHitNumbers().size(); ip = ip + 2) { + for (size_t ip = 0; ip < edm4hep_trk.getSubdetectorHitNumbers().size(); ip = ip + 2) { detID++; - debug() << " det id " << detID << " , nhits in track = " << lcio_trk.getSubdetectorHitNumbers()[ip] - << " , nhits in fit = " << lcio_trk.getSubdetectorHitNumbers()[ip + 1] << endmsg; - if (lcio_trk.getSubdetectorHitNumbers()[ip] > 0) + debug() << " det id " << detID << " , nhits in track = " << edm4hep_trk.getSubdetectorHitNumbers()[ip] + << " , nhits in fit = " << edm4hep_trk.getSubdetectorHitNumbers()[ip + 1] << endmsg; + if (edm4hep_trk.getSubdetectorHitNumbers()[ip] > 0) // TODO: is detID - 1 correct? - lcio_trk.setType(lcio_trk.getType() | (1 << detID)); + edm4hep_trk.setType(edm4hep_trk.getType() | (1 << detID)); } - // TODO: - // if (input_rel_col) { - // auto mcParticleVec = relation->getRelatedToObjects(track); - // auto weightVec = relation->getRelatedToWeights(track); - // for (size_t i = 0; i < mcParticleVec.size(); ++i) { - // LCRelationImpl* relationTrack = new LCRelationImpl(lcioTrkPtr, mcParticleVec[i], weightVec[i]); - // trackRelationCollection->addElement(relationTrack); - // } - // } - - trackVec.push_back(lcio_trk); - } // for loop to the tracks - - // TODO: - // if (input_rel_col) { - // evt->addCollection(trackRelationCollection, _output_track_rel_name); - // } - return std::make_tuple(std::move(trackVec), std::move(trackRelationCollection)); -} - -int RefitFinal::FitInit2(const edm4hep::Track& track, GaudiDDKalTestTrack& marlinTrk) const { - edm4hep::TrackState trackState; + // if required apply the ReducedChi2 cut + if (m_ReducedChi2Cut > 0.0 && edm4hep_trk.getChi2() / edm4hep_trk.getNdf() > m_ReducedChi2Cut) { + debug() << "Track Discarded due to ReducedChi2 cut: Chi2 = " << edm4hep_trk.getChi2() + << " Ndf = " << edm4hep_trk.getNdf() << " ReducedChi2 = " << edm4hep_trk.getChi2() / edm4hep_trk.getNdf() + << endmsg; + ++counter; + continue; + } - size_t refPoint = m_refPoint.value() == -1 ? 0 : static_cast(m_refPoint.value()); - if (refPoint >= track.getTrackStates().size()) { - error() << "Cannot find trackstate for " << m_refPoint << endmsg; - return 1; - } - trackState = track.getTrackStates()[refPoint]; + //trackVec.push_back(lcio_trk); + // create the Track-MCParticle relation + edm4hep::MutableTrackMCParticleLink relation = trackRelationCollection.create(); + relation.setFrom(edm4hep_trk); + relation.setTo(trackIndexToMCParticle.at(static_cast(iTrack))); + } // for loop to the tracks - const bool direction = m_extrapolateForward ? true : false; - marlinTrk.initialise(trackState, direction); + debug() << "Final number of Tracks after refit = " << trackVec.size() + << " Number of discarded Tracks = " << counter << endmsg; - return 0; + return std::make_tuple(std::move(trackVec), std::move(trackRelationCollection)); } From d0916c7611795ac5a46bc58c9fef5fc1892eb462 Mon Sep 17 00:00:00 2001 From: eat4978 Date: Wed, 10 Dec 2025 16:27:07 -0800 Subject: [PATCH 11/14] Move EventHolder definition inside of namespace --- k4Reco/Overlay/components/OverlayTiming.cpp | 2 +- k4Reco/Overlay/components/OverlayTiming.h | 41 ++++++++-------- .../components/OverlayTimingRandomMix.cpp | 2 +- .../components/OverlayTimingRandomMix.h | 49 ++++++++++--------- 4 files changed, 48 insertions(+), 46 deletions(-) diff --git a/k4Reco/Overlay/components/OverlayTiming.cpp b/k4Reco/Overlay/components/OverlayTiming.cpp index a38486e..5bc76aa 100644 --- a/k4Reco/Overlay/components/OverlayTiming.cpp +++ b/k4Reco/Overlay/components/OverlayTiming.cpp @@ -69,7 +69,7 @@ StatusCode OverlayTiming::initialize() { // TODO:: shuffle input files // std::shuffle(inputFiles.begin(), inputFiles.end(), rng_engine); - m_bkgEvents = make_unique(inputFiles); + m_bkgEvents = make_unique(inputFiles); for (auto& val : m_bkgEvents->m_totalNumberOfEvents) { if (val == 0) { std::string err = "No events found in the background files"; diff --git a/k4Reco/Overlay/components/OverlayTiming.h b/k4Reco/Overlay/components/OverlayTiming.h index b705ef1..9e1c385 100644 --- a/k4Reco/Overlay/components/OverlayTiming.h +++ b/k4Reco/Overlay/components/OverlayTiming.h @@ -52,28 +52,29 @@ #include #include -struct EventHolder { - std::vector> m_fileNames; - std::vector m_rootFileReaders; - std::vector m_totalNumberOfEvents; - std::map m_events; - - std::vector m_nextEntry; - - EventHolder(const std::vector>& fileNames) : m_fileNames(fileNames) { - for (auto& names : m_fileNames) { - m_rootFileReaders.emplace_back(podio::makeReader(names)); - m_totalNumberOfEvents.push_back(m_rootFileReaders.back().getEntries("events")); +namespace OverlayTimingNS { + struct EventHolder { + std::vector> m_fileNames; + std::vector m_rootFileReaders; + std::vector m_totalNumberOfEvents; + std::map m_events; + std::vector m_nextEntry; + + EventHolder(const std::vector>& fileNames) : m_fileNames(fileNames) { + for (auto& names : m_fileNames) { + m_rootFileReaders.emplace_back(podio::makeReader(names)); + m_totalNumberOfEvents.push_back(m_rootFileReaders.back().getEntries("events")); + } + m_nextEntry.resize(m_fileNames.size(), 0); } - m_nextEntry.resize(m_fileNames.size(), 0); - } - EventHolder() = default; + EventHolder() = default; - // TODO: Cache functionality - // podio::Frame& read + // TODO: Cache functionality + // podio::Frame& read - size_t size() const { return m_fileNames.size(); } -}; + size_t size() const { return m_fileNames.size(); } + }; +} using retType = std::tuple, @@ -138,7 +139,7 @@ struct OverlayTiming : public k4FWCore::MultiTransformer m_deltaT{this, "Delta_t", float(0.5), "Time difference between BXs in the BXtrain"}; - mutable std::unique_ptr m_bkgEvents{}; + mutable std::unique_ptr m_bkgEvents{}; Gaudi::Property>> m_timeWindows{ this, "TimeWindows", std::map>(), "Time windows for the different collections"}; diff --git a/k4Reco/Overlay/components/OverlayTimingRandomMix.cpp b/k4Reco/Overlay/components/OverlayTimingRandomMix.cpp index d0fbd92..6341b5a 100644 --- a/k4Reco/Overlay/components/OverlayTimingRandomMix.cpp +++ b/k4Reco/Overlay/components/OverlayTimingRandomMix.cpp @@ -91,7 +91,7 @@ StatusCode OverlayTimingRandomMix::initialize() { } } - m_bkgEvents = make_unique(inputFiles); + m_bkgEvents = make_unique(inputFiles); for (auto& valset : m_bkgEvents->m_totalNumberOfEvents) { for (auto& val : valset) { if (val == 0) { diff --git a/k4Reco/Overlay/components/OverlayTimingRandomMix.h b/k4Reco/Overlay/components/OverlayTimingRandomMix.h index afab73b..3fff3a1 100644 --- a/k4Reco/Overlay/components/OverlayTimingRandomMix.h +++ b/k4Reco/Overlay/components/OverlayTimingRandomMix.h @@ -53,34 +53,35 @@ #include #include -struct EventHolder { - std::vector> m_fileNames; - std::vector> m_totalNumberOfEvents; - - std::vector> m_nextEntry; - - EventHolder(const std::vector>& fileNames) : m_fileNames(fileNames) { - m_totalNumberOfEvents.resize(m_fileNames.size()); - m_nextEntry.resize(m_fileNames.size()); - - for (int group = 0; group < m_fileNames.size(); group++) { - m_nextEntry[group].resize(m_fileNames[group].size()); - for (auto& name : m_fileNames[group]) { - m_totalNumberOfEvents[group].push_back(1);//m_rootFileReaders[group].back().getEntries("events")); +namespace OverlayTimingRandomMixNS { + struct EventHolder { + std::vector> m_fileNames; + std::vector> m_totalNumberOfEvents; + std::vector> m_nextEntry; + + EventHolder(const std::vector>& fileNames) : m_fileNames(fileNames) { + m_totalNumberOfEvents.resize(m_fileNames.size()); + m_nextEntry.resize(m_fileNames.size()); + + for (int group = 0; group < m_fileNames.size(); group++) { + m_nextEntry[group].resize(m_fileNames[group].size()); + for (auto& name : m_fileNames[group]) { + m_totalNumberOfEvents[group].push_back(1);//m_rootFileReaders[group].back().getEntries("events")); + } } } - } - EventHolder() = default; + EventHolder() = default; - podio::Reader open(int groupIndex, int index) { - return podio::makeReader(m_fileNames[groupIndex][index]); - } + podio::Reader open(int groupIndex, int index) { + return podio::makeReader(m_fileNames[groupIndex][index]); + } - // TODO: Cache functionality - // podio::Frame& read + // TODO: Cache functionality + // podio::Frame& read - size_t size() const { return m_fileNames.size(); } -}; + size_t size() const { return m_fileNames.size(); } + }; +} using retType = std::tuple, @@ -135,7 +136,7 @@ struct OverlayTimingRandomMix : public k4FWCore::MultiTransformer m_deltaT{this, "Delta_t", float(0.5), "Time difference between BXs in the BXtrain"}; - mutable std::unique_ptr m_bkgEvents{}; + mutable std::unique_ptr m_bkgEvents{}; Gaudi::Property>> m_timeWindows{this, "TimeWindows", std::map>(), "Time windows for the different collections"}; Gaudi::Property m_allowReusingBackgroundFiles{this, "AllowReusingBackgroundFiles", false, "If true the same background file can be used for the same event"}; From 8dd7dda93749bf7c296d0977400a2a8c68b94687 Mon Sep 17 00:00:00 2001 From: samf25 Date: Mon, 25 Aug 2025 19:23:33 -0400 Subject: [PATCH 12/14] FastJet --- CMakeLists.txt | 4 +- cmake/FindFastJet.cmake | 29 +++ k4Reco/CMakeLists.txt | 5 + k4Reco/FastJet/include/EClusterMode.h | 35 +++ k4Reco/FastJet/include/FastJetAlg.hxx | 135 +++++++++++ k4Reco/FastJet/src/FastJetAlg.cpp | 336 ++++++++++++++++++++++++++ 6 files changed, 543 insertions(+), 1 deletion(-) create mode 100644 cmake/FindFastJet.cmake create mode 100644 k4Reco/FastJet/include/EClusterMode.h create mode 100644 k4Reco/FastJet/include/FastJetAlg.hxx create mode 100644 k4Reco/FastJet/src/FastJetAlg.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 63faae9..2241c56 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -26,6 +26,7 @@ set(${PROJECT_NAME}_VERSION_MINOR 3) set(${PROJECT_NAME}_VERSION_PATCH 0) set(${PROJECT_NAME}_VERSION "${${PROJECT_NAME}_VERSION_MAJOR}.${${PROJECT_NAME}_VERSION_MINOR}.${${PROJECT_NAME}_VERSION_PATCH}") +set(CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake ${CMAKE_MODULE_PATH}) option(BUILD_TRACKING "Build conformal tracking" ON) @@ -38,7 +39,7 @@ find_package(ROOT REQUIRED COMPONENTS RIO Tree) find_package(k4SimGeant4 REQUIRED) find_package(k4geo REQUIRED) find_package(CLHEP REQUIRED) - +find_package(FastJet REQUIRED) include(cmake/Key4hepConfig.cmake) include(GNUInstallDirs) @@ -52,3 +53,4 @@ endif() add_subdirectory(python) include(cmake/CreateProjectConfig.cmake) +include(cmake/FindFastJet.cmake) diff --git a/cmake/FindFastJet.cmake b/cmake/FindFastJet.cmake new file mode 100644 index 0000000..dedebbe --- /dev/null +++ b/cmake/FindFastJet.cmake @@ -0,0 +1,29 @@ +# - Locate FastJet library +# Defines: +# +# FASTJET_FOUND +# FASTJET_INCLUDE_DIR +# FASTJET_INCLUDE_DIRS (not cached) +# FASTJET_LIBRARY +# FASTJET_LIBRARIES (not cached) +# FASTJET_LIBRARY_DIRS (not cached) + +find_path(FASTJET_INCLUDE_DIR fastjet/version.hh + HINTS $ENV{FASTJET_ROOT_DIR}/include ${FASTJET_ROOT_DIR}/include) + +find_library(FASTJET_LIBRARY NAMES fastjet + HINTS $ENV{FASTJET_ROOT_DIR}/lib ${FASTJET_ROOT_DIR}/lib) + +find_library(FASTJETPLUGINS_LIBRARY NAMES fastjetplugins + HINTS $ENV{FASTJET_ROOT_DIR}/lib ${FASTJET_ROOT_DIR}/lib) + +# handle the QUIETLY and REQUIRED arguments and set FASTJET_FOUND to TRUE if +# all listed variables are TRUE +INCLUDE(FindPackageHandleStandardArgs) +FIND_PACKAGE_HANDLE_STANDARD_ARGS(FastJet DEFAULT_MSG FASTJET_INCLUDE_DIR FASTJET_LIBRARY) + +mark_as_advanced(FASTJET_FOUND FASTJET_INCLUDE_DIR FASTJET_LIBRARY) + +set(FASTJET_INCLUDE_DIRS ${FASTJET_INCLUDE_DIR}) +set(FASTJET_LIBRARIES ${FASTJET_LIBRARY} ${FASTJETPLUGINS_LIBRARY}) +get_filename_component(FASTJET_LIBRARY_DIRS ${FASTJET_LIBRARY} PATH) diff --git a/k4Reco/CMakeLists.txt b/k4Reco/CMakeLists.txt index a34cb7e..56be941 100644 --- a/k4Reco/CMakeLists.txt +++ b/k4Reco/CMakeLists.txt @@ -18,6 +18,7 @@ limitations under the License. ]] set(sources + FastJet/src/FastJetAlg.cpp DDPlanarDigi/components/DDPlanarDigi.cpp Overlay/components/OverlayTiming.cpp Overlay/components/OverlayTimingRandomMix.cpp @@ -53,6 +54,7 @@ set(libs ROOT::MathCore ROOT::Tree CLHEP::Random + ${FASTJET_LIBRARIES} ) if(BUILD_TRACKING) @@ -124,6 +126,9 @@ gaudi_add_module(k4RecoPlugins # Common include directories set(k4RecoPlugins_public_includes $ + $ + $ + ${FASTJET_INCLUDE_DIRS} ) target_include_directories(k4RecoPlugins PUBLIC ${k4RecoPlugins_public_includes}) diff --git a/k4Reco/FastJet/include/EClusterMode.h b/k4Reco/FastJet/include/EClusterMode.h new file mode 100644 index 0000000..75f2350 --- /dev/null +++ b/k4Reco/FastJet/include/EClusterMode.h @@ -0,0 +1,35 @@ +#ifndef EClusterMode_h +#define EClusterMode_h 1 + +#include + + +// The enum, name and value of the enum for the Cluster Mode +enum EClusterMode { + NONE = 0, + FJ_exclusive_yCut = 1, // exclusive clustering mode implemented in FastJet + FJ_exclusive_nJets = 2, // exclusive clustering mode implemented in FastJet + FJ_inclusive = 4, // inclusive "-" + OWN_inclusiveIteration = 8 // use FJ inclusive Clustering, but iterate until we have the desired number of jets +}; + +namespace { + +std::ostream& operator<<(std::ostream& out, EClusterMode& m) { + switch (m) { + case OWN_inclusiveIteration: + out << "InclusiveIterativeNJets"; break; + case FJ_inclusive: + out << "Inclusive"; break; + case FJ_exclusive_nJets: + out << "ExclusiveNJets"; break; + case FJ_exclusive_yCut: + out << "ExclusiveYCut"; break; + default: + out << "unknown"; break; + } + return out; +} + +} +#endif // EClusterMode_h \ No newline at end of file diff --git a/k4Reco/FastJet/include/FastJetAlg.hxx b/k4Reco/FastJet/include/FastJetAlg.hxx new file mode 100644 index 0000000..9e80872 --- /dev/null +++ b/k4Reco/FastJet/include/FastJetAlg.hxx @@ -0,0 +1,135 @@ +/* + * FastJetAlg.h + * + * Created on: 25.05.2010 + * Author: Lars Weuste (MPP Munich) - weuste@mpp.mpg.de + * iterative inclusive algorithm based on design by Marco Battaglia (CERN) - Marco.Battaglia@cern.ch + * Converted to Gaudi on: 25.08.2025 + * Conversion: Samuel Ferraro - samuel.rowles.ferraro@cern.ch + */ + +#ifndef FASTJETALG_H_ +#define FASTJETALG_H_ + +#include "EClusterMode.h" + +#include +#include + +//FastJet +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#define ITERATIVE_INCLUSIVE_MAX_ITERATIONS 20 + +//Forward declaration +typedef std::vector< fastjet::PseudoJet > PseudoJetList; + +class SkippedFixedNrJetException: public std::runtime_error { +public: + SkippedFixedNrJetException():std::runtime_error("") {} +}; + +class SkippedMaxIterationException: public std::runtime_error { +public: + SkippedMaxIterationException(PseudoJetList& jets) :std::runtime_error(""), m_jets(jets) {} + PseudoJetList m_jets; +}; + +struct FastJetAlg : k4FWCore::MultiTransformer< + std::tuple(const edm4hep::ReconstructedParticleCollection&)> { +public: + FastJetAlg(const std::string& name, ISvcLocator* svcLoc); + + /** Called at the begin of the job before anything is read. + * Use to initialize the processor, e.g. book histograms. + */ + StatusCode initialize(); + + /** Called for every run. + */ + std::tuple operator()( + const edm4hep::ReconstructedParticleCollection& inputCollection) const; + +private: + std::vector defaultJetAlgoNameAndParams{"kt_algorithm", "0.7"}; + std::vector defaultClusterMode{"Inclusive", "0.0"}; + + Gaudi::Property m_jetRecoSchemeName{this, "recombinationScheme", std::string("E_scheme"), "The recombination scheme used when merging 2 particles. Usually there is no need to use anything else than 4-Vector addition: E_scheme."}; + Gaudi::Property> m_jetAlgoNameAndParams{this, "algorithm", defaultJetAlgoNameAndParams, "Selects the algorithm and its parameters. E.g. 'kt_algorithm 0.7' or 'ee_kt_algorithm'. For a full list of supported algorithms, see the logfile after execution."}; + Gaudi::Property> m_clusterModeNameAndParam{this, "clusteringMode", defaultClusterMode, "One of 'Inclusive ', 'InclusiveIterativeNJets ', 'ExclusiveNJets ', 'ExclusiveYCut '. Note: not all modes are available for all algorithms."}; + + // jet algorithm + std::string m_jetAlgoName; + fastjet::JetDefinition* m_jetAlgo; + fastjet::JetAlgorithm m_jetAlgoType; + + // clustering mode + std::string m_clusterModeName; + EClusterMode m_clusterMode; + + // jet reco scheme + fastjet::RecombinationScheme m_jetRecoScheme; + + // jet strategy + std::string m_strategyName; + fastjet::Strategy m_strategy; + + // parameters + unsigned m_requestedNumberOfJets; + double m_yCut; + double m_minPt; + double m_minE; + +private: + bool isJetAlgo(std::string algo, int nrParams, int supportedModes) const; +}; + +DECLARE_COMPONENT(FastJetAlg) + +bool FastJetAlg::isJetAlgo(std::string algo, int nrParams, int supportedModes) const { + info() << " " << algo; + + // check if the chosen algorithm is the same as it was passed + if (m_jetAlgoName.compare(algo) != 0) { + return false; + } + info() << "*"; // mark the current algorithm as the selected + // one, even before we did our checks on nr of + // parameters etc. + + // check if we have enough number of parameters + if ((int)m_jetAlgoNameAndParams.size() - 1 != nrParams) { + error() << std::endl + << "Wrong numbers of parameters for algorithm: " << algo << std::endl + << "We need " << nrParams << " params, but we got " << m_jetAlgoNameAndParams.size() - 1 << std::endl; + throw GaudiException("You have insufficient number of parameters for this algorithm! See log for more details.", name(), StatusCode::FAILURE); + } + + // check if the mode is supported via a binary AND + if ((supportedModes & m_clusterMode) != m_clusterMode) { + error() << std::endl + << "This algorithm is not capable of running in this clustering mode (" + << m_clusterMode << "). Sorry!" << std::endl; + throw GaudiException("This algorithm is not capable of running in this mode", name(), StatusCode::FAILURE); + } + + return true; +} + + +#endif /* FASTJETALG_H_ */ \ No newline at end of file diff --git a/k4Reco/FastJet/src/FastJetAlg.cpp b/k4Reco/FastJet/src/FastJetAlg.cpp new file mode 100644 index 0000000..3c5db90 --- /dev/null +++ b/k4Reco/FastJet/src/FastJetAlg.cpp @@ -0,0 +1,336 @@ +/* + * FastJetAlg.cpp + * + * Created on: 25.05.2010 + * Author: Lars Weuste (MPP Munich) - weuste@mpp.mpg.de + * iterative inclusive algorithm based on design by Marco Battaglia (CERN) - Marco.Battaglia@cern.ch + * Converted to Gaudi on: 25.08.2025 + * Conversion: Samuel Ferraro - samuel.rowles.ferraro@cern.ch + */ + +#include "FastJetAlg.hxx" + +#include +#include +#include +#include +#include + +#include + +FastJetAlg::FastJetAlg(const std::string& name, ISvcLocator* svcLoc) : MultiTransformer(name, svcLoc, + { KeyValues("recParticleIn", {"MCParticle"}) }, + { KeyValues("jetOut", {"JetOut"}), + KeyValues("recParticleOut", {"Constituents"}) }), + m_jetAlgoName(""), + m_jetAlgo(NULL), + m_jetAlgoType(), + m_clusterModeName(""), + m_clusterMode( NONE ), + m_jetRecoScheme(), + m_strategyName(""), + m_strategy(), + m_requestedNumberOfJets(0), + m_yCut(0.0), + m_minPt(0.0), + m_minE(0.0) {} + +StatusCode FastJetAlg::initialize() { + // ------------------ Init Strategy ------------------ + m_strategy = fastjet::Best; + m_strategyName = "Best"; + info() << "Strategy: " << m_strategyName << endmsg; + + // ------------------ Init Reco Scheme ------------------ + if (m_jetRecoSchemeName.value().compare("E_scheme") == 0) + m_jetRecoScheme = fastjet::E_scheme; + else if (m_jetRecoSchemeName.value().compare("pt_scheme") == 0) + m_jetRecoScheme = fastjet::pt_scheme; + else if (m_jetRecoSchemeName.value().compare("pt2_scheme") == 0) + m_jetRecoScheme = fastjet::pt2_scheme; + else if (m_jetRecoSchemeName.value().compare("Et_scheme") == 0) + m_jetRecoScheme = fastjet::Et_scheme; + else if (m_jetRecoSchemeName.value().compare("Et2_scheme") == 0) + m_jetRecoScheme = fastjet::Et2_scheme; + else if (m_jetRecoSchemeName.value().compare("BIpt_scheme") == 0) + m_jetRecoScheme = fastjet::BIpt_scheme; + else if (m_jetRecoSchemeName.value().compare("BIpt2_scheme") == 0) + m_jetRecoScheme = fastjet::BIpt2_scheme; + else { + error() << "Unknown recombination scheme: " << m_jetRecoSchemeName << endmsg; + throw GaudiException("Unknown FastJet recombination scheme! See log for more details.", name(), StatusCode::FAILURE); + } + info() << "recombination scheme: " << m_jetRecoSchemeName << endmsg; + + // ------------------ Init Cluster Mode ------------------ + // at least a name has to be given + if (m_clusterModeNameAndParam.size() == 0) + throw GaudiException("Cluster mode not specified", name(), StatusCode::FAILURE); + // save the name of the cluster mode + m_clusterModeName = m_clusterModeNameAndParam[0]; + m_clusterMode = NONE; + // check the different cluster mode possibilities, and check if the number of parameters are correct + if (m_clusterModeName.compare("Inclusive") == 0) { + if (m_clusterModeNameAndParam.size() != 2) { + error() << "Wrong number of values for parameter clusteringMode 'Inclusive': missing minPt" << endmsg; + throw GaudiException("Wrong Parameter(s) for Clustering Mode. Expected:\n Inclusive ", name(), StatusCode::FAILURE); + } + m_minPt = atof(m_clusterModeNameAndParam[1].c_str()); + m_clusterMode = FJ_inclusive; + } else if (m_clusterModeName.compare("InclusiveIterativeNJets") == 0) { + if (m_clusterModeNameAndParam.size() != 3) { + throw GaudiException("Wrong Parameter(s) for Clustering Mode. Expected:\n InclusiveIterativeNJets ", name(), StatusCode::FAILURE); + } + m_requestedNumberOfJets = atoi(m_clusterModeNameAndParam[1].c_str()); + m_minE = atoi(m_clusterModeNameAndParam[2].c_str()); + m_clusterMode = OWN_inclusiveIteration; + } else if (m_clusterModeName.compare("ExclusiveNJets") == 0) { + if (m_clusterModeNameAndParam.size() != 2) { + throw GaudiException("Wrong Parameter(s) for Clustering Mode. Expected:\n ExclusiveNJets ", name(), StatusCode::FAILURE); + } + m_requestedNumberOfJets = atoi(m_clusterModeNameAndParam[1].c_str()); + m_clusterMode = FJ_exclusive_nJets; + } else if (m_clusterModeName.compare("ExclusiveYCut") == 0) { + if (m_clusterModeNameAndParam.size() != 2) { + throw GaudiException("Wrong Parameter(s) for Clustering Mode. Expected:\n ExclusiveYCut ", name(), StatusCode::FAILURE); + } + m_yCut = atof(m_clusterModeNameAndParam[1].c_str()); + m_clusterMode = FJ_exclusive_yCut; + } else { + throw GaudiException("Unknown cluster mode.", name(), StatusCode::FAILURE); + } + info() << "Cluster mode: " << m_clusterMode << endmsg; + + // ------------------ Init Jet Algorithm ------------------ + // sanity check + if (m_jetAlgoNameAndParams.size() < 1) + throw GaudiException("No Jet algorithm provided!", name(), StatusCode::FAILURE); + // save the name + m_jetAlgoName = m_jetAlgoNameAndParams[0]; + // check all supported algorithms and create the appropriate FJ instance + m_jetAlgo = NULL; + info() << "Algorithms: "; // the isJetAlgo function will write to streamlog_out(MESSAGE), so that we get a list of available algorithms in the log + // example: kt_algorithm, needs 1 parameter, supports inclusive, inclusiveIterative, exlusiveNJets and exlusiveYCut clustering + if (isJetAlgo("kt_algorithm", 1, FJ_inclusive | FJ_exclusive_nJets | FJ_exclusive_yCut | OWN_inclusiveIteration)) { + m_jetAlgoType = fastjet::kt_algorithm; + m_jetAlgo = new fastjet::JetDefinition( + m_jetAlgoType, atof(m_jetAlgoNameAndParams[1].c_str()), m_jetRecoScheme, m_strategy); + } + if (isJetAlgo("cambridge_algorithm", 1, FJ_inclusive | FJ_exclusive_nJets | FJ_exclusive_yCut | OWN_inclusiveIteration)) { + m_jetAlgoType = fastjet::cambridge_algorithm; + m_jetAlgo = new fastjet::JetDefinition( + m_jetAlgoType, atof(m_jetAlgoNameAndParams[1].c_str()), m_jetRecoScheme, m_strategy); + } + if (isJetAlgo("antikt_algorithm", 1, FJ_inclusive | OWN_inclusiveIteration)) { + m_jetAlgoType = fastjet::antikt_algorithm; + m_jetAlgo = new fastjet::JetDefinition( + m_jetAlgoType, atof(m_jetAlgoNameAndParams[1].c_str()), m_jetRecoScheme, m_strategy); + } + if (isJetAlgo("genkt_algorithm", 2, FJ_inclusive | OWN_inclusiveIteration | FJ_exclusive_nJets | FJ_exclusive_yCut)) { + m_jetAlgoType = fastjet::genkt_algorithm; + m_jetAlgo = new fastjet::JetDefinition( + m_jetAlgoType, atof(m_jetAlgoNameAndParams[1].c_str()), atof(m_jetAlgoNameAndParams[2].c_str()), m_jetRecoScheme, m_strategy); + } + if (isJetAlgo("cambridge_for_passive_algorithm", 1, FJ_inclusive | OWN_inclusiveIteration | FJ_exclusive_nJets | FJ_exclusive_yCut)) { + m_jetAlgoType = fastjet::cambridge_for_passive_algorithm; + m_jetAlgo = new fastjet::JetDefinition( + m_jetAlgoType, atof(m_jetAlgoNameAndParams[1].c_str()), m_jetRecoScheme, m_strategy); + } + if (isJetAlgo("genkt_for_passive_algorithm", 1, FJ_inclusive | OWN_inclusiveIteration)) { + m_jetAlgoType = fastjet::genkt_for_passive_algorithm; + m_jetAlgo = new fastjet::JetDefinition( + m_jetAlgoType, atof(m_jetAlgoNameAndParams[1].c_str()), m_jetRecoScheme, m_strategy); + } + if (isJetAlgo("ee_kt_algorithm", 0, FJ_exclusive_nJets | FJ_exclusive_yCut)) { + m_jetAlgoType = fastjet::ee_kt_algorithm; + m_jetAlgo = new fastjet::JetDefinition( + m_jetAlgoType, m_jetRecoScheme, m_strategy); + } + // backwards compatibility for using 1 parameter only assuming exponent to be 1. + bool commentOnAlgo = false; + if ((m_jetAlgoNameAndParams[0]=="ee_genkt_algorithm") && ((int)m_jetAlgoNameAndParams.size() == 2)){ + m_jetAlgoNameAndParams.value().push_back("1."); + commentOnAlgo = true; + } + if (isJetAlgo("ee_genkt_algorithm", 2, FJ_inclusive | FJ_exclusive_nJets | FJ_exclusive_yCut)) { + m_jetAlgoType = fastjet::ee_genkt_algorithm; + m_jetAlgo = new fastjet::JetDefinition( + m_jetAlgoType, atof(m_jetAlgoNameAndParams[1].c_str()), atof(m_jetAlgoNameAndParams[2].c_str()), m_jetRecoScheme, m_strategy); + } + if (isJetAlgo("SISConePlugin", 2, FJ_inclusive | OWN_inclusiveIteration)) { + fastjet::SISConePlugin* pl; + pl = new fastjet::SISConePlugin( + atof(m_jetAlgoNameAndParams[1].c_str()), + atof(m_jetAlgoNameAndParams[2].c_str()) + ); + m_jetAlgo = new fastjet::JetDefinition(pl); + m_jetAlgo->delete_plugin_when_unused(); + } + if (isJetAlgo("SISConeSphericalPlugin", 2, FJ_inclusive | OWN_inclusiveIteration)) { + fastjet::SISConeSphericalPlugin* pl; + pl = new fastjet::SISConeSphericalPlugin( + atof(m_jetAlgoNameAndParams[1].c_str()), + atof(m_jetAlgoNameAndParams[2].c_str()) + ); + m_jetAlgo = new fastjet::JetDefinition(pl); + m_jetAlgo->delete_plugin_when_unused(); + } + if (isJetAlgo("ValenciaPlugin", 3, FJ_exclusive_nJets | FJ_exclusive_yCut)) { + fastjet::contrib::ValenciaPlugin* pl; + pl = new fastjet::contrib::ValenciaPlugin( + atof(m_jetAlgoNameAndParams[1].c_str()), // R value + atof(m_jetAlgoNameAndParams[2].c_str()), // beta value + atof(m_jetAlgoNameAndParams[3].c_str()) // gamma value + ); + m_jetAlgo = new fastjet::JetDefinition(pl); + m_jetAlgo->delete_plugin_when_unused(); + } + info() << endmsg; // end of list of available algorithms + //ee_genkt_algorithm + if (commentOnAlgo) {info() << "When only 1 parameter is provided for ee_genkt_algorithm it is assumed to be R, and the exponent p is assumed to be equal to 1" << endmsg;} + if (!m_jetAlgo) { + error() << "The given algorithm \"" << m_jetAlgoName << "\" is unknown to me!" << endmsg; + throw GaudiException("Unknown FastJet algorithm.", name(), StatusCode::FAILURE); + } + info() << "jet algorithm: " << m_jetAlgo->description() << endmsg; + + return StatusCode::SUCCESS; +} + +std::tuple FastJetAlg::operator()( + const edm4hep::ReconstructedParticleCollection& inputCollection) const { + edm4hep::ReconstructedParticleCollection outputCollection; + outputCollection->setSubsetCollection(true); + edm4hep::ReconstructedParticleCollection jetCollection; + + PseudoJetList jets; + PseudoJetList pjList; + for (int i = 0; i < inputCollection.size(); ++i) { + edm4hep::ReconstructedParticle par = inputCollection.at(i); + pjList.push_back( fastjet::PseudoJet( + par.getMomentum().x, + par.getMomentum().y, + par.getMomentum().z, + par.getEnergy() ) ); + pjList.back().set_user_index(i); // save the id of this recParticle + } + + fastjet::ClusterSequence cs = fastjet::ClusterSequence(pjList, *m_jetAlgo); + try { + if (m_clusterMode == FJ_inclusive) { + jets = cs.inclusive_jets(m_minPt); + } else if (m_clusterMode == FJ_exclusive_yCut) { + jets = cs.exclusive_jets_ycut(m_yCut); + } else if (m_clusterMode == FJ_exclusive_nJets) { + // sanity check: if we have not enough particles, FJ will cause an assert + if (inputCollection.size() < (int)m_requestedNumberOfJets) { + warning() << "Not enough elements in the input collection to create " << m_requestedNumberOfJets << " jets." << endmsg; + throw SkippedFixedNrJetException(); + } else { + jets = cs.exclusive_jets((int)(m_requestedNumberOfJets)); + } + } else if (m_clusterMode == OWN_inclusiveIteration) { + // sanity check: if we have not enough particles, FJ will cause an assert + if (inputCollection.size() < (int)m_requestedNumberOfJets) { + warning() << "Not enough elements in the input collection to create " << m_requestedNumberOfJets << " jets." << endmsg; + throw SkippedFixedNrJetException(); + } else { + // lets do a iterative procedure until we found the correct number of jets + // for that we will do inclusive clustering, modifying the R parameter in some kind of minimization + // this is based on Marco Battaglia's FastJetClustering + double R = M_PI_4; // maximum of R is Pi/2, minimum is 0. So we start hat Pi/4 + double RDiff = R / 2; // the step size we modify the R parameter at each iteration. Its size for the n-th step is R/(2n), i.e. starts with R/2 + PseudoJetList jets_it; + unsigned nJets; + int iIter = 0; // nr of current iteration + // these variables are only used if the SisCone(Spherical)Plugin is selected + // This is necessary, as these are plugins and hence use a different constructor than + // the built in fastjet algorithms + // here we save pointer to the plugins, so that if they are created, we can delete them again after usage + fastjet::SISConePlugin* pluginSisCone = NULL; + fastjet::SISConeSphericalPlugin* pluginSisConeSph = NULL; + // check if we use the siscones + bool useSisCone = m_jetAlgoName.compare("SISConePlugin") == 0; + bool useSisConeSph = m_jetAlgoName.compare("SISConeSphericalPlugin") == 0; + // save the 2nd parameter of the SisCones + double sisConeOverlapThreshold = 0; + if (useSisCone || useSisConeSph) + sisConeOverlapThreshold = atof(m_jetAlgoNameAndParams[2].c_str()); + // do a maximum of N iterations + for (iIter=0; iIter m_minE) + jets.push_back(jets_it[j]); + nJets = jets.size(); + debug() << iIter << " " << R << " " << jets_it.size() << " " << nJets << endmsg; + if (nJets == m_requestedNumberOfJets) { // if the number of jets is correct: success! + delete pluginSisCone; pluginSisCone = NULL; + delete pluginSisConeSph; pluginSisConeSph = NULL; + delete jetDefinition; + break; + } else if (nJets < m_requestedNumberOfJets) { + // if number of jets is too small: we need a smaller Radius per jet (so + // that we get more jets) + R -= RDiff; + } else if (nJets > m_requestedNumberOfJets) { // if the number of jets is too + // high: increase the Radius + R += RDiff; + } + RDiff /= 2; + // clean up + delete pluginSisCone; pluginSisCone = NULL; + delete pluginSisConeSph; pluginSisConeSph = NULL; + delete jetDefinition; + } + if (iIter == ITERATIVE_INCLUSIVE_MAX_ITERATIONS) { + warning() << "Maximum number of iterations reached. Canceling" << endmsg; + throw SkippedMaxIterationException( jets ); + // Currently we will return the latest results, independent if the number is actually matched + // jets.clear(); + } + } + } + } catch(const SkippedFixedNrJetException& e ) { + } catch(const SkippedMaxIterationException& e ) { + jets = e.m_jets; + } + + PseudoJetList::iterator it; + for (it=jets.begin(); it != jets.end(); it++) { + // create a reconstructed particle for this jet, and add all the containing particles to it + edm4hep::MutableReconstructedParticle rec = jetCollection->create(); + rec.setEnergy( (*it).E() ); + rec.setMass( (*it).m() ); + edm4hep::Vector3f mom((*it).px(), (*it).py(), (*it).pz()); + rec.setMomentum(mom); + for (unsigned int n = 0; n < cs.constituents(*it).size(); ++n) { + rec.addToParticles(inputCollection.at(cs.constituents(*it)[n].user_index())); + } + + // add jet constituents to output collection + for (unsigned int n = 0; n < cs.constituents(*it).size(); ++n) { + edm4hep::ReconstructedParticle p = inputCollection.at((cs.constituents(*it))[n].user_index()); + outputCollection.push_back( p ); + } + } + + return std::make_tuple(std::move(outputCollection), std::move(jetCollection)); +} \ No newline at end of file From 8140955d45cbb2fa706f1c58516eac755f7b5ef6 Mon Sep 17 00:00:00 2001 From: samf25 Date: Tue, 26 Aug 2025 11:56:46 -0400 Subject: [PATCH 13/14] Remove Valencia so that it compiles --- k4Reco/CMakeLists.txt | 6 ++++++ k4Reco/FastJet/include/FastJetAlg.hxx | 2 +- k4Reco/FastJet/src/FastJetAlg.cpp | 21 +++++++++++---------- 3 files changed, 18 insertions(+), 11 deletions(-) diff --git a/k4Reco/CMakeLists.txt b/k4Reco/CMakeLists.txt index 56be941..345bf87 100644 --- a/k4Reco/CMakeLists.txt +++ b/k4Reco/CMakeLists.txt @@ -42,6 +42,12 @@ set(sources CaloDigi/src/RealisticCaloReco.cc ) +if(DEFINED FASTJETPLUGINS_LIBRARY) + message(STATUS "FASTJETPLUGINS_LIBRARY: ${FASTJETPLUGINS_LIBRARY}") +else() + message(WARNING "FASTJETPLUGINS_LIBRARY is not defined") +endif() + set(libs podio::podioIO EDM4HEP::edm4hep diff --git a/k4Reco/FastJet/include/FastJetAlg.hxx b/k4Reco/FastJet/include/FastJetAlg.hxx index 9e80872..1dcd149 100644 --- a/k4Reco/FastJet/include/FastJetAlg.hxx +++ b/k4Reco/FastJet/include/FastJetAlg.hxx @@ -27,7 +27,7 @@ #include #include #include -#include +// #include #include #include diff --git a/k4Reco/FastJet/src/FastJetAlg.cpp b/k4Reco/FastJet/src/FastJetAlg.cpp index 3c5db90..c659ec9 100644 --- a/k4Reco/FastJet/src/FastJetAlg.cpp +++ b/k4Reco/FastJet/src/FastJetAlg.cpp @@ -175,16 +175,17 @@ StatusCode FastJetAlg::initialize() { m_jetAlgo = new fastjet::JetDefinition(pl); m_jetAlgo->delete_plugin_when_unused(); } - if (isJetAlgo("ValenciaPlugin", 3, FJ_exclusive_nJets | FJ_exclusive_yCut)) { - fastjet::contrib::ValenciaPlugin* pl; - pl = new fastjet::contrib::ValenciaPlugin( - atof(m_jetAlgoNameAndParams[1].c_str()), // R value - atof(m_jetAlgoNameAndParams[2].c_str()), // beta value - atof(m_jetAlgoNameAndParams[3].c_str()) // gamma value - ); - m_jetAlgo = new fastjet::JetDefinition(pl); - m_jetAlgo->delete_plugin_when_unused(); - } + // *********** VALENCIA WAS NOT IN THE FASTJET INSTALLATION SO IT IS NOT ACTIVE FOR NOW *********** + // if (isJetAlgo("ValenciaPlugin", 3, FJ_exclusive_nJets | FJ_exclusive_yCut)) { + // fastjet::contrib::ValenciaPlugin* pl; + // pl = new fastjet::contrib::ValenciaPlugin( + // atof(m_jetAlgoNameAndParams[1].c_str()), // R value + // atof(m_jetAlgoNameAndParams[2].c_str()), // beta value + // atof(m_jetAlgoNameAndParams[3].c_str()) // gamma value + // ); + // m_jetAlgo = new fastjet::JetDefinition(pl); + // m_jetAlgo->delete_plugin_when_unused(); + // } info() << endmsg; // end of list of available algorithms //ee_genkt_algorithm if (commentOnAlgo) {info() << "When only 1 parameter is provided for ee_genkt_algorithm it is assumed to be R, and the exponent p is assumed to be equal to 1" << endmsg;} From 7cd9f43cc0362c8ce7334ab531c94b6e7719e175 Mon Sep 17 00:00:00 2001 From: Andrew David Loeliger Date: Mon, 2 Mar 2026 23:24:54 +0100 Subject: [PATCH 14/14] Fix compilation errors --- k4Reco/FastJet/src/FastJetAlg.cpp | 6 +++--- .../Overlay/components/OverlayTimingRandomMix.cpp | 14 +++++++------- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/k4Reco/FastJet/src/FastJetAlg.cpp b/k4Reco/FastJet/src/FastJetAlg.cpp index c659ec9..d665483 100644 --- a/k4Reco/FastJet/src/FastJetAlg.cpp +++ b/k4Reco/FastJet/src/FastJetAlg.cpp @@ -201,7 +201,7 @@ StatusCode FastJetAlg::initialize() { std::tuple FastJetAlg::operator()( const edm4hep::ReconstructedParticleCollection& inputCollection) const { edm4hep::ReconstructedParticleCollection outputCollection; - outputCollection->setSubsetCollection(true); + outputCollection.setSubsetCollection(true); edm4hep::ReconstructedParticleCollection jetCollection; PseudoJetList jets; @@ -317,7 +317,7 @@ std::tuplecreate(); + edm4hep::MutableReconstructedParticle rec = jetCollection.create(); rec.setEnergy( (*it).E() ); rec.setMass( (*it).m() ); edm4hep::Vector3f mom((*it).px(), (*it).py(), (*it).pz()); @@ -334,4 +334,4 @@ std::tuplepush_back(part.clone(false)); + oparticles.push_back(part.clone(false)); } // Fix relations to point to the new particles for (size_t i = 0; i < particles.size(); ++i) { @@ -171,7 +171,7 @@ retType OverlayTimingRandomMix::operator()( auto nhit = simTrackerHit.clone(false); if (simTrackerHit.getParticle().getObjectID().index != -1) nhit.setParticle(oparticles[simTrackerHit.getParticle().getObjectID().index]); - ocoll->push_back(nhit); + ocoll.push_back(nhit); } } osimTrackerHits.emplace_back(std::move(ocoll)); @@ -291,7 +291,7 @@ retType OverlayTimingRandomMix::operator()( npart.setTime(bgParticles[i].getTime() + timeOffset); npart.setOverlay(true); - oparticles->push_back(npart); + oparticles.push_back(npart); for (const auto& parent : bgParticles[i].getParents()) { parentDaughterMap[j].first.push_back(parent.getObjectID().index); } @@ -354,7 +354,7 @@ retType OverlayTimingRandomMix::operator()( nhit.setMomentum({(float)mom.x, (float)mom.y, (float)mom.z}); } } - ocoll->push_back(nhit); + ocoll.push_back(nhit); } } @@ -429,17 +429,17 @@ retType OverlayTimingRandomMix::operator()( for (const auto& [index, calHitMap] : cellIDsMap) { auto ocoll = edm4hep::SimCalorimeterHitCollection(); for (const auto& [cellID, hit] : calHitMap) { - ocoll->push_back(std::move(hit)); + ocoll.push_back(std::move(hit)); } osimCaloHits.emplace_back(std::move(ocoll)); } debug() << "\n\t\tCollection\t\t|\t\tPre BIB\t\t|\t\tPost BIB\t\t\n--------------------------------------------------------------------------\n"; for (int trkCol = 0; trkCol < simTrackerHits.size(); trkCol++) { - debug() << "\tTrackerHits " << trkCol << "\t|\t\t" << simTrackerHits[trkCol]->size() << "\t\t|\t\t" << osimTrackerHits[trkCol]->size() << "\n"; + debug() << "\tTrackerHits " << trkCol << "\t|\t\t" << simTrackerHits[trkCol]->size() << "\t\t|\t\t" << osimTrackerHits[trkCol].size() << "\n"; } for (int calCol = 0; calCol < simCaloHits.size(); calCol++) { - debug() << "\tCaloHits " << calCol << "\t|\t\t" << simCaloHits[calCol]->size() << "\t\t|\t\t" << osimCaloHits[calCol]->size() << "\n"; + debug() << "\tCaloHits " << calCol << "\t|\t\t" << simCaloHits[calCol]->size() << "\t\t|\t\t" << osimCaloHits[calCol].size() << "\n"; } debug() << endmsg;