RomanPotsReconstruction_factory initial version (#536)

This PR finally provides a working Roman Pots reconstruction, for the
simplest, static matrix case (ideal for DVCS reconstruction at one beam
energy).

This will be rapidly improved in the next PRs, but we needed a working
starting first.

The current code is full of commented out portions - this will all be
cleaned up with the next PR, which will aim to improve the baseline.

src/detectors/RPOTS/RPOTS.cc

@@ -1,35 +1,22 @@
-// Copyright 2022, David Lawrence
+// Copyright 2023, Alex Jentsch
 // Subject to the terms in the LICENSE file found in the top-level directory.
 //
 //
 
 #include <JANA/JApplication.h>
+#include <JANA/JFactoryGenerator.h>
 
 #include <extensions/jana/JChainFactoryGeneratorT.h>
 
-#include <global/digi/SiliconTrackerDigi_factory.h>
-#include <global/tracking/TrackerHitReconstruction_factory.h>
+#include "RomanPotsReconstruction_factory.h"
 
-#include <algorithms/digi/SiliconTrackerDigiConfig.h>
-#include <algorithms/tracking/TrackerHitReconstructionConfig.h>
 
 extern "C" {
 void InitPlugin(JApplication *app) {
     InitJANAPlugin(app);
     using namespace eicrecon;
 
-    // Digitization
-    SiliconTrackerDigiConfig digi_default_cfg;
-    digi_default_cfg.threshold = 0;
-    digi_default_cfg.timeResolution = 8;
-    app->Add(new JChainFactoryGeneratorT<SiliconTrackerDigi_factory>({"ForwardRomanPotHits"}, "ForwardRomanPotDigiHits", digi_default_cfg));
 
-    // Convert raw digitized hits into hits with geometry info (ready for tracking)
-    TrackerHitReconstructionConfig hit_reco_cfg;
-    hit_reco_cfg.time_resolution = 8;
-    app->Add(new JChainFactoryGeneratorT<TrackerHitReconstruction_factory>(
-            {"ForwardRomanPotDigiHits"},     // Input data collection tags
-            "ForwardRomanPotRecHits",          // Output data tag
-            hit_reco_cfg));           // Hit reco default config for factories
+    app->Add(new JFactoryGeneratorT<RomanPotsReconstruction_factory>());
 }
 }

src/detectors/RPOTS/RomanPotsReconstruction_factory.cc

@@ -0,0 +1,222 @@
+// Created by Alex Jentsch to do Roman Pots reconstruction
+// Subject to the terms in the LICENSE file found in the top-level directory.
+//
+
+#include <JANA/JEvent.h>
+#include "RomanPotsReconstruction_factory.h"
+#include "services/log/Log_service.h"
+#include "extensions/spdlog/SpdlogExtensions.h"
+#include "extensions/string/StringHelpers.h"
+
+namespace eicrecon {
+
+
+    RomanPotsReconstruction_factory::RomanPotsReconstruction_factory(){ SetTag("ForwardRomanPotRecParticles"); }
+
+    void RomanPotsReconstruction_factory::Init() {
+
+        auto app = GetApplication();
+
+	m_log = app->GetService<Log_service>()->logger("ForwardRomanPotRecParticles");
+
+	/*
+        auto id_spec = detector->readout(m_readout).idSpec();
+        try {
+            id_dec = id_spec.decoder();
+            if (!m_sectorField.empty()) {
+                sector_idx = id_dec->index(m_sectorField);
+                m_log->info("Find sector field {}, index = {}", m_sectorField, sector_idx);
+            }
+            if (!m_layerField.empty()) {
+                layer_idx = id_dec->index(m_layerField);
+                m_log->info("Find layer field {}, index = {}", m_layerField, layer_idx);
+            }
+        } catch (...) {
+            m_log->error("Failed to load ID decoder for {}", m_readout);
+            throw JException("Failed to load ID decoder");
+        }
+
+        // local detector name has higher priority
+        if (!m_localDetElement.empty()) {
+            try {
+                local = detector->detector(m_localDetElement);
+                m_log->info("Local coordinate system from DetElement {}", m_localDetElement);
+            } catch (...) {
+                m_log->error("Failed to locate local coordinate system from DetElement {}", m_localDetElement);
+                throw JException("Failed to locate local coordinate system");
+            }
+            // or get from fields
+        } else {
+            std::vector<std::pair<std::string, int>> fields;
+            for (auto &f: u_localDetFields) {
+                fields.emplace_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) {
+            m_log->error("Reco matrix determinant = 0! Matrix cannot be inverted! Double-check matrix!");
+            throw JException("Reco matrix determinant = 0! Matrix cannot be inverted! Double-check matrix!");
+        }
+
+        aXRPinv[0][0] = aXRP[1][1] / det;
+        aXRPinv[0][1] = -aXRP[0][1] / det;
+        aXRPinv[1][0] = -aXRP[1][0] / det;
+        aXRPinv[1][1] = aXRP[0][0] / det;
+
+        det = aYRP[0][0] * aYRP[1][1] - aYRP[0][1] * aYRP[1][0];
+        aYRPinv[0][0] = aYRP[1][1] / det;
+        aYRPinv[0][1] = -aYRP[0][1] / det;
+        aYRPinv[1][0] = -aYRP[1][0] / det;
+        aYRPinv[1][1] = aYRP[0][0] / det;
+
+
+    }
+
+
+    void RomanPotsReconstruction_factory::ChangeRun(const std::shared_ptr<const JEvent> &event) {
+    	// Nothing to do here
+    }
+
+    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.
+
+
+	auto rawhits =  event->Get<edm4hep::SimTrackerHit>("ForwardRomanPotHits");
+
+
+        //---- begin Roman Pot Reconstruction code ----
+
+        int eventReset = 0; // counter for IDing at least one hit per layer
+        std::vector<double> hitx;
+        std::vector<double> hity;
+        std::vector<double> hitz;
+
+	double goodHitX[2] = {0.0, 0.0};
+	double goodHitY[2] = {0.0, 0.0};
+	double goodHitZ[2] = {0.0, 0.0};
+
+	bool goodHit1 = false;
+	bool goodHit2 = false;
+
+        for (const auto h: rawhits) {
+
+            // auto cellID = h->getCellID();
+            // The actual hit position in Global Coordinates
+            auto pos0 = h->getPosition();
+
+	    //TODO: this code is for the local coordinates conversion -- next PR
+            //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
+
+	    if(!goodHit2 && pos0.z > 27099.0 && pos0.z < 28022.0){
+
+		goodHitX[1] = pos0.x;
+		goodHitY[1] = pos0.y;
+		goodHitZ[1] = pos0.z;
+	    	goodHit2 = true;
+
+	    }
+	    if(!goodHit1 && pos0.z > 25099.0 && pos0.z < 26022.0){
+
+		goodHitX[0] = pos0.x;
+		goodHitY[0] = pos0.y;
+		goodHitZ[0] = pos0.z;
+		goodHit1 = true;
+
+	    }
+
+	}// end loop over hits
+
+        // NB:
+        // This is a "dumb" algorithm - I am just checking the basic thing works with a simple single-proton test.
+        // Will need to update and modify for generic # of hits for more complicated final-states.
+
+        if (goodHit1 && goodHit2) {
+
+            // 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 YL[2] = {goodHitY[0], goodHitY[1]};
+
+            double base = goodHitZ[1] - goodHitZ[0];
+
+            if (base == 0) {
+                m_log->info("Detector separation = 0! Cannot calculate slope!");
+                throw JException("Detector separation = 0! Cannot calculate slope!");
+            }
+
+            double Xip[2] = {0.0, 0.0};
+            double Xrp[2] = {XL[1], (1000 * (XL[1] - XL[0]) / (base)) - local_x_slope_offset}; //- _SX0RP_;
+            double Yip[2] = {0.0, 0.0};
+            double Yrp[2] = {YL[1], (1000 * (YL[1] - YL[0]) / (base)) - local_y_slope_offset}; //- _SY0RP_;
+
+            // use the hit information and calculated slope at the RP + the transfer matrix inverse to calculate the
+            // Polar Angle and deltaP at the IP
+
+            for (unsigned i0 = 0; i0 < 2; i0++) {
+                for (unsigned i1 = 0; i1 < 2; i1++) {
+                    Xip[i0] += aXRPinv[i0][i1] * Xrp[i1];
+                    Yip[i0] += aYRPinv[i0][i1] * Yrp[i1];
+                }
+            }
+
+            // convert polar angles to radians
+            double rsx = Xip[1] / 1000.;
+            double rsy = Yip[1] / 1000.;
+
+            // calculate momentum magnitude from measured deltaP – using thin lens optics.
+            double p = nomMomentum * (1 + 0.01 * Xip[0]);
+            double norm = std::sqrt(1.0 + rsx * rsx + rsy * rsy);
+
+            float prec[3] = {static_cast<float>((p * rsx) / norm), static_cast<float>((p * rsy) / norm),
+                                   static_cast<float>(p / norm)};
+
+	    float refPoint[3] = {static_cast<float>(goodHitX[0]), static_cast<float>(goodHitY[0]), static_cast<float>(goodHitZ[0])};
+
+            //----- end RP reconstruction code ------
+
+            edm4eic::MutableReconstructedParticle rpTrack;
+            rpTrack.setType(0);
+            rpTrack.setMomentum({prec});
+            rpTrack.setEnergy(std::hypot(edm4eic::magnitude(rpTrack.getMomentum()), .938272));
+            rpTrack.setReferencePoint({refPoint});
+            rpTrack.setCharge(1);
+            rpTrack.setMass(.938272);
+            rpTrack.setGoodnessOfPID(1.);
+            rpTrack.setPDG(2212);
+            //rpTrack.covMatrix(); // @TODO: Errors
+            outputRPTracks.push_back(new edm4eic::ReconstructedParticle(rpTrack));
+
+
+        } // END matrix reco
+
+	Set(outputRPTracks);
+
+    }
+
+}

src/detectors/RPOTS/RomanPotsReconstruction_factory.h

@@ -0,0 +1,86 @@
+// Created by Alex Jentsch
+// Subject to the terms in the LICENSE file found in the top-level directory.
+//
+
+#pragma once
+
+#include <algorithm>
+#include <cmath>
+#include <fmt/format.h>
+#include <spdlog/spdlog.h>
+
+#include <DDRec/CellIDPositionConverter.h>
+#include <DDRec/Surface.h>
+#include <DDRec/SurfaceManager.h>
+
+
+// Event Model related classes
+#include <edm4eic/MutableReconstructedParticle.h>
+#include <edm4eic/ReconstructedParticle.h>
+#include <edm4eic/TrackerHit.h>
+#include <edm4eic/vector_utils.h>
+#include <edm4hep/SimTrackerHit.h>
+
+#include <extensions/jana/JChainFactoryT.h>
+#include <extensions/spdlog/SpdlogMixin.h>
+#include <spdlog/logger.h>
+
+namespace eicrecon {
+
+    class RomanPotsReconstruction_factory : public JFactoryT<edm4eic::ReconstructedParticle>{
+
+    public:
+
+	RomanPotsReconstruction_factory(); //constructer
+
+        /** One time initialization **/
+        void Init() override;
+
+        /** On run change preparations **/
+        void ChangeRun(const std::shared_ptr<const JEvent> &event) override;
+
+        /** Event by event processing **/
+        void Process(const std::shared_ptr<const JEvent> &event) override;
+
+	//----- 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_slope_offset = -0.00622147;
+	const double local_y_slope_offset = -0.0451035;
+	const double crossingAngle = -0.025;
+	const double nomMomentum = 275.0;
+
+	std::string m_readout;
+	std::string m_layerField;
+	std::string m_sectorField;
+
+	dd4hep::BitFieldCoder *id_dec = nullptr;
+	size_t sector_idx{0}, layer_idx{0};
+
+	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}};
+	const double aYRP[2][2] = {{0.0000159900, 3.94082098},
+	                           {0.0000079946, -0.1402995}};
+
+	double aXRPinv[2][2] = {{0.0, 0.0},
+	                        {0.0, 0.0}};
+	double aYRPinv[2][2] = {{0.0, 0.0},
+	                        {0.0, 0.0}};
+
+
+
+	private:
+		std::shared_ptr<spdlog::logger> m_log;              /// Logger for this factory
+
+    };
+
+} // eicrecon

src/services/io/podio/JEventProcessorPODIO.cc

@@ -94,7 +94,7 @@ JEventProcessorPODIO::JEventProcessorPODIO() {
             "B0TrackerRecHits",
 
             //
-            "ForwardRomanPotParticles",
+            "ForwardRomanPotRecParticles",
             "SmearedFarForwardParticles",
 
             // Reconstructed data