GPS energy vs. transferred momentum

Particles, Track, Event, Run and Biasing
1

Dear experts,

I have a question about the General Particle Source, where I do not understand the mismatch between the requested energy via “/gps/ene/mono” and the transferred momentum of the primary particle. The transferred momentum is equal to the energy for massless particles, but for massive particles the momentum is larger than the energy.

Using the Geant4 “hadrontherapy” example, I can reproduce this by setting the verbose levels accordingly:

/gps/verbose 2
/event/verbose 2
I get the following output for one example event:

…
Particle name: proton
Energy: 61.8414
…
Primary Vetex generated !

G4EventManager::ProcessOneEvent()

1 vertices passed from G4Event.
Primary particle (proton) — Transfered with momentum (346.226,-0,0)
…

Thank you for your help!
Best
Johannes

2

Hi jerdmann,

I consider as intermediate but can slightly help you with gps selection and momentum storing. I am not getting what exactly you doing with gps and your geometry. In simple terms, gps/ene/mono use when defines the profile of the beam. I mean if you choosing rectangular beam or Gaussian or etc. Better you change and try with simple case
/gps/particle proton
/gps/number 1
/gps/energy 61.8414 keV

Regarding momentum I suggest you to call in ProcessHits or Hits class (in case sensitive detectors) defining accumulation of G4step with momentum in EventAction for Event by Event.

3

Thanks for your answer. Indeed, I can reproduce the mismatch if I use /gps/energy instead and create a mono-energetic beam.

To be clear: my question is independent of geometry, i.e. I wonder if the transferred momentum if the beam particle is correctly calculated from the energy provided to gps for a massive particle.

4

Hi Jerdmann,

confine gps with physical volume of material G4_Galactic with PV name as Source.
/gps/pos/confine Source

if parentID == 0 && stepVolume == Source
{
G4double pX_part = prePoint->GetMomentum().x();
G4double pY_part = prePoint->GetMomentum().y();
G4double pZ_part = prePoint->GetMomentum().z();
}

I hope it will work for sure.

VRS

5

Hi,

thanks a lot for following up!

  1. I now defined a volume that is definitely large enough to contain the source:

G4VSolid* sourceBox = new G4Box(“SourceBox”, 10m/2., 10m/2., 10*m/2.);

G4LogicalVolume* sourceLV = new G4LogicalVolume(sourceBox, G4Material::GetMaterial(“Galactic”), “SourceLV”);

//G4LogicalVolume* sourceLV = new G4LogicalVolume(sourceBox, G4Material::GetMaterial(“G4_Pb”), “SourceLV”);

G4VPhysicalVolume* sourcePV = new G4PVPlacement(0, G4ThreeVector(0,0,0), sourceLV, “SourcePV”, worldLV, false, 0, true);

  1. I also added "/gps/pos/confine SourcePV” to the run file.

  2. With “/gps/ene/mono 10000 MeV” I still see the debug output:

Primary particle (proton) — Transfered with momentum (0,0,10898)

  1. I also implemented a UserSteppingAction, which however does not contain any steps for material “Galactic” (as expected) but contains steps if I change the material of the SourcePV to G4_Pb (just for testing purposes). This is the output I get:

=====================================

G4EventManager::ProcessOneEvent()

=====================================

1 vertices passed from G4Event.

Primary particle (proton) — Transfered with momentum (0,0,10898)

A new track 0x14c978ac0 (trackID 1, parentID 0) is passed to G4StackManager.

1 primaries are passed from G4EventTransformer.

!!! Now start processing an event !!!

###pop requested out of 1 stacked tracks.

Track 0x14c978ac0 (trackID 1, parentID 0) is passed to G4TrackingManager.

parentID = 0, trackVolume = SourcePV, prePoint x/y/z-momentum: 0/0/10898

parentID = 0, trackVolume = SourcePV, prePoint x/y/z-momentum: 20.2309/8.21415/10889.5

parentID = 0, trackVolume = SourcePV, prePoint x/y/z-momentum: 17.0722/15.7584/10884.5

…

So, it seems that this does not fix the issue, i.e. the z-momentum is still > 10 GeV for the proton. Did I miss something?

Best

Johannes

6

In case if the direction of the particle is ‘x’ below is a code that can help you to understand. You need to define ‘prePoint’ using G4StepPoint*.
After you can call it with different methods under ProcessHits. e.g.
G4cout << “K.E. PDMass " << partM << " is " << prePoint->GetKineticEnergy() << " MeV, and Momentum: " << prePoint->GetMomentum().x() << " MeV/c”<< std::endl;

Result for e- of energy 2 MeV
K.E. PDMass 0.510999 is 2 MeV, and Momentum: 2.45845 MeV/c

Mathematical Analysis to prove momentum calculation by geant4 is 100% correct.
we know: mc^2=0.511MeV
pc=sqrt(K^2+2Kmc^2) ​= sqrt((2.00MeV)^2+2(2.00MeV)(0.511MeV))
This yields p = 2.46MeV/c

Now your question: your first post, the answer is below;
K.E. PDMass 938.272 is 61.8414 MeV, and Momentum: 346.226 MeV/c

:slight_smile:

I hope you caught me. Good Luck.

7

I was missing that /gps/ene/mono defines the kinetic energy, not the total energy. Thanks for your help!

8

Any particle that moves consist of kinetic energy considered as its total energy.