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

Fixed RP code to use local coordinates for matrix reconstruction. #545

Merged
merged 5 commits into from
Mar 16, 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
75 changes: 42 additions & 33 deletions src/detectors/RPOTS/RomanPotsReconstruction_factory.cc
Original file line number Diff line number Diff line change
Expand Up @@ -11,16 +11,25 @@
namespace eicrecon {


RomanPotsReconstruction_factory::RomanPotsReconstruction_factory(){ SetTag("ForwardRomanPotRecParticles"); }
RomanPotsReconstruction_factory::RomanPotsReconstruction_factory(){ SetTag(m_output_tag); }

void RomanPotsReconstruction_factory::Init() {

auto app = GetApplication();
std::string plugin_name = eicrecon::str::ReplaceAll(GetPluginName(), ".so", "");
std::string param_prefix = plugin_name + ":" + m_input_tag + ":";

m_log = app->GetService<Log_service>()->logger("ForwardRomanPotRecParticles");
auto app = GetApplication();

/*
auto id_spec = detector->readout(m_readout).idSpec();
m_log = app->GetService<Log_service>()->logger(m_output_tag);

m_readout = m_input_tag;

app->SetDefaultParameter(param_prefix+"readoutClass", m_readout);
m_geoSvc = app->GetService<JDD4hep_service>();

if(m_readout.empty()){ m_log->error("READOUT IS EMPTY!"); return; }

auto id_spec = m_geoSvc->detector()->readout(m_readout).idSpec();
try {
id_dec = id_spec.decoder();
if (!m_sectorField.empty()) {
Expand All @@ -36,6 +45,9 @@ namespace eicrecon {
throw JException("Failed to load ID decoder");
}


m_log->info("RP Decoding complete...");

// local detector name has higher priority
if (!m_localDetElement.empty()) {
try {
Expand All @@ -60,7 +72,7 @@ namespace eicrecon {
// fmt::join(fields, ", "))
// << endmsg;
}
*/

double det = aXRP[0][0] * aXRP[1][1] - aXRP[0][1] * aXRP[1][0];

if (det == 0) {
Expand Down Expand Up @@ -89,15 +101,10 @@ namespace eicrecon {

void RomanPotsReconstruction_factory::Process(const std::shared_ptr<const JEvent> &event) {

std::vector<edm4eic::ReconstructedParticle*> outputRPTracks;

// 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.
std::vector<edm4eic::ReconstructedParticle*> outputRPTracks;

auto converter = m_geoSvc->cellIDPositionConverter();

auto rawhits = event->Get<edm4hep::SimTrackerHit>("ForwardRomanPotHits");

Expand All @@ -118,33 +125,35 @@ namespace eicrecon {

for (const auto h: rawhits) {

// auto cellID = h->getCellID();
// The actual hit position in Global Coordinates
auto pos0 = h->getPosition();
auto cellID = h->getCellID();

//global --> local begins here -----

//TODO: this code is for the local coordinates conversion -- next PR
//auto gpos = converter->position(cellID);
// local positions
auto gpos = converter->position(cellID);

// local positions
//if (m_localDetElement.empty()) {
// auto volman = 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 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

//information is stored in cm, we need mm - divide by dd4hep::mm

if(!goodHit2 && pos0.z > 27099.0 && pos0.z < 28022.0){
if(!goodHit2 && gpos.z()/dd4hep::mm > 27099.0 && gpos.z()/dd4hep::mm < 28022.0){

goodHitX[1] = pos0.x;
goodHitY[1] = pos0.y;
goodHitZ[1] = pos0.z;
goodHitX[1] = pos0.x()/dd4hep::mm;
goodHitY[1] = pos0.y()/dd4hep::mm;
goodHitZ[1] = gpos.z()/dd4hep::mm;
goodHit2 = true;

}
if(!goodHit1 && pos0.z > 25099.0 && pos0.z < 26022.0){
if(!goodHit1 && gpos.z()/dd4hep::mm > 25099.0 && gpos.z()/dd4hep::mm < 26022.0){

goodHitX[0] = pos0.x;
goodHitY[0] = pos0.y;
goodHitZ[0] = pos0.z;
goodHitX[0] = pos0.x()/dd4hep::mm;
goodHitY[0] = pos0.y()/dd4hep::mm;
goodHitZ[0] = gpos.z()/dd4hep::mm;
goodHit1 = true;

}
Expand All @@ -160,7 +169,7 @@ namespace eicrecon {
// extract hit, subtract orbit offset – this is to get the hits in the coordinate system of the orbit
// trajectory -- should eventually be in local coordinates.
//
double XL[2] = {goodHitX[0] - local_x_offset_station_1, goodHitX[1] - local_x_offset_station_2};
double XL[2] = {goodHitX[0], goodHitX[1]};
double YL[2] = {goodHitY[0], goodHitY[1]};

double base = goodHitZ[1] - goodHitZ[0];
Expand Down
12 changes: 8 additions & 4 deletions src/detectors/RPOTS/RomanPotsReconstruction_factory.h
Original file line number Diff line number Diff line change
Expand Up @@ -9,10 +9,11 @@
#include <fmt/format.h>
#include <spdlog/spdlog.h>

#include "DD4hep/Detector.h"
#include <DDRec/CellIDPositionConverter.h>
#include <DDRec/Surface.h>
#include <DDRec/SurfaceManager.h>

#include <services/geometry/dd4hep/JDD4hep_service.h>

// Event Model related classes
#include <edm4eic/MutableReconstructedParticle.h>
Expand Down Expand Up @@ -44,8 +45,8 @@ namespace eicrecon {

//----- Define constants here ------

const double local_x_offset_station_1 = -833.3878326;
const double local_x_offset_station_2 = -924.342804;
//const double local_x_offset_station_1 = -833.3878326;
//const double local_x_offset_station_2 = -924.342804;
const double local_x_slope_offset = -0.00622147;
const double local_y_slope_offset = -0.0451035;
const double crossingAngle = -0.025;
Expand All @@ -54,17 +55,18 @@ namespace eicrecon {
std::string m_readout;
std::string m_layerField;
std::string m_sectorField;
std::string m_geoSvcName;
ajentsch marked this conversation as resolved.
Show resolved Hide resolved

dd4hep::BitFieldCoder *id_dec = nullptr;
size_t sector_idx{0}, layer_idx{0};

std::shared_ptr<JDD4hep_service> m_geoSvc;
std::string m_localDetElement;
std::vector<std::string> u_localDetFields;

dd4hep::DetElement local;
size_t local_mask = ~0;
dd4hep::Detector *detector = nullptr;
std::shared_ptr<const dd4hep::rec::CellIDPositionConverter> m_cellid_converter = nullptr;

const double aXRP[2][2] = {{2.102403743, 29.11067626},
{0.186640381, 0.192604619}};
Expand All @@ -80,6 +82,8 @@ namespace eicrecon {

private:
std::shared_ptr<spdlog::logger> m_log; /// Logger for this factory
std::string m_input_tag = "ForwardRomanPotHits";
std::string m_output_tag = "ForwardRomanPotRecParticles";

};

Expand Down