Rebase of muon collider fork

CMakeLists.txt

@@ -26,17 +26,20 @@ 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)
 
-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)
 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)
+find_package(FastJet REQUIRED)
 include(cmake/Key4hepConfig.cmake)
 
 include(GNUInstallDirs)
@@ -50,3 +53,4 @@ endif()
 add_subdirectory(python)
 
 include(cmake/CreateProjectConfig.cmake)
+include(cmake/FindFastJet.cmake)

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)

k4Reco/CMakeLists.txt

@@ -18,8 +18,14 @@ limitations under the License.
 ]]
 
 set(sources
+  FastJet/src/FastJetAlg.cpp
   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
@@ -32,8 +38,16 @@ set(sources
   GaudiLumiCalClusterer/src/LCCluster.cpp
   GaudiLumiCalClusterer/src/SuperTrueClusterWeights.cpp
   GaudiLumiCalClusterer/src/ProjectionInfo.cpp
+  CaloDigi/src/RealisticCaloDigi.cc
+  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
@@ -45,6 +59,8 @@ set(libs
         ROOT::Core
         ROOT::MathCore
         ROOT::Tree
+        CLHEP::Random
+        ${FASTJET_LIBRARIES}
 )
 
 if(BUILD_TRACKING)
@@ -97,6 +113,7 @@ if(BUILD_TRACKING)
        Tracking/src/ClonesAndSplitTracksFinder.cpp
        Tracking/src/RefitFinal.cpp
        Tracking/src/TruthTrackFinder.cpp
+       Tracking/src/FilterDoubleLayerHits.cc
   )
 
   list(APPEND libs
@@ -115,13 +132,17 @@ gaudi_add_module(k4RecoPlugins
 # Common include directories
 set(k4RecoPlugins_public_includes
   $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>
+  $<BUILD_INTERFACE:${CMAKE_CURRENT_LIST_DIR}/CaloDigi/include>
+  $<BUILD_INTERFACE:${CMAKE_CURRENT_LIST_DIR}/FastJet/include>
+  ${FASTJET_INCLUDE_DIRS}
 )
 
 target_include_directories(k4RecoPlugins PUBLIC ${k4RecoPlugins_public_includes})
 target_include_directories(k4RecoPlugins PUBLIC
   $<BUILD_INTERFACE:${CMAKE_CURRENT_LIST_DIR}/ConformalTracking/include>
   $<BUILD_INTERFACE:${CMAKE_CURRENT_LIST_DIR}/GaudiLumiCalClusterer/include>
   $<BUILD_INTERFACE:${CMAKE_CURRENT_LIST_DIR}/Tracking/include>
+  $<BUILD_INTERFACE:${CMAKE_CURRENT_LIST_DIR}/CaloDigi/include>
   $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>)
 
 if(BUILD_TRACKING)

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 <iostream>
+
+/** Helper class for decoding/encoding lcio::CalorimeterHit types for the ILD
+ *  detector. The encoding is: caloType + 10 * caloID + 1000 * layout + 10000 * layerNum <br>
+ *  (see enums: CaloType, CaloID and Layout for possible values).<br>
+ *  Example usage: <br>
+ *  <pre>
+ *     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 ;
+ *
+ *  </pre>
+ *
+ *  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

k4Reco/CaloDigi/include/RealisticCaloDigi.h

@@ -0,0 +1,118 @@
+#ifndef REALISTICCALODIGI_H
+#define REALISTICCALODIGI_H 1
+
+#include <k4FWCore/Transformer.h>
+#include <edm4hep/SimCalorimeterHit.h>
+#include <edm4hep/SimCalorimeterHitCollection.h>
+#include <edm4hep/CalorimeterHitCollection.h>
+#include <edm4hep/EventHeaderCollection.h>
+#include <edm4hep/CaloHitSimCaloHitLinkCollection.h>
+
+#include "k4Interface/IGeoSvc.h"
+#include "k4Interface/IUniqueIDGenSvc.h"
+
+#include "TRandom2.h"
+
+#include <string>
+#include <tuple>
+#include <vector>
+#include <functional>
+#include <optional>
+
+
+/** === RealisticCaloDigi Processor === <br>
+    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<edm4hep::CalorimeterHitCollection, 
+    edm4hep::CaloHitSimCaloHitLinkCollection>(
+      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<edm4hep::CalorimeterHitCollection, 
+               edm4hep::CaloHitSimCaloHitLinkCollection> 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<float,float>;
+    using integr_res_opt = std::optional<integr_res>;
+    using integr_function = std::function<integr_res_opt(const edm4hep::SimCalorimeterHit*)>;
+
+   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<int> m_time_apply{this, "timingCut", 0, "Use hit times"};
+   Gaudi::Property<int> m_time_correctForPropagation{this, "timingCorrectForPropagation", 0, "Correct hit times for propagation: radial distance/c"};
+   Gaudi::Property<float> m_time_windowMin{this, "timingWindowMin", -10.0f, "Time Window minimum time in ns"};
+   Gaudi::Property<float> m_time_windowMax{this, "timingWindowMax", 100.0f, "Time Window maximum time in ns"};
+   Gaudi::Property<std::string> m_integration_method{this, "integrationMethod", "Standard", "Energy integration and time calculation method. Options: Standard, ROC"};
+   Gaudi::Property<float> m_fast_shaper{this, "fastShaper", 0.f, "Fast shaper value. Unit in ns"};
+   Gaudi::Property<float> m_slow_shaper{this, "slowShaper", 0.f, "Slow shaper value. Unit in ns"};
+   Gaudi::Property<float> m_time_resol{this, "timingResolution", 0.f, "Time resolution to apply (gaussian smearing). Unit in ns"};
+   // additional digi effects
+   Gaudi::Property<float> m_calib_mip{this, "calibration_mip", 1.0e-4f, "Average G4 deposited energy by MIP for calibration"};
+   Gaudi::Property<float> m_misCalib_uncorrel{this, "miscalibration_uncorrel", 0.0f, "Uncorrelated random Gaussian miscalibration (as a fraction: 1.0 = 100%)"};
+   Gaudi::Property<float> m_misCalib_correl{this, "miscalibration_correl", 0.0f, "Correlated random Gaussian miscalibration (as a fraction: 1.0 = 100%)"};
+   Gaudi::Property<float> m_deadCell_fraction{this, "deadCell_fraction", 0.0f, "Random dead cell fraction (as a fraction: 0->1)"};
+   // simple model of electronics properties
+   Gaudi::Property<float> m_elec_noiseMip{this, "elec_noise_mip", 0.0f, "Typical electronics noise (in MIP units)"};
+   Gaudi::Property<float> 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<std::string> m_encodingStringVariable{this, "EncodingStringParameterName", "GlobalCalorimeterReadoutID", "The name of the DD4hep constant that contains the Encoding string for tracking detectors"};
+   // energy threshold
+   Gaudi::Property<float> m_threshold_value{this, "threshold", 0.5f, "Threshold for Hit"};
+   Gaudi::Property<std::string> 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<std::string> m_calo_type {this, "CaloType", "em", "Calorimeter Type: em, had, mu"};
+   Gaudi::Property<std::string> m_calo_id {this, "CaloID", "ecal", "Calorimeter ID: ecal, hcal, yoke, lcal, lhcal, bcal"};
+   Gaudi::Property<std::string> m_calo_layout {this, "CaloLayout", "barrel", "Calorimeter Layout: barrel, endcap, ring, plug"};
+
+
+   int m_threshold_iunit{};  
+   inline static thread_local TRandom2 m_engine;
+   SmartIF<IGeoSvc>                    m_geoSvc;
+   SmartIF<IUniqueIDGenSvc>            m_uidSvc;
+
+
+   integr_function m_integr_function{};
+
+} ;
+
+#endif
+
+
+