From 018a844f3438a73d3a59038ae361d1fb070a5d47 Mon Sep 17 00:00:00 2001
From: Dave Casper <dcasper@uci.edu>
Date: Wed, 22 Jun 2022 20:56:53 -0700
Subject: [PATCH] Update SimHitExample
---
.../SimHitExample/CMakeLists.txt | 3 +-
.../python/SimHitExampleConfig.py | 67 +++++++++++++++++++
.../share/SimHitExample_jobOptions.py | 4 ++
.../SimHitExample/src/SimHitAlg.cxx | 32 +++++----
.../SimHitExample/src/SimHitAlg.h | 8 +++
.../NeutrinoSimEvent/NeutrinoHit.h | 2 +-
Neutrino/NeutrinoSimEvent/src/NeutrinoHit.cxx | 13 +++-
Tracker/TrackerSimEvent/src/FaserSiHit.cxx | 11 +++
8 files changed, 125 insertions(+), 15 deletions(-)
create mode 100644 Control/CalypsoExample/SimHitExample/python/SimHitExampleConfig.py
diff --git a/Control/CalypsoExample/SimHitExample/CMakeLists.txt b/Control/CalypsoExample/SimHitExample/CMakeLists.txt
index 19aa729b..c4ff6c38 100644
--- a/Control/CalypsoExample/SimHitExample/CMakeLists.txt
+++ b/Control/CalypsoExample/SimHitExample/CMakeLists.txt
@@ -6,4 +6,5 @@ atlas_add_component( SimHitExample
LINK_LIBRARIES AthenaBaseComps GeneratorObjects TrackerSimEvent ScintSimEvent FaserCaloSimEvent NeutrinoSimEvent
)
-atlas_install_joboptions( share/*.py )
\ No newline at end of file
+atlas_install_joboptions( share/*.py )
+atlas_install_python_modules( python/*.py )
\ No newline at end of file
diff --git a/Control/CalypsoExample/SimHitExample/python/SimHitExampleConfig.py b/Control/CalypsoExample/SimHitExample/python/SimHitExampleConfig.py
new file mode 100644
index 00000000..5e71b860
--- /dev/null
+++ b/Control/CalypsoExample/SimHitExample/python/SimHitExampleConfig.py
@@ -0,0 +1,67 @@
+# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS and FASER collaborations
+
+#!/usr/bin/env python
+import sys
+from AthenaCommon.Constants import VERBOSE, INFO
+from AthenaConfiguration.ComponentFactory import CompFactory
+
+def SimHitAlgCfg(flags, name="SimHitAlg", **kwargs):
+
+ # Initialize GeoModel
+ from FaserGeoModel.FaserGeoModelConfig import FaserGeometryCfg
+ a = FaserGeometryCfg(flags)
+
+ # Configure the algorithm itself
+ SimHitAlg = CompFactory.SimHitAlg
+ a.addEventAlgo(SimHitAlg(name, **kwargs))
+
+ # Set up histogramming
+ thistSvc = CompFactory.THistSvc()
+ thistSvc.Output += ["HIST DATAFILE='simHitAlg.root' OPT='RECREATE'"]
+ a.addService(thistSvc)
+
+ return a
+
+if __name__ == "__main__":
+ from AthenaCommon.Logging import log#, logging
+ from AthenaCommon.Configurable import Configurable
+ from CalypsoConfiguration.AllConfigFlags import ConfigFlags
+
+ Configurable.configurableRun3Behavior = True
+
+# Flags for this job
+ ConfigFlags.Input.Files = ["my.HITS.pool.root"] # input file(s)
+ ConfigFlags.Input.isMC = True # Needed to bypass autoconfig
+ ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-02" # Always needed; must match FaserVersion
+ ConfigFlags.GeoModel.FaserVersion = "FASERNU-03" # Default FASER geometry
+ ConfigFlags.Detector.GeometryEmulsion = True
+ ConfigFlags.Detector.GeometryTrench = True
+ ConfigFlags.lock()
+
+# Configure components
+# Core framework
+ from CalypsoConfiguration.MainServicesConfig import MainServicesCfg
+ acc = MainServicesCfg(ConfigFlags)
+
+# Data input
+ from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
+ acc.merge(PoolReadCfg(ConfigFlags))
+
+# Algorithm
+ acc.merge(SimHitAlgCfg(ConfigFlags, McEventCollection = "TruthEvent",
+ PrintNeutrinoHits = False,
+ PrintTrackerHits = False,
+ PrintScintillatorHits = False,
+ PrintCalorimeterHits = False))
+
+# Configure verbosity
+ msgSvc = acc.getService("MessageSvc")
+ msgSvc.Format = "% F%30W%S%7W%R%T %0W%M"
+ # ConfigFlags.dump()
+ # logging.getLogger('forcomps').setLevel(VERBOSE)
+ #acc.foreach_component("*").OutputLevel = VERBOSE
+ #acc.foreach_component("*ClassID*").OutputLevel = INFO
+ #log.setLevel(VERBOSE)
+
+# Execute and finish
+ sys.exit(int(acc.run(maxEvents=-1).isFailure()))
diff --git a/Control/CalypsoExample/SimHitExample/share/SimHitExample_jobOptions.py b/Control/CalypsoExample/SimHitExample/share/SimHitExample_jobOptions.py
index 85285711..ac6bc487 100644
--- a/Control/CalypsoExample/SimHitExample/share/SimHitExample_jobOptions.py
+++ b/Control/CalypsoExample/SimHitExample/share/SimHitExample_jobOptions.py
@@ -1,3 +1,7 @@
+
+# Note: this job configuration method is obsolete
+# Please see python/SimHitAlgConfig.py for the correct method to use
+
from AthenaCommon.GlobalFlags import globalflags
globalflags.InputFormat.set_Value_and_Lock('pool')
diff --git a/Control/CalypsoExample/SimHitExample/src/SimHitAlg.cxx b/Control/CalypsoExample/SimHitExample/src/SimHitAlg.cxx
index 100b0eea..7fc37147 100644
--- a/Control/CalypsoExample/SimHitExample/src/SimHitAlg.cxx
+++ b/Control/CalypsoExample/SimHitExample/src/SimHitAlg.cxx
@@ -18,9 +18,9 @@ StatusCode SimHitAlg::initialize()
ATH_CHECK(histSvc()->regHist("/HIST/modulesSide1", m_moduleSide1));
ATH_CHECK(histSvc()->regHist("/HIST/modulesSide2", m_moduleSide2));
- m_plate_preshower = new TH2D("plate", "Scint Hit Plate", 3, -1, 1, 4, 0, 1 );
- m_plate_trigger = new TH2D("plate", "Scint Hit Plate", 3, -1, 1, 4, 0, 1 );
- m_plate_veto = new TH2D("plate", "Scint Hit Plate", 3, -1, 1, 4, 0, 1 );
+ m_plate_preshower = new TH2D("preshowerplate", "Preshower Hit", 3, -1, 1, 4, 0, 1 );
+ m_plate_trigger = new TH2D("triggerplate", "Trigger Hit", 3, -1, 1, 4, 0, 1 );
+ m_plate_veto = new TH2D("vetoplate", "Veto Hit", 3, -1, 1, 4, 0, 1 );
ATH_CHECK(histSvc()->regHist("/HIST/plate_preshower", m_plate_preshower));
ATH_CHECK(histSvc()->regHist("/HIST/plate_trigger", m_plate_trigger));
ATH_CHECK(histSvc()->regHist("/HIST/plate_veto", m_plate_veto));
@@ -28,6 +28,9 @@ StatusCode SimHitAlg::initialize()
m_ecalEnergy = new TH1D("ecalEnergy", "Ecal Energy Fraction", 100, 0.0, 0.20);
ATH_CHECK(histSvc()->regHist("/HIST/ecal_energy", m_ecalEnergy));
+ m_emulsionPDG = new TH1D("emulsionPDG", "Energy-weighted PDG code of emulsion hits", 4425, -2212.5, 2212.5);
+ ATH_CHECK(histSvc()->regHist("/HIST/emulsionPDG", m_emulsionPDG));
+
// initialize data handle keys
ATH_CHECK( m_mcEventKey.initialize() );
ATH_CHECK( m_faserSiHitKey.initialize() );
@@ -72,13 +75,13 @@ StatusCode SimHitAlg::execute()
}
SG::ReadHandle<ScintHitCollection> h_preshowerHits(m_preshowerHitKey);
- ATH_MSG_INFO("Read ScintHitCollection/Preshower with " << h_preshowerHits->size() << " hits");
+ ATH_MSG_INFO("Read ScintHitCollection/PreshowerHits with " << h_preshowerHits->size() << " hits");
SG::ReadHandle<ScintHitCollection> h_triggerHits(m_triggerHitKey);
- ATH_MSG_INFO("Read ScintHitCollection/Trigger with " << h_triggerHits->size() << " hits");
+ ATH_MSG_INFO("Read ScintHitCollection/TriggerHits with " << h_triggerHits->size() << " hits");
SG::ReadHandle<ScintHitCollection> h_vetoHits(m_vetoHitKey);
- ATH_MSG_INFO("Read ScintHitCollectionVeto with " << h_vetoHits->size() << " hits");
-
+ ATH_MSG_INFO("Read ScintHitCollection/VetoHits with " << h_vetoHits->size() << " hits");
SG::ReadHandle<CaloHitCollection> h_ecalHits(m_ecalHitKey);
+ ATH_MSG_INFO("Read CaloHitCollection with " << h_ecalHits->size() << " hits");
// Since we have no pile-up, there should always be a single GenEvent in the container
const HepMC::GenEvent* ev = (*h_mcEvents)[0];
@@ -101,7 +104,11 @@ StatusCode SimHitAlg::execute()
//Loop over all hits; print
for (const NeutrinoHit& hit : *h_emulsionHits)
{
- hit.print();
+ if (m_printNeutrino) hit.print();
+ if (hit.particleLink().isValid())
+ {
+ m_emulsionPDG->Fill(hit.particleLink()->pdg_id(), hit.energyLoss());
+ }
}
}
@@ -112,7 +119,7 @@ StatusCode SimHitAlg::execute()
// Loop over all hits; print and fill histogram
for (const FaserSiHit& hit : *h_siHits)
{
- hit.print();
+ if (m_printTracker) hit.print();
m_hist->Fill( hit.energyLoss() );
m_module->Fill( hit.getModule(), hit.getRow() );
if (hit.getSensor() == 0)
@@ -130,7 +137,7 @@ StatusCode SimHitAlg::execute()
if (h_preshowerHits->size()!=0){
for (const ScintHit& hit : *h_preshowerHits)
{
- hit.print();
+ if (m_printScintillator) hit.print();
m_hist->Fill(hit.energyLoss());
m_plate_preshower->Fill(hit.getStation(),hit.getPlate(),hit.energyLoss());
@@ -140,7 +147,7 @@ StatusCode SimHitAlg::execute()
if (h_triggerHits->size()!=0){
for (const ScintHit& hit : *h_triggerHits)
{
- hit.print();
+ if (m_printScintillator) hit.print();
m_hist->Fill(hit.energyLoss());
m_plate_trigger->Fill(hit.getStation(),hit.getPlate(),hit.energyLoss());
@@ -150,7 +157,7 @@ StatusCode SimHitAlg::execute()
if (h_vetoHits->size()!=0){
for (const ScintHit& hit : *h_vetoHits)
{
- hit.print();
+ if (m_printScintillator) hit.print();
m_hist->Fill(hit.energyLoss());
m_plate_veto->Fill(hit.getStation(),hit.getPlate(),hit.energyLoss());
@@ -162,6 +169,7 @@ StatusCode SimHitAlg::execute()
double ecalTotal = 0.0;
for (const CaloHit& hit : *h_ecalHits)
{
+ if (m_printCalorimeter) hit.print();
ecalTotal += hit.energyLoss();
}
if (ePrimary > 0) m_ecalEnergy->Fill(ecalTotal/ePrimary);
diff --git a/Control/CalypsoExample/SimHitExample/src/SimHitAlg.h b/Control/CalypsoExample/SimHitExample/src/SimHitAlg.h
index 476f2999..b385c87e 100644
--- a/Control/CalypsoExample/SimHitExample/src/SimHitAlg.h
+++ b/Control/CalypsoExample/SimHitExample/src/SimHitAlg.h
@@ -34,6 +34,9 @@ class SimHitAlg : public AthHistogramAlgorithm
// Ecal histogram
TH1* m_ecalEnergy;
+ // Emulsion PDG_ID
+ TH1* m_emulsionPDG;
+
// Read handle keys for data containers
// Any other event data can be accessed identically
// Note the key names ("BeamTruthEvent" or "SCT_Hits") are Gaudi properties and can be configured at run-time
@@ -54,4 +57,9 @@ class SimHitAlg : public AthHistogramAlgorithm
//EcalHits
SG::ReadHandleKey<CaloHitCollection> m_ecalHitKey { this, "CaloHitCollection", "EcalHits" };
+
+ BooleanProperty m_printNeutrino { this, "PrintNeutrinoHits", false };
+ BooleanProperty m_printTracker { this, "PrintTrackerHits", false };
+ BooleanProperty m_printScintillator { this, "PrintScintillatorHits", false };
+ BooleanProperty m_printCalorimeter { this, "PrintCalorimeterHits", false };
};
\ No newline at end of file
diff --git a/Neutrino/NeutrinoSimEvent/NeutrinoSimEvent/NeutrinoHit.h b/Neutrino/NeutrinoSimEvent/NeutrinoSimEvent/NeutrinoHit.h
index 1d28bd3e..2b331c5a 100644
--- a/Neutrino/NeutrinoSimEvent/NeutrinoSimEvent/NeutrinoHit.h
+++ b/Neutrino/NeutrinoSimEvent/NeutrinoSimEvent/NeutrinoHit.h
@@ -180,4 +180,4 @@ inline float hitTime(const NeutrinoHit& hit)
return (float) hit.meanTime();
}
-#endif // CALOSIMEVENT_CALOHIT_H
+#endif // NEUTRINOSIMEVENT_NEUTRINOHIT_H
diff --git a/Neutrino/NeutrinoSimEvent/src/NeutrinoHit.cxx b/Neutrino/NeutrinoSimEvent/src/NeutrinoHit.cxx
index 0101b412..2fb0fe88 100644
--- a/Neutrino/NeutrinoSimEvent/src/NeutrinoHit.cxx
+++ b/Neutrino/NeutrinoSimEvent/src/NeutrinoHit.cxx
@@ -149,8 +149,19 @@ int NeutrinoHit::getModule() const {
void NeutrinoHit::print() const {
std::cout << "*** Neutrino Hit" << std::endl;
std::cout << " Module Number " << getModule() << std::endl;
- std::cout << " Base Number " << getBase() << std::endl;
+ std::cout << " Base Number " << getBase() << std::endl;
std::cout << " Film Number " << getFilm() << std::endl;
+ std::cout << " Energy (MeV) " << energyLoss() << std::endl;
+ if (particleLink().isValid())
+ {
+ std::cout << " Barcode " << particleLink().barcode() << std::endl;
+ const HepMC::GenParticle* particle = (particleLink());
+ std::cout << " PDG ID " << particle->pdg_id() << std::endl;
+ }
+ else
+ {
+ std::cout << "Barcode/PDG ID Unknown " << std::endl;
+ }
}
int NeutrinoHit::trackNumber() const
diff --git a/Tracker/TrackerSimEvent/src/FaserSiHit.cxx b/Tracker/TrackerSimEvent/src/FaserSiHit.cxx
index 80b632ff..9fa64b3d 100644
--- a/Tracker/TrackerSimEvent/src/FaserSiHit.cxx
+++ b/Tracker/TrackerSimEvent/src/FaserSiHit.cxx
@@ -167,6 +167,17 @@ void FaserSiHit::print() const {
std::cout << " Row Number " << getRow() << std::endl;
std::cout << " Module Number " << getModule() << std::endl;
std::cout << " Sensor " << getSensor() << std::endl;
+ std::cout << " Energy (MeV) " << energyLoss() << std::endl;
+ if (particleLink().isValid())
+ {
+ std::cout << " Barcode " << particleLink().barcode() << std::endl;
+ const HepMC::GenParticle* particle = (particleLink());
+ std::cout << " PDG ID " << particle->pdg_id() << std::endl;
+ }
+ else
+ {
+ std::cout << "Barcode/PDG ID Unknown " << std::endl;
+ }
}
int FaserSiHit::trackNumber() const
--
GitLab