Added distinction between light and heavy nuclei in PDG ID plots

LDMX-Software · Nov 30, 2023 · c060401 · c060401
1 parent a32b864
commit c060401
Show file tree

Hide file tree

Showing 3 changed files with 23 additions and 10 deletions.
diff --git a/DQM/python/dqm.py b/DQM/python/dqm.py
@@ -600,20 +600,20 @@ def __init__(self, name='SampleValidation') :
         super().__init__(name, 'dqm::SampleValidation', 'DQM')
 
         # primary histograms
-        self.build1DHistogram("pdgid_primaries", "PDG ID of primary particles", 19, 0, 19)
+        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", 19, 0, 19)
+        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", 19, 0, 19)
+        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", 19, 0, 19)
+        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)
 
 

diff --git a/DQM/src/DQM/SampleValidation.cxx b/DQM/src/DQM/SampleValidation.cxx
@@ -21,7 +21,7 @@ namespace dqm {
 
     std::vector<int> primary_daughters;
 
-    double hard_thresh = 2500;
+    double hard_thresh;
 
     //Loop over all SimParticles
     for (auto const& it : particle_map) {
@@ -87,7 +87,7 @@ namespace dqm {
   }
 
   int SampleValidation::pdgid_label(const int pdgid) {
-    int label = 18; // initially assign label as "anything else"/overflow value, only change if the pdg id is something of interest
+    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-
@@ -116,11 +116,24 @@ namespace dqm {
 
     if (pdgid == 310) label = 14; // K-Short
 
-    if (pdgid == 3122 || pdgid == 3222 || pdgid == 3212 || pdgid == 3112 || pdgid == 3322 || pdgid == 3312) label = 16; // strange baryon
+    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 > 10000) label = 15; //nuclei
+    if (pdgid > 1000000000) { //nuclei
+      if (((pdgid / 10) % 1000) <= 4) {
+	label = 15; // light nuclei
+      }
+      else {
+	label = 16; // heavy nuclei
+      }
+    }
 
-    if (pdgid == 622) label = 17; // dark photon, need pdg id for other models like ALPs and SIMPs
+    if (pdgid == 622) label = 18; // dark photon, need pdg id for other models like ALPs and SIMPs
 
     return label;
   }

diff --git a/Validation/src/Validation/simparticles.py b/Validation/src/Validation/simparticles.py
@@ -7,7 +7,7 @@
 @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', 'nucleus', 'strange baryon', "A\'", 'something else'] # finish later 
+    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,