New algorithm: TelescopeAnalysis for telescope reference plots

src/algorithms/TelescopeAnalysis/CMakeLists.txt

+# Define module and return the generated name as ALGORITHM_NAME
+CORRYVRECKAN_ALGORITHM(ALGORITHM_NAME)
+
+# Add source files to library
+CORRYVRECKAN_ALGORITHM_SOURCES(${ALGORITHM_NAME}
+    TelescopeAnalysis.cpp
+    # ADD SOURCE FILES HERE...
+)
+
+# Provide standard install target
+CORRYVRECKAN_ALGORITHM_INSTALL(${ALGORITHM_NAME})

src/algorithms/TelescopeAnalysis/README.md

+## TelescopeAnalysis
+**Maintainer**: Simon Spannagel (<simon.spannagel@cern.ch>)
+**Status**: Work in progress
+
+#### Description
+This algorithm produces reference histograms for the telescope performance used for tracking. It produces local and global biased residuals for each of the telescope planes, and if Monte Carlo information is available, calculates the residuals between track and Monte Carlo particle.
+
+Furthermore, the telescope resolution at the position of the DUT detector is plotted of Monte Carlo information is available. The Monte Carlo particle position is compared with the track interception with the DUT.
+
+#### Parameters
+* `chi2ndofCut`: Chi2 over number of degrees of freedom for the track to be taken into account. Tracks with a larger value are discarded. Default value is `3`.
+
+#### Plots produced
+* Telescope resolution at position of DUT
+
+For each detector participating in tracking, the following plots are produced:
+* Biased local track residual, for X and Y;
+* Biased global track residual, for X and Y;
+* Local residuals with track and Monte Carlo particle, for X and Y;
+* Global residuals with track and Monte Carlo particle, for X and Y;
+
+#### Usage
+```toml
+[NAME]
+chi2ndofCut = 3
+```
+Parameters to be used in multiple algorithms can also be defined globally at the top of the configuration file. This is highly encouraged for parameters such as `DUT` and `reference`.

src/algorithms/TelescopeAnalysis/TelescopeAnalysis.cpp

+#include "TelescopeAnalysis.h"
+#include "objects/Cluster.h"
+#include "objects/MCParticle.h"
+
+using namespace corryvreckan;
+using namespace std;
+
+TelescopeAnalysis::TelescopeAnalysis(Configuration config, std::vector<Detector*> detectors)
+    : Algorithm(std::move(config), std::move(detectors)) {
+
+    chi2ndofCut = m_config.get<double>("chi2ndofCut", 3.);
+}
+
+void TelescopeAnalysis::initialise() {
+
+    telescopeResolution = new TH1F("telescopeResolution", "telescopeResolution", 600, -0.2, 0.2);
+
+    // Initialise histograms per device
+    for(auto& detector : get_detectors()) {
+        std::string name = "residualX_local_" + detector->name();
+        telescopeResidualsLocalX[detector->name()] = new TH1F(name.c_str(), name.c_str(), 400, -0.2, 0.2);
+        name = "residualY_local_" + detector->name();
+        telescopeResidualsLocalY[detector->name()] = new TH1F(name.c_str(), name.c_str(), 400, -0.2, 0.2);
+        name = "residualX_global_" + detector->name();
+        telescopeResidualsX[detector->name()] = new TH1F(name.c_str(), name.c_str(), 400, -0.2, 0.2);
+        name = "residualY_global_" + detector->name();
+        telescopeResidualsY[detector->name()] = new TH1F(name.c_str(), name.c_str(), 400, -0.2, 0.2);
+
+        name = "residualX_MC_local_" + detector->name();
+        telescopeMCresidualsLocalX[detector->name()] = new TH1F(name.c_str(), name.c_str(), 400, -0.2, 0.2);
+        name = "residualY_MC_local_" + detector->name();
+        telescopeMCresidualsLocalY[detector->name()] = new TH1F(name.c_str(), name.c_str(), 400, -0.2, 0.2);
+        name = "residualX_MC_global_" + detector->name();
+        telescopeMCresidualsX[detector->name()] = new TH1F(name.c_str(), name.c_str(), 400, -0.2, 0.2);
+        name = "residualY_MC_global_" + detector->name();
+        telescopeMCresidualsY[detector->name()] = new TH1F(name.c_str(), name.c_str(), 400, -0.2, 0.2);
+    }
+}
+
+ROOT::Math::XYZPoint TelescopeAnalysis::closestApproach(ROOT::Math::XYZPoint position, MCParticles* particles) {
+    // Find the closest MC particle
+    double smallestDistance(DBL_MAX);
+    ROOT::Math::XYZPoint particlePosition;
+    for(auto& particle : (*particles)) {
+        ROOT::Math::XYZPoint entry = particle->getLocalStart();
+        ROOT::Math::XYZPoint exit = particle->getLocalEnd();
+        ROOT::Math::XYZPoint centre((entry.X() + exit.X()) / 2., (entry.Y() + exit.Y()) / 2., (entry.Z() + exit.Z()) / 2.);
+        double distance = sqrt((centre.X() - position.X()) * (centre.X() - position.X()) +
+                               (centre.Y() - position.Y()) * (centre.Y() - position.Y()));
+        if(distance < smallestDistance) {
+            particlePosition.SetXYZ(centre.X(), centre.Y(), centre.Z());
+        }
+    }
+    return particlePosition;
+}
+
+StatusCode TelescopeAnalysis::run(Clipboard* clipboard) {
+
+    // Get the tracks from the clipboard
+    Tracks* tracks = (Tracks*)clipboard->get("tracks");
+    if(tracks == NULL) {
+        LOG(DEBUG) << "No tracks on the clipboard";
+        return Success;
+    }
+
+    for(auto& track : (*tracks)) {
+
+        // Cut on the chi2/ndof
+        if(track->chi2ndof() > chi2ndofCut) {
+            continue;
+        }
+
+        // Loop over clusters of the track:
+        for(auto& cluster : track->clusters()) {
+            auto detector = get_detector(cluster->detectorID());
+            auto name = detector->name();
+
+            // Ignore DUT
+            if(name == m_config.get<std::string>("DUT")) {
+                continue;
+            }
+
+            // Get the MC particles from the clipboard
+            MCParticles* mcParticles = (MCParticles*)clipboard->get(name, "mcparticles");
+            if(mcParticles == NULL) {
+                continue;
+            }
+
+            ROOT::Math::XYZPoint intercept = track->intercept(cluster->globalZ());
+            auto interceptLocal = detector->globalToLocal(intercept);
+            ROOT::Math::XYZPoint particlePosition = closestApproach(cluster->local(), mcParticles);
+
+            telescopeResidualsLocalX[name]->Fill(cluster->localX() - interceptLocal.X());
+            telescopeResidualsLocalY[name]->Fill(cluster->localY() - interceptLocal.Y());
+            telescopeResidualsX[name]->Fill(cluster->globalX() - intercept.X());
+            telescopeResidualsY[name]->Fill(cluster->globalY() - intercept.Y());
+
+            telescopeMCresidualsLocalX[name]->Fill(cluster->localX() + detector->size().X() / 2 - particlePosition.X());
+            telescopeMCresidualsLocalY[name]->Fill(cluster->localY() + detector->size().Y() / 2 - particlePosition.Y());
+            telescopeMCresidualsX[name]->Fill(interceptLocal.X() + detector->size().X() / 2 - particlePosition.X());
+            telescopeMCresidualsY[name]->Fill(interceptLocal.Y() + detector->size().Y() / 2 - particlePosition.Y());
+        }
+
+        // Calculate telescope resolution at DUT
+        auto detector = get_detector(m_config.get<std::string>("DUT"));
+
+        // Get the MC particles from the clipboard
+        MCParticles* mcParticles = (MCParticles*)clipboard->get(detector->name(), "mcparticles");
+        if(mcParticles == NULL) {
+            continue;
+        }
+
+        auto intercept = detector->getIntercept(track);
+        auto interceptLocal = detector->globalToLocal(intercept);
+        auto particlePosition = closestApproach(interceptLocal, mcParticles);
+
+        telescopeResolution->Fill(interceptLocal.X() + detector->size().X() / 2 - particlePosition.X());
+    }
+
+    // Return value telling analysis to keep running
+    return Success;
+}

src/algorithms/TelescopeAnalysis/TelescopeAnalysis.h

+#ifndef TelescopeAnalysis_H
+#define TelescopeAnalysis_H 1
+
+#include <iostream>
+#include "TH1F.h"
+#include "core/algorithm/Algorithm.h"
+#include "objects/MCParticle.h"
+#include "objects/Track.h"
+
+namespace corryvreckan {
+    /** @ingroup Algorithms
+     */
+    class TelescopeAnalysis : public Algorithm {
+
+    public:
+        // Constructors and destructors
+        TelescopeAnalysis(Configuration config, std::vector<Detector*> detectors);
+        ~TelescopeAnalysis() {}
+
+        // Functions
+        void initialise();
+        StatusCode run(Clipboard* clipboard);
+        void finalise(){};
+
+    private:
+        ROOT::Math::XYZPoint closestApproach(ROOT::Math::XYZPoint position, MCParticles* particles);
+
+        // Histograms for each of the devices
+        std::map<std::string, TH1F*> telescopeMCresidualsLocalX;
+        std::map<std::string, TH1F*> telescopeMCresidualsLocalY;
+        std::map<std::string, TH1F*> telescopeMCresidualsX;
+        std::map<std::string, TH1F*> telescopeMCresidualsY;
+
+        std::map<std::string, TH1F*> telescopeResidualsLocalX;
+        std::map<std::string, TH1F*> telescopeResidualsLocalY;
+        std::map<std::string, TH1F*> telescopeResidualsX;
+        std::map<std::string, TH1F*> telescopeResidualsY;
+
+        // Histograms at the position of the DUT
+        TH1F* telescopeResolution;
+
+        // Parameters
+        double chi2ndofCut;
+    };
+} // namespace corryvreckan
+#endif // TelescopeAnalysis_H