Merge branch 'gblNew_tracking_merged' into 'master'

Add GeneralBrokenLines as new Track Model

See merge request !230

src/core/detector/Detector.cpp

@@ -310,19 +310,25 @@ Configuration Detector::getConfiguration() const {
 // Function to get global intercept with a track
 PositionVector3D<Cartesian3D<double>> Detector::getIntercept(const Track* track) const {
 
-    // Get the distance from the plane to the track initial state
-    double distance = (m_origin.X() - track->state(m_detectorName).X()) * m_normal.X();
-    distance += (m_origin.Y() - track->state(m_detectorName).Y()) * m_normal.Y();
-    distance += (m_origin.Z() - track->state(m_detectorName).Z()) * m_normal.Z();
-    distance /= (track->direction(m_detectorName).X() * m_normal.X() + track->direction(m_detectorName).Y() * m_normal.Y() +
-                 track->direction(m_detectorName).Z() * m_normal.Z());
-
-    // Propagate the track
-    PositionVector3D<Cartesian3D<double>> globalIntercept(
-        track->state(m_detectorName).X() + distance * track->direction(m_detectorName).X(),
-        track->state(m_detectorName).Y() + distance * track->direction(m_detectorName).Y(),
-        track->state(m_detectorName).Z() + distance * track->direction(m_detectorName).Z());
-    return globalIntercept;
+    // FIXME: this is else statement can only be temporary
+    if(track->getType() == "gbl") {
+        return track->state(name());
+    } else {
+        // Get the distance from the plane to the track initial state
+        double distance = (m_origin.X() - track->state(m_detectorName).X()) * m_normal.X();
+        distance += (m_origin.Y() - track->state(m_detectorName).Y()) * m_normal.Y();
+        distance += (m_origin.Z() - track->state(m_detectorName).Z()) * m_normal.Z();
+        distance /=
+            (track->direction(m_detectorName).X() * m_normal.X() + track->direction(m_detectorName).Y() * m_normal.Y() +
+             track->direction(m_detectorName).Z() * m_normal.Z());
+
+        // Propagate the track
+        PositionVector3D<Cartesian3D<double>> globalIntercept(
+            track->state(m_detectorName).X() + distance * track->direction(m_detectorName).X(),
+            track->state(m_detectorName).Y() + distance * track->direction(m_detectorName).Y(),
+            track->state(m_detectorName).Z() + distance * track->direction(m_detectorName).Z());
+        return globalIntercept;
+    }
 }
 
 PositionVector3D<Cartesian3D<double>> Detector::getLocalIntercept(const Track* track) const {

src/modules/Tracking4D/README.md

@@ -17,7 +17,10 @@ Clusters from the first plane in Z (named the seed plane) are related to cluster
 * `exclude_dut`: Boolean to choose if the DUT plane is included in the track finding. Default value is `true`.
 * `require_detectors`: Names of detectors which are required to have a cluster on the track. If a track does not have a cluster from all detectors listed here, it is rejected. If empty, no detector is required. Default is empty.
 * `timestamp_from`: Defines the detector which provides the track timestamp. This detector also needs to be set as `required_detector`. If empty, the average timestamp of all clusters on the track will be used. Empty by default.
-* `track_model`: Select the track model used for reconstruction. Defaults to the only supported fit  `straightline`.
+* `track_model`: Select the track model used for reconstruction. A simple line fit ignoring scattering (`straightline`) and a General-Broken-Lines (`gbl`) are currently supported. Defaults to  `straightline`.
+* `momentum`: Set the beam momentum. Defaults to 5 GeV
+* `volume_scattering`: Select if volume scattering will be taken into account - defaults to false
+* `volume_radiation_length`: Define the radiation length of the volume around the telescope. Defaults to dry air with a radiation length of`304.2 m`
 
 ### Plots produced
 

src/modules/Tracking4D/Tracking4D.cpp

@@ -33,6 +33,16 @@ Tracking4D::Tracking4D(Configuration config, std::vector<std::shared_ptr<Detecto
     requireDetectors = m_config.getArray<std::string>("require_detectors", {""});
     timestampFrom = m_config.get<std::string>("timestamp_from", {});
     trackModel = m_config.get<std::string>("track_model", "straightline");
+    momentum = m_config.get<double>("momentum", Units::get<double>(5, "GeV"));
+    volumeRadiationLength = m_config.get<double>("volume_radiation_length", Units::get<double>(304.2, "m"));
+    useVolumeScatterer = m_config.get<bool>("volume_scattering", false);
+
+    // print a warning if volumeScatterer are used as this causes fit failures
+    // that are still not understood
+    if(useVolumeScatterer) {
+        LOG_ONCE(WARNING) << "Taking volume scattering effects into account is still WIP and causes the GBL to fail - these "
+                             "tracks are rejected";
+    }
     // spatial cut, relative (x * spatial_resolution) or absolute:
     if(m_config.count({"spatial_cut_rel", "spatial_cut_abs"}) > 1) {
         throw InvalidCombinationError(
@@ -71,11 +81,6 @@ void Tracking4D::initialise() {
     for(auto& detector : get_detectors()) {
         auto detectorID = detector->name();
 
-        // Do not create plots for detectors not participating in the tracking:
-        if(excludeDUT && detector->isDUT()) {
-            continue;
-        }
-
         // Do not created plots for auxiliary detectors:
         if(detector->isAuxiliary()) {
             continue;
@@ -83,29 +88,44 @@ void Tracking4D::initialise() {
 
         TDirectory* directory = getROOTDirectory();
         TDirectory* local_directory = directory->mkdir(detectorID.c_str());
+
         if(local_directory == nullptr) {
             throw RuntimeError("Cannot create or access local ROOT directory for module " + this->getUniqueName());
         }
         local_directory->cd();
 
+        title = detectorID + " kink X;kink [rad];events";
+        kinkX[detectorID] = new TH1F("kinkX", title.c_str(), 500, -0.01, -0.01);
+        title = detectorID + " kinkY ;kink [rad];events";
+        kinkY[detectorID] = new TH1F("kinkY", title.c_str(), 500, -0.01, -0.01);
+
+        // Do not create plots for detectors not participating in the tracking:
+        if(excludeDUT && detector->isDUT()) {
+            continue;
+        }
+
         title = detectorID + " Residual X;x_{track}-x [mm];events";
         residualsX[detectorID] = new TH1F("residualsX", title.c_str(), 500, -0.1, 0.1);
-        title = detectorID + " Residual X, size 1;x_{track}-x [mm];events";
+        title = detectorID + " Residual X, cluster column width 1;x_{track}-x [mm];events";
         residualsXwidth1[detectorID] = new TH1F("residualsXwidth1", title.c_str(), 500, -0.1, 0.1);
-        title = detectorID + " Residual X, size 2;x_{track}-x [mm];events";
+        title = detectorID + " Residual X, cluster column width  2;x_{track}-x [mm];events";
         residualsXwidth2[detectorID] = new TH1F("residualsXwidth2", title.c_str(), 500, -0.1, 0.1);
-        title = detectorID + " Residual X, size 3;x_{track}-x [mm];events";
+        title = detectorID + " Residual X, cluster column width  3;x_{track}-x [mm];events";
         residualsXwidth3[detectorID] = new TH1F("residualsXwidth3", title.c_str(), 500, -0.1, 0.1);
         title = detectorID + " Residual Y;y_{track}-y [mm];events";
         residualsY[detectorID] = new TH1F("residualsY", title.c_str(), 500, -0.1, 0.1);
-        title = detectorID + " Residual Y, size 1;y_{track}-y [mm];events";
+        title = detectorID + " Residual Y, cluster row width 1;y_{track}-y [mm];events";
         residualsYwidth1[detectorID] = new TH1F("residualsYwidth1", title.c_str(), 500, -0.1, 0.1);
-        title = detectorID + " Residual Y, size 2;y_{track}-y [mm];events";
+        title = detectorID + " Residual Y, cluster row width 2;y_{track}-y [mm];events";
         residualsYwidth2[detectorID] = new TH1F("residualsYwidth2", title.c_str(), 500, -0.1, 0.1);
-        title = detectorID + " Residual Y, size 3;y_{track}-y [mm];events";
+        title = detectorID + " Residual Y, cluster row width 3;y_{track}-y [mm];events";
         residualsYwidth3[detectorID] = new TH1F("residualsYwidth3", title.c_str(), 500, -0.1, 0.1);
 
-        directory->cd();
+        title = detectorID + " Pull X;x_{track}-x/resolution;events";
+        pullX[detectorID] = new TH1F("pullX", title.c_str(), 500, -5, 5);
+
+        title = detectorID + " Pull Y;y_{track}-y/resolution;events";
+        pullY[detectorID] = new TH1F("pully", title.c_str(), 500, -5, 5);
     }
 }
 
@@ -141,6 +161,9 @@ StatusCode Tracking4D::run(std::shared_ptr<Clipboard> clipboard) {
             KDTree* clusterTree = new KDTree();
             clusterTree->buildTimeTree(*tempClusters);
             trees[detectorID] = clusterTree;
+        }
+        // the detector always needs to be listed as we would like to add the material budget information
+        if(!detector->isAuxiliary()) {
             detectors.push_back(detectorID);
         }
     }
@@ -166,15 +189,25 @@ StatusCode Tracking4D::run(std::shared_ptr<Clipboard> clipboard) {
         LOG(DEBUG) << "Looking at next seed cluster";
 
         auto track = Track::Factory(trackModel);
+        // The track finding is based on a straight line. Therefore a refTrack to extrapolate to the next plane is used here
+        auto refTrack = new StraightLineTrack();
         // Add the cluster to the track
         track->addCluster(cluster);
         track->setTimestamp(cluster->timestamp());
+        track->setVolumeScatter(volumeRadiationLength);
+        track->useVolumeScatterer(useVolumeScatterer);
+        refTrack->addCluster(cluster);
+        refTrack->setTimestamp(cluster->timestamp());
 
+        track->setParticleMomentum(momentum);
         // Loop over each subsequent plane and look for a cluster within the timing cuts
         for(auto& detectorID : detectors) {
             // Get the detector
             auto det = get_detector(detectorID);
 
+            // always add the material budget:
+            track->addMaterial(detectorID, det->materialBudget(), det->displacement().z());
+            LOG(TRACE) << "added material budget for " << detectorID << " at z = " << det->displacement().z();
             if(trees.count(detectorID) == 0) {
                 LOG(TRACE) << "Skipping detector " << det->name() << " as it has 0 clusters.";
                 continue;
@@ -214,10 +247,10 @@ StatusCode Tracking4D::run(std::shared_ptr<Clipboard> clipboard) {
 
             // Now look for the spatially closest cluster on the next plane
             double interceptX, interceptY;
-            if(track->nClusters() > 1) {
-                track->fit();
+            if(refTrack->nClusters() > 1) {
+                refTrack->fit(); // fixme: this is not really a nice way to get the details
 
-                PositionVector3D<Cartesian3D<double>> interceptPoint = det->getIntercept(track);
+                PositionVector3D<Cartesian3D<double>> interceptPoint = det->getIntercept(refTrack);
                 interceptX = interceptPoint.X();
                 interceptY = interceptPoint.Y();
             } else {
@@ -265,7 +298,8 @@ StatusCode Tracking4D::run(std::shared_ptr<Clipboard> clipboard) {
             // Add the cluster to the track
             LOG(DEBUG) << "- added cluster to track";
             track->addCluster(closestCluster);
-        } //*/
+            refTrack->addCluster(closestCluster);
+        }
 
         // check if track has required detector(s):
         auto foundRequiredDetector = [this](Track* t) {
@@ -289,37 +323,55 @@ StatusCode Tracking4D::run(std::shared_ptr<Clipboard> clipboard) {
             delete track;
             continue;
         }
-
+        delete refTrack;
         // Fit the track and save it
         track->fit();
-        tracks->push_back(track);
-
+        if(track->isFitted()) {
+            tracks->push_back(track);
+        } else {
+            LOG_N(WARNING, 100) << "Rejected a track due to failure in fitting";
+            delete track;
+            continue;
+        }
         // Fill histograms
         trackChi2->Fill(track->chi2());
         clustersPerTrack->Fill(static_cast<double>(track->nClusters()));
         trackChi2ndof->Fill(track->chi2ndof());
-        trackAngleX->Fill(atan(track->direction(track->clusters().front()->detectorID()).X()));
-        trackAngleY->Fill(atan(track->direction(track->clusters().front()->detectorID()).Y()));
-
+        if(!(trackModel == "gbl")) {
+            trackAngleX->Fill(atan(track->direction(track->clusters().front()->detectorID()).X()));
+            trackAngleY->Fill(atan(track->direction(track->clusters().front()->detectorID()).Y()));
+        }
         // Make residuals
         auto trackClusters = track->clusters();
         for(auto& trackCluster : trackClusters) {
             string detectorID = trackCluster->detectorID();
-            ROOT::Math::XYZPoint intercept = track->intercept(trackCluster->global().z());
-            residualsX[detectorID]->Fill(intercept.X() - trackCluster->global().x());
+            residualsX[detectorID]->Fill(track->residual(detectorID).X());
+            pullX[detectorID]->Fill(track->residual(detectorID).X() / track->clusters().front()->errorX());
+            pullY[detectorID]->Fill(track->residual(detectorID).Y() / track->clusters().front()->errorY());
             if(trackCluster->columnWidth() == 1)
-                residualsXwidth1[detectorID]->Fill(intercept.X() - trackCluster->global().x());
+                residualsXwidth1[detectorID]->Fill(track->residual(detectorID).X());
             if(trackCluster->columnWidth() == 2)
-                residualsXwidth2[detectorID]->Fill(intercept.X() - trackCluster->global().x());
+                residualsXwidth2[detectorID]->Fill(track->residual(detectorID).X());
             if(trackCluster->columnWidth() == 3)
-                residualsXwidth3[detectorID]->Fill(intercept.X() - trackCluster->global().x());
-            residualsY[detectorID]->Fill(intercept.Y() - trackCluster->global().y());
+                residualsXwidth3[detectorID]->Fill(track->residual(detectorID).X());
+            residualsY[detectorID]->Fill(track->residual(detectorID).Y());
             if(trackCluster->rowWidth() == 1)
-                residualsYwidth1[detectorID]->Fill(intercept.Y() - trackCluster->global().y());
+                residualsYwidth1[detectorID]->Fill(track->residual(detectorID).Y());
             if(trackCluster->rowWidth() == 2)
-                residualsYwidth2[detectorID]->Fill(intercept.Y() - trackCluster->global().y());
+                residualsYwidth2[detectorID]->Fill(track->residual(detectorID).Y());
             if(trackCluster->rowWidth() == 3)
-                residualsYwidth3[detectorID]->Fill(intercept.Y() - trackCluster->global().y());
+                residualsYwidth3[detectorID]->Fill(track->residual(detectorID).Y());
+        }
+
+        for(auto& det : detectors) {
+            if(!kinkX.count(det)) {
+                LOG(WARNING) << "Skipping writing kinks due to missing init of histograms for  " << det;
+                continue;
+            }
+
+            XYPoint kink = track->kink(det);
+            kinkX.at(det)->Fill(kink.x());
+            kinkY.at(det)->Fill(kink.y());
         }
 
         if(timestampFrom.empty()) {

src/modules/Tracking4D/Tracking4D.h

@@ -41,10 +41,17 @@ namespace corryvreckan {
         std::map<std::string, TH1F*> residualsYwidth2;
         std::map<std::string, TH1F*> residualsYwidth3;
 
+        std::map<std::string, TH1F*> kinkX;
+        std::map<std::string, TH1F*> kinkY;
+        std::map<std::string, TH1F*> pullX;
+        std::map<std::string, TH1F*> pullY;
         // Cuts for tracking
+        double momentum{};
+        double volumeRadiationLength{};
         double time_cut_reference_;
         size_t minHitsOnTrack;
         bool excludeDUT;
+        bool useVolumeScatterer{};
         std::vector<std::string> requireDetectors;
         std::map<std::shared_ptr<Detector>, double> time_cuts_;
         std::map<std::shared_ptr<Detector>, XYVector> spatial_cuts_;

src/objects/CMakeLists.txt

@@ -25,6 +25,7 @@ ROOT_GENERATE_DICTIONARY(CorryvreckanObjectsDictionary
     ${CMAKE_CURRENT_SOURCE_DIR}/SpidrSignal.hpp
     ${CMAKE_CURRENT_SOURCE_DIR}/Track.hpp
     ${CMAKE_CURRENT_SOURCE_DIR}/StraightLineTrack.hpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/GblTrack.hpp
     ${CMAKE_CURRENT_SOURCE_DIR}/MCParticle.hpp
     ${CMAKE_CURRENT_SOURCE_DIR}/Event.hpp
     LINKDEF
@@ -56,6 +57,7 @@ ADD_LIBRARY(CorryvreckanObjects SHARED
     Cluster.cpp
     Track.cpp
     StraightLineTrack.cpp
+    GblTrack.cpp
     KDTree.cpp
     Event.cpp
     ${CMAKE_CURRENT_BINARY_DIR}/CorryvreckanObjectsDictionary.cxx.o

src/objects/GblTrack.cpp

+#include "GblTrack.hpp"
+#include "Track.hpp"
+#include "exceptions.h"
+
+#include "GblPoint.h"
+#include "GblTrajectory.h"
+
+using namespace corryvreckan;
+using namespace gbl;
+
+GblTrack::GblTrack() : Track() {}
+
+GblTrack::GblTrack(const GblTrack& track) : Track(track) {
+    if(track.getType() != this->getType())
+        throw TrackModelChanged(typeid(*this), track.getType(), this->getType());
+    m_planes = track.m_planes;
+}
+
+void GblTrack::fit() {
+
+    // Fitting with 2 clusters or less is pointless
+    if(m_trackClusters.size() < 2) {
+        throw TrackError(typeid(GblTrack), " attempting to fit a track with less than 2 clusters");
+    }
+    // store the used clusters in a map for easy access:
+    std::map<std::string, Cluster*> clusters;
+    for(auto& c : m_trackClusters) {
+        auto cluster = dynamic_cast<Cluster*>(c.GetObject());
+        clusters[cluster->detectorID()] = cluster;
+    }
+    // create a list of planes and sort it, also calculate the material budget:
+    m_planes.clear();
+    double total_material = 0;
+    for(auto l : m_materialBudget) {
+        total_material += l.second.first;
+        Plane current(l.second.second, l.second.first, l.first, clusters.count(l.first));
+        if(current.hasCluster()) {
+            current.setPosition(clusters.at(current.name())->global().z());
+            current.setCluster(clusters.at(current.name()));
+        }
+        m_planes.push_back(current);
+    }
+    std::sort(m_planes.begin(), m_planes.end());
+    // add volume scattering length - we ignore for now the material thickness while considering air
+    total_material += (m_planes.end()->postion() - m_planes.front().postion()) / m_scattering_length_volume;
+
+    std::vector<GblPoint> points;
+    // get the seedcluster for the fit - simply the first one in the list
+    auto seedcluster = m_planes.front().cluster();
+    double prevPos = 0;
+
+    // lambda to calculate the scattering theta
+    auto scatteringTheta = [this](double mbCurrent, double mbTotal) -> double {
+        return (13.6 / m_momentum * sqrt(mbCurrent) * (1 + 0.038 * log(mbTotal)));
+    };
+    // lambda to add measurement (if existing) and scattering from single plane
+    auto addScattertoGblPoint = [&total_material, &scatteringTheta](GblPoint& point, double material) {
+        Eigen::Vector2d scatterWidth;
+        scatterWidth(0) = 1 / (scatteringTheta(material, total_material) * scatteringTheta(material, total_material));
+        scatterWidth(1) = scatterWidth(0);
+        point.addScatterer(Eigen::Vector2d::Zero(), scatterWidth);
+
+    };
+    auto JacToNext = [](double val) {
+        Matrix5d Jac;
+        Jac = Matrix5d::Identity();
+        Jac(3, 1) = val;
+        Jac(4, 2) = Jac(3, 1);
+        return Jac;
+    };
+    auto addMeasurementtoGblPoint = [&seedcluster](GblPoint& point, std::vector<Plane>::iterator& p) {
+        auto cluster = p->cluster();
+        Eigen::Vector2d initialResidual;
+        initialResidual(0) = cluster->global().x() - seedcluster->global().x();
+        initialResidual(1) = cluster->global().y() - seedcluster->global().y();
+        // FIXME uncertainty of single hit - rotations ignored
+        Eigen::Matrix2d covv = Eigen::Matrix2d::Identity();
+        covv(0, 0) = 1. / cluster->errorX() / cluster->errorX();
+        covv(1, 1) = 1. / cluster->errorY() / cluster->errorY();
+        point.addMeasurement(initialResidual, covv);
+    };
+    // lambda to add plane (not the first one) and air scatterers //FIXME: Where to put them?
+    auto addPlane = [&JacToNext, &prevPos, &addMeasurementtoGblPoint, &addScattertoGblPoint, &points, this](
+        std::vector<Plane>::iterator& plane) {
+        double dist = plane->postion() - prevPos;
+        double frac1 = 0.21, frac2 = 0.58;
+        // Current layout
+        // |        |        |       |
+        // |  frac1 | frac2  | frac1 |
+        // |        |        |       |
+
+        // first air scatterer:
+        auto point = GblPoint(JacToNext(frac1 * dist));
+        addScattertoGblPoint(point, dist / 2. / m_scattering_length_volume);
+        if(m_use_volume_scatter) {
+            points.push_back(point);
+            point = GblPoint(JacToNext(frac2 * dist));
+            addScattertoGblPoint(point, dist / 2. / m_scattering_length_volume);
+            points.push_back(point);
+        } else {
+            frac1 = 1;
+        }
+        point = GblPoint(JacToNext(frac1 * dist));
+        addScattertoGblPoint(point, plane->materialbudget());
+        if(plane->hasCluster()) {
+            addMeasurementtoGblPoint(point, plane);
+        }
+        prevPos = plane->postion();
+        points.push_back(point);
+        plane->setGblPos(unsigned(points.size())); // gbl starts counting at 1
+    };
+
+    // First GblPoint
+    std::vector<Plane>::iterator pl = m_planes.begin();
+    auto point = GblPoint(Matrix5d::Identity());
+    addScattertoGblPoint(point, pl->materialbudget());
+    addMeasurementtoGblPoint(point, pl);
+    points.push_back(point);
+    pl->setGblPos(1);
+    pl++;
+    // add all other points
+    for(; pl != m_planes.end(); ++pl) {
+        addPlane(pl);
+    }
+
+    // Make sure we missed nothing
+    if((points.size() != ((m_materialBudget.size() * 3) - 2) && m_use_volume_scatter) ||
+       (points.size() != m_materialBudget.size() && !m_use_volume_scatter)) {
+        throw TrackError(typeid(GblTrack),
+                         "Number of planes " + std::to_string(m_materialBudget.size()) +
+                             " doesn't match number of GBL points on trajectory " + std::to_string(points.size()));
+    }
+    // perform fit
+    GblTrajectory traj(points, false); // false = no magnetic field
+    double lostWeight = 0;
+    int ndf = 0;
+    // fit it
+    auto fitReturnValue = traj.fit(m_chi2, ndf, lostWeight);
+    if(fitReturnValue != 0) { // Is this a good ieda? should we discard track candidates that fail?
+        fitReturnValue = traj.fit(m_chi2, ndf, lostWeight);
+        if(fitReturnValue != 0) { // Is this a good ieda? should we discard track candidates that fail?
+            return;
+            //            throw TrackFitError(typeid(GblTrack), "internal GBL Error " + std::to_string(fitReturnValue));
+        }
+    }
+
+    // copy the results
+    m_ndof = double(ndf);
+    m_chi2ndof = (m_ndof < 0.0) ? -1 : (m_chi2 / m_ndof);
+    Eigen::VectorXd localPar(5);
+    Eigen::MatrixXd localCov(5, 5);
+    Eigen::VectorXd gblCorrection(5);
+    Eigen::MatrixXd gblCovariance(5, 5);
+    Eigen::VectorXd gblResiduals(2);
+    Eigen::VectorXd gblErrorsMeasurements(2);
+    Eigen::VectorXd gblErrorsResiduals(2);
+    Eigen::VectorXd gblDownWeights(2);
+    unsigned int numData = 2;
+    for(auto plane : m_planes) {
+        traj.getScatResults(
+            plane.gblPos(), numData, gblResiduals, gblErrorsMeasurements, gblErrorsResiduals, gblDownWeights);
+        m_kink[plane.name()] = ROOT::Math::XYPoint(gblResiduals(0), gblResiduals(1));
+
+        traj.getResults(int(plane.gblPos()), localPar, localCov);
+        m_corrections[plane.name()] = ROOT::Math::XYZPoint(localPar(3), localPar(4), 0);
+        if(plane.hasCluster()) {
+            traj.getMeasResults(
+                plane.gblPos(), numData, gblResiduals, gblErrorsMeasurements, gblErrorsResiduals, gblDownWeights);
+            m_residual[plane.name()] = ROOT::Math::XYPoint(gblResiduals(0), gblResiduals(1));
+        }
+    }
+    m_isFitted = true;
+}
+
+ROOT::Math::XYZPoint GblTrack::intercept(double z) const {
+    // find the detector with largest z-positon <= z, assumes detectors sorted by z position
+    std::string layer = "";
+    bool found = false;
+
+    if(!m_isFitted) {
+        throw TrackError(typeid(GblTrack), "An interception is requested befor the track is fitted");
+    }
+    for(auto l : m_planes) {
+        layer = l.name();
+        if(l.postion() >= z) {
+            found = true;
+            break;
+        }
+    }
+    if(!found) {
+        throw TrackError(typeid(GblTrack), "Z-Position of " + std::to_string(z) + " is outside the telescopes z-coverage");
+    }
+    return (state(layer) + direction(layer) * (z - state(layer).z()));
+}
+
+ROOT::Math::XYZPoint GblTrack::state(std::string detectorID) const {
+    // The track state at any plane is the seed (always first cluster for now) plus the correction for the plane
+    // And as rotations are ignored, the z position is simply the detectors z postion
+    // Let's check first if the data is fitted and all components are there
+
+    if(!m_isFitted)
+        throw TrackError(typeid(GblTrack), " detector " + detectorID + " state is not defined before fitting");
+    if(m_materialBudget.count(detectorID) != 1) {
+        std::string list;
+        for(auto l : m_materialBudget)
+            list.append("\n " + l.first + " with <materialBudget, z position>(" + std::to_string(l.second.first) + ", " +
+                        std::to_string(l.second.second) + ")");
+        throw TrackError(typeid(GblTrack),
+                         " detector " + detectorID + " is not appearing in the material budget map" + list);
+    }
+    if(m_corrections.count(detectorID) != 1)
+        throw TrackError(typeid(GblTrack), " detector " + detectorID + " is not appearing in the corrections map");
+
+    // Using the global detector position here is of course not correct, it works for small/no rotations
+    // For larger rotations it is an issue
+    return ROOT::Math::XYZPoint(clusters().at(0)->global().x() + correction(detectorID).x(),
+                                clusters().at(0)->global().y() + correction(detectorID).y(),
+                                m_materialBudget.at(detectorID).second);
+}
+
+ROOT::Math::XYZVector GblTrack::direction(std::string detectorID) const {
+
+    // Defining the direction following the particle results in the direction
+    // beeing definded from the requested plane onwards to the next one
+    ROOT::Math::XYZPoint point = state(detectorID);
+
+    // searching for the next detector layer
+    bool found = false;
+    std::string nextLayer = "";
+    for(auto& layer : m_materialBudget) {
+        if(found) {
+            nextLayer = layer.first;
+            break;
+        } else if(layer.first == detectorID) {
+            found = true;
+        }
+    }
+    if(nextLayer == "")
+        throw TrackError(typeid(GblTrack), "Direction after the last telescope plane not defined");
+    ROOT::Math::XYZPoint pointAfter = state(nextLayer);
+    return ((pointAfter - point) / (pointAfter.z() - point.z()));
+}
+
+void GblTrack::print(std::ostream& out) const {
+    out << "GblTrack with nhits = " << m_trackClusters.size() << " and nscatterers = " << m_materialBudget.size()
+        << ", chi2 = " << m_chi2 << ", ndf = " << m_ndof;
+}

src/objects/GblTrack.hpp

+#ifndef GblTrack_H
+#define GblTrack_H 1
+
+#include "Track.hpp"
+namespace corryvreckan {
+    /**
+     * @ingroup Objects
+     * @brief GblTrack object
+     *
+     * This class is a general broken line track which knows how to fit itself. It is dervied from Track
+     */
+
+    class GblTrack : public Track {
+
+    public:
+        // Constructors and destructors
+        GblTrack();
+
+        // copy constructor
+        GblTrack(const GblTrack& track);
+
+        virtual GblTrack* clone() const override { return new GblTrack(*this); }
+
+        void print(std::ostream& out) const override;
+
+        /**
+         * @brief The fitting routine
+         */
+        void fit() override;
+
+        /**
+         * @brief Get the track position for a certain z position