Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

EnergyInCaloLayers: Using UserDataCollections for the output #51

Merged
merged 1 commit into from
Oct 6, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
65 changes: 31 additions & 34 deletions Detector/DetStudies/src/components/EnergyInCaloLayers.cpp
Original file line number Diff line number Diff line change
@@ -1,11 +1,8 @@
#include "EnergyInCaloLayers.h"

// Key4HEP
#include "k4Interface/IGeoSvc.h"

// EDM4HEP
#include "edm4hep/CalorimeterHitCollection.h"
#include "edm4hep/MCParticleCollection.h"
// std
#include <cstddef>
#include <podio/UserDataCollection.h>

// DD4hep
#include "DD4hep/Detector.h"
Expand Down Expand Up @@ -67,12 +64,12 @@ StatusCode EnergyInCaloLayers::execute() {
auto decoder = m_geoSvc->getDetector()->readout(m_readoutName).idSpec().decoder();

// Initialize output variables
std::vector<double>* energyInLayer = m_energyInLayer.createAndPut();
energyInLayer->assign(m_numLayers, 0.);
std::vector<double>* energyInCryo = m_energyInCryo.createAndPut();
energyInCryo->assign(6, 0.);
std::vector<double>* particleVec = m_particleVec.createAndPut();
particleVec->assign(4, 0.);
auto* energyInLayerColl = m_energyInLayer.createAndPut();
energyInLayerColl->vec().assign(m_numLayers, 0.);
auto* energyInCryoColl = m_energyInCryo.createAndPut();
energyInCryoColl->vec().assign(6, 0.);
auto* particleVecColl = m_particleVec.createAndPut();
particleVecColl->vec().assign(4, 0.);

// Fill particle vector
const auto particles = m_particle.get();
Expand All @@ -83,10 +80,10 @@ StatusCode EnergyInCaloLayers::execute() {
if (particles->size() > 1) {
warning() << "Found more than one initial particle!" << endmsg;
}
particleVec->at(0) = particles->at(0).getMass();
particleVec->at(1) = particles->at(0).getMomentum().x;
particleVec->at(2) = particles->at(0).getMomentum().y;
particleVec->at(3) = particles->at(0).getMomentum().z;
particleVecColl->vec().at(0) = particles->at(0).getMass();
particleVecColl->vec().at(1) = particles->at(0).getMomentum().x;
particleVecColl->vec().at(2) = particles->at(0).getMomentum().y;
particleVecColl->vec().at(3) = particles->at(0).getMomentum().z;

// Calculate particle theta
double rxy = std::sqrt(std::pow(particles->at(0).getMomentum().x, 2) + std::pow(particles->at(0).getMomentum().y, 2));
Expand All @@ -100,48 +97,48 @@ StatusCode EnergyInCaloLayers::execute() {
size_t cryoID = decoder->get(cellID, "cryo");
if (cryoID == 0) {
size_t layerID = decoder->get(cellID, "layer");
energyInLayer->at(layerID) += hit.getEnergy();
energyInLayerColl->vec().at(layerID) += hit.getEnergy();
} else {
size_t cryoTypeID = decoder->get(cellID, "type");
energyInCryo->at(0) += hit.getEnergy();
energyInCryo->at(cryoTypeID) += hit.getEnergy();
energyInCryoColl->vec().at(0) += hit.getEnergy();
energyInCryoColl->vec().at(cryoTypeID) += hit.getEnergy();
}
}

// Calibrate energy in the calorimeter layers
for (size_t i = 0; i < energyInLayer->size(); ++i) {
energyInLayer->at(i) /= m_samplingFractions[i];
for (size_t i = 0; i < energyInLayerColl->size(); ++i) {
energyInLayerColl->vec().at(i) /= m_samplingFractions[i];
}

// Energy deposited in the whole calorimeter
double energyInCalo = 0.;
for (size_t i = 0; i < energyInLayer->size(); ++i) {
energyInCalo += energyInLayer->at(i);
for (size_t i = 0; i < energyInLayerColl->size(); ++i) {
energyInCalo += energyInLayerColl->vec().at(i);
}

// Printouts
verbose() << "********************************************************************" << endmsg;
verbose() << "Particle theta: " << particleTheta << " deg" << endmsg << endmsg;

verbose() << "Energy in layers:" << endmsg;
for (size_t i = 0; i < energyInLayer->size(); ++i) {
verbose() << " * layer " << i << ": " << energyInLayer->at(i) << " GeV" << endmsg;
for (size_t i = 0; i < energyInLayerColl->size(); ++i) {
verbose() << " * layer " << i << ": " << energyInLayerColl->vec().at(i) << " GeV" << endmsg;
}
verbose() << endmsg;

verbose() << "Energy in calorimeter: " << energyInCalo << " GeV" << endmsg;
verbose() << "Energy in cryostat: " << energyInCryo->at(0) << " GeV" << endmsg;
verbose() << "Energy in cryostat front: " << energyInCryo->at(1) << " GeV" << endmsg;
verbose() << "Energy in cryostat back: " << energyInCryo->at(2) << " GeV" << endmsg;
verbose() << "Energy in cryostat sides: " << energyInCryo->at(3) << " GeV" << endmsg;
verbose() << "Energy in cryostat LAr bath front: " << energyInCryo->at(4) << " GeV" << endmsg;
verbose() << "Energy in cryostat LAr bath back: " << energyInCryo->at(5) << " GeV" << endmsg << endmsg;
verbose() << "Energy in cryostat: " << energyInCryoColl->vec().at(0) << " GeV" << endmsg;
verbose() << "Energy in cryostat front: " << energyInCryoColl->vec().at(1) << " GeV" << endmsg;
verbose() << "Energy in cryostat back: " << energyInCryoColl->vec().at(2) << " GeV" << endmsg;
verbose() << "Energy in cryostat sides: " << energyInCryoColl->vec().at(3) << " GeV" << endmsg;
verbose() << "Energy in cryostat LAr bath front: " << energyInCryoColl->vec().at(4) << " GeV" << endmsg;
verbose() << "Energy in cryostat LAr bath back: " << energyInCryoColl->vec().at(5) << " GeV" << endmsg << endmsg;

verbose() << "Energy in calorimeter and in the cryostat front: "
<< energyInCalo + energyInCryo->at(1) + energyInCryo->at(4) << " GeV" << endmsg;
<< energyInCalo + energyInCryoColl->vec().at(1) + energyInCryoColl->vec().at(4) << " GeV" << endmsg;
verbose() << "Energy in calorimeter and in the cryostat back: "
<< energyInCalo + energyInCryo->at(2) + energyInCryo->at(5) << " GeV" << endmsg;
verbose() << "Energy in calorimeter and in the cryostat: " << energyInCalo + energyInCryo->at(0) << " GeV" << endmsg;
<< energyInCalo + energyInCryoColl->vec().at(2) + energyInCryoColl->vec().at(5) << " GeV" << endmsg;
verbose() << "Energy in calorimeter and in the cryostat: " << energyInCalo + energyInCryoColl->vec().at(0) << " GeV" << endmsg;

return StatusCode::SUCCESS;
}
Expand Down
18 changes: 9 additions & 9 deletions Detector/DetStudies/src/components/EnergyInCaloLayers.h
Original file line number Diff line number Diff line change
Expand Up @@ -5,14 +5,14 @@
#include "GaudiAlg/GaudiAlgorithm.h"

// Key4HEP
#include "k4Interface/IGeoSvc.h"
#include "k4FWCore/DataHandle.h"
class IGeoSvc;

// EDM4HEP
namespace edm4hep {
class CalorimeterHitCollection;
class MCParticleCollection;
}
// EDM4hep & Podio
#include "edm4hep/CalorimeterHitCollection.h"
#include "edm4hep/MCParticleCollection.h"
#include "podio/UserDataCollection.h"


/** @class EnergyInCaloLayers EnergyInCaloLayers.h
*
Expand Down Expand Up @@ -51,11 +51,11 @@ class EnergyInCaloLayers : public GaudiAlgorithm {
/// Handle for the particle
DataHandle<edm4hep::MCParticleCollection> m_particle{"det/particles", Gaudi::DataHandle::Reader, this};
/// Handle for vector with energy deposited in every layer
DataHandle<std::vector<double>> m_energyInLayer {"energyInLayer", Gaudi::DataHandle::Writer, this};
DataHandle<podio::UserDataCollection<double>> m_energyInLayer {"energyInLayer", Gaudi::DataHandle::Writer, this};
/// Handle for vector with energy deposited in cryostat and in its parts
DataHandle<std::vector<double>> m_energyInCryo {"energyInCryo", Gaudi::DataHandle::Writer, this};
DataHandle<podio::UserDataCollection<double>> m_energyInCryo {"energyInCryo", Gaudi::DataHandle::Writer, this};
/// Handle for initial particle vector
DataHandle<std::vector<double>> m_particleVec {"particleVec", Gaudi::DataHandle::Writer, this};
DataHandle<podio::UserDataCollection<double>> m_particleVec {"particleVec", Gaudi::DataHandle::Writer, this};

/// Pointer to the geometry service
ServiceHandle<IGeoSvc> m_geoSvc;
Expand Down
Loading