Using TGenPhaseSpace to calculate nuclear reaction

Leave a comment

In CERN Root analysis software,  there is a class called TGenPhaseSpace. This class can generate all possible solution of a nuclear reaction ( elastic, inelastic, 2-body, 3-body, decay, etc. ) base on isotropic distribution and balancing four-momenta.

The unit of mass, energy, and momentum are GeV/c^2, GeV, and GeV/c, respectively. The typical usage is

TLorentzVector target; 
TLorentzVector beam; 

Double_t masses[3] = { m1, m2, m3}; // what particles after reaction

TGenPhaseSpace event;
event.SetDecay(target + beam, 3, masses);

for( Int_t n = 0; n < numEvent; n++){
   Double_t xsec = event.Generate();       // cross-section
   TLorentzVector *p1 = event.GetDecay(0); // get the 4-momentum of m1
   TLorentzVector *p2 = event.GetDecay(1);
   TLorentzVector *p3 = event.GetDecay(2); 

   // fill histogram;


In the proton-proton elastic reaction at 300 MeV,  the calculation is very nice.

However, in the 23F(p,2p)22O inverse knockout reaction at 300A MeV , the Fermi-momentum of the bound proton, which is given by

P_F - P_O = P_k ,

where P_F is the 4-momentum of the 23F, P_F is the 4-momentum of the 22O, can be very large.

The Fermi-energy of a bound nucleon has maximum 40 MeV, or the momentum can be at most 400 MeV/c. But the TGenPhaseSpace only calculate all possible solution that balancing the 4-momenta, the Fermi-momentum can be more than 1 GeV/c.


This gives an un-realistic result. For example, the proton KE can be as large as 1200 MeV, while only 300A MeV of 23F KE. In realistic situation, the maximum KE of the scattered proton is at most 300 MeV at zero degree.

Now, if we restricted the Fermi-momentum to be 100 MeV/c +- 10 MeV/c, we have,


This is more realistic result.

If we can control the distribution of Fermi-momentum, and also understand the estimation of the cross-section, we can use it to simulate many reactions.





a GEANT4 Simulation


The GEANT4 program structure was borrow from the example B5. I found that the most confusing part is the Action. Before that, let me start with the main().

GEANT4 is a set of library and toolkit, so, to use it, basically, you add a alot GEANT4 header files on your c++ programs. And, every c++ program started with main(). I write the tree diagram for my simplified exampleB5,

 +- DetectorConstruction.hh
    +- Construct()
       +- HodoscopeSD.hh     // SD = sensitive detector
          +- HodoscopeHit.hh //information for a hit
          +- ProcessHits()   //save the hit information
       +- ConstructSDandField() //define which solid is SD
       +- ConstructMaterials()  //define material
+- PhysicsList.hh  // use FTFP_BERT, a default physics for high energy physics
+- ActionInitialization.hh
   +- PrimaryGeneratorAction.hh // define the particle gun
   +- RunAction.hh // define what to do when a Run start, like define root tree
   +- Analysis.h  // call for g4root.h
   +- EventAction.hh //fill the tree and show some information during an event


A GEANT4 program contains 3 parts: Detector Construction, Physics, and Action. The detector construction is very straight forward. GEANT4 manual is very good place to start. The physics is a kind of mystery for me. The Action part can be complicated, because there could be many things to do, like getting the 2ndary beam, the particle type, the reaction channel, energy deposit, etc…

Anyway, I managed to do a simple scattering simulation, 6Li(2mm) + 22Ne(60A MeV) scattering in vacuum. A 100 events were drawn. The detector is a 2 layers plastic hodoscope, 1 mm for dE detector, 5 mm for E detector. I generated 1 million events. The default color code is Blue for positive charge, Green for neutral, Red for negative charge. So, the green rays could be gamma or neutron. The red rays could be positron, because it passed through the dE detector.

Screenshot from 2016-01-30 00:34:34.png

The histogram for the dE-TOF isScreenshot from 2016-01-29 23:26:34.png

We can see a tiny spot on (3.15,140), this is the elastics scattered beam, which is 20Ne. We can see 11 loci, started from Na on the top, to H at the very bottom.

The histogram of dE-E

Screenshot from 2016-01-29 23:30:38.png

For Mass identification, I follow the paper T. Shimoda et al., Nucl. Instrum. Methods 165, 261 (1979).

Screenshot from 2016-01-30 00:06:02.png

I counted the 20Na from 0.1 billion beam, the cross section is 2.15 barn.


GEANT4 installation – Concept

Leave a comment

for the moment, i think GEANT4 is a set of libraries (of physical process, meterial) and a tool (Monte Carlo method).

As you can see, Visualization is not included in GEANT4. I mean, GEANT4 can output some files for other Visualization programs, but those programs have to be install separately.

A complete environment = GEANT4 + Visualization 

//——————– GEANT4 installation

To install GEANT4, you can go to to download to source code.

The installation method can be found in

basically, unzip the source code.

tar -zxvf geant4-source
mkdir geant4-build
cd geant4-build
cmake -DGEANT4_INSTALL_DATA=ON -DGEANT4_USE_NETWORKVRML=ON -DCMAKE_INSTALL_PREFIX=/path/to/geant4-install  geant4-source
make install


There are two options, 1 is the GEANT4_INSTALL_DATA, this will download and install the physical process for radiation-matter interaction. 2 is the GEANT4_USE_NETWORKVRML, this will enable network VRML visuallization. This is a local install, all the installed files are located at /path/to/geant4-install. After the installation, better to write this in .bashrc

cd ~/geant4.10.02-install/bin
cd ~/


The only difficulty is the CMAKE 3.0 or above. But all errors i met can be easily found in google.

//————— Visualization program.

I skipped the openGL, as it is well known for unfriendly for NVidia display card.

I tried DAWN and freeWRL. DAWN is a static visualization. It read *.prim files and output an *.eps file for fixed angle. freeWRL is a VRML display using java-script, an alternative is openVRML. But I cannot install openVRML.

DAWN can be installed from, I have zero problem.

FreeWRL is a bit tricky. I can only install the version, not earlier, not later. I don’t know why. It require javascript engine, apt-get install libmozjs185-dev. when it said any thing missing, just look for libXXXX.

//——————————– Example.

I took the example form the geant4-source/example/basic/B1, copy it to ~/geant4Works/B1.

mkdir B1-build
cd B1-build
cmake -DGeant4_DIR=/path/to/geant4-install/lib/Geant4-10.2.0  B1

to run, you simple ./exampleB1

the macro files are *.mac. vis.mac is visualization macro, and will be called by init_vis.mac. You should go to vis.mac and modify it.

To run,

Idle> /control/execute run1.mac


I modified the last line be /run/beamOn 100, the result is

Screenshot from 2016-01-27 13:46:48.png

you can see there are 100 protons passing through. I don’t understand the example B1, I’m just showing you what is a proper result.