A model for the interaction of high-energy particles in straight and bent crystals implemented in Geant4

A model for the simulation of orientational effects in straight and bent periodic atomic structures is presented. The continuum potential approximation has been adopted.The model allows the manipulation of particle trajectories by means of straight and bent crystals and the scaling of the cross sections of hadronic and electromagnetic processes for channeled particles. Based on such a model, an extension of the Geant4 toolkit has been developed. The code has been validated against data from channeling experiments carried out at CERN.


Introduction
The interaction of either charged or neutral particles with crystals is an area of science under development. Coherent effects of ultra-relativistic particles in crystals allow the manipulation of particle trajectories thanks to the strong electric field generated between crystal planes and axes [1,2,3]. Important examples of the interaction of neutral particles in crystals include production of electron-positron pairs and birefringence of high energy gamma quanta [4,5]. Radiation emission due to curved trajectories in bent crystals has been seen to enhance photon production through bremsstrahlung, channeling radiation, PXR, undulators [6,7,8,9,10] and recently through volume reflection and multiple volume reflection [11,12]. The inelastic nuclear interaction rate is known to be modified by channeling and volume reflection [13].
Various applications of orientational phenomena with crystals have been proposed and investigated such as • beam steering [14] • extraction and collimation in circular accelerators [15,16,17,18] • splitting and focusing of external beams [19].
On the strength of recent optimizations in manufacturing techniques [20,21] and crystal holders [22], bent crystals have also been proposed as beam collimators [23] and extractors [24,25] for the LHC. Indeed, the reduction of the nuclear interaction rate for channeled positive particles in bent crystals [13] has been proven to lower beam losses throughout the SPS synchrotron [26]. In addition, crystals have been proven to work in the same way for Pb ions [27]. Beam extraction with bent crystals obviates the need for electric current (and its associated maintenance costs) to deflect a particle beam, especially at high energy.
Coherent effects for the interaction of particles with aligned structures have always exploited opportunities furnished by the most advanced computers and computational methods of the current period. In 1963 Oen and Robinson [28] investigated particle behavior in crystal lattices using an IBM 7090 and noticed that particles moving near planes and axes had anomalously long ranges. Their Monte Carlo code was based on a binary collision model and was capable of predicting the experimental results observed in 1963 [29]. The binary collision model allows the determination of the trajectory of a low-energy particle in a crystal with high precision, but it is computationally expensive due to the need to solve the equation of motion of a particle with an integration step smaller than the cell distance between two neighboring atoms, which is typically less than 1Å.
In 1965 Lindhard showed [30] that the motion of relativistic charged particles under channeling is well approximated with classical physics equations and, as a consequence, also the potential and its related quantities. In particular, Lindhard proposed the idea of approximating the interaction of swift particles with aligned atoms by the interaction with a single string of atoms, i.e. the so-called continuous potential approximation. The same idea holds for interactions with a crystalline plane. By adopting the continuous approximation, the equation of motion can be solved in one dimension for planar channeling with an integration step of up to 1 µm for GeV particles [31,32,33,34,35,36], with a high computational cost for each particle due to the necessity of integrating over the full particle trajectory. As an example of the capability of such a method, in 1987 Vorobiev and Taratin predicted the volume reflection phenomenon in bent crystals [3] which was first observed in 2006 by the H8RD22 collaboration [37] .
In the last few years experiments at CERN external beam lines [38,39,40,41,42,43,44,45,46,47,48,13,49,50,51,52,53] have made available a large amount of high-resolution data thanks to the ability to track trajectories of ultra-relativistic particles with µrad resolution [54]. Such improvements in the experimental knowledge of orientational effects has led to the development of Monte Carlo simulations based on the experimental cross sections of orientational phenomena [55,51]. With this model, very high throughput is achieved but the scaling of dechanneling models is inaccurate due to the lack of a dedicated campaign of measurement. Moreover such an approach is not suitable to describe the cross section variation of physical phenomena for channeled particles.
In 2013 a model for the simulation of planar channeling of positive particles in bent crystals was proposed [56] for Fluka [57]. The model relies on the continuous potential approximation. Channeling particles follow the bend of the crystal until they either exit the crystal or are dechanneled. Single scattering is treated via a microscopic approach and its effects on transverse energy are cumulative. The number of events caused by nuclear interactions is modified for channeled particles due to the variation of the relevant cross sections.
Nowadays Monte Carlo simulations of the interaction of particles with matter are usually done with download-able toolkits such as Geant4 [58] and Fluka [57]. Such Monte Carlo codes are continuously expanded and improved thanks to the collaborative effort of scientists from around the world. Geant4, an objectoriented toolkit, has seen a large expansion of its user community in recent years. As an example, applications simulated by Geant4 range from particle transportation in the ATLAS detector [59] to calculations of dose distribution curves for a typical proton therapy beam line [60], and from radiation analysis for space instruments [61] to early biological damage induced by ionizing radiation at the DNA scale [62].
In November 2012 a version of Geant4 with the first implementation of a physical process in a crystal was released. In fact, an extension was developed to manage properties of periodic structures in Geant4. The process of phonon propagation within Ge crystals was added to the toolkit [63,64], but no orientational effects for charged particles were developed at that time. All physical processes implemented in Geant4 require the interaction length of the process to be calculated. Thus, the cross section of the process must be known or computed in order for the process to be added to the toolkit. In addition, the concurrent presence of many physical processes forces the use of an integration step greater than a µm to limit the computational time. As a result, the full solution of the equation of motion is not suitable. An alternative approach would be to simulate orientational effects using experimental data, but such data (channeling of negative particles in bent crystals, for example) do not currently exist.
In this paper we present a general model for the simulation of orientational effects in straight and bent crystals for high energy charged particles. The model is based on the continuous potential approximation but does not rely on the full integration of particle motion. The model has been implemented in Geant4, and validated against experimental data.

Model
In this section the models for channeling and volume reflection are presented. Since they are based on the continuous potential approximation, a resume of the Lindhard work and its range of applicability is presented here.

Continuous approximation
The continuous approximation was developed by Lindhard to describe channeling and its related phenomena, but can be extended to all orientational phenomena because the same approximations hold. Coherent effects are primary phenomena; they govern the paths of primary particles, and not secondary ones, which are determined by the path. Thus, four basic assumptions can be introduced for particles under orientational conditions: • scattering angles may be assumed to be small. Indeed, scattering at large angles implies complete loss of the original direction.
• Because the particle moves at small angles with respect to an aligned pattern of atoms and because collisions with atoms in a crystal demand proximity, correlations between collisions occur.
• Since the coherence length l of a scattering process (l = 2E/q 2 , where E is the particle energy and q the transferred momentum) is larger than the lattice constant, a classical picture can be adopted.
• The idealized case of a perfect lattice may be used as a first approximation.
Under these assumptions, the continuous approximation can be inferred, and the potential of a plane of atoms U (x) can be computed by taking the average of the detailed potential along the direction of motion of the particle.
where d p is the interplanar distance, N is the atomic density and V (r) is the potential of a particle-atom interaction. By using the screened Coulomb potential approximation for V (r), the interplanar potential becomes where Z 1 and Z 2 are the atomic numbers of particle and ion, respectively, e is the elementary charge, a T F is the Thomas-Fermi radius and U max is the maximum of the potential.
The basis of the continuous approximation relies on the qualitative assumption that many consecutive atoms contribute to the deflection of a particle trajectory. Thus, for relativistic particles, the time of collision ∆t ≈ ∆z/c multiplied by the momentum component parallel to the plane of atoms, p z ∼pcosθ, must be large compared to the distance d z between atoms along the particle direction, where p is particle momentum, p z is the momentum component along the particle direction and θ is the angle between the particle direction and the crystal plane orientation. Since the collision time is approximately ∼r min /(vsinθ), where r min is the minimal distance of approach, the condition for the continuous approximation to hold is In the most restrictive form r min is determined by the condition that the transverse kinetic energy cannot exceed the transverse potential energy at r min .
Therefore, from previous equations, a condition can be derived for which the continuous approximation is still valid: Two terms appear in this condition. One refers to the Lindhard angle of channeling θ L = 2U max /(pβ), which determines the maximum angle for channeling. The other is more interesting because it implies that θ < a T F /d z ∼ 0.5Å/1Å ∼ 0.5 is very large compared to θ L at high energy. Thus, the continuous approximation is still valid for angles greater than θ L as long as that particle does not approach closer than r min to a nucleus.
The continuous potential approximation can be extended to regions closer than r min to the atomic position by treating in more detail atomic displacement in the structure. In fact, since the crystal temperature is usually higher than 0 K degree, atoms vibrate around their center of mass. By averaging the thermal vibration amplitude over space and time, the probability density function for the position of atoms can be derived. Thus, the continuous approximation can be extended to regions closer to the center of vibration of atoms. Because the averaging is due to thermal fluctuations, such an approximation is not valid at very low temperatures and the limits of the continuous approximation must be kept in mind.

Channeling
When a charged particle hits a crystal aligned with an atomic plane it can be trapped by the strong electromagnetic field between two planes, thus undergoing planar channeling. Channeled particles follow the direction of the crystal plane, oscillating between or across planes if the particle charge is positive or negative, as shown in Fig. 1.a. Under channeling conditions positive particles penetrate deeper into the crystal relative to the un-aligned orientation because the trajectory is repelled from the nuclei. On the other hand, negative particles interact more frequently because of their attraction to zones with high densities of nuclei.
The continuous interplanar potential for main planes in crystals [65] can be approximately described by a harmonic potential well for positive particles, as shown in Fig. 1.c. However for negative particles, being attracted by nuclei, the interplanar potential must be reversed and becomes non-harmonic with a minimum in the middle of the potential well, as shown in Fig. 1 Because the trajectory is strongly affected by such a potential, positive and negative particles under channeling trace different shapes in phase space (see Fig. 1 Channeling holds for particles with transverse energy E x,θ lower than the maximum of the potential well depth U 0 , i.e., E x,θ < U 0 . Such particles follow the channeling plane or axes until they exit the crystal or are dechanneled. The dechanneling mechanism behaves the same for both straight and bent crystals. If all the processes which lead to dechanneling are disabled, a particle remains under channeling for the entire crystal length as long as E xin,θin < U 0 , where x in and θ in are the impact position and incoming angle with respect to channeling plane. Thus, conservation of transverse energy allows the treatment of channeling through a knowledge of the initial impact position on a crystal channel x in and the angle with respect to the crystal plane θ in : Solving the equation of motion requires point-by-point knowledge of the transverse position and transverse momentum of an oscillating particle. However, by choosing a crystal which extends along the beam for more than one oscillation period of a channeled particle, the energy level occupied by a particle in the electrostatic potential well generated between atomic planes or axes is the only physical quantity to link initial to final parameters in a real-case study. By imposing a continuous and uniform distribution in position x out for a channeled particle of energy E t , an outgoing angle θ out is generated by evaluating Therefore, information regarding x in and θ in can be condensed into a single variable E t , which determines the occupied energy level of a channeled particle and allows the computation of the outgoing distribution of channeled particles by means of E t and the continuous potential.
For bent crystals the model is still valid. The sole difference relies on the modified potential in the non-inertial reference frame orthogonal to the crystal plane or axis. In fact, centrifugal force acting on the particle in this frame pulls down the potential barrier resulting in a shallower potential well. Thus, the condition for channeling holds with a modified maximum potential and transverse energy related to the non-inertial reference system.
The presence of torsion in a crystal spoils channeling efficiency in bent crystals [22]. Indeed, the orientation of the channeling angle with respect to the beam direction changes with the impact position on the crystal surface. Since a beam has a finite size, two particles with the same direction and the same impact position on the potential well but different impact positions on the crystal surface have different transverse energies. This effect is introduced in the simulation by changing the plane direction with respect to the impact position on the crystal surface.
Another important parameter for channeling in bent crystals is the miscut [66], which is the angle between the lateral surface of a crystal and the atomic planes. Only the trajectories of particles channeled near a crystal edge are affected by the presence of the miscut, because it modifies the total length of the bent plane of channeling. This effect is introduced by defining the plane orientation independently of the crystal volume.

Dechanneling and volume capture
Particles which no longer satisfy the channeling condition have suffered dechanneling. Un-channeled particles which enter the channeling state undergo volume capture. Dechanneling and volume capture take place when particles interact incoherently with nuclei or electrons. Indeed, a channeled particle can acquire enough transverse energy to leave the channeling state by exceeding the maximum of the potential well, or an un-channeled particle can lose energy and decrease its transverse energy by passing under the maximum of the potential well.
All physical phenomena occurring for a channeled particle are strongly affected by the occupied energy levels. As shown in Fig. 3, the average density of material seen by a particle traversing a crystal aligned with its planes is strongly affected by the transverse energy of the particle. Thus, the probability of interaction with nuclei and electrons has to be weighted as a function of the transverse energy. Kitagawa and Ohtsuki [67] demonstrated that there is a linear dependence between the incoherent interaction rate and the material density. Therefore, the modified cross section σ(E t ) of each phenomenon is where σ am is the cross section in amorphous material, n am is the density of the amorphous material and n(E t ) is the modified density. Since each incoherent phenomenon produces a variation of transverse energy ∆E t , the maximum distance traveled by a channeled particle is limited by the probability of overcoming or under-passing the potential barrier during each step. As an example, the modified rms σ is (E t ) for incoherent scattering on nuclei depends approximately on the square root of the traversed material density [68]. Thus, the step can be limited by the condition where σ is,am ∼ Es pβ z X0 is the rms of the incoherent scattering in the amorphous material, E s = 13.6 MeV, β is the particle velocity in units of the speed of light and X 0 is the radiation length of the material. Thus, the step ∆z is Such an approach can be applied to all the concurrent incoherent processes to determine the maximum step size by comparing the contributions. The model is still valid for positive and negative particles since no restriction has been applied.
The same response for the interaction probability is not obtained by averaging the density over an oscillation period. Therefore, this model must be adapted for crystal with lengths along the beam greater than one oscillation period. On the contrary, by integrating the particle trajectory it is possible to determine the interaction probability for each step depending on the position in the channel. Thus, the peculiar characteristic of channeling in the first layers of a crystal is lost by the averaging.

Volume reflection
When charged particles cross a bent crystal tangent to its planes they are "reflected" in the direction opposite to the bending curvature. This is called volume reflection. In fact, the particle is deflected by the continuous potential barrier of one plane, but immediately leaves the channel because the barrier of the opposite plane is lowered due to bending, and thus the particle cannot be trapped under channeling. Therefore, the condition for volume reflection holds when the projection of the particle momentum on the direction of a plane changes sign. Volume reflection and related phenomena limit the maximum allowed step length. Indeed, particles can be captured into a channeling state if they lose enough transverse energy to fulfill the channeling condition E x < U 0 . Thus, the step length must be comparable to the oscillation period near the turning point. The distance of a particle to the tangency point in a bent crystal must be evaluated at each step to set the step size at the proximity of the interesting region.
For a slightly bent crystal, the mean deflection angle for volume reflection is approximately ∼ 1.4θ L for positive particles [39] and ∼ 0.8θ L for negative particles [69], where θ L = (2U 0 pβ) 1/2 is the Lindhard angle [30]. In fact, a positive particle spends more time in a zone within which the particle has higher transverse velocity, while the opposite is true for a negative particle, as shown in Fig 2.a. By decreasing the radius of curvature, the mean deflection angle decreases. Indeed, the deflection angle of volume reflection depends on the transverse energy at the turning point. The more the crystal is bent, the larger is the angular spread [39]. In Fig. 2.b the potential in the non-inertial reference frame orthogonal to the crystal plane is shown. The maximum energy difference delimits the reflection region. The potential shapes the trajectories of particles with transverse energies a bit above the maximum of the potential. These are the so-called over-barrier particles. By approximating the potential in the reflection region with a linear function, the volume reflection angle becomes proportional to the position of the turning point. Thus, the deflection angle can be generated by adopting a continuous and uniform distribution proportional to the reflection region.

Average density
The average density seen by a particle undergoing orientiational effects is a very important parameter for the model proposed in this paper. As shown in Fig. 3 the average density is strongly affected by the transverse energy for channeled and over-barrier particles. The computation is made with the DYNECHARM++ code [36] in which all incoherent processes are disabled. The DYNECHARM++ code is based on the full solution of the equation of motion in the continuous potential and allows the computation of electric characteristics of the crystal through the ECHARM calculation method [70,71]. Therefore, the density as a function of transverse energy for complex atomic structures and for many planes and axes can be computed.
The calculation of average density by DYNECHARM++ is very accurate, but, at the same time, can be very slow. Thus, a fast version has been developed to compute the average density. By integrating the density of all the possible states in which a particle with transverse energy E t can exist, the approximate average density ρ(E t ) can be computed.
Since no approximation was imposed on U (x) and ρ(x), this approximation is still valid for any potential and average density function. Experiments with orientational effects rely mostly on the use of crystal planes with low Miller indexes. For these the averaged nuclei density ρ(x) is analytically derived starting with the density of nuclei function averaged over thermal and space fluctuations where u t is the thermal vibration amplitude. Thus, the average density is where E t ≤ U 0 and U (x) = U 0 at x = d p /2. Fast calculation models do not take into consideration the time spent in a particular region by a particle. The model can be applied to compute density for both positive and negative particles. In Fig. 3 the average density ratio of electrons and nuclei is shown for both the fast model and DYNECHARM++.

Geant4 implementation
The Geant4 toolkit allows new physical processes to be added to the standard ones it already provides. The new process must provide its mean interaction length and how particle properties are affected by the interaction. The model proposed in this paper has been implemented by a class describing the orientational process, and two wrapper classes that modify the material density in existing processes. In addition, classes for the fast calculation of crystal characteristics have been developed to allow the simulation of orientational processes with no need for external software.

ProcessChanneling
The ProcessChanneling class used for the implementation of orientational processes inherits from the G4VDiscreteProcess virtual class. The channeling process is valid only in a volume with a crystal lattice. A crystal lattice is present if a XVPhysicalLattice class is attached to the physical volume. The XVPhysicalLattice class collects all the crystal data, such as unit cell and lattice type.
When a particle crosses the boundary between two geometrical volumes, one with and one without a crystal lattice, the channeling process limits the step of the particle and checks if the particle is subject to orientational effects. A uniformly distributed random number is generated to determine the impact position of the particle on the crystal channel x in and, consequently, to compute the initial potential energy U 0 . The particle momentum is projected on the channeling plane to evaluate the transverse momentum. The initial transverse energy E xin,θin is computed through Eq. 6. Thus, E xin,θin is adapted to find the modified density ρ(E t ). If the particle satisfies the channeling condition, E xin,θin < U 0 , the channeling process proposes to the Geant4 core an alignment of the particle momentum with the direction of the channeling plane. The condition for channeling is recomputed until the particle exits the volume with the crystal lattice.
Volume reflection occurs only for bent crystals under the condition defined in section 1.4. Under volume reflection, the particle momentum vector is rotated by the volume reflection angle around the axis orthogonal to the channeling plane.

Crystal
The XVPhysicalLattice class was introduced into Geant4 to attach a crystal lattice to a physical volume. This class was extended to collect all the crystal data, such as unit cell and lattice type. Specifically, two classes, XLogicalBase and XUnitCell, were added. The first contains the kind and disposition of the atoms.
The second groups the unit cell information, i.e., the sizes and the angles of the cell, and holds a vector of pointers to as many XLogicalBase objects as needed. In the XVPhysicalLattice class a pointer to an XUnitCell object is stored. The information stored in a XUnitCell object may be used to compute electrical characteristics under the continuous approximation of the channeling processes. Thus, the XVCrystalCharacteristic class was developed.

Wrappers
At each step in a crystal, the particle momentum can be modified by any of the Geant4 processes. Such modifications vary the transverse energy of a particle and may cause dechanneling, that is, the overcoming of the potential well maximum. As stated in section 1.3, the average densities of nuclei and electrons change as a function of the transverse energy of the particle. Thus, these densities should be recomputed at each step and used to modify the cross section of the physics processes which depend on the traversed quantity of matter (see section 1.5 ).
In order to modify the cross-section of existing processes and to preserve code reusability for future releases of Geant4, wrapper classes for the G4VDiscreteProcess and G4VEnergyLoss classes were developed. For both these classes, the interaction length of discrete processes is resized proportionally to the modified material density. For the continuous energy loss of the G4VEnergyLoss processes, the traversed length is resized proportionally to the modified average density. For each wrapped process a wrapper object must be instantiated. The wrappers need only the average density to recompute the process cross section. Thus, in principle, it may work independently of the ProcessChanneling class.

Examples of calculation
Model validation has been completed by comparison with published experimental data. Experiments studying the efficiency of channeling vs. incoming angle [44], the rate of inelastic nuclear interaction under channeling [13], and the channeling efficiency dependence on radius of curvature for bent crystals [72], were simulated for positive particles. For negative particles, simulations of the dechanneling length for high-energy pions [52] was performed. Comparing simulations to experiments allowed both the precision of the model and the quality of the Geant4 implementation to be checked.
A bent crystal was modeled as a small fraction of a toroid with a bending radius on the order of a meter and a length on the order of a mm along the beam direction, matching the dimensions used in the experiment. Though torsion can be simulated, it was set to zero for all the current simulations. This has no effect on the agreement of simulation with data, even though the experimental data have been corrected for torsion. In addition, the miscut value has no influence because only particles impinging far from crystal edges have been used in the analyses.
As in the experimental setups, three silicon detectors were inserted into the simulation along the beam direction to track the particle. For measurement of the rate of inelastic nuclear interaction, two scintillators were added to reproduce the experimental setup of Ref. [13]. To speed up simulation, volumes other than crystal and detectors have been filled with galactic vacuum (G4 Galactic material).

Positive particles
In Fig. 4, the channeling efficiency as a function of incoming angle is superimposed on experimental results (Fig. 3 of Ref. [44]) and a Monte Carlo simulation with complete integration of the trajectories. The maximum efficiency for channeling in Geant4 is in good agreement with experimental data as well as efficiency in the tails. However, fair agreement is obtained in the region between maximum efficiency and tail, with ∼ 5% deviation in efficiency. In this region the model lacks accuracy because the trajectories were not completely integrated. Thus, such behavior is to be ascribed to the shape of the interplanar potential used in simulation for both the models. Fig. 5 should be compared with Fig. 5 of Ref. [13]. The rate of secondary particles as a function of the beam angular spread is shown normalized to an amorphous condition. The standard Geant4 release without the channeling extension has been used for simulations with amorphous Si and with no crystal. Simulations are in agreement with experimental data. The channeling extension allows the correct modification of the cross sections of incoherent phenomena, reducing the rate with respect to amorphous materials. Discrepancies are observed for small angles and the slope of the two curves are different. Table 1 presents the deflection efficiency for channeling vs. radius of curvature. Experimental data and DYNECHARM++ simulations are taken from Ref. [72]. As the critical radius R = R c is approached, the discrepancy between experimental data and simulation increases. Such behavior is also recorded for DYNECHARM++ simulations. As stated in Ref. [72], the discrepancy must be ascribed to the lack of knowledge of the exact density distribution between atomic planes. An important feature of Geant4 is its capability to evaluate the number of particles which suffer nuclear interaction or are scattered at large angles. In Table 1 the fraction of "lost" particles , i.e., which do not hit the last detector, is reported.

Negative particles
In Ref. [52] the interaction of 150 GeV/c negative pions with a bent Si crystal has been studied in order to measure the dechanneling length for negative particles. A dechanneling length of 1.54 ± 0.05 mm was obtained by Geant4 simulation with the density computed by DYNECHARM++ code, compared to 0.71 ± 0.05 mm with the density computed by the new Geant4 model (see Fig. 6.a). The dechanneling rate is increased due to the stronger incoherent scattering with nuclei and electrons. Thus, the model for negative particles is very sensitive to the interaction rate in one oscillation period, since a big discrepancy between the two simulations exists. Indeed, the discrepancy of the dechanneling lengths becomes large for the channeling efficiency, which goes from 26.8 ± 0.5% to 6.2 ± 0.5%. As a consequence, by computing accurately the average density experienced by a particle the model is able to output the measured dechanneling length. The same configuration was used to simulate channeling of 150 GeV/c π + . The comparison between positive and negative pions is shown in Fig. 6.b. The deflection efficiency for π + is ∼ 70%, which is greater than than for π − . This result demonstrates that the channeling model developed for Geant4 allows positive and negative particles to be managed differently thanks to the wrapper classes.

Computation time
The Geant4 code has been compared to the DYNECHARM++ code in order to evaluate advantages of the approach proposed in this paper in terms of computation. The same initial conditions have been used as in Ref. [36] : a 400 GeV/c proton beam interacting with a 1.94 mm thick (1 1 0) Si bent crystal with a 38 m radius of curvature. The Geant4 single-threaded version 10.00b has been adopted and only a discrete single scattering model [73] has been added to its list of physics processes. The computer was the same as that used for DYNECHARM++ test, i.e. a personal computer with 8 GB of RAM and an Intel(R) Core(TM) i7-2600K CPU running at 3.40GHz. Computation time was approximately 14 ms per particle in Geant4 vs. 38 ms per particle in DYNECHARM++, in spite of the greater complexity of the Geant4 code. This result is explained by considering the number of steps required by the two models adopted for the simulation. Full integration of trajectories requires step sizes much smaller than the oscillation period in the potential well. On the contrary, the Geant4-based model allows the use of a step size comparable to the oscillation period.

Conclusions
The exploitation of orientational processes in crystals to manipulate particle trajectories is currently a topic of intense interest in physical research, with possible applications for the LHC for beam collimation [23] and extraction [24,25]. A physical model suitable for the Monte Carlo simulation of such processes has been developed. This model relies on the continuous potential approximation. The model makes use of the transverse energy in the non-inertial reference frame orthogonal to the channeling plane in order to discriminate between channeled and un-channeled particles. The average density experienced by a channeled particle is evaluated in order to compute the modification of the cross section for hadronic and electromagnetic processes. The model represents an extension of the Geant4 toolkit. The code has been validated against data collected by experiments at CERN. It demonstrates that Geant4 is able to compute the deflection efficiency for channeling and the variation of the rate of inelastic interactions under channeling.     GeV/c π − passed through a 1.91 mm long silicon crystal bent at R = 19.2 m along (110) planes. Only particles hitting the crystal within an angle of 5 µrad are selected. The average density experienced by a channeled particle has been computed by DYNECHARM++ (Geant4 a.) and an algorithm implemented in Geant4 (Geant4 b.). Experimental data was published in Ref. [52]. The dechanneling length was evaluated with the method proposed in the same reference. b) Geant4 simulation for the same crystal with 150 GeV/c π − (a. and b.) and 150 GeV/c π + .