From d9e494ec81e8627307ae398c74302565e32736f3 Mon Sep 17 00:00:00 2001
From: "s.zharko@gsi.de" <s.zharko@gsi.de>
Date: Wed, 26 Oct 2022 20:13:20 +0200
Subject: [PATCH] L1: added writer for reconstructable MC tracks
---
macro/run/run_reco.C | 2 +
reco/L1/CbmL1.cxx | 13 +++++-
reco/L1/CbmL1.h | 38 ++++++++++++++---
reco/L1/CbmL1Performance.cxx | 83 +++++++++++++++++++++++++++++++++++-
4 files changed, 128 insertions(+), 8 deletions(-)
diff --git a/macro/run/run_reco.C b/macro/run/run_reco.C
index bbb6675c92..5c3a90adc3 100644
--- a/macro/run/run_reco.C
+++ b/macro/run/run_reco.C
@@ -381,7 +381,9 @@ void run_reco(TString input = "", Int_t nTimeSlices = -1, Int_t firstTimeSlice =
// L1 tracking
auto l1 = (debugWithMC) ? new CbmL1("L1", 2, 3, 0) : new CbmL1("L1", 0);
+ // L1 configuration file (optional)
//l1->SetInputConfigName(TString(gSystem->Getenv("VMCWORKDIR")) + "/reco/L1/L1Algo/L1ConfigExample.yaml");
+ if (debugWithMC) { l1->SetOutputMcTracksNtupleFilename("mcTracksNtuple"); }
// --- Material budget file names
TString mvdGeoTag;
diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index 96179436bb..ca0d359f2e 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -49,6 +49,7 @@
#include "TGeoManager.h"
#include "TGeoNode.h"
#include "TMatrixD.h"
+#include "TNtuple.h"
#include "TProfile2D.h"
#include "TROOT.h"
#include "TRandom3.h"
@@ -1180,6 +1181,7 @@ void CbmL1::Reconstruct(CbmEvent* event)
EfficienciesPerformance();
HistoPerformance();
TrackFitPerformance();
+ if (fsMcTracksOutputFilename.size()) { DumpMCTracksToNtuple(); }
// TimeHist();
/// WriteSIMDKFData();
}
@@ -1228,6 +1230,13 @@ void CbmL1::Finish()
outfile->Delete();
}
+ if (fpMcTracksOutFile) {
+ fpMcTracksOutFile->cd();
+ fpMcTracksTree->Write();
+ fpMcTracksOutFile->Close();
+ fpMcTracksOutFile->Delete();
+ }
+
gFile = currentFile;
gDirectory = curr;
fpAlgo->Finish();
@@ -1500,6 +1509,7 @@ void CbmL1::WriteSIMDKFData()
}
}
+// ---------------------------------------------------------------------------------------------------------------------
// NOTE: this function should be called before fInitManager.SendParameters(fpAlgo)
std::vector<L1Material> CbmL1::ReadMaterialBudget(L1DetectorID detectorID)
{
@@ -1544,7 +1554,8 @@ std::vector<L1Material> CbmL1::ReadMaterialBudget(L1DetectorID detectorID)
return result;
}
-
+// ---------------------------------------------------------------------------------------------------------------------
+//
double CbmL1::boundedGaus(double sigma)
{
assert(sigma > 0. && std::isfinite(sigma));
diff --git a/reco/L1/CbmL1.h b/reco/L1/CbmL1.h
index 1bccc2849e..f3b5edc187 100644
--- a/reco/L1/CbmL1.h
+++ b/reco/L1/CbmL1.h
@@ -69,6 +69,7 @@ class CbmMCDataObject;
class CbmEvent;
class TProfile2D;
+class TNtuple;
/// TODO: SZh 21.09.2022: Replace instances of this class with L1Hit
class CbmL1HitStore {
@@ -266,6 +267,11 @@ public:
/// \param filename Name of the input tracking configuration file
void SetOutputConfigName(const char* filename) { fInitManager.SetOutputConfigName(std::string(filename)); }
+ /// \brief Sets output file for MC tracks ntuple
+ /// If the filename is empty string, ntuple is not filled
+ /// \param filename Name of the output file name
+ void SetOutputMcTracksNtupleFilename(const char* filename) { fsMcTracksOutputFilename = std::string(filename); }
+
/// Sets flag: to correct input hits on MC or not
/// \param flag: true - hits will be corrected on MC information
void SetCorrectHitsOnMC(bool flag) { fIfCorrectHitsOnMC = flag; }
@@ -399,11 +405,29 @@ private:
void InputPerformance(); // Build histograms about input data, like hit pulls, etc.
void TimeHist();
- /// Reconstruction Performance
- void TrackMatch(); // Procedure for match Reconstructed and MC Tracks. Should be called before Performances
- void EfficienciesPerformance(); // calculate efficiencies
- void TrackFitPerformance(); // pulls & residuals. Can be called only after Performance()
- void HistoPerformance(); // fill some histograms and calculate efficiencies
+ // ********************************
+ // ** Reconstruction Performance **
+ // ********************************
+
+
+ /// Matches reconstructed and MC tracks
+ /// \note Should be called before Performances
+ void TrackMatch();
+
+ /// Calculates tracking efficiencies (counters)
+ void EfficienciesPerformance();
+
+ /// Builds pulls and residuals
+ /// \note Should be called only after CbmL1::Performance()
+ void TrackFitPerformance();
+
+ /// Fills performance histograms
+ void HistoPerformance();
+
+ /// Writes MC tracks to ntuple
+ /// \note Executed only if the filename for MC tracks ntuple output is defined
+ void DumpMCTracksToNtuple();
+
// ** STandAlone Package service-functions **
@@ -636,13 +660,15 @@ private:
static const int fNGhostHistos = 9;
TH1F* fGhostHisto[fNGhostHistos] {nullptr};
+ TFile* fpMcTracksOutFile = nullptr; ///< File to save MC-tracks ntuple
+ TNtuple* fpMcTracksTree = nullptr; ///< Ntuple to save MC-tracks
+ std::string fsMcTracksOutputFilename = ""; ///< Name of file to save MC-tracks ntuple
int fFindParticlesMode {0}; // 0 - don't run FindParticles
// 1 - run, all MC particle is reco-able
// 2 - run, MC particle is reco-able if created from reco-able tracks
// 3 - run, MC particle is reco-able if created from reconstructed tracks
-
std::unordered_map<L1DetectorID, TString>
fMatBudgetFileName {}; ///< Map for material budget file names vs. detectorID
diff --git a/reco/L1/CbmL1Performance.cxx b/reco/L1/CbmL1Performance.cxx
index ed597abe04..3bfd392054 100644
--- a/reco/L1/CbmL1Performance.cxx
+++ b/reco/L1/CbmL1Performance.cxx
@@ -37,13 +37,17 @@
#include "CbmTrdHit.h"
#include "CbmTrdPoint.h"
+#include "FairRunAna.h"
#include "FairTrackParam.h" // for vertex pulls
+#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "TMath.h"
+#include "TNtuple.h"
#include "TProfile.h"
-#include <TFile.h>
+
+#include <boost/filesystem.hpp>
#include <iostream>
#include <list>
@@ -2244,3 +2248,80 @@ void CbmL1::InputPerformance()
} // void CbmL1::InputPerformance()
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void CbmL1::DumpMCTracksToNtuple()
+{
+ if (!fpMcTracksTree) {
+ auto* currentDir = gDirectory;
+ auto* currentFile = gFile;
+
+
+ // Get prefix and directory
+ boost::filesystem::path p = (FairRunAna::Instance()->GetUserOutputFileName()).Data();
+ std::string dir = p.parent_path().string();
+ if (dir.empty()) dir = ".";
+ std::string filename = dir + "/" + fsMcTracksOutputFilename + "." + p.filename().string();
+ std::cout << "\033[1;31mSAVING A TREE: " << filename << "\033[0m\n";
+
+ fpMcTracksOutFile = new TFile(filename.c_str(), "RECREATE");
+ fpMcTracksOutFile->cd();
+ fpMcTracksTree = new TNtuple("t", "McTracks", "motherId:pdg:p:q:zv:s:x0:y0:z0:x1:y1:z1:x2:y2:z2");
+ // motherID - id of mother particle (< 0, if primary)
+ // pdg - PDG code of particle
+ // p - absolute value of track momentum [GeV/c]
+ // q - charge of track [e]
+ // zv - z-component of track vertex [cm]
+ // s - global index of station
+ // x0, y0, z0 - position of the left MC point in a triplet [cm]
+ // x1, y1, z1 - position of the middle MC point in a triplet [cm]
+ // x2, y2, z2 - position of the right MC point in a triplet [cm]
+
+ gFile = currentFile;
+ gDirectory = currentDir;
+ }
+
+ struct ReducedMcPoint {
+ int s; ///< global active index of tracking station
+ float x; ///< x-component of point position [cm]
+ float y; ///< y-component of point position [cm]
+ float z; ///< z-component of point position [cm]
+ bool operator<(const ReducedMcPoint& other) { return s < other.s; }
+ };
+
+ for (const auto& tr : fvMCTracks) {
+ // Use only reconstructable tracks
+ if (!tr.IsReconstructable()) { continue; }
+
+ std::vector<ReducedMcPoint> vPoints;
+ vPoints.reserve(tr.Points.size());
+ for (unsigned int iP = 0; iP < tr.Points.size(); ++iP) {
+ const auto& point = fvMCPoints[tr.Points[iP]];
+ vPoints.emplace_back(ReducedMcPoint {point.iStation, float(point.x), float(point.y), float(point.z)});
+ }
+
+ std::sort(vPoints.begin(), vPoints.end());
+ for (unsigned int i = 0; i + 2 < vPoints.size(); ++i) {
+ // Condition to collect only triplets without gaps in stations
+ // TODO: SZh 20.10.2022 Add cases for jump iterations
+ if (vPoints[i + 1].s == vPoints[i].s + 1 && vPoints[i + 2].s == vPoints[i].s + 2) {
+ fpMcTracksTree->Fill(tr.mother_ID, ///< index of mother particle
+ tr.pdg, ///< PDG code
+ tr.p, ///< absolute value of track momentum [GeV/c]
+ tr.q, ///< charge of track [e]
+ tr.z, ///< z-position of track vertex [cm]
+ vPoints[i].s, ///< global index of active tracking station
+ vPoints[i].x, ///< x-component of the left MC point in a triplet [cm]
+ vPoints[i].y, ///< y-component of the left MC point in a triplet [cm]
+ vPoints[i].z, ///< z-component of the left MC point in a triplet [cm]
+ vPoints[i + 1].x, ///< x-component of the middle MC point in a triplet [cm]
+ vPoints[i + 1].y, ///< y-component of the middle MC point in a triplet [cm]
+ vPoints[i + 1].z, ///< z-component of the middle MC point in a triplet [cm]
+ vPoints[i + 2].x, ///< x-component of the right MC point in a triplet [cm]
+ vPoints[i + 2].y, ///< y-component of the right MC point in a triplet [cm]
+ vPoints[i + 2].z); ///< z-component of the right MC point in a triplet [cm]
+ }
+ }
+ }
+}
--
GitLab