From 54c0674d4afa5b0662a9967ff2669664a8f9c488 Mon Sep 17 00:00:00 2001 From: Giovanni Marchiori <39376142+giovannimarchiori@users.noreply.github.com> Date: Mon, 9 Sep 2024 18:06:29 +0200 Subject: [PATCH] expanded version of CreateCaloCells that also add proper cell positions (#100) * expanded version of CreateCaloCells that also add proper cell positions * GaudiAlgorithm -> Gaudi::Algorithm * improve comment * Follow Juraj's suggestion and fix a bug * remove commented includes * replace CreateCaloCells+CellPositions.. with CreatePositionedCaloCells. Include also hcal endcap in test, and move to new hcal positioning tool for hcal barrel) * remove commented includes * enable code to cluster hcal endcap cells * fix a bug in the shell script for automatic tests, allow option for running over particle gun samples * fix bug when crosstalk is off (readCrosstalkMap undefined) * move code to RecCalorimeter, remove unused pointer to geoSvc * remove multiple spaces --- .../components/CreatePositionedCaloCells.cpp | 189 +++++++++ .../components/CreatePositionedCaloCells.h | 97 +++++ .../components/CreateCaloCellPositionsFCCee.h | 2 +- .../tests/options/ALLEGRO_o1_v03_digi_reco.py | 359 +++++++++--------- 4 files changed, 460 insertions(+), 187 deletions(-) create mode 100644 RecCalorimeter/src/components/CreatePositionedCaloCells.cpp create mode 100644 RecCalorimeter/src/components/CreatePositionedCaloCells.h diff --git a/RecCalorimeter/src/components/CreatePositionedCaloCells.cpp b/RecCalorimeter/src/components/CreatePositionedCaloCells.cpp new file mode 100644 index 00000000..3c741d49 --- /dev/null +++ b/RecCalorimeter/src/components/CreatePositionedCaloCells.cpp @@ -0,0 +1,189 @@ +#include "CreatePositionedCaloCells.h" + +// k4geo +#include "detectorCommon/DetUtils_k4geo.h" + +// edm4hep +#include "edm4hep/CalorimeterHit.h" + +DECLARE_COMPONENT(CreatePositionedCaloCells) + +CreatePositionedCaloCells::CreatePositionedCaloCells(const std::string& name, ISvcLocator* svcLoc) : +Gaudi::Algorithm(name, svcLoc) { + declareProperty("hits", m_hits, "Hits from which to create cells (input)"); + declareProperty("cells", m_cells, "The created calorimeter cells (output)"); + + declareProperty("positionsTool", m_cellPositionsTool, "Handle for cell positions tool"); + declareProperty("crosstalkTool", m_crosstalkTool, "Handle for the cell crosstalk tool"); + declareProperty("calibTool", m_calibTool, "Handle for tool to calibrate Geant4 energy to EM scale tool"); + declareProperty("noiseTool", m_noiseTool, "Handle for the calorimeter cells noise tool"); + declareProperty("geometryTool", m_geoTool, "Handle for the geometry tool"); +} + +StatusCode CreatePositionedCaloCells::initialize() { + StatusCode sc = Gaudi::Algorithm::initialize(); + if (sc.isFailure()) return sc; + + info() << "CreatePositionedCaloCells initialized" << endmsg; + info() << "do calibration : " << m_doCellCalibration << endmsg; + info() << "add cell noise : " << m_addCellNoise << endmsg; + info() << "remove cells below threshold : " << m_filterCellNoise << endmsg; + info() << "emulate crosstalk : " << m_addCrosstalk << endmsg; + + // Initialization of tools + + // Cell position tool + if (!m_cellPositionsTool.retrieve()) { + error() << "Unable to retrieve the cell positions tool!!!" << endmsg; + return StatusCode::FAILURE; + } + // Cell crosstalk tool + if (m_addCrosstalk) { + if (!m_crosstalkTool.retrieve()) { + error() << "Unable to retrieve the cell crosstalk tool!!!" << endmsg; + return StatusCode::FAILURE; + } + } + // Calibrate Geant4 energy to EM scale tool + if (m_doCellCalibration) { + if (!m_calibTool.retrieve()) { + error() << "Unable to retrieve the calo cells calibration tool!!!" << endmsg; + return StatusCode::FAILURE; + } + } + // Cell noise tool + if (m_addCellNoise || m_filterCellNoise) { + if (!m_noiseTool.retrieve()) { + error() << "Unable to retrieve the calo cells noise tool!!!" << endmsg; + return StatusCode::FAILURE; + } + // Geometry settings + if (!m_geoTool.retrieve()) { + error() << "Unable to retrieve the geometry tool!!!" << endmsg; + return StatusCode::FAILURE; + } + // Prepare map of all existing cells in calorimeter to add noise to all + StatusCode sc_prepareCells = m_geoTool->prepareEmptyCells(m_cellsMap); + if (sc_prepareCells.isFailure()) { + error() << "Unable to create empty cells!" << endmsg; + return StatusCode::FAILURE; + } + } + + // Copy over the CellIDEncoding string from the input collection to the output collection + auto hitsEncoding = m_hitsCellIDEncoding.get_optional(); + if (!hitsEncoding.has_value()) { + error () << "Missing cellID encoding for input collection" << endmsg; + return StatusCode::FAILURE; + } + m_cellsCellIDEncoding.put(hitsEncoding.value()); + + return StatusCode::SUCCESS; +} + +StatusCode CreatePositionedCaloCells::execute(const EventContext&) const { + // Get the input collection with Geant4 hits + const edm4hep::SimCalorimeterHitCollection* hits = m_hits.get(); + debug() << "Input Hit collection size: " << hits->size() << endmsg; + + // 0. Clear all cells + if (m_addCellNoise) { + std::for_each(m_cellsMap.begin(), m_cellsMap.end(), [](std::pair& p) { p.second = 0; }); + } else { + m_cellsMap.clear(); + } + + + // 1. Merge energy deposits into cells + // If running with noise, map was already prepared in initialize(). + // Otherwise it is being created below + for (const auto& hit : *hits) { + verbose() << "CellID : " << hit.getCellID() << endmsg; + m_cellsMap[hit.getCellID()] += hit.getEnergy(); + } + debug() << "Number of calorimeter cells after merging of hits: " << m_cellsMap.size() << endmsg; + + // 2. Emulate cross-talk (if asked) + if (m_addCrosstalk) { + // Derive the cross-talk contributions without affecting yet the nominal energy + // (one has to emulate crosstalk based on cells free from any cross-talk contributions) + m_crosstalkCellsMap.clear(); // this is a temporary map to hold energy exchange due to cross-talk, without affecting yet the nominal energy + // loop over cells with nominal energies + for (const auto& this_cell : m_cellsMap) { + uint64_t this_cellId = this_cell.first; + auto vec_neighbours = m_crosstalkTool->getNeighbours(this_cellId); // a vector of neighbour IDs + auto vec_crosstalks = m_crosstalkTool->getCrosstalks(this_cellId); // a vector of crosstalk coefficients + // loop over crosstalk neighbours of the cell under study + for (unsigned int i_cell=0; i_cellcalibrate(m_cellsMap); + } + + // 4. Add noise to all cells + if (m_addCellNoise) { + m_noiseTool->addRandomCellNoise(m_cellsMap); + } + + // 5. Filter cells + if (m_filterCellNoise) { + m_noiseTool->filterCellNoise(m_cellsMap); + } + + // 6. Copy information to CaloHitCollection + edm4hep::CalorimeterHitCollection* edmCellsCollection = new edm4hep::CalorimeterHitCollection(); + for (const auto& cell : m_cellsMap) { + if (m_addCellNoise || (!m_addCellNoise && cell.second != 0)) { + auto newCell = edmCellsCollection->create(); + newCell.setEnergy(cell.second); + uint64_t cellid = cell.first; + newCell.setCellID(cellid); + + // add cell position + auto cached_pos = m_positions_cache.find(cellid); + if(cached_pos == m_positions_cache.end()) { + // retrieve position from tool + dd4hep::Position posCell = m_cellPositionsTool->xyzPosition(cellid); + edm4hep::Vector3f edmPos; + edmPos.x = posCell.x() / dd4hep::mm; + edmPos.y = posCell.y() / dd4hep::mm; + edmPos.z = posCell.z() / dd4hep::mm; + m_positions_cache[cellid] = edmPos; + newCell.setPosition(edmPos); + } + else { + newCell.setPosition(cached_pos->second); + } + + debug() << "Cell energy (GeV) : " << newCell.getEnergy() << "\tcellID " << newCell.getCellID() << endmsg; + debug() << "Position of cell (mm) : \t" << newCell.getPosition().x/dd4hep::mm + << "\t" << newCell.getPosition().y/dd4hep::mm + << "\t" << newCell.getPosition().z/dd4hep::mm << endmsg; + } + } + + // push the CaloHitCollection to event store + m_cells.put(edmCellsCollection); + + debug() << "Output Cell collection size: " << edmCellsCollection->size() << endmsg; + + return StatusCode::SUCCESS; +} + +StatusCode CreatePositionedCaloCells::finalize() { return Gaudi::Algorithm::finalize(); } diff --git a/RecCalorimeter/src/components/CreatePositionedCaloCells.h b/RecCalorimeter/src/components/CreatePositionedCaloCells.h new file mode 100644 index 00000000..0aae29d4 --- /dev/null +++ b/RecCalorimeter/src/components/CreatePositionedCaloCells.h @@ -0,0 +1,97 @@ +#ifndef RECCALORIMETER_CREATEPOSITIONEDCALOCELLS_H +#define RECCALORIMETER_CREATEPOSITIONEDCALOCELLS_H + +// k4FWCore +#include "k4FWCore/DataHandle.h" +#include "k4FWCore/MetaDataHandle.h" +#include "k4Interface/ICellPositionsTool.h" +#include "k4Interface/ICalibrateCaloHitsTool.h" +#include "k4Interface/ICalorimeterTool.h" +#include "k4Interface/INoiseCaloCellsTool.h" +#include "k4Interface/ICaloReadCrosstalkMap.h" + +// Gaudi +#include "Gaudi/Algorithm.h" +#include "GaudiKernel/ToolHandle.h" + +// edm4hep +#include "edm4hep/CalorimeterHitCollection.h" +#include "edm4hep/SimCalorimeterHitCollection.h" + +/** @class CreatePositionedCaloCells + * + * Algorithm for creating positioned calorimeter cells from Geant4 hits. + * Based on CreateCaloCells, adding a positioning tool so that output + * cells from the digitisation contain the correct position in 3D space. + * + * Flow of the program: + * 1/ Merge Geant4 energy deposits with same cellID + * 2/ Emulate cross-talk (if switched on) + * 3/ Calibrate to electromagnetic scale (if calibration switched on) + * 4/ Add random noise to each cell (if noise switched on) + * 5/ Filter cells and remove those with energy below threshold (if noise + + * filtering switched on) + * 6/ Add cell positions + * + * Tools called: + * - CalibrateCaloHitsTool + * - NoiseCaloCellsTool + * - CaloReadCrosstalkMap + * - CalorimeterTool (for geometry) + * - CellPositionsTool + * + * @author Giovanni Marchiori + * @date 2024-07 + * + */ + +class CreatePositionedCaloCells : public Gaudi::Algorithm { + +public: + CreatePositionedCaloCells(const std::string& name, ISvcLocator* svcLoc); + + StatusCode initialize(); + + StatusCode execute(const EventContext&) const; + + StatusCode finalize(); + +private: + /// Handle for tool to get cells positions + ToolHandle m_cellPositionsTool{"CellPositionsTool", this}; + /// Handle for the calorimeter cells crosstalk tool + mutable ToolHandle m_crosstalkTool{"ReadCaloCrosstalkMap", this}; + /// Handle for tool to calibrate Geant4 energy to EM scale tool + mutable ToolHandle m_calibTool{"CalibrateCaloHitsTool", this}; + /// Handle for the calorimeter cells noise tool + mutable ToolHandle m_noiseTool{"NoiseCaloCellsFlatTool", this}; + /// Handle for the geometry tool + ToolHandle m_geoTool{"TubeLayerPhiEtaCaloTool", this}; + + /// Add crosstalk to cells? + Gaudi::Property m_addCrosstalk{this, "addCrosstalk", false, "Add crosstalk effect?"}; + /// Calibrate to EM scale? + Gaudi::Property m_doCellCalibration{this, "doCellCalibration", true, "Calibrate to EM scale?"}; + /// Add noise to cells? + Gaudi::Property m_addCellNoise{this, "addCellNoise", true, "Add noise to cells?"}; + /// Save only cells with energy above threshold? + Gaudi::Property m_filterCellNoise{this, "filterCellNoise", false, + "Save only cells with energy above threshold?"}; + + /// Handle for calo hits (input collection) + mutable DataHandle m_hits{"hits", Gaudi::DataHandle::Reader, this}; + /// Handle for the cellID encoding string of the input collection + MetaDataHandle m_hitsCellIDEncoding{m_hits, edm4hep::labels::CellIDEncoding, Gaudi::DataHandle::Reader}; + /// Handle for calo cells (output collection) + mutable DataHandle m_cells{"cells", Gaudi::DataHandle::Writer, this}; + MetaDataHandle m_cellsCellIDEncoding{m_cells, edm4hep::labels::CellIDEncoding, Gaudi::DataHandle::Writer}; + /// Maps of cell IDs (corresponding to DD4hep IDs) on final energies to be used for clustering + mutable std::unordered_map m_cellsMap; + /// Maps of cell IDs (corresponding to DD4hep IDs) on transfer of signals due to crosstalk + mutable std::unordered_map m_crosstalkCellsMap; + + /// Cache position vs cellID + mutable std::unordered_map m_positions_cache{}; +}; + +#endif /* RECCALORIMETER_CREATEPOSITIONEDCALOCELLS_H */ diff --git a/RecFCCeeCalorimeter/src/components/CreateCaloCellPositionsFCCee.h b/RecFCCeeCalorimeter/src/components/CreateCaloCellPositionsFCCee.h index 22890915..675c955b 100644 --- a/RecFCCeeCalorimeter/src/components/CreateCaloCellPositionsFCCee.h +++ b/RecFCCeeCalorimeter/src/components/CreateCaloCellPositionsFCCee.h @@ -47,7 +47,7 @@ class CreateCaloCellPositionsFCCee : public Gaudi::Algorithm { StatusCode finalize(); private: - /// Handle for tool to get positions width + /// Handle for tool to get positions ToolHandle m_cellPositionsTool{}; /// Input collection mutable DataHandle m_hits{"hits/hits", Gaudi::DataHandle::Reader, this}; diff --git a/RecFCCeeCalorimeter/tests/options/ALLEGRO_o1_v03_digi_reco.py b/RecFCCeeCalorimeter/tests/options/ALLEGRO_o1_v03_digi_reco.py index c55207b0..2aa04472 100644 --- a/RecFCCeeCalorimeter/tests/options/ALLEGRO_o1_v03_digi_reco.py +++ b/RecFCCeeCalorimeter/tests/options/ALLEGRO_o1_v03_digi_reco.py @@ -23,11 +23,9 @@ # - what to save in output file # -# always drop unpositioned / uncalibrated cells, excepts for tests and debugging +# always drop uncalibrated cells, except for tests and debugging # dropUncalibratedCells = True -# dropUnpositionedCells = True dropUncalibratedCells = False -dropUnpositionedCells = False # 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 @@ -105,18 +103,20 @@ # 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) +ecalBarrelReadoutName = "ECalBarrelModuleThetaMerged" # barrel, original segmentation (baseline) +ecalBarrelReadoutName2 = "ECalBarrelModuleThetaMerged2" # barrel, after re-segmentation (for optimisation studies) +ecalEndcapReadoutName = "ECalEndcapTurbine" # endcap, turbine-like (baseline) # - HCAL readouts if runHCal: - hcalBarrelReadoutName = "HCalBarrelReadout" # barrel, original segmentation (row-phi) - hcalBarrelReadoutName2 = "BarHCal_Readout_phitheta" # barrel, groups together cells of different row within same theta slice - hcalEndcapReadoutName = "HCalEndcapReadout" # endcap, original segmentation + hcalBarrelReadoutName = "HCalBarrelReadout" # barrel, original segmentation (row-phi) + hcalBarrelReadoutName2 = "BarHCal_Readout_phitheta" # barrel, groups together cells of different row within same theta slice + hcalEndcapReadoutName = "HCalEndcapReadout" # endcap, original segmentation + hcalEndcapReadoutName2 = "HCalEndcapReadout_phitheta" # endcap, groups together cells of different row within same theta slice else: hcalBarrelReadoutName = "" hcalBarrelReadoutName2 = "" hcalEndcapReadoutName = "" + hcalEndcapReadoutName2 = "" # - EM scale calibration (sampling fraction) from Configurables import CalibrateInLayersTool @@ -140,10 +140,57 @@ calibHCalEndcap = CalibrateCaloHitsTool( "CalibrateHCalEndcap", invSamplingFraction="29.4202") # FIXME: to be updated for ddsim -from Configurables import CreateCaloCells -from Configurables import CreateCaloCellPositionsFCCee +# - cell positioning tools from Configurables import CellPositionsECalBarrelModuleThetaSegTool +cellPositionEcalBarrelTool = CellPositionsECalBarrelModuleThetaSegTool( + "CellPositionsECalBarrel", + readoutName=ecalBarrelReadoutName, + OutputLevel=INFO +) +if resegmentECalBarrel: + cellPositionEcalBarrelTool2 = CellPositionsECalBarrelModuleThetaSegTool( + "CellPositionsECalBarrel2", + readoutName=ecalBarrelReadoutName2, + OutputLevel=INFO + ) + from Configurables import CellPositionsECalEndcapTurbineSegTool +cellPositionEcalEndcapTool = CellPositionsECalEndcapTurbineSegTool( + "CellPositionsECalEndcap", + readoutName=ecalEndcapReadoutName, + OutputLevel=INFO +) + +if runHCal: + from Configurables import CellPositionsHCalPhiThetaSegTool + cellPositionHCalBarrelTool = CellPositionsHCalPhiThetaSegTool( + "CellPositionsHCalBarrel", + readoutName=hcalBarrelReadoutName, + detectorName="HCalBarrel", + OutputLevel=INFO + ) + cellPositionHCalBarrelTool2 = CellPositionsHCalPhiThetaSegTool( + "CellPositionsHCalBarrel2", + readoutName=hcalBarrelReadoutName2, + detectorName="HCalBarrel", + OutputLevel=INFO + ) + cellPositionHCalEndcapTool = CellPositionsHCalPhiThetaSegTool( + "CellPositionsHCalEndcap", + readoutName=hcalEndcapReadoutName, + detectorName="HCalThreePartsEndcap", + numLayersHCalThreeParts=[6, 9, 22], + OutputLevel=INFO + ) + cellPositionHCalEndcapTool2 = CellPositionsHCalPhiThetaSegTool( + "CellPositionsHCalEndcap2", + readoutName=hcalEndcapReadoutName2, + detectorName="HCalThreePartsEndcap", + numLayersHCalThreeParts=[6, 9, 22], + OutputLevel=INFO + ) + +# - crosstalk tool if addCrosstalk: from Configurables import ReadCaloCrosstalkMap # read the crosstalk map @@ -153,44 +200,28 @@ else: readCrosstalkMap = None -# 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" -createEcalBarrelCells = CreateCaloCells("CreateECalBarrelCells", - doCellCalibration=True, - calibTool=calibEcalBarrel, - crosstalksTool=readCrosstalkMap, - addCrosstalk=addCrosstalk, - 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 -) +# Create cells in ECal barrel (calibrated and positioned - optionally with xtalk and noise added) +# from uncalibrated cells (+cellID info) from ddsim ecalBarrelPositionedCellsName = ecalBarrelReadoutName + "Positioned" -createEcalBarrelPositionedCells = CreateCaloCellPositionsFCCee( - "CreateECalBarrelPositionedCells", - OutputLevel=INFO -) -createEcalBarrelPositionedCells.positionsTool = cellPositionEcalBarrelTool -createEcalBarrelPositionedCells.hits.Path = ecalBarrelCellsName -createEcalBarrelPositionedCells.positionedHits.Path = ecalBarrelPositionedCellsName +from Configurables import CreatePositionedCaloCells +createEcalBarrelCells = CreatePositionedCaloCells("CreatePositionedECalBarrelCells", + doCellCalibration=True, + positionsTool=cellPositionEcalBarrelTool, + calibTool=calibEcalBarrel, + crosstalkTool=readCrosstalkMap, + addCrosstalk=addCrosstalk, + addCellNoise=False, + filterCellNoise=False, + OutputLevel=INFO, + hits=ecalBarrelReadoutName, + cells=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 + # 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 + # Step a: compute new cellID of cells based on new readout # (merged module-theta segmentation with variable merging vs layer) from Configurables import RedoSegmentation resegmentEcalBarrelTool = RedoSegmentation("ReSegmentationEcal", @@ -202,64 +233,40 @@ newReadoutName=ecalBarrelReadoutName2, OutputLevel=INFO, debugPrint=200, - inhits=ecalBarrelCellsName, + inhits=ecalBarrelPositionedCellsName, outhits="ECalBarrelCellsMerged") - # Step 2b: merge new cells with same cellID together + # Step b: 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 - ) + # noise and xtalk off assuming they were applied earlier ecalBarrelPositionedCellsName2 = ecalBarrelReadoutName2 + "Positioned" - createEcalBarrelPositionedCells2 = CreateCaloCellPositionsFCCee( - "CreateECalBarrelPositionedCells2", - OutputLevel=INFO - ) - createEcalBarrelPositionedCells2.positionsTool = cellPositionEcalBarrelTool2 - createEcalBarrelPositionedCells2.hits.Path = ecalBarrelCellsName2 - createEcalBarrelPositionedCells2.positionedHits.Path = ecalBarrelPositionedCellsName2 - + createEcalBarrelCells2 = CreatePositionedCaloCells("CreatePositionedECalBarrelCells2", + doCellCalibration=False, + positionsTool=cellPositionEcalBarrelTool2, + calibTool=None, + crosstalkTool=None, + addCrosstalk=False, + addCellNoise=False, + filterCellNoise=False, + OutputLevel=INFO, + hits="ECalBarrelCellsMerged", + cells=ecalBarrelPositionedCellsName2) # 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 - +createEcalEndcapCells = CreatePositionedCaloCells("CreatePositionedECalEndcapCells", + doCellCalibration=True, + positionsTool=cellPositionEcalEndcapTool, + calibTool=calibEcalEndcap, + crosstalkTool=None, + addCrosstalk=False, + addCellNoise=False, + filterCellNoise=False, + OutputLevel=INFO, + hits=ecalEndcapReadoutName, + cells=ecalEndcapPositionedCellsName) if addNoise: ecalBarrelNoisePath = "elecNoise_ecalBarrelFCCee_theta.root" @@ -312,46 +319,28 @@ if runHCal: - # Create cells in HCal barrel - # 1 - merge hits into cells with the default readout - hcalBarrelCellsName = "HCalBarrelCells" - createHCalBarrelCells = CreateCaloCells("CreateHCalBarrelCells", - doCellCalibration=True, - calibTool=calibHCalBarrel, - addCellNoise=False, - filterCellNoise=False, - addPosition=True, - hits=hcalBarrelReadoutName, - cells=hcalBarrelCellsName, - OutputLevel=INFO) - - # 2 - attach positions to the cells (cell positions needed for RedoSegmentation!) - from Configurables import CellPositionsHCalBarrelPhiThetaSegTool - cellPositionHCalBarrelTool = CellPositionsHCalBarrelPhiThetaSegTool( - "CellPositionsHCalBarrel", - readoutName=hcalBarrelReadoutName, - OutputLevel=INFO - ) - hcalBarrelPositionedCellsName = "HCalBarrelPositionedCells" - createHCalBarrelPositionedCells = CreateCaloCellPositionsFCCee( - "CreateHCalBarrelPositionedCells", - OutputLevel=INFO - ) - createHCalBarrelPositionedCells.positionsTool = cellPositionHCalBarrelTool - createHCalBarrelPositionedCells.hits.Path = hcalBarrelCellsName - createHCalBarrelPositionedCells.positionedHits.Path = hcalBarrelPositionedCellsName - - # 3 - compute new cellID of cells based on new readout - removing row information + # Apply calibration and positioning to cells in HCal barrel + hcalBarrelPositionedCellsName = hcalBarrelReadoutName + "Positioned" + createHCalBarrelCells = CreatePositionedCaloCells("CreateHCalBarrelCells", + doCellCalibration=True, + calibTool=calibHCalBarrel, + positionsTool=cellPositionHCalBarrelTool, + addCellNoise=False, + filterCellNoise=False, + hits=hcalBarrelReadoutName, + cells=hcalBarrelPositionedCellsName, + OutputLevel=INFO) + + # Compute new cellID of cells based on new readout - removing row information # We use a RedoSegmentation. Using a RewriteBitField with removeIds=["row"], # wont work because there are tiles with same layer/theta/phi but different row # as a consequence there will be multiple cells with same cellID in the output collection # and this will screw up the SW clustering - hcalBarrelCellsName2 = "HCalBarrelCells2" # first we create new hits with the readout without the row information - # and then merge them into new cells + # and then merge them into new cells, wihotut applying the calibration again from Configurables import RedoSegmentation - rewriteHCalBarrel = RedoSegmentation("ReSegmentationHCal", + rewriteHCalBarrel = RedoSegmentation("ReSegmentationHCalBarrel", # old bitfield (readout) oldReadoutName=hcalBarrelReadoutName, # specify which fields are going to be altered (deleted/rewritten) @@ -363,47 +352,59 @@ inhits=hcalBarrelPositionedCellsName, outhits="HCalBarrelCellsWithoutRow") - createHCalBarrelCells2 = CreateCaloCells("CreateHCalBarrelCells2", - doCellCalibration=False, - addCellNoise=False, - filterCellNoise=False, - OutputLevel=INFO, - hits=rewriteHCalBarrel.outhits.Path, - cells=hcalBarrelCellsName2) - - # 4 - attach positions to the new cells - from Configurables import CellPositionsHCalBarrelPhiThetaSegTool hcalBarrelPositionedCellsName2 = hcalBarrelReadoutName2 + "Positioned" - cellPositionHCalBarrelTool2 = CellPositionsHCalBarrelPhiThetaSegTool( - "CellPositionsHCalBarrel2", - readoutName=hcalBarrelReadoutName2, - OutputLevel=INFO - ) - createHCalBarrelPositionedCells2 = CreateCaloCellPositionsFCCee( - "CreateHCalBarrelPositionedCells2", - OutputLevel=INFO - ) - createHCalBarrelPositionedCells2.positionsTool = cellPositionHCalBarrelTool2 - createHCalBarrelPositionedCells2.hits.Path = hcalBarrelCellsName2 - createHCalBarrelPositionedCells2.positionedHits.Path = hcalBarrelPositionedCellsName2 + createHCalBarrelCells2 = CreatePositionedCaloCells("CreateHCalBarrelCells2", + doCellCalibration=False, + positionsTool=cellPositionHCalBarrelTool2, + addCellNoise=False, + filterCellNoise=False, + OutputLevel=INFO, + hits=rewriteHCalBarrel.outhits.Path, + cells=hcalBarrelPositionedCellsName2) # Create cells in HCal endcap - # createHCalEndcapCells = CreateCaloCells("CreateHCalEndcapCaloCells", - # doCellCalibration=True, - # calibTool=calibHCalEndcap, - # addCellNoise=False, - # filterCellNoise=False, - # OutputLevel=INFO) - # createHCalEndcapCells.hits.Path="HCalEndcapHits" - # createHCalEndcapCells.cells.Path="HCalEndcapCells" + hcalEndcapPositionedCellsName = hcalEndcapReadoutName + "Positioned" + createHCalEndcapCells = CreatePositionedCaloCells("CreateHCalEndcapCells", + doCellCalibration=True, + calibTool=calibHCalEndcap, + addCellNoise=False, + filterCellNoise=False, + positionsTool=cellPositionHCalEndcapTool, + OutputLevel=INFO, + hits=hcalEndcapReadoutName, + cells=hcalEndcapPositionedCellsName) + + rewriteHCalEndcap = RedoSegmentation("ReSegmentationHCalEndcap", + # old bitfield (readout) + oldReadoutName=hcalEndcapReadoutName, + # specify which fields are going to be altered (deleted/rewritten) + oldSegmentationIds=["row", "theta", "phi"], + # new bitfield (readout), with new segmentation (theta-phi grid) + newReadoutName=hcalEndcapReadoutName2, + OutputLevel=INFO, + debugPrint=200, + inhits=hcalEndcapPositionedCellsName, + outhits="HCalEndcapCellsWithoutRow") + + hcalEndcapPositionedCellsName2 = hcalEndcapReadoutName2 + "Positioned" + createHCalEndcapCells2 = CreatePositionedCaloCells("CreateHCalEndcapCells2", + doCellCalibration=False, + positionsTool=cellPositionHCalEndcapTool2, + addCellNoise=False, + filterCellNoise=False, + OutputLevel=INFO, + hits=rewriteHCalEndcap.outhits.Path, + cells=hcalEndcapPositionedCellsName2) else: - hcalBarrelCellsName = "emptyCaloCells" hcalBarrelPositionedCellsName = "emptyCaloCells" - hcalBarrelCellsName2 = "emptyCaloCells" hcalBarrelPositionedCellsName2 = "emptyCaloCells" + hcalEndcapPositionedCellsName = "emptyCaloCells" + hcalEndcapPositionedCellsName2 = "emptyCaloCells" cellPositionHCalBarrelTool = None cellPositionHCalBarrelTool2 = None + cellPositionHCalEndcapTool = None + cellPositionHCalEndcapTool2 = None # Empty cells for parts of calorimeter not implemented yet from Configurables import CreateEmptyCaloCellsCollection @@ -579,7 +580,7 @@ ecalFwdReadoutName="", hcalBarrelReadoutName=hcalBarrelReadoutName2, hcalExtBarrelReadoutName="", - hcalEndcapReadoutName="", + hcalEndcapReadoutName=hcalEndcapReadoutName2, hcalFwdReadoutName="", OutputLevel=INFO) towers.ecalBarrelCells.Path = ecalBarrelPositionedCellsName @@ -587,7 +588,7 @@ towers.ecalFwdCells.Path = "emptyCaloCells" towers.hcalBarrelCells.Path = hcalBarrelPositionedCellsName2 towers.hcalExtBarrelCells.Path = "emptyCaloCells" - towers.hcalEndcapCells.Path = "emptyCaloCells" + towers.hcalEndcapCells.Path = hcalEndcapPositionedCellsName2 towers.hcalFwdCells.Path = "emptyCaloCells" createClusters = CreateCaloClustersSlidingWindowFCCee("CreateCaloClusters", @@ -822,27 +823,17 @@ out.outputCommands.append("drop %s" % ecalEndcapReadoutName) if runHCal: out.outputCommands.append("drop %s" % hcalBarrelReadoutName) - # out.outputCommands.append("drop %s" % hcalEndcapReadoutName) - -# drop the unpositioned ECal and HCal barrel and endcap cells -if dropUnpositionedCells: - if runHCal: - out.outputCommands += ["drop %s" % ecalBarrelCellsName, - "drop %s" % ecalEndcapCellsName, - "drop %s" % hcalBarrelCellsName, - "drop %s" % hcalBarrelPositionedCellsName, - "drop %s" % hcalBarrelCellsName2, - "drop %s" % rewriteHCalBarrel.outhits.Path] - + out.outputCommands.append("drop %s" % hcalEndcapReadoutName) else: - out.outputCommands += ["drop HCal*", - "drop %s" % ecalBarrelCellsName, - "drop %s" % ecalEndcapCellsName] + out.outputCommands += ["drop HCal*"] # 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 the intermediate hcal barrel cells before resegmentation + if runHCal: + out.outputCommands.append("drop %s" % hcalBarrelPositionedCellsName) + out.outputCommands.append("drop %s" % hcalEndcapPositionedCellsName) # drop lumi, vertex, DCH, Muons (unless want to keep for event display) out.outputCommands.append("drop Lumi*") @@ -862,6 +853,7 @@ out.outputCommands.append("drop %s" % ecalBarrelPositionedCellsName2) if runHCal: out.outputCommands.append("drop %s" % hcalBarrelPositionedCellsName2) + out.outputCommands.append("drop %s" % hcalEndcapPositionedCellsName2) if not saveClusterCells: out.outputCommands.append("drop Calo*ClusterCells*") @@ -886,39 +878,34 @@ TopAlg = [ input_reader, createEcalBarrelCells, - createEcalBarrelPositionedCells, createEcalEndcapCells, - createEcalEndcapPositionedCells, ] createEcalBarrelCells.AuditExecute = True -createEcalBarrelPositionedCells.AuditExecute = True createEcalEndcapCells.AuditExecute = True -createEcalEndcapPositionedCells.AuditExecute = True if resegmentECalBarrel: TopAlg += [ resegmentEcalBarrelTool, createEcalBarrelCells2, - createEcalBarrelPositionedCells2, ] resegmentEcalBarrelTool.AuditExecute = True createEcalBarrelCells2.AuditExecute = True - createEcalBarrelPositionedCells2.AuditExecute = True if runHCal: TopAlg += [ createHCalBarrelCells, - createHCalBarrelPositionedCells, rewriteHCalBarrel, createHCalBarrelCells2, - createHCalBarrelPositionedCells2, - # createHCalEndcapCells + createHCalEndcapCells, + rewriteHCalEndcap, + createHCalEndcapCells2, ] createHCalBarrelCells.AuditExecute = True - createHCalBarrelPositionedCells.AuditExecute = True rewriteHCalBarrel.AuditExecute = True createHCalBarrelCells2.AuditExecute = True - createHCalBarrelPositionedCells2.AuditExecute = True + createHCalEndcapCells.AuditExecute = True + rewriteHCalEndcap.AuditExecute = True + createHCalEndcapCells2.AuditExecute = True if doSWClustering or doTopoClustering: TopAlg += [createemptycells]