Validation of 8GeV Signal

.github/validation_samples/signal/.gitignore

@@ -0,0 +1,2 @@
+env.sh
+scratch

.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)

.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

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()
         ])

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<std::string>& 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<std::string, int> 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<int, int> known_elements_ = {
+    {1, 1},
+    {6, 2},
+    {8, 3},
+    {11, 4},
+    {14, 5},
+    {20, 6},
+    {29, 7},
+    {74, 8}
+  };
+};
+
+}

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