Commit ab57c5af authored by Whitney Armstrong's avatar Whitney Armstrong

Fixing units in crystal digi/reco and others

- Note that units do not need to be used internally. Only when values
are needed outside of gaudi.
- The default units are (GeV, mm, ns, rad)
- Fixed units in algorithms
parent 8e69a864
Pipeline #4972 failed with stages
in 5 minutes and 51 seconds
......@@ -48,17 +48,18 @@ namespace Jug {
}
StatusCode execute() override {
// Note the energy is in units of GeV from dd4hep
// input collections
const dd4pod::CalorimeterHitCollection* simhits = m_inputHitCollection.get();
// Create output collections
auto rawhits = m_outputHitCollection.createAndPut();
eic::RawCalorimeterHitCollection* rawHitCollection = new eic::RawCalorimeterHitCollection();
for (const auto& ahit : *simhits) {
double res = m_gaussDist()/sqrt(ahit.energyDeposit()/Gaudi::Units::GeV);
double res = m_gaussDist()/sqrt(ahit.energyDeposit());
eic::RawCalorimeterHit rawhit(
(long long) ahit.cellID(),
(long long) ahit.energyDeposit() * (1. + res)/Gaudi::Units::MeV * 100.0,
(double) ahit.truth().time/Gaudi::Units::ns);
(long long) ahit.energyDeposit() * (1. + res)*1.0e6, // convert to keV integer
(double) ahit.truth().time);
rawhits->push_back(rawhit);
}
return StatusCode::SUCCESS;
......
......@@ -126,7 +126,7 @@ private:
for(auto &hit : group)
{
// not a qualified center
if(hit.energy() < m_minClusterCenterEdep) {
if(hit.energy() < m_minClusterCenterEdep/GeV) {
continue;
}
......
......@@ -9,8 +9,8 @@
#include "GaudiKernel/PhysicalConstants.h"
#include "DDRec/CellIDPositionConverter.h"
#include "DDRec/SurfaceManager.h"
#include "DDRec/Surface.h"
#include "DDRec/SurfaceManager.h"
// FCCSW
#include "JugBase/DataHandle.h"
......@@ -23,74 +23,65 @@
using namespace Gaudi::Units;
namespace Jug::Reco {
class CrystalEndcapsReco : public GaudiAlgorithm
{
public:
Gaudi::Property<double> m_minModuleEdep{this, "minModuleEdep", 0.5*MeV};
DataHandle<eic::RawCalorimeterHitCollection>
m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this};
DataHandle<eic::CalorimeterHitCollection>
m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer, this};
class CrystalEndcapsReco : public GaudiAlgorithm {
public:
Gaudi::Property<double> m_minModuleEdep{this, "minModuleEdep", 0.5 * MeV};
DataHandle<eic::RawCalorimeterHitCollection> m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader,
this};
DataHandle<eic::CalorimeterHitCollection> m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer,
this};
/// Pointer to the geometry service
SmartIF<IGeoSvc> m_geoSvc;
// ill-formed: using GaudiAlgorithm::GaudiAlgorithm;
CrystalEndcapsReco(const std::string& name, ISvcLocator* svcLoc)
: GaudiAlgorithm(name, svcLoc)
CrystalEndcapsReco(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc)
{
declareProperty("inputHitCollection", m_inputHitCollection, "");
declareProperty("outputHitCollection", m_outputHitCollection, "");
declareProperty("inputHitCollection", m_inputHitCollection, "");
declareProperty("outputHitCollection", m_outputHitCollection, "");
}
StatusCode initialize() override
{
if (GaudiAlgorithm::initialize().isFailure()) {
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;
}
return StatusCode::SUCCESS;
if (GaudiAlgorithm::initialize().isFailure()) {
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;
}
return StatusCode::SUCCESS;
}
StatusCode execute() override
{
// input collections
const auto &rawhits = *m_inputHitCollection.get();
// Create output collections
auto &hits = *m_outputHitCollection.createAndPut();
// input collections
const auto& rawhits = *m_inputHitCollection.get();
// Create output collections
auto& hits = *m_outputHitCollection.createAndPut();
// energy time reconstruction
for (auto &rh : rawhits) {
float energy = rh.amplitude()/100.*MeV;
if (energy >= m_minModuleEdep) {
float time = rh.timeStamp()*ns;
auto id = rh.cellID();
// global positions
auto gpos = m_geoSvc->cellIDPositionConverter()->position(id);
// local positions
auto pos = m_geoSvc->cellIDPositionConverter()->findContext(id)->volumePlacement().position();
// cell dimension
auto dim = m_geoSvc->cellIDPositionConverter()->cellDimensions(id);
hits.push_back(eic::CalorimeterHit{
id, energy, time,
{gpos.x(), gpos.y(), gpos.z()},
{pos.x(), pos.y(), pos.z()},
{dim[0], dim[1], 0.0},
0
});
}
// energy time reconstruction
for (auto& rh : rawhits) {
float energy = rh.amplitude() / 1.0e6; // convert keV -> GeV
if (energy >= (m_minModuleEdep / GeV)) {
float time = rh.timeStamp();
auto id = rh.cellID();
// global positions
auto gpos = m_geoSvc->cellIDPositionConverter()->position(id);
// local positions
auto pos = m_geoSvc->cellIDPositionConverter()->findContext(id)->volumePlacement().position();
// cell dimension
auto dim = m_geoSvc->cellIDPositionConverter()->cellDimensions(id);
hits.push_back(eic::CalorimeterHit{
id, energy, time, {gpos.x(), gpos.y(), gpos.z()}, {pos.x(), pos.y(), pos.z()}, {dim[0], dim[1], 0.0}, 0});
}
}
return StatusCode::SUCCESS;
return StatusCode::SUCCESS;
}
};
};
DECLARE_COMPONENT(CrystalEndcapsReco)
DECLARE_COMPONENT(CrystalEndcapsReco)
} // namespace Jug::Reco
......@@ -119,28 +119,28 @@ namespace Jug {
//m_err_eT_fit = sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime));
}
auto tsize = traj.trajectory().first.size();
debug() << "# fitted parameters : " << tsize << endmsg;
if(tsize == 0 ) continue;
traj.trajectory().second.visitBackwards(tsize-1, [&](auto&& trackstate) {
//debug() << trackstate.hasPredicted() << endmsg;
//debug() << trackstate.predicted() << endmsg;
auto params = trackstate.predicted() ;//<< endmsg;
double p0 = (1.0 / params[Acts::eBoundQOverP]) / Acts::UnitConstants::GeV;
debug() << "track predicted p = " << p0 << " GeV" << endmsg;
if ( std::abs(p0) > 500) {
debug() << "skipping" << endmsg;
return;
}
eic::Particle p({params[Acts::eBoundPhi], params[Acts::eBoundTheta], 1.0 / std::abs(params[Acts::eBoundQOverP]), 0.000511},
{0.0, 0.0, 0.0, params[Acts::eBoundTime]},
(long long)11 * params[Acts::eBoundQOverP] / std::abs(params[Acts::eBoundQOverP]), 0);
//debug() << p << endmsg;
rec_parts->push_back(p);
});
auto tsize = traj.trajectory().first.size();
debug() << "# fitted parameters : " << tsize << endmsg;
if(tsize == 0 ) continue;
traj.trajectory().second.visitBackwards(tsize-1, [&](auto&& trackstate) {
//debug() << trackstate.hasPredicted() << endmsg;
//debug() << trackstate.predicted() << endmsg;
auto params = trackstate.predicted() ;//<< endmsg;
double p0 = (1.0 / params[Acts::eBoundQOverP]) / Acts::UnitConstants::GeV;
debug() << "track predicted p = " << p0 << " GeV" << endmsg;
if ( std::abs(p0) > 500) {
debug() << "skipping" << endmsg;
return;
}
eic::Particle p({params[Acts::eBoundPhi], params[Acts::eBoundTheta], 1.0 / std::abs(params[Acts::eBoundQOverP]), 0.000511},
{0.0, 0.0, 0.0, params[Acts::eBoundTime]},
(long long)11 * params[Acts::eBoundQOverP] / std::abs(params[Acts::eBoundQOverP]), 0);
//debug() << p << endmsg;
rec_parts->push_back(p);
});
}
return StatusCode::SUCCESS;
......
......@@ -79,7 +79,7 @@ public:
for (auto &rh : rawhits) {
float npe = (rh.amplitude() - m_pedMean)/m_speMean;
if (npe >= m_minNpe) {
float time = rh.timeStamp()*m_timeStep;
float time = rh.timeStamp()*(m_timeStep/ns);
auto id = rh.cellID();
// global positions
auto gpos = m_geoSvc->cellIDPositionConverter()->position(id);
......
......@@ -69,7 +69,8 @@ namespace Jug::Reco {
std::vector<eic::CalorimeterHit> hits_B;
std::vector<eic::CalorimeterHit>& the_hits = hits_A;
std::vector<eic::CalorimeterHit>& remaining_hits = hits_B;
double max_dist = m_maxDistance.value();
double max_dist = m_maxDistance.value() / mm;
eic::CalorimeterHit ref_hit;
ref_hit.energy(0.0);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment