From 7f3d719f95bf8f56f936ecf85362ea04e8a95912 Mon Sep 17 00:00:00 2001 From: tomeichlersmith Date: Fri, 29 Sep 2023 12:50:24 -0500 Subject: [PATCH 01/10] merge old and new DQM dark brem stuff replace old ntuplizer with one that puts data in event bus better energy calculation in old, but the new adds the info to the event bus and extracts vertex material information for us --- DQM/python/dqm.py | 5 + DQM/src/DQM/DarkBremInteraction.cxx | 215 ++++++++++++++++++++ DQM/src/DQM/NtuplizeDarkBremInteraction.cxx | 151 -------------- 3 files changed, 220 insertions(+), 151 deletions(-) create mode 100644 DQM/src/DQM/DarkBremInteraction.cxx delete mode 100644 DQM/src/DQM/NtuplizeDarkBremInteraction.cxx diff --git a/DQM/python/dqm.py b/DQM/python/dqm.py index 7e5320c98..90a3a8dbd 100644 --- a/DQM/python/dqm.py +++ b/DQM/python/dqm.py @@ -293,6 +293,11 @@ def __init__(self,name='sim_dqm',sim_pass='') : self.sim_pass = sim_pass +class DarkBremInteraction(ldmxcfg.Producer) : + def __init__(self) : + super().__init__('db_kinematics','dqm::DarkBremInteraction','DQM') + + class HCalRawDigi(ldmxcfg.Analyzer) : def __init__(self, input_name) : super().__init__('hcal_pedestals','dqm::HCalRawDigi','DQM') diff --git a/DQM/src/DQM/DarkBremInteraction.cxx b/DQM/src/DQM/DarkBremInteraction.cxx new file mode 100644 index 000000000..e6f9c2c19 --- /dev/null +++ b/DQM/src/DQM/DarkBremInteraction.cxx @@ -0,0 +1,215 @@ +#include "Framework/EventProcessor.h" +#include "SimCore/Event/SimParticle.h" + +namespace dqm { + +/** + * @class DarkBremInteraction + * + * Go through the particle map and find the dark brem products, + * storing their vertex and the dark brem outgoing kinematics + * for further study. + * + * ## Products + * APrime{Px,Py,Pz} - 3-vector momentum of A' at dark brem + * APrimeEnergy - energy of A' at dark brem + * Recoil{Px,Py,Pz} - 3-vector momentum of electron recoiling from dark brem + * RecoilEnergy - energy of recoil at dark brem + * Incident{Px,Py,Pz} - 3-vector momentum of electron incident to dark brem + * IncidentEnergy - energy of incident electron at dark brem + * APrimeParentID - TrackID of A' parent + * DarkBremVertexMaterial - integer corresponding to index of known_materials + * parameter OR -1 if not found in known_materials + * DarkBremVertexMaterialZ - elemental Z value for element chosen by random from + * the elements in the material + * DarkBrem{X,Y,Z} - physical space location where dark brem occurred + */ +class DarkBremInteraction : public framework::Producer { + public: + DarkBremInteraction(const std::string& n, framework::Process& p) + : framework::Producer(n,p) {} + virtual void produce(framework::Event& e) final override; + private: + /** + * the list of known materials assiging them to material ID numbers + * + * During the simulation, we can store the name of the logical volume + * that the particle originated in. There can be many copies of logical + * volumes in different places but they all will be the same material + * by construction of how we designed our GDML. In the ecal GDML, the + * beginning the 'volume' tags list the logical volumes and you can + * see there which materials they all are in. + * + * We go through this list on each event, checking if any of these entries + * match a substring of the logical volume name stored. If we don't find any, + * the integer ID is set to -1. + * + * The inverse LUT that can be used on the plotting side is + * + * material_lut = { + * -1 : 'Unknown', + * 1 : 'C', + * 2 : 'PCB', + * 3 : 'Glue', + * 4 : 'Si', + * 5 : 'Al', + * 6 : 'W' + * } + * + * This is kind of lazy, we could instead do a full LUT where we list all known + * logical volume names and their associated materials but this analysis isn't + * as important so I haven't invested that much time in it yet. + */ + std::map known_materials_ = { + { "Carbon", 1 }, + { "PCB", 2 }, // in v12, the motherboards were simple rectangles with 'PCB' in the name + { "Glue", 3 }, + { "Si", 4 }, + { "Al", 5 }, + { "W" , 6 }, + { "strongback" , 5 }, // strongback is made of aluminum + { "motherboard" , 2 }, // motherboards are PCB + { "support" , 5 }, // support box is aluminum + { "CFMix" , 3 }, // in v12, we called the Glue layers CFMix + { "C_volume" , 1 } // in v12, we called the carbon cooling planes C but this is too general for substr matching + }; +}; + +/** + * calculate total energy from 3-momentum and mass + * + * @param[in] p 3-momentum + * @param[in] m mass + * @return total energy + */ +static double energy(const std::vector& p, const double& m) { + return sqrt(p.at(0)*p.at(0)+ p.at(1)*p.at(1)+ p.at(2)*p.at(2)+ m*m); +} + +/** + * extract the kinematics of the dark brem interaction from the SimParticles + * + * Sometimes the electron that undergoes the dark brem is not in a region + * where it should be saved (i.e. it is a shower electron inside of the ECal). + * In this case, we need to reconstruct the incident momentum from the outgoing + * products (the recoil electron and the dark photon) which should be saved by + * the biasing filter used during the simulation. + * + * Since the dark brem model does not include a nucleus, it only is able to + * conserve momentum, so we need to reconstruct the incident particle's 3-momentum + * and then use the electron mass to calculate its total energy. + */ +void DarkBremInteraction::produce(framework::Event& event) { + const auto& particle_map{event.getMap("SimParticles")}; + const ldmx::SimParticle *recoil{nullptr}, *aprime{nullptr}, *beam{nullptr}; + for (const auto& [track_id, particle] : particle_map) { + if (track_id == 1) beam = &particle; + if (particle.getProcessType() == ldmx::SimParticle::ProcessType::eDarkBrem) { + if (particle.getPdgID() == 622) { + if (aprime != nullptr) { + EXCEPTION_RAISE("BadEvent", "Found multiple A' in event."); + } + aprime = &particle; + } else { + recoil = &particle; + } + } + } + + if (recoil == nullptr and aprime == nullptr) { + /* dark brem did not occur during the simulation + * IF PROPERLY CONFIGURED, this occurs because the simulation + * exhausted the maximum number of tries to get a dark brem + * to occur. We just leave early so that the entries in the + * ntuple are the unphysical numeric minimum. + * + * This can also happen during development, so I leave a debug + * printout here to be uncommented when developing the dark + * brem simulation. + std::cout << "Event " << e.getEventNumber() + << " did not have a dark brem occur within it." << std::endl; + */ + return; + } + + if (recoil == nullptr or aprime == nullptr or beam == nullptr) { + // we are going to end processing so let's take our time to + // construct a nice error message + std::stringstream err_msg; + err_msg + << "Unable to final all necessary particles for DarkBrem interaction." + << " Missing: [ " + << (recoil == nullptr ? "recoil " : "") + << (aprime == nullptr ? "aprime " : "") + << (beam == nullptr ? "beam " : "") + << "]" << std::endl; + EXCEPTION_RAISE("BadEvent", err_msg.str()); + return; + } + + const auto& recoil_p = recoil->getMomentum(); + const auto& aprime_p = aprime->getMomentum(); + + std::vector incident_p = recoil_p; + for (std::size_t i{0}; i < recoil_p.size(); ++i) incident_p[i] += aprime_p.at(i); + + double incident_energy = energy(incident_p, recoil->getMass()); + double recoil_energy = energy(recoil_p, recoil->getMass()); + double visible_energy = (beam->getEnergy() - incident_energy) + recoil_energy; + + std::vector ap_vertex{aprime->getVertex()}; + std::string ap_vertex_volume{aprime->getVertexVolume()}; + auto ap_vertex_material_it = std::find_if( + known_materials_.begin(), known_materials_.end(), + [&](const auto& mat_pair) { + return ap_vertex_volume.find(mat_pair.first) != std::string::npos; + } + ); + int ap_vertex_material = (ap_vertex_material_it != known_materials_.end()) ? + ap_vertex_material_it->second : -1; + + int ap_parent_id{-1}; + if (aprime->getParents().size() > 0) { + ap_parent_id = aprime->getParents().at(0); + } else { + ldmx_log(error) << "Found A' without a parent ID!"; + } + + float aprime_energy = energy(aprime_p, aprime->getMass()); + int aprime_genstatus = aprime->getGenStatus(); + double aprime_px{aprime_p.at(0)}, aprime_py{aprime_p.at(1)}, aprime_pz{aprime_p.at(2)}; + event.add("APrimeEnergy", aprime_energy); + event.add("APrimePx", aprime_px); + event.add("APrimePy", aprime_py); + event.add("APrimePz", aprime_pz); + event.add("APrimeParentID", ap_parent_id); + event.add("APrimeGenStatus", aprime_genstatus); + + int recoil_genstatus = recoil->getGenStatus(); + double recoil_px{recoil_p.at(0)}, recoil_py{recoil_p.at(1)}, recoil_pz{recoil_p.at(2)}; + event.add("RecoilEnergy", recoil_energy); + event.add("RecoilPx", recoil_px); + event.add("RecoilPy", recoil_py); + event.add("RecoilPz", recoil_pz); + event.add("RecoilGenStatus", recoil_genstatus); + + event.add("IncidentEnergy", incident_energy); + double incident_px{incident_p.at(0)}, incident_py{incident_p.at(1)}, incident_pz{incident_p.at(2)}; + event.add("IncidentPx", incident_px); + event.add("IncidentPy", incident_py); + event.add("IncidentPz", incident_pz); + + double vtx_x{aprime->getVertex().at(0)}, + vtx_y{aprime->getVertex().at(1)}, + vtx_z{aprime->getVertex().at(2)}; + event.add("DarkBremX", vtx_x); + event.add("DarkBremY", vtx_y); + event.add("DarkBremZ", vtx_z); + event.add("DarkBremVertexMaterial", ap_vertex_material); + float db_material_z = event.getEventHeader().getFloatParameter("db_material_z"); + event.add("DarkBremVertexMaterialZ", db_material_z); +} + +} + +DECLARE_ANALYZER_NS(dqm,DarkBremInteraction); diff --git a/DQM/src/DQM/NtuplizeDarkBremInteraction.cxx b/DQM/src/DQM/NtuplizeDarkBremInteraction.cxx deleted file mode 100644 index 15c99bdd7..000000000 --- a/DQM/src/DQM/NtuplizeDarkBremInteraction.cxx +++ /dev/null @@ -1,151 +0,0 @@ -#include "Framework/EventProcessor.h" -#include "SimCore/Event/SimParticle.h" - -namespace dqm { -class NtuplizeDarkBremInteraction : public framework::Analyzer { - public: - NtuplizeDarkBremInteraction(const std::string& n, framework::Process& p) - : framework::Analyzer(n,p) {} - virtual void onProcessStart() final override; - virtual void analyze(const framework::Event& e) final override; -}; - -void NtuplizeDarkBremInteraction::onProcessStart() { - getHistoDirectory(); - ntuple_.create("dbint"); - ntuple_.addVar("dbint","x"); - ntuple_.addVar("dbint","y"); - ntuple_.addVar("dbint","z"); - ntuple_.addVar("dbint","weight"); - ntuple_.addVar("dbint","incident_pdg"); - ntuple_.addVar("dbint","incident_genstatus"); - ntuple_.addVar("dbint","incident_mass"); - ntuple_.addVar("dbint","incident_energy"); - ntuple_.addVar("dbint","incident_px"); - ntuple_.addVar("dbint","incident_py"); - ntuple_.addVar("dbint","incident_pz"); - ntuple_.addVar("dbint","recoil_pdg"); - ntuple_.addVar("dbint","recoil_genstatus"); - ntuple_.addVar("dbint","recoil_mass"); - ntuple_.addVar("dbint","recoil_energy"); - ntuple_.addVar("dbint","recoil_px"); - ntuple_.addVar("dbint","recoil_py"); - ntuple_.addVar("dbint","recoil_pz"); - ntuple_.addVar("dbint","aprime_pdg"); - ntuple_.addVar("dbint","aprime_genstatus"); - ntuple_.addVar("dbint","aprime_mass"); - ntuple_.addVar("dbint","aprime_energy"); - ntuple_.addVar("dbint","aprime_px"); - ntuple_.addVar("dbint","aprime_py"); - ntuple_.addVar("dbint","aprime_pz"); - ntuple_.addVar("dbint","visible_energy"); - ntuple_.addVar("dbint","beam_energy"); -} - -static double energy(const std::vector& p, const double& m) { - return sqrt(p.at(0)*p.at(0)+ p.at(1)*p.at(1)+ p.at(2)*p.at(2)+ m*m); -} - -/** - * extract the kinematics of the dark brem interaction from the SimParticles - * - * Sometimes the electron that undergoes the dark brem is not in a region - * where it should be saved (i.e. it is a shower electron inside of the ECal). - * In this case, we need to reconstruct the incident momentum from the outgoing - * products (the recoil electron and the dark photon) which should be saved by - * the biasing filter used during the simulation. - * - * Since the dark brem model does not include a nucleus, it only is able to - * conserve momentum, so we need to reconstruct the incident particle's 3-momentum - * and then use the electron mass to calculate its total energy. - */ -void NtuplizeDarkBremInteraction::analyze(const framework::Event& e) { - const auto& particle_map{e.getMap("SimParticles")}; - const ldmx::SimParticle *recoil{nullptr}, *aprime{nullptr}, *beam{nullptr}; - for (const auto& [track_id, particle] : particle_map) { - if (track_id == 1) beam = &particle; - if (particle.getProcessType() == ldmx::SimParticle::ProcessType::eDarkBrem) { - if (particle.getPdgID() == 62) { - if (aprime != nullptr) { - EXCEPTION_RAISE("BadEvent", "Found multiple A' in event."); - } - aprime = &particle; - } else { - recoil = &particle; - } - } - } - - if (recoil == nullptr and aprime == nullptr) { - /* dark brem did not occur during the simulation - * IF PROPERLY CONFIGURED, this occurs because the simulation - * exhausted the maximum number of tries to get a dark brem - * to occur. We just leave early so that the entries in the - * ntuple are the unphysical numeric minimum. - * - * This can also happen during development, so I leave a debug - * printout here to be uncommented when developing the dark - * brem simulation. - std::cout << "Event " << e.getEventNumber() - << " did not have a dark brem occur within it." << std::endl; - */ - return; - } - - if (recoil == nullptr or aprime == nullptr or beam == nullptr) { - // we are going to end processing so let's take our time to - // construct a nice error message - std::stringstream err_msg; - err_msg - << "Unable to final all necessary particles for DarkBrem interaction." - << " Missing: [ " - << (recoil == nullptr ? "recoil " : "") - << (aprime == nullptr ? "aprime " : "") - << (beam == nullptr ? "beam " : "") - << "]" << std::endl; - EXCEPTION_RAISE("BadEvent", err_msg.str()); - return; - } - - const auto& recoil_p = recoil->getMomentum(); - const auto& aprime_p = aprime->getMomentum(); - - std::vector incident_p = recoil_p; - for (std::size_t i{0}; i < recoil_p.size(); ++i) incident_p[i] += aprime_p.at(i); - - double incident_energy = energy(incident_p, recoil->getMass()); - double recoil_energy = energy(recoil_p, recoil->getMass()); - double visible_energy = (beam->getEnergy() - incident_energy) + recoil_energy; - - ntuple_.setVar("x", aprime->getVertex().at(0)); - ntuple_.setVar("y", aprime->getVertex().at(1)); - ntuple_.setVar("z", aprime->getVertex().at(2)); - ntuple_.setVar("weight", e.getEventWeight()); - ntuple_.setVar("incident_pdg", recoil->getPdgID()); - ntuple_.setVar("incident_genstatus", -1); - ntuple_.setVar("incident_mass", recoil->getMass()); - ntuple_.setVar("incident_energy", incident_energy); - ntuple_.setVar("incident_px", incident_p.at(0)); - ntuple_.setVar("incident_py", incident_p.at(1)); - ntuple_.setVar("incident_pz", incident_p.at(2)); - ntuple_.setVar("recoil_pdg", recoil->getPdgID()); - ntuple_.setVar("recoil_genstatus", recoil->getGenStatus()); - ntuple_.setVar("recoil_mass", recoil->getMass()); - ntuple_.setVar("recoil_energy", recoil_energy); - ntuple_.setVar("recoil_px", recoil_p.at(0)); - ntuple_.setVar("recoil_py", recoil_p.at(1)); - ntuple_.setVar("recoil_pz", recoil_p.at(2)); - ntuple_.setVar("aprime_pdg", aprime->getPdgID()); - ntuple_.setVar("aprime_genstatus", aprime->getGenStatus()); - ntuple_.setVar("aprime_mass", aprime->getMass()); - ntuple_.setVar("aprime_energy", energy(aprime_p,aprime->getMass())); - ntuple_.setVar("aprime_px", aprime_p.at(0)); - ntuple_.setVar("aprime_py", aprime_p.at(1)); - ntuple_.setVar("aprime_pz", aprime_p.at(2)); - ntuple_.setVar("beam_energy", beam->getEnergy()); - ntuple_.setVar("visible_energy", visible_energy); -} - -} - -DECLARE_ANALYZER_NS(dqm,NtuplizeDarkBremInteraction); From ad744a22aaa235b66ee3877f50f6e05e0d28a1b8 Mon Sep 17 00:00:00 2001 From: tomeichlersmith Date: Mon, 2 Oct 2023 12:24:07 -0500 Subject: [PATCH 02/10] update energy thresholds for target dark brem to 8GeV --- Biasing/python/target.py | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/Biasing/python/target.py b/Biasing/python/target.py index a39a36d98..833468685 100644 --- a/Biasing/python/target.py +++ b/Biasing/python/target.py @@ -177,7 +177,7 @@ def gamma_mumu( detector, generator ) : return sim -def dark_brem( ap_mass , lhe, detector ) : +def dark_brem( ap_mass , lhe, detector) : """Example configuration for producing dark brem interactions in the target. This configures the sim to fire a 4 GeV electron upstream of the @@ -214,26 +214,29 @@ def dark_brem( ap_mass , lhe, detector ) : sim.description = "One e- fired far upstream with Dark Brem turned on and biased up in target" sim.setDetector( detector , True ) - sim.generators.append( generators.single_4gev_e_upstream_tagger() ) + sim.generators.append( generators.single_8gev_e_upstream_tagger() ) sim.beamSpotSmear = [ 20., 80., 0. ] #mm #Activiate dark bremming with a certain A' mass and LHE library from LDMX.SimCore import dark_brem db_model = dark_brem.G4DarkBreMModel(lhe) - db_model.threshold = 2. #GeV - minimum energy electron needs to have to dark brem + db_model.threshold = 4. #GeV - minimum energy electron needs to have to dark brem db_model.epsilon = 0.01 #decrease epsilon from one to help with Geant4 biasing calculations sim.dark_brem.activate( ap_mass , db_model ) + import math + mass_power = max(math.log10(sim.dark_brem.ap_mass), 2.) + #Biasing dark brem up inside of the target sim.biasing_operators = [ - bias_operators.DarkBrem.target(sim.dark_brem.ap_mass**2 / db_model.epsilon**2) + bias_operators.DarkBrem.target(sim.dark_brem.ap_mass**mass_power / db_model.epsilon**2) ] sim.actions.extend([ #make sure electron reaches target with 3.5GeV - filters.TaggerVetoFilter(3500.), + filters.TaggerVetoFilter(7000.), #make sure dark brem occurs in the target where A' has at least 2GeV - filters.TargetDarkBremFilter(2000.), + filters.TargetDarkBremFilter(4000.), #keep all prodcuts of dark brem(A' and recoil electron) util.TrackProcessFilter.dark_brem() ]) From 5d4aed3bb47f91128df0f08b433488aac607240d Mon Sep 17 00:00:00 2001 From: tomeichlersmith Date: Mon, 2 Oct 2023 12:39:21 -0500 Subject: [PATCH 03/10] init a signal validation sample use mA=10MeV and use the template simulation from Biasing.target module use the init.sh script infrastructure to generate the reference lbirary before attempting the simulation --- .github/validation_samples/signal/.gitignore | 2 + .github/validation_samples/signal/config.py | 71 ++++++++++++++++++++ .github/validation_samples/signal/init.sh | 24 +++++++ 3 files changed, 97 insertions(+) create mode 100644 .github/validation_samples/signal/.gitignore create mode 100644 .github/validation_samples/signal/config.py create mode 100644 .github/validation_samples/signal/init.sh diff --git a/.github/validation_samples/signal/.gitignore b/.github/validation_samples/signal/.gitignore new file mode 100644 index 000000000..443daedf4 --- /dev/null +++ b/.github/validation_samples/signal/.gitignore @@ -0,0 +1,2 @@ +env.sh +scratch diff --git a/.github/validation_samples/signal/config.py b/.github/validation_samples/signal/config.py new file mode 100644 index 000000000..140f95fed --- /dev/null +++ b/.github/validation_samples/signal/config.py @@ -0,0 +1,71 @@ +from LDMX.Framework import ldmxcfg +p = ldmxcfg.Process('test') + +p.maxTriesPerEvent = 10000 + +from LDMX.Biasing import target +mySim = target.dark_brem( + #A' mass in MeV - set in init.sh to same value in GeV + 10.0, + # library path is uniquely determined by arguments given to `dbgen run` in init.sh + # easiest way to find this path out is by running `. init.sh` locally to see what + # is produced + 'electron_tungsten_MaxE_8.0_MinE_4.0_RelEStep_0.1_UndecayedAP_mA_0.01_run_1', + 'ldmx-det-v14-8gev' +) + +p.sequence = [ mySim ] + +################################################################## +# Below should be the same for all sim scenarios + +import os +import sys + +p.maxEvents = int(os.environ['LDMX_NUM_EVENTS']) +p.run = int(os.environ['LDMX_RUN_NUMBER']) + +p.histogramFile = f'hist.root' +p.outputFiles = [f'events.root'] + +import LDMX.Ecal.EcalGeometry +import LDMX.Ecal.ecal_hardcoded_conditions +import LDMX.Hcal.HcalGeometry +import LDMX.Hcal.hcal_hardcoded_conditions +import LDMX.Ecal.digi as ecal_digi +import LDMX.Ecal.vetos as ecal_vetos +import LDMX.Hcal.digi as hcal_digi + +from LDMX.TrigScint.trigScint import TrigScintDigiProducer +from LDMX.TrigScint.trigScint import TrigScintClusterProducer +from LDMX.TrigScint.trigScint import trigScintTrack +ts_digis = [ + TrigScintDigiProducer.pad1(), + TrigScintDigiProducer.pad2(), + TrigScintDigiProducer.pad3(), + ] +for d in ts_digis : + d.randomSeed = 1 + +from LDMX.Recon.electronCounter import ElectronCounter +from LDMX.Recon.simpleTrigger import TriggerProcessor + +count = ElectronCounter(1,'ElectronCounter') +count.input_pass_name = '' + +from LDMX.DQM import dqm + +p.sequence.extend([ + ecal_digi.EcalDigiProducer(), + ecal_digi.EcalRecProducer(), + ecal_vetos.EcalVetoProcessor(), + hcal_digi.HcalDigiProducer(), + hcal_digi.HcalRecProducer(), + *ts_digis, + TrigScintClusterProducer.pad1(), + TrigScintClusterProducer.pad2(), + TrigScintClusterProducer.pad3(), + trigScintTrack, + count, TriggerProcessor('trigger'), + dqm.DarkBremInteraction() + ] + dqm.all_dqm) diff --git a/.github/validation_samples/signal/init.sh b/.github/validation_samples/signal/init.sh new file mode 100644 index 000000000..e66a10d9d --- /dev/null +++ b/.github/validation_samples/signal/init.sh @@ -0,0 +1,24 @@ +#!/bin/bash + +############################################################################### +# init.sh +# Pre-validation initializing for dark brem signal samples +# +# We need to produce the dark brem event library. +############################################################################### + +start_group Produce Dark Brem Library +wget https://raw.githubusercontent.com/LDMX-Software/dark-brem-lib-gen/main/env.sh +source env.sh +# commented out lines are dbgen's defaults for reference +#dbgen use latest +#dbgen cache ${HOME} <- only matters for apptainer/singularity +#dbgen work /tmp +#dbgen dest $PWD +mkdir scratch +dbgen work scratch +dbgen run \ + --run 1 \ + --max_energy 8.0 \ + --apmass 0.01 +end_group From f9b8e1336930a534a23ef68e6fa1590d1a80cdca Mon Sep 17 00:00:00 2001 From: tomeichlersmith Date: Mon, 2 Oct 2023 15:36:15 -0500 Subject: [PATCH 04/10] add most prominent histograms for PR Validations --- DQM/python/dqm.py | 27 ++++++ DQM/src/DQM/DarkBremInteraction.cxx | 124 +++++++++++++++++++++++++++- 2 files changed, 148 insertions(+), 3 deletions(-) diff --git a/DQM/python/dqm.py b/DQM/python/dqm.py index 90a3a8dbd..f65b5803d 100644 --- a/DQM/python/dqm.py +++ b/DQM/python/dqm.py @@ -297,6 +297,33 @@ class DarkBremInteraction(ldmxcfg.Producer) : def __init__(self) : super().__init__('db_kinematics','dqm::DarkBremInteraction','DQM') + self.build1DHistogram('aprime_energy', + 'Dark Photon Energy [MeV]',101,0,8080) + self.build1DHistogram('aprime_pt', + 'Dark Photon pT [MeV]',100,0,2000) + + self.build1DHistogram('recoil_energy', + 'Recoil Electron Energy [MeV]',101,0,8080) + self.build1DHistogram('recoil_pt', + 'Recoil Electron pT [MeV]',100,0,2000) + + self.build1DHistogram('incident_energy', + 'Incident Electron Energy [MeV]',101,0,8080) + self.build1DHistogram('incident_pt', + 'Incident Electron pT [MeV]',100,0,2000) + + # weird binning so we can see the target and trigger pads + self.build1DHistogram('dark_brem_z', + 'Z Location of Dark Brem [mm]', + [-5.0, -4.6752, -3.5502, -2.4252, -1.3002, -0.1752, 0.1752, 1.]) + # elements are hydrogen and carbon (for trigger pads) and tungsten target + self.build1DHistogram('dark_brem_element', + 'Element in which Dark Brem Occurred', + 10, 0, 10) + self.build1DHistogram('dark_brem_material', + 'Material in which Dark Brem Occurred', + 8, 0, 8) + class HCalRawDigi(ldmxcfg.Analyzer) : def __init__(self, input_name) : diff --git a/DQM/src/DQM/DarkBremInteraction.cxx b/DQM/src/DQM/DarkBremInteraction.cxx index e6f9c2c19..c278b8d3b 100644 --- a/DQM/src/DQM/DarkBremInteraction.cxx +++ b/DQM/src/DQM/DarkBremInteraction.cxx @@ -28,8 +28,13 @@ class DarkBremInteraction : public framework::Producer { public: DarkBremInteraction(const std::string& n, framework::Process& p) : framework::Producer(n,p) {} + virtual void onProcessStart() final override; virtual void produce(framework::Event& e) final override; private: + /** + * Set the labels of the histogram of the input name with the input labels + */ + void setHistLabels(const std::string& name, const std::vector& labels); /** * the list of known materials assiging them to material ID numbers * @@ -47,13 +52,14 @@ class DarkBremInteraction : public framework::Producer { * The inverse LUT that can be used on the plotting side is * * material_lut = { - * -1 : 'Unknown', + * 0 : 'Unknown', * 1 : 'C', * 2 : 'PCB', * 3 : 'Glue', * 4 : 'Si', * 5 : 'Al', - * 6 : 'W' + * 6 : 'W', + * 7 : 'PVT' * } * * This is kind of lazy, we could instead do a full LUT where we list all known @@ -67,12 +73,49 @@ class DarkBremInteraction : public framework::Producer { { "Si", 4 }, { "Al", 5 }, { "W" , 6 }, + { "target", 6 }, + { "trigger_pad", 7 }, { "strongback" , 5 }, // strongback is made of aluminum { "motherboard" , 2 }, // motherboards are PCB { "support" , 5 }, // support box is aluminum { "CFMix" , 3 }, // in v12, we called the Glue layers CFMix { "C_volume" , 1 } // in v12, we called the carbon cooling planes C but this is too general for substr matching }; + + /** + * The list of known elements assigning them to the bins that we are putting them into. + * + * There are two failure modes for this: + * 1. The dark brem didn't happen, in which case, the element reported by the event header + * will be -1. We give this an ID of 0. + * 2. The dark brem occurred within an element not listed here, in which case we give it + * the last bin. + * + * The inverset LUT that can be used if studying the output tree is + * + * element_lut = { + * 0 : 'did_not_happen', + * 1 : 'H 1', + * 2 : 'C 6', + * 3 : 'O 8', + * 4 : 'Na 11', + * 5 : 'Si 14', + * 6 : 'Ca 20', + * 7 : 'Cu 29', + * 8 : 'W 74', + * 9 : 'unlisted' + * } + */ + std::map known_elements_ = { + {1, 1}, + {6, 2}, + {8, 3}, + {11, 4}, + {14, 5}, + {20, 6}, + {29, 7}, + {74, 8} + }; }; /** @@ -86,6 +129,58 @@ static double energy(const std::vector& p, const double& m) { return sqrt(p.at(0)*p.at(0)+ p.at(1)*p.at(1)+ p.at(2)*p.at(2)+ m*m); } +/** + * calculate the sum in quadrature of the passed list of doubles + */ +static double quadsum(const std::initializer_list& list) { + double sum{0}; + for (const double& elem : list) sum += elem*elem; + return sqrt(sum); +} + + +void DarkBremInteraction::setHistLabels( + const std::string& name, + const std::vector& labels) { + auto h{histograms_.get(name)}; + for (std::size_t ibin{1}; ibin <= labels.size(); ibin++) { + h->GetXaxis()->SetBinLabel(ibin, labels[ibin-1].c_str()); + } +} + +/** + * update the labels of some categorical histograms + */ +void DarkBremInteraction::onProcessStart() { + setHistLabels( + "dark_brem_material", + { + "Unknown", + "C", + "PCB", + "Glue", + "Si", + "Al", + "W", + "PVT" + }); + + setHistLabels( + "dark_brem_element", + { + "did not happen", + "H 1", + "C 6", + "O 8", + "Na 11", + "Si 14", + "Ca 20", + "Cu 29", + "W 74", + "unlisted" + }); +} + /** * extract the kinematics of the dark brem interaction from the SimParticles * @@ -166,7 +261,7 @@ void DarkBremInteraction::produce(framework::Event& event) { } ); int ap_vertex_material = (ap_vertex_material_it != known_materials_.end()) ? - ap_vertex_material_it->second : -1; + ap_vertex_material_it->second : 0; int ap_parent_id{-1}; if (aprime->getParents().size() > 0) { @@ -185,6 +280,9 @@ void DarkBremInteraction::produce(framework::Event& event) { event.add("APrimeParentID", ap_parent_id); event.add("APrimeGenStatus", aprime_genstatus); + histograms_.fill("aprime_energy", aprime_energy); + histograms_.fill("aprime_pt", quadsum({aprime_px, aprime_py})); + int recoil_genstatus = recoil->getGenStatus(); double recoil_px{recoil_p.at(0)}, recoil_py{recoil_p.at(1)}, recoil_pz{recoil_p.at(2)}; event.add("RecoilEnergy", recoil_energy); @@ -193,12 +291,18 @@ void DarkBremInteraction::produce(framework::Event& event) { event.add("RecoilPz", recoil_pz); event.add("RecoilGenStatus", recoil_genstatus); + histograms_.fill("recoil_energy", recoil_energy); + histograms_.fill("recoil_pt", quadsum({recoil_px, recoil_py})); + event.add("IncidentEnergy", incident_energy); double incident_px{incident_p.at(0)}, incident_py{incident_p.at(1)}, incident_pz{incident_p.at(2)}; event.add("IncidentPx", incident_px); event.add("IncidentPy", incident_py); event.add("IncidentPz", incident_pz); + histograms_.fill("incident_energy", incident_energy); + histograms_.fill("incident_pt", quadsum({incident_px, incident_py})); + double vtx_x{aprime->getVertex().at(0)}, vtx_y{aprime->getVertex().at(1)}, vtx_z{aprime->getVertex().at(2)}; @@ -208,6 +312,20 @@ void DarkBremInteraction::produce(framework::Event& event) { event.add("DarkBremVertexMaterial", ap_vertex_material); float db_material_z = event.getEventHeader().getFloatParameter("db_material_z"); event.add("DarkBremVertexMaterialZ", db_material_z); + + histograms_.fill("dark_brem_z", vtx_z); + + int i_element = 0; + if (db_material_z > 0) { + if (known_elements_.find(static_cast(db_material_z)) == known_elements_.end()) { + i_element = known_elements_.size(); + } else { + i_element = known_elements_.at(static_cast(db_material_z)); + } + } + + histograms_.fill("dark_brem_element", i_element); + histograms_.fill("dark_brem_material", ap_vertex_material); } } From 1f61df88fa99f0a40a7b7c3ad10c9760dec0167e Mon Sep 17 00:00:00 2001 From: tomeichlersmith Date: Tue, 3 Oct 2023 10:15:41 -0500 Subject: [PATCH 05/10] add plotting of dark brem DQM histograms not updating labels for material/elements but that's how it goes --- Validation/src/Validation/__init__.py | 1 + Validation/src/Validation/dark_brem.py | 32 ++++++++++++++++++++++++++ 2 files changed, 33 insertions(+) create mode 100644 Validation/src/Validation/dark_brem.py diff --git a/Validation/src/Validation/__init__.py b/Validation/src/Validation/__init__.py index 4ada6e544..e75aa91c1 100644 --- a/Validation/src/Validation/__init__.py +++ b/Validation/src/Validation/__init__.py @@ -6,3 +6,4 @@ from . import trigscint from . import hcal from . import photonuclear +from . import dark_brem diff --git a/Validation/src/Validation/dark_brem.py b/Validation/src/Validation/dark_brem.py new file mode 100644 index 000000000..5f7b5d579 --- /dev/null +++ b/Validation/src/Validation/dark_brem.py @@ -0,0 +1,32 @@ +"""Plotting of ECal-related validation plots""" + +from ._differ import Differ +from ._plotter import plotter +import logging + +log = logging.getLogger('dark_brem') + +@plotter(hist=True,event=False) +def kinematics(d : Differ, out_dir = None) : + """Plot Dark Brem interaction histograms + + Parameters + ---------- + d : Differ + Differ containing files that are not event files (presumably histogram files) + """ + + features = [ + ('aprime_energy', 'Dark Photon Energy [MeV]'), + ('aprime_pt', 'Dark Photon $p_T$ [MeV]'), + ('recoil_energy', 'Recoil Energy [MeV]'), + ('recoil_pt', 'Recoil $p_T$ [MeV]'), + ('incident_energy', 'Incident Energy [MeV]'), + ('incident_pt', 'Incident $p_T$ [MeV]'), + ('dark_brem_z', 'Dark Brem Z Location [mm]'), + ('dark_brem_element', 'Element in which Dark Brem Occurred'), + ('dark_brem_material', 'Material in which Dark Brem Occurred') + ] + for h, name in features : + log.info(f'plotting {h}') + d.plot1d(f'db_kinematics/db_kinematics_{h}', name, out_dir = out_dir) From a7a199b7f7c68ddf6b27e5c08596758b616ed19d Mon Sep 17 00:00:00 2001 From: tomeichlersmith Date: Tue, 3 Oct 2023 10:28:05 -0500 Subject: [PATCH 06/10] add tick labels to element/material dark brem histograms - move tick labels to bin centers so the number of labels is equal to the number of bins (and not the number of bin edges) - update hcal and photonuclear modules to conform to this change, removing the trailing empty label present in those tick label lists --- Validation/src/Validation/_differ.py | 2 +- Validation/src/Validation/dark_brem.py | 38 +++++++++++++++++++++-- Validation/src/Validation/hcal.py | 2 +- Validation/src/Validation/photonuclear.py | 6 ++-- 4 files changed, 41 insertions(+), 7 deletions(-) diff --git a/Validation/src/Validation/_differ.py b/Validation/src/Validation/_differ.py index 54c39a013..fe6ca1bd8 100644 --- a/Validation/src/Validation/_differ.py +++ b/Validation/src/Validation/_differ.py @@ -121,7 +121,7 @@ def plot1d(self, column, xlabel, legend_kw['title'] = self.grp_name if tick_labels is not None: - ax.set_xticks(bins) + ax.set_xticks((bins[1:]+bins[:-1])/2) ax.set_xticklabels(tick_labels) ax.tick_params(axis='x', rotation=90) diff --git a/Validation/src/Validation/dark_brem.py b/Validation/src/Validation/dark_brem.py index 5f7b5d579..9ba6e7364 100644 --- a/Validation/src/Validation/dark_brem.py +++ b/Validation/src/Validation/dark_brem.py @@ -24,9 +24,43 @@ def kinematics(d : Differ, out_dir = None) : ('incident_energy', 'Incident Energy [MeV]'), ('incident_pt', 'Incident $p_T$ [MeV]'), ('dark_brem_z', 'Dark Brem Z Location [mm]'), - ('dark_brem_element', 'Element in which Dark Brem Occurred'), - ('dark_brem_material', 'Material in which Dark Brem Occurred') ] for h, name in features : log.info(f'plotting {h}') d.plot1d(f'db_kinematics/db_kinematics_{h}', name, out_dir = out_dir) + + log.info('plotting dark_brem_element') + d.plot1d( + 'db_kinematics/db_kinematics_dark_brem_element', + 'Element in which Dark Brem Occurred', + out_dir = out_dir, + tick_labels = [ + "did not happen", + "H 1", + "C 6", + "O 8", + "Na 11", + "Si 14", + "Ca 20", + "Cu 29", + "W 74", + "unlisted", + ] + ) + + log.info('plotting dark_brem_material') + d.plot1d( + 'db_kinematics/db_kinematics_dark_brem_material', + 'Material in which Dark Brem Occurred', + out_dir = out_dir, + tick_labels = [ + "Unknown", + "C", + "PCB", + "Glue", + "Si", + "Al", + "W", + "PVT" + ] + ) diff --git a/Validation/src/Validation/hcal.py b/Validation/src/Validation/hcal.py index 8ab24733e..6c9deb7d8 100644 --- a/Validation/src/Validation/hcal.py +++ b/Validation/src/Validation/hcal.py @@ -42,7 +42,7 @@ def dqm(d: Differ, out_dir=None): tick_labels=['', 'Back', 'Top', 'Bottom', 'Right', 'Left', 'Any', 'Both', 'Back only', 'Side only', - 'Neither', '', ''], + 'Neither', ''], out_dir=out_dir, yscale='linear', ) diff --git a/Validation/src/Validation/photonuclear.py b/Validation/src/Validation/photonuclear.py index b43e612cf..d8f97bc00 100644 --- a/Validation/src/Validation/photonuclear.py +++ b/Validation/src/Validation/photonuclear.py @@ -11,10 +11,10 @@ def pndqm(d: Differ, out_dir=None): event_type_labels = ['', 'Nothing hard', 'n', 'nn', '≥ 3n', 'π', 'ππ', 'π₀', 'πA', 'π2A', 'ππA', 'π₀A', 'π₀2A', 'π₀πA', 'p', 'pp', 'pn', 'K_L⁰X', 'KX', - 'K_S⁰X', 'exotics', 'multi-body', '', '', ''] + 'K_S⁰X', 'exotics', 'multi-body', '', ''] - compact_event_type_labels = ['', 'n', 'K^{±}X', 'K⁰', 'nn', 'soft', 'other', '',''] - neutron_event_type_labels = ['', '', 'nn', 'pn', 'π^+n', 'π⁰n', '', ''] + compact_event_type_labels = ['', 'n', 'K^{±}X', 'K⁰', 'nn', 'soft', 'other', ''] + neutron_event_type_labels = ['', '', 'nn', 'pn', 'π^+n', 'π⁰n', ''] d.plot1d("PN/PN_event_type", "Event category (200 MeV cut)", tick_labels=event_type_labels, From c16889a9f53626cfab83c326348549391cb52662 Mon Sep 17 00:00:00 2001 From: tomeichlersmith Date: Tue, 3 Oct 2023 11:19:36 -0500 Subject: [PATCH 07/10] make sure to weight the histograms by the event weight this does not affect the distribution shapes greatly since the weights are extremely close to exactly 1/B within each sample --- DQM/src/DQM/DarkBremInteraction.cxx | 1 + 1 file changed, 1 insertion(+) diff --git a/DQM/src/DQM/DarkBremInteraction.cxx b/DQM/src/DQM/DarkBremInteraction.cxx index c278b8d3b..47e785f9a 100644 --- a/DQM/src/DQM/DarkBremInteraction.cxx +++ b/DQM/src/DQM/DarkBremInteraction.cxx @@ -195,6 +195,7 @@ void DarkBremInteraction::onProcessStart() { * and then use the electron mass to calculate its total energy. */ void DarkBremInteraction::produce(framework::Event& event) { + histograms_.setWeight(event.getEventHeader().getWeight()); const auto& particle_map{event.getMap("SimParticles")}; const ldmx::SimParticle *recoil{nullptr}, *aprime{nullptr}, *beam{nullptr}; for (const auto& [track_id, particle] : particle_map) { From ad64e869c4d05f05dae97ca512448062d762a477 Mon Sep 17 00:00:00 2001 From: tomeichlersmith Date: Mon, 9 Oct 2023 10:31:00 -0500 Subject: [PATCH 08/10] move declaration of DarkBremInteraction into its own header --- DQM/include/DQM/DarkBremInteraction.h | 145 ++++++++++++++++++++++++++ DQM/src/DQM/DarkBremInteraction.cxx | 144 ++----------------------- 2 files changed, 154 insertions(+), 135 deletions(-) create mode 100644 DQM/include/DQM/DarkBremInteraction.h diff --git a/DQM/include/DQM/DarkBremInteraction.h b/DQM/include/DQM/DarkBremInteraction.h new file mode 100644 index 000000000..c5ba3d255 --- /dev/null +++ b/DQM/include/DQM/DarkBremInteraction.h @@ -0,0 +1,145 @@ +#include "Framework/EventProcessor.h" +#include "SimCore/Event/SimParticle.h" + +namespace dqm { +/** + * @class DarkBremInteraction + * + * Go through the particle map and find the dark brem products, + * storing their vertex and the dark brem outgoing kinematics + * for further study. + * + * While histograms are filled to be automatically validated and plotted, + * we also put these values into the event tree so users can look at the + * variables related to the dark brem in detail. + * + * ## Products + * APrime{Px,Py,Pz} - 3-vector momentum of A' at dark brem + * APrimeEnergy - energy of A' at dark brem + * Recoil{Px,Py,Pz} - 3-vector momentum of electron recoiling from dark brem + * RecoilEnergy - energy of recoil at dark brem + * Incident{Px,Py,Pz} - 3-vector momentum of electron incident to dark brem + * IncidentEnergy - energy of incident electron at dark brem + * APrimeParentID - TrackID of A' parent + * DarkBremVertexMaterial - integer corresponding to index of known_materials + * parameter OR -1 if not found in known_materials + * DarkBremVertexMaterialZ - elemental Z value for element chosen by random from + * the elements in the material + * DarkBrem{X,Y,Z} - physical space location where dark brem occurred + */ +class DarkBremInteraction : public framework::Producer { + public: + DarkBremInteraction(const std::string& n, framework::Process& p) + : framework::Producer(n,p) {} + /** + * update the labels of some categorial histograms + * + * This is helpful for downstream viewers of the histograms + * so that ROOT will display the bins properly. + */ + virtual void onProcessStart() final override; + + /** + * extract the kinematics of the dark brem interaction from the SimParticles + * + * Sometimes the electron that undergoes the dark brem is not in a region + * where it should be saved (i.e. it is a shower electron inside of the ECal). + * In this case, we need to reconstruct the incident momentum from the outgoing + * products (the recoil electron and the dark photon) which should be saved by + * the biasing filter used during the simulation. + * + * Since the dark brem model does not include a nucleus, it only is able to + * conserve momentum, so we need to reconstruct the incident particle's 3-momentum + * and then use the electron mass to calculate its total energy. + */ + virtual void produce(framework::Event& e) final override; + private: + /** + * Set the labels of the histogram of the input name with the input labels + */ + void setHistLabels(const std::string& name, const std::vector& labels); + + /** + * the list of known materials assigning them to material ID numbers + * + * During the simulation, we can store the name of the logical volume + * that the particle originated in. There can be many copies of logical + * volumes in different places but they all will be the same material + * by construction of how we designed our GDML. In the ecal GDML, the + * beginning the 'volume' tags list the logical volumes and you can + * see there which materials they all are in. + * + * We go through this list on each event, checking if any of these entries + * match a substring of the logical volume name stored. If we don't find any, + * the integer ID is set to -1. + * + * The inverse LUT that can be used on the plotting side is + * + * material_lut = { + * 0 : 'Unknown', + * 1 : 'C', + * 2 : 'PCB', + * 3 : 'Glue', + * 4 : 'Si', + * 5 : 'Al', + * 6 : 'W', + * 7 : 'PVT' + * } + * + * This is kind of lazy, we could instead do a full LUT where we list all known + * logical volume names and their associated materials but this analysis isn't + * as important so I haven't invested that much time in it yet. + */ + std::map known_materials_ = { + { "Carbon", 1 }, + { "PCB", 2 }, // in v12, the motherboards were simple rectangles with 'PCB' in the name + { "Glue", 3 }, + { "Si", 4 }, + { "Al", 5 }, + { "W" , 6 }, + { "target", 6 }, + { "trigger_pad", 7 }, + { "strongback" , 5 }, // strongback is made of aluminum + { "motherboard" , 2 }, // motherboards are PCB + { "support" , 5 }, // support box is aluminum + { "CFMix" , 3 }, // in v12, we called the Glue layers CFMix + { "C_volume" , 1 } // in v12, we called the carbon cooling planes C but this is too general for substr matching + }; + + /** + * The list of known elements assigning them to the bins that we are putting them into. + * + * There are two failure modes for this: + * 1. The dark brem didn't happen, in which case, the element reported by the event header + * will be -1. We give this an ID of 0. + * 2. The dark brem occurred within an element not listed here, in which case we give it + * the last bin. + * + * The inverset LUT that can be used if studying the output tree is + * + * element_lut = { + * 0 : 'did_not_happen', + * 1 : 'H 1', + * 2 : 'C 6', + * 3 : 'O 8', + * 4 : 'Na 11', + * 5 : 'Si 14', + * 6 : 'Ca 20', + * 7 : 'Cu 29', + * 8 : 'W 74', + * 9 : 'unlisted' + * } + */ + std::map known_elements_ = { + {1, 1}, + {6, 2}, + {8, 3}, + {11, 4}, + {14, 5}, + {20, 6}, + {29, 7}, + {74, 8} + }; +}; + +} diff --git a/DQM/src/DQM/DarkBremInteraction.cxx b/DQM/src/DQM/DarkBremInteraction.cxx index 47e785f9a..d10368240 100644 --- a/DQM/src/DQM/DarkBremInteraction.cxx +++ b/DQM/src/DQM/DarkBremInteraction.cxx @@ -1,123 +1,7 @@ -#include "Framework/EventProcessor.h" -#include "SimCore/Event/SimParticle.h" +#include "DQM/DarkBremInteraction.h" namespace dqm { -/** - * @class DarkBremInteraction - * - * Go through the particle map and find the dark brem products, - * storing their vertex and the dark brem outgoing kinematics - * for further study. - * - * ## Products - * APrime{Px,Py,Pz} - 3-vector momentum of A' at dark brem - * APrimeEnergy - energy of A' at dark brem - * Recoil{Px,Py,Pz} - 3-vector momentum of electron recoiling from dark brem - * RecoilEnergy - energy of recoil at dark brem - * Incident{Px,Py,Pz} - 3-vector momentum of electron incident to dark brem - * IncidentEnergy - energy of incident electron at dark brem - * APrimeParentID - TrackID of A' parent - * DarkBremVertexMaterial - integer corresponding to index of known_materials - * parameter OR -1 if not found in known_materials - * DarkBremVertexMaterialZ - elemental Z value for element chosen by random from - * the elements in the material - * DarkBrem{X,Y,Z} - physical space location where dark brem occurred - */ -class DarkBremInteraction : public framework::Producer { - public: - DarkBremInteraction(const std::string& n, framework::Process& p) - : framework::Producer(n,p) {} - virtual void onProcessStart() final override; - virtual void produce(framework::Event& e) final override; - private: - /** - * Set the labels of the histogram of the input name with the input labels - */ - void setHistLabels(const std::string& name, const std::vector& labels); - /** - * the list of known materials assiging them to material ID numbers - * - * During the simulation, we can store the name of the logical volume - * that the particle originated in. There can be many copies of logical - * volumes in different places but they all will be the same material - * by construction of how we designed our GDML. In the ecal GDML, the - * beginning the 'volume' tags list the logical volumes and you can - * see there which materials they all are in. - * - * We go through this list on each event, checking if any of these entries - * match a substring of the logical volume name stored. If we don't find any, - * the integer ID is set to -1. - * - * The inverse LUT that can be used on the plotting side is - * - * material_lut = { - * 0 : 'Unknown', - * 1 : 'C', - * 2 : 'PCB', - * 3 : 'Glue', - * 4 : 'Si', - * 5 : 'Al', - * 6 : 'W', - * 7 : 'PVT' - * } - * - * This is kind of lazy, we could instead do a full LUT where we list all known - * logical volume names and their associated materials but this analysis isn't - * as important so I haven't invested that much time in it yet. - */ - std::map known_materials_ = { - { "Carbon", 1 }, - { "PCB", 2 }, // in v12, the motherboards were simple rectangles with 'PCB' in the name - { "Glue", 3 }, - { "Si", 4 }, - { "Al", 5 }, - { "W" , 6 }, - { "target", 6 }, - { "trigger_pad", 7 }, - { "strongback" , 5 }, // strongback is made of aluminum - { "motherboard" , 2 }, // motherboards are PCB - { "support" , 5 }, // support box is aluminum - { "CFMix" , 3 }, // in v12, we called the Glue layers CFMix - { "C_volume" , 1 } // in v12, we called the carbon cooling planes C but this is too general for substr matching - }; - - /** - * The list of known elements assigning them to the bins that we are putting them into. - * - * There are two failure modes for this: - * 1. The dark brem didn't happen, in which case, the element reported by the event header - * will be -1. We give this an ID of 0. - * 2. The dark brem occurred within an element not listed here, in which case we give it - * the last bin. - * - * The inverset LUT that can be used if studying the output tree is - * - * element_lut = { - * 0 : 'did_not_happen', - * 1 : 'H 1', - * 2 : 'C 6', - * 3 : 'O 8', - * 4 : 'Na 11', - * 5 : 'Si 14', - * 6 : 'Ca 20', - * 7 : 'Cu 29', - * 8 : 'W 74', - * 9 : 'unlisted' - * } - */ - std::map known_elements_ = { - {1, 1}, - {6, 2}, - {8, 3}, - {11, 4}, - {14, 5}, - {20, 6}, - {29, 7}, - {74, 8} - }; -}; - /** * calculate total energy from 3-momentum and mass * @@ -131,6 +15,8 @@ static double energy(const std::vector& p, const double& m) { /** * calculate the sum in quadrature of the passed list of doubles + * + * @param[in] list `{`-bracket-enclosed list of doubles to square, sum, and square-root */ static double quadsum(const std::initializer_list& list) { double sum{0}; @@ -138,19 +24,20 @@ static double quadsum(const std::initializer_list& list) { return sqrt(sum); } - void DarkBremInteraction::setHistLabels( const std::string& name, const std::vector& labels) { + /** + * We could probably move this into Framework since it is a relatively + * common task. I could even imagine a way of constructing a StrCategory + * histogram. + */ auto h{histograms_.get(name)}; for (std::size_t ibin{1}; ibin <= labels.size(); ibin++) { h->GetXaxis()->SetBinLabel(ibin, labels[ibin-1].c_str()); } } -/** - * update the labels of some categorical histograms - */ void DarkBremInteraction::onProcessStart() { setHistLabels( "dark_brem_material", @@ -181,19 +68,6 @@ void DarkBremInteraction::onProcessStart() { }); } -/** - * extract the kinematics of the dark brem interaction from the SimParticles - * - * Sometimes the electron that undergoes the dark brem is not in a region - * where it should be saved (i.e. it is a shower electron inside of the ECal). - * In this case, we need to reconstruct the incident momentum from the outgoing - * products (the recoil electron and the dark photon) which should be saved by - * the biasing filter used during the simulation. - * - * Since the dark brem model does not include a nucleus, it only is able to - * conserve momentum, so we need to reconstruct the incident particle's 3-momentum - * and then use the electron mass to calculate its total energy. - */ void DarkBremInteraction::produce(framework::Event& event) { histograms_.setWeight(event.getEventHeader().getWeight()); const auto& particle_map{event.getMap("SimParticles")}; @@ -233,7 +107,7 @@ void DarkBremInteraction::produce(framework::Event& event) { // construct a nice error message std::stringstream err_msg; err_msg - << "Unable to final all necessary particles for DarkBrem interaction." + << "Unable to find all necessary particles for DarkBrem interaction." << " Missing: [ " << (recoil == nullptr ? "recoil " : "") << (aprime == nullptr ? "aprime " : "") From 158057592b1d3029ae9415b16b71a5fa39859428 Mon Sep 17 00:00:00 2001 From: tomeichlersmith Date: Mon, 9 Oct 2023 10:54:14 -0500 Subject: [PATCH 09/10] use density since the histogram's weight scales are so different --- Validation/src/Validation/dark_brem.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/Validation/src/Validation/dark_brem.py b/Validation/src/Validation/dark_brem.py index 9ba6e7364..9726e7116 100644 --- a/Validation/src/Validation/dark_brem.py +++ b/Validation/src/Validation/dark_brem.py @@ -27,7 +27,7 @@ def kinematics(d : Differ, out_dir = None) : ] for h, name in features : log.info(f'plotting {h}') - d.plot1d(f'db_kinematics/db_kinematics_{h}', name, out_dir = out_dir) + d.plot1d(f'db_kinematics/db_kinematics_{h}', name, out_dir = out_dir, density=True, ylabel='Weighted Fraction') log.info('plotting dark_brem_element') d.plot1d( @@ -45,7 +45,9 @@ def kinematics(d : Differ, out_dir = None) : "Cu 29", "W 74", "unlisted", - ] + ], + density=True, + ylabel='Weighted Fraction' ) log.info('plotting dark_brem_material') @@ -62,5 +64,7 @@ def kinematics(d : Differ, out_dir = None) : "Al", "W", "PVT" - ] + ], + density=True, + ylabel='Weighted Fraction' ) From 1de89db49c9cec55eca02904c9078556e9204c95 Mon Sep 17 00:00:00 2001 From: Tom Eichlersmith <31970302+tomeichlersmith@users.noreply.github.com> Date: Tue, 10 Oct 2023 11:49:41 -0500 Subject: [PATCH 10/10] add rational for existence of energy function explain in doxygen comment why this function exists basically copying the rational written for the produce function as well. --- DQM/src/DQM/DarkBremInteraction.cxx | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/DQM/src/DQM/DarkBremInteraction.cxx b/DQM/src/DQM/DarkBremInteraction.cxx index d10368240..ef8ede592 100644 --- a/DQM/src/DQM/DarkBremInteraction.cxx +++ b/DQM/src/DQM/DarkBremInteraction.cxx @@ -5,6 +5,10 @@ namespace dqm { /** * calculate total energy from 3-momentum and mass * + * Since the dark brem model does not include a nucleus, it only is able to + * conserve momentum, so we need to reconstruct the incident particle's 3-momentum + * and then use the known particle mass to calculate its total energy. + * * @param[in] p 3-momentum * @param[in] m mass * @return total energy