Simulation code

Structure of the program and relevant subroutines

Geometry description

In testDetectorConstruction.hh and testDetectorConstruction.cc, the geometry is defined. The detector is a compact cylinder of radius 2.5 cm and a height of 5 cm, with its axis parallel to the Z axis, and separated 10 cm from the origin of coordinates. A simulation was done with a height of 10 cm too. The material was set to Ge, and later on, to CsI. The cylinder is placed in the vacuum.

// Detector as a Germanium cylinder
G4double minRadius  = 0.0*cm;
G4double maxRadius  = 2.5*cm;
G4double startPhi   = 0.*deg;
G4double deltaPhi   = 360.*deg;
G4double halfLength = 0.25*cm;

G4double distance              = 10*cm + halfLength;
G4ThreeVector positionDetector = G4ThreeVector(0.*cm, 0.*cm, distance);

G4Tubs* solidDetector            = new G4Tubs("Detector", minRadius, maxRadius, halfLength, startPhi, deltaPhi);
G4LogicalVolume* logicDetector   = new G4LogicalVolume(solidDetector, Ge, "Detector");
G4VPhysicalVolume* physiDetector = new G4PVPlacement(0, positionDetector, "Detector", logicDetector, physiWorld, false, 0);

Cylindrical geometry setup Detector front view The figures above show the cylinder geometry and position.

Event generator

In testPrimaryGenerationAction.hh and testPrimaryGenerationAction.cc, the primary generator is set to a gamma of 1 MeV, which is launched from the origin with a fixed direction, fixed by the angular coordinates θ and φ, whose value is generated randomly. The minimum value of θ is 0 degrees, and the maximum is fixed by the system's aperture, which in this case is 14 degrees. The φ angle goes from 0 to 360 degrees. With these values, we create a unitary vector in spherical coordinates that sets the direction of the gammas momentum.

Calculation of the azimuthal angle. Knowing the radius and the distance from the origin, we can calculate the azimuthal angle.

void testPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) {

	// Generation of primary particles

	// source energy
	particleGun->SetParticleEnergy(1.0*MeV);

	// source position
	G4double x0 = 0.*cm;
	G4double y0 = 0.*cm;
	G4double z0 = 0.*cm;
	particleGun->SetParticlePosition(G4ThreeVector(x0,y0,z0));

	// random direction within aperture
	G4double aperture = 14.; // in degrees

	if(aperture > 0.) {

		G4double theta = acos(1. + (cos(aperture/180.*pi) -1.)*G4UniformRand());
		G4double phi   = G4UniformRand()*twopi;

		G4double ux = sin(theta)*cos(phi);
		G4double uy = sin(theta)*sin(phi);
		G4double uz = cos(theta);

		particleGun->SetParticleMomentumDirection(G4ThreeVector(ux,uy,uz));

	} else {
		particleGun->SetParticleMomentumDirection(G4ThreeVector(0., 0., 1.));
	}

	particleGun->GeneratePrimaryVertex(anEvent);
}

Particles and processes included

In testPhysicsList.hh and testPhysicsList.cc, gammas, electrons and positrons are created (&gamma, e^- and e⁺). The physical processes are different for each particle. For gammas, we ste discrete processes like Compton, photoelectric and pair production, and for electrons/positrons, multiple collisions, ionization and bremsstrahlung. For positrons, we add anihillation as well.

if (particleName == "gamma") {     // gamma

	pmanager->AddDiscreteProcess(new G4PhotoElectricEffect);
	pmanager->AddDiscreteProcess(new G4ComptonScattering);
	pmanager->AddDiscreteProcess(new G4GammaConversion);

} else if (particleName == "e-") { //electron

	pmanager->AddProcess(new G4MultipleScattering, -1,  1, 1);
	pmanager->AddProcess(new G4eIonisation,        -1,  2, 2);
	pmanager->AddProcess(new G4eBremsstrahlung,    -1, -1, 3);

} else if (particleName == "e+") { //positron

	pmanager->AddProcess(new G4MultipleScattering, -1,  1, 1);
	pmanager->AddProcess(new G4eIonisation,        -1,  2, 2);
	pmanager->AddProcess(new G4eBremsstrahlung,    -1, -1, 3);
	pmanager->AddProcess(new G4eplusAnnihilation,   0, -1, 4);
}

Bremsstrahlung and ionization mix discrete and continuous processes, so cuts are defined for the energy. Over this energy cuts, we consider the processes as discrete, below it, as continuous. In this way, if the energy loss is very low we can consider it approximately continuous, and if it is very big, we don't follow the particle.

testPhysicsList::testPhysicsList() : G4VUserPhysicsList() {

	// 10 mu default cut value
	defaultCutValue = 0.01*mm;
}

Gathering information

According to Geant4 terminology, a run is a collection of events. Through testRunAction and testEventAction, we have control over the start and end of a run and an event, and in testSteppingAction we have control over the end of each step.

void testSteppingAction::UserSteppingAction(const G4Step* aStep) {

	// Collect the energy deposited in the sensitive volume
	const G4String currentPhysicalName = aStep->GetPreStepPoint()->GetPhysicalVolume()->GetName();

	if (currentPhysicalName == "Detector") {

		G4double EdepStep = aStep->GetTotalEnergyDeposit();

		if (EdepStep > 0.)
			eventAction->addEdep(EdepStep);
	 }
}

void testEventAction::EndOfEventAction(const G4Event* evt) {

	if (TotalEnergyDeposit > 0.) {

		// fill histograms
		runAction-> fillHistogram(TotalEnergyDeposit);

		// add total counts
		runAction-> addTotalCount();
	}
}

void testRunAction::EndOfRunAction(const G4Run*) {

	//write efficiency
	ofstream efficiencyOutputFile("test.eff");

	efficiencyOutputFile << " Efficiency for test detector" << G4endl;
	efficiencyOutputFile << " Total efficiency = " << TotalCounts << G4endl;

	//write histograms
	ofstream histogramOutputFile("test.vec");

	if (NumberChannel > 0) {
		for(G4int ind=0; ind < NumberChannel; ind++) {
			histogramOutputFile << setw(6) << ind+1 <<
			setw(10) << Histogram[ind] << G4endl;
		}
	}
}

In testSteppingAction.hh and testSteppingAction.cc, we collect the energy deposited in each step of an event, and in testEventAction.hh and testEventAction.cc we collect the total energy deposited. Later, in testRunAction.hh and testRunAction.cc the histograms and efficiency are generated and saved in test.vec (a file with two columns, one for the channel, and another for the content of the channel) and test.eff (which contains the value of the sum of the total counts or the peak counts).

void testEventAction::addEdep(G4double Edep) {
	TotalEnergyDeposit += Edep;
}

The batch file hplot.cc is used to plot the histogram using the ROOT command interpreter.

Visualization of the geometry

The macro vis.mac allows us to visualize and manipulate our geometric configuration using DAWN. If we modify the runs, the type of particle and the energy, we can try different experimental configurations.

# Visualization and execution macro
#    /control/execute vis.mac

# Visualization commands
/vis/scene/create
/vis/open DAWNFILE

/vis/viewer/reset
/vis/viewer/set/style w
/vis/viewer/set/viewpointThetaPhi 40 40
/vis/viewer/zoom 2

/vis/scene/add/axes 0 0 0 0.05
/vis/scene/add/volume

###################################
# Uncomment these for tracks visualization

/tracking/storeTrajectory 1
/vis/scene/add/trajectories
/vis/scene/endOfEventAction accumulate

## Particle type and energy
#/gun/particle e+
#/gun/energy 1.0 MeV

## Number of events to run
/run/beamOn 10

# End of uncomment
###################################

# Refresh view
/vis/viewer/flush

Tracks inside the cylinder, snapshot 1 The figures above show some processes happening already, like pair creation for example.

Tracks inside the cylinder, snapshot 3 Examples with a lot of runs. Now we can see more processes like anihillation, pair creation, etc. The second image shows the cylinder with a height of 10 cm.

Detailed information about the tracks

The macro print-track.mac obtains information on what's happening to the photons, electrons and positrons in real time, and whether they leave the material, which of the processes we defined are they suffering, which secondary particles are produced in those interactions, how much energy are they losing, etc. We can see the processes as they happen and which particles they generate run by run for each event.

# Execution macro to print track information
#    /control/execute print-track.mac

# Information level
/run/verbose 0
/event/verbose 0
/tracking/verbose 1

## Particle type and energy
#/gun/particle gamma
#/gun/energy 1.0 MeV

## Number of events to run (only few!)
/run/beamOn 2

Deposited-energy spectra

Using the macro spec-eff.mac we can obtain the plots of the energy deposited in the detector for the particles defined in the testPhysicsList, and loading the batch file hplot.cc in ROOT we can visualize the histograms, as we did before.

# Execution macro
#    /control/execute spec-eff.mac

## Particle type and energy
/gun/particle gamma
#/gun/energy 5.0 MeV

## Number of events to run
/run/beamOn 10000

We can get the theta efficiency from the file test.eff.

Efficiency for test detector
 Total efficiency = 1367

Histogram plotting the energy vs the counts Gammas at 1 MeV in CsI.

Histogram plotting the energy vs the counts Electrons at 1 MeV in CsI.

Histogram plotting the energy vs the counts Positrons at 1 MeV in CsI.

Want some resources?

Then, check this list out:

Alex Bielajew, Fundamentals of the Monte Carlo method for neutral and charged particle transport, http://www-personal.engin.umich.edu/~bielajew/MCBook/book.pdf
The Geant4 official site http://geant4.web.cern.ch
The ROOT official site https://root.cern.ch

Published: 2006-04-29 22:28:31