Weihe Zeng
Tsinghua University
Dec. 11, 2017
Content
- Physics models
- electromagnetic/hadronic model inventory
- reference
physicslist
and use cases - build your own physics model
- Geant4 Simulation Control
- Initialization, idle and "beam on"
- Basic structure:
run -> event -> track -> step
- User action(hook) for information extraction
Interaction of particles (e.g. $\gamma$)
Physics processes in example B1
/process/list
Output of this UI command:
Transportation, msc, hIoni, ionIoni
hBrems, hPairProd, CoulombScat, eIoni
eBrem, annihil, phot, compt
conv, muIoni, muBrems, muPairProd
photonNuclear, electronNuclear, positronNuclear, muonNuclear
Decay, hadElastic,hFritiofCaptureAtRest,hBertiniCaptureAtRest
muMinusCaptureAtRest, dInelastic, tInelastic, He3Inelastic
alphaInelastic, ionInelastic, anti_He3Inelastic,anti_alphaInelastic
anti_deuteronInelastic,anti_lambdaInelastic,anti_neutronInelastic,anti_omega-Inelastic
anti_protonInelastic,anti_sigma+Inelastic,anti_sigma-Inelastic,anti_tritonInelastic
anti_xi-Inelastic, anti_xi0Inelastic, kaon+Inelastic, kaon-Inelastic
kaon0LInelastic, kaon0SInelastic, lambdaInelastic, neutronInelastic
nCapture, omega-Inelastic, pi+Inelastic, pi-Inelastic
protonInelastic, sigma+Inelastic, sigma-Inelastic, sigma0Inelastic
xi-Inelastic, xi0Inelastic, nKiller
Electromagnetic models
Hadronic models
GEANT4 packaged physics lists
See header files in source/physics_lists/lists/include
// use built-in physics lists
FTFP_BERT* physlist = new FTFP_BERT(verbose);
runManager->SetUserInitialization(physlist);
Reference Physics Lists:
- FTFPBERT, FTFPBERT_HP
- QGSPBERT, QGSPBERT_HP
- QGSP_BIC
- QGSPFTFPBERT
Recommended use cases
| Use case | Physics list | | ---------------------------------------- | ---------------------------------------- | | High energy physics calorimetry | FTFPBERT, QGSPBERT | | Linear collider neutron fluxes | QGSPBERTHP, QGSPBICHP | | Shielding applications | QGSPBERTHP, QGSPBICHP, QGSPINCLXX, Shielding | | Low energy dosimetric applications | QGSPBERTHP, QGSPBICHP | | Medical and industrial neutron applications | QGSPBERTHP, QGSPBICHP, QGSPINCLXX | | Underground physics | Shielding |
ref: use case
Build Your Own Physics List
Abstract Class for Physics list is G4VUserPhysicsList
. In main()
function, we set our concrete implementation of physics list as initialization for G4RunManager
using the following method:
void SetUserInitialization(G4VUserPhysicsList* userInit);
To build our own physics list, i.e. implement G4VUserPhysicsList
, we can derive a concrete class and complete all the (pure) virtual of it, or use another higher level abstract class G4VModularPhysicsList
.
Virtual functions in G4VUserPhysicsList
// mandatory
virtual void ConstructParticle() = 0;
virtual void ConstructProcess() = 0;
//optional
// mistake in Textbook 6.1.2.4
virtual void SetCuts();
Build physics list from scratch
examples/extended/optical/OpNovice
void OpNovicePhysicsList::ConstructParticle()
{
G4BosonConstructor bConstructor;
bConstructor.ConstructParticle();
G4LeptonConstructor lConstructor;
lConstructor.ConstructParticle();
G4MesonConstructor mConstructor;
mConstructor.ConstructParticle();
G4BaryonConstructor rConstructor;
rConstructor.ConstructParticle();
G4IonConstructor iConstructor;
iConstructor.ConstructParticle();
}
void OpNovicePhysicsList::ConstructProcess(){
AddTransportation(); // implement in G4VUserPhysicsList
ConstructDecay(); // new methods in this class
ConstructEM();
ConstructOp();
}
void OpNovicePhysicsList::ConstructEM(){
auto particleIterator=GetParticleIterator(); // new in 10.3
particleIterator->reset();
while( (*particleIterator)() ){
G4ParticleDefinition* particle = particleIterator->value();
G4ProcessManager* pmanager = particle->GetProcessManager();
G4String particleName = particle->GetParticleName();
if (particleName == "gamma") {
// Construct processes for gamma
pmanager->AddDiscreteProcess(new G4GammaConversion());
pmanager->AddDiscreteProcess(new G4ComptonScattering());
pmanager->AddDiscreteProcess(new G4PhotoElectricEffect());
} else if (particleName == "e-") {
pmanager->AddProcess(new G4eMultipleScattering(),-1, 1, 1);
pmanager->AddProcess(new G4eIonisation(), -1, 2, 2);
pmanager->AddProcess(new G4eBremsstrahlung(), -1, 3, 3);
/* ... */
}
Energy Cut v.s. Range Cut (Chapter 5.4)
To avoid infrared divergence, some electromagnetic processes require a threshold below which no secondary will be generated. Because of this requirement, gammas, electrons and positrons require production threshold. This threshold should be defined as a distance, or range cut-off, which is internally converted to an energy for individual materials.
void OpNovicePhysicsList::SetCuts(){
if (value<0.0) return;
defaultCutValue = 1.0*mm;
// set cut values for gamma at first and for e- and e+
SetCutValue(defaultCutValue, "gamma");
SetCutValue(defaultCutValue, "e-");
SetCutValue(defaultCutValue, "e+");
SetCutValue(defaultCutValue, "proton");
}
Build Physics List using Modular Builder
example/basic/B3
class B3PhysicsList: public G4VModularPhysicsList{
public:
B3PhysicsList();
virtual ~B3PhysicsList();
virtual void SetCuts();
};
B3PhysicsList::B3PhysicsList() : G4VModularPhysicsList(){
SetVerboseLevel(1);
// Default physics
RegisterPhysics(new G4DecayPhysics());
// Radioactive decay
RegisterPhysics(new G4RadioactiveDecayPhysics());
// EM physics
RegisterPhysics(new G4EmStandardPhysics());
}
You can register physics in source/physics_lists/constructors/
. Only one physics constructor can be registered for each "physics_type".
Some built-in Physics Lists are also constructed using G4VModularPhysicsList
QBBC::QBBC( G4int ver, const G4String&){
G4DataQuestionaire it(photon, neutronxs);
G4cout << "<<< Reference Physics List QBBC "
<<G4endl;
defaultCutValue = 0.7*mm;
SetVerboseLevel(ver);
// EM Physics
RegisterPhysics( new G4EmStandardPhysics(ver) );
// Synchroton Radiation & GN Physics
RegisterPhysics( new G4EmExtraPhysics(ver) );
// Decays
RegisterPhysics( new G4DecayPhysics(ver) );
// Hadron Physics
RegisterPhysics( new G4HadronElasticPhysicsXS(ver) );
RegisterPhysics( new G4StoppingPhysics(ver) );
RegisterPhysics( new G4IonPhysics(ver) );
RegisterPhysics( new G4HadronInelasticQBBC(ver));
// Neutron tracking cut
RegisterPhysics( new G4NeutronTrackingCut(ver) );
}
Geant4 Simulation Management
Basic Simulation Structure
Run
- As an analogy of the real experiment, a run of Geant4 starts with “Beam On”.
- Within a run, the user cannot change
- detector setup
- settings of physics processes
- Conceptually, a run is a collection of events which share the same detector and physics conditions.
G4RunManager
class manages processing a run, a run is represented by G4Run class or a user-defined class derived from G4Run.- A run class may have a summary results of the run.
G4UserRunAction
is the optional user hook.
Event
- An event is the basic unit of simulation in Geant4.
- At beginning of processing, primary tracks are generated. These primary tracks are pushed into a stack.
- A track is popped up from the stack one by one and “tracked”. Resulting secondary tracks are pushed into the stack.
- When the stack becomes empty, processing of one event is over.
G4Event
class represents an event.G4EventManager
class manages processing an event.G4UserEventAction
is the optional user hook.
Stack (data structure)
Stacking mechanism in GEANT4
Use G4UserStackingAction
to manage tracks before they're pushed into stack.
G4ClassificationOfNewTrack ClassifyNewTrack(const G4Track*)
return values: fUrgent, fWaiting, fPostpone, fKill
Track
- Track is a snapshot of a particle.
- It has physical quantities of current instance only. It does not record previous quantities.
- Step is a “delta” information to a track. Track is not a collection of steps. Instead, a track is being updated by steps.
- Track object is deleted when
- it goes out of the world volume,
- it disappears (by e.g. decay, inelastic scattering),
- it goes down to zero kinetic energy
- the user decides to kill it artificially.
- No track object persists at the end of event.
G4TrackingManager
manages processing a track, a track is represented byG4Track
class.G4UserTrackingAction
is the optional user hook.
Step
- Step has two points and also “delta” information of a particle (energy loss on the step, time-of-flight spent by the step, etc.).
- Each point knows the volume (and material). In case a step is limited by a volume boundary, the end point physically stands on the boundary, and it logically belongs to the next volume.
G4SteppingManager
class manages processing a step, a step is represented byG4Step
class.G4UserSteppingAction
is the optional user hook.
Getting information from simulation (in reverse order)
G4UserSteppingAction
void B1SteppingAction::UserSteppingAction(const G4Step* step){
// get logical volume of the current step
G4LogicalVolume* volume = step->GetPreStepPoint()
->GetTouchableHandle()->GetVolume()->GetLogicalVolume();
// collect information of this step
G4double edepStep = step->GetTotalEnergyDeposit();
// record information in Event/Tracking Action
fEventAction->AddEdep(edepStep);
}
Recall that Step is a “delta” information of a track. We should pass the information (energy deposition in this case) to upper level (tracking/event) and record it (sum up in this case).
Some other information in step that would be useful:
step->GetPostStepPoint()->GetPosition().x();
step->GetPostStepPoint()->GetPosition().y();
step->GetPostStepPoint()->GetPosition().z();
step->GetPreStepPoint()->GetPosition().z();
step->GetPreStepPoint()->GetGlobalTime();
// get current track information
G4Track* track = step->GetTrack();
/* some examples of track information */
track->GetParticleDefinition()->GetPDGEncoding();
track->GetTrackID();
track->GetParentID();
track->GetCurrentStepNumber();
track->GetKineticEnergy();
G4UserEventAction
// function that get information from stepping action
void AddEdep(G4double edep) { fEdep += edep; }
// reset total energy deposition to 0
void B1EventAction::BeginOfEventAction(const G4Event*){
fEdep = 0.;
}
// passing energy deposition of one event to Run action
void B1EventAction::EndOfEventAction(const G4Event*){
// accumulate statistics in run action
fRunAction->AddEdep(fEdep);
}
In this example (examples/basic/B1
), user passes recorded energy deposition to the Run level.
G4UserRunAction
void B1RunAction::AddEdep(G4double edep){
fEdep += edep, fEdep2 += edep*edep;}
void B1RunAction::BeginOfRunAction(const G4Run*){
fEdep = 0.0, fEdep2 = 0.0;}
void B1RunAction::EndOfRunAction(const G4Run* run){
G4int nofEvents = run->GetNumberOfEvent();
// Compute dose = total energy deposit in a run and its variance
G4double rms = fEdep2 - fEdep*fEdep/nofEvents;
if (rms > 0.) rms = std::sqrt(rms); else rms = 0.;
// get mass from detector logical volume
G4double dose = edep/mass;
G4double rmsDose = rms/mass;
// output results to standard output
G4cout << " The run consists of " << nofEvents << G4endl
<< " Cumulated dose per run, in scoring volume : "
<< G4BestUnit(dose,"Dose") << " rms = "
<< G4BestUnit(rmsDose,"Dose") << G4endl}
Here we write average dose (energy deposition) to standard output. For other simulation purposes, we can also record energy deposition for each event in an array/vector
and write to a file. (Next week)
Homeworks
- Read remaining code in
examples/basic/B1
and understand the information passing mechanism. - (optional) Install ROOT on your computer.
Resources and references
- Tutorial handouts @ MIT
- Built-in physicslist use cases
- Textbook: Chapter 5 & 6, page 168-243