Sample scripts for validating future samples

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

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", 19, 0, 19)
+        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("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("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("startZ_hardbremdaughters", "Start z position of hard brem daughters", 200, -1000, 1000)
 
 
 ecal_dqm = [

DQM/src/DQM/SampleValidation.cxx

@@ -0,0 +1,130 @@
+#include "DQM/SampleValidation.h"
+#include "SimCore/Event/SimParticle.h"
+#include "SimCore/Event/SimTrackerHit.h"
+#include "Framework/NtupleManager.h"
+#include <iostream>
+#include <fstream>
+#include <algorithm>
+
+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<ldmx::SimTrackerHit>("TargetScoringPlaneHits"));
+    auto particle_map{event.getMap<int, ldmx::SimParticle>("SimParticles")};
+
+    std::vector<int> primary_daughters;
+
+    double hard_thresh = 2500;
+
+    //Loop over all SimParticles
+    for (auto const& it : particle_map) {
+      ldmx::SimParticle p = it.second;
+      int pdgid = p.getPdgID();
+      std::vector<double> vertex = p.getVertex();
+      double energy = p.getEnergy();
+      std::vector<int> parents_track_ids = p.getParents();
+      std::vector<int> 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<std::vector<int>> 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<int> &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 = 18; // 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 = 16; // strange baryon
+
+    if (pdgid > 10000) label = 15; //nuclei
+
+    if (pdgid == 622) label = 17; // dark photon, need pdg id for other models like ALPs and SIMPs
+
+    return label;
+  }
+
+}
+
+DECLARE_ANALYZER_NS(dqm, SampleValidation)

Validation/src/Validation/__init__.py

@@ -7,3 +7,4 @@
 from . import hcal
 from . import photonuclear
 from . import dark_brem
+from . import simparticles

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', '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)