Merge branch 'master' into documentation

src/core/config/Configuration.hpp

@@ -165,6 +165,13 @@ namespace corryvreckan {
         // TODO [doc] Provide second template parameter to specify the vector type to return it in
         template <typename T> void setArray(const std::string& key, const std::vector<T>& val);
 
+        /**
+         * @brief Set matrix of values for a key in a given type
+         * @param key Key to set values of
+         * @param val List of values to assign to the key
+         */
+        template <typename T> void setMatrix(const std::string& key, const Matrix<T>& val);
+
         /**
          * @brief Set default value for a key only if it is not defined yet
          * @param key Key to possible set value of

src/core/config/Configuration.tpp

@@ -90,6 +90,10 @@ namespace corryvreckan {
             Matrix<T> matrix;
             auto node = parse_value(str);
             for(auto& child : node->children) {
+                if(child->children.empty()) {
+                    throw std::invalid_argument("matrix has less than two dimensions");
+                }
+
                 std::vector<T> array;
                 // Create subarray of matrix
                 for(auto& subchild : child->children) {
@@ -99,10 +103,6 @@ namespace corryvreckan {
                         throw InvalidKeyError(key, getName(), subchild->value, typeid(T), e.what());
                     }
                 }
-                if(!child->value.empty()) {
-                    throw std::invalid_argument("matrix has less than two dimensions");
-                }
-
                 matrix.push_back(array);
             }
             return matrix;
@@ -154,6 +154,27 @@ namespace corryvreckan {
         config_[key] = str;
     }
 
+    template <typename T> void Configuration::setMatrix(const std::string& key, const Matrix<T>& val) {
+        // NOTE: not the most elegant way to support arrays
+        if(val.empty()) {
+            return;
+        }
+
+        std::string str = "[";
+        for(auto& col : val) {
+            str += "[";
+            for(auto& el : col) {
+                str += corryvreckan::to_string(el);
+                str += ",";
+            }
+            str.pop_back();
+            str += "],";
+        }
+        str.pop_back();
+        str += "]";
+        config_[key] = str;
+    }
+
     template <typename T> void Configuration::setDefault(const std::string& key, const T& val) {
         if(!has(key)) {
             set<T>(key, val);

src/core/detector/Detector.cpp

@@ -44,6 +44,7 @@ Detector::Detector(const Configuration& config) {
     m_nPixelsX = npixels.x();
     m_nPixelsY = npixels.y();
     m_timingOffset = config.get<double>("time_offset", 0.0);
+    m_roi = config.getMatrix<int>("roi", std::vector<std::vector<int>>());
 
     this->initialise();
 
@@ -172,11 +173,13 @@ Configuration Detector::getConfiguration() {
         config.set("mask_file", m_maskfile_name);
     }
 
+    config.setMatrix("roi", m_roi);
+
     return config;
 }
 
 // Function to get global intercept with a track
-PositionVector3D<Cartesian3D<double>> Detector::getIntercept(Track* track) {
+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();
@@ -192,8 +195,12 @@ PositionVector3D<Cartesian3D<double>> Detector::getIntercept(Track* track) {
     return globalIntercept;
 }
 
+PositionVector3D<Cartesian3D<double>> Detector::getLocalIntercept(const Track* track) {
+    return globalToLocal(getIntercept(track));
+}
+
 // Function to check if a track intercepts with a plane
-bool Detector::hasIntercept(Track* track, double pixelTolerance) {
+bool Detector::hasIntercept(const Track* track, double pixelTolerance) {
 
     // First, get the track intercept in global co-ordinates with the plane
     PositionVector3D<Cartesian3D<double>> globalIntercept = this->getIntercept(track);
@@ -240,11 +247,11 @@ bool Detector::hitMasked(Track* track, int tolerance) {
 }
 
 // Functions to get row and column from local position
-double Detector::getRow(PositionVector3D<Cartesian3D<double>> localPosition) {
+double Detector::getRow(const PositionVector3D<Cartesian3D<double>> localPosition) {
     double row = ((localPosition.Y() + m_pitch.Y() / 2.) / m_pitch.Y()) + m_nPixelsY / 2.;
     return row;
 }
-double Detector::getColumn(PositionVector3D<Cartesian3D<double>> localPosition) {
+double Detector::getColumn(const PositionVector3D<Cartesian3D<double>> localPosition) {
     double column = ((localPosition.X() + m_pitch.X() / 2.) / m_pitch.X()) + m_nPixelsX / 2.;
     return column;
 }
@@ -257,13 +264,91 @@ PositionVector3D<Cartesian3D<double>> Detector::getLocalPosition(double row, dou
 }
 
 // Function to get in-pixel position
-double Detector::inPixelX(PositionVector3D<Cartesian3D<double>> localPosition) {
+double Detector::inPixelX(const PositionVector3D<Cartesian3D<double>> localPosition) {
     double column = getColumn(localPosition);
     double inPixelX = m_pitch.X() * (column - floor(column));
     return inPixelX;
 }
-double Detector::inPixelY(PositionVector3D<Cartesian3D<double>> localPosition) {
+double Detector::inPixelY(const PositionVector3D<Cartesian3D<double>> localPosition) {
     double row = getRow(localPosition);
     double inPixelY = m_pitch.Y() * (row - floor(row));
     return inPixelY;
 }
+
+// Check if track position is within ROI:
+bool Detector::isWithinROI(const Track* track) {
+
+    // Check that track is within region of interest using winding number algorithm
+    auto localIntercept = this->getLocalIntercept(track);
+    auto coordinates = std::make_pair(this->getColumn(localIntercept), this->getRow(localIntercept));
+    if(winding_number(coordinates, m_roi) != 0) {
+        return true;
+    }
+
+    // Outside ROI:
+    return false;
+}
+
+// Check if cluster is within ROI and/or touches ROI border:
+bool Detector::isWithinROI(const Cluster* cluster) {}
+
+/* isLeft(): tests if a point is Left|On|Right of an infinite line.
+ * via: http://geomalgorithms.com/a03-_inclusion.html
+ *    Input:  three points P0, P1, and P2
+ *    Return: >0 for P2 left of the line through P0 and P1
+ *            =0 for P2  on the line
+ *            <0 for P2  right of the line
+ *    See: Algorithm 1 "Area of Triangles and Polygons"
+ */
+int Detector::isLeft(std::pair<int, int> pt0, std::pair<int, int> pt1, std::pair<int, int> pt2) {
+    return ((pt1.first - pt0.first) * (pt2.second - pt0.second) - (pt2.first - pt0.first) * (pt1.second - pt0.second));
+}
+
+/* Winding number test for a point in a polygon
+ * via: http://geomalgorithms.com/a03-_inclusion.html
+ *      Input:   x, y = a point,
+ *               polygon = vector of vertex points of a polygon V[n+1] with V[n]=V[0]
+ *      Return:  wn = the winding number (=0 only when P is outside)
+ */
+int Detector::winding_number(std::pair<int, int> probe, std::vector<std::vector<int>> polygon) {
+    // Two points don't make an area
+    if(polygon.size() < 3) {
+        LOG(DEBUG) << "No ROI given.";
+        return 0;
+    }
+
+    int wn = 0; // the  winding number counter
+
+    // loop through all edges of the polygon
+
+    // edge from V[i] to  V[i+1]
+    for(int i = 0; i < polygon.size(); i++) {
+        auto point_this = std::make_pair(polygon.at(i).at(0), polygon.at(i).at(1));
+        auto point_next = (i + 1 < polygon.size() ? std::make_pair(polygon.at(i + 1).at(0), polygon.at(i + 1).at(1))
+                                                  : std::make_pair(polygon.at(0).at(0), polygon.at(0).at(1)));
+
+        // start y <= P.y
+        if(point_this.second <= probe.second) {
+            // an upward crossing
+            if(point_next.second > probe.second) {
+                // P left of  edge
+                if(isLeft(point_this, point_next, probe) > 0) {
+                    // have  a valid up intersect
+                    ++wn;
+                }
+            }
+        } else {
+            // start y > P.y (no test needed)
+
+            // a downward crossing
+            if(point_next.second <= probe.second) {
+                // P right of  edge
+                if(isLeft(point_this, point_next, probe) < 0) {
+                    // have  a valid down intersect
+                    --wn;
+                }
+            }
+        }
+    }
+    return wn;
+}

src/core/detector/Detector.hpp

@@ -79,10 +79,12 @@ namespace corryvreckan {
         void update();
 
         // Function to get global intercept with a track
-        PositionVector3D<Cartesian3D<double>> getIntercept(Track* track);
+        PositionVector3D<Cartesian3D<double>> getIntercept(const Track* track);
+        // Function to get local intercept with a track
+        PositionVector3D<Cartesian3D<double>> getLocalIntercept(const Track* track);
 
         // Function to check if a track intercepts with a plane
-        bool hasIntercept(Track* track, double pixelTolerance = 0.);
+        bool hasIntercept(const Track* track, double pixelTolerance = 0.);
 
         // Function to check if a track goes through/near a masked pixel
         bool hitMasked(Track* track, int tolerance = 0.);
@@ -101,6 +103,9 @@ namespace corryvreckan {
         ROOT::Math::XYZPoint localToGlobal(ROOT::Math::XYZPoint local) { return m_localToGlobal * local; };
         ROOT::Math::XYZPoint globalToLocal(ROOT::Math::XYZPoint global) { return m_globalToLocal * global; };
 
+        bool isWithinROI(const Track* track);
+        bool isWithinROI(const Cluster* cluster);
+
     private:
         // Member variables
         // Detector information
@@ -111,6 +116,10 @@ namespace corryvreckan {
         int m_nPixelsY;
         double m_timingOffset;
 
+        std::vector<std::vector<int>> m_roi;
+        static int winding_number(std::pair<int, int> probe, std::vector<std::vector<int>> polygon);
+        inline static int isLeft(std::pair<int, int> pt0, std::pair<int, int> pt1, std::pair<int, int> pt2);
+
         // Displacement and rotation in x,y,z
         ROOT::Math::XYZPoint m_displacement;
         ROOT::Math::XYZVector m_orientation;

src/modules/CLICpix2Analysis/CLICpix2Analysis.cpp

@@ -12,6 +12,7 @@ CLICpix2Analysis::CLICpix2Analysis(Configuration config, std::vector<Detector*>
     m_DUT = m_config.get<std::string>("DUT");
 
     m_timeCutFrameEdge = m_config.get<double>("timeCutFrameEdge", Units::convert(20, "ns"));
+    m_roi = m_config.getMatrix<int>("roi", std::vector<std::vector<int>>());
 
     spatialCut = m_config.get<double>("spatialCut", Units::convert(50, "um"));
     chi2ndofCut = m_config.get<double>("chi2ndofCut", 3.);
@@ -49,6 +50,8 @@ void CLICpix2Analysis::initialise() {
     // Per-pixel histograms
     hHitMapAssoc =
         new TH2F("hitMapAssoc", "hitMapAssoc", det->nPixelsX(), 0, det->nPixelsX(), det->nPixelsY(), 0, det->nPixelsY());
+    hHitMapROI =
+        new TH2F("hitMapROI", "hitMapROI", det->nPixelsX(), 0, det->nPixelsX(), det->nPixelsY(), 0, det->nPixelsY());
     hPixelToTAssoc = new TH1F("pixelToTAssoc", "pixelToTAssoc", 32, 0, 31);
     hPixelToTMapAssoc = new TProfile2D("pixelToTMapAssoc",
                                        "pixelToTMapAssoc",
@@ -71,19 +74,50 @@ void CLICpix2Analysis::initialise() {
     residualsY2pix = new TH1F("residualsY2pix", "residualsY2pix", 400, -0.2, 0.2);
 
     clusterTotAssoc = new TH1F("clusterTotAssociated", "clusterTotAssociated", 10000, 0, 10000);
+    clusterTotAssocNorm = new TH1F("clusterTotAssociatedNormalized", "clusterTotAssociatedNormalized", 10000, 0, 10000);
     clusterSizeAssoc = new TH1F("clusterSizeAssociated", "clusterSizeAssociated", 30, 0, 30);
 
+    // In-pixel studies:
+    auto pitch_x = Units::convert(det->pitch().X(), "um");
+    auto pitch_y = Units::convert(det->pitch().Y(), "um");
+    auto mod_axes = "x_{track} mod " + std::to_string(pitch_x) + "#mum;y_{track} mod " + std::to_string(pitch_y) + "#mum;";
+
+    std::string title = "DUT x resolution;" + mod_axes + "MAD(#Deltax) [#mum]";
+    rmsxvsxmym = new TProfile2D("rmsxvsxmym", title.c_str(), pitch_x, 0, pitch_x, pitch_y, 0, pitch_y);
+
+    title = "DUT y resolution;" + mod_axes + "MAD(#Deltay) [#mum]";
+    rmsyvsxmym = new TProfile2D("rmsyvsxmym", title.c_str(), pitch_x, 0, pitch_x, pitch_y, 0, pitch_y);
+
+    title = "DUT resolution;" + mod_axes + "MAD(#sqrt{#Deltax^{2}+#Deltay^{2}}) [#mum]";
+    rmsxyvsxmym = new TProfile2D("rmsxyvsxmym", title.c_str(), pitch_x, 0, pitch_x, pitch_y, 0, pitch_y);
+
+    title = "DUT cluster charge map;" + mod_axes + "<cluster charge> [ke]";
+    qvsxmym = new TProfile2D("qvsxmym", title.c_str(), pitch_x, 0, pitch_x, pitch_y, 0, pitch_y, 0, 250);
+
+    title = "DUT cluster charge map, Moyal approx;" + mod_axes + "cluster charge MPV [ke]";
+    qMoyalvsxmym = new TProfile2D("qMoyalvsxmym", title.c_str(), pitch_x, 0, pitch_x, pitch_y, 0, pitch_y, 0, 250);
+
+    title = "DUT seed pixel charge map;" + mod_axes + "<seed pixel charge> [ke]";
+    pxqvsxmym = new TProfile2D("pxqvsxmym", title.c_str(), pitch_x, 0, pitch_x, pitch_y, 0, pitch_y, 0, 250);
+
+    title = "DUT cluster size map;" + mod_axes + "<pixels/cluster>";
+    npxvsxmym = new TProfile2D("npxvsxmym", title.c_str(), pitch_x, 0, pitch_x, pitch_y, 0, pitch_y, 0, 4.5);
+
+    title = "DUT 1-pixel cluster map;" + mod_axes + "clusters";
+    npx1vsxmym = new TH2F("npx1vsxmym", title.c_str(), pitch_x, 0, pitch_x, pitch_y, 0, pitch_y);
+
+    title = "DUT 2-pixel cluster map;" + mod_axes + "clusters";
+    npx2vsxmym = new TH2F("npx2vsxmym", title.c_str(), pitch_x, 0, pitch_x, pitch_y, 0, pitch_y);
+
+    title = "DUT 3-pixel cluster map;" + mod_axes + "clusters";
+    npx3vsxmym = new TH2F("npx3vsxmym", title.c_str(), pitch_x, 0, pitch_x, pitch_y, 0, pitch_y);
+
+    title = "DUT 4-pixel cluster map;" + mod_axes + "clusters";
+    npx4vsxmym = new TH2F("npx4vsxmym", title.c_str(), pitch_x, 0, pitch_x, pitch_y, 0, pitch_y);
+
     // Efficiency maps
-    hPixelEfficiencyMap = new TProfile2D("hPixelEfficiencyMap",
-                                         "hPixelEfficiencyMap",
-                                         Units::convert(det->pitch().X(), "um"),
-                                         0,
-                                         Units::convert(det->pitch().X(), "um"),
-                                         Units::convert(det->pitch().Y(), "um"),
-                                         0,
-                                         Units::convert(det->pitch().Y(), "um"),
-                                         0,
-                                         1);
+    hPixelEfficiencyMap =
+        new TProfile2D("hPixelEfficiencyMap", "hPixelEfficiencyMap", pitch_x, 0, pitch_x, pitch_y, 0, pitch_y, 0, 1);
     hChipEfficiencyMap = new TProfile2D("hChipEfficiencyMap",
                                         "hChipEfficiencyMap",
                                         det->nPixelsX(),
@@ -134,25 +168,39 @@ StatusCode CLICpix2Analysis::run(Clipboard* clipboard) {
 
     // Loop over all tracks
     for(auto& track : (*tracks)) {
-        // Flag for efficiency calculation
-        bool cluster_associated = false;
+        // Flags to select clusters and tracks
+        bool has_associated_cluster = false;
+        bool is_within_roi = true;
 
         LOG(DEBUG) << "Looking at next track";
 
         // Cut on the chi2/ndof
         if(track->chi2ndof() > chi2ndofCut) {
+            LOG(DEBUG) << " - track discarded due to Chi2/ndof";
             continue;
         }
 
         // Check if it intercepts the DUT
-        PositionVector3D<Cartesian3D<double>> globalIntercept = detector->getIntercept(track);
+        auto globalIntercept = detector->getIntercept(track);
+        auto localIntercept = detector->globalToLocal(globalIntercept);
+
         if(!detector->hasIntercept(track, 1.)) {
             LOG(DEBUG) << " - track outside DUT area";
             continue;
         }
 
+        // Check that track is within region of interest using winding number algorithm
+        if(!track->is_within_roi()) {
+            is_within_roi = false;
+        }
+
+        // Check that it doesn't go through/near a masked pixel
+        if(detector->hitMasked(track, 1.)) {
+            LOG(DEBUG) << " - track close to masked pixel";
+            continue;
+        }
+
         // Discard tracks which are very close to the frame edges
-        //
         if(fabs(track->timestamp() - clipboard->get_persistent("currentTime")) < m_timeCutFrameEdge) {
             // Late edge - the currentTime has been updated by Timexpi3EventLoader to point to the end of the frame`
             LOG(DEBUG) << " - track close to end of readout frame: "
@@ -170,13 +218,9 @@ StatusCode CLICpix2Analysis::run(Clipboard* clipboard) {
             continue;
         }
 
-        // Check that it doesn't go through/near a masked pixel
-        if(detector->hitMasked(track, 1.)) {
-            LOG(DEBUG) << " - track close to masked pixel";
-            continue;
-        }
-
-        // FIXME check that track is within fiducial area
+        // Calculate in-pixel position of track in microns
+        double xmod = Units::convert(detector->inPixelX(localIntercept), "um");
+        double ymod = Units::convert(detector->inPixelY(localIntercept), "um");
 
         // Get the DUT clusters from the clipboard
         Clusters* clusters = (Clusters*)clipboard->get(m_DUT, "clusters");
@@ -210,7 +254,7 @@ StatusCode CLICpix2Analysis::run(Clipboard* clipboard) {
                            << Units::display(yabsdistance, {"um"}) << ")";
 
                 // We now have an associated cluster
-                cluster_associated = true;
+                has_associated_cluster = true;
                 // FIXME need to understand local coord of clusters - why shifted? what's normal?
                 auto clusterLocal = detector->globalToLocal(cluster->global());
                 hClusterMapAssoc->Fill(detector->getColumn(clusterLocal), detector->getRow(clusterLocal));
@@ -220,9 +264,17 @@ StatusCode CLICpix2Analysis::run(Clipboard* clipboard) {
 
                 clusterTotAssoc->Fill(cluster->tot());
 
+                // Cluster charge normalized to path length in sensor:
+                double norm = 1; // FIXME fabs(cos( turn*wt )) * fabs(cos( tilt*wt ));
+                auto normalized_charge = cluster->tot() * norm;
+                clusterTotAssocNorm->Fill(normalized_charge);
+
                 // Fill per-pixel histograms
                 for(auto& pixel : (*cluster->pixels())) {
                     hHitMapAssoc->Fill(pixel->column(), pixel->row());
+                    if(is_within_roi) {
+                        hHitMapROI->Fill(pixel->column(), pixel->row());
+                    }
                     hPixelToTAssoc->Fill(pixel->tot());
                     hPixelToTMapAssoc->Fill(pixel->column(), pixel->row(), pixel->tot());
                 }
@@ -244,6 +296,31 @@ StatusCode CLICpix2Analysis::run(Clipboard* clipboard) {
 
                 clusterSizeAssoc->Fill(cluster->size());
 
+                // Fill in-pixel plots: (all as function of track position within pixel cell)
+                if(is_within_roi) {
+                    qvsxmym->Fill(xmod, ymod, cluster->tot());                     // cluster charge profile
+                    qMoyalvsxmym->Fill(xmod, ymod, exp(-normalized_charge / 3.5)); // norm. cluster charge profile
+
+                    // mean charge of cluster seed
+                    pxqvsxmym->Fill(xmod, ymod, cluster->getSeedPixel()->charge());
+
+                    // mean cluster size
+                    npxvsxmym->Fill(xmod, ymod, cluster->size());
+                    if(cluster->size() == 1)
+                        npx1vsxmym->Fill(xmod, ymod);
+                    if(cluster->size() == 2)
+                        npx2vsxmym->Fill(xmod, ymod);
+                    if(cluster->size() == 3)
+                        npx3vsxmym->Fill(xmod, ymod);
+                    if(cluster->size() == 4)
+                        npx4vsxmym->Fill(xmod, ymod);
+
+                    // residual MAD x, y, combined (sqrt(x*x + y*y))
+                    rmsxvsxmym->Fill(xmod, ymod, xabsdistance);
+                    rmsyvsxmym->Fill(xmod, ymod, yabsdistance);
+                    rmsxyvsxmym->Fill(xmod, ymod, fabs(sqrt(xdistance * xdistance + ydistance * ydistance)));
+                }
+
                 track->addAssociatedCluster(cluster);
                 hAssociatedTracksGlobalPosition->Fill(globalIntercept.X(), globalIntercept.Y());
 
@@ -253,13 +330,14 @@ StatusCode CLICpix2Analysis::run(Clipboard* clipboard) {
         }
 
         // Efficiency plots:
-        hGlobalEfficiencyMap->Fill(globalIntercept.X(), globalIntercept.Y(), cluster_associated);
+        hGlobalEfficiencyMap->Fill(globalIntercept.X(), globalIntercept.Y(), has_associated_cluster);
+        hChipEfficiencyMap->Fill(
+            detector->getColumn(localIntercept), detector->getRow(localIntercept), has_associated_cluster);
 
-        auto localIntercept = detector->globalToLocal(globalIntercept);
-        hChipEfficiencyMap->Fill(detector->getColumn(localIntercept), detector->getRow(localIntercept), cluster_associated);
-        hPixelEfficiencyMap->Fill(Units::convert(detector->inPixelX(localIntercept), "um"),
-                                  Units::convert(detector->inPixelY(localIntercept), "um"),
-                                  cluster_associated);
+        // For pixels, only look at the ROI:
+        if(is_within_roi) {
+            hPixelEfficiencyMap->Fill(xmod, ymod, has_associated_cluster);
+        }
     }
 
     // Return value telling analysis to keep running

src/modules/CLICpix2Analysis/CLICpix2Analysis.h

@@ -24,7 +24,7 @@ namespace corryvreckan {
 
     private:
         // Histograms
-        TH2F *hClusterMapAssoc, *hHitMapAssoc;
+        TH2F *hClusterMapAssoc, *hHitMapAssoc, *hHitMapROI;
         TProfile2D *hClusterSizeMapAssoc, *hClusterToTMapAssoc;
 
         TH1F* hPixelToTAssoc;
@@ -36,9 +36,14 @@ namespace corryvreckan {
         TH1F *residualsX1pix, *residualsY1pix;
         TH1F *residualsX2pix, *residualsY2pix;
 
-        TH1F* clusterTotAssoc;
+        TH1F *clusterTotAssoc, *clusterTotAssocNorm;
         TH1F* clusterSizeAssoc;
 
+        TProfile2D *rmsxvsxmym, *rmsyvsxmym, *rmsxyvsxmym;
+        TProfile2D *qvsxmym, *qMoyalvsxmym, *pxqvsxmym;
+        TProfile2D* npxvsxmym;
+        TH2F *npx1vsxmym, *npx2vsxmym, *npx3vsxmym, *npx4vsxmym;
+
         TProfile2D* hPixelEfficiencyMap;
         TProfile2D* hChipEfficiencyMap;
         TProfile2D* hGlobalEfficiencyMap;
@@ -57,6 +62,7 @@ namespace corryvreckan {
         std::string m_DUT;
         double spatialCut, m_timeCutFrameEdge;
         double chi2ndofCut;
+        std::vector<std::vector<int>> m_roi;
     };
 } // namespace corryvreckan
 

src/modules/CLICpix2Analysis/README.md

@@ -15,6 +15,7 @@ Analysis module for CLICpix2 prototypes. This module is still work in progress,
 * 2D Map of cluster sizes for associated clusters
 * 2D Map of cluster ToT values from associated clusters
 * 2D Map of associated hits
+* 2D Map of associated hits within the defined region-of-interest
 * Distribution of pixel ToT values from associated clusters
 * 2D Map of pixel ToT values from associated clusters
 * Track residuals in X and Y

src/modules/Clicpix2EventLoader/Clicpix2EventLoader.cpp

@@ -28,9 +28,9 @@ void Clicpix2EventLoader::initialise() {
     // Open the root directory
     DIR* directory = opendir(inputDirectory.c_str());
     if(directory == NULL) {
-        LOG(ERROR) << "Directory " << inputDirectory << " does not exist";
-        return;
+        throw ModuleError("Directory " + inputDirectory + " does not exist");
     }
+
     dirent* entry;
     dirent* file;
 
@@ -52,8 +52,7 @@ void Clicpix2EventLoader::initialise() {
     }
 
     if(m_matrix.empty()) {
-        LOG(ERROR) << "No matrix configuration file found in " << inputDirectory;
-        return;
+        throw ModuleError("No matrix configuration file found in " + inputDirectory);
     }
 
     // Read the matrix configuration:

src/modules/ClicpixAnalysis/ClicpixAnalysis.cpp

@@ -409,8 +409,9 @@ StatusCode ClicpixAnalysis::run(Clipboard* clipboard) {
 
           // Get the pixel global position
           LOG(TRACE) <<"New pixel, row = "<<(*itPixel)->row()<<", column = "<<(*itPixel)->column();
-          PositionVector3D<Cartesian3D<double> > pixelPositionLocal = get_detector(dutID)->getLocalPosition((*itPixel)->row(),(*itPixel)->column());
-          PositionVector3D<Cartesian3D<double> > pixelPositionGlobal = *(get_detector(dutID)->m_localToGlobal) * pixelPositionLocal;
+          PositionVector3D<Cartesian3D<double> > pixelPositionLocal =
+      get_detector(dutID)->getLocalPosition((*itPixel)->row(),(*itPixel)->column()); PositionVector3D<Cartesian3D<double> >
+      pixelPositionGlobal = *(get_detector(dutID)->m_localToGlobal) * pixelPositionLocal;
 
           // Check if it is close to the track
           if( fabs( pixelPositionGlobal.X() - trackIntercept.X() ) > m_associationCut ||

src/modules/Millepede/Millepede.h

@@ -12,7 +12,7 @@ namespace corryvreckan {
      *
      *  @author Christoph Hombach
      *  @date   2012-06-19
-    */
+     */
     class Millepede : public Module {
     public:
         /// Constructor

src/objects/Cluster.h

@@ -78,6 +78,19 @@ namespace corryvreckan {
         double rowWidth() { return m_rowWidth; }
         Pixels* pixels() { return (&m_pixels); }
 
+        // Retrieve the seed pixel of the cluster, defined as the one with the highest charge:
+        Pixel* getSeedPixel() {
+            Pixel* seed;
+            double maxcharge = -1;
+            for(auto& px : m_pixels) {
+                if(px->charge() > maxcharge) {
+                    maxcharge = px->charge();
+                    seed = px;
+                }
+            }
+            return seed;
+        }
+
         // Set cluster parameters
         void setRow(double row) { m_row = row; }
         void setColumn(double col) { m_column = col; }
@@ -116,7 +129,7 @@ namespace corryvreckan {
         std::map<int, bool> m_columnHits;
 
         // ROOT I/O class definition - update version number when you change this class!
-        ClassDef(Cluster, 5)
+        ClassDef(Cluster, 6)
     };
 
     // Vector type declaration

src/objects/Pixel.h

@@ -13,7 +13,7 @@ namespace corryvreckan {
         virtual ~Pixel() {}
         Pixel(std::string detectorID, int row, int col, int tot) : Pixel(detectorID, row, col, tot, 0.) {}
         Pixel(std::string detectorID, int row, int col, int tot, double timestamp)
-            : Object(detectorID, timestamp), m_row(row), m_column(col), m_adc(tot) {}
+            : Object(detectorID, timestamp), m_row(row), m_column(col), m_adc(tot), m_charge(tot) {}
 
         int row() const { return m_row; }
         int column() const { return m_column; }
@@ -40,6 +40,7 @@ namespace corryvreckan {
         int m_row;
         int m_column;