The Monte Carlo Race 2011

We saw earlier how different cars (and a box of Lego) could fit the description of various Monte Carlo particle transport codes. Every year there is the prestigious Monaco Grand Prix race (according to Wikipedia, I must stress, I am absolutely not into cars!), and this race passes the famous Monte Carlo region with its casino.

What you may not be aware of is the little known, sister event called the Monte Carlo Grand Prix, the year long publishing competition of Monte Carlo codes. I’m going to give you a run down of the latest results. Get a seat in the grandstand and plugin your ear plugs!


Monte Carlo codes are surging in radiotherapy treatment planning. For a course at the Aarhus University Hospital I had to give [footnote1] a little review of MC codes in radiotherapy. Inspired by a talk given by Emiliano Spezi I did a little study done using the ISI web of knowledge which clearly shows this trend:

Number of publications on “Monte Carlo Treatment Planning”. Data for 2011 are still incomplete for obvious reasons.

For a real race we of course need several competitors. In addition to the Geant4, FLUKA and SHIELD-HIT cars, we will have a few new participants joining the race!

From the country of The United States of America there are the two participants, MCNP and MCNPX developed in the secret labs of the Bombs-R-Us factory (a.k.a. Los Alamos), now released to the public (with a few exclusions) against a modest fee. MCNP helps physicists transporting neutrons as well as photon and electrons. The grand old MNCP brings along its more recent offspring MCNPX which also can transport heavy ions!

Just of north of the US, the Canadian National Research Council funded the development of EGS4/EGSnrc and as a courtesy to our medical minded people BEAMnrc was developed which is particularly suited for simulations of linear accelerators for radiation therapy. Rumours say that this nifty product is based on US technology from University of Wisconsin, ‘traitors’ some might argue, I’d say this is a remarkable example of pan American friendship!

The Japanese people, known for their industrial know-how, diligence and incomprehensible machine translated manuals contribute with this gemstone of technology: the PHITS multi purpose particle transport vehicle... sorry... code, maintained by the Japanese Atomic Energy Agency! Jovial attitude to technology breeds wonderful gadgets, often decades ahead of what eventually will appear in Europe.

Finally, and being among us for a while are the diligent efforts of our European friends in University of Barcelona in Spain. Hola PENELOPE!  PENELOPE can transport photons and electrons, evil tongues may even say much faster and smoother than our Canadian opponents. Let the masses judge how well they perform on the curved streets of Monaco.

Anyway, the race started many years ago and is still going on, let us see the current results:

Number of publications on various MC codes (not restricted to therapy planning)

Clearly MCNP has been roaring through the Mediterranean streets long before computers became mainstream and when 1 MByte was still an unbelievable amount of memory. It would have been a dull race, if EGS4 had not slowly appeared on the horizon. Geant4 was not conceived of when FLUKA already had an early start. Both Geant4 and FLUKA are developed at CERN, which makes their performance particularly exciting. Geant3 does not even appear here, it had a trauma from an earlier rally - it might recover, but bones tend to heal slowly at that age.

The fiery tempered Spanish PENELOPE and the latest ace from the United States MCNPX entered the rally roughly at the same time, and here the race gets really exciting! Let us take a closer look at the last 10 years. All codes are surging, so I will plot it relative to the total amount of publications per year:

The Monte Carlo Particle Transportation Grand Prix, during last 10 years.

It’s clearly a close race between the elderly and not-so-elderly codes. MCNP is clearly losing against its younger offspring MCPNX which also offers heavy ion transport. However, it is difficult to hide that MCNPX is loosing towards Geant4. The good news for MCNP is that it will be merging the MCNPX code base so it will get a boost of sorts! The success of Geant4 is truly amazing! Even if Geant4 cars are complicated, DIY affairs, the possibility to have total freedom with modding your car, giving it personal style, seems to attract a solid customer base.

FLUKA clearly shows a remarkably constant relative user base, going steady, going strong. But woo-hoo... whats happening with EGS4? The medical physicists came out of their closet declaring themselves as BEAMnrc users, and EGS4 is superseded by EGSnrc... 2011 could be one of the last years we hear of EGS4, but no-one will really miss it, we’ve got EGSnrc and BEAMnrc.

Trailing still is PHITS. PHITS users have still not published this year (2011) at the time writing these lines, but obviously I would not be surprised if it suddenly might come into play. Hopefully Japanese physicists will find time to catchup after they fix their issues with their nuclear power plants.

Where are the Russians? I can see smoke coming up on the horizon somewhere behind, I hope it is not a fire, but just a demonstration of the new engine. I am confident the Russian mechanics (with some aid from their friends in Aarhus :-) are working on it, and we soon will see fast SHIELD-HIT Niwas with whining tires and roaring physics engines in the curved streets of Monte Carlo.

The race isn’t over yet, and frankly I am not sure if it is supposed to end at all... but what can be said, some codes will prevail, others may lag and eventually disappear from the race.

See you next year at the Monte Carlo rally 2012 - readers are welcome to place their bets for 2012 below !



[footnote1]
“You have to give a 5 ECTS course in dosimetry.” Yeah right, I thought, it took me almost 4 months to prepare and run it, do all the exams (400 pages + external reviewer), no time for research whatever, looking forward for the day when it was all over. What followed? Research? Guess again: more teaching!

The "ideal" Monte Carlo user interface

If you have read many of our posts, you probably know that radiation Monte Carlo software plays a large role in both Niels' and my research efforts. Niels wrote an insightful and entertaining overview of the available Monte Carlo engines relevant for particle therapy. I'm going to talk about one of the things that most MC users have surely thought about or at least been frustrated by: user interfaces.

I'll start by making a disclaimer. Monte Carlo for radiation / particle transport is not a simple problem. It takes many person years worth of effort to produce a codebase that gives reasonable results. This post is mostly my musings about some ideal user interface. I am the proverbial beggar who also wants to be a chooser :)

Most interaction with MC programs is performed by creating input files of some sort (in Geant4 these are called macro files and can function in largely the same way). Usually these input files allow you to set the majority of important parameters relevant to your simulation (e.g. beam parameters, target geometry, detector geometry, scoring parameters, etc). While initially cryptic (see Fluka and MCNP input files), these text input files are extremely flexible and can be generated programatically or with a GUI. The problem of course, is that your typical input file format is not very intuitive, designed to be machine parsible, rather than human friendly. And as all new MC users know, input files and their syntax are one of the major hurdles in getting up and running.

A Fluka input file.

So what would the ideal MC user interface look like? Radiation physics users of MC engines want to irradiate something and score some quantity. For medical physics users, that usually means irradiating a human or water phantom and scoring dose or fluence. To me the obvious interface for MC codes would be identical to 3D CAD (computer-aided design).

The open source FreeCAD CAD program.

A 3D CAD style interface would put you directly in the geometry of the simulation world. Build your target from simple shapes or import existing geometries (e.g. DICOM files), graphically designate your beam source, type, and direction, and set up your detector geometries. More importantly, you would be able to manage your simulation end-to-end in the interface.

It can be argued that 3D CAD is as hard or harder to learn than a given MC engine. My approach for a user interface would be to expose the minimum useful controls, making advanced options discoverable through menus and configurable with shortcuts.

The follow-up disclaimer is that 3D CAD is also not an easy problem, so we are unlikely to see this soon. In fact many of the MC programs can import geometries from CAD programs (see SimpleGeo), but I'm unaware of any that have attempted to fully integrate a CAD-style GUI as a primary user interface.

What's your ideal Monte Carlo user interface? Leave us a comment and let us know.

Notes from the 2011 Geant4 Winter Course Tutorial

In January I traveled to Texas A&M University in College Station, TX, USA to attend the Geant4 Winter Course Tutorial to brush up on my Geant4 skills. Here are some of my impressions from the tutorial. The tutorial was from 10-14 January and, despite being in Texas, was on the cold side. Though hosted by the Texas A&M Dept. of Nuclear Engineering, the tutorial itself was held in the main auditorium of the Texas A&M Institute for Preclinical Studies (TIPS). There were about 100 attendees from all over the world, with a sizable fraction from A&M.

So what goes on at one of these tutorials? As described in Niels' earlier post about the different Monte Carlo "codes", Geant4 is not a program, per se, but rather a toolkit, primarily consisting of a set of C++ libraries and data files. The tutorial was aimed at bringing Geant4 newbies and post-newbies up to speed and slightly beyond. To accomplish this, we attended classes for five days in the stadium style auditorium with lectures covering a myriad of topics. Each day had one or two "hands on" sessions, in which we'd work through examples and have our questions answered. In many ways, the hands on sessions were the most useful, because of the one-on-one help. The actual exercises may not have been terribly useful, but having the time to ask questions with laptop-at-the-ready was when things came together. I certainly asked a few questions that were completely unrelated to the tutorial (but not unrelated to my research!) (Thanks again Sebastian!). The other thing I found extremely helpful was simply having time to discuss things with the Geant4 developers and other users. On the last day there were break-out sessions for medical physics, high energy physics, and on DNA damage. I attended the medical physics session with about 30 other people.

Other random notes:

• The attendees were a good mix of medical physics, nuclear applications, high energy physics, and others.
• I met a number of people involved with proton therapy projects as well as a post-doc from Wayne State working on their re-commissioned fast neutron therapy project.
• To keep small animals outside of TIPS, they had several poisonous plastic rocks outside the building (yes, really).
• We got a nice tour of TIPS, which included seeing the "most powerful" PET/CT scanner in existence.
• Catering by Jason's Deli all week
• Possibility of using cloud computing for distributed Geant4 was mentioned in a lecture by Asai. This is a subject of one of my projects (see arXiv:1009,5282v1 [physics.med-ph]).
• Looked a lot like a web surfing conference.

Overall it was a good tutorial. We learned that we were taking some wrong (or at least more difficult) paths in our code. Certainly Geant4 is a vast topic and a tutorial like this can be very helpful, if only to meet the right people for when you need to ask questions.

One Stopping Power Library to Rule Them All: libdEdx

Stopping powers describe the energy loss of charged particles traversing matter. In particle therapy, stopping powers are an essential ingredient for calculating the dose distribution of ion beams.

A direct way of calculating the stopping power is using the Bethe equation. However, this equation requires a good knowledge of the ionization potential (in particle therapy jargon: "the I-value") and for compounds this value is pretty ill researched. The I-value can be found experimentally, but there are only few experimental data available for compounds relevant for particle therapy, such as various tissue types and even for water. Our postdoc. Armin Lühr has recently submitted a paper where he takes a closer look on these issues which is available on the archive: http://arxiv.org/abs/1010.5356

Several programs provide stopping powers, either as analytic calculations or inter/extrapolated experimental data. The list below may be incomplete, but these are the codes I have been in touch with so far:

• ESTAR: for electrons. (Fortran77)
• PSTAR: for protons. (Fortran77)
• ASTAR: for alpha particles (Fortran77)
• MSTAR: for alpha particles and heavier ions up to Z=18. Basically MSTAR is scaling ASTAR data, where the scaling factors are fitted to experimental data. (Fortran77)
• ATIMA: code developed by GSI. (Java wrapping a Fortran core)
• TRIM/SRIM: Application which can simulate any ion on any compound.

And various tables provided by the ICRU:

• ICRU 49: proton data and alphas, equivalent with PSTAR/ASTAR.
• ICRU 73: ions heavier than helium, old version
• ICRU 73: new version after errata, quite similar to MSTAR, even if calculated differently.

and in summary they don't agree too well at lower energies (say below 10 MeV/u).

Nonetheless: from an application developer point of view this leaves you in a dilemma of choosing the proper code for your application, and most likely we will see more updates, since the issues are not solved yet with the erratum of the ICRU73.

In addition, from a technical point of view most stopping power codes are not written in convenient ways which allows a clean integration into your application. This is a common thing which happens when physicists develop code: the code will always be optimized to work nicely on the developers computer and fulfill their specific needs. Platform portability, an API, and install scripts following de-facto standards are mostly absent.

This is why we started the development of libdEdx.

libdEdx is meant to be a platform independent stopping power library which contains multiple tables to choose from, and can be extended with additional algorithms/tables. Currently libdEdx includes ESTAR,PSTAR,ASTAR,MSTAR and ICRU73 data as well as an implementation of the Bethe-equation akin to that found in SHIELD-HIT.
The material database is based on the extensive ICRU set found in the ESTAR table, and using Braggs additivity rule libdEdx can extend the MSTAR/ICRU73 data set to cover the ESTAR material list.

libdEdx comes along with an installer based on cmake, and the first official release should be able to run on Windoze just as well as Unix/Posix/Linux compliant platforms.

Surely, Geant4 also offers a range of stopping power tables to choose from, but you do not want to install entire Geant4 just to access the stopping power values.

Our bachelor student Jakob Toftegaard (who very conveniently has a background in physics and computer science) did most of the coding. A very early experimental pre-release is available on sourceforge, but if you want to experiment with it and even contribute you can also grab the most recent version from the SVN repository.

P.S.: Armin and David will present libdEdx on the MC2010 conference in Stockholm on the 9th to the 12th November. Posters can be found here.

Monte Carlo Programs in Particle Therapy Research

- a note on software design.

In my research, I use Monte Carlo particle transport codes a lot in order to simulate the interactions of a particle beams with matter.  There are several of these particle transport codes which are available for free, but if I have to simulate ions heavier than protons, then there are only few codes available which are capable of simulating the physics processes at clinical relevant energies. The four most common codes are FLUKA, Geant4, MCNPX and SHIELD-HIT.

Each of these transport codes have their scope of applicability, advantages and disadvantages. Here, in our research group APTG we work with FLUKA, Geant4 and SHIELD-HIT. In fact, we are actively developing the SHIELD-HIT code (visit our SHIELD-HIT developer page). Recently, I had a discussion with a colleague (and non-programmer) where I wanted to clarify why we put so much effort into SHIELD-HIT, now that basically all functions are more or less available in FLUKA or can be build in Geant4. I found it difficult to explain him the reasons.

That is why I came up with the idea of comparing the particle transport codes with cars. I am not particularly fond of cars, on the contrary, but in the spirit of particle transportation (pun intended) I could not resist the temptation.

A FLUKA car, Fortran 77 style.

A FLUKA car would be a regular car which does the job it is supposed to do. Imagine VW Passat or whatever. If you enter such a car for the first time, it will take you less than 10 seconds to orient yourself, since things are as expected: clutch, throttle, brakes etc.. If you forget to fasten your seat belt, you will get an error message which tells you exactly what is wrong. If you need some extras which goes beyond the standard equipment, you can add them yourselves (by the FLUKA user routines), although the developers limit what you can attach to your car. You are not allowed to touch the engine at all, the developers tune it for you to best possible performance.
FLUKA is closed source, and black-box like. If something fails with the engine (the physics models inside FLUKA), and this happens very rarely, you should take the car to repair ( = notifying the developers). Only the developers have access to the engine interior, and can fix it. While they do this, you realize that the interior may look rather old fashioned (the source is written in Fortran 77), and imagine it might be difficult to maintain, but the developers are experts, and have worked with it for many years. No reason to change this, as long as the engine runs smoothly and you still can find the special gasoline for the car (= a Fortran 77 compiler, deprecated on many newer Linux distributions).
If you want to publish benchmarks against other cars (other codes or experimental data), you must first talk back to the developers, and they want to assure that you have been operating the car right. (Unlike removing the brakes, then publishing how the car crashed.) The FLUKA car is increasingly popular, since it has a steep learning curve, is easy to run and the driver does not need to know how the engine is working.

A Geant4 car. Geant4 provides you all the pieces you need to build a car. ANY car.

Geant4 is not a car. It is a large box of Lego where you - in principle - can build your own car. We are not talking about normal standard Lego bricks, but the most fancy kind of them, Lego Technic C++ style bricks! If you ever have programmed in C++ and played with Lego, you know that equally to Lego you can also attach these code bricks in ways which they are not really supposed to (weak typing). Installing the Lego bricks in your laboratory may require some expert knowledge.
No-one has ever build a Lego based car similar to the FLUKA car, where you just step in, start the engine and go places. But there is no doubt, given infinite resources, you can build yourself a Ferrari, if you want. You can build a diamond inlaid SUV if such things turns you on. In principle anything is possible, but you need a lot of good developers and plenty of time. No research group has access to indefinite money supplies (or time), so instead research groups using Geant4, focuses on their specific needs. For instance, a famous research group in Boston have developed an engine (doing protons), which works well for their specific needs. In Japan a group is working on a special vehicle, lets say a Caterpillar (think gMocren is a part of this), which eventually can do specialized tasks (treatment planning with ions). That is fine, but if you want to adapt it to your own needs which goes beyond the application originally was designed for, it can again be quite some task and requires a good deal of programming knowledge. (Assuming you get the source code at all. If you get it, you still need to understand what is going on).
Perhaps some group has already developed higher level parts such as a carburetor and a light generator which are available to other researchers, yet they are not obliged to give these parts away. Geant4 is not GPL.

A SHIELD-HIT type car. Note, the driver found the light switch.

SHIELD-HIT is a russian Niva. Well, not precisely a regular Niva. At first glance you think, this is a normal car, not unlike a FLUKA car. You have to get a special contract with the code owners at INR before you may access the car.
Once you got it, you feel confident you can run this car effortlessly. However, as soon as you get inside, you realize something is very different. First of all you need three different keys to start the engine. The clutch is mounted on the steering wheel, the light switch is hidden under the seat, and if you need to do left turns, you must configure the car to do so before you start it. You need an English speaking Russian to tell you all this, because essential parts of the manual is written in Russian, and the manual itself is incomplete, only covering the light switch part.
If you look on the engine, you again realize it has lot in common with the FLUKA engine (SHIELD-HIT is also written in Fortran 77, and the geometry parser is also CG).
The key difference to the FLUKA car is, if something is broken, you are allowed (or even encouraged) to repair it yourself. Clear error messages or other indications of what is wrong, are seldom. You may simply have operated your Niva wrong, of there may be a real bug in the engine. But you do have access to the source code - and this enables you to do all the modification you want in the code. If you know exactly what you want, this is actually a big advantage. Imagine yourself stranded in a village in the middle of Siberia, you will be happy you drive a Niva. When your Niva is fixed, the car runs smoothly - just as smooth as a FLUKA car or a custom Geant4 based "forward going vehicle" would do. And the Niva is tolerant to various gasoline types (you can compile with g77, gfortran, ifort etc...)

So, our MSc student David has invested a lot of time in building the next generation Niva, the Niva version 2.0 (currently also known as SHIELD-HIT10A). It is supposed still to resemble a regular car but with improved user friendliness . The clutch and light switch are moved to more intuitive positions. Only one key should be necessary for starting the engine. And, as I mentioned in an earlier blog-post I am currently preparing an English manual for the new Niva 2.0.

David is also benchmarking the Niva 2.0 and overhauling the engine (meaning, better parameters for the physics models matching the recently released experimental data by E. Hättner). Finally, the Niva 2.0 will feature a lot of new features which are not readily available in any other Monte Carlo particle transport vehicle, such as air conditioning and cup holders (one for each passenger).

SHIELD-HIT10A, aka Niva 2.0. Fortran 77 style artwork is still visible despite of the upgrade. (APTG developer impression).

Now, let us assume, when David finishes his MSc February 2010, he would like to do a PhD project involving a pair of windscreen wipers and a light generator.
He has three codes to choose from. Of course, we all like playing with Lego. Lego is fun, and someone has already developed a Lego light generator ... so we only need to build the windscreen wipers, and off we go, doing a lot of research. (I really get carried away now.)
Alternatively David could choose to continue to work with the Niva 2.0. As a side effect the Niva will be upgraded with windscreen wipers which may benefit the continuation of the SHIELD-HIT project. Not as much fun as working with Lego, but probably as useful. Personally, I'd hit it, but working with Fortran77 is really demotivating. Especially if the code is full of GOTO statements.
Finally we have the choice of working with FLUKA. Perhaps this is the easiest, but in my opinion a bit dull since we are not contributing with much new on the developer side of the code.

There will be more of the Niva 2.0 at the MC2010 conference in Stockholm, where David will present his work (either poster or talk, we don't know yet).

Come and meet us!