# 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 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);

G4LogicalVolume* logicDetector   = new G4LogicalVolume(solidDetector, Ge, "Detector");
G4VPhysicalVolume* physiDetector = new G4PVPlacement(0, positionDetector, "Detector", logicDetector, physiWorld, false, 0);``````  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. 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

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

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

}``````

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.)
}
}``````
``````void testEventAction::EndOfEventAction(const G4Event* evt) {

if (TotalEnergyDeposit > 0.) {

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

}
}``````
``````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

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

/tracking/storeTrajectory 1
/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``````  The figures above show some processes happening already, like pair creation for example.  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`````` Gammas at 1 MeV in CsI. Electrons at 1 MeV in CsI. Positrons at 1 MeV in CsI.

## Want some resources?

Then, check this list out:

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