dqm::DarkBremInteraction Class Reference

Go through the particle map and find the dark brem products, storing their vertex and the dark brem outgoing kinematics for further study. More...

#include <DarkBremInteraction.h>

Public Member Functions
	DarkBremInteraction (const std::string &n, framework::Process &p)
virtual void	produce (framework::Event &e) override
	extract the kinematics of the dark brem interaction from the SimParticles
void	configure (framework::config::Parameters &parameters) override
	Callback for the EventProcessor to configure itself from the given set of parameters.
Public Member Functions inherited from framework::Producer
	Producer (const std::string &name, Process &process)
	Class constructor.
virtual void	process (Event &event) final
	Processing an event for a Producer is calling produce.
Public Member Functions inherited from framework::EventProcessor
	DECLARE_FACTORY (EventProcessor, EventProcessor *, const std::string &, Process &)
	declare that we have a factory for this class
	EventProcessor (const std::string &name, Process &process)
	Class constructor.
virtual	~EventProcessor ()=default
	Class destructor.
virtual void	beforeNewRun (ldmx::RunHeader &run_header)
	Callback for Producers to add parameters to the run header before conditions are initialized.
virtual void	onNewRun (const ldmx::RunHeader &run_header)
	Callback for the EventProcessor to take any necessary action when the run being processed changes.
virtual void	onFileOpen (EventFile &event_file)
	Callback for the EventProcessor to take any necessary action when a new event input ROOT file is opened.
virtual void	onFileClose (EventFile &event_file)
	Callback for the EventProcessor to take any necessary action when a event input ROOT file is closed.
virtual void	onProcessStart ()
	Callback for the EventProcessor to take any necessary action when the processing of events starts, such as creating histograms.
virtual void	onProcessEnd ()
	Callback for the EventProcessor to take any necessary action when the processing of events finishes, such as calculating job-summary quantities.
template<class T >
const T &	getCondition (const std::string &condition_name)
	Access a conditions object for the current event.
TDirectory *	getHistoDirectory ()
	Access/create a directory in the histogram file for this event processor to create histograms and analysis tuples.
void	setStorageHint (framework::StorageControl::Hint hint)
	Mark the current event as having the given storage control hint from this module_.
void	setStorageHint (framework::StorageControl::Hint hint, const std::string &purposeString)
	Mark the current event as having the given storage control hint from this module and the given purpose string.
int	getLogFrequency () const
	Get the current logging frequency from the process.
int	getRunNumber () const
	Get the run number from the process.
std::string	getName () const
	Get the processor name.
void	createHistograms (const std::vector< framework::config::Parameters > &histos)
	Internal function which is used to create histograms passed from the python configuration @parma histos vector of Parameters that configure histograms to create.

Private Attributes
std::map< std::string, int >	known_materials_
	the list of known materials assigning them to material ID numbers
std::map< int, int >	known_elements_
	The list of known elements assigning them to the bins that we are putting them into.
std::string	particle_coll_name_
std::string	particle_passname_

Additional Inherited Members
Protected Member Functions inherited from framework::EventProcessor
void	abortEvent ()
	Abort the event immediately.
Protected Attributes inherited from framework::EventProcessor
HistogramPool	histograms_
	helper object for making and filling histograms
NtupleManager &	ntuple_ {NtupleManager::getInstance()}
	Manager for any ntuples.
logging::logger	the_log_
	The logger for this EventProcessor.

Detailed Description

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

Definition at line 29 of file DarkBremInteraction.h.

Constructor & Destructor Documentation

DarkBremInteraction()

dqm::DarkBremInteraction::DarkBremInteraction ( const std::string & n, framework::Process & p )

inline

Definition at line 31 of file DarkBremInteraction.h.

      : framework::Producer(n, p) {}

Member Function Documentation

configure()

void dqm::DarkBremInteraction::configure ( framework::config::Parameters & parameters )

overridevirtual

Callback for the EventProcessor to configure itself from the given set of parameters.

The parameters a processor has access to are the member variables of the python class in the sequence that has class_name equal to the EventProcessor class name.

For an example, look at MyProcessor.

Parameters

parameters

Parameters for configuration.

Reimplemented from framework::EventProcessor.

Definition at line 8 of file DarkBremInteraction.cxx.

                                                                         {
  particle_coll_name_ = parameters.get<std::string>("particle_coll_name");
  particle_passname_ = parameters.get<std::string>("particle_passname");
}

References framework::config::Parameters::get().

produce()

void dqm::DarkBremInteraction::produce ( framework::Event & e )

overridevirtual

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.

Implements framework::Producer.

Definition at line 41 of file DarkBremInteraction.cxx.

                                                       {
  histograms_.setWeight(event.getEventHeader().getWeight());
  const auto& particle_map{event.getMap<int, ldmx::SimParticle>(
      particle_coll_name_, particle_passname_)};
  const ldmx::SimParticle *recoil{nullptr}, *aprime{nullptr}, *beam{nullptr};
  for (const auto& [track_id, particle] : particle_map) {
    if (track_id == 1) beam = &particle;
    if (particle.getProcessType() ==
        ldmx::SimParticle::ProcessType::eDarkBrem) {
      if (particle.getPdgID() == 622) {
        if (aprime != nullptr) {
          EXCEPTION_RAISE("BadEvent", "Found multiple A' in event.");
        }
        aprime = &particle;
      } else {
        recoil = &particle;
      }
    }
  }
 
  if (recoil == nullptr and aprime == nullptr) {
    /* dark brem did not occur during the simulation
     *    IF PROPERLY CONFIGURED, this occurs because the simulation
     *    exhausted the maximum number of tries to get a dark brem
     *    to occur. We just leave early so that the entries in the
     *    ntuple are the unphysical numeric minimum.
     */
    ldmx_log(error) << " No dark brem occured in this event";
    return;
  }
 
  if (recoil == nullptr or aprime == nullptr or beam == nullptr) {
    // we are going to end processing so let's take our time to
    // construct a nice error message
    ldmx_log(fatal)
        << "Unable to find all necessary particles for DarkBrem interaction."
        << " Missing: [ " << (recoil == nullptr ? " recoil " : "")
        << (aprime == nullptr ? " aprime " : "")
        << (beam == nullptr ? " beam " : "") << "]";
    EXCEPTION_RAISE(
        "BadEvent",
        "Unable to find all necessary particles for DarkBrem interaction.");
    return;
  }
 
  const auto& recoil_p = recoil->getMomentum();
  const auto& aprime_p = aprime->getMomentum();
  ROOT::Math::XYZVector recoil_pvec(recoil_p[0], recoil_p[1], recoil_p[2]);
  ROOT::Math::XYZVector aprime_pvec(aprime_p[0], aprime_p[1], aprime_p[2]);
 
  std::vector<double> incident_p = recoil_p;
  for (std::size_t i{0}; i < recoil_p.size(); ++i)
    incident_p[i] += aprime_p.at(i);
 
  double incident_energy = energy(incident_p, recoil->getMass());
  double recoil_energy = energy(recoil_p, recoil->getMass());
 
  std::vector<double> ap_vertex{aprime->getVertex()};
  std::string ap_vertex_volume{aprime->getVertexVolume()};
  auto ap_vertex_material_it = std::find_if(
      known_materials_.begin(), known_materials_.end(),
      [&](const auto& mat_pair) {
        return ap_vertex_volume.find(mat_pair.first) != std::string::npos;
      });
  int ap_vertex_material = (ap_vertex_material_it != known_materials_.end())
                               ? ap_vertex_material_it->second
                               : 0;
 
  if (ap_vertex_material == 0) {
    ldmx_log(warn) << "Dark brem interaction occurred in an unknown material: "
                   << ap_vertex_volume;
  }
 
  int ap_parent_id{-1};
  if (aprime->getParents().size() > 0) {
    ap_parent_id = aprime->getParents().at(0);
  } else {
    ldmx_log(error) << "Found A' without a parent ID!";
  }
 
  float aprime_energy = energy(aprime_p, aprime->getMass());
  int aprime_genstatus = aprime->getGenStatus();
  double aprime_px{aprime_p.at(0)}, aprime_py{aprime_p.at(1)},
      aprime_pz{aprime_p.at(2)};
  event.add("APrimeEnergy", aprime_energy);
  event.add("APrimePx", aprime_px);
  event.add("APrimePy", aprime_py);
  event.add("APrimePz", aprime_pz);
  event.add("APrimeParentID", ap_parent_id);
  event.add("APrimeGenStatus", aprime_genstatus);
 
  histograms_.fill("aprime_energy", aprime_energy);
  histograms_.fill("aprime_pt", quadsum({aprime_px, aprime_py}));
  histograms_.fill("aprime_theta", aprime_pvec.Theta() * (180 / 3.14159));
 
  int recoil_genstatus = recoil->getGenStatus();
  double recoil_px{recoil_p.at(0)}, recoil_py{recoil_p.at(1)},
      recoil_pz{recoil_p.at(2)};
  event.add("RecoilEnergy", recoil_energy);
  event.add("RecoilPx", recoil_px);
  event.add("RecoilPy", recoil_py);
  event.add("RecoilPz", recoil_pz);
  event.add("RecoilGenStatus", recoil_genstatus);
 
  histograms_.fill("recoil_energy", recoil_energy);
  histograms_.fill("recoil_pt", quadsum({recoil_px, recoil_py}));
  histograms_.fill("recoil_theta", recoil_pvec.Theta() * (180 / 3.14159));
 
  event.add("IncidentEnergy", incident_energy);
  double incident_px{incident_p.at(0)}, incident_py{incident_p.at(1)},
      incident_pz{incident_p.at(2)};
  event.add("IncidentPx", incident_px);
  event.add("IncidentPy", incident_py);
  event.add("IncidentPz", incident_pz);
 
  histograms_.fill("incident_energy", incident_energy);
  histograms_.fill("incident_pt", quadsum({incident_px, incident_py}));
 
  double vtx_x{aprime->getVertex().at(0)}, vtx_y{aprime->getVertex().at(1)},
      vtx_z{aprime->getVertex().at(2)};
  event.add("DarkBremX", vtx_x);
  event.add("DarkBremY", vtx_y);
  event.add("DarkBremZ", vtx_z);
  event.add("DarkBremVertexMaterial", ap_vertex_material);
  float db_material_z =
      event.getEventHeader().getFloatParameter("db_material_z");
  event.add("DarkBremVertexMaterialZ", db_material_z);
  float aprime_conversion_material_z =
      event.getEventHeader().getFloatParameter("aprime_conversion_material_z");
 
  histograms_.fill("dark_brem_z", vtx_z);
 
  int i_element = 0;
  if (db_material_z > 0) {
    if (known_elements_.find(static_cast<int>(db_material_z)) ==
        known_elements_.end()) {
      i_element = known_elements_.size();
      ldmx_log(warn)
          << "Dark brem interaction occurred in an unknown element with Z = "
          << db_material_z << ". Using index " << i_element
          << " for this element.";
    } else {
      i_element = known_elements_.at(static_cast<int>(db_material_z));
    }
  }
 
  histograms_.fill("dark_brem_element", i_element + 0.5);
  histograms_.fill("dark_brem_material", ap_vertex_material + 0.5);
 
  // Get the daughters of the A' if it decayed within the simulation
  std::vector<int> aprime_daughters = aprime->getDaughters();
  int n_ap_daughters = aprime_daughters.size();
  ldmx_log(debug) << "A' with energy " << aprime->getEnergy() << " and momentum"
                  << " (" << aprime_px << ", " << aprime_py << ", " << aprime_pz
                  << ") GeV " << " has " << n_ap_daughters << " daughters";
  if (n_ap_daughters == 0) {
    histograms_.fill("aprime_daughter_pdgid", 0);
    histograms_.fill("aprime_daughter_material", 0.5);
    histograms_.fill("aprime_daughter_element", 0.5);
  } else {
    // Loop again on the particles to find the daughters of the A'
    for (const auto& [track_id, daughter_particle] : particle_map) {
      for (const auto& primary_daughter : aprime_daughters) {
        if (track_id == primary_daughter) {
          auto const& daughter_p = daughter_particle.getMomentum();
          double daughter_px{daughter_p.at(0)}, daughter_py{daughter_p.at(1)},
              daughter_pz{daughter_p.at(2)};
 
          ldmx_log(debug) << "  Daughter track ID " << track_id
                          << " with PDG ID " << daughter_particle.getPdgID()
                          << " and energy " << daughter_particle.getEnergy()
                          << " and charge " << daughter_particle.getCharge()
                          << " and mass " << daughter_particle.getMass()
                          << " GeV" << " and momentum (" << daughter_px << ", "
                          << daughter_py << ", " << daughter_pz << ") GeV";
          histograms_.fill("aprime_daughter_energy",
                           daughter_particle.getEnergy());
          histograms_.fill("aprime_daughter_pt",
                           quadsum({daughter_px, daughter_py}));
 
          // Fill histogram for daughter creation vertex Z position
          double daughter_start_z = daughter_particle.getVertex().at(2);
          histograms_.fill("aprime_daughter_start_z", daughter_start_z);
 
          // Fill histogram for material where A' daughter was created.
          // Prefer the explicit interaction material; fall back to vertex
          // volume.
          std::string daughter_material_name =
              daughter_particle.getInteractionMaterial();
          std::string daughter_vertex_volume =
              daughter_particle.getVertexVolume();
          int daughter_material = 0;
 
          if (daughter_material_name.find("Carbon") != std::string::npos) {
            daughter_material = 1;
          } else if (daughter_material_name.find("FR4") != std::string::npos ||
                     daughter_material_name.find("PCB") != std::string::npos ||
                     daughter_vertex_volume.find("motherboard") !=
                         std::string::npos ||
                     daughter_vertex_volume.find("PCB") != std::string::npos) {
            daughter_material = 2;
          } else if (daughter_material_name.find("Glue") != std::string::npos ||
                     daughter_vertex_volume.find("Glue") != std::string::npos ||
                     daughter_vertex_volume.find("CFMix") !=
                         std::string::npos) {
            daughter_material = 3;
          } else if (daughter_material_name.find("Silicon") !=
                         std::string::npos ||
                     daughter_material_name.find("Si") != std::string::npos ||
                     daughter_vertex_volume.find("Si") != std::string::npos ||
                     daughter_vertex_volume.find("Sensor") !=
                         std::string::npos ||
                     daughter_vertex_volume.find("sensor") !=
                         std::string::npos) {
            daughter_material = 4;
          } else if (daughter_material_name.find("Al") != std::string::npos ||
                     daughter_material_name.find("Aluminum") !=
                         std::string::npos ||
                     daughter_vertex_volume.find("strongback") !=
                         std::string::npos ||
                     daughter_vertex_volume.find("support") !=
                         std::string::npos) {
            daughter_material = 5;
          } else if (daughter_material_name.find("W") != std::string::npos ||
                     daughter_material_name.find("Tungsten") !=
                         std::string::npos ||
                     daughter_vertex_volume.find("target") !=
                         std::string::npos ||
                     daughter_vertex_volume.find("W_front_volume") !=
                         std::string::npos ||
                     daughter_vertex_volume.find("W_cooling") !=
                         std::string::npos) {
            daughter_material = 6;
          } else if (daughter_material_name.find("Polyvinyltoluene") !=
                         std::string::npos ||
                     daughter_material_name.find("PVT") != std::string::npos ||
                     daughter_vertex_volume.find("trigger_pad") !=
                         std::string::npos) {
            daughter_material = 7;
          } else if (daughter_material_name.find("Air") != std::string::npos ||
                     daughter_vertex_volume.find("Air") != std::string::npos) {
            daughter_material = 8;
          } else {
            ldmx_log(warn) << "Daughter particle track ID " << track_id
                           << " created in unknown material: "
                           << daughter_material_name
                           << " and vertex volume: " << daughter_vertex_volume;
          }
          histograms_.fill("aprime_daughter_material", daughter_material + 0.5);
 
          // Fill histogram for element where A' conversion happened.
          // This comes from the conversion process selecting an element in
          // material.
          int daughter_element = 0;
          if (aprime_conversion_material_z > 0) {
            if (known_elements_.find(static_cast<int>(
                    aprime_conversion_material_z)) == known_elements_.end()) {
              daughter_element = known_elements_.size();
            } else {
              daughter_element = known_elements_.at(
                  static_cast<int>(aprime_conversion_material_z));
            }
          }
          histograms_.fill("aprime_daughter_element", daughter_element + 0.5);
 
          if (daughter_particle.getPdgID() == 11) {
            histograms_.fill("aprime_daughter_pdgid", 1.5);
          } else if (daughter_particle.getPdgID() == -11) {
            histograms_.fill("aprime_daughter_pdgid", 2.5);
          } else if (daughter_particle.getPdgID() == 13) {
            histograms_.fill("aprime_daughter_pdgid", 3.5);
          } else if (daughter_particle.getPdgID() == -13) {
            histograms_.fill("aprime_daughter_pdgid", 4.5);
          } else if (daughter_particle.getPdgID() == 17) {
            histograms_.fill("aprime_daughter_pdgid", 5.5);
          } else if (daughter_particle.getPdgID() == -17) {
            histograms_.fill("aprime_daughter_pdgid", 6.5);
          } else if (daughter_particle.getPdgID() == 211) {
            histograms_.fill("aprime_daughter_pdgid", 7.5);
          } else if (daughter_particle.getPdgID() == -211) {
            histograms_.fill("aprime_daughter_pdgid", 8.5);
          } else {
            histograms_.fill("aprime_daughter_pdgid", 9.5);
          }
        }  // end if track_id matches primary daughter
      }  // end loop over A' daughters
    }  // end loop over particles
  }  // end if n_ap_daughters > 0
 
  // Get recoil electron daughters if it underwent bremsstrahlung
  std::vector<int> recoil_daughters = recoil->getDaughters();
  int n_recoil_brem_daughters = 0;
  // Loop again on the particles to find the daughters of the recoil electron
  for (const auto& [track_id, daughter_particle] : particle_map) {
    for (const auto& primary_daughter : recoil_daughters) {
      if (track_id == primary_daughter) {
        if (daughter_particle.getEnergy() > (0.2 * recoil->getEnergy()) &&
            (daughter_particle.getPdgID() == 22)) {
          n_recoil_brem_daughters++;
          histograms_.fill("recoil_brem_daughter_energy",
                           daughter_particle.getEnergy());
          histograms_.fill("recoil_brem_daughter_energy_ratio",
                           daughter_particle.getEnergy() / recoil->getEnergy());
          ldmx_log(debug) << "  Recoil electron daughter track ID " << track_id
                          << " with PDG ID " << daughter_particle.getPdgID()
                          << " and energy " << daughter_particle.getEnergy();
        }
      }
    }  // end loop over recoil daughters
  }  // end loop over particles
  histograms_.fill("recoil_brem_daughter_num", n_recoil_brem_daughters);
}  // end of produce

References framework::HistogramPool::fill(), framework::Event::getEventHeader(), ldmx::SimParticle::getMomentum(), ldmx::EventHeader::getWeight(), framework::EventProcessor::histograms_, known_elements_, known_materials_, and framework::HistogramPool::setWeight().

Member Data Documentation

known_elements_

std::map<int, int> dqm::DarkBremInteraction::known_elements_

private

Initial value:

= {{1, 1},  {6, 2},  {8, 3},  {11, 4},
                                        {14, 5}, {20, 6}, {29, 7}, {39, 8},
                                        {71, 9}, {74, 10}}

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 : 'Y 39', 9 : 'Lu 71', 10 : 'W 74', 11 : 'unlisted' }

Definition at line 130 of file DarkBremInteraction.h.

                                     {{1, 1},  {6, 2},  {8, 3},  {11, 4},
                                        {14, 5}, {20, 6}, {29, 7}, {39, 8},
                                        {71, 9}, {74, 10}};

Referenced by produce().

known_materials_

std::map<std::string, int> dqm::DarkBremInteraction::known_materials_

private

Initial value:

= {
      {"Carbon", 1},
      {"PCB", 2},  
                   
      {"Glue", 3},
      {"Si", 4},
      {"Al", 5},
      {"W", 6},
      {"target", 6},
      {"trigger_pad", 7},
      {"strongback", 5},   
                           
                           
      {"motherboard", 2},  
      {"support", 5},      
      {"CFMix", 3},        
      {"C_volume", 1},     
                           
      {"Air", 8}}

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', 8 : 'Air' }

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.

Definition at line 83 of file DarkBremInteraction.h.

                                            {
      {"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
      {"Air", 8}};

Referenced by produce().

particle_coll_name_

std::string dqm::DarkBremInteraction::particle_coll_name_

private

Definition at line 133 of file DarkBremInteraction.h.

particle_passname_

std::string dqm::DarkBremInteraction::particle_passname_

private

Definition at line 134 of file DarkBremInteraction.h.

---

The documentation for this class was generated from the following files:

• DQM/include/DQM/DarkBremInteraction.h
• DQM/src/DQM/DarkBremInteraction.cxx