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 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() ]) 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/python/dqm.py b/DQM/python/dqm.py index 7e5320c98..f65b5803d 100644 --- a/DQM/python/dqm.py +++ b/DQM/python/dqm.py @@ -293,6 +293,38 @@ 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') + + 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) : 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..ef8ede592 --- /dev/null +++ b/DQM/src/DQM/DarkBremInteraction.cxx @@ -0,0 +1,212 @@ +#include "DQM/DarkBremInteraction.h" + +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 + */ +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 + * + * @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}; + for (const double& elem : list) sum += elem*elem; + 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()); + } +} + +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" + }); +} + +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) { + 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 find 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 : 0; + + 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); + + 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); + event.add("RecoilPx", recoil_px); + event.add("RecoilPy", recoil_py); + 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)}; + 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); + + 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); +} + +} + +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); 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/_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 new file mode 100644 index 000000000..9726e7116 --- /dev/null +++ b/Validation/src/Validation/dark_brem.py @@ -0,0 +1,70 @@ +"""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]'), + ] + for h, name in features : + log.info(f'plotting {h}') + 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( + '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", + ], + density=True, + ylabel='Weighted Fraction' + ) + + 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" + ], + density=True, + ylabel='Weighted Fraction' + ) 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,