diff --git a/CMakeLists.txt b/CMakeLists.txt index c1b02b9..c1fbae7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -38,6 +38,7 @@ find_package(PandoraSDK) find_package(LCContent) find_package(Gaudi) find_package(DD4hep) +find_package(CLHEP REQUIRED) #--------------------------------------------------------------- @@ -53,6 +54,37 @@ endif() include(CTest) + +#---RPATH options------------------------------------------------------------------------------- +option(K4GAUDIPANDORA_SET_RPATH "Link libraries with built-in RPATH (run-time search path)" ON) + +# When building, don't use the install RPATH already (but later on when installing) +set(CMAKE_SKIP_BUILD_RPATH FALSE) # don't skip the full RPATH for the build tree +set(CMAKE_BUILD_WITH_INSTALL_RPATH FALSE) # use always the build RPATH for the build tree +set(CMAKE_MACOSX_RPATH TRUE) # use RPATH for MacOSX +set(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE) # point to directories outside the build tree to the install RPATH + +# Check whether to add RPATH to the installation (the build tree always has the RPATH enabled) +if(APPLE) + set(CMAKE_INSTALL_NAME_DIR "@rpath") + set(CMAKE_INSTALL_RPATH "@loader_path/../lib") # self relative LIBDIR + # the RPATH to be used when installing, but only if it's not a system directory + list(FIND CMAKE_PLATFORM_IMPLICIT_LINK_DIRECTORIES "${CMAKE_INSTALL_PREFIX}/lib" isSystemDir) + if("${isSystemDir}" STREQUAL "-1") + set(CMAKE_INSTALL_RPATH "@loader_path/../lib") + endif("${isSystemDir}" STREQUAL "-1") +elseif(K4GAUDIPANDORA_SET_RPATH) + set(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/lib") # install LIBDIR + # the RPATH to be used when installing, but only if it's not a system directory + list(FIND CMAKE_PLATFORM_IMPLICIT_LINK_DIRECTORIES "${CMAKE_INSTALL_PREFIX}/lib" isSystemDir) + if("${isSystemDir}" STREQUAL "-1") + set(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/lib") + endif("${isSystemDir}" STREQUAL "-1") +else() + set(CMAKE_SKIP_INSTALL_RPATH TRUE) # skip the full RPATH for the install tree +endif() + + add_subdirectory(k4GaudiPandora) include(cmake/CreateProjectConfig.cmake) diff --git a/k4GaudiPandora/CMakeLists.txt b/k4GaudiPandora/CMakeLists.txt index c60834d..3d87864 100644 --- a/k4GaudiPandora/CMakeLists.txt +++ b/k4GaudiPandora/CMakeLists.txt @@ -18,21 +18,10 @@ limitations under the License. ]] set(_plugin_sources ${CMAKE_CURRENT_LIST_DIR}/ - # src/components/*.cpp src/CalorimeterHitType.cc - # src/DDBFieldPlugin.cc - # src/DDCaloDigi.cc - # src/DDCaloHitCreator.cc - # src/DDExternalClusteringAlgorithm.cc - # src/DDGeometryCreator.cc - # src/DDMCParticleCreator.cc - # src/DDPandoraPFANewProcessor.cc - # src/DDPfoCreator.cc - # src/DDScintillatorPpdDigi.cc + src/DDCaloDigi.cc + src/DDScintillatorPpdDigi.cc src/DDSimpleMuonDigi.cc - # src/DDTrackCreatorBase.cc - # src/DDTrackCreatorCLIC.cc - # src/DDTrackCreatorILD.cc ) gaudi_add_module(k4GaudiPandoraPlugins @@ -42,6 +31,7 @@ gaudi_add_module(k4GaudiPandoraPlugins EDM4HEP::edm4hep DD4hep::DDRec DD4hep::DDCore + CLHEP::CLHEP ) target_include_directories(k4GaudiPandoraPlugins PRIVATE include) install(TARGETS k4GaudiPandoraPlugins diff --git a/k4GaudiPandora/include/DDCaloDigi.h b/k4GaudiPandora/include/DDCaloDigi.h index 20a6173..2ced46e 100644 --- a/k4GaudiPandora/include/DDCaloDigi.h +++ b/k4GaudiPandora/include/DDCaloDigi.h @@ -16,22 +16,26 @@ * See the License for the specific language governing permissions and * limitations under the License. */ -#ifndef DDCCALODIGI_H -#define DDCCALODIGI_H 1 +#ifndef DDCALODIGI_H +#define DDCALODIGI_H 1 +#include "Gaudi/Accumulators/RootHistogram.h" #include "Gaudi/Property.h" -#include "GaudiAlg/Transformer.h" +#include "edm4hep/CaloHitSimCaloHitLinkCollection.h" #include "edm4hep/CalorimeterHitCollection.h" #include "edm4hep/EventHeaderCollection.h" -#include "edm4hep/MCRecoCaloAssociationCollection.h" #include "edm4hep/SimCalorimeterHitCollection.h" -#include "k4FWCore/BaseClass.h" +#include "k4FWCore/DataHandle.h" +#include "k4FWCore/MetaDataHandle.h" +#include "k4FWCore/Transformer.h" #include "k4Interface/IGeoSvc.h" #include "k4Interface/IUniqueIDGenSvc.h" -#include "CLHEP/Random/MTwistEngine.h" #include "CalorimeterHitType.h" #include "DDScintillatorPpdDigi.h" + +#include "CLHEP/Random/MTwistEngine.h" +#include "CLHEP/Random/RandGauss.h" #include "TFile.h" #include "TH1.h" #include "TH2.h" @@ -97,445 +101,341 @@ const int MAX_LAYERS = 200; const int MAX_STAVES = 16; -using retType = std::tuple, - edm4hep::MCRecoCaloAssociationCollection>; //FIXME: Does this need to be a map as well??? +using retType = std::tuple; + +using DDCaloDigi_t = k4FWCore::MultiTransformer; -struct DDCaloDigi final : k4FWCore ::MultiTransformer << retType > - (const std::map, const edm4hep::EventHeaderCollection&) > -{ +struct DDCaloDigi final : DDCaloDigi_t { DDCaloDigi(const std::string& name, ISvcLocator* svcLoc); StatusCode initialize() override; StatusCode finalize() override; - std::tuple operator()( - const edm4hep::SimCalorimeterHitCollection& simCaloHits, - const edm4hep::EventHeaderCollection& headers) const override; + retType operator()(const edm4hep::SimCalorimeterHitCollection& simCaloHits, + const edm4hep::EventHeaderCollection& headers) const; private: - Gaudi::Property _thresholdEcal{this, "ECALThreshold", {5.0e-5}, "Threshold for ECAL Hits in GeV"}; - Gaudi::Property _unitThresholdEcal{ + // check if input collection is ECAL or HCAL + Gaudi::Property m_inputColIsECAL{ + this, "InputColIsECAL", {true}, "Input SimCalorimeterHit collection is ECAL? If false, then is HCAL."}; + + // digitazing parameters for ECAL and HCAL + Gaudi::Property m_thresholdEcal{this, "ECALThreshold", {5.0e-5}, "Threshold for ECAL Hits in GeV"}; + Gaudi::Property m_unitThresholdEcal{ this, "ECALThresholdUnit", {"GeV"}, "Unit for ECAL Threshold. Can be \"GeV\", \"MIP\" or \"px\". MIP and px need properly set calibration constants"}; - Gaudi::Property _thresholdHcal{ - this, - "HCALThreshold", - {0.00004}, - "Unit for ECAL Threshold. Can be \"GeV\", \"MIP\" or \"px\". MIP and px need properly set calibration constants"}; - Gaudi::Property _unitThresholdHcal{ + Gaudi::Property> m_thresholdHcal{ + this, "HCALThreshold", {0.00004}, "Threshold for ECAL Hits in GeV"}; + Gaudi::Property m_unitThresholdHcal{ this, "HCALThresholdUnit", {"GeV"}, "Unit for HCAL Threshold. Can be \"GeV\", \"MIP\" or \"px\". MIP and px need properly set calibration constants"}; - std::vector ecalLayers; - ecalLayers.push_back(20); - ecalLayers.push_back(100); - Gaudi::Property> _ecalLayer{this, "ECALLayers", {ecalLayers}, "Index of ECal Layers"}; - std::vector hcalLayers; - hcalLayers.push_back(100); - Gaudi::Property> _hcalLayer{this, "ECALLayers", {hcalLayers}, "Index of HCal Layers"}; - std::vector calibrEcal; - calibrEcal.push_back(40.91); - calibrEcal.push_back(81.81); - Gaudi::Property> _calibrCoeffEcal{ - this, "CalibrECAL", {calibrEcal}, "Calibration coefficients for ECAL"}; - std::vector calibrHcalBarrel; - calibrHcalBarrel.push_back(0.); - Gaudi::Property> _calibrCoeffHcalBarrel{ - this, "CalibrHCALBarrel", {calibrHcalBarrel}, "Calibration coefficients for Barrel HCAL"}; - std::vector calibrHcalEndCap; - calibrHcalEndCap.push_back(0.); - Gaudi::Property> _calibrCoeffHcalEndCap{ - this, "CalibrHCALEndCap", {calibrHcalEndCap}, "Calibration coefficients for EndCap HCAL"}; - std::vector calibrHcalOther; - calibrHcalOther.push_back(0.); - Gaudi::Property> _calibrCoeffHcalOther{ - this, "CalibrHCALOther", {calibrHcalOther}, "Calibration coefficients for Other HCAL"}; - Gaudi::Property _digitalEcal{this, "IfDigitalEcal", {0}, "Digital Ecal"}; - Gaudi::Property _mapsEcalCorrection{ + Gaudi::Property> m_ecalLayers{this, "ECALLayers", {20, 100}, "Index of ECAL Layers"}; + Gaudi::Property> m_hcalLayers{this, "HCALLayers", {100}, "Index of HCAL Layers"}; + Gaudi::Property> m_calibrCoeffEcal{ + this, "CalibrECAL", {40.91, 81.81}, "Calibration coefficients for ECAL"}; + Gaudi::Property> m_calibrCoeffHcalBarrel{ + this, "CalibrHCALBarrel", {0.0}, "Calibration coefficients for Barrel HCAL"}; + Gaudi::Property> m_calibrCoeffHcalEndcap{ + this, "CalibrHCALEndcap", {0.0}, "Calibration coefficients for Endcap HCAL"}; + Gaudi::Property> m_calibrCoeffHcalOther{ + this, "CalibrHCALOther", {0.0}, "Calibration coefficients for Other HCAL"}; + Gaudi::Property m_digitalEcal{this, "IfDigitalEcal", {0}, "Digital ECAL"}; + Gaudi::Property m_mapsEcalCorrection{ this, "MapsEcalCorrection", {0}, "Ecal correction for theta dependency of calibration for MAPS"}; - Gaudi::Property _digitalHcal{this, "IfDigitalHcal", {0}, "Digital Hcal"}; - Gaudi::Property _ecalGapCorrection{ - this, "ECAconst float slop = 0.25; // (mm)LGapCorrection", {1}, "Correct for ECAL gaps"}; - Gaudi::Property _hcalGapCorrection{this, "HCALGapCorrection", {1}, "Correct for HCAL gaps"}; - Gaudi::Property _ecalEndcapCorrectionFactor{ + Gaudi::Property m_digitalHcal{this, "IfDigitalHcal", {0}, "Digital Hcal"}; + Gaudi::Property m_ecalGapCorrection{this, "ECALGapCorrection", {1}, "Correct for ECAL gaps"}; + Gaudi::Property m_hcalGapCorrection{this, "HCALGapCorrection", {1}, "Correct for HCAL gaps"}; + Gaudi::Property m_ecalEndcapCorrectionFactor{ this, "ECALEndcapCorrectionFactor", {1.025}, "Energy correction for ECAL endcap"}; - Gaudi::Property _hcalEndcapCorrectionFactor{ + Gaudi::Property m_hcalEndcapCorrectionFactor{ this, "HCALEndcapCorrectionFactor", {1.025}, "Energy correction for HCAL endcap"}; - Gaudi::Property _ecalGapCorrectionFactor{ + Gaudi::Property m_ecalGapCorrectionFactor{ this, "ECALGapCorrectionFactor", {1.0}, "Factor applied to gap correction"}; - Gaudi::Property _ecalModuleGapCorrectionFactor{ + Gaudi::Property m_ecalModuleGapCorrectionFactor{ this, "ECALModuleGapCorrectionFactor", {0.5}, "Factor applied to module gap correction ECAL"}; - Gaudi::Property _hcalModuleGapCorrectionFactor{ + Gaudi::Property m_hcalModuleGapCorrectionFactor{ this, "HCALModuleGapCorrectionFactor", {0.5}, "Factor applied to module gap correction HCAL"}; - Gaudi::Property _histograms{this, "Histograms", {0}, "Hit times histograms"}; - Gaudi::Property _useEcalTiming{this, "UseEcalTiming", {0}, "Use ECAL hit times"}; - Gaudi::Property _ecalCorrectTimesForPropagation{ + + // histograms that will be filled --- <1> or <2> is the number of dimensions of the histogram (1D or 2D) + mutable Gaudi::Accumulators::WeightedHistogram<1> fEcal{this, "fEcal", "Ecal time profile", {1000, 0., 1000.0}}; + mutable Gaudi::Accumulators::WeightedHistogram<1> fHcal{this, "fHcal", "Hcal time profile", {1000, 0., 1000.0}}; + + mutable Gaudi::Accumulators::WeightedHistogram<1> fEcalC{this, "fEcalC", "Ecal time profile cor", {1000, 0., 1000.0}}; + mutable Gaudi::Accumulators::WeightedHistogram<1> fHcalC{this, "fHcalC", "Hcal time profile cor", {1000, 0., 1000.0}}; + + mutable Gaudi::Accumulators::WeightedHistogram<1> fEcalC1{ + this, "fEcalC1", "Ecal time profile cor", {100, 0., 1000.0}}; + mutable Gaudi::Accumulators::WeightedHistogram<1> fHcalC1{ + this, "fHcalC1", "Hcal time profile cor", {100, 0., 1000.0}}; + + mutable Gaudi::Accumulators::WeightedHistogram<1> fEcalC2{this, "fEcalC2", "Ecal time profile cor", {10, 0., 1000.0}}; + mutable Gaudi::Accumulators::WeightedHistogram<1> fHcalC2{this, "fHcalC2", "Hcal time profile cor", {10, 0., 1000.0}}; + + mutable Gaudi::Accumulators::RootHistogram<2> fHcalCvsE{ + this, "fHcalCvsE", "Hcal time profile cor", {{100, 0., 500.0}, {100, 0., 10.}}}; + + mutable Gaudi::Accumulators::RootHistogram<2> fHcalLayer1{ + this, "fHcalLayer1", "Hcal layer 1 map", {{300, -4500., 4500.0}, {300, -4500, 4500.}}}; + mutable Gaudi::Accumulators::RootHistogram<2> fHcalLayer11{ + this, "fHcalLayer11", "Hcal layer 11 map", {{300, -4500., 4500.0}, {300, -4500, 4500.}}}; + mutable Gaudi::Accumulators::RootHistogram<2> fHcalLayer21{ + this, "fHcalLayer21", "Hcal layer 21 map", {{300, -4500., 4500.0}, {300, -4500, 4500.}}}; + mutable Gaudi::Accumulators::RootHistogram<2> fHcalLayer31{ + this, "fHcalLayer31", "Hcal layer 31 map", {{300, -4500., 4500.0}, {300, -4500, 4500.}}}; + mutable Gaudi::Accumulators::RootHistogram<2> fHcalLayer41{ + this, "fHcalLayer41", "Hcal layer 41 map", {{300, -4500., 4500.0}, {300, -4500, 4500.}}}; + mutable Gaudi::Accumulators::RootHistogram<2> fHcalLayer51{ + this, "fHcalLayer51", "Hcal layer 51 map", {{300, -4500., 4500.0}, {300, -4500, 4500.}}}; + mutable Gaudi::Accumulators::RootHistogram<2> fHcalLayer61{ + this, "fHcalLayer61", "Hcal layer 61 map", {{300, -4500., 4500.0}, {300, -4500, 4500.}}}; + mutable Gaudi::Accumulators::RootHistogram<2> fHcalLayer71{ + this, "fHcalLayer71", "Hcal layer 71 map", {{300, -4500., 4500.0}, {300, -4500, 4500.}}}; + + mutable Gaudi::Accumulators::RootHistogram<1> fHcalRLayer1{this, "fHcalRLayer1", "Hcal R layer 1", {50, 0., 5000.0}}; + mutable Gaudi::Accumulators::RootHistogram<1> fHcalRLayer11{ + this, "fHcalRLayer11", "Hcal R layer 11", {50, 0., 5000.0}}; + mutable Gaudi::Accumulators::RootHistogram<1> fHcalRLayer21{ + this, "fHcalRLayer21", "Hcal R layer 21", {50, 0., 5000.0}}; + mutable Gaudi::Accumulators::RootHistogram<1> fHcalRLayer31{ + this, "fHcalRLayer31", "Hcal R layer 31", {50, 0., 5000.0}}; + mutable Gaudi::Accumulators::RootHistogram<1> fHcalRLayer41{ + this, "fHcalRLayer41", "Hcal R layer 41", {50, 0., 5000.0}}; + mutable Gaudi::Accumulators::RootHistogram<1> fHcalRLayer51{ + this, "fHcalRLayer51", "Hcal R layer 51", {50, 0., 5000.0}}; + mutable Gaudi::Accumulators::RootHistogram<1> fHcalRLayer61{ + this, "fHcalRLayer61", "Hcal R layer 61", {50, 0., 5000.0}}; + mutable Gaudi::Accumulators::RootHistogram<1> fHcalRLayer71{ + this, "fHcalRLayer71", "Hcal R layer 71", {50, 0., 5000.0}}; + mutable Gaudi::Accumulators::WeightedHistogram<1> fHcalRLayerNorm{ + this, "fHcalRLayerNorm", "Hcal R layer Norm", {50, 0., 5000.0}}; + + mutable Gaudi::Accumulators::RootHistogram<2> fEcalLayer1{ + this, "fEcalLayer1", "Ecal layer 1 map", {{1800, -4500., 4500.0}, {1800, -4500, 4500.}}}; + mutable Gaudi::Accumulators::RootHistogram<2> fEcalLayer11{ + this, "fEcalLayer11", "Ecal layer 11 map", {{1800, -4500., 4500.0}, {1800, -4500, 4500.}}}; + mutable Gaudi::Accumulators::RootHistogram<2> fEcalLayer21{ + this, "fEcalLayer21", "Ecal layer 21 map", {{1800, -4500., 4500.0}, {1800, -4500, 4500.}}}; + + mutable Gaudi::Accumulators::RootHistogram<1> fEcalRLayer1{this, "fEcalRLayer1", "Ecal R layer 1", {100, 0., 5000.0}}; + mutable Gaudi::Accumulators::RootHistogram<1> fEcalRLayer11{ + this, "fEcalRLayer11", "Ecal R layer 11", {100, 0., 5000.0}}; + mutable Gaudi::Accumulators::RootHistogram<1> fEcalRLayer21{ + this, "fEcalRLayer21", "Ecal R layer 21", {100, 0., 5000.0}}; + mutable Gaudi::Accumulators::WeightedHistogram<1> fEcalRLayerNorm{ + this, "fEcalRLayerNorm", "Ecal R layer Norm", {100, 0., 5000.0}}; + + // timing parameters for ECAL + Gaudi::Property m_useEcalTiming{this, "UseEcalTiming", {0}, "Use ECAL hit times"}; + Gaudi::Property m_ecalCorrectTimesForPropagation{ this, "ECALCorrectTimesForPropagation", {0}, "Correct ECAL hit times for propagation: radial distance/c"}; - Gaudi::Property _ecalTimeWindowMin{this, "ECALTimeWindowMin", {-10.0}, "ECAL Time Window minimum time in ns"}; - Gaudi::Property _ecalEndcapTimeWindowMax{ - this, "ECALEndcapTimeWindowMax", {100}, "ECAL Endcap Time Window maximum time in ns"}; - Gaudi::Property _ecalBarrelTimeWindowMax{ - this, "ECALBarrelTimeWindowMax", {100}, "ECAL Barrel Time Window maximum time in ns"}; - Gaudi::Property _ecalDeltaTimeHitResolution{ - this, "ECALDeltaTimeHitResolution", {10}, "ECAL Minimum Delta Time in ns for resolving two hits"}; - Gaudi::Property _ecalTimeResolution{ - this, "ECALTimeResolution", {10}, "ECAL Time Resolution used to smear hit times"}; - Gaudi::Property _ecalSimpleTimingCut{ + Gaudi::Property m_ecalTimeWindowMin{this, "ECALTimeWindowMin", {-10.0}, "ECAL Time Window minimum time in ns"}; + Gaudi::Property m_ecalEndcapTimeWindowMax{ + this, "ECALEndcapTimeWindowMax", {100.0}, "ECAL Endcap Time Window maximum time in ns"}; + Gaudi::Property m_ecalBarrelTimeWindowMax{ + this, "ECALBarrelTimeWindowMax", {100.0}, "ECAL Barrel Time Window maximum time in ns"}; + Gaudi::Property m_ecalDeltaTimeHitResolution{ + this, "ECALDeltaTimeHitResolution", {10.0}, "ECAL Minimum Delta Time in ns for resolving two hits"}; + Gaudi::Property m_ecalTimeResolution{ + this, "ECALTimeResolution", {10.0}, "ECAL Time Resolution used to smear hit times"}; + Gaudi::Property m_ecalSimpleTimingCut{ this, "ECALSimpleTimingCut", {true}, "Use simple time window cut on hit times? If false: use original hit-time clustering algorithm. If true: use " "time window defined by ECALBarrelTimeWindowMin and ECALBarrelTimeWindowMax"}; - Gaudi::Property _useHcalTiming{this, "UseHcalTiming", {1}, "Use HCAL hit times"}; - Gaudi::Property _hcalCorrectTimesForPropagation{ + + // timing parameters for HCAL + Gaudi::Property m_useHcalTiming{this, "UseHcalTiming", {1}, "Use HCAL hit times"}; + Gaudi::Property m_hcalCorrectTimesForPropagation{ this, "HCALCorrectTimesForPropagation", {0}, "Correct HCAL hit times for propagation: radial distance/c"}; - Gaudi::Property _hcalTimeWindowMin{this, "HCALTimeWindowMin", {-10.0}, "HCAL Time Window minimum time in ns"}; - Gaudi::Property _hcalEndcapTimeWindowMax{ - this, "HCALEndcapTimeWindowMax", {100}, "HCAL Endcap Time Window maximum time in ns"}; - Gaudi::Property _hcalBarrelTimeWindowMax{ - this, "HCALBarrelTimeWindowMax", {100}, "HCAL Barrel Time Window maximum time in ns"}; - Gaudi::Property _hcalDeltaTimeHitResolution{ - this, "HCALDeltaTimeHitResolution", {10}, "HCAL Minimum Delta Time in ns for resolving two hits"}; - Gaudi::Property _hcalTimeResolution{ - this, "HCALTimeResolution", {10}, "HCAL Time Resolution used to smear hit times"}; - Gaudi::Property _hcalSimpleTimingCut{ + Gaudi::Property m_hcalTimeWindowMin{this, "HCALTimeWindowMin", {-10.0}, "HCAL Time Window minimum time in ns"}; + Gaudi::Property m_hcalEndcapTimeWindowMax{ + this, "HCALEndcapTimeWindowMax", {100.0}, "HCAL Endcap Time Window maximum time in ns"}; + Gaudi::Property m_hcalBarrelTimeWindowMax{ + this, "HCALBarrelTimeWindowMax", {100.0}, "HCAL Barrel Time Window maximum time in ns"}; + Gaudi::Property m_hcalDeltaTimeHitResolution{ + this, "HCALDeltaTimeHitResolution", {10.0}, "HCAL Minimum Delta Time in ns for resolving two hits"}; + Gaudi::Property m_hcalTimeResolution{ + this, "HCALTimeResolution", {10.0}, "HCAL Time Resolution used to smear hit times"}; + Gaudi::Property m_hcalSimpleTimingCut{ this, "HCALSimpleTimingCut", {true}, "Use simple time window cut on hit times? If false: use original hit-time clustering algorithm. If true: use " "time window defined by HCALBarrelTimeWindowMin and HCALBarrelTimeWindowMax"}; - Gaudi::Property _calibEcalMip{ - this, "CalibECALMIP", {1.0e-4}, "calibration to convert ECAL deposited energy to MIPs"}; - Gaudi::Property _applyEcalDigi{this, - "ECAL_apply_realistic_digi", - {0}, - "apply realistic digitisation to ECAL hits? (0=none, 1=silicon, 2=scintillator)"}; - Gaudi::Property _ecal_PPD_pe_per_mip{ - this, "ECAL_PPD_PE_per_MIP", {7.0}, "# Photo-electrons per MIP (scintillator): used to poisson smear #PEs if >0"}; - Gaudi::Property _ecal_PPD_n_pixels{ + + // parameters for extra ECAL digitization effects + Gaudi::Property m_calibEcalMip{ + this, "CalibECALMIP", {1.0e-4}, "Calibration to convert ECAL deposited energy to MIPs"}; + Gaudi::Property m_applyEcalDigi{ + this, + "ECALApplyRealisticDigi", + {0}, + "Apply realistic digitisation to ECAL hits? (0=none, 1=silicon, 2=scintillator)"}; + Gaudi::Property m_ecal_PPD_pe_per_mip{ + this, "ECAL_PPD_PE_per_MIP", {7.0}, "# photoelectrons per MIP (scintillator): used to poisson smear #PEs if > 0"}; + Gaudi::Property m_ecal_PPD_n_pixels{ this, "ECAL_PPD_N_Pixels", {10000}, "ECAL total number of MPPC/SiPM pixels for implementation of saturation effect"}; - Gaudi::Property _ecal_misCalibNpix{ + Gaudi::Property m_ecal_misCalibNpix{ this, "ECAL_PPD_N_Pixels_uncertainty", {0.05}, "ECAL fractional uncertainty of effective total number of MPPC/SiPM pixels"}; - Gaudi::Property _misCalibEcal_uncorrel{ + Gaudi::Property m_misCalibEcal_uncorrel{ this, "ECAL_miscalibration_uncorrel", {0.0}, - "uncorrelated ECAL random gaussian miscalibration (as a fraction: 1.0 = 100%"}; - Gaudi::Property _misCalibEcal_uncorrel_keep{ + "Uncorrelated ECAL random gaussian miscalibration (as a fraction: 1.0 = 100%)"}; + Gaudi::Property m_misCalibEcal_uncorrel_keep{ this, "ECAL_miscalibration_uncorrel_memorise", {false}, - "store oncorrelated ECAL miscalbrations in memory? (WARNING: can take a lot of memory if used..."}; - Gaudi::Property _misCalibEcal_correl{ + "Store uncorrelated ECAL miscalbrations in memory? (WARNING: Can take a lot of memory if used...)"}; + Gaudi::Property m_misCalibEcal_correl{ this, "ECAL_miscalibration_correl", {0.0}, - "correlated ECAL random gaussian miscalibration (as a fraction: 1.0 = 100%"}; - Gaudi::Property _deadCellFractionEcal{ - this, "ECAL_deadCellRate", {0.0}, "ECAL random dead cell fraction (as a fraction: 0->1)"}; - Gaudi::Property _deadCellEcal_keep{ + "Correlated ECAL random gaussian miscalibration (as a fraction: 1.0 = 100%)"}; + Gaudi::Property m_deadCellFractionEcal{ + this, "ECAL_deadCellRate", {0.0}, "ECAL random dead cell fraction (as a fraction: [0;1])"}; + Gaudi::Property m_deadCellEcal_keep{ this, "ECAL_deadCell_memorise", {false}, - "store dead ECAL cells in memory? (WARNING: can take a lot of memory if used..."}; - Gaudi::Property _strip_abs_length{ - this, "ECAL_strip_absorbtionLength", {1000000.0}, "length scale for absorbtion along scintillator strip (mm)"}; - Gaudi::Property _ecal_pixSpread{ - this, "ECAL_pixel_spread", {0.05}, "variation of mppc/sipm pixels capacitance in ECAL (as a fraction: 0.01=1%)"}; - Gaudi::Property _ecal_elec_noise{ - this, "ECAL_elec_noise_mips", {0.}, "typical electronics noise (ECAL, in MIP units)"}; - Gaudi::Property _ehEnergy{ - this, "energyPerEHpair", {3.6}, "energy required to create e-h pair in silicon (in eV)"}; - Gaudi::Property _ecalMaxDynMip{ - this, "ECAL_maxDynamicRange_MIP", {2500.}, "maximum of dynamic range for ECAL (in MIPs)"}; - Gaudi::Property _ecalStrip_default_nVirt{ - this, "StripEcal_default_nVirtualCells", {9}, "default number of virtual cells (used if not found in gear file)"}; - Gaudi::Property _ecal_deafult_layer_config{ + "Store dead ECAL cells in memory? (WARNING: Can take a lot of memory if used...)"}; + Gaudi::Property m_strip_abs_length{ + this, "ECAL_strip_absorbtionLength", {1000000.0}, "Length scale for absorbtion along scintillator strip (mm)"}; + Gaudi::Property m_ecal_pixSpread{ + this, "ECAL_pixel_spread", {0.05}, "Variation of MPPC/SiPM pixels capacitance in ECAL (as a fraction: 0.01=1%)"}; + Gaudi::Property m_ecal_elec_noise{ + this, "ECAL_elec_noise_mips", {0.0}, "Typical electronics noise for ECAL (in MIP units)"}; + Gaudi::Property m_ehEnergy{ + this, "energyPerEHpair", {3.6}, "Energy required to create e-h pair in silicon (in eV)"}; + Gaudi::Property m_ecalMaxDynMip{ + this, "ECAL_maxDynamicRange_MIP", {2500.0}, "Maximum of electronis dynamic range for ECAL (in MIPs)"}; + Gaudi::Property m_ecalStrip_default_nVirt{ + this, "StripEcal_default_nVirtualCells", {9}, "Default number of virtual cells (used if not found in gear file)"}; + Gaudi::Property m_ecal_deafult_layer_config{ this, "ECAL_default_layerConfig", {"000000000000000"}, - "default ECAL layer configuration (used if not found in gear file"}; - Gaudi::Property _calibHcalMip{ - this, "CalibHCALMIP", {1.0e-4}, "calibration to convert HCAL deposited energy to MIPs"}; - Gaudi::Property _applyHcalDigi{this, - "HCAL_apply_realistic_digi", - {0}, - "apply realistic digitisation to HCAL hits? (0=none, 1=silicon, 2=scintillator)"}; - Gaudi::Property _hcal_PPD_pe_per_mip{ - this, "HCAL_PPD_PE_per_MIP", {7.0}, "# Photo-electrons per MIP (scintillator): used to poisson smear #PEs if >0"}; - Gaudi::Property _hcal_PPD_n_pixels{ + "Default ECAL layer configuration (used if not found in gear file"}; + + // parameters for extra HCAL digitization effects + Gaudi::Property m_calibHcalMip{ + this, "CalibHCALMIP", {1.0e-4}, "Calibration to convert HCAL deposited energy to MIPs"}; + Gaudi::Property m_applyHcalDigi{ + this, + "HCALApplyRealisticDigi", + {0}, + "Apply realistic digitisation to HCAL hits? (0=none, 1=silicon, 2=scintillator)"}; + Gaudi::Property m_hcal_PPD_pe_per_mip{ + this, + "HCAL_PPD_PE_per_MIP", + {7.0}, + "# photo-electrons per MIP (scintillator): used to poisson smear #PEs if > 0"}; + Gaudi::Property m_hcal_PPD_n_pixels{ this, "HCAL_PPD_N_Pixels", {10000}, "HCAL total number of MPPC/SiPM pixels for implementation of saturation effect"}; - Gaudi::Property _hcal_misCalibNpix{ + Gaudi::Property m_hcal_misCalibNpix{ this, "HCAL_PPD_N_Pixels_uncertainty", {0.05}, "HCAL fractional uncertainty of effective total number of MPPC/SiPM pixels"}; - Gaudi::Property _misCalibHcal_uncorrel{ + Gaudi::Property m_misCalibHcal_uncorrel{ this, "HCAL_miscalibration_uncorrel", {0.0}, - "uncorrelated HCAL random gaussian miscalibration (as a fraction: 1.0 = 100%"}; - Gaudi::Property _misCalibHcal_uncorrel_keep{ + "Uncorrelated HCAL random gaussian miscalibration (as a fraction: 1.0 = 100%)"}; + Gaudi::Property m_misCalibHcal_uncorrel_keep{ this, "HCAL_miscalibration_uncorrel_memorise", {false}, - "store oncorrelated HCAL miscalbrations in memory? (WARNING: can take a lot of memory if used..."}; - Gaudi::Property _misCalibHcal_correl{ + "Store oncorrelated HCAL miscalbrations in memory? (WARNING: Can take a lot of memory if used...)"}; + Gaudi::Property m_misCalibHcal_correl{ this, "HCAL_miscalibration_correl", {0.0}, - "correlated HCAL random gaussian miscalibration (as a fraction: 1.0 = 100%"}; - Gaudi::Property _deadCellFractionHcal{ - this, "HCAL_deadCellRate", {0.0}, "HCAL random dead cell fraction (as a fraction: 0->1)"}; - Gaudi::Property _deadCellHcal_keep{ + "Correlated HCAL random gaussian miscalibration (as a fraction: 1.0 = 100%)"}; + Gaudi::Property m_deadCellFractionHcal{ + this, "HCAL_deadCellRate", {0.0}, "HCAL random dead cell fraction (as a fraction: [0;1])"}; + Gaudi::Property m_deadCellHcal_keep{ this, "HCAL_deadCell_memorise", {false}, - "store dead HCAL cells in memory? (WARNING: can take a lot of memory if used..."}; - Gaudi::Property _hcal_pixSpread{ - this, "HCAL_pixel_spread", {0.05}, "variation of mppc/sipm pixels capacitance in HCAL (as a fraction: 0.01=1%)"}; - Gaudi::Property _hcal_elec_noise{ - this, "HCAL_elec_noise_mips", {0.}, "typical electronics noise (ECAL, in MIP units)"}; - Gaudi::Property _hcalMaxDynMip{ - this, "HCAL_maxDynamicRange_MIP", {2500.}, "maximum of dynamic range for HCAL (in MIPs)"}; - //Gaudi::Property _ecalStrip_default_nVirt{this, "StripEcal_default_nVirtualCells", {9}, "default number of virtual cells (used if not found in gear file)"}; - //Gaudi::Property _ecal_deafult_layer_config{this, "ECAL_default_layerConfig", {"000000000000000"}, "default ECAL layer configuration (used if not found in gear file"}; - - const float slop = 0.25; // (mm) - //DDCaloDigi() ; - // DDCaloDigi(const DDCaloDigi&) = delete; - // DDCaloDigi& operator=(const DDCaloDigi&) = delete; - - virtual void fillECALGaps(); - - float digitalHcalCalibCoeff(CHT::Layout, float energy); - - float analogueHcalCalibCoeff(CHT::Layout, int layer); - - float digitalEcalCalibCoeff(int layer); - - float analogueEcalCalibCoeff(int layer); - - float ecalEnergyDigi(float energy, int id); - float ahcalEnergyDigi(float energy, int id); - - float siliconDigi(float energy); - float scintillatorDigi(float energy, bool isEcal); - auto combineVirtualStripCells(auto col, bool isBarrel, int stripOrientation); - - int getNumberOfVirtualCells(); - std::vector>& getLayerConfig(); - void checkConsistency(std::string colName, int layer); - std::pair getLayerProperties(std::string colName, int layer); - int getStripOrientationFromColName(std::string colName); - - int _nRun = 0; - int _nEvt = 0; - - //LCFlagImpl _flag{}; - - std::string _outputRelCollection = ""; - - float _thresholdEcal = 5.0e-5; - std::string _unitThresholdEcal = "GeV"; - std::vector _thresholdHcal{}; - std::string _unitThresholdHcal = "GeV"; - - int _digitalEcal = 0; - int _mapsEcalCorrection = 0; - int _digitalHcal = 0; - - //bool _ECAL_stripHits; - - std::vector _calibrCoeffEcal{}; - std::vector _calibrCoeffHcalBarrel{}; - std::vector _calibrCoeffHcalEndCap{}; - std::vector _calibrCoeffHcalOther{}; - - std::vector _ecalLayers{}; - std::vector _hcalLayers{}; - - int _ecalGapCorrection = 1; - float _ecalGapCorrectionFactor = 1; - float _ecalModuleGapCorrectionFactor = 0.5; - float _ecalEndcapCorrectionFactor = 1.025; - float _hcalEndcapCorrectionFactor = 1.025; - int _hcalGapCorrection = 1; - float _hcalModuleGapCorrectionFactor = 0.5; - - std::vector _calHitsByStaveLayer[MAX_STAVES][MAX_LAYERS]; - std::vector _calHitsByStaveLayerModule[MAX_STAVES][MAX_LAYERS]; - - float _zOfEcalEndcap = 0.0; - float _barrelPixelSizeT[MAX_LAYERS]; - float _barrelPixelSizeZ[MAX_LAYERS]; - float _endcapPixelSizeX[MAX_LAYERS]; - float _endcapPixelSizeY[MAX_LAYERS]; - float _barrelStaveDir[MAX_STAVES][2]; - - int _histograms = 0; - - // timing - int _useEcalTiming = 0; - int _ecalCorrectTimesForPropagation = 0; - float _ecalTimeWindowMin = -10.0; - float _ecalBarrelTimeWindowMax = 100.0; - float _ecalEndcapTimeWindowMax = 100.0; - float _ecalDeltaTimeHitResolution = 10.0; - float _ecalTimeResolution = 10.0; - bool _ecalSimpleTimingCut = true; - - int _useHcalTiming = 1; - int _hcalCorrectTimesForPropagation = 0; - float _hcalTimeWindowMin = -10.0; - float _hcalBarrelTimeWindowMax = 100.0; - float _hcalEndcapTimeWindowMax = 100.0; - float _hcalDeltaTimeHitResolution = 10.0; - float _hcalTimeResolution = 10.0; - bool _hcalSimpleTimingCut = true; - - std::unique_ptr _scEcalDigi{}; - std::unique_ptr _scHcalDigi{}; + "Store dead HCAL cells in memory? (WARNING: Can take a lot of memory if used...)"}; + Gaudi::Property m_hcal_pixSpread{ + this, "HCAL_pixel_spread", {0.05}, "Variation of MPPC/SiPM pixels capacitance in HCAL (as a fraction: 0.01=1%)"}; + Gaudi::Property m_hcal_elec_noise{ + this, "HCAL_elec_noise_mips", {0.}, "Typical electronics noise for HCAL (in MIP units)"}; + Gaudi::Property m_hcalMaxDynMip{ + this, "HCAL_maxDynamicRange_MIP", {2500.}, "Maximum of electronis dynamic range for HCAL (in MIPs)"}; + //Gaudi::Property m_hcalStrip_default_nVirt{this, "StripHcal_default_nVirtualCells", {9}, "Default number of virtual cells (used if not found in gear file)"}; + //Gaudi::Property m_hcal_deafult_layer_config{this, "HCAL_default_layerConfig", {"000000000000000"}, "Default HCAL layer configuration (used if not found in gear file"}; + //Gaudi::Property m_encodingStringVariable{this, "EncodingStringParameterName", "GlobalTrackerReadoutID", + //"The name of the DD4hep constant that contains the Encoding string for tracking detectors"}; - // parameters for extra ECAL digitization effects - float _calibEcalMip = 1.0e-4; // MIP calibration factor - int _applyEcalDigi = 0; // which realistic calib to apply - float _ecal_PPD_pe_per_mip = 7; // # photoelectrons/MIP for MPPC - int _ecal_PPD_n_pixels = 10000; // # pixels in MPPC - float _ehEnergy = 3.6; // energy to create e-h pair in silicon - float _ecal_misCalibNpix = 0.05; // miscalibration of # MPPC pixels - - float _misCalibEcal_uncorrel = 0.0; // general ECAL miscalibration (uncorrelated between channels) - bool _misCalibEcal_uncorrel_keep = - false; // if true, use the same ECAL cell miscalibs in each event (requires more memory) - float _misCalibEcal_correl = 0.0; // general ECAL miscalibration (100% uncorrelated between channels) - - float _deadCellFractionEcal = 0.0; // fraction of random dead channels - bool _deadCellEcal_keep = false; // keep same cells dead between events? - - float _strip_abs_length = 1000000; // absorption length along strip for non-uniformity modeling - float _ecal_pixSpread = 0.05; // relative spread of MPPC pixel signal - float _ecal_elec_noise = 0; // electronics noise (as fraction of MIP) - float _ecalMaxDynMip = 2500; // electronics dynamic range (in terms of MIPs) - int _ecalStrip_default_nVirt = - 9; // # virtual cells used in Mokka simulation of strips (if available, this is taken from gear file) - std::string _ecal_deafult_layer_config = - "000000000000000"; // ECAL layer configuration (if available, this is taken from gear file) - - // parameters for extra AHCAL digitization effects - float _calibHcalMip = 1.0e-4; // MIP calibration factor - int _applyHcalDigi = 0; // which realistic calib to apply - float _hcal_PPD_pe_per_mip = 10; // # photoelectrons/MIP for MPPC - int _hcal_PPD_n_pixels = 400; // # pixels in MPPC - float _hcal_misCalibNpix = 0.05; // miscalibration of # MPPC pixels - - float _misCalibHcal_uncorrel = 0.0; // general ECAL miscalibration (uncorrelated between channels) - bool _misCalibHcal_uncorrel_keep = - false; // if true, use the same AHCAL cell miscalibs in each event (requires more memory) - float _misCalibHcal_correl = 0.0; // general ECAL miscalibration (100% uncorrelated between channels) - - float _deadCellFractionHcal = 0.0; // fraction of random dead channels - bool _deadCellHcal_keep = false; // keep same cells dead between events? - float _hcal_pixSpread = 0.0; // relative spread of MPPC pixel signal - float _hcal_elec_noise = 0.0; // electronics noise (as fraction of MIP) - float _hcalMaxDynMip = 200; // electronics dynamic range (in terms of MIPs) + SmartIF m_geoSvc; + SmartIF m_uidSvc; - // internal variables - std::vector> _layerTypes{}; - int _strip_virt_cells = 999; - int _countWarnings = 0; - std::string _ecalLayout = ""; + //const float slop = 0.25; // (mm) - float _event_correl_miscalib_ecal = 0.0; - float _event_correl_miscalib_hcal = 0.0; + virtual void fillECALGaps(std::vector m_calHitsByStaveLayer[MAX_STAVES][MAX_LAYERS], + std::vector m_calHitsByStaveLayerModule[MAX_STAVES][MAX_LAYERS]) const; - CLHEP::MTwistEngine* _randomEngineDeadCellEcal = NULL; - CLHEP::MTwistEngine* _randomEngineDeadCellHcal = NULL; + float digitalEcalCalibCoeff(int layer) const; + float analogueEcalCalibCoeff(int layer) const; - std::map, float> _ECAL_cell_miscalibs{}; - std::map, bool> _ECAL_cell_dead{}; - std::map, float> _HCAL_cell_miscalibs{}; - std::map, bool> _HCAL_cell_dead{}; + float digitalHcalCalibCoeff(CHT::Layout, float energy) const; + float analogueHcalCalibCoeff(CHT::Layout, int layer) const; - enum { SQUARE, STRIP_ALIGN_ALONG_SLAB, STRIP_ALIGN_ACROSS_SLAB, SIECAL = 0, SCECAL }; + float ecalEnergyDigi(float energy, int id) const; + float hcalEnergyDigi(float energy, int id) const; + + float siliconDigi(float energy) const; + float scintillatorDigi(float energy, bool isEcal) const; + + std::vector> getLayerConfig() const; + void checkConsistency(std::string colName, int layer) const; + std::pair getLayerProperties(std::string const& colName, int layer) const; + int getStripOrientationFromColName(std::string const& colName) const; + + float m_zOfEcalEndcap = 0.0; + float m_barrelPixelSizeT[MAX_LAYERS]; + float m_barrelPixelSizeZ[MAX_LAYERS]; + float m_endcapPixelSizeX[MAX_LAYERS]; + float m_endcapPixelSizeY[MAX_LAYERS]; + float m_barrelStaveDir[MAX_STAVES][2]; - TH1F* fEcal = NULL; - TH1F* fHcal = NULL; - TH1F* fEcalC = NULL; - TH1F* fHcalC = NULL; - TH1F* fEcalC1 = NULL; - TH1F* fHcalC1 = NULL; - TH1F* fEcalC2 = NULL; - TH1F* fHcalC2 = NULL; - TH2F* fHcalCvsE = NULL; - TH2F* fHcalLayer1 = NULL; - TH2F* fHcalLayer11 = NULL; - TH2F* fHcalLayer21 = NULL; - TH2F* fHcalLayer31 = NULL; - TH2F* fHcalLayer41 = NULL; - TH2F* fHcalLayer51 = NULL; - TH2F* fHcalLayer61 = NULL; - TH2F* fHcalLayer71 = NULL; - TH1F* fHcalRLayer1 = NULL; - TH1F* fHcalRLayer11 = NULL; - TH1F* fHcalRLayer21 = NULL; - TH1F* fHcalRLayer31 = NULL; - TH1F* fHcalRLayer41 = NULL; - TH1F* fHcalRLayer51 = NULL; - TH1F* fHcalRLayer61 = NULL; - TH1F* fHcalRLayer71 = NULL; - TH1F* fHcalRLayerNorm = NULL; - - TH1F* fEcalRLayerNorm = NULL; - TH2F* fEcalLayer1 = NULL; - TH2F* fEcalLayer11 = NULL; - TH2F* fEcalLayer21 = NULL; - TH1F* fEcalRLayer1 = NULL; - TH1F* fEcalRLayer11 = NULL; - TH1F* fEcalRLayer21 = NULL; + std::unique_ptr m_scEcalDigi{}; + std::unique_ptr m_scHcalDigi{}; + + // internal variables + mutable int m_countWarnings = 0; + std::string m_ecalLayout = ""; + + CLHEP::MTwistEngine* m_randomEngineDeadCellEcal = NULL; + CLHEP::MTwistEngine* m_randomEngineDeadCellHcal = NULL; + float m_event_correl_miscalib_ecal = CLHEP::RandGauss::shoot(1.0, m_misCalibEcal_correl); + float m_event_correl_miscalib_hcal = CLHEP::RandGauss::shoot(1.0, m_misCalibHcal_correl); + + std::map m_ECAL_cell_miscalibs{}; + std::map m_ECAL_cell_dead{}; + std::map m_HCAL_cell_miscalibs{}; + std::map m_HCAL_cell_dead{}; + + enum { SQUARE, STRIP_ALIGN_ALONG_SLAB, STRIP_ALIGN_ACROSS_SLAB, SIECAL = 0, SCECAL }; }; + DECLARE_COMPONENT(DDCaloDigi) -#endif - -// ecalCollections.push_back(std::string("EcalRingCollection")); -// std::vector _ecalCollections; // this is for silicon -// _ecalCollections.push_back(std::string("EcalBarrelCollection")); -// _ecalCollections.push_back(std::string("EcalEndcapCollection")); -//Gaudi::Property> _ecalCollections{this, "ECALCollections", {_ecalCollections}, "ECAL Collection Names"}; - -// std::vector _hcalCollections; -// _hcalCollections.push_back(std::string("HcalBarrelRegCollection")); -// _hcalCollections.push_back(std::string("HcalEndcapRingsCollection")); -// _hcalCollections.push_back(std::string("HcalEndcapsCollection")); -//Gaudi::Property> _hcalCollections{this, "HCALCollections", {_hcalCollections}, "HCAL Collection Names"}; - -// std::vector _outputEcalCollections{}; -// _outputEcalCollections.push_back(std::string("ECALBarrel")); -// _outputEcalCollections.push_back(std::string("ECALEndcap")); -// _outputEcalCollections.push_back(std::string("ECALOther")); - -// std::vector _outputHcalCollections{}; -// _outputHcalCollections.push_back(std::string("HCALBarrel")); -// _outputHcalCollections.push_back(std::string("HCALEndcap")); -// _outputHcalCollections.push_back(std::string("HCALOther")); - -// Gaudi::Property _outputEcalCollections[0]{this, "output_ECALCollections0", {"ECALBarrel"}, "ECAL Collection of real Hits in Barrel"}; -// Gaudi::Property _outputEcalCollections[1]{this, "output_ECALCollections1", {"ECALEndcap"}, "ECAL Collection of real Hits in EndCap"}; -// Gaudi::Property _outputEcalCollections[2]{this, "output_ECALCollections2", {"ECALOther"}, "ECAL Collection of real Hits other"}; - -// Gaudi::Property _outputHcalCollections[0]{this, "output_HCALCollections0", {"HCALBarrel"}, "HCAL Collection of real Hits in Barrel"}; -// Gaudi::Property _outputHcalCollections[1]{this, "output_HCALCollections1", {"HCALEndcap"}, "HCAL Collection of real Hits in EndCap"}; -// Gaudi::Property _outputHcalCollections[2]{this, "output_HCALCollections2", {"HCALOther"}, "HCAL Collection of real Hits other"}; -// Gaudi::Property _outputRelCollection{this, "RelationOutputCollection", {"RelationCaloHit"}, "CaloHit Relation Collection"}; +#endif \ No newline at end of file diff --git a/k4GaudiPandora/include/DDScintillatorPpdDigi.h b/k4GaudiPandora/include/DDScintillatorPpdDigi.h index 0a9e851..d5fd87b 100644 --- a/k4GaudiPandora/include/DDScintillatorPpdDigi.h +++ b/k4GaudiPandora/include/DDScintillatorPpdDigi.h @@ -16,8 +16,8 @@ * See the License for the specific language governing permissions and * limitations under the License. */ -#ifndef _DDScintillatorPpdDigi_h_ -#define _DDScintillatorPpdDigi_h_ +#ifndef DDSCINTILLATORPPDDIGI_H +#define DDSCINTILLATORPPDDIGI_H class DDScintillatorPpdDigi { public: @@ -25,38 +25,38 @@ class DDScintillatorPpdDigi { ~DDScintillatorPpdDigi() {} // expected # photoelectrons / MIP - void setPEperMIP(float x) { _pe_per_mip = x; } + void setPEperMIP(float x) { m_PEperMIP = x; } // calibration factor from input hit energy to MIPs - void setCalibMIP(float x) { _calib_mip = x; } + void setCalibMIP(float x) { m_calibMIP = x; } - // #pixels of PPD - void setNPix(int x) { _npix = x; } + // # pixels of PPD + void setNPix(int x) { m_Npix = x; } - // random miscalibration of total #pixels (as a fraction of pixel number: 0.05 = 5% miscalibration) - void setRandomMisCalibNPix(float x) { _misCalibNpix = x; } + // random miscalibration of total # pixels (as a fraction of pixel number: 0.05 = 5% miscalibration) + void setRandomMisCalibNPix(float x) { m_misCalibNpix = x; } // spread in pixel capacitance (as a fraction: 0.05 = 5% spread) - void setPixSpread(float x) { _pixSpread = x; } + void setPixSpread(float x) { m_pixSpread = x; } // electronics noise (in MIP units) - void setElecNoise(float x) { _elecNoise = x; } + void setElecNoise(float x) { m_elecNoise = x; } // electronics dynamic range (in MIP units) - void setElecRange(float x) { _elecMaxDynRange_MIP = x; } + void setElecRange(float x) { m_elecMaxDynRange_MIP = x; } float getDigitisedEnergy(float energy); void printParameters(); private: - float _pe_per_mip = -99; - float _calib_mip = -99; - float _npix = -99; - float _misCalibNpix = 0; - float _pixSpread = 0; - float _elecNoise = 0; - float _elecMaxDynRange_MIP = 0; + float m_PEperMIP = -99; + float m_calibMIP = -99; + float m_Npix = -99; + float m_misCalibNpix = 0; + float m_pixSpread = 0; + float m_elecNoise = 0; + float m_elecMaxDynRange_MIP = 0; }; #endif diff --git a/k4GaudiPandora/options/runDDCaloDigi.py b/k4GaudiPandora/options/runDDCaloDigi.py new file mode 100644 index 0000000..d30b447 --- /dev/null +++ b/k4GaudiPandora/options/runDDCaloDigi.py @@ -0,0 +1,167 @@ +# +# Copyright (c) 2020-2024 Key4hep-Project. +# +# This file is part of Key4hep. +# See https://key4hep.github.io/key4hep-doc/ for further info. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# +from Gaudi.Configuration import INFO +from k4FWCore import ApplicationMgr, IOSvc +from Configurables import EventDataSvc +from Configurables import DDCaloDigi + +from Configurables import GeoSvc +from Configurables import UniqueIDGenSvc +from Configurables import RootHistSvc +from Configurables import Gaudi__Histograming__Sink__Root as RootHistoSink +import os + +id_service = UniqueIDGenSvc("UniqueIDGenSvc") + +geoservice = GeoSvc("GeoSvc") +geoservice.detectors = [os.environ["K4GEO"]+"/FCCee/CLD/compact/CLD_o2_v06/CLD_o2_v06.xml"] +geoservice.OutputLevel = INFO +geoservice.EnableGeant4Geo = False + +calodigi = [DDCaloDigi("ECALBarrelDigi"), + DDCaloDigi("ECALEndcapDigi"), + DDCaloDigi("HCALBarrelDigi"), + DDCaloDigi("HCALEndcapDigi"), + DDCaloDigi("HCALRingDigi")] + +ECALorHCAL = [True, True, False, False, False] + +inputcollections = [["ECalBarrelCollection"], + ["ECalEndcapCollection"], + ["HCalBarrelCollection"], + ["HCalEndcapCollection"], + ["HCalRingCollection"]] + +outputcollections = [["ECALBarrel"], + ["ECALEndcap"], + ["HCALBarrel"], + ["HCALEndcap"], + ["HCALRing"]] + +relcollections = [["RelationCaloHitECALBarrel"], + ["RelationCaloHitECALEndcap"], + ["RelationCaloHitHCALBarrel"], + ["RelationCaloHitHCALEndcap"], + ["RelationCaloHitHCALRing"]] + +#set properties +for calodigicol, ecalorhcal, inputcol, outputcol, relcol in zip(calodigi, ECALorHCAL, inputcollections, outputcollections, relcollections): + + calodigicol.InputColIsECAL = ecalorhcal # True -- ECAL // False -- HCAL + calodigicol.InputCaloHitCollection = inputcol # "ECalBarrelCollection","ECalEndcapCollection" + # "HCalBarrelCollection","HCalEndcapCollection","HCalRingCollection" + calodigicol.OutputCaloHitCollection = outputcol + calodigicol.RelCollection = relcol + + # digitazing parameters for ECAL and HCAL + calodigicol.ECALThreshold = 5.0e-5 + calodigicol.ECALThresholdUnit = "GeV" + calodigicol.HCALThreshold = [0.00025] + calodigicol.HCALThresholdUnit = "GeV" + calodigicol.ECALLayers = [41,100] + calodigicol.HCALLayers = [100] + calodigicol.CalibrECAL = [37.5227197175, 37.5227197175] + calodigicol.CalibrHCALBarrel = [45.9956826061] + calodigicol.CalibrHCALEndcap = [46.9252540291] + calodigicol.CalibrHCALOther = [57.4588011802] + calodigicol.IfDigitalEcal = 0 + calodigicol.MapsEcalCorrection = 0 + calodigicol.IfDigitalHcal = 0 + calodigicol.ECALGapCorrection = 1 + calodigicol.HCALGapCorrection = 1 + calodigicol.ECALEndcapCorrectionFactor = 1.03245503522 + calodigicol.HCALEndcapCorrectionFactor = 1.000 + calodigicol.ECALGapCorrectionFactor = 1.0 + calodigicol.ECALModuleGapCorrectionFactor = 0.0 + calodigicol.HCALModuleGapCorrectionFactor = 0.5 + + # timing parameters for ECAL + calodigicol.UseEcalTiming = 1 + calodigicol.ECALCorrectTimesForPropagation = 1 + calodigicol.ECALTimeWindowMin = -1.0 + calodigicol.ECALEndcapTimeWindowMax = 10.0 + calodigicol.ECALBarrelTimeWindowMax = 10.0 + calodigicol.ECALDeltaTimeHitResolution = 10.0 + calodigicol.ECALTimeResolution = 10.0 + calodigicol.ECALSimpleTimingCut = True + + # timing parameters for HCAL + calodigicol.UseHcalTiming = 1 + calodigicol.HCALCorrectTimesForPropagation = 1 + calodigicol.HCALTimeWindowMin = -1.0 + calodigicol.HCALEndcapTimeWindowMax = 10.0 + calodigicol.HCALBarrelTimeWindowMax = 10.0 + calodigicol.HCALDeltaTimeHitResolution = 10.0 + calodigicol.HCALTimeResolution = 10.0 + calodigicol.HCALSimpleTimingCut = True + + # parameters for extra ECAL digitization effects + calodigicol.CalibECALMIP = 1.0e-4 + calodigicol.ECALApplyRealisticDigi = 0 + calodigicol.ECAL_PPD_PE_per_MIP = 7.0 + calodigicol.ECAL_PPD_N_Pixels = 10000 + calodigicol.ECAL_PPD_N_Pixels_uncertainty = 0.05 + calodigicol.ECAL_miscalibration_uncorrel = 0.0 + calodigicol.ECAL_miscalibration_uncorrel_memorise = False + calodigicol.ECAL_miscalibration_correl = 0.0 + calodigicol.ECAL_deadCellRate = 0.0 + calodigicol.ECAL_deadCell_memorise = False + calodigicol.ECAL_strip_absorbtionLength = 1.0e6 + calodigicol.ECAL_pixel_spread = 0.05 + calodigicol.ECAL_elec_noise_mips = 0.0 + calodigicol.energyPerEHpair = 3.6 + calodigicol.ECAL_maxDynamicRange_MIP = 2500.0 + calodigicol.StripEcal_default_nVirtualCells = 9 + calodigicol.ECAL_default_layerConfig = "000000000000000" + + # parameters for extra HCAL digitization effects + calodigicol.CalibHCALMIP = 1.0e-4 + calodigicol.HCALApplyRealisticDigi = 0 + calodigicol.HCAL_PPD_PE_per_MIP = 10.0 + calodigicol.HCAL_PPD_N_Pixels = 400 + calodigicol.HCAL_PPD_N_Pixels_uncertainty = 0.05 + calodigicol.HCAL_miscalibration_uncorrel = 0.0 + calodigicol.HCAL_miscalibration_uncorrel_memorise = False + calodigicol.HCAL_miscalibration_correl = 0.0 + calodigicol.HCAL_deadCellRate = 0.0 + calodigicol.HCAL_deadCell_memorise = False + calodigicol.HCAL_pixel_spread = 0.0 + calodigicol.HCAL_elec_noise_mips = 0.0 + calodigicol.HCAL_maxDynamicRange_MIP = 200.0 + + +# +iosvc = IOSvc() +iosvc.Input = "../simulation/sim_partgun_1000.root" +iosvc.Output = "../outputfiles/DDCaloDigi/outputCaloDigi_Gaudi.root" + +hps = RootHistSvc("HistogramPersistencySvc") +root_hist_svc = RootHistoSink("RootHistoSink") +root_hist_svc.FileName = "../outputfiles/DDCaloDigi/ddcalodigi_hist.root" + +ApplicationMgr(TopAlg=[calodigi[0], + calodigi[1], + calodigi[2], + calodigi[3], + calodigi[4]], + EvtSel="NONE", + EvtMax=-1, + ExtSvc=[EventDataSvc("EventDataSvc"), root_hist_svc], + OutputLevel=INFO, + ) diff --git a/k4GaudiPandora/src/DDCaloDigi.cc b/k4GaudiPandora/src/DDCaloDigi.cc index 36d8b8a..00387b5 100644 --- a/k4GaudiPandora/src/DDCaloDigi.cc +++ b/k4GaudiPandora/src/DDCaloDigi.cc @@ -25,106 +25,68 @@ #include #include #include "DD4hep/DD4hepUnits.h" +#include "DD4hep/DetType.h" #include "DD4hep/Detector.h" +#include "DD4hep/DetectorSelector.h" +#include "DD4hep/Factories.h" #include "DDRec/DetectorData.h" +#include "DDRec/MaterialManager.h" #include "GaudiKernel/MsgStream.h" #include "edm4hep/CalorimeterHit.h" #include "edm4hep/Constants.h" -#include "k4FWCore/BaseClass.h" #include "CLHEP/Random/RandFlat.h" #include "CLHEP/Random/RandGauss.h" #include "CLHEP/Random/RandPoisson.h" -#include "DD4hep/DD4hepUnits.h" -#include "DD4hep/DetType.h" -#include "DD4hep/Factories.h" -#include "DDRec/DetectorData.h" -#include "DDRec/MaterialManager.h" -using namespace std; using namespace dd4hep; using namespace DDSegmentation; -// protect against rounding errors -// will not find caps smaller than this -const float slop = 0.25; // (mm) -//const float pi = acos(-1.0); ///FIXME -//const float twopi = 2.0*pi; ///FIXME: DD4HEP INTERFERES WITH THESE +// Forward Declaration, gets linked in from DDPandoraPFANewProcessor --- FIXME +// now declare here, to be linked from DDPandoraPFANewProcessor +dd4hep::rec::LayeredCalorimeterData* getExtension(unsigned int includeFlag, unsigned int excludeFlag = 0) { + dd4hep::rec::LayeredCalorimeterData* theExtension = 0; + + dd4hep::Detector& mainDetector = dd4hep::Detector::getInstance(); + const std::vector& theDetectors = + dd4hep::DetectorSelector(mainDetector).detectors(includeFlag, excludeFlag); + + std::cout << " getExtension : includeFlag: " << dd4hep::DetType(includeFlag) + << " excludeFlag: " << dd4hep::DetType(excludeFlag) << " found : " << theDetectors.size() + << " - first det: " << theDetectors.at(0).name() << std::endl; + + if (theDetectors.size() != 1) { + std::stringstream es; + es << " getExtension: selection is not unique (or empty) includeFlag: " << dd4hep::DetType(includeFlag) + << " excludeFlag: " << dd4hep::DetType(excludeFlag) << " --- found detectors : "; + for (unsigned i = 0, N = theDetectors.size(); i < N; ++i) { + es << theDetectors.at(i).name() << ", "; + } + throw std::runtime_error(es.str()); + } + + theExtension = theDetectors.at(0).extension(); -DECLARE_COMPONENT(DDCaloDigi) + return theExtension; +}; DDCaloDigi::DDCaloDigi(const std::string& aName, ISvcLocator* aSvcLoc) - : MultiTransformer( - aName, aSvcLoc, - { - KeyValues("SimCaloHitCollections", {"SimCalorimeterHitCollection1", "SimCalorimeterHitCollection2", - "SimCalorimeterHitCollection3"}), - KeyValues("HeaderName", {"EventHeader"}), - }, - {KeyValues("output_ECALCollections", {"ECalorimeterHit1", "ECalorimeterHit2", "ECalorimeterHit3"}), - KeyValues("output_HCALCollections", {"HCalorimeterHit1", "HCalorimeterHit2", "HCalorimeterHit3"}), - KeyValues("RelationOutputCollection", {"RelationSimCaloHit"})}) { + : MultiTransformer(aName, aSvcLoc, + { + KeyValues("InputCaloHitCollection", {"ECalBarrelCollection"}), + KeyValues("HeaderName", {"EventHeader"}), + }, + {KeyValues("OutputCaloHitCollection", {"ECalorimeterHit1"}), + KeyValues("RelCollection", {"RelationCaloHit"})}) { m_uidSvc = service("UniqueIDGenSvc", true); if (!m_uidSvc) { error() << "Unable to get UniqueIDGenSvc" << endmsg; } + m_geoSvc = serviceLocator()->service("GeoSvc"); // important to initialize m_geoSvc } -//forward declarations. See in DDPandoraPFANewProcessor.cc -dd4hep::rec::LayeredCalorimeterData* getExtension(unsigned int includeFlag, unsigned int excludeFlag = 0); - -// // helper struct for string comparision -// struct XToLower{ -// int operator() ( int ch ) { -// return std::tolower ( ch ); -// } -// } - -if (_histograms) { - fEcal = new TH1F("fEcal", "Ecal time profile", 1000, 0., 1000.0); - fHcal = new TH1F("fHcal", "Hcal time profile", 1000, 0., 1000.0); - - fEcalC = new TH1F("fEcalC", "Ecal time profile cor", 1000, 0., 1000.0); - fHcalC = new TH1F("fHcalC", "Hcal time profile cor", 1000, 0., 1000.0); - - fEcalC1 = new TH1F("fEcalC1", "Ecal time profile cor", 100, 0., 1000.0); - fHcalC1 = new TH1F("fHcalC1", "Hcal time profile cor", 100, 0., 1000.0); - - fEcalC2 = new TH1F("fEcalC2", "Ecal time profile cor", 10, 0., 1000.0); - fHcalC2 = new TH1F("fHcalC2", "Hcal time profile cor", 10, 0., 1000.0); - - fHcalLayer1 = new TH2F("fHcalLayer1", "Hcal layer 1 map", 300, -4500., 4500.0, 300, -4500, 4500.); - fHcalLayer11 = new TH2F("fHcalLayer11", "Hcal layer 11 map", 300, -4500., 4500.0, 300, -4500, 4500.); - fHcalLayer21 = new TH2F("fHcalLayer21", "Hcal layer 21 map", 300, -4500., 4500.0, 300, -4500, 4500.); - fHcalLayer31 = new TH2F("fHcalLayer31", "Hcal layer 31 map", 300, -4500., 4500.0, 300, -4500, 4500.); - fHcalLayer41 = new TH2F("fHcalLayer41", "Hcal layer 41 map", 300, -4500., 4500.0, 300, -4500, 4500.); - fHcalLayer51 = new TH2F("fHcalLayer51", "Hcal layer 51 map", 300, -4500., 4500.0, 300, -4500, 4500.); - fHcalLayer61 = new TH2F("fHcalLayer61", "Hcal layer 61 map", 300, -4500., 4500.0, 300, -4500, 4500.); - fHcalLayer71 = new TH2F("fHcalLayer71", "Hcal layer 71 map", 300, -4500., 4500.0, 300, -4500, 4500.); - fHcalRLayer1 = new TH1F("fHcalRLayer1", "Hcal R layer 1", 50, 0., 5000.0); - fHcalRLayer11 = new TH1F("fHcalRLayer11", "Hcal R layer 11", 50, 0., 5000.0); - fHcalRLayer21 = new TH1F("fHcalRLayer21", "Hcal R layer 21", 50, 0., 5000.0); - fHcalRLayer31 = new TH1F("fHcalRLayer31", "Hcal R layer 31", 50, 0., 5000.0); - fHcalRLayer41 = new TH1F("fHcalRLayer41", "Hcal R layer 41", 50, 0., 5000.0); - fHcalRLayer51 = new TH1F("fHcalRLayer51", "Hcal R layer 51", 50, 0., 5000.0); - fHcalRLayer61 = new TH1F("fHcalRLayer61", "Hcal R layer 61", 50, 0., 5000.0); - fHcalRLayer71 = new TH1F("fHcalRLayer71", "Hcal R layer 71", 50, 0., 5000.0); - fHcalRLayerNorm = new TH1F("fHcalRLayerNorm", "Hcal R layer Norm", 50, 0., 5000.0); - - fEcalLayer1 = new TH2F("fEcalLayer1", "Ecal layer 1 map", 1800, -4500., 4500.0, 1800, -4500, 4500.); - fEcalLayer11 = new TH2F("fEcalLayer11", "Ecal layer 11 map", 1800, -4500., 4500.0, 1800, -4500, 4500.); - fEcalLayer21 = new TH2F("fEcalLayer21", "Ecal layer 21 map", 1800, -4500., 4500.0, 1800, -4500, 4500.); - fEcalRLayer1 = new TH1F("fEcalRLayer1", "Ecal R layer 1", 100, 0., 5000.0); - fEcalRLayer11 = new TH1F("fEcalRLayer11", "Ecal R layer 11", 100, 0., 5000.0); - fEcalRLayer21 = new TH1F("fEcalRLayer21", "Ecal R layer 21", 100, 0., 5000.0); - fEcalRLayerNorm = new TH1F("fEcalRLayerNorm", "Ecal R layer Norm", 100, 0., 5000.0); - m_geoSvc = serviceLocator()->service(m_geoSvcName); -}; - StatusCode DDCaloDigi::initialize() { - //fHcalCvsE = new TH2F("fHcalCvsE", "Hcal time profile cor",100, 0., 500.0,100,0.,10.) - _strip_virt_cells = -999; - _countWarnings = 0; + m_countWarnings = 0; try { dd4hep::rec::LayeredCalorimeterData* ecalBarrelData = @@ -138,422 +100,391 @@ StatusCode DDCaloDigi::initialize() { const std::vector& ecalBarrelLayers = ecalBarrelData->layers; const std::vector& ecalEndcapLayers = ecalEndcapData->layers; - // determine geometry of ECAL - int symmetry = ecalBarrelData->inner_symmetry; - _zOfEcalEndcap = ecalEndcapData->extent[2] / dd4hep::mm; + // Determine geometry of ECAL + int symmetry = ecalBarrelData->inner_symmetry; + m_zOfEcalEndcap = ecalEndcapData->extent[2] / dd4hep::mm; // Determine ECAL polygon angles - // Store radial vectors perpendicular to stave layers in _ecalBarrelStaveDir + // Store radial vectors perpendicular to stave layers in m_ecalBarrelStaveDir // ASSUMES Mokka Stave numbering 0 = top, then numbering increases anti-clockwise if (symmetry > 1) { float nFoldSymmetry = static_cast(symmetry); float phi0 = ecalBarrelData->phi0 / dd4hep::rad; for (int i = 0; i < symmetry; ++i) { - float phi = phi0 + i * dd4hep::twopi / nFoldSymmetry; - _barrelStaveDir[i][0] = cos(phi); - _barrelStaveDir[i][1] = sin(phi); + float phi = phi0 + i * dd4hep::twopi / nFoldSymmetry; + m_barrelStaveDir[i][0] = cos(phi); + m_barrelStaveDir[i][1] = sin(phi); } } for (unsigned int i = 0; i < ecalBarrelLayers.size(); ++i) { - _barrelPixelSizeT[i] = ecalBarrelLayers[i].cellSize0; - _barrelPixelSizeZ[i] = ecalBarrelLayers[i].cellSize1; - streamlog_out(DEBUG) << "barrel pixel size " << i << " " << _barrelPixelSizeT[i] << " " << _barrelPixelSizeZ[i] - << endl; + m_barrelPixelSizeT[i] = ecalBarrelLayers[i].cellSize0; + m_barrelPixelSizeZ[i] = ecalBarrelLayers[i].cellSize1; + debug() << "barrel pixel size " << i << " " << m_barrelPixelSizeT[i] << " " << m_barrelPixelSizeZ[i] << std::endl; } for (unsigned int i = 0; i < ecalEndcapLayers.size(); ++i) { - _endcapPixelSizeX[i] = ecalEndcapLayers[i].cellSize0; - _endcapPixelSizeY[i] = ecalEndcapLayers[i].cellSize1; - streamlog_out(DEBUG) << "endcap pixel size " << i << " " << _endcapPixelSizeX[i] << " " << _endcapPixelSizeY[i] - << endl; + m_endcapPixelSizeX[i] = ecalEndcapLayers[i].cellSize0; + m_endcapPixelSizeY[i] = ecalEndcapLayers[i].cellSize1; + debug() << "endcap pixel size " << i << " " << m_endcapPixelSizeX[i] << " " << m_endcapPixelSizeY[i] << std::endl; } - _strip_virt_cells = _ecalStrip_default_nVirt; - streamlog_out(WARNING) << "taking number of virtual cells from steering file (FIXME!): " << _strip_virt_cells - << endl; - _ecalLayout = _ecal_deafult_layer_config; - streamlog_out(WARNING) << "taking layer layout from steering file (FIXME): " << _ecalLayout << endl; + m_ecalLayout = m_ecal_deafult_layer_config; + warning() << "taking layer layout from steering file (FIXME): " << m_ecalLayout << std::endl; } catch (std::exception& e) { - streamlog_out(ERROR) << "Could not get ECAL parameters from DD4hep!" << endl; + error() << "Could not get ECAL parameters from DD4hep!" << endmsg; } - //convert ECAL thresholds to GeV units - if (_unitThresholdEcal.compare("GeV") == 0) { - //ECAL threshold unit is GeV, do nothing - } else if (_unitThresholdEcal.compare("MIP") == 0) { - //ECAL threshold unit is MIP, convert via MIP2GeV - _thresholdEcal *= _calibEcalMip; - } else if (_unitThresholdEcal.compare("px") == 0) { - //ECAL threshold unit is pixels, convert via MIP2GeV and lightyield - _thresholdEcal *= _ecal_PPD_pe_per_mip * _calibEcalMip; + // Convert ECAL thresholds to GeV units + if (m_unitThresholdEcal.value().compare("GeV") == 0) { + // ECAL threshold unit is GeV, do nothing + } else if (m_unitThresholdEcal.value().compare("MIP") == 0) { + // ECAL threshold unit is MIP, convert via MIP2GeV + m_thresholdEcal.value() *= m_calibEcalMip.value(); + } else if (m_unitThresholdEcal.value().compare("px") == 0) { + // ECAL threshold unit is pixels, convert via MIP2GeV and lightyield + m_thresholdEcal.value() *= m_ecal_PPD_pe_per_mip.value() * m_calibEcalMip.value(); } else { - streamlog_out(ERROR) << "could not identify ECAL threshold unit. Please use \"GeV\", \"MIP\" or \"px\"! Aborting." - << std::endl; + error() << "Could not identify ECAL threshold unit. Please use \"GeV\", \"MIP\" or \"px\"! Aborting." << endmsg; assert(0); } - //convert HCAL thresholds to GeV units - if (_unitThresholdHcal.compare("GeV") == 0) { - //HCAL threshold unit is GeV, do nothing - } else if (_unitThresholdHcal.compare("MIP") == 0) { - //HCAL threshold unit is MIP, convert via MIP2GeV - for (unsigned int i = 0; i < _thresholdHcal.size(); i++) { - _thresholdHcal.at(i) *= _calibHcalMip; + // Convert HCAL thresholds to GeV units + if (m_unitThresholdHcal.value().compare("GeV") == 0) { + // HCAL threshold unit is GeV, do nothing + } else if (m_unitThresholdHcal.value().compare("MIP") == 0) { + // HCAL threshold unit is MIP, convert via MIP2GeV + for (unsigned int i = 0; i < m_thresholdHcal.value().size(); i++) { + m_thresholdHcal.value()[i] *= m_calibHcalMip.value(); } - } else if (_unitThresholdHcal.compare("px") == 0) { - //HCAL threshold unit is pixels, convert via MIP2GeV and lightyield - for (unsigned int i = 0; i < _thresholdHcal.size(); i++) { - _thresholdHcal.at(i) *= _hcal_PPD_pe_per_mip * _calibHcalMip; + } else if (m_unitThresholdHcal.value().compare("px") == 0) { + // HCAL threshold unit is pixels, convert via MIP2GeV and lightyield + for (unsigned int i = 0; i < m_thresholdHcal.size(); i++) { + m_thresholdHcal[i] *= m_hcal_PPD_pe_per_mip.value() * m_calibHcalMip.value(); } } else { - streamlog_out(ERROR) << "could not identify HCAL threshold unit. Please use \"GeV\", \"MIP\" or \"px\"! Aborting." - << std::endl; + error() << "Could not identify HCAL threshold unit. Please use \"GeV\", \"MIP\" or \"px\"! Aborting." << endmsg; assert(0); } - // set up the scintillator/MPPC digitiser - _scEcalDigi = std::unique_ptr(new DDScintillatorPpdDigi()); - _scEcalDigi->setPEperMIP(_ecal_PPD_pe_per_mip); - _scEcalDigi->setCalibMIP(_calibEcalMip); - _scEcalDigi->setNPix(_ecal_PPD_n_pixels); - _scEcalDigi->setRandomMisCalibNPix(_ecal_misCalibNpix); - _scEcalDigi->setPixSpread(_ecal_pixSpread); - _scEcalDigi->setElecNoise(_ecal_elec_noise); - _scEcalDigi->setElecRange(_ecalMaxDynMip); - cout << "ECAL sc digi:" << endl; - _scEcalDigi->printParameters(); - - _scHcalDigi = std::unique_ptr(new DDScintillatorPpdDigi()); - _scHcalDigi->setPEperMIP(_hcal_PPD_pe_per_mip); - _scHcalDigi->setCalibMIP(_calibHcalMip); - _scHcalDigi->setNPix(_hcal_PPD_n_pixels); - _scHcalDigi->setRandomMisCalibNPix(_hcal_misCalibNpix); - _scHcalDigi->setPixSpread(_hcal_pixSpread); - _scHcalDigi->setElecNoise(_hcal_elec_noise); - _scHcalDigi->setElecRange(_hcalMaxDynMip); - cout << "HCAL sc digi:" << endl; - _scHcalDigi->printParameters(); - - //set up the random engines for ecal and hcal dead cells: (could use a steering parameter though) - if (_deadCellEcal_keep) { - _randomEngineDeadCellEcal = new CLHEP::MTwistEngine(0, 0); + // Set up the scintillator/MPPC digitizer + m_scEcalDigi = std::unique_ptr(new DDScintillatorPpdDigi()); + m_scEcalDigi->setPEperMIP(m_ecal_PPD_pe_per_mip); + m_scEcalDigi->setCalibMIP(m_calibEcalMip); + m_scEcalDigi->setNPix(m_ecal_PPD_n_pixels); + m_scEcalDigi->setRandomMisCalibNPix(m_ecal_misCalibNpix); + m_scEcalDigi->setPixSpread(m_ecal_pixSpread); + m_scEcalDigi->setElecNoise(m_ecal_elec_noise); + m_scEcalDigi->setElecRange(m_ecalMaxDynMip); + std::cout << "ECAL sc digi:" << std::endl; + m_scEcalDigi->printParameters(); + + m_scHcalDigi = std::unique_ptr(new DDScintillatorPpdDigi()); + m_scHcalDigi->setPEperMIP(m_hcal_PPD_pe_per_mip); + m_scHcalDigi->setCalibMIP(m_calibHcalMip); + m_scHcalDigi->setNPix(m_hcal_PPD_n_pixels); + m_scHcalDigi->setRandomMisCalibNPix(m_hcal_misCalibNpix); + m_scHcalDigi->setPixSpread(m_hcal_pixSpread); + m_scHcalDigi->setElecNoise(m_hcal_elec_noise); + m_scHcalDigi->setElecRange(m_hcalMaxDynMip); + std::cout << "HCAL sc digi:" << std::endl; + m_scHcalDigi->printParameters(); + + // Set up the random engines for ECAL and HCAL dead cells: (could use a steering parameter though) + if (m_deadCellEcal_keep) { + m_randomEngineDeadCellEcal = new CLHEP::MTwistEngine(0, 0); } else { - _randomEngineDeadCellEcal = 0; + m_randomEngineDeadCellEcal = 0; } - - if (_deadCellHcal_keep) { - _randomEngineDeadCellHcal = new CLHEP::MTwistEngine(0, 0); + if (m_deadCellHcal_keep) { + m_randomEngineDeadCellHcal = new CLHEP::MTwistEngine(0, 0); } else { - _randomEngineDeadCellHcal = 0; + m_randomEngineDeadCellHcal = 0; } + + return StatusCode::SUCCESS; } -retType DDCaloDigi::operator()( - const std::map& simCaloHitCollections, - const edm4hep::EventHeaderCollection& headers) const { - debug() << " process event : " << headers[0].getEventNumber() << " - run " << headers[0].getRunNumber() +retType DDCaloDigi::operator()(const edm4hep::SimCalorimeterHitCollection& simCaloHitCollection, + const edm4hep::EventHeaderCollection& headers) const { + debug() << " process event : " << headers[0].getEventNumber() << " - run " << headers[0].getRunNumber() << endmsg; - << endmsg; // headers[0].getRunNumber(),headers[0].getEventNumber() + auto col = edm4hep::CalorimeterHitCollection(); // create output CalorimeterHit collection + auto Relcol = + edm4hep::CaloHitSimCaloHitLinkCollection(); // create relation collection CalorimeterHit-SimCalorimeterHit - //auto chschcol = edm4hep::CalorimeterHitCollection(); //collection - std::map outputCollections; - auto Relcol = edm4hep::MCRecoCaloAssociationCollection(); // Relation collection CalorimeterHit, SimCalorimeterHit + // decide on this event's correlated miscalibration + //float m_event_correl_miscalib_ecal = CLHEP::RandGauss::shoot(1.0, m_misCalibEcal_correl.value()); + //float m_event_correl_miscalib_hcal = CLHEP::RandGauss::shoot(1.0, m_misCalibHcal_correl.value()); - // copy the flags from the input collection - //_flag.setBit(LCIO::CHBIT_LONG); - //_flag.setBit(LCIO::RCHBIT_TIME); //store timing on output hits. + std::string colName = m_inputLocations[0][0].key(); // take input collection name + std::cout << "looking for collection: " << colName << std::endl; - // decide on this event's correlated miscalibration - _event_correl_miscalib_ecal = CLHEP::RandGauss::shoot(1.0, _misCalibEcal_correl); - _event_correl_miscalib_hcal = CLHEP::RandGauss::shoot(1.0, _misCalibHcal_correl); + if (colName.find("dummy") != std::string::npos) { + debug() << "Ignoring input collection name (looks like dummy name)" << colName << endmsg; + } + + CHT::Layout caloLayout = layoutFromString(colName); + + //FIXME: take this from the metadata !!! + //std::string initString = m_geoSvc->constantAsString(m_encodingStringVariable.value()); + std::string initString = "system:5,side:2,module:8,stave:4,layer:9,submodule:4,x:32:-16,y:-16"; + + dd4hep::DDSegmentation::BitFieldCoder bitFieldCoder(initString); // check! + // check if decoder contains "layer" + + std::vector m_calHitsByStaveLayer[MAX_STAVES][MAX_LAYERS]; + std::vector m_calHitsByStaveLayerModule[MAX_STAVES][MAX_LAYERS]; + + int orientation = getStripOrientationFromColName(colName); // // * Reading Collections of ECAL Simulated Hits * // + if (m_inputColIsECAL) { + //loop over all SimCalorimeterHits in this collection + for (const auto& hit : simCaloHitCollection) { + const int cellID = hit.getCellID(); + float energy = hit.getEnergy(); + + // apply threshold cut + if (energy > m_thresholdEcal) { + unsigned int layer = bitFieldCoder.get(cellID, "layer"); + unsigned int stave = bitFieldCoder.get(cellID, "stave"); + unsigned int module = bitFieldCoder.get(cellID, "module"); + + // check that layer and assumed layer type are compatible + checkConsistency(colName, layer); + + float x = hit.getPosition()[0]; + float y = hit.getPosition()[1]; + float z = hit.getPosition()[2]; + float r = sqrt(x * x + y * y + z * z); + float rxy = sqrt(x * x + y * y); + float cost = fabs(z) / r; + + // fill ECAL Layer histograms + if (z > 0) { + if (layer == 1) + ++fEcalLayer1[{x, y}]; + if (layer == 11) + ++fEcalLayer11[{x, y}]; + if (layer == 21) + ++fEcalLayer21[{x, y}]; + if (layer == 1) + ++fEcalRLayer1[rxy]; + if (layer == 11) + ++fEcalRLayer11[rxy]; + if (layer == 21) + ++fEcalRLayer21[rxy]; + } - for (auto const& inputPair : simCaloHitCollections) { - std::string colName = inputPair.first; - auto& ecalCol = outputCollections["new" + colName]; - auto const& inputCaloHitCollection = inputPair.second; // (or is it ->second ??) - streamlog_out(DEBUG) << "looking for collection: " << colName << endl; - - if (colName.find("dummy") != string::npos) { - streamlog_out(DEBUG) << "ignoring input ECAL collection name (looks like dummy name)" << colName << endl; - continue; - } - - //fg: need to establish the subdetetcor part here - // use collection name as cellID does not seem to have that information - CHT::Layout caloLayout = layoutFromString(colName); - - try { - // auto col = evt->getCollection( colName.c_str() ) ; - const auto initString = cellIDHandle.get(); - //std::string initString = m_geoSvc->constantAsString(m_encodingStringVariable.value()); - dd4hep::DDSegmentation::BitFieldCoder bitFieldCoder(initString); // check! - // check if decoder contains "layer" - //CellIDDecoder idDecoder( col ); - - // create new collection - //LCCollectionVec *ecalcol = new LCCollectionVec(LCIO::CALORIMETERHIT); - //auto ecalcol = edm4hep::CalorimeterHit(); - //ecalcol->setFlag(_flag.getFlag()); - - // if making gap corrections clear the vectors holding pointers to calhits - if (_ecalGapCorrection != 0) { - for (int is = 0; is < MAX_STAVES; is++) { - for (int il = 0; il < MAX_LAYERS; il++) { - _calHitsByStaveLayer[is][il].clear(); - _calHitsByStaveLayerModule[is][il].clear(); + // evaluate the ECAL calibration coefficient + float calibr_coeff = 1.; + if (m_digitalEcal) { + calibr_coeff = this->digitalEcalCalibCoeff(layer); + if (m_mapsEcalCorrection) { + if (caloLayout == CHT::barrel) { + float correction = 1.1387 - 0.068 * cost - 0.191 * cost * cost; + calibr_coeff /= correction; + } else { + float correction = 0.592 + 0.590 * cost; + calibr_coeff /= correction; + } } + } else { + calibr_coeff = this->analogueEcalCalibCoeff(layer); + } + // if(fabs(hit->getPosition()[2])>=_zOfEcalEndcap)calibr_coeff *= m_ecalEndcapCorrectionFactor; + if (caloLayout != CHT::endcap) { + calibr_coeff *= m_ecalEndcapCorrectionFactor; // more robust } - } - - // deal with strips split into virtual cells - // if this collection is a strip which has been split into virtual cells, they need to be recombined - int orientation = getStripOrientationFromColName(colName); - if (orientation != SQUARE && getNumberOfVirtualCells() > 1) { - auto col = combineVirtualStripCells(ecalCol, caloLayout == CHT::barrel, orientation); - } - // for (int j(0); j < numElements; ++j) { - for (const auto& hit : inputCaloHitCollection) { - // SimCalorimeterHit * hit = dynamic_cast( col->getElementAt( j ) ) ; - float energy = hit.getEnergy(); - - // apply threshold cut - if (energy > _thresholdEcal) { - int cellid = hit.getCellID(); - int layer = idDecoder(hit)["layer"]; - int stave = idDecoder(hit)["stave"]; - int module = idDecoder(hit)["module"]; - - // check that layer and assumed layer type are compatible - checkConsistency(colName, layer); - - // save hits by module/stave/layer if required later - float calibr_coeff(1.); - float x = hit.getPosition()[0]; - float y = hit.getPosition()[1]; - float z = hit.getPosition()[2]; - float r = sqrt(x * x + y * y + z * z); - float rxy = sqrt(x * x + y * y); - float cost = fabs(z) / r; - - if (z > 0 && _histograms) { - if (layer == 1) - fEcalLayer1->Fill(x, y); - if (layer == 11) - fEcalLayer11->Fill(x, y); - if (layer == 21) - fEcalLayer21->Fill(x, y); - if (layer == 1) - fEcalRLayer1->Fill(rxy); - if (layer == 11) - fEcalRLayer11->Fill(rxy); - if (layer == 21) - fEcalRLayer21->Fill(rxy); + // apply timing cut for ECAL + if (m_useEcalTiming) { + float ecalTimeWindowMax; + if (caloLayout == CHT::barrel) { // current SimHit is in barrel, use barrel timing cut + ecalTimeWindowMax = m_ecalBarrelTimeWindowMax; + } else { // current simhit is not in barrel, use endcap timing cut + ecalTimeWindowMax = m_ecalEndcapTimeWindowMax; } - if (_digitalEcal) { - calibr_coeff = this->digitalEcalCalibCoeff(layer); - if (_mapsEcalCorrection) { - if (caloLayout == CHT::barrel) { - float correction = 1.1387 - 0.068 * cost - 0.191 * cost * cost; - calibr_coeff /= correction; - } else { - float correction = 0.592 + 0.590 * cost; - calibr_coeff /= correction; - } + + float dt = r / 300. - 0.1; + auto ecalSingleHit = hit.getContributions(); + int count = 0; + float eCellInTime = 0.; + float eCellOutput = 0.; + const unsigned int n = ecalSingleHit.size(); + std::vector used(n, false); + + for (unsigned int i_t = 0; i_t < n; i_t++) { // loop over all subhits + float timei = ecalSingleHit[i_t].getTime(); // absolute hit timing of current subhit + float energyi = ecalSingleHit[i_t].getEnergy(); // energy of current subhit + float energySum = 0; + + float deltat = 0; + if (m_ecalCorrectTimesForPropagation) { + deltat = dt; //deltat now carries hit timing correction. + } + if (timei - deltat > m_ecalTimeWindowMin.value() && timei - deltat < ecalTimeWindowMax) { + float ecor = energyi * calibr_coeff; + eCellInTime += ecor; } - } else { - calibr_coeff = this->analogueEcalCalibCoeff(layer); - } - // if(fabs(hit->getPosition()[2])>=_zOfEcalEndcap)calibr_coeff *= _ecalEndcapCorrectionFactor; - if (caloLayout != CHT::barrel) - calibr_coeff *= _ecalEndcapCorrectionFactor; // more robust - - // if you want to understand the timing cut code, please refer to the hcal timing cut further below. it is functionally identical, but has comments, explanations and excuses. - if (_useEcalTiming) { - float ecalTimeWindowMax = _ecalEndcapTimeWindowMax; - if (caloLayout == CHT::barrel) - ecalTimeWindowMax = _ecalBarrelTimeWindowMax; - float dt = r / 300. - 0.1; - - std::vector used(n, false); - //for(unsigned int i =0; i _ecalTimeWindowMin && timei - deltat < ecalTimeWindowMax) { - float ecor = energyi * calibr_coeff; - eCellInTime += ecor; - } - if (!used[i_t]) { - // merge with other hits? - used[i_t] = true; - for (unsigned int j_t = i_t + 1; j_t < n; j_t++) { - if (!used[j_t]) { - float timej = EsingleHit[j_t].getTime(); - float energyj = EsingleHit[j_t].getEnergy(); - if (_ecalSimpleTimingCut) { - float deltat_ij = _ecalCorrectTimesForPropagation ? dt : 0; - if (timej - deltat_ij > _ecalTimeWindowMin && timej - deltat_ij < ecalTimeWindowMax) { - energySum += energyj; - if (timej < timei) { - timei = timej; - } + if (!used + [i_t]) { //if current subhit has not been merged with previous hits already, take current hit as starting point to merge hits + // merge with other hits? + used[i_t] = true; + for (unsigned int j_t = i_t + 1; j_t < n; j_t++) { //loop through all hits after current hit + if (!used[j_t]) { + float timej = ecalSingleHit[j_t].getTime(); + float energyj = ecalSingleHit[j_t].getEnergy(); + if (m_ecalSimpleTimingCut) { + float deltat_ij = m_ecalCorrectTimesForPropagation ? dt : 0; + if (timej - deltat_ij > m_ecalTimeWindowMin && timej - deltat_ij < ecalTimeWindowMax) { + energySum += energyj; + if (timej < timei) { + timei = timej; //use earliest hit time for simpletimingcut } - } else { - float deltat_ij = fabs(timei - timej); - if (deltat_ij < _ecalDeltaTimeHitResolution) { - if (energyj > energyi) - timei = timej; - energyi += energyj; - used[j_t] = true; + } + } else { + float deltat_ij = fabs(timei - timej); + //if this subhit is close to current subhit, add this hit's energy to + if (deltat_ij < m_ecalDeltaTimeHitResolution) { + if (energyj > energyi) { + timei = timej; } + energyi += energyj; + used[j_t] = true; } } } - if (_ecalSimpleTimingCut) { - used = vector(n, true); // mark everything as used to terminate for loop on next run - energyi += energySum; //fill energySum back into energyi to have rest of loop behave the same. - } - if (_digitalEcal) { - calibr_coeff = this->digitalEcalCalibCoeff(layer); - if (_mapsEcalCorrection) { - if (caloLayout == CHT::barrel) { - float correction = 1.1387 - 0.068 * cost - 0.191 * cost * cost; - calibr_coeff /= correction; - } else { - float correction = 0.592 + 0.590 * cost; - calibr_coeff /= correction; - } + } + if (m_ecalSimpleTimingCut) { + used = std::vector(n, true); // mark everything as used to terminate for loop on next run + energyi += energySum; // fill energySum back into energyi to have rest of loop behave the same. + } + + // variables and their behaviour at this point: + // if SimpleTimingCut == false + // energyi carries the sum of subhit energies within +- one hcal time resolution - the timecluster energy. + // timei carries something vaguely similar to the central hit time of the merged subhits + + // if SimpleTimingCut == true + // energyi carries the sum of subhit energies within timeWindowMin and timeWindowMax + // timei carries the time of the earliest hit within this window + + if (m_digitalEcal) { + calibr_coeff = this->digitalEcalCalibCoeff(layer); + if (m_mapsEcalCorrection) { + if (caloLayout == CHT::barrel) { + float correction = 1.1387 - 0.068 * cost - 0.191 * cost * cost; + calibr_coeff /= correction; + } else { + float correction = 0.592 + 0.590 * cost; + calibr_coeff /= correction; } - } else { - calibr_coeff = this->analogueEcalCalibCoeff(layer); - } - // if(fabs(hit->getPosition()[2])>=_zOfEcalEndcap)calibr_coeff *= _ecalEndcapCorrectionFactor; - if (caloLayout != CHT::barrel) - calibr_coeff *= _ecalEndcapCorrectionFactor; // more robust - - if (_histograms) { - fEcal->Fill(timei, energyi * calibr_coeff); - fEcalC->Fill(timei - dt, energyi * calibr_coeff); - fEcalC1->Fill(timei - dt, energyi * calibr_coeff); - fEcalC2->Fill(timei - dt, energyi * calibr_coeff); } + } else { + calibr_coeff = this->analogueEcalCalibCoeff(layer); + } + // if(fabs(hit->getPosition()[2])>=_zOfEcalEndcap)calibr_coeff *= _ecalEndcapCorrectionFactor; + if (caloLayout != CHT::barrel) { + calibr_coeff *= m_ecalEndcapCorrectionFactor; // more robust + } - // apply extra energy digitisation effects - energyi = ecalEnergyDigi(energyi, cellID); - - if (energyi > _thresholdEcal) { - float timeCor = 0; - if (_ecalCorrectTimesForPropagation) - timeCor = dt; - timei = timei - timeCor; - if (timei > _ecalTimeWindowMin && timei < ecalTimeWindowMax) { - count++; - // CalorimeterHitImpl * calhit = new CalorimeterHitImpl(); - edm4hep::MutableCalorimeterHit calHit = ecalCol.create(); - if (_ecalGapCorrection != 0) { - _calHitsByStaveLayer[stave][layer].push_back(calhit); - _calHitsByStaveLayerModule[stave][layer].push_back(module); - } - calhit.setCellID(cellid); - - if (_digitalEcal) { - calhit.setEnergy(calibr_coeff); - } else { - calhit.setEnergy(calibr_coeff * energyi); - // calhit->setEnergy(energyi); - } - - eCellOutput += energyi * calibr_coeff; + // fill the ECAL time profile histograms (weighted by the energy) + fEcal[timei] += energyi * calibr_coeff; + fEcalC[timei - dt] += energyi * calibr_coeff; + fEcalC1[timei - dt] += energyi * calibr_coeff; + fEcalC2[timei - dt] += energyi * calibr_coeff; - calhit.setTime(timei); - calhit.setPosition(hit.getPosition()); - calhit.setType(CHT(CHT::em, CHT::ecal, caloLayout, layer)); - calhit.setRawHit(hit); + // apply extra energy digitisation effects + energyi = ecalEnergyDigi(energyi, cellID); // this only uses the current subhit "timecluster"! - // Set relation with LCRelationNavigator - auto caloRel = Relcol.create(); - caloRel.setRec(calHit); - caloRel.setSim(hit); - // calohitNav.addRelation(calhit, hit, 1.0); + if (energyi > m_thresholdEcal) { // now would be the correct time to do threshold comparison + float timeCor = 0; + if (m_ecalCorrectTimesForPropagation) { + timeCor = dt; + } + timei = timei - timeCor; + if (timei > m_ecalTimeWindowMin && + timei < ecalTimeWindowMax) { // if current subhit timecluster is within specified timing window, + // create new CalorimeterHit and add to collections etc. + count++; + edm4hep::MutableCalorimeterHit calHit = col.create(); + if (m_ecalGapCorrection != 0) { + m_calHitsByStaveLayer[stave][layer].push_back(&calHit); + m_calHitsByStaveLayerModule[stave][layer].push_back(module); + } + calHit.setCellID(cellID); + if (m_digitalEcal) { + calHit.setEnergy(calibr_coeff); } else { - // if(caloLayout==CHT::barrel)std::cout << " Drop ECAL Barrel hit : " << timei << " " << calibr_coeff*energyi << std::endl; + calHit.setEnergy(calibr_coeff * energyi); } + + eCellOutput += energyi * calibr_coeff; + + calHit.setTime(timei); + calHit.setPosition(hit.getPosition()); + calHit.setType(CHT(CHT::em, CHT::ecal, caloLayout, layer)); + + // set relation with CaloHitSimCaloHitLinkCollection + auto ecaloRel = Relcol.create(); + ecaloRel.setFrom(calHit); + ecaloRel.setTo(hit); + + } else { + // if(caloLayout==CHT::barrel) std::cout << " Drop ECAL Barrel hit : " << timei << " " << calibr_coeff*energyi << std::endl; } } } - } else { // don't use timing - // CalorimeterHitImpl * calhit = new CalorimeterHitImpl(); - edm4hep::MutableCalorimeterHit calHit = ecalCol.create(); - if (_ecalGapCorrection != 0) { - _calHitsByStaveLayer[stave][layer].push_back(calhit); - _calHitsByStaveLayerModule[stave][layer].push_back(module); - } - float energyi = hit.getEnergy(); + } + } else { // if timing cut is not used + edm4hep::MutableCalorimeterHit calHit = col.create(); + if (m_ecalGapCorrection != 0) { + m_calHitsByStaveLayer[stave][layer].push_back(&calHit); + m_calHitsByStaveLayerModule[stave][layer].push_back(module); + } + float energyi = hit.getEnergy(); - // apply extra energy digitisation effects - energyi = ecalEnergyDigi(energyi, cellid); + // apply extra energy digitisation effects + energyi = ecalEnergyDigi(energyi, cellID); - calhit.setCellID(cellid); - if (_digitalEcal) { - calhit.setEnergy(calibr_coeff); - } else { - calhit.setEnergy(calibr_coeff * energyi); - } - calhit.setTime(0); - calhit.setPosition(hit.getPosition()); - calhit.setType(CHT(CHT::em, CHT::ecal, caloLayout, layer)); - calhit.setRawHit(hit); - - // Set relation with LCRelationNavigator - auto caloRel = Relcol.create(); - caloRel.setRec(calHit); - caloRel.setSim(hit); - //calohitNav.addRelation(calhit, hit, 1.0); - } // timing if...else - - // std::cout << hit->getTimeCont(0) << " count = " << count << " E ECAL = " << energyCal << " - " << eCellInTime << " - " << eCellOutput << std::endl; - } // energy threshold - } - // if requested apply gap corrections in ECAL ? - if (_ecalGapCorrection != 0) - this->fillECALGaps(); - // add ECAL collection to event - // ecalcol->parameters().setValue(LCIO::CellIDEncoding,initString); - return std::make_tuple(std::move(outputCollections), std::move(caloRel)); - //evt->addCollection(ecalcol,_outputEcalCollections[i].c_str()); - } catch (DataNotAvailableException& e) { - streamlog_out(DEBUG) << "could not find input ECAL collection " << colName << std::endl; + calHit.setCellID(cellID); + if (m_digitalEcal) { + calHit.setEnergy(calibr_coeff); + } else { + calHit.setEnergy(calibr_coeff * energyi); + } + calHit.setTime(0); + calHit.setPosition(hit.getPosition()); + calHit.setType(CHT(CHT::em, CHT::ecal, caloLayout, layer)); + + // set relation with CaloHitSimCaloHitLinkCollection + auto ecaloRel = Relcol.create(); + ecaloRel.setFrom(calHit); + ecaloRel.setTo(hit); + } // timing if...else end + + //std::cout << hit->getTimeCont(0) << " count = " << count << " E ECAL = " << energyCal << " - " << eCellInTime << " - " << eCellOutput << std::endl; + } // energy threshold end + } // hits loop end + + // if requested apply gap corrections in ECAL ? + if (m_ecalGapCorrection != 0) { + this->fillECALGaps(m_calHitsByStaveLayer, m_calHitsByStaveLayerModule); } - } - //} - if (_histograms) { // fill normalisation of HCAL occupancy plots for (float x = 15; x < 3000; x += 30) { for (float y = 15; y < 3000; y += 30) { if (x > 430 || y > 430) { float r = sqrt(x * x + y * y); - fHcalRLayerNorm->Fill(r, 4.); + fHcalRLayerNorm[r] += 4.; } } } @@ -562,284 +493,211 @@ retType DDCaloDigi::operator()( for (float x = 2.5; x < 3000; x += 5) { for (float y = 2.5; y < 3000; y += 5) { float r = sqrt(x * x + y * y); - if (r > 235) - fEcalRLayerNorm->Fill(r, 4.); + if (r > 235) { + fEcalRLayerNorm[r] += 4.; + } } } - } + } // end of ECAL digitization // // * Reading HCAL Collections of Simulated Hits * // - - for (auto const& inputPair : simCaloHitCollections) { - std::string colName = inputPair.first; - auto& hcalCol = outputCollections["new" + colName]; - auto const& inputCaloHitCollection = inputPair.second; // (or is it ->second ??) - streamlog_out(DEBUG) << "looking for collection: " << colName << endl; - - if (colName.find("dummy") != string::npos) { - streamlog_out(DEBUG) << "ignoring input HCAL collection name (looks like dummy name)" << colName << endl; - continue; - } - - //fg: need to establish the subdetetcor part here - // use collection name as cellID does not seem to have that information - CHT::Layout caloLayout = layoutFromString(colName); - - try { - //LCCollection * col = evt->getCollection( _hcalCollections[i].c_str() ) ; - std::string initString = m_geoSvc->constantAsString(m_encodingStringVariable.value()); - dd4hep::DDSegmentation::BitFieldCoder bitFieldCoder(initString); // check! - - hcalcol->setFlag(_flag.getFlag()); - // for (int j(0); j < numElements; ++j) { //loop over all SimCalorimeterHits in this collection - for (const auto& hit : inputCaloHitCollection) { - float energy = hit.getEnergy(); - - if (energy > - _thresholdHcal[0] / - 2) { //preselect for SimHits with energySum>threshold. Doubtful at least, as lower energy hit might fluctuate up and still be counted - int cellid = hit.getCellID(); - float calibr_coeff(1.); - int layer = idDecoder(hit)["layer"]; - // NOTE : for a digital HCAL this does not allow for varying layer thickness - // with depth - would need a simple mod to make response proportional to layer thickness - if (_digitalHcal) { - calibr_coeff = this->digitalHcalCalibCoeff(caloLayout, energy); - } else { - calibr_coeff = this->analogueHcalCalibCoeff(caloLayout, layer); + else { + //loop over all SimCalorimeterHits in this collection + for (const auto& hit : simCaloHitCollection) { + float energy = hit.getEnergy(); + //preselect for SimHits with energySum>threshold. Doubtful at least, as lower energy hit might fluctuate up and still be counted + if (energy > m_thresholdHcal[0] / 2) { + int cellID = hit.getCellID(); + float calibr_coeff = 1.0; + unsigned int layer = bitFieldCoder.get(cellID, "layer"); + + // NOTE : for a digital HCAL this does not allow for varying layer thickness + // with depth - would need a simple mod to make response proportional to layer thickness + if (m_digitalHcal) { + calibr_coeff = this->digitalHcalCalibCoeff(caloLayout, energy); + } else { + calibr_coeff = this->analogueHcalCalibCoeff(caloLayout, layer); + } + // if(fabs(hit->getPosition()[2])>=_zOfEcalEndcap)calibr_coeff*=_hcalEndcapCorrectionFactor; + if (caloLayout != CHT::endcap) + calibr_coeff *= m_hcalEndcapCorrectionFactor; // more robust, is applied to ALL hits outside of barrel. + + //float energyCal = energy*calibr_coeff + float x = hit.getPosition()[0]; + float y = hit.getPosition()[1]; + float z = hit.getPosition()[2]; + //float r = sqrt(x*x+y*y); + if (m_useHcalTiming) { + float hcalTimeWindowMax; + if (caloLayout == CHT::barrel) { //current SimHit is in barrel, use barrel timing cut + hcalTimeWindowMax = m_hcalBarrelTimeWindowMax; + } else { //current SimHit is not in barrel, use endcap timing cut + hcalTimeWindowMax = m_hcalEndcapTimeWindowMax; } - // if(fabs(hit->getPosition()[2])>=_zOfEcalEndcap)calibr_coeff*=_hcalEndcapCorrectionFactor; - if (caloLayout != CHT::barrel) - calibr_coeff *= _hcalEndcapCorrectionFactor; // more robust, is applied to ALL hits outside of barrel. - - //float energyCal = energy*calibr_coeff - float x = hit.getPosition()[0]; - float y = hit.getPosition()[1]; - float z = hit.getPosition()[2]; - // float r = sqrt(x*x+y*y); - if (_useHcalTiming) { - float hcalTimeWindowMax; - if (caloLayout == CHT::barrel) { //current SimHit is in barrel, use barrel timing cut - hcalTimeWindowMax = _hcalBarrelTimeWindowMax; - } else { //current simhit is not in barrel, use endcap timing cut - hcalTimeWindowMax = _hcalEndcapTimeWindowMax; - } - float r = sqrt( - x * x + y * y + - z * z); //this is a crude approximation. assumes initial particle originated at the very center of the detector. - float dt = r / 300 - 0.1; //magic numbers! ~ - auto HsingleHit = hit.getContributions(); - const unsigned int n = HsingleHit.size(); - ; //number of subhits of this SimHit - - std::vector used(n, false); - - int count = 0; - - for (unsigned int i_t = 0; i_t < n; i_t++) { // loop over all subhits - float timei = HsingleHit[i_t].getTime(); //absolute hit timing of current subhit - float energyi = HsingleHit[i_t].getEnergy(); //energy of current subhit - float energySum = 0; - //float deltat = 0; - //if(_hcalCorrectTimesForPropagation)deltat=dt; - //deltat now carries hit timing correction. - //std::cout <<"outer:" << i << " " << n << std::endl; - - //idea of the following section: - //if simpletimingcut == false - //sum up hit energies which lie within one calo timing resolution to "timecluster" of current subhit - //then treat each single timecluster as one hit over threshold and digitise separately. this means there can be more than one CalorimeterHit with the same cellIDs, but different hit times (!) - // - //if simpletimingcut == true - //i'm very sorry. this is the worst code you will ever see. - //sum up hit energies within timeWindowMin and timeWindowMax, use earliest subhit in this window as hit time for resulting calohit. - //only one calorimeterhit will be generated from this. - - if (!used - [i_t]) { //if current subhit has not been merged with previous hits already, take current hit as starting point to merge hits - // merge with other hits? - used[i_t] = true; - for (unsigned int j_t = i_t; j_t < n; j_t++) { //loop through all hits after current hit - //std::cout << "inner:" << i << " " << j << " " << n << std::endl; - if (!used[j_t]) { - float timej = HsingleHit[j_t].getTime(); - float energyj = HsingleHit[j_t].getEnergy(); - // std::cout << " HCAL deltat_ij : " << deltat_ij << std::endl; - if (_hcalSimpleTimingCut) { - float deltat_ij = _hcalCorrectTimesForPropagation ? dt : 0; - if (timej - deltat_ij > _hcalTimeWindowMin && timej - deltat_ij < hcalTimeWindowMax) { - energySum += energyj; - if (timej < timei) { - timei = timej; //use earliest hit time for simpletimingcut - } + float r = sqrt( + x * x + y * y + + z * z); //this is a crude approximation. assumes initial particle originated at the very center of the detector. + float dt = r / 300 - 0.1; //magic numbers! ~ + auto hcalSingleHit = hit.getContributions(); + const unsigned int n = hcalSingleHit.size(); + + std::vector used(n, false); + + int count = 0; + + for (unsigned int i_t = 0; i_t < n; i_t++) { // loop over all subhits + float timei = hcalSingleHit[i_t].getTime(); //absolute hit timing of current subhit + float energyi = hcalSingleHit[i_t].getEnergy(); //energy of current subhit + float energySum = 0; + //float deltat = 0; ??? + //if(_hcalCorrectTimesForPropagation)deltat=dt;??? + //deltat now carries hit timing correction. ???? + //std::cout <<"outer:" << i << " " << n << std::endl; ??? + + //idea of the following section: + //if simpletimingcut == false + //sum up hit energies which lie within one calo timing resolution to "timecluster" of current subhit + //then treat each single timecluster as one hit over threshold and digitise separately. this means there can be more than one CalorimeterHit with the same cellIDs, but different hit times (!) + // + //if simpletimingcut == true + //i'm very sorry. this is the worst code you will ever see. + //sum up hit energies within timeWindowMin and timeWindowMax, use earliest subhit in this window as hit time for resulting calohit. + //only one calorimeterhit will be generated from this. + + if (!used + [i_t]) { //if current subhit has not been merged with previous hits already, take current hit as starting point to merge hits + // merge with other hits? + used[i_t] = true; + for (unsigned int j_t = i_t; j_t < n; j_t++) { //loop through all hits after current hit + //std::cout << "inner:" << i << " " << j << " " << n << std::endl; + if (!used[j_t]) { + float timej = hcalSingleHit[j_t].getTime(); + float energyj = hcalSingleHit[j_t].getEnergy(); + // std::cout << " HCAL deltat_ij : " << deltat_ij << std::endl; + if (m_hcalSimpleTimingCut) { + float deltat_ij = m_hcalCorrectTimesForPropagation ? dt : 0; + if (timej - deltat_ij > m_hcalTimeWindowMin.value() && timej - deltat_ij < hcalTimeWindowMax) { + energySum += energyj; + if (timej < timei) { + timei = timej; //use earliest hit time for simpletimingcut } - } else { - float deltat_ij = fabs(timei - timej); - if (deltat_ij < - _hcalDeltaTimeHitResolution) { //if this subhit is close to current subhit, add this hit's energy to timecluster - if (energyj > energyi) - timei = - timej; //this is probably not what was intended. i guess this should find the largest hit of one timecluster and use its hittime for the cluster, but instead it compares the current hit energy to the sum of already found hit energies - //std::cout << timei << " - " << timej << std::endl; - //std::cout << energyi << " - " << energyj << std::endl; - energyi += energyj; - used[j_t] = true; - //std::cout << timei << " " << energyi << std::endl; + } + } else { + float deltat_ij = fabs(timei - timej); + //if this subhit is close to current subhit, add this hit's energy to timecluster + if (deltat_ij < m_hcalDeltaTimeHitResolution) { + if (energyj > energyi) { + //this is probably not what was intended. i guess this should find the largest hit of one timecluster and use its hittime for the cluster, but instead it compares the current hit energy to the sum of already found hit energies + timei = timej; } + //std::cout << timei << " - " << timej << std::endl; + //std::cout << energyi << " - " << energyj << std::endl; + energyi += energyj; + used[j_t] = true; + //std::cout << timei << " " << energyi << std::endl; } } } - if (_hcalSimpleTimingCut) { - used = vector( - n, true); //mark everything as used to terminate loop. this is worse than goto. i'm sorry. - energyi += energySum; //fill energySum back into energyi to have rest of loop behave the same. - } - - //variables and their behaviour at this point: - //if SimpleTimingCut == false - //energyi carries the sum of subhit energies within +- one hcal time resolution - the timecluster energy. - //timei carries something vaguely similar to the central hit time of the merged subhits - - //if SimpleTimingCut == true - //energyi carries the sum of subhit energies within timeWindowMin and timeWindowMax - //timei carries the time of the earliest hit within this window - - // apply extra energy digitisation effects - energyi = ahcalEnergyDigi(energyi, cellid); //this only uses the current subhit "timecluster"! - - if (energyi > _thresholdHcal[0]) { //now would be the correct time to do threshold comparison - float timeCor = 0; - if (_hcalCorrectTimesForPropagation) - timeCor = dt; - timei = timei - timeCor; - if (timei > _hcalTimeWindowMin && - timei < - hcalTimeWindowMax) { //if current subhit timecluster is within specified timing window, create new CalorimeterHit and add to collections etc. - count++; - edm4hep::MutableCalorimeterHit calHit = hcalCol.create(); - calhit.setCellID(cellid); - if (_digitalHcal) { - calhit.setEnergy(calibr_coeff); - //eCellOutput+= calibr_coeff; - } else { - calhit.setEnergy(calibr_coeff * energyi); - //eCellOutput+= energyi*calibr_coeff; - } - calhit.setTime(timei); - calhit.setPosition(hit.getPosition()); - calhit.setType(CHT(CHT::had, CHT::hcal, caloLayout, layer)); - // Set relation with LCRelationNavigator - auto hcaloRel = Relcol.create(); - hcaloRel.setRec(calHit); - hcaloRel.setSim(hit); - //calohitNav.addRelation(calhit, hit, 1.0); + } + if (m_hcalSimpleTimingCut) { + used = std::vector( + n, true); //mark everything as used to terminate loop. this is worse than goto. i'm sorry. + energyi += energySum; //fill energySum back into energyi to have rest of loop behave the same. + } + //variables and their behaviour at this point: + //if SimpleTimingCut == false + //energyi carries the sum of subhit energies within +- one hcal time resolution - the timecluster energy. + //timei carries something vaguely similar to the central hit time of the merged subhits + + //if SimpleTimingCut == true + //energyi carries the sum of subhit energies within timeWindowMin and timeWindowMax + //timei carries the time of the earliest hit within this window + + // apply extra energy digitisation effects + energyi = hcalEnergyDigi(energyi, cellID); //this only uses the current subhit "timecluster"! + + if (energyi > m_thresholdHcal[0]) { //now would be the correct time to do threshold comparison + float timeCor = 0; + if (m_hcalCorrectTimesForPropagation) + timeCor = dt; + timei = timei - timeCor; + if (timei > m_hcalTimeWindowMin && + timei < + hcalTimeWindowMax) { //if current subhit timecluster is within specified timing window, create new CalorimeterHit and add to collections etc. + count++; + edm4hep::MutableCalorimeterHit calHit = col.create(); + calHit.setCellID(cellID); + if (m_digitalHcal) { + calHit.setEnergy(calibr_coeff); + //eCellOutput+= calibr_coeff; } else { - // std::cout << "Drop HCAL hit : " << timei << " " << calibr_coeff*energyi << std::endl; + calHit.setEnergy(calibr_coeff * energyi); + //eCellOutput+= energyi*calibr_coeff; } + calHit.setTime(timei); + calHit.setPosition(hit.getPosition()); + calHit.setType(CHT(CHT::had, CHT::hcal, caloLayout, layer)); + // Set relation with CaloHitSimCaloHitLinkCollection + auto hcaloRel = Relcol.create(); + hcaloRel.setFrom(calHit); + hcaloRel.setTo(hit); + } else { + //std::cout << "Drop HCAL hit : " << timei << " " << calibr_coeff*energyi << std::endl; ??? } } - //std::cout <<"hello" << std::endl; } - } else { // don't use timing - // CalorimeterHitImpl * calhit = new CalorimeterHitImpl(); - edm4hep::MutableCalorimeterHit calHit = hcalCol.create(); - calhit->setCellID(cellid); - float energyi = hit.getEnergy(); + } // end loop over all subhits + } else { // if timing cuo is not used + edm4hep::MutableCalorimeterHit calHit = col.create(); + calHit.setCellID(cellID); + float energyi = hit.getEnergy(); - // apply realistic digitisation - energyi = ahcalEnergyDigi(energyi, cellid); + // apply realistic digitisation + energyi = hcalEnergyDigi(energyi, cellID); - if (_digitalHcal) { - calhit.setEnergy(calibr_coeff); - } else { - calhit.setEnergy(calibr_coeff * energyi); - } - //eCellOutput+= energyi*calibr_coeff; - calhit.setTime(0); - calhit.setPosition(hit.getPosition()); - calhit.setType(CHT(CHT::had, CHT::hcal, caloLayout, layer)); - calhit.setRawHit(hit); - - auto hcaloRel = Relcol.create(); - hcaloRel.setRec(calHit); - hcaloRel.setSim(hit); + if (m_digitalHcal) { + calHit.setEnergy(calibr_coeff); + } else { + calHit.setEnergy(calibr_coeff * energyi); } - - // std::cout << hit->getTimeCont(0) << " count = " << count << " EHCAL = " << energyCal << " - " << eCellInTime << " - " << eCellOutput << std::endl; + //eCellOutput+= energyi*calibr_coeff; + calHit.setTime(0); + calHit.setPosition(hit.getPosition()); + calHit.setType(CHT(CHT::had, CHT::hcal, caloLayout, layer)); + + auto hcaloRel = Relcol.create(); + hcaloRel.setFrom(calHit); + hcaloRel.setTo(hit); } + // std::cout << hit->getTimeCont(0) << " count = " << count << " EHCAL = " << energyCal << " - " << eCellInTime << " - " << eCellOutput << std::endl; } - // add HCAL collection to event - // hcalcolparameters().setValue(LCIO::CellIDEncoding,initString); - // evt->addCollection(hcalcol,_outputHcalCollections[i].c_str()); - return std::make_tuple(std::move(outputCollections), std::move(hcaloRel)); - } catch (DataNotAvailableException& e) { - streamlog_out(DEBUG) << "could not find input HCAL collection " << colName << std::endl; } - } - - // Create and add relation collection for ECAL/HCAL to event - //chschcol = calohitNav.createLCCollection(); - // evt->addCollection(chschcol,_outputRelCollection.c_str()); + } // end of HCAL digitization - // _nEvt++; + return std::make_tuple(std::move(col), std::move(Relcol)); } -StatusCode DDCaloDigi::finalize() { return StatusCode::SUCCESS; } //fix - -if (_histograms) { - TFile* hfile = new TFile("calTiming.root", "recreate"); - fEcal->TH1F::Write(); - fHcal->TH1F::Write(); - fEcalC->TH1F::Write(); - fHcalC->TH1F::Write(); - fEcalC1->TH1F::Write(); - fHcalC1->TH1F::Write(); - fEcalC2->TH1F::Write(); - fHcalC2->TH1F::Write(); - //fHcalCvsE->TH2F::Write(); - fHcalLayer1->TH2F::Write(); - fHcalLayer11->TH2F::Write(); - fHcalLayer21->TH2F::Write(); - fHcalLayer31->TH2F::Write(); - fHcalLayer41->TH2F::Write(); - fHcalLayer51->TH2F::Write(); - fHcalLayer61->TH2F::Write(); - fHcalLayer71->TH2F::Write(); - fHcalRLayer1->TH1F::Write(); - fHcalRLayer11->TH1F::Write(); - fHcalRLayer21->TH1F::Write(); - fHcalRLayer31->TH1F::Write(); - fHcalRLayer41->TH1F::Write(); - fHcalRLayer51->TH1F::Write(); - fHcalRLayer61->TH1F::Write(); - fHcalRLayer71->TH1F::Write(); - fHcalRLayerNorm->TH1F::Write(); - fEcalLayer1->TH2F::Write(); - fEcalLayer11->TH2F::Write(); - fEcalLayer21->TH2F::Write(); - fEcalRLayer1->TH1F::Write(); - fEcalRLayer11->TH1F::Write(); - fEcalRLayer21->TH1F::Write(); - fEcalRLayerNorm->TH1F::Write(); - - hfile->Close(); - delete hfile; -} +StatusCode DDCaloDigi::finalize() { + //delete randomengines if needed + if (m_randomEngineDeadCellHcal != 0) { + delete m_randomEngineDeadCellHcal; + } -//delete randomengines if needed -if (_randomEngineDeadCellHcal != 0) { - delete _randomEngineDeadCellHcal; -} + if (m_randomEngineDeadCellEcal != 0) { + delete m_randomEngineDeadCellEcal; + } -if (_randomEngineDeadCellEcal != 0) { - delete _randomEngineDeadCellEcal; + return StatusCode::SUCCESS; } -void DDCaloDigi::fillECALGaps() { +void DDCaloDigi::fillECALGaps( + std::vector m_calHitsByStaveLayer[MAX_STAVES][MAX_LAYERS], + std::vector m_calHitsByStaveLayerModule[MAX_STAVES][MAX_LAYERS]) const { + const float slop = 0.25; // (mm) // Loop over hits in the Barrel // For each layer calculated differences in hit positions // Look for gaps based on expected separation of adjacent hits @@ -847,25 +705,25 @@ void DDCaloDigi::fillECALGaps() { for (int is = 0; is < MAX_STAVES; ++is) { for (int il = 0; il < MAX_LAYERS; ++il) { - if (_calHitsByStaveLayer[is][il].size() > 1) { + if (m_calHitsByStaveLayer[is][il].size() > 1) { // compare all pairs of hits just once (j>i) - for (unsigned int i = 0; i < _calHitsByStaveLayer[is][il].size() - 1; ++i) { - edm4hep::MutableCalorimeterHit hiti = _calHitsByStaveLayer[is][il][i]; - int modulei = _calHitsByStaveLayerModule[is][il][i]; - float xi = hiti.getPosition()[0]; - float yi = hiti.getPosition()[1]; - float zi = hiti.getPosition()[2]; - - for (unsigned int j = i + 1; j < _calHitsByStaveLayer[is][il].size(); ++j) { - edm4hep::MutableCalorimeterHit hitj = _calHitsByStaveLayer[is][il][j]; - int modulej = _calHitsByStaveLayerModule[is][il][j]; - float xj = hitj.getPosition()[0]; - float yj = hitj.getPosition()[1]; - float zj = hitj.getPosition()[2]; - float dz = fabs(zi - zj); + for (unsigned int i = 0; i < m_calHitsByStaveLayer[is][il].size() - 1; ++i) { + edm4hep::MutableCalorimeterHit* hiti = m_calHitsByStaveLayer[is][il][i]; + int modulei = m_calHitsByStaveLayerModule[is][il][i]; + float xi = hiti->getPosition()[0]; + float yi = hiti->getPosition()[1]; + float zi = hiti->getPosition()[2]; + + for (unsigned int j = i + 1; j < m_calHitsByStaveLayer[is][il].size(); ++j) { + edm4hep::MutableCalorimeterHit* hitj = m_calHitsByStaveLayer[is][il][j]; + int modulej = m_calHitsByStaveLayerModule[is][il][j]; + float xj = hitj->getPosition()[0]; + float yj = hitj->getPosition()[1]; + float zj = hitj->getPosition()[2]; + float dz = fabs(zi - zj); // *** BARREL CORRECTION *** - if (fabs(zi) < _zOfEcalEndcap && fabs(zj) < _zOfEcalEndcap) { + if (fabs(zi) < m_zOfEcalEndcap && fabs(zj) < m_zOfEcalEndcap) { // account for stave directions using normals // calculate difference in hit postions in z and along stave float dx = xi - xj; @@ -879,12 +737,12 @@ void DDCaloDigi::fillECALGaps() { bool mgap = false; // gaps between ECAL modules // criteria gaps in the z and t direction - float zminm = 1.0 * _barrelPixelSizeZ[il] - slop; - float zmin = 1.0 * _barrelPixelSizeZ[il] + slop; - float zmax = 2.0 * _barrelPixelSizeZ[il] - slop; - float tminm = 1.0 * _barrelPixelSizeT[il] - slop; - float tmin = 1.0 * _barrelPixelSizeT[il] + slop; - float tmax = 2.0 * _barrelPixelSizeT[il] - slop; + float zminm = 1.0 * m_barrelPixelSizeZ[il] - slop; + float zmin = 1.0 * m_barrelPixelSizeZ[il] + slop; + float zmax = 2.0 * m_barrelPixelSizeZ[il] - slop; + float tminm = 1.0 * m_barrelPixelSizeT[il] - slop; + float tmin = 1.0 * m_barrelPixelSizeT[il] + slop; + float tmax = 2.0 * m_barrelPixelSizeT[il] - slop; // criteria for gaps // WOULD BE BETTER TO USE GEAR TO CHECK GAPS ARE OF EXPECTED SIZE @@ -896,31 +754,31 @@ void DDCaloDigi::fillECALGaps() { ztgap = true; if (modulei != modulej) { - if (dz > zmin && dz < 3.0 * _barrelPixelSizeZ[il] - slop && dt < tmin) + if (dz > zmin && dz < 3.0 * m_barrelPixelSizeZ[il] - slop && dt < tmin) mgap = true; } // found a gap now apply a correction based on area of gap/area of pixel if (zgap || tgap || ztgap || mgap) { float ecor = 1.; - float f = _ecalGapCorrectionFactor; // fudge + float f = m_ecalGapCorrectionFactor; // fudge if (mgap) - f = _ecalModuleGapCorrectionFactor; + f = m_ecalModuleGapCorrectionFactor; if (zgap || mgap) - ecor = 1. + f * (dz - _barrelPixelSizeZ[il]) / 2. / _barrelPixelSizeZ[il]; + ecor = 1. + f * (dz - m_barrelPixelSizeZ[il]) / 2. / m_barrelPixelSizeZ[il]; if (tgap) - ecor = 1. + f * (dt - _barrelPixelSizeT[il]) / 2. / _barrelPixelSizeT[il]; + ecor = 1. + f * (dt - m_barrelPixelSizeT[il]) / 2. / m_barrelPixelSizeT[il]; if (ztgap) - ecor = 1. + f * (dt - _barrelPixelSizeT[il]) * (dz - _barrelPixelSizeZ[il]) / 4. / - _barrelPixelSizeT[il] / _barrelPixelSizeZ[il]; - float ei = hiti.getEnergy() * ecor; - float ej = hitj.getEnergy() * ecor; - hiti.setEnergy(ei); - hitj.setEnergy(ej); + ecor = 1. + f * (dt - m_barrelPixelSizeT[il]) * (dz - m_barrelPixelSizeZ[il]) / 4. / + m_barrelPixelSizeT[il] / m_barrelPixelSizeZ[il]; + float ei = hiti->getEnergy() * ecor; + float ej = hitj->getEnergy() * ecor; + hiti->setEnergy(ei); + hitj->setEnergy(ej); } // *** ENDCAP CORRECTION *** - } else if (fabs(zi) > _zOfEcalEndcap && fabs(zj) > _zOfEcalEndcap && dz < 100) { + } else if (fabs(zi) > m_zOfEcalEndcap && fabs(zj) > m_zOfEcalEndcap && dz < 100) { float dx = fabs(xi - xj); float dy = fabs(yi - yj); bool xgap = false; @@ -931,11 +789,11 @@ void DDCaloDigi::fillECALGaps() { // x and y need to be swapped in different staves of endcap. float pixsizex, pixsizey; if (is % 2 == 1) { - pixsizex = _endcapPixelSizeY[il]; - pixsizey = _endcapPixelSizeX[il]; + pixsizex = m_endcapPixelSizeY[il]; + pixsizey = m_endcapPixelSizeX[il]; } else { - pixsizex = _endcapPixelSizeX[il]; - pixsizey = _endcapPixelSizeY[il]; + pixsizex = m_endcapPixelSizeX[il]; + pixsizey = m_endcapPixelSizeY[il]; } float xmin = 1.0 * pixsizex + slop; @@ -953,14 +811,14 @@ void DDCaloDigi::fillECALGaps() { xygap = true; if (xgap || ygap || xygap) { - // cout <<"NewLDCCaloDigi found endcap gap, adjusting energy! " << xgap << " " << ygap << " " << xygap << " , " << il << endl; - // cout << "stave " << is << " layer " << il << endl; - // cout << " dx, dy " << dx<< " " << dy << " , sizes = " << pixsizex << " " << pixsizey << endl; - // cout << " xmin... " << xmin << " " << xminm << " " << xmax << " ymin... " << ymin << " " << yminm << " " << ymax << endl; + // std::cout <<"NewLDCCaloDigi found endcap gap, adjusting energy! " << xgap << " " << ygap << " " << xygap << " , " << il << std::endl; + // std::cout << "stave " << is << " layer " << il << std::endl; + // std::cout << " dx, dy " << dx<< " " << dy << " , sizes = " << pixsizex << " " << pixsizey << std::endl; + // std::cout << " xmin... " << xmin << " " << xminm << " " << xmax << " ymin... " << ymin << " " << yminm << " " << ymax << std::endl; // found a gap make correction float ecor = 1.; - float f = _ecalGapCorrectionFactor; // fudge + float f = m_ecalGapCorrectionFactor; // fudge if (xgap) ecor = 1. + f * (dx - pixsizex) / 2. / pixsizex; if (ygap) @@ -968,10 +826,10 @@ void DDCaloDigi::fillECALGaps() { if (xygap) ecor = 1. + f * (dx - pixsizex) * (dy - pixsizey) / 4. / pixsizex / pixsizey; - // cout << "correction factor = " << ecor << endl; + // std::cout << "correction factor = " << ecor << std::endl; - hiti.setEnergy(hiti.getEnergy() * ecor); - hitj.setEnergy(hitj.getEnergy() * ecor); + hiti->setEnergy(hiti->getEnergy() * ecor); + hitj->setEnergy(hitj->getEnergy() * ecor); } } } @@ -983,76 +841,73 @@ void DDCaloDigi::fillECALGaps() { return; } -float DDCaloDigi::digitalHcalCalibCoeff(CHT::Layout caloLayout, float energy) { +float DDCaloDigi::digitalHcalCalibCoeff(CHT::Layout caloLayout, float energy) const { float calib_coeff = 0; unsigned int ilevel = 0; - for (unsigned int ithresh = 1; ithresh < _thresholdHcal.size(); ithresh++) { + for (unsigned int ithresh = 1; ithresh < m_thresholdHcal.size(); ithresh++) { // Assume!!! hit energies are stored as floats, i.e. 1, 2 or 3 - if (energy > _thresholdHcal[ithresh]) + if (energy > m_thresholdHcal[ithresh]) ilevel = ithresh; // ilevel = 0 , 1, 2 } switch (caloLayout) { case CHT::barrel: - if (ilevel > _calibrCoeffHcalBarrel.size() - 1) { - streamlog_out(ERROR) << " Semi-digital level " << ilevel - << " greater than number of HCAL Calibration Constants (" << _calibrCoeffHcalBarrel.size() - << ")" << std::endl; + if (ilevel > m_calibrCoeffHcalBarrel.value().size() - 1) { + error() << " Semi-digital level " << ilevel << " greater than number of HCAL Calibration Constants (" + << m_calibrCoeffHcalBarrel.value().size() << ")" << std::endl; } else { - calib_coeff = _calibrCoeffHcalBarrel[ilevel]; + calib_coeff = m_calibrCoeffHcalBarrel.value()[ilevel]; } break; case CHT::endcap: - if (ilevel > _calibrCoeffHcalEndCap.size() - 1) { - streamlog_out(ERROR) << " Semi-digital level " << ilevel - << " greater than number of HCAL Calibration Constants (" << _calibrCoeffHcalEndCap.size() - << ")" << std::endl; + if (ilevel > m_calibrCoeffHcalEndcap.value().size() - 1) { + error() << " Semi-digital level " << ilevel << " greater than number of HCAL Calibration Constants (" + << m_calibrCoeffHcalEndcap.value().size() << ")" << std::endl; } else { - calib_coeff = _calibrCoeffHcalEndCap[ilevel]; + calib_coeff = m_calibrCoeffHcalEndcap.value()[ilevel]; } break; case CHT::plug: - if (ilevel > _calibrCoeffHcalOther.size() - 1) { - streamlog_out(ERROR) << " Semi-digital level " << ilevel - << " greater than number of HCAL Calibration Constants (" << _calibrCoeffHcalOther.size() - << ")" << std::endl; + if (ilevel > m_calibrCoeffHcalOther.value().size() - 1) { + error() << " Semi-digital level " << ilevel << " greater than number of HCAL Calibration Constants (" + << m_calibrCoeffHcalOther.value().size() << ")" << std::endl; } else { - calib_coeff = _calibrCoeffHcalOther[ilevel]; + calib_coeff = m_calibrCoeffHcalOther.value()[ilevel]; } break; default: - streamlog_out(ERROR) << " Unknown HCAL Hit Type " << std::endl; + error() << " Unknown HCAL Hit Type " << std::endl; break; } return calib_coeff; } -float DDCaloDigi::analogueHcalCalibCoeff(CHT::Layout caloLayout, int layer) { +float DDCaloDigi::analogueHcalCalibCoeff(CHT::Layout caloLayout, int layer) const { float calib_coeff = 0; - for (unsigned int k(0); k < _hcalLayers.size(); ++k) { + for (unsigned int k(0); k < m_hcalLayers.size(); ++k) { int min, max; if (k == 0) min = 0; else - min = _hcalLayers[k - 1]; + min = m_hcalLayers[k - 1]; - max = _hcalLayers[k]; + max = m_hcalLayers[k]; if (layer >= min && layer < max) { switch (caloLayout) { case CHT::barrel: - calib_coeff = _calibrCoeffHcalBarrel[k]; + calib_coeff = m_calibrCoeffHcalBarrel[k]; break; case CHT::endcap: - calib_coeff = _calibrCoeffHcalEndCap[k]; + calib_coeff = m_calibrCoeffHcalEndcap[k]; break; case CHT::plug: case CHT::ring: - calib_coeff = _calibrCoeffHcalOther[k]; + calib_coeff = m_calibrCoeffHcalOther[k]; break; default: - streamlog_out(ERROR) << " Unknown HCAL Hit Type " << std::endl; + error() << " Unknown HCAL Hit Type " << std::endl; break; } } @@ -1061,103 +916,103 @@ float DDCaloDigi::analogueHcalCalibCoeff(CHT::Layout caloLayout, int layer) { return calib_coeff; } -float DDCaloDigi::digitalEcalCalibCoeff(int layer) { - float calib_coeff = 0; +float DDCaloDigi::digitalEcalCalibCoeff(int layer) const { + float calib_coeff = 0.; - for (unsigned int k(0); k < _ecalLayers.size(); ++k) { + for (unsigned int k(0); k < m_ecalLayers.size(); ++k) { int min, max; if (k == 0) min = 0; else - min = _ecalLayers[k - 1]; + min = m_ecalLayers[k - 1]; - max = _ecalLayers[k]; + max = m_ecalLayers[k]; if (layer >= min && layer < max) { - calib_coeff = _calibrCoeffEcal[k]; + calib_coeff = m_calibrCoeffEcal[k]; break; } } - return calib_coeff; } -float DDCaloDigi::analogueEcalCalibCoeff(int layer) { +float DDCaloDigi::analogueEcalCalibCoeff(int layer) const { float calib_coeff = 0; // retrieve calibration constants - for (unsigned int k(0); k < _ecalLayers.size(); ++k) { + for (unsigned int k(0); k < m_ecalLayers.size(); ++k) { int min, max; if (k == 0) { min = 0; } else { - min = _ecalLayers[k - 1]; + min = m_ecalLayers[k - 1]; } - max = _ecalLayers[k]; + max = m_ecalLayers[k]; if (layer >= min && layer < max) { - calib_coeff = _calibrCoeffEcal[k]; + calib_coeff = m_calibrCoeffEcal[k]; break; } } - return calib_coeff; } -float DDCaloDigi::ecalEnergyDigi(float energy, int id) { +float DDCaloDigi::ecalEnergyDigi(float energy, int id) const { // some extra digi effects (daniel) - // controlled by _applyEcalDigi = 0 (none), 1 (silicon), 2 (scintillator) + // controlled by m_applyEcalDigi = 0 (none), 1 (silicon), 2 (scintillator) // small update for time-constant uncorrelated miscalibrations. DJ, Jan 2015 - float e_out(energy); - if (_applyEcalDigi == 1) + float e_out = energy; + if (m_applyEcalDigi == 1) { e_out = siliconDigi(energy); // silicon digi - else if (_applyEcalDigi == 2) + } else if (m_applyEcalDigi == 2) { e_out = scintillatorDigi(energy, true); // scintillator digi - + } // add electronics dynamic range // Sept 2015: Daniel moved this to the ScintillatorDigi part, so it is applied before unfolding of sipm response - // if (_ecalMaxDynMip>0) e_out = min (e_out, _ecalMaxDynMip*_calibEcalMip); + // if (_ecalMaxDynMip>0) e_out = min (e_out, _ecalMaxDynMip*_calibEcalMip); ??? // random miscalib - if (_misCalibEcal_uncorrel > 0) { + if (m_misCalibEcal_uncorrel > 0) { float miscal(0); - if (_misCalibEcal_uncorrel_keep) { - int id; - if (_ECAL_cell_miscalibs.find(id) != - _ECAL_cell_miscalibs.end()) { // this cell was previously seen, and a miscalib stored - miscal = _ECAL_cell_miscalibs[id]; + if (m_misCalibEcal_uncorrel_keep) { + int id{0}; + if (m_ECAL_cell_miscalibs.find(id) != + m_ECAL_cell_miscalibs.end()) { // this cell was previously seen, and a miscalib stored + miscal = m_ECAL_cell_miscalibs.at(id); } else { // we haven't seen this one yet, get a miscalib for it - miscal = CLHEP::RandGauss::shoot(1.0, _misCalibEcal_uncorrel); - _ECAL_cell_miscalibs[id] = miscal; + miscal = CLHEP::RandGauss::shoot(1.0, m_misCalibEcal_uncorrel); + // FIXME: this is storing miscalibration globally for a run??? + // FIXME _ECAL_cell_miscalibs[id] = miscal; ??? } } else { - miscal = CLHEP::RandGauss::shoot(1.0, _misCalibEcal_uncorrel); + miscal = CLHEP::RandGauss::shoot(1.0, m_misCalibEcal_uncorrel); } e_out *= miscal; } - if (_misCalibEcal_correl > 0) - e_out *= _event_correl_miscalib_ecal; + if (m_misCalibEcal_correl > 0) + e_out *= m_event_correl_miscalib_ecal; // random cell kill - if (_deadCellFractionEcal > 0) { - if (_deadCellEcal_keep == true) { + if (m_deadCellFractionEcal > 0) { + if (m_deadCellEcal_keep == true) { int id; - if (_ECAL_cell_dead.find(id) != _ECAL_cell_dead.end()) { // this cell was previously seen - if (_ECAL_cell_dead[id] == true) { + if (m_ECAL_cell_dead.find(id) != m_ECAL_cell_dead.end()) { // this cell was previously seen + if (m_ECAL_cell_dead.at(id) == true) { e_out = 0; } } else { // we haven't seen this one yet, get a miscalib for it - bool thisDead = (CLHEP::RandFlat::shoot(_randomEngineDeadCellEcal, .0, 1.0) < _deadCellFractionEcal); - _ECAL_cell_dead[id] = thisDead; + bool thisDead = (CLHEP::RandFlat::shoot(m_randomEngineDeadCellEcal, .0, 1.0) < m_deadCellFractionEcal); + // FIXME global map ??? + //_ECAL_cell_dead[id] = thisDead; ??? if (thisDead == true) { e_out = 0; } } } else { - if (CLHEP::RandFlat::shoot(0.0, 1.0) < _deadCellFractionEcal) + if (CLHEP::RandFlat::shoot(0.0, 1.0) < m_deadCellFractionEcal) e_out = 0; } } @@ -1165,14 +1020,14 @@ float DDCaloDigi::ecalEnergyDigi(float energy, int id) { return e_out; } -float DDCaloDigi::ahcalEnergyDigi(float energy, int id) { +float DDCaloDigi::hcalEnergyDigi(float energy, int id) const { // some extra digi effects (daniel) // controlled by _applyHcalDigi = 0 (none), 1 (scintillator/SiPM) // small update for time-constant uncorrelated miscalibrations. DJ, Jan 2015 float e_out(energy); - if (_applyHcalDigi == 1) + if (m_applyHcalDigi == 1) e_out = scintillatorDigi(energy, false); // scintillator digi // add electronics dynamic range @@ -1181,72 +1036,74 @@ float DDCaloDigi::ahcalEnergyDigi(float energy, int id) { // random miscalib // if (_misCalibHcal_uncorrel>0) e_out*=CLHEP::RandGauss::shoot( 1.0, _misCalibHcal_uncorrel ); - if (_misCalibHcal_uncorrel > 0) { + if (m_misCalibHcal_uncorrel > 0) { float miscal(0); - if (_misCalibHcal_uncorrel_keep) { + if (m_misCalibHcal_uncorrel_keep) { int id; - if (_HCAL_cell_miscalibs.find(id) != - _HCAL_cell_miscalibs.end()) { // this cell was previously seen, and a miscalib stored - miscal = _HCAL_cell_miscalibs[id]; + if (m_HCAL_cell_miscalibs.find(id) != + m_HCAL_cell_miscalibs.end()) { // this cell was previously seen, and a miscalib stored + miscal = m_HCAL_cell_miscalibs.at(id); } else { // we haven't seen this one yet, get a miscalib for it - miscal = CLHEP::RandGauss::shoot(1.0, _misCalibHcal_uncorrel); - _HCAL_cell_miscalibs[id] = miscal; + miscal = CLHEP::RandGauss::shoot(1.0, m_misCalibHcal_uncorrel); + // FIXME: same as above + //_HCAL_cell_miscalibs[id] = miscal; ??? } } else { - miscal = CLHEP::RandGauss::shoot(1.0, _misCalibHcal_uncorrel); + miscal = CLHEP::RandGauss::shoot(1.0, m_misCalibHcal_uncorrel); } e_out *= miscal; } - if (_misCalibHcal_correl > 0) - e_out *= _event_correl_miscalib_hcal; + if (m_misCalibHcal_correl > 0) + e_out *= m_event_correl_miscalib_hcal; // random cell kill - if (_deadCellFractionHcal > 0) { - if (_deadCellHcal_keep == true) { + if (m_deadCellFractionHcal > 0) { + if (m_deadCellHcal_keep == true) { int id; - if (_HCAL_cell_dead.find(id) != _HCAL_cell_dead.end()) { // this cell was previously seen - if (_HCAL_cell_dead[id] == true) { + if (m_HCAL_cell_dead.find(id) != m_HCAL_cell_dead.end()) { // this cell was previously seen + if (m_HCAL_cell_dead.at(id) == true) { e_out = 0; } } else { // we haven't seen this one yet, get a miscalib for it - bool thisDead = (CLHEP::RandFlat::shoot(_randomEngineDeadCellHcal, .0, 1.0) < _deadCellFractionHcal); - _HCAL_cell_dead[id] = thisDead; + bool thisDead = (CLHEP::RandFlat::shoot(m_randomEngineDeadCellHcal, .0, 1.0) < m_deadCellFractionHcal); + // FIXME globally dead cell map??? + //FIXME _HCAL_cell_dead[id] = thisDead; ??? if (thisDead == true) { e_out = 0; } } } else { - if (CLHEP::RandFlat::shoot(0.0, 1.0) < _deadCellFractionHcal) + if (CLHEP::RandFlat::shoot(0.0, 1.0) < m_deadCellFractionHcal) e_out = 0; } } return e_out; } -float DDCaloDigi::siliconDigi(float energy) { +float DDCaloDigi::siliconDigi(float energy) const { // applies extra digitisation to silicon hits // calculate #e-h pairs - float nehpairs = 1e9 * energy / _ehEnergy; // check units of energy! _ehEnergy is in eV, energy in GeV + float nehpairs = 1.0e9 * energy / m_ehEnergy; // check units of energy! m_ehEnergy is in eV, energy in GeV // fluctuate it by Poisson float smeared_energy = energy * CLHEP::RandPoisson::shoot(nehpairs) / nehpairs; // limited electronics dynamic range // Daniel moved electronics dyn range to here - if (_ecalMaxDynMip > 0) - smeared_energy = std::min(smeared_energy, _ecalMaxDynMip * _calibEcalMip); + if (m_ecalMaxDynMip > 0) + smeared_energy = std::min(smeared_energy, m_ecalMaxDynMip * m_calibEcalMip); // add electronics noise - if (_ecal_elec_noise > 0) - smeared_energy += CLHEP::RandGauss::shoot(0, _ecal_elec_noise * _calibEcalMip); + if (m_ecal_elec_noise > 0) + smeared_energy += CLHEP::RandGauss::shoot(0, m_ecal_elec_noise * m_calibEcalMip); return smeared_energy; } -float DDCaloDigi::scintillatorDigi(float energy, bool isEcal) { +float DDCaloDigi::scintillatorDigi(float energy, bool isEcal) const { // this applies some extra digitisation to scintillator+PPD hits (PPD=SiPM, MPPC) // - poisson fluctuates the number of photo-electrons according to #PEs/MIP // - applies PPD saturation according to #pixels @@ -1254,481 +1111,142 @@ float DDCaloDigi::scintillatorDigi(float energy, bool isEcal) { float digiEn(0); if (isEcal) { - digiEn = _scEcalDigi->getDigitisedEnergy(energy); //CHECK!!! + digiEn = m_scEcalDigi->getDigitisedEnergy(energy); //CHECK!!! } else { - digiEn = _scHcalDigi->getDigitisedEnergy(energy); + digiEn = m_scHcalDigi->getDigitisedEnergy(energy); } return digiEn; } -auto DDCaloDigi::combineVirtualStripCells(auto col, bool isBarrel, int stripOrientation) { - // combines the virtual cells in a strip - // input collection is for virtual cells - // returns collection of strips - - // daniel jeans, jan/feb 2014 - - // sanity check - if (stripOrientation == SQUARE) { - streamlog_out(ERROR) << "DDCaloDigi::combineVirtualStripCells trying to deal with silicon strips??? I do not know " - "how to do that, refusing to do anything!" - << std::endl; - return NULL; - } - - //CellIDDecoder idDecoder( col ); - // This variable is for tracking detectors only, not for calorimeters - // std::string initString = m_geoSvc->constantAsString(m_encodingStringVariable.value()); - - // get the encoding string from the original collection - const auto initString = cellIDHandle.get(); - dd4hep::DDSegmentation::BitFieldCoder bitFieldCoder(initString); - - // make the new output collection - //edm4hep::MutableSimCalorimeterHit tempstripcol = inputCaloHitCollection.create();; - auto stripcol = edm4hep::SimCalorimeterHitCollection(); - stripcol.setFlag(_flag.getFlag()); - //FIXME: add initString/ cellIDEncoding to stripcol - - // link between strip CellID and hits storage for strip hits - std::map newhits; //SORT - - // loop over input collection - int numElements = col.size(); - - // // sum energy for check - float tempenergysum(0); - for (const auto& hit : inputCaloHitCollection) { - tempenergysum += hit.getEnergy(); - } - - float scTVirtLengthBar(-99); - float scLVirtLengthBar(-99); - float scTVirtLengthEnd(-99); - float scLVirtLengthEnd(-99); - - for (const auto& hit : inputCaloHitCollection) { // loop over virtual cells - - // for now, ignore preshower layer. - // in "usual" (non-preshower) ecal collections, km1 = 0 is the first layer AFTER the preshower - int km1 = bitFieldCoder(hit)["layer"]; // FIXME - - int gearlayer = km1 + 1; // in "gearlayer", preshower is 0, first layer after absorber is 1 - - // after fix of MOKKA, this should be not needed - // int gearlayer(-99); // this is the layer index which the GEAR geometry knows about. depends on driver version... - // if ( _ecal_driver_version == 0 ) { // SEcal04 driver - // } else if ( _ecal_driver_version == 1 ) { // SEcal05 driver - // gearlayer = km1+1; - // } else { - // streamlog_out ( ERROR ) << "ERROR unknown value for ecal driver version, _ecal_driver_version=" << _ecal_driver_version << endl; - // } - // assert (gearlayer>=0); - - int m = bitFieldCoder(hit)["module"]; - int s = bitFieldCoder(hit)["stave"]; - int i_index = bitFieldCoder(hit)["x"]; - int j_index = bitFieldCoder(hit)["y"]; - - //cout << "ORIG: " << km1 << " " << m << " " << s << " " << i_index << " " << j_index << " : " << - // hit->getPosition()[0] << " "<< hit->getPosition()[1] << " "<< hit->getPosition()[2] << endl; - - // should we combine virtual cells in the I or J direction? - // - // in barrel: - // i is along slab (phi in barrel) - // j is across slab (z in barrel) - // increasing j = increasing z; increasing i = decreasing phi - // - // in endcap: - // i is across slab - // j is along slab - // different staves are rotated around beam axis. - // in +ve z endcap (m==6), - // stave 0 is at -ve x and +ve y. in this stave, i is parallel to -x, j to +y - // stave 3 is at +ve x and +ve y. in this stave, i is parallel to +y, j to +x - // -ve z endcap (m==0) is same, just everything flipped in x. - - bool collateAlongI = - isBarrel ? stripOrientation == STRIP_ALIGN_ALONG_SLAB : // I is the index along the slab (R-phi in barrel) - stripOrientation == STRIP_ALIGN_ACROSS_SLAB; // for some reason seems to be reversed in endcap...!! - - streamlog_out(DEBUG) << "isBarrel = " << isBarrel << " : stripOrientation= " << stripOrientation - << " gear layer = " << gearlayer << endl; - - // let's get the length of the virtual cell in the direction of the strip - // fixed rare crashes, and streamlined code to get virtualCellSizeAlongStrip. Sept 2014 - float virtualCellSizeAlongStrip(0); - - float* layerpixsize(0); - float* savedPixSize(0); - if (isBarrel) { - if (stripOrientation == STRIP_ALIGN_ACROSS_SLAB) { - layerpixsize = _barrelPixelSizeT; - savedPixSize = &scTVirtLengthBar; - } else { - layerpixsize = _barrelPixelSizeZ; - savedPixSize = &scLVirtLengthBar; - } - } else { - if (stripOrientation == STRIP_ALIGN_ACROSS_SLAB) { - layerpixsize = _endcapPixelSizeX; - savedPixSize = &scTVirtLengthEnd; - } else { - layerpixsize = _endcapPixelSizeY; - savedPixSize = &scLVirtLengthEnd; - } - } - - virtualCellSizeAlongStrip = layerpixsize[gearlayer]; - if (virtualCellSizeAlongStrip > 0) { // looks OK - if (*savedPixSize < 0) { // no saved size yet, keep this one - *savedPixSize = virtualCellSizeAlongStrip; - } - } else { // no valid info in gear file for this layer - if (*savedPixSize > 0) { // use saved value from other layer - virtualCellSizeAlongStrip = *savedPixSize; - } else { // look at previous layers for one with same orientation (extra check, sept 2014) - std::vector> layerTypes = getLayerConfig(); - - streamlog_out(DEBUG) << "could not get valid info from gear file..." << endl - << "looking through previous layers to get a matching orientation" << endl - << "this gear layer " << gearlayer << " type: " << layerTypes[gearlayer].first << " " - << layerTypes[gearlayer].second << endl; - - for (int il = gearlayer - 1; il >= 0; il--) { - // layer types include the preshower at posn "0" - // gearlayer has preshower as layer 0 - if (layerTypes[il] == layerTypes[gearlayer]) { // found layer with same setup - streamlog_out(DEBUG) << "found a match! " << il << " " << layerTypes[il].first << " " - << layerTypes[il].second << " : " << layerpixsize[il] << endl; - virtualCellSizeAlongStrip = layerpixsize[il]; - *savedPixSize = virtualCellSizeAlongStrip; - break; - } - } - } - } - streamlog_out(DEBUG) << "virtualCellSizeAlongStrip = " << virtualCellSizeAlongStrip << endl; - - // calculate the new strip's I,J coordinates - int i_new = collateAlongI ? i_index / getNumberOfVirtualCells() : i_index; - int j_new = collateAlongI ? j_index : j_index / getNumberOfVirtualCells(); - - // apply position-dependent energy correction - // - // relative virtual cell position within strip (0 = centre) - // distance between centre of virtual cell and centre of strip - int relativeIndx = collateAlongI ? i_index : j_index; - float relativePos = - relativeIndx % getNumberOfVirtualCells(); // this is index within strip wrt strip end (0->_strip_virt_cells-1) - relativePos -= (getNumberOfVirtualCells() - 1) / 2.0; // index wrt strip centre - relativePos *= virtualCellSizeAlongStrip; // distance wrt strip centre - - // effect of absorbtion length within scintillator - // TODO: should check the polarity is consistent with mppc position, to make sure larger response nearer to mppc.... - float energy_new = hit.getEnergy(); - - float energyNonuniformityScaling(1.); - if (_strip_abs_length > 0) { - energyNonuniformityScaling = exp(-relativePos / _strip_abs_length); - energy_new *= energyNonuniformityScaling; - } - - // calculate ID of the new Strip // FIXME??? - bitFieldCoder["system"] = m; // FIXME ???? - bitFieldCoder["module"] = m; - bitFieldCoder["stave"] = s; - bitFieldCoder["layer"] = km1; - bitFieldCoder["x"] = i_new; - bitFieldCoder["y"] = j_new; - int cellID = bitFieldCoder.getCellID; // FIXME - - // calculate position of centre of this new strip - // it depends if we are in the barrel or endcap, and in which stave - // (n.b. we could do this later, after checking for duplicates. however, keeping it here allows us to - // check that we haven't supplied wrong nVirtualCells parameter) - - float stripPos[3] = {0}; - if (isBarrel) { - if (stripOrientation == STRIP_ALIGN_ACROSS_SLAB) { // it's along z - stripPos[0] = hit.getPosition()[0]; - stripPos[1] = hit.getPosition()[1]; - stripPos[2] = hit.getPosition()[2] - relativePos; // increasing j = increasing z - } else { // it's along t - float phi = atan2(_barrelStaveDir[s][1], _barrelStaveDir[s][0]); - stripPos[0] = hit.getPosition()[0] - relativePos * cos(phi); // negative: increasing I = decreasing phi - stripPos[1] = hit.getPosition()[1] - relativePos * sin(phi); - stripPos[2] = hit.getPosition()[2]; - } - } else { // endcap - stripPos[0] = hit.getPosition()[0]; - stripPos[1] = hit.getPosition()[1]; - stripPos[2] = hit.getPosition()[2]; // z never changes - // there's almost certainly a neater way to get these different polarities in here...DJ - // ....but this is, I believe, correct! - int xsign = m == 6 ? -1 : +1; // endcaps are mirrored in x, not in y - - switch (s) { - case 0: - if (collateAlongI) - stripPos[0] += -xsign * relativePos; - else - stripPos[1] += -relativePos; - break; - case 1: - if (collateAlongI) - stripPos[1] += +relativePos; - else - stripPos[0] += -xsign * relativePos; - break; - case 2: - if (collateAlongI) - stripPos[0] += +xsign * relativePos; - else - stripPos[1] += +relativePos; - break; - case 3: - if (collateAlongI) - stripPos[1] += -relativePos; - else - stripPos[0] += +xsign * relativePos; - break; - } - - } // endcap - - // cout << "new strip pos: " << stripPos[0] << " " << stripPos[1] << " " << stripPos[2] << endl; - - edm4hep::MutableSimCalorimeterHit* newhit = nullptr; - - // check if we have already seen a virtual cell from this strip - if (newhits.find(cellid) != newhits.end()) { // already have one from this strip - - edm4hep::MutableSimCalorimeterHit htt = newhits.find(cellid)->second; // previously made hit in this stirp - // check that calculated positions are same! - bool isOK = true; - for (int i = 0; i < 3; i++) { - if (fabs(htt.getPosition()[i] - stripPos[i]) > 1) { // should be identical, but allow a little slop - isOK = false; - break; - } - } - if (!isOK) { - streamlog_out(ERROR) << "strip virtual cell merging: inconsistent position calc!" << std::endl - << "please check that the parameter ECAL_strip_nVirtualCells is correctly set..." - << getNumberOfVirtualCells() << std::endl - << " layer = (k-1) " << km1 << " (gear) " << gearlayer << endl; - - for (int i = 0; i < 3; i++) { - streamlog_out(DEBUG) << stripPos[i] << " "; - } - streamlog_out(DEBUG) << endl; - - for (int i = 0; i < 3; i++) { - streamlog_out(DEBUG) << htt.getPosition()[i] << " "; - } - streamlog_out(DEBUG) << endl; - - // std::cout << "K-1 = " << km1 ; - // std::cout << "; GEARlay = " << gearlayer; - // std::cout << "; m = " << m; - // std::cout << "; s = " << s; - // std::cout << "; i = " << i_index; - // std::cout << "; j = " << j_index << std::endl; - - assert(0); - } - - // set it to the existing one - newhit = htt; - - } else { // it's the first hit from this strip - // create a new hit for the strip - newhit = &stripcol.create(); - ; - newhits[cellid] = newhit; - newHit.setCellID(cellID); - newhit.setPosition(stripPos); - newhit.setEnergy(0); // this is added when we add MC contributions - } - - // add the MC contriutions - float eadd(0); - auto singleHit = hit.getContributions(); - - for (int ij = 0; ij < singleHit.size(); ij++) { - auto contrib = stripcol.create(); - contrib.setParticle(hit.getParticle(ij)); - contrib.setEnergy(hit.getEnergy(ij) * energyNonuniformityScaling); - contrib.setTime(hit.getTime(ij)); - contrib.setPDG(hit.getPDG(ij)); - newhit.addToContributions(contrib); - eadd += hit.getEnergy(ij) * energyNonuniformityScaling; - } - - float esum(0); - auto singlenewHit = newhit.getContributions(); - for (int ij = 0; ij < singlenewHit.size(); ij++) { - esum += newhit.getEnergy(ij); - } - } // loop over hits - - // move the hits from the temporary storage to the output collection - for (std::map::iterator itt = newhits.begin(); itt != newhits.end(); itt++) { - stripcol->addElement(itt->second); - } - - // return the new collection - return stripcol; -} - -int DDCaloDigi::getNumberOfVirtualCells() { - // if ( _strip_virt_cells < 0 ) { // not yet set, try to get from gear file - // // number of virtual cells in scintillator strip - // try { - // const gear::GearParameters& pMokka = Global::GEAR->getGearParameters("MokkaParameters"); - // std::string nVirtualMokkaS = pMokka.getStringVal("Ecal_Sc_number_of_virtual_cells"); - // _strip_virt_cells = std::atoi( nVirtualMokkaS.c_str() ); - // streamlog_out (MESSAGE) << "INFO: got number of virtual cells from gear file: " << _strip_virt_cells << endl; - // } catch(gear::UnknownParameterException &e) { // not found in gear file, get from steering file parameter... - // streamlog_out (MESSAGE) << "INFO: taking number of virtual cells from steering file (not found in gear file): " << _ecalStrip_default_nVirt << endl; - // _strip_virt_cells = _ecalStrip_default_nVirt; - // } - // } - return _strip_virt_cells; -} - -std::vector>& DDCaloDigi::getLayerConfig() { +//FIXME should return a reference to somewhere else? or do this every event? +std::vector> DDCaloDigi::getLayerConfig() const { // get the layer layout (silicon, scintillator) // first element of layerTypes is the preshower - - if (_layerTypes.size() == 0) { - for (std::string::size_type i = 0; i < _ecalLayout.size(); ++i) { + std::vector> m_layerTypes{}; + if (m_layerTypes.size() == 0) { + for (std::string::size_type i = 0; i < m_ecalLayout.size(); ++i) { // convert each element of string to integer // int type = std::atoi( &ccdd ); // this is not well done (must be null-terminated) - int etype = _ecalLayout[i] - '0'; // this is less obvious, but works... + int etype = m_ecalLayout[i] - '0'; // this is less obvious, but works... switch (etype) { // these originally defined in Mokka driver SEcalSD04 case 0: - _layerTypes.push_back(std::pair(SIECAL, SQUARE)); - _layerTypes.push_back(std::pair(SIECAL, SQUARE)); + m_layerTypes.push_back(std::pair(SIECAL, SQUARE)); + m_layerTypes.push_back(std::pair(SIECAL, SQUARE)); break; case 1: - _layerTypes.push_back(std::pair(SCECAL, STRIP_ALIGN_ALONG_SLAB)); - _layerTypes.push_back(std::pair(SCECAL, STRIP_ALIGN_ALONG_SLAB)); + m_layerTypes.push_back(std::pair(SCECAL, STRIP_ALIGN_ALONG_SLAB)); + m_layerTypes.push_back(std::pair(SCECAL, STRIP_ALIGN_ALONG_SLAB)); break; case 2: - _layerTypes.push_back(std::pair(SCECAL, STRIP_ALIGN_ACROSS_SLAB)); - _layerTypes.push_back(std::pair(SCECAL, STRIP_ALIGN_ACROSS_SLAB)); + m_layerTypes.push_back(std::pair(SCECAL, STRIP_ALIGN_ACROSS_SLAB)); + m_layerTypes.push_back(std::pair(SCECAL, STRIP_ALIGN_ACROSS_SLAB)); break; case 3: - _layerTypes.push_back(std::pair(SCECAL, STRIP_ALIGN_ALONG_SLAB)); - _layerTypes.push_back(std::pair(SCECAL, STRIP_ALIGN_ACROSS_SLAB)); + m_layerTypes.push_back(std::pair(SCECAL, STRIP_ALIGN_ALONG_SLAB)); + m_layerTypes.push_back(std::pair(SCECAL, STRIP_ALIGN_ACROSS_SLAB)); break; case 4: - _layerTypes.push_back(std::pair(SCECAL, STRIP_ALIGN_ACROSS_SLAB)); - _layerTypes.push_back(std::pair(SCECAL, STRIP_ALIGN_ALONG_SLAB)); + m_layerTypes.push_back(std::pair(SCECAL, STRIP_ALIGN_ACROSS_SLAB)); + m_layerTypes.push_back(std::pair(SCECAL, STRIP_ALIGN_ALONG_SLAB)); break; case 5: - _layerTypes.push_back(std::pair(SIECAL, SQUARE)); - _layerTypes.push_back(std::pair(SCECAL, STRIP_ALIGN_ALONG_SLAB)); + m_layerTypes.push_back(std::pair(SIECAL, SQUARE)); + m_layerTypes.push_back(std::pair(SCECAL, STRIP_ALIGN_ALONG_SLAB)); break; case 6: - _layerTypes.push_back(std::pair(SIECAL, SQUARE)); - _layerTypes.push_back(std::pair(SCECAL, STRIP_ALIGN_ACROSS_SLAB)); + m_layerTypes.push_back(std::pair(SIECAL, SQUARE)); + m_layerTypes.push_back(std::pair(SCECAL, STRIP_ALIGN_ACROSS_SLAB)); break; case 7: - _layerTypes.push_back(std::pair(SCECAL, STRIP_ALIGN_ALONG_SLAB)); - _layerTypes.push_back(std::pair(SIECAL, SQUARE)); + m_layerTypes.push_back(std::pair(SCECAL, STRIP_ALIGN_ALONG_SLAB)); + m_layerTypes.push_back(std::pair(SIECAL, SQUARE)); break; case 8: - _layerTypes.push_back(std::pair(SCECAL, STRIP_ALIGN_ACROSS_SLAB)); - _layerTypes.push_back(std::pair(SIECAL, SQUARE)); + m_layerTypes.push_back(std::pair(SCECAL, STRIP_ALIGN_ACROSS_SLAB)); + m_layerTypes.push_back(std::pair(SIECAL, SQUARE)); break; default: - streamlog_out(ERROR) << "ERROR, unknown layer type " << etype << endl; + error() << "ERROR, unknown layer type " << etype << std::endl; } } } - return _layerTypes; + return m_layerTypes; } -void DDCaloDigi::checkConsistency(std::string colName, int layer) { - if (_applyEcalDigi == 0 || _countWarnings > 20) +void DDCaloDigi::checkConsistency(std::string colName, int layer) const { + if (m_applyEcalDigi == 0 || m_countWarnings > 20) return; std::pair thislayersetup = getLayerProperties(colName, layer); - if (_applyEcalDigi == 1 && thislayersetup.first != SIECAL) { - streamlog_out(ERROR) << "collection: " << colName << endl; - streamlog_out(ERROR) << "you seem to be trying to apply ECAL silicon digitisation to scintillator? Refusing!" - << endl; - streamlog_out(ERROR) << "check setting of ECAL_apply_realistic_digi: " << _applyEcalDigi << endl; + if (m_applyEcalDigi == 1 && thislayersetup.first != SIECAL) { + error() << "collection: " << colName << endmsg; + error() << "You seem to be trying to apply ECAL silicon digitisation to scintillator? Refusing!" << endmsg; + error() << "Check setting of ECAL_apply_realistic_digi: " << m_applyEcalDigi << endmsg; assert(0); - _countWarnings++; + m_countWarnings++; } - if (_applyEcalDigi == 2 && thislayersetup.first != SCECAL) { - streamlog_out(ERROR) << "collection: " << colName << endl; - streamlog_out(ERROR) << "you seem to be trying to apply ECAL scintillator digitisation to silicon? Refusing!" - << endl; - streamlog_out(ERROR) << "check setting of ECAL_apply_realistic_digi: " << _applyEcalDigi << endl; + if (m_applyEcalDigi == 2 && thislayersetup.first != SCECAL) { + error() << "collection: " << colName << endmsg; + error() << "You seem to be trying to apply ECAL scintillator digitisation to silicon? Refusing!" << endmsg; + error() << "Check setting of ECAL_apply_realistic_digi: " << m_applyEcalDigi << endmsg; assert(0); - _countWarnings++; + m_countWarnings++; } if (thislayersetup.second != getStripOrientationFromColName(colName)) { - streamlog_out(ERROR) << "collection: " << colName << endl; - streamlog_out(ERROR) << "some inconsistency in strip orientation?" << endl; - streamlog_out(ERROR) << " from collection name: " << getStripOrientationFromColName(colName) << endl; - streamlog_out(ERROR) << " from layer config string: " << thislayersetup.second << endl; - _countWarnings++; + error() << "collection: " << colName << endmsg; + error() << "Some inconsistency in strip orientation?" << endmsg; + error() << " from collection name: " << getStripOrientationFromColName(colName) << endmsg; + error() << " from layer config string: " << thislayersetup.second << endmsg; + m_countWarnings++; } return; } -std::pair DDCaloDigi::getLayerProperties(std::string colName, int layer) { +std::pair DDCaloDigi::getLayerProperties(std::string const& colName, int layer) const { std::pair thislayersetup(-99, -99); std::string colNameLow(colName); std::transform(colNameLow.begin(), colNameLow.end(), colNameLow.begin(), ::tolower); - if (colNameLow.find("presh") != string::npos) { // preshower + if (colNameLow.find("presh") != std::string::npos) { // preshower if (layer != 0) { - streamlog_out(WARNING) << "preshower layer with layer index = " << layer << " ??? " << endl; + warning() << "preshower layer with layer index = " << layer << " ??? " << endmsg; } else { thislayersetup = getLayerConfig()[layer]; } } else if (colNameLow.find("ring") != - string::npos) { // ecal ring (has no preshower), and is actually always all silicon + std::string::npos) { // ecal ring (has no preshower), and is actually always all silicon if (layer < int(getLayerConfig().size())) { // thislayersetup = getLayerConfig()[layer]; thislayersetup = std::pair(SIECAL, SQUARE); } else { - streamlog_out(WARNING) << "unphysical layer number? " << layer << " " << getLayerConfig().size() << endl; + warning() << "unphysical layer number? " << layer << " " << getLayerConfig().size() << endmsg; } } else { // endcap, barrel if (layer + 1 < int(getLayerConfig().size())) { thislayersetup = getLayerConfig()[layer + 1]; } else { - streamlog_out(WARNING) << "unphysical layer number? " << layer << " " << getLayerConfig().size() << endl; + warning() << "unphysical layer number? " << layer << " " << getLayerConfig().size() << endmsg; } } return thislayersetup; } -int DDCaloDigi::getStripOrientationFromColName(std::string colName) { +int DDCaloDigi::getStripOrientationFromColName(std::string const& colName) const { int orientation(-99); std::string colNameLow(colName); std::transform(colNameLow.begin(), colNameLow.end(), colNameLow.begin(), ::tolower); - if (colNameLow.find("trans") != string::npos) { + if (colNameLow.find("trans") != std::string::npos) { orientation = STRIP_ALIGN_ACROSS_SLAB; - } else if (colNameLow.find("long") != string::npos) { + } else if (colNameLow.find("long") != std::string::npos) { orientation = STRIP_ALIGN_ALONG_SLAB; } else { // assume square... orientation = SQUARE; - // cout << "WARNING, cannot guess strip orientation! for collection " << colName << endl; + //std::cout << "WARNING, cannot guess strip orientation! for collection " << colName << std::endl; } return orientation; -} +} \ No newline at end of file diff --git a/k4GaudiPandora/src/DDScintillatorPpdDigi.cc b/k4GaudiPandora/src/DDScintillatorPpdDigi.cc index f2b778b..9fea288 100644 --- a/k4GaudiPandora/src/DDScintillatorPpdDigi.cc +++ b/k4GaudiPandora/src/DDScintillatorPpdDigi.cc @@ -36,13 +36,13 @@ DDScintillatorPpdDigi::DDScintillatorPpdDigi() {} void DDScintillatorPpdDigi::printParameters() { cout << "--------------------------------" << endl; cout << "DDScintillatorPpdDigi parameters" << endl; - cout << " pe_per_mip = " << _pe_per_mip << endl; - cout << " calib_mip = " << _calib_mip << endl; - cout << " npix = " << _npix << endl; - cout << " misCalibNpix = " << _misCalibNpix << endl; - cout << " pixSpread = " << _pixSpread << endl; - cout << " elecDynRange = " << _elecMaxDynRange_MIP << endl; - cout << " elecNoise = " << _elecNoise << endl; + cout << " PEperMIP = " << m_PEperMIP << endl; + cout << " calibMIP = " << m_calibMIP << endl; + cout << " Npix = " << m_Npix << endl; + cout << " misCalibNpix = " << m_misCalibNpix << endl; + cout << " pixSpread = " << m_pixSpread << endl; + cout << " elecDynRange = " << m_elecMaxDynRange_MIP << endl; + cout << " elecNoise = " << m_elecNoise << endl; cout << "--------------------------------" << endl; return; } @@ -50,75 +50,74 @@ void DDScintillatorPpdDigi::printParameters() { float DDScintillatorPpdDigi::getDigitisedEnergy(float energy) { float correctedEnergy(energy); - if (_pe_per_mip <= 0 || _calib_mip <= 0 || _npix <= 0) { - cout << "ERROR, crazy parameters for DDScintillatorPpdDigi: PE/MIP=" << _pe_per_mip << ", MIP calib=" << _calib_mip - << ", #pixels=" << _npix << endl; - cout << "you must specify at least the #PE/MIP, MIP calibration, and #pixels for realistic scintillator " + if (m_PEperMIP <= 0 || m_calibMIP <= 0 || m_Npix <= 0) { + cout << "ERROR, crazy parameters for DDScintillatorPpdDigi: PE/MIP = " << m_PEperMIP + << ", MIP calibration = " << m_calibMIP << ", # pixels = " << m_Npix << endl; + cout << "You must specify at least the #PE/MIP, MIP calibration, and #pixels for realistic scintillator " "digitisation!!" << endl; - cout << "refusing to proceed!" << endl; + cout << "Refusing to proceed!" << endl; assert(0); - } - - // 1. convert energy to expected # photoelectrons (via MIPs) - float npe = _pe_per_mip * energy / _calib_mip; + } else { + // 1. convert energy to expected # photoelectrons (via MIPs) + float Npe = m_PEperMIP * energy / m_calibMIP; - //oh: commented out Daniel's digitisation model. (npe -> poisson -> saturation -> stoykov smearing). - // Stoykov smearing used with Gaussian shape for lack of better model. - /* - // 2. smear according to poisson (PE statistics) - npe = CLHEP::RandPoisson::shoot( npe ); + //oh: commented out Daniel's digitisation model. (npe -> poisson -> saturation -> stoykov smearing). + // Stoykov smearing used with Gaussian shape for lack of better model. + /* + // 2. smear according to poisson (PE statistics) + npe = CLHEP::RandPoisson::shoot( npe ); - if (_npix>0) { - // 3. apply average sipm saturation behaviour - npe = _npix*(1.0 - exp( -npe/_npix ) ); + if (_npix>0) { + // 3. apply average sipm saturation behaviour + npe = _npix*(1.0 - exp( -npe/_npix ) ); - // 4. apply intrinsic sipm fluctuations (see e.g arXiv:0706.0746 [physics.ins-det]) - float alpha = npe/_npix; // frac of hit pixels - float width = sqrt( _npix*exp(-alpha)*( 1. - (1.+alpha)*exp(-alpha) ) ); + // 4. apply intrinsic sipm fluctuations (see e.g arXiv:0706.0746 [physics.ins-det]) + float alpha = npe/_npix; // frac of hit pixels + float width = sqrt( _npix*exp(-alpha)*( 1. - (1.+alpha)*exp(-alpha) ) ); - // make an integer correction - int dpix = int( floor( CLHEP::RandGauss::shoot(0, width) + 0.5 ) ); + // make an integer correction + int dpix = int( floor( CLHEP::RandGauss::shoot(0, width) + 0.5 ) ); - npe += dpix; - } -*/ + npe += dpix; + } + */ - //AHCAL TB style digitisation: npe -> saturation -> binomial smearing - //shown to be mathematically equivalent to Daniel's model above, but slightly faster and automatically generates correct shape instead of Gaussian approximation + //AHCAL TB style digitisation: npe -> saturation -> binomial smearing + //shown to be mathematically equivalent to Daniel's model above, but slightly faster and automatically generates correct shape instead of Gaussian approximation - if (_npix > 0) { - // apply average sipm saturation behaviour - npe = _npix * (1.0 - exp(-npe / _npix)); + if (m_Npix > 0) { + // apply average SiPM saturation behaviour + Npe = m_Npix * (1.0 - exp(-Npe / m_Npix)); - //apply binomial smearing - float p = npe / _npix; // fraction of hit pixels on SiPM - npe = CLHEP::RandBinomial::shoot(_npix, p); //npe now quantised to integer pixels - } + //apply binomial smearing + float p = Npe / m_Npix; // fraction of hit pixels on SiPM + Npe = CLHEP::RandBinomial::shoot(m_Npix, p); // # photoelectrons now quantised to integer pixels + } - if (_pixSpread > 0) { - // variations in pixel capacitance - npe *= CLHEP::RandGauss::shoot(1, _pixSpread / sqrt(npe)); - } + if (m_pixSpread > 0) { + // variations in pixel capacitance + Npe *= CLHEP::RandGauss::shoot(1., m_pixSpread / sqrt(Npe)); + } - if (_elecMaxDynRange_MIP > 0) { - // limited dynamic range of readout electronics - // Daniel moved this here, before the unfolding of saturation (September 2015) - npe = std::min(npe, _elecMaxDynRange_MIP * _pe_per_mip); - } + if (m_elecMaxDynRange_MIP > 0) { + // limited dynamic range of readout electronics + // Daniel moved this here, before the unfolding of saturation (September 2015) + Npe = std::min(Npe, m_elecMaxDynRange_MIP * m_PEperMIP); + } - if (_elecNoise > 0) { - // add electronics noise - npe += CLHEP::RandGauss::shoot(0, _elecNoise * _pe_per_mip); - } + if (m_elecNoise > 0) { + // add electronics noise + Npe += CLHEP::RandGauss::shoot(0., m_elecNoise * m_PEperMIP); + } - if (_npix > 0) { - // 4. unfold the saturation - // - miscalibration of npix - float smearedNpix = _misCalibNpix > 0 ? _npix * CLHEP::RandGauss::shoot(1.0, _misCalibNpix) : _npix; + if (m_Npix > 0) { + // 4. unfold the saturation + // - miscalibration of Npix + float smearedNpix = m_misCalibNpix > 0 ? m_Npix * CLHEP::RandGauss::shoot(1.0, m_misCalibNpix) : m_Npix; - //oh: commented out daniel's implmentation of dealing with hits>smearedNpix. using linearisation of saturation-reconstruction for high amplitude hits instead. - /* + //oh: commented out daniel's implmentation of dealing with hits>smearedNpix. using linearisation of saturation-reconstruction for high amplitude hits instead. + /* // - this is to deal with case when #pe is larger than #pixels (would mean taking log of negative number) float epsilon=1; // any small number... if ( npe>smearedNpix-epsilon ) npe=smearedNpix-epsilon; @@ -126,18 +125,18 @@ float DDScintillatorPpdDigi::getDigitisedEnergy(float energy) { npe = -smearedNpix * std::log ( 1. - ( npe / smearedNpix ) ); */ - const float r = - 0.95; //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; //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. - if (npe < r * smearedNpix) { //current hit below linearisation threshold, reconstruct energy normally: - npe = -smearedNpix * std::log(1. - (npe / smearedNpix)); - } else { //current hit is aove linearisation threshold, reconstruct using linear continuation function: - npe = 1 / (1 - r) * (npe - r * smearedNpix) - smearedNpix * std::log(1 - r); + if (Npe < r * smearedNpix) { //current hit below linearisation threshold, reconstruct energy normally: + Npe = -smearedNpix * std::log(1. - (Npe / smearedNpix)); + } else { //current hit is aove linearisation threshold, reconstruct using linear continuation function: + Npe = 1 / (1 - r) * (Npe - r * smearedNpix) - smearedNpix * std::log(1 - r); + } } - } - - // convert back to energy - correctedEnergy = _calib_mip * npe / _pe_per_mip; + // convert back to energy + correctedEnergy = m_calibMIP * Npe / m_PEperMIP; + } return correctedEnergy; }