Abstract
Even though a computer-aided design (CAD)-based geometry can be directly implemented in Geant4 as polygon-mesh using the G4TessellatedSolid class, the computation speed becomes very slow, especially when the geometry is composed of a large number of facets. To address this problem, in the present study, a new Geant4 solid class, named DagSolid, was developed based on the direct accelerated geometry for the Monte Carlo (DAGMC) library which provides the ray-tracing acceleration algorithm functions. To develop the DagSolid class, the new solid class was derived from the G4VSolid class, and its ray-tracing functions were linked to the corresponding functions of the DAGMC library. The results of this study show that the use of the DagSolid class drastically improves the computation speed. The improvement was more significant when there were more facets, meaning that the DagSolid class can be used more effectively for complicated geometries with many facets than for simple geometries. The maximum difference of computation speed was 1562 and 680 times for Geantino and ChargedGeantino, respectively. For real particles (gammas, electrons, neutrons, and protons), the difference of computation speed was less significant, but still was within the range of 53–685 times depending on the type of beam particles simulated.
Export citation and abstract BibTeX RIS
General scientific summary Currently, there is focus on direct implementation of CAD-based geometries to various general-purpose Monte Carlo particle transport codes for sophisticated geometry modeling in various applications such as nuclear, medical, and space science. However, a considerable amount of computation time for CAD-based geometries is required due to CAD-based geometries being generally composed of a large number of facets. It is, therefore, important to study how to reduce the computation time for CAD-based geometries. In the present study, we developed a new Geant4 class, named DagSolid used for fast Monte Carlo simulation with CAD-based geometries in Geant4 Monte Carlo code. The DagSolid class was developed based on the DAGMC library which provides an accelerated ray-tracing algorithm. The present study shows that the DagSolid class can improve computation speeds with the range of 53–685 times depending on particle-types giving correct dose values as compared with an existing class in Geant4.
1. Introduction
Computer-aided design (CAD) is popularly used in industries to design and manufacture physical products such as buildings, bridges, roads, aircrafts, ships, and cars. Several investigators (Schwarz et al 2005, Badal et al 2009, Tautges et al 2009, Wu and Team 2009, Constantin et al 2010, Sato et al 2010) have implemented CAD-based geometries in general-purpose Monte Carlo codes such as Geant4 (Agostinelli et al 2003), MCNP/MCNPX (Briesmeister 2000), and PENELOPE (Salvat et al 2006).
The polygon-mesh is typically used in 3D computer graphics and can be used to model highly complex geometries such as human organs, medical equipment, and brachytherapy applicator modeling to a great detail. Although a polygon-mesh geometry can contribute to very accurate modeling of geometry, it usually requires a large number of facets, resulting in very long computation time. For instance, Kim et al (2011) developed a very detailed polygon-mesh based human phantom from a voxel based phantom for more precise dose calculation. Their study, however, reported that the polygon-mesh phantom was 70–150 times slower than the original voxel based phantom in Monte Carlo dose calculations.
The slowdown of computation speed is mainly due to the repeated use of the computational expensive ray-tracing functions. The ray-tracing functions are repeatedly called for each of the facets around the particle whenever a particle is to move, and dominate the total computation time. This repetitive work proportionally increases with the number of the facets used in the geometry.
Recently, a research team at the University of Wisconsin developed a ray-tracing acceleration technique called direct accelerated geometry for Monte Carlo (DAGMC), and they integrated the technique into MCNP5 (Tautges et al 2009, Wilson et al 2010, Smith and Wilson 2011). The DAGMC is based on a hierarchical structure of oriented bounding box (OBB) tree for rapid interference detection. Badal et al (2009) developed a modular program named penMesh to implement triangle mesh geometries in PENELOPE. The penMesh also includes a ray-tracing acceleration algorithm which is based on Octree hierarchical structure.
The ray-tracing acceleration technique, however, has not been implemented in Geant4, and the objective of the present study is to develop a new Geant4 solid class based on the ray-tracing acceleration algorithm functions of the DAGMC library to improve computation speed for polygon-mesh geometries. The performance and accuracy of the developed DagSolid class were evaluated by comparing the computation time and calculation result with those of the conventional G4TessellatedSolid class for identical simulation cases.
2. Materials and methods
2.1. Direct accelerated geometry for Monte Carlo (DAGMC)
To improve the computation speed in a polygon-mesh geometry in Geant4, this study developed a new class, named DagSolid, based on the ray-tracing functions of the DAGMC library. The DAGMC library is again based on the mesh-oriented datABase (MOAB) library for more generic functionality such as tree construction and searching, and the common geometry module (CGM) library for importing and evaluating polygon-mesh geometries (Tautges et al 2004, Tautges 2001). The DAGMC is designed to be integrated with Monte Carlo particle transport codes. DAG-MCNP5, integrating the DAGMC with the MCNP5 code, has been developed and used to analyze modified international thermonuclear experimental reactor first wall/shield designs. Lately, DAG-Tripoli4, integrating the DAGMC with the Tripoli4 code, has been also developed.
The main idea of the DAGMC library for ray-tracing in a polygon-mesh geometry is to use the OBB tree structure for a very fast hierarchical search of the facet that is intersected by the projected particle motion, which is modeled by a straight ray. This kind of acceleration technique is frequently used in computer graphics for collision detection.
Figure 1 shows the schematic diagram of a hierarchical OBB tree. First, the entire surface, represented by a set of facets, is surrounded by a bounding box called root bounding box (figure 1(a)). The set of facets is then geometrically divided into two subsets of facets, and the subsets are individually surrounded by bounding boxes (figures 1(b), (c)). Bounding boxes are continuously formed until there is a bounding box for a single facet (figure 1(d)). As a result, a binary OBB tree is built, which is then used for hierarchically search of the bounding box that includes the facet collided by extension of the particle track. This hierarchical search process is much faster than the regular ray-tracing process in which all of the facets are individually checked if it is collided by the particle track. Use of the binary OBB tree theoretically decreases the number of repeated calculations to find the collided facet from N (for the regular ray-tracing algorithm) to log2N, where N is the number of the facets around the particle.
Figure 1. Schematic diagram of hierarchical OBBs ((a)–(d)) in DAGMC and example of finding a facet collided with particle track (e).
Download figure:
Standard image High-resolution image2.2. Solid class in Geant4
In Geant4, all solid classes inherit from the G4VSolid class or one of its derived classes. The G4TessellatedSolid class used for polygon-mesh geometries also inherits from the G4VSolid class. The G4VSolid class has mandatory member functions and some functions are used for ray-tracing; that is, Inside, DistanceToIn, DistanceToOut, SurfaceNormal, and CalculateExtent functions, which are used to determine the particle step length as follows:
- (1)Geant4 finds the solid in which the particle is positioned by using the Inside function of the solid class. Then, the '/physical step-length' is calculated by the GetPhysicalInteraction-Length (GPIL) function based on the material information of the solid.
- (2)The calculated 'physical-step-length' is compared with the isotropic distance to the nearest boundary in geometry called 'safety'. The 'safety' is calculated by the DistanceToIn (p) and DistanceToOut (p) function of the solid classes for the mother and daughter volumes.
- (3)If the 'physical-step-length' is shorter than the 'safety,' the "physical-step-length' is selected as a particle step length without any geometric calculations.
- (4)In the opposite case, the DistanceToIn (p, v) or DistanceToOut (p, v) function calculates the 'geometric-step-length,' which is the distance, along the particle direction, to a facet from the point p. Then, the smaller of the 'physical-step-length' and the 'geometric-step-length' is selected as a particle step length.
The SurfaceNormal and CalculateExtent functions are selectively used as necessary during the aforementioned procedure. The SurfaceNormal function returns a unit normal vector of the facet closest to the particle, and the CalculateExtent function calculates the minimum bounding box for the solid. The procedure is repeated whenever a particle is to move, resulting in a significant of computation time, typically occupying most computation time during Monte Carlo simulation.
2.3. Development of DagSolid class
In the present study, a new Geant4 solid class, named DagSolid, was developed based on the G4VSolid class and the ray-tracing acceleration algorithms of the DAGMC library. The DAGMC library was linked to Geant4 for the code to use the DAGMC functions in this solid's implementation.
The DAGMC library has functions corresponding to the previously explained functions for determination of particle step length in the G4VSolid class except for CalculateExtent function corresponding to get_coords function included in MOAB library. Figure 2 shows the corresponding relationship of functions between the DAGMC library and the G4VSolid class. To develop the DagSolid class, the new Geant4 solid class was derived from the G4VSolid class and then the ray-tracing functions of the class were linked to the ray-tracing functions of the DAGMC library.
Figure 2. Corresponding relationship between functions of DAGMC and G4VSolid class.
Download figure:
Standard image High-resolution image2.4. Monte Carlo simulation set-up
We performed Monte Carlo simulations to evaluate the performance of the developed DagSolid class. The computation times of the DagSolid class and the G4TessellatedSolid class were compared under identical simulation conditions. For this, two types of geometries were considered; (1) a sphere polygon-mesh model, representing a simple geometry and (2) a polygon mesh model of a small intestine, representing a complex geometry.
The sphere polygon-mesh model with a diameter of 1 m was constructed by Rapidform® XOS/SCAN™. The number of polygon facets of the sphere was 5000, 20 000, 80 000, and 320 000 as shown in figure 3. The polygon-mesh model of the small intestine was extracted from the PSRK-Man developed by Kim et al (2011). The size of bounding box of the small intestine was 18.97 cm × 11.16 cm × 25.59 cm. The numbers of polygon facets of the small intestine were also equal to those of the sphere.
Figure 3. Polygon mesh sphere model with different number of facets, (a) 5000 facets, (b) 20 000 facets, (c) 80 000 facets, (d) 320 000 facets.
Download figure:
Standard image High-resolution imageFirst, in order to emphasize the effect of the ray-tracing acceleration algorithm on particle transport process, the Geant4 virtual particles, Geantino and ChargedGeantino, were simulated. Then, real particles (gamma, electron, neutron, and proton) were simulated to see the effect in real radiation transport situations.
The Geantino and ChargedGeantino are virtual particles and do not interact with the medium. The Geantino has no charge, but the ChargedGeantino has charge and therefore interacts with the magnetic field. The movement of the ChargedGeantino depends on its energy and the direction and strength of the magnetic field. In the present study, the energy of ChargedGeantinos and the strength of the magnetic field were assumed to be 10 MeV and 3 T, respectively.
The particles were generated randomly and isotropically within the bounding box of the solids. The number of the simulated primary particles was 104. Figure 4 shows the movement of Geantinos and ChargedGeantinos in the sphere and small intestine models.
Figure 4. Tracks of Geantinos and ChargedGeantino in sphere and small intestine: (a) Geantinos in sphere, (b) Geantinos in small intestine, (c) ChargedGeantinos with sphere, and (d) ChargedGeantinos in small intestine.
Download figure:
Standard image High-resolution imageTo implement the polygon-mesh geometries in Geant4 with the DagSolid class, the HDF5 (*.h5m, binary file) file format must be used. Because the sphere and small intestine polygon mesh models were first obtained as a 'obj' file, which is one of the popular formats for polygon mesh geometries, a class, named Obj2h5m, was developed separately to convert the 'obj' files to the 'h5m' files.
The computation time of each ray-tracing function and total computation time were measured by Intel® VTune™ Amplifier XE 2011. The simulation in the present study was performed on a single core of the AMD Opteron™ 6176 @ 2.3 GHz and 64 GB memory with Geant4 9.5 patch 01, MOAB 4.0, DAGMC 0.99, and CGM 10.2.3.
3. Results and discussion
3.1. Performance evaluation for virtual particles
The performance of the developed DagSolid class was evaluated by comparing the computation time with that of the conventional G4TessellatedSolid class for identical simulation cases. For this, in the present study, we measured the total computation time and the class computation time which is used by the DagSolid or G4TessellatedSolid class.
Figure 5 shows the total and class computation times as function of the number of facets for Geantino and ChargedGeantino in the sphere and small intestine models. For the cases of the G4TessellatedSolid class, we can hardly see the difference between the total computation time and the class computation time. The class computation time approaches the total computation time, meaning that, during simulation, the geometric calculation to determine the particle step length occupies most of the computation time. The computation time also rapidly increases with the number of facets which is mainly due to the increased usage of the ray-tracing functions in the G4TessellatedSolid class. Note that the computation time is almost proportional to the number of the facets.
Figure 5. Total and class computation times as a function of the number of polygon facets: (a) sphere with Geantino, (b) small intestine with Geantino, (c) sphere with ChargedGeantino, (b) small intestine with ChargedGeantino.
Download figure:
Standard image High-resolution imageFor the cases of the DagSolid class, when compared to the G4TessellatedSolid class, the computation time is significantly reduced for all of the cases for both the total and class computation times. Furthermore, in contrast to the G4TessellatedSolid class cases, the computation time does not increases rapidly as the number of facets increases, which is mainly due to the use of the ray-tracing acceleration algorithm functions in the DagSolid class. The increase of the class computation time approximately follow the theoretically expected log2N relationship, where N is the number of the facets around the particle. Accordingly, the difference of computation time increases with the number of the facets, and the maximum difference of computation time was 1562 and 680 times for 320 000 facets for Geantino and ChargedGeantino, respectively. This result shows that the use of the DagSolid class will effectively reduce computation time for complicated polygon-mesh models using many facets.
Tables 1 and 2 show the detailed computation time information of the ray-tracing functions for Geantino and ChargedGeantino, respectively. The result of the G4TessellatedSolid class shows that the Inside function takes the largest computation time in Geantino simulation, but the other functions take the more time in ChargedGeantino simulation. The DistanceToIn and DistanceToOut functions occupy the most computation time in ChargedGeantino simulation. This result can be explained by the fact that the ChargedGeantino, which moves helically in the magnetic field, creates more steps and thus more repeatedly calls the DistanceToIn, DistanceToOut functions, and SurfaceNormal functions.
Table 1. Computation time of ray-tracing functions in G4TessellatedSolid and DagSolid classes for Geantino.
| Inside/point_in_volume | DistanceToIn + DistanceToOut/closest_to_location + ray_fire | SurfaceNormal/get_angle | |||||
|---|---|---|---|---|---|---|---|
| Geometry | Number of facets | G4TessellatedSolid (s) | DagSolid (s) | G4TessellatedSolid (s) | DagSolid (s) | G4TessellatedSolid (s) | DagSolid (s) |
| Sphere | 5000 | 19.45 | 0.4 | 9.86 | 0.92 | 0.52 | 0.11 |
| 20 000 | 194.47 | 0.56 | 95.89 | 1.5 | 5.64 | 0.12 | |
| 80 000 | 1146.24 | 1.24 | 591.6 | 2.91 | 49.82 | 0.17 | |
| 320 000 | 3643.36 | 1.88 | 1901.32 | 5.3 | 158.38 | 0.21 | |
| Small intestine | 5000 | 33.15 | 0.31 | 25.42 | 1.14 | 2.55 | 0.26 |
| 20 000 | 177.24 | 0.39 | 152.84 | 2.05 | 16.4 | 0.57 | |
| 80 000 | 1017.25 | 0.48 | 1007.34 | 2.38 | 134.05 | 0.59 | |
| 320 000 | 3656.45 | 0.64 | 3499.94 | 3.62 | 466.01 | 0.62 | |
Table 2. Computation time of ray-tracing functions in G4TessellatedSolid and DagSolid classes for ChargedGeantino.
| Inside/point_in_volume | DistanceToIn + DistanceToOut/closest_to_location + ray_fire | SurfaceNormal/get_angle | |||||
|---|---|---|---|---|---|---|---|
| Geometry | Number of facets | G4TessellatedSolid (s) | DagSolid (s) | G4TessellatedSolid (s) | DagSolid (s) | G4TessellatedSolid (s) | DagSolid (s) |
| Sphere | 5000 | 82.42 | 0.24 | 662.33 | 9.14 | 33.60 | 2.07 |
| 20 000 | 348.55 | 0.63 | 3338.34 | 14.75 | 189.84 | 3.82 | |
| 80 000 | 1248.09 | 1.19 | 12 456.26 | 20.90 | 841.25 | 6.35 | |
| 320 000 | 3818.72 | 1.93 | 37 719.40 | 36.58 | 2871.52 | 26.79 | |
| Small intestine | 5000 | 113.37 | 0.27 | 1492.82 | 55.57 | 179.21 | 16.77 |
| 20 000 | 543.60 | 0.43 | 8332.11 | 86.77 | 1220.28 | 27.53 | |
| 80 000 | 1183.36 | 0.54 | 19 044.76 | 93.75 | 2977.50 | 34.80 | |
| 320 000 | 6745.34 | 0.69 | 108 023.60 | 134.45 | 18 309.61 | 65.19 | |
In general, the results show that all of the ray-tracing functions of the DagSolid class are much faster than the corresponding functions in the G4TessellatedSolid class for all of the cases considered in the present study. This result assures that the DagSolid class will be always faster than the G4TessellatedSolid class, regardless of the complexity of geometry, or the number of facets used in the model, and the kind of particles (charged or uncharged) simulated in Monte Carlo simulation.
3.2. Performance evaluation for real particles
The performance of the DagSolid class was evaluated for real particles (gamma, electron, neutron, and proton) again by comparing the computation time with that of the G4TessellatedSolid class. For this study, a 1 m diameter water sphere and a small intestine, which were composed of 80 000 facets, were irradiated by the 10 MeV pencil beam. To avoid the memory overflow of the profiler used for the computation time calculation, the number of primary particles emitted from the disk source was limited to 100. The physics library used for gamma and electron simulations was G4EmLivermorePhysics. For neutron and proton simulations, this study used the physics libraries of G4EmLivermorePhysics, HadronPhysicsQGSP_BIC_HP, G4EmExtraPhysics, G4Hadron-ElasticPhysics, G4QStoppingPhysics, and G4IonBinaryCascadePhysics. The default value of the range cut value was used for all particles.
Table 3 compares the computing times of the DagSolid and G4TessellatedSolid classes for real particles. The table also shows the ratio of the G4TessellatedSolid computation time and the DagSolid computation time. The DagSolid class is always much faster than the G4TessellatedSolid class. The difference of computation time was about 55–2395 times depending on transported particles for class computation time. For total computation time, the difference was less significant, but still 53–685 times. It is also shown that the improvement of the computation speed tends to be more effective for charged particles than for uncharged particles.
Table 3. Comparison of computation times of G4TessellatedSolid and DagSolid classes for real particle (gammas, electrons, neutrons, and protons).
| Class computation timea | Total computation time | ||||||
|---|---|---|---|---|---|---|---|
| Geometry | Particle | G4TessellatedSolid (s) | DagSolid (s) | Ratio | Using G4TessellatedSolid (s) | Using DagSolid (s) | Ratio |
| Sphere | Gamma | 12 138.56 (±1275.37) | 218.54 (±9.75) | 55.54 (±6.34) | 12 155.07 (±1275.07) | 229.34 (±9.88) | 53.00 (±6.01) |
| Electron | 13 650.85 (±736.94) | 12.04 (±0.30) | 1133.79 (±67.43) | 13 667.30 (±736.56) | 19.95 (±0.33) | 685.08 (±38.58) | |
| Neutron | 12 527.23 (±1293.54) | 51.67 (±3.30) | 242.43 (±29.44) | 12 564.68 (±1293.73) | 79.48 (±3.34) | 158.08 (±17.58) | |
| Proton | 18 898.92 (±1829.96) | 7.89 (±0.06) | 2395.30 (±232.65) | 18 943.00 (±1831.05) | 35.87 (±0.04) | 528.10 (±51.05) | |
| Gamma | 866.30 (±46.79) | 3.16 (±0.11) | 274.05 (±17.79) | 873.19 (±46.85) | 9.70 (±0.13) | 90.20 (±4.98) | |
| Small intestine | Electron | 7433.65 (±378.09) | 15.86 (±0.28) | 468.70 (±25.22) | 7444.74 (±378.22) | 23.86 (±0.13) | 312.02 (±16.49) |
| Neutron | 2394.46 (±124.83) | 5.13 (±0.29) | 466.75 (±35.96) | 2422.47 (±124.83) | 32.49 (±0.35) | 74.56 (±3.99) | |
| Proton | 11 888.01 (±1318.30) | 16.73 (±0.31) | 710.58 (±79.91) | 11 930.21 (±1318.32) | 46.09 (±0.44) | 258.84 (±28.71) | |
aClass computation time includes polygon mesh model's data loading time.
Firstly, in the case of the sphere, it was found that the improvement of class computation speeds for charged particles were much better than those of uncharged particles due to the 'point_in_volume' function of DagSolid class. In the present study, it was observed that the computation speed of point_in_volume function can be faster when the interaction occurs near facets of geometry. In principle, charged particles have relatively a lot of interactions near the facets after entering the geometry as compared with uncharged particles, resulting that the performance of the point_in_volume function is more effective for charged particles. On the other hand, in the case of the small intestine, it was found that it is no longer necessary to find the trend shown in the case of the sphere, because in the complex geometry, complexity of geometry has a significant effect on possibility of interactions near the facets regardless of particle types.
Secondly, It is difficult to find a direct correlation in speed improvement between the class and total computation times. In other words, in case of sphere geometry, while the improvement of class computation time was most significant for proton (2395 times), the improvement of total computation time for proton (528 times) was not even close to those of class computation time. Otherwise, the improvement of class computation time and total computation time for gamma (55 and 53 times) were almost same. This is due to the fact that the fraction of the computation time for the geometric calculation is different for different particles. The majority of computation time is generally used to calculate 'geometric-step-length' and 'physical-step-length,' which are calculated by the Solid class and GetPhysicalInteraction-Length (GPIL) function, respectively. The computation time of the GPIL function is different for different particles. In general, the GPIL function uses more computation time for hadrons (neutrons and protons) than gammas or electrons mainly due to the fact that hadrons have more types of interaction, for each of which the physical interaction length must be sampled in Geant4, than gammas and electrons. In addition, the data structure and physics models tend to be more complicated for hadrons. Figure 6 shows the fraction of the computation times which is one of the simulation results for the DagSolid class and the physics models called by the GPIL function. As a result, for gamma, the fraction of the DagSolid computation time is ∼95% of the total computation time, which explains the reason why the improvement of the computation speed for total computation time is ∼95%. Otherwise, for proton, even though the fraction of the DagSolid class computation time is ∼77%, the improvement of the computation speed for total computation time is ∼22%.
Figure 6. Fraction of computation times for DagSolid class and physics-related classes in total computation time for gammas, electrons, neutrons, and protons.
Download figure:
Standard image High-resolution image3.3. Accuracy evaluation
The accuracy of the DagSolid class was also evaluated by comparing the calculation results of the DagSolid class with those of the G4TessellatedSolid class. For this study, the polygon-mesh computational human phantom named PSRK-Man, developed by Kim et al (2011), was irradiated by a broad parallel photon beam at the anterior-posterior direction. Then, organ dose conversion coefficients (DCCs), defined as organ absorbed dose per air kerma, were calculated and compared. In the calculation, a 90 cm radius source disk was defined to generate the broad parallel photon beam. The simulated photon energies were from 0.015 to 10 MeV. The number of primary photons emitted from the source disk was from 107 to 108 considering the statistical uncertainties of the calculated values. For precise dose calculations, the Livermore electromagnetic physics library was used and the secondary production cut values for photons and electrons were set to 1 µm.
Figure 7 compares the calculated DCC values for six selected organs: bone, brain, large intestine, liver, lungs, and skin. It can be shown that the DagSolid and G4TessellatedSolid classes produce the same results within statistical uncertainties, validating the accuracy of the DagSolid class.
Figure 7. Comparison of organ DCCs of the PSRK-Man when using G4TessellatedSolid and DagSolid classes. Error bars are included on each data point.
Download figure:
Standard image High-resolution image
4. Conclusions
In this work, a new Geant4 solid class named DagSolid was developed based on the ray-tracing acceleration algorithms of the DAGMC library. The performance of the developed class was then evaluated for both virtual and real particles. The use of the developed DagSolid class significantly improved the computation speed compared to using the existing G4TessellatedSolid class for all of the cases considered in the present study. For real particles (gammas, electrons, neutrons, and protons), the improvement of computation speed was within the range of 53–685 times depending on simulated particles. The improvement of the computation speed was more significant when there were more facets, meaning that the DagSolid class can be used more effectively for complicated geometries with many facets than for simple geometries. The improvement was also more significant for charged particles than for neutral particles. The detailed analysis of the ray-tracing functions assures that the DagSolid class will be always faster than the G4TessellatedSolid class, regardless of the complexity or the number of facets used in the geometry and the kind of particles (charged or uncharged) transported in the simulation. In the present study, the developed DagSolid class was used with the G4NormalNavigation class which is not optimized for the DagSolid class. We believe that the computation speed of the DagSolid class can be further improved by future development of a dedicated navigation system with optimization.
Acknowledgments
This research was supported by the National Nuclear R and D Program, Basic Science Research Program, and Global PhD Fellowship Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (no. 2012-K001146, 2011-0025496, 2011-0030970, 2011-0007318). The authors thank for the valuable comments from Dr Makoto Asai (SLAC) on this research.