Skip to content
This repository was archived by the owner on May 6, 2024. It is now read-only.

Commit 2e40938

Browse files
EinarElenEinarElen
authored
Merge pull request #88 from LDMX-Software/iss86PNModels
Iss86 pn models
2 parents 9a28a9e + f0a9b64 commit 2e40938

23 files changed

Lines changed: 914 additions & 42 deletions

‎CMakeLists.txt‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -40,6 +40,10 @@ setup_geant4_target()
4040
# include dark brem simulation
4141
add_subdirectory(G4DarkBreM)
4242

43+
setup_library(module SimCore
44+
name PhotoNuclearModels
45+
dependencies SimCore::SimCore)
46+
4347
# Set the LHE Reading Library
4448
setup_library(module SimCore
4549
name LHE

‎include/SimCore/Event/SimParticle.h‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -63,6 +63,12 @@ class SimParticle {
6363
*/
6464
double getEnergy() const { return energy_; }
6565

66+
/**
67+
* Get the kinetic energy of this particle [MeV].
68+
*
69+
* @return The kinetic energy of this particle.
70+
*/
71+
double getKineticEnergy() const { return energy_ - mass_; }
6672
/**
6773
* Get the PDG ID of this particle.
6874
*

‎include/SimCore/GammaPhysics.h‎

Lines changed: 30 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -15,14 +15,17 @@
1515
#include "G4ProcessManager.hh"
1616
#include "G4VPhysicsConstructor.hh"
1717
#include "G4VProcess.hh"
18-
18+
#include "SimCore/PhotoNuclearModel.h"
1919
namespace simcore {
2020

2121
/**
22-
* @class GammaPhysics
23-
* @brief Adds extra gamma particle physics for simulation
22+
* @class GammaPhysics @brief Adds extra gamma particle physics for simulation
23+
* and sets up the photonuclear model to use from the configuration
2424
*
2525
* @note
26+
*
27+
* Is responsible for selecting the photonuclear model from the python
28+
* configuration.
2629
* Currently adds gamma -> mumu reaction using the
2730
* <i>G4GammaConversionToMuons</i> process. Also changes ordering of
2831
* gamma processes such that photonNuclear and GammaToMuMu are called first.
@@ -32,7 +35,13 @@ class GammaPhysics : public G4VPhysicsConstructor {
3235
/**
3336
* Class constructor.
3437
* @param name The name of the physics.
38+
* @param parameters The python configuration
3539
*/
40+
GammaPhysics(const G4String& name,
41+
const framework::config::Parameters& parameters)
42+
: G4VPhysicsConstructor(name),
43+
modelParameters{parameters.getParameter<framework::config::Parameters>(
44+
"photonuclear_model")} {}
3645
GammaPhysics(const G4String& name = "GammaPhysics");
3746

3847
/**
@@ -50,11 +59,29 @@ class GammaPhysics : public G4VPhysicsConstructor {
5059
*/
5160
void ConstructProcess();
5261

62+
/**
63+
* Search the list for the process "photoNuclear". When found,
64+
* change the calling order so photonNuclear is called before
65+
* EM processes. The biasing operator needs the photonNuclear
66+
* process to be called first because the cross-section is
67+
* needed to bias down other process.
68+
*/
69+
void SetPhotonNuclearAsFirstProcess() const;
70+
5371
private:
72+
/*
73+
* Returns the process manager object for the G4Gamma class from the list of
74+
* Geant4 particles.
75+
*/
76+
G4ProcessManager* GetGammaProcessManager() const;
5477
/**
5578
* The gamma to muons process.
5679
*/
5780
G4GammaConversionToMuons gammaConvProcess;
81+
/**
82+
* Parameters from the configuration to pass along to the photonuclear model.
83+
*/
84+
framework::config::Parameters modelParameters;
5885
};
5986

6087
} // namespace simcore

‎include/SimCore/PhotoNuclearModel.h‎

Lines changed: 98 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,98 @@
1+
#ifndef SIMCORE_PHOTONUCLEAR_MODEL_H
2+
#define SIMCORE_PHOTONUCLEAR_MODEL_H
3+
#include <G4CrossSectionDataSetRegistry.hh>
4+
#include <G4HadronInelasticProcess.hh>
5+
#include <G4HadronicInteraction.hh>
6+
#include <G4PhotoNuclearCrossSection.hh>
7+
#include <G4ProcessManager.hh>
8+
#include <string>
9+
#include <utility>
10+
11+
#include "Framework/Configure/Parameters.h"
12+
#include "SimCore/Factory.h"
13+
/*
14+
** Dynamically loadable photonuclear models either from SimCore or external
15+
** libraries implementing this interface. For example implementations, see the
16+
** BertiniNothingHard model in SimCore.
17+
**
18+
** Allows for replacing the default Bertini model from Geant4 with any other
19+
** G4HadronicInteraction process. The library is used from within the
20+
** GammaPhysics module in SimCore which ensures that the removeExistingModel and
21+
** ConstructModel functions are called in the right order and that the
22+
** photonNuclear process is located in the right part of the G4Gamma process
23+
** list.
24+
*/
25+
namespace simcore {
26+
27+
class PhotoNuclearModel {
28+
public:
29+
/**
30+
* Base class does not take any parameters or do anything in particular, but
31+
* any derived class may.
32+
*
33+
* @param[in] name unique instance name for this model
34+
* @param[in] parameters python configuration
35+
*/
36+
PhotoNuclearModel(const std::string& name,
37+
const framework::config::Parameters& parameters){};
38+
virtual ~PhotoNuclearModel() = default;
39+
40+
/**
41+
* The primary part of the model interface, responsible for adding the desired
42+
* G4HadronicInteraction to the process manager for the G4Gamma class.
43+
*
44+
* @param[in] processManager the process manager for the G4Gamma class, passed
45+
* in automatically by the GammaPhysics module.
46+
*/
47+
virtual void ConstructGammaProcess(G4ProcessManager* processManager) = 0;
48+
49+
/**
50+
* The factory for PhotoNuclearModels.
51+
*/
52+
using Factory =
53+
::simcore::Factory<PhotoNuclearModel, std::shared_ptr<PhotoNuclearModel>,
54+
const std::string&,
55+
const framework::config::Parameters&>;
56+
57+
/**
58+
* Removes any existing photonNuclear process from the process manager of the
59+
* G4Gamma class. Should in general not be overridden for most models other
60+
* than the default Bertini model (which just retains the default
61+
* interaction).
62+
*
63+
* @param[in] processManager the process manager for the G4Gamma class, passed
64+
* in automatically by the GammaPhysics module.
65+
*/
66+
virtual void removeExistingModel(G4ProcessManager* processManager);
67+
68+
/**
69+
* Default implementation for adding XS data for the process.
70+
*
71+
* The default implementation is adapted from G4PhotoNuclearProcess.hh but can
72+
* be overridden to add cross sections in another way (e.g. from FLUKA). If no
73+
* cross-section is added to the process, the simulation process will halt
74+
* when attempting to calculate the mean free path so there is no way of
75+
* accidentally forgetting to enable the XS data.
76+
*
77+
* Typically called during `ConstructGammaProcess`.
78+
*
79+
*/
80+
virtual void addPNCrossSectionData(G4HadronInelasticProcess* process) const;
81+
};
82+
} // namespace simcore
83+
84+
/**
85+
* @macro DECLARE_PHOTONUCLEAR_MODEL
86+
*
87+
* Defines a builder for the declared derived photonuclear model class and
88+
* registers it as a possible photonuclear model. If you are implementing your
89+
* own photonuclear model class, make sure to invoke this macro in your
90+
* implementation file.
91+
*
92+
* See e.g. SimCore/src/SimCore/PhotoNuclearModels/BertiniModel.cxx
93+
*/
94+
#define DECLARE_PHOTONUCLEAR_MODEL(CLASS) \
95+
namespace { \
96+
auto v = ::simcore::PhotoNuclearModel::Factory::get().declare<CLASS>(); \
97+
}
98+
#endif /* SIMCORE_PHOTONUCLEAR_MODEL_H */

‎include/SimCore/PhotoNuclearModels/BertiniAtLeastNProductsModel.h‎

Lines changed: 63 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,63 @@
1+
#ifndef SIMCORE_BERTINI_AT_LEAST_N_PRODUCTS_MODEL_H
2+
#define SIMCORE_BERTINI_AT_LEAST_N_PRODUCTS_MODEL_H
3+
#include <G4CrossSectionDataSetRegistry.hh>
4+
#include <G4Gamma.hh>
5+
#include <G4HadProjectile.hh>
6+
#include <G4HadronInelasticProcess.hh>
7+
#include <G4Nucleus.hh>
8+
#include <G4ProcessManager.hh>
9+
10+
#include "Framework/Configure/Parameters.h"
11+
#include "SimCore/PhotoNuclearModel.h"
12+
#include "SimCore/PhotoNuclearModels/BertiniEventTopologyProcess.h" /* */
13+
namespace simcore {
14+
class BertiniAtLeastNProductsProcess : public BertiniEventTopologyProcess {
15+
public:
16+
BertiniAtLeastNProductsProcess(double threshold, int Zmin, double Emin,
17+
std::vector<int> pdg_ids, int min_products)
18+
: BertiniEventTopologyProcess{},
19+
threshold_{threshold},
20+
Zmin_{Zmin},
21+
Emin_{Emin},
22+
pdg_ids_{pdg_ids},
23+
min_products_{min_products} {}
24+
virtual ~BertiniAtLeastNProductsProcess() = default;
25+
bool acceptProjectile(const G4HadProjectile& projectile) const override {
26+
return projectile.GetKineticEnergy() >= Emin_;
27+
}
28+
bool acceptTarget(const G4Nucleus& targetNucleus) const override {
29+
return targetNucleus.GetZ_asInt() >= Zmin_;
30+
}
31+
bool acceptEvent() const override;
32+
33+
private:
34+
double threshold_;
35+
int Zmin_;
36+
double Emin_;
37+
std::vector<int> pdg_ids_;
38+
int min_products_;
39+
};
40+
41+
class BertiniAtLeastNProductsModel : public PhotoNuclearModel {
42+
public:
43+
BertiniAtLeastNProductsModel(const std::string& name,
44+
const framework::config::Parameters& parameters)
45+
: PhotoNuclearModel{name, parameters},
46+
threshold_{parameters.getParameter<double>("hard_particle_threshold")},
47+
Zmin_{parameters.getParameter<int>("zmin")},
48+
Emin_{parameters.getParameter<double>("emin")},
49+
pdg_ids_{parameters.getParameter<std::vector<int>>("pdg_ids")},
50+
min_products_{parameters.getParameter<int>("min_products")} {}
51+
virtual ~BertiniAtLeastNProductsModel() = default;
52+
void ConstructGammaProcess(G4ProcessManager* processManager) override;
53+
54+
private:
55+
double threshold_;
56+
int Zmin_;
57+
double Emin_;
58+
std::vector<int> pdg_ids_;
59+
int min_products_;
60+
};
61+
62+
} // namespace simcore
63+
#endif /* SIMCORE_BERTINI_AT_LEAST_N_PRODUCTS_MODEL_H */

‎include/SimCore/PhotoNuclearModels/BertiniEventTopologyProcess.h‎

Lines changed: 137 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,137 @@
1+
#ifndef SIMCORE_BERTINI_EVENTTOPOLOGY_PROCESS_H
2+
#define SIMCORE_BERTINI_EVENTTOPOLOGY_PROCESS_H
3+
4+
#include <G4CascadeInterface.hh>
5+
#include <G4EventManager.hh>
6+
#include <G4HadFinalState.hh>
7+
#include <G4HadProjectile.hh>
8+
#include <G4HadronicInteraction.hh>
9+
#include <G4Nucleus.hh>
10+
#include <iostream>
11+
12+
#include "SimCore/PhotoNuclearModel.h"
13+
#include "SimCore/UserEventInformation.h"
14+
namespace simcore {
15+
16+
/*
17+
** A wrapper interface around the Bertini cascade process (G4CascadeInterface)
18+
** which reruns the event generator until a particular condition is met for the
19+
** products. The decision of whether or not to rerun the event generator is
20+
** handled by the acceptEvent virtual function.
21+
**
22+
** For example of a derived class, see
23+
** SimCore/include/SimCore/PhotoNuclearModels/BertiniNothingHardModel.h
24+
**
25+
** Note: You almost certainly want to use these types of models together with a
26+
** photonuclear topology filter.
27+
**
28+
** Note: When performing N attempts, this will increment the event weight in the
29+
** UserEventInformation by 1/N. To change this behaviour, override the
30+
** incrementEventWeight function.
31+
*/
32+
33+
class BertiniEventTopologyProcess : public G4CascadeInterface {
34+
public:
35+
BertiniEventTopologyProcess(bool count_light_ions = true)
36+
: G4CascadeInterface{}, count_light_ions_{count_light_ions} {}
37+
38+
/*
39+
* The primary function for derived classes to customize. After each call to
40+
* the Bertini cascade, this function will be called to see whether or not to
41+
* keep the event. The products can be accessed from the `theParticleChange`
42+
* member inherited from G4CascadeInterface (i.e. the Bertini cascade).
43+
*/
44+
virtual bool acceptEvent() const = 0;
45+
46+
/*
47+
* Is the projectile of interest?
48+
*
49+
* If false, the cascade will only run once. Example use includes only
50+
* applying repeated simulations for particular energy ranges. Is called
51+
* automatically during `ApplyYourself`.
52+
*
53+
**/
54+
virtual bool acceptProjectile(const G4HadProjectile& projectile) const = 0;
55+
56+
/*
57+
* Is the target nucleus of interest?
58+
*
59+
* If false, the cascade will only run once. Example use includes only
60+
* applying repeated simulations for high Z materials. Is called
61+
* automatically during `ApplyYourself`.
62+
*
63+
**/
64+
virtual bool acceptTarget(const G4Nucleus& targetNucleus) const = 0;
65+
66+
/*
67+
* Run the Bertini cascade until the condition given by `acceptEvent` is
68+
* matched by the reaction products. Increments the event weight by 1/n where
69+
* n is the number of attempts needed to produce a match.
70+
*
71+
**/
72+
G4HadFinalState* ApplyYourself(const G4HadProjectile& projectile,
73+
G4Nucleus& targetNucleus) override;
74+
75+
/*
76+
* Geant4 assumes that secondaries produced from the bertini cascade are owned
77+
* by some other part of the code. Since we are re-running the cascade until
78+
* we get something matching our condition in `acceptEvent`, we have to make
79+
* sure to clean up any secondaries from failed attempts.
80+
*
81+
*/
82+
void cleanupSecondaries();
83+
84+
/**
85+
* Check if the PDG code corresponds to a light ion nucleus.
86+
*
87+
* Nuclear PDG codes are given by ±10LZZZAAAI So to find the atomic
88+
* number, we first divide by 10 (to lose the I-component) and then
89+
* take the modulo with 1000.
90+
*
91+
*/
92+
constexpr bool isLightIon(const int pdgCode) const {
93+
//
94+
// TODO: Is the < check necessary?
95+
if (pdgCode > 1000000000 && pdgCode < 10000000000) {
96+
// Check if the atomic number is less than or equal to 4
97+
return ((pdgCode / 10) % 1000) <= 4;
98+
}
99+
return false;
100+
}
101+
102+
/**
103+
* Whether or not to include a particular particle type in any counting.
104+
* Unless \ref count_light_ions_ is set, we don't count anything with a
105+
* nuclear PDG code. This is consistent with the counting behaviour used in
106+
* the PhotoNuclearDQM.
107+
*
108+
* If \ref count_light_ions_ is set, we also match PDG codes for nuclei with
109+
* atomic number < 4. 
110+
*
111+
* @see isLightIon
112+
*
113+
*/
114+
constexpr bool skipCountingParticle(const int pdgcode) const {
115+
return !(pdgcode < 10000 || count_light_ions_ && isLightIon(pdgcode));
116+
}
117+
118+
/*
119+
* Update the event weight depending on the number of attempts made.
120+
*
121+
* @param N The number of attempts used for a successful topology
122+
* production.
123+
*
124+
**/
125+
126+
virtual void incrementEventWeight(int N) {
127+
auto event_info{static_cast<UserEventInformation*>(
128+
G4EventManager::GetEventManager()->GetUserInformation())};
129+
event_info->incWeight(1. / N);
130+
}
131+
132+
private:
133+
bool count_light_ions_;
134+
};
135+
} // namespace simcore
136+
137+
#endif /* SIMCORE_BERTINI_EVENTTOPOLOGY_PROCESS_H */

0 commit comments

Comments
 (0)