From 346ec1fb05f8c098e4a4b4f49cbb57da723e538f Mon Sep 17 00:00:00 2001
From: Deion Elgin Fellers <dfellers@lxplus791.cern.ch>
Date: Fri, 24 Feb 2023 04:46:28 +0100
Subject: [PATCH 1/2] update ntuple blinding default, add truth positions to DP
daughters, add raw calibrated energy
---
.../NtupleDumper/src/NtupleDumperAlg.cxx | 200 +++++++++---------
.../NtupleDumper/src/NtupleDumperAlg.h | 4 +-
2 files changed, 104 insertions(+), 100 deletions(-)
diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
index c06051d2a..995b69ca6 100644
--- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
+++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
@@ -170,7 +170,9 @@ StatusCode NtupleDumperAlg::initialize()
addBranch("Calo_total_nMIP", &m_calo_total_nMIP);
addBranch("Calo_total_E_dep", &m_calo_total_E_dep);
addBranch("Calo_total_E_EM", &m_calo_total_E_EM);
+ addBranch("Calo_total_raw_E_EM", &m_calo_total_raw_E_EM);
addBranch("Calo_total_E_EM_fudged", &m_calo_total_E_EM_fudged);
+ addBranch("Calo_total_raw_E_EM_fudged", &m_calo_total_raw_E_EM_fudged);
addCalibratedBranches("Preshower",2,12);
addBranch("Preshower_total_nMIP", &m_preshower_total_nMIP);
@@ -582,119 +584,114 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
if (isMC) { // if simulation find MC cross section and primary lepton
// Work out effective cross section for MC
- if (m_useFlukaWeights)
- {
- double flukaWeight = truthEventContainer->at(0)->weights()[0];
- ATH_MSG_ALWAYS("Found fluka weight = " << flukaWeight);
- m_crossSection = m_baseEventCrossSection * flukaWeight;
- }
- else if (m_useGenieWeights)
- {
- m_crossSection = m_baseEventCrossSection;
- }
- else
- {
+ if (m_useFlukaWeights) {
+ double flukaWeight = truthEventContainer->at(0)->weights()[0];
+ ATH_MSG_ALWAYS("Found fluka weight = " << flukaWeight);
+ m_crossSection = m_baseEventCrossSection * flukaWeight;
+ } else if (m_useGenieWeights) {
+ m_crossSection = m_baseEventCrossSection;
+ } else {
//ATH_MSG_WARNING("Monte carlo event with no weighting scheme specified. Setting crossSection (weight) to " << m_baseEventCrossSection << " fb.");
- m_crossSection = m_baseEventCrossSection;
+ m_crossSection = m_baseEventCrossSection;
}
- // Find the M d0 and d1 truth information
+ // Find truth particle information
SG::ReadHandle<xAOD::TruthParticleContainer> truthParticleContainer { m_truthParticleContainer, ctx };
- if (truthParticleContainer.isValid() && truthParticleContainer->size() > 0)
- {
-
+ if (truthParticleContainer.isValid() && truthParticleContainer->size() > 0) {
int ipart(0);
- for (auto particle : *truthParticleContainer)
- {
+ for (auto particle : *truthParticleContainer) { // loop over first 10 truth particles (for non A' samples)
+ if (ipart++ > 9) break;
- // loop over first 10 truth particles (for non A' samples)
-
- if (ipart++ < 10) {
-
- m_truth_P.push_back(particle->p4().P());
- m_truth_px.push_back(particle->p4().X());
- m_truth_py.push_back(particle->p4().Y());
- m_truth_pz.push_back(particle->p4().Z());
- m_truth_m.push_back(particle->m());
- m_truth_pdg.push_back(particle->pdgId());
+ m_truth_P.push_back(particle->p4().P());
+ m_truth_px.push_back(particle->p4().X());
+ m_truth_py.push_back(particle->p4().Y());
+ m_truth_pz.push_back(particle->p4().Z());
+ m_truth_m.push_back(particle->m());
+ m_truth_pdg.push_back(particle->pdgId());
if ( particle->hasProdVtx()) {
- m_truth_prod_x.push_back(particle->prodVtx()->x());
- m_truth_prod_y.push_back(particle->prodVtx()->y());
- m_truth_prod_z.push_back(particle->prodVtx()->z());
- } else {
- m_truth_prod_x.push_back(999999);
- m_truth_prod_y.push_back(999999);
- m_truth_prod_z.push_back(999999);
- }
+ m_truth_prod_x.push_back(particle->prodVtx()->x());
+ m_truth_prod_y.push_back(particle->prodVtx()->y());
+ m_truth_prod_z.push_back(particle->prodVtx()->z());
+ } else {
+ m_truth_prod_x.push_back(999999);
+ m_truth_prod_y.push_back(999999);
+ m_truth_prod_z.push_back(999999);
+ }
if ( particle->hasDecayVtx()) {
m_truth_dec_x.push_back(particle->decayVtx()->x());
m_truth_dec_y.push_back(particle->decayVtx()->y());
m_truth_dec_z.push_back(particle->decayVtx()->z());
- } else {
- m_truth_dec_x.push_back(999999);
- m_truth_dec_y.push_back(999999);
- m_truth_dec_z.push_back(999999);
+ } else {
+ m_truth_dec_x.push_back(999999);
+ m_truth_dec_y.push_back(999999);
+ m_truth_dec_z.push_back(999999);
}
- }
-
-
- if ( particle->barcode() == 1 || particle->barcode() == 2 || particle->barcode() == 3 ) {
-
- if ( particle->pdgId() == 32) { // mother particle (A')
-
- m_truthM_P.push_back(particle->p4().P());
- m_truthM_px.push_back(particle->p4().X());
- m_truthM_py.push_back(particle->p4().Y());
- m_truthM_pz.push_back(particle->p4().Z());
-
- if ( particle->hasDecayVtx()) { // decay vertex for A' particle
- m_truthM_x.push_back(particle->decayVtx()->x());
- m_truthM_y.push_back(particle->decayVtx()->y());
- m_truthM_z.push_back(particle->decayVtx()->z());
- } else {
- m_truthM_x.push_back(999999);
- m_truthM_y.push_back(999999);
- m_truthM_z.push_back(999999);
- }
- }
-
- if ( particle->pdgId() == 11) // daughter particle (positron)
- {
- m_truthd0_P.push_back(particle->p4().P());
- m_truthd0_px.push_back(particle->p4().X());
- m_truthd0_py.push_back(particle->p4().Y());
- m_truthd0_pz.push_back(particle->p4().Z());
-
- if ( particle->hasProdVtx()) {
- m_truthd0_x.push_back(particle->prodVtx()->x());
- m_truthd0_y.push_back(particle->prodVtx()->y());
- m_truthd0_z.push_back(particle->prodVtx()->z());
- } else {
- m_truthd0_x.push_back(999999);
- m_truthd0_y.push_back(999999);
- m_truthd0_z.push_back(999999);
- }
- }
- if ( particle->pdgId() == -11) // daughter particle (electron)
- {
- m_truthd1_P.push_back(particle->p4().P());
- m_truthd1_px.push_back(particle->p4().X());
- m_truthd1_py.push_back(particle->p4().Y());
- m_truthd1_pz.push_back(particle->p4().Z());
-
- if ( particle->hasProdVtx()) {
- m_truthd1_x.push_back(particle->prodVtx()->x());
- m_truthd1_y.push_back(particle->prodVtx()->y());
- m_truthd1_z.push_back(particle->prodVtx()->z());
- } else {
- m_truthd1_x.push_back(999999);
- m_truthd1_y.push_back(999999);
- m_truthd1_z.push_back(999999);
- }
- }
- }
+
+ // Find the M d0 and d1 truth information for dark photon
+ if ( particle->barcode() == 1 || particle->barcode() == 2 || particle->barcode() == 3 ) {
+ auto positions = m_fiducialParticleTool->getTruthPositions(particle->barcode());
+ if ( particle->pdgId() == 32) { // mother particle (A')
+ m_truthM_P.push_back(particle->p4().P());
+ m_truthM_px.push_back(particle->p4().X());
+ m_truthM_py.push_back(particle->p4().Y());
+ m_truthM_pz.push_back(particle->p4().Z());
+
+ if ( particle->hasDecayVtx()) { // decay vertex for A' particle
+ m_truthM_x.push_back(particle->decayVtx()->x());
+ m_truthM_y.push_back(particle->decayVtx()->y());
+ m_truthM_z.push_back(particle->decayVtx()->z());
+ } else {
+ m_truthM_x.push_back(999999);
+ m_truthM_y.push_back(999999);
+ m_truthM_z.push_back(999999);
+ }
+ }
+
+ if ( particle->pdgId() == 11) { // daughter particle (electron)
+ m_truthd0_P.push_back(particle->p4().P());
+ m_truthd0_px.push_back(particle->p4().X());
+ m_truthd0_py.push_back(particle->p4().Y());
+ m_truthd0_pz.push_back(particle->p4().Z());
+
+ if ( particle->hasProdVtx()) {
+ m_truthd0_x.push_back(particle->prodVtx()->x());
+ m_truthd0_y.push_back(particle->prodVtx()->y());
+ m_truthd0_z.push_back(particle->prodVtx()->z());
+ } else {
+ m_truthd0_x.push_back(999999);
+ m_truthd0_y.push_back(999999);
+ m_truthd0_z.push_back(999999);
+ }
+ for (int station = 1; station < 4; ++station) {
+ m_truthd0_x.push_back(positions[station].x());
+ m_truthd0_y.push_back(positions[station].y());
+ m_truthd0_z.push_back(positions[station].z());
+ }
+ }
+ if ( particle->pdgId() == -11) { // daughter particle (positron)
+ m_truthd1_P.push_back(particle->p4().P());
+ m_truthd1_px.push_back(particle->p4().X());
+ m_truthd1_py.push_back(particle->p4().Y());
+ m_truthd1_pz.push_back(particle->p4().Z());
+
+ if ( particle->hasProdVtx()) {
+ m_truthd1_x.push_back(particle->prodVtx()->x());
+ m_truthd1_y.push_back(particle->prodVtx()->y());
+ m_truthd1_z.push_back(particle->prodVtx()->z());
+ } else {
+ m_truthd1_x.push_back(999999);
+ m_truthd1_y.push_back(999999);
+ m_truthd1_z.push_back(999999);
+ }
+ for (int station = 1; station < 4; ++station) {
+ m_truthd1_x.push_back(positions[station].x());
+ m_truthd1_y.push_back(positions[station].y());
+ m_truthd1_z.push_back(positions[station].z());
+ }
+ }
+ }
}
}
}
@@ -711,6 +708,7 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
m_calo_total_nMIP += hit->Nmip();
m_calo_total_E_dep += hit->E_dep();
m_calo_total_E_EM += hit->E_EM();
+ m_calo_total_raw_E_EM += (hit->E_EM()*hit->fit_to_raw_ratio());
ATH_MSG_DEBUG("Calibrated calo: 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());
@@ -728,8 +726,10 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
// add a fudged energy variable that corrects the MC to agree with the testbeam results
m_calo_total_E_EM_fudged = m_calo_total_E_EM;
+ m_calo_total_raw_E_EM_fudged = m_calo_total_raw_E_EM;
if (isMC) {
m_calo_total_E_EM_fudged = m_calo_total_E_EM_fudged * m_caloMC_FudgeFactor;
+ m_calo_total_raw_E_EM_fudged = m_calo_total_raw_E_EM_fudged * m_caloMC_FudgeFactor;
}
// load in calibrated preshower container
@@ -1211,7 +1211,9 @@ NtupleDumperAlg::clearTree() const
m_calo_total_nMIP=0;
m_calo_total_E_dep=0;
m_calo_total_E_EM=0;
+ m_calo_total_raw_E_EM=0;
m_calo_total_E_EM_fudged=0;
+ m_calo_total_raw_E_EM_fudged=0;
m_preshower_total_nMIP=0;
m_preshower_total_E_dep=0;
diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
index 9c4a90240..322e38c90 100644
--- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
+++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
@@ -117,7 +117,7 @@ private:
DoubleProperty m_genieLuminosity { this, "GenieLuminosity", 150.0, "Genie luminosity in inverse fb." };
DoubleProperty m_minMomentum { this, "MinMomentum", 50000.0, "Write out all truth particles with a momentum larger MinMomentum"};
- DoubleProperty m_blindingCaloThreshold {this, "BlindingCaloThreshold", 25000.0, "Blind events with Ecal energy above threshold (in MeV)"};
+ DoubleProperty m_blindingCaloThreshold {this, "BlindingCaloThreshold", 100000.0, "Blind events with Ecal energy above threshold (in MeV)"};
DoubleProperty m_caloMC_FudgeFactor {this, "CaloFudgeFactorMC", 1.088, "Correct the MC energy calibration by this fudge factor"};
double m_baseEventCrossSection {1.0};
@@ -167,8 +167,10 @@ private:
mutable float m_calo_total_nMIP;
mutable float m_calo_total_E_dep;
mutable float m_calo_total_E_EM;
+ mutable float m_calo_total_raw_E_EM;
mutable float m_calo_total_E_EM_fudged;
+ mutable float m_calo_total_raw_E_EM_fudged;
mutable float m_preshower_total_nMIP;
mutable float m_preshower_total_E_dep;
--
GitLab
From b5863048c8432fd4b87a4e499a060bd216ac0f31 Mon Sep 17 00:00:00 2001
From: Deion Elgin Fellers <dfellers@lxplus791.cern.ch>
Date: Fri, 24 Feb 2023 04:57:06 +0100
Subject: [PATCH 2/2] cleanup: fix annoying tab spacing
---
.../NtupleDumper/src/NtupleDumperAlg.cxx | 56 +++++++++----------
1 file changed, 28 insertions(+), 28 deletions(-)
diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
index 995b69ca6..8b6aba734 100644
--- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
+++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
@@ -603,40 +603,40 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
if (ipart++ > 9) break;
m_truth_P.push_back(particle->p4().P());
- m_truth_px.push_back(particle->p4().X());
- m_truth_py.push_back(particle->p4().Y());
- m_truth_pz.push_back(particle->p4().Z());
- m_truth_m.push_back(particle->m());
- m_truth_pdg.push_back(particle->pdgId());
-
- if ( particle->hasProdVtx()) {
- m_truth_prod_x.push_back(particle->prodVtx()->x());
- m_truth_prod_y.push_back(particle->prodVtx()->y());
- m_truth_prod_z.push_back(particle->prodVtx()->z());
- } else {
- m_truth_prod_x.push_back(999999);
- m_truth_prod_y.push_back(999999);
- m_truth_prod_z.push_back(999999);
- }
-
- if ( particle->hasDecayVtx()) {
- m_truth_dec_x.push_back(particle->decayVtx()->x());
- m_truth_dec_y.push_back(particle->decayVtx()->y());
- m_truth_dec_z.push_back(particle->decayVtx()->z());
- } else {
- m_truth_dec_x.push_back(999999);
- m_truth_dec_y.push_back(999999);
- m_truth_dec_z.push_back(999999);
- }
+ m_truth_px.push_back(particle->p4().X());
+ m_truth_py.push_back(particle->p4().Y());
+ m_truth_pz.push_back(particle->p4().Z());
+ m_truth_m.push_back(particle->m());
+ m_truth_pdg.push_back(particle->pdgId());
+
+ if ( particle->hasProdVtx()) {
+ m_truth_prod_x.push_back(particle->prodVtx()->x());
+ m_truth_prod_y.push_back(particle->prodVtx()->y());
+ m_truth_prod_z.push_back(particle->prodVtx()->z());
+ } else {
+ m_truth_prod_x.push_back(999999);
+ m_truth_prod_y.push_back(999999);
+ m_truth_prod_z.push_back(999999);
+ }
+
+ if ( particle->hasDecayVtx()) {
+ m_truth_dec_x.push_back(particle->decayVtx()->x());
+ m_truth_dec_y.push_back(particle->decayVtx()->y());
+ m_truth_dec_z.push_back(particle->decayVtx()->z());
+ } else {
+ m_truth_dec_x.push_back(999999);
+ m_truth_dec_y.push_back(999999);
+ m_truth_dec_z.push_back(999999);
+ }
// Find the M d0 and d1 truth information for dark photon
if ( particle->barcode() == 1 || particle->barcode() == 2 || particle->barcode() == 3 ) {
auto positions = m_fiducialParticleTool->getTruthPositions(particle->barcode());
if ( particle->pdgId() == 32) { // mother particle (A')
- m_truthM_P.push_back(particle->p4().P());
+ m_truthM_P.push_back(particle->p4().P());
m_truthM_px.push_back(particle->p4().X());
- m_truthM_py.push_back(particle->p4().Y());
- m_truthM_pz.push_back(particle->p4().Z());
+ m_truthM_py.push_back(particle->p4().Y());
+ m_truthM_pz.push_back(particle->p4().Z());
if ( particle->hasDecayVtx()) { // decay vertex for A' particle
m_truthM_x.push_back(particle->decayVtx()->x());
--
GitLab