Skip to content

Track based distortions update #4318

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

Merged
osbornjd merged 9 commits into from
Jun 25, 2026
Merged

Track based distortions update #4318

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
54554dd
min pt cut set to 0.2 by default.
hupereir Jun 16, 2026
a547e64
added negative and positive charge-separated matrix containers.
hupereir Jun 22, 2026
913c0b2
added accessor to all entries in container.
hupereir Jun 23, 2026
54e9407
Added charge-dependent matrix containers.
hupereir Jun 23, 2026
efe5881
added some debug output
hupereir Jun 23, 2026
c7d6d62
fixed writing charge dependent matrices to output.
hupereir Jun 23, 2026
008bfe7
restored original name
hupereir Jun 23, 2026
4d747f2
add option to ignore edge clusters, as per Christof's suggestion.
hupereir Jun 25, 2026
81cccb7
clang-tidy
hupereir Jun 25, 2026
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
112 changes: 78 additions & 34 deletions offline/packages/tpccalib/PHTpcResiduals.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -128,6 +128,8 @@ namespace
PHTpcResiduals::PHTpcResiduals(const std::string& name)
: SubsysReco(name)
, m_matrix_container(new TpcSpaceChargeMatrixContainerv2)
, m_matrix_container_pos(new TpcSpaceChargeMatrixContainerv2)
, m_matrix_container_neg(new TpcSpaceChargeMatrixContainerv2)
{
}

Expand All @@ -139,6 +141,7 @@ int PHTpcResiduals::Init(PHCompositeNode* /*topNode*/)
std::cout << "PHTpcResiduals::Init - m_maxTBeta: " << m_maxTBeta << std::endl;
std::cout << "PHTpcResiduals::Init - m_maxResidualDrphi: " << m_maxResidualDrphi << " cm" << std::endl;
std::cout << "PHTpcResiduals::Init - m_maxResidualDz: " << m_maxResidualDz << " cm" << std::endl;
std::cout << "PHTpcResiduals::Init - m_ignoreEdgeClusters: " << m_ignoreEdgeClusters << std::endl;
std::cout << "PHTpcResiduals::Init - m_minRPhiErr: " << m_minRPhiErr << " cm" << std::endl;
std::cout << "PHTpcResiduals::Init - m_minZErr: " << m_minZErr << " cm" << std::endl;
std::cout << "PHTpcResiduals::Init - m_minPt: " << m_minPt << " GeV/c" << std::endl;
Expand Down Expand Up @@ -186,11 +189,25 @@ int PHTpcResiduals::End(PHCompositeNode* /*topNode*/)
std::cout << "PHTpcResiduals::End - writing matrices to " << m_outputfile << std::endl;

// save matrix container in output file
if (m_matrix_container)
if (m_matrix_container || m_matrix_container_pos || m_matrix_container_neg)
{
std::unique_ptr<TFile> outputfile(TFile::Open(m_outputfile.c_str(), "RECREATE"));
outputfile->cd();
Comment on lines 194 to 195

@coderabbitai coderabbitai Bot Jun 25, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🩺 Stability & Availability | 🟠 Major | ⚡ Quick win

Guard against a failed TFile::Open before dereferencing.

TFile::Open returns nullptr if the output file cannot be created (bad path, permissions, disk full). The subsequent outputfile->cd() would then dereference a null pointer and crash the job at End. Note that add_from_file already guards its TFile::Open with a null check, so this path is inconsistent.

🛡️ Proposed guard 
     std::unique_ptr<TFile> outputfile(TFile::Open(m_outputfile.c_str(), "RECREATE"));
+    if (!outputfile)
+    {
+      std::cout << "PHTpcResiduals::End - could not create output file " << m_outputfile << std::endl;
+      return Fun4AllReturnCodes::EVENT_OK;
+    }
     outputfile->cd();
📝 Committable suggestion

‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.

Suggested change
std::unique_ptr<TFile> outputfile(TFile::Open(m_outputfile.c_str(), "RECREATE"));
outputfile->cd();
std::unique_ptr<TFile> outputfile(TFile::Open(m_outputfile.c_str(), "RECREATE"));
if (!outputfile)
{
std::cout << "PHTpcResiduals::End - could not create output file " << m_outputfile << std::endl;
return Fun4AllReturnCodes::EVENT_OK;
}
outputfile->cd();

m_matrix_container->Write("TpcSpaceChargeMatrixContainer");

if( m_matrix_container ) {
std::cout << "PHTpcResiduals::End - writing TpcSpaceChargeMatrixContainer object. entries: " << m_matrix_container->get_entries() << std::endl;
m_matrix_container->Write("TpcSpaceChargeMatrixContainer");
}

if( m_matrix_container_pos ) {
std::cout << "PHTpcResiduals::End - writing TpcSpaceChargeMatrixContainer_pos object. entries: " << m_matrix_container_pos->get_entries() << std::endl;
m_matrix_container_pos->Write("TpcSpaceChargeMatrixContainer_pos");
}

if( m_matrix_container_neg ) {
std::cout << "PHTpcResiduals::End - writing TpcSpaceChargeMatrixContainer_neg object. entries: " << m_matrix_container_neg->get_entries() << std::endl;
m_matrix_container_neg->Write("TpcSpaceChargeMatrixContainer_neg");
}
}

// print counters
Expand Down Expand Up @@ -456,8 +473,6 @@ void PHTpcResiduals::processTrack(SvtxTrack* track)

for (const auto& cluskey : get_cluster_keys(track))
{
// increment counter
++m_total_clusters;

// make sure cluster is from TPC
const auto detId = TrkrDefs::getTrkrId(cluskey);
Expand All @@ -466,6 +481,10 @@ void PHTpcResiduals::processTrack(SvtxTrack* track)
continue;
}

// increment counter
/* only counts TPC clusters */
++m_total_clusters;

// find matching track state
const auto stateiter = std::find_if( track->begin_states(), track->end_states(), [&cluskey]( const auto& state_pair )
{ return state_pair.second->get_cluskey() == cluskey; } );
Expand All @@ -478,6 +497,14 @@ void PHTpcResiduals::processTrack(SvtxTrack* track)

// calculate residuals with respect to cluster
auto* const cluster = m_clusterContainer->findCluster(cluskey);

// check cluster
if( !cluster ) { continue; }

// check if cluster is an edge
if( m_ignoreEdgeClusters && cluster->getEdge() > 0 && cluster->getEdge() < std::numeric_limits<char>::max() )
{ continue; }

const auto globClusPos = m_globalPositionWrapper.getGlobalPositionDistortionCorrected(cluskey, cluster, crossing);
const double clusR = get_r(globClusPos(0), globClusPos(1));
const double clusPhi = std::atan2(globClusPos(1), globClusPos(0));
Expand Down Expand Up @@ -650,43 +677,58 @@ void PHTpcResiduals::processTrack(SvtxTrack* track)
continue;
}

// Fill distortion matrices
m_matrix_container->add_to_lhs(index, 0, 0, square(clusR) / erp);
m_matrix_container->add_to_lhs(index, 0, 1, 0);
m_matrix_container->add_to_lhs(index, 0, 2, clusR * trackAlpha / erp);
std::vector<TpcSpaceChargeMatrixContainer*> containers;
containers.emplace_back( m_matrix_container.get() );
if( track->get_positive_charge() ) {
containers.emplace_back( m_matrix_container_pos.get() );
} else {
containers.emplace_back( m_matrix_container_neg.get() );
}


m_matrix_container->add_to_lhs(index, 1, 0, 0);
m_matrix_container->add_to_lhs(index, 1, 1, 1. / ez);
m_matrix_container->add_to_lhs(index, 1, 2, trackBeta / ez);
for( auto& container:containers )
{

if( !container ) { continue; }

// Fill distortion matrices
container->add_to_lhs(index, 0, 0, square(clusR) / erp);
container->add_to_lhs(index, 0, 1, 0);
container->add_to_lhs(index, 0, 2, clusR * trackAlpha / erp);

m_matrix_container->add_to_lhs(index, 2, 0, clusR * trackAlpha / erp);
m_matrix_container->add_to_lhs(index, 2, 1, trackBeta / ez);
m_matrix_container->add_to_lhs(index, 2, 2, square(trackAlpha) / erp + square(trackBeta) / ez);
container->add_to_lhs(index, 1, 0, 0);
container->add_to_lhs(index, 1, 1, 1. / ez);
container->add_to_lhs(index, 1, 2, trackBeta / ez);

m_matrix_container->add_to_rhs(index, 0, clusR * drphi / erp);
m_matrix_container->add_to_rhs(index, 1, dz / ez);
m_matrix_container->add_to_rhs(index, 2, trackAlpha * drphi / erp + trackBeta * dz / ez);
container->add_to_lhs(index, 2, 0, clusR * trackAlpha / erp);
container->add_to_lhs(index, 2, 1, trackBeta / ez);
container->add_to_lhs(index, 2, 2, square(trackAlpha) / erp + square(trackBeta) / ez);

// also update rphi reduced matrices
m_matrix_container->add_to_lhs_rphi(index, 0, 0, square(clusR) / erp);
m_matrix_container->add_to_lhs_rphi(index, 0, 1, clusR * trackAlpha / erp);
m_matrix_container->add_to_lhs_rphi(index, 1, 0, clusR * trackAlpha / erp);
m_matrix_container->add_to_lhs_rphi(index, 1, 1, square(trackAlpha) / erp);
container->add_to_rhs(index, 0, clusR * drphi / erp);
container->add_to_rhs(index, 1, dz / ez);
container->add_to_rhs(index, 2, trackAlpha * drphi / erp + trackBeta * dz / ez);

m_matrix_container->add_to_rhs_rphi(index, 0, clusR * drphi / erp);
m_matrix_container->add_to_rhs_rphi(index, 1, trackAlpha * drphi / erp);
// also update rphi reduced matrices
container->add_to_lhs_rphi(index, 0, 0, square(clusR) / erp);
container->add_to_lhs_rphi(index, 0, 1, clusR * trackAlpha / erp);
container->add_to_lhs_rphi(index, 1, 0, clusR * trackAlpha / erp);
container->add_to_lhs_rphi(index, 1, 1, square(trackAlpha) / erp);

// also update z reduced matrices
m_matrix_container->add_to_lhs_z(index, 0, 0, 1. / ez);
m_matrix_container->add_to_lhs_z(index, 0, 1, trackBeta / ez);
m_matrix_container->add_to_lhs_z(index, 1, 0, trackBeta / ez);
m_matrix_container->add_to_lhs_z(index, 1, 1, square(trackBeta) / ez);
container->add_to_rhs_rphi(index, 0, clusR * drphi / erp);
container->add_to_rhs_rphi(index, 1, trackAlpha * drphi / erp);

m_matrix_container->add_to_rhs_z(index, 0, dz / ez);
m_matrix_container->add_to_rhs_z(index, 1, trackBeta * dz / ez);
// also update z reduced matrices
container->add_to_lhs_z(index, 0, 0, 1. / ez);
container->add_to_lhs_z(index, 0, 1, trackBeta / ez);
container->add_to_lhs_z(index, 1, 0, trackBeta / ez);
container->add_to_lhs_z(index, 1, 1, square(trackBeta) / ez);

// update entries in cell
m_matrix_container->add_to_entries(index);
container->add_to_rhs_z(index, 0, dz / ez);
container->add_to_rhs_z(index, 1, trackBeta * dz / ez);

// update entries in cell
container->add_to_entries(index);
}

// increment number of accepted clusters
++m_accepted_clusters;
Expand Down Expand Up @@ -776,5 +818,7 @@ int PHTpcResiduals::getNodes(PHCompositeNode* topNode)
//____________________________________________________________________________
void PHTpcResiduals::setGridDimensions(const int phiBins, const int rBins, const int zBins)
{
m_matrix_container->set_grid_dimensions(phiBins, rBins, zBins);
if( m_matrix_container ) { m_matrix_container->set_grid_dimensions(phiBins, rBins, zBins); }
if( m_matrix_container_pos ) { m_matrix_container_pos->set_grid_dimensions(phiBins, rBins, zBins); }
if( m_matrix_container_neg ) { m_matrix_container_neg->set_grid_dimensions(phiBins, rBins, zBins); }
}
18 changes: 18 additions & 0 deletions offline/packages/tpccalib/PHTpcResiduals.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -63,11 +63,20 @@ class PHTpcResiduals : public SubsysReco
}
//@}

/// true to remove "edge" clusters
/** these are clusters that touch the edge of a detector and are considered pathological */
void setIgnoreEdgeClusters( bool value )
{ m_ignoreEdgeClusters = value; }

/// minimum value for RPhi error.
/** a too small rphi error is usually a sign of pathological TPC cluster */
void setMinRPhiErr(float minRPhiErr)
{
m_minRPhiErr = minRPhiErr;
}

/// minimum value for z error.
/** a too small z error is usually a sign of pathological TPC cluster */
void setMinZErr(float minZErr)
{
m_minZErr = minZErr;
Expand Down Expand Up @@ -171,6 +180,9 @@ class PHTpcResiduals : public SubsysReco
float m_maxTBeta = 1.5;
float m_maxResidualDz = 0.5; // cm

/// ignore edge clusters
bool m_ignoreEdgeClusters = false;

float m_minRPhiErr = 0.005; // 0.005cm -- 50um
float m_minZErr = 0.01; // 0.01cm -- 100um

Expand All @@ -193,6 +205,12 @@ class PHTpcResiduals : public SubsysReco
/// matrix container
std::unique_ptr<TpcSpaceChargeMatrixContainer> m_matrix_container;

/// matrix container positive charges only
std::unique_ptr<TpcSpaceChargeMatrixContainer> m_matrix_container_pos;

/// matrix container negative charges only
std::unique_ptr<TpcSpaceChargeMatrixContainer> m_matrix_container_neg;

// TODO: check if needed
int m_event = 0;

Expand Down
4 changes: 4 additions & 0 deletions offline/packages/tpccalib/TpcSpaceChargeMatrixContainer.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,6 +45,10 @@ class TpcSpaceChargeMatrixContainer : public PHObject
virtual int get_cell_index( int /*iphibin*/, int /*irbin*/, int /*izbin*/ ) const
{ return -1; }

/// get all entries
virtual int get_entries() const
{ return 0; }

/// get entries for a given cell
virtual int get_entries( int /*cell_index*/ ) const
{ return 0; }
Expand Down
6 changes: 6 additions & 0 deletions offline/packages/tpccalib/TpcSpaceChargeMatrixContainerv2.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,8 @@

#include "TpcSpaceChargeMatrixContainerv2.h"

#include <numeric>

//___________________________________________________________
TpcSpaceChargeMatrixContainerv2::TpcSpaceChargeMatrixContainerv2()
{
Expand Down Expand Up @@ -56,6 +58,10 @@ int TpcSpaceChargeMatrixContainerv2::get_cell_index(int iphi, int ir, int iz) co
return iz + m_zbins * (ir + m_rbins * iphi);
}

//___________________________________________________________
int TpcSpaceChargeMatrixContainerv2::get_entries() const
{ return std::accumulate( m_entries.begin(), m_entries.end(), 0); }

//___________________________________________________________
int TpcSpaceChargeMatrixContainerv2::get_entries(int cell_index) const
{
Expand Down
3 changes: 3 additions & 0 deletions offline/packages/tpccalib/TpcSpaceChargeMatrixContainerv2.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,9 @@ class TpcSpaceChargeMatrixContainerv2 : public TpcSpaceChargeMatrixContainer
/// get grid index for given sub-indexes
int get_cell_index(int iphibin, int irbin, int izbin) const override;

/// get all entries
int get_entries() const override;

/// get entries for a given cell
int get_entries(int cell_index) const override;

Expand Down
10 changes: 10 additions & 0 deletions offline/packages/tpccalib/TpcSpaceChargeMatrixInversion.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -165,6 +165,14 @@ bool TpcSpaceChargeMatrixInversion::add_from_file(const std::string& shortfilena
return false;
}

if( Verbosity() ) {
std::cout << "TpcSpaceChargeMatrixInversion::add_from_file -"
<< " file: " << filename
<< " objectname: " << objectname
<< " entries: " << source->get_entries()
<< std::endl;
}

// add object
return add(*source);
}
Expand Down Expand Up @@ -200,6 +208,8 @@ void TpcSpaceChargeMatrixInversion::calculate_distortion_corrections(const Inver
exit(1);
}

std::cout << "TpcSpaceChargeMatrixInversion::calculate_distortion_corrections - entries: " << m_matrix_container->get_entries() << std::endl;

// get grid dimensions from matrix container
int phibins = 0;
int rbins = 0;
Expand Down
6 changes: 3 additions & 3 deletions offline/packages/trackreco/PHMicromegasTpcTrackMatching.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,7 @@ class PHMicromegasTpcTrackMatching : public SubsysReco
int InitRun(PHCompositeNode* topNode) override;
int process_event(PHCompositeNode*) override;
int End(PHCompositeNode*) override;

// deprecated calls
inline void set_sc_calib_mode(const bool) {}
inline void set_collision_rate(const double) {}
Expand Down Expand Up @@ -77,7 +77,7 @@ class PHMicromegasTpcTrackMatching : public SubsysReco
unsigned int _max_tpc_layer = 55;

// pt cut for field-on data
float _pt_cut = 0.5;
float _pt_cut = 0.2;

// delta_phi window between the last cluster in the tracklet and the projection
float _dphi_cut = 0.9;
Expand All @@ -90,7 +90,7 @@ class PHMicromegasTpcTrackMatching : public SubsysReco
TrackSeedContainer* _si_track_map{nullptr};

std::string _clustermap_name = "TRKR_CLUSTER";

//! default rphi search window for each layer
std::array<double, _n_mm_layers> _rphi_search_win{0.25, 13.0};

Expand Down