Mark ; Thomas
email: thomam@gosh.nhs.uk

Great

M L Thomas

D Roebuck

R F Howard

The effect of volume of local anaesthetic on the anatomical spread of caudal block in children age 1-7 years

Keywords: child < Age, infant < Age, regional Anaesthetic techniques, regional, caudal, paediatric, anatomy

Objectives

To examine the anatomical spread of caudal local anaesthetic solution in children age 1-7 years.

Aim

To determine whether incremental increases in the volume of caudal injections of 0.5, 0.75 and 1.0 ml/kg result in reliable (>90%) and potentially clinically significant increases in the number of vertebral segments reached.

3

Results

There were 15 children in each group, they were similar in terms of age, height and weight. Spread was observed between the 5 th Lumbar (L5) and 12 th Thoracic (T12) vertebral levels. 1 ml/Kg results in a small but significantly greater spread of solution than 0.5 ml/Kg (p<0.05), but there was no difference between 0.5 and 0.75ml or between 0.75 and 1.0ml. No volume reliably reached a level higher than the second lumbar vertebra (L2).

Conclusions

Incrementally increasing the volume of injectate between 0.5 and 1.0 results in a modest increase in spread of the caudal solution. It is unlikely that volumes of less than 1ml will reliably reach a vertebral level that is higher than L2.

Background

Caudal block is one of the most frequently performed paediatric regional analgesic techniques. Traditional formulae suggest that changes in the volume of caudal injectate in the range 0.5-1.0 ml/kg would have clinically useful effects.

Methods

In a single blind design, 45 children age 1-7 years undergoing caudal block received one of three pre-determined volumes (0.5, 0.75 and 1 ml/Kg) of local anaesthetic solution containing radio-opaque contrast under controlled conditions. Following X-ray examination, the anatomical spread of the block was reported by a radiologist blinded to the volume of solution received.

Introduction

Caudal epidural block is one of the most popular regional analgesic techniques performed in children [START_REF] Cook | A national census of central neuraxial block in the UK: results of the snapshot phase of the Third National Audit Project of the Royal College of Anaesthetists[END_REF][START_REF] Armitage | Local anaesthetic techniques for prevention of postoperative pain[END_REF]. The extent of analgesia following epidural blockade with local anaesthesia is dependant on the anatomical spread of solution within the epidural space, this in turn depends on the volume injected [START_REF] Hong | A comparison of high volume/low concentration and low volume/high concentration ropivacaine in caudal analgesia for pediatric orchiopexy[END_REF][START_REF] Mcgown | Caudal analgesia in children. Five hundred cases for procedures below the diaphragm[END_REF]. Physical and developmental factors also influence the spread of a given volume of solution and because of this a number of predictive formulae of varying complexity have been devised in order to predict the required volume for different circumstances; it has been suggested that analgesia reaching to T10 or above can be obtained with only 0.7-1.0 ml/Kg by some authors, but much higher volumes have been recommended by others [START_REF] Hong | A comparison of high volume/low concentration and low volume/high concentration ropivacaine in caudal analgesia for pediatric orchiopexy[END_REF][START_REF] Schulte-Steinberg | Spread of extradural analgesia following caudal injection in children. A statistical study[END_REF](6)(7)(8). However, we are not aware of any studies directly examining the anatomical spread of volumes less than 1ml/kg. surgery with a caudal block as part of their normal anaesthetic plan were recruited. Sample size calculation was based on detecting a 2 vertebral level difference between incremental volumes with an 80% power and a significance level of P<0.05. The study was approved by the local ethics committee and the study was performed according to the declaration of Helsinki. Parental consent was obtained and child assent sought and obtained when appropriate. 15 children were allocated to each of the three caudal volume groups. In general, for perineal/lower limb surgery 0.5 ml/kg was used, for groin surgery 0.75 ml/Kg, and those children scheduled for higher abdominal surgery received 1 ml/Kg: children were allocated to the appropriate group based on the operation for which they were scheduled. The slight variability in the allocation of surgery type to each group reflects difference of opinion that exists between the relationship of the height of the block to the volume of caudal used.

The local anaesthetic solution was made up just prior to administration by adding 1 ml Omnipaque  300 mg I/ml X-ray contrast medium to every 9 ml of 0.25% bupivacaine. This solution was well mixed prior to administration. This resulted in a 0.225% solution of bupivacaine containing 30 mg I/ml of Omnipaque . Omnipaque (300 mg I/ml) has been shown to be chemically compatible with bupivacaine 0.25% and 0.125% at room temperature [START_REF] Van Asten | Compatibility of bupivacaine and iohexol in two mixtures for paediatric regional anaesthesia[END_REF]. The pH of the 0.25% solution was 7.1 and the osmolality 370 mOsm/kg. Thus the A single lateral X-ray was taken upon positioning the patient on the operating table. The gonadal region of the patient was shielded so as to minimise the Xray exposure (the dose of radiation required represents less than half the expected annual exposure to a member of the public, parents were explicitly made aware of this risk during the consent process).

Results

Demographic variables between children in the three groups were similar in terms of age, weight, height and the time from caudal injection to X-ray, as shown in Table 1. A typical X-ray obtained following caudal injection containing radio-opaque dye with the dark column of dye clearly visible in the epidural space is shown in Figure 1.

The number of blocks reaching a given vertebral height for each volume group is shown in Table 2, the proportion for each group is given in parenthesis. The highest vertebral level achieved in at least 90% of subjects in each of the 3 volume groups is shown in Figure 2; this was L2 for 1ml/kg and L3 for both 0.75 and 0.5ml/kg respectively. The difference between 0.5ml and 1.0ml/kg was statistically significant (P<0.01).

Table 3 shows the types of surgery allocated, on clinical grounds, to each of the volume groups.

Table 4 shows the median ± SEM height of spread of caudal solution. There was no statistically significant difference between incremental groups but a statistically significant difference (p<0.01) between the 0.5 and the 1 ml/kg groups.

Discussion

This study shows that increasing the volume of caudal injectate in the range 0.5-1.0ml/kg (100% increase) has only a modest effect on the anatomical spread of the solution in the epidural space. This small increase may also be of little clinical significance as even when the volume is doubled from 0.5 ml/kg to 1 ml/kg the anatomical spread can only be reliably increased by 1 vertebral level, from L3 to L2 in the group of children we studied, well below the 10 th thoracic (T10) vertebral level supplying the dermatome that includes the umbilicus. Although a simple and memorable formula based on weight has obvious advantages, a disadvantage may be that when applied to a general clinical population other important variables are ignored, thereby reducing it's predictive accuracy. It is know that a number of factors, including weight, influence the spread of caudal solutions and mathematical formulae using weight, age and height have been devised to predict the correct volume.

Satayoshi and Kamiyama suggest that the distance from C7 to the sacral hiatus in cm (D) is related to the desired volume (V) to clinically block to T4 by the formula V=D-13 [START_REF] Satoyoshi | Caudal anaesthesia for upper abdominal surgery in infants and children: a simple calculation of the volume of local anaesthetic[END_REF]. Weight was found to be a greater predictor of spread of analgesia in very young patients having caudals whereas age was a better guide in older children in the study by Busoni and Andreuccetti on 763 children(8). Schulte-Steinberg et. al. proposed a predictive formula in ml per dermatome to be blocked on the basis of age, this group found their formula to hold true for younger children aged 0-4 years as well as for older children up to the age of 12; although predicted volumes were relatively small using this method the authors suggested that on the basis of an efficacy model clinical analgesia several segments higher than the actual level of solution might be expected [START_REF] Schulte-Steinberg | Spread of extradural analgesia following caudal injection in children. A statistical study[END_REF].

Allocation of patients to their respective volume group by type of surgery rather than by chance should not introduce any bias. Had we set out to determine the clinical efficacy of the block as an outcome measure then clearly one might unwittingly assign more painful surgeries to the higher volume groups, for example. However we were simply interested in the physical spread of the caudal solution as our primary outcome and since this was measured before surgery commenced, should be independent of the intended operation.

The relationship between anatomical height of caudal solution and the extent of analgesia is not fully understood. We did not attempt to measure clinical efficacy of the caudal block in this study as it is notoriously difficult to accurately measure the extent of a regional sensory block in anaesthetised children. However, in clinical practice the aim is always to achieve reliable spread of the solution to the vertebral level corresponding to the required nerve root supply. Indeed, studies that have compared the efficacy low volume high concentration caudal injections with high volume low concentration techniques have tended to favour higher volumes, this implies that spread of the solution to the required level is desirable [START_REF] Hong | A comparison of high volume/low concentration and low volume/high concentration ropivacaine in caudal analgesia for pediatric orchiopexy[END_REF][START_REF] Verghese | Caudal anesthesia in children: effect of volume versus concentration of bupivacaine on blocking spermatic cord traction response during orchidopexy[END_REF].

There are several reasons why our results might reflect a lower height of block spread than the traditional formulae would predict. At approximately 10 minutes post injection (when the X-ray image was obtained), the caudal level. Also, there may be 'pharmacological spread' that is not visible to the naked eye on X-ray examination, representing significant numbers of local anaesthetic molecules reaching higher vertebral levels.

Hong et al [START_REF] Hong | A comparison of high volume/low concentration and low volume/high concentration ropivacaine in caudal analgesia for pediatric orchiopexy[END_REF] provides an interesting comparison to our study in that their 1 ml/kg group reached a median height of T11, only 1 vertebral level higher than our median for the 1ml/kg group, with their range being T8-L2. This is not incompatible with our findings. What we are stressing is that to choose a mean height block as the intended dermatomal cover is to condemn half the children (by definition of the mean) in that group to a block lower than the mean and thus potentially to a painful experience.

The 90 th centile is at L2 in our study and thus 2 vertebral levels lower than the median. Essentially we are suggesting that when it comes to a potentially painful dermatome it is better to be 90% sure of reaching it than 50%.

Whilst there may be a statistically significant difference between the 0.5 and 1 ml/kg groups, a difference of only one vertebral level is unlikely to be clinically significant unless the surgical site happens to span those dermatomes.

In conclusion, we found that in order to be confident that the L2 vertebral level is reached in 90% of cases in a group of children aged from 1-7 years, then a volume of at least 1ml/kg is required.

  aged between 1 and 7 years scheduled for elective

  the requirements for epidurography, furthermore Omnipaque has a long clinical history of safe use in studies on both caudal and epidural injection.The caudal was administered via a 22G abbocath intravenous catheter with the child in the left lateral position. The catheter was placed by a Consultant Paediatric Anaesthetist in the conventional manner under asceptic conditions once the child was anaesthetised. An Alaris Asena infusion pump delivered the pre-determined dose of caudal injectate at a standard rate of 600 ml/hr (i.e. 1 ml/6sec).

  A single Consultant Radiologist (DR), blinded to the volume administered, reported the extent of spread as the highest whole vertebral level visualised (fig 1). Patient weight, height, age and the time from caudal injection to X-ray were recorded. Results were analysed using GraphPad Prism 5 ® for windows. Parametric data were analysed using t-tests ± ANOVA as appropriate. Kruskal-Wallis with Dunn's Multiple comparison test was used to compare the 3 volume groups.

  be creeping cranially and not yet have reached its highest

  Payne K, Heydenrych JJ, Martins M, et al. Caudal block for analgesia after paediatric inguinal surgery. S Afr Med J 1987; 72: 629-630.

Figure 1 .

 1 Figure 1. A typical X-ray obtained with the radio-opaque dye clearly visible

Figure 1 FFigure 2 .

 12 Figure 1

Figure 2 1

 2 Figure 2

Table 1

 1 Demographics of each group by height, weight, age and time from injection to X-ray.

	Table 1
		Mean height	Mean weight	Mean age	Time (mins)
		(cms) ± SD	(Kg) ± SD	(Years) ± SD	from Caudal
	Volume
	(ml/Kg)				injection to X-
					ray
	0.5 (n=15)	102 ± 11.47	14.60 ± 3.36	3.51 ± 1.27	11.47 ± 6.9
	0.75 (n=15)	95 ± 14.50	14.53 ± 4.36	2.74 ± 1.33	11.00 ± 4.2
	1.0 (n=15)	104 ± 13.30	15.73 ± 3.73	3.55 ± 1.80	11.87 ± 3.9

Table 2 .

 2 The number of blocks and proportion in each group achieving vertebral levels in the range L4-T12.

	Table 2
	Volume	0.5 ml/kg	0.75 ml/kg	1 ml/kg
		(n=15)	(n=15)	(n=15)
	Vertebral level
	T12	2 (13%)	2 (13%)	5 (33%)
	L1	2 (13%)	6 (40%)	8 (53%)
	L2	4 (27%)	10 (67%)	14 (93%)
	L3	14 (93%)	14 (93%)	15 (100%)
	L4	15 (100%)	15 (100%)	15 (100%)

Table 3

 3 The different types of surgery allocated, clinically, to each of the Table4shows the median ± SEM height of anatomical spread for each caudal volume.

	caudal volume groups
	Table 3
	Volume of Caudal (ml/kg)	Type of surgery
	0.5	Wire removal (foot), circumcision (x8),
		hypospadias (x2), ureteric re-implant,
		excision glans cyst, anal dilatation,
		injection bladder neck.
	0.75	Inguinal hernia (x2), orchidopexy (x6),
		hypospadias (x4), club foot,
		circumcision, release buried penis.
	1	Pyeloplasty (x3), hypospadias (x5),
		orchidopexy (x5), inguinal hernia, left
		hip metalwork removal.

Acknowledgements

Our thanks to Professor Tim Cole in the Institute of Child Health for help with statistical advice and support.


Loïc Lestand

Paul Force

Gérard Montarou

Nicoleta Pauna

PCCF RI

Keywords: parameters available in Geant4 . . . . . . . . . . . . . . . . . Hadrontherapy, in beam TOF-PET, β + emitter, prompt γ, QMD, Fermi BreakUp, Evaporation, Geant4

The different steps of the G4QMDReaction::ApplyYourself method in Geant4 . 3.2 α, β, γ, ρ 0 , L and E sym parameters available in Geant4 . . . . . . . . . . . . . 3.3 Total effective potential parameters available in Geant4 . . . . . . . . . . . . 3.4 Ground state

List of Figures

Introduction

The discovery of the X and γ rays respectively in 1895 by Wilhem Röntgen, and in 1896 by Henri Becquerel, initiate a new field of activity which are ionizing radiation medical applications. Physics of ionizing particles has been widely applied for medical imaging as well as for therapy.

Different strategies are currently used for the cancer treatment : surgery, radiotherapy and chemotherapy. They could be used separatly or combined to adjust the best treatment to each patient. Among all the radiotherapy strategies, one can distinguish external radiotherapy by means of photon beams and/or electron beams, from internal radiotherapy currently named brachytherapy. The availability of new radiopharmaceutic products made of new isotopes are now more and more used for medical diagnosis and also for therapeutic purpose. For instance internal targeted radiotherapy (ITR), which consists in irradiation of small disseminated tumour lesions using injected radiopharmaceuticals has been for a long time successfully used for differentiated thyroid carcinoma treatment. New efficient radiopharmaceutical products for tumour targeting become available, making feasible such successful ITR treatment for other endocrine tumours.

The field of external radiotherapy has been deeply developed for the last 20 years. The aim of each improvement is the conformation to the treatment volume to optimize the dose deposition and to improve the healthy tissue sparing. Among all the techniques available today, the most remarkable are Intensity Modulated RadioTherapy (IMRT), storeotactic radiosurgery or Volumetric Modulated Arc Therapy (VMAT).

However, despite all these improvements, some cancers remain uncurable by conventional radiotherapy. The first type of cancers is the deep seated tumors in close vecinity to organs at risqs (OAR). The main difficulty in this case is the geometrical conformation of the dose deposition ; it requires a very precise balistic since OAR are very sensitive to a high dose deposition. Otherwise the surounding OAR could be over irradiated and the treatment would fail. This type of tumours requires the dose deposition to be well segmented.

11

The second type of uncurable cancers comes from radioresistant tumors. The treatment efficiency relies on the fact that after irradiation healthy tissue repair is better than malignant tissue repair. By fractionnating the whole treatment, this repair behavior induces a cell survival differential effect. This differential effect is the basis of the treatment effectiveness. In the case of radioresistant tumors this differential effect is unfavorable to healthy tissues. This type of tumors requires a more effective radiotherapy technique.

In the middle of the 40s Robert R.Wilson published an article called : "Radiological Use of Fast Protons" [START_REF] Wilson | Radiological use of fast protons[END_REF]]. Point out two arguments of Wilson's paper :

The proton proceeds through the tissue in very nearly a straight line, and the tissue is ionized at the expense of the energy of the proton until the proton is stopped. The dosage is proportional to the ionization per centimeter of path, or specific ionization, and this varies almost inversely with the energy of the proton. Thus the specific ionization or dose is many times less where the proton enters the tissue at high energy than it is in the last centimeter of the path where the ion is brought to rest.

The way the protons loose their energy in the interacting middle is greatly favorable for the treatment of deep seated tumors in close vecinity of OAR. The second developped argument is :

It is well known that the biological damage depends not only on the number of ions produced in a cell, but also upon the density of ionization. Thus the biological effects near the end of the range will be considerably enhanced due to greater specific ionization, the degree of enhancement depending critically upon the type of cell irradiated.

So the ionizing radiation efficiency depends on the density of ionization which increases at the end of the ion range. So the treatment of radioresistant cells could be achieved by means of ions. Ion therapy or hadrontherapy opens the way to the treatment of deep seated tumors and/or radioresistant cancers.

The in situ dose delivery monitoring is a crucial issue of these treatments quality control since ion range and therefore dosimetry is very sensitive to anatomical changes during the whole treatment period as well to the patient positionning and to the stoechiometric calibration of the CT images used as virtual patient for the dosimetric calculations. Figure 1: Geant4 v 9.2 simulation : dose profile from 12 6 C 6+ at 200 MeV/u in water equivalent material for different densities. Even a small variation in the density induces a significant modification of the range.

The dose deposition is induced by the electromagnetic processes. In addition, a lot of secondary particles and fragments are produced through nuclear collisions of the primary ions beam. Among the whole type of fragments, one can distinguish β + emitters which are essentially 10 C, 11 C and 15 O. They are short lived radioisotopes: T10 C = 20s, T11 C = 20min and T15 O = 2min. During irradiation, the β + emitters concentration reach an equilibrium. The induced β + activity measurement using PET technique could be one of the possible hadrontherapy treatment control principle.

However, a large amount of γ from nuclear deexcitation are superimposed over the 511 keV annihilation γ. Moreover, those nuclear γ are emitted within few keV and few MeV energy range. Their emission probability is deeply correlated to the primary ions beam temporal structure. Those prompt γ could induce a lot of random coincidences which are considered as background noise regarding the PET technique. To reduce this prompt background noise an acquisition out of ion spill out would be suitable. This technique CONTENTS was used by [START_REF] Alexandre | Optimization of In-Beam Positron Emission Tomography for Monitoring Heavy Ion Tumor Therapy[END_REF]] at GSI. The PET acquisition was performed between each 2s spill extraction of the GSI synchrotron beam structure. However, the number of true coincidences is reduced too.

Recent experiments performed at Ganil in Caen with a cyclotron, have demonstrated, given a 80 ns period accelerator signal, that acquisition with a veto applied on the prompt events allows the extraction of the 511 keV annihilation γ from background noise. As a consequence, the conception of a PET device dedicated to the hadrontherapy requires to find a compromise between different parameters such as the beam temporal structure and the detection device performances.

The acquisition of true event will be improved thanks to the design of fast electronics with a high bandwidth and without dead time despite the high rate of prompt γ.

Detectors and dedicated electronics specifications rely on accurate simulations of count rate and bandwidth assessment. Simulations have to reproduce β + emitters yields (true signal) and prompt γ yields (background noise).

Sections 1 and 2 will describe the physics of nuclear collisions at hadrontherapy energy range. More specifically, section 1 will consider nuclear dynamics whereas section 2 will be about the statistical deexcitation of fragments. Section 3 will deal with the implementation in Geant4 of both aspects of the physics of nuclear collisions. Section 4 will present β + emitters yield simulations and comparison with experimental data.

Finally, section 5 will report on the prompt γ yield study. Up to now, Geant4 reproduces the β + production rate quite well but overestimates the prompt γ rate up to a factor of 12 [START_REF] Foulher | Monte carlo simulations of prompt γ emission during carbon ion irradiation[END_REF]].

Chapter 1

Nuclear collisions dynamics

It is commonly assumed that nuclear collisions can be split into three different phases [I. Pshenichnov et al. (2009)].

During the first phase of the reaction, depending on the incident energy of the projectile and the impact parameter, a more or less large amount of nucleons of the projectile and the target (participant nucleons) will interact through binary collisions and/or mean field interactions. The theoritical description of this step is not trivial. At high incident energy (> 500 MeV/u typically), De Broglie wavelength and free mean path of any nucleon involded in the collision is short enough to avoid any quantum mechanics description. For instance, the Pauli exclusion principle is negligeable given that the free mean path is very short. During this phase of the collision, it is assumed that nucleons interact mainly through binary collisions which are governed by kinematics parameters such as the minimal distance of interaction d min = σ nn /π, where σ nn is the nucleon-nucleon interaction cross section. At lower incident energy (< 100 MeV/u typically) quantum mechanics cannot be neglected. Thus nucleons interact only through mean field interactions. At intermediate energies both descriptions are required.

The second stage of the collision is associated with the pre-equilibrium dynamics. During this phase of the reaction, the highly excited system is split into an intermediate set of excited fragments. As the dynamical stage leaves the nuclear matter very excited, energetic nucleon interact through collisions that could lead to particle-hole creation. Then the nuclear matter condensate which means that the nucleons which are close enough to each others are considered belonging to the same fragment [START_REF] Machner | Pre-equilibrium decay into complex particle channles analysys of double differential cross sections[END_REF]] while the others are considered as free nucleons.

The third stage of the reaction is the deexcitation of the resulting fragments. During this stage, all the excited fragments produced in the preequilibrium stage will deexcite via different deexcitation channels. At the end of this stage, all the remaining fragments are supposed to be in ground state.

The complete theoretical description of heavy ions collisions is not an easy task, because of the great number of implied particles. To mitigate this disadvantage, one generally use simplifying assumptions which make possible to define less complex ideal models. These models can be classified summarily in two main categories : the macroscopic and the microscopic models.

Macroscopic models : general concepts

The hydrodynamic models represent a first example of macroscopic models. They are based on the following assumption: the mean free path of the nucleons, λ, is very small in comparison with the dimensions R of the system. Nuclei are then identified with fluids, and the collision is treated by solving the fluid transport equation, where the nuclear matter equation of state is explicitely involved.

The thermodynamic models are another example of "macroscopic" models. The validity of these models corresponds to the following inequality R > λ > r where r is the range of the nuclear forces. This inequality means that the system can be treated like a perfect gas of nucleons in a state of thermodynamic equilibrium. The "fireball" model was one of the first model of this kind : it is considered that all the participating nucleons contained in the intersection volume of the two reacting nucleus (defined by a geometrical cut of the projectile in the target) is in a state of thermodynamic equilibrium. Moreover the impulsion distribution of the nucleons emitted in the the centre of mass of the "fireball" follows the Boltzmann equation of an ideal gas.

The use of statistical thermodynamics concepts to describe the fireball implies the establishment of a chemical equilibrium. That is a chemical balance between the various components of the fireball. This means that the creation rate of the fragments is balanced by their destruction rate until the chemical freezing. The chemical freezing is the begining of the evaporation phase, where the system is very diluted. At that time of the reaction, the number of interactions between the various components of the fireball is then sufficiently small so that the concentration of the various types of fragments evolve no more.

Microscopic models : Quantum Molecular Dynamics approach

The Quantum Molecular Dynamics (QMD) is one of the most used microscopic model that could describe heavy ion collision at intermediate energy. It is based on molecular dynamics models with some quantum mechanics considerations such as the Pauli exclusion principle. Such types of models are dealing with the time evolution of any N body system. The N individual components of the system mutally interact all together. QMD models take into account binary collisions between nucleons as well as mean field interactions.

Physical concepts

In this section we briefly present the QMD models main concepts.

The first concept is that any nucleon is represented by a gaussian wave packet which depends on two parameters (x α , p α ) [J. Aichelin (1993)] :

φ α (x, t) = 2L π 3/4 e -L(x-xα(t)) 2 e ixpα(t) e -ip 2 α (t)t/2m , ∀α ∈ [1, ..., A p + A t ] (1.1)
The parameter L, which is related to the extension of the wave packet in phase space, is fixed all along the reaction [START_REF] Hartnack | Modelling the many-body dynamics of heavy ion collisions : Present status and future perspective[END_REF]].

The second concept is that, neglecting the Pauli exclusion principle, the total wave function is the product of the wave packet of each nucleon :

φ = α φ α (x, x α , p α , t) , ∀α ∈ [1, ..., A p + A t ]
The initial value of the parameters are chosen in such a way that the ensemble of the A p + A t nucleons gives an appropriate density distribution as well as an appropriate momentum distribution of the projectile and target nucleus. The evolution of the system is calculated by means of a generalized variational principle starting out from the action.

S = t 2 t 1 L[φ, φ * ]dt (1.2)
with the Lagrange functional L :

L = φ ih d dt -H φ (1.3)
The total time derivative includes the derivation with respect to the parameters. The evolution is obtained by the requirement that the action is stationnary under the allowed variation of the wave function.

δS = δ t 2 t 1 L[φ, φ * ]dt = 0 (1.4)
Performing the variation with the total N body wave function we obtain for each parameter λ an Euler Lagrange equation :

d dt ∂L ∂ λ - ∂L ∂λ (1.5)
For the coherent states and an Hamiltonian of the form H = i T i + 1 2 ij V ij (with T i and V ij refer to the kinetic and potential energies respectively), the lagrangian and the Euler Lagrange function can be easily calculated :

L = α   -ẋα p α + β V αβ + 3 8Lm   (1.6)
This leads to the two following equations of motion

ẋα = p α m + ∇p α β V αβ (1.7) ṗα = -∇x α β V αβ (1.8)
Thus the variational approach has reduced the N body Schrödinger equations to a set of 6N differential equations for the parameters which can be solved on present days through computational procedures.

Description of the Hamiltonian

The Hamilonian of the system has to deal with nuclear interactions which is the crucial ingredient to reach a good description of the nuclear collision dynamics. It can be written as the following equation :

H = i T i + 1 2 ij V ij = i p 2 i 2m + 1 2 ij V ij (1.9)
During the time evolution the nucleons interact via nuclear interactions (Yukawa potential), charged particles interactions (Coulomb potential) and mean field interactions (Skyrme potential) [J. Aichelin (1993)].

The potential is commonly written as :

V ij = V Coul + V Y uk + V (2) loc + V (3) loc + V M DI
(1.10)

Let's describe each term :

• The Coulomb potential is V Coul (charged particles interactions)

V Coul = i j =i V Coul ij (1.11) with V Coul ij = e 2 r ij erf r ij 2 √ L (1.12) • The Yukawa potential is V Y uk (nuclear interactions) V Y uk = t3 e -|r 1 -r 2 |/m |r 1 -r 2 | /m (1.13)
• The two body local Skyrme interaction is V

(2) loc :

V (2) loc = t 1 δ(r 1 -r 2 ) (1.14)
• The three body local Skyrme interaction is V

(3) loc :

V (3) loc = t 2 δ(r 1 -r 2 )δ(r 1 -r 3 ) (1.15)
• V M DI is an additionnal momentum dependent interaction, which is function of two parameters : t 4 and t 5 . They are determined by a fit to the measured optical potential.

The parameters t 1 and t 2 are determined from the infinite nuclear matter state which is defined by its average binding energy E/A = -15.75M eV and the normal nuclear matter density ρ = ρ 0 = 0.17 f m -3 . As a consequence, the local interaction is a function of the density only.

V loc (α, β, γ) = α ρ ρ 0 + β ρ ρ 0 γ (1.16)
The parameter γ is introduced to take into account the nuclear matter compressibility. The choice of those three parameters is depending on the nuclear matter equation state (EOS) [J. Aichelin (1993),J. Aichelin (1991), [START_REF] Polanski | Development of methods for simulation of electronuclear process[END_REF]] :

E A = 3p 2 F (ρ) 10m + α 2 ρ ρ 0 + β 1 + γ ρ ρ 0 γ (1.17)
and on the nuclear matter compressibility : All four EOS's give the same ground state binding energy i.e., E/A=-15.75 MeV at ρ = ρ 0 , but they differ drastically at high densities. Here the hard EOS leads to much more compressional energy than the soft EOS at the same density. For infinite nuclear matter at rest the inclusion of the momentum dependent interactions does not change the compressional energy as shown in figure 1.2. We see no difference between the H and HM EOS. Considering heavy ion collisions, the additionnal repulsion due to the initial separation of projectile and target in momentum space shifts the curve for SM and HM interactions to higher energies.

K = 9ρ 2 ∂ 2 E/A ∂ρ 2 ρ=ρ 0 (1.18)

Time evolution

The time evolution of the system in the space phase proceed in two different ways. The nuclear matter equation of state : the density dependence of the energy per particle in nuclear matter at temperature T=0 is displayed for 4 different sets of parameters [J. Aichelin (1993)]

The first one consists of the propagation of the centroid of each gaussian wave packet associated to each nucleon by solving the equations of motion :

ṙ α = {p α , H} ṗα = {r α , H} (1.19)
The second way consists of two body collisions. A pair of nucleon suffers a collision if the distance between the centroids becomes smaller or equal to d min = σ nn /π, where σ nn is the nucleon-nucleon interaction cross section. Thanks to this binary collision process, it is possible to introduce a Pauli blocking which allows to take the fermionic nature of the nucleons into account. The collision between a pair of nucleon is allowed if the new states in the phase space are not already occupied. Otherwise, the collision is blocked and the nucleons propagate away through the effective potential. The lower the energy, the more collisions are blocked. Thus for very low energy reactions the collisions are not so important. Anyway, QMD remains a phenomenological approach and has to be treated as such.

The differential equations describing the evolution of the centroid of each gaussian wave packet associated to each nucleon are solved by a discretization of the time into small ∆t intervals. The integration procedure determines the precision of the calculation. The procedure which is often chosen is the two step Euler method [START_REF] Polanski | Development of methods for simulation of electronuclear process[END_REF]] :

f isrt step :        r α n + 1 2 = r α (n) + ∆t 2 p α (n) m + ∇V p α n + 1 2 = p α (n) - ∆t 2 ∇V
(1.20)

second step :            r α (n + 1) = r α (n) + ∆t     p α n + 1 2 m + ∇V     p α (n + 1) = p α (n) -∆t∇V (1.21)
where r α (n) or p α (n) means position or impulsion at time n∆t. Depending on the initial conditions (incident energy, impact parameter, system size), the time propagation is stopped after 100-200 fm/c.

Clustering and deexcitation

At the end of the propagation procedure, a cluster reasearch is performed [START_REF] Polanski | Development of methods for simulation of electronuclear process[END_REF]]. The nucleons are spread out in the phase space but two nucleons belong to the same cluster if their relative distance in positions space is lower than a certain distance R coal . R coal depends on the fragments formation coalescence time. Its value is different if the coalescence time is fixed at the chemical freezing time of the reaction (15-30 fm/c [J. [START_REF] Kinet | Pion production in central high energy nuclear collisions[END_REF]]), or in the final time of the collision (300 fm/c in our case). Such a theory of coalescence for the fragment formation also required that two nucleons belong to the same cluster if their relative momentum |∆ -

→ p | = | - → p 1 -- → p 2 | do not
exceed a maximal relative momentum p 0 . This parameter was evaluated to 160 MeV/c [START_REF] Westfall | Light particle spectra from 35 mev/nucleon 12 c-induced reactions on 197 au[END_REF]]. One can approximate R coal value as a function of the coalescence time t through the following relation :

R coal = r gc + p 0 m (t -t gc ) (1.22)
where, m is the mass of the nucleon (m = 938.3M eV ), t gc is the time of the chemical freezing and r gc is the chemical freezing corresponding coalescence radius (r gc = 3 fm).

For each produced fragment, one determines the mass, total momentum and energy E int , defined as the sum of energies of each nucleon. The total momentum of fragment, P f rag is given by the vectorial sum of the momenta of the nucleons which compose this fragment :

- → P f rag = i - → p i (1.23)
Each fragment has thus an excitation energy E * given by

E * = E int -E = N i p 2 i + m 2 -E (1.24)
where E is the total energy of the fragment. This energy is calculated starting from the total momentum of the fragment, P f rag , and the mass of this fragment estimated by a traditionnal formula of mass, p i is the momentum of the i th nucleon of the fragment and m is the nucleon mass. Some of the fragments produced by coalescence do not exist in ground state, or are too far away from the line of stability of the table of isotopes (too large number of neutrons and protons). These fragments will desintegrate in lighter stable nuclei.

Some other fragments, because of the momenta of the nuclei which compose them, are affected a too large energy which makes them unstable. They will evacuate their excitation energy through emission of photons or light particules.

Each fragment produced by coalescence is thus likely to undergo a process of desintegration or evaporation to form new fragments in their ground state. This define the third step of the whole collision process which is described in the next section.

Chapter 2

Statistical deexcitation

In the previous section, we have briefly described, how excited secondary fragments are produced from a collision of projectile and target ions. In this section we are going to discuss about the models which can be used to deal with the fragments deexcitation.

For excitation energy of about 1 MeV/u, the main deexcitation mechanisms are evaporation of particles and fission for heavier and slightly excited fragments. At these excitation energies, the nuclear density is close to the saturation density ρ 0 ∼ 0.17 fm -3 and the global nuclear properties are well described by the liquid drop model.

For larger excitation energies, ∼3 MeV/u, quasi continuous emission of nucleon or light clusters is expected.

At higher excitation energies (5)(6)(7)(8) MeV/u) a more explosive deexcitation mechanism is more relevant to describe the equilibrium return [START_REF] Bondorf | Statistical multifragmentation of nuclei[END_REF]].

There are many models to describe the deexcitation of excited fragments and the emission of intermediate mass fragments (IMF). The statistical multifragmentation model is based on a statistical description of those phenomena. J Bondorf et al define clearly in a Physics Report [START_REF] Bondorf | Statistical multifragmentation of nuclei[END_REF]] the main idea of the statistical multifragmentation: the statistical approach which despite of simplicity is very suitable for describing highly excited nuclear systems. Its main idea is that at high excitation energy a very large number of degrees of freedom is involved in the process and the probabilities of different decay channels are mainly determined by statistical weights rather than by the detailed dynamics of the process. This approach makes possible a full description of the highly excited system of hadrons and nucleon clusters in an economic way. It gives wide opportunities to simply imple-25 ment many specific features of nuclear systems. But the price for this simplicity is the assumption of statistical equilibrium.

Physical concepts

Consider now only the mechanisms involved in the fragments deexcitation, and assume that each excited fragment is characterized by a volume V, an excitation energy E * , a charge Z 0 and a mass A 0 . Due the excitation energy, the internal pressure is increasing and the volume of the excited fragment will expand. As long as the fragment volume increases, the distance between each nucleon is raising so that the attractive force is not able to maintain any more the cohesion of the fragment. This is followed by the breakup of the fragment into n smaller fragments. For excitation energies considered ( ∼ 1 -10 MeV/u), the number of decay channels is very large. A statistical description of all possible decay channels is relevant. Each channel is weighted by its probability of occurrence. The probabilities of each fragment deexcitation channel are calculated and randomly chosen.

Depending on the excited system size, different models are available : Fermi BreakUp model for light fragments and Evaporation model for heavier systems.

The Fermi Break-Up model

For the light clusters (A≤16) even a small excitation energy is comparable with the binding energy [START_REF] Bondorf | Statistical simulation of the break-up of highly excited nuclei[END_REF]]. Then, the main break-up mechanism is assumed to be an explosive like mechanism. The probability of each channel can be written as

W ∝ E 3n -5 2 avai (2.1)
where E avai is total energy available for the set of fragments at the break-up time. This energy is expressed as

E avai = E * + M (A, Z)c 2 - n i=1 m i c 2 + ǫ i -B Coul (2.2)
E * + M (A, Z)c 2 refers to the total energy of the fragment in the center of mass (E * is the excitation energy and M (A, Z)c 2 is the mass energy) .

m i c 2 + ǫ i refers to the total energy of each subfragment in the center of mass (ǫ i is the excitation energy and m i c 2 is the mass energy). B Coul is the Coulomb barrier of the splitting system.

B Coul = 3 5 e 2 r 0 1 + V f r V 0 -1/3 Z 2 A 1/3 - n i=1 Z 2 i A 1/3 i (2.3)
V f r corresponds to the freeze out volume that may be larger than the volume of the nuclear matter at the saturation point. The final state fragments are supposed to be in their ground or quantized excited state.

As a matter of fact, by forcing the excitation energy of the resulting fragments to be quantized, the whole number of channel is significantly reduced. This is a way to build nuclear level during the deexcitation process.

The evaporation model

Particle evaporation

For heavier fragments (A>16 typically), the evaporation process through the emission of lighter clusters is assumed to be the main deexcitation mechanism. This model is based on the Weisskopf and Ewing evaporation scheme [START_REF] Bondorf | Statistical multifragmentation of nuclei[END_REF], [START_REF] Weisskopf | On the yield of nuclear reactions with heavy elements[END_REF]]. It is relying on the following equality :

P p→f = P f →p (2.4)
The probability of a particle p emits a particle l c to become a particle f is equal to the probability of a particle f captures the same particle l c to become the particule p. P f →p is proportional to the cross section σ inv of the inverse reaction (A ′ , Z ′ )+l c = (A, Z). Therefore, the probability that a parent nucleus p with an excitation energy E * emits a particle l c in its ground state with kinetic energy ǫ kin is [START_REF] Weisskopf | On the yield of nuclear reactions with heavy elements[END_REF]]:

P lc ∝ σ inv ρ f (E max -ǫ kin ) ρ p (E * ) (2.5)
where ρ p (E * ) is the excited level density of the evaporating nucleus, ρ f (E maxǫ kin ) that of the daugther (residual) nucleus after emission of a fragment l c and E max is the maximum energy that can be carried by the ejectile.

The available evaporating channels are : proton channel, neutron channel, deuteron channel, triton channel, 3 He channel, α channel. This list can be extended to heavier nuclei, e.g. 18 O, and to photons.

Photon evaporation

The photon emission processes can be added to the deexcitation procedure. Those photons are called statistical photons since their emission is weighted by probabilities. As an exemple, they can come from giant resonnances, resulting from nucleon collective motion. The nuclear matter undergoes compression and expansion mechanisms which set in motion the nucleons which lead to collective motions.

The main assumption considers these photons coming from a particular type of giant resonances the so called giant dipole resonances (GDR) [START_REF] Iljinov | Phenomenological statistical analysis of level densities decay of lifetimes of excited nuclei[END_REF]].

The probability to evaporate γ in the energy interval (ǫ γ , ǫ γ + dǫ γ ) is :

W γ (ǫ γ ) = 1 π 2 (hc) 3 σ γ (ǫ γ ) ρ (E * -ǫ γ ) ρ (E * ) ǫ 2 γ (2.6)
where σ γ is the γ absorption cross section and ρ is the nucleus level density. σ γ is defined by :

σ γ = σ 0 ǫ 2 γ Γ 2 GDP ǫ 2 γ -Γ 2 GDP 2 + (ǫ γ Γ GDP ) 2 (2.7)
with σ 0 = 2.5A mb, Γ GDR = 0.3E GDP and E GDP = 40.3A -1/5 M eV . Then the total probability is calculated as follows

W γ = 3 π 2 (hc) 3 E * 0 σ γ (ǫ γ ) ρ (E * -ǫ γ ) ρ (E * ) ǫ 2 γ dǫ γ (2.8)

Nuclear radiative transitions

There is an other photon emission phenomenon which is the nuclear radiative transitions.

Just after the collision dynamical stage, we can imagine that nuclear levels do not exist yet. While the nuclear matter is cooling down, nuclear levels appear. For exemple, particles coming from the Fermi Break Up process are in their ground state or in a quantized exited nuclear state. Therefore, the nuclear transition is a low excitation energy process. So we can distinguish photons coming from giant resonances which are treated as a statistical deexcitation channel and photons coming from nuclear transitions which are treated independently from tabulated data and supposed to occur at the end of the deexcitation procedure.

Discrete photons production is governed by the nuclear level transition rules. Each excited level is characterized by a mean lifetime τ and the transition probability per unit time is inversly porportional to τ . Moreover, to reproduce the nuclear magic numbers, we need to introduce a strong spin orbit coupling for each nucleon. The existence of this interaction lead to several transition rules since there are energy and angular momenta exchanges for each transition. Those rules are essentially governed by the quantum mechanics and deal with the addition of angular momenta. The figure 2.1 shows schematically the different type of photons which are emitted through nuclear radiative transitions and photon evaporation. In conclusion, the Fermi Break-Up model is a good way to treat the main part of fragments. Nuclear transitions model can be used as additionnal model since excitation energy of the fragments coming from Fermi Break-Up splitting procedure may coincide with nuclear level.

Ground state S 0 S 1 S 2 S 3 Continuum Binding energy E ∞ S i -→ S j discrete radiative transition. C i -→ S j continuous radiative transition.

Chapter 3

Geant4 implementation

Geant4 is a Monte Carlo toolkit dedicated to the simulation of the passage of particles through matter. Its domain of application include high energy, nuclear and accelerator physics. It has been extensively used for high energy detector conception and validation. More recently Geant4 is applied for studies in medical and space sciences.

Such a tool has a lot of advantages in comparison with other Monte Carlo code: it is an open source that is also very well documented. The two main reference papers for Geant4 are published in Nuclear Instruments and Methods in Physics Research [J. Sulkimoah et al. (2003)], and IEEE Transactions on Nuclear Science [J. Allison et al. (2006)]. A large community worldwide uses the code and induces so a continuous process of validation. Regular upgrades of the code are distributed (twice a year) through the Geant4 collaboration and its website (http://geant4.cern.ch/). There are many courses and tutorial available as well an efficient online user support. This section is dedicated to the description of the implementation of the physical models for ions collision in the Geant4 toolkit. Hereafter we only consider versions 9.2 till 9.4 of the code.

QMD package description

One of the most important ingredient for Geant4 simulation is the "physics list". It defines all the particles and physical processes that will be handled by the code in the user defined application. Geant4 offers a large variety of nuclear processes embedded in several packages. Two different packages are useful for hadrontherapy simulations: Binary Cascade and QMD packages. The term package refers to the association of a collision and deexcitation models. QMD package contains QMD and statistical deexcitation models.

QMD

QMD implementation

The table 3.1 describes the QMD model in Geant4.

The collision impact parameter

The collision is performed given an impact parameter b. The maximum impact parameter b max is calculated from the Shen cross section. The general Shen cross section are given by :

σ R = πR 2 1 - B E CM (3.1)
where R is the interaction radius and B is the nuclear-nuclear interaction barrier. The Shen cross section gives better agreement at energies below 30 MeV/u [START_REF] Shen | Total reaction cross section for heavy-ion collisions and its relation to the neutron excess degree of freedom[END_REF]]. The maximum impact parameter is calculated from cross section :

b max = σ R /π (3.2)
Then, the impact parameter is randomly chosen according to a uniform distribution:

b reaction = √ rb max (3.3)
G4QMDReaction::ApplyYourself(projectile,target)

• Maximum Impact parameter b max calculated from Shen Cross Section [START_REF] Shen | Total reaction cross section for heavy-ion collisions and its relation to the neutron excess degree of freedom[END_REF]]

• Impact parameter of the reaction chosen randomly, b reaction ∈ [0, b max ]

• Space phase initialization :

-G4QMDGroundStateNucleus(G4int z,G4int a) : (r α , p α ) initialization -G4QMDMeanField meanFiel->SetSystem(projectile+target) : parametrizable mean fiel construction -G4QMDCollision collision->SetMeanField ( meanField )

• for ( G4int i = 0 ; i < 100 fm/c ; i++ ) meanField->DoPropagation(1 fm/c) Space phase temporal evolution collision->CalKinematicsOfBinaryCollisions(1 fm/c) Binary collisions + Pauli blocking

• Clustering -R coal = 4 fm • A set of excited fragments {A,Z,E * } Table 3.1:
The different steps of the G4QMDReaction::ApplyYourself method in Geant4

Space phase initialazation

The second step consists in initializing the centroid of each wave packet (r α , p α ) to build the initial ground state. Moreover, since QMD model is based on N body interactions through an effective potential, the definition of this effective potential is a crucial ingredient because it rules the dynamics of the system. The total effective potential (cf. section 1) is :

V ij = V Coul + V Y uk + V (2) loc + V (3) loc + V M DI (3.4)
The Geant4 class which handles the main parameters of the effective potential is G4QMDParameters. In the section 1, some parameters were introduced : t 1 and t 2 the two and three body Skyrme interactions, t 3 for the nuclear interaction and t 4 and t 5 for the additionnal momentum dependent interaction.

The parameters handled by the G4QMDParameters class for the effective potential definition are called differently. In fact, they are called c 0 , c 3 , c s and c l . The two first refer to the two and three body interaction contribution. c s , refers to the nuclear potential and c l refers to the Coulomb interaction. In the Geant4 toolkit, there is no momentum dependent interaction. Each parameter is defined as followed :

c 0 = α 2ρ 0 (4πL) 3/2 c 3 = β (γ + 1) (4πL) 3γ/2 2ρ γ 0 c s = E sym 2ρ 0 (4πL) 3/2 c l = c Coul 2 (3.5)
Those parameters are in fact function of more physical parameters which are given in the table 3.2. Finally, ground states of the projectile and the target are determined using the four paramters c 0p , c 3p , c sp ,c lp .

ρ 0 L E sym α β γ 0.168 f m -3 2 f m 25 M eV -219.418 M eV 165.327 M eV 1.33333
c 0p = 2c 0 c 3p = c 3 (γ + 1) c sp = 2c s c lp = 2c l (3.6)
The table 3.4 sums up the ground state parameters available in Geant4.

Time evolution

The propagation procedure is handled by the G4QMDMeanField::DoPropagation method and is performed following a two step Euler method. The total evolution time is set to 100 fm/c and the time step is set to 1 fm/c.

Two body collision and Pauli blockling

The collision procedure is managed by the CalKinematicsOfBinaryCollisions method from the G4QMDCollision class. It is based on the G4Scatterer Geant4 class. More precisely the method theScatterer->Scatter(kt1 ,kt2) returns a kinetictrackvector non dummy if the two nucleons, refered by kt1 and kt2 have scattered through elastic process. A Pauli blocking procedure is performed if the elastic collision is allowed.

Clustering

After 100 fm/c, the propagation procedure is stopped. The G4QMDMeanField:: DoClus-terJudgment method checks if two nucleons are close enough in the phase space. to be considered to belong to the same cluster. Two conditions are requiered : one distance condition (r min = 4 f m) and one momentum condition (p min = 0.0966 GeV ). p min is deduced from r min according to the Heisenberg uncertainty principle :

r min × p min 2 < hc (3.7)
The 1/2 multiplication factor is introduced to take the nucleons fermionic nature into account.

Statistical deexcitation implementation

General overview

The figure 3.2 illustrates the deexcitation processes as implemented in the QMD package of the Geant4 code. A>4 and A<17 and Z<9 G4FermiBreakUp::

BreakItUp(A ′ , Z ′ , E * ′ ) E * >3 MeV/u G4StatMF::BreakItUp(A ′ ,Z ′ ,E * ′ ) A>17 or Z>9 G4Evaporation::BreakItUp(A ′ ,Z ′ ,E * ′ ) A set of intermediate excited fragments {A ′′ ,Z ′′ ,E * ′′ } G4PhotonEvaporation::BreakItUp(A ′′ ,Z ′′ ,E * ′′ ) if (no G4PhotonEvaporation::BreakItUp(A ′′ ,Z ′′ ,E * ′′ ) forced γ emission, E γ =E * ′′ .
A set of stable fragments

Fermi BreakUp model implementation

The G4FermiBreakUp::BreakItUp(A,Z,E * ) method is the main method which handles the Fermi BreakUp deexcitation process. Given an excited fragment, all the break-up channels are evaluated, which means, for each break-up configuration the following quantity is calculated

E avai = E + M (A, Z)c 2 - n i=1 m i c 2 + ǫ i -B Coul (3.8)
where B Coul is the Coulomb barrier of the splitting system.

B Coul = 3 5 e 2 r 0 1 + V f r V 0 -1/3 Z 2 A 1/3 - n i=1 Z 2 i A 1/3 i (3.9)
If E avai is strictely positive the break-up probability is calculated. Then each channel probability is sorted and the break-up channel is randomly chosen among all available channels. Note that specific channels are also considered:

• 5 3 Li -→ α + p • 5 2 He -→ α + n • 8 4 Be -→ α + α • 9 4 Be -→ α + α + p
The secondary set of fragments are either in ground state or with some quantized excitation level which allows discrete γ transitions.

We will not describe the G4StatMF::BreakItUp since it is dedicated for very excited and heavy fragments. Those fragments are very rare in nuclear collisions involded in the hadrontherapy energy range.

Evaporation model implementation

Particles evaporation

The evaporation model used by default in the QMD package is the generalized evaporation model. This model handles emission of fragments heavier than α particles. Here are the particles which could be emitted through an evaporation process : p, n, d, t, 3 2 He, 4 2 He, and all the stable nuclei till 28 12 M g. In addition competitive fission and photon emission channels are considered.

The main method which handles the evaporation process is BreakItUp(A,Z,E * ) from the class G4Evaporation. This method calculates the probability of each channel according to the Weisskopf and Ewing evaporation scheme. For all the channels available, the probabilities are sorted and the channel is randomly chosen using all those probabilities.

Finally it is the BreakUp(A,Z,E * ) method which proceeds to the deexcitation procedure given an evaporation channel.

Unlike Fermi BreakUp model, fragment excitation energies after evaporation are not quantized since excitation energy is completly determined by the kinetic energy emission of evaporated fragment. This energy is of equal probability distributed between the Coulomb barrier and the maximal kinetic energy that can be carried by fragment. Thus secondary fragment excitation energy is continuously distributed between the emitted fragments.

Photon evaporation

If the photon evaporation channel is chosen, the BreakUp(A,Z,E * ) method from the G4PhotonEvaporation class handles the process. The photon evaporation is divided into two components : "discrete" evaporation through nuclear radiative transitions and "continuous" evaporation through the giant dipole resonance.

Discrete levels are tabulated from the first excited level till the last discrete level just before the continuum. If the excitation energy is lower than the last discrete level, refered as the maximum level, the method proceeds to evaporation of a "discrete" photon. Otherwise evaporation of "continuous" photon is performed.

If E * > E max
level , the probability to evaporate γ in the energy interval (ǫ γ , ǫ γ + dǫ γ ) is :

W γ (ǫ γ ) = 1 π 2 (hc) 3 σ γ (ǫ γ ) ρ (E * -ǫ γ ) ρ (E * ) ǫ 2 γ (3.10)
σ γ is the γ absorption cross section and ρ is the nucleus level density. σ γ is defined :

σ γ = σ 0 ǫ 2 γ Γ 2 GDP ǫ 2 γ -Γ 2 GDP 2 + (ǫ γ Γ GDP ) 2 (3.11)
with σ 0 = 2.5A mb, Γ GDR = 0.3E GDP and E GDP = 40.3A -1/5 M eV . Then the total probability is calculated following the equation :

W γ = 3 π 2 (hc) 3 E * 0 σ γ (ǫ γ ) ρ (E * -ǫ γ ) ρ (E * ) ǫ 2 γ dǫ γ (3.12)
If "discrete" evaporation is chosen i.e. E * < E max level , nuclear transitions are performed according tabulated data. If the excitation level is lower than E min level /2, the excitation is set to 0.

In each case, i.e. "continuous" and "discrete" evaporation, only one photon is evaporated. There is no a cascade.

Photon emission and forced γ emission

The G4PhotonEvaporation::BreakItUp(A,Z,E * ) method handles the last deexcitation procedure. It acts exactly like G4PhotonEvaporation::BreakUp(A,Z,E * ) method, however, photon emission is performed while fragments are excited. In case there is no nuclear level found in the nuclear data base, the deexcitation is forced through an artificial γ emission procedure. This procedure emits a γ with an energy equal to the residual excitation energy.

Chapter 4 β + emitter production

In the previous sections, the Quantum Molecular Dynamics model and its implementation in the Geant4 toolkit has been described. PET technique is based on the detection of a pair of 511 keV annihilation γ in a narrow time interval (coincidence). Such a technique could be used to the hadrontherapy treatment control with some in beam specific adaptations. Nevertheless, the design of such a device requires that one is able to accuratly evaluate the β + induced count rate. However, PET technique for hadrontherapy application will be used in conditions where, in the spill extraction, a lot of random coincidences coming from the prompt γ are superimposing over the 511 keV coincident γ. This requires to modify the trigger regarding the PET acquisition. Thus, one should also be able to predict such a background noise. Anyway, this section is dealing with the comparison between experimental and Geant4 simulated β + emitters production rates. Comparisons are performed in terms of integrated production rates and depth profiles shapes.

Description of the experimental data

Two different types of measurement are available : production rates and production rates versus depth.

Production rates

The production rate given a β + emitters nucleus i is defined as following :

Y i = N i N ion (4.1)
where N i is the number of produced β + emitters given N ion primary ions. These experimental data are reported in [START_REF] Parodi | On the feasibility of dose quantification with in-beam PET data in radiotherapy with 12 C and proton beams[END_REF]][I. Pshenichnov et al. (2006)]. One consider carbon ions interacting with a 300 mm × 90 mm × 90 mm PMMA target (PMMA is an equivalent tissue medium). The main resulting β + emitters nuclei are 10 C (T 10C = 2 s), 11 C (T 11C = 20 min) and 15 O (T 15O = 2 min). 10 C, 11 C come from target and projectile fragmentation whereas 15 O comes from target fragmentation only. Measurements where performed for three different energies : 212,12 MeV/u, 259.5 MeV/u and 343,46 MeV/u.

MeV/u

Nuclei

Experiment 10 C (8 ± 3) × 10 -3 11 C (10.5 ± 1.3) × 10 -2 15 O (2.1 ± 0.3) × 10 -2 259.5 MeV/u Nuclei Experiment 10 C (1.2 ± 0.3) × 10 -2 11 C (14.7 ± 1.6) × 10 -2 15 O (3.1 ± 0.4) × 10 -2

MeV/u

Nuclei Experiment 10 C

(1.5 ± 0.3) × 10 -2 11 C (19.9 ± 2.4) × 10 -2 15 O

(5 ± 0.4) × 10 -2

Table 4.1: β + emitters yields in PMMA

The table 4.1 gives an overview of the experimental production rates of β + emitters nuclei induced in a PMMA target and for different incident energies.

Depth profiles

Depth profiles measurement principles are reported in [START_REF] Priegnitz | A novel approach for predicting the positron emitter distributions produced during therapeutic ion irradiation[END_REF]]. They are obtained from reconstructed images of a PMMA target irradiation by carbon ions beam for different energies. Only three energies are reported here : 215 MeV/u, 260 MeV/u and 337 MeV/u. The reconstruction algorithm used is filtered back projection including random coincidences correction. Detection efficiency, solid angle and attenuation corrections are also taken into account. These data are shared by the wp6 ENVISION europeen project.

The experimental depth profiles are illustrated by the figure 4.1. Depth profiles show the spatial distribution of the main β + emitters nuclei. Moreover, by integrating the profiles one can determine the total production rates over the whole target. Once different experimental data sets have been introduced, we are going to present some comparisons with simulations. They were performed with QMD package since it gives better results for hadrontherapy applications [START_REF] Böhlen | Benchmarking nuclear models of fluka and geant4 for carbon ion therapy[END_REF]]. We have divided this comparison study into two steps. The first step is consisting in the simulation of the integrated production rates and their comparison with experimental data. We have considered the measured production rate introduced in the subsection 4.1.1 and the calculated production rates from the depth profiles introduced in the current subsection. The second step will report comparison between simulated and measured production rates versus depth in terms of shape (see section 4.3).

Total production rates : simulations versus measurements

Total production rates, according to the equation 4.1, were simulated with Geant4. A 90×90×300 mm 3 PMMA phantom was chosen as interacting target. Results are summarized in the table 4.2.

MeV/u

Y 1 Y 2 QMD ǫ Y1 ǫ Y2 10 C (8 ± 3) × 10 -3 (1.1 ± 0.3) × 10 -2
1.0 × 10 -2 25% 9% 11 C (10.5 ± 1.3) × 10 -2 (13.9 ± 0.2) × 10 -2 9.0 × 10 -2 15% 35% 15 O

(2.1 ± 0.3) × 10 -2 (2.8 ± 0.4) × 10 -2 2.0 × 10 -2 5% 28%

MeV/u

Y 1 Y 2 QMD ǫ Y1 ǫ Y2 10 C (1.2 ± 0.3) × 10 -2 (1.7 ± 0.2) × 10 -2 1.4 × 10 -2 16% 17% 11 C (14.7 ± 1.6) × 10 -2 (24.1 ± 0.3) × 10 -2 13.0 × 10 -2 12% 46% 15 O (3.1 ± 0.4) × 10 -2 (5.1 ± 0.5) × 10 -2 3.3 × 10 -2 6% 35% 343.46 MeV/u Y 1 Y 2 QMD ǫ Y1 ǫ Y2 10 C (1.5 ± 0.3) × 10 -2 (2.2 ± 0.3) × 10 -2 2.3 × 10 -2 53% 4% 11 C (19.9 ± 2.4) × 10 -2 (30 ± 0.4)× 10 -2 20.9 × 10 -2 5% 30% 15 O (5 ± 0.4) × 10 -2 (7.8 ± 0.5) × 10 -2 5.9 × 10 -2 18% 24%
Table 4.2: β + emitter yields in PMMA. Y 1 refers to production rate given by [I. Pshenichnov et al. (2006)] [ [START_REF] Parodi | On the feasibility of dose quantification with in-beam PET data in radiotherapy with 12 C and proton beams[END_REF]]. Y 2 refers to production calculated from depth profiles. Simulations were performed using Geant4 v.9.4.b01 using QMD package. ǫ Y 1 and ǫ Y 2 are the relative deviations between QMD and Y 1 and Y 2 respectively.

Y 1 column contains the total production rates given by [I. Pshenichnov et al. (2006)] [ [START_REF] Parodi | On the feasibility of dose quantification with in-beam PET data in radiotherapy with 12 C and proton beams[END_REF]]. Y 2 column gives the total production rates calculated from experimental depth profiles. In each case the relative diviations ǫ Y 1 and ǫ Y 2 between QMD and each yield were calculated.

As we can see, in each case, results given by simulations are in good agreement with both data sets. We can notice that 11 C and 15 O simulated yields are in better agreement with Y 1 whereas 10 C simulated yields are in better agreement with Y 2 . Moreover it seems that Y 2 are systematically larger than Y 1 . Further investigations are required to conclude regarding those two experimental data sets.

The total β + emitters yields induced by carbone ions in PMMA target for a [200 MeV/u;350 MeV/u] energy range are well reproduced by simulations given both experimental data. Such observable are intergrated production rate over the whole interacting medium. An other important physical observable is the β + emitters distribution over the target. In the section bellow comparisons with experimental depth profiles are presented.

Depth profiles : simulations versus measurements

In the previous section we saw that simulated 10 C production rates are in good agreement with Y 2 whereas 11 C and 15 O simulated yields are in better agreement with Y 1 . For this study we are focusing only on the shapes of the profiles. Both experimental data proviously introduced were not acquired following the same procedure. Y 1 experimental yields were obtained from the radioctive decay curves measured after the target irradiation. Y 2 depth profiles were measured from reconstructed images. But, depth profiles which represent the way the β + activity is distributed in the target is much more important than the production rates themselves. Since β + activity is deeply correlated to the dose distribution, the activity distribution is crucial.

Since 10 C simulated yields are in good agreement with one of the two experimental data sets introduced, whereas 11 C and 15 O production rates are in better agreement with the other, we have chosen to compare experimental depth profiles to the simulated ones after renormalizing each profile (simulated and experimental) to unity. The shapes of the simulated profiles are in quite good agreement with experimental data. On can notice on the 15 O simulated profiles there are some little discrepencies near the Bragg pic depth. This point is more obvious from the figures 4.2 and 4.3. Anyway, the shapes are well reproduced which is very promising.

Discussion

On the one hand β + simulated production rates are in good agreement with experimental data. But 10 C simulated yields are in good agreement with one of the two experimental data sets introduced, whereas 11 C and 15 O production rates are in better agreement with the other. Moreover, experimental yields calculated from depth profiles are systematically higher than those measured (Y 2 ∼ 1.5 × Y 1 ).

On the other hand simulated depth profiles are in rather good agreement with experimental measurements in terms of shape. Depth profiles shape was the only criterion used to perform comparisons because we probably missed something concerning the understanding of the experimental depth profiles.

In conclusion, Geant4 QMD package gives encouraging results in terms of β + emitters nuclei production rates. The shapes are rather well reproduced. Concerning the total production rates, further investigations are however required. If we only consider the experimental depth profiles (see section 4.1.2, 11 C and 15 O simulated profiles have to be weighted by a factor of 1.5 approximately (Y 2 /Y 1 , cf. table 4.2). In case we consider production yields (see section 4.1.1), 10 C simulated profiles have to be weighted by a factor of 0.6 approximately (Y 1 /Y 2 , cf. table 4.2). If we focus on both experimental data it is impossible to conclude.

Chapter 5

Prompt γ yield

During irradiation, a part of all γ resulting from positron annihilation are produced simultaneously with a lot of secondary particles. Among them prompt γ can be considered as physical background event regarding the PET principle detection. In order, to design a dedicated device able to control the treatment, it is very important to have an accurate estimation of secondary prompt γ yield.

Prompt γ depth profile : present status

Prompt γ imaging for the dose profile monitoring has been already investigated by Testa et al. [START_REF] Testa | Monitoring the bragg peak location of 73 mev/u carbon ions by means of prompt γ-ray measurements[END_REF]], [START_REF] Testa | Dose profile monitoring with carbon ions by means of prompt-gamma measurements[END_REF]]. First experiments performed at Ganil in Caen (France) paved the way of real time monitoring carbon therapy treatments. Prompt γ depth profiles simulations have already been performed by Le Foulher et al [START_REF] Foulher | Monte carlo simulations of prompt γ emission during carbon ion irradiation[END_REF]]. Two different experiments were done by the same group. The first experiment was held at Ganil facility in Caen to measure the prompt γ depth profile induced by a 95 MeV/u carbon ions beam in a PMMA target. The second experiment was held at GSI facility in Darmstadt (Germany) to measure the the prompt γ depth profile induced by a 310 MeV/u carbon ions beam in a water target.

Each experimental set-up was designed to detect prompt γ emitted at 90 degrees of the bean direction. Only prompt gammas that deposited more than 1 or 2 MeV in the detector, in a few ns Time of Flight (TOF) window were recorded [START_REF] Foulher | Monte carlo simulations of prompt γ emission during carbon ion irradiation[END_REF]]. Since a lot of secondary neutrons are induced by nuclear collisions, this TOF technique was introduced to discriminate these neutrons from the prompt γ. [START_REF] Foulher | Monte carlo simulations of prompt γ emission during carbon ion irradiation[END_REF]]

The physics list used was composed of binary cascade model to handle nuclear collisions, precompound model used as a smooth transition, Evaporation model to deexcite fragments and Photon emission model to finish the deexcitation chain. Geant4 v.9. In each case Exp refers to the measurements [START_REF] Foulher | Monte carlo simulations of prompt γ emission during carbon ion irradiation[END_REF]].

As illustrate in the figure 5.2 simulations overestimates largely the experimental prompt γ rate. Moreover, the shapes of the profiles are not well reproduced. In order to understand the nature of those discrepencies, we have to focus on the deexcitation methods since they handle the prompt γ emission. For our study one used the version 9.2 of the Geant4 toolkit.

Excited fragments spectrum

We have performed preliminary studies that consisted in characterize the excited fragments spectrum in terms of mass A, charge Z and excitation energy E * produced during 12 C -12 C for 75 MeV/u and 200 MeV/u incident energies. Processes which are supposed to be involded in the deexcitation procedures have been thus clearly identified. The figures 5.3 and 5.4 present Z versus A and A versus E * spectra obtained using the Geant4 v.9.2 QMD package. As shown in the figures 5.3 and 5.4, the excited fragments population is mainly composed of light fragments with Z≤7 and A≤14. Those fragments have to be deexcited by Fermi BreakUp process. The heavier fragments have to be treated according to Evaporation process.

Deexcitation chain

The G4ExcitationHandler class handles all the deexcitation models. The default class of the Geant4 toolkit has been described in the section 3. Here some results obtained with this default class are presented.

G4ExcitationHandler default class

We consider the prompt γ energy spectrum and the prompt γ emission rate during deexcitation. Photon emission is treated exclusively during deexcitation. One can distringuish 5 different modes of γ depending on the way they are produced.

The two first production modes consist of photons emitted through evaporation process. There are discrete emitted photons which are coming from quantized transitions and also continuous emitted photons from dipole transition consequently to giant dipole resonances.

The two second production modes are compound of discrete and continuous emitted photons too. But, unlike the two first types, those photons are produced at the end of the deexcitation chain, after Fermi BreakUp or Evaporation processes.

The last production mode of photons is the "forced" γ. An artificial final photon emission is performed if photon evaporation can't deexcite fragments. The table 5.1 sums up the contribution of each type of photons in the total spectrum. The main component is the forced γ component. This point is confirmed by the figures 5.5 and 5.6. This can be explained as followed. On the one hand, light clusters (A≤4) are not deexcited by Fermi BreakUp. They are only deexcited by PhotonEvaporation. But any nuclear transition is impossible since the database does not exist for such light nuclei. Indeed, nuclear transitions are tabulated from the PhotonEvaporation2.0 database. So, every light excited clusters are deexcited by "forced" emission. Moreover, their contribution is not negligeable as shown in figures 5.3, 5.4. On the other hand, excitation energy of fragments produced by Fermi BreakUp is supposed to be quantized. The Geant4 method which manages this quantization is the G4FermiBreakUp:: G4FermiFragmentPool. A look at this class shows that energy units are keV not MeV. Moreover, there is a mismatch between values from the Pho-tonEvaporation2.0 database and those which are given by the G4FermiFragmentPool method. This have some consequences because, some fragments atfer Fermi BreakUp can be only deexcited by "forced" emission since nuclear level do not coincide with the one listed in the Geant4 database.

Those two points show that are some modifications have to be done in order to reduce the simulated prompt γ emission rate, just because some γ are coming from unphysical processes. In the next section we are going to present in detail modifications introduced in the code and the consequences on the γ emission rate.

New handler algorithm

The new algorithm is divided into two different steps. One step handles the deexcitation through Fermi BreakUp and/or Evaporation as shown in figure 5.7.

At the end of this first step, one can distinguish two different types of fragments : unexcited fragments which are stored as the final products of the reaction and excited fragments which are going to be deexcite through PhotonEvaporation process. The algorithm which is used to handle the different step PhotonEvaporation process is illustrated on the figure 5.8. The table 5.2, and the figures 5.9 and 5.10 show energy spectum obtained with this new handler and the contribution of each photon depending on their production mode. One can note some important changes.

The first change is that there is no more "forced" γ. The second change we can observe is that photons coming from nuclear transitions are the main photons of the total spectrum since Fermi BreakUp is the main deexciting process. The energy range of the emitted γ is quite reduced too. Since Fermi BreakUp process is prefered given the mass spectrum of the primary fragments, nuclear transitions reaction is the natural process which follows. But, despite the discrete photon relative contribution coming from PhotonEvaporation is larger with this new algorithm, the number (normalized by the primary ion number) of this type of photon is not larger. By way of constrast the continuous PhotonEvaporation photons yields are significantly reduced.

We can give an explanation from the figure 5.11. This figure shows the excitation energy versus mass distribution of fragments just before the PhotonEvaporation process for two different incident energies (75 MeV/u and 200 MeV/u) and for the two different handler.

The default handler gives heavy high excited fragments which will be deexcited through continuous PhotonEvaporation. The new handler algorithm lead to the deexcitation of those heavy excited fragments. That's why continuous PhotonEvaporation γ yield is much more important in case of default handler. The figure 5.12 shows the total prompt γ yield for two incident energies obtained with the two handler. The new handler reduces by a factor 5 the total simulated prompt γ yield.

The figure 5.13 presents the secondary particle yields (p,n,d,t, 3 He, 4 He, 10 C, 11 C) produced by 12 C -12 C collisions simulated with Geant4 v.9.2. This new handler produces a little bit more light clusters such as p, n, d and 4 He. In fact, Fermi BreakUp model is the main deexcitation process and it is a more or less explosive like process. This leads to split fragments via the emission of light clusters. However, yields are not very different which means this new algorithm impacts significantly on photon yield only.

To complete this study one can consider the excitation state of fragments at the end of the deexcitation chain. The figure 5.14 shows the excitation state of the final simulated fragments ((a) for 75 MeV/u incident energy and (b) for 200 MeV/u incident energy). Some fragments are not completly deexcited. Those fragments account for ∼ 15% of the total secondary particles produced by nuclear collisions. QMD and statistical deexcitation models do not take into account nuclear spin orbit neither pairing couplings between nucleons.

However such kind of interaction are very important since they govern the nuclear

∆ ≈ 12A -1/2 M eV (5.1)
For light clusters, pairing gaps account for few MeV which is comparable to the residual excitation energy. On the other hand, ground state calculated are not the true ground state since they are calculated in order to fit nuclear saturation point (ρ = ρ 0 and E = -15.75 M eV ). In conclusion given all these arguments this is not a real problem to have some fragments with a small residual excitation energy.

Until now were presented some modifications introduced to reduce prompt γ production. Comparisons with experimental data are required to validate those modifications. In the following section, simulations were performed with the last released Geant4 version (9.4). This report were initiated before this version was available. Some modifications were introduced (compared to the previous versions) in the G4ExcitationHandler class and G4PhotonEvaporation classes. However, the new G4ExcitationHandler::BreakIt Up() method still remain different from our own modified method. In fact, despite Fermi BreakUp model is extended to the light excited clusters such as 4 He or t, the Evaporation model is called in case Fermi BreakUp can not deexcite fragments. In those precise cases, we have noticed that the photon evaporation channel is always chosen. Moreover, the method used to perform the photon evaporation process produces only one photon. As a consequence, there is no photon emitted when fragments are processed by the last photon emission loop. That's why we choose to implement our own G4ExcitationHandler::BreakItUp() method in the last released Geant4 version.

Simulations vs measurements

Preliminaries

Before discussing comparison between simulations ang Ganil experimental results explain the nuclear level selection procedure as implemented in the G4PhotonEvaporation classes. This selection procedure is performed by the NearestLevel method from the GEANT4 class G4NuclearLevel. This method tackes two arugments. The first one is the excitation energy of the fragment whereas the second argument is an artificial tolerance.

Basically, a radiative transition is allowed when the difference between the excitation energy and a level found in the nuclear database is bellow the tolerance.

It is possible to tune this artificial tolerance. In case there is a large tolerance, the closer tabulated level is chosen even if the difference between the excitation energy and the tabulated one is huge, e.g. 2 MeV. On the other hand, if the tolerance is very low, only fragments which the excitation is striclty quantized will be deexcited. By default, this tolerance is very large (10 TeV), which means there is no constraint on the transition selection.

Ganil experiment

The figure 5.15 shows comparison between simulations and the Ganil experiment results. Simulations are in rather good agreement with experimental data for a 200 keV tolerance value. Shape is well reproduced despite a small discrepancie near the Bragg peak depth. But this deviation is no larger than few dozen percent which is in much more better agreement than previous results.

By tuning the nuclear level selection procedure, we are able to to reproduce Ganil experimental data. The prompt γ rays at 95 MeV/u are well reproduced. However, we do not know if this energy level selection procedure tolerance value will give good agreements at higher energies and for different targets and projectiles.

GSI experiment

Consider experimental results obtained at GSI in Darmstadt for carbon ions at 310 MeV/u interacting with a water phantom. The figure 5.17 shows measured and simulated depth profile. Simulations were performed with the QMD package, and the modified handler. Two tolerance values were considered : a 200 keV tolerance (red plot) and the default one (black plot). .17: Measured versus simulated prompt γ depth profile : 310 MeV/u 12 6 C 6+ on water target at the GSI At this energy range, we can observe that the prompt γ detection yield is less sensitive to the level tolerance value than at low energy (∼ 100 MeV/u, cf. previous section). Simulations with a 200 keV tolerance and with a 10 TeV tolerance are very close. This point is very interesting. Since results are very similar whatever the tolerance, is it reasonnable to choose a 200 keV tolerance value whatever the incident energy of the primary ions. In conclusion, the new handler gives quite good results. The prompt γ yield is significantly reduced. However, we have to mitigate one point. We have introduced a parameter in order to controle the prompt γ production rate. This is quite dubious since this tolerance create band gaps. In fact, any fragment with an excitation energy out of this energy band will not be deexcited.

On the other hand, this tolerance is a way to control the photon emission. If any constraint is added to the selection procedure, any fragment with a little excitation energy will be deexcited. This leads to produce too much prompt γ. This tolerance is a kind of trick but not a real physical solution. The correct solution should be to modify the nuclear collision models to avoid the formation of such fragments. This observation requires further measurements and simulations.

Conclusion

The section 1 has presented different type of nuclear collision models. More precisly, QMD model was developped and each ingredient of the model were introduced. QMD was prefered to other nuclear models since it gives for now the best results for hadrontherapy applications. Moreover, it is a theoritically able to simulate low energy nuclear collisions because it takes mutual interactions through Skyrme potential into account.

The section 2 was dealing with the different models used to describe the fragment deexcitation. Fermi Break model is suitable to deexcite light fragments (A<16). Evaporation model seems more appropriate to the treatment of heavier fragments. The different channels can be divided into three components : particles evaporation (p, n, d, t, 3 He, 4 He and heavier cluster such as 29 M g), photon evaporation essentially governed by the giant dipole resonances, and the competitive fission. In order to complete the deexcitation, an additionnal photon emission based on nuclear radiative transition is also suitable. Fermi BreakUp and nuclear transition seem dedicated to the deexcitation of fragments induced by nuclear collisions between biological element (O,C,H,N).

The section 3 was dedicated to the description of the implementation of the previously presented models in the Geant4 monte carlo toolkit.

The section 4 has reported β + emitter production yield simulation and comparison with experimental data. QMD package was used for the simulations. Simulations provide good results in comparisons with experimental data. However, we can not compare simulations with data acquired at low energy (< 100 MeV/u) because of the lack of measurements. This point will be investigated later since an experiment is planned in the current year to provide such data for 95 MeV/u carbon ions beam at Ganil in Caen.

The section 5 we have investigated the present status on prompt γ ray. We have investigated on the the prompt γ yield reduction and provide some comparisons with measured data. Simulations give good result at 95 MeV/u and and at 310 MeV/u but further measurements especially at 310 MeV/u are required to conclude properly.

Results presented in the sections 4 and 5 are relatively consistent regarding the experimental data. Annihilation γ and prompt γ are well reproduced in terms of spatial distribution and production yields. The proper estimation of those two types of particles are very important for the conception of a PET based device to the control of hadrontherapy treatment.

This document was initiated in order to summarize all the work we have done on Monte Carlo simulations dedicated to hadrontherapy applications. We are concerned by R&D on TOF PET detectors. The proper estimation of secondariy particles produced counting rates is crucial. Only QMD model was investigated. Despite the other nuclear collision models available in Geant4, such as binary cascade model, QMD model is the more precise for particle therapy energy range. However, it is clear that those models were not concieved to study carbon ion fragmentation at about a hundred MeV/u. Many groups in Europe are working on Monte Carlo simulations for particle therapy purposes and they need accurate simulation tools. Geant4 gives good results, but many improvements are required. The modifications we have introduced are not sufficient. The use of a tolerance value to control prompt γ yield is just a palliative. On the one hand, the only available data are prompt γ detection yield. It is hard to provide a nuclear model validation from those data sets. On the other hand, prompt γ are produced during fragment deexcitation. But we can not check if the primary fragment distributions are correctly predicted by the simulation code. Obviously, we have to work on the nuclear collision models.

Finally, we are able to correctly reproduce some existing experimental data. This point is paving the way of accurate simulations to the design of an in beam TOF PET dedicated to the hadrontherapy treatment.

Figure 1

 1 Figure1illustrates clearly that even a small variation in the density induces significant modifications of the range.

  Figure1.2: The nuclear matter equation of state : the density dependence of the energy per particle in nuclear matter at temperature T=0 is displayed for 4 different sets of parameters [J.Aichelin (1993)]

Figure 2 . 1 :

 21 Figure 2.1: Diagram of the nuclear radiative transition.

Figure 3

 3 Figure3.1: QMD package description. The different steps of the QMD package are shown from top to bottom. the choice of the projectile and target configurations, then the initial dynamic phase that finish with a set of excited fragments. the "handler" that deals with all the deexcitation channels and produces at the end a set of stable particles and fragments.

  and Z<9 G4FermiBreakUp::BreakItUp(A,Z,E * ) E * >3 MeV/u G4StatMF::BreakItUp(A,Z,E * ) A>17 or Z>9 G4Evaporation::BreakItUp(A,Z,E * ) 2 nd step A set of intermediate excited fragments {A ′ ,Z ′ ,E * ′ }

Figure 3 . 2 :

 32 Figure 3.2: Deexcitation processes implemented in the QMD package of the Geant4 code

Figure 4

 4 Figure 4.1: β + emitter experimental profiles versus depth for 215 MeV/u, 260 MeV/u and 337 MeV/u carbon ions.

Figure 4 Figure 4 Figure 4

 444 Figure 4.2: β + emitter profile versus depth for 215 MeV/u carbon ions. data refers to normalized experimental depth profiles. QMD represents normalized simulated depth profiles using QMD package(v.9.4.b01)

Figure 5 . 1 :

 51 Figure 5.1: Experimental set-up dedicated to the prompt γ depth profiles measurement. (a) Ganil set-up : 95 MeV/u carbon ions on PMMA target. (b) GSI set-up : 310 MeV/u carbon ions on water target.[START_REF] Foulher | Monte carlo simulations of prompt γ emission during carbon ion irradiation[END_REF]

Figure 5 . 2 :

 52 Figure 5.2: Measured versus simulated prompt γ depth profile. (a) Ganil experiment. (b) GSI experiment.In each case Exp refers to the measurements[START_REF] Foulher | Monte carlo simulations of prompt γ emission during carbon ion irradiation[END_REF]].

Figure 5 . 3 :Figure 5

 535 Figure 5.3: Z versus A spectrum of the QMD Geant4 v.9.2 simulated fragments population produced during 12 C -12 C collisions at 75 MeV/u (a) and 200 MeV/u (b)

Figure 5 . 5 :

 55 Figure 5.5: Simulation of prompt γ energy spectrum in 12 C -12 C collisions at 75 MeV/u

Figure 5 . 6 :

 56 Figure 5.6: Simulation of prompt γ energy spectrum in 12 C -12 C collisions at 200 MeV/u

Figure 5 . 9 :

 59 Figure 5.9: Geant4 v.9.2 simulated prompt γ energy spectrum for 12 C -12 C collisions at 200 MeV/u. The new handler algorithm is used to handle the deexcitation chain

Figure 5 .Figure 5 . 11 :Figure 5 .

 55115 Figure 5.10: Geant4 v.9.2 simulated prompt γ energy spectrum for 12 C -12 C collisions at 200 MeV/u. The new handler algorithm is used to handle the deexcitation chain

Figure 5 . 13 :Figure 5 . 14 :

 513514 Figure 5.13: Secondary particle Geant4 v.9.2 simulated yield for 12 C -12 C collisions at 75 MeV/u (a) and at 200 MeV/u (b) with the two different handler

Figure 5 . 15 :Figure 5 . 16 :

 515516 Figure 5.15: Measured versus Geant4 v.9.4.b01 + new handler simulated prompt γ depth profile: 95 MeV/u 12 6 C 6+ on PMMA target at the GANILThe blue plot refers to the measured prompt γ yield as a function of depth in the PMMA target whereas the red plot refers to the simulated depth profile using QMD package. The

  Figure 5.17: Measured versus simulated prompt γ depth profile : 310 MeV/u 12 6 C 6+ on water target at the GSI

1

  Geant4 v 9.2 simulation : dose profile from 12 6 C 6+ at 200 MeV/u in water equivalent material for different densities. Even a small variation in the density induces a significant modification of the range. . . . . . . . . . . . . . . . 1.1 Parameters of the potential for different type of EOS. . . . . . . . . . . . . . . Deexcitation processes implemented in the QMD package of the Geant4 code 4.1 β + emitter experimental profiles versus depth for 215 MeV/u, 260 MeV/u and 337 MeV/u carbon ions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 β + emitter profile versus depth for 215 MeV/u carbon ions. data refers to normalized experimental depth profiles. QMD represents normalized simulated depth profiles using QMD package (v.9.4.b01) . . . . . . . . . . . . . . . . . . 4.3 β + emitter profile versus depth for 260 MeV/u carbon ions. data refers to normalized experimental depth profiles. QMD represents normalized simulated depth profiles using QMD package (v.9.4.b01) . . . . . . . . . . . . . . . . . . 4.4 β + emitter profile versus depth for 337 MeV/u carbon ions. data refers to normalized experimental depth profiles. QMD represents normalized simulated depth profiles using QMD package (v.9.4.b01) . . . . . . . . . . . . . . . . . .

1.2

The nuclear matter equation of state : the density dependence of the energy per particle in nuclear matter at temperature T=0 is displayed for 4 different sets of parameters [J.

Aichelin (1993)

] . . . . . . . . . . . . . . . . . . . . . . . 2.1 Diagram of the nuclear radiative transition. . . . . . . . . . . . . . . . . . . . . 3.1 QMD package description. The different steps of the QMD package are shown from top to bottom. the choice of the projectile and target configurations, then the initial dynamic phase that finish with a set of excited fragments. the "handler" that deals with all the deexcitation channels and produces at the end a set of stable particles and fragments. . . . . . . . . . . . . . . . . . . . 3.2

Table 3 .

 3 2: α, β, γ, ρ 0 , L and E sym parameters available in Geant4The table 3.3 gives an overview of the parameters of the effective potential.

	c 0	c 3	c s	c l	c Coul
	-0.00518291 0.00121007 0.000590529 0.000719883 0.001439767

Table 3 . 3 :

 33 Total effective potential parameters available in Geant4

Table 3 . 4 :

 34 Ground state parameters available in Geant4

75 MeV/u

  Table 5.1: Prompt γ spectrum description . Geant4 v.9.2 simulated 12 C -12 C collisions

	Continuous	Discret	Continuous	Discret	Forced
	Evaporation Evaporation PhotonEvaporation PhotonEvaporation Emission
	4.10 -4	2.10 -3	7.10 -2	8.10 -2	12.10 -2
			200 MeV/u
	Continuous	Discret	Continuous	Discret	Forced
	Evaporation Evaporation PhotonEvaporation PhotonEvaporation Emission
	3.10 -4	1.10 -3	9.10 -2	33.10 -2	80.10 -2

  First step of the deexcitation. Only Fermi BreakUp and Evaporation are considered. Second step of the deexcitation chain. Only PhotonEvaporation is considered.

				E * = 0 ?
			γ	A,Z,E * prompt energy spectrum	final products
				γ	prompt
	list size>1 ? dN/dE (1/ion/200 keV)	-5 -4 -2 -3 10 10 10 10	E * > 0 ? (A < 17 && Z < 9) Fermi BreakUp else Evaporation Intermediate fragment list list size=1 ? photonEvap list prompt γ forced Discrete PhotonEvaporation Continuum PhotonEvaporation prompt γ prompt γ prompt γ Discrete Evaporation prompt γ Continuum Evaporation
		-6 10
	E * = 0 ? PhotonEvaporation E * > 0 ? Intermediate fragment list list size=1 ? final products final products Figure 5.8: Energy (MeV) list size>1 ? Figure 5.7: photonEvap list 0 2 4 6 8 10 12 14 16	18	20

Table 5 .

 5 2: Total spectrum obtained with the new handler

			75 MeV/u
	Continuous	Discret	Continuous	Discret	Forced
	Evaporation Evaporation PhotonEvaporation PhotonEvaporation Emission
	1.10 -3	5.10 -3	1.10 -3	5.10 -2	0
			200 MeV/u
	Continuous	Discret	Continuous	Discret	Forced
	Evaporation Evaporation PhotonEvaporation PhotonEvaporation Emission
	1.10 -3	4.10 -3	1.10 -2	33.10 -2	0