The unit of getting Total Kinetic Energy/wrong GetTotalEnergy()

Particles, Track, Event, Run and Biasing
1

Hello,
I am a beginner in GEANT4, i wrote the following code in my detector.cc file. (The Headers are given in the detector.hh file):

#include "detector.hh"

sensitivedetector::sensitivedetector(G4String name) : G4VSensitiveDetector(name){}
sensitivedetector::~sensitivedetector(){}


G4bool sensitivedetector:: ProcessHits(G4Step *aStep, G4TouchableHistory *ROhist){

    G4Track* track = aStep->GetTrack();
    G4ParticleDefinition* particle = track->GetDefinition();
 
    if(particle==G4Alpha::Definition()){
        // G4double edep= aStep->GetTotalEnergyDeposit();
        G4double incene = track->GetTotalEnergy();
        // G4double time = aStep->GetPreStepPoint()->GetGlobalTime();
        std::ofstream file("output.txt", std::ios::app);
        file.seekp(0, std::ios::end);
        file <<incene<<G4endl;
        file.close();
    }  

}

I am running a simulation of radioactive decay of Am^241, i am using a simple cylindrical detector to measure the energy of the produced alpha particles. theoretically the obtained energy should fall in range of 5.4MeV but what i am getting in output.txt file is:

3727.38
3727.38
3727.38
3727.38
3727.38
3727.38
3727.38
3727.38
3727.38
3727.38
3727.38
3727.38
3727.38
3727.38
3727.38
3727.38

I believe that the output’s unit is MeV, if So, how is the Energy coming in GeV? Am i doing anything wrong? Please let me know if theres any fix for this.

My detector Construction is :

#include "construction.hh"

detectorconstruction::detectorconstruction(){}

detectorconstruction::~detectorconstruction(){}


G4VPhysicalVolume *detectorconstruction::Construct(){

    // Define materials
    G4NistManager* nistManager = G4NistManager::Instance();
    G4Material* air = nistManager->FindOrBuildMaterial("G4_AIR");
    G4Material* silicon = nistManager->FindOrBuildMaterial("G4_Si");


    G4VisAttributes* WW = new G4VisAttributes(G4Colour(0.85,0.85,1,0));
    WW->SetForceSolid(true);

    G4VisAttributes* detectorvis = new G4VisAttributes(G4Colour(0.85,0.85,0.45,0.5));
    detectorvis->SetForceSolid(true);


    G4bool checkoverlap = true;

    //Defining World Volume
    G4Box *solidworld = new G4Box("world",2.5*m,2.5*m,2.5*m);
    G4LogicalVolume *logicworld = new G4LogicalVolume(solidworld, air,"logicalworld");
    G4VPhysicalVolume *physicalworld = new G4PVPlacement(nullptr,G4ThreeVector(0,0,0*km),logicworld, "physicalworld", nullptr, false,0,checkoverlap);
    logicworld->SetVisAttributes(WW);




    G4Tubs* alpha_detector = new G4Tubs("Alpha detector", 0.01*mm,4*mm,1*cm,0,360*deg);
    logicdetector = new G4LogicalVolume(alpha_detector,silicon, "alpha_detector");
    new G4PVPlacement(nullptr, G4ThreeVector(0.,0.,1.5*cm), logicdetector, "alpha_detector", logicworld, false, 0,checkoverlap);


    logicdetector->SetVisAttributes(detectorvis);
    return physicalworld;



}



void detectorconstruction::ConstructSDandField(){
    sensitivedetector *sensdet = new sensitivedetector("SD");
    logicdetector->SetSensitiveDetector(sensdet);
}

My Generator is defined as:

#include "generator.hh"
#include <cmath>




primarygenerator::primarygenerator(){

    fparticlegun = new G4ParticleGun(1);

    G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
    G4ParticleDefinition* particle = particleTable->FindParticle("geantino"); 




    fparticlegun->SetParticleDefinition(G4MuonMinus::Definition());  
    G4ThreeVector pos(0.*mm,0.0*mm,0.0*m);
    G4ThreeVector mom(0,0,0);
    fparticlegun->SetParticlePosition(pos);
    fparticlegun->SetParticleMomentumDirection(mom);
    fparticlegun->SetParticleMomentum(0*GeV);
     
}

primarygenerator::~primarygenerator(){
    delete fparticlegun;
}


void primarygenerator::GeneratePrimaries(G4Event *anEvent)
{
    G4ParticleDefinition* particle = fparticlegun->GetParticleDefinition();
    if(particle==G4Geantino::Geantino());
    {
        G4int Z = 95;
        G4int A = 241;
        G4double Charge   = 0.*eplus;
        G4double Energy  = 0.*keV;
        G4ParticleDefinition* ion = G4IonTable::GetIonTable()->GetIon(Z,A,Energy);
        fparticlegun->SetParticleDefinition(ion);
        fparticlegun->SetParticleCharge(Charge);

    }
    

    fparticlegun->GeneratePrimaryVertex(anEvent);
}

Almost all the codes that are written by me are taken from this YouTube Channel Tutorial on Geant4 by Physics Matters by a guy named Mustafa.

Geant4 Version: 11.1.2
Operating System: Linux Mint 21 Cinnamon
Compiler/Version: g++ & gcc 11.4
CMake Version: 3.22.1

2

Total energy includes mass. The mass of an alpha particle is 3727 MeV/c^2

You want kinetic energy

3

Yeah, I changed the part where it reads the Total Energy to Kinetic energy like:

#include "detector.hh"

sensitivedetector::sensitivedetector(G4String name) : G4VSensitiveDetector(name){}
sensitivedetector::~sensitivedetector(){}


G4bool sensitivedetector:: ProcessHits(G4Step *aStep, G4TouchableHistory *ROhist){

    G4Track* track = aStep->GetTrack();
    G4ParticleDefinition* particle = track->GetDefinition();
 
    if(particle==G4Alpha::Definition()){
        // G4double edep= aStep->GetTotalEnergyDeposit();
        G4double incene = track->GetKineticEnergy();
        // G4double time = aStep->GetPreStepPoint()->GetGlobalTime();
        std::ofstream file("output.txt", std::ios::app);
        file.seekp(0, std::ios::end);
        file <<incene<<G4endl;
        file.close();
    }  

}

But Now, in my output file, the printed out values are all 0, and rarely, some of the entries are of few keV (Theorectically, majority of it should be around 5.5 MeV). Any idea why is that?

4

You are tracking every step of every alpha particle. So even steps in which the kinetic energy may be 0 will be printed. Also, the default units for geant4 are MeV. If you want to print keV, you need to have

file << incene*keV << G4endl;

If you want to look at the energy of the alpha particles deposited in your detector, you need to use GettotalEnergyDeposit()

5

Not exactly. If you want to print values in keV, you divide by the unit:

G4cout << incene/keV << " keV" << G4endl;
6

This topic was automatically closed 7 days after the last reply. New replies are no longer allowed.