|
|
Flattening proton dose distributions is an important method for large-field radiotherapy. (see: A.M. Koehler, R.J. Schneider and J.M. Sisterson: Flattening of proton dose distributions for large-field radiotherapy. Med Phys 4: p.297-301 (1977). A simple example of how to do this with the help of Turtle's Monte-Carlo capabilities has been presented here. When trying to do more complicated examples like those of Koehler et al, then ring-shaped scatterers have to be used. An example of modelling a ring-shaped scatterer with Turtle is given here.
Recent simulations of PSI's SINQ beam line have shown some discrepancies between the results computed with the Windows version (TurtleNT) and the one for Linux (Turtlin). A comparison of the corresponding source codes showed some missing code fragments in the Linux (and Mac OS-X11 = 100% Linux compatible) version. These differences have now been eliminated.
In case an arbitrary matrix which should influence the momentum spread (coordinate 6) was inserted into the input deck, then the momentum spread was not updated at all. For this case Turtle histogram output did not produce data corresponding to the program Transport. This bug (it was discovered by Pierre Roussel of ipno.in2p3.fr) has now been corrected.
The focal-length expression f in the used formula for computing the transfer matrix elements for the electrostatic accelerator was changed from the non-relativistic to the relativistic expression {f = -2*p0*v0/(E1-E0), was -4*T0/(E1-E0) , see Compendium of Transport Enhancements}. A longtime nuisance for some users was the fact, that imported Turtle input files from the Windows environment did not work under Linux/MacOSX-X11 (due to a non-feature of the GNU g77 FORTRAN compiler) without removing the CR-character (ASCII-value 13) before the LF-character at the end of each line. This problem has now been resolved by having the (Fortran-) code replacing each CR-character by a blank-character during the read-in process.
When trying to run the 'Turtle Histogram Option' (gratur) with some of the newest Linux distributions (e.g. SuSE 9.2) you got the following error message: ./gratur: relocation error: ./gratur: symbol errno, version GLIBC_2.0 not defined in file libc.so.6 with link time referenceThis was caused by the fact, that all sources are still compiled and linked with the gcc-version of RedHat 6.0 (year 2000). It was assumed that this would guarantee the upward-compatibility with all newer Linux distributions. By substituting the library routine 'errno' through an own dummy routine this complication could be removed. It has been observed, that the statistical data (sums and RMS-widths) of histograms were getting inaccurate or even wrong if more than typically 100000 rays were run through a beam line. The main reason for this was the usage of only single precision real arrays for the accumulators. By computing the critical statistical data in double precision for each histogram at the end of a TURTLE run, this bug has now been corrected without affecting the overall performance of the program. (Why running more than 100000 particles through a system ? This has to do with poor statistics when using targets, scattering foils, slits, collimators or degraders.)
This type code only worked for collimators (a series of slits, which is the default for the type code 6 elements). If the slit functionality is inverted with the instruction line '16. 160. ;' to target-disc instead of slit, then the type code 29 was not applicable anymore. This has now been corrected. The new feature may be called 'Jump over Target'. Of course the 'Jump over Collimator' still works. The effect of the instruction line '16. 160. ;' may be reverted back to default with the instruction line '16. 161. ;'. Example for 'Jump over Target-disc': 16. 160. /invs/ ; (enable inverse slit = target disc) 29. 0.004125 2.5 2.5 2.5 2.5 5. 1100. /Jump/ ; 9. 1100. ; 6. 1. 2.5 3. 2.5 /Disc/ ; 1. 3.352 6.33 50.48 .506 119.6 0.0155 0.0377 0.0 /W/ ; (tungsten) 3. 0.00000375 ; 9. 0. ;Addendum to type code 29 (Jump over Slit). (added in December 2008) Because no checks what so ever are made over which type codes the particle has to jump, the "target" may also be replaced by a normal slit. Therefore if the slit is preceded by a "jump outside over" sequence and right afterwards the normal slit-scattering is enabled again, then we have nothing else than a ring-scatterer (if the whole thing is embedded in a repetition sequence). The following example demonstrates this with a tungsten ring of 1 cm thickness and an inner diameter of 10 mm and an outer diameter of 20 mm. The requested 2-dim histogram (r over p) at the exit of the ring proves the desired interaction of the protons with the "tungsten geometry". /Test Ring scatterer/ 100000 15. 1.0 /MM/ 0.1 ; 15. 11.0 /MeVc/ 0.001 ; 15. 12.0 /MeV/ 0.001 ; 1.0 1.0 10. 1.0 10. 0.0 0.1 369.1 /BEAM/ ; (Ekin = 70 MeV beam) 16. 3.0 1836.15 /m0/ ; 16. 165. /MULT/ ; (enable multiple scattering) 16. 190. 0. 100. /Meta/ ; 3. 2.0 ; 9. 100. ; (100 * 0.0001 m = 0.01 m = 1 cm) 16. 160. /invs/ ; (enable outside jumping) 29. 0.0001 20. 20. 20. 20. 4. 1. /Jump/ ; (jump over next 4 cards when outside 20 mm) 16. 161. /norm/ ; (enable normal slit scattering) 6. 1. 10. 3. 10. /Scat/ ; (slit of inner diameter 10 mm) 1. 3.123 9.74 5.05 .51 24.7 0.082 0.207 0.0 /W/ ; (tungsten scattering material) 3. 0.0001 ; 9. 0. ; 51. 11. 200. 400. 10. ; (x-coord. = 11 = "p") 52. 21. 0. 40. 2. ; (y-coord. = 21 = "r") SENTINEL; SENTINEL;
If the aperture constraint instruction 13. 10. + id was present, then no particles were passing through an octupole because the aperture data were not properly stored. But this bug did not always show up, because it was masked in case a 16. 100. a. b. ; instruction (Quadrupole apertures) was present. This bug has now been corrected.
Because of the relatively long fringe field of solenoids, the tracking method with type code 28 may cause overlapping problems, because with 2 nearby placed solenoids the first type code 28 sequence may not yet be finished when the second one should start. This problem has been solved with some extensions for the type code 28 element. For the description of these new features, a sample input-code and the corresponding output-histograms click here. This port is 99.9% compatible with the Linux-version. So all features, which are running under Linux are also running on the Mac. It has been tested with Mac OS X Version 10.3.2 and X11 1.0 - XFree86 4.3.0. Besides X11 (from Apple) tcltk and xemacs (e.g. from http:// fink.sourceforge.net/index.php) are needed for fully functioning.
The type code 19 for tracking particles through solenoids with Turtle has proven to be not very satisfactory, because of missing codes for higher order effects on the beam. Therefore an additional new type code has been implemented, which allows through segmentation to include the parameterization of the field map of a solenoid. For the description of this type code, a sample input-code and the corresponding output-histograms click here. When tracking particles through a bending magnet (type code 4) with a zero B-field (e.g. a switch-yard magnet), then all particles were disappearing inside this magnet. This bug has been eliminated; for B=0 the bend is treated as a drift space.
In order to assure that the compiled and linked Turtle program modules run with all Linux versions of the last few years, they are compiled and linked with 'gcc' under RedHat 6.0. But when running these modules with newer Linux distributions, an exception may show up because a needed shared object (libncurses.so.4) is usually missing and has painfully to be searched for in older distributions and put into /usr/lib in order to get the images 'turtle' and 'muscat' running. In order to avoid this, the needed code of 'libncurses' is now statically linked into the 2 mentioned program modules.
The multiple scattering code in Turtle was delivering results between inaccurate and wrong, if two or more considerably degrading devices of very different materials were put in series. This was due to the usage of the parameter PBEAM (original central momentum) instead of the parameter P0 (design momentum) in some parts of the code. This has now been corrected, but requires from the Turtle user, that he has to keep track by himself of the changing design momentum P0 along the beam line. This is also required, because bending magnets and the expression dp/p depend on this parameter P0 as well. Furthermore, the first energy on the dE/dx-card (type-code 30, see Compendium of Turtle Enhancements) has now to correspond (within +/- 2%) to the actual P0 (and not to the initial central momentum PBEAM as before). How this has to be done now is shown in the following Turtle input sample code: /Degrading of 250 MeV protons/ 100000 15. 11.0 /MeVc/ .001 ; 15. 12.0 /MeV/ .001 ; 15. 6.0 /PM/ .1 ; 15. 1.0 /MM/ .1 ; 13. 10. ; 16. 190. 0. 100. /FILE/ ; 1.0 1. 2. 1. 2. 0. 0.1 729.1 /BEAM/ ; (E=250 MeV, PBEAM=P0=729.1 MeV/c) (initial design momentum P0 set internally to the original central momentum PBEAM) 16. 3.0 1836.7 /MASS/ ; 16. 165.0 /MULT/ ; 16. 198. 100.0 ; . . (Beam line data) . . (Graphite degrader: 250 MeV --> 70 MeV) 16. 160. /INVS/ ; 30. 250. 3.55 200. 4.05 100. 6.60 10.0 41.5 1.0 240.0 0.0 /dEdx/; (E0=250 MeV) 9. 895. ; (loop over 895 degrader segments) 6. 1.0 50.0 3.0 50.0 /DEGR/ ; 1. 0.879 9.30 295.5 .630 144.2 0.0315 0.1261 0.0 /C/ ; 3. 0.0002 ; (895 * 0.0002 m = 0.179 m) 9. 0. ; (end of loop) 16. 161. /NORM/ ; 16. 11. 369.1 /P0_2/ ; (new design momentum for Ekin=70 MeV) 50. 6. -50. 50. 2.5 ; (dp/p histogram) 50. 11. 360. 380. .5 ; (Momentum p histogram --> P0_2 may be extracted) . . (Beam line data) . . (Tantalum scattering foil of 0.7 mm thicknes) 16. 160. /INVS/ ; 30. 70. 4.748 50. 6.036 20. 11.61 10. 18.65 1. 67.79 0.0 /dEdx/ ; (E0=70 MeV) 9. 76. ; 6. 1.0 50.0 3.0 50.0 /Foil/ ; 1. 5.058 7.29 58.22 0.506 114.8 0.0330 0.0818 0.0 /Ta/ ; 3. 0.00001 ; (76 * 0.01 mm = 0.76 mm) 9. 0. ; 16. 161. /NORM/ ; 16. 11. 342. /P0_3/ ; (new design momentum) 50. 6. -50. 50. 2.5 ; (dp/p histogram) 50. 11. 335. 355. .5 ; (Momentum p histogram --> P0_3 may be extracted) . . (Beam line data) . . SENTINEL SENTINEL When displaying Lego- and Grid-plots, the phi and theta angle may be varied simultaneously while dragging the corners of the plot's base rectangle with the mouse pointer. The printouts of the Graphic Histogram option have now the same aspect ratio for both the landscape and the portrait format (scaled down by a factor of 1.4).
The to small font size and the to small left margin of text printout have been corrected. The empty trailing page has been removed. Note: Printing of the content of text windows may be started from its system menu. The Info text of the GUI-editor (help file) may now also be printed via its system menu.
In order to generate a smooth particle distribution in phase space at z = 0 with lesser particles, a different module than the standard random generator may be enabled with the 16. 31. ; code for being used in the subroutine emit. This alternative routine produces so-called quasi-random numbers, i.e. numbers that are generated by number-theoretical methods, and which tend to distribute themselves (in one or several dimensions) in a self-avoiding manner. These numbers combine some virtues of both randomly and regularly distributed numbers. In order to avoid correlations with different parameters, for all other statistical processes used in Turtle (multiple scattering, RMS-add to beam, decay and magnet-jitter) only the standard random generator is always being used. If the "seed" (code 16. 30. iseed; iseed > 0) is frozen in, then the 6 used Halton-series start with each run at the same numbers. Otherwise these starting numbers vary statistically from run to run (default). For a demonstration of this new feature use the file test.dat as input (FOR001.DAT).
Aperture constraints observation may now be switched on (13. 10.+id ;) or switched off (13. 11.+id ;) along the beam line. The default is still no observation of aperture constraints. (Before this modification there was only the choice between on or off for the whole beam line.) This new feature may be helpful for 'what-happens-if-decisions'.
Due to an error in the program logic, the exit fringe field was not computed, if the quadrupole was divided into 2 or more sections. Example: 2. 1000. ; 5. L/2. B. a. ; 5. L/2. B. a. ; 2. 1000. ; (had no effect) 2. 1000. ; 5. L. B. a. ; 2. 1000. ; (was computed)This error has now been eliminated !!
While dragging with depressed left mouse and control buttons at one of the corners at the base of a Lego- or grid-plot (started e.g. with taskbar item 'Show Histograms'), the 3D-histogram is rotated. In order to avoid the initially existing flickering while dragging, two measures against it were undertaken:
The new type code '16. 191. unit count ;' writes in third order Turtle runs tables of the 3rd order contributions for the first 'count' rays along the beam line onto a file with the name FOR0xy.DAT (with unit=xy). Each term Uijkl is already multiplied with the corresponding ray-coordinates j, k and l (in order to see immediately the magnitude of each component). On top and bottom of each table the {x,x',y,y'}-ray vector before and after the addition of the 3rd order components are displayed. So far only *QUAD* and *QUAD FRINGE* elements are displayed. Example: 16. 191. 33. 100. ; produces an output to FOR033.DAT with the first 100 rays. Sample output (only first few lines):
Turtle third order info:
========================
*QUAD FRINGE* with ray 1 at z= 0.450 m
x = -7.789 cm, x' =-183.945 mr, y = 5.929 cm, y' = 144.498 mr
U1111= -0.035 , U2111= -0.402 , U3333= -0.014 , U4333= 0.162
U1133= -0.057 , U2133= -0.657 , U3311= -0.077 , U4311= 0.890
, U2112= 2.402 , U4334= 1.063
, U2134= -2.879 , U4312= -3.547
, U2233= 1.309 , U4411= 1.950
x = -7.881 cm, x' =-184.172 mr, y = 5.838 cm, y' = 145.016 mr
*QUAD* Q1 with ray 1 at z= 0.850 m
x = -8.932 cm, x' = 137.683 mr, y = 17.460 cm, y' = 500.256 mr
U1111= 0.046 , U2111= 5.574 , U3333= 0.114 , U4333= 1.708
U1112= -0.028 , U2112= -1.090 , U3334= 0.440 , U4334= 13.071
U1122= 0.126 , U2122= 0.828 , U3344= 0.405 , U4344= 15.666
U1133= 0.104 , U2133= 15.496 , U3113= -0.051 , U4113= -0.703
U1134= 0.146 , U2134= 23.360 , U3123= -0.187 , U4123= -17.139
U1144= 0.120 , U2144= 10.745 , U3223= 0.116 , U4223= 1.318
U1222= 0.063 , U2222= 1.904 , U3444= 0.097 , U4444= 4.458
U1233= 0.136 , U2233= 3.240 , U3114= 0.189 , U4114= 7.119
U1234= 0.088 , U2234= 5.079 , U3124= -0.076 , U4124= -10.456
U1244= 0.053 , U2244= 3.397 , U3224= 0.026 , U4224= -2.524
x = -8.077 cm, x' = 206.216 mr, y = 18.533 cm, y' = 512.773 mr
.
.
.
The matrix elements R21 and R43 were too large by a factor of 2. The 3rd order contributions depending on the fringing field integral I4 were also wrong by a factor depending on the particle momentum. These errors have now been corrected. The degree of higher order calculation along the beam line may be changed. This feature is useful for testing the higher order effects of individual elements. The codes are: 17. 1. /1st/ ; (first order, default) 17. 2. /2nd/ ; (second order, 17. /SECO/ ; is also for second order) 17. 3. /3rd/ ; (third order)Example: /Test, second order run with 2 quads in third order/ 100000 17. /SECO/ ; . . 17. 3. /3rd/ ; 5. 0.5 4.11 15. /Q1/ ; 3. 0.3 ; 5. 0.5 -3.22 15. /Q2/ ; 17. 2. /2nd/ ; . . SENTINEL SENTINEL For studiing the flattening of gaussian distributions with the help of an annulus and a following second scatterer it is usefull to have some code producing a cut-off of the particles outside a radius r1 and inside a radius r2 (>r1) [ring]. The following code is doing this: 6. 1. r1. 1. r2. ; (both coordinates have to be 1 and r2 > r1) The codes 16. 160. ; for inversion and 16. 177. dr1. dr2. ; for incrementing/decrementing r1 and r2 in code 9. n. ; - loops as well as the add to beam feature (multiple scattering) are also working in connection with this ring-coding.
In order to make the installation of Turtle for Windows easier and more transparent, the self extracting file (turwin32.exe) has been compressed with PKZIP's PKSFX option, which allows it to specify the unpacking directory and to start setup.exe in a dialog window. During the design phase of a new beam line there is always the question of how stable the currents generated by the power supplies for the quads and bends have to be. For this reason the pole tip field of the quads and the design momentum p0 inside the bending magnets may be varied randomly by a small amount. (The reason why p0 has to be varied inside a bending magnet is because changing the B-field of a bending magnet does only affect the focusing and not the direction of the outgoing beam, which is the main concern caused by unstable bending magnet supplies.) The randomly added values may have a cosine or a linear distribution (see below). If several elements in a beam line are varied together, then the randomly added values for each element are uncorrelated. (Exception: If quads or bends are divided into several consecutive fragments, then all fragments which belong together {quads having the same field value or bends being between the same pole face rotation elements} are changed for each ray by the same amount.) For a bending magnet this 'randomly add to p0' feature is only working if it has a preceding entrance and a trailing exit pole face rotation element. The fringing field elements of a quad are not taken into account, because they represent only small contributions to the effect of the whole quad. Newly created type codes: 16. 35. DQAmp. ; (random amplitude in units of the pole tip field [def=kG]) 16. 36. dpp0. ; (random amplitude in units of dp/p [def=%]) (if DQAmp or dpp0 are > 0.0, then a cosine, otherwise (<0.0) a linear distribution will be selected.) Examples: 16. 35. 0.003 ; (random amplitude is 3 Gauss) 5. 0.51 -3.1 25. /Q1/ ; 3. 0.20 ; 5. 0.51 2.5 25. /Q2/ ; 16. 35. 0.0 ; (disables random addition for following quads); 16. 36. 0.01 ; (random amplitude is 0.01 %, equal to 1 Gauss in B) 2. 0. ; 4. 1.0 10.0 0.0 /B1/ ; 2. 0. ; 16. 36. 0.0 ; (disables random addition to p0 for following bends);
You have now the possibility to get connected to the PSI Turtle home page via the Web by either clicking at the "Help" menu item or the WWW button of the 3rd "about" box. New feature: In case the "short" button in the "View in List-box" section is pressed, then the type code 13, 15 and 16 element's functions are described in a few words depending on their meaning they have in Turtle. The implementation of the electrostatic deflecting device has been done as close as possible to the corresponding magnetic device: 16. 94. ; (enable electrostatic deflection) 4. L. alpha. pq. /EPRI/ ; 16. 93. ; (back to magnetic bending)The meaning of the parameters is the following:
L = central length of prism in meters. [deflection radius re0 = (180.0 / pi) * L / alpha]
alpha = deflection angle in degrees.
pq = 0.0 for cylindrical shape.
1.0 for spherical shape.
all other values: toroidal shape with pq = re0 / Re0. [Re0 = vertical curvature radius]
- Only first order matrix elements are included.- No type code 2. for fringing field is available. - Use 16. 4. d/2. ; for horizontal aperture extent of the device. - Use 16. 5. h/2. ; for vertical aperture extent of the device. The implementation of the fringing field for electrostatic quadrupoles up to 3rd order has been done the same way as for magnetic quadrupoles. (see Compendium of Turtle Enhancements): -17. 3. /3rd/ ; (may be enabled) 16. 25. I1. ; (fringe field integral I1) 16. 26. I2. ; (fringe field integral I2) 16. 27. I3. ; (fringe field integral I3) 16. 28. I4. ; (fringe field integral I4) 16. 96. ; (switch to electrostatic QUADs) 2. 1000. ; (compute entrance fringe field) 5.0 Leff. Ep. a. /EQ/ ; (Ep is the electric field [kV/cm] at the pole tip) 2. 1000. ; (compute exit fringe field) 16. 95. ; (switch back to magnetic QUADs)Note: If third order computation is enabled, then also the inner region of the electrostatic quadrupole (i.e. the whole device) is computed up to 3rd order. All input parameters are now displayed in the input dialog box. When you exit the program the last input values are stored on a file and taken as default for the next session. The positioning of the different windows has been improved (works much better under Windows NT4/2000 than under Windows 9x). For the 'graphical histogram display option' viewing some statistical info (like mean value and rms-width) together with the histogram has been added. This feature may be toggled on/off (default=on).
For the 2 'graphical histogram display options' viewing some statistical info (like mean value and rms-width) together with the histogram has been added. This feature may be switched on or off (default=on). The auxiliary program Gratur32.exe has been split into 2 threads: The first thread is a C++ class build around a dialog box to handle the input of the many parameters. The second thread is the initial Fortran program being able to display the histograms stored on the metafile FOR100.DAT with Cernlib's Hplot/higz. Both threads are communicating with each other via a small data file. The big advantage of the new interface is to bypass the complicated and time consuming sequential input via standard Fortran-I/O. Traditionally Lahey Fortran with its TNT DOS extender was used for the Windows 9x (95,98,Me) operating systems. But with the apparent limitation to only 8 char directory names it has been dropped. Now for all Fortran programs the Digital (Compac) Fortran compiler is used for all Windows operating systems. (Today there is no need anymore for using a DOS extender.) The auxiliary program Muscat has been split into 2 threads: The first thread is a C++ class build around dialog boxes to handle the input of parameters. The second thread is the initial Fortran program being able to compute the parameters for multiple scattering. Both threads are communicating with each other via a small data file. For each input parameter some help is now available and the input procedure is more comfortable. When you install Turtle for Windows, you are asked into which directory it should be installed. Then - besides some DLL's, an editor and a font - all files (executables and data files) are put there. By default all newly created data files are also stored there. If you want to change this behavior - that means you want the data file access redirected to a different location - then you have to change the 1st command line parameter with which Ttmenu32.exe (yellow/black/red turtle icon) is called (default= C:\turtle\). You may access this parameter in the properties/shortcut/target entry. By default you would find there: C:\turtle\ttmenu32.exe C:\turtle\ C C (program name with 3 parameters on the command line). The 1st parameter you have to replace with the new directory name. (Important: If you forget the backslash at the end of the directory name it will be added!!! but please do not choose names with blanks in it!!!). What do the other 3 parameters mean ? The first C stands for colored histogram screen graphics (nC would be for only black and white). The second C means the same for printer outputs (nC for only black and white [no dithering] printer output). For low energy beam lines (Tandems, Ion Sources) it is desirable to have the possibility to include electrostatic elements instead of magnetic ones into the beam lines. The electrostatic QUAD has now been included the following way into Turtle (very close to magnetic QUADS): 16. 96. ; (switch to electrostatic QUADs) 5.vvv Leff. Ep. a. /EQ/ ; (Ep is the electric field [kV/cm] at the pole tip) 16. 95. ; (switch back to magnetic QUADs)- The used formula for kQ is (see Reference [9] p.99ff): kQ = sqrt( Ep * CHARGE / (a * RiE) ) with RiE = Ekin *(gamma+1)/gamma. (Ekin=particle energy in electron volts, CHARGE=particle charge in units of electrons) The transfer matrix element R11 e.g. is computed according to the following formula: R11 = cos(kQ * Leff). - The electric field gradient E is assumed to rise linearly from 0 at the axis to Ep = 2 * (Vp / a) at the pole tip with Vp to be the voltage at the pole tip with the distance a from the axis. - The dimension of Ep = 2 * (Vp / a) is fixed to kV/cm and may not be changed. - At the moment fringe field computations for this type of QUAD are not available like for magnetic QUADS. The dimension of this space has been increased from 1200 to 10000. The reason for this is the way multiple scattering in absorbers, collimators and degraders is implemented in Turtle. The thinner the slices of the above mentioned devices are, the more apertures are needed for a certain thickness (1 aperture per slice). Recent studies with a graphite degrader for 250 MeV protons have shown that the thickness of 1 slice has to be 0.1 mm in order to properly describe the energy straggling (range = 230 mm --> 2300 apertures). In order to reduce the amount of output in FOR002.DAT in the 'Summary of particles stopped by apertures' section a concatenation factor has been introduced:
16. 198. nnn. ; ( nnn consecutive slit losses are added to 1 value. )
The charge q of the particle (in electron units , default = 1) may be defined. The formula given on page 16 of reference [1] rho0 = p0 / (q * B0) ) was only implemented with q for one electron charge unit. Example: For He++ beams the line 16. 29. 2. ; should be inserted at the place of the stripper. If you start with a He- beam, do not forget to place a 16. 29. 1. ; line after the beam-card. The following type codes are influenced by q: 2, 4, 5, 11, 18, 19, 21, 26 and 27.
|
).