Skip to content

Commit

Permalink
Converting all remaining algorithms to EDM4HEP
Browse files Browse the repository at this point in the history
  • Loading branch information
@kjvbrt @vvolkl
kjvbrt authored and vvolkl committed Jun 28, 2021
1 parent 4c1533e commit 5cfb866
Show file tree
Hide file tree
Showing 12 changed files with 287 additions and 277 deletions.
30 changes: 15 additions & 15 deletions RecCalorimeter/src/components/ConeSelection.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
#include "ConeSelection.h"

// FCCSW
// FCC Detectors
#include "DetCommon/DetUtils.h"

// DD4hep
Expand All @@ -11,10 +11,10 @@
#include "TVector3.h"
#include "TMath.h"

// our EDM
#include "datamodel/CaloHit.h"
#include "datamodel/CaloHitCollection.h"
#include "datamodel/MCParticleCollection.h"
// EDM4HEP
#include "edm4hep/CalorimeterHit.h"
#include "edm4hep/CalorimeterHitCollection.h"
#include "edm4hep/MCParticleCollection.h"

DECLARE_COMPONENT(ConeSelection)

Expand Down Expand Up @@ -47,39 +47,39 @@ StatusCode ConeSelection::execute() {
m_cellsMap.clear();

// Get the input collection with Geant4 hits
const fcc::CaloHitCollection* cells = m_cells.get();
const edm4hep::CalorimeterHitCollection* cells = m_cells.get();
debug() << "Input Cell collection size: " << cells->size() << endmsg;
// Get the input particle collection
const fcc::MCParticleCollection* particles = m_particles.get();
const edm4hep::MCParticleCollection* particles = m_particles.get();
debug() << "Input Particle collection size: " << particles->size() << endmsg;

fcc::CaloHitCollection* edmCellsCollection = new fcc::CaloHitCollection();
edm4hep::CalorimeterHitCollection* edmCellsCollection = new edm4hep::CalorimeterHitCollection();
// Loop over all generated particles
for (const auto& part : *particles) {
TVector3 genVec(part.p4().px,part.p4().py,part.p4().pz);
TVector3 genVec(part.getMomentum().x,part.getMomentum().y,part.getMomentum().z);
auto genEta = genVec.Eta();
auto genPhi = genVec.Phi();

debug() << "Particle direction eta= " << genEta << ", phi= " << genPhi << endmsg;
// Select cells within cone around particle direction
for (const auto& cell : *cells) {
auto posCell = m_cellPositionsTool->xyzPosition(cell.cellId());
auto posCell = m_cellPositionsTool->xyzPosition(cell.getCellID());
auto eta = posCell.Eta();
auto phi = posCell.Phi();
auto circPhi = TVector2::Phi_mpi_pi(phi-genPhi);
double deltaR = double(sqrt(pow(circPhi,2)+pow((eta-genEta),2)));
if (deltaR < m_r){
//debug() << "Found a cell in cone: " << cell.cellId() << endmsg;
m_cellsMap[cell.cellId()] = cell.energy();
//debug() << "Found a cell in cone: " << cell.getCellID() << endmsg;
m_cellsMap[cell.getCellID()] = cell.getEnergy();
}
}
debug() << "Number of selected cells: " << m_cellsMap.size() << endmsg;
}

for (const auto& cell : m_cellsMap) {
fcc::CaloHit newCell = edmCellsCollection->create();
newCell.core().energy = cell.second;
newCell.core().cellId = cell.first;
edm4hep::CalorimeterHit newCell = edmCellsCollection->create();
newCell.setEnergy(cell.second);
newCell.setCellID(cell.first);
}

// push the CaloHitCollection to event store
Expand Down
21 changes: 11 additions & 10 deletions RecCalorimeter/src/components/ConeSelection.h
View file
Original file line number Diff line number Diff line change
@@ -1,16 +1,17 @@
#ifndef RECCALORIMETER_CONESELECTION_H
#define RECCALORIMETER_CONESELECTION_H

// FCCSW
#include "k4FWCore/DataHandle.h"
#include "k4Interface/ICellPositionsTool.h"

// Gaudi
#include "GaudiAlg/GaudiAlgorithm.h"

#include "datamodel/CaloHit.h"
#include "datamodel/CaloHitCollection.h"
#include "datamodel/MCParticleCollection.h"
// Key4HEP
#include "k4FWCore/DataHandle.h"
#include "k4Interface/ICellPositionsTool.h"

// EDM4HEP
#include "edm4hep/CalorimeterHit.h"
#include "edm4hep/CalorimeterHitCollection.h"
#include "edm4hep/MCParticleCollection.h"

class IGeoSvc;

Expand Down Expand Up @@ -39,11 +40,11 @@ class ConeSelection : public GaudiAlgorithm {
ToolHandle<ICellPositionsTool> m_cellPositionsTool{"CellPositionsTool", this};

/// Handle for calo hits (input collection)
DataHandle<fcc::CaloHitCollection> m_cells{"cells", Gaudi::DataHandle::Reader, this};
DataHandle<edm4hep::CalorimeterHitCollection> m_cells{"cells", Gaudi::DataHandle::Reader, this};
/// Handle for calo hits (input collection)
DataHandle<fcc::MCParticleCollection> m_particles{"particles", Gaudi::DataHandle::Reader, this};
DataHandle<edm4hep::MCParticleCollection> m_particles{"particles", Gaudi::DataHandle::Reader, this};
/// Handle for calo cells (output collection)
DataHandle<fcc::CaloHitCollection> m_selCells{"selCells", Gaudi::DataHandle::Writer, this};
DataHandle<edm4hep::CalorimeterHitCollection> m_selCells{"selCells", Gaudi::DataHandle::Writer, this};
/// Map of cell IDs (corresponding to DD4hep IDs) and energy
std::unordered_map<uint64_t, double> m_cellsMap;

Expand Down
81 changes: 39 additions & 42 deletions RecCalorimeter/src/components/CorrectECalBarrelSliWinCluster.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,12 +10,12 @@
#include "DD4hep/Detector.h"
#include "DDSegmentation/MultiSegmentation.h"

// our EDM
#include "datamodel/CaloCluster.h"
#include "datamodel/CaloClusterCollection.h"
#include "datamodel/CaloHitCollection.h"
#include "datamodel/MCParticleCollection.h"
#include "datamodel/GenVertexCollection.h"
// our EDM4HEP
#include "edm4hep/Cluster.h"
#include "edm4hep/ClusterCollection.h"
#include "edm4hep/CalorimeterHitCollection.h"
#include "edm4hep/MCParticleCollection.h"
#include "edm4hep/VertexCollection.h"

// Root
#include "TFile.h"
Expand Down Expand Up @@ -203,8 +203,8 @@ StatusCode CorrectECalBarrelSliWinCluster::initialize() {

StatusCode CorrectECalBarrelSliWinCluster::execute() {
// Get the input collection with clusters
const fcc::CaloClusterCollection* inClusters = m_inClusters.get();
fcc::CaloClusterCollection* correctedClusters = m_correctedClusters.createAndPut();
const edm4hep::ClusterCollection* inClusters = m_inClusters.get();
edm4hep::ClusterCollection* correctedClusters = m_correctedClusters.createAndPut();

// for single particle events compare with truth particles
TVector3 momentum;
Expand All @@ -216,10 +216,10 @@ StatusCode CorrectECalBarrelSliWinCluster::execute() {
const auto vertex = m_vertex.get();
if (particle->size() == 1 && vertex->size() == 1) {
for(const auto& part : *particle) {
momentum = TVector3(part.core().p4.px, part.core().p4.py, part.core().p4.pz);
momentum = TVector3(part.getMomentum().x, part.getMomentum().y, part.getMomentum().z);
etaVertex = momentum.Eta();
phiVertex = momentum.Phi();
zVertex = vertex->begin()->position().z;
zVertex = vertex->begin()->getPosition().z;
thetaVertex = 2 * atan( exp( - etaVertex ) );
verbose() << " vertex eta " << etaVertex << " phi = " << phiVertex << " theta = " << thetaVertex << " z = " << zVertex << endmsg;
}
Expand All @@ -236,14 +236,14 @@ StatusCode CorrectECalBarrelSliWinCluster::execute() {
std::vector<TLorentzVector> clustersMassInvScaled;
for (const auto& cluster : *inClusters) {
double oldEnergy = 0;
TVector3 pos(cluster.core().position.x, cluster.core().position.y, cluster.core().position.z);
TVector3 pos(cluster.getPosition().x, cluster.getPosition().y, cluster.getPosition().z);
double oldEta = pos.Eta();
double oldPhi = pos.Phi();
for (auto cell = cluster.hits_begin(); cell != cluster.hits_end(); cell++) {
oldEnergy += cell->core().energy;
oldEnergy += cell->getEnergy();
}
verbose() << " OLD ENERGY = " << oldEnergy << " from " << cluster.hits_size() << " cells" << endmsg;
verbose() << " OLD CLUSTER ENERGY = " << cluster.core().energy << endmsg;
verbose() << " OLD CLUSTER ENERGY = " << cluster.getEnergy() << endmsg;

// Do everything only using the first defined calorimeter (default: Ecal barrel)
double oldEtaId = -1;
Expand All @@ -254,38 +254,36 @@ StatusCode CorrectECalBarrelSliWinCluster::execute() {
}

// 0. Create new cluster, copy information from input
fcc::CaloCluster newCluster = correctedClusters->create();
edm4hep::Cluster newCluster = correctedClusters->create();
double energy = 0;
newCluster.core().position.x = cluster.core().position.x;
newCluster.core().position.y = cluster.core().position.y;
newCluster.core().position.z = cluster.core().position.z;
newCluster.setPosition(cluster.getPosition());
for (auto cell = cluster.hits_begin(); cell != cluster.hits_end(); cell++) {
if (m_segmentationMulti[systemId] != nullptr) {
segmentation = dynamic_cast<const dd4hep::DDSegmentation::FCCSWGridPhiEta*>(&m_segmentationMulti[systemId]->subsegmentation(cell->core().cellId));
segmentation = dynamic_cast<const dd4hep::DDSegmentation::FCCSWGridPhiEta*>(&m_segmentationMulti[systemId]->subsegmentation(cell->getCellID()));
oldEtaId = int(floor((oldEta + 0.5 * segmentation->gridSizeEta() - segmentation->offsetEta()) / segmentation->gridSizeEta()));
oldPhiId = int(floor((oldPhi + 0.5 * segmentation->gridSizePhi() - segmentation->offsetPhi()) / segmentation->gridSizePhi()));
}
if (m_decoder[systemId]->get(cell->core().cellId, "system") == systemId) {
uint layerId = m_decoder[systemId]->get(cell->core().cellId, "layer");
if (m_decoder[systemId]->get(cell->getCellID(), "system") == systemId) {
uint layerId = m_decoder[systemId]->get(cell->getCellID(), "layer");
if(m_nPhiFinal[layerId] > 0 && m_nEtaFinal[layerId] > 0) {
uint etaId = m_decoder[systemId]->get(cell->core().cellId, "eta");
uint phiId = m_decoder[systemId]->get(cell->core().cellId, "phi");
uint etaId = m_decoder[systemId]->get(cell->getCellID(), "eta");
uint phiId = m_decoder[systemId]->get(cell->getCellID(), "phi");
if ( etaId >= (oldEtaId - m_halfEtaFin[layerId]) && etaId <= (oldEtaId + m_halfEtaFin[layerId]) &&
phiId >= phiNeighbour((oldPhiId - m_halfPhiFin[layerId]), segmentation->phiBins()) && phiId <= phiNeighbour((oldPhiId + m_halfPhiFin[layerId]), segmentation->phiBins()) ) {
if (m_ellipseFinalCluster) {
if ( pow( (etaId - oldEtaId) / (m_nEtaFinal[layerId] / 2.), 2) + pow( (phiId - oldPhiId) / (m_nPhiFinal[layerId] / 2.), 2) < 1) {
newCluster.addhits(*cell);
energy += cell->core().energy;
newCluster.addToHits(*cell);
energy += cell->getEnergy();
}
} else {
newCluster.addhits(*cell);
energy += cell->core().energy;
newCluster.addToHits(*cell);
energy += cell->getEnergy();
}
}
}
}
}
newCluster.core().energy = energy;
newCluster.setEnergy(energy);

// 1. Correct eta position with log-weighting
double sumEnFirstLayer = 0;
Expand All @@ -300,9 +298,9 @@ StatusCode CorrectECalBarrelSliWinCluster::execute() {
sumWeightLayer.assign(m_numLayers, 0);
// first check the energy deposited in each layer
for (auto cell = newCluster.hits_begin(); cell != newCluster.hits_end(); cell++) {
dd4hep::DDSegmentation::CellID cID = cell->core().cellId;
dd4hep::DDSegmentation::CellID cID = cell->getCellID();
uint layer = m_decoder[systemId]->get(cID, m_layerFieldName) + m_firstLayerId;
sumEnLayer[layer] += cell->core().energy;
sumEnLayer[layer] += cell->getEnergy();
}
// sort energy to check value of 2nd highest, 3rd highest etc
for (uint iLayer = 0; iLayer < m_numLayers; iLayer++) {
Expand All @@ -313,12 +311,12 @@ StatusCode CorrectECalBarrelSliWinCluster::execute() {
// repeat but calculating eta barycentre in each layer
for (auto cell = newCluster.hits_begin(); cell != newCluster.hits_end(); cell++) {
if (m_segmentationMulti[systemId] != nullptr) {
segmentation = dynamic_cast<const dd4hep::DDSegmentation::FCCSWGridPhiEta*>(&m_segmentationMulti[systemId]->subsegmentation(cell->core().cellId));
segmentation = dynamic_cast<const dd4hep::DDSegmentation::FCCSWGridPhiEta*>(&m_segmentationMulti[systemId]->subsegmentation(cell->getCellID()));
}
dd4hep::DDSegmentation::CellID cID = cell->core().cellId;
dd4hep::DDSegmentation::CellID cID = cell->getCellID();
uint layer = m_decoder[systemId]->get(cID, m_layerFieldName) + m_firstLayerId;
double weightLog = std::max(0., m_etaRecalcLayerWeights[layer] + log(cell->core().energy / sumEnLayer[layer]));
double eta = segmentation->eta(cell->core().cellId);
double weightLog = std::max(0., m_etaRecalcLayerWeights[layer] + log(cell->getEnergy() / sumEnLayer[layer]));
double eta = segmentation->eta(cell->getCellID());
sumEtaLayer[layer] += (weightLog * eta);
sumWeightLayer[layer] += weightLog;
}
Expand Down Expand Up @@ -349,9 +347,8 @@ StatusCode CorrectECalBarrelSliWinCluster::execute() {
// alter Cartesian position of a cluster using new eta position
double radius = pos.Perp();
double phi = pos.Phi();
newCluster.core().position.x = radius * cos(phi);
newCluster.core().position.y = radius * sin(phi);
newCluster.core().position.z = radius * sinh(newEta);
auto newClusterPosition = edm4hep::Vector3f(radius * cos(phi), radius * sin(phi), radius * sinh(newEta));
newCluster.setPosition(newClusterPosition);

// 2. Correct energy for pileup noise
uint numCells = newCluster.hits_size();
Expand All @@ -363,7 +360,7 @@ StatusCode CorrectECalBarrelSliWinCluster::execute() {
} else {
noise = m_constPileupNoise * m_gauss.shoot() * std::sqrt(static_cast<int>(m_mu));
}
newCluster.core().energy += noise;
newCluster.setEnergy(newCluster.getEnergy() + noise);
m_hPileupEnergy->Fill(noise);

// 3. Correct for energy upstream
Expand All @@ -384,16 +381,16 @@ StatusCode CorrectECalBarrelSliWinCluster::execute() {
warning() << "cluster eta = " << newEta << " is larger than last defined eta values." << endmsg;
//return StatusCode::FAILURE;
}
double presamplerShift = P00 + P01 * cluster.core().energy;
double presamplerScale = P10 + P11 * sqrt(cluster.core().energy);
double presamplerShift = P00 + P01 * cluster.getEnergy();
double presamplerScale = P10 + P11 * sqrt(cluster.getEnergy());
double energyFront = presamplerShift + presamplerScale * sumEnFirstLayer * m_samplingFraction[0];
m_hUpstreamEnergy->Fill(energyFront);
newCluster.core().energy += energyFront;
newCluster.setEnergy(newCluster.getEnergy() + energyFront);

// Fill histograms
m_hEnergyPreAnyCorrections->Fill(oldEnergy);
m_hEnergyPostAllCorrections->Fill(newCluster.core().energy);
m_hEnergyPostAllCorrectionsAndScaling->Fill(newCluster.core().energy / m_response);
m_hEnergyPostAllCorrections->Fill(newCluster.getEnergy());
m_hEnergyPostAllCorrectionsAndScaling->Fill(newCluster.getEnergy() / m_response);

// Position resolution
m_hEta->Fill(newEta);
Expand Down
21 changes: 11 additions & 10 deletions RecCalorimeter/src/components/CorrectECalBarrelSliWinCluster.h
View file
Original file line number Diff line number Diff line change
@@ -1,22 +1,23 @@
#ifndef RECCALORIMETER_CORRECTECALBARRELSLIWINCLUSTER_H
#define RECCALORIMETER_CORRECTECALBARRELSLIWINCLUSTER_H

// FCCSW
// Key4HEP
#include "k4FWCore/DataHandle.h"
#include "GaudiKernel/RndmGenerators.h"
class IGeoSvc;
class IRndmGenSvc;
class ITHistSvc;

// Gaudi
#include "GaudiAlg/GaudiAlgorithm.h"
#include "GaudiKernel/ToolHandle.h"
#include "GaudiKernel/RndmGenerators.h"

namespace fcc {
class CaloClusterCollection;
class CaloHitCollection;
// EDM4HEP
namespace edm4hep {
class ClusterCollection;
class CalorimeterHitCollection;
class MCParticleCollection;
class GenVertexCollection;
class VertexCollection;
}

namespace dd4hep {
Expand Down Expand Up @@ -87,13 +88,13 @@ class CorrectECalBarrelSliWinCluster : public GaudiAlgorithm {
*/
double getNoiseConstantPerCluster(double aEta, uint numCells);
/// Handle for clusters (input collection)
DataHandle<fcc::CaloClusterCollection> m_inClusters{"clusters", Gaudi::DataHandle::Reader, this};
DataHandle<edm4hep::ClusterCollection> m_inClusters{"clusters", Gaudi::DataHandle::Reader, this};
/// Handle for corrected clusters (output collection)
DataHandle<fcc::CaloClusterCollection> m_correctedClusters{"correctedClusters", Gaudi::DataHandle::Writer, this};
DataHandle<edm4hep::ClusterCollection> m_correctedClusters{"correctedClusters", Gaudi::DataHandle::Writer, this};
/// Handle for particles with truth position and energy information: for SINGLE PARTICLE EVENTS (input collection)
DataHandle<fcc::MCParticleCollection> m_particle{"particles", Gaudi::DataHandle::Reader, this};
DataHandle<edm4hep::MCParticleCollection> m_particle{"particles", Gaudi::DataHandle::Reader, this};
/// Handle for the genvertices to read vertex position information
DataHandle<fcc::GenVertexCollection> m_vertex{"vertices", Gaudi::DataHandle::Reader, this};
DataHandle<edm4hep::VertexCollection> m_vertex{"vertices", Gaudi::DataHandle::Reader, this};
/// Pointer to the interface of histogram service
ServiceHandle<ITHistSvc> m_histSvc;
/// Pointer to the geometry service
Expand Down
Loading

0 comments on commit 5cfb866

Please sign in to comment.