Roy and Niels

Roy and Niels

Thursday, March 6, 2014

A Stopping Power App for Android



We (as in APTG) have just released the dE/dx Android app on Google Play!

Screnshot of the dE/dx app, ver. 1.0 from my LG Nexus 4 phone.

The app is written by student Casper Cristensen in collaboration with the Dept. of Physics and Astronomy and the School of Engineering at Aarhus University. The app is based on the stopping power library libdEdx-1.2.1  which was presented in an earlier blog entry here.

The cool thing about the app is that all tables are stored locally on phone (this is why it needs access to USB memory), so no network access is necessary - useful if you want to look up stopping power while working in a shielded area with no network.

In summary the features are:
  • ICRU 49 and ICRU 73 (revised version) electronic stopping power tables
  • MSTAR tables provided by Helmut Paul
  • BETHE_EXT00, a generic algorithm as used in SHIELD-HIT12A. Applicable for any ion to all 278 ICRU compounds (beware, use with care, can be off at lower energies), .
  • no network access needed after installation
  • calculation of CSDA range of particle in target material
  • inverse CSDA lookups possible (if you need to figure out what specific energy is needed to produce a certain range in a material.

Disclaimer: of course we do not claim correct reproduction of any of the data, use at own responsibility!

Thanks to:

  • Casper Christensen (programming this app)
  • Jesper Rosholm Tørresø (being Casper's supervisor and our contact person to the engineering school)
  • Jakob Toftegaard (writing most of the libdEdx backend and dedx.au.dk webpage)


Enjoy! :-)

Sunday, October 20, 2013

I Just Bought a Dosimeter

I recently bought a cheap personal dosimeter by the Ukrainian company SOEKS. The best time to buy these is, when there has *not* recently been a nuclear disaster - prices are more reasonable then. I went for their "Defender" model, since it is capable of measuring accumulated dose, and I do not need fast counting provided by those devices with two GM-tubes.

The SOEKS Defender dosimeter. www.soeks.ru
Ok, I said "dose", but in fact I am not sure what it measures, it should be the operational quantity "personal dose equivalent", but it could also be something completely different (see ICRP 103). In fact I would not be surprised if they do some sort of "dose equivalent" calibration following ICRP 26 protocol from 1971. Many still confuse this.). Just for now, I'll just call it "dose", well knowingly that this is wrong.

It is the first time I order stuff from Ukraine, but shipment was really professional. They even email you a picture of your parcel with your address on, before they ship it. Contrary to some of the reports, I had no trouble with customs either. It was directly delivered to my institute address with no fuzz. Good thing!

Unboxing it, I found a nifty small gadget, which fits in a pocket, It comes with rechargeable AAA batteries (SOEKS brand, these people are really dedicated!), which can be recharged via the USB connector. I might have preferred AA batteries, as these simply perform much better than AAAs. But still, nice that you can always switch to normal batteries if you are in the field, run out of power and have no time to recharge.

One thing I really miss though, is a loop for attaching a lanyard. :-/

Switching it on, I immediately measured the normal background radiation levels here in Aarhus, around 0.12 µSv/hour.

One of the first things I tried out, was to take it along a trip to Turkey (yeah, the NUFRA2013 conference), to see what dose I'd receive during the flight.

I left it turned on while my cabin luggage was X-rayed. I did this three times, and each of these measurements gave 0.8 µSv.

On board the airplane I could measure a background level of roughly 2 µSv/hour, triggering the "Dangerous radiation" warning of the dosimeter. The threshold for the warning can be adjusted by the user. The 2 µSv/hour is probably a very conservative measurement anyway, as the device does not count the contribution from neutrons.

Video documentation of that flight (at roughly 10 km altitude):


Dose rate increased rapidly when we came above 8 km altitude. Under 8 km, it was pretty much just the background radiation. For the 4-5 hour trip from Antalya to Amsterdam I measured a total dose of 8.19 µSv. (This is 10 times the amount my luggage received when it was X-rayed. Very interesting.)

Arrived at Shiphol, Amsterdam from Antalya. Flight time 4 hours 47 minutes. Accumulated dose 8.19 µSv.
At our institute I can test the dosimeter at several sources. None of them are calibrated in any way, it is merely to see some signal. The strongest sources we have at the institute are couple of Americium-Beryllium sources for neutron generation. Whenever you have neutron fields you also have plenty of gammas, and if the sources is packed in a hydrogen rich moderator material there are also energetic protons (recoiling hydrogen nuclei).

The AmBe source sown below is locked up in a safe, but even when the safe is closed, you can see an elevated photon background in a meters distance:

Elevated background radiation near the neutron cabinet.

Inside the safe, the AmBe source is located in a paraffin moderator. The source is used for various activation experiments, where you can lower the materials to be activated into the two plastic tubes.

The AmBe neutron source is hidden in the paraffin block. Through the tubes you can immerse various materials for neutron activation analysis.

Here, I recorded the highest dose rate recorded so far, 320 µSv/hour. Again, this is probably not a very correct measure, since that area is flooded with neutrons as well. How neutron interact depends strongly on the target material.

320 µSv/h. Occupational limit is 20 mSv/year, public limit is 1 mSv/year, so let's quickly close the door again.

We have several dosimeters at the institute, so I tried to compare the elevated background with a calibrated device, the quite common and more pricey "Gamma scout" dosimeter:
Comparison of the SOEKS Defender against a Gamma Scout. In fact this was the best agreement I ever saw between the two devices.
They both measure the same background to approximately 0.7 µSv/hour, but the agreement varied a lot depending where in the room I was. Mostly, I noticed a considerable underestimation of the dose reported by the SOEKS device, sometimes up to a factor of 2. Probably this is just the precision you can get, it's good enough for my (non-work related) purposes. If you want high-precision monitoring, you need professional equipment which is well calibrated, and that is certainly a different price tag.

I also have a little Japanese silver medallion. I activated it by leaving it in on top of the AmBe paraffin block. Neutrons generate Ag-110 and Ag-108 isotopes with T_1/2 of 26 sec and 2.4 min, respectively. I could follow the cooling down of the medal with the dosimeter:

Japanese silver medallion activated by neutrons. After 15-30 minutes it cools down to normal background levels.

Various sources for student exercises.

Anyway, I still hope to find time to realize my Chernobyl visit, will certainly bring this handy device with me. :)


-------------------

P.s.: Not really related, I recently came across a very old video of how I in my young ages zapped a CCD detector using a medical linear accelerator. The chip suffering in the video is an engineering grade CCD47-20 by EEV (E2V / Marconi), which was intended to fly on a satellite mission. I later also irradiated one with protons. (Sorry, video is a polyglot mix of German and Danish.)

Wednesday, November 21, 2012

SHIELD-HIT12A demo version released

Here is a little sneak-preview of the upcoming release of SHIELD-HIT12A.

Pimp my Niva... (artist's impression of SHIELD-HIT12A).

SHIELD-HIT12A is a Monte Carlo particle transport code capable of transporting heavy ions through arbitrary media. The -A fork was made in 2010 from SHIELD-HIT08, and since then we have added plenty of new features, solved many bugs, increased calculation speed and optimized the nuclear models to new data on carbon-12 fragmentation.

A free demo version where the random seed and statistics are fixed to 10.000 particles can be downloaded from the project development page. There are builds for Linux and Windows systems (32- and 64 bit). It is a beta-release, but we would like to bring this demo version to a broader community, and thereby hopefully also fix some more bugs before we release the full version.

The new features in SHIELD-HIT12A are:
  • New simplified material and beam parameter parser in free format and extensible without breaking downward compability.
  • Including 279 ICRU default materials and elements, it has never been so easy to specify a material.
  • New beam model: divergence and focus distance can now be specified (thanks to Uli Weber from Marburg).
  • Arbitrary starting beam directions now possible.
  • New routine for Vavilov straggling, 5-6 times faster than the original one by Rotondi and Montagna which was used in Geant3.21. In total, this alone means a speed improvement of roughly 30-40%.
  • Ripple filter has two modes of operation, Monte Carlo type or Modulus type.
  • Logarithmic energy binning in SPC files for TRiP
  • Full howto for generating DDD files for TRiP
  • Now only three input files are needed to setup a run.
  • Improved documentation.
  • Scoring by zones using detect.dat (complementary to Cartesian mesh and cylindrical scoring)
  • Alanine response model included, so SHIELD-HIT12A can directly calculate the dose equivalent response in alanine.
  • Flat circular and square beams can be defined
  • Neutron data for natural Argon was added (needed for detailed simulations of air)
  • Another ton of bug fixes.
Of course SHIELD-HIT12A includes the features from SHIELD-HIT10A (which was never released):
  • Totally new (parallelizable) scoring system:
    • Arbitrary Mesh and Cylindrical scoring
    • Lots of detectors such as energy, fluence, dose-averaged LET, track-averaged LET, average velocity (beta), dose to medium (where medium can be changed if you want to calculate stopping power ratios) etc....
  • finally SHIELD-HIT10A is parallelizable
  • New random number generator, which gives a massive performance boost
  • New adjusted inelastic cross sections for carbon ions based on recent data
  • Fine tuning of the fermi-breakup parameters
  • SHIELD-HIT10A can be configured without accessing the source code anymore, so no programming knowledge required to use SHIELD-HIT10A.
  • SHIELD-HIT10A is installable
  • Runs on linux again, even when compiling with code optimizations, ok with GNU gfortran, Intel and Portland compilers.
  • Many bug fixes
Enjoy! And please drop me a line when you find bugs in the software and errors in the manual, they are there, but we hid them well. :o)


----

Friday, November 16, 2012

Medical physics papers on the arXiv: Some stats

Last year I made a post about the state of open access in medical physics and the use of arXiv.org to make medical physics papers freely available to everyone. Today I've decided to follow that up by taking a look at what is going on over at arXiv.org and sorting through some of their data.

The arXiv is "an openly accessible, moderated repository for scholarly articles in specific scientific disciplines". Beyond simply allowing all comers to read and download the articles they host, the arXiv also makes information about its vast collection easily available via a set of APIs. I decided to make use of the API to download meta-data (i.e. author, title, abstract, etc) for all articles in the medical physics category. Let's take a look!

First some raw numbers. As of mid-November 2012, the physics.med-ph category had (about) 980 articles. Those 980 articles where co-listed in 76 of the possible 126 other arXiv categories along with the medical physics category. In Figure 1 I've plotted the number of submissions by year (with a partial count for 2012). It's clear that the submission rate to physics.med-ph has greatly increased. On the plot I've fit a logistic growth curve (Gompertz), as I'm assuming that the submission number will saturate at some point. You can see faint exponential and linear fits as well. The logistic model predicts 174 submissions for 2012, but 140-145 seems more likely this year.


Figure 1. Submissions per year to the physics.med-ph category on arXiv.org.

Another item of interest is how the medical physics category fits into arXiv. If you've browsed physics.med-ph before, you know that it contains a broader range of topics than the popular "medical physics" journals, such as Physics in Medicine and Biology or Medical Physics. As mentioned above, many of the 980 articles were listed in multiple categories (663 to be exact). Figure 2 shows the most popular co-categories for articles in physics.med-ph. As might be expected, biophysics (physics.bio-ph) was the most popular co-category, along with topics that seem well aligned, such as instrumentation and detectors (physics.ins-det), organs and tissues (q-bio.TO), and computational physics (physics.comp-ph). But less obvious categories were also co-submitted as well, such as chaotic dynamics (nlin.CD).

Figure 2. "Co-categories" of papers submitted to physics.med-ph.
To try to get a better idea of how these categories interplay with one another, I made some simple network graph visualizations. Figures 3 and 4 show the connections between the co-categories. All of the papers are clearly in physics.med-ph and all other categories are linked to that. The other lines in the network plots represent when a paper was simultaneously in more than two categories (e.g. medical physics, biophysics, and physics - data analysis).

Figure 3. Network graph of categories for papers submitted to physics.med-ph. Click for larger version.
Figure 4. Network graph of categories for papers submitted to physics.med-ph. Click to view larger version with category labels on the nodes.
In Figure 4 the region near the center of the graph are the categories that are the most likely to be co-listed with other categories. As you might expect, the categories with a dense set of lines connecting them tend to also be the most frequently occurring categories, as seen in Figure 2. This is more easily seen in the larger version of Figure 4 with category names. (click the above figure).

As you can see the arXiv physics.med-ph category as a wide ranging and growing category for open access articles related to medical physics. It will be interesting to see how it fits in with the wider trends of funder mandates for open access and the general growing acceptance and demand for open access in out community.

N.B. This arXiv meta-data is relatively easily pulled down and processed with Python tools using the arXiv API, going from XML to JSON to your computer screen at home!

Thursday, June 21, 2012

libdEdx 1.2.0 released - stopping powers for the masses!


Our master's student Jakob Toftegaard has been very busy lately, and we're now ready to release a new version of the open-source stopping power library libdEdx - version 1.2.0.

Changes to the last version 1.0 are
  • first of all: a new, nice and clean API. This breaks compatibility with 1.0, but we will do our best to avoid this happening again in the future. (We now have extensive use of structs, which can be extended with new members.)
  • generic ICRU table included, which combines ICRU49 and the revised ICRU73.
  • all four calculation modes in MSTAR are now supported, the default is that recommended by Helmut Paul.
  • aggregate state can be specified, following ICRU recommendations
  • I-values can be overridden for analytical functions (BETHE_EXT)
  • provides a bunch of new functions
    • calculate CSDA range
    • inverse range lookup - a given range will return required particle energy in CSDA approximation. (Yes, this is the feature you have been waiting for all your life!)
    • inverse dEdx lookup - a given stopping power will yield an energy (either high or low value, depending on what the user requested)
  • version string of libdEdx can be accessed
  • memory leak fixes
  • typo fixes in material lists
  • code should now be thread-safe
And were just getting started!

A real gem is our new web based frontend at http://dedx.au.dk to libdEdx. Here you can lookup stopping power functions using various tables and energies, use it as a supplement to those from NIST. The website includes a nice plotting feature as well, where you can add multiple plots for comparison.

The web frontend is still in beta-testing phase and may reside in this state for a long time. Any feedback is appreciated.

We have a lot of plans regarding how to continue with this. In the next version of libdEdx we plan to include more features such as
  • more Bethe-based stopping power functions such as
  • Bethe-Bloch
  • Bethe-Bloch-Barkas
We also would like to include additional stopping power programs such as ESTAR, ATIMA and SRIM (just the stopping power part, of course), yet the outcome will depend on the willingness of the respective authors to contribute.
  • algorithms for nuclear stopping power
  • … and a surprise which we won’t reveal yet … :-)
Other contributors will be most welcome, the project is available on sourceforge for inspection.

We greatly acknowledge our hero Helmut Paul for contributing to the development with very fruitful discussions and suggestions. We also acknowledge the permission from the ICRU to use their stopping power tables in libdEdx.

Yet, we do not claim that the produced results are correct in any way, so any use of the data are on own risk. Nonetheless, if you DO find discrepancies, errors, misbehaviour of code, we would really appreciate if you tell us.

Enjoy!

Friday, June 1, 2012

Don't Mention the War

Last year the Danish government announced that it will finance a national particle therapy facility together with some private funding agencies. Two applications were submitted for evaluation:

and of course there is a strong interest from each side to host this project. Naturally, since I am part of the Aarhus Particle Therapy Group, I will right out admit I am naturally in favour of our project in Aarhus.

The plan is now that an international expert commission is being setup by the ministry of health. That commission will then evaluate those two applications and give a suggestion where to place such a facility.

If you compare those two applications, you will see that they both apply for a 3-gantry treatment facility. Copenhagen lock themselves on a cyclotron solution, where we in Aarhus are also open to a synchrotron based solution.

Basically the differences can be summarized here:

Copenhagen Aarhus
3 treatment rooms with 2 gantries 3 treatment rooms with 2 gantries,
+1 research room
3rd gantry will be installed later 3rd gantry will be installed later,
a fourth treatment room with gantry can be added
ProtonsProtons
CyclotronCyclotron or Synchrotron
1000 patients per year1000 patients per year


So basically what we apply for is rather similar.

What about the costs? If you compare those budgets, you will note that the prices of the equipment and the building are comparable too.

Copenhagen Aarhus
Equipment450 mio. DKK 475 mio. DKK
Building382 mio. DKK 295 mio. DKK
Supply systems and lines195 mio. DKK (incl. in building)
Purchase and demolition of Rockefeller Building147 mio. DKK (none)
Total1174 mio. DKK 770 mio. DKK
Optional 3rd gantry(not mentioned) 75 mio. DKK

Copenhagen need to build in the middle of the city, whereas we in Aarhus have a pristine green field adjacent to the New University Hospital. This means that the expenses for the Copenhagen proposal are significantly higher than the Aarhus solution, since a building has to be demolished first (they have simply a lack of space in Copenhagen). The "supply lines" cover, as far as I heard, establishing extra power for the facility, please correct me if I am wrong here. The annual operational costs are also similar.

There has been a few articles in the Danish press about the two applications. None of the Journalists managed to actually compare the two budgets, so basically the message was that Aarhus offeres particle therapy at dumping prices. Yeah right.

In addition, I hear a lot of weird arguments, e.g. such a facility MUST be build in Copenhagen, since it is the capital city, and it has a proximity to the airport.

Huh?

How many of the existing facilities are located in capital cities? Lets pick some of the most famous ones: HIT is in Heidelberg (147.000 inhabitants, 1hours drive to next airport). PSI is in Villingen. The Swedish facility is built in Uppsala, and not Stockholm. Why is that so? These research institutions simply had the longest track record in terms of particle therapy! I don't get the point of the closeness of an airport either, since you can reach any point within a few hours in DK by car, and sorry guys... particle therapy is not that acute.

Now, today this conflict took a very curious turn, as Copenhagen announced in Dagens Medicin that they will cut the price 75 % (login required to view link). What is going on?
It seems that Copenhagen suddenly decided to go for a single gantry facility by Mevion, formerly known as StillRiver. To me this looks like a very poor decision. The Mevion facility has been controversial for a long list of reasons. In a scientific article by Schippers and Lomax they list some of the issues with this facility:

  • the very compact design of the cyclotron, may lead to large beam losses and resulting activation. 
  • repetition frequency of synchrocyclotrons are typically 1kHz which limits applicability of the different pencil beam scanning methods.
  • There is no beam analysis and no magnets in the beam path, in other words an energy selection system is missing. This means that the beam will have a poor distal edge, always looking as if it traversed 25-30 cm of material, e.g. when treating head an neck cancers. This limits possible treatment planning techniques, such as patch fields.
  • Neutron contamination may be an issue, due to the proximity of the degrader, modulator and collimator(s). 
To this, I might add that the construction of the prototype is still awaiting FDA approval, and we have not seen any data on how this machine operates. We are just in the middle of a never ending scandal about a couple of diesel IC4 trains ordered from AnsaldoBreda a decade ago, where the trains still are not functional and with a significant risk that they never will be. Billions of taxpayer money are lost here. Most Danes are therefore allergic to order unproven technology abroad.

Update: Mevion just got their FDA approval.

Another update: it's a "Premarket Notification" and not a "New Device Approval". Read it here. Thanks to Klaus Seiersen for pointing this out.

The recent article in Dagens Medicin does not even mention that cutting 75% of the price leads to a single-gantry facility, which means the patient numbers needs to be adjusted down, perhaps by 2/3rds.

Finally, if you need three treatment rooms, then measured in costs per treatment room, a conventional solution is cheaper than the Mevion solution, according to Schippers.

I am very astonished, it seems that Copenhagen just scored an own goal. However, I fear this may delay the decision process even more. In the end this is about patients, and if the money is there, we should not be satisfied with giving our patients the most inferior kind of proton therapy, and rely on unproven technology.

Sunday, January 22, 2012

The Quick and Dirty Guide for Parallelizing FLUKA

(Single PC version)

Imagine you got a desktop or laptop PC with 4 or perhaps even 8 CPU cores available, and you want to run the Monte Carlo particle transport program  FLUKA on it using all CPU cores.
The FLUKA execution script rfluka however was designed to run in "serial" mode. That is, if you request to repeat your simulation a lot of times (say, 100) issuing the command rfluka -N0 -M100 example, each process is launched serially, instead of utilizing all available cores on your PC.

A solution can be to use a job queuing system and a scheduler. Here, I'll present one way to do it on a Debian based Linux system. Ubuntu might work just as well, since Ubuntu is very similar to Debian. A feature of the method presented here, is that it can easily be extended to cover several PCs on your network, so you can use the computing power of your colleagues when they do not use their PCs (e.g. at night). However, this post will try to make it very simple, namely set it just on your own PC. In less than 10 minutes you'll have it up and running...

The idea is to use TORQUE in a very minimal configuration. There will be no fuzz with Maui or similar schedulers, we will only use packages we can get from the Debian/Ubuntu software repositories.
In order to be friendly to all the Ubuntu users out there, all commands issued as root are here prefixed with the "sudo" command. As a Debian user you can become root using the "su" command first.

First install these packages:

$ sudo apt-get install torque-server torque-scheduler 
$ sudo apt-get install torque-common torque-mom libtorque2
and either
$ sudo apt-get install torque-client
or
$ sudo apt-get install torque-client-x11

after installation we need to setup torque properly. I here assume that your PC hostname cannot be resolved by DNS, which is quite common on small local networks. You can test whether your hostname can be resolved by the "host" command. Assuming your PC has the name "kepler", you may get an answer like:

$ host $HOSTNAME
Host kepler not found: 3(NXDOMAIN)

this means you may need to edit the /etc/hosts file, so your PC can associate an IP number with your hostname. Debian like distros may have a propensity to assign the hostname to 127.0.1.1 which will not work with torque. Instead I looked up my IP number (which in my case is pretty static) using /sbin/ifconfig, and edited the /etc/hosts accordingly, using your favourite text editor (emacs, gedit, vi...)
My /etc/hosts file ended up looking like this:

127.0.0.1 localhost
#127.0.1.1 kepler.lan kepler
192.168.1.108   kepler

If your hostname of your PC can be resolved, you can ommit the last line, but under all circumstances you must comment out the line starting with 127.0.1.1.


Once this is done, execute the following commands to configure torque:
$ sudo echo $HOSTNAME > /etc/torque/server_name
$ sudo echo $HOSTNAME > /var/spool/torque/server_name
$ sudo pbs_server -t create
$ sudo echo $HOSTNAME np=`grep proc /proc/cpuinfo | wc -l` > /var/spool/torque/server_priv/nodes 
$ sudo qterm
$ sudo pbs_server
$ sudo pbs_mom

(Update: If qterm fails, you probably have a problem with your /etc/hosts file. You can still kill the server with $killall -r "pbs_*".)

Now let's  see if things are running as expected:
$ pbsnodes -a
kepler
     state = free
     np = 4
     ntype = cluster
     status = rectime=1326926041,varattr=,jobs=,state=free,netload=3304768553,gres=,loadave=0.09,ncpus=4,physmem=3988892kb,availmem=6643852kb,totmem=7876584kb,idletime=2518,nusers=2,nsessions=8,sessions=1183 1760 2170 2271 2513 15794 16067 16607,uname=Linux kepler 3.1.0-1-amd64 #1 SMP Tue Jan 10 05:01:58 UTC 2012 x86_64,opsys=linux

and also
$sudo momctl -d 0 -h $HOSTNAME

Host: kepler/kepler   Version: 2.4.16   PID: 16835
Server[0]: kepler (192.168.1.108:15001)
  Last Msg From Server:   279 seconds (CLUSTER_ADDRS)
  Last Msg To Server:     9 seconds
HomeDirectory:          /var/spool/torque/mom_priv
MOM active:             280 seconds
LogLevel:               0 (use SIGUSR1/SIGUSR2 to adjust)
NOTE:  no local jobs detected

Now setup a queue, which here is called "batch".
$ sudo qmgr -c 'create queue batch'
$ sudo qmgr -c 'set queue batch queue_type = Execution'
$ sudo qmgr -c 'set queue batch resources_default.nodes = 1'
$ sudo qmgr -c 'set queue batch resources_default.walltime = 01:00:00'
$ sudo qmgr -c 'set queue batch enabled = True'
$ sudo qmgr -c 'set queue batch started = True'
$ sudo qmgr -c 'set server default_queue = batch'
$ sudo qmgr -c 'set server scheduling = True'

[update: you may want to increase walltime to 10:00:00 so jobs dont stop after 1 hour]

and start the scheduler:
$ sudo pbs_sched

The rest of the commands can be issued as a normal user (i.e. non-root).

Let's see if all servers are running:
$ ps -e | grep pbs
 1286 ?        00:00:00 pbs_mom
 1293 ?        00:00:00 pbs_server
 2174 ?        00:00:00 pbs_sched

Anything in the queue?
$ qstat
$ 
Nope, it's empty.

Lets try to submit a simple job
echo "sleep 20" | qsub

and within the next 20 seconds you can test, if its in the queue:
$ qstat
Job id                    Name             User            Time Use S Queue
------------------------- ---------------- --------------- -------- - -----
0.kepler                 STDIN            bassler                0 R batch


Great, now were ready to rock 'n roll! This is really a minimalistic setup, which just works. For more bells and whistles, check the torque manual.

All we need, is a simple FLUKA job submission script: rtfluka.sh
#!/bin/bash
#
# how to use this
# change to directory with the files you want to run
# and enter:
# $ qsub -V -t 0-9 -d . rtfluka.sh
#
#PBS -N FLUKA_JOB
#
start="$PBS_ARRAYID"
let stop="$start+1"
stop_pad=`printf "%03i\n" $stop`
#
# Init new random number sequence for each calculation. 
# This may be a poor solution.
cp $FLUPRO/random.dat ranexample$stop_pad
sed -i '/RANDOMIZE        1.0/c\RANDOMIZE        1.0 '"${RANDOM}"'.0 \' example.inp
$FLUPRO/flutil/rfluka -N$start -M$stop example -e flukadpm3

Update: Note that your RANDOMIZE card in your own .inp file must match the sed regular expression above, else you may repeat the exact same simulation over and over again...


Let's submit 10 jobs:
$ qsub -V -t 0-9 -d . rtfluka.sh

And watch the blinkenlichts.
$ qstat
Job id                    Name             User            Time Use S Queue
------------------------- ---------------- --------------- -------- - -----
15-0.kepler               FLUKA_JOB-0      bassler                0 R batch          
15-1.kepler               FLUKA_JOB-1      bassler                0 R batch          
15-2.kepler               FLUKA_JOB-2      bassler                0 R batch          
15-3.kepler               FLUKA_JOB-3      bassler                0 R batch          
15-4.kepler               FLUKA_JOB-4      bassler                0 Q batch          
15-5.kepler               FLUKA_JOB-5      bassler                0 Q batch          
15-6.kepler               FLUKA_JOB-6      bassler                0 Q batch          
15-7.kepler               FLUKA_JOB-7      bassler                0 Q batch          
15-8.kepler               FLUKA_JOB-8      bassler                0 Q batch          
15-9.kepler               FLUKA_JOB-9      bassler                0 Q batch 

Surely, this can be improved a lot, suggestions are most welcome in the comments below. One problem is for instance, that the random number seed is limited to a 16 bit integer, which only covers a very small fraction of the possible seeds for the RANDOMIZE card.
Update: There is also a very small risk that the same seed occasionally is used twice (or more often). Alternatively one could just add a random number to a starting seed after each run. (Any MC random number experts out there?)

Output data can be processed in regular ways, using flair
Alternatively you may use some of the scripts in the auflukatools package, which for instance can do the merging of USRBIN output with a single command. Auflukatools also includes rtfluka.sh as well as a CONDOR job submission script rcfluka.py, which is better suited for heterogenous clusters.

Finally, here is a job script for SHIELD_HITxxA, (which is even shorter):

#!/bin/bash
#
# how to use
# change to directory you want to run
# $ qsub -V -t 0-9 -d . rtshield.sh
#
#PBS -N SHIELD_JOB
shield_exe  -N$PBS_ARRAYID

Enjoy!

Totally unrelated: englishrussia.com just posted some nice pics from the Budker institute for Nuclear Physics in Novosibirsk, Russia. Certainly worth visiting, have a look at:
http://englishrussia.com/2012/01/21/the-budker-institute-of-nuclear-physics/
 :-) Heaps of pioneering accelerator technology was developed there, such as electron cooling, the first collider, lithium lenses (e.g. for capturing antiprotons), and they supplied the conventional magnets for the beam transfer lines to the LHC at CERN. I visited the center many years ago but my pics are not as good. :-/ The German wiki about Budker himself, is also worth reading.