Merge branch 'f/fakeRate' into 'master'

AnalysisEfficiency: Added fake rate analysis and plots.

See merge request corryvreckan/corryvreckan!610

src/modules/AnalysisEfficiency/AnalysisEfficiency.cpp

@@ -27,6 +27,10 @@ AnalysisEfficiency::AnalysisEfficiency(Configuration& config, std::shared_ptr<De
     config_.setDefault<XYVector>("inpixel_cut_edge", {Units::get(5.0, "um"), Units::get(5.0, "um")});
     config_.setDefault<double>("masked_pixel_distance_cut", 1.);
     config_.setDefault<double>("spatial_cut_sensoredge", 1.);
+    config_.setDefault<FakeRateMethod>("fake_rate_method", FakeRateMethod::RADIUS);
+    config_.setDefault<double>("fake_rate_distance", 2.);
+    config_.setDefault<int>("n_charge_bins", 1000);
+    config_.setDefault<double>("charge_histo_range", 1000.0);
 
     m_timeCutFrameEdge = config_.get<double>("time_cut_frameedge");
     m_chi2ndofCut = config_.get<double>("chi2ndof_cut");
@@ -35,6 +39,10 @@ AnalysisEfficiency::AnalysisEfficiency(Configuration& config, std::shared_ptr<De
     m_inpixelEdgeCut = config_.get<XYVector>("inpixel_cut_edge");
     m_maskedPixelDistanceCut = config_.get<int>("masked_pixel_distance_cut");
     spatial_cut_sensoredge = config_.get<double>("spatial_cut_sensoredge");
+    m_fake_rate_method = config_.get<FakeRateMethod>("fake_rate_method");
+    m_fake_rate_distance = config_.get<double>("fake_rate_distance");
+    m_n_charge_bins = config_.get<int>("n_charge_bins");
+    m_charge_histo_range = config_.get<double>("charge_histo_range");
 }
 void AnalysisEfficiency::initialize() {
 
@@ -294,6 +302,48 @@ void AnalysisEfficiency::initialize() {
     auto nCols = static_cast<size_t>(m_detector->nPixels().X());
     std::vector<double> v_row(nRows, 0.); // create vector will zeros of length <nRows>
     prev_hit_ts.assign(nCols, v_row);     // use vector v_row to construct matrix
+
+    createFakeRatePlots();
+}
+
+void AnalysisEfficiency::createFakeRatePlots() {
+    TDirectory* directory = getROOTDirectory();
+    TDirectory* fake_rate_directory = directory->mkdir("fake_rate");
+    if(fake_rate_directory == nullptr) {
+        throw RuntimeError("Cannot create or access fake rate ROOT directory for module " + this->getUniqueName());
+    }
+    fake_rate_directory->cd();
+
+    std::string title = m_detector->getName() + " number of fake hits per event; hits; events";
+    hFakePixelPerEvent = new TH1D("hFakePixelPerEvent", title.c_str(), 25, 0 - 0.5, 25 - 0.5);
+
+    title = m_detector->getName() + " pixel fake hits per event;x [px];y [px]; hits";
+    fakePixelPerEventMap = new TH2D("fakePixelPerEventMap",
+                                    title.c_str(),
+                                    m_detector->nPixels().X(),
+                                    -0.5,
+                                    m_detector->nPixels().X() - 0.5,
+                                    m_detector->nPixels().Y(),
+                                    -0.5,
+                                    m_detector->nPixels().Y() - 0.5);
+
+    title = m_detector->getName() + " pixel fake hits per event vs. time; time [s]; hits";
+    fakePixelPerEventVsTime = new TProfile("fakePixelPerEventVsTime", title.c_str(), 3000, 0, 3000);
+
+    title = m_detector->getName() + " pixel fake hits per event vs. time; time [s]; hits";
+    fakePixelPerEventVsTimeLong = new TProfile("efficiencyVsTimeLong", title.c_str(), 3000, 0, 30000);
+
+    title = m_detector->getName() + " charge distribution for fake pixels; charge [a.u.]; entries";
+    hFakePixelCharge = new TH1D("hFakePixelCharge", title.c_str(), m_n_charge_bins, 0.0, m_charge_histo_range);
+
+    title = m_detector->getName() + " charge distribution for fake clusters; charge [a.u.]; entries";
+    hFakeClusterCharge = new TH1D("hFakeClusterCharge", title.c_str(), m_n_charge_bins, 0.0, m_charge_histo_range);
+
+    title = m_detector->getName() + " number of fake clusters per event; clusters; events";
+    hFakeClusterPerEvent = new TH1D("hFakeClusterPerEvent", title.c_str(), 25, 0 - 0.5, 25 - 0.5);
+
+    title = m_detector->getName() + " cluster size of fake clusters; cluster size; events";
+    hFakeClusterSize = new TH1D("hFakeClusterSize", title.c_str(), 25, 0 - 0.5, 25 - 0.5);
 }
 
 StatusCode AnalysisEfficiency::run(const std::shared_ptr<Clipboard>& clipboard) {
@@ -516,6 +566,104 @@ StatusCode AnalysisEfficiency::run(const std::shared_ptr<Clipboard>& clipboard)
         }
     }
 
+    // fake rate analysis
+    if(m_fake_rate_method == FakeRateMethod::RADIUS) {
+        LOG_ONCE(STATUS) << "Estimating fake rate based on radial cut around pixels (RADIUS method).";
+
+        int fake_hits = 0;
+        int fake_clusters = 0;
+
+        // get and iterate dut clusters from clipboard
+        auto clusters = clipboard->getData<Cluster>(m_detector->getName());
+        for(auto& cluster : clusters) {
+
+            bool track_too_close = false;
+            // iterate the tracks from the clipboard
+            for(auto& track : tracks) {
+
+                // discard tracks without intercept, using a tolerance defined
+                // by the radial cut which we will use later.
+                if(m_detector->hasIntercept(track.get(), -m_fake_rate_distance)) {
+                    continue;
+                }
+
+                // now check if the track is too close to the considered cluster,
+                // using the distance in units of the pitch.
+                auto interceptLocal = m_detector->getLocalIntercept(track.get());
+                double norm_xdistance = (interceptLocal.X() - cluster->local().x()) / m_detector->getPitch().X();
+                double norm_ydistance = (interceptLocal.Y() - cluster->local().y()) / m_detector->getPitch().Y();
+                double norm_rdistance = sqrt(norm_xdistance * norm_xdistance + norm_ydistance * norm_ydistance);
+                if(norm_rdistance < m_fake_rate_distance) {
+                    track_too_close = true;
+                    break;
+                }
+            }
+
+            // now study cluster and pixel properties, which seems to be a fake
+            // cluster, as there is no track nearby.
+            if(!track_too_close) {
+
+                fake_clusters++;
+                hFakeClusterCharge->Fill(cluster->charge());
+                hFakeClusterSize->Fill(static_cast<double>(cluster->size()));
+
+                for(auto& pixel : cluster->pixels()) {
+                    fake_hits++;
+                    hFakePixelCharge->Fill(pixel->charge());
+                    fakePixelPerEventMap->Fill(pixel->column(), pixel->row(), 1);
+                }
+            }
+        }
+
+        hFakePixelPerEvent->Fill(fake_hits);
+        fakePixelPerEventVsTime->Fill(static_cast<double>(Units::convert(event->start(), "s")), fake_hits);
+        fakePixelPerEventVsTimeLong->Fill(static_cast<double>(Units::convert(event->start(), "s")), fake_hits);
+        hFakeClusterPerEvent->Fill(fake_clusters);
+    }
+    if(m_fake_rate_method == FakeRateMethod::EDGE) {
+        LOG_ONCE(STATUS) << "Estimating fake rate based on events without DUT intercepting tracks (EDGE method).";
+
+        bool track_in_active = false;
+        int fake_hits = 0;
+        int fake_clusters = 0;
+
+        // iterate the tracks from the clipboard
+        for(auto& track : tracks) {
+            // check if one of the tracks intercepts the dut
+            // A positive value of m_fake_rate_radius means that we want to
+            // check an area that is larger than the sensor matrix.
+            if(m_detector->hasIntercept(track.get(), -m_fake_rate_distance)) {
+                track_in_active = true;
+                break;
+            }
+        }
+
+        // Study the dut response if there is no track pointing to the active
+        // area of the dut. There might still be particles, that we failed to
+        // track!
+        if(!track_in_active) {
+
+            // iterate the dut pixels from clipboard
+            for(auto& pixel : pixels) {
+                fake_hits++;
+                hFakePixelCharge->Fill(pixel->charge());
+                fakePixelPerEventMap->Fill(pixel->column(), pixel->row(), 1);
+            }
+            hFakePixelPerEvent->Fill(fake_hits);
+            fakePixelPerEventVsTime->Fill(static_cast<double>(Units::convert(event->start(), "s")), fake_hits);
+            fakePixelPerEventVsTimeLong->Fill(static_cast<double>(Units::convert(event->start(), "s")), fake_hits);
+
+            // get and iterate dut clusters from clipboard
+            auto clusters = clipboard->getData<Cluster>(m_detector->getName());
+            for(auto& cluster : clusters) {
+                fake_clusters++;
+                hFakeClusterCharge->Fill(cluster->charge());
+                hFakeClusterSize->Fill(static_cast<double>(cluster->size()));
+            }
+            hFakeClusterPerEvent->Fill(fake_clusters);
+        }
+    }
+
     return StatusCode::Success;
 }
 
@@ -553,4 +701,8 @@ void AnalysisEfficiency::finalize(const std::shared_ptr<ReadonlyClipboard>&) {
             }
         }
     }
+    // normalize fake rate map, if it exists
+    if(hFakeClusterPerEvent->GetEntries() > 0) {
+        fakePixelPerEventMap->Scale(1. / hFakePixelPerEvent->GetEntries());
+    }
 }

src/modules/AnalysisEfficiency/AnalysisEfficiency.h

@@ -17,6 +17,7 @@
 
 #include "core/module/Module.hpp"
 
+#include <TDirectory.h>
 #include "TEfficiency.h"
 #include "TH2D.h"
 #include "TNamed.h"
@@ -84,7 +85,25 @@ namespace corryvreckan {
         TH2D* hPosDiffPrevTrack_noAssocCluster;
         TH2D* hDistanceCluster_track;
         TH2D* htimeRes_cluster_size;
-        double m_chi2ndofCut, m_timeCutFrameEdge, m_inpixelBinSize, spatial_cut_sensoredge;
+
+        // fake rate plots
+        TH1D* hFakePixelPerEvent;
+        TH2D* fakePixelPerEventMap;
+        TProfile* fakePixelPerEventVsTime;
+        TProfile* fakePixelPerEventVsTimeLong;
+        TH1D* hFakePixelCharge;
+        TH1D* hFakeClusterPerEvent;
+        TH1D* hFakeClusterCharge;
+        TH1D* hFakeClusterSize;
+
+        enum class FakeRateMethod {
+            RADIUS,
+            EDGE,
+        } m_fake_rate_method;
+
+        double m_chi2ndofCut, m_timeCutFrameEdge, m_inpixelBinSize, spatial_cut_sensoredge, m_fake_rate_distance,
+            m_charge_histo_range;
+        int m_n_charge_bins;
         XYVector m_inpixelEdgeCut;
         int m_maskedPixelDistanceCut = 1;
         int total_tracks = 0;
@@ -97,6 +116,8 @@ namespace corryvreckan {
         std::vector<std::string> require_associated_cluster_on_;
 
         Matrix<double> prev_hit_ts; // matrix containing previous hit timestamp for every pixel
+
+        void createFakeRatePlots();
     };
 
 } // namespace corryvreckan

src/modules/AnalysisEfficiency/README.md

@@ -3,9 +3,9 @@
 # SPDX-License-Identifier: CC-BY-4.0 OR MIT
 ---
 # AnalysisEfficiency
-**Maintainer**: Simon Spannagel (<simon.spannagel@cern.ch>), Jens Kroeger (<jens.kroeger@cern.ch>)  
-**Module Type**: *DUT*  
-**Detector Type**: *all*  
+**Maintainer**: Simon Spannagel (<simon.spannagel@cern.ch>), Jens Kroeger (<jens.kroeger@cern.ch>)
+**Module Type**: *DUT*
+**Detector Type**: *all*
 **Status**: Functional
 
 ### Description
@@ -27,6 +27,10 @@ More information can be found in the ROOT `TEfficiency` class reference, section
 * `masked_pixel_distance_cut`: Distance (in pixels) to exclude tracks passing close to masked pixel. Defaults to `1`.
 * `require_associated_cluster_on`: Names of detectors which are required to have an associated cluster to the telescope tracks. Detectors listed here must be marked as `role = DUT` in the detector configuration file. Only tracks satisfying this requirement are accepted for the efficiency measurement. If empty, no detector is required. Default is empty.
 * `spatial_cut_sensoredge`: Parameter to discard telescope tracks at the sensor edges in fractions of pixel pitch. Defaults to `1`.
+* `fake_rate_distance`: Distance cut used in the fake rate estimate. Given in units of the pixel pitch. Defaults to `2`.
+* `fake_rate_method`: Method used in the fake rate estimate. The main idea is to look for clusters and hits that are far away from reconstructed tracks. Those are defined as fake, which neglects effects of tracking in-efficiency. There are two slightly different methods. The `RADIUS` method considers DUT clusters without reconstructed track in a radius of `fake_rate_distance` pixel pitches around them as fake. The `EDGE` method looks for events without tracks intercepting the DUT within its active area plus `fake_rate_distance` pixel pitches in each direction, and considers all DUT activity in these events as fake. The former method is intended for large DUTs, the latter for small ones. Defaults to `RADIUS`.
+* `n_charge_bins`: Number of bins for pixel and cluster charge distributions. Defaults to `1000`.
+* `charge_histo_range`: Maximum value for pixel and cluster charge distributions. Defaults to `1000`.
 
 ### Plots produced
 
@@ -46,6 +50,12 @@ For the DUT, the following plots are produced:
   * Histograms of the time difference of a matched (non-matched) cluster to a previous hit (not matter if noise or track)
   * Distribution of cluster-track distances
 
+* Fake rate plots (if enabled):
+  * Number of fake pixels per event as histogram, map, and as function of time.
+  * Number of fake clusters per event as histogram.
+  * Pixel and cluster charge distributions for fake pixels and clusters.
+  * Cluster size of fake clusters as histogram.
+
 * Other:
   * Value of total efficiency as `TEfficiency` including (asymmetric) error bars (total and restricted to in-pixel ROI)
   * Efficiency as function of column and row, and vs. time