diff --git a/DQM/include/DQM/SampleValidation.h b/DQM/include/DQM/SampleValidation.h new file mode 100644 index 000000000..b8ef8e26f --- /dev/null +++ b/DQM/include/DQM/SampleValidation.h @@ -0,0 +1,31 @@ +#ifndef DQM_SAMPLEVALIDATION_H +#define DQM_SAMPLEVALIDATION_H + +//LDMX Framework +#include "Framework/Configure/Parameters.h" +#include "Framework/EventProcessor.h" + +namespace dqm { + + /** + * @class SampleValidation + * @brief + */ + + class SampleValidation : public framework::Analyzer { + public: + + SampleValidation(const std::string& name, framework::Process& process) : Analyzer(name, process) {} + + virtual void configure(framework::config::Parameters& ps); + + virtual void analyze(const framework::Event& event); + + int pdgid_label(const int pdgid); + private: + + }; + +} + +#endif diff --git a/DQM/python/dqm.py b/DQM/python/dqm.py index ecb71a0b7..240bfb860 100644 --- a/DQM/python/dqm.py +++ b/DQM/python/dqm.py @@ -594,6 +594,27 @@ def __init__(self,name='Trigger',coll='Trigger') : self.trigger_name = coll self.trigger_pass = '' + +class SampleValidation(ldmxcfg.Analyzer) : + def __init__(self, name='SampleValidation') : + super().__init__(name, 'dqm::SampleValidation', 'DQM') + + # primary histograms + self.build1DHistogram("pdgid_primaries", "PDG ID of primary particles", 20, 0, 20) + self.build1DHistogram("energy_primaries", "Energy of primary particles", 90, 0, 9000) # range applicable for 4 GeV beam + self.build2DHistogram("beam_smear", "x", 30, -150, 150, "y", 30, -150, 150) + self.build1DHistogram("pdgid_primarydaughters", "PDG ID of primary daughtesr", 20, 0, 20) + self.build1DHistogram("energy_daughterphoton", "Energy spectrum of all photons from primary", 170, 0, 8500) + + # primary daughter of interest (brem / dark brem) histograms + self.build1DHistogram("pdgid_harddaughters", "PDG ID of primary daughters", 20, 0, 20) + self.build1DHistogram("startZ_hardbrem", "Start z position of hard primary daughter", 100, -500, 500) + self.build1DHistogram("endZ_hardbrem", "End z position of hard primary daughter", 100, -500, 500) + self.build1DHistogram("energy_hardbrem", "Energy spectrum of hard primary daughter", 130, 2000, 8500) + + # daughters of hard brem histograms + self.build1DHistogram("pdgid_hardbremdaughters", "PDG ID of hard brem daughters", 20, 0, 20) + self.build1DHistogram("startZ_hardbremdaughters", "Start z position of hard brem daughters", 200, -1000, 1000) ecal_dqm = [ diff --git a/DQM/src/DQM/SampleValidation.cxx b/DQM/src/DQM/SampleValidation.cxx new file mode 100644 index 000000000..a03652c66 --- /dev/null +++ b/DQM/src/DQM/SampleValidation.cxx @@ -0,0 +1,143 @@ +#include "DQM/SampleValidation.h" +#include "SimCore/Event/SimParticle.h" +#include "SimCore/Event/SimTrackerHit.h" +#include "Framework/NtupleManager.h" +#include +#include +#include + +namespace dqm { + + void SampleValidation::configure(framework::config::Parameters& ps) { + + return; + } + + void SampleValidation::analyze(const framework::Event& event) { + + //Grab the SimParticle Map and Target Scoring Plane Hits + auto targetSPHits(event.getCollection("TargetScoringPlaneHits")); + auto particle_map{event.getMap("SimParticles")}; + + std::vector primary_daughters; + + double hard_thresh; + + //Loop over all SimParticles + for (auto const& it : particle_map) { + ldmx::SimParticle p = it.second; + int pdgid = p.getPdgID(); + std::vector vertex = p.getVertex(); + double energy = p.getEnergy(); + std::vector parents_track_ids = p.getParents(); + std::vector daughters = p.getDaughters(); + + for (auto const& parent_track_id: parents_track_ids) { + if (parent_track_id == 0) { + histograms_.fill("pdgid_primaries", pdgid_label(pdgid)); + histograms_.fill("energy_primaries", energy); + hard_thresh = (2500/4000)*energy; + primary_daughters = daughters; + for (const ldmx::SimTrackerHit &sphit : targetSPHits) { + if (sphit.getTrackID() == it.first && sphit.getPosition()[2] < 0) { + histograms_.fill("beam_smear", vertex[0], vertex[1]); + } + } + } + } + } + + std::vector> hardbrem_daughters; + + for (auto const& it : particle_map) { + int trackid = it.first; + ldmx::SimParticle p = it.second; + for (auto const& primary_daughter : primary_daughters) { + if (trackid == primary_daughter) { + histograms_.fill("pdgid_primarydaughters", pdgid_label(p.getPdgID())); + if (p.getPdgID() == 22) { + histograms_.fill("energy_daughterphoton", p.getEnergy()); + } + if (p.getEnergy() >= hard_thresh) { + histograms_.fill("pdgid_harddaughters", pdgid_label(p.getPdgID())); + histograms_.fill("startZ_hardbrem", p.getVertex()[2]); + histograms_.fill("endZ_hardbrem", p.getEndPoint()[2]); + histograms_.fill("energy_hardbrem", p.getEnergy()); + hardbrem_daughters.push_back(p.getDaughters()); + } + } + } + } + + for (auto const& it : particle_map) { + int trackid = it.first; + ldmx::SimParticle p = it.second; + for (const std::vector &daughter_track_id : hardbrem_daughters){ + for (const int &daughter_id : daughter_track_id) { + if (trackid == daughter_id) { + histograms_.fill("pdgid_hardbremdaughters", pdgid_label(p.getPdgID())); + histograms_.fill("startZ_hardbremdaughters", p.getVertex()[2]); + } + } + } + } + + + return; + } + + int SampleValidation::pdgid_label(const int pdgid) { + int label = 19; // initially assign label as "anything else"/overflow value, only change if the pdg id is something of interest + if (pdgid == -11) label = 1; // e+ + + if (pdgid == 11) label = 2; // e- + + if (pdgid == -13) label = 3; // μ+ + + if (pdgid == 13) label = 4; // μ- + + if (pdgid == 22) label = 5; // γ + + if (pdgid == 2212) label = 6; // proton + + if (pdgid == 2112) label = 7; // neutron + + if (pdgid == 211) label = 8; //π+ + + if (pdgid == -211) label = 9; //π- + + if (pdgid == 111) label = 10; //π0 + + if (pdgid == 321) label = 11; // K+ + + if (pdgid == -321) label = 12; // K- + + if (pdgid == 130) label = 13; // K-Long + + if (pdgid == 310) label = 14; // K-Short + + if (pdgid == 3122 || pdgid == 3222 || pdgid == 3212 || pdgid == 3112 || pdgid == 3322 || pdgid == 3312) label = 17; // strange baryon + + /* + * Nuclear PDG codes are given by ±10LZZZAAAI so to find the atomic + * number, we divide by 10 (to lose I) and then take the modulo + * with 1000. + */ + + if (pdgid > 1000000000) { //nuclei + if (((pdgid / 10) % 1000) <= 4) { + label = 15; // light nuclei + } + else { + label = 16; // heavy nuclei + } + } + + if (pdgid == 622) label = 18; // dark photon, need pdg id for other models like ALPs and SIMPs + + return label; + } + +} + +DECLARE_ANALYZER_NS(dqm, SampleValidation) diff --git a/Validation/src/Validation/__init__.py b/Validation/src/Validation/__init__.py index e75aa91c1..76d6a9e63 100644 --- a/Validation/src/Validation/__init__.py +++ b/Validation/src/Validation/__init__.py @@ -7,3 +7,4 @@ from . import hcal from . import photonuclear from . import dark_brem +from . import simparticles diff --git a/Validation/src/Validation/simparticles.py b/Validation/src/Validation/simparticles.py new file mode 100644 index 000000000..af15bb65b --- /dev/null +++ b/Validation/src/Validation/simparticles.py @@ -0,0 +1,53 @@ +from ._differ import Differ +from ._plotter import plotter +import logging + +log = logging.getLogger('8GeV') + +@plotter(hist=True, event=False) +def beamenergy_comp(d: Differ, out_dir=None): + + pdgid_labels = ['', 'e+', 'e-', 'μ+', 'μ-', 'γ', 'p', 'n', 'π+', 'π-', 'π0', 'K+', 'K-', 'K-L', 'K-S', 'light nucleus', 'heavy nucleus', 'strange baryon', "A\'", 'something else'] # finish later + + d.plot1d("SampleValidation/SampleValidation_pdgid_primaries", "PDG ID, primaries", + tick_labels=pdgid_labels, + out_dir=out_dir, + density=True) + + d.plot1d("SampleValidation/SampleValidation_energy_primaries", "Energy of primaries (MeV)", + out_dir=out_dir, + density=True) + + d.plot1d("SampleValidation/SampleValidation_pdgid_primarydaughters", "PDG ID, primary daughters", + tick_labels=pdgid_labels, + out_dir=out_dir, + density=True) + d.plot1d("SampleValidation/SampleValidation_energy_daughterphoton", "Energy spectrum of all photons from primary", + out_dir=out_dir, + density=True) + + d.plot1d("SampleValidation/SampleValidation_pdgid_harddaughters", "PDG ID of hard primary daughter", + tick_labels=pdgid_labels, + out_dir=out_dir, + density=True) + + d.plot1d("SampleValidation/SampleValidation_startZ_hardbrem", "Start z position of hard primary daughter", + out_dir=out_dir, + density=True) + + d.plot1d("SampleValidation/SampleValidation_endZ_hardbrem", "End z position of hard primary daughter", + out_dir=out_dir, + density=True) + + d.plot1d("SampleValidation/SampleValidation_energy_hardbrem", "Energy spectrum of hard primary daughter", + out_dir=out_dir, + density=True) + + d.plot1d("SampleValidation/SampleValidation_pdgid_hardbremdaughters", "PDG ID, hard brem daughters", + tick_labels=pdgid_labels, + out_dir=out_dir, + density=True) + + d.plot1d("SampleValidation/SampleValidation_startZ_hardbremdaughters", "Start z position of hard brem daughters", + out_dir=out_dir, + density=True)