Commit 03f9fc25 authored by Simon Spannagel's avatar Simon Spannagel
Browse files

Track: move functions into cpp file

parent 1e09d645
......@@ -225,16 +225,16 @@ Configuration Detector::getConfiguration() {
PositionVector3D<Cartesian3D<double>> Detector::getIntercept(const Track* track) {
// Get the distance from the plane to the track initial state
double distance = (m_origin.X() - track->m_state.X()) * m_normal.X();
distance += (m_origin.Y() - track->m_state.Y()) * m_normal.Y();
distance += (m_origin.Z() - track->m_state.Z()) * m_normal.Z();
distance /= (track->m_direction.X() * m_normal.X() + track->m_direction.Y() * m_normal.Y() +
track->m_direction.Z() * m_normal.Z());
double distance = (m_origin.X() - track->state().X()) * m_normal.X();
distance += (m_origin.Y() - track->state().Y()) * m_normal.Y();
distance += (m_origin.Z() - track->state().Z()) * m_normal.Z();
distance /= (track->direction().X() * m_normal.X() + track->direction().Y() * m_normal.Y() +
track->direction().Z() * m_normal.Z());
// Propagate the track
PositionVector3D<Cartesian3D<double>> globalIntercept(track->m_state.X() + distance * track->m_direction.X(),
track->m_state.Y() + distance * track->m_direction.Y(),
track->m_state.Z() + distance * track->m_direction.Z());
PositionVector3D<Cartesian3D<double>> globalIntercept(track->state().X() + distance * track->direction().X(),
track->state().Y() + distance * track->direction().Y(),
track->state().Z() + distance * track->direction().Z());
return globalIntercept;
}
......
......@@ -160,8 +160,8 @@ void AlignmentDUTResidual::MinimiseResiduals(Int_t&, Double_t*, Double_t& result
// Loop over all tracks
for(auto& track : globalTracks) {
LOG(TRACE) << "track has chi2 " << track->chi2();
LOG(TRACE) << "- track has gradient x " << track->m_direction.X();
LOG(TRACE) << "- track has gradient y " << track->m_direction.Y();
LOG(TRACE) << "- track has gradient x " << track->direction().X();
LOG(TRACE) << "- track has gradient y " << track->direction().Y();
// Find the cluster that needs to have its position recalculated
for(auto& associatedCluster : track->associatedClusters()) {
......
......@@ -273,8 +273,8 @@ bool AlignmentMillepede::putTrack(Track* track, const size_t nPlanes) {
/// Refit the track for the reference states.
track->fit();
const double tx = track->m_state.X();
const double ty = track->m_state.Y();
const double tx = track->state().X();
const double ty = track->state().Y();
// Iterate over each cluster on the track.
for(auto& cluster : track->clusters()) {
......
......@@ -205,8 +205,8 @@ StatusCode Tracking4D::run(std::shared_ptr<Clipboard> clipboard) {
trackChi2->Fill(track->chi2());
clustersPerTrack->Fill(static_cast<double>(track->nClusters()));
trackChi2ndof->Fill(track->chi2ndof());
trackAngleX->Fill(atan(track->m_direction.X()));
trackAngleY->Fill(atan(track->m_direction.Y()));
trackAngleX->Fill(atan(track->direction().X()));
trackAngleY->Fill(atan(track->direction().Y()));
// Make residuals
Clusters trackClusters = track->clusters();
......
......@@ -178,8 +178,8 @@ StatusCode TrackingSpatial::run(std::shared_ptr<Clipboard> clipboard) {
trackChi2->Fill(track->chi2());
clustersPerTrack->Fill(static_cast<double>(track->nClusters()));
trackChi2ndof->Fill(track->chi2ndof());
trackAngleX->Fill(atan(track->m_direction.X()));
trackAngleY->Fill(atan(track->m_direction.Y()));
trackAngleX->Fill(atan(track->direction().X()));
trackAngleY->Fill(atan(track->direction().Y()));
// Make residuals
Clusters trackClusters = track->clusters();
......
......@@ -22,7 +22,7 @@ ROOT_GENERATE_DICTIONARY(CorryvreckanObjectsDictionary
${CMAKE_CURRENT_SOURCE_DIR}/Pixel.h
${CMAKE_CURRENT_SOURCE_DIR}/SpidrSignal.h
${CMAKE_CURRENT_SOURCE_DIR}/Object.hpp
${CMAKE_CURRENT_SOURCE_DIR}/Track.h
${CMAKE_CURRENT_SOURCE_DIR}/Track.hpp
${CMAKE_CURRENT_SOURCE_DIR}/MCParticle.h
LINKDEF
${CMAKE_CURRENT_SOURCE_DIR}/Linkdef.h
......@@ -50,6 +50,7 @@ ADD_CUSTOM_COMMAND(OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/CorryvreckanObjectsDiction
ADD_LIBRARY(CorryvreckanObjects SHARED
Object.cpp
Pixel.cpp
Track.cpp
${CMAKE_CURRENT_BINARY_DIR}/CorryvreckanObjectsDictionary.cxx.o
)
......
#include "Track.hpp"
using namespace corryvreckan;
Track::Track() : m_direction(0, 0, 1.), m_state(0, 0, 0.) {}
Track::Track(Track* track) {
Clusters trackClusters = track->clusters();
for(auto& track_cluster : trackClusters) {
Cluster* cluster = new Cluster(track_cluster);
m_trackClusters.push_back(cluster);
}
Clusters associatedClusters = track->associatedClusters();
for(auto& assoc_cluster : associatedClusters) {
Cluster* cluster = new Cluster(assoc_cluster);
m_associatedClusters.push_back(cluster);
}
m_state = track->m_state;
m_direction = track->m_direction;
}
void Track::addCluster(Cluster* cluster) {
m_trackClusters.push_back(cluster);
}
void Track::addAssociatedCluster(Cluster* cluster) {
m_associatedClusters.push_back(cluster);
}
double Track::distance2(Cluster* cluster) const {
// Get the track X and Y at the cluster z position
double trackX = m_state.X() + m_direction.X() * cluster->globalZ();
double trackY = m_state.Y() + m_direction.Y() * cluster->globalZ();
// Calculate the 1D residuals
double dx = (trackX - cluster->globalX());
double dy = (trackY - cluster->globalY());
// Return the distance^2
return (dx * dx + dy * dy);
}
void Track::calculateChi2() {
// Get the number of clusters
m_ndof = static_cast<double>(m_trackClusters.size()) - 2.;
m_chi2 = 0.;
m_chi2ndof = 0.;
// Loop over all clusters
for(auto& cluster : m_trackClusters) {
// Get the distance^2 and the error^2
double error2 = cluster->error() * cluster->error();
m_chi2 += (this->distance2(cluster) / error2);
}
// Store also the chi2/degrees of freedom
m_chi2ndof = m_chi2 / m_ndof;
}
double Track::operator()(const double* parameters) {
// Update the track gradient and intercept
this->m_direction.SetX(parameters[0]);
this->m_state.SetX(parameters[1]);
this->m_direction.SetY(parameters[2]);
this->m_state.SetY(parameters[3]);
// Calculate the chi2
this->calculateChi2();
// Return this to minuit
return m_chi2;
}
void Track::fit() {
double vecx[2] = {0., 0.};
double vecy[2] = {0., 0.};
double matx[2][2] = {{0., 0.}, {0., 0.}};
double maty[2][2] = {{0., 0.}, {0., 0.}};
// Loop over all clusters and fill the matrices
for(auto& cluster : m_trackClusters) {
// Get the global point details
double x = cluster->globalX();
double y = cluster->globalY();
double z = cluster->globalZ();
double er2 = cluster->error() * cluster->error();
// Fill the matrices
vecx[0] += x / er2;
vecx[1] += x * z / er2;
vecy[0] += y / er2;
vecy[1] += y * z / er2;
matx[0][0] += 1. / er2;
matx[1][0] += z / er2;
matx[1][1] += z * z / er2;
maty[0][0] += 1. / er2;
maty[1][0] += z / er2;
maty[1][1] += z * z / er2;
}
// Invert the matrices
double detx = matx[0][0] * matx[1][1] - matx[1][0] * matx[1][0];
double dety = maty[0][0] * maty[1][1] - maty[1][0] * maty[1][0];
// Check for singularities.
if(detx == 0. || dety == 0.)
return;
// Get the track parameters
double slopex = (vecx[1] * matx[0][0] - vecx[0] * matx[1][0]) / detx;
double slopey = (vecy[1] * maty[0][0] - vecy[0] * maty[1][0]) / dety;
double interceptx = (vecx[0] * matx[1][1] - vecx[1] * matx[1][0]) / detx;
double intercepty = (vecy[0] * maty[1][1] - vecy[1] * maty[1][0]) / dety;
// Set the track parameters
m_state.SetX(interceptx);
m_state.SetY(intercepty);
m_state.SetZ(0.);
m_direction.SetX(slopex);
m_direction.SetY(slopey);
m_direction.SetZ(1.);
// Calculate the chi2
this->calculateChi2();
}
bool Track::isAssociated(Cluster* cluster) const {
if(std::find(m_associatedClusters.begin(), m_associatedClusters.end(), cluster) != m_associatedClusters.end()) {
return true;
}
return false;
}
ROOT::Math::XYZPoint Track::intercept(double z) const {
return m_state + m_direction * z;
}
......@@ -19,163 +19,43 @@ namespace corryvreckan {
public:
// Constructors and destructors
Track() {
m_direction.SetZ(1.);
m_state.SetZ(0.);
}
Track();
// Copy constructor (also copies clusters from the original track)
Track(Track* track) {
Clusters trackClusters = track->clusters();
for(auto& track_cluster : trackClusters) {
Cluster* cluster = new Cluster(track_cluster);
m_trackClusters.push_back(cluster);
}
Clusters associatedClusters = track->associatedClusters();
for(auto& assoc_cluster : associatedClusters) {
Cluster* cluster = new Cluster(assoc_cluster);
m_associatedClusters.push_back(cluster);
}
m_state = track->m_state;
m_direction = track->m_direction;
}
// -----------
// Functions
// -----------
Track(Track* track);
// Add a new cluster to the track
void addCluster(Cluster* cluster) { m_trackClusters.push_back(cluster); }
void addCluster(Cluster* cluster);
// Add a new cluster to the track (which will not be in the fit)
void addAssociatedCluster(Cluster* cluster) { m_associatedClusters.push_back(cluster); }
void addAssociatedCluster(Cluster* cluster);
// Calculate the 2D distance^2 between the fitted track and a cluster
double distance2(Cluster* cluster) {
double distance2(Cluster* cluster) const;
// Get the track X and Y at the cluster z position
double trackX = m_state.X() + m_direction.X() * cluster->globalZ();
double trackY = m_state.Y() + m_direction.Y() * cluster->globalZ();
// Minimisation operator used by Minuit. Minuit passes the current iteration of the parameters and checks if the chi2
// is better or worse
double operator()(const double* parameters);
// Calculate the 1D residuals
double dx = (trackX - cluster->globalX());
double dy = (trackY - cluster->globalY());
// Fit the track (linear regression)
void fit();
// Return the distance^2
return (dx * dx + dy * dy);
}
// Retrieve track parameters
double chi2() const { return m_chi2; }
double chi2ndof() const { return m_chi2ndof; }
double ndof() const { return m_ndof; }
Clusters clusters() const { return m_trackClusters; }
Clusters associatedClusters() const { return m_associatedClusters; }
bool isAssociated(Cluster* cluster) const;
// Calculate the chi2 of the track
void calculateChi2() {
// Get the number of clusters
m_ndof = static_cast<double>(m_trackClusters.size()) - 2.;
m_chi2 = 0.;
m_chi2ndof = 0.;
// Loop over all clusters
for(auto& cluster : m_trackClusters) {
// Get the distance^2 and the error^2
double error2 = cluster->error() * cluster->error();
m_chi2 += (this->distance2(cluster) / error2);
}
// Store also the chi2/degrees of freedom
m_chi2ndof = m_chi2 / m_ndof;
}
// Minimisation operator used by Minuit. Minuit passes the current iteration
// of
// the parameters and checks if the chi2 is better or worse
double operator()(const double* parameters) {
// Update the track gradient and intercept
this->m_direction.SetX(parameters[0]);
this->m_state.SetX(parameters[1]);
this->m_direction.SetY(parameters[2]);
this->m_state.SetY(parameters[3]);
// Calculate the chi2
this->calculateChi2();
// Return this to minuit
return m_chi2;
}
size_t nClusters() const { return m_trackClusters.size(); }
// Fit the track (linear regression)
void fit() {
double vecx[2] = {0., 0.};
double vecy[2] = {0., 0.};
double matx[2][2] = {{0., 0.}, {0., 0.}};
double maty[2][2] = {{0., 0.}, {0., 0.}};
// Loop over all clusters and fill the matrices
for(auto& cluster : m_trackClusters) {
// Get the global point details
double x = cluster->globalX();
double y = cluster->globalY();
double z = cluster->globalZ();
double er2 = cluster->error() * cluster->error();
// Fill the matrices
vecx[0] += x / er2;
vecx[1] += x * z / er2;
vecy[0] += y / er2;
vecy[1] += y * z / er2;
matx[0][0] += 1. / er2;
matx[1][0] += z / er2;
matx[1][1] += z * z / er2;
maty[0][0] += 1. / er2;
maty[1][0] += z / er2;
maty[1][1] += z * z / er2;
}
// Invert the matrices
double detx = matx[0][0] * matx[1][1] - matx[1][0] * matx[1][0];
double dety = maty[0][0] * maty[1][1] - maty[1][0] * maty[1][0];
// Check for singularities.
if(detx == 0. || dety == 0.)
return;
// Get the track parameters
double slopex = (vecx[1] * matx[0][0] - vecx[0] * matx[1][0]) / detx;
double slopey = (vecy[1] * maty[0][0] - vecy[0] * maty[1][0]) / dety;
double interceptx = (vecx[0] * matx[1][1] - vecx[1] * matx[1][0]) / detx;
double intercepty = (vecy[0] * maty[1][1] - vecy[1] * maty[1][0]) / dety;
// Set the track parameters
m_state.SetX(interceptx);
m_state.SetY(intercepty);
m_state.SetZ(0.);
m_direction.SetX(slopex);
m_direction.SetY(slopey);
m_direction.SetZ(1.);
// Calculate the chi2
this->calculateChi2();
}
ROOT::Math::XYZPoint intercept(double z) const;
ROOT::Math::XYZPoint state() const { return m_state; }
ROOT::Math::XYZVector direction() const { return m_direction; }
// Retrieve track parameters
double chi2() { return m_chi2; }
double chi2ndof() { return m_chi2ndof; }
double ndof() { return m_ndof; }
Clusters clusters() { return m_trackClusters; }
Clusters associatedClusters() { return m_associatedClusters; }
bool isAssociated(Cluster* cluster) {
if(std::find(m_associatedClusters.begin(), m_associatedClusters.end(), cluster) != m_associatedClusters.end()) {
return true;
}
return false;
}
size_t nClusters() { return m_trackClusters.size(); }
ROOT::Math::XYZPoint intercept(double z) {
ROOT::Math::XYZPoint point = m_state + m_direction * z;
return point;
}
private:
// Calculate the chi2 of the track
void calculateChi2();
// Member variables
Clusters m_trackClusters;
......@@ -183,11 +63,11 @@ namespace corryvreckan {
double m_chi2;
double m_ndof;
double m_chi2ndof;
ROOT::Math::XYZPoint m_state;
ROOT::Math::XYZVector m_direction;
ROOT::Math::XYZPoint m_state;
// ROOT I/O class definition - update version number when you change this class!
ClassDef(Track, 2)
ClassDef(Track, 3)
};
// Vector type declaration
......
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