Skip to content

diamond track reconstruction with local coordinates #100

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.

Sign up for GitHub

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

Draft
grzanka wants to merge 8 commits into
base: master
Choose a base branch
from
Draft

diamond track reconstruction with local coordinates #100

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
57bbde9
diamond track reconstruction with local coordinates
grzankatest Jun 18, 2025
3df80f7
Merged pps_DQMdiamondRecoLocalCoord from repository CTPPS with cms-me…
grzankatest Jul 22, 2025
78d68a3
Enhance CTPPSDiamondLocalTrackFitter to support local and global trac…
grzankatest Jul 23, 2025
b5d43e1
Refactor track analysis logic in CTPPSDiamondDQMSource for improved c…
grzankatest Jul 23, 2025
af0156d
Refactor CTPPSDiamondDQMSource to incorporate rotation angle in track…
grzankatest Jul 23, 2025
dd6a123
Adjust channel cords to defined rotation
kkawula Dec 18, 2025
097f58d
Update debug logs
kkawula Dec 18, 2025
c24c17c
Add rotation to `ProtonProducer`
kkawula Dec 19, 2025
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
78 changes: 59 additions & 19 deletions DQM/CTPPS/plugins/CTPPSDiamondDQMSource.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,19 +42,56 @@
//----------------------------------------------------------------------------------------------------

// Utility for efficiency computations
// NOTE: "CTPPSDiamondLocalTrack" refers to "local reconstruction", NOT coordinate system!
// - track.x0() is in GLOBAL coordinates (output from CTPPSDiamondLocalTrackFitter)
// - track.x0Sigma() is in LOCAL coordinates (track width)
// - det->translation() is in GLOBAL coordinates (channel center)
// - det->getDiamondDimensions() are in LOCAL coordinates
bool channelAlignedWithTrack(const CTPPSGeometry* geom,
const CTPPSDiamondDetId& detid,
const CTPPSDiamondLocalTrack& localTrack,
const CTPPSDiamondLocalTrack& track,
const float tolerance = 1) {
const DetGeomDesc* det = geom->sensor(detid);
const float x_pos = det->translation().x(),
x_width = 2.0 * det->getDiamondDimensions().xHalfWidth; // parameters stand for half the size
return ((x_pos + 0.5 * x_width > localTrack.x0() - localTrack.x0Sigma() - tolerance &&
x_pos + 0.5 * x_width < localTrack.x0() + localTrack.x0Sigma() + tolerance) ||
(x_pos - 0.5 * x_width > localTrack.x0() - localTrack.x0Sigma() - tolerance &&
x_pos - 0.5 * x_width < localTrack.x0() + localTrack.x0Sigma() + tolerance) ||
(x_pos - 0.5 * x_width < localTrack.x0() - localTrack.x0Sigma() - tolerance &&
x_pos + 0.5 * x_width > localTrack.x0() + localTrack.x0Sigma() + tolerance));
const auto& rot = det->rotation();
const auto& trans = det->translation();

// Channel dimensions in LOCAL coordinates
const float x_half_width_local = det->getDiamondDimensions().xHalfWidth;
const float channel_x_min_local = -x_half_width_local;
const float channel_x_max_local = +x_half_width_local;

// Transform track position from GLOBAL to LOCAL coordinates of this channel
auto trackGlobal = CTPPSGeometry::Vector(track.x0(), track.y0(), track.z0());
auto trackLocal = rot.Inverse() * (trackGlobal - trans);

// Now track position is in local coordinates
const float track_x_local = trackLocal.x();
const float track_x_sigma = track.x0Sigma(); // already in local coordinates
const float track_x_min_local = track_x_local - track_x_sigma - tolerance;
const float track_x_max_local = track_x_local + track_x_sigma + tolerance;

// Check overlap in LOCAL coordinates (new method)
bool overlaps = (channel_x_min_local <= track_x_max_local) && (track_x_min_local <= channel_x_max_local);

// double rot_xx, rot_xy, rot_xz, rot_yx, rot_yy, rot_yz, rot_zx, rot_zy, rot_zz;
// rot.GetComponents(rot_xx, rot_xy, rot_xz, rot_yx, rot_yy, rot_yz, rot_zx, rot_zy, rot_zz);
// double angle_z = std::atan2(rot_yx, rot_xx) * 180.0 / M_PI;

// std::cout << "=== channelAlignedWithTrack ===" << std::endl;
// std::cout << "DetId: arm=" << detid.arm() << " station=" << detid.station() << " rp=" << detid.rp() << " plane=" << detid.plane() << " channel=" << detid.channel() << std::endl;
// std::cout << "Channel position (GLOBAL): " << trans.x() << " mm" << std::endl;
// std::cout << "Channel half-width (LOCAL): " << x_half_width_local << " mm" << std::endl;
// std::cout << "Rotation around Z axis: " << angle_z << " degrees" << std::endl;
// std::cout << "Track GLOBAL: x=" << track.x0() << ", x0sigma=" << track.x0Sigma() << std::endl;
// std::cout << "Track LOCAL (after transform): x=" << track_x_local << std::endl;
// std::cout << "Channel LOCAL X range: [" << channel_x_min_local << ", " << channel_x_max_local << "]" << std::endl;
// std::cout << "Track LOCAL X range: [" << track_x_min_local << ", " << track_x_max_local << "]" << std::endl;
// std::cout << "Track X range OLD: [" << track.x0() - track.x0Sigma() - tolerance << ", " << track.x0() + track.x0Sigma() + tolerance << "]" << std::endl;
// std::cout << "Channel X range OLD: [" << trans.x() - x_half_width_local << ", " << trans.x() + x_half_width_local << "]" << std::endl;
// std::cout << "NEW method result: " << (overlaps ? "TRUE" : "FALSE") << std::endl;
// std::cout << "==========================================" << std::endl;

return overlaps;
}

namespace dds {
Expand Down Expand Up @@ -812,6 +849,11 @@ void CTPPSDiamondDQMSource::analyze(const edm::Event& event, const edm::EventSet
// to preprare correlation plot, we need to select tracks from nr_hr pot in 220m station
const CTPPSDiamondDetId detId_220nr_hr(tracks_220nr_hr.detId());

// extract the rotation angle
auto cosRotAngle_220nr_hr = diamRotations_.at(detId_220nr_hr) * CTPPSGeometry::Vector(1, 0, 0);



// selecting only tracks from 220nr station, realised as skipping tracks from 220cyl station
if ((detId_220nr_hr.rp() == CTPPS_DIAMOND_CYL_RP_ID) &&
(detId_220nr_hr.station() == CTPPS_DIAMOND_CYL_STATION_ID))
Expand All @@ -826,12 +868,13 @@ void CTPPSDiamondDQMSource::analyze(const edm::Event& event, const edm::EventSet
continue;

// get the bins from per-pot plots
int startBin_220nr_hr = trackHistoInTimeTmp->FindBin(
track_220nr_hr.x0() - diamShifts_[detId_220nr_hr.rpId()].global - track_220nr_hr.x0Sigma());
auto localTrackX_220nr_hr = (track_220nr_hr.x0() - diamTranslations_.at(detId_220nr_hr).x() + diamHalfWidths_.at(detId_220nr_hr) - track_220nr_hr.x0Sigma()) / cosRotAngle_220nr_hr.x();
int startBin_220nr_hr = trackHistoInTimeTmp->FindBin(localTrackX_220nr_hr);
int numOfBins_220nr_hr = 2 * track_220nr_hr.x0Sigma() * INV_DISPLAY_RESOLUTION_FOR_HITS_MM;

for (const auto& tracks_220cyl : *diamondLocalTracks) {
CTPPSDiamondDetId detId_220cyl(tracks_220cyl.detId());
auto cosRotAngle_220cyl = diamRotations_.at(detId_220cyl) * CTPPSGeometry::Vector(1, 0, 0);

// select tracks in the same arm, but belonging to the cylindrical pot
// that means skipping tracks from the opposite arm (skip if 220nr_hr.arm != 220cyl.arm)
Expand All @@ -849,8 +892,8 @@ void CTPPSDiamondDQMSource::analyze(const edm::Event& event, const edm::EventSet
for (const auto& track_220cyl : tracks_220cyl) {
if (!track_220cyl.isValid())
continue;
int startBin_220cyl = trackHistoTmpYAxis->FindBin(
track_220cyl.x0() - diamShifts_[detId_220cyl.rpId()].global - track_220cyl.x0Sigma());
auto localTrackX_220cyl = (track_220cyl.x0() - diamTranslations_.at(detId_220cyl).x() + diamHalfWidths_.at(detId_220cyl) - track_220cyl.x0Sigma()) / cosRotAngle_220cyl.x();
int startBin_220cyl = trackHistoTmpYAxis->FindBin(localTrackX_220cyl);
int numOfBins_220cyl = 2 * track_220cyl.x0Sigma() * INV_DISPLAY_RESOLUTION_FOR_HITS_MM;

// fill the correlation plot
Expand Down Expand Up @@ -1028,7 +1071,6 @@ void CTPPSDiamondDQMSource::analyze(const edm::Event& event, const edm::EventSet
// Using CTPPSDiamondLocalTrack
for (const auto& tracks : *diamondLocalTracks) {
const CTPPSDiamondDetId detId(tracks.detId()), detId_pot(detId.rpId());
const auto& x_shift = diamShifts_.at(detId_pot);
// extract the rotation angle for the diamond pot
auto cosRotAngle = diamRotations_.at(detId_pot) * CTPPSGeometry::Vector(1, 0, 0);

Expand All @@ -1042,9 +1084,11 @@ void CTPPSDiamondDQMSource::analyze(const edm::Event& event, const edm::EventSet
if (potPlots_.count(detId_pot) == 0)
continue;

auto localTrackX = (track.x0() - diamTranslations_.at(detId_pot).x() + diamHalfWidths_.at(detId_pot) - track.x0Sigma()) / cosRotAngle.x();

TH2F* trackHistoOOTTmp = potPlots_[detId_pot].trackDistributionOOT->getTH2F();
TAxis* trackHistoOOTTmpYAxis = trackHistoOOTTmp->GetYaxis();
int startBin = trackHistoOOTTmpYAxis->FindBin(track.x0() - x_shift.global - track.x0Sigma());
int startBin = trackHistoOOTTmpYAxis->FindBin((localTrackX));
int numOfBins = 2 * track.x0Sigma() * INV_DISPLAY_RESOLUTION_FOR_HITS_MM;
for (int i = 0; i < numOfBins; ++i)
trackHistoOOTTmp->Fill(track.ootIndex(), trackHistoOOTTmpYAxis->GetBinCenter(startBin + i));
Expand All @@ -1053,10 +1097,6 @@ void CTPPSDiamondDQMSource::analyze(const edm::Event& event, const edm::EventSet
TH1F* trackHistoInTimeTmp = potPlots_[detId_pot].trackDistribution->getTH1F();

// X coordinate of the left edge of the track in the local coordinate system
auto localTrackX =
(track.x0() - diamTranslations_.at(detId_pot).x() + diamHalfWidths_.at(detId_pot) - track.x0Sigma()) /
cosRotAngle.x();

int startBin = trackHistoInTimeTmp->FindBin((localTrackX));
int numOfBins = 2 * track.x0Sigma() * INV_DISPLAY_RESOLUTION_FOR_HITS_MM;
for (int i = 0; i < numOfBins; ++i)
Expand Down
76 changes: 69 additions & 7 deletions RecoPPS/Local/plugins/CTPPSDiamondLocalTrackFitter.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,8 @@
#include "DataFormats/CTPPSReco/interface/CTPPSDiamondLocalTrack.h"

#include "RecoPPS/Local/interface/CTPPSDiamondTrackRecognition.h"
#include "Geometry/Records/interface/VeryForwardRealGeometryRecord.h"
#include "Geometry/VeryForwardGeometryBuilder/interface/CTPPSGeometry.h"

class CTPPSDiamondLocalTrackFitter : public edm::stream::EDProducer<> {
public:
Expand All @@ -37,22 +39,27 @@ class CTPPSDiamondLocalTrackFitter : public edm::stream::EDProducer<> {
edm::EDGetTokenT<edm::DetSetVector<CTPPSDiamondRecHit> > recHitsToken_;
const edm::ParameterSet trk_algo_params_;
std::unordered_map<CTPPSDetId, std::unique_ptr<CTPPSDiamondTrackRecognition> > trk_algo_;
edm::ESGetToken<CTPPSGeometry, VeryForwardRealGeometryRecord> ctppsGeometryEventToken_;
};

CTPPSDiamondLocalTrackFitter::CTPPSDiamondLocalTrackFitter(const edm::ParameterSet& iConfig)
: recHitsToken_(
consumes<edm::DetSetVector<CTPPSDiamondRecHit> >(iConfig.getParameter<edm::InputTag>("recHitsTag"))),
trk_algo_params_(iConfig.getParameter<edm::ParameterSet>("trackingAlgorithmParams")) {
trk_algo_params_(iConfig.getParameter<edm::ParameterSet>("trackingAlgorithmParams")),
ctppsGeometryEventToken_(esConsumes<CTPPSGeometry, VeryForwardRealGeometryRecord>()){
produces<edm::DetSetVector<CTPPSDiamondLocalTrack> >();
}

void CTPPSDiamondLocalTrackFitter::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
// prepare the output
auto pOut = std::make_unique<edm::DetSetVector<CTPPSDiamondLocalTrack> >();
auto pOutLocal = std::make_unique<edm::DetSetVector<CTPPSDiamondLocalTrack> >();
auto pOutGlobal = std::make_unique<edm::DetSetVector<CTPPSDiamondLocalTrack> >();

edm::Handle<edm::DetSetVector<CTPPSDiamondRecHit> > recHits;
iEvent.getByToken(recHitsToken_, recHits);

const CTPPSGeometry* ctppsGeometry = &iSetup.getData(ctppsGeometryEventToken_);

// clear all hits possibly inherited from previous event
for (auto& algo_vs_id : trk_algo_)
algo_vs_id.second->clear();
Expand All @@ -65,17 +72,72 @@ void CTPPSDiamondLocalTrackFitter::produce(edm::Event& iEvent, const edm::EventS
trk_algo_[detid] = std::make_unique<CTPPSDiamondTrackRecognition>(trk_algo_params_);
for (const auto& hit : vec)
// skip hits without a leading edge
if (hit.ootIndex() != CTPPSDiamondRecHit::TIMESLICE_WITHOUT_LEADING)
trk_algo_[detid]->addHit(hit);
if (hit.ootIndex() != CTPPSDiamondRecHit::TIMESLICE_WITHOUT_LEADING){

// create a copy of the hit to transform it to the local coordinates in X and Y
CTPPSDiamondRecHit hitCopy(hit);
auto localVector = CTPPSGeometry::Vector(hit.x(), hit.y(), hit.z());
const auto diam = ctppsGeometry->sensor(detid);
// do the global to local transformation
// Global = Rotation * Local + Translation
// Local = Rotation^-1 * (Global - Translation)
localVector -= diam->translation();
localVector = diam->rotation().Inverse() * localVector;
hitCopy.setX(localVector.x());
hitCopy.setY(localVector.y());
// print X,Y,Z coordinates before and after transformation
// std::cout << "Hit before transformation: X=" << hit.x() << ", Y=" << hit.y() << ", Z=" << hit.z() << std::endl;
// std::cout << "Hit after transformation: X=" << localVector.x() << std::endl;

trk_algo_[detid]->addHit(hitCopy);
}

}

// build the tracks for all stations
// build the local tracks for all stations
for (auto& algo_vs_id : trk_algo_) {
auto& tracks = pOut->find_or_insert(algo_vs_id.first);
auto& tracks = pOutLocal->find_or_insert(algo_vs_id.first);
algo_vs_id.second->produceTracks(tracks);
}

iEvent.put(std::move(pOut));
for (const auto& vec : *pOutLocal) {
const CTPPSDiamondDetId raw_detid(vec.detId()), detid(raw_detid.arm(), raw_detid.station(), raw_detid.rp());
const auto diam = ctppsGeometry->sensor(detid);
pOutGlobal->find_or_insert(detid);
for (const auto& track : vec) {
auto trackCopy = track; // make a copy of the track to transform it
auto globalVector = CTPPSGeometry::Vector(track.x0(), track.y0(), track.z0());

// do the global to local transformation, doing first the rotation and then the translation
globalVector = diam->rotation() * globalVector;
globalVector += diam->translation();
auto globalPoint = math::XYZPoint(globalVector.x(), globalVector.y(), track.z0());

trackCopy.setPosition(globalPoint);
pOutGlobal->operator[](detid).push_back(trackCopy);
}
}

// // print the local tracks
// std::cout << "Local tracks for event " << iEvent.id() << ":" << std::endl;
// for (const auto& vec : *pOutLocal) {
// const CTPPSDiamondDetId raw_detid(vec.detId()), detid(raw_detid.arm(), raw_detid.station(), raw_detid.rp());
// std::cout << " DetId: " << detid << std::endl;
// for (const auto& track : vec) {
// std::cout << " Track: x0=" << track.x0() << ", y0=" << track.y0() << ", z0=" << track.z0() << std::endl;
// }
// }
// // print the global tracks
// std::cout << "Global tracks for event " << iEvent.id() << ":" << std::endl;
// for (const auto& vec : *pOutGlobal) {
// const CTPPSDiamondDetId raw_detid(vec.detId()), detid(raw_detid.arm(), raw_detid.station(), raw_detid.rp());
// std::cout << " DetId: " << detid << std::endl;
// for (const auto& track : vec) {
// std::cout << " Track: x0=" << track.x0() << ", y0=" << track.y0() << ", z0=" << track.z0() << std::endl;
// }
// }

iEvent.put(std::move(pOutGlobal));
}

void CTPPSDiamondLocalTrackFitter::fillDescriptions(edm::ConfigurationDescriptions& descr) {
Expand Down
43 changes: 35 additions & 8 deletions RecoPPS/ProtonReconstruction/plugins/CTPPSProtonProducer.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -372,23 +372,50 @@ void CTPPSProtonProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSe
if (timing_RP_track_multiplicity[tr_ti.rpId()] > max_n_timing_tracks_)
continue;

// interpolation from tracking RPs
// Interpolation in GLOBAL coordinates from tracking RPs (pixels)
const double z_ti = hGeometry->rpTranslation(tr_ti.rpId()).z();
const double f_i = (z_ti - z_j) / (z_i - z_j), f_j = (z_i - z_ti) / (z_i - z_j);
const double x_inter = f_i * tr_i.x() + f_j * tr_j.x();
const double x_inter_global = f_i * tr_i.x() + f_j * tr_j.x();
const double y_inter_global = f_i * tr_i.y() + f_j * tr_j.y();
const double x_inter_unc_sq =
f_i * f_i * tr_i.xUnc() * tr_i.xUnc() + f_j * f_j * tr_j.xUnc() * tr_j.xUnc();

const double de_x = tr_ti.x() - x_inter;
const double de_x_unc = sqrt(tr_ti.xUnc() * tr_ti.xUnc() + x_inter_unc_sq);
// Get timing detector geometry for coordinate transformation
const DetGeomDesc* det = hGeometry->sensor(tr_ti.rpId());
const auto& rot = det->rotation();
const auto& trans = det->translation();

// Transform interpolated point from GLOBAL to LOCAL coordinates of timing (diamond) detector
auto interGlobal = CTPPSGeometry::Vector(x_inter_global, y_inter_global, z_ti);
auto interLocal = rot.Inverse() * (interGlobal - trans);
const double x_inter_local = interLocal.x();

// Transform timing track from GLOBAL to LOCAL coordinates
auto trackGlobal = CTPPSGeometry::Vector(tr_ti.x(), tr_ti.y(), z_ti);
auto trackLocal = rot.Inverse() * (trackGlobal - trans);
const double x_ti_local = trackLocal.x();

// Compare in LOCAL coordinates (both now in same rotated frame)
const double de_x = x_ti_local - x_inter_local;

// The track uncertainty tr_ti.xUnc() is already in LOCAL coordinates
const double x_inter_unc_local = sqrt(x_inter_unc_sq);
const double de_x_unc = sqrt(tr_ti.xUnc() * tr_ti.xUnc() + x_inter_unc_local * x_inter_unc_local);
const double r = (de_x_unc > 0.) ? de_x / de_x_unc : 1E100;

const bool matching = (ac.getTiTrMin() < r && r < ac.getTiTrMax());

if (verbosity_)
ssLog << "ti=" << ti << ", i=" << i << ", j=" << j << " | z_ti=" << z_ti << ", z_i=" << z_i
<< ", z_j=" << z_j << " | x_ti=" << tr_ti.x() << ", x_inter=" << x_inter << ", de_x=" << de_x
<< ", de_x_unc=" << de_x_unc << ", matching=" << matching << std::endl;
if (verbosity_) {
double rot_xx, rot_xy, rot_xz, rot_yx, rot_yy, rot_yz, rot_zx, rot_zy, rot_zz;
rot.GetComponents(rot_xx, rot_xy, rot_xz, rot_yx, rot_yy, rot_yz, rot_zx, rot_zy, rot_zz);
double angle_z = std::atan2(rot_yx, rot_xx) * 180.0 / M_PI;
ssLog << "ti=" << ti << ", i=" << i << ", j=" << j
<< " | angle_z=" << angle_z << " deg"
<< " | x_ti_global=" << tr_ti.x() << ", x_ti_local=" << x_ti_local
<< " | x_inter_global=" << x_inter_global << ", x_inter_local=" << x_inter_local
<< " | de_x=" << de_x << ", de_x_unc=" << de_x_unc
<< ", r=" << r << ", matching=" << matching << std::endl;
}

if (!matching)
continue;
Expand Down