Commit dfb1f48a authored by Alex Jentsch's avatar Alex Jentsch Committed by Sylvester Joosten
Browse files

Resolve "Update RP/OMD coordinates to local"

parent 2e4368e2
......@@ -8,7 +8,12 @@
#include "GaudiAlg/Transformer.h"
#include "GaudiKernel/RndmGenerators.h"
#include "DDRec/CellIDPositionConverter.h"
#include "DDRec/Surface.h"
#include "DDRec/SurfaceManager.h"
#include "JugBase/DataHandle.h"
#include "JugBase/IGeoSvc.h"
#include "JugBase/UniqueID.h"
// Event Model related classes
......@@ -31,6 +36,21 @@ class FarForwardParticles : public GaudiAlgorithm, AlgorithmIDMixin<> {
Gaudi::Property<double> crossingAngle{this, "crossingAngle", -0.025};
Gaudi::Property<double> nomMomentum{this, "beamMomentum", 275.0};
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", ""};
Gaudi::Property<std::string> m_sectorField{this, "sectorField", ""};
SmartIF<IGeoSvc> m_geoSvc;
dd4hep::BitFieldCoder* id_dec = nullptr;
size_t sector_idx, layer_idx;
Gaudi::Property<std::string> m_localDetElement{this, "localDetElement", ""};
Gaudi::Property<std::vector<std::string>> u_localDetFields{this, "localDetFields", {}};
dd4hep::DetElement local;
size_t local_mask = ~0;
const double aXRP[2][2] = {{2.102403743, 29.11067626}, {0.186640381, 0.192604619}};
const double aYRP[2][2] = {{0.0000159900, 3.94082098}, {0.0000079946, -0.1402995}};
......@@ -44,10 +64,73 @@ public:
declareProperty("outputCollection", m_outputParticles, "ReconstructedParticles");
}
// See Wouter's example to extract local coordinates CalorimeterHitReco.cpp
//includes DDRec/CellIDPositionConverter.here
//See tutorial
// auto converter = m_GeoSvc ....
//https://eicweb.phy.anl.gov/EIC/juggler/-/blob/master/JugReco/src/components/CalorimeterHitReco.cpp - line 200
//include the Eigen libraries, used in ACTS, for the linear algebra.
StatusCode initialize() override {
if (GaudiAlgorithm::initialize().isFailure())
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;
}
// do not get the layer/sector ID if no readout class provided
if (m_readout.value().empty()) {
return StatusCode::SUCCESS;
}
auto id_spec = m_geoSvc->detector()->readout(m_readout).idSpec();
try {
id_dec = id_spec.decoder();
if (m_sectorField.value().size()) {
sector_idx = id_dec->index(m_sectorField);
info() << "Find sector field " << m_sectorField.value() << ", index = " << sector_idx << endmsg;
}
if (m_layerField.value().size()) {
layer_idx = id_dec->index(m_layerField);
info() << "Find layer field " << m_layerField.value() << ", index = " << sector_idx << endmsg;
}
} catch (...) {
error() << "Failed to load ID decoder for " << m_readout << endmsg;
return StatusCode::FAILURE;
}
// local detector name has higher priority
if (m_localDetElement.value().size()) {
try {
local = m_geoSvc->detector()->detector(m_localDetElement.value());
info() << "Local coordinate system from DetElement " << m_localDetElement.value()
<< endmsg;
} catch (...) {
error() << "Failed to locate local coordinate system from DetElement "
<< m_localDetElement.value() << endmsg;
return StatusCode::FAILURE;
}
// or get from fields
} else {
std::vector<std::pair<std::string, int>> fields;
for (auto& f : u_localDetFields.value()) {
fields.push_back({f, 0});
}
local_mask = id_spec.get_mask(fields);
// use all fields if nothing provided
if (fields.empty()) {
local_mask = ~0;
}
//info() << fmt::format("Local DetElement mask {:#064b} from fields [{}]", local_mask,
// fmt::join(fields, ", "))
// << endmsg;
}
double det = aXRP[0][0] * aXRP[1][1] - aXRP[0][1] * aXRP[1][0];
if (det == 0) {
......@@ -74,6 +157,8 @@ public:
const eic::TrackerHitCollection* rawhits = m_inputHitCollection.get();
auto& rc = *(m_outputParticles.createAndPut());
auto converter = m_geoSvc->cellIDPositionConverter();
// for (const auto& part : mc) {
// if (part.genStatus() > 1) {
// if (msgLevel(MSG::DEBUG)) {
......@@ -92,25 +177,38 @@ public:
int32_t idx = 0;
for (const auto& h : *rawhits) {
// The actual hit position:
auto pos0 = h.position();
auto cellID = h.cellID();
// The actual hit position in Global Coordinates
//auto pos0 = h.position();
auto gpos = converter->position(cellID);
// local positions
if (m_localDetElement.value().empty()) {
auto volman = m_geoSvc->detector()->volumeManager();
local = volman.lookupDetElement(cellID);
}
auto pos0 = local.nominal().worldToLocal(dd4hep::Position(gpos.x(), gpos.y(), gpos.z())); //hit position in local coordinates
// auto mom0 = h.momentum;
// auto pidCode = h.g4ID;
auto eDep = h.edep();
if (eDep < 0.00001) {
continue;
}
if (eventReset < 2) {
hitx.push_back(pos0.x - local_x_offset_station_2);
hitx.push_back(pos0.x());// - local_x_offset_station_2);
} // use station 2 for both offsets since it is used for the reference orbit
else {
hitx.push_back(pos0.x - local_x_offset_station_2);
hitx.push_back(pos0.x()); // - local_x_offset_station_2);
}
hity.push_back(pos0.y);
hitz.push_back(pos0.z);
hity.push_back(pos0.y());
hitz.push_back(pos0.z());
eventReset++;
}
......
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