Transition to using G4DarkBreM

[tomeichlersmith]

This was relatively simple and so the C++ updates were changed quickly. I think this will work pretty well since the validation repo is modeled after this one, but I want to test and make sure.

To Do

• [x] Make sure running is operational
• [x] Test target and ecal dark brem results are physical (small scale)
• [x] Document breaking changes that will need to be applied downstream in Biasing

[tomeichlersmith]

In order to get this to compile (no run test yet), I needed to update two files in the Biasing module of ldmx-sw

diff --git a/Biasing/src/Biasing/EcalDarkBremFilter.cxx b/Biasing/src/Biasing/EcalDarkBremFilter.cxx
index 58157ff..78276b9 100644
--- a/Biasing/src/Biasing/EcalDarkBremFilter.cxx
+++ b/Biasing/src/Biasing/EcalDarkBremFilter.cxx
@@ -10,8 +10,8 @@
 #include "Biasing/EcalDarkBremFilter.h"

 #include "G4LogicalVolumeStore.hh"      //for the store
-#include "SimCore/DarkBrem/G4APrime.h"  //checking if particles match A'
-#include "SimCore/DarkBrem/G4eDarkBremsstrahlung.h"  //checking for dark brem secondaries
+#include "G4DarkBreM/G4APrime.h"  //checking if particles match A'
+#include "G4DarkBreM/G4DarkBremsstrahlung.h"  //checking for dark brem secondaries
 #include "SimCore/UserTrackInformation.h"            //make sure A' is saved

 namespace biasing {
@@ -60,7 +60,7 @@ void EcalDarkBremFilter::BeginOfEventAction(const G4Event*) {

 G4ClassificationOfNewTrack EcalDarkBremFilter::ClassifyNewTrack(
     const G4Track* aTrack, const G4ClassificationOfNewTrack& cl) {
-  if (aTrack->GetParticleDefinition() == simcore::darkbrem::G4APrime::APrime()) {
+  if (aTrack->GetParticleDefinition() == G4APrime::APrime()) {
     // there is an A'! Yay!
     /* Debug message
     std::cout << "[ EcalDarkBremFilter ]: "
@@ -98,12 +98,12 @@ void EcalDarkBremFilter::PostUserTrackingAction(const G4Track* track) {

   const G4VProcess* creator = track->GetCreatorProcess();
   if (creator and creator->GetProcessName().contains(
-                      simcore::darkbrem::G4eDarkBremsstrahlung::PROCESS_NAME)) {
+                      G4DarkBremsstrahlung::PROCESS_NAME)) {
     // make sure all secondaries of dark brem process are saved
     simcore::UserTrackInformation* userInfo = simcore::UserTrackInformation::get(track);
     // make sure A' is persisted into output file
     userInfo->setSaveFlag(true);
-    if (track->GetParticleDefinition() == simcore::darkbrem::G4APrime::APrime()) {
+    if (track->GetParticleDefinition() == G4APrime::APrime()) {
       // check if A' was made in the desired volume and has the minimum energy
       if (not inDesiredVolume(track)) {
         AbortEvent("A' wasn't produced inside of the requested volume.");
diff --git a/Biasing/src/Biasing/TargetDarkBremFilter.cxx b/Biasing/src/Biasing/TargetDarkBremFilter.cxx
index 35fc16f..07e362d 100644
--- a/Biasing/src/Biasing/TargetDarkBremFilter.cxx
+++ b/Biasing/src/Biasing/TargetDarkBremFilter.cxx
@@ -10,7 +10,7 @@
 #include "Biasing/TargetDarkBremFilter.h"

 #include "G4Electron.hh"                   //to check if track is electron
-#include "SimCore/DarkBrem/G4APrime.h"     //checking if particles match A'
+#include "G4DarkBreM/G4APrime.h"     //checking if particles match A'
 #include "SimCore/UserTrackInformation.h"  //make sure A' is saved

 namespace biasing {
@@ -32,7 +32,7 @@ void TargetDarkBremFilter::stepping(const G4Step* step) {
   // Leave if track is not an electron
   auto particle_def{track->GetParticleDefinition()};
   if (particle_def != G4Electron::Electron()) {
-    if (particle_def == simcore::darkbrem::G4APrime::APrime() and
+    if (particle_def == G4APrime::APrime() and
         track->GetCurrentStepNumber() == 1) {
       /** check on first step of A' to see if it originated in correct volume
        * this needs to be here because sometimes the
@@ -71,8 +71,7 @@ void TargetDarkBremFilter::stepping(const G4Step* step) {
     } else {
       // check secondaries to see if we made a dark brem
       for (auto& secondary_track : *secondaries) {
-        if (secondary_track->GetParticleDefinition() ==
-            simcore::darkbrem::G4APrime::APrime()) {
+        if (secondary_track->GetParticleDefinition() == G4APrime::APrime()) {
           // we found an A', woo-hoo!

           if (secondary_track->GetTotalEnergy() < threshold_) {

[tomeichlersmith]

The updated xsec calculation avoids the need for the more complicated biasing factor

diff --git a/Biasing/python/eat.py b/Biasing/python/eat.py
index 3f193f3..278dad5 100644
--- a/Biasing/python/eat.py
+++ b/Biasing/python/eat.py
@@ -129,13 +129,9 @@ def dark_brem( ap_mass , lhe, detector ) :
     sim.dark_brem.activate( ap_mass , db_model )

     #Biasing dark brem up inside of the ecal volumes
-    from math import log10
-    #need a higher power for the higher mass A'
-    mass_power = max(log10(sim.dark_brem.ap_mass),2.)
-
     from LDMX.SimCore import bias_operators
     sim.biasing_operators = [ 
-            bias_operators.DarkBrem.ecal(sim.dark_brem.ap_mass**mass_power / db_model.epsilon**2)
+            bias_operators.DarkBrem.ecal(sim.dark_brem.ap_mass**2 / db_model.epsilon**2)
             ]

     sim.actions = [
diff --git a/Biasing/python/target.py b/Biasing/python/target.py
index 9bff37f..2515276 100644
--- a/Biasing/python/target.py
+++ b/Biasing/python/target.py
@@ -222,11 +222,8 @@ def dark_brem( ap_mass , lhe, detector ) :
     sim.dark_brem.activate( ap_mass , db_model )

     #Biasing dark brem up inside of the target
-    #need to bias up high mass A' by more than 2 so that they can actually happen
-    from math import log10
-    mass_power = max(log10(sim.dark_brem.ap_mass),2.)
     sim.biasing_operators = [
-            bias_operators.DarkBrem.target(sim.dark_brem.ap_mass**mass_power / db_model.epsilon**2)
+            bias_operators.DarkBrem.target(sim.dark_brem.ap_mass**2 / db_model.epsilon**2)
             ]

     sim.actions.extend([

[tomeichlersmith]

Should rename VertextLibrary to G4DarkBreM in Biasing python modules as well for future folks to use as example.

[tomeichlersmith]

Validation of Transition to using G4DarkBreM

The main method of validation in this case is to check the weight and vertex position distributions. The kinematics are validated within G4DarkBreM itself and the associated paper ArXiV:2211.03873.

In order to do this validation, I am using the dark brem event library shipped with G4DarkBreM (4 GeV electrons on tungsten with $m_{A'} = 100$ MeV) to simulate 10k events within LDMX. Both target and ECal samples are generated.

The transition also required some patches to areas outside of SimCore, so in order to re-run this test you will need to checkout a specific branch of ldmx-sw and SimCore.

cd ldmx-s
git fetch
git checkout SimCore-83-external-g4db
cd SimCore
git fetch
git checkout 83-external-g4db

%matplotlib inline
import matplotlib.pyplot as plt
import mplhep
plt.style.use(mplhep.style.ROOT)

import uproot
import pandas as pd

def plt_hist(x, *,
             xlabel = None, ylabel = None,
             yscale = 'log',
             **hist_kwargs) :
    if 'histtype' not in hist_kwargs :
        hist_kwargs['histtype'] = 'step'
    if 'linewidth' not in hist_kwargs :
        hist_kwargs['linewidth'] = 2
    plt.hist(x, **hist_kwargs)
    plt.xlabel(xlabel)
    if ylabel is None :
        ylabel = 'Event Count'
        if 'weights' in hist_kwargs :
            ylabel = 'Weighted '+ylabel
    plt.ylabel(ylabel)
    plt.yscale(yscale)
    plt.show()

ECal Sample

The ECal sample is a "thick target" scenario, so we want a falling exponential distribution for the event weights and distributions in position coordinates that roughly follow the shower development.

This sample was generated with the config copied here.

import os
from LDMX.Framework import ldmxcfg
p = ldmxcfg.Process('ecal_db')
from LDMX.Biasing import eat
from LDMX.Ecal import EcalGeometry
from LDMX.Hcal import HcalGeometry
p.sequence = [
    eat.dark_brem( 
        100., #MeV - mass of A'
        os.environ['LDMX_BASE']+'/ldmx-sw/SimCore/G4DarkBreM/data/electron_tungsten_MaxE_4.0_MinE_0.2_RelEStep_0.1_UndecayedAP_mA_0.1_run_3000.csv.gz',
        'ldmx-det-v14' , #name of geometry to use
        ),
    ldmxcfg.Analyzer('dbint','dqm::NtuplizeDarkBremInteraction','DQM')
    ]
p.maxTriesPerEvent = 1000
p.maxEvents = 10000
p.outputFiles = ['/dev/null']
p.histogramFile = 'ecal_db.root'

with uproot.open('ecal_db.root') as f :
    ecal_db = f['dbint/dbint'].arrays(library='pd')

plt_hist(
    ecal_db['weight']*(100/0.01)**2, bins='auto',
    xlabel='Event Weight * Biasing Factor',
)

plt_hist(
    ecal_db['z'], weights=ecal_db['weight'], bins = 100,
    xlabel = 'Dark Brem z [mm]',
    ylabel = 'Weighted Event Count'
)

plt_hist(
    ecal_db['x'], weights = ecal_db['weight'], bins = 100,
    xlabel = 'Dark Brem x [mm]',
)

plt_hist(
    ecal_db['y'], weights = ecal_db['weight'], bins = 100,
    xlabel = 'Dark Brem y [mm]',
)

Target Sample

The tungsten target is a "thin target" scenario so we expect a sharply peaked weight distribution and a roughly uniform distribution in the position coordinates.

This sample was generated with

import os
from LDMX.Framework import ldmxcfg
p = ldmxcfg.Process('target_db')
from LDMX.Biasing import target
import LDMX.Ecal.EcalGeometry
import LDMX.Hcal.HcalGeometry
p.sequence = [
    target.dark_brem( 
        100., #MeV - mass of A'
        os.environ['LDMX_BASE']+'/ldmx-sw/SimCore/G4DarkBreM/data/electron_tungsten_MaxE_4.0_MinE_0.2_RelEStep_0.1_UndecayedAP_mA_0.1_run_3000.csv.gz',
        'ldmx-det-v14', #name of geometry to use
        ),
    ldmxcfg.Analyzer('dbint','dqm::NtuplizeDarkBremInteraction','DQM')
    ]
p.maxEvents = 10000
p.maxTriesPerEvent = 10000
p.outputFiles = [ '/dev/null' ]
p.histogramFile = 'target_db.root'

with uproot.open('target_db.root') as f :
    target_db = f['dbint/dbint'].arrays(library='pd')
# the target doesn't always successfully dark brem since it is so small,
# we just need to drop the events that didn't dark brem
# this can be avoided by increasing the maximum number of tries per event
target_db.drop(target_db[(target_db.weight < 1e-100)].index, inplace=True)

plt_hist(
    target_db['weight']*(100/0.01)**2, bins='auto',
    xlabel='Event Weight * Biasing Factor',
)

plt_hist(
    target_db['z'], weights=target_db['weight'], bins = 100,
    xlabel = 'Dark Brem z [mm]',
    ylabel = 'Weighted Event Count'
)

plt_hist(
    target_db['x'], weights = target_db['weight'], bins = 100,
    xlabel = 'Dark Brem x [mm]',
)

plt_hist(
    target_db['y'], weights = target_db['weight'], bins = 100,
    xlabel = 'Dark Brem y [mm]',
)