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$)

online image

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 by G4Track 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 by G4Step 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