Newer
Older
// It converts digitized ADC/TDC values to energy/time, and looks for geometrical information of the
// readout pixels Author: Chao Peng Date: 06/02/2021
#include <algorithm>
#include <bitset>
#include "Gaudi/Property.h"
#include "GaudiAlg/GaudiAlgorithm.h"
#include "GaudiAlg/GaudiTool.h"
#include "GaudiAlg/Transformer.h"
#include "GaudiKernel/PhysicalConstants.h"
#include "GaudiKernel/RndmGenerators.h"
#include "GaudiKernel/ToolHandle.h"
#include "DDRec/CellIDPositionConverter.h"
#include "DDRec/Surface.h"
#include "DDRec/SurfaceManager.h"
// FCCSW
#include "JugBase/DataHandle.h"
#include "JugBase/IGeoSvc.h"
// Event Model related classes
#include "eicd/ImagingPixelCollection.h"
#include "eicd/RawCalorimeterHitCollection.h"
using namespace Gaudi::Units;
namespace Jug::Reco {
class ImagingPixelReco : public GaudiAlgorithm {
public:
Gaudi::Property<std::string> m_geoSvcName{this, "geoServiceName", "GeoSvc"};
Gaudi::Property<std::string> m_readout{this, "readoutClass", ""};
Gaudi::Property<std::string> m_layerField{this, "layerField", "layer"};
Gaudi::Property<std::string> m_sectorField{this, "sectorField", "sector"};
Gaudi::Property<double> m_lUnit{this, "lengthUnit", dd4hep::mm};
Gaudi::Property<int> m_capADC{this, "capacityADC", 8096};
Gaudi::Property<int> m_pedMeanADC{this, "pedestalMean", 400};
Gaudi::Property<double> m_dyRangeADC{this, "dynamicRangeADC", 100 * MeV};
Gaudi::Property<double> m_pedSigmaADC{this, "pedestalSigma", 3.2};
Gaudi::Property<double> m_thresholdADC{this, "thresholdFactor", 3.0};
DataHandle<eic::RawCalorimeterHitCollection> m_inputHitCollection{
"inputHitCollection", Gaudi::DataHandle::Reader, this};
DataHandle<eic::ImagingPixelCollection> m_outputHitCollection{"outputHitCollection",
Gaudi::DataHandle::Writer, this};
// Pointer to the geometry service
SmartIF<IGeoSvc> m_geoSvc;
// visit readout fields
dd4hep::BitFieldCoder* id_dec;
size_t sector_idx, layer_idx;
ImagingPixelReco(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc)
{
declareProperty("inputHitCollection", m_inputHitCollection, "");
declareProperty("outputHitCollection", m_outputHitCollection, "");
if (GaudiAlgorithm::initialize().isFailure()) {
return StatusCode::FAILURE;
}
m_geoSvc = service(m_geoSvcName);
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_readout.value().empty()) {
error() << "readoutClass is not provided, it is needed to know the fields in readout ids"
<< endmsg;
return StatusCode::FAILURE;
}
try {
id_dec = m_geoSvc->detector()->readout(m_readout).idSpec().decoder();
sector_idx = id_dec->index(m_sectorField);
layer_idx = id_dec->index(m_layerField);
} catch (...) {
error() << "Failed to load ID decoder for " << m_readout << endmsg;
return StatusCode::FAILURE;
}
return StatusCode::SUCCESS;
// input collections
const auto& rawhits = *m_inputHitCollection.get();
// Create output collections
auto& hits = *m_outputHitCollection.createAndPut();
// energy time reconstruction
for (const auto& rh : rawhits) {
// did not pass the threshold
if ((rh.amplitude() - m_pedMeanADC) < m_thresholdADC * m_pedSigmaADC) {
continue;
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
double edep = (rh.amplitude() - m_pedMeanADC) / (double)m_capADC *
m_dyRangeADC; // convert ADC -> energy
double time = rh.timeStamp(); // ns
auto id = rh.cellID();
int lid = (int)id_dec->get(id, layer_idx);
int sid = (int)id_dec->get(id, sector_idx);
// global positions
auto gpos = m_geoSvc->cellIDPositionConverter()->position(id);
// local positions
auto volman = m_geoSvc->detector()->volumeManager();
auto alignment = volman.lookupDetElement(id).nominal();
auto pos = alignment.worldToLocal(dd4hep::Position(gpos.x(), gpos.y(), gpos.z()));
// polar coordinates
double r = std::sqrt(gpos.x() * gpos.x() + gpos.y() * gpos.y() + gpos.z() * gpos.z());
double th = std::acos(gpos.z() / r);
double eta = -std::log(std::tan(th / 2.));
double phi = std::atan2(gpos.y(), gpos.x());
hits.push_back(eic::ImagingPixel{
-1,
lid,
sid,
-1, // cluster id, layer id, sector id, hit id
edep,
time,
eta, // edep, time, pseudo-rapidity
{pos.x() / m_lUnit, pos.y() / m_lUnit, pos.z() / m_lUnit}, // local pos
{gpos.x() / m_lUnit, gpos.y() / m_lUnit, gpos.z() / m_lUnit}, // global pos
{r, th, phi} // polar global pos
});
}
return StatusCode::SUCCESS;