diff --git a/RecCalorimeter/src/components/ReadNoiseFromFileTool.cpp b/RecCalorimeter/src/components/ReadNoiseFromFileTool.cpp index 032b88b9..d7c3f6cc 100644 --- a/RecCalorimeter/src/components/ReadNoiseFromFileTool.cpp +++ b/RecCalorimeter/src/components/ReadNoiseFromFileTool.cpp @@ -19,6 +19,7 @@ DECLARE_COMPONENT(ReadNoiseFromFileTool) ReadNoiseFromFileTool::ReadNoiseFromFileTool(const std::string& type, const std::string& name, const IInterface* parent) : GaudiTool(type, name, parent) { declareInterface(this); + declareProperty("cellPositionsTool", m_cellPositionsTool, "Handle for tool to retrieve cell positions"); } StatusCode ReadNoiseFromFileTool::initialize() { @@ -29,10 +30,24 @@ StatusCode ReadNoiseFromFileTool::initialize() { << "Make sure you have GeoSvc and SimSvc in the right order in the configuration." << endmsg; return StatusCode::FAILURE; } - m_segmentation = dynamic_cast(m_geoSvc->getDetector()->readout(m_readoutName).segmentation().segmentation()); - if (m_segmentation == nullptr) { - error() << "There is no phi-eta segmentation!!!!" << endmsg; - return StatusCode::FAILURE; + + // Check if cell position tool available if m_useSeg==false; if tool not + // available, try using segmentation instead + if (!m_cellPositionsTool.retrieve() and !m_useSeg) { + info() << "Unable to retrieve cell positions tool, try eta-phi segmentation." << endmsg; + m_useSeg = true; + } + + // Get PhiEta segmentation + if (m_useSeg) { + m_segmentation = dynamic_cast( + m_geoSvc->getDetector()->readout(m_readoutName).segmentation().segmentation()); + if (m_segmentation == nullptr) { + error() << "There is no phi-eta segmentation." << endmsg; + return StatusCode::FAILURE; + } + else + info() << "Found phi-eta segmentation." << endmsg; } // open and check file, read the histograms with noise constants @@ -43,6 +58,10 @@ StatusCode ReadNoiseFromFileTool::initialize() { // Take readout bitfield decoder from GeoSvc m_decoder = m_geoSvc->getDetector()->readout(m_readoutName).idSpec().decoder(); + if (m_decoder == nullptr) { + error() << "Cannot create decore for readout " << m_readoutName << endmsg; + return StatusCode::FAILURE; + } StatusCode sc = GaudiTool::initialize(); if (sc.isFailure()) return sc; @@ -134,8 +153,11 @@ double ReadNoiseFromFileTool::getNoiseConstantPerCell(uint64_t aCellId) { // Get cell coordinates: eta and radial layer dd4hep::DDSegmentation::CellID cID = aCellId; - double cellEta = m_segmentation->eta(cID); - + double cellEta; + if (m_useSeg) + cellEta = m_segmentation->eta(aCellId); + else + cellEta = m_cellPositionsTool->xyzPosition(aCellId).Eta(); unsigned cellLayer = m_decoder->get(cID, m_activeFieldName); // All histograms have same binning, all bins with same size @@ -143,12 +165,16 @@ double ReadNoiseFromFileTool::getNoiseConstantPerCell(uint64_t aCellId) { unsigned index = 0; if (m_histoElecNoiseConst.size() != 0) { int Nbins = m_histoElecNoiseConst.at(index).GetNbinsX(); + // GM: basically the same as FindBin ... + /* double deltaEtaBin = (m_histoElecNoiseConst.at(index).GetBinLowEdge(Nbins) + m_histoElecNoiseConst.at(index).GetBinWidth(Nbins) - m_histoElecNoiseConst.at(index).GetBinLowEdge(1)) / Nbins; // find the eta bin for the cell int ibin = floor(fabs(cellEta) / deltaEtaBin) + 1; + */ + int ibin = m_histoElecNoiseConst.at(index).FindBin(fabs(cellEta)); if (ibin > Nbins) { error() << "eta outside range of the histograms! Cell eta: " << cellEta << " Nbins in histogram: " << Nbins << endmsg; @@ -189,7 +215,11 @@ double ReadNoiseFromFileTool::getNoiseOffsetPerCell(uint64_t aCellId) { // Get cell coordinates: eta and radial layer dd4hep::DDSegmentation::CellID cID = aCellId; - double cellEta = m_segmentation->eta(cID); + double cellEta; + if (m_useSeg) + cellEta = m_segmentation->eta(aCellId); + else + cellEta = m_cellPositionsTool->xyzPosition(aCellId).Eta(); unsigned cellLayer = m_decoder->get(cID, m_activeFieldName); // All histograms have same binning, all bins with same size @@ -197,12 +227,16 @@ double ReadNoiseFromFileTool::getNoiseOffsetPerCell(uint64_t aCellId) { unsigned index = 0; if (m_histoElecNoiseOffset.size() != 0) { int Nbins = m_histoElecNoiseOffset.at(index).GetNbinsX(); - double deltaEtaBin = + // find the eta bin for the cell + // GM: basically the same as FindBin ... + /* + double deltaEtaBin = (m_histoElecNoiseOffset.at(index).GetBinLowEdge(Nbins) + m_histoElecNoiseOffset.at(index).GetBinWidth(Nbins) - m_histoElecNoiseOffset.at(index).GetBinLowEdge(1)) / Nbins; - // find the eta bin for the cell - int ibin = floor(fabs(cellEta) / deltaEtaBin) + 1; + int ibin = floor(fabs(cellEta) / deltaEtaBin) + 1; + */ + int ibin = m_histoElecNoiseOffset.at(index).FindBin(fabs(cellEta)); if (ibin > Nbins) { error() << "eta outside range of the histograms! Cell eta: " << cellEta << " Nbins in histogram: " << Nbins << endmsg; diff --git a/RecCalorimeter/src/components/ReadNoiseFromFileTool.h b/RecCalorimeter/src/components/ReadNoiseFromFileTool.h index fbbc55b0..369dfb38 100644 --- a/RecCalorimeter/src/components/ReadNoiseFromFileTool.h +++ b/RecCalorimeter/src/components/ReadNoiseFromFileTool.h @@ -7,11 +7,12 @@ // FCCSW #include "k4FWCore/DataHandle.h" -#include "DetSegmentation/FCCSWGridPhiEta.h" #include "k4Interface/ICalorimeterTool.h" #include "k4Interface/INoiseConstTool.h" #include "k4Interface/ICellPositionsTool.h" +#include "DetSegmentation/FCCSWGridPhiEta.h" + class IGeoSvc; // Root @@ -43,15 +44,18 @@ class ReadNoiseFromFileTool : public GaudiTool, virtual public INoiseConstTool { double getNoiseOffsetPerCell(uint64_t aCellID); private: + /// Handle for tool to get cell positions + ToolHandle m_cellPositionsTool; + /// use segmentation (eta-phi only so far!) in case no cell position tool is defined + Gaudi::Property m_useSeg{this, "useSegmentation", true, "Specify if segmentation is used to determine cell position."}; /// Add pileup contribution to the electronics noise? (only if read from file) Gaudi::Property m_addPileup{this, "addPileup", true, "Add pileup contribution to the electronics noise? (only if read from file)"}; /// Noise offset, if false, mean is set to 0 Gaudi::Property m_setNoiseOffset{this, "setNoiseOffset", true, "Set a noise offset per cell"}; - /// Name of the file with noise constants Gaudi::Property m_noiseFileName{this, "noiseFileName", "", "Name of the file with noise constants"}; - /// Name of the detector readout + /// Name of the detector readout (needed to get the decoder to extract e.g. layer information) Gaudi::Property m_readoutName{this, "readoutName", "ECalHitsPhiEta", "Name of the detector readout"}; /// Name of active layers for sampling calorimeter Gaudi::Property m_activeFieldName{this, "activeFieldName", "active_layer", @@ -66,28 +70,23 @@ class ReadNoiseFromFileTool : public GaudiTool, virtual public INoiseConstTool { "Name of electronics noise offset histogram"}; /// Name of pileup offset histogram Gaudi::Property m_pileupOffsetHistoName{this, "pileupOffsetHistoName", "h_pileup_layer", "Name of pileup offset histogram"}; - /// Number of radial layers Gaudi::Property m_numRadialLayers{this, "numRadialLayers", 3, "Number of radial layers"}; - /// Factor to apply to the noise values to get them in GeV if e.g. they were produced in MeV Gaudi::Property m_scaleFactor{this, "scaleFactor", 1, "Factor to apply to the noise values"}; - /// Histograms with pileup constants (index in array - radial layer) std::vector m_histoPileupConst; /// Histograms with electronics noise constants (index in array - radial layer) std::vector m_histoElecNoiseConst; - /// Histograms with pileup offset (index in array - radial layer) std::vector m_histoPileupOffset; /// Histograms with electronics noise offset (index in array - radial layer) std::vector m_histoElecNoiseOffset; - /// Pointer to the geometry service SmartIF m_geoSvc; /// PhiEta segmentation dd4hep::DDSegmentation::FCCSWGridPhiEta* m_segmentation; - // Decoder + // Decoder for ECal layers dd4hep::DDSegmentation::BitFieldCoder* m_decoder; }; diff --git a/RecFCCeeCalorimeter/CMakeLists.txt b/RecFCCeeCalorimeter/CMakeLists.txt index fd7560a5..0c7f74ec 100644 --- a/RecFCCeeCalorimeter/CMakeLists.txt +++ b/RecFCCeeCalorimeter/CMakeLists.txt @@ -12,6 +12,7 @@ gaudi_add_module(k4RecFCCeeCalorimeterPlugins DD4hep::DDCore EDM4HEP::edm4hep FCCDetectors::DetSegmentation + FCCDetectors::DetCommon DD4hep::DDG4 ROOT::Core ROOT::Hist diff --git a/RecFCCeeCalorimeter/src/components/CaloTopoClusterFCCee.cpp b/RecFCCeeCalorimeter/src/components/CaloTopoClusterFCCee.cpp index 646db5f6..7e4669d5 100644 --- a/RecFCCeeCalorimeter/src/components/CaloTopoClusterFCCee.cpp +++ b/RecFCCeeCalorimeter/src/components/CaloTopoClusterFCCee.cpp @@ -117,8 +117,15 @@ StatusCode CaloTopoClusterFCCee::execute() { }); std::map>> preClusterCollection; - CaloTopoClusterFCCee::buildingProtoCluster(m_neighbourSigma, m_lastNeighbourSigma, firstSeeds, m_allCells, - preClusterCollection); + StatusCode sc_buildProtoClusters = CaloTopoClusterFCCee::buildingProtoCluster(m_neighbourSigma, + m_lastNeighbourSigma, + firstSeeds, m_allCells, + preClusterCollection); + if (sc_buildProtoClusters.isFailure()) { + error() << "Unable to build the protoclusters!" << endmsg; + return StatusCode::FAILURE; + } + // Build Clusters in edm debug() << "Building " << preClusterCollection.size() << " cluster." << endmsg; double checkTotEnergy = 0.; @@ -132,10 +139,10 @@ StatusCode CaloTopoClusterFCCee::execute() { double energy = 0.; double deltaR = 0.; std::vector posPhi (i.second.size()); - std::vector posEta (i.second.size()); + std::vector posTheta (i.second.size()); std::vector vecEnergy (i.second.size()); double sumPhi = 0.; - double sumEta = 0.; + double sumTheta = 0.; std::map system; for (auto pair : i.second) { @@ -178,10 +185,10 @@ StatusCode CaloTopoClusterFCCee::execute() { posY += posCell.Y() * newCell.getEnergy(); posZ += posCell.Z() * newCell.getEnergy(); posPhi.push_back(posCell.Phi()); - posEta.push_back(posCell.Eta()); + posTheta.push_back(posCell.Theta()); vecEnergy.push_back(newCell.getEnergy()); sumPhi += posCell.Phi() * newCell.getEnergy(); - sumEta += posCell.Eta() * newCell.getEnergy(); + sumTheta += posCell.Theta() * newCell.getEnergy(); cluster.addToHits(newCell); auto er = m_allCells.erase(cID); @@ -193,10 +200,10 @@ StatusCode CaloTopoClusterFCCee::execute() { cluster.setPosition(edm4hep::Vector3f(posX / energy, posY / energy, posZ / energy)); // store deltaR of cluster in time for the moment.. sumPhi = sumPhi / energy; - sumEta = sumEta / energy; + sumTheta = sumTheta / energy; int counter = 0; - for (auto entryEta : posEta){ - deltaR += sqrt(pow(entryEta-sumEta,2) + pow(posEta[counter]-sumPhi,2)) * vecEnergy[counter]; + for (auto entryTheta : posTheta){ + deltaR += sqrt(pow(entryTheta-sumTheta,2) + pow(posPhi[counter]-sumPhi,2)) * vecEnergy[counter]; counter++; } cluster.addToShapeParameters(deltaR / energy); @@ -208,7 +215,7 @@ StatusCode CaloTopoClusterFCCee::execute() { clusterWithMixedCells++; posPhi.clear(); - posEta.clear(); + posTheta.clear(); vecEnergy.clear(); } @@ -225,12 +232,19 @@ void CaloTopoClusterFCCee::findingSeeds(const std::unordered_map>& aSeeds) { for (const auto& cell : aCells) { // retrieve the noise const and offset assigned to cell - // // first try to use the cache int system = m_decoder->get(cell.first, m_index_system); if (system == 4) { //ECal barrel int layer = m_decoder_ecal->get(cell.first, m_index_layer_ecal); + double min_threshold = m_min_offset[layer] + m_min_noise[layer] * aNumSigma; + + debug() << "m_min_offset[layer] = " << m_min_offset[layer] << endmsg; + debug() << "m_min_noise[layer] = " << m_min_noise[layer] << endmsg; + debug() << "aNumSigma = " << aNumSigma << endmsg; + debug() << "min_threshold = " << min_threshold << endmsg; + debug() << "abs(cell.second) = " << abs(cell.second) << endmsg; + if (abs(cell.second) < min_threshold) { // if we are below the minimum threshold for the full layer, no need to retrieve the exact value continue; @@ -240,9 +254,10 @@ void CaloTopoClusterFCCee::findingSeeds(const std::unordered_mapnoiseOffset(cell.first) + ( m_noiseTool->noiseRMS(cell.first) * aNumSigma); if (msgLevel() <= MSG::VERBOSE){ - verbose() << "noise offset = " << m_noiseTool->noiseOffset(cell.first) << "GeV " << endmsg; - verbose() << "noise rms = " << m_noiseTool->noiseRMS(cell.first) << "GeV " << endmsg; - verbose() << "seed threshold = " << threshold << "GeV " << endmsg; + debug() << "noise offset = " << m_noiseTool->noiseOffset(cell.first) << "GeV " << endmsg; + debug() << "noise rms = " << m_noiseTool->noiseRMS(cell.first) << "GeV " << endmsg; + debug() << "seed threshold = " << threshold << "GeV " << endmsg; + debug() << "======================================" << endmsg; } if (abs(cell.second) > threshold) { aSeeds.emplace_back(cell.first, cell.second); diff --git a/RecFCCeeCalorimeter/src/components/CaloTopoClusterFCCee.h b/RecFCCeeCalorimeter/src/components/CaloTopoClusterFCCee.h index a524b55d..fa46358b 100644 --- a/RecFCCeeCalorimeter/src/components/CaloTopoClusterFCCee.h +++ b/RecFCCeeCalorimeter/src/components/CaloTopoClusterFCCee.h @@ -126,7 +126,7 @@ class CaloTopoClusterFCCee : public GaudiAlgorithm { //ToolHandle m_cellPositionsHFwdTool{"CellPositionsCaloDiscsTool", this}; /// no segmentation used in HCal - Gaudi::Property m_noSegmentationHCalUsed{this, "noSegmentationHCal", true, "HCal Barrel readout without DD4hep eta-phi segmentation used."}; + Gaudi::Property m_noSegmentationHCalUsed{this, "noSegmentationHCal", true, "HCal Barrel readout without DD4hep theta-module segmentation used."}; /// Seed threshold in sigma Gaudi::Property m_seedSigma{this, "seedSigma", 4, "number of sigma in noise threshold"}; /// Neighbour threshold in sigma @@ -134,7 +134,7 @@ class CaloTopoClusterFCCee : public GaudiAlgorithm { /// Last neighbour threshold in sigma Gaudi::Property m_lastNeighbourSigma{this, "lastNeighbourSigma", 0, "number of sigma in noise threshold"}; /// Name of the electromagnetic calorimeter readout - Gaudi::Property m_readoutName{this, "readoutName", "ECalBarrelPhiEta"}; + Gaudi::Property m_readoutName{this, "readoutName", "ECalBarrelModuleThetaMerged"}; /// General decoder to encode the calorimeter sub-system to determine which positions tool to use dd4hep::DDSegmentation::BitFieldCoder* m_decoder = new dd4hep::DDSegmentation::BitFieldCoder("system:4"); diff --git a/RecFCCeeCalorimeter/src/components/CreateFCCeeCaloNeighbours.cpp b/RecFCCeeCalorimeter/src/components/CreateFCCeeCaloNeighbours.cpp new file mode 100644 index 00000000..fa25906b --- /dev/null +++ b/RecFCCeeCalorimeter/src/components/CreateFCCeeCaloNeighbours.cpp @@ -0,0 +1,489 @@ +#include "CreateFCCeeCaloNeighbours.h" + +#include "DD4hep/Detector.h" +#include "DetCommon/DetUtils.h" +#include "k4Interface/IGeoSvc.h" + +#include "TFile.h" +#include "TTree.h" + +DECLARE_COMPONENT(CreateFCCeeCaloNeighbours) + +CreateFCCeeCaloNeighbours::CreateFCCeeCaloNeighbours(const std::string& aName, ISvcLocator* aSL) + : base_class(aName, aSL) { + declareProperty( "outputFileName", m_outputFileName, "Name of the output file"); +} + +CreateFCCeeCaloNeighbours::~CreateFCCeeCaloNeighbours() {} + +StatusCode CreateFCCeeCaloNeighbours::initialize() { + // Initialize necessary Gaudi components + if (Service::initialize().isFailure()) { + error() << "Unable to initialize Service()" << endmsg; + return StatusCode::FAILURE; + } + m_geoSvc = service("GeoSvc"); + 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; + } + std::unordered_map> map; + + // will be used for volume connecting + int eCalLastLayer; + std::pair extremaECalLastLayerModule; + std::pair extremaECalLastLayerTheta; + double eCalThetaOffset = 0; + double eCalThetaSize = 0; + //double eCalPhiOffset = 0; + double eCalModuleSize = 0; + double hCalThetaOffset = 0; + double hCalThetaSize = 0; + double hCalPhiOffset = 0; + dd4hep::DDSegmentation::BitFieldCoder* decoderECalBarrel = nullptr; + // will be used for volume connecting + std::pair extremaHCalFirstLayerPhi; + std::pair extremaHCalFirstLayerTheta; + std::pair extremaHCalFirstLayerZ; + dd4hep::DDSegmentation::BitFieldCoder* decoderHCalBarrel = nullptr; + + //////////////////////////////////// + /// SEGMENTED THETA-MODULE VOLUMES /// + //////////////////////////////////// + + for (uint iSys = 0; iSys < m_readoutNamesSegmented.size(); iSys++) { + // Check if readouts exist + info() << "Readout: " << m_readoutNamesSegmented[iSys] << endmsg; + if (m_geoSvc->getDetector()->readouts().find(m_readoutNamesSegmented[iSys]) == m_geoSvc->getDetector()->readouts().end()) { + error() << "Readout <<" << m_readoutNamesSegmented[iSys] << ">> does not exist." << endmsg; + return StatusCode::FAILURE; + } + + // get Theta-Module Merged segmentation + dd4hep::DDSegmentation::FCCSWGridModuleThetaMerged* segmentation; + segmentation = dynamic_cast( + m_geoSvc->getDetector()->readout(m_readoutNamesSegmented[iSys]).segmentation().segmentation()); + if (segmentation == nullptr) { + error() << "There is no Theta-Module Merged segmentation!!!!" << endmsg; + return StatusCode::FAILURE; + } + info() << "FCCSWGridModuleThetaMerged: size in Theta " << segmentation->gridSizeTheta() << " , bins in Module " << segmentation->nModules() + << endmsg; + info() << "FCCSWGridModuleThetaMerged: offset in Theta " << segmentation->offsetTheta() << endmsg; + + // retrieve decoders and other info needed for volume (ECal-HCal) connection + auto decoder = m_geoSvc->getDetector()->readout(m_readoutNamesSegmented[iSys]).idSpec().decoder(); + if (m_fieldNamesSegmented[iSys] == "system" && m_fieldValuesSegmented[iSys] == 5) { + decoderECalBarrel = decoder; + eCalThetaSize = segmentation->gridSizeTheta(); + eCalModuleSize = 2 * M_PI / segmentation->nModules(); + //eCalModuleSize = 2 * M_PI / segmentation->phiBins(); + eCalThetaOffset = segmentation->offsetTheta(); + //eCalPhiOffset = segmentation->offsetPhi(); + } + if (m_fieldNamesSegmented[iSys] == "system" && m_fieldValuesSegmented[iSys] == 8) { + decoderHCalBarrel = decoder; + hCalThetaSize = segmentation->gridSizeTheta(); + hCalThetaOffset = segmentation->offsetTheta(); + //hCalPhiOffset = segmentation->offsetPhi(); + } + + // Loop over all cells in the calorimeter and retrieve existing cellIDs + // Loop over active layers + std::vector> extrema; + // extrema[0]: min layer, n layers + extrema.push_back(std::make_pair(0, m_activeVolumesNumbersSegmented[iSys] - 1)); + extrema.push_back(std::make_pair(0, 0)); + extrema.push_back(std::make_pair(0, 0)); + for (unsigned int ilayer = 0; ilayer < m_activeVolumesNumbersSegmented[iSys]; ilayer++) { + dd4hep::DDSegmentation::CellID volumeId = 0; + // Get VolumeID + // m_fieldValuesSegmented: in .py systemValuesModuleTheta = [4] + // m_activeFieldNamesSegmented: in .py activeFieldNamesModuleTheta = ["layer"] + (*decoder)[m_fieldNamesSegmented[iSys]].set(volumeId, m_fieldValuesSegmented[iSys]); + (*decoder)[m_activeFieldNamesSegmented[iSys]].set(volumeId, ilayer); + (*decoder)["theta"].set(volumeId, 0); + (*decoder)["module"].set(volumeId, 0); + //(*decoder)["phi"].set(volumeId, 0); + // Get number of segmentation cells within the active volume + // numberOfCells: return Array of the number of cells in (module, theta) and the minimum theta ID. + auto numCells = det::utils::numberOfCells(volumeId, *segmentation); + // extrema 1: min module number (0), max module number + extrema[1] = std::make_pair(0, (numCells[0] - 1)*segmentation->mergedModules(ilayer)); + // extrema[2]: min theta ID, n (merged) theta cells + extrema[2] = std::make_pair(numCells[2], numCells[2] + (numCells[1] - 1)*segmentation->mergedThetaCells(ilayer)); + + // for layer N-1 of ECal barrel, will be used for volume connecting + // should 5 be systemValuesModuleTheta instead? + if (ilayer == (m_activeVolumesNumbersSegmented[iSys] - 1) && m_fieldNamesSegmented[iSys] == "system" && + m_fieldValuesSegmented[iSys] == 5) { + eCalLastLayer = m_activeVolumesNumbersSegmented[iSys] - 1; + extremaECalLastLayerModule = std::make_pair(0, numCells[0] - 1); + extremaECalLastLayerTheta = std::make_pair(numCells[2], numCells[1] + numCells[2] - 1); + } + else if(m_fieldNamesSegmented[iSys] == "system" && + m_fieldValuesSegmented[iSys] == 8 && m_readoutNamesSegmented[iSys]=="BarHCal_Readout_phitheta"){ +//// + uint cellsTheta = ceil(( 2*m_activeVolumesTheta[ilayer] - segmentation->gridSizeTheta() ) / 2 / segmentation->gridSizeTheta()) * 2 + 1; //ceil( 2*m_activeVolumesRadii[ilayer] / segmentation->gridSizeTheta()); + uint minThetaID = int(floor(( - m_activeVolumesTheta[ilayer] + 0.5 * segmentation->gridSizeTheta() - segmentation->offsetTheta()) / segmentation->gridSizeTheta())); + numCells[1]=cellsTheta; + numCells[2]=minThetaID; + // for layer 0 of HCal barrel, will be used for volume connecting + if (ilayer == 0){ + extremaHCalFirstLayerPhi = std::make_pair(0, numCells[0] - 1); + extremaHCalFirstLayerTheta = std::make_pair(numCells[2], numCells[1] + numCells[2] - 1); + } + } + debug() << "Layer: " << ilayer << endmsg; + debug() << "Extrema[0]: " << extrema[0].first << " , " << extrema[0].second << endmsg; + debug() << "Extrema[1]: " << extrema[1].first << " , " << extrema[1].second << endmsg; + debug() << "Extrema[2]: " << extrema[2].first << " , " << extrema[2].second << endmsg; + debug() << "Number of segmentation cells in (module,theta): " << numCells << endmsg; + // Loop over segmentation cells + for (int imodule = extrema[1].first; imodule <= extrema[1].second; imodule += segmentation->mergedModules(ilayer)) { + for (int itheta = extrema[2].first; itheta <= extrema[2].second; itheta += segmentation->mergedThetaCells(ilayer)) { + dd4hep::DDSegmentation::CellID cellId = volumeId; + decoder->set(cellId, "module", imodule); + decoder->set(cellId, "theta", itheta); // start from the minimum existing theta cell in this layer + uint64_t id = cellId; + map.insert(std::pair>( + id, + det::utils::neighbours_ModuleThetaMerged( + *segmentation, + *decoder, + {m_activeFieldNamesSegmented[iSys], + "module", "theta"}, + extrema, + id, + {false, true, false}, + false))); + } + } + } + if (msgLevel() <= MSG::DEBUG) { + std::vector counter; + counter.assign(40, 0); + for (const auto& item : map) { + counter[item.second.size()]++; + } + for (uint iCount = 0; iCount < counter.size(); iCount++) { + if (counter[iCount] != 0) { + info() << counter[iCount] << " cells have " << iCount << " neighbours" << endmsg; + } + } + } + info() << "total number of cells: " << map.size() << endmsg; + } + + ////////////////////////////////// + /// NESTED VOLUMES /// + ////////////////////////////////// + + for (uint iSys = 0; iSys < m_readoutNamesNested.size(); iSys++) { + // Sanity check + if (m_activeFieldNamesNested.size() != 3) { + error() << "Property activeFieldNamesNested requires 3 names." << endmsg; + return StatusCode::FAILURE; + } + // Check if readouts exist + info() << "Readout: " << m_readoutNamesNested[iSys] << endmsg; + if (m_geoSvc->getDetector()->readouts().find(m_readoutNamesNested[iSys]) == m_geoSvc->getDetector()->readouts().end()) { + error() << "Readout <<" << m_readoutNamesNested[iSys] << ">> does not exist." << endmsg; + return StatusCode::FAILURE; + } + auto decoder = m_geoSvc->getDetector()->readout(m_readoutNamesNested[iSys]).idSpec().decoder(); + // will be used for volume connecting + if (m_fieldNameNested == "system" && m_fieldValuesNested[iSys] == 8) { + decoderHCalBarrel = decoder; + } + //hCalPhiOffset = m_hCalPhiOffset; + // Get VolumeID + dd4hep::DDSegmentation::CellID volumeId = 0; + decoder->set(volumeId, m_fieldNameNested, m_fieldValuesNested[iSys]); + // Get the total number of given hierarchy of active volumes + auto highestVol = gGeoManager->GetTopVolume(); + std::vector numVolumes; + numVolumes.reserve(m_activeVolumeNamesNested.size()); + for (const auto& volName : m_activeVolumeNamesNested) { + numVolumes.push_back(det::utils::countPlacedVolumes(highestVol, volName)); + info() << "Number of active volumes named " << volName << " is " << numVolumes.back() << endmsg; + if (numVolumes.back() == 0) { + error() << "Volume name " << volName << " not found! Check naming in detector description." << endmsg; + return StatusCode::FAILURE; + } + } + // First sort to figure out which volume is inside which one + std::vector> numVolumesMap; + for (unsigned int it = 0; it < m_activeVolumeNamesNested.size(); it++) { + // names of volumes (m_activeVolumeNamesNested) not needed anymore, only corresponding bitfield names are used + // (m_activeFieldNamesNested) + numVolumesMap.push_back(std::pair(m_activeFieldNamesNested[it], numVolumes[it])); + } + std::sort( + numVolumesMap.begin(), numVolumesMap.end(), + [](std::pair vol1, std::pair vol2) { return vol1.second < vol2.second; }); + // now recompute how many volumes exist within the larger volume + for (unsigned int it = numVolumesMap.size() - 1; it > 0; it--) { + if (numVolumesMap[it].second % numVolumesMap[it - 1].second != 0) { + error() << "Given volumes are not nested in each other!" << endmsg; + return StatusCode::FAILURE; + } + numVolumesMap[it].second /= numVolumesMap[it - 1].second; + } + // Debug calculation of total number of cells + if (msgLevel() <= MSG::DEBUG) { + unsigned int checkTotal = 1; + for (const auto& vol : numVolumesMap) { + debug() << "Number of active volumes named " << vol.first << " is " << vol.second << endmsg; + checkTotal *= vol.second; + } + debug() << "Total number of cells ( " << numVolumesMap.back().first << " ) is " << checkTotal << endmsg; + } + // make sure the ordering above is as in property activeFieldNamesNested + std::map activeVolumesNumbersNested; + for (const auto& name : m_activeFieldNamesNested) { + for (const auto& number : numVolumesMap) { + if (name == number.first) { + activeVolumesNumbersNested.insert(std::make_pair(number.first, number.second)); + } + } + } + + // Loop over all cells in the calorimeter and retrieve existing cellIDs + // Loop over active layers + std::vector> extrema; + extrema.push_back(std::make_pair(0, activeVolumesNumbersNested.find(m_activeFieldNamesNested[0])->second - 1)); + extrema.push_back(std::make_pair(0, activeVolumesNumbersNested.find(m_activeFieldNamesNested[1])->second - 1)); + extrema.push_back(std::make_pair(0, activeVolumesNumbersNested.find(m_activeFieldNamesNested[2])->second - 1)); + // for layer 0 of HCal barrel + if (m_fieldNameNested == "system" && m_fieldValuesNested[iSys] == 8) { + extremaHCalFirstLayerPhi = + std::make_pair(0, activeVolumesNumbersNested.find(m_activeFieldNamesNested[1])->second - 1); + extremaHCalFirstLayerZ = + std::make_pair(0, activeVolumesNumbersNested.find(m_activeFieldNamesNested[2])->second - 1); + } + for (unsigned int ilayer = 0; ilayer < activeVolumesNumbersNested.find(m_activeFieldNamesNested[0])->second; + ilayer++) { + for (unsigned int iphi = 0; iphi < activeVolumesNumbersNested.find(m_activeFieldNamesNested[1])->second; iphi++) { + for (unsigned int iz = 0; iz < activeVolumesNumbersNested.find(m_activeFieldNamesNested[2])->second; iz++) { + + dd4hep::DDSegmentation::CellID cID = volumeId; + decoder->set(cID, m_activeFieldNamesNested[0], ilayer); + decoder->set(cID, m_activeFieldNamesNested[1], iphi); + decoder->set(cID, m_activeFieldNamesNested[2], iz); + + map.insert(std::pair>( + cID, det::utils::neighbours(*decoder, {m_activeFieldNamesNested[0], m_activeFieldNamesNested[1], + m_activeFieldNamesNested[2]}, + extrema, cID, {false, true, false}, true))); + } + } + } + if (msgLevel() <= MSG::DEBUG) { + std::vector counter; + counter.assign(40, 0); + for (const auto& item : map) { + counter[item.second.size()]++; + } + for (uint iCount = 0; iCount < counter.size(); iCount++) { + if (counter[iCount] != 0) { + info() << counter[iCount] << " cells have " << iCount << " neighbours" << endmsg; + } + } + } + } + + ////////////////////////////////////////////////// + /// BARREL: connection ECAL + HCAL /// + ///////////////////////////////////////////////// + int count=0; +/* + if (m_connectBarrels) { + // first check if ECAL barrel (system==5) and HCal barrel (system==8) are configured + if (decoderECalBarrel == nullptr || decoderHCalBarrel == nullptr) { + error() << "ECAL barrel and/or HCal barrel are not configured correctly!" << endmsg; + return StatusCode::FAILURE; + } + // print how many cells in each dimensions will be matched + if(m_readoutNamesNested.size()!=0){ + info() << "ECAL layer " << eCalLastLayer << " is a neighbour of HCAL layer 0." << endmsg; + info() << "ECAL phi cells " << extremaECalLastLayerModule.first << " - " << extremaECalLastLayerModule.second + << " will be matched to HCAL " << m_activeFieldNamesNested[1] << "(s) " << extremaHCalFirstLayerPhi.first + << " - " << extremaHCalFirstLayerPhi.second << endmsg; + info() << "ECAL theta cells " << extremaECalLastLayerTheta.first << " - " << extremaECalLastLayerTheta.second + << " will be matched to HCAL " << m_activeFieldNamesNested[2] << "(s) " << extremaHCalFirstLayerZ.first + << " - " << extremaHCalFirstLayerZ.second << endmsg; + } + else{ + info() << "ECAL layer " << eCalLastLayer << " is a neighbour of HCAL layer 0." << endmsg; + info() << "ECAL phi cells " << extremaECalLastLayerModule.first << " - " << extremaECalLastLayerModule.second + << " will be matched to HCAL cells " << extremaHCalFirstLayerPhi.first + << " - " << extremaHCalFirstLayerPhi.second << endmsg; + info() << "ECAL theta cells " << extremaECalLastLayerTheta.first << " - " << extremaECalLastLayerTheta.second + << " will be matched to HCAL " << extremaHCalFirstLayerTheta.first + << " - " << extremaHCalFirstLayerTheta.second << endmsg; + } + + std::unordered_map> thetaNeighbours; + std::unordered_map> phiNeighbours; + double hCalPhiSize = 2 * M_PI / (extremaHCalFirstLayerPhi.second - extremaHCalFirstLayerPhi.first + 1); + // loop over z and find which theta cells to add + if (m_readoutNamesNested.size()!=0){ + for (int iZ = 0; iZ < extremaHCalFirstLayerZ.second + 1; iZ++) { + double lowZ = m_hCalZOffset + iZ * m_hCalZSize; + double highZ = m_hCalZOffset + (iZ + 1) * m_hCalZSize; + double lowTheta = 0, highTheta = 0; + if (fabs(lowZ) < 1e-3) { + lowTheta = 0; + } else { + lowTheta = +////TODO + lowZ / fabs(lowZ) * atan(m_hCalRinner / lowZ); + //lowZ / fabs(lowZ) * (-log(fabs(tan(atan(m_hCalRinner / lowZ) / 2.)))); // theta = atan(m_hCalRinner / lowZ) + } + if (fabs(highZ) < 1e-3) { + highTheta = 0; + } else { + highTheta = highZ / fabs(highZ) * (-log(fabs(tan(atan(m_hCalRinner / highZ) / 2.)))); + } + debug() << "HCal z id : " << iZ << endmsg; + debug() << "HCal theta range : " << lowTheta << " - " << highTheta << endmsg; + int lowId = floor((lowTheta - 0.5 * eCalThetaSize - eCalThetaOffset) / eCalThetaSize); + int highId = floor((highTheta + 0.5 * eCalThetaSize - eCalThetaOffset) / eCalThetaSize); + debug() << "ECal theta range : " << lowId * eCalThetaSize + eCalThetaOffset << " - " + << highId * eCalThetaSize + eCalThetaOffset << endmsg; + std::vector neighbours; + for (int idThetaToAdd = lowId; idThetaToAdd <= highId; idThetaToAdd++) { + if (idThetaToAdd >= extremaECalLastLayerTheta.first && idThetaToAdd <= extremaECalLastLayerTheta.second) { + neighbours.push_back(idThetaToAdd); + } + } + debug() << "HCal z id : " << iZ << endmsg; + debug() << "Found ECal Neighbours in theta : " << neighbours.size() << endmsg; + for (auto id : neighbours) { + debug() << "ECal Neighbours id : " << id << endmsg; + } + thetaNeighbours.insert(std::pair>(iZ, neighbours)); + } + } + else{ // loop over theta cells to match in theta + for (int iTheta = extremaHCalFirstLayerTheta.first; iTheta < extremaHCalFirstLayerTheta.second + 1; iTheta++) { + double lowTheta = hCalThetaOffset + iTheta * hCalThetaSize; + double highTheta = hCalThetaOffset + (iTheta + 1) * hCalThetaSize; + debug() << "HCal theta range : " << lowTheta << " - " << highTheta << endmsg; + int lowId = floor((lowTheta - 0.5 * eCalThetaSize - eCalThetaOffset) / eCalThetaSize); + int highId = floor((highTheta + 0.5 * eCalThetaSize - eCalThetaOffset) / eCalThetaSize); + debug() << "ECal theta range : " << lowId * eCalThetaSize + eCalThetaOffset << " - " + << highId * eCalThetaSize + eCalThetaOffset << endmsg; + std::vector neighbours; + for (int idThetaToAdd = lowId; idThetaToAdd <= highId; idThetaToAdd++) { + neighbours.push_back(det::utils::cyclicNeighbour(idThetaToAdd, extremaECalLastLayerTheta)); + } + debug() << "HCal theta id : " << iTheta << endmsg; + debug() << "Found ECal Neighbours in theta : " << neighbours.size() << endmsg; + for (auto id : neighbours) { + debug() << "ECal Neighbours id : " << id << endmsg; + } + thetaNeighbours.insert(std::pair>(iTheta, neighbours)); + } + } + // loop over phi and find which phi cells to add + for (int iPhi = 0; iPhi < extremaHCalFirstLayerPhi.second +1; iPhi++) { + double lowPhi = hCalPhiOffset + iPhi * hCalPhiSize; + double highPhi = hCalPhiOffset + (iPhi + 1) * hCalPhiSize; + debug() << "HCal phi range : " << lowPhi << " - " << highPhi << endmsg; + int lowId = floor((lowPhi - 0.5 * eCalModuleSize - eCalPhiOffset) / eCalModuleSize); + int highId = floor((highPhi + 0.5 * eCalModuleSize - eCalPhiOffset) / eCalModuleSize); + debug() << "ECal phi range : " << lowId * eCalModuleSize + eCalPhiOffset << " - " + << highId * eCalModuleSize + eCalPhiOffset << endmsg; + std::vector neighbours; + for (int idPhiToAdd = lowId; idPhiToAdd <= highId; idPhiToAdd++) { + neighbours.push_back(det::utils::cyclicNeighbour(idPhiToAdd, extremaECalLastLayerModule)); + } + debug() << "HCal phi id : " << iPhi << endmsg; + debug() << "Found ECal Neighbours in phi : " << neighbours.size() << endmsg; + for (auto id : neighbours) { + debug() << "ECal Neighbours id : " << id << endmsg; + } + phiNeighbours.insert(std::pair>(iPhi, neighbours)); + } + // add neighbours to both ecal cell and hcal cells + dd4hep::DDSegmentation::CellID ecalCellId = 0; + dd4hep::DDSegmentation::CellID hcalCellId = 0; + (*decoderECalBarrel)["system"].set(ecalCellId, 5); + (*decoderECalBarrel)[m_activeFieldNamesSegmented[0]].set(ecalCellId, eCalLastLayer); + (*decoderHCalBarrel)["system"].set(hcalCellId, 8); + // loop over nested hcal cells + if (m_readoutNamesNested.size()!=0){ + (*decoderHCalBarrel)[m_activeFieldNamesNested[0]].set(hcalCellId, 0); + for (auto iZ : thetaNeighbours) { + (*decoderHCalBarrel)[m_activeFieldNamesNested[2]].set(hcalCellId, iZ.first); + for (auto iMod : phiNeighbours) { + (*decoderHCalBarrel)[m_activeFieldNamesNested[1]].set(hcalCellId, iMod.first); + for (auto iTheta : iZ.second) { + (*decoderECalBarrel)["theta"].set(ecalCellId, iTheta); + for (auto iPhi : iMod.second) { + (*decoderECalBarrel)["phi"].set(ecalCellId, iPhi); + map.find(hcalCellId)->second.push_back(ecalCellId); + map.find(ecalCellId)->second.push_back(hcalCellId); + count++; + } + } + } + } + } + // loop over segmented hcal cells + else { + (*decoderHCalBarrel)[m_activeFieldNamesSegmented[1]].set(hcalCellId, 0); + for (auto iThetaHCal : thetaNeighbours) { + (*decoderHCalBarrel)["theta"].set(hcalCellId, iThetaHCal.first); + for (auto iPhiHCal : phiNeighbours) { + (*decoderHCalBarrel)["phi"].set(hcalCellId, iPhiHCal.first); + for (auto iTheta : iThetaHCal.second) { + (*decoderECalBarrel)["theta"].set(ecalCellId, iTheta); + for (auto iPhi : iPhiHCal.second) { + (*decoderECalBarrel)["phi"].set(ecalCellId, iPhi); + map.find(hcalCellId)->second.push_back(ecalCellId); + map.find(ecalCellId)->second.push_back(hcalCellId); + count ++; + } + } + } + } + } + } +*/ + if (msgLevel() <= MSG::DEBUG) { + std::vector counter; + counter.assign(40, 0); + for (const auto& item : map) { + counter[item.second.size()]++; + } + for (uint iCount = 0; iCount < counter.size(); iCount++) { + if (counter[iCount] != 0) { + debug() << counter[iCount] << " cells have " << iCount << " neighbours" << endmsg; + } + } + } + debug() << "cells with neighbours across Calo boundaries: " << count << endmsg; + + std::unique_ptr file(TFile::Open(m_outputFileName.c_str(), "RECREATE")); + file->cd(); + TTree tree("neighbours", "Tree with map of neighbours"); + uint64_t saveCellId; + std::vector saveNeighbours; + tree.Branch("cellId", &saveCellId, "cellId/l"); + tree.Branch("neighbours", &saveNeighbours); + for (const auto& item : map) { + saveCellId = item.first; + saveNeighbours = item.second; + tree.Fill(); + } + file->Write(); + file->Close(); + + return StatusCode::SUCCESS; +} + +StatusCode CreateFCCeeCaloNeighbours::finalize() { return Service::finalize(); } diff --git a/RecFCCeeCalorimeter/src/components/CreateFCCeeCaloNeighbours.h b/RecFCCeeCalorimeter/src/components/CreateFCCeeCaloNeighbours.h new file mode 100644 index 00000000..199816fa --- /dev/null +++ b/RecFCCeeCalorimeter/src/components/CreateFCCeeCaloNeighbours.h @@ -0,0 +1,93 @@ +#ifndef RECCALORIMETER_CREATEFCCEECALONEIGHBOURS_H +#define RECCALORIMETER_CREATEFCCEECALONEIGHBOURS_H + +// Gaudi +#include "GaudiKernel/Service.h" +#include "k4Interface/ICaloCreateMap.h" + +#include "DetCommon/DetUtils.h" +#include "k4Interface/IGeoSvc.h" +#include "DetSegmentation/FCCSWGridPhiEta.h" +#include "DetSegmentation/FCCSWGridPhiTheta.h" +#include "DetSegmentation/FCCSWGridModuleThetaMerged.h" +#include "k4FWCore/DataHandle.h" +#include "k4Interface/ICellPositionsTool.h" + +#include "DD4hep/Readout.h" +#include "DD4hep/Volumes.h" +#include "DDSegmentation/Segmentation.h" +#include "TGeoManager.h" + +class IGeoSvc; + +/** @class CreateFCCeeCaloNeighbours + * + * Service building a map of neighbours for all existing cells in the geometry. + * The volumes for which the neighbour map is created can be either segmented in theta-module (e.g. ECal inclined), + * or can contain nested volumes (e.g. HCal barrel). + * + * @author Anna Zaborowska + */ + +class CreateFCCeeCaloNeighbours : public extends1 { +public: + /// Standard constructor + explicit CreateFCCeeCaloNeighbours(const std::string& aName, ISvcLocator* aSL); + /// Standard destructor + virtual ~CreateFCCeeCaloNeighbours(); + /** Initialize the map creator service. + * @return status code + */ + virtual StatusCode initialize() final; + /** Finalize the map creator service. + * @return status code + */ + virtual StatusCode finalize() final; + +private: + /// Pointer to the geometry service + SmartIF m_geoSvc; + + /// Names of the detector readout for volumes with theta-module segmentation + Gaudi::Property> m_readoutNamesSegmented{this, "readoutNamesModuleTheta", {"ECalBarrelModuleThetaMerged"}}; + /// Name of the fields describing the segmented volume + Gaudi::Property> m_fieldNamesSegmented{this, "systemNamesModuleTheta", {"system"}}; + /// Values of the fields describing the segmented volume + Gaudi::Property> m_fieldValuesSegmented{this, "systemValuesModuleTheta", {5}}; + /// Names of the active volume in geometry along radial axis (e.g. layer), the others are "module", "theta" + Gaudi::Property> m_activeFieldNamesSegmented{this, "activeFieldNamesModuleTheta", {"layer"}}; + /// Number of layers in the segmented volume + Gaudi::Property> m_activeVolumesNumbersSegmented{this, "activeVolumesNumbers", {12}}; + // Radii of layers in the segmented volume + Gaudi::Property> m_activeVolumesTheta{this, "activeVolumesTheta"}; + + /// Names of the detector readout for volumes with nested volume structure and no segmentation + Gaudi::Property> m_readoutNamesNested{this, "readoutNamesVolumes"}; + /// Name of the field describing the nested volume + Gaudi::Property m_fieldNameNested{this, "systemNameNested"}; + /// Values of the fields describing the nested volume + Gaudi::Property> m_fieldValuesNested{this, "systemValuesNested"}; + /// Names of the active volume in geometry: along radial axis, azimuthal angle, and along z axis + Gaudi::Property> m_activeFieldNamesNested{ + this, "activeFieldNamesNested"}; + /// Names of the nested volumes - to retrieve the number of active volumes, need to correspond to m_activeFieldNamesNested + Gaudi::Property> m_activeVolumeNamesNested{ + this, + "activeVolumeNamesNested", + {"layerVolume", "moduleVolume", "wedgeVolume"}}; // to find out number of volumes + /// Name of output file + std::string m_outputFileName; + + // For combination of barrels: flag if ECal and HCal barrels should be merged + Gaudi::Property m_connectBarrels{this, "connectBarrels", true}; + // For combination of barrels: size of HCal cell along z axis + Gaudi::Property m_hCalZSize{this, "hCalZsize", 18}; + // For combination of barrels: offset of HCal detector in z (lower edge) + Gaudi::Property m_hCalZOffset{this, "hCalZoffset", -4590}; + // For combination of barrels: HCal inner radius for calculation of eta from z ??? + Gaudi::Property m_hCalRinner{this, "hCalRinner", 2850}; + // For combination of barrels: offset of HCal modules in phi (lower edge) + Gaudi::Property m_hCalPhiOffset{this, "hCalPhiOffset"}; +}; + +#endif /* RECALORIMETER_CREATEFCCHHCALONEIGHBOURS_H */ diff --git a/RecFCCeeCalorimeter/src/components/CreateFCCeeCaloNoiseLevelMap.cpp b/RecFCCeeCalorimeter/src/components/CreateFCCeeCaloNoiseLevelMap.cpp new file mode 100644 index 00000000..fcbe940d --- /dev/null +++ b/RecFCCeeCalorimeter/src/components/CreateFCCeeCaloNoiseLevelMap.cpp @@ -0,0 +1,231 @@ +#include "CreateFCCeeCaloNoiseLevelMap.h" + +#include "DD4hep/Detector.h" +#include "DetCommon/DetUtils.h" +#include "k4Interface/IGeoSvc.h" +#include "DetSegmentation/FCCSWGridModuleThetaMerged.h" + +#include "TFile.h" +#include "TTree.h" + +DECLARE_COMPONENT(CreateFCCeeCaloNoiseLevelMap) + +CreateFCCeeCaloNoiseLevelMap::CreateFCCeeCaloNoiseLevelMap(const std::string& aName, ISvcLocator* aSL) + : base_class(aName, aSL) { + declareProperty("ECalBarrelNoiseTool", m_ecalBarrelNoiseTool, "Handle for the cells noise tool of Barrel ECal"); + declareProperty("HCalBarrelNoiseTool", m_hcalBarrelNoiseTool, "Handle for the cells noise tool of Barrel HCal"); + declareProperty( "outputFileName", m_outputFileName, "Name of the output file"); +} + +CreateFCCeeCaloNoiseLevelMap::~CreateFCCeeCaloNoiseLevelMap() {} + +StatusCode CreateFCCeeCaloNoiseLevelMap::initialize() { + // Initialize necessary Gaudi components + if (Service::initialize().isFailure()) { + error() << "Unable to initialize Service()" << endmsg; + return StatusCode::FAILURE; + } + m_geoSvc = service("GeoSvc"); + 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; + } + + //m_ecalBarrelNoiseTool.retrieve(); + //m_hcalBarrelNoiseTool.retrieve(); + + std::unordered_map> map; + + /////////////////////////////////////// + /// SEGMENTED MODULE-THETA VOLUMES /// + /////////////////////////////////////// + + for (uint iSys = 0; iSys < m_readoutNamesSegmented.size(); iSys++) { + // Check if readouts exist + info() << "Readout: " << m_readoutNamesSegmented[iSys] << endmsg; + if (m_geoSvc->getDetector()->readouts().find(m_readoutNamesSegmented[iSys]) == m_geoSvc->getDetector()->readouts().end()) { + error() << "Readout <<" << m_readoutNamesSegmented[iSys] << ">> does not exist." << endmsg; + return StatusCode::FAILURE; + } + // get segmentation + dd4hep::DDSegmentation::FCCSWGridModuleThetaMerged* segmentation; + segmentation = dynamic_cast( + m_geoSvc->getDetector()->readout(m_readoutNamesSegmented[iSys]).segmentation().segmentation()); + if (segmentation == nullptr) { + error() << "There is no Module-Theta segmentation!!!!" << endmsg; + return StatusCode::FAILURE; + } + + info() << "FCCSWGridModuleThetaMerged: size in Theta " << segmentation->gridSizeTheta() << " , bins in Module " << segmentation->nModules() + << endmsg; + info() << "FCCSWGridModuleThetaMerged: offset in Theta " << segmentation->offsetTheta() << endmsg; + + auto decoder = m_geoSvc->getDetector()->readout(m_readoutNamesSegmented[iSys]).idSpec().decoder(); + // Loop over all cells in the calorimeter and retrieve existing cellIDs + // Loop over active layers + std::vector> extrema; + extrema.push_back(std::make_pair(0, m_activeVolumesNumbersSegmented[iSys] - 1)); + extrema.push_back(std::make_pair(0, 0)); + extrema.push_back(std::make_pair(0, 0)); + for (unsigned int ilayer = 0; ilayer < m_activeVolumesNumbersSegmented[iSys]; ilayer++) { + dd4hep::DDSegmentation::CellID volumeId = 0; + // Get VolumeID + (*decoder)[m_fieldNamesSegmented[iSys]].set(volumeId, m_fieldValuesSegmented[iSys]); + (*decoder)[m_activeFieldNamesSegmented[iSys]].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, *segmentation); + // Range of module ID + extrema[1] = std::make_pair(0, (numCells[0] - 1)*segmentation->mergedModules(ilayer)); + // Range and min of theta ID + // for HCAL use alternative calculation + if(m_fieldNamesSegmented[iSys] == "system" && + m_fieldValuesSegmented[iSys] == m_hcalBarrelSysId){ + uint cellsTheta = ceil(( 2*m_activeVolumesTheta[ilayer] - segmentation->gridSizeTheta() ) / 2 / segmentation->gridSizeTheta()) * 2 + 1; //ceil( 2*m_activeVolumesRadii[ilayer] / segmentation->gridSizeTheta()); + uint minThetaID = int(floor(( - m_activeVolumesTheta[ilayer] + 0.5 * segmentation->gridSizeTheta() - segmentation->offsetTheta()) / segmentation->gridSizeTheta())); + numCells[1]=cellsTheta; + numCells[2]=minThetaID; + } + extrema[2] = std::make_pair(numCells[2], numCells[2] + (numCells[1] - 1)*segmentation->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*segmentation->mergedModules(ilayer)); + decoder->set(cellId, "theta", numCells[2] + itheta*segmentation->mergedThetaCells(ilayer)); // start from the minimum existing theta cell in this layer + uint64_t id = cellId; + double noise = 0.; + double noiseOffset = 0.; + if (m_fieldValuesSegmented[iSys] == m_hcalBarrelSysId){ + noise = m_hcalBarrelNoiseTool->getNoiseConstantPerCell(id); + noiseOffset = m_hcalBarrelNoiseTool->getNoiseOffsetPerCell(id); + } else if (m_fieldValuesSegmented[iSys] == m_ecalBarrelSysId){ + noise = m_ecalBarrelNoiseTool->getNoiseConstantPerCell(id); + noiseOffset = m_ecalBarrelNoiseTool->getNoiseOffsetPerCell(id); + } + map.insert( std::pair >(id, std::make_pair(noise, noiseOffset) ) ); + } + } + } + } // end of SEGMENTED MODULE-THETA VOLUMES + + ////////////////////////////////// + /// NESTED VOLUMES /// + ////////////////////////////////// + + for (uint iSys = 0; iSys < m_readoutNamesNested.size(); iSys++) { + // Sanity check + if (m_activeFieldNamesNested.size() != 3) { + error() << "Property activeFieldNamesNested requires 3 names." << endmsg; + return StatusCode::FAILURE; + } + // Check if readouts exist + info() << "Readout: " << m_readoutNamesNested[iSys] << endmsg; + if (m_geoSvc->getDetector()->readouts().find(m_readoutNamesNested[iSys]) == m_geoSvc->getDetector()->readouts().end()) { + error() << "Readout <<" << m_readoutNamesNested[iSys] << ">> does not exist." << endmsg; + return StatusCode::FAILURE; + } + auto decoder = m_geoSvc->getDetector()->readout(m_readoutNamesNested[iSys]).idSpec().decoder(); + // Get VolumeID + dd4hep::DDSegmentation::CellID volumeId = 0; + decoder->set(volumeId, m_fieldNameNested, m_fieldValuesNested[iSys]); + // Get the total number of given hierarchy of active volumes + auto highestVol = gGeoManager->GetTopVolume(); + std::vector numVolumes; + numVolumes.reserve(m_activeVolumeNamesNested.size()); + for (const auto& volName : m_activeVolumeNamesNested) { + numVolumes.push_back(det::utils::countPlacedVolumes(highestVol, volName)); + info() << "Number of active volumes named " << volName << " is " << numVolumes.back() << endmsg; + if (numVolumes.back() == 0) { + error() << "Volume name " << volName << " not found! Check naming in detector description." << endmsg; + return StatusCode::FAILURE; + } + } + // First sort to figure out which volume is inside which one + std::vector> numVolumesMap; + for (unsigned int it = 0; it < m_activeVolumeNamesNested.size(); it++) { + // names of volumes (m_activeVolumeNamesNested) not needed anymore, only corresponding bitfield names are used + // (m_activeFieldNamesNested) + numVolumesMap.push_back(std::pair(m_activeFieldNamesNested[it], numVolumes[it])); + } + std::sort( + numVolumesMap.begin(), numVolumesMap.end(), + [](std::pair vol1, std::pair vol2) { return vol1.second < vol2.second; }); + // now recompute how many volumes exist within the larger volume + for (unsigned int it = numVolumesMap.size() - 1; it > 0; it--) { + if (numVolumesMap[it].second % numVolumesMap[it - 1].second != 0) { + error() << "Given volumes are not nested in each other!" << endmsg; + return StatusCode::FAILURE; + } + numVolumesMap[it].second /= numVolumesMap[it - 1].second; + } + // Debug calculation of total number of cells + if (msgLevel() <= MSG::DEBUG) { + unsigned int checkTotal = 1; + for (const auto& vol : numVolumesMap) { + debug() << "Number of active volumes named " << vol.first << " is " << vol.second << endmsg; + checkTotal *= vol.second; + } + debug() << "Total number of cells ( " << numVolumesMap.back().first << " ) is " << checkTotal << endmsg; + } + // make sure the ordering above is as in property activeFieldNamesNested + std::map activeVolumesNumbersNested; + for (const auto& name : m_activeFieldNamesNested) { + for (const auto& number : numVolumesMap) { + if (name == number.first) { + activeVolumesNumbersNested.insert(std::make_pair(number.first, number.second)); + } + } + } + + // Loop over all cells in the calorimeter and retrieve existing cellIDs + // Loop over active layers + std::vector> extrema; + extrema.push_back(std::make_pair(0, activeVolumesNumbersNested.find(m_activeFieldNamesNested[0])->second - 1)); + extrema.push_back(std::make_pair(0, activeVolumesNumbersNested.find(m_activeFieldNamesNested[1])->second - 1)); + extrema.push_back(std::make_pair(0, activeVolumesNumbersNested.find(m_activeFieldNamesNested[2])->second - 1)); + for (unsigned int ilayer = 0; ilayer < activeVolumesNumbersNested.find(m_activeFieldNamesNested[0])->second; + ilayer++) { + for (unsigned int iphi = 0; iphi < activeVolumesNumbersNested.find(m_activeFieldNamesNested[1])->second; iphi++) { + for (unsigned int iz = 0; iz < activeVolumesNumbersNested.find(m_activeFieldNamesNested[2])->second; iz++) { + + dd4hep::DDSegmentation::CellID cID = volumeId; + decoder->set(cID, m_activeFieldNamesNested[0], ilayer); + decoder->set(cID, m_activeFieldNamesNested[1], iphi); + decoder->set(cID, m_activeFieldNamesNested[2], iz); + + double noise = m_hcalBarrelNoiseTool->getNoiseConstantPerCell(cID); + double noiseOffset = m_hcalBarrelNoiseTool->getNoiseOffsetPerCell(cID); + + map.insert( std::pair >(cID, std::make_pair(noise, noiseOffset) ) ); + } + } + } + } // end of NESTED VOLUMES + + std::unique_ptr file(TFile::Open(m_outputFileName.c_str(), "RECREATE")); + file->cd(); + TTree tree("noisyCells", "Tree with map of noise per cell"); + uint64_t saveCellId; + double saveNoiseLevel; + double saveNoiseOffset; + tree.Branch("cellId", &saveCellId, "cellId/l"); + tree.Branch("noiseLevel", &saveNoiseLevel); + tree.Branch("noiseOffset", &saveNoiseOffset); + for (const auto& item : map) { + saveCellId = item.first; + saveNoiseLevel = item.second.first; + saveNoiseOffset = item.second.second; + tree.Fill(); + } + file->Write(); + file->Close(); + + return StatusCode::SUCCESS; +} + +StatusCode CreateFCCeeCaloNoiseLevelMap::finalize() { return Service::finalize(); } diff --git a/RecFCCeeCalorimeter/src/components/CreateFCCeeCaloNoiseLevelMap.h b/RecFCCeeCalorimeter/src/components/CreateFCCeeCaloNoiseLevelMap.h new file mode 100644 index 00000000..da65dbfb --- /dev/null +++ b/RecFCCeeCalorimeter/src/components/CreateFCCeeCaloNoiseLevelMap.h @@ -0,0 +1,80 @@ +#ifndef RECCALORIMETER_CREATEFCCEECALONOISELEVELMAP_H +#define RECCALORIMETER_CREATEFCCEECALONOISELEVELMAP_H + +// Gaudi +#include "GaudiKernel/Service.h" +#include "k4Interface/ICaloCreateMap.h" +#include "k4Interface/INoiseConstTool.h" +#include "k4Interface/ICellPositionsTool.h" + +class IGeoSvc; + +/** @class CreateFCCeeCaloNoiseLevelMap + * + * Service building a map from cellIds to noise level per cell. + * The volumes for which the neighbour map is created can be either segmented in Module-Theta (e.g. ECal inclined), + * or can contain nested volumes (e.g. HCal barrel). + * + * @author Coralie Neubueser + */ + +class CreateFCCeeCaloNoiseLevelMap : public extends1 { +public: + /// Standard constructor + explicit CreateFCCeeCaloNoiseLevelMap(const std::string& aName, ISvcLocator* aSL); + /// Standard destructor + virtual ~CreateFCCeeCaloNoiseLevelMap(); + /** Initialize the map creator service. + * @return status code + */ + virtual StatusCode initialize() final; + /** Finalize the map creator service. + * @return status code + */ + virtual StatusCode finalize() final; + +private: + /// Pointer to the geometry service + SmartIF m_geoSvc; + + /// Handle for the cells noise tool in ECal + ToolHandle m_ecalBarrelNoiseTool{"ReadNoiseFromFileTool", this}; + Gaudi::Property m_ecalBarrelSysId{this, "ecalBarrelSysId", 5}; + /// Handle for the cells noise tool in HCal + ToolHandle m_hcalBarrelNoiseTool{"ReadNoiseFromFileTool", this}; + Gaudi::Property m_hcalBarrelSysId{this, "hcalBarrelSysId", 8}; + + /// Names of the detector readout for volumes with Module-Theta segmentation + Gaudi::Property> m_readoutNamesSegmented{this, "readoutNamesModuleTheta", {"ECalBarrelModuleThetaMerged"}}; + /// Name of the fields describing the segmented volume + Gaudi::Property> m_fieldNamesSegmented{this, "systemNamesModuleTheta", {"system"}}; + /// Values of the fields describing the segmented volume + Gaudi::Property> m_fieldValuesSegmented{this, "systemValuesModuleTheta", {5}}; + /// Names of the active volume in geometry along radial axis (e.g. layer), the others are "module", "theta" + Gaudi::Property> m_activeFieldNamesSegmented{this, "activeFieldNamesModuleTheta", {"layer"}}; + /// Number of layers in the segmented volume + Gaudi::Property> m_activeVolumesNumbersSegmented{this, "activeVolumesNumbers", {8}}; + // Radii of layers in the segmented volume + Gaudi::Property> m_activeVolumesTheta{this, "activeVolumesTheta"}; + + /// Names of the detector readout for volumes with nested volume structure and no segmentation + Gaudi::Property> m_readoutNamesNested{this, "readoutNamesVolumes", {"HCalBarrelReadout"}}; + /// Name of the field describing the nested volume + Gaudi::Property m_fieldNameNested{this, "systemNameNested", "system"}; + /// Values of the fields describing the nested volume + Gaudi::Property> m_fieldValuesNested{this, "systemValuesNested", {8}}; + /// Names of the active volume in geometry: along radial axis, azimuthal angle, and along z axis + Gaudi::Property> m_activeFieldNamesNested{ + this, "activeFieldNamesNested", {"layer", "module", "row"}}; + /// Names of the nested volumes - to retrieve the number of active volumes, need to correspond to + /// m_activeFieldNamesNested + Gaudi::Property> m_activeVolumeNamesNested{ + this, + "activeVolumeNamesNested", + {"layerVolume", "moduleVolume", "wedgeVolume"}}; // to find out number of volumes + + /// Name of output file + std::string m_outputFileName; +}; + +#endif /* RECALORIMETER_CREATEFCCEECALONOISELEVELMAP_H */