Merge branch 'master' of ssh://gitlab.cern.ch:7999/corryvreckan/corryvreckan into fix_magnus

doc/usermanual/references.bib

@@ -301,3 +301,9 @@ keywords = "Simulation, Silicon detectors, Geant4, TCAD, Drift–diffusion"
       url           = "https://cds.cern.ch/record/2649493",
       note          = {accessed 11~1019}
 }
+@online{root-tefficiency-class-ref,
+    title = {ROOT TEfficiency Class Reference},
+    author = {},
+    url = {https://root.cern.ch/doc/master/classTEfficiency.html#ae80c3189bac22b7ad15f57a1476ef75b},
+    note = {accessed 01~2020},
+    }

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/AnalysisEfficiency/README.md

@@ -8,6 +8,13 @@
 This module measures the efficiency of the DUT by comparing its cluster positions with the interpolated track position at the DUT.
 It also comprises a range of histograms to investigate where inefficiencies might come from.
 
+The efficiency is calculated as the fraction of tracks with associated clusters on the DUT over the the total number of tracks intersecting the DUT (or region-of-interest, if defined).
+It is stored in a ROOT `TEfficiency` object (see below).
+Its uncertainty is calculated using the default ROOT `TEfficiency` method which is applying a Clopper-Pearson confidence interval of one sigma.
+Analogue to a Gaussian sigma, this corresponds to the central 68.3% of a binomial distribution for the given efficiency but taking into account a lower limit of 0 and an upper limit of 1.
+This method is recommended by the Particle Data Group.
+More information can be found in the ROOT `TEfficiency` class reference, section `ClopperPearson()` @root-tefficiency-class-ref.
+
 ### Parameters
 * `time_cut_frameedge`: Parameter to discard telescope tracks at the frame edges (start and end of the current event window). Defaults to `20ns`.
 * `chi2ndof_cut`: Acceptance criterion for telescope tracks, defaults to a value of `3`.
@@ -22,7 +29,7 @@ For the DUT, the following plots are produced:
   * 2D Maps of chip efficiency in local and global coordinates, filled at the position of the track intercept point or at the position of the associated cluster center
   * 2D Maps of the position difference of a track with and without associated cluster to the previous track
   * 2D Map of the distance between track intersection and associated cluster
-  
+
 * 1D histograms:
   * Histogram of all single-pixel efficiencies
   * Histograms of time difference of the matched and non-matched track time to the previous track
@@ -37,3 +44,4 @@ For the DUT, the following plots are produced:
 [AnalysisEfficiency]
 chi2ndof_cut = 5
 ```
+[@root-tefficiency-class-ref]: https://root.cern.ch/doc/master/classTEfficiency.html#ae80c3189bac22b7ad15f57a1476ef75b

src/modules/EventLoaderEUDAQ2/EventLoaderEUDAQ2.cpp

@@ -120,7 +120,7 @@ void EventLoaderEUDAQ2::initialise() {
                          2.1e4,
                          -10,
                          200);
-            std::string histTitle = "hPixelTriggerTimeResidualOverTime_0;time [us];pixel_ts - trigger_ts [us];# entries";
+            std::string histTitle = "hPixelTriggerTimeResidualOverTime_0;time [s];pixel_ts - trigger_ts [us];# entries";
             hPixelTriggerTimeResidualOverTime =
                 new TH2D("hPixelTriggerTimeResidualOverTime_0", histTitle.c_str(), 3e3, 0, 3e3, 1e4, -50, 50);
         }

src/modules/EventLoaderTimepix3/EventLoaderTimepix3.cpp

@@ -486,7 +486,7 @@ bool EventLoaderTimepix3::loadData(std::shared_ptr<Clipboard> clipboard,
         if(header == 0x6) {
             const UChar_t header2 = ((pixdata & 0x0F00000000000000) >> 56) & 0xF;
             if(header2 == 0xF) {
-                unsigned long long int stamp = (pixdata & 0x1E0) >> 5;
+                unsigned int stamp = (pixdata & 0x1E0) >> 5;
                 long long int timestamp_raw = static_cast<long long int>(pixdata & 0xFFFFFFFFE00) >> 9;
                 long long int timestamp = 0;
                 // int triggerNumber = ((pixdata & 0xFFF00000000000) >> 44);
@@ -502,7 +502,7 @@ bool EventLoaderTimepix3::loadData(std::shared_ptr<Clipboard> clipboard,
                 timestamp = timestamp_raw + (static_cast<long long int>(m_TDCoverflowCounter) << 35);
 
                 double triggerTime =
-                    static_cast<double>(timestamp + static_cast<long long int>(stamp) / 12) / (8. * 0.04); // 320 MHz clock
+                    (static_cast<double>(timestamp) + static_cast<double>(stamp) / 12) / (8. * 0.04); // 320 MHz clock
                 SpidrSignal* triggerSignal = new SpidrSignal("trigger", triggerTime);
                 spidrData->push_back(triggerSignal);
             }

src/modules/OnlineMonitor/README.md

@@ -5,7 +5,7 @@
 
 ### Description
 This module opens a GUI to monitor the progress of the reconstruction.
-Since Linux allows concurrent (reading) file access, this can already be started while a run is recorded to disk and thus serves as online monitoring tool during data taking.
+Since Linux allows concurrent (reading) file access, this can already be started while a run is recorded to disk and thus serves as an online monitoring tool during data taking.
 A set of canvases is available to display a variety of information ranging from hitmaps and basic correlation plots to more advances results such as tracking quality and track angles.
 The plots on each of the canvases contain real time data, automatically updated every `update` events.
 
@@ -13,7 +13,7 @@ The displayed plots and their source can be configured via the framework configu
 Here, each canvas is configured via a matrix containing the path of the plot and its plotting options in each row, e.g.
 
 ```toml
-DUTPlots =  [["EventLoaderEUDAQ2/%DUT%/hitmap", "colz"],
+dut_plots =  [["EventLoaderEUDAQ2/%DUT%/hitmap", "colz"],
              ["EventLoaderEUDAQ2/%DUT%/hPixelRawValues", "log"]]
 ```
 
@@ -23,7 +23,7 @@ In addition, the `log` keyword is supported, which switches the Y axis of the re
 Several keywords can be used in the plot path, which are parsed and interpreted by the OnlineMonitor module:
 
 * `%DETECTOR%`: If this keyword is found, the plot is added for each of the available detectors by replacing the keyword with the respective detector name.
-* `%DUT%`: This keyword is replaced by the vale of the `DUT` configuration key of the framework.
+* `%DUT%`: This keyword is replaced by the value of the `DUT` configuration key of the framework.
 * `%REFERENCE%`: This keyword is replaced by the vale of the `reference` configuration key of the framework.
 
 The "corryvreckan" namespace is not required to be added to the plot path.
@@ -35,8 +35,8 @@ The "corryvreckan" namespace is not required to be added to the plot path.
 
 
 * `overview`: List of plots to be placed on the "Overview" canvas of the online monitor. The list of plots created in the default configuration is listed below.
-* `dut_plots`: List of plots to be placed on the "DUTPlots" canvas of the online monitor. By default, this canvas contains plots collected from the `EventLoaderEUDAQ2` as well as the `AnalysisDUT` modules for the each configured DUT. This canvas should be customized for the respective DUT.
-* `hitmaps`: List of plots to be placed on the "HitMaps" canvas of the online monitor. By default, this canvas displays `Correlations/%DETECTOR%/hitmap` for all detectors.
+* `dut_plots`: List of plots to be placed on the "DUTs" canvas of the online monitor. By default, this canvas contains plots collected from the `EventLoaderEUDAQ2` as well as the `AnalysisDUT` modules for the each configured DUT. This canvas should be customized for the respective DUT.
+* `hitmaps`: List of plots to be placed on the "Hitmaps" canvas of the online monitor. By default, this canvas displays `Correlations/%DETECTOR%/hitmap` for all detectors.
 * `tracking`: List of plots to be placed on the "Tracking" canvas of the online monitor. The list of plots created in the default configuration is listed below.
 * `residuals`: List of plots to be placed on the "Residuals" canvas of the online monitor. By default, this canvas displays `Tracking4D/%DETECTOR%/residualsX` for all detectors.
 * `correlation_x`: List of plots to be placed on the "CorrelationX" canvas of the online monitor.  By default, this canvas displays `Correlations/%DETECTOR%/correlationX` for all detectors.

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;
+}