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.




Compiling Fortran-77 code in Ubuntu-16

Leave a comment

Fortran-77 is a very old code, who lives in 32-bit computer.

In Ubuntu-16, the compiler g++, gcc, or gfortran are “basically the same” (as far as I understand, correct me if I am wrong.) that they only support fortran-95.

In order to compile Fortran-77 code, I tried many way, but the only way is install g77 from external source, and add -m32 for the compiling flag.

The g77 compiler can be downloaded in here (I download from the web, If I violated some copy right, please let me know):

or, people can search in google by



people need to extract, change the mod of to be executable.

 tar -xzvf g77_x64_debian_and_ubuntu.tar.gz
 cd g77_x64_debian_and_ubuntu
 chmod +x ./


Somehow, you may face an error in apt-get, saying

Errors were encountered while processing:

you can remove that by

cd /var/lib/dpkg/info
sudo rm g77-3.4-doc*
sudo dpkg --remove --force-remove-reinstreq g77-3.4-doc

Thanks (here)

Hope it help. :)


Algorithm of Wavelet Transform (with Qt class)

Leave a comment

There are many kind of wavelet transform, and I think the names are quite confusing.

For instance, there are continuous and discrete wavelet transforms, in which, the “continuous” and “discrete” are for the wavelet parameters, not for the “data” itself. Therefore, for discrete data, there are “continuous” and “discrete” wavelet transforms, and for function, there are also “continuous” and “discrete” wavelet transforms.

In here, we will focus on discrete wavelet transform for function first. This discrete wavelet transform is also called as wavelet series, which express a compact support function into series of wavelet.

For simplicity, we also focus on orthonormal wavelet.

As the wavelet span the entire space, any compact function can be expressed as

\displaystyle f(t) = \sum_{j,k} \left<f(t)|\psi_{j,k}(t)\right> \psi_{j,k}(t)

\psi_{j,k}(t) = 2^{j/2} \psi(2^j t - k)

where j, k are integer.

Now, we move to discrete data discrete wavelet transform. The data is discrete, we can imagine only t_n = t_0 + n \Delta points are known with finite n .

\displaystyle f_n = f(t_n) = \sum_{j,k} \left<f_n|\psi_{j,k}(t_n) \right> \psi_{j,k}(t_n)

the integration becomes a finite sum.

Without loss of generality, we can set t_0 = 0, \Delta = 1, and then the time axis becomes an integer number axis. We found that j  \leq 0 as the wavelet can only be expand, not shrink. Because there are finite number of data point, i.e. n < \infty, -Log_2(n) < j \leq 0 .

However, this double summation for each f_n is very time consuming. There is a Fast Discrete Wavelet Transform. Before we continuous, we must study the wavelet.

From the last post, we know that the scaling function that generate a MRA must be:

\displaystyle \phi(t) = \sum_{k} g_0(k) \phi(2t-k)

\left<\phi(t-k) | \phi(t-k') \right> = \delta_{kk'}

, where k are integer. The set of shifted scaling function span a space V_0 . For the wavelet,

\displaystyle \psi(t) = \sum_{k} g_1(k) \psi(2t-k)

\left<\psi(t-k) | \psi(t-k') \right> = \delta_{kk'}

The set of shifted wavelet span a space W_0, so that W_0 \perp V_0, so that

\left<\phi(t-k)|\psi(t-k') \right> = 0

Since the wavelet is generated from the scaling function, we expect the coefficient of g_0(k) and g_1(k) are related. In fact, the relationship for orthonormal scaling function and wavelet is

g_1(k) = (-1)^k g_0(1-k)

For discrete data x_i , it can be decomposed into the MRA space. We start by the largest V_0 space, where the wavelet is most shrunken.

\displaystyle x_i = \sum_{k} v_{0,k} \phi(i-k)

to decompose to the V_{-1} and W_{-1} space. We can use the nested property of the MRA space, \phi(2t) can be decomposed into \phi(t-k) and \psi(t-k) ,

\displaystyle \psi(2t-l) = \sum_{k} h_0(2k-l) \phi(t-k) + h_1(2k-l) \psi(t-k)

where (given that \phi(t) and $\latex \psi(t)$ are orthonormal ),

h_0(2k-l) = \left< \phi(2t-l) | \phi(t-k) \right>

h_1(2k-l) = \left< \phi(2t-l) | \psi(t-k) \right>

Therefore, using the coefficient of h_0 and h_1, the wavelet coefficient v_{0,k} can be decomposed to

\displaystyle v_{s-1,k} = \sum_{l} h_0(2k-l) v_{s,l}  

\displaystyle w_{s-1,k} = \sum_{l} h_1(2k-l) v_{s,l}

in graphic representation


This is a fast discrete wavelet transform.

Due to the nested space of MRA, we also expect that the coefficient h_0 and h_1 are related to g_0 . For orthonormal wavelet,

\displaystyle h_0(k) = \frac{1}{2} g_0(-k)

\displaystyle h_1(k) = \frac{1}{2} (-1)^{k} g_0 (k+1)

Since the g_0 is finite, the g_1, h_0, h_1 are all finite. That greatly reduce the computation cost of the discrete wavelet transform.

To reconstruct the discrete data x_i, we don’t need to use

\displaystyle v_{s+1,l} = \sum_{k} v_{s,k} \phi(l - k) + w_{s,k} \psi(l-k)

using the nested space of MRA, \psi(t) = \sum_{k} g_1(k) \psi(2t-k) ,

\displaystyle v_{s+1,l} = \sum_{k} g_0(l-2k) v_{s,k} + g_1(l-2k) w_{s,k}

in graphical representation,


I attached the wavelet transfrom class for Qt, feel free to modify.

in the code, the data did not transform to MRA space. The code treats the data already in the MRA space. Some people said this is a “crime”. But for the seek of “speed”, it is no need to map the original discrete data into MRA space. But i agree, for continuous function, we must map to MRA space.



Drawing energy level with Latex

Leave a comment

First, make sure you have the tikz, so that you can


without any error.

in linux,

sudo apt-get install texlive-pictures pgf

with the tikz package, we can use \draw to draw line or arrow, \node to draw a text. I will demonstrate to draw a simple energy levels scheme. The advantage of using latex is that the energy levels can be drawn accurately and once the template is set, it is very easy.

I define new command

\draw [level](#3,#1) -- (#3+\len,#1) ;
\draw (#3+\len +0.1,#1) -- (#3+\len +0.3, #1 + #4) -- (#3+\len +0.6, #1 + #4);
\node[right] at (#3+\len +0.6,#1+#4) {#1, #2};

with usage,


with this, we can draw something like this.


I attached the template in here.


Maximum Likelihood

Leave a comment

In data analysis, especially the number of data is small, in order to found out the parameter of the distribution, which fit the data the best, maximum likelihood method is a mathematical tool to do so.

The ideal can be found in Wikipedia. For illustration, I generate 100 data points from a Gaussian distribution with mean = 1, and sigma = 2.

In Mathematica,

Data = RandomVariate[NormalDistribution[1, 2], 100]
MaxLikeliHood = Flatten[Table[{
 Log[Product[PDF[NormalDistribution[mean, sigma], Data[[i]]] // N, {i, 1,100}]],
 {mean, -3, 3, 1}, {sigma, 0.5, 3.5, 0.5}], 1]

This calculate the a table of mean form -3 to 3, step 1, sigma from 0.5 to 3.5 step 0.5. To find the maximum of the LogProduct in the table:

maxL=MaxLikeliHood[[1 ;; -1, 3]];
maxN = Position[maxL[[1 ;; -1, 3]], %]

The result is


which is the correct mean and sigma.


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.


Older Entries