I am a newbie with Geant4. I want to use it for dumping some proton cross sections to the terminal. I came up with the MWE seen below.
My question is whether it’s possible to simplify this example further. In other words, can I skip the world generation with “air”? Do I really need to define a particle gun or is there a simpler approach? Shouldn’t the cross sections just depend on the active physics list and the particle / material, and not on the “world” and “run” ?
I was thinking that something like G4CrossSectionDataStore cds; G4DynamicParticle dp(particle, 230 * MeV, G4ThreeVector(1.0,0,0)); cds.GetCrossSection(&dp, elm, mat); could be enough, but it crashes, probably because it’s missing the Datasets that it should look for. Any help in simplifying this example is greatly appreciated.
#include "G4HadronicProcessStore.hh"
#include "G4NistManager.hh"
#include "G4SystemOfUnits.hh"
#include "G4RunManager.hh"
#include "G4PhysListFactory.hh"
#include "G4VUserDetectorConstruction.hh"
#include "G4NistManager.hh"
#include "G4Box.hh"
#include "G4PVPlacement.hh"
#include "G4LogicalVolume.hh"
#include "G4VUserPrimaryGeneratorAction.hh"
#include "G4ParticleGun.hh"
#include "G4Proton.hh"
G4String elementName = "O";
class BasicDetectorConstruction : public G4VUserDetectorConstruction{
public:
G4VPhysicalVolume* Construct()
{
//Obtain pointer to NIST material manager
G4NistManager* nist = G4NistManager::Instance();
//Build materials for world and target
G4Material* mat = nist->FindOrBuildMaterial("G4_"+elementName);
G4Material* air = nist->FindOrBuildMaterial("G4_AIR");
//Define the world solid. The world will be 20m x 20m x20m
G4double world_sizeX=10*m;
G4double world_sizeY=10*m;
G4double world_sizeZ=10*m;
G4Box* solidWorld =
new G4Box("World",world_sizeX,world_sizeY,world_sizeZ);
//Fill the world with air
G4LogicalVolume* logicWorld =
new G4LogicalVolume(solidWorld, air, "myWorld");
//Create the world physical volume. The world is the only
//physical volume with no mother volume.
G4VPhysicalVolume* physWorld =
new G4PVPlacement(0, //no rotation
G4ThreeVector(), //at (0,0,0)
logicWorld, //its logical volume
"World", //its name
0, //its mother volume
false, //no boolean operation
0, //copy number
true); //overlaps checking
//Create the shape of a target to fire particles at
G4double target_sizeX=1*m;
G4double target_sizeY=1*m;
G4double target_sizeZ=1*m;
G4Box* solidTarget =
new G4Box("Target",target_sizeX, target_sizeY, target_sizeZ);
//Create the target logical volume by
//assigning the material of the target to be mat
G4LogicalVolume* logicTarget =
new G4LogicalVolume(solidTarget, mat, "myTarget");
//Create the target physical volume by placing it in the
//"logicWorld" logical volume.
G4VPhysicalVolume* physTarget =
new G4PVPlacement(0, //no rotation
G4ThreeVector(), //at (0,0,0)
logicTarget, //its logical volume
"World", //its name
logicWorld, //its mother volume
false, //no boolean operation
0, //copy number
true); //overlaps checking
return physWorld;
}
};
class BasicPrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction{
public:
void GeneratePrimaries(G4Event* anEvent){
//Create particle gun
G4ParticleGun* myGun = new G4ParticleGun();
//Specify particle to be emitted
myGun->SetParticleDefinition(G4Proton::ProtonDefinition());
//Set particle kinetic energy and direction of travel
myGun->SetParticleEnergy(1.*MeV);
myGun->SetParticleMomentumDirection(G4ThreeVector(1.0,0,0));
//Generate one instance of specified particle
myGun->GeneratePrimaryVertex(anEvent);
}
};
int main(int argc,char** argv)
{
G4RunManager * runManager = new G4RunManager;
runManager->SetUserInitialization(new BasicDetectorConstruction());
// Physics list
G4PhysListFactory factory;
G4VModularPhysicsList* physicsList = factory.GetReferencePhysList("QGSP_BIC_AllHP");// FTFP_BERT
runManager->SetUserInitialization(physicsList);
// Primary generator action
runManager->SetUserAction(new BasicPrimaryGeneratorAction());
// Initialize G4 kernel
runManager->Initialize();
//Cause the run manager to generate a single event using the
//primary generator action registered above.
runManager->BeamOn(1);
G4HadronicProcessStore* store = G4HadronicProcessStore::Instance();
const G4ParticleDefinition* particle =
G4ParticleTable::GetParticleTable()->FindParticle("proton");
G4cout << "-------------------------------------------------------------" << G4endl;
G4cout << " N E(MeV) Elastic(b) Inelastic(b)";
G4cout << " Total(b)" << G4endl;
G4cout << "-------------------------------------------------------------" << G4endl;
const G4Element* elm = G4NistManager::Instance()->FindOrBuildElement(elementName);
const G4Material* mat = G4NistManager::Instance()->FindOrBuildMaterial("G4_" + elementName);
G4double e, p, xs, xtot;
for (auto i = 1; i < 230; i++) {
e = i * MeV;
G4cout << std::setw(5) << i << std::setw(12) << e;
xs = store->GetElasticCrossSectionPerAtom(particle, e, elm, mat);
xtot = xs;
G4cout << std::setw(12) << xs / barn;
xs = store->GetInelasticCrossSectionPerAtom(particle, e, elm, mat);
xtot += xs;
G4cout << " " << std::setw(12) << xs / barn;
G4cout << " " << std::setw(12) << xtot / barn << G4endl;
}
//After the run is complete, free up the memory used by run manager and return
delete runManager;
return 0;
}
Geant4 Version: 4.11.03
Operating System: AlmaLinux 9
Compiler/Version: g++ 11.5
CMake Version: 3.26