From 261b11235c4f397fa51f1f54d5cc8a8d8cf7f655 Mon Sep 17 00:00:00 2001
From: FASER Reco <faserrec@lxplus928.cern.ch>
Date: Sun, 3 Nov 2024 00:52:43 +0100
Subject: [PATCH] Fixes to ntuple
---
.../Reconstruction/scripts/faser_reco.py | 4 +-
.../scripts/faser_ntuple_maker.py | 12 +-
.../NtupleDumper/src/NtupleDumperAlg.cxx | 109 +++++++++++++++---
.../NtupleDumper/src/NtupleDumperAlg.h | 3 +
4 files changed, 105 insertions(+), 23 deletions(-)
diff --git a/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py b/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
index f2d32291..c5fc6d5f 100755
--- a/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
+++ b/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
@@ -188,7 +188,9 @@ if len(args.cond):
if useCKF:
# Enable ACTS material corrections, this crashes testbeam geometries
- configFlags.TrackingGeometry.MaterialSource = "/cvmfs/faser.cern.ch/repo/sw/database/DBRelease/current/acts/material-maps.json"
+ # Use alma-9 specific file
+ print(">>> Using alma9 specific ACTS material file! <<<")
+ configFlags.TrackingGeometry.MaterialSource = "/cvmfs/faser.cern.ch/repo/sw/database/DBRelease/current/acts/material-maps-alma9.json"
# Must use original input string here, as pathlib mangles double // in path names
configFlags.Input.Files = [ args.file_path ]
diff --git a/PhysicsAnalysis/NtupleDumper/scripts/faser_ntuple_maker.py b/PhysicsAnalysis/NtupleDumper/scripts/faser_ntuple_maker.py
index 83d6e54c..f137bc27 100755
--- a/PhysicsAnalysis/NtupleDumper/scripts/faser_ntuple_maker.py
+++ b/PhysicsAnalysis/NtupleDumper/scripts/faser_ntuple_maker.py
@@ -348,16 +348,19 @@ if args.isMC:
else:
acc.merge(NtupleDumperAlgCfg(configFlags, outfile, DoBlinding=(not args.unblind), OnlyBlinded=args.onlyblind, DoScintFilter = args.scintFilt, DoTrackFilter = (not args.NoTrackFilt), DoTrigFilter = args.trigFilt, StableOnly = (not args.no_stable), DoRandomFilter = args.randomTrigFilt, DoRandomOnlyFilter=args.randomOnlyTrigFilt ,**grl_kwargs))
+from AthenaConfiguration.ComponentFactory import CompFactory
+AthenaEventLoopMgr = CompFactory.AthenaEventLoopMgr()
+
if not args.verbose:
- from AthenaConfiguration.ComponentFactory import CompFactory
- AthenaEventLoopMgr = CompFactory.AthenaEventLoopMgr()
AthenaEventLoopMgr.EventPrintoutInterval=1000
- acc.addService(AthenaEventLoopMgr)
else:
nd = acc.getEventAlgo("NtupleDumperAlg")
nd.OutputLevel = VERBOSE
+ AthenaEventLoopMgr.EventPrintoutInterval=1
+acc.addService(AthenaEventLoopMgr)
+
# Hack to avoid problem with our use of MC databases when isMC = False
if not args.isMC:
replicaSvc = acc.getService("DBReplicaSvc")
@@ -368,11 +371,12 @@ if not args.isMC:
if args.verbose:
log.setLevel(VERBOSE)
+ acc.getService("MessageSvc").debugLimit = 100000
acc.printConfig(withDetails=True)
configFlags.dump()
else:
log.setLevel(INFO)
-
+
acc.getService("MessageSvc").Format = "% F%40W%S%7W%R%T %0W%M"
# Execute and finish
diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
index d964b1d2..02c4a266 100644
--- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
+++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
@@ -47,25 +47,32 @@ void NtupleDumperAlg::addBranch(const std::string &name,
m_tree->Branch(name.c_str(),var,(name+"/I").c_str());
}
-// Have to declare this cost to call during execute
-void NtupleDumperAlg::defineWaveBranches() const {
+std::map<std::string, std::list<int>>&
+NtupleDumperAlg::getChannelMap() const {
- ATH_MSG_DEBUG("defineWaveBranches called");
-
// Use waveform map to find all defined waveform channels
auto mapping = m_mappingTool->getCableMapping();
ATH_MSG_DEBUG("Cable mapping contains " << mapping.size() << " entries");
// mapping is std::map<int, std::pair<std::string, Identifier> >
// use this to fill my own map of channel lists keyed by type
-
- std::map<std::string, std::list<int>> wave_map;
+ static std::map<std::string, std::list<int>> wave_map;
for (const auto& [key, value] : mapping) {
wave_map[value.first].push_back(key);
ATH_MSG_DEBUG("Found mapping " << value.first << " chan " << key);
}
+ return wave_map;
+}
+
+// Have to declare this cost to call during execute
+void NtupleDumperAlg::defineWaveBranches() const {
+
+ ATH_MSG_DEBUG("defineWaveBranches called");
+
+ auto wave_map = getChannelMap();
+
// Now go through found types and define ntuple entries
// Keys are defined in cable map and used by RawWaveformDecoder
for (const auto& [key, value] : wave_map) {
@@ -574,6 +581,10 @@ StatusCode NtupleDumperAlg::initialize()
ATH_MSG_INFO("Blinding will NOT be enforced for real data.");
}
+ // Make sure charge histogram pointers are empty
+ for (unsigned int chan = 0; chan<max_chan; chan++)
+ m_HistRandomCharge[chan] = NULL;
+
return StatusCode::SUCCESS;
}
@@ -593,7 +604,6 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
if (truthEventContainer.isValid() && truthEventContainer->size() > 0)
{
isMC = true;
-
}
SG::ReadHandle<McEventCollection> mcEventContainer {m_mcEventContainer, ctx};
@@ -732,8 +742,10 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
// Fill histograms
for (unsigned int chan = 0; chan<max_chan; chan++) {
// Only fill histograms that have been defined
- if (m_HistRandomCharge[chan])
+ if (m_HistRandomCharge[chan]) {
+ ATH_MSG_DEBUG("Filling random charge histogram chan= " << chan);
m_HistRandomCharge[chan]->Fill(m_wave_raw_charge[chan]);
+ }
}
}
@@ -770,18 +782,81 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
}
+
if (m_doScintFilter) {
+ // Get channel mapping
+ auto wave_map = getChannelMap();
+
// filter events, but use waveform peak cuts instead of triggers, as triggers could miss signals slightly out of time
// digitizer channels described here: https://twiki.cern.ch/twiki/bin/viewauth/FASER/CaloScintillator
- bool calo_trig = (m_wave_raw_peak[0] > 3.0 || m_wave_raw_peak[1] > 3.0 || m_wave_raw_peak[2] > 3.0 || m_wave_raw_peak[3] > 3.0);
- bool vetoNu_trig = (m_wave_raw_peak[4] > 25.0 && m_wave_raw_peak[5] > 25.0);
- bool vetoSt2_trig = (m_wave_raw_peak[6] > 25.0 && m_wave_raw_peak[7] > 25.0);
- bool timingScint_trig = ((m_wave_raw_peak[8] > 25.0 && m_wave_raw_peak[9] > 25.0) || (m_wave_raw_peak[10] > 25.0 && m_wave_raw_peak[11] > 25.0));
- bool preshower_trig = (m_wave_raw_peak[12] > 3.0 && m_wave_raw_peak[13] > 3.0);
- bool vetoSt1_trig = (m_wave_raw_peak[14] > 25.0);
+ // Make calorimeter trigger, use CaloHi if it exists
+ bool calo_trig = false;
+ if (wave_map.contains(std::string("calo2"))) {
+ for (auto chan : wave_map[std::string("calo2")]) {
+ calo_trig |= m_wave_raw_peak[chan] > 3.0;
+ }
+ }
+ else if (wave_map.contains(std::string("calo"))) {
+ for (auto chan : wave_map[std::string("calo")]) {
+ calo_trig |= m_wave_raw_peak[chan] > 3.0;
+ }
+ }
- bool veto_OR_trig = (vetoNu_trig || vetoSt1_trig || vetoSt2_trig);
+ // VetoNu trigger
+ bool vetoNu_trig = true;
+ if (wave_map.contains(std::string("vetonu"))) {
+ for (auto chan : wave_map[std::string("vetonu")]) {
+ vetoNu_trig &= m_wave_raw_peak[chan] > 25.0;
+ }
+ } else {
+ vetoNu_trig = false;
+ }
+ // Veto Station 2
+ bool vetoSt1_trig = true;
+ bool vetoSt2_trig = true;
+ if (wave_map.contains(std::string("veto"))) {
+ int ncount = 0;
+ for (auto chan : wave_map[std::string("veto")]) {
+ if (++ncount < 3) {
+ vetoSt2_trig &= m_wave_raw_peak[chan] > 25.0;
+ } else {
+ vetoSt1_trig &= m_wave_raw_peak[chan] > 25.0;
+ }
+ }
+ } else {
+ vetoSt1_trig = false;
+ vetoSt2_trig = false;
+ }
+ bool veto_OR_trig = vetoNu_trig || vetoSt1_trig || vetoSt2_trig;
+
+ // Timing layer (either top or bottom)
+ bool timingScint_trig = true;
+ bool timingScint_trig1 = true;
+ bool timingScint_trig2 = true;
+ if (wave_map.contains(std::string("trigger"))) {
+ int ncount = 0;
+ for (auto chan : wave_map[std::string("trigger")]) {
+ if (++ncount < 3) {
+ timingScint_trig1 &= m_wave_raw_peak[chan] > 25.0;
+ } else {
+ timingScint_trig2 &= m_wave_raw_peak[chan] > 25.0;
+ }
+ }
+ timingScint_trig = timingScint_trig1 || timingScint_trig2;
+ } else {
+ timingScint_trig = false;
+ }
+
+ bool preshower_trig = true;
+ if (wave_map.contains(std::string("preshower"))) {
+ for (auto chan : wave_map[std::string("preshower")]) {
+ preshower_trig &= m_wave_raw_peak[chan] > 3.0;
+ }
+ } else {
+ preshower_trig = false;
+ }
+
bool passes_ScintFilter = false;
if (calo_trig) {
passes_ScintFilter = true;
@@ -1052,7 +1127,6 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
m_calo_total_E_EM_fudged = m_calo_total_E_EM_fudged * m_caloMC_FudgeFactor;
}
-
// load in calibrated preshower container
SG::ReadHandle<xAOD::CalorimeterHitContainer> preshowerCalibratedContainer { m_preshowerCalibratedContainer, ctx };
ATH_CHECK(preshowerCalibratedContainer.isValid());
@@ -1066,8 +1140,7 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
ATH_MSG_DEBUG("Calibrated preshower: ch is " << ch << ", edep is " << hit->E_dep() << ", E_EM is " << hit->E_EM() << ", Nmip is " << hit->Nmip() << ", fit_to_raw_ratio is " << hit->fit_to_raw_ratio());
}
-
-
+
// enforce blinding such that events that miss the veto layers and have a large calo signal are skipped and not in the output root file
// Only blind colliding BCIDs (+/- 1)
bool blinded_event = false;
diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
index c5051ae7..85d36d67 100644
--- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
+++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
@@ -71,6 +71,9 @@ private:
void addCalibratedBranches(const std::string &name, int nchannels, int first);
double radius(const Acts::Vector3 &position) const;
+ // Return map with string as key and list<int> as value
+ std::map<std::string, std::list<int>>& getChannelMap() const;
+
ServiceHandle <ITHistSvc> m_histSvc;
SG::ReadHandleKey<xAOD::TruthEventContainer> m_truthEventContainer { this, "EventContainer", "TruthEvents", "Truth event container name." };
--
GitLab