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..52a07d44 100644 --- a/RecFCCeeCalorimeter/src/components/CaloTopoClusterFCCee.h +++ b/RecFCCeeCalorimeter/src/components/CaloTopoClusterFCCee.h @@ -31,8 +31,7 @@ class Segmentation; } } -/** @class CaloTopoClusterFCCeeAlgorithm Reconstruction/RecFCCeeCalorimeter/src/components/CaloTopoClusterFCCee.h - * CombinedCaloTopoCluster.h +/** @class CaloTopoClusterFCCee k4RecCalorimeter/RecFCCeeCalorimeter/src/components/CaloTopoClusterFCCee.h * * Algorithm building the topological clusters for the energy reconstruction, following ATLAS note * ATL-LARG-PUB-2008-002. @@ -126,7 +125,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 +133,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..8e8eabbe --- /dev/null +++ b/RecFCCeeCalorimeter/src/components/CreateFCCeeCaloNeighbours.cpp @@ -0,0 +1,488 @@ +#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); + // 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, + m_includeDiagonalCells + ))); + } + } + } + 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..07bed908 --- /dev/null +++ b/RecFCCeeCalorimeter/src/components/CreateFCCeeCaloNeighbours.h @@ -0,0 +1,95 @@ +#ifndef RECFCCEECALORIMETER_CREATEFCCEECALONEIGHBOURS_H +#define RECFCCEECALORIMETER_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 Giovanni Marchiori + */ + +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"}; + /// Whether to create the geant4 geometry or not + Gaudi::Property m_includeDiagonalCells{this, "includeDiagonalCells", false, "If True will consider also diagonal neighbours in volumes with theta-module segmentation"}; + + /// 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 /* RECFCCEECALORIMETER_CREATEFCCEECALONEIGHBOURS_H */ diff --git a/RecFCCeeCalorimeter/src/components/CreateFCCeeCaloNoiseLevelMap.cpp b/RecFCCeeCalorimeter/src/components/CreateFCCeeCaloNoiseLevelMap.cpp new file mode 100644 index 00000000..deaf2261 --- /dev/null +++ b/RecFCCeeCalorimeter/src/components/CreateFCCeeCaloNoiseLevelMap.cpp @@ -0,0 +1,238 @@ +#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; + } + + if (!m_ecalBarrelNoiseTool.retrieve()) { + error() << "Unable to retrieve the ECAL noise tool!!!" << endmsg; + return StatusCode::FAILURE; + } + + if (!m_hcalBarrelNoiseTool.retrieve()) { + error() << "Unable to retrieve the HCAL noise tool!!!" << endmsg; + return StatusCode::FAILURE; + } + + 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 */ diff --git a/RecFCCeeCalorimeter/src/components/ReadNoiseVsThetaFromFileTool.cpp b/RecFCCeeCalorimeter/src/components/ReadNoiseVsThetaFromFileTool.cpp new file mode 100644 index 00000000..3ae2c2e7 --- /dev/null +++ b/RecFCCeeCalorimeter/src/components/ReadNoiseVsThetaFromFileTool.cpp @@ -0,0 +1,259 @@ +#include "ReadNoiseVsThetaFromFileTool.h" + +// FCCSW +#include "DetCommon/DetUtils.h" +#include "k4Interface/IGeoSvc.h" +#include "DDSegmentation/Segmentation.h" + +// DD4hep +#include "DD4hep/Detector.h" +#include "DD4hep/Readout.h" + +// Root +#include "TFile.h" +#include "TH1F.h" + +DECLARE_COMPONENT(ReadNoiseVsThetaFromFileTool) + +ReadNoiseVsThetaFromFileTool::ReadNoiseVsThetaFromFileTool(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 ReadNoiseVsThetaFromFileTool::initialize() { + // Get GeoSvc + 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; + } + + // Check if cell position tool available if m_useSeg==false; if tool not + // available, try using segmentation instead + if (!m_useSeg){ + if (!m_cellPositionsTool.retrieve()) { + error() << "Unable to retrieve cell positions tool, and useSegmentation is false." << endmsg; + return StatusCode::FAILURE; + } + } + + // Get segmentation + if (m_useSeg) { + m_segmentation = dynamic_cast( + m_geoSvc->getDetector()->readout(m_readoutName).segmentation().segmentation()); + if (m_segmentation == nullptr) { + error() << "There is no module-theta segmentation." << endmsg; + return StatusCode::FAILURE; + } + else + info() << "Found module-theta segmentation." << endmsg; + } + + // open and check file, read the histograms with noise constants + if (ReadNoiseVsThetaFromFileTool::initNoiseFromFile().isFailure()) { + error() << "Couldn't open file with noise constants!!!" << endmsg; + return StatusCode::FAILURE; + } + + // 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; + + return sc; +} + +StatusCode ReadNoiseVsThetaFromFileTool::finalize() { + StatusCode sc = GaudiTool::finalize(); + return sc; +} + +StatusCode ReadNoiseVsThetaFromFileTool::initNoiseFromFile() { + // check if file exists + if (m_noiseFileName.empty()) { + error() << "Name of the file with noise values not set" << endmsg; + return StatusCode::FAILURE; + } + std::unique_ptr file(TFile::Open(m_noiseFileName.value().c_str(), "READ")); + if (file->IsZombie()) { + error() << "Couldn't open the file with noise constants" << endmsg; + return StatusCode::FAILURE; + } else { + info() << "Opening the file with noise constants: " << m_noiseFileName << endmsg; + } + + std::string elecNoiseLayerHistoName, pileupLayerHistoName; + std::string elecNoiseOffsetLayerHistoName, pileupOffsetLayerHistoName; + // Read the histograms with electronics noise and pileup from the file + for (unsigned i = 0; i < m_numRadialLayers; i++) { + elecNoiseLayerHistoName = m_elecNoiseHistoName + std::to_string(i + 1); + debug() << "Getting histogram with a name " << elecNoiseLayerHistoName << endmsg; + m_histoElecNoiseConst.push_back(*dynamic_cast(file->Get(elecNoiseLayerHistoName.c_str()))); + if (m_histoElecNoiseConst.at(i).GetNbinsX() < 1) { + error() << "Histogram " << elecNoiseLayerHistoName + << " has 0 bins! check the file with noise and the name of the histogram!" << endmsg; + return StatusCode::FAILURE; + } + if (m_setNoiseOffset){ + elecNoiseOffsetLayerHistoName = m_elecNoiseOffsetHistoName + std::to_string(i + 1); + debug() << "Getting histogram with a name " << elecNoiseOffsetLayerHistoName << endmsg; + m_histoElecNoiseOffset.push_back(*dynamic_cast(file->Get(elecNoiseOffsetLayerHistoName.c_str()))); + if (m_histoElecNoiseOffset.at(i).GetNbinsX() < 1) { + error() << "Histogram " << elecNoiseOffsetLayerHistoName + << " has 0 bins! check the file with noise and the name of the histogram!" << endmsg; + return StatusCode::FAILURE; + } + } + if (m_addPileup) { + pileupLayerHistoName = m_pileupHistoName + std::to_string(i + 1); + debug() << "Getting histogram with a name " << pileupLayerHistoName << endmsg; + m_histoPileupConst.push_back(*dynamic_cast(file->Get(pileupLayerHistoName.c_str()))); + if (m_histoPileupConst.at(i).GetNbinsX() < 1) { + error() << "Histogram " << pileupLayerHistoName + << " has 0 bins! check the file with noise and the name of the histogram!" << endmsg; + return StatusCode::FAILURE; + } + if (m_setNoiseOffset == true){ + pileupOffsetLayerHistoName = m_pileupOffsetHistoName + std::to_string(i + 1); + debug() << "Getting histogram with a name " << pileupOffsetLayerHistoName << endmsg; + m_histoPileupOffset.push_back(*dynamic_cast(file->Get(pileupOffsetLayerHistoName.c_str()))); + if (m_histoElecNoiseOffset.at(i).GetNbinsX() < 1) { + error() << "Histogram " << pileupOffsetLayerHistoName + << " has 0 bins! check the file with noise and the name of the histogram!" << endmsg; + return StatusCode::FAILURE; + } + } + } + } + // Check if we have same number of histograms (all layers) for pileup and electronics noise + if (m_histoElecNoiseConst.size() == 0 ) { + error() << "No histograms with noise found!!!!" << endmsg; + return StatusCode::FAILURE; + } + if (m_addPileup) { + if (m_histoElecNoiseConst.size() != m_histoPileupConst.size()) { + error() << "Missing histograms! Different number of histograms for electronics noise and pileup!!!!" << endmsg; + return StatusCode::FAILURE; + } + } + + return StatusCode::SUCCESS; +} + +double ReadNoiseVsThetaFromFileTool::getNoiseConstantPerCell(uint64_t aCellId) { + + double elecNoise = 0.; + double pileupNoise = 0.; + + // Get cell coordinates: eta/theta and radial layer + dd4hep::DDSegmentation::CellID cID = aCellId; + double cellTheta; + // checked that for baseline theta-module merged segmentation + // the two approaches give identical theta. + // however, code based on positioning tool is more general + // since it can be run for any segmentation class without the need + // to change the interface of this tool.. + if (m_useSeg) + cellTheta = m_segmentation->theta(aCellId); + else + cellTheta = m_cellPositionsTool->xyzPosition(aCellId).Theta(); + + unsigned cellLayer = m_decoder->get(cID, m_activeFieldName); + + // All histograms have same binning, all bins with same size + // Using the histogram in the first layer to get the bin size + unsigned index = 0; + if (m_histoElecNoiseConst.size() != 0) { + int Nbins = m_histoElecNoiseConst.at(index).GetNbinsX(); + int ibin = m_histoElecNoiseConst.at(index).FindBin(cellTheta); + if (ibin > Nbins) { + error() << "theta outside range of the histograms! Cell theta: " << cellTheta << " Nbins in histogram: " << Nbins + << endmsg; + ibin = Nbins; + } + // Check that there are not more layers than the constants are provided for + if (cellLayer < m_histoElecNoiseConst.size()) { + elecNoise = m_histoElecNoiseConst.at(cellLayer).GetBinContent(ibin); + if (m_addPileup) { + pileupNoise = m_histoPileupConst.at(cellLayer).GetBinContent(ibin); + } + } else { + error() + << "More radial layers than we have noise for!!!! Using the last layer for all histograms outside the range." + << endmsg; + } + } else { + error() << "No histograms with noise constants!!!!! " << endmsg; + } + + // Total noise: electronics noise + pileup + double totalNoise = sqrt(pow(elecNoise, 2) + pow(pileupNoise, 2)) * m_scaleFactor; + + if (totalNoise < 1e-6) { + warning() << "Zero noise: cell theta " << cellTheta << " layer " << cellLayer << " noise " << totalNoise << endmsg; + } + + return totalNoise; +} + +double ReadNoiseVsThetaFromFileTool::getNoiseOffsetPerCell(uint64_t aCellId) { + + if (!m_setNoiseOffset) + return 0.; + else { + double elecNoise = 0.; + double pileupNoise = 0.; + + // Get cell coordinates: eta and radial layer + dd4hep::DDSegmentation::CellID cID = aCellId; + double cellTheta; + if (m_useSeg) + cellTheta = m_segmentation->theta(aCellId); + else + cellTheta = m_cellPositionsTool->xyzPosition(aCellId).Theta(); + unsigned cellLayer = m_decoder->get(cID, m_activeFieldName); + + // All histograms have same binning, all bins with same size + // Using the histogram in the first layer to get the bin size + unsigned index = 0; + if (m_histoElecNoiseOffset.size() != 0) { + int Nbins = m_histoElecNoiseOffset.at(index).GetNbinsX(); + int ibin = m_histoElecNoiseOffset.at(index).FindBin(cellTheta); + if (ibin > Nbins) { + error() << "theta outside range of the histograms! Cell theta: " << cellTheta << " Nbins in histogram: " << Nbins + << endmsg; + ibin = Nbins; + } + + // Check that there are not more layers than the constants are provided for + if (cellLayer < m_histoElecNoiseOffset.size()) { + elecNoise = m_histoElecNoiseOffset.at(cellLayer).GetBinContent(ibin); + if (m_addPileup) { + pileupNoise = m_histoPileupOffset.at(cellLayer).GetBinContent(ibin); + } + } else { + error() + << "More radial layers than we have noise for!!!! Using the last layer for all histograms outside the range." + << endmsg; + } + } else { + error() << "No histograms with noise offset!!!!! " << endmsg; + } + + // Total noise: electronics noise + pileup + double totalNoise = sqrt(pow(elecNoise, 2) + pow(pileupNoise, 2)) * m_scaleFactor; + + if (totalNoise < 1e-6) { + warning() << "Zero noise: cell theta " << cellTheta << " layer " << cellLayer << " noise " << totalNoise << endmsg; + } + + return totalNoise; + } +} diff --git a/RecFCCeeCalorimeter/src/components/ReadNoiseVsThetaFromFileTool.h b/RecFCCeeCalorimeter/src/components/ReadNoiseVsThetaFromFileTool.h new file mode 100644 index 00000000..c39c4cb5 --- /dev/null +++ b/RecFCCeeCalorimeter/src/components/ReadNoiseVsThetaFromFileTool.h @@ -0,0 +1,99 @@ +#ifndef RECFCCEECALORIMETER_READNOISEVSTHETAFROMFILETOOL_H +#define RECFCCEECALORIMETER_READNOISEVSTHETAFROMFILETOOL_H + +// from Gaudi +#include "GaudiAlg/GaudiTool.h" +#include "GaudiKernel/ToolHandle.h" + +// FCCSW +//#include "k4FWCore/DataHandle.h" +//#include "k4Interface/ICalorimeterTool.h" +#include "k4Interface/INoiseConstTool.h" +#include "k4Interface/ICellPositionsTool.h" + +#include "DetSegmentation/FCCSWGridModuleThetaMerged.h" + +class IGeoSvc; + +// Root +class TH1F; + +/** @class ReadNoiseVsThetaFromFileTool + * + * Tool to read the stored noise constant per cell in the calorimeters + * Access noise constants from TH1F histogram (noise vs. theta) + * Cell positioning can be done via segmentation or via a cell + * positioning tool + * + * @author Tong Li, Giovanni Marchiori + * @date 2023-09 + * + */ + +class ReadNoiseVsThetaFromFileTool : public GaudiTool, virtual public INoiseConstTool { +public: + ReadNoiseVsThetaFromFileTool(const std::string& type, const std::string& name, const IInterface* parent); + virtual ~ReadNoiseVsThetaFromFileTool() = default; + + virtual StatusCode initialize() final; + + virtual StatusCode finalize() final; + + /// Open file and read noise histograms in the memory + StatusCode initNoiseFromFile(); + /// Find the appropriate noise constant from the histogram + double getNoiseConstantPerCell(uint64_t aCellID); + double getNoiseOffsetPerCell(uint64_t aCellID); + +private: + /// Handle for tool to get cell positions (optional - not needed if segmentation class is used) + /// Provided as a way to use the tool for other theta-based segmentations (using the positioning + /// tool instead of the segmentation to compute the theta of a cell) + /// (maybe it suffices to use a GridTheta segmentation base class?) + ToolHandle m_cellPositionsTool; + /// use segmentation (theta-module 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 (needed to get the decoder to extract e.g. layer information) + Gaudi::Property m_readoutName{this, "readoutName", "ECalBarrelThetaModuleMerged", "Name of the detector readout"}; + /// Name of active layers for sampling calorimeter + Gaudi::Property m_activeFieldName{this, "activeFieldName", "active_layer", + "Name of active layers for sampling calorimeter"}; + /// Name of pileup histogram + Gaudi::Property m_pileupHistoName{this, "pileupHistoName", "h_pileup_layer", "Name of pileup histogram"}; + /// Name of electronics noise histogram + Gaudi::Property m_elecNoiseHistoName{this, "elecNoiseHistoName", "h_elecNoise_layer", + "Name of electronics noise histogram"}; + /// Name of electronics noise offset histogram + Gaudi::Property m_elecNoiseOffsetHistoName{this, "elecNoiseOffsetHistoName", "h_mean_pileup_layer", + "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; + /// Theta-Module segmentation + dd4hep::DDSegmentation::FCCSWGridModuleThetaMerged* m_segmentation; + // Decoder for ECal layers + dd4hep::DDSegmentation::BitFieldCoder* m_decoder; + +}; + +#endif /* RECFCCEECALORIMETER_READNOISEVSTHETAFROMFILETOOL_H */