AnalysisEfficiency: merged .cpp from jens_dev (next: update README). No old functionality is lost

src/modules/AnalysisEfficiency/AnalysisEfficiency.cpp

@@ -30,6 +30,8 @@ void AnalysisEfficiency::initialise() {
     auto pitch_x = static_cast<double>(Units::convert(m_detector->pitch().X(), "um"));
     auto pitch_y = static_cast<double>(Units::convert(m_detector->pitch().Y(), "um"));
 
+    hPixelEfficiency = new TH1D("hPixelEfficiency", "hPixelEfficiency", 220, 0, 1.1);
+
     std::string title = m_detector->name() + " Pixel efficiency map;x_{track} mod " + std::to_string(pitch_x) +
                         "#mum;y_{track} mod " + std::to_string(pitch_y) + "#mum;efficiency";
     hPixelEfficiencyMap_trackPos = new TProfile2D("pixelEfficiencyMap_trackPos",
@@ -86,6 +88,78 @@ void AnalysisEfficiency::initialise() {
                                                    1.5 * m_detector->size().Y(),
                                                    0,
                                                    1);
+    eTotalEfficiency = new TEfficiency("eTotalEfficiency", "totalEfficiency", 1, 0, 1);
+    totalEfficiency = new TNamed("totalEffiency", "totalEffiency");
+
+    hTrackTimeToPrevHit_matched =
+        new TH1D("trackTimeToPrevHit_matched", "trackTimeToPrevHit_matched;time to prev hit [us];# events", 1e6, 0, 1e6);
+    hTrackTimeToPrevHit_notmatched = new TH1D(
+        "trackTimeToPrevHit_notmatched", "trackTimeToPrevHit_notmatched;time to prev hit [us];# events", 1e6, 0, 1e6);
+
+    title = m_detector->name() + "time difference to previous track (if this has assoc cluster)";
+    hTimeDiffPrevTrack_assocCluster = new TH1D("timeDiffPrevTrack_assocCluster", title.c_str(), 11000, -1000, 10000);
+    hTimeDiffPrevTrack_assocCluster->GetXaxis()->SetTitle("time diff [#mus]");
+    hTimeDiffPrevTrack_assocCluster->GetYaxis()->SetTitle("events");
+    title = m_detector->name() + "time difference to previous track (if this has no assoc cluster)";
+    hTimeDiffPrevTrack_noAssocCluster = new TH1D("timeDiffPrevTrack_noAssocCluster", title.c_str(), 11000, -1000, 10000);
+    hTimeDiffPrevTrack_noAssocCluster->GetXaxis()->SetTitle("time diff [#mus]");
+    hTimeDiffPrevTrack_noAssocCluster->GetYaxis()->SetTitle("events");
+
+    hRowDiffPrevTrack_assocCluster = new TH1D("rowDiffPrevTrack_assocCluster",
+                                              "rowDiffPrevTrack_assocCluster",
+                                              2 * m_detector->nPixels().Y(),
+                                              -m_detector->nPixels().Y(),
+                                              m_detector->nPixels().Y());
+    hRowDiffPrevTrack_assocCluster->GetXaxis()->SetTitle("row difference");
+    hRowDiffPrevTrack_assocCluster->GetYaxis()->SetTitle("events");
+    hColDiffPrevTrack_assocCluster = new TH1D("colDiffPrevTrack_assocCluster",
+                                              "colDiffPrevTrack_assocCluster",
+                                              2 * m_detector->nPixels().X(),
+                                              -m_detector->nPixels().X(),
+                                              m_detector->nPixels().X());
+    hColDiffPrevTrack_assocCluster->GetXaxis()->SetTitle("column difference");
+    hColDiffPrevTrack_assocCluster->GetYaxis()->SetTitle("events");
+
+    hRowDiffPrevTrack_noAssocCluster = new TH1D("rowDiffPrevTrack_noAssocCluster",
+                                                "rowDiffPrevTrack_assocCluster",
+                                                2 * m_detector->nPixels().Y(),
+                                                -m_detector->nPixels().Y(),
+                                                m_detector->nPixels().Y());
+    hRowDiffPrevTrack_noAssocCluster->GetXaxis()->SetTitle("row difference");
+    hRowDiffPrevTrack_noAssocCluster->GetYaxis()->SetTitle("events");
+    hColDiffPrevTrack_noAssocCluster = new TH1D("colDiffPrevTrack_noAassocCluster",
+                                                "colDiffPrevTrack_assocCluster",
+                                                2 * m_detector->nPixels().X(),
+                                                -m_detector->nPixels().X(),
+                                                m_detector->nPixels().X());
+    hColDiffPrevTrack_noAssocCluster->GetXaxis()->SetTitle("column difference");
+    hColDiffPrevTrack_noAssocCluster->GetYaxis()->SetTitle("events");
+
+    hPosDiffPrevTrack_assocCluster = new TH2D("posDiffPrevTrack_assocCluster",
+                                              "posDiffPrevTrack_assocCluster",
+                                              2 * m_detector->nPixels().X(),
+                                              -m_detector->nPixels().X(),
+                                              m_detector->nPixels().X(),
+                                              2 * m_detector->nPixels().Y(),
+                                              -m_detector->nPixels().Y(),
+                                              m_detector->nPixels().Y());
+    hPosDiffPrevTrack_noAssocCluster = new TH2D("posDiffPrevTrack_noAssocCluster",
+                                                "posDiffPrevTrack_noAssocCluster",
+                                                2 * m_detector->nPixels().X(),
+                                                -m_detector->nPixels().X(),
+                                                m_detector->nPixels().X(),
+                                                2 * m_detector->nPixels().Y(),
+                                                -m_detector->nPixels().Y(),
+                                                m_detector->nPixels().Y());
+
+    // initialize matrix with hit timestamps to all 0:
+    std::vector<double> v_row;
+    for(int irow = 0; irow < m_detector->nPixels().Y(); irow++) { // create vector will zeros of length <nRows>
+        v_row.push_back(0);
+    }
+    for(int icol = 0; icol < m_detector->nPixels().X(); icol++) { // push back to prev_hits_ts <nCols> times
+        prev_hit_ts.push_back(v_row);
+    }
 }
 
 StatusCode AnalysisEfficiency::run(std::shared_ptr<Clipboard> clipboard) {
@@ -122,6 +196,7 @@ StatusCode AnalysisEfficiency::run(std::shared_ptr<Clipboard> clipboard) {
         if(!m_detector->isWithinROI(track)) {
             LOG(DEBUG) << " - track outside ROI";
             is_within_roi = false;
+            // here we don't continue because only some particular histograms shall be effected
         }
 
         // Check that it doesn't go through/near a masked pixel
@@ -186,21 +261,79 @@ StatusCode AnalysisEfficiency::run(std::shared_ptr<Clipboard> clipboard) {
         // For pixels, only look at the ROI:
         if(is_within_roi) {
             hPixelEfficiencyMap_trackPos->Fill(xmod, ymod, has_associated_cluster);
+            eTotalEfficiency->Fill(has_associated_cluster, 0); // use 0th bin for total efficiency
         }
-        if(has_associated_cluster == false) {
+
+        size_t intercept_col = static_cast<size_t>(m_detector->getColumn(localIntercept));
+        size_t intercept_row = static_cast<size_t>(m_detector->getRow(localIntercept));
+
+        if(has_associated_cluster) {
+            hTimeDiffPrevTrack_assocCluster->Fill(
+                static_cast<double>(Units::convert(track->timestamp() - last_track_timestamp, "us")));
+            hRowDiffPrevTrack_assocCluster->Fill(m_detector->getRow(localIntercept) - last_track_row);
+            hColDiffPrevTrack_assocCluster->Fill(m_detector->getColumn(localIntercept) - last_track_col);
+            hPosDiffPrevTrack_assocCluster->Fill(m_detector->getColumn(localIntercept) - last_track_col,
+                                                 m_detector->getRow(localIntercept) - last_track_row);
+            if((prev_hit_ts.at(intercept_col)).at(intercept_row) != 0) {
+                hTrackTimeToPrevHit_matched->Fill(static_cast<double>(
+                    Units::convert(track->timestamp() - prev_hit_ts.at(intercept_col).at(intercept_row), "us")));
+            }
+        } else {
             hGlobalEfficiencyMap_clustPos->Fill(globalIntercept.X(), globalIntercept.Y(), has_associated_cluster);
             hChipEfficiencyMap_clustPos->Fill(
                 m_detector->getColumn(localIntercept), m_detector->getRow(localIntercept), has_associated_cluster);
+
+            hTimeDiffPrevTrack_noAssocCluster->Fill(
+                static_cast<double>(Units::convert(track->timestamp() - last_track_timestamp, "us")));
+            hRowDiffPrevTrack_noAssocCluster->Fill(m_detector->getRow(localIntercept) - last_track_row);
+            hColDiffPrevTrack_noAssocCluster->Fill(m_detector->getColumn(localIntercept) - last_track_col);
+            hPosDiffPrevTrack_noAssocCluster->Fill(m_detector->getColumn(localIntercept) - last_track_col,
+                                                   m_detector->getRow(localIntercept) - last_track_row);
+            if((prev_hit_ts.at(intercept_col)).at(intercept_row) != 0) {
+                LOG(DEBUG) << "Found a time difference of "
+                           << static_cast<double>(
+                                  Units::convert(track->timestamp() - prev_hit_ts.at(intercept_col).at(intercept_row), "us"))
+                           << " us";
+                hTrackTimeToPrevHit_notmatched->Fill(static_cast<double>(
+                    Units::convert(track->timestamp() - prev_hit_ts.at(intercept_col).at(intercept_row), "us")));
+            }
         }
+        last_track_timestamp = track->timestamp();
+        last_track_col = m_detector->getColumn(localIntercept);
+        last_track_row = m_detector->getRow(localIntercept);
+    } // end loop over tracks
+
+    // Before going to the next event, loop over all pixels (all hits incl. noise)
+    // and fill matrix with timestamps of previous pixels.
+    Pixels* pixels = reinterpret_cast<Pixels*>(clipboard->get(m_detector->name(), "pixels"));
+    if(pixels == nullptr) {
+        LOG(DEBUG) << "Detector " << m_detector->name() << " does not have any pixels on the clipboard";
+        return StatusCode::Success;
+    }
+    for(auto& pixel : (*pixels)) {
+        prev_hit_ts.at(static_cast<size_t>(pixel->column())).at(static_cast<size_t>(pixel->row())) = pixel->timestamp();
     }
 
     return StatusCode::Success;
 }
 
 void AnalysisEfficiency::finalise() {
-    LOG(INFO) << "No. matched tracks=" << matched_tracks;
-    LOG(INFO) << "Total no. tracks=" << total_tracks;
-    LOG(STATUS) << "Total efficiency of detector " << m_detector->name() << ": "
-                << (100 * matched_tracks / (total_tracks > 0 ? total_tracks : 1)) << "%, measured with " << total_tracks
-                << " tracks";
+    double totalEff =
+        100 * static_cast<double>(matched_tracks) / (total_tracks > 0 ? static_cast<double>(total_tracks) : 1.);
+    LOG(STATUS) << "Total efficiency of detector " << m_detector->name() << ": " << totalEff << "%, measured with "
+                << matched_tracks << "/" << total_tracks << " matched/total tracks";
+
+    totalEfficiency->SetName((to_string(totalEff) + " %").c_str());
+    totalEfficiency->Write();
+
+    for(int icol = 1; icol < m_detector->nPixels().X() + 1; icol++) {
+        for(int irow = 1; irow < m_detector->nPixels().Y() + 1; irow++) {
+            // calculate total efficiency: (just to double check the other calculation)
+            double eff = hChipEfficiencyMap_trackPos->GetBinContent(icol, irow);
+            if(eff > 0) {
+                LOG(DEBUG) << "col/row = " << icol << "/" << irow << ", binContent = " << eff;
+                hPixelEfficiency->Fill(hChipEfficiencyMap_trackPos->GetBinContent(icol, irow));
+            }
+        }
+    }
 }