ldmx-software/simcore

Integration of Geant4 simulation framework into a more generalized event processing framework.

simcore's Issues

Add an option to check if a track is a descendant of another track in the trackmap

See LDMX-Software/ldmx-sw#1273 :)

See if kaon properties can be configured at runtime

One of the changes in https://github.com/LDMX-Software/geant4/tree/LDMX.upKaons_mod which would be useful to be able to make use of even without the kaon enhancement features of that branch is the changes to the decay modes and lifetimes of the charged kaons (Geant4/geant4@25228b8). If it is straight-forward to implement this kind of customization as a part of LDMX-sw rather than direct changes to Geant4, this could make it easier to do other studies with challenging kaon events.

I think this actually should be quite straight-forward to do, you should be able to access and configuring the decay table when constructing the physics.

Dark Brem Process Enhancements

In order to accomodate the ECal as Target (EaT) project, several changes need to be made to the dark brem process. I am in the process of validating these changes for the EaT project and then I can test to make sure that the dark brem process can be used to generate signal in the target case.

Still TODO

• Update documentation note on DocDB
• Label used/not-used vertices for the Forward Only running mode
• Write functionality test : #3
• Include a dark brem vertex library that won't work for production but will help people get up and running for small samples

Update to make use of Modern CMake

Organize primary generators into subdirectory

This will help clean up the module a little more. This issue will also give me an opportunity to have the primary generators be passed the run header and add their configuration themselves.

Transition to using external G4DarkBreM

Our dark brem simulation method has been moved to its own repository https://github.com/tomeichlersmith/G4DarkBreM in order to be shareable with the outside world. We should update our simulation to pull in this module as a submodule and use it.

An example use case very close to how we use it is https://github.com/tomeichlersmith/ldmx-sim-technique

Update for move of Exception to Framework

Visualize Detector using Geant4

A lot of people have been wanting to be able to visualize the detector like we used to, using Geant4, so that we can verify that the gdml is correctly defining the geometry. I think we can work this into our current method using two steps:

1. Define a C++ program that links with Geant4 and its visualization libraries. This would take a path to the gdml as an input. Potential issues:

   • Are the Geant4 visualization libraries built in the container installation? I think so, it looks like they are built and already required.
   • Can we open the Geant4 visualization from the container with a specific mounting procedure? Probably, although this will change depending on what the host OS is.
   • How independent can this program be? Maybe we can make it a configuration (cmake) option so people can skip it if it is requiring problems.

2. Define a bash alias in the ldmx-env.sh script that will do the following. These changes will need to be in the central ldmx-sw repository.

   • Go into the container with the correct mounting procedure
   • Run the visualization command with the input gdml

@Hongyin-debug I am "assigning" this to you because we talked about it and I think you are working on it. If you stop working on it, please un-assign yourself. If you are working on it, please put helpful notes/links here so we can work together on sharing this with everyone.

Vertex Smearing is not being Seeded Correctly

@bryngemark found this while comparing inclusive overlay files and the sim hits in the tagger trigger pad. The tagger trigger pad is really close to the primary vertex, so that is why we were able to observe this bug there.

The primary generator action has its own random number generator that it seeds when it is constructed:

SimCore/src/SimCore/PrimaryGeneratorAction.cxx

Lines 38 to 39 in 0a05e1d

	random_ = std::make_unique<TRandom3>();
	random_->SetSeed(CLHEP::HepRandom::getTheSeed());

The run manager creates the primary generator action in RunManager::Initialize:

SimCore/src/SimCore/RunManager.cxx

Lines 124 to 125 in 0a05e1d

	auto primaryGeneratorAction{new PrimaryGeneratorAction(parameters_)};
	SetUserAction(primaryGeneratorAction);

which is called in Simulator::onProcessStart:

SimCore/src/SimCore/Simulator.cxx

Lines 327 to 328 in 0a05e1d

	// initialize run
	runManager_->Initialize();

But the random seed is not set until Simulator::onNewRun:

SimCore/src/SimCore/Simulator.cxx

Lines 283 to 291 in 0a05e1d

	void Simulator::onNewRun(const ldmx::RunHeader&) {
	const framework::RandomNumberSeedService& rseed =
	getCondition<framework::RandomNumberSeedService>(
	framework::RandomNumberSeedService::CONDITIONS_OBJECT_NAME);
	std::vector<int> seeds;
	seeds.push_back(rseed.getSeed("Simulator[0]"));
	seeds.push_back(rseed.getSeed("Simulator[1]"));
	setSeeds(seeds);
	}

So I think that the primary generator action (and only the primary generator action) is not being seeded correctly. This is a bug that should be "easy" to fix because we just need to move the seeding of the PGA to after the seeding of Geant4 (since the PGA uses G4's CLHEP seed). The reason we haven't seen this before is because this only affects primary vertex smearing between runs.

A simple solution is to use Geant4's internal random generator (which uses the seed after it is set in onNewRun). Basically, we would remove the random_ member variable of the Primary Generator Action, and instead use G4UniformRand() (which still needs to be scaled). i.e. replace

SimCore/src/SimCore/PrimaryGeneratorAction.cxx

Lines 140 to 142 in 0a05e1d

	double x0_f = random_->Uniform(x0_i - IPWidthX, x0_i + IPWidthX);
	double y0_f = random_->Uniform(y0_i - IPWidthY, y0_i + IPWidthY);
	double z0_f = random_->Uniform(z0_i - IPWidthZ, z0_i + IPWidthZ);

with

        double x0_f = G4UniformRand()*2*IPWidthX - IPWidthX + x0_i;
        double y0_f = G4UniformRand()*2*IPWidthY - IPWidthY + y0_i;
        double z0_f = G4UniformRand()*2*IPWidthZ - IPWidthZ + z0_i;

Clean-Up Sensitive Detector Implementation

There are a lot of new and delete floating around when we use Geant4's sensitive detectors, hit collections, and hits. This includes custom memory allocation!! To me, the process of going "through" G4-style hits and hit collections just to have them translated to our hits and hit collections is a waste of time and space. We already control the senstive detector implementation, so we could skip the G4-style nonsense (and custom memory allocation) and go straight to our STL collections which can be added directly to the event bus.

Downside: This would effectively remove the possibility of running this repository independent of the Framework repository. I think that makes sense - we are specializing Geant4 to our software engine, but I don't know how much an independent G4 simulation is desired.

Resim from ECal scoring plane is being done incorrectly.

Currently, the resim code loops through all particles in the ECal scoring plane collection and re-fires them. Looking at the collection more closely, it contains multiple hits from the same particle and also has the incident electron and photon. The code needs to be updated to check for duplicates and remove and incident particles.

Add Z-Axis to Dark Brem Reference Library

The current dark brem simulation uses a reference library of MadGraph-generated events to simulate the outgoing kinematics after Geant4 determines that a dark brem should happen within a given material. Currently, this reference library has only a single axis - the incident energy of the electron. This has proven successful for simulating a wide range of incident electron energies due to various upstream interactions occuring. Adding another axis to the library corresponding to the atomic number (Z) of the nucleus that the dark brem was generated within MG with would allow this implementation of dark brem to happen realistically in a material that has a wider variety of atoms in it (most importantly, a potential active target made of LYSO and materials within the ECal like PCB).

Expected Changes

• Update LHE reader to extract the atomic number alongside the incident energy, recoil pT and recoil KE fraction
• Introduce a method for random selection of atom when a composite material is chosen to dark brem
  • probably relying on xsec-weighted random selection
• Choose outgoing kinematics based on nearest incident energy and nearest Z
  • materials close-in-Z pretty realistically model each other according to MG, so we don't need to have a library entry for each unique atom

RunNumber cleanup

RunNumber cleanup, see LDMX-Software/ldmx-sw#1045

By running a quick ag runNumber, the following cases in SimCore show up

• test/basic.py
• python/examples.py.in

Increase the number of processes in the process map

By default, we are only including entries that are interesting for electrons and photons. This leads to most entries in the particle tables to be marked as "unkown" which has caused confusion for some users. When working on #106, I needed to be able to check that the daughters of a primary kaon were from its decay or not so I needed an entry for at least the Geant4 "Decay" process.

Adding more entries to the process map will increase the runtime for SimParticle::findProcessType which will have to search a larger map. I would argue that if we are running into performance issues because of this, we should look at the datastructure and algorithm we are using (map + string comparison) first.

We should think through what particle + process combinations we are especially interested in, off the top of my head

• Kaons
• Neutrons

Dark Brem Process Testing

Not necessarily testing that it is physically realistic, more that it works the way we expect. Not even going to put it into a simulation run, just testing the public members of the class that they work in the way they are used within Geant4 and our simulation.

• Import LHE vertices
• Scale vertices depending on configuration
• Produce A' and electron secondaries from original electron

Switch to using Hcal Geometry condition for scintillator orientations

See LDMX-Software/ldmx-sw#1135 and LDMX-Software/ldmx-sw#1117 for details. We have added more detailed information into the HCal geometry condition, so the need to use layerIsHorizontal which only was useful for the back HCal isn't needed any more other than when decoding older geometry copy-numbers. At the same time, for the side hcal we can now handle 3D readout the same way which hopefully can simplify some things in here.

Trajectory Handling

Currently, the track -> trajectory -> SimParticle handling in the simulation is very hap-hazard and complicated. I think we would benefit from cleaning it up and we can look into removing the trajectory step all together. Moreover, I have a feeling that the intricacies in how the tracking action decides which particles are saved are the cause for some of the issues @cmantill and @AminaJLi are seeing with kaon-less kaon events.

Hcal SimHit updates

We need to record some additional information in the SimHits for LDMX-Software/Hcal#12 by @ehrlich-uva. In particular Ralf needs

• pre- and post-step positions and times so that I can calculate the path through the counter.
• The true path length is useful too, because the true path length may be longer than the distance between both points.
  Can be obtained from G4Step::GetStepLength()
• I also need the charge of the particle
  Already accessible from the PDG code
• The visible deposited energy
  In the Hcal, we already apply Birk's law to the deposited energy. We should discuss whether or not this implementation is what we want to use in the end, but for now it should be fine.
• The speed/beta of the particle (or some quantities that allow me to calculate it) for the Cerenkov light
  Can we obtained either from the track or we could take the average between the pre/post step points

For SimHits with just one contrib, like what we do in Hcal, this would be straight-forward to add to just record as part of the SimCalorimeterHit. I don't think recording all of this for individual contributions would make sense for

• Path length
• Pre/post step t,x,y,z
• Velocity
  we will need to decide on what to do within the EcalSD. They could be left defaulted (either 0 or some signal value), we could record the values for one of the contribs (if so, which to pick? Time-wise earliest?). Latter would bloat the add/update contrib interface even further. @tomeichlersmith, do you have any thoughts here? We discussed this briefly in Slack, but it would be better to have things recorded here.

I'm not all that familiar with what happens underneath the hood with ROOT, does this revision require us to bump the ClassDef number?

@ehrlich-uva what units do you want these quantities in?

(Minor bug) wrong unit in simulation's energy printout

(Minor bug so I just skipped the formality)

We set the particle gun's energy using GeV in the config file. However, the printout when we fire the config appears to be in MeV:

RunHeader { run: 1, detectorName: ldmx-det-v12, description: Single electron gun
  intParameters: 
    Compress ECal Hit Contribs = 1
    Included Scoring Planes = 1
    RandomNumberMasterSeed[v12] = 1
    Save ECal Hit Contribs = 1
    Use Random Seed from Event Header = 0
  floatParameters: 
    Gen 1 Direction X = 0
    Gen 1 Direction Y = 0
    Gen 1 Direction Z = 1
    Gen 1 Energy [MeV] = 4
    Gen 1 Position [mm] X = 0
    Gen 1 Position [mm] Y = 0
    Gen 1 Position [mm] Z = -1.2
    Gen 1 Time [ns] = 0
    Smear Beam Spot [mm] X = 20
    Smear Beam Spot [mm] Y = 80
    Smear Beam Spot [mm] Z = 0
  stringParameters: 
    Geant4 revision = a50bbe0a781fa634a2811f70ed6ec56cc89e101f
    Gen 1 Class = ldmx::ParticleGun
    Gen 1 Particle = e-
    ldmx-sw revision = 1c87e34ce12d06b65284e0d0834f5e603799676f
}[ Simulator ] : Started 100 events to produce 100 events.

Here is the code that prints the energy:

SimCore/src/Simulator.cxx

Line 193 in 908687e

header.setFloatParameter( genID + " Energy [MeV]" , gen.getParameter<double>("energy") );

Load biasing operators dynamically (like UserActions, PrimaryGenerators, etc...)

The base class XsecBiasingOperator already does most of the work interfacing with Geant4, so we can have our derived classes instead be modified to be "biasing calculators" or something similar that are dynamically loaded into the XsecBiasingOperator and then specific functions can be called.

Add a field to event info for tracking whether the primary particle entered the tagger

See comments under LDMX-Software/ldmx-sw#1146

Update LHE reader to enable placing A' decay vertex at a specific position

This is required when looking at long lived visibly decaying dark photons.

Use HcalGeometry condition for layer parity

See LDMX-Software/ldmx-sw#1063 for details. Add the HcalGeometry condition to HcalSD to be able to get the horizontal parity for different layers since this no longer is the same for all geometries.

Allow saving of auxiliary photonuclear info

Specifically, we want to tag and save the photonuclear secondaries directly associated with the photon to their own collection. Whether this collection is created should be configurable.

Allow for rotation of PrimaryGeneratorVertices

Alternative PN models

As discussed in the Hcal meeting yesterday, there are cases where it would be useful to be able to exchange the model used for PN reactions or to extend it.

Examples:

1. Replacing PN products with known complicated PN topologies for comparing different detector geometries
2. Replacing Bertini with alternative models (e.g. what I've done with FLUKA)
3. Instrumenting a model to extend it in some way, e.g. to run the event generator until you produce an event topology that you are interested in, but could in principle be used to do anything you need (e.g. debugging)

This is relatively straight-forward to do and would mostly end up concerning the GammaPhysics code, which would get which PN-model to use as a parameter. By default, do what we already do (don't touch the model but move its place in the process list) so no regular uses would be affected. If the model parameter is something else, replace the photonNuclear process with a different one and place it in the correct location. I have examples of how to do both 2 and 3 above that hopefully should be straight-forward to use for anyone wishing to do something similar.

Currently, I'm doing this directly inside of GammaPhysics but maybe it would be good to be able to load such models dynamically. For example, any model based on using FLUKA might need to be in a separate/private repository. I haven't played around with how to do this nicely yet, but I'm guessing the existing plugin-style code would be a good place for inspiration.

This doesn't really let you do the same thing as #49 but might be an alternative solution for some cases. For example, to look at the difference between "Nothing hard" category events with FLUKA, I've been running a version of the FLUKA hadronics PN model where we re-run the event generator until we produce such an event. This means that we could get away with storing/simulating only 1e5-1e7 events and is much faster than running a full event each time.

LHEReader doesn't work with all madgraph outputs

I am analyzing some LHE files I got from our theory colleagues. In them there is auxiliary event-level information after each particle is listed. An example can be found here:

 7      1 +2.9536814e+04 1.71289800e+03 7.81860800e-03 8.21987300e-02
       11 -1    0    0    0    0 +0.0000000000e+00 +0.0000000000e+00 +3.9999999674e+01 4.0000000000e+01 5.1099890000e-03 0.0000e+00 -1.0000e+00
      623 -1    0    0    0    0 -0.0000000000e+00 -0.0000000000e+00 -5.9685589804e-12 1.7124600000e+03 1.7124600000e+03 0.0000e+00 1.0000e+00
       11  1    1    2    0    0 +1.9043651668e+00 -6.9159080595e-01 +8.7837770258e+00 9.0144145391e+00 5.1099890000e-03 0.0000e+00 -1.0000e+00
      623  1    1    2    0    0 +1.1962853951e-03 +3.8271048071e-01 +3.4922114418e-01 1.7124600784e+03 1.7124600000e+03 0.0000e+00 1.0000e+00
      622  2    1    2    0    0 -1.9055614522e+00 +3.0888032524e-01 +3.0867001504e+01 3.0985507087e+01 1.8982350000e+00 0.0000e+00 0.0000e+00
      -11  1    5    5    0    0 -2.8644920854e-01 -6.8086220641e-01 +5.8139920221e+00 5.8607299749e+00 5.1099890000e-03 0.0000e+00 -1.0000e+00
       11  1    5    5    0    0 -1.6191122511e+00 +9.8974253286e-01 +2.5053009602e+01 2.5124777233e+01 5.1099890000e-03 0.0000e+00 1.0000e+00
<mgrwt>
<rscale>  0 0.17128984E+04</rscale>
<asrwt>0</asrwt>
<pdfrwt beam="1">  1       11 0.10000000E+01 0.17128984E+04</pdfrwt>
<pdfrwt beam="2">  1      623 0.10000000E+01 0.17128984E+04</pdfrwt>
<totfact> 0.10000000E+01</totfact>
</mgrwt>
</event>

This breaks the LHEReader because the code expects that every line in between and (except the first?) is particle data. However, a simple fix can bypass this information and keep the code compatible with the usual LHE format we have: just change this to:

if (line == "</event>" or line=="<mgrwt>") {

This will short circuit the event block analysis and totally ignore weights. If these are useful, I guess someone will need to implement a class for interpreting and persisting this data, but for now, it's not useful to us, as far as I can tell.

Update namespace to simcore to reflect the module name

Simulation Functionalities

• Create and Configure UserActions
• Create and Configure PrimaryGenerators
• Run with each individual generator
• Run with mutliple generators WON'T DO
• Filtering events using event fileters (Biasing Module)
• Saving tracks from certain processes (Biasing Module)
• Turn scoring planes on and off
• Biasing operations can be configured and effect the simulation

Sensitive detector update doesn't register prototype volumes

The for check for isSensDet in

SimCore/include/SimCore/SDs/HcalSD.h

Line 41 in 7d101f9

bool isSensDet(G4LogicalVolume* volume) const final override {

depends on the name of the volume in the GDML sources containing ScintBox. Since that was spelled scint_box in the prototype GDML files, well, they aren't registered. This could be fixed either by retroactively renaming the variables in the prototype geometry or by changing the check to work with either. I'm not sure which would be best

8GeV generator not centering the beam

The 8GeV electron generator is not centering the beam on the target and now that I mention it, the 4GeV beam isn't doing that either although they are both within ~5mm of the center in the horizontal direction and seemingly perfectly centered vertically.

These distributions were made using the 4gev and 8gev particle gun generators currently in SimCore and then a filter on the total energy of the primary electron as it reaches the ECal was applied. This was the sample I had readily available for my work with EaT, it'll probably be a good idea to do an unfiltered sample to make sure my selection isn't biasing how the electron reaches the target.

Add parent ID to Calorimeter hit contributions.

Is your feature request related to a problem? Please describe.
Currently, there is no easy way to associate the MC particle contribution to an ECal sim hit to its parent particle. This makes it difficult to associate showers to different incident particles.

Describe the solution you'd like
Add parentID to the contributions.

Clean up biasing framework

We can clean up the biasing framework so that each operator is python-configured with its own parameters (just like EventProcessors, UserActions, and PrimaryGenerators). This might open the door to biasing more than one process at once, but no promises.

Re-Simulation Processor

We already have two generators that require an input file formatted the same as our event files and these generators attempt to restore the configuration of the random number generator in Geant4 at the time of the first simulation.

I've been trying to use these generators to do some re-sim and it isn't working.

I think this is because the Simulator processor is designed to be the first Producer in a production chain. I think we can get around this issue by writing another processor that is specifically focused on re-simulation. This processor would be designed to provide the current event bus to other simulation parts before starting the simulation, so that Geant4 can be reconfigured the same. Then we can run the ReSimulator processor as the first Producer in an analysis chain.

This will probably require some modifications to other parts of our simulation infrastructure. Another potential solution is to try to have the Simulator work for both running modes; however, I ran into issues trying a few "easy" fixes, so I think separating them will be easier.

Add support for loading of ILC LCDD detectors

The Heavy Photon Search experiment makes use of the LCDD detector description language to describe its detectors. For HPS to use the Simulator, support for LCDD is needed.

Addition of more detailed information to RunHeader during simulation

Currently, the RunHeader has very sparse information: The name of the geometry and a hash from Geant. I think we should look in to adding more specific information automatically. For example, the names of all of the plugins loaded into Geant.

I'm not very familiar with the simulation code; however, so I'm wondering how feasible this is. Could we access some of the settings in the macro? If we can, what should we save to the RunHeader?

My opinion is that we should save as many details as possible; then users can avoid re-simulating any runs that they forgot the settings for.

Patch persistency manager for no-SD case

In a simpler running mode where we want to validate our simulation, we may be running using a detector construction that has no sensitive detectors. This means the list of Hit Collections in an event are never initialized, so the following lines produce a seg-fault.

SimCore/src/SimCore/Persist/RootPersistencyManager.cxx

Lines 126 to 127 in e5e154c

	G4HCofThisEvent *hce = anEvent->GetHCofThisEvent();
	int nColl = hce->GetNumberOfCollections();

This is a simple patch, we just need to check if there is a set of hit collections before trying to see how many there are.