From af83315a749bf118e9b499b75cbb3904f6d07ef6 Mon Sep 17 00:00:00 2001 From: Zhibo Wu <91186438+zwu0922@users.noreply.github.com> Date: Fri, 13 Sep 2024 13:39:49 +0200 Subject: [PATCH] Add tools for the implmentation of noise for ALLEGRO v3 (#107) * Add tools for the implmentation of noise for ALLEGRO v3 * Steering file for addNoise * Conflict: remove ALLEGRO tests * Temporarily add the test for ALLEGRO addNoise (to be moved to ddsim+k4run) * Try to get the new test working --------- Co-authored-by: Zhibo Wu Co-authored-by: Zhibo Wu --- .../src/components/CreateCaloCells.h | 8 +- .../TubeLayerModuleThetaMergedCaloTool.cpp | 99 +++ .../TubeLayerModuleThetaMergedCaloTool.h | 59 ++ RecFCCeeCalorimeter/CMakeLists.txt | 7 + ...unEcalBarrel_addNoise_thetamodulemerged.py | 593 ++++++++++++++++++ 5 files changed, 762 insertions(+), 4 deletions(-) create mode 100644 RecCalorimeter/src/components/TubeLayerModuleThetaMergedCaloTool.cpp create mode 100644 RecCalorimeter/src/components/TubeLayerModuleThetaMergedCaloTool.h create mode 100644 RecFCCeeCalorimeter/tests/options/runEcalBarrel_addNoise_thetamodulemerged.py diff --git a/RecCalorimeter/src/components/CreateCaloCells.h b/RecCalorimeter/src/components/CreateCaloCells.h index 70e4b0ce..c4c1de05 100644 --- a/RecCalorimeter/src/components/CreateCaloCells.h +++ b/RecCalorimeter/src/components/CreateCaloCells.h @@ -28,7 +28,7 @@ class IGeoSvc; /** @class CreateCaloCells * * Algorithm for creating calorimeter cells from Geant4 hits. - * Tube geometry with PhiEta segmentation expected. + * Tube geometry with ModuleThetaMerged segmentation expected. * * Flow of the program: * 1/ Merge Geant4 energy deposits with same cellID @@ -70,7 +70,7 @@ class CreateCaloCells : public Gaudi::Algorithm { /// Handle for the calorimeter cells noise tool mutable ToolHandle m_noiseTool{"NoiseCaloCellsFlatTool", this}; /// Handle for the geometry tool - ToolHandle m_geoTool{"TubeLayerPhiEtaCaloTool", this}; + ToolHandle m_geoTool{"TubeLayerModuleThetaMergedCaloTool", this}; /// Add crosstalk to cells? Gaudi::Property m_addCrosstalk{this, "addCrosstalk", false, "Add crosstalk effect?"}; @@ -92,7 +92,7 @@ class CreateCaloCells : public Gaudi::Algorithm { mutable DataHandle m_cells{"cells", Gaudi::DataHandle::Writer, this}; MetaDataHandle m_cellsCellIDEncoding{m_cells, edm4hep::labels::CellIDEncoding, Gaudi::DataHandle::Writer}; /// Name of the detector readout - Gaudi::Property m_readoutName{this, "readoutName", "ECalBarrelPhiEta", "Name of the detector readout"}; + Gaudi::Property m_readoutName{this, "readoutName", "ECalBarrelModuleThetaMerged", "Name of the detector readout"}; /// Name of active volumes Gaudi::Property m_activeVolumeName{this, "activeVolumeName", "_sensitive", "Name of the active volumes"}; /// Name of active layers for sampling calorimeter @@ -116,7 +116,7 @@ class CreateCaloCells : public Gaudi::Algorithm { * This property won't be needed anymore. */ unsigned int m_activeVolumesNumber; - /// Use only volume ID? If false, using PhiEtaSegmentation + /// Use only volume ID? If false, using ModuleThetaMergedSegmentation bool m_useVolumeIdOnly; /// Pointer to the geometry service diff --git a/RecCalorimeter/src/components/TubeLayerModuleThetaMergedCaloTool.cpp b/RecCalorimeter/src/components/TubeLayerModuleThetaMergedCaloTool.cpp new file mode 100644 index 00000000..17bef729 --- /dev/null +++ b/RecCalorimeter/src/components/TubeLayerModuleThetaMergedCaloTool.cpp @@ -0,0 +1,99 @@ +#include "TubeLayerModuleThetaMergedCaloTool.h" + +// segm +#include "DD4hep/Detector.h" +#include "DD4hep/MultiSegmentation.h" +#include "detectorCommon/DetUtils_k4geo.h" +#include "k4Interface/IGeoSvc.h" + +DECLARE_COMPONENT(TubeLayerModuleThetaMergedCaloTool) + +TubeLayerModuleThetaMergedCaloTool::TubeLayerModuleThetaMergedCaloTool(const std::string& type, const std::string& name, + const IInterface* parent) + : AlgTool(type, name, parent), m_geoSvc("GeoSvc", name) { + declareInterface(this); +} + +StatusCode TubeLayerModuleThetaMergedCaloTool::initialize() { + StatusCode sc = AlgTool::initialize(); + if (sc.isFailure()) return sc; + + if (!m_geoSvc) { + error() << "Unable to locate Geometry Service. " + << "Make sure you have GeoSvc and SimSvc in the right order in the configuration." << endmsg; + return StatusCode::FAILURE; + } + if (m_readoutName != "") { + // Check if readouts exist + info() << "Readout: " << m_readoutName << endmsg; + if (m_geoSvc->getDetector()->readouts().find(m_readoutName) == m_geoSvc->getDetector()->readouts().end()) { + error() << "Readout <<" << m_readoutName << ">> does not exist." << endmsg; + return StatusCode::FAILURE; + } + } + return sc; +} + +StatusCode TubeLayerModuleThetaMergedCaloTool::finalize() { return AlgTool::finalize(); } + +StatusCode TubeLayerModuleThetaMergedCaloTool::prepareEmptyCells(std::unordered_map& aCells) { + // Get the total number of active volumes in the geometry + auto highestVol = gGeoManager->GetTopVolume(); + unsigned int numLayers; + if (!m_activeVolumesNumber) { + numLayers = det::utils::countPlacedVolumes(highestVol, m_activeVolumeName); + } else { + // used when MergeLayers tool is used. To be removed once MergeLayer gets replaced by RedoSegmentation. + numLayers = m_activeVolumesNumber; + } + info() << "Number of active layers " << numLayers << endmsg; + + // get segmentation + dd4hep::DDSegmentation::Segmentation *aSegmentation = m_geoSvc->getDetector()->readout(m_readoutName).segmentation().segmentation(); + std::string segmentationType = aSegmentation->type(); + info() << "Segmentation type : " << segmentationType << endmsg; + + dd4hep::DDSegmentation::FCCSWGridModuleThetaMerged_k4geo *moduleThetaSegmentation = nullptr; + if (segmentationType == "FCCSWGridModuleThetaMerged_k4geo") { + moduleThetaSegmentation = dynamic_cast(aSegmentation); + info() << "Segmentation: bins in Module " << moduleThetaSegmentation->nModules() << endmsg; + } + else { + error() << "Unable to cast segmentation pointer!!!! Tool only applicable to FCCSWGridModuleThetaMerged_k4geo segmentation." << endmsg; + return StatusCode::FAILURE; + } + + // get decoder and extrema + auto decoder = m_geoSvc->getDetector()->readout(m_readoutName).idSpec().decoder(); + std::vector> extrema; + extrema.push_back(std::make_pair(0, numLayers - 1)); + extrema.push_back(std::make_pair(0, 0)); + extrema.push_back(std::make_pair(0, 0)); + + for (unsigned int ilayer = 0; ilayer < numLayers; ilayer++) { + dd4hep::DDSegmentation::CellID volumeId = 0; + (*decoder)["system"].set(volumeId, 4); + (*decoder)["layer"].set(volumeId, ilayer); + (*decoder)["theta"].set(volumeId, 0); + (*decoder)["module"].set(volumeId, 0); + // Get number of segmentation cells within the active volume, for given layer + auto numCells = det::utils::numberOfCells(volumeId, *moduleThetaSegmentation); + // extrema[1]: Range of module ID (0, max module ID) + extrema[1] = std::make_pair(0, (numCells[0] - 1) * moduleThetaSegmentation->mergedModules(ilayer)); + // extrema[2]: Range of theta ID (min theta ID, max theta ID + extrema[2] = std::make_pair(numCells[2], numCells[2] + (numCells[1] - 1) * moduleThetaSegmentation->mergedThetaCells(ilayer)); + debug() << "Layer: " << ilayer << endmsg; + debug() << "Number of segmentation cells in (module, theta): " << numCells << endmsg; + // Loop over segmentation cells + for (unsigned int imodule = 0; imodule < numCells[0]; imodule++) { + for (unsigned int itheta = 0; itheta < numCells[1]; itheta++) { + dd4hep::DDSegmentation::CellID cellId = volumeId; + decoder->set(cellId, "module", imodule * moduleThetaSegmentation->mergedModules(ilayer)); + decoder->set(cellId, "theta", numCells[2] + itheta * moduleThetaSegmentation->mergedThetaCells(ilayer)); // start from the minimum existing theta cell in this layer + aCells.insert(std::pair(cellId, 0)); + } + } + } + + return StatusCode::SUCCESS; +} diff --git a/RecCalorimeter/src/components/TubeLayerModuleThetaMergedCaloTool.h b/RecCalorimeter/src/components/TubeLayerModuleThetaMergedCaloTool.h new file mode 100644 index 00000000..b414da0e --- /dev/null +++ b/RecCalorimeter/src/components/TubeLayerModuleThetaMergedCaloTool.h @@ -0,0 +1,59 @@ +#ifndef RECCALORIMETER_TUBELAYERMODULETHETAMERGEDCALOTOOL_H +#define RECCALORIMETER_TUBELAYERMODULETHETAMERGEDCALOTOOL_H + +// from Gaudi +#include "GaudiKernel/AlgTool.h" + +// k4FWCore +#include "k4Interface/ICalorimeterTool.h" + +class IGeoSvc; + +/** @class TubeLayerModuleThetaMergedCaloTool Reconstruction/RecCalorimeter/src/components/TubeLayerModuleThetaMergedCaloTool.h + *TubeLayerModuleThetaMergedCaloTool.h + * + * Tool for geometry-dependent settings of the digitisation. + * It assumes cylindrical geometry (layers) and phi-theta segmentation. + * + * @author Anna Zaborowska + * @author Zhibo Wu + */ + +class TubeLayerModuleThetaMergedCaloTool : public AlgTool, virtual public ICalorimeterTool { +public: + TubeLayerModuleThetaMergedCaloTool(const std::string& type, const std::string& name, const IInterface* parent); + virtual ~TubeLayerModuleThetaMergedCaloTool() = default; + virtual StatusCode initialize() final; + virtual StatusCode finalize() final; + /** Prepare a map of all existing cells in current geometry. + * Active layers (cylindrical tubes) are looked in the geometry manager by name ('\b activeVolumeName'). + * Corresponding bitfield name is given in '\b activeFieldName'. + * If users wants to limit the number of active layers, it is possible by setting '\b activeVolumesNumber'. + * The total number of cells N = n_layer * n_theta * n_phi, where + * n_layer is number of layers (taken from geometry if activeVolumesNumber not set), + * n_theta is number of eta bins in that layer, + * n_phi is number of phi bins (the same for each layer). + * For more explanation please [see reconstruction documentation](@ref md_reconstruction_doc_reccalorimeter). + * @param[out] aCells map of existing cells (and deposited energy, set to 0) + * return Status code. + */ + virtual StatusCode prepareEmptyCells(std::unordered_map& aCells) final; + +private: + /// Pointer to the geometry service + ServiceHandle m_geoSvc; + /// Name of the detector readout + Gaudi::Property m_readoutName{this, "readoutName", ""}; + /// Name of active volumes + Gaudi::Property m_activeVolumeName{this, "activeVolumeName", "LAr_sensitive"}; + /// Name of active layers for sampling calorimeter + Gaudi::Property m_activeFieldName{this, "activeFieldName", "active_layer"}; + /// Name of the fields describing the segmented volume + Gaudi::Property> m_fieldNames{this, "fieldNames"}; + /// Values of the fields describing the segmented volume + Gaudi::Property> m_fieldValues{this, "fieldValues"}; + /// Temporary: for use with MergeLayer tool + Gaudi::Property m_activeVolumesNumber{this, "activeVolumesNumber", 0}; +}; + +#endif /* RECCALORIMETER_TUBELAYERMODULETHETAMERGEDCALOTOOL_H */ diff --git a/RecFCCeeCalorimeter/CMakeLists.txt b/RecFCCeeCalorimeter/CMakeLists.txt index 69686b19..3f2c04e5 100644 --- a/RecFCCeeCalorimeter/CMakeLists.txt +++ b/RecFCCeeCalorimeter/CMakeLists.txt @@ -50,3 +50,10 @@ add_test(NAME FCCeeLAr_benchmarkCorrection WORKING_DIRECTORY ${PROJECT_SOURCE_DIR} ) set_test_env(FCCeeLAr_benchmarkCorrection) + +# test for ALLEGRO addNoise, to be moved ddsim+k4run +add_test(NAME FCCeeLAr_ecalNoise + COMMAND k4run RecFCCeeCalorimeter/tests/options/runEcalBarrel_addNoise_thetamodulemerged.py + WORKING_DIRECTORY ${PROJECT_SOURCE_DIR} +) +set_test_env(FCCeeLAr_ecalNoise) diff --git a/RecFCCeeCalorimeter/tests/options/runEcalBarrel_addNoise_thetamodulemerged.py b/RecFCCeeCalorimeter/tests/options/runEcalBarrel_addNoise_thetamodulemerged.py new file mode 100644 index 00000000..4f0ba9bd --- /dev/null +++ b/RecFCCeeCalorimeter/tests/options/runEcalBarrel_addNoise_thetamodulemerged.py @@ -0,0 +1,593 @@ +# +# IMPORTS +# +from Configurables import ApplicationMgr +from Configurables import AuditorSvc, ChronoAuditor +# Input/output +from Configurables import k4DataSvc, PodioInput +from Configurables import PodioOutput +# Geometry +from Configurables import GeoSvc +# Create cells +from Configurables import CreateCaloCells +from Configurables import CreateEmptyCaloCellsCollection +# Cell positioning tools +from Configurables import CreateCaloCellPositionsFCCee +from Configurables import CellPositionsECalBarrelModuleThetaSegTool +from Configurables import CellPositionsECalEndcapTurbineSegTool +# Redo segmentation for ECAL and HCAL +from Configurables import RedoSegmentation +# Read noise values from file and generate noise in cells +from Configurables import NoiseCaloCellsVsThetaFromFileTool +# Apply sampling fraction corrections +from Configurables import CalibrateCaloHitsTool +from Configurables import CalibrateInLayersTool +# Up/down stream correction +from Configurables import CorrectCaloClusters +# SW clustering +from Configurables import CaloTowerToolFCCee +from Configurables import CreateCaloClustersSlidingWindowFCCee +# Topo clustering +from Configurables import CaloTopoClusterInputTool +from Configurables import TopoCaloNeighbours +from Configurables import TopoCaloNoisyCells +from Configurables import CaloTopoClusterFCCee +# Decorate clusters with shower shape parameters +from Configurables import AugmentClustersFCCee +# Logger +from Gaudi.Configuration import INFO, VERBOSE, DEBUG +# units and physical constants +from GaudiKernel.SystemOfUnits import GeV, tesla, mm +from GaudiKernel.PhysicalConstants import pi, halfpi, twopi +# python libraries +import os +from math import cos, sin, tan + +# +# SETTINGS +# + +# - general settings +# +inputfile = "https://fccsw.web.cern.ch/fccsw/filesForSimDigiReco/ALLEGRO/ALLEGRO_o1_v03/ALLEGRO_sim.root" # input file produced with ddsim +Nevts = 2 # -1 means all events +addNoise = True # add noise or not to the cell energy +dumpGDML = False # create GDML file of detector model + +# - what to save in output file +# +# for big productions, save significant space removing hits and cells +# however, hits and cluster cells might be wanted for small productions for detailed event displays +# cluster cells are not needed for the training of the MVA energy regression nor the photon ID since needed quantities are stored in cluster shapeParameters +# saveHits = False +# saveCells = False +# saveClusterCells = False +saveHits = True +saveCells = True +saveClusterCells = True + +# ECAL barrel parameters for digitisation +samplingFraction=[0.3800493723322256] * 1 + [0.13494147915064658] * 1 + [0.142866851721152] * 1 + [0.14839315921940666] * 1 + [0.15298362570665006] * 1 + [0.15709704561942747] * 1 + [0.16063717490147533] * 1 + [0.1641723795419055] * 1 + [0.16845490287689746] * 1 + [0.17111520115997653] * 1 + [0.1730605163148862] * 1 +upstreamParameters = [[0.028158491043365624, -1.564259408365951, -76.52312805346982, 0.7442903558010191, -34.894692961350195, -74.19340877431723]] +downstreamParameters = [[0.00010587711361028165, 0.0052371999097777355, 0.69906696456064, -0.9348243433360095, -0.0364714212117143, 8.360401126995626]] + +ecalBarrelLayers = len(samplingFraction) +resegmentECalBarrel = False + +# - parameters for clustering +# +doSWClustering = True +doTopoClustering = True + +# cluster energy corrections +# simple parametrisations of up/downstream losses +# not to be applied for ECAL+HCAL clustering +applyUpDownstreamCorrections = False + +# calculate cluster energy and barycenter per layer and save it as extra parameters +addShapeParameters = True +ecalBarrelThetaWeights = [-1, 3.0, 3.0, 3.0, 4.25, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0] # to be recalculated for V03, separately for topo and calo clusters... + +# +# ALGORITHMS AND SERVICES SETUP +# + +# Input: load the output of the SIM step +podioevent = k4DataSvc('EventDataSvc') +podioevent.input = inputfile +input_reader = PodioInput('InputReader') + + +# Detector geometry +# prefix all xmls with path_to_detector +# if K4GEO is empty, this should use relative path to working directory +geoservice = GeoSvc("GeoSvc") +path_to_detector = os.environ.get("K4GEO", "") +detectors_to_use = [ + 'FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ALLEGRO_o1_v03.xml' +] +geoservice.detectors = [ + os.path.join(path_to_detector, _det) for _det in detectors_to_use +] +geoservice.OutputLevel = INFO + +# GDML dump of detector model +if dumpGDML: + from Configurables import GeoToGdmlDumpSvc + gdmldumpservice = GeoToGdmlDumpSvc("GeoToGdmlDumpSvc") + +# Digitisation (merging hits into cells, EM scale calibration via sampling fractions) + +# - ECAL readouts +ecalBarrelReadoutName = "ECalBarrelModuleThetaMerged" # barrel, original segmentation (baseline) +ecalBarrelReadoutName2 = "ECalBarrelModuleThetaMerged2" # barrel, after re-segmentation (for optimisation studies) +ecalEndcapReadoutName = "ECalEndcapTurbine" # endcap, turbine-like (baseline) +# - HCAL readouts +hcalBarrelReadoutName = "" +hcalBarrelReadoutName2 = "" +hcalEndcapReadoutName = "" + +# - EM scale calibration (sampling fraction) +# * ECAL barrel +calibEcalBarrel = CalibrateInLayersTool("CalibrateECalBarrel", + samplingFraction=samplingFraction, + readoutName=ecalBarrelReadoutName, + layerFieldName="layer") +# * ECAL endcap +calibEcalEndcap = CalibrateInLayersTool("CalibrateECalEndcap", + samplingFraction = [0.16419] * 1 + [0.192898] * 1 + [0.18783] * 1 + [0.193203] * 1 + [0.193928] * 1 + [0.192286] * 1 + [0.199959] * 1 + [0.200153] * 1 + [0.212635] * 1 + [0.180345] * 1 + [0.18488] * 1 + [0.194762] * 1 + [0.197775] * 1 + [0.200504] * 1 + [0.205555] * 1 + [0.203601] * 1 + [0.210877] * 1 + [0.208376] * 1 + [0.216345] * 1 + [0.201452] * 1 + [0.202134] * 1 + [0.207566] * 1 + [0.208152] * 1 + [0.209889] * 1 + [0.211743] * 1 + [0.213188] * 1 + [0.215864] * 1 + [0.22972] * 1 + [0.192515] * 1 + [0.0103233] * 1, + readoutName=ecalEndcapReadoutName, + layerFieldName="layer") + +# Create cells in ECal barrel (needed if one wants to apply cell calibration, +# which is not performed by ddsim) +# - merge hits into cells according to initial segmentation +ecalBarrelCellsName = "ECalBarrelCells" +if not addNoise: # Create ECAL barrel cells without adding noise + createEcalBarrelCells = CreateCaloCells("CreateECalBarrelCells", + doCellCalibration=True, + calibTool=calibEcalBarrel, + addCellNoise=False, + filterCellNoise=False, + addPosition=True, + OutputLevel=INFO, + hits=ecalBarrelReadoutName, + cells=ecalBarrelCellsName) + +# - add to Ecal barrel cells the position information +# (good for physics, all coordinates set properly) +cellPositionEcalBarrelTool = CellPositionsECalBarrelModuleThetaSegTool( + "CellPositionsECalBarrel", + readoutName=ecalBarrelReadoutName, + OutputLevel=INFO +) +ecalBarrelPositionedCellsName = ecalBarrelReadoutName + "Positioned" +createEcalBarrelPositionedCells = CreateCaloCellPositionsFCCee( + "CreateECalBarrelPositionedCells", + OutputLevel=INFO +) +createEcalBarrelPositionedCells.positionsTool = cellPositionEcalBarrelTool +createEcalBarrelPositionedCells.hits.Path = ecalBarrelCellsName +createEcalBarrelPositionedCells.positionedHits.Path = ecalBarrelPositionedCellsName + +# - now, if we want to also save cells with coarser granularity: +if resegmentECalBarrel: + # 2. step - rewrite the cellId using the merged theta-module segmentation + # (merging several modules and severla theta readout cells). + # Add noise at this step if you derived the noise already assuming merged cells + # Step 2a: compute new cellID of cells based on new readout + # (merged module-theta segmentation with variable merging vs layer) + resegmentEcalBarrelTool = RedoSegmentation("ReSegmentationEcal", + # old bitfield (readout) + oldReadoutName=ecalBarrelReadoutName, + # specify which fields are going to be altered (deleted/rewritten) + oldSegmentationIds=["module", "theta"], + # new bitfield (readout), with new segmentation (merged modules and theta cells) + newReadoutName=ecalBarrelReadoutName2, + OutputLevel=INFO, + debugPrint=200, + inhits=ecalBarrelCellsName, + outhits="ECalBarrelCellsMerged") + + # Step 2b: merge new cells with same cellID together + # do not apply cell calibration again since cells were already + # calibrated in Step 1 + ecalBarrelCellsName2 = "ECalBarrelCells2" + createEcalBarrelCells2 = CreateCaloCells("CreateECalBarrelCells2", + doCellCalibration=False, + addCellNoise=False, + filterCellNoise=False, + OutputLevel=INFO, + hits="ECalBarrelCellsMerged", + cells=ecalBarrelCellsName2) + + + cellPositionEcalBarrelTool2 = CellPositionsECalBarrelModuleThetaSegTool( + "CellPositionsECalBarrel2", + readoutName=ecalBarrelReadoutName2, + OutputLevel=INFO + ) + createEcalBarrelPositionedCells2 = CreateCaloCellPositionsFCCee( + "CreateECalBarrelPositionedCells2", + OutputLevel=INFO + ) + createEcalBarrelPositionedCells2.positionsTool = cellPositionEcalBarrelTool2 + createEcalBarrelPositionedCells2.hits.Path = ecalBarrelCellsName2 + createEcalBarrelPositionedCells2.positionedHits.Path = ecalBarrelReadoutName2 + "Positioned" + + +# Create cells in ECal endcap (needed if one wants to apply cell calibration, +# which is not performed by ddsim) +ecalEndcapCellsName = "ECalEndcapCells" +createEcalEndcapCells = CreateCaloCells("CreateEcalEndcapCaloCells", + doCellCalibration=True, + calibTool=calibEcalEndcap, + addCellNoise=False, + filterCellNoise=False, + OutputLevel=INFO, + hits=ecalEndcapReadoutName, + cells=ecalEndcapCellsName) + +# Add to Ecal endcap cells the position information +# (good for physics, all coordinates set properly) +cellPositionEcalEndcapTool = CellPositionsECalEndcapTurbineSegTool( + "CellPositionsECalEndcap", + readoutName=ecalEndcapReadoutName, + OutputLevel=INFO +) +ecalEndcapPositionedCellsName = ecalEndcapReadoutName + "Positioned" +createEcalEndcapPositionedCells = CreateCaloCellPositionsFCCee( + "CreateECalEndcapPositionedCells", + OutputLevel=INFO +) +createEcalEndcapPositionedCells.positionsTool = cellPositionEcalEndcapTool +createEcalEndcapPositionedCells.hits.Path = ecalEndcapCellsName +createEcalEndcapPositionedCells.positionedHits.Path = ecalEndcapPositionedCellsName + + +if addNoise: + # FIXME Input histograms of ECAL noise needs to be updated to match the latest ALLEGRO geometry + ecalBarrelNoisePath = "https://fccsw.web.cern.ch/fccsw/filesForSimDigiReco/ALLEGRO/ALLEGRO_o1_v03/elecNoise_ecalBarrelFCCee_theta.root" + ecalBarrelNoiseRMSHistName = "h_elecNoise_fcc_" + from Configurables import NoiseCaloCellsVsThetaFromFileTool + noiseBarrel = NoiseCaloCellsVsThetaFromFileTool("NoiseBarrel", + cellPositionsTool=cellPositionEcalBarrelTool, + readoutName=ecalBarrelReadoutName, + noiseFileName=ecalBarrelNoisePath, + elecNoiseRMSHistoName=ecalBarrelNoiseRMSHistName, + setNoiseOffset=False, + activeFieldName="layer", + addPileup=False, + filterNoiseThreshold=0, + numRadialLayers=11, + scaleFactor=1 / 1000., # MeV to GeV + OutputLevel=INFO) + + # barrel geometry tool is migrated to match ALLEGRO v3 segmentation + from Configurables import TubeLayerModuleThetaMergedCaloTool + barrelGeometry = TubeLayerModuleThetaMergedCaloTool("EcalBarrelGeo", + readoutName=ecalBarrelReadoutName, + activeVolumeName="LAr_sensitive", + activeFieldName="layer", + activeVolumesNumber=11, + fieldNames=["system"], + fieldValues=[4]) + + # cells with noise not filtered + createEcalBarrelCells = CreateCaloCells("CreateECalBarrelCellsNoise", + doCellCalibration=True, + calibTool=calibEcalBarrel, + addCellNoise=True, + filterCellNoise=False, + noiseTool=noiseBarrel, + addPosition=True, + geometryTool=barrelGeometry, + OutputLevel=INFO, + hits=ecalBarrelReadoutName, + cells=ecalBarrelCellsName) + + # cells with noise filtered + # createEcalBarrelCellsNoise = CreateCaloCells("CreateECalBarrelCellsNoise_filtered", + # doCellCalibration=False, + # addCellNoise=True, + # filterCellNoise=True, + # OutputLevel=INFO, + # hits="ECalBarrelCellsStep2", + # noiseTool=noiseBarrel, + # geometryTool=barrelGeometry, + # cells=EcalBarrelCellsName) + +hcalBarrelCellsName = "emptyCaloCells" +hcalBarrelPositionedCellsName = "emptyCaloCells" +hcalBarrelCellsName2 = "emptyCaloCells" +hcalBarrelPositionedCellsName2 = "emptyCaloCells" +cellPositionHcalBarrelTool = None +cellPositionHcalBarrelTool2 = None + +# Empty cells for parts of calorimeter not implemented yet +createemptycells = CreateEmptyCaloCellsCollection("CreateEmptyCaloCells") +createemptycells.cells.Path = "emptyCaloCells" + +# Produce sliding window clusters +if doSWClustering: + towers = CaloTowerToolFCCee("towers", + deltaThetaTower=4 * 0.009817477 / 4, deltaPhiTower=2 * 2 * pi / 1536., + ecalBarrelReadoutName=ecalBarrelReadoutName, + ecalEndcapReadoutName=ecalEndcapReadoutName, + ecalFwdReadoutName="", + hcalBarrelReadoutName=hcalBarrelReadoutName2, + hcalExtBarrelReadoutName="", + hcalEndcapReadoutName="", + hcalFwdReadoutName="", + OutputLevel=INFO) + towers.ecalBarrelCells.Path = ecalBarrelPositionedCellsName + towers.ecalEndcapCells.Path = "emptyCaloCells" + towers.ecalFwdCells.Path = "emptyCaloCells" + towers.hcalBarrelCells.Path = hcalBarrelPositionedCellsName2 + towers.hcalExtBarrelCells.Path = "emptyCaloCells" + towers.hcalEndcapCells.Path = "emptyCaloCells" + towers.hcalFwdCells.Path = "emptyCaloCells" + + # Cluster variables + windT = 9 + windP = 17 + posT = 5 + posP = 11 + dupT = 7 + dupP = 13 + finT = 9 + finP = 17 + # Minimal energy to create a cluster in GeV (FCC-ee detectors have to reconstruct low energy particles) + threshold = 0.040 + + createClusters = CreateCaloClustersSlidingWindowFCCee("CreateClusters", + towerTool=towers, + nThetaWindow=windT, nPhiWindow=windP, + nThetaPosition=posT, nPhiPosition=posP, + nThetaDuplicates=dupT, nPhiDuplicates=dupP, + nThetaFinal=finT, nPhiFinal=finP, + energyThreshold=threshold, + energySharingCorrection=False, + attachCells=True, + OutputLevel=INFO + ) + createClusters.clusters.Path = "CaloClusters" + createClusters.clusterCells.Path = "CaloClusterCells" + + if applyUpDownstreamCorrections: + correctCaloClusters = CorrectCaloClusters("CorrectCaloClusters", + inClusters=createClusters.clusters.Path, + outClusters="Corrected" + createClusters.clusters.Path, + systemIDs=[4], + numLayers=[ecalBarrelLayers], + firstLayerIDs=[0], + lastLayerIDs=[ecalBarrelLayers-1], + readoutNames=[ecalBarrelReadoutName], + upstreamParameters=upstreamParameters, + upstreamFormulas=[ + ['[0]+[1]/(x-[2])', '[0]+[1]/(x-[2])']], + downstreamParameters=downstreamParameters, + downstreamFormulas=[ + ['[0]+[1]*x', '[0]+[1]/sqrt(x)', '[0]+[1]/x']], + OutputLevel=INFO + ) + + if addShapeParameters: + augmentCaloClusters = AugmentClustersFCCee("augmentCaloClusters", + inClusters=createClusters.clusters.Path, + outClusters="Augmented" + createClusters.clusters.Path, + systemIDs=[4], + systemNames=["EMB"], + numLayers=[ecalBarrelLayers], + readoutNames=[ecalBarrelReadoutName], + layerFieldNames=["layer"], + thetaRecalcWeights=[ecalBarrelThetaWeights], + do_photon_shapeVar=False, + OutputLevel=INFO + ) + +if doTopoClustering: + # Produce topoclusters (ECAL only or ECAL+HCAL) + createTopoInput = CaloTopoClusterInputTool("CreateTopoInput", + ecalBarrelReadoutName=ecalBarrelReadoutName, + ecalEndcapReadoutName="", + ecalFwdReadoutName="", + hcalBarrelReadoutName=hcalBarrelReadoutName2, + hcalExtBarrelReadoutName="", + hcalEndcapReadoutName="", + hcalFwdReadoutName="", + OutputLevel=INFO) + + createTopoInput.ecalBarrelCells.Path = ecalBarrelPositionedCellsName + createTopoInput.ecalEndcapCells.Path = "emptyCaloCells" + createTopoInput.ecalFwdCells.Path = "emptyCaloCells" + createTopoInput.hcalBarrelCells.Path = hcalBarrelPositionedCellsName2 + createTopoInput.hcalExtBarrelCells.Path = "emptyCaloCells" + createTopoInput.hcalEndcapCells.Path = "emptyCaloCells" + createTopoInput.hcalFwdCells.Path = "emptyCaloCells" + cellPositionHcalBarrelNoSegTool = None + cellPositionHcalExtBarrelTool = None + + neighboursMap = "https://fccsw.web.cern.ch/fccsw/filesForSimDigiReco/ALLEGRO/ALLEGRO_o1_v03/neighbours_map_ecalB_thetamodulemerged.root" + noiseMap = "https://fccsw.web.cern.ch/fccsw/filesForSimDigiReco/ALLEGRO/ALLEGRO_o1_v03/cellNoise_map_electronicsNoiseLevel_ecalB_thetamodulemerged.root" + readNeighboursMap = TopoCaloNeighbours("ReadNeighboursMap", + fileName=neighboursMap, + OutputLevel=INFO) + + # Noise levels per cell + readNoisyCellsMap = TopoCaloNoisyCells("ReadNoisyCellsMap", + fileName=noiseMap, + OutputLevel=INFO) + + createTopoClusters = CaloTopoClusterFCCee("CreateTopoClusters", + TopoClusterInput=createTopoInput, + # expects neighbours map from cellid->vec < neighbourIds > + neigboursTool=readNeighboursMap, + # tool to get noise level per cellid + noiseTool=readNoisyCellsMap, + # cell positions tools for all sub - systems + positionsECalBarrelTool=cellPositionEcalBarrelTool, + positionsHCalBarrelTool=cellPositionHcalBarrelTool2, + # positionsHCalBarrelNoSegTool=cellPositionHcalBarrelNoSegTool, + # positionsHCalExtBarrelTool=cellPositionHcalExtBarrelTool, + # positionsHCalExtBarrelTool = HCalExtBcells, + # positionsEMECTool = EMECcells, + # positionsHECTool = HECcells, + # positionsEMFwdTool = ECalFwdcells, + # positionsHFwdTool = HCalFwdcells, + noSegmentationHCal=False, + seedSigma=4, + neighbourSigma=2, + lastNeighbourSigma=0, + OutputLevel=INFO) + createTopoClusters.clusters.Path = "CaloTopoClusters" + createTopoClusters.clusterCells.Path = "CaloTopoClusterCells" + + + # Correction below is for EM-only clusters + # Need something different for EM+HCAL + if applyUpDownstreamCorrections: + correctCaloTopoClusters = CorrectCaloClusters( + "CorrectCaloTopoClusters", + inClusters=createTopoClusters.clusters.Path, + outClusters="Corrected" + createTopoClusters.clusters.Path, + systemIDs=[4], + numLayers=[ecalBarrelLayers], + firstLayerIDs=[0], + lastLayerIDs=[ecalBarrelLayers-1], + readoutNames=[ecalBarrelReadoutName], + # do not split the following line or it will break scripts that update the values of the corrections + upstreamParameters=upstreamParameters, + upstreamFormulas=[['[0]+[1]/(x-[2])', '[0]+[1]/(x-[2])']], + # do not split the following line or it will break scripts that update the values of the corrections + downstreamParameters=downstreamParameters, + downstreamFormulas=[['[0]+[1]*x', '[0]+[1]/sqrt(x)', '[0]+[1]/x']], + OutputLevel=INFO + ) + + if addShapeParameters: + augmentCaloTopoClusters = AugmentClustersFCCee("augmentCaloTopoClusters", + inClusters=createTopoClusters.clusters.Path, + outClusters="Augmented" + createTopoClusters.clusters.Path, + systemIDs=[4], + systemNames=["EMB"], + numLayers=[ecalBarrelLayers], + readoutNames=[ecalBarrelReadoutName], + layerFieldNames=["layer"], + thetaRecalcWeights=[ecalBarrelThetaWeights], + do_photon_shapeVar=False, + OutputLevel=INFO) + +# Output +out = PodioOutput("out", + OutputLevel=INFO) +out.filename = "ALLEGRO_sim_digi_reco.root" + +# drop the unpositioned ECal and HCal barrel and endcap cells +out.outputCommands = ["keep *", + "drop emptyCaloCells", + "drop %s" % ecalBarrelCellsName, + "drop %s" % ecalEndcapCellsName] + +# drop the uncalibrated cells +out.outputCommands.append("drop %s" % ecalBarrelReadoutName) +out.outputCommands.append("drop %s" % ecalBarrelReadoutName2) +out.outputCommands.append("drop %s" % ecalEndcapReadoutName) + +# drop the intermediate ecal barrel cells in case of a resegmentation +if resegmentECalBarrel: + out.outputCommands.append("drop ECalBarrelCellsMerged") + out.outputCommands.append("drop %s" % ecalBarrelCellsName2) + +# drop lumi, vertex, DCH, Muons (unless want to keep for event display) +out.outputCommands.append("drop Lumi*") +# out.outputCommands.append("drop Vertex*") +# out.outputCommands.append("drop DriftChamber_simHits*") +out.outputCommands.append("drop MuonTagger*") + +# drop hits/positioned cells/cluster cells if desired +if not saveHits: + out.outputCommands.append("drop *%sContributions" % ecalBarrelReadoutName) + out.outputCommands.append("drop *%sContributions" % ecalBarrelReadoutName2) + out.outputCommands.append("drop *%sContributions" % ecalEndcapReadoutName) +if not saveCells: + out.outputCommands.append("drop %s" % ecalBarrelPositionedCellsName) + out.outputCommands.append("drop %s" % ecalEndcapPositionedCellsName) + if resegmentECalBarrel: + out.outputCommands.append("drop %s" % ecalBarrelPositionedCellsName2) +if not saveClusterCells: + out.outputCommands.append("drop Calo*ClusterCells*") + +# if we decorate the clusters, we can drop the non-decorated ones +if addShapeParameters: + out.outputCommands.append("drop %s" % augmentCaloClusters.inClusters) + +# CPU information +chra = ChronoAuditor() +audsvc = AuditorSvc() +audsvc.Auditors = [chra] +out.AuditExecute = True + +# Configure list of external services +ExtSvc = [geoservice, podioevent, audsvc] +if dumpGDML: + ExtSvc += [gdmldumpservice] + +# Setup alg sequence +TopAlg = [ + input_reader, + createEcalBarrelCells, + createEcalBarrelPositionedCells, +] +createEcalBarrelCells.AuditExecute = True +createEcalBarrelPositionedCells.AuditExecute = True + +if resegmentECalBarrel: + TopAlg += [ + resegmentEcalBarrelTool, + createEcalBarrelCells2, + createEcalBarrelPositionedCells2, + ] + resegmentEcalBarrelTool.AuditExecute = True + createEcalBarrelCells2.AuditExecute = True + createEcalBarrelPositionedCells2.AuditExecute = True + +if doSWClustering or doTopoClustering: + TopAlg += [createemptycells] + createemptycells.AuditExecute = True + + if doSWClustering: + TopAlg += [createClusters] + createClusters.AuditExecute = True + + if applyUpDownstreamCorrections: + TopAlg += [correctCaloClusters] + correctCaloClusters.AuditExecute = True + + if addShapeParameters: + TopAlg += [augmentCaloClusters] + augmentCaloClusters.AuditExecute = True + + if doTopoClustering: + TopAlg += [createTopoClusters] + createTopoClusters.AuditExecute = True + + if applyUpDownstreamCorrections: + TopAlg += [correctCaloTopoClusters] + correctCaloTopoClusters.AuditExecute = True + + if addShapeParameters: + TopAlg += [augmentCaloTopoClusters] + augmentCaloTopoClusters.AuditExecute = True + +TopAlg += [ + out +] + +ApplicationMgr( + TopAlg=TopAlg, + EvtSel='NONE', + EvtMax=Nevts, + ExtSvc=ExtSvc, + StopOnSignal=True, +)