Changeset 1313 for trunk/examples/advanced/hadrontherapy

trunk/examples/advanced/hadrontherapy/GNUmakefile

-                      r1230
+                      r1313
 # $Id: GNUmakefile,v 1.12 2009/08/13 20:48:04 kaitanie Exp $
+# $Id: GNUmakefile,v 1.17 2010/05/18 11:14:39 cirrone Exp $
 # --------------------------------------------------------------
 # GNUmakefile for examples module.  Gabriele Cosmo, 06/04/98.
 …
 G4EXLIB := true
+# Debug info
+#CPPFLAGS += -g
 ifndef G4INSTALL
   G4INSTALL = ../../..
+  G4INSTALL = ../..
 endif
 …
 include $(G4INSTALL)/config/architecture.gmk
+ifdef G4ANALYSIS_USE
+CPPFLAGS += -DANALYSIS_USE
+endif
+ifndef G4ANALYSIS_USE      # If we don't have AIDA
+ifdef G4ANALYSIS_USE_ROOT   # And we have ROOT
+CPPFLAGS += -DANALYSIS_USE -DG4ANALYSIS_USE_ROOT
+ifdef G4ANALYSIS_USE_ROOT   # If we have ROOT
+CPPFLAGS += -DG4ANALYSIS_USE_ROOT
 CPPFLAGS += $(shell root-config --cflags)
 LDFLAGS  += $(shell root-config --glibs)
-endif
 endif
 include $(G4INSTALL)/config/binmake.gmk
-ifdef G4ANALYSIS_USE
-endif

trunk/examples/advanced/hadrontherapy/Hadrontherapy.cc

-                      r1230
+                      r1313
 // ********************************************************************
 //
+// Hadrontherapy.cc
+//
+// Main of the Hadrontherapy example;
+// Released with the Geant4 9.3 version (December 2009)
+//
+// Last modified: G.A.P.Cirrone
+//
+// This is the *basic* version of Hadrontherapy, a Geant4-based application
 // See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
 //
+// To obtain the full version visit the pages: http://sites.google.com/site/hadrontherapy/
 // ----------------------------------------------------------------------------
 //                 GEANT 4 - Hadrontherapy example
 // ----------------------------------------------------------------------------
 // Code developed by:
+// Main Authors:
 //
 // G.A.P. Cirrone(a)Â°, G.Cuttone(a), F.Di Rosa(a), E.Mazzaglia(a), F.Romano(a)
 //
 // Contributor authors:
 // P.Kaitaniemi(d), A.Heikkinen(d), Gillis Danielsen (d)
+// P.Kaitaniemi(d), A.Heikkinen(d), G.Danielsen (d)
 //
 // Past authors:
 …
 #endif
+#ifdef G4UI_USE
+#include "G4UIExecutive.hh"
+#endif
 //////////////////////////////////////////////////////////////////////////////////////////////
 int main(int argc ,char ** argv)
+{
+  // Set the Random engine
+  CLHEP::HepRandom::setTheEngine(new CLHEP::RanecuEngine());
+  G4RunManager* runManager = new G4RunManager;
+  //Initialize possible analysis needs, needs to come early in order to pick up metadata
+#ifdef ANALYSIS_USE
+  HadrontherapyAnalysisManager* analysis = HadrontherapyAnalysisManager::getInstance();
+  analysis -> book();
+#endif
+  // Geometry controller is responsible for instantiating the
+  // geometries. All geometry specific setup tasks are now in class
+  // HadrontherapyGeometryController.
+  HadrontherapyGeometryController *geometryController = new HadrontherapyGeometryController();
+  // Connect the geometry controller to the G4 user interface
+  HadrontherapyGeometryMessenger *geometryMessenger = new HadrontherapyGeometryMessenger(geometryController);
+  G4ScoringManager *scoringManager = G4ScoringManager::GetScoringManager();
+  scoringManager->SetVerboseLevel(1);
+  scoringManager->SetScoreWriter(new IAEAScoreWriter());
+  // Initialize the default Hadrontherapy geometry
+  geometryController->SetGeometry("default");
+  // Initialize command based scoring
+  G4ScoringManager::GetScoringManager();
+  // Initialize the physics
+  runManager -> SetUserInitialization(new HadrontherapyPhysicsList());
+  // Initialize the primary particles
+  HadrontherapyPrimaryGeneratorAction *pPrimaryGenerator = new HadrontherapyPrimaryGeneratorAction();
+  runManager -> SetUserAction(pPrimaryGenerator);
+  // Optional UserActions: run, event, stepping
+  HadrontherapyRunAction* pRunAction = new HadrontherapyRunAction();
+  runManager -> SetUserAction(pRunAction);
+  HadrontherapyEventAction* pEventAction = new HadrontherapyEventAction();
+  runManager -> SetUserAction(pEventAction);
+  HadrontherapySteppingAction* steppingAction = new HadrontherapySteppingAction(pRunAction);
+  runManager -> SetUserAction(steppingAction);
+  // Interaction data: stopping powers
+  HadrontherapyInteractionParameters* pInteraction = new HadrontherapyInteractionParameters();
+                // Set the Random engine
+        CLHEP::HepRandom::setTheEngine(new CLHEP::RanecuEngine());
+        // Only if an initial random seed is needed
+        // G4int seed = time(0);
+        // CLHEP::HepRandom::setTheSeed(seed);
+        G4RunManager* runManager = new G4RunManager;
+        // Geometry controller is responsible for instantiating the
+        // geometries. All geometry specific setup tasks are now in class
+        // HadrontherapyGeometryController.
+        HadrontherapyGeometryController *geometryController = new HadrontherapyGeometryController();
+        // Connect the geometry controller to the G4 user interface
+        HadrontherapyGeometryMessenger *geometryMessenger = new HadrontherapyGeometryMessenger(geometryController);
+        G4ScoringManager *scoringManager = G4ScoringManager::GetScoringManager();
+        scoringManager->SetVerboseLevel(1);
+        scoringManager->SetScoreWriter(new IAEAScoreWriter());
+                // Initialize the default Hadrontherapy geometry
+        geometryController->SetGeometry("default");
+                // Initialize command based scoring
+        G4ScoringManager::GetScoringManager();
+                // Initialize the physics
+        runManager -> SetUserInitialization(new HadrontherapyPhysicsList());
+                // Initialize the primary particles
+        HadrontherapyPrimaryGeneratorAction *pPrimaryGenerator = new HadrontherapyPrimaryGeneratorAction();
+        runManager -> SetUserAction(pPrimaryGenerator);
+                // Optional UserActions: run, event, stepping
+        HadrontherapyRunAction* pRunAction = new HadrontherapyRunAction();
+        runManager -> SetUserAction(pRunAction);
+        HadrontherapyEventAction* pEventAction = new HadrontherapyEventAction();
+        runManager -> SetUserAction(pEventAction);
+        HadrontherapySteppingAction* steppingAction = new HadrontherapySteppingAction(pRunAction);
+        runManager -> SetUserAction(steppingAction);
+                // Interaction data: stopping powers
+        HadrontherapyInteractionParameters* pInteraction = new HadrontherapyInteractionParameters(true);
 #ifdef G4VIS_USE
   // Visualization manager
   G4VisManager* visManager = new G4VisExecutive;
   visManager -> Initialize();
+                // Visualization manager
+        G4VisManager* visManager = new G4VisExecutive;
+        visManager -> Initialize();
 #endif
+G4UImanager* UI = G4UImanager::GetUIpointer();
+ if (argc!=1)   // batch mode
+   {
+     G4String command = "/control/execute ";
+     G4String fileName = argv[1];
+     UI->ApplyCommand(command+fileName);
+   }
+ else  // interactive mode : define visualization UI terminal
+   {
+     G4UIsession* session = 0;
+     // If the enviroment variable for the TCSH terminal is active, it is used and the
+     // defaultMacro.mac file is executed
+        G4UImanager* UI = G4UImanager::GetUIpointer();
+        // Batch mode: in this configuration, i.e. when the User explicitly call a macro file,
+        // the program run and, at the end, authomatically call the exit state.
+        // With a such configuration, it is possible to send in background a simulation
+        // and correctly get the output file.
+        if (argc!=1)
+        {
+                G4UIsession* session = 0;
+                G4String command = "/control/execute ";
+                G4String fileName = argv[1];
+                UI -> ApplyCommand(command+fileName);
 #if defined(G4UI_USE_TCSH)
+     session = new G4UIterminal(new G4UItcsh);
+     UI->ApplyCommand("/control/execute defaultMacro.mac");
+     // Alternatively (if G4UI_USE_TCSH is not defined)  the program search for the
+     // G$UI_USE_QT variable. It starts a graphical user interface based on the QT libraries
+     // In the following case the GUI.mac file is also executed
+     //
+                session = new G4UIterminal(new G4UItcsh);
+#else
+                session = new G4UIterminal();
+#endif
+                // Only if the User Interface is needed!
+                //session -> SessionStart();
+                //delete session;
+        }
+        else  // interactive mode : define visualization UI terminal
+        {
+                G4UIsession* session = 0;
+                        // If the enviroment variable for the TCSH terminal is active, it is used and the
+                        // defaultMacro.mac file is executed
+#if defined(G4UI_USE_TCSH)
+                session = new G4UIterminal(new G4UItcsh);
+                UI->ApplyCommand("/control/execute defaultMacro.mac");
+                session -> SessionStart();
+                // Alternatively (if G4UI_USE_TCSH is not defined)  the program search for the
+                // G4UI_USE_QT variable. It starts a graphical user interface based on the QT libraries
+                // In the following case the GUI.mac file is also executed
 #elif defined(G4UI_USE_QT)
      session = new G4UIQt(argc,argv);
      UI->ApplyCommand("/control/execute macro/GUI.mac");
      // As final option, the simpler user interface terminal is opened
+                session = new G4UIQt(argc,argv);
+                UI->ApplyCommand("/control/execute macro/GUI.mac");
+                session -> SessionStart();
+                // As final option, the basic user interface terminal is opened
 #else
+    session = new G4UIterminal();
+    UI->ApplyCommand("/control/execute defaultMacro.mac");
+#endif
+    session->SessionStart();
+    delete session;
+   }
+    HadrontherapyMatrix* matrix = HadrontherapyMatrix::getInstance();
+    if (matrix) matrix -> TotalEnergyDeposit();
+#ifdef ANALYSIS_USE
+ analysis -> finish();
+#endif
+ // Job termination
+#ifdef G4VIS_USE
+ delete visManager;
+                session = new G4UIterminal();
+                UI -> ApplyCommand("/control/execute defaultMacro.mac");
+                session -> SessionStart();
+#endif
+                delete session;
+        }
+                // Job termination
+#ifdef G4ANALYSIS_USE_ROOT
+    HadrontherapyAnalysisManager::GetInstance() -> flush();     // Finalize the root file
+#endif
+#ifdef G4UI_USE
+      G4UIExecutive * ui = new G4UIExecutive(argc,argv);
+      ui->SessionStart();
+      delete ui;
 #endif
 …
   delete runManager;
   return 0;
+}

trunk/examples/advanced/hadrontherapy/History

-                      r1230
+                      r1313
              History file of the Hadrontherapy application
          ====================================================
+.11.2009 S.E.Mazzaglia & F.Romano; Tag:Hadrontherapy-V09-02-40
+.06.2010 J.Perl; Tag: Hadrontherapy-V09-03-04
+           - Updated vis usage
+.06.2010: G.A.P.Cirrone; Tag: harontherapy-V09-03-03
+           - Added G4UIExecutive
+           - Main file updated
+.05.2010 S.E.Mazzaglia; Tag: hadrontherapy-V09-03-02
+           - Updated ROOT scripts in folder RootScripts/proton/BraggPeak
+.05.2010 S.E.Mazzaglia; Tag: hadrontherapy-V09-03-01
+           - Solved some bugs related to the preprocessor variable G4ANALYSIS_USE_ROOT.
+.05.2010 S.E.Mazzaglia; Tag: hadrontherapy-V09-03-00
+           - Removed the possibility to calculate Let.
+           - Removed AIDA and variables: G4_ANALYSIS_USE and ANALYSIS_USE.
+           - Improved the commands to modify phantom/detector geometry
+           - Added a command to set cuts just only in the detector: /physic/setDetectorCuts
+           - Added the possibility to store dose & fluence for all secondaries
+             produced in the same ASCII file. This feature can be switched on by the macro
+             command: /analysis/secondaries true.
+.12.2009 S.E.Mazzaglia; Tag: hadrontherapy-V09-02-42
+           - Added the possibility to calculate the Let too, using an ASCII file to store it.
+           - Added the new class HadrontherapyLet.
+           - Minor revisions of the code.
+.11.2009 S.E.Mazzaglia; Tag: hadrontherapy-V09-02-41
+           - Added the possibility to store dose and fluence for every particle/ion
+             produced during the simulation
+.11.2009 S.E.Mazzaglia & F.Romano; Tag: hadrontherapy-V09-02-40
            - Corrected a bug in HadrontherapyDetectorConstruction class
            - Added G4RadiactiveDecayPhysics class to the Physics List.
 .11.2009 S.E.Mazzaglia; Tag: Hadrontherapy-V09-02-39
+           - Added G4RadiactivedecayPhysics class to the Physics List.
+.11.2009 S.E.Mazzaglia; Tag: hadrontherapy-V09-02-39
            - Correction in the initialization of the passiveProtonBeamLine class.

trunk/examples/advanced/hadrontherapy/README

-                      r807
+                      r1313
+     =========================================================
+                  Geant4 - Hadrontherapy example
+     =========================================================
+                                README file
+                          ----------------------
+                                   AUTHORS
+G.A.P. CIRRONE(a *), G.CUTTONE(a), F. DI ROSA(a), G.RUSSO(a)
+a. Laboratori Nazionali del Sud - Istituto Nazionale di Fisica Nucleare
+Catania, Italy
+* e-mail:cirrone@lns.infn.it
+M.G. PIA(b)
+b. Istituto Nazionale di Fisica Nucleare, Sezione di Genova Via Dodecaneso, 33
+, Genova, Italy
+A. LECHNER (c)
+c. CERN, Switzerland
+More informations on the Hadrontherapy example can be found in the
+Hadrontherapy Documentation available at http://www.ge.infn.it/geant4/examples/
+Alternatevely send an e-mail to cirrone@lns.infn.it.
+---->0. INTRODUCTION.
+The hadrontherapy example simulates a hadron therapy beam line.
+In particular the example models the specific proton therapy beam line
+installed at Laboratori Nazionali del Sud (INFN) in Catania, Sicily (Italy).
+For more information on the proton therapy center of Catania
+or/and proton/hadron therapy in general, please visit the
+pages:
+http://www.lns.infn.it/catanaweb/catana/
+---->1. GEOMETRY SET-UP.
+             =========================================================
+                     Text version of the Hadrontherapy README file
+             =========================================================
+Last revsion: G.A.P.Cirrone, November 2009;
+Released with the Geant4 9.3 version (December 2009)
+------------------------------------------------------------------------------------------------
+ADVERTISEMENT: this is the text version of the README file of the hadrontherapy application,
+as it has been released in the official Geant4 9.3 release
+A more complete and updated version of this file is published inside the web pages of Hadrontherapy:
+http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
+Please refer to the web-published version to know the last and future developments of hadrontherapy
+-------------------------------------------------------------------------------------------------
+             =========================================================
+                                HADRONTHERAPY
+             =========================================================
+ Code developed by:
+ G.A.P. Cirrone(a)°, G.Cuttone(a), F.Di Rosa(a), S.E.Mazzaglia(a), F.Romano(a)
+The elements simulated are:
+ Contributor authors:
+ P.Kaitaniemi(d), A.Heikkinen(d), G.Danielsen (d)
+ Past authors:
+ M.G.Pia(b), S.Guatelli(c), G.Russo(a), M.Russo(a), A.Lechner(e)
+ (a) Laboratori Nazionali del Sud
+     of the INFN, Catania, Italy
+ (b) INFN Section of Genova, Italy
+ (c) University of Wallongong, Australia
+ (d) Helsinki Institute of Physics, Helsinki, Finland
+ (e) CERN, (CH)
+  *Corresponding author, email to cirrone@lns.infn.it
+-------------------------------------------------------------------------------------------------
+HADRONTHERAPY:
+WHAT IT IS, WHAT IT DOES AND WHAT IT WILL PROVIDE
+Hadrontherapy is a Geant4-based application specifically developed to address typical needs related to the proton and ion therapy.
+It first release was in 2004. At that time Hadrontherapy was only capable to simulate a well specified proton therapy facility: the passive transport beam line installed at  Laboratori Nazionali del Sud (INFN) in Catania, Italy.
+Today Hadrontontherapy, except that it is in continuous development, is more flexible and show many additional capabilities as respect the past.
+Its geometrical set-up, for example, is now completely interchangeable permitting a simple switch between different geometrical configurations.
+In the actual version two geometrical configuration are available: the 'passive beam line'
+and the 'IAEA Benchmark' geometry. See the paragraph Geometry set-up for more information.
+Folder structure of Hadrontherapy
+Hadrontherapy distribution contain different sub-folders:
+\src: where source .cc files are stored
+\include: where header .hh files are stored
+\macro: where a set of ready-to-use macro files are provided
+\experimentalData: in this director a set of reference (both experimental and analithycal) data are stored. These data are then used to perform a direct comparison with simulation results that are stored in the simulationResults folder. Data stored are better described in the README file contained inside.
+\SimulationOutputs: when one of the .mac file contained in the macro folder is used, simulation results are directly stored in this directory.
+\RootScripts: if the ROOT program is installed the User can use the scripts contained in this directory to compare directly results from the his/her simulation with reference data provided inside the experimentalData folder.
+Currently this folders structure is in development and reference data as well as ROOT scripts will be added in the meanwhile new features and capabilities will be added. Moreover some ROOT script can be missed. Apologize for this and contact author if you need more information, clarification or useful discussion.
+Description of the \macro folder
+In the example directory, inside the "macro" folder two macro files are actually provided for the use of hadrontherapy with proton and carbon beams: proton_therapy.mac and ion_therapy.mac.
+The proton_therapy.mac permits to run a simulation with the whole passive beam line installed in Catania.
+The carbon_therapy.mac excludes all the elements (moving the origin of the ion beam close to the water phantom) and reproduce a simple passive beam line for the use with carbon beams.
+The macro iaea.mac uses an alternative geometry that was created for the IAEA benchmark. It features a geometry with water target, aluminum beam window, and a particle detector behind the target.
+DOWNLOAD AND INSTALLATION
+Hadrontherapy source code is actually released inside the official distribution of the Geant4 toolkit in the $G4INSTALL/examples/AdvancedExamples folder.
+To run Hadrontherapy you must first install the Geant4 package. Once Geant4 is installed the example must be first compiled (with the command gmake inside the
+../Hadrontherapy folder). When compilation is completed the program can be executed.
+A complete guide for the Geant4 installation in different operating systems can be found inside the official installation Geant4 pages.
+If you have troubles with the Geant4 installation please send an e-mail to us.
+SISTEM SET-UP: enviroment variables
+A standard Geant4 example GNUmakefile is provided
+The following section reports the environment variables that are necessary for the run of Hadrontherapy.
+#-------------------------------------
+#    SET UP LINUX  GCC
+#-------------------------------------
+VERSION="geant4-09-03"
+# Path to the directory in wich you have put data files and CLHEP
+LIBPATH=$HOME/Geant4Library
+export G4SYSTEM=Linux-g++
+# Path to the directory in wich you put your Geant4 installation
+export G4INSTALL=$HOME/${VERSION}
+export G4LIB=$G4INSTALL/lib
+export G4WORKDIR=$G4INSTALL/workdir
+export G4EXE=$G4WORKDIR/bin/$G4SYSTEM
+export CLHEP_BASE_DIR=$LIBPATH/CLHEP2.0.4.5
+export G4LEDATA=$LIBPATH/G4EMLOW6.9
+export G4LEVELGAMMADATA=$LIBPATH/PhotonEvaporation2.0
+export G4NEUTRONHPDATA=$LIBPATH/G4NDL3.13
+export G4RADIOACTIVEDATA=$LIBPATH/RadioactiveDecay3.2
+export G4ABLADATA=$LIBPATH/G4ABLA3.0
+export LD_LIBRARY_PATH=$CLHEP_BASE_DIR/lib:$LD_LIBRARY_PATH
+# For the generation .root file directly using the ROOT (if ROOT is
+# instaled in you machine)
+export G4ANALYSIS_USE_ROOT=1
+export LD_LIBRARY_PATH=$ROOTSYS/lib:$LD_LIBRARY_PATH
+#-------------------------------------
+#    SET UP    VRML VIEW
+#-------------------------------------
+export G4VIS_BUILD_VRML_DRIVER=1
+export G4VIS_USE_VRML=1
+export G4VIS_USE_VRMLFILE=1
+export G4VRMLFILE_MAX_FILE_NUM=100
+export G4VRMLFILE_VIEWER=vrmlview    #if installed
+# Add path to your VRML installation
+export PATH=$PATH:~/VRML
+#-------------------------------------
+#    SET UP    OpenGL o Mesa
+#-------------------------------------
+export G4VIS_BUILD_OPENGLX_DRIVER=1
+export G4VIS_USE_OPENGLX=1
+# Add path to your OpenGL installation
+#export OGLHOME=/usr/lib
+#-------------------------------------
+#    SET UP    DAWN (if installed)
+#-------------------------------------
+export G4VIS_BUILD_DAWN_DRIVER=1
+export G4VIS_BUILD_DOWNFILE_DRIVER=1
+export G4VIS_USE_DAWN=1
+export G4VIS_USE_DAWNFILE=1
+# Add path to your DAWN installation
+# export PATH=$PATH:~/dawn_3_86a
+# VARIOUS USER INTERFACES
+export G4UI_USE_XM=1
+export G4UI_USE_TCSH=1
+export G4UI_BUILD_QT_SESSION=1
+export G4UI_USE_QT=1
+# VARIOUS GRPHICAL USER INTERFACES
+export G4VIS_BUILD_QT_SESSION=1
+export G4VIS_BUILD_OPENGLQT_DRIVER=1
+export G4VIS_USE_OPENGLQT=1
+# If the QT libraries want be used for the User interfaces than the
+# correct path must be addressed
+export QTHOME=/usr/lib/qt-3.3
+export PATH=$PATH:/usr/lib/qt-3.3/include/
+export PATH=$PATH:/usr/lib/qt-3.3/
+GEOMETRICAL SET-UP
+The idea of Hadrontherapy is to provide a tool useful for Users interested in the field of proton and ion therapy. These can include the simple calculation of dose distribution curves in water or other materials, the derivation of important transport parameters (stopping powers, ranges, etc.) in different geometrical set-ups and for different materials, up to the complete simulation of a real transport beam line for therapy.
+The main component of the simulation is the phantom, a box that can be filled with different material and where the score of different information (at moment  the dose deposited and fluence in voxels) can be performed. A more complete description of the phantom is given in the next subsection.
+At moment the Hadrontherapy include the simulation of the proton beam line for eye-treatments installed at the INFN-LNS facility in Catania. This is a passive beam line and it is simulated in the file PassiveProtonBeamLine.cc.
+In the next future an ActiveProtonBeamLine.cc will be provided for the simulation of the active scanning treatment modality.
+Moreover the possibility to add a very simple set-up (a beam, a phantom where collect the informations and some simple component) will be also provided.
+All these configuration will be setted by macro commands.
+There is also a facility that allows the user to make a choice between alternative geometry set-ups. This can be done by using command:
+/geometrySetup/selectGeometry <name>
+where <name> is either "default" for the standard hadrontherapy geometry or "IAEA" for the IAEA benchmark geometry. Both geometries are described below. By default the standard "default" geometry is used.
+The water phantom to collect informations
+At the end of the beam line a phantom (a box of uniform material) is reproduced. Inside it, a user-defined region is divided (via the ROGeomtry classes of Geant4) in cubic and identical voxels. The voxels size can be varied as well as the voxelized region.
+At the end of a simulation run the dose deposited by primaries and secondaries in each voxel is collected. This information is available as an .out file or as a .root (if the G4ANALYSIS_USE variable is defined and the AIDA interface is activated).
+The default sizes of the active voxelized region are 40x40x40 mm and actually the default voxel configuration is 200 x 1 x 1, which means 200 slices with 0.2 mm of thickness.
+Of course this default can be modified in order to obtain, for example, a matrix of 80x80x80 cubic voxels each with a lateral dimension of 0.5 mm.
+As concern the cut and stepMax values, the default configuration implies a cut value of 0.01 mm in the whole  world (use the command /physic/setCuts <length>  in order to set the cut for all, and the command /physic/setDetectorCuts <length> to set the cut for the detector only)  and a stepMax of 0.01 mm just in the phantom (use the command /Step/waterPhantomStepMax 0.01 mm).
+In any case it is strongly recommended to use a stepMax value not bigger than 5% of the dose slice thickness.
+The Proton passive beam line class file
+The following is the description of the elements of the passive proton beam line of the Laboratori Nazionali del Sud in Catania (I). This line is completely simulated inside this class.
+The main elements are:
+* The SCATTERING SYSTEM: to transversally enlarge the original beam
+* The COLLIMATORS: placed along the beam line to collimate the beam;
+* The RANGE SHIFTERS: to decrease the energy of the primary proton beam to a specific value;
+* The MODULATOR WHEEL: to modulate the energy of the primary and mono-energetic beam in to a wide spectrum. The energy modulation is necessary to homogeneously irradiate a tumour volume that can extends in depth up to 20 mm;
+* The MONITOR CHAMBERS: very thin ionisation chamber that permit the dose monitoring during the patient irradiation;
+* The MOPI detector: microstrips, air free detector utilised for the check of the beam symmetry during the treatment;
+* The PATIENT COLLIMATOR: a brass, tumour-shaped collimator able to  confine the proton irradiation field in order to irradiate just the tumour mass in the transverse direction;
+The user has the possibility to vary, via messenger, almost all the geometrical characteristics of the beam line elements (i.e. their position along the beam line, their thickness, etc.).
+The elements simulated in the PassiveBeamLine.cc file are:
 . A scattering system, to spread geometrically the beam;
 …
 . A system of collimators, to avoid the scattering radiation;
+. A modulation system that spreads the beam in energy and
+   produces the so-called spread out bragg peak;
+   It is constituted by a rotating wheel of different thichnesses.
+   The wheel rotates around is axis (parallel to the proton
+   beam axis) and its movement can be obtained by means of a
+   messenger between runs.
+. A set of monitor chambers (special transmission ionisation
+   chambers used to control hadron flux during the irradiation);
+. A "nozzle" and a final collimator defining the final shape
+   of the beam before reaching the patient.
+. A water phantom: it is a box of water where the energy deposit is
+   calculated.
+   The use of  the water phantom is required by the international protocol
+   on the measure of dose in the case of proton and ion beams (IAEA 398, 2000).
+---->2. EXPERIMENTAL SET-UP.
+The application simulates the proton therapy beam line
+installed at Laboratori Nazionali del Sud.
+The default beam line is a typical treatment line composed by several elements all
+devoted to create the so-called "terapeutical beam", i.e. a beam ideal
+for a radiotherapeutic treatment.
+   The main elements are:
+** The COLLIMATORS: placed along the beam line to collimate the beam;
+** The RANGE SHIFTERS: to decrease the energy of the primary proton beam
+   to a specific value;
+** The MODULATOR WHEEL: to modulate the energy of the primary and monoenergetic
+   beam in to a wide spectrum. The energy modulation is necessary to
+   homogeneusly irradiate a tumour volume that can extends in depth
+   up to 20 mm;
+** The MONITOR CHAMBERS: very thin ionisation chamber that permit the
+   dose monitoring during the patient irradiation;
+** The PATIENT COLLIMATOR: a brass, tumour-shaped collimator able to
+   confine the proton irradiation field in order to irradiate just the tumour mass
+   in the trasverse direction;
+The user has the possibility to vary, via messenger, almost all the geometrical
+characteristics of the beam line elements (i.e. their position along the beam line,
+their thickness, etc.). More details on the available user messengers can be
+found in the Hadronterapy Documentation (http://www.ge.infn.it/geant4/examples/).
+At the end of the beam line, a typical water phantom is reproduced.
+A user-defined region of the phantom is divided (via the ROGeomtry class) in
+cubic and identical voxels. The voxels size can be varied. At the end of the simulation
+the energy deposited by primary protons, and secondaries in each voxel
+is collected. This information is available as an .hbk file (if the
+G4ANALYSIS_USE variable is defined).
+The default sizes of the active voxelized region are 40x40x40 mm corresponding
+to a matric of 80x80x80 cubic voxels each with a lateral dimension of 0.5 mm.
+---->3. SET-UP
+- a standard Geant4 example GNUmakefile is provided
+setup with:
+compiler = gcc-3.2.3
+G4SYSTEM = linux-g++
+The following section reports the necessary environment variables
+necessary for the run of Hadrontherapy.
+---->3.1  ENVIROMENT VARIABLES
+ - G4SYSTEM = Linux-g++
+ - G4INSTALL              points to the installation directory of GEANT4;
+ - G4LIB                  point to the compiled libraries of GEANT4;
+ - G4WORKDIR              points to the work directory;
+ - CLHEP_BASE_DIR         points to the installation directory of CHLEP;
+ - G4LEVELGAMMADATA       points to the photoevaporation library;
+ - NeutronHPCrossSections points to the neutron data files;
+ - G4RADIOACTIVEDATA      points to the libraries for radio-active decay
+                          hadronic processes;
+ - G4LEDATA               points to the low energy electromagnetic libraries
+ - LD_LIBRARY_PATH = $CLHEP_BASE_DIR/lib
+---->3.2  VISUALISATION
+The user can visualise the experimental set-up with OpenGL, DAWN and vrml
+---->4. HOW TO RUN THE EXAMPLE
+In interactive mode:
+. A modulation system that spreads the beam in energy and produces the so-called spread out Bragg peak; It is constituted by a rotating wheel of different thicknesses. The wheel  rotates around its axis (parallel to the proton beam axis) and its movement can be obtained by means of a messenger between runs.
+. A set of monitor chambers (special transmission ionisation chambers used to control the particle flux during the irradiation);
+. A final long collimator and a patient collimator defining the final shape of the beam before reaching the patient.
+. A water phantom: it is a box of water where the dose deposit is calculated. The use of  the water phantom is required by the international protocol on the measure of dose in the case of proton and ion beams (IAEA 398, 2000).
+Geometry for the IAEA benchmark
+Simple geometry for benchmark purposes contains water phantom (thickness can be set using a macro command), an aluminum beam window and Plexi-glass. Behind the phantom we have a detector that records outcoming particles. The IAEA geometry can be activated by using the command:
+/geometrySetup/selectGeometry IAEA
+Example IAEA benchmark run can be done as follows:
+run:        Hadrontherapy macro/iaea.mac
+analysis: root RootScripts/iaeaBenchmark/fragmentEnergy.C
+PHYSICS PROCESSES AND PHYSICS MODELS IMPLEMENTATION
+A particular care is addressed to the simulation of the physic processes.
+Three different approaches can be used for the choose of the physic models.
+Approach 1:
+Using the macro command:
+/physic/addPhysics/<physics List name>.
+In this case the models (for electromagnetic, hadronic elastic and hadronic inelastic) can be activated directly calling the name of the Physics Lists that are available inside the
+Geant4 kernel in the directory:
+$G4INSTALL/source/physics_lists/builders/include
+An example of the use of the Physics List can be found in the macro files:
+proton_therapy.mac and ion_therapy.mac
+Approach 2:
+A set of built-in physic models are also contained inside the Hadrontherapy directory. These are called Local*.cc and Local*.hh and can be activated using the macro command:
+/physic/addPhysics/<name>.
+Approach 3:
+We developed this approach in order to simplify the choice of the physic models to
+be used in the application.
+With this approach the user must only insert a command line in his/her .mac file using the: /physics/addPackage <PACKAGE_NAME>
+This permits to switch-on an already build physic package.
+Various packages are already present in the Geant4 tree: they are in the directory: geant4/source/physics_lists/lists/include
+In this case hadronic inelastic models are directly activated for every particle.
+NOTE: we do not recommend the use of local physics lists while we recommend the use of the Physics Lists or of the Reference Physics Lists (Approach 1 or 3)
+--------------------------------------------------------------------------------
+SUGGESTED PHYSIC FOR ION BEAMS IN THE RANGE 0 - 400 MeV
+--------------------------------------------------------------------------------
+Two macro files (proton_therapy.mac and ion_therapy.mac) can be used for proton and ion
+simulations.
+Also the QGSP_BIC_EMY package can be used if the *Approach 3* is preferred.
+INTERACTIVE COMMANDS
+How to change Phantom and Detector geometries
+In order to let the end user to change phantom and detector geometries and voxelization, some interactive commands have been provided. All parameters are mandatory, except those inside square brackets.
+Detector geometry
+The user can change:
+(1) The detector (box) size.
+(2) The voxels sizes. Changing this parameters, and/or the detector sizes, user should choose values in order to be divisors of the detector correspondent sizes.
+For both above commands, zero or negative values mean << don't change it >>
+(3) The displacement between the phantom and the detector.  Displacement parameters refer to the lower left corner of the detector respect to that of the phantom, by the point of view of the beam. In this case zero or positive values are allowed, while the negatives ones mean: << don't change it>>.
+Command synopsis: the <default values> follow "#"
+. /changeDetector/size <dimX> <dimY> <dimZ> <[unit]> # 4 4 4 cm
+. /changeDetector/voxelSize <dimX> <dimY> <dimZ> <[unit]> # 0.2 40 40 mm
+. /changeDetector/displacement <dispX> <dispY> <dispZ> <[unit]> # 0 18 18 cm
+where the X dimension is that along the beam direction
+Phantom geometry
+(1) The phantom size. As usually, zero or negatives values mean: <<don't change it>>.
+(2) The phantom position respect to the world. In this case specified values refer to the three components of the position of the phantom's center respect to the world's.
+Command synopsis:
+. /changePhantom/size <dimX> <dimY> <dimZ> <[unit]> # 40 40 40 cm
+. /changePhantom/position <posX> <posY> <posZ> <[unit]> # 20 0 0 cm
+All these commands must be followed by the command /changePhantom/update
+in order to check and apply them to the real geometry.
+Moreover must be issued between runs (so where you want but after the /run/initialize initialization command, or the G4State_Idle Geant4 state machine).
+Obviously all the previous sizes must be set in order to maintain the detector inside the phantom, otherwise system complains.
+ Some examples follow:
+/changeDetector/size 40 0 0 cm
+# Will extend detector X size to cover in full the phantom X size
+/changeDetector/size 0 4.5 0 cm
+# Will extend the Y size to 4.5 cm (note that voxel size Y is automatically
+#  rounded to 4.5 cm because the default value along Y is 4 cm)
+/changePhantom/update
+# Remember to always update the geometry before the beamOn command!!
+/changeDetector/size 0 8 0 cm
+# Will extend the Y size to 8 cm. In this case voxel size Y doesn't change, but
+# the number of voxel along Y doubles.
+/changePhantom/update
+/changeDetector/voxelSize 100 0 0 um
+# 100 um should be a divisor of detector size X
+# Will change only slabs X size to 100 um, without affecting the other.
+/changePhantom/update
+/changeDetector/displacement 0 0 0 # default unit mm
+# Will place the detector in the left lower corner (from the point of view of the beam) of #the  phantom.
+/changePhantom/update
+Stopping powers calculation
+It is possible for the end-user to calculate, via macro command, stopping powers only for those materials inserted into G4NistMaterialBuilder class (about 300).
+To get stopping powers user must provide this command line on the idle interactive terminal (or into a macro file) :
+/parameter/getstopping <G4_material> <Emin> <Emax> <nPoints> <[particle]> <[output_filename]>
+All parameters are mandatory except those inside square brackets [].
+Default values for parameters inside square brackets are respectively proton and standard output (usually the user console terminal).
+Parameters are respectively:
+* The material (NIST) name (something like G4_..., the complete list of elements and materials is available into the G4NistMaterialBuilder class and can be printed  to the terminal screen via the macro command: /parameter/nist )
+* Kinetic energy range in MeV and the number of data points to be retrieved (in a logarithmically uniform space)
+* The particle name (proton, e+, e-, He3, neutron,... a full list can be gotten via the macro command: /particle/list).
+         Only for ions, user must firstly give them to the particle gun, for example issuing the  macro commands:
+a. /gun/particle ion
+b. /gun/ion <Z> <A> <[charge]>
+* The output filename: if users leave this blank then the standard output is used.
+Below is an example in order to calculate the stopping power for alphas into Hydrogen between 1 keV to 150 MeV for 15 points:
+/parameter/getstopping G4_H 0.001 150 15 alpha
+# and for C12 ion:
+/gun/particle ion
+/gun/ion 6 12 6
+/parameter/getstopping G4_H 0.001 150 15 C12[0.0]
+# Value inside square brackets is the excitation energy of the ion (ground state).
+HOW RUN HADRONTHERAPY
+Run the example in interactive mode
 > $G4WORDIR/bin/Linux-g++/Hadrontherapy
+The defaultMacro.mac is executed
+The primary particle beam parameter are:
+Radiation:                proton beam;
+Energy distribution:      gaussian;
+Mean energy:              63.4 MeV;
+Energy spread:            300 keV;
+Beam spot size:           1 mm;
+Beam angular spread:      0.057 deg;
+The modulator wheel can be rotated via the messenger:
+Idle>/modulator/angle/xx deg
+To produce a Spread Out Bragg Peak using the modulator a macro
+(modulatorMacro.mac) is provided. With this macro the modulator is
+rotated of 360 degree at 1 deg steps. In each run 1000 protons are
+generated as primary particles. Obviously a bigger resolution can be obtained
+with smaller angles or increasing the protons number in each run.
+Modulator wheel can be omitted setting its material air.
+run $G4WORKDIR/bin/Linux-g++/Hadrontherapy visualisationMacro.mac
+to visualise the experimental set-up with OpenGL
+---->5. PHYSICS
+Both electromagnetic and hadronic physic processes are activated for
+the particles of the experimental set-up.
+Different physics models can be activated by the user interactively.
+Examples of activation are provided in the macro files starting
+with the string "physics":
+All possible physics options are summarized in the file
+physicsAllOptions.mac.
+Different options concerning electromagnetic interactions of protons
+and neutrons can be tested with the files:
+physicsElectromagneticICRU49.mac
+physicsElectromagneticZiegler77.mac
+physicsElectromagneticZiegler85.mac
+Different options concerning hadronic interactions of protons and
+neutrons can be tested with the files:
+physicsHadronicBertini.mac
+physicsHadronicBinary.mac
+physicsHadronicLEP.mac
+physicsHadronicPrecompound.mac
+NOTE: Apart from the different models for protons and neutrons, a user
+can also select among several interaction models for particles like
+electrons or photons. All possible options are listed in the file
+physicsAllOptions.mac.
+---->6. SIMULATION OUTPUT
+The output is an .hbk file (hadrontherapy.hbk) produced
+if the variable G4ANALYSIS_USE is set to 1 and the analysis tool (AIDA
+interface) correctly installed.
+The file contains an histogram and an n-tuple.
+The histogram contains the Bragg curve: energy deposited
+by the proton beam (in MeV) versus the depth in water (in mm).
+The n-tuple contains the total 3D energy deposit in the phantom; the information
+is energy deposit in each voxel with respect to the position of the voxel.
+Setup for analysis: AIDA 3.2.1
+Users can download the analysis tools from:
+http://aida.freehep.org/
+Note that the same information can be stored in any different format.
+Please contact cirrone@lns.infn.it if you want store the information in
+a different format.
+---------------------------------------------------------------------------
+for comments, advices, doubts and questions please contact:
+cirrone@lns.infn.it, giorgiorusso@lns.infn.it
+last modified: A. Lechner, 16/11/2007
+In this case the main file (Hadrontherapy.cc) performs different operations depending on which environment variable is activated;
+For example, if the environment variable G4UI_USE_TCSH is activated, Hadrontherapy will start with the TCSH User Interface that has many useful functionalities. On the other hand, if this first variables is not defined, the program will continue searching for the G4UI_USE_QT variable and, finally, will open the standard G4UITerminal.
+Run the example using macro files
+Hadrontherapy can be launched using a macro file:
+> $G4WORDIR/bin/Linux-g++/Hadrontherapy macroFile.mac
+The defaultMacro.mac file is contained in the main directory of Hadrontherapy and is automatically readed in case the user launch the executable without a parameter. A big number of other macro are inside the /macro folder.
+Example use of command based scoring
+In the IAEA geometry it is possible to collect the energy deposition data using the Geant4 command based scoring feature. This allows the users to define scorers interactively in the user interface without writing a single line of C++. Below is listed an example usage of command based scoring:
+/score/create/boxMesh boxMesh_1
+#Box size is the radius of the box ie 20x20x20 gives 40x40x40 outer dimensions
+/score/mesh/boxSize 13.95 16. 16. cm
+/score/mesh/translate/xyz 69 0.0 0.0 cm
+/score/mesh/nBin 400 1 1      # 400 bins in x-direction, 1 in y and z directions
+SIMULATION OUTPUT
+ASCII file
+A .out file is generated at the end of each run, Dose.out is its default name that can be changed in the HadrontherapyMatrix.cc file.
+The file can contain four columns; the first three columns represent the voxel indexes (that univocally identify the voxel volume), while the last column represent the dose deposited in the given voxel.
+Alternatively, if the macro command /analysis/secondaries true is given, ordinated dose and fluence, for every secondary produced, is added to the file.
+Setting the name of the ROOT output file
+By default the name of the ROOT output file is DoseDistribution.root. The name of the file can be set by using command:
+/analysis/setAnalysisFile <filename>
+It is also possible to create multiple new output files in the same simulation session. For example:
+/beam/energy/meanEnergy 4800 MeV
+/analysis/setAnalysisFile firstRun.root
+/run/beamOn 1000
+/beam/energy/meanEnergy 3000 MeV
+/analysis/setAnalysisFile secondRun.root
+/run/beamOn 1000
+Use of the ROOT analysis
+It is possible to use ROOT data analysis package directly for the production of output files.
+In the last version, anyway, this functionality must be implemented by User. This can be accomplished by setting an ad-hoc environment variable (i.e. G4ANALYSIS_USE_ROOT) to 1, adding in the code lines to create outputs with the ROOT libraries and recompiling the application.
+In this case you must have the ROOT framework installed in your machine.
+FUTURE CHALLENGES AND USERS' REQUESTS
+This is a list of future components that will be added in Hadrontherapy and of main Users requests that we hope to fulfill in the next future.
+What is in progress
+. A module for the simulation of an active beam line will be provided.
+The Korean Group of the Proton therapy center, National Cancer Center is developing this.
+. Modules for LET and RBE (Relative Biological Effectiveness) calculation. The Catania Group in Collaboration with the Turin one is working on this.
+. Hadrontherapy will permit the simulation of a 'basic' experiment usefull in hadrontherapy applications. User will be able to simulate thin/thick target configuration to calculate quantities like double differential cross sections of secondaries produces, or fluence and yelds of primary and secondaries in a thick block. This work is maily discussed with P. Kaitaniemi and the Group from the Helsinki Institute of Physics.
+What is requested from Users
+. Dicom interface
+Please contact cirrone@lns.infn.it for more details or suggestions and feedbacks on this document

trunk/examples/advanced/hadrontherapy/defaultMacro.mac

-                      r1230
+                      r1313
 # G.A.P.Cirrone
+ # G.A.P.Cirrone
+#
 # Default macro file. It is called if no argument is provided at run
 …
 # Set of the physic models
+#
 /physic/addPhysics emstandard_opt3                     # Electromagnetic model
+/physic/addPhysics standard_opt3                       # Electromagnetic model
 /physic/addPhysics elastic                             # Hadronic elastic model
 /physic/addPhysics binary                              # Hadronic inelastic model
 /physic/addPhysics local_ion_ion_inelastic             # Hadronic inelastic model for ions (local physic list)
 ##########################
 …
 # Visualisation
+#
+/vis/scene/create
+#/vis/open OGLIQt # only if QT library are installed
+/vis/open OGLIX
+/vis/open OGL
 /vis/viewer/flush
 /vis/viewer/set/viewpointThetaPhi 30 140 deg
 /vis/viewer/zoom 1
 /vis/viewer/pan -10  0 cm
+/tracking/storeTrajectory 1
+#/vis/scene/endOfEventAction accumulate
+/vis/scene/endOfEventAction accumulate -1
+/vis/viewer/update
+/vis/scene/add/trajectories
+/vis/scene/endOfEventAction accumulate
 ##########################
 …
 # Start of the run
+#
+/run/beamOn 10
+# If secondaries dose/fluence are needed
+#/analysis/secondaries true
+/run/beamOn 100

trunk/examples/advanced/hadrontherapy/include/CVS/Entries

-                      r1230
+                      r1313
 /Decay.hh/1.2/Wed Nov 18 14:37:51 2009//Tgeant4-09-03
 /HadrontherapyAnalysisFileMessenger.hh/1.4/Fri Jan  8 13:32:05 2010//Tgeant4-09-03
 /HadrontherapyAnalysisManager.hh/1.16/Wed Nov 18 14:37:51 2009//Tgeant4-09-03
 /HadrontherapyDetectorConstruction.hh/1.24/Fri Jan  8 13:32:05 2010//Tgeant4-09-03
 /HadrontherapyDetectorHit.hh/1.2/Wed Nov 18 14:37:51 2009//Tgeant4-09-03
 /HadrontherapyDetectorMessenger.hh/1.15/Fri Jan  8 13:32:05 2010//Tgeant4-09-03
 /HadrontherapyDetectorROGeometry.hh/1.3/Wed Nov 18 14:37:51 2009//Tgeant4-09-03
 /HadrontherapyDetectorSD.hh/1.2/Wed Nov 18 14:37:51 2009//Tgeant4-09-03
 /HadrontherapyDummySD.hh/1.6/Wed Nov 18 14:37:51 2009//Tgeant4-09-03
 /HadrontherapyEventAction.hh/1.15/Wed Nov 18 14:37:51 2009//Tgeant4-09-03
 /HadrontherapyEventActionMessenger.hh/1.2/Wed Nov 18 14:37:51 2009//Tgeant4-09-03
 /HadrontherapyGeometryController.hh/1.2/Wed Nov 18 14:37:51 2009//Tgeant4-09-03
 /HadrontherapyGeometryMessenger.hh/1.2/Wed Nov 18 14:37:51 2009//Tgeant4-09-03
 /HadrontherapyInteractionParameters.hh/1.5/Fri Jan  8 13:32:05 2010//Tgeant4-09-03
 /HadrontherapyMatrix.hh/1.6/Fri Jan  8 13:32:05 2010//Tgeant4-09-03
 /HadrontherapyModulator.hh/1.5/Thu Jun 29 15:58:08 2006//Tgeant4-09-03
 /HadrontherapyParameterMessenger.hh/1.5/Fri Jan  8 13:32:05 2010//Tgeant4-09-03
 /HadrontherapyParticles.hh/1.7/Wed Nov 18 14:37:51 2009//Tgeant4-09-03
 /HadrontherapyPhysicsList.hh/1.21/Fri Jan  8 13:32:05 2010//Tgeant4-09-03
 /HadrontherapyPhysicsListMessenger.hh/1.9/Wed Nov 18 14:37:51 2009//Tgeant4-09-03
 /HadrontherapyPrimaryGeneratorAction.hh/1.17/Wed Nov 18 14:37:51 2009//Tgeant4-09-03
 /HadrontherapyPrimaryGeneratorMessenger.hh/1.7/Wed Nov 18 14:37:51 2009//Tgeant4-09-03
 /HadrontherapyRunAction.hh/1.14/Wed Nov 18 14:37:51 2009//Tgeant4-09-03
 /HadrontherapyStepMax.hh/1.2/Wed Nov 18 14:37:51 2009//Tgeant4-09-03
 /HadrontherapyStepMaxMessenger.hh/1.3/Wed Nov 18 14:37:51 2009//Tgeant4-09-03
 /HadrontherapySteppingAction.hh/1.16/Wed Nov 18 14:37:51 2009//Tgeant4-09-03
 /IAEADetectorConstruction.hh/1.5/Wed Nov 18 14:37:51 2009//Tgeant4-09-03
 /IAEADetectorMessenger.hh/1.2/Wed Nov 18 14:37:51 2009//Tgeant4-09-03
 /IAEAScoreWriter.hh/1.3/Wed Nov 18 14:37:51 2009//Tgeant4-09-03
 /LocalINCLIonIonInelasticPhysic.hh/1.2/Wed Nov 18 14:37:51 2009//Tgeant4-09-03
 /LocalIonIonInelasticPhysic.hh/1.2/Wed Nov 18 14:37:51 2009//Tgeant4-09-03
 /PassiveProtonBeamLine.hh/1.5/Wed Nov 18 14:37:51 2009//Tgeant4-09-03
 /PassiveProtonBeamLineMessenger.hh/1.2/Wed Nov 18 14:37:51 2009//Tgeant4-09-03
+/Decay.hh/1.2/Wed Nov 18 14:37:51 2009//Tgeant4-09-04-beta-cand-00
+/HadrontherapyAnalysisFileMessenger.hh/1.6/Mon Jun 14 13:47:41 2010//Tgeant4-09-04-beta-cand-00
+/HadrontherapyAnalysisManager.hh/1.20/Mon Jun 14 13:47:41 2010//Tgeant4-09-04-beta-cand-00
+/HadrontherapyDetectorConstruction.hh/1.28/Mon Jun 14 13:47:41 2010//Tgeant4-09-04-beta-cand-00
+/HadrontherapyDetectorHit.hh/1.2/Wed Nov 18 14:37:51 2009//Tgeant4-09-04-beta-cand-00
+/HadrontherapyDetectorMessenger.hh/1.18/Mon Jun 14 13:47:41 2010//Tgeant4-09-04-beta-cand-00
+/HadrontherapyDetectorROGeometry.hh/1.3/Wed Nov 18 14:37:51 2009//Tgeant4-09-04-beta-cand-00
+/HadrontherapyDetectorSD.hh/1.3/Mon Jun 14 13:47:41 2010//Tgeant4-09-04-beta-cand-00
+/HadrontherapyDummySD.hh/1.6/Wed Nov 18 14:37:51 2009//Tgeant4-09-04-beta-cand-00
+/HadrontherapyEventAction.hh/1.16/Mon Jun 14 13:47:41 2010//Tgeant4-09-04-beta-cand-00
+/HadrontherapyEventActionMessenger.hh/1.2/Wed Nov 18 14:37:51 2009//Tgeant4-09-04-beta-cand-00
+/HadrontherapyGeometryController.hh/1.2/Wed Nov 18 14:37:51 2009//Tgeant4-09-04-beta-cand-00
+/HadrontherapyGeometryMessenger.hh/1.2/Wed Nov 18 14:37:51 2009//Tgeant4-09-04-beta-cand-00
+/HadrontherapyInteractionParameters.hh/1.8/Mon Jun 14 13:47:41 2010//Tgeant4-09-04-beta-cand-00
+/HadrontherapyMatrix.hh/1.12/Mon Jun 14 13:47:41 2010//Tgeant4-09-04-beta-cand-00
+/HadrontherapyModulator.hh/1.5/Thu Jun 29 15:58:08 2006//Tgeant4-09-04-beta-cand-00
+/HadrontherapyParameterMessenger.hh/1.5/Fri Jan  8 13:32:05 2010//Tgeant4-09-04-beta-cand-00
+/HadrontherapyParticles.hh/1.7/Wed Nov 18 14:37:51 2009//Tgeant4-09-04-beta-cand-00
+/HadrontherapyPhysicsList.hh/1.22/Mon Jun 14 13:47:41 2010//Tgeant4-09-04-beta-cand-00
+/HadrontherapyPhysicsListMessenger.hh/1.10/Mon Jun 14 13:47:41 2010//Tgeant4-09-04-beta-cand-00
+/HadrontherapyPrimaryGeneratorAction.hh/1.18/Mon Jun 14 13:47:41 2010//Tgeant4-09-04-beta-cand-00
+/HadrontherapyPrimaryGeneratorMessenger.hh/1.7/Wed Nov 18 14:37:51 2009//Tgeant4-09-04-beta-cand-00
+/HadrontherapyRunAction.hh/1.16/Mon Jun 14 13:47:41 2010//Tgeant4-09-04-beta-cand-00
+/HadrontherapyStepMax.hh/1.2/Wed Nov 18 14:37:51 2009//Tgeant4-09-04-beta-cand-00
+/HadrontherapyStepMaxMessenger.hh/1.3/Wed Nov 18 14:37:51 2009//Tgeant4-09-04-beta-cand-00
+/HadrontherapySteppingAction.hh/1.16/Wed Nov 18 14:37:51 2009//Tgeant4-09-04-beta-cand-00
+/IAEADetectorConstruction.hh/1.5/Wed Nov 18 14:37:51 2009//Tgeant4-09-04-beta-cand-00
+/IAEADetectorMessenger.hh/1.2/Wed Nov 18 14:37:51 2009//Tgeant4-09-04-beta-cand-00
+/IAEAScoreWriter.hh/1.3/Wed Nov 18 14:37:51 2009//Tgeant4-09-04-beta-cand-00
+/LocalINCLIonIonInelasticPhysic.hh/1.2/Wed Nov 18 14:37:51 2009//Tgeant4-09-04-beta-cand-00
+/LocalIonIonInelasticPhysic.hh/1.2/Wed Nov 18 14:37:51 2009//Tgeant4-09-04-beta-cand-00
+/PassiveProtonBeamLine.hh/1.6/Mon Jun 14 13:47:41 2010//Tgeant4-09-04-beta-cand-00
+/PassiveProtonBeamLineMessenger.hh/1.3/Mon Jun 14 13:47:41 2010//Tgeant4-09-04-beta-cand-00
+D

trunk/examples/advanced/hadrontherapy/include/CVS/Tag

r1230	r1313
1		Ngeant4-09-03
	1	Ngeant4-09-04-beta-cand-00

trunk/examples/advanced/hadrontherapy/include/HadrontherapyAnalysisFileMessenger.hh

-                      r1230
+                      r1313
 #define HadrontherapyAnalysisFileMessenger_h 1
-#ifdef ANALYSIS_USE
 #include "G4UImessenger.hh"
 …
 class HadrontherapyAnalysisManager; ///< Provides SetanalysisFileName()
+class G4UIcmdWithAString; ///< Provides possibility to add commands
+class G4UIcmdWithAString;
+class G4UIcmdWithABool;
 /**
 …
         * Constructor requires command name and messenger class(this).
         */
+    G4UIcmdWithAString* FileNameCmd;
+    G4UIcmdWithABool *secondariesCmd;
+#ifdef G4ANALYSIS_USE_ROOT
+    G4UIcmdWithAString *FileNameCmd;
+#endif
 };
 #endif
-#endif

trunk/examples/advanced/hadrontherapy/include/HadrontherapyAnalysisManager.hh

-                      r1230
+                      r1313
 //
 // ********************************************************************
 // * License and Disclaimer                                           *
 // *                                                                  *
 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
 // * conditions of the Geant4 Software License,  included in the file *
 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
 // * include a list of copyright holders.                             *
 // *                                                                  *
 // * Neither the authors of this software system, nor their employing *
 // * institutes,nor the agencies providing financial support for this *
 // * work  make  any representation or  warranty, express or implied, *
 // * regarding  this  software system or assume any liability for its *
 // * use.  Please see the license in the file  LICENSE  and URL above *
 // * for the full disclaimer and the limitation of liability.         *
 // *                                                                  *
 // * This  code  implementation is the result of  the  scientific and *
 // * technical work of the GEANT4 collaboration.                      *
 // * By using,  copying,  modifying or  distributing the software (or *
 // * any work based  on the software)  you  agree  to acknowledge its *
 // * use  in  resulting  scientific  publications,  and indicate your *
 // * acceptance of all terms of the Geant4 Software license.          *
 // ********************************************************************
 //
 // HadrontherapyAnalysisManager.hh; May 2005
 // See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
+        //
+        // ********************************************************************
+        // * License and Disclaimer                                           *
+        // *                                                                  *
+        // * The  Geant4 software  is  copyright of the Copyright Holders  of *
+        // * the Geant4 Collaboration.  It is provided  under  the terms  and *
+        // * conditions of the Geant4 Software License,  included in the file *
+        // * LICENSE and available at  http://cern.ch/geant4/license .  These *
+        // * include a list of copyright holders.                             *
+        // *                                                                  *
+        // * Neither the authors of this software system, nor their employing *
+        // * institutes,nor the agencies providing financial support for this *
+        // * work  make  any representation or  warranty, express or implied, *
+        // * regarding  this  software system or assume any liability for its *
+        // * use.  Please see the license in the file  LICENSE  and URL above *
+        // * for the full disclaimer and the limitation of liability.         *
+        // *                                                                  *
+        // * This  code  implementation is the result of  the  scientific and *
+        // * technical work of the GEANT4 collaboration.                      *
+        // * By using,  copying,  modifying or  distributing the software (or *
+        // * any work based  on the software)  you  agree  to acknowledge its *
+        // * use  in  resulting  scientific  publications,  and indicate your *
+        // * acceptance of all terms of the Geant4 Software license.          *
+        // ********************************************************************
+        //
+        // HadrontherapyAnalysisManager.hh; May 2005
+        // See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
 #ifndef HADRONTHERAPYANALYSISMANAGER_HH
 …
 #include "globals.hh"
-#ifdef ANALYSIS_USE ///< If we use analysis
-#ifdef G4ANALYSIS_USE ///< If analysis is done via AIDA
-#include <AIDA/AIDA.h>
-namespace AIDA{
-  class ITree;
-  class IAnalysisFactory;
-  class ITreeFactory;
+}
-#endif
 #ifdef G4ANALYSIS_USE_ROOT ///< If analysis is done directly with ROOT
 …
 #include "TH1F.h"
 #endif
 /**
  * Messenger class for analysis-settings for HadronTherapyAnalysisManager
 …
         /**
          * Analysis manager is a singleton object (there is only one instance).
          * The pointer to this object is available through the use of the method getInstance();
+         * The pointer to this object is available through the use of the method GetInstance();
+         *
          * @see getInstance
          */
   HadrontherapyAnalysisManager();
+         * @see GetInstance
+         */
+        HadrontherapyAnalysisManager();
 public:
+  ~HadrontherapyAnalysisManager();
+  /**
+   * Get the pointer to the analysis manager.
+   */
+  static HadrontherapyAnalysisManager* getInstance();
+  /**
+  * Book the histograms and ntuples in an AIDA or ROOT file.
+  */
+  void book();
+  /**
+   * Set name for the analysis file .root (used by macro)
+   */
+  void SetAnalysisFileName(G4String);
+  /**
+  * Fill the ntuple with the energy deposit in the phantom
+  */
+  void FillEnergyDeposit(G4int voxelXId, G4int voxelYId, G4int voxelZId,
+                         G4double energyDeposit);
+  void BraggPeak(G4int, G4double); ///< Fill 1D histogram with the Bragg peak in the phantom
+  void SecondaryProtonEnergyDeposit(G4int slice, G4double energy);
+  ///< Fill 1D histogram with the energy deposit of secondary protons
+   void SecondaryNeutronEnergyDeposit(G4int slice, G4double energy);
+  ///< Fill 1D histogram with the energy deposit of secondary neutrons
+  void SecondaryAlphaEnergyDeposit(G4int slice, G4double energy);
+  ///< Fill 1D histogram with the energy deposit of secondary alpha particles
+  void SecondaryGammaEnergyDeposit(G4int slice, G4double energy);
+  ///< Fill 1D histogram with the energy deposit of secondary gamma
+  void SecondaryElectronEnergyDeposit(G4int slice, G4double energy);
+  ///< Fill 1D histogram with the energy deposit of secondary electrons
+  void SecondaryTritonEnergyDeposit(G4int slice, G4double energy);
+  ///< Fill 1D histogram with the energy deposit of secondary tritons
+  void SecondaryDeuteronEnergyDeposit(G4int slice, G4double energy);
+  ///< Fill 1D histogram with the energy deposit of secondary deuterons
+  void SecondaryPionEnergyDeposit(G4int slice, G4double energy);
+  ///< Fill 1D histogram with the energy deposit of secondary pions
+  void electronEnergyDistribution(G4double secondaryParticleKineticEnergy);
+  ///< Energy distribution of secondary electrons originated in the phantom
+  void gammaEnergyDistribution(G4double secondaryParticleKineticEnergy);
+  ///< Energy distribution of secondary gamma originated in the phantom
+  void deuteronEnergyDistribution(G4double secondaryParticleKineticEnergy);
+  ///< Energy distribution of secondary deuterons originated in the phantom
+  void tritonEnergyDistribution(G4double secondaryParticleKineticEnergy);
+  ///< Energy distribution of secondary tritons originated in the phantom
+  void alphaEnergyDistribution(G4double secondaryParticleKineticEnergy);
+  ///< Energy distribution of secondary alpha originated in the phantom
+  void heliumEnergy(G4double secondaryParticleKineticEnergy);
+  ///< Energy distribution of the helium (He3 and alpha) particles after the phantom
+  void hydrogenEnergy(G4double secondaryParticleKineticEnergy);
+  ///< Energy distribution of the hydrogen (proton, d, t) particles after the phantom
+  void fillFragmentTuple(G4int A, G4double Z, G4double energy, G4double posX, G4double posY, G4double posZ);
+  ///< Energy ntuple
+  void genericIonInformation(G4int, G4double, G4int, G4double);
+  void ThintargetBeamDisp(G4double,G4double);
+  void startNewEvent();
+  ///< Tell the analysis manager that a new event is starting
+  void setGeometryMetaData(G4double, G4double, G4double);
+  ///< from the detector construction information about the geometry can be written as metadata
+  void setBeamMetaData(G4double, G4double);
+  ///< metadata about the beam can be written this way
+  void finish();
+  ///< Close the .hbk file with the histograms and the ntuples
+void flush();
+#ifdef G4ANALYSIS_USE_ROOT
+        ~HadrontherapyAnalysisManager();
+        /**
+         * Get the pointer to the analysis manager.
+         */
+        static HadrontherapyAnalysisManager* GetInstance();
+#ifdef G4ANALYSIS_USE_ROOT
+         /**
+         * Clear analysis manager heap.
+         */
+        void Clear();
+        /**
+         * Book the histograms and ntuples in an AIDA or ROOT file.
+         */
+        void book();
+        /**
+         * Set name for the analysis file .root (used by macro)
+         */
+        void SetAnalysisFileName(G4String);
+        /**
+         * Fill the ntuple with the energy deposit in the phantom
+         */
+        void FillEnergyDeposit(G4int voxelXId, G4int voxelYId, G4int voxelZId,
+                                                   G4double energyDeposit);
+        void BraggPeak(G4int, G4double); ///< Fill 1D histogram with the Bragg peak in the phantom
+        void SecondaryProtonEnergyDeposit(G4int slice, G4double energy);
+                ///< Fill 1D histogram with the energy deposit of secondary protons
+        void SecondaryNeutronEnergyDeposit(G4int slice, G4double energy);
+                ///< Fill 1D histogram with the energy deposit of secondary neutrons
+        void SecondaryAlphaEnergyDeposit(G4int slice, G4double energy);
+                ///< Fill 1D histogram with the energy deposit of secondary alpha particles
+        void SecondaryGammaEnergyDeposit(G4int slice, G4double energy);
+                ///< Fill 1D histogram with the energy deposit of secondary gamma
+        void SecondaryElectronEnergyDeposit(G4int slice, G4double energy);
+                ///< Fill 1D histogram with the energy deposit of secondary electrons
+        void SecondaryTritonEnergyDeposit(G4int slice, G4double energy);
+                ///< Fill 1D histogram with the energy deposit of secondary tritons
+        void SecondaryDeuteronEnergyDeposit(G4int slice, G4double energy);
+                ///< Fill 1D histogram with the energy deposit of secondary deuterons
+        void SecondaryPionEnergyDeposit(G4int slice, G4double energy);
+                ///< Fill 1D histogram with the energy deposit of secondary pions
+        void electronEnergyDistribution(G4double secondaryParticleKineticEnergy);
+                ///< Energy distribution of secondary electrons originated in the phantom
+        void gammaEnergyDistribution(G4double secondaryParticleKineticEnergy);
+                ///< Energy distribution of secondary gamma originated in the phantom
+        void deuteronEnergyDistribution(G4double secondaryParticleKineticEnergy);
+                ///< Energy distribution of secondary deuterons originated in the phantom
+        void tritonEnergyDistribution(G4double secondaryParticleKineticEnergy);
+                ///< Energy distribution of secondary tritons originated in the phantom
+        void alphaEnergyDistribution(G4double secondaryParticleKineticEnergy);
+                ///< Energy distribution of secondary alpha originated in the phantom
+        void heliumEnergy(G4double secondaryParticleKineticEnergy);
+                ///< Energy distribution of the helium (He3 and alpha) particles after the phantom
+        void hydrogenEnergy(G4double secondaryParticleKineticEnergy);
+                ///< Energy distribution of the hydrogen (proton, d, t) particles after the phantom
+                //Kinetic energy by voxel, mass number A and atomic number Z.
+        void FillKineticFragmentTuple(G4int i, G4int j, G4int k, G4int A, G4double Z, G4double kinEnergy);
+                //Kinetic energy by voxel, mass number A and atomic number Z of only primary particles
+        void FillKineticEnergyPrimaryNTuple(G4int i, G4int j, G4int k, G4double kinEnergy);
+                ///< Energy by voxel, mass number A and atomic number Z.
+        void FillVoxelFragmentTuple(G4int i, G4int j, G4int k, G4int A, G4double Z, G4double energy, G4double fluence);
+        void FillFragmentTuple(G4int A, G4double Z, G4double energy, G4double posX, G4double posY, G4double posZ);
+                ///< Energy ntuple
+        void FillLetFragmentTuple(G4int i, G4int j, G4int k, G4int A, G4double Z, G4double letT, G4double letD);
+                ///< let ntuple
+        void genericIonInformation(G4int, G4double, G4int, G4double);
+        void ThintargetBeamDisp(G4double,G4double);
+        void startNewEvent();
+                ///< Tell the analysis manager that a new event is starting
+        void setGeometryMetaData(G4double, G4double, G4double);
+                ///< from the detector construction information about the geometry can be written as metadata
+        void setBeamMetaData(G4double, G4double);
+                ///< metadata about the beam can be written this way
+        void flush();
+                ///< Close the .hbk file with the histograms and the ntuples
 private:
+  TH1F *createHistogram1D(const TString name, const TString title, int bins, double xmin, double xmax) {
+    TH1F *histo = new TH1F(name, title, bins, xmin, xmax);
+    histo->SetLineWidth(2);
+    return histo;
+  }
+#endif
+        TH1F *createHistogram1D(const TString name, const TString title, int bins, double xmin, double xmax) {
+                TH1F *histo = new TH1F(name, title, bins, xmin, xmax);
+                histo->SetLineWidth(2);
+                return histo;
+        }
 private:
+  static HadrontherapyAnalysisManager* instance;
+  HadrontherapyAnalysisFileMessenger* fMess;
+  G4String analysisFileName;
+#ifdef G4ANALYSIS_USE
+  AIDA::IAnalysisFactory* aFact;
+  AIDA::ITree* theTree;
+  AIDA::IHistogramFactory *histFact;
+  AIDA::ITupleFactory *tupFact;
+  AIDA::IHistogram1D *h1;
+  AIDA::IHistogram1D *h2;
+  AIDA::IHistogram1D *h3;
+  AIDA::IHistogram1D *h4;
+  AIDA::IHistogram1D *h5;
+  AIDA::IHistogram1D *h6;
+  AIDA::IHistogram1D *h7;
+  AIDA::IHistogram1D *h8;
+  AIDA::IHistogram1D *h9;
+  AIDA::IHistogram1D *h10;
+  AIDA::IHistogram1D *h11;
+  AIDA::IHistogram1D *h12;
+  AIDA::IHistogram1D *h13;
+  AIDA::IHistogram1D *h14;
+  AIDA::IHistogram1D *h15;
+  AIDA::IHistogram1D *h16;
+  AIDA::ITuple *ntuple;
+  AIDA::ITuple *ionTuple;
+  AIDA::ITuple *fragmentTuple;
+#endif
+#ifdef G4ANALYSIS_USE_ROOT
+  TFile *theTFile;
+  TH1F *histo1;
+  TH1F *histo2;
+  TH1F *histo3;
+  TH1F *histo4;
+  TH1F *histo5;
+  TH1F *histo6;
+  TH1F *histo7;
+  TH1F *histo8;
+  TH1F *histo9;
+  TH1F *histo10;
+  TH1F *histo11;
+  TH1F *histo12;
+  TH1F *histo13;
+  TH1F *histo14;
+  TH1F *histo15;
+  TH1F *histo16;
+  TNtuple *theROOTNtuple;
+  TNtuple *theROOTIonTuple;
+  TNtuple *fragmentNtuple; // fragments
+  TNtuple *metaData;
+#endif
+  G4long eventCounter;      // Simulation metadata
+  G4double detectorDistance;
+  G4double phantomDepth;
+  G4double beamEnergy;
+  G4double energyError;
+  G4double phantomCenterDistance;
+#endif
+        static HadrontherapyAnalysisManager* instance;
+        HadrontherapyAnalysisFileMessenger* fMess;
+#ifdef G4ANALYSIS_USE_ROOT
+        G4String analysisFileName;
+        TFile *theTFile;
+        TH1F *histo1;
+        TH1F *histo2;
+        TH1F *histo3;
+        TH1F *histo4;
+        TH1F *histo5;
+        TH1F *histo6;
+        TH1F *histo7;
+        TH1F *histo8;
+        TH1F *histo9;
+        TH1F *histo10;
+        TH1F *histo11;
+        TH1F *histo12;
+        TH1F *histo13;
+        TH1F *histo14;
+        TH1F *histo15;
+        TH1F *histo16;
+        TNtuple *kinFragNtuple;
+        TNtuple *kineticEnergyPrimaryNtuple;
+        // ntuple containing the fluence of all the particle in any voxel
+        TNtuple *doseFragNtuple;
+        // ntuple containing the fluence of all the particle in any voxel
+        TNtuple *fluenceFragNtuple;
+        // ntuple containing the fluence of all the particle in any voxel
+        TNtuple *letFragNtuple;
+        TNtuple *theROOTNtuple;
+        TNtuple *theROOTIonTuple;
+        TNtuple *fragmentNtuple; // fragments
+        TNtuple *metaData;
+        G4long eventCounter;      // Simulation metadata
+        G4double detectorDistance;
+        G4double phantomDepth;
+        G4double beamEnergy;
+        G4double energyError;
+        G4double phantomCenterDistance;
+#endif
 };
 #endif
+#endif

trunk/examples/advanced/hadrontherapy/include/HadrontherapyDetectorConstruction.hh

-                      r1230
+                      r1313
   void ConstructPhantom();
   void ConstructDetector();
+  void ConstructSensitiveDetector(G4ThreeVector position_respect_to_WORLD);
+  void ConstructSensitiveDetector(G4ThreeVector positionToWORLD);
+  void ParametersCheck();
 public:
 // Get detector position relative to WORLD
 …
+  }
 /////////////////////////////////////////////////////////////////////////////
 // Get displacement between phantom and detector by detector position, phantom and detector sizes
+// Get displacement between phantom and detector by detector position (center of), phantom (center of) and detector sizes
 inline G4ThreeVector GetDetectorToPhantomPosition()
+{
     return G4ThreeVector(phantomSizeX - detectorSizeX + detectorPosition.getX(),
                          phantomSizeY - detectorSizeY + detectorPosition.getY(),
                          phantomSizeZ - detectorSizeZ + detectorPosition.getZ()
+    return G4ThreeVector(phantomSizeX/2 - detectorSizeX/2 + detectorPosition.getX(),
+                         phantomSizeY/2 - detectorSizeY/2 + detectorPosition.getY(),
+                         phantomSizeZ/2 - detectorSizeZ/2 + detectorPosition.getZ()
                           );
+}
 …
+  {
           // Adjust detector position
           detectorPosition.setX(detectorToPhantomPosition.getX() - phantomSizeX + detectorSizeX);
           detectorPosition.setY(detectorToPhantomPosition.getY() - phantomSizeY + detectorSizeY);
           detectorPosition.setZ(detectorToPhantomPosition.getZ() - phantomSizeZ + detectorSizeZ);
+          detectorPosition.setX(detectorToPhantomPosition.getX() - phantomSizeX/2 + detectorSizeX/2);
+          detectorPosition.setY(detectorToPhantomPosition.getY() - phantomSizeY/2 + detectorSizeY/2);
+          detectorPosition.setZ(detectorToPhantomPosition.getZ() - phantomSizeZ/2 + detectorSizeZ/2);
+      if (detectorPhysicalVolume) detectorPhysicalVolume -> SetTranslation(detectorPosition);
+    //G4cout << "*************** DetectorToPhantomPosition " << detectorToPhantomPosition/cm << "\n";
+    //G4cout << "*************** DetectorPosition " << detectorPosition/cm << "\n";
+  }
 /////////////////////////////////////////////////////////////////////////////
 // Check whether detector is inside phantom
 inline bool IsInside(G4double detectorHalfX,
                      G4double detectorHalfY,
                      G4double detectorHalfZ,
                      G4double phantomHalfX,
                      G4double phantomHalfY,
                      G4double phantomHalfZ,
+inline bool IsInside(G4double detectorX,
+                     G4double detectorY,
+                     G4double detectorZ,
+                     G4double phantomX,
+                     G4double phantomY,
+                     G4double phantomZ,
                      G4ThreeVector detectorToPhantomPosition)
+{
 // Dimensions check... X Y and Z
 // Firstly check what dimension we are modifying
-        if (detectorHalfX > 0. && phantomHalfX > 0. && detectorToPhantomPosition.getX() >=0.)
+        {
             if (detectorHalfX > phantomHalfX)
+            if (detectorX > phantomX)
+                 {
                     G4cout << "Error: Detector X dimension must be smaller or equal to the corrispondent of the phantom" << G4endl;
                     return false;
+                 }
             if ( 2*(phantomHalfX - detectorHalfX) < detectorToPhantomPosition.getX())
+            if ( (phantomX - detectorX) < detectorToPhantomPosition.getX())
+                 {
                     G4cout << "Error: X dimension doesn't fit with detector to phantom relative position" << G4endl;
 …
+        }
-        if (detectorHalfY > 0. && phantomHalfY > 0.&& detectorToPhantomPosition.getY() >=0.)
+        {
             if (detectorHalfY > phantomHalfY)
+            if (detectorY > phantomY)
+                 {
                     G4cout << "Error: Detector Y dimension must be smaller or equal to the corrispondent of the phantom" << G4endl;
                     return false;
+                 }
             if ( 2*(phantomHalfY - detectorHalfY) < detectorToPhantomPosition.getY())
+            if ( (phantomY - detectorY) < detectorToPhantomPosition.getY())
+             {
                    G4cout << "Error: Y dimension doesn't fit with detector to phantom relative position" << G4endl;
 …
+        }
-        if (detectorHalfZ > 0. && phantomHalfZ > 0.&& detectorToPhantomPosition.getZ() >=0.)
+        {
             if (detectorHalfZ > phantomHalfZ)
+            if (detectorZ > phantomZ)
+                 {
                     G4cout << "Error: Detector Z dimension must be smaller or equal to the corrispondent of the phantom" << G4endl;
                     return false;
+                 }
             if ( 2*(phantomHalfZ - detectorHalfZ) < detectorToPhantomPosition.getZ())
+            if ( (phantomZ - detectorZ) < detectorToPhantomPosition.getZ())
+             {
                    G4cout << "Error: Z dimension doesn't fit with detector to phantom relative position" << G4endl;
 …
+             }
+        }
+/*
+    G4cout << "Displacement between Phantom and Detector is: ";
+    G4cout << "DX= "<< G4BestUnit(detectorToPhantomPosition.getX(),"Length") <<
+              "DY= "<< G4BestUnit(detectorToPhantomPosition.getY(),"Length") <<
+              "DZ= "<< G4BestUnit(detectorToPhantomPosition.getZ(),"Length") << G4endl;
+*/
         return true;
+}
 /////////////////////////////////////////////////////////////////////////////
+  G4bool SetNumberOfVoxelBySize(G4double sizeX, G4double sizeY, G4double sizeZ);
+  G4bool SetDetectorSize(G4double sizeX, G4double sizeY, G4double sizeZ);
+  G4bool SetPhantomSize(G4double sizeX, G4double sizeY, G4double sizeZ);
+  G4bool SetPhantomPosition(G4ThreeVector);
+  G4bool SetDetectorToPhantomPosition(G4ThreeVector DetectorToPhantomPosition);
+  G4bool  SetPhantomMaterial(G4String material);
+  void SetVoxelSize(G4double sizeX, G4double sizeY, G4double sizeZ);
+  void SetDetectorSize(G4double sizeX, G4double sizeY, G4double sizeZ);
+  void SetPhantomSize(G4double sizeX, G4double sizeY, G4double sizeZ);
+  void SetPhantomPosition(G4ThreeVector);
+  void SetDetectorToPhantomPosition(G4ThreeVector DetectorToPhantomPosition);
+  void UpdateGeometry();
+  void PrintParameters();
   G4LogicalVolume* GetDetectorLogicalVolume(){ return detectorLogicalVolume;}
 …
   HadrontherapyMatrix*             matrix;
+  G4VPhysicalVolume* phantomPhysicalVolume;
+  G4LogicalVolume*   phantomLogicalVolume;
+  G4LogicalVolume*   detectorLogicalVolume;
+  G4VPhysicalVolume* detectorPhysicalVolume;
+  G4Box *phantom , *detector;
+  G4LogicalVolume *phantomLogicalVolume, *detectorLogicalVolume;
+  G4VPhysicalVolume *phantomPhysicalVolume, *detectorPhysicalVolume;
   G4double phantomSizeX;
 …
   G4int numberOfVoxelsAlongZ;
+  G4Box* phantom;
+  G4Box* detector;
+  G4double volumeOfVoxel, massOfVoxel;
+  G4Material *phantomMaterial, *detectorMaterial;
+  G4Region* aRegion;
 };
 #endif

trunk/examples/advanced/hadrontherapy/include/HadrontherapyDetectorMessenger.hh

-                      r1230
+                      r1313
 class HadrontherapyDetectorConstruction;
 class G4UIdirectory;
-class G4UIcmdWithADoubleAndUnit;
-class G4UIcmdWithAString;
 class G4UIcmdWith3VectorAndUnit;
+class G4UIcmdWithoutParameter;
+class G4UIcmdWithAString;
 class HadrontherapyDetectorMessenger: public G4UImessenger
 …
   G4UIdirectory *changeThePhantomDir,  *changeTheDetectorDir;
+  G4UIcmdWithoutParameter   *updateCmd;
+  G4UIcmdWithAString        *changeThePhantomMaterialCmd;
   G4UIcmdWith3VectorAndUnit *changeThePhantomSizeCmd,
                             *changeThePhantomPositionCmd,

trunk/examples/advanced/hadrontherapy/include/HadrontherapyDetectorSD.hh

r1230	r1313
43	43	~HadrontherapyDetectorSD();
44	44
	45	std::ofstream ofs;
45	46	void Initialize(G4HCofThisEvent*);
46	47

trunk/examples/advanced/hadrontherapy/include/HadrontherapyEventAction.hh

r1230	r1313
59	59	G4String drawFlag; //Visualisation flag
60	60	G4int hitsCollectionID;
61		HadrontherapyMatrix *matrix;
	61	//HadrontherapyMatrix *matrix;
62	62	G4int printModulo;
63	63	HadrontherapyEventActionMessenger* pointerEventMessenger;

trunk/examples/advanced/hadrontherapy/include/HadrontherapyInteractionParameters.hh

-                      r1230
+                      r1313
 #include "G4NistElementBuilder.hh"
+#ifdef G4ANALYSIS_USE_ROOT
+#include "TROOT.h"
+#include "TCanvas.h"
+#include "TFile.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TGraph.h"
+#include "TLegend.h"
+#include "TLegendEntry.h"
+#include "TStyle.h"
+#endif
 class HadrontherapyDetectorConstruction;
 class HadrontherapyParameterMessenger;
+class G4ParticleDefinition;
+class G4Material;
 class HadrontherapyInteractionParameters : public G4EmCalculator
+{
 public:
     HadrontherapyInteractionParameters();
         ~HadrontherapyInteractionParameters();
+    HadrontherapyInteractionParameters(G4bool);
+    ~HadrontherapyInteractionParameters();
 // Get data for Mass SP (MeV*cm2/g)
+// Get data for Mass SP
 // G4NistMaterialBuilder class materials
 // User must provide: material kinetic energy lower limit, kinetic energy upper limit, number of points to retrieve,
 // [particle], [output filename].
+    bool GetStoppingTable (const G4String& vararg);
+    G4bool GetStoppingTable (const G4String& vararg);
+    G4double GetStopping (G4double energy,
+                          const G4ParticleDefinition*,
+                          const G4Material*,
+                          G4double density = 0.);
+#ifdef G4ANALYSIS_USE_ROOT
+    void PlotStopping(const G4String&);
+#endif
     void ListOfNistMaterials (const G4String& vararg);
     void BeamOn();
 …
 private:
     G4Material* GetNistMaterial(G4String material);
     G4NistElementBuilder* nistEle;
     G4NistMaterialBuilder* nistMat;
-    G4double kinEmin, kinEmax, npoints;
-    G4String particle, material, filename;
     std::ofstream outfile;
     std::ostream data;
     G4Material* Pmaterial;
-    G4double density;
-    G4EmCalculator* emCal;
     HadrontherapyParameterMessenger* pMessenger;
     bool beamFlag;
+#ifdef G4ANALYSIS_USE_ROOT
+    TCanvas *theRootCanvas;
+    TGraph *theRootGraph;
+    TAxis *axisX, *axisY;
+#endif
+    G4double kinEmin, kinEmax, npoints;
+    G4String particle, material, filename;
+    G4double dedxtot, density;
+    std::vector<G4double> energy;
+    std::vector<G4double> massDedx;
 };
 #endif

trunk/examples/advanced/hadrontherapy/include/HadrontherapyMatrix.hh

-                      r1230
+                      r1313
 #ifndef HadrontherapyMatrix_H
 #define HadrontherapyMatrix_H 1
+#include <G4ParticleDefinition.hh>
 #include "globals.hh"
 #include <vector>
+#include <fstream>
 // The information: energy deposit and position in the phantom
 // is stored in a matrix
+// type struct useful to store nucludes data
+struct ion
+ {
+     G4bool isPrimary;   // true if particle is primary
+     G4int PDGencoding;  // Particle data group id for the particle
+     //G4String extName; //  AZ[excitation energy]: like He3[1277.4], He4[0.0], Li7[231.4], ...
+     G4String name;      // simple name without excitation energy: He3, He4, Li7, ...
+     G4int len;          // name length
+     G4int Z;            // atomic number
+     G4int A;            // mass number
+     G4double *dose;     // pointer to dose matrix
+     unsigned int    *fluence;  // pointer to fluence matrix
+     //friend bool operator<(const ion& a, const ion& b) {return (a.Z == b.Z) ? b.A < a.A : b.Z < a.Z ;}
+     G4bool operator<(const ion& a) const{return (this->Z == a.Z) ? this-> A < a.A : this->Z < a.Z ;}
+ };
 class HadrontherapyMatrix
+{
 private:
+  HadrontherapyMatrix(G4int voxelX, G4int voxelY, G4int voxelZ); //<--- this is supposed to be a singleton
+  HadrontherapyMatrix(G4int numberOfVoxelAlongX,
+                      G4int numberOfVoxelAlongY,
+                      G4int numberOfVoxelAlongZ,
+                      G4double massOfVoxel); //< this is supposed to be a singleton
 public:
   ~HadrontherapyMatrix();
+// Get object instance only
+  static HadrontherapyMatrix* getInstance();
+// Make & Get instance
+  static HadrontherapyMatrix* getInstance(G4int nX, G4int nY, G4int nZ);
+  void flush();
+  // Get object instance only
+  static HadrontherapyMatrix* GetInstance();
+  // Make & Get instance
+  static HadrontherapyMatrix* GetInstance(G4int nX, G4int nY, G4int nZ, G4double mass);
+  static G4bool secondaries;
+  // Full list of generated nuclides
+  void PrintNuclides();
+  // Hit array marker (useful to avoid multiple counts of fluence)
+  void ClearHitTrack();
+  G4int* GetHitTrack(G4int i, G4int j, G4int k);
+  // All the elements of the matrix are initialised to zero
   void Initialize();
+  // All the elements of the matrix are initialised to zero
+  void Clear();
+  // Fill DOSE/fluence matrix for particle:
+  // if fluence parameter is true then fluence at voxel (i, j, k) is increased
+  // else energyDeposit fill the dose matrix for voxel (i,j,k)
+  G4bool Fill(G4int, G4ParticleDefinition* particleDef, G4int i, G4int j, G4int k, G4double energyDeposit, G4bool fluence=false);
+  // Fill TOTAL DOSE matrix for primary particles only
   void Fill(G4int i, G4int j, G4int k, G4double energyDeposit);
   // The matrix is filled with the energy deposit
   // in the element corresponding to the voxel of the phantom where
   // the energy deposit was registered
+  void TotalEnergyDeposit();
   // Store the information of the matrix in a ntuple and in
   // a 1D Histogram
+  inline G4int Index(G4int i, G4int j, G4int k){ return (i * numberVoxelY + j) * numberVoxelZ + k; }
+  // Get a unique index from a three dimensional voxel information
+  void TotalEnergyDeposit();
+  // Store single matrix data to filename
+  void StoreMatrix(G4String file, void* data,size_t psize);
+  // Store all fluence data to filenames
+  void StoreFluenceData();
+  // Store all dose data to filenames
+  void StoreDoseData();
+  // Store all data (except the total dose) to ONE filename
+  void StoreDoseFluenceAscii();
+#ifdef G4ANALYSIS_USE_ROOT
+  void StoreDoseFluenceRoot();
+#endif
+  inline G4int Index(G4int i, G4int j, G4int k) { return (i * numberOfVoxelAlongY + j) * numberOfVoxelAlongZ + k; }
+  // Get a unique index from  a three dimensional one
+  G4double * GetMatrix(){return matrix;}
+  G4int GetNvoxel(){return numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ;}
+  // Total number of voxels read only access
+  G4int GetNumberOfVoxelAlongX(){return numberOfVoxelAlongX;}
+  G4int GetNumberOfVoxelAlongY(){return numberOfVoxelAlongY;}
+  G4int GetNumberOfVoxelAlongZ(){return numberOfVoxelAlongZ;}
 private:
   static HadrontherapyMatrix* instance;
+  G4int numberVoxelX;
+  G4int numberVoxelY;
+  G4int numberVoxelZ;
+  G4int numberOfVoxelAlongX;
+  G4int numberOfVoxelAlongY;
+  G4int numberOfVoxelAlongZ;
+  G4double massOfVoxel;
   G4double* matrix;
+  G4int* hitTrack;
+  G4String filename;
+  std::ofstream ofs;
+  // Dose&fluence data store
+  std::vector <ion> ionStore;
+  // want secondaries?
+  G4double doseUnit;
 };
 #endif

trunk/examples/advanced/hadrontherapy/include/HadrontherapyPhysicsList.hh

r1230	r1313
51	51	void SetCutForElectron(G4double);
52	52	void SetCutForPositron(G4double);
53
	53	void SetDetectorCut(G4double cut);
54	54	void AddPhysicsList(const G4String& name);
55	55	void ConstructProcess();

trunk/examples/advanced/hadrontherapy/include/HadrontherapyPhysicsListMessenger.hh

r1230	r1313
58	58	G4UIcmdWithADoubleAndUnit* protoCutCmd;
59	59	G4UIcmdWithADoubleAndUnit* allCutCmd;
	60	G4UIcmdWithADoubleAndUnit* allDetectorCmd;
60	61	G4UIcmdWithAString* pListCmd;
61	62	G4UIcmdWithAString* packageListCmd;

trunk/examples/advanced/hadrontherapy/include/HadrontherapyPrimaryGeneratorAction.hh

-                      r1230
+                      r1313
   void SetsigmaMomentumZ(G4double);
   G4double GetmeanKineticEnergy(void);
+  G4ParticleGun *GetParticleGun(void){return particleGun;}
 private:
 …
 private:
   G4ParticleGun*                particleGun;
+  G4ParticleGun*                          particleGun;
   HadrontherapyPrimaryGeneratorMessenger* gunMessenger;
   G4double sigmaX;
+  G4double                                sigmaX;
 };

trunk/examples/advanced/hadrontherapy/include/PassiveProtonBeamLine.hh

-                      r1230
+                      r1313
 private:
   void SetDimensions(); //passive proton line dimensions
+//passive proton line dimensions
+  void SetDefaultDimensions();
   void ConstructPassiveProtonBeamLine();

trunk/examples/advanced/hadrontherapy/include/PassiveProtonBeamLineMessenger.hh

-                      r1230
+                      r1313
 // ********************************************************************
 //
 // HadrontherapyDetectorMessenger.hh;
+// PassiveProtonBeamLineMessenger.hh;
 // See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
 #ifndef HadrontherapyDetectorMessenger_h
 #define HadrontherapyDetectorMessenger_h 1
+#ifndef PassiveProtonBeamLineMessenger_h
+#define PassiveProtonBeamLineMessenger_h 1
 #include "globals.hh"

trunk/examples/advanced/hadrontherapy/src/HadrontherapyAnalysisManager.cc

-                      r1230
+                      r1313
 //
 // ********************************************************************
 // * License and Disclaimer                                           *
 // *                                                                  *
 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
 // * conditions of the Geant4 Software License,  included in the file *
 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
 // * include a list of copyright holders.                             *
 // *                                                                  *
 // * Neither the authors of this software system, nor their employing *
 // * institutes,nor the agencies providing financial support for this *
 // * work  make  any representation or  warranty, express or implied, *
 // * regarding  this  software system or assume any liability for its *
 // * use.  Please see the license in the file  LICENSE  and URL above *
 // * for the full disclaimer and the limitation of liability.         *
 // *                                                                  *
 // * This  code  implementation is the result of  the  scientific and *
 // * technical work of the GEANT4 collaboration.                      *
 // * By using,  copying,  modifying or  distributing the software (or *
 // * any work based  on the software)  you  agree  to acknowledge its *
 // * use  in  resulting  scientific  publications,  and indicate your *
 // * acceptance of all terms of the Geant4 Software license.          *
 // ********************************************************************
 //
 // $Id: HadrontherapyAnalisysManager.cc;
 // See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
+        //
+        // ********************************************************************
+        // * License and Disclaimer                                           *
+        // *                                                                  *
+        // * The  Geant4 software  is  copyright of the Copyright Holders  of *
+        // * the Geant4 Collaboration.  It is provided  under  the terms  and *
+        // * conditions of the Geant4 Software License,  included in the file *
+        // * LICENSE and available at  http://cern.ch/geant4/license .  These *
+        // * include a list of copyright holders.                             *
+        // *                                                                  *
+        // * Neither the authors of this software system, nor their employing *
+        // * institutes,nor the agencies providing financial support for this *
+        // * work  make  any representation or  warranty, express or implied, *
+        // * regarding  this  software system or assume any liability for its *
+        // * use.  Please see the license in the file  LICENSE  and URL above *
+        // * for the full disclaimer and the limitation of liability.         *
+        // *                                                                  *
+        // * This  code  implementation is the result of  the  scientific and *
+        // * technical work of the GEANT4 collaboration.                      *
+        // * By using,  copying,  modifying or  distributing the software (or *
+        // * any work based  on the software)  you  agree  to acknowledge its *
+        // * use  in  resulting  scientific  publications,  and indicate your *
+        // * acceptance of all terms of the Geant4 Software license.          *
+        // ********************************************************************
+        //
+        // $Id: HadrontherapyAnalisysManager.cc;
+        // See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
 #include "HadrontherapyAnalysisManager.hh"
 …
 #include "HadrontherapyAnalysisFileMessenger.hh"
 #include <time.h>
+#ifdef ANALYSIS_USE
 HadrontherapyAnalysisManager* HadrontherapyAnalysisManager::instance = 0;
+#ifdef G4ANALYSIS_USE_ROOT
+#undef G4ANALYSIS_USE
+#endif
+/////////////////////////////////////////////////////////////////////////////
+#ifdef G4ANALYSIS_USE
+HadrontherapyAnalysisManager::HadrontherapyAnalysisManager() :
+  analysisFileName("DoseDistribution.root"), aFact(0), theTree(0), histFact(0), tupFact(0), h1(0), h2(0), h3(0),
+  h4(0), h5(0), h6(0), h7(0), h8(0), h9(0), h10(0), h11(0), h12(0), h13(0), h14(0), h15(0), h16(0), ntuple(0),
+  ionTuple(0),
+  fragmentTuple(0),
+  eventCounter(0)
+{
+        fMess = new HadrontherapyAnalysisFileMessenger(this);
+}
+#endif
+#ifdef G4ANALYSIS_USE_ROOT
+HadrontherapyAnalysisManager::HadrontherapyAnalysisManager() :
+  analysisFileName("DoseDistribution.root"),theTFile(0), histo1(0), histo2(0), histo3(0),
+  histo4(0), histo5(0), histo6(0), histo7(0), histo8(0), histo9(0), histo10(0), histo11(0), histo12(0), histo13(0), histo14(0), histo15(0), histo16(0),
+  theROOTNtuple(0),
+  theROOTIonTuple(0),
+  fragmentNtuple(0),
+  metaData(0),
+  eventCounter(0)
+{
+  fMess = new HadrontherapyAnalysisFileMessenger(this);
+HadrontherapyAnalysisManager::HadrontherapyAnalysisManager():
+#ifdef G4ANALYSIS_USE_ROOT
+analysisFileName("DoseDistribution.root"),theTFile(0), histo1(0), histo2(0), histo3(0),
+histo4(0), histo5(0), histo6(0), histo7(0), histo8(0), histo9(0), histo10(0), histo11(0), histo12(0), histo13(0), histo14(0), histo15(0), histo16(0),
+kinFragNtuple(0),
+kineticEnergyPrimaryNtuple(0),
+doseFragNtuple(0),
+fluenceFragNtuple(0),
+letFragNtuple(0),
+theROOTNtuple(0),
+theROOTIonTuple(0),
+fragmentNtuple(0),
+metaData(0),
+eventCounter(0)
+{
+    fMess = new HadrontherapyAnalysisFileMessenger(this);
+}
 #endif
 /////////////////////////////////////////////////////////////////////////////
 HadrontherapyAnalysisManager::~HadrontherapyAnalysisManager()
+{
+delete(fMess); //kill the messenger
+#ifdef G4ANALYSIS_USE
+  delete fragmentTuple;
+  fragmentTuple = 0;
+  delete ionTuple;
+  ionTuple = 0;
+  delete ntuple;
+  ntuple = 0;
+  delete h16;
+  h16 = 0;
+  delete h15;
+  h15 = 0;
+  delete h14;
+  h14 = 0;
+  delete h13;
+  h13 = 0;
+  delete h12;
+  h12 = 0;
+  delete h11;
+  h11 = 0;
+  delete h10;
+  h10 = 0;
+  delete h9;
+  h9 = 0;
+  delete h8;
+  h8 = 0;
+  delete h7;
+  h7 = 0;
+  delete h6;
+  h6 = 0;
+  delete h5;
+  h5 = 0;
+  delete h4;
+  h4 = 0;
+  delete h3;
+  h3 = 0;
+  delete h2;
+  h2 = 0;
+  delete h1;
+  h1 = 0;
+  delete tupFact;
+  tupFact = 0;
+  delete histFact;
+  histFact = 0;
+  delete theTree;
+  theTree = 0;
+  delete aFact;
+  aFact = 0;
+    delete(fMess);
+#ifdef G4ANALYSIS_USE_ROOT
+    Clear();
 #endif
+}
+HadrontherapyAnalysisManager* HadrontherapyAnalysisManager::GetInstance()
+{
+        if (instance == 0) instance = new HadrontherapyAnalysisManager;
+        return instance;
+}
 #ifdef G4ANALYSIS_USE_ROOT
+  delete metaData;
+  metaData = 0;
+  delete fragmentNtuple;
+  fragmentNtuple = 0;
+  delete theROOTIonTuple;
+  theROOTIonTuple = 0;
+  delete theROOTNtuple;
+  theROOTNtuple = 0;
+  delete histo16;
+  histo14 = 0;
+  delete histo15;
+  histo14 = 0;
+  delete histo14;
+  histo14 = 0;
+  delete histo13;
+  histo13 = 0;
+  delete histo12;
+  histo12 = 0;
+  delete histo11;
+  histo11 = 0;
+  delete histo10;
+  histo10 = 0;
+  delete histo9;
+  histo9 = 0;
+  delete histo8;
+  histo8 = 0;
+  delete histo7;
+  histo7 = 0;
+  delete histo6;
+  histo6 = 0;
+  delete histo5;
+  histo5 = 0;
+  delete histo4;
+  histo4 = 0;
+  delete histo3;
+  histo3 = 0;
+  delete histo2;
+  histo2 = 0;
+  delete histo1;
+  histo1 = 0;
+void HadrontherapyAnalysisManager::Clear()
+{
+    // XXX delete ALL variables!!
+    // but this seems to be automatically done by
+    // theTFile->Close() command!
+    // so we get a *double free or corruption*
+    delete metaData;
+    metaData = 0;
+    delete fragmentNtuple;
+    fragmentNtuple = 0;
+    delete theROOTIonTuple;
+    theROOTIonTuple = 0;
+    delete theROOTNtuple;
+    theROOTNtuple = 0;
+    delete histo16;
+    histo16 = 0;
+    delete histo15;
+    histo15 = 0;
+    delete histo14;
+    histo14 = 0;
+    delete histo13;
+    histo13 = 0;
+    delete histo12;
+    histo12 = 0;
+    delete histo11;
+    histo11 = 0;
+    delete histo10;
+    histo10 = 0;
+    delete histo9;
+    histo9 = 0;
+    delete histo8;
+    histo8 = 0;
+    delete histo7;
+    histo7 = 0;
+    delete histo6;
+    histo6 = 0;
+    delete histo5;
+    histo5 = 0;
+    delete histo4;
+    histo4 = 0;
+    delete histo3;
+    histo3 = 0;
+    delete histo2;
+    histo2 = 0;
+    delete histo1;
+    histo1 = 0;
+}
+/////////////////////////////////////////////////////////////////////////////
+void HadrontherapyAnalysisManager::SetAnalysisFileName(G4String aFileName)
+{
+        this->analysisFileName = aFileName;
+}
+        /////////////////////////////////////////////////////////////////////////////
+void HadrontherapyAnalysisManager::book()
+{
+        delete theTFile; // this is similar to theTFile->Close() => delete all associated variables created via new, moreover it delete itself.
+        theTFile = new TFile(analysisFileName, "RECREATE");
+        // Create the histograms with the energy deposit along the X axis
+        histo1 = createHistogram1D("braggPeak","slice, energy", 400, 0., 80); //<different waterthicknesses are accoutned for in ROOT-analysis stage
+        histo2 = createHistogram1D("h20","Secondary protons - slice, energy", 400, 0., 400.);
+        histo3 = createHistogram1D("h30","Secondary neutrons - slice, energy", 400, 0., 400.);
+        histo4 = createHistogram1D("h40","Secondary alpha - slice, energy", 400, 0., 400.);
+        histo5 = createHistogram1D("h50","Secondary gamma - slice, energy", 400, 0., 400.);
+        histo6 = createHistogram1D("h60","Secondary electron - slice, energy", 400, 0., 400.);
+        histo7 = createHistogram1D("h70","Secondary triton - slice, energy", 400, 0., 400.);
+        histo8 = createHistogram1D("h80","Secondary deuteron - slice, energy", 400, 0., 400.);
+        histo9 = createHistogram1D("h90","Secondary pion - slice, energy", 400, 0., 400.);
+        histo10 = createHistogram1D("h100","Energy distribution of secondary electrons", 70, 0., 70.);
+        histo11 = createHistogram1D("h110","Energy distribution of secondary photons", 70, 0., 70.);
+        histo12 = createHistogram1D("h120","Energy distribution of secondary deuterons", 70, 0., 70.);
+        histo13 = createHistogram1D("h130","Energy distribution of secondary tritons", 70, 0., 70.);
+        histo14 = createHistogram1D("h140","Energy distribution of secondary alpha particles", 70, 0., 70.);
+        histo15 = createHistogram1D("heliumEnergyAfterPhantom","Energy distribution of secondary helium fragments after the phantom",
+, 0., 500.);
+        histo16 = createHistogram1D("hydrogenEnergyAfterPhantom","Energy distribution of secondary helium fragments after the phantom",
+, 0., 500.);
+        kinFragNtuple  = new TNtuple("kinFragNtuple",
+                "Kinetic energy by voxel & fragment",
+                "i:j:k:A:Z:kineticEnergy");
+        kineticEnergyPrimaryNtuple= new TNtuple("kineticEnergyPrimaryNtuple",
+                "Kinetic energy by voxel of primary",
+                "i:j:k:kineticEnergy");
+        doseFragNtuple = new TNtuple("doseFragNtuple",
+                "Energy deposit by voxel & fragment",
+                "i:j:k:A:Z:energy");
+        fluenceFragNtuple = new TNtuple("fluenceFragNtuple",
+                "Fluence by voxel & fragment",
+                "i:j:k:A:Z:fluence");
+        letFragNtuple = new TNtuple("letFragNtuple",
+                "Let by voxel & fragment",
+                "i:j:k:A:Z:letT:letD");
+        theROOTNtuple =   new TNtuple("theROOTNtuple",
+                "Energy deposit by slice",
+                "i:j:k:energy");
+        theROOTIonTuple = new TNtuple("theROOTIonTuple",
+                "Generic ion information",
+                "a:z:occupancy:energy");
+        fragmentNtuple =  new TNtuple("fragmentNtuple",
+                "Fragments",
+                "A:Z:energy:posX:posY:posZ");
+        metaData =        new TNtuple("metaData",
+                "Metadata",
+                "events:detectorDistance:waterThickness:beamEnergy:energyError:phantomCenterDistance");
+}
+        /////////////////////////////////////////////////////////////////////////////
+void HadrontherapyAnalysisManager::FillEnergyDeposit(G4int i,
+                                                                                                         G4int j,
+                                                                                                         G4int k,
+                                                                                                         G4double energy)
+{
+        if (theROOTNtuple) {
+                theROOTNtuple->Fill(i, j, k, energy);
+        }
+}
+        /////////////////////////////////////////////////////////////////////////////
+void HadrontherapyAnalysisManager::BraggPeak(G4int slice, G4double energy)
+{
+        histo1->SetBinContent(slice, energy); //This uses setbincontent instead of fill to get labels correct
+}
+        /////////////////////////////////////////////////////////////////////////////
+void HadrontherapyAnalysisManager::SecondaryProtonEnergyDeposit(G4int slice, G4double energy)
+{
+        histo2->Fill(slice, energy);
+}
+        /////////////////////////////////////////////////////////////////////////////
+void HadrontherapyAnalysisManager::SecondaryNeutronEnergyDeposit(G4int slice, G4double energy)
+{
+        histo3->Fill(slice, energy);
+}
+        /////////////////////////////////////////////////////////////////////////////
+void HadrontherapyAnalysisManager::SecondaryAlphaEnergyDeposit(G4int slice, G4double energy)
+{
+        histo4->Fill(slice, energy);
+}
+        /////////////////////////////////////////////////////////////////////////////
+void HadrontherapyAnalysisManager::SecondaryGammaEnergyDeposit(G4int slice, G4double energy)
+{
+        histo5->Fill(slice, energy);
+}
+        /////////////////////////////////////////////////////////////////////////////
+void HadrontherapyAnalysisManager::SecondaryElectronEnergyDeposit(G4int slice, G4double energy)
+{
+        histo6->Fill(slice, energy);
+}
+        /////////////////////////////////////////////////////////////////////////////
+void HadrontherapyAnalysisManager::SecondaryTritonEnergyDeposit(G4int slice, G4double energy)
+{
+        histo7->Fill(slice, energy);
+}
+        /////////////////////////////////////////////////////////////////////////////
+void HadrontherapyAnalysisManager::SecondaryDeuteronEnergyDeposit(G4int slice, G4double energy)
+{
+        histo8->Fill(slice, energy);
+}
+        /////////////////////////////////////////////////////////////////////////////
+void HadrontherapyAnalysisManager::SecondaryPionEnergyDeposit(G4int slice, G4double energy)
+{
+        histo9->Fill(slice, energy);
+}
+        /////////////////////////////////////////////////////////////////////////////
+void HadrontherapyAnalysisManager::electronEnergyDistribution(G4double energy)
+{
+        histo10->Fill(energy);
+}
+        /////////////////////////////////////////////////////////////////////////////
+void HadrontherapyAnalysisManager::gammaEnergyDistribution(G4double energy)
+{
+        histo11->Fill(energy);
+}
+        /////////////////////////////////////////////////////////////////////////////
+void HadrontherapyAnalysisManager::deuteronEnergyDistribution(G4double energy)
+{
+        histo12->Fill(energy);
+}
+        /////////////////////////////////////////////////////////////////////////////
+void HadrontherapyAnalysisManager::tritonEnergyDistribution(G4double energy)
+{
+        histo13->Fill(energy);
+}
+        /////////////////////////////////////////////////////////////////////////////
+void HadrontherapyAnalysisManager::alphaEnergyDistribution(G4double energy)
+{
+        histo14->Fill(energy);
+}
+        /////////////////////////////////////////////////////////////////////////////
+void HadrontherapyAnalysisManager::heliumEnergy(G4double secondaryParticleKineticEnergy)
+{
+        histo15->Fill(secondaryParticleKineticEnergy);
+}
+        /////////////////////////////////////////////////////////////////////////////
+void HadrontherapyAnalysisManager::hydrogenEnergy(G4double secondaryParticleKineticEnergy)
+{
+        histo16->Fill(secondaryParticleKineticEnergy);
+}
+        /////////////////////////////////////////////////////////////////////////////
+        // FillKineticFragmentTuple create an ntuple where the voxel indexs, the atomic number and mass and the kinetic
+        // energy of all the particles interacting with the phantom, are stored
+void HadrontherapyAnalysisManager::FillKineticFragmentTuple(G4int i, G4int j, G4int k, G4int A, G4double Z, G4double kinEnergy)
+{
+    kinFragNtuple -> Fill(i, j, k, A, Z, kinEnergy);
+}
+        /////////////////////////////////////////////////////////////////////////////
+        // FillKineticEnergyPrimaryNTuple creates a ntuple where the voxel indexs and the kinetic
+        // energies of ONLY primary particles interacting with the phantom, are stored
+void HadrontherapyAnalysisManager::FillKineticEnergyPrimaryNTuple(G4int i, G4int j, G4int k, G4double kinEnergy)
+{
+    kineticEnergyPrimaryNtuple -> Fill(i, j, k, kinEnergy);
+}
+        /////////////////////////////////////////////////////////////////////////////
+        // This function is called only if ROOT is activated.
+        // It is called by the HadrontherapyMatric.cc class file and it is used to create two ntuples containing
+        // the total energy deposited and the fluence values, in each voxel and per any particle (primary beam
+        // and secondaries )
+void HadrontherapyAnalysisManager::FillVoxelFragmentTuple(G4int i, G4int j, G4int k, G4int A, G4double Z, G4double energy, G4double fluence)
+{
+                // Fill the ntuple containing the voxel, mass and atomic number and the energy deposited
+    doseFragNtuple ->    Fill( i, j, k, A, Z, energy );
+                // Fill the ntuple containing the voxel, mass and atomic number and the fluence
+        if (i==1 && Z==1) {
+                fluenceFragNtuple -> Fill( i, j, k, A, Z, fluence );
+        }
+}
+void HadrontherapyAnalysisManager::FillLetFragmentTuple(G4int i, G4int j, G4int k, G4int A, G4double Z, G4double letT, G4double letD)
+{
+        letFragNtuple -> Fill( i, j, k, A, Z, letT, letD);
+}
+        /////////////////////////////////////////////////////////////////////////////
+void HadrontherapyAnalysisManager::FillFragmentTuple(G4int A, G4double Z, G4double energy, G4double posX, G4double posY, G4double posZ)
+{
+        fragmentNtuple->Fill(A, Z, energy, posX, posY, posZ);
+}
+        /////////////////////////////////////////////////////////////////////////////
+void HadrontherapyAnalysisManager::genericIonInformation(G4int a,
+                                                                                                                 G4double z,
+                                                                                                                 G4int electronOccupancy,
+                                                                                                                 G4double energy)
+{
+        if (theROOTIonTuple) {
+                theROOTIonTuple->Fill(a, z, electronOccupancy, energy);
+        }
+}
+        /////////////////////////////////////////////////////////////////////////////
+void HadrontherapyAnalysisManager::startNewEvent()
+{
+        eventCounter++;
+}
+        /////////////////////////////////////////////////////////////////////////////
+void HadrontherapyAnalysisManager::setGeometryMetaData(G4double endDetectorPosition, G4double waterThickness, G4double phantomCenter)
+{
+        this->detectorDistance = endDetectorPosition;
+        this->phantomDepth = waterThickness;
+        this->phantomCenterDistance = phantomCenter;
+}
+void HadrontherapyAnalysisManager::setBeamMetaData(G4double meanKineticEnergy,G4double sigmaEnergy)
+{
+        this->beamEnergy = meanKineticEnergy;
+        this->energyError = sigmaEnergy;
+}
+/////////////////////////////////////////////////////////////////////////////
+// Flush data & close the file
+void HadrontherapyAnalysisManager::flush()
+{
+    if (theTFile)
+    {
+        theTFile -> Write();
+        theTFile -> Close();
+    }
+    eventCounter = 0;
+}
 #endif
+}
-/////////////////////////////////////////////////////////////////////////////
-HadrontherapyAnalysisManager* HadrontherapyAnalysisManager::getInstance()
+{
-  if (instance == 0) instance = new HadrontherapyAnalysisManager;
-  return instance;
+}
-/////////////////////////////////////////////////////////////////////////////
-void HadrontherapyAnalysisManager::SetAnalysisFileName(G4String aFileName)
+{
-  this->analysisFileName = aFileName;
+}
-/////////////////////////////////////////////////////////////////////////////
-void HadrontherapyAnalysisManager::book()
+{
-#ifdef G4ANALYSIS_USE
-  // Build up  the  analysis factory
-  aFact = AIDA_createAnalysisFactory();
-  AIDA::ITreeFactory* treeFact = aFact -> createTreeFactory();
-  // Create the .hbk or the .root file
-  G4String fileName = "DoseDistribution.hbk";
-  std::string opts = "export=root";
-  theTree = treeFact -> create(fileName,"hbook",false,true);
-  theTree = treeFact -> create(analysisFileName,"ROOT",false,true,opts);
-  // Factories are not "managed" by an AIDA analysis system.
-  // They must be deleted by the AIDA user code.
-  delete treeFact;
-  // Create the histogram and the ntuple factory
-  histFact = aFact -> createHistogramFactory(*theTree);
-  tupFact = aFact -> createTupleFactory(*theTree);
-  // Create the histograms with the enrgy deposit along the X axis
-  h1 = histFact -> createHistogram1D("10","slice, energy", 400, 0., 400. );
-  h2 = histFact -> createHistogram1D("20","Secondary protons - slice, energy", 400, 0., 400. );
-  h3 = histFact -> createHistogram1D("30","Secondary neutrons - slice, energy", 400, 0., 400. );
-  h4 = histFact -> createHistogram1D("40","Secondary alpha - slice, energy", 400, 0., 400. );
-  h5 = histFact -> createHistogram1D("50","Secondary gamma - slice, energy", 400, 0., 400. );
-  h6 = histFact -> createHistogram1D("60","Secondary electron - slice, energy", 400, 0., 400. );
-  h7 = histFact -> createHistogram1D("70","Secondary triton - slice, energy", 400, 0., 400. );
-  h8 = histFact -> createHistogram1D("80","Secondary deuteron - slice, energy", 400, 0., 400. );
-  h9 = histFact -> createHistogram1D("90","Secondary pion - slice, energy", 400, 0., 400. );
-  h10 = histFact -> createHistogram1D("100","Energy distribution of secondary electrons", 70, 0., 70. );
-  h11 = histFact -> createHistogram1D("110","Energy distribution of secondary photons", 70, 0., 70. );
-  h12 = histFact -> createHistogram1D("120","Energy distribution of secondary deuterons", 70, 0., 70. );
-  h13 = histFact -> createHistogram1D("130","Energy distribution of secondary tritons", 70, 0., 70. );
-  h14 = histFact -> createHistogram1D("140","Energy distribution of secondary alpha particles", 70, 0., 70. );
-  h15 = histFact -> createHistogram1D("150","Energy distribution of helium fragments after the phantom", 70, 0., 500.);
-  h16 = histFact -> createHistogram1D("160","Energy distribution of hydrogen fragments after the phantom", 70, 0., 500.);
-  // Create the ntuple
-  G4String columnNames = "int i; int j; int k; double energy;";
-  G4String options = "";
-  if (tupFact) ntuple = tupFact -> create("1","1",columnNames, options);
-  // Create the ntuple
-  G4String columnNames2 = "int a; double z;  int occupancy; double energy;";
-  G4String options2 = "";
-  if (tupFact) ionTuple = tupFact -> create("2","2", columnNames2, options2);
-  // Create the fragment ntuple
-  G4String columnNames3 = "int a; double z; double energy; double posX; double posY; double posZ;";
-  G4String options3 = "";
-  if (tupFact) fragmentTuple = tupFact -> create("3","3", columnNames3, options3);
-#endif
-#ifdef G4ANALYSIS_USE_ROOT
-  // Use ROOT
-  theTFile = new TFile(analysisFileName, "RECREATE");
-  // Create the histograms with the energy deposit along the X axis
-  histo1 = createHistogram1D("braggPeak","slice, energy", 400, 0., 27.9); //<different waterthicknesses are accoutned for in ROOT-analysis stage
-  histo2 = createHistogram1D("h20","Secondary protons - slice, energy", 400, 0., 400.);
-  histo3 = createHistogram1D("h30","Secondary neutrons - slice, energy", 400, 0., 400.);
-  histo4 = createHistogram1D("h40","Secondary alpha - slice, energy", 400, 0., 400.);
-  histo5 = createHistogram1D("h50","Secondary gamma - slice, energy", 400, 0., 400.);
-  histo6 = createHistogram1D("h60","Secondary electron - slice, energy", 400, 0., 400.);
-  histo7 = createHistogram1D("h70","Secondary triton - slice, energy", 400, 0., 400.);
-  histo8 = createHistogram1D("h80","Secondary deuteron - slice, energy", 400, 0., 400.);
-  histo9 = createHistogram1D("h90","Secondary pion - slice, energy", 400, 0., 400.);
-  histo10 = createHistogram1D("h100","Energy distribution of secondary electrons", 70, 0., 70.);
-  histo11 = createHistogram1D("h110","Energy distribution of secondary photons", 70, 0., 70.);
-  histo12 = createHistogram1D("h120","Energy distribution of secondary deuterons", 70, 0., 70.);
-  histo13 = createHistogram1D("h130","Energy distribution of secondary tritons", 70, 0., 70.);
-  histo14 = createHistogram1D("h140","Energy distribution of secondary alpha particles", 70, 0., 70.);
-  histo15 = createHistogram1D("heliumEnergyAfterPhantom","Energy distribution of secondary helium fragments after the phantom",
-, 0., 500.);
-  histo16 = createHistogram1D("hydrogenEnergyAfterPhantom","Energy distribution of secondary helium fragments after the phantom",
-, 0., 500.);
-  theROOTNtuple = new TNtuple("theROOTNtuple", "Energy deposit by slice", "i:j:k:energy");
-  theROOTIonTuple = new TNtuple("theROOTIonTuple", "Generic ion information", "a:z:occupancy:energy");
-  fragmentNtuple = new TNtuple("fragmentNtuple", "Fragments", "A:Z:energy:posX:posY:posZ");
-  metaData = new TNtuple("metaData", "Metadata", "events:detectorDistance:waterThickness:beamEnergy:energyError:phantomCenterDistance");
-#endif
+}
-/////////////////////////////////////////////////////////////////////////////
-void HadrontherapyAnalysisManager::FillEnergyDeposit(G4int i,
-                                                     G4int j,
-                                                     G4int k,
-                                                     G4double energy)
+{
-#ifdef G4ANALYSIS_USE
- if (ntuple)     {
-       G4int iSlice = ntuple -> findColumn("i");
-       G4int jSlice = ntuple -> findColumn("j");
-       G4int kSlice = ntuple -> findColumn("k");
-       G4int iEnergy = ntuple -> findColumn("energy");
-       ntuple -> fill(iSlice,i);
-       ntuple -> fill(jSlice,j);
-       ntuple -> fill(kSlice,k);
-       ntuple -> fill(iEnergy, energy);     }
-  ntuple -> addRow();
-#endif
-#ifdef G4ANALYSIS_USE_ROOT
-  if (theROOTNtuple) {
-    theROOTNtuple->Fill(i, j, k, energy);
+  }
-#endif
+}
-/////////////////////////////////////////////////////////////////////////////
-void HadrontherapyAnalysisManager::BraggPeak(G4int slice, G4double energy)
+{
-#ifdef G4ANALYSIS_USE
-  h1 -> fill(slice,energy);
-#endif
-#ifdef G4ANALYSIS_USE_ROOT
-  histo1->SetBinContent(slice, energy); //This uses setbincontent instead of fill to get labels correct
-#endif
+}
-/////////////////////////////////////////////////////////////////////////////
-void HadrontherapyAnalysisManager::SecondaryProtonEnergyDeposit(G4int slice, G4double energy)
+{
-#ifdef G4ANALYSIS_USE
-  h2 -> fill(slice,energy);
-#endif
-#ifdef G4ANALYSIS_USE_ROOT
-  histo2->Fill(slice, energy);
-#endif
+}
-/////////////////////////////////////////////////////////////////////////////
-void HadrontherapyAnalysisManager::SecondaryNeutronEnergyDeposit(G4int slice, G4double energy)
+{
-#ifdef G4ANALYSIS_USE
-  h3 -> fill(slice,energy);
-#endif
-#ifdef G4ANALYSIS_USE_ROOT
-  histo3->Fill(slice, energy);
-#endif
+}
-/////////////////////////////////////////////////////////////////////////////
-void HadrontherapyAnalysisManager::SecondaryAlphaEnergyDeposit(G4int slice, G4double energy)
+{
-#ifdef G4ANALYSIS_USE
-  h4 -> fill(slice,energy);
-#endif
-#ifdef G4ANALYSIS_USE_ROOT
-  histo4->Fill(slice, energy);
-#endif
+}
-/////////////////////////////////////////////////////////////////////////////
-void HadrontherapyAnalysisManager::SecondaryGammaEnergyDeposit(G4int slice, G4double energy)
+{
-#ifdef G4ANALYSIS_USE
-  h5 -> fill(slice,energy);
-#endif
-#ifdef G4ANALYSIS_USE_ROOT
-  histo5->Fill(slice, energy);
-#endif
+}
-/////////////////////////////////////////////////////////////////////////////
-void HadrontherapyAnalysisManager::SecondaryElectronEnergyDeposit(G4int slice, G4double energy)
+{
-#ifdef G4ANALYSIS_USE
-  h6 -> fill(slice,energy);
-#endif
-#ifdef G4ANALYSIS_USE_ROOT
-  histo6->Fill(slice, energy);
-#endif
+}
-/////////////////////////////////////////////////////////////////////////////
-void HadrontherapyAnalysisManager::SecondaryTritonEnergyDeposit(G4int slice, G4double energy)
+{
-#ifdef G4ANALYSIS_USE
-  h7 -> fill(slice,energy);
-#endif
-#ifdef G4ANALYSIS_USE_ROOT
-  histo7->Fill(slice, energy);
-#endif
+}
-/////////////////////////////////////////////////////////////////////////////
-void HadrontherapyAnalysisManager::SecondaryDeuteronEnergyDeposit(G4int slice, G4double energy)
+{
-#ifdef G4ANALYSIS_USE
-  h8 -> fill(slice,energy);
-#endif
-#ifdef G4ANALYSIS_USE_ROOT
-  histo8->Fill(slice, energy);
-#endif
+}
-/////////////////////////////////////////////////////////////////////////////
-void HadrontherapyAnalysisManager::SecondaryPionEnergyDeposit(G4int slice, G4double energy)
+{
-#ifdef G4ANALYSIS_USE
-  h9 -> fill(slice,energy);
-#endif
-#ifdef G4ANALYSIS_USE_ROOT
-  histo9->Fill(slice, energy);
-#endif
+}
-/////////////////////////////////////////////////////////////////////////////
-void HadrontherapyAnalysisManager::electronEnergyDistribution(G4double energy)
+{
-#ifdef G4ANALYSIS_USE
-  h10 -> fill(energy);
-#endif
-#ifdef G4ANALYSIS_USE_ROOT
-  histo10->Fill(energy);
-#endif
+}
-/////////////////////////////////////////////////////////////////////////////
-void HadrontherapyAnalysisManager::gammaEnergyDistribution(G4double energy)
+{
-#ifdef G4ANALYSIS_USE
-  h11 -> fill(energy);
-#endif
-#ifdef G4ANALYSIS_USE_ROOT
-  histo11->Fill(energy);
-#endif
+}
-/////////////////////////////////////////////////////////////////////////////
-void HadrontherapyAnalysisManager::deuteronEnergyDistribution(G4double energy)
+{
-#ifdef G4ANALYSIS_USE
-  h12 -> fill(energy);
-#endif
-#ifdef G4ANALYSIS_USE_ROOT
-  histo12->Fill(energy);
-#endif
+}
-/////////////////////////////////////////////////////////////////////////////
-void HadrontherapyAnalysisManager::tritonEnergyDistribution(G4double energy)
+{
-#ifdef G4ANALYSIS_USE
-  h13 -> fill(energy);
-#endif
-#ifdef G4ANALYSIS_USE_ROOT
-  histo13->Fill(energy);
-#endif
+}
-/////////////////////////////////////////////////////////////////////////////
-void HadrontherapyAnalysisManager::alphaEnergyDistribution(G4double energy)
+{
-#ifdef G4ANALYSIS_USE
-  h14 -> fill(energy);
-#endif
-#ifdef G4ANALYSIS_USE_ROOT
-  histo14->Fill(energy);
-#endif
+}
-/////////////////////////////////////////////////////////////////////////////
-void HadrontherapyAnalysisManager::heliumEnergy(G4double secondaryParticleKineticEnergy)
+{
-#ifdef G4ANALYSIS_USE
-  h15->fill(secondaryParticleKineticEnergy);
-#endif
-#ifdef G4ANALYSIS_USE_ROOT
-  histo15->Fill(secondaryParticleKineticEnergy);
-#endif
+}
-/////////////////////////////////////////////////////////////////////////////
-void HadrontherapyAnalysisManager::hydrogenEnergy(G4double secondaryParticleKineticEnergy)
+{
-#ifdef G4ANALYSIS_USE
-  h16->fill(secondaryParticleKineticEnergy);
-#endif
-#ifdef G4ANALYSIS_USE_ROOT
-  histo16->Fill(secondaryParticleKineticEnergy);
-#endif
+}
-void HadrontherapyAnalysisManager::fillFragmentTuple(G4int A, G4double Z, G4double energy, G4double posX, G4double posY, G4double posZ)
+{
-#ifdef G4ANALYSIS_USE
-  if (fragmentTuple)    {
-       G4int aIndex = fragmentTuple -> findColumn("a");
-       G4int zIndex = fragmentTuple -> findColumn("z");
-       G4int energyIndex = fragmentTuple -> findColumn("energy");
-       G4int posXIndex = fragmentTuple -> findColumn("posX");
-       G4int posYIndex = fragmentTuple -> findColumn("posY");
-       G4int posZIndex = fragmentTuple -> findColumn("posZ");
-       fragmentTuple -> fill(aIndex,A);
-       fragmentTuple -> fill(zIndex,Z);
-       fragmentTuple -> fill(energyIndex, energy);
-       fragmentTuple -> fill(posXIndex, posX);
-       fragmentTuple -> fill(posYIndex, posY);
-       fragmentTuple -> fill(posZIndex, posZ);
-       fragmentTuple -> addRow();
+  }
-#endif
-#ifdef G4ANALYSIS_USE_ROOT
-  //G4cout <<" A = " << A << "  Z = " << Z << " energy = " << energy << G4endl;
-  fragmentNtuple->Fill(A, Z, energy, posX, posY, posZ);
-#endif
+}
-/////////////////////////////////////////////////////////////////////////////
-void HadrontherapyAnalysisManager::genericIonInformation(G4int a,
-                                                         G4double z,
-                                                         G4int electronOccupancy,
-                                                         G4double energy)
+{
-#ifdef G4ANALYSIS_USE
-  if (ionTuple)    {
-       G4int aIndex = ionTuple -> findColumn("a");
-       G4int zIndex = ionTuple -> findColumn("z");
-       G4int electronIndex = ionTuple -> findColumn("occupancy");
-       G4int energyIndex = ionTuple -> findColumn("energy");
-       ionTuple -> fill(aIndex,a);
-      ionTuple -> fill(zIndex,z);
-      ionTuple -> fill(electronIndex, electronOccupancy);
-       ionTuple -> fill(energyIndex, energy);
+     }
-   ionTuple -> addRow();
-#endif
-#ifdef G4ANALYSIS_USE_ROOT
-   if (theROOTIonTuple) {
-     theROOTIonTuple->Fill(a, z, electronOccupancy, energy);
+   }
-#endif
+}
-/////////////////////////////////////////////////////////////////////////////
-void HadrontherapyAnalysisManager::startNewEvent()
+{
-  eventCounter++;
+}
-/////////////////////////////////////////////////////////////////////////////
-void HadrontherapyAnalysisManager::setGeometryMetaData(G4double endDetectorPosition, G4double waterThickness, G4double phantomCenter)
+{
-  this->detectorDistance = endDetectorPosition;
-  this->phantomDepth = waterThickness;
-  this->phantomCenterDistance = phantomCenter;
+}
-void HadrontherapyAnalysisManager::setBeamMetaData(G4double meanKineticEnergy,G4double sigmaEnergy)
+{
-  this->beamEnergy = meanKineticEnergy;
-  this->energyError = sigmaEnergy;
+}
-/////////////////////////////////////////////////////////////////////////////
-void HadrontherapyAnalysisManager::flush()
+{
-  HadrontherapyMatrix* matrix = HadrontherapyMatrix::getInstance();
-  matrix->TotalEnergyDeposit();
-#ifdef G4ANALYSIS_USE
-  theTree -> commit();
-  theTree ->close();
-#endif
-#ifdef G4ANALYSIS_USE_ROOT
-  metaData->Fill((Float_t) eventCounter,(Float_t) detectorDistance, (Float_t) phantomDepth, (Float_t) beamEnergy,(Float_t) energyError, (Float_t) phantomCenterDistance);
-  metaData->Write();
-  theROOTNtuple->Write();
-  theROOTIonTuple->Write();
-  fragmentNtuple->Write();
-  theTFile->Write();
-  //  theTFile->Clear();
-  theTFile->Close();
-#endif
-  eventCounter = 0;
-  matrix->flush();
+}
-/////////////////////////////////////////////////////////////////////////////
-void HadrontherapyAnalysisManager::finish()
+{
-#ifdef G4ANALYSIS_USE
-  // Write all histograms to file
-  theTree -> commit();
-  // Close (will again commit)
-  theTree ->close();
-#endif
-#ifdef G4ANALYSIS_USE_ROOT
-  metaData->Fill((Float_t) eventCounter,(Float_t) detectorDistance, (Float_t) phantomDepth, (Float_t) beamEnergy,(Float_t) energyError, (Float_t) phantomCenterDistance);
-  metaData->Write();
-  theROOTNtuple->Write();
-  theROOTIonTuple->Write();
-  fragmentNtuple->Write();
-  theTFile->Write();
-  theTFile->Close();
-#endif
-  eventCounter = 0;
+}
-#endif

trunk/examples/advanced/hadrontherapy/src/HadrontherapyDetectorConstruction.cc

-                      r1230
+                      r1313
 // See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
 #include "G4SDManager.hh"
 #include "G4RunManager.hh"
+#include "G4GeometryManager.hh"
+#include "G4SolidStore.hh"
+#include "G4PhysicalVolumeStore.hh"
+#include "G4LogicalVolumeStore.hh"
 #include "G4Box.hh"
 #include "G4LogicalVolume.hh"
 …
 #include "G4VisAttributes.hh"
 #include "G4NistManager.hh"
+#include "HadrontherapyDetectorConstruction.hh"
 #include "HadrontherapyDetectorROGeometry.hh"
 #include "HadrontherapyDetectorMessenger.hh"
 #include "HadrontherapyDetectorSD.hh"
-#include "HadrontherapyDetectorConstruction.hh"
 #include "HadrontherapyMatrix.hh"
+#include "HadrontherapyAnalysisManager.hh"
+#include <cmath>
 /////////////////////////////////////////////////////////////////////////////
 HadrontherapyDetectorConstruction::HadrontherapyDetectorConstruction(G4VPhysicalVolume* physicalTreatmentRoom)
   : motherPhys(physicalTreatmentRoom),
+  : motherPhys(physicalTreatmentRoom), // pointer to WORLD volume
     detectorSD(0), detectorROGeometry(0), matrix(0),
     phantomPhysicalVolume(0),
     detectorLogicalVolume(0), detectorPhysicalVolume(0),
     phantomSizeX(20.*cm), phantomSizeY(20.*cm), phantomSizeZ(20.*cm), // Default half dimensions
     detectorSizeX(2.*cm), detectorSizeY(2.*cm), detectorSizeZ(2.*cm),
+    phantomPosition(20.*cm, 0.*cm, 0.*cm),
     detectorToPhantomPosition(0.*cm,18.*cm,18.*cm)// Default displacement of the detector respect to the phantom
+{
+    phantom(0), detector(0),
+    phantomLogicalVolume(0), detectorLogicalVolume(0),
+    phantomPhysicalVolume(0), detectorPhysicalVolume(0),
+    aRegion(0)
+{
+  HadrontherapyAnalysisManager::GetInstance();
   // NOTE! that the HadrontherapyDetectorConstruction class
   // does NOT inherit from G4VUserDetectorConstruction G4 class
 …
   // Default detector voxels size
   // 200 slabs along the beam direction (X)
+  sizeOfVoxelAlongX = 200 *um;
+  sizeOfVoxelAlongY = 2 * detectorSizeY;
+  sizeOfVoxelAlongZ = 2 * detectorSizeZ;
+  // Calculate (and eventually set) detector position by displacement, phantom size and detector size
+  SetDetectorPosition();
+  // Build phantom and associated detector
+  ConstructPhantom();
+  ConstructDetector();
+  // Set number of the detector voxels along X Y and Z directions.
+  // This will construct also the sensitive detector, the ROGeometry
+  // and the matrix where the energy deposited is collected!
+  SetNumberOfVoxelBySize(sizeOfVoxelAlongX, sizeOfVoxelAlongY, sizeOfVoxelAlongZ);
+  sizeOfVoxelAlongX = 200 *um;
+  sizeOfVoxelAlongY = 4 *cm;
+  sizeOfVoxelAlongZ = 4 *cm;
+  // Define here the material of the water phantom and of the detector
+  SetPhantomMaterial("G4_WATER");
+  // Construct geometry (messenger commands)
+  SetDetectorSize(4.*cm, 4.*cm, 4.*cm);
+  SetPhantomSize(40. *cm, 40. *cm, 40. *cm);
+  SetPhantomPosition(G4ThreeVector(20. *cm, 0. *cm, 0. *cm));
+  SetDetectorToPhantomPosition(G4ThreeVector(0. *cm, 18. *cm, 18. *cm));
+  // Write virtual parameters to the real ones and check for consistency
+  UpdateGeometry();
+}
 …
 HadrontherapyDetectorConstruction::~HadrontherapyDetectorConstruction()
+{
     delete detectorROGeometry;// This should be safe in C++ even if the argument is a NULL pointer
+    delete detectorROGeometry;
     delete matrix;
     delete detectorMessenger;
+}
+/////////////////////////////////////////////////////////////////////////////
+// ConstructPhantom() is the method that reconstuct a water box (called phantom
+// (or water phantom) in the usual Medical physicists slang).
+// A water phantom can be considered a good
+// approximation of a an human body.
 void HadrontherapyDetectorConstruction::ConstructPhantom()
+{
+  //----------------------------------------
+  // Phantom:
+  // A box used to approximate tissues
+  //----------------------------------------
+    G4bool isotopes =  false;
+    G4Material* waterNist = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER", isotopes);
+    phantom = new G4Box("Phantom",phantomSizeX, phantomSizeY, phantomSizeZ);
+    // Definition of the solid volume of the Phantom
+    phantom = new G4Box("Phantom",
+                        phantomSizeX/2,
+                        phantomSizeY/2,
+                        phantomSizeZ/2);
+// Definition of the logical volume of the Phantom
     phantomLogicalVolume = new G4LogicalVolume(phantom,
                                              waterNist,
+                                             phantomMaterial,
                                              "phantomLog", 0, 0, 0);
+    // Definition of the physics volume of the Phantom
     phantomPhysicalVolume = new G4PVPlacement(0,
                                             phantomPosition,
 …
     red -> SetVisibility(true);
     red -> SetForceSolid(true);
 //red -> SetForceWireframe(true);
+    //red -> SetForceWireframe(true);
     phantomLogicalVolume -> SetVisAttributes(red);
+}
 /////////////////////////////////////////////////////////////////////////////
+// ConstructDetector() it the method the reconstruct a detector region
+// inside the water phantom. It is a volume, located inside the water phantom
+// and with two coincident faces:
+//
+//           **************************
+//           *   water phantom        *
+//           *                        *
+//           *                        *
+//           *---------------         *
+//  Beam     *              -         *
+//  ----->   * detector     -         *
+//           *              -         *
+//           *---------------         *
+//           *                        *
+//           *                        *
+//           *                        *
+//           **************************
+//
+// The detector is the volume that can be dived in slices or voxelized
+// and in it we can collect a number of usefull information:
+// dose distribution, fluence distribution, LET and so on
 void HadrontherapyDetectorConstruction::ConstructDetector()
+{
+  //-----------
+  // Detector
+  //-----------
+    G4bool isotopes =  false;
+    G4Material* waterNist = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER", isotopes);
+    detector = new G4Box("Detector",detectorSizeX,detectorSizeY,detectorSizeZ);
+    // Definition of the solid volume of the Detector
+    detector = new G4Box("Detector",
+                         detectorSizeX/2,
+                         detectorSizeY/2,
+                         detectorSizeZ/2);
+    // Definition of the logic volume of the Phantom
     detectorLogicalVolume = new G4LogicalVolume(detector,
                                                 waterNist,
+                                                detectorMaterial,
                                                 "DetectorLog",
 ,0,0);
 // Detector is attached by default to the phantom face directly exposed to the beam
     detectorPhysicalVolume = new G4PVPlacement(0,
                                              detectorPosition, // Setted by displacement
                                             "DetectorPhys",
                                              detectorLogicalVolume,
                                              phantomPhysicalVolume,
                                              false,0);
+// Definition of the physical volume of the Phantom
+    detectorPhysicalVolume = new G4PVPlacement(0,
+                                               detectorPosition, // Setted by displacement
+                                               "DetectorPhys",
+                                               detectorLogicalVolume,
+                                               phantomPhysicalVolume,
+                                               false,0);
 // Visualisation attributes of the detector
     skyBlue = new G4VisAttributes( G4Colour(135/255. , 206/255. ,  235/255. ));
     skyBlue -> SetVisibility(true);
     skyBlue -> SetForceSolid(true);
 //skyBlue -> SetForceWireframe(true);
+    //skyBlue -> SetForceWireframe(true);
     detectorLogicalVolume -> SetVisAttributes(skyBlue);
+  // **************
+  // Cut per Region
+  // **************
+}
+/////////////////////////////////////////////////////////////////////////////
+  // A smaller cut is fixed in the phantom to calculate the energy deposit with the
+  // required accuracy
+    if (!aRegion)
+    {
+        aRegion = new G4Region("DetectorLog");
+        detectorLogicalVolume -> SetRegion(aRegion);
+        aRegion -> AddRootLogicalVolume(detectorLogicalVolume);
+    }
+}
+/////////////////////////////////////////////////////////////////////////////
 void  HadrontherapyDetectorConstruction::ConstructSensitiveDetector(G4ThreeVector detectorToWorldPosition)
+{
     // Install new Sensitive Detector and ROGeometry
     delete detectorROGeometry; // this should be safe in C++ also if we have a NULL pointer
+    //if (detectorSD) detectorSD->PrintAll();
+    //delete detectorSD;
     // Sensitive Detector and ReadOut geometry definition
     G4SDManager* sensitiveDetectorManager = G4SDManager::GetSDMpointer();
     G4String sensitiveDetectorName = "Detector";
+    static G4String sensitiveDetectorName = "Detector";
     if (!detectorSD)
+        {
 …
+        }
     // The Read Out Geometry is instantiated
     G4String ROGeometryName = "DetectorROGeometry";
+    static G4String ROGeometryName = "DetectorROGeometry";
     detectorROGeometry = new HadrontherapyDetectorROGeometry(ROGeometryName,
                                                             detectorToWorldPosition,
                                                             detectorSizeX,
                                                             detectorSizeY,
                                                             detectorSizeZ,
+                                                            detectorSizeX/2,
+                                                            detectorSizeY/2,
+                                                            detectorSizeZ/2,
                                                             numberOfVoxelsAlongX,
                                                             numberOfVoxelsAlongY,
 …
+        }
+}
+void  HadrontherapyDetectorConstruction::ParametersCheck()
+{
+    // Check phantom/detector sizes & relative position
+    if (!IsInside(detectorSizeX,
+                detectorSizeY,
+                detectorSizeZ,
+                phantomSizeX,
+                phantomSizeY,
+                phantomSizeZ,
+                detectorToPhantomPosition
+                ))
+        G4Exception("Error at HadrontherapyDetectorConstruction::ParametersCheck(). Detector is not fully inside Phantom!");
+    // Check Detector sizes respect to the voxel ones
+    if ( detectorSizeX < sizeOfVoxelAlongX) {
+        G4Exception("Error at HadrontherapyDetectorConstruction::ParametersCheck(). Detector X size must be bigger or equal than that of Voxel X");
+    }
+    if ( detectorSizeY < sizeOfVoxelAlongY) {
+        G4Exception("Error at HadrontherapyDetectorConstruction::ParametersCheck(). Detector Y size must be bigger or equal than that of Voxel Y");
+    }
+    if ( detectorSizeZ < sizeOfVoxelAlongZ) {
+        G4Exception("Error at HadrontherapyDetectorConstruction::ParametersCheck(). Detector Z size must be bigger or equal than that of Voxel Z");
+    }
+}
 /////////////////
 // MESSENGERS //
 ////////////////
+G4bool HadrontherapyDetectorConstruction::SetNumberOfVoxelBySize(G4double sizeX, G4double sizeY, G4double sizeZ)
+{
+    // Only change positive dimensions
+    // XXX numberOfVoxels must be an integer, warn the user
+    if (sizeX > 0)
+G4bool HadrontherapyDetectorConstruction::SetPhantomMaterial(G4String material)
+{
+    if (G4Material* pMat = G4NistManager::Instance()->FindOrBuildMaterial(material, false) )
+    {
+        if (sizeX > 2*detectorSizeX)
+        phantomMaterial  = pMat;
+        detectorMaterial = pMat;
+        if (detectorLogicalVolume && phantomLogicalVolume)
+        {
+            G4cout << "WARNING: Voxel X size must be smaller or equal than that of detector X" << G4endl;
+            return false;
+            detectorLogicalVolume -> SetMaterial(pMat);
+            phantomLogicalVolume ->  SetMaterial(pMat);
+            G4RunManager::GetRunManager() -> PhysicsHasBeenModified();
+            G4RunManager::GetRunManager() -> GeometryHasBeenModified();
+            G4cout << "The material of Phantom/Detector has been changed to " << material << G4endl;
+        }
+        // Round to the nearest integer
+        numberOfVoxelsAlongX = lrint(2 * detectorSizeX / sizeX);
+        sizeOfVoxelAlongX = (2 * detectorSizeX / numberOfVoxelsAlongX );
+        if(sizeOfVoxelAlongX!=sizeX) G4cout << "Rounding " <<
+                                                G4BestUnit(sizeX, "Length") << " to " <<
+                                                G4BestUnit(sizeOfVoxelAlongX, "Length") << G4endl;
+    }
+    if (sizeY > 0)
+    }
+    else
+    {
+        if (sizeY > 2*detectorSizeY)
+        {
+            G4cout << "WARNING: Voxel Y size must be smaller or equal than that of detector Y" << G4endl;
+            return false;
+        }
+        numberOfVoxelsAlongY = lrint(2 * detectorSizeY / sizeY);
+        sizeOfVoxelAlongY = (2 * detectorSizeY / numberOfVoxelsAlongY );
+        if(sizeOfVoxelAlongY!=sizeY) G4cout << "Rounding " <<
+                                                G4BestUnit(sizeY, "Length") << " to " <<
+                                                G4BestUnit(sizeOfVoxelAlongY, "Length") << G4endl;
+    }
+    if (sizeZ > 0)
+        G4cout << "WARNING: material \"" << material << "\" doesn't exist in NIST elements/materials"
+            " table [located in $G4INSTALL/source/materials/src/G4NistMaterialBuilder.cc]" << G4endl;
+        G4cout << "Use command \"/parameter/nist\" to see full materials list!" << G4endl;
+        return false;
+    }
+    return true;
+}
+/////////////////////////////////////////////////////////////////////////////
+void HadrontherapyDetectorConstruction::SetPhantomSize(G4double sizeX, G4double sizeY, G4double sizeZ)
+{
+    if (sizeX > 0.) phantomSizeX = sizeX;
+    if (sizeY > 0.) phantomSizeY = sizeY;
+    if (sizeZ > 0.) phantomSizeZ = sizeZ;
+}
+/////////////////////////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////////////////////////
+void HadrontherapyDetectorConstruction::SetDetectorSize(G4double sizeX, G4double sizeY, G4double sizeZ)
+{
+    if (sizeX > 0.) {detectorSizeX = sizeX;}
+    if (sizeY > 0.) {detectorSizeY = sizeY;}
+    if (sizeZ > 0.) {detectorSizeZ = sizeZ;}
+    SetVoxelSize(sizeOfVoxelAlongX, sizeOfVoxelAlongY, sizeOfVoxelAlongZ);
+}
+/////////////////////////////////////////////////////////////////////////////
+void HadrontherapyDetectorConstruction::SetVoxelSize(G4double sizeX, G4double sizeY, G4double sizeZ)
+{
+    if (sizeX > 0.) {sizeOfVoxelAlongX = sizeX;}
+    if (sizeY > 0.) {sizeOfVoxelAlongY = sizeY;}
+    if (sizeZ > 0.) {sizeOfVoxelAlongZ = sizeZ;}
+}
+void HadrontherapyDetectorConstruction::SetPhantomPosition(G4ThreeVector pos)
+{
+    phantomPosition = pos;
+}
+/////////////////////////////////////////////////////////////////////////////
+void HadrontherapyDetectorConstruction::SetDetectorToPhantomPosition(G4ThreeVector displ)
+{
+    detectorToPhantomPosition = displ;
+}
+/////////////////////////////////////////////////////////////////////////////
+void HadrontherapyDetectorConstruction::UpdateGeometry()
+{
+    ParametersCheck();
+    //G4RunManager::GetRunManager() -> PhysicsHasBeenModified();
+    G4GeometryManager::GetInstance() -> OpenGeometry();
+    if (phantom)
+    {
+        if (sizeZ > 2*detectorSizeZ)
+        {
+            G4cout << "WARNING: Voxel Z size must be smaller or equal than that of detector Z" << G4endl;
+            return false;
+        }
+        numberOfVoxelsAlongZ = lrint(2 * detectorSizeZ / sizeZ);
+        sizeOfVoxelAlongZ = (2 * detectorSizeZ / numberOfVoxelsAlongZ );
+        if(sizeOfVoxelAlongZ!=sizeZ) G4cout << "Rounding " <<
+                                                G4BestUnit(sizeZ, "Length") << " to " <<
+                                                G4BestUnit(sizeOfVoxelAlongZ, "Length") << G4endl;
+    }
+    G4cout << "The (X, Y, Z) sizes of the Voxels are: (" <<
+                G4BestUnit(sizeOfVoxelAlongX, "Length")  << ", " <<
+                G4BestUnit(sizeOfVoxelAlongY, "Length")  << ", " <<
+                G4BestUnit(sizeOfVoxelAlongZ, "Length") << ')' << G4endl;
+            phantom -> SetXHalfLength(phantomSizeX/2);
+            phantom -> SetYHalfLength(phantomSizeY/2);
+            phantom -> SetZHalfLength(phantomSizeZ/2);
+            phantomPhysicalVolume -> SetTranslation(phantomPosition);
+    }
+    else   ConstructPhantom();
+    // Get the center of the detector
+    SetDetectorPosition();
+    if (detector)
+    {
+                        detector -> SetXHalfLength(detectorSizeX/2);
+                        detector -> SetYHalfLength(detectorSizeY/2);
+                        detector -> SetZHalfLength(detectorSizeZ/2);
+                        detectorPhysicalVolume -> SetTranslation(detectorPosition);
+    }
+    else    ConstructDetector();
+// Round to nearest integer number of voxel
+        numberOfVoxelsAlongX = lrint(detectorSizeX / sizeOfVoxelAlongX);
+        sizeOfVoxelAlongX = ( detectorSizeX / numberOfVoxelsAlongX );
+        numberOfVoxelsAlongY = lrint(detectorSizeY / sizeOfVoxelAlongY);
+        sizeOfVoxelAlongY = ( detectorSizeY / numberOfVoxelsAlongY );
+        numberOfVoxelsAlongZ = lrint(detectorSizeZ / sizeOfVoxelAlongZ);
+        sizeOfVoxelAlongZ = ( detectorSizeZ / numberOfVoxelsAlongZ );
+    //G4cout << "*************** DetectorToWorldPosition " << GetDetectorToWorldPosition()/cm << "\n";
+    ConstructSensitiveDetector(GetDetectorToWorldPosition());
+    volumeOfVoxel = sizeOfVoxelAlongX * sizeOfVoxelAlongY * sizeOfVoxelAlongZ;
+    massOfVoxel = detectorMaterial -> GetDensity() * volumeOfVoxel;
+    //  This will clear the existing matrix (together with all data inside it)!
+    matrix = HadrontherapyMatrix::GetInstance(numberOfVoxelsAlongX,
+                                              numberOfVoxelsAlongY,
+                                              numberOfVoxelsAlongZ,
+                                              massOfVoxel);
+    // Inform the kernel about the new geometry
+    G4RunManager::GetRunManager() -> GeometryHasBeenModified();
+    G4RunManager::GetRunManager() -> PhysicsHasBeenModified();
+    PrintParameters();
+}
+void HadrontherapyDetectorConstruction::PrintParameters()
+{
+    G4cout << "The (X,Y,Z) dimensions of the phantom are : (" <<
+        G4BestUnit( phantom -> GetXHalfLength()*2., "Length") << ',' <<
+        G4BestUnit( phantom -> GetYHalfLength()*2., "Length") << ',' <<
+        G4BestUnit( phantom -> GetZHalfLength()*2., "Length") << ')' << G4endl;
+    G4cout << "The (X,Y,Z) dimensions of the detector are : (" <<
+        G4BestUnit( detector -> GetXHalfLength()*2., "Length") << ',' <<
+        G4BestUnit( detector -> GetYHalfLength()*2., "Length") << ',' <<
+        G4BestUnit( detector -> GetZHalfLength()*2., "Length") << ')' << G4endl;
+    G4cout << "Displacement between Phantom and World is: ";
+    G4cout << "DX= "<< G4BestUnit(phantomPosition.getX(),"Length") <<
+        "DY= "<< G4BestUnit(phantomPosition.getY(),"Length") <<
+        "DZ= "<< G4BestUnit(phantomPosition.getZ(),"Length") << G4endl;
+    G4cout << "The (X,Y,Z) sizes of the Voxels are: (" <<
+        G4BestUnit(sizeOfVoxelAlongX, "Length")  << ',' <<
+        G4BestUnit(sizeOfVoxelAlongY, "Length")  << ',' <<
+        G4BestUnit(sizeOfVoxelAlongZ, "Length") << ')' << G4endl;
     G4cout << "The number of Voxels along (X,Y,Z) is: (" <<
+                numberOfVoxelsAlongX  << ", " <<
+                numberOfVoxelsAlongY  << ", "  <<
+                numberOfVoxelsAlongZ  << ')' << G4endl;
+    //  This will clear the existing matrix (together with data inside it)!
+    matrix = HadrontherapyMatrix::getInstance(numberOfVoxelsAlongX,
+                                              numberOfVoxelsAlongY,
+                                              numberOfVoxelsAlongZ);
+    // Here construct the Sensitive Detector and Read Out Geometry
+    ConstructSensitiveDetector(GetDetectorToWorldPosition());
+    G4RunManager::GetRunManager() -> GeometryHasBeenModified();
+    return true;
+}
+/////////////////////////////////////////////////////////////////////////////
+G4bool HadrontherapyDetectorConstruction::SetDetectorSize(G4double sizeX, G4double sizeY, G4double sizeZ)
+{
+// Check that the detector stay inside the phantom
+        if (sizeX > 0 && sizeX < sizeOfVoxelAlongX) {G4cout << "WARNING: Detector X size must be bigger than that of Voxel X" << G4endl; return false;}
+        if (sizeY > 0 && sizeY < sizeOfVoxelAlongY) {G4cout << "WARNING: Detector Y size must be bigger than that of Voxel Y" << G4endl; return false;}
+        if (sizeZ > 0 && sizeZ < sizeOfVoxelAlongZ) {G4cout << "WARNING: Detector Z size must be bigger than that of Voxel Z" << G4endl; return false;}
+        if (!IsInside(sizeX/2,
+                      sizeY/2,
+                      sizeZ/2,
+                      phantomSizeX,
+                      phantomSizeY,
+                      phantomSizeZ,
+                      detectorToPhantomPosition))
+        {return false;}
+// Negative or null values mean don't change it!
+    if (sizeX > 0) {
+                        detectorSizeX = sizeX/2;
+                        detector -> SetXHalfLength(detectorSizeX);
+                   }
+    if (sizeY > 0) {
+                        detectorSizeY = sizeY/2;
+                        detector -> SetYHalfLength(detectorSizeY);
+                   }
+    if (sizeZ > 0) {
+                        detectorSizeZ = sizeZ/2;
+                        detector -> SetZHalfLength(detectorSizeZ);
+                    }
+    G4cout << "The (X, Y, Z) dimensions of the detector are : (" <<
+                  G4BestUnit( detector -> GetXHalfLength()*2., "Length") << ", " <<
+                  G4BestUnit( detector -> GetYHalfLength()*2., "Length") << ", " <<
+                  G4BestUnit( detector -> GetZHalfLength()*2., "Length") << ')' << G4endl;
+// Adjust detector position
+    SetDetectorPosition();
+// Adjust voxels number accordingly to new detector geometry
+// Matrix will be re-instantiated!
+// Voxels and ROGeometry must follow the detector!
+    SetNumberOfVoxelBySize(sizeOfVoxelAlongX, sizeOfVoxelAlongY, sizeOfVoxelAlongZ);
+    G4RunManager::GetRunManager() -> GeometryHasBeenModified();
+    return true;
+}
+/////////////////////////////////////////////////////////////////////////////
+G4bool HadrontherapyDetectorConstruction::SetPhantomSize(G4double sizeX, G4double sizeY, G4double sizeZ)
+{
+        if (!IsInside(detectorSizeX,
+                                  detectorSizeY,
+                                  detectorSizeZ,
+                                  sizeX/2,//method parameters
+                                  sizeY/2,
+                                  sizeZ/2,
+                                  detectorToPhantomPosition
+                                  ))
+        return false;
+// Only change positive dimensions
+    if (sizeX > 0) {
+                     phantomSizeX = sizeX/2;
+                     phantom -> SetXHalfLength(phantomSizeX);
+                   }
+    if (sizeY > 0) {
+                     phantomSizeY = sizeY/2;
+                     phantom -> SetYHalfLength(phantomSizeY);
+                   }
+    if (sizeZ > 0) {
+                     phantomSizeZ = sizeZ/2;
+                     phantom -> SetZHalfLength(phantomSizeZ);
+                   }
+    G4cout << "The (X, Y, Z) dimensions of the phantom are : (" <<
+                  G4BestUnit( phantom -> GetXHalfLength()*2., "Length") << ", " <<
+                  G4BestUnit( phantom -> GetYHalfLength()*2., "Length") << ", " <<
+                  G4BestUnit( phantom -> GetZHalfLength()*2., "Length") << ')' << G4endl;
+//G4cout << '\n' << "Coordinate volume: " << phantomPhysicalVolume -> GetTranslation() << G4endl;
+// Adjust detector position inside phantom
+    SetDetectorPosition();
+    ConstructSensitiveDetector(GetDetectorToWorldPosition());
+    G4RunManager::GetRunManager() -> GeometryHasBeenModified();
+    return true;
+}
+/////////////////////////////////////////////////////////////////////////////
+G4bool HadrontherapyDetectorConstruction::SetPhantomPosition(G4ThreeVector displacement)
+{
+// Set Phantom position respect to the World
+// TODO check for overlap!
+    phantomPosition = displacement;
+    if (phantomPhysicalVolume)
+        {
+            phantomPhysicalVolume -> SetTranslation(phantomPosition);
+            G4cout << "Displacement between Phantom and World is: ";
+            G4cout << "DX= "<< G4BestUnit(phantomPosition.getX(),"Length") << ", " <<
+                      "DY= "<< G4BestUnit(phantomPosition.getY(),"Length") << ", " <<
+                      "DZ= "<< G4BestUnit(phantomPosition.getZ(),"Length") << G4endl;
+// Redraw ROGeometry!
+            ConstructSensitiveDetector(GetDetectorToWorldPosition());
+            G4RunManager::GetRunManager() -> GeometryHasBeenModified();
+        }
+    return true;
+}
+/////////////////////////////////////////////////////////////////////////////
+G4bool HadrontherapyDetectorConstruction::SetDetectorToPhantomPosition(G4ThreeVector displacement)
+{
+// Ignore negative values
+    if (displacement[0] < 0.) displacement[0] = detectorToPhantomPosition[0];
+    if (displacement[1] < 0.) displacement[1] = detectorToPhantomPosition[1];
+    if (displacement[2] < 0.) displacement[2] = detectorToPhantomPosition[2];
+    if (!IsInside(detectorSizeX,
+                  detectorSizeY,
+                  detectorSizeZ,
+                  phantomSizeX,
+                  phantomSizeY,
+                  phantomSizeZ,
+                  displacement // method parameter!
+                  ))
+    {return false;}
+    detectorToPhantomPosition = displacement;
+// Adjust detector position inside phantom
+    SetDetectorPosition();
+    ConstructSensitiveDetector(GetDetectorToWorldPosition());
+    G4RunManager::GetRunManager() -> GeometryHasBeenModified();
+    return true;
+}
+        numberOfVoxelsAlongX  << ',' <<
+        numberOfVoxelsAlongY  <<','  <<
+        numberOfVoxelsAlongZ  << ')' << G4endl;
+}

trunk/examples/advanced/hadrontherapy/src/HadrontherapyDetectorMessenger.cc

-                      r1230
+                      r1313
 #include "HadrontherapyDetectorConstruction.hh"
 #include "G4UIdirectory.hh"
+#include "G4UIcmdWithADoubleAndUnit.hh"
+#include "G4UIcmdWith3VectorAndUnit.hh"
+#include "G4UIcmdWithoutParameter.hh"
 #include "G4UIcmdWithAString.hh"
-#include "G4UIcmdWith3VectorAndUnit.hh"
 /////////////////////////////////////////////////////////////////////////////
 …
                                                 "PhantomSizeAlongZ", false);
     changeThePhantomSizeCmd -> SetDefaultUnit("mm");
     changeThePhantomSizeCmd -> SetUnitCandidates("um mm cm");
+    changeThePhantomSizeCmd -> SetUnitCandidates("nm um mm cm");
     changeThePhantomSizeCmd -> AvailableForStates(G4State_Idle);
+    // Change Phantom material
+    changeThePhantomMaterialCmd = new G4UIcmdWithAString("/changePhantom/material", this);
+    changeThePhantomMaterialCmd -> SetGuidance("Change the Phantom and the detector material");
+    changeThePhantomMaterialCmd -> SetParameterName("PhantomMaterial", false);
+    changeThePhantomMaterialCmd -> SetDefaultValue("G4_WATER");
+    changeThePhantomMaterialCmd -> AvailableForStates(G4State_Idle);
     // Change Phantom position
 …
                                                     "PositionAlongZ", false);
     changeThePhantomPositionCmd -> SetDefaultUnit("mm");
     changeThePhantomPositionCmd -> SetUnitCandidates("mm cm m");
+    changeThePhantomPositionCmd -> SetUnitCandidates("um mm cm m");
     changeThePhantomPositionCmd -> AvailableForStates(G4State_Idle);
+    updateCmd = new G4UIcmdWithoutParameter("/changePhantom/update",this);
+    updateCmd->SetGuidance("Update Phantom/Detector geometry.");
+    updateCmd->SetGuidance("This command MUST be applied before \"beamOn\" ");
+    updateCmd->SetGuidance("if you changed geometrical value(s).");
+    updateCmd->AvailableForStates(G4State_Idle);
     //  Change detector size
 …
     changeTheDetectorSizeCmd -> SetParameterName("DetectorSizeAlongX", "DetectorSizeAlongY", "DetectorSizeAlongZ", false);
     changeTheDetectorSizeCmd -> SetDefaultUnit("mm");
     changeTheDetectorSizeCmd -> SetUnitCandidates("um mm cm");
+    changeTheDetectorSizeCmd -> SetUnitCandidates("nm um mm cm");
     changeTheDetectorSizeCmd -> AvailableForStates(G4State_Idle);
 …
                                                               "DisplacementAlongZ", false);
     changeTheDetectorToPhantomPositionCmd -> SetDefaultUnit("mm");
     changeTheDetectorToPhantomPositionCmd -> SetUnitCandidates("um mm cm");
+    changeTheDetectorToPhantomPositionCmd -> SetUnitCandidates("nm um mm cm");
     changeTheDetectorToPhantomPositionCmd -> AvailableForStates(G4State_Idle);
 …
     changeTheDetectorVoxelCmd -> SetParameterName("VoxelSizeAlongX", "VoxelSizeAlongY", "VoxelSizeAlongZ", false);
     changeTheDetectorVoxelCmd -> SetDefaultUnit("mm");
     changeTheDetectorVoxelCmd -> SetUnitCandidates("um mm cm");
+    changeTheDetectorVoxelCmd -> SetUnitCandidates("nm um mm cm");
     changeTheDetectorVoxelCmd -> AvailableForStates(G4State_Idle);
+   }
 …
     delete changeThePhantomSizeCmd;
     delete changeThePhantomPositionCmd;
+    delete changeThePhantomMaterialCmd;
+    delete updateCmd;
     delete changeTheDetectorDir;
     delete changeTheDetectorSizeCmd;
 …
          hadrontherapyDetector -> SetPhantomPosition(size);
+  }
+  else if (command == changeThePhantomMaterialCmd)
+  {
+      hadrontherapyDetector -> SetPhantomMaterial(newValue);
+  }
   else if (command == changeTheDetectorSizeCmd)
+  {
 …
+  {
         G4ThreeVector size = changeTheDetectorVoxelCmd  -> GetNew3VectorValue(newValue);
+        hadrontherapyDetector -> SetNumberOfVoxelBySize(size.getX(),size.getY(),size.getZ());
+        hadrontherapyDetector -> SetVoxelSize(size.getX(),size.getY(),size.getZ());
+  }
+  else if (command == updateCmd)
+  {
+      hadrontherapyDetector -> UpdateGeometry();
+  }
+}

trunk/examples/advanced/hadrontherapy/src/HadrontherapyEventAction.cc

-                      r1230
+                      r1313
 // ********************************************************************
 //
+// $Id: HadrontherapyEventAction.cc;
+// See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
+// $Id: HadrontherapyEventAction.cc;
+//
+// See more at: http://workgroup.lngs.infn.it/geant4lns
+//
+// ----------------------------------------------------------------------------
+//                 GEANT 4 - Hadrontherapy example
+// ----------------------------------------------------------------------------
+// Code developed by:
+//
+// G.A.P. Cirrone(a)*, G. Candiano, F. Di Rosa(a), S. Guatelli(b), G. Russo(a)
+//
+// (a) Laboratori Nazionali del Sud
+//     of the National Institute for Nuclear Physics, Catania, Italy
+// (b) National Institute for Nuclear Physics Section of Genova, genova, Italy
+//
+// * cirrone@lns.infn.it
+// --------------------------------------------------------------
 #include "G4Event.hh"
 #include "G4EventManager.hh"
 #include "G4HCofThisEvent.hh"
 #include "G4VHitsCollection.hh"
-#include "G4TrajectoryContainer.hh"
-#include "G4Trajectory.hh"
-#include "G4VVisManager.hh"
 #include "G4SDManager.hh"
 #include "G4VVisManager.hh"
 #include "HadrontherapyEventAction.hh"
 #include "HadrontherapyDetectorHit.hh"
 …
   //printing survey
   if (evtNb%printModulo == 0)
     G4cout << "\n---> Begin of Event: " << evtNb << G4endl;
+     G4cout << "\n---> Begin of Event: " << evtNb << G4endl;
   G4SDManager* pSDManager = G4SDManager::GetSDMpointer();
   if(hitsCollectionID == -1)
     hitsCollectionID = pSDManager -> GetCollectionID("HadrontherapyDetectorHitsCollection");
+}
 /////////////////////////////////////////////////////////////////////////////
 void HadrontherapyEventAction::EndOfEventAction(const G4Event* evt)
+{
+{
   if(hitsCollectionID < 0)
   return;
   G4HCofThisEvent* HCE = evt -> GetHCofThisEvent();
+  // Clear voxels hit list
+  HadrontherapyMatrix* matrix = HadrontherapyMatrix::GetInstance();
+  if (matrix) matrix -> ClearHitTrack();
   if(HCE)
+  {
 …
     if(CHC)
+     {
-           matrix = HadrontherapyMatrix::getInstance();
        if(matrix)
+          {
 …
+    }
+  }
-  // Extract the trajectories and draw them in the visualisation
-  if (G4VVisManager::GetConcreteInstance())
+    {
-      G4TrajectoryContainer * trajectoryContainer = evt -> GetTrajectoryContainer();
-      G4int n_trajectories = 0;
-      if (trajectoryContainer) n_trajectories = trajectoryContainer -> entries();
-      for (G4int i = 0; i < n_trajectories; i++)
+        {
-          G4Trajectory* trj = (G4Trajectory*)
-            ((*(evt -> GetTrajectoryContainer()))[i]);
-          if(drawFlag == "all") trj -> DrawTrajectory(50);
-          else if((drawFlag == "charged")&&(trj -> GetCharge() != 0.))
-            trj -> DrawTrajectory(50);
-          else if ((drawFlag == "neutral")&&(trj -> GetCharge() == 0.))
-            trj -> DrawTrajectory(50);
+        }
+    }
+}

trunk/examples/advanced/hadrontherapy/src/HadrontherapyMatrix.cc

-                      r1230
+                      r1313
 //
 // ********************************************************************
 // * License and Disclaimer                                           *
 // *                                                                  *
 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
 // * conditions of the Geant4 Software License,  included in the file *
 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
 // * include a list of copyright holders.                             *
 // *                                                                  *
 // * Neither the authors of this software system, nor their employing *
 // * institutes,nor the agencies providing financial support for this *
 // * work  make  any representation or  warranty, express or implied, *
 // * regarding  this  software system or assume any liability for its *
 // * use.  Please see the license in the file  LICENSE  and URL above *
 // * for the full disclaimer and the limitation of liability.         *
 // *                                                                  *
 // * This  code  implementation is the result of  the  scientific and *
 // * technical work of the GEANT4 collaboration.                      *
 // * By using,  copying,  modifying or  distributing the software (or *
 // * any work based  on the software)  you  agree  to acknowledge its *
 // * use  in  resulting  scientific  publications,  and indicate your *
 // * acceptance of all terms of the Geant4 Software license.          *
 // ********************************************************************
 //
 // $Id: HadrontherapyMatrix.cc;
 // See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy//
+        //
+        // ********************************************************************
+        // * License and Disclaimer                                           *
+        // *                                                                  *
+        // * The  Geant4 software  is  copyright of the Copyright Holders  of *
+        // * the Geant4 Collaboration.  It is provided  under  the terms  and *
+        // * conditions of the Geant4 Software License,  included in the file*
+        // * LICENSE and available at  http://cern.ch/geant4/license .  These *
+        // * include a list of copyright holders.                             *
+        // *                                                                  *
+        // * Neither the authors of this software system, nor their employing *
+        // * institutes,nor the agencies providing financial support for this *
+        // * work  make  any representation or  warranty, express or implied, *
+        // * regarding  this  software system or assume any liability for its *
+        // * use.  Please see the license in the file  LICENSE  and URL above *
+        // * for the full disclaimer and the limitation of liability.         *
+        // *                                                                  *
+        // * This  code  implementation is the result of  the  scientific and *
+        // * technical work of the GEANT4 collaboration.                      *
+        // * By using,  copying,  modifying or  distributing the software (or *
+        // * any work based  on the software)  you  agree  to acknowledge its *
+        // * use  in  resulting  scientific  publications,  and indicate your *
+        // * acceptance of all terms of the Geant4 Software license.          *
+        // ********************************************************************
+        //
+        // $Id: HadrontherapyMatrix.cc;
+        // See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy//
 #include "HadrontherapyMatrix.hh"
 #include "HadrontherapyAnalysisManager.hh"
+#include "G4RunManager.hh"
+#include "HadrontherapyPrimaryGeneratorAction.hh"
+#include "G4ParticleGun.hh"
+#include <fstream>
+#include <unistd.h>
+#include <iostream>
+#include <sstream>
+#include <iomanip>
 #include "globals.hh"
+#include <fstream>
+// Units definition: CLHEP/Units/SystemOfUnits.h
+//
 HadrontherapyMatrix* HadrontherapyMatrix::instance = NULL;
+// Static method that only return a pointer to the matrix object
+HadrontherapyMatrix* HadrontherapyMatrix::getInstance()
+{
+  return instance;
+}
+// This STATIC method delete (!) the old matrix and rewrite a new object returning a pointer to it
+HadrontherapyMatrix* HadrontherapyMatrix::getInstance(G4int voxelX, G4int voxelY, G4int voxelZ)
+{
+    if (instance) delete instance;
+    instance = new HadrontherapyMatrix(voxelX, voxelY, voxelZ);
+    instance -> Initialize();
+    return instance;
+}
+HadrontherapyMatrix::HadrontherapyMatrix(G4int voxelX, G4int voxelY, G4int voxelZ):
+    matrix(0)
+G4bool HadrontherapyMatrix::secondaries = false;
+        // Only return a pointer to matrix
+HadrontherapyMatrix* HadrontherapyMatrix::GetInstance()
+{
+        return instance;
+}
+        // This STATIC method delete (!) the old matrix and rewrite a new object returning a pointer to it
+        // TODO A check on the parameters is required!
+HadrontherapyMatrix* HadrontherapyMatrix::GetInstance(G4int voxelX, G4int voxelY, G4int voxelZ, G4double mass)
+{
+        if (instance) delete instance;
+        instance = new HadrontherapyMatrix(voxelX, voxelY, voxelZ, mass);
+        instance -> Initialize();// XXX
+        return instance;
+}
+HadrontherapyMatrix::HadrontherapyMatrix(G4int voxelX, G4int voxelY, G4int voxelZ, G4double mass):
+    filename("Dose.out"),
+    doseUnit(MeV/g)
+{
+// Number of the voxels of the phantom
+// For Y = Z = 1 the phantom is divided in slices (and not in voxels)
+// orthogonal to the beam axis
+  numberVoxelX = voxelX;
+  numberVoxelY = voxelY;
+  numberVoxelZ = voxelZ;
+  // Create the dose matrix
+  matrix = new G4double[numberVoxelX*numberVoxelY*numberVoxelZ];
+  if (matrix) G4cout << "Matrix: Memory space to store physical dose into " <<
+                        numberVoxelX*numberVoxelY*numberVoxelZ <<
+                        " voxels has been allocated " << G4endl;
+  else G4Exception("Can't allocate memory to store physical dose!");
+}
+        // Number of the voxels of the phantom
+        // For Y = Z = 1 the phantom is divided in slices (and not in voxels)
+        // orthogonal to the beam axis
+        numberOfVoxelAlongX = voxelX;
+        numberOfVoxelAlongY = voxelY;
+        numberOfVoxelAlongZ = voxelZ;
+        massOfVoxel = mass;
+        // Create the dose matrix
+        matrix = new G4double[numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ];
+        if (matrix)
+        {
+                G4cout << "HadrontherapyMatrix: Memory space to store physical dose into " <<
+                numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ <<
+                " voxels has been allocated " << G4endl;
+        }
+        else G4Exception(" HadrontherapyMatrix::HadrontherapyMatrix. Can't allocate memory to store physical dose!");
+                // Hit voxel (TrackID) marker
+                // This array mark the status of voxel, if a hit occur, with the trackID of the particle
+                // Must be initialized
+        hitTrack = new G4int[numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ];
+        ClearHitTrack();
+}
+/////////////////////////////////////////////////////////////////////////////
 HadrontherapyMatrix::~HadrontherapyMatrix()
+{
     delete[] matrix;
+}
+void HadrontherapyMatrix::flush()
+{
+        if(matrix)
+        for(int i=0; i<numberVoxelX*numberVoxelY*numberVoxelZ; i++)
+        {
+          matrix[i] = 0;
+        }
+}
+    delete[] hitTrack;
+                // clear fluences/dose data
+    Clear();
+}
+/////////////////////////////////////////////////////////////////////////////
+void HadrontherapyMatrix::Clear()
+{
+    for (size_t i=0; i<ionStore.size(); i++)
+    {
+        delete[] ionStore[i].dose;
+        delete[] ionStore[i].fluence;
+    }
+    ionStore.clear();
+}
+/////////////////////////////////////////////////////////////////////////////
+// Initialise the elements of the matrix to zero
 void HadrontherapyMatrix::Initialize()
+{
+// Initialise the elements of the matrix to zero
+  for(G4int i = 0; i < numberVoxelX; i++)
+    {
+      for(G4int j = 0; j < numberVoxelY; j++)
+           {
+              for(G4int k = 0; k < numberVoxelZ; k++)
+              matrix[Index(i,j,k)] = 0.;
+           }
+    }
+}
+    Clear();
+    for(int i=0;i<numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ;i++)
+    {
+                matrix[i] = 0;
+    }
+}
+        /////////////////////////////////////////////////////////////////////////////
+        /////////////////////////////////////////////////////////////////////////////
+        // Print generated nuclides list
+void HadrontherapyMatrix::PrintNuclides()
+{
+    for (size_t i=0; i<ionStore.size(); i++)
+    {
+                G4cout << ionStore[i].name << G4endl;
+    }
+}
+        /////////////////////////////////////////////////////////////////////////////
+        // Clear Hit voxel (TrackID) markers
+void HadrontherapyMatrix::ClearHitTrack()
+{
+        for(G4int i=0; i<numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ; i++) hitTrack[i] = 0;
+}
+        // Return Hit status
+G4int* HadrontherapyMatrix::GetHitTrack(G4int i, G4int j, G4int k)
+{
+        return &(hitTrack[Index(i,j,k)]);
+}
+/////////////////////////////////////////////////////////////////////////////
+// Dose methods...
+// Fill DOSE/fluence matrix for particle/ion:
+// If fluence parameter is true (default value is FALSE) then fluence at voxel (i, j, k) is increased.
+// The energyDeposit parameter fill the dose matrix for voxel (i,j,k)
+/////////////////////////////////////////////////////////////////////////////
+G4bool HadrontherapyMatrix::Fill(G4int trackID,
+                                 G4ParticleDefinition* particleDef,
+                                 G4int i, G4int j, G4int k,
+                                 G4double energyDeposit,
+                                 G4bool fluence)
+{
+    if (energyDeposit <=0. && !fluence || !secondaries) return false;
+    // Get Particle Data Group particle ID
+    G4int PDGencoding = particleDef -> GetPDGEncoding();
+    PDGencoding -= PDGencoding%10;
+    // Search for already allocated data...
+    for (size_t l=0; l < ionStore.size(); l++)
+    {
+                if (ionStore[l].PDGencoding == PDGencoding )
+                {   // Is it a primary or a secondary particle?
+                        if ( trackID ==1 && ionStore[l].isPrimary || trackID !=1 && !ionStore[l].isPrimary)
+                        {
+                                if (energyDeposit > 0.) ionStore[l].dose[Index(i, j, k)] += energyDeposit/massOfVoxel;
+                                        // Fill a matrix per each ion with the fluence
+                                if (fluence) ionStore[l].fluence[Index(i, j, k)]++;
+                                return true;
+                        }
+                }
+    }
+    G4int Z = particleDef-> GetAtomicNumber();
+    G4int A = particleDef-> GetAtomicMass();
+    G4String fullName = particleDef -> GetParticleName();
+    G4String name = fullName.substr (0, fullName.find("[") ); // cut excitation energy
+  // Let's put a new particle in our store...
+    ion newIon =
+    {
+                (trackID == 1) ? true:false,
+                PDGencoding,
+                name,
+                name.length(),
+                Z,
+                A,
+                new G4double[numberOfVoxelAlongX * numberOfVoxelAlongY * numberOfVoxelAlongZ],
+                new unsigned int[numberOfVoxelAlongX * numberOfVoxelAlongY * numberOfVoxelAlongZ]
+    };
+                // Initialize data
+    if (newIon.dose && newIon.fluence)
+    {
+                for(G4int m=0; m<numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ; m++)
+                {
+                        newIon.dose[m] = 0.;
+                        newIon.fluence[m] = 0;
+                }
+                if (energyDeposit > 0.) newIon.dose[Index(i, j, k)] += energyDeposit/massOfVoxel;
+                if (fluence) newIon.fluence[Index(i, j, k)]++;
+                ionStore.push_back(newIon);
+                // TODO Put some verbosity check
+                /*
+                 G4cout << "Memory space to store the DOSE/FLUENCE into " <<
+                 numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ <<
+                 " voxels has been allocated for the nuclide " << newIon.name <<
+                 " (Z = " << Z << ", A = " << A << ")" << G4endl ;
+                 */
+                return true;
+    }
+    else // XXX Out of memory! XXX
+    {
+                return false;
+    }
+}
+        /////////////////////////////////////////////////////////////////////////////
+        /////////////////////////////////////////////////////////////////////////////
+        // Methods to store data to filenames...
+        ////////////////////////////////////////////////////////////////////////////
+        ////////////////////////////////////////////////////////////////////////////
+        //
+        // General method to store matrix data to filename
+void HadrontherapyMatrix::StoreMatrix(G4String file, void* data, size_t psize)
+{
+    if (data)
+    {
+                ofs.open(file, std::ios::out);
+                if (ofs.is_open())
+                {
+                        for(G4int i = 0; i < numberOfVoxelAlongX; i++)
+                                for(G4int j = 0; j < numberOfVoxelAlongY; j++)
+                                        for(G4int k = 0; k < numberOfVoxelAlongZ; k++)
+                                        {
+                                                G4int n = Index(i, j, k);
+                                                        // Check for data type: u_int, G4double, XXX
+                                                if (psize == sizeof(unsigned int))
+                                                {
+                                                        unsigned int* pdata = (unsigned int*)data;
+                                                        if (pdata[n]) ofs << i << '\t' << j << '\t' <<
+                                                                k << '\t' << pdata[n] << G4endl;
+                                                }
+                                                else if (psize == sizeof(G4double))
+                                                {
+                                                        G4double* pdata = (G4double*)data;
+                                                        if (pdata[n]) ofs << i << '\t' << j << '\t' <<
+                                                                k << '\t' << pdata[n] << G4endl;
+                                                }
+                                        }
+                        ofs.close();
+                }
+    }
+}
+        // Store fluence per single ion in multiple files
+void HadrontherapyMatrix::StoreFluenceData()
+{
+    for (size_t i=0; i < ionStore.size(); i++){
+                StoreMatrix(ionStore[i].name + "_Fluence.out", ionStore[i].fluence, sizeof(unsigned int));
+    }
+}
+        // Store dose per single ion in multiple files
+void HadrontherapyMatrix::StoreDoseData()
+{
+    for (size_t i=0; i < ionStore.size(); i++){
+                StoreMatrix(ionStore[i].name + "_Dose.out", ionStore[i].dose, sizeof(G4double));
+    }
+}
+        /////////////////////////////////////////////////////////////////////////
+        // Store dose for all ions into a single file and into ntuples.
+        // Please note that this function is called via messenger commands
+        // defined in the HadrontherapyAnalysisFileMessenger.cc class file
+void HadrontherapyMatrix::StoreDoseFluenceAscii()
+{
+    // Sort like periodic table
+    std::sort(ionStore.begin(), ionStore.end());
+#define width 15L
+    ofs.open(filename, std::ios::out);
+    if (ofs.is_open())
+    {
+        //
+        // Write the voxels index and the list of particles/ions
+        ofs << std::setprecision(6) << std::left <<
+            "i\tj\tk\t";
+/*
+            G4RunManager *runManager = G4RunManager::GetRunManager();
+            HadrontherapyPrimaryGeneratorAction *pPGA = (HadrontherapyPrimaryGeneratorAction*)runManager -> GetUserPrimaryGeneratorAction();
+            G4String name = pPGA -> GetParticleGun() -> GetParticleDefinition() -> GetParticleName();
+*/
+        // Total dose
+        ofs << std::setw(width) << "Dose";
+        for (size_t l=0; l < ionStore.size(); l++)
+        {
+            G4String a = (ionStore[l].isPrimary) ? "_1":""; // is it a primary?
+            ofs << std::setw(width) << ionStore[l].name + a <<
+                std::setw(width) << ionStore[l].name  + a;
+        }
+        ofs << G4endl;
+/*
+        // PDGencoding
+        ofs << std::setprecision(6) << std::left <<
+        "0\t0\t0\t";
+        // Total dose
+        ofs << std::setw(width) << '0';
+        for (size_t l=0; l < ionStore.size(); l++)
+        {
+        ofs << std::setw(width) << ionStore[l].PDGencoding  <<
+        std::setw(width) << ionStore[l].PDGencoding;
+        }
+        ofs << G4endl;
+*/
+        // Write data
+        for(G4int i = 0; i < numberOfVoxelAlongX; i++)
+            for(G4int j = 0; j < numberOfVoxelAlongY; j++)
+                for(G4int k = 0; k < numberOfVoxelAlongZ; k++)
+                {
+                    G4int n = Index(i, j, k);
+                    // Write only not identically null data lines
+                    if (matrix[n])
+                    {
+                        ofs << G4endl;
+                        ofs << i << '\t' << j << '\t' << k << '\t';
+                        // Total dose
+                        ofs << std::setw(width) << matrix[n]/(doseUnit);
+                        {
+                            for (size_t l=0; l < ionStore.size(); l++)
+                            {
+                                // Fill ASCII file rows
+                                ofs << std::setw(width) << ionStore[l].dose[n]/(doseUnit) <<
+                                    std::setw(width) << ionStore[l].fluence[n];
+                            }
+                        }
+                    }
+                }
+        ofs.close();
+    }
+}
+/////////////////////////////////////////////////////////////////////////////
+#ifdef G4ANALYSIS_USE_ROOT
+void HadrontherapyMatrix::StoreDoseFluenceRoot()
+{
+    HadrontherapyAnalysisManager* analysis = HadrontherapyAnalysisManager::GetInstance();
+    for(G4int i = 0; i < numberOfVoxelAlongX; i++)
+        for(G4int j = 0; j < numberOfVoxelAlongY; j++)
+            for(G4int k = 0; k < numberOfVoxelAlongZ; k++)
+            {
+                G4int n = Index(i, j, k);
+                for (size_t l=0; l < ionStore.size(); l++)
+                {
+                    // Do the same work for .root file: fill dose/fluence ntuple
+                    analysis -> FillVoxelFragmentTuple( i, j, k,
+                                                        ionStore[l].A,
+                                                        ionStore[l].Z,
+                                                        ionStore[l].dose[n]/(doseUnit),
+                                                        ionStore[l].fluence[n] );
+                }
+            }
+}
+#endif
 void HadrontherapyMatrix::Fill(G4int i, G4int j, G4int k,
                                G4double energyDeposit)
+{
+  if (matrix)
+    matrix[Index(i,j,k)] += energyDeposit;
+  // Store the energy deposit in the matrix element corresponding
+  // to the phantom voxel  i, j, k
+}
+    if (matrix)
+                matrix[Index(i,j,k)] += energyDeposit;
+    // Store the energy deposit in the matrix element corresponding
+    // to the phantom voxel
+}
 void HadrontherapyMatrix::TotalEnergyDeposit()
+{
+  // Convert energy deposited to dose.
   // Store the information of the matrix in a ntuple and in
   // a 1D Histogram
+  G4int k;
+  G4int j;
+  G4int i;
+  if (matrix)
+#ifdef G4ANALYSIS_USE_ROOT
+    HadrontherapyAnalysisManager* analysis = HadrontherapyAnalysisManager::GetInstance();
+#endif
+    if (matrix)
+    {
+      std::ofstream ofs;
+      ofs.open("DoseDistribution.out");
+      for(G4int l = 0; l < numberVoxelZ; l++)
+        {
+          k = l;
+          for(G4int m = 0; m < numberVoxelY; m++)
+            {
+              j = m * numberVoxelZ + k;
+                for(G4int n = 0; n <  numberVoxelX; n++)
+                  {
+                    i =  n* numberVoxelZ * numberVoxelY + j;
+                    if(matrix[i] != 0)
+                      {
+                        ofs << n << '\t' << m << '\t' <<
+                          k << '\t' << matrix[i] << G4endl;
+#ifdef ANALYSIS_USE
+                        HadrontherapyAnalysisManager* analysis =
+                        HadrontherapyAnalysisManager::getInstance();
+                        analysis -> FillEnergyDeposit(n, m, k, matrix[i]);
+                        analysis -> BraggPeak(n, matrix[i]);
+        for(G4int i = 0; i < numberOfVoxelAlongX; i++)
+            for(G4int j = 0; j < numberOfVoxelAlongY; j++)
+                for(G4int k = 0; k < numberOfVoxelAlongZ; k++)
+                {
+                    G4int n = Index(i,j,k);
+                    matrix[n] = matrix[n] / massOfVoxel;
+#ifdef G4ANALYSIS_USE_ROOT
+                    analysis -> FillEnergyDeposit(i, j, k, matrix[n]/doseUnit);
+                    analysis -> BraggPeak(i, matrix[n]/doseUnit);
 #endif
+                      }
+}
+              }
+          }
+    }
+}
+                }
+    }
+}

trunk/examples/advanced/hadrontherapy/src/HadrontherapyPhysicsList.cc

-                      r1230
+                      r1313
 //                        and local physic for ion-ion enelastic processes)
+#include "G4RunManager.hh"
+#include "G4Region.hh"
+#include "G4RegionStore.hh"
 #include "HadrontherapyPhysicsList.hh"
 #include "HadrontherapyPhysicsListMessenger.hh"
 …
 #include "LocalINCLIonIonInelasticPhysic.hh"         // Physic dedicated to the ion-ion inelastic processes using INCL/ABLA
 // Physic lists (contained inside the Geant4 distribution)
+// Physic lists (contained inside the Geant4 source code, in the 'physicslists folder')
 #include "G4EmStandardPhysics_option3.hh"
 #include "G4EmLivermorePhysics.hh"
 …
+{
   // transportation
-  //
   AddTransportation();
   // electromagnetic physics list
-  //
   emPhysicsList->ConstructProcess();
   em_config.AddModels();
   // decay physics list
-  //
   decPhysicsList->ConstructProcess();
 …
   //   ELECTROMAGNETIC MODELS
   /////////////////////////////////////////////////////////////////////////////
+  if (name == "standard_opt3") {
+    if (name == "standard_opt3") {
     emName = name;
     delete emPhysicsList;
     emPhysicsList = new G4EmStandardPhysics_option3();
+    G4RunManager::GetRunManager() -> PhysicsHasBeenModified();
     G4cout << "THE FOLLOWING ELECTROMAGNETIC PHYSICS LIST HAS BEEN ACTIVATED: G4EmStandardPhysics_option3" << G4endl;
 …
     delete emPhysicsList;
     emPhysicsList = new G4EmLivermorePhysics();
+    G4RunManager::GetRunManager()-> PhysicsHasBeenModified();
     G4cout << "THE FOLLOWING ELECTROMAGNETIC PHYSICS LIST HAS BEEN ACTIVATED: G4EmLivermorePhysics" << G4endl;
 …
   SetCutValue(cutForPositron, "e+");
+  // Set cuts for detector
+  SetDetectorCut(defaultCutValue);
   if (verboseLevel>0) DumpCutValuesTable();
+}
 …
   SetParticleCuts(cutForPositron, G4Positron::Positron());
+}
+void HadrontherapyPhysicsList::SetDetectorCut(G4double cut)
+{
+  G4String regionName = "DetectorLog";
+  G4Region* region = G4RegionStore::GetInstance()->GetRegion(regionName);
+  G4ProductionCuts* cuts = new G4ProductionCuts ;
+  cuts -> SetProductionCut(cut,G4ProductionCuts::GetIndex("gamma"));
+  cuts -> SetProductionCut(cut,G4ProductionCuts::GetIndex("e-"));
+  cuts -> SetProductionCut(cut,G4ProductionCuts::GetIndex("e+"));
+  region -> SetProductionCuts(cuts);
+}

trunk/examples/advanced/hadrontherapy/src/HadrontherapyPhysicsListMessenger.cc

-                      r1230
+                      r1313
   allCutCmd->AvailableForStates(G4State_PreInit,G4State_Idle);
+  allDetectorCmd = new G4UIcmdWithADoubleAndUnit("/physic/setDetectorCuts",this);
+  allDetectorCmd->SetGuidance("Set cut for all. into Detector");
+  allDetectorCmd->SetParameterName("cut",false);
+  allDetectorCmd->SetUnitCategory("Length");
+  allDetectorCmd->SetRange("cut>0.0");
+  allDetectorCmd->AvailableForStates(G4State_PreInit,G4State_Idle);
   pListCmd = new G4UIcmdWithAString("/physic/addPhysics",this);
   pListCmd->SetGuidance("Add physics list.");
   pListCmd->SetParameterName("PList",false);
   pListCmd->AvailableForStates(G4State_PreInit);
+  pListCmd->AvailableForStates(G4State_PreInit, G4State_Idle);
   packageListCmd = new G4UIcmdWithAString("/physic/addPackage",this);
   packageListCmd->SetGuidance("Add physics package.");
   packageListCmd->SetParameterName("package",false);
   packageListCmd->AvailableForStates(G4State_PreInit);
+  packageListCmd->AvailableForStates(G4State_PreInit, G4State_Idle);
+}
 …
   delete protoCutCmd;
   delete allCutCmd;
+  delete allDetectorCmd;
   delete pListCmd;
   delete physDir;
 …
    { pPhysicsList->SetCutForGamma(gammaCutCmd->GetNewDoubleValue(newValue));}
   if( command == electCutCmd )
+  else if( command == electCutCmd )
    { pPhysicsList->SetCutForElectron(electCutCmd->GetNewDoubleValue(newValue));}
   if( command == protoCutCmd )
+  else if( command == protoCutCmd )
    { pPhysicsList->SetCutForPositron(protoCutCmd->GetNewDoubleValue(newValue));}
   if( command == allCutCmd )
+  else if( command == allCutCmd )
+    {
       G4double cut = allCutCmd->GetNewDoubleValue(newValue);
 …
       pPhysicsList->SetCutForPositron(cut);
+    }
+  if( command == pListCmd )
+  else if( command == allDetectorCmd)
+  {
+      G4double cut = allDetectorCmd -> GetNewDoubleValue(newValue);
+      pPhysicsList -> SetDetectorCut(cut);
+  }
+  else if( command == pListCmd )
    { pPhysicsList->AddPhysicsList(newValue);}
   if( command == packageListCmd )
+  else if( command == packageListCmd )
    { pPhysicsList->AddPackage(newValue);}

trunk/examples/advanced/hadrontherapy/src/HadrontherapyPrimaryGeneratorAction.cc

-                      r1230
+                      r1313
   sigmaEnergy = defaultsigmaEnergy;
 #ifdef ANALYSIS_USE
+#ifdef G4ANALYSIS_USE_ROOT
   // Write these values into the analysis if needed. Have to be written separately on change.
   HadrontherapyAnalysisManager::getInstance()->setBeamMetaData(meanKineticEnergy, sigmaEnergy);
+  HadrontherapyAnalysisManager::GetInstance()->setBeamMetaData(meanKineticEnergy, sigmaEnergy);
 #endif
 …
 void HadrontherapyPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
+{
 #ifdef ANALYSIS_USE
+#ifdef G4ANALYSIS_USE_ROOT
   // Increment the event counter
   HadrontherapyAnalysisManager::getInstance()->startNewEvent();
+  HadrontherapyAnalysisManager::GetInstance()->startNewEvent();
 #endif
 …
+{
         meanKineticEnergy = val;
 #ifdef ANALYSIS_USE
+#ifdef G4ANALYSIS_USE_ROOT
   // Update the beam-data in the analysis manager
   HadrontherapyAnalysisManager::getInstance()->setBeamMetaData(meanKineticEnergy, sigmaEnergy);
+  HadrontherapyAnalysisManager::GetInstance()->setBeamMetaData(meanKineticEnergy, sigmaEnergy);
 #endif
 …
+{
         sigmaEnergy = val;
 #ifdef ANALYSIS_USE
+#ifdef G4ANALYSIS_USE_ROOT
   // Update the sigmaenergy in the metadata.
   HadrontherapyAnalysisManager::getInstance()->setBeamMetaData(meanKineticEnergy, sigmaEnergy);
+  HadrontherapyAnalysisManager::GetInstance()->setBeamMetaData(meanKineticEnergy, sigmaEnergy);
 #endif
+}
 …
 G4double HadrontherapyPrimaryGeneratorAction::GetmeanKineticEnergy(void)
 { return meanKineticEnergy;}

trunk/examples/advanced/hadrontherapy/src/HadrontherapyRunAction.cc

-                      r1230
+                      r1313
 // See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
 #include "HadrontherapyRunAction.hh"
 #include "HadrontherapyEventAction.hh"
 …
 #include "HadrontherapyRunAction.hh"
+#include "HadrontherapyAnalysisManager.hh"
+#include "HadrontherapyMatrix.hh"
 HadrontherapyRunAction::HadrontherapyRunAction()
+{
 …
 void HadrontherapyRunAction::BeginOfRunAction(const G4Run* aRun)
+{
    G4RunManager::GetRunManager()->SetRandomNumberStore(true);
    G4cout << "Run " << aRun -> GetRunID() << " starts ..." << G4endl;
+    G4RunManager::GetRunManager()-> SetRandomNumberStore(true);
+    G4cout << "Run " << aRun -> GetRunID() << " starts ..." << G4endl;
+   electromagnetic = 0;
+   hadronic = 0;
+    // Warning! any beamOn will reset all data (dose, fluence, histograms, etc)!
+    //
+    // Initialize matrix with energy of primaries clearing data inside
+    if (HadrontherapyMatrix::GetInstance())
+    {
+        HadrontherapyMatrix::GetInstance() -> Initialize();
+    }
+#ifdef G4ANALYSIS_USE_ROOT
+    HadrontherapyAnalysisManager::GetInstance() -> flush();     // Finalize the root file
+    // Initialize root analysis ----> book
+    HadrontherapyAnalysisManager::GetInstance() -> book();
+#endif
+    electromagnetic = 0;
+    hadronic = 0;
+}
 void HadrontherapyRunAction::EndOfRunAction(const G4Run*)
+{
+  //   G4cout << " Summary of Run " << aRun -> GetRunID() <<" :"<< G4endl;
+  //G4cout << "Number of electromagnetic processes of primary particles in the phantom:"
+  //     << electromagnetic << G4endl;
+  //G4cout << "Number of hadronic processes of primary particles in the phantom:"
+  //     << hadronic << G4endl;
+    // Store dose & fluence data to ASCII & ROOT files
+    if ( HadrontherapyMatrix::GetInstance() )
+    {
+        HadrontherapyMatrix::GetInstance() -> TotalEnergyDeposit();
+        HadrontherapyMatrix::GetInstance() -> StoreDoseFluenceAscii();
+#ifdef G4ANALYSIS_USE_ROOT
+        if (HadrontherapyAnalysisManager::GetInstance())
+        {
+            HadrontherapyMatrix::GetInstance() -> StoreDoseFluenceRoot();
+        }
+#endif
+    }
+    //G4cout << " Summary of Run " << aRun -> GetRunID() <<" :"<< G4endl;
+    //G4cout << "Number of electromagnetic processes of primary particles in the phantom:"
+    //     << electromagnetic << G4endl;
+    //G4cout << "Number of hadronic processes of primary particles in the phantom:"
+    //     << hadronic << G4endl;
+}
 void HadrontherapyRunAction::AddEMProcess()
 …
+}

trunk/examples/advanced/hadrontherapy/src/HadrontherapySteppingAction.cc

-                      r1230
+                      r1313
 #include "G4ParticleDefinition.hh"
 #include "G4ParticleTypes.hh"
+#include "G4UserEventAction.hh"
 #include "HadrontherapyAnalysisManager.hh"
 …
 HadrontherapySteppingAction::HadrontherapySteppingAction( HadrontherapyRunAction *run)
+{
   runAction = run;
+    runAction = run;
+}
 …
 void HadrontherapySteppingAction::UserSteppingAction(const G4Step* aStep)
+{
+    /*
+    // USEFULL METHODS TO RETRIEVE INFORMATION DURING THE STEPS
+    if( (aStep->GetTrack()->GetVolume()->GetName() == "DetectorPhys")
+    && aStep->GetTrack()->GetDefinition()->GetParticleName() == "proton")
+    //G4int evtNb = G4RunManager::GetRunManager()->GetCurrentEvent() -> GetEventID();
+    {
+    G4cout << "ENERGIA:   " << aStep->GetTrack()->GetKineticEnergy()
+    << "   VOLUME  " << aStep->GetTrack()->GetVolume()->GetName()
+    << "   MATERIALE    " <<  aStep -> GetTrack() -> GetMaterial() -> GetName()
+    << "   EVENTO     " << G4RunManager::GetRunManager()->GetCurrentEvent() -> GetEventID()
+    << "   POS     " << aStep->GetTrack()->GetPosition().x()
+    << G4endl;
+    }
+    */
     if( aStep->GetTrack()->GetVolume()->GetName() == "NewDetectorPhys"){
 #ifdef ANALYSIS_USE
       G4ParticleDefinition *def = aStep->GetTrack()->GetDefinition();
       G4double secondaryParticleKineticEnergy =  aStep->GetTrack()->GetKineticEnergy();
       G4String particleType = def->GetParticleType(); // particle type = nucleus for d, t, He3, alpha, and heavier nuclei
       G4String particleName = def->GetParticleName(); // e.g. for alpha: the name = "alpha" and type = "nucleus"
       if(particleType == "nucleus") {
         G4int A = def->GetBaryonNumber();
         G4double Z = def->GetPDGCharge();
         G4double posX = aStep->GetTrack()->GetPosition().x() / cm;
         G4double posY = aStep->GetTrack()->GetPosition().y() / cm;
         G4double posZ = aStep->GetTrack()->GetPosition().z() / cm;
         G4double energy = secondaryParticleKineticEnergy / A / MeV;
         HadrontherapyAnalysisManager* analysisMgr =  HadrontherapyAnalysisManager::getInstance();
         analysisMgr->fillFragmentTuple(A, Z, energy, posX, posY, posZ);
       } else if(particleName == "proton") {   // proton (hydrogen-1) is a special case
         G4double posX = aStep->GetTrack()->GetPosition().x() / cm ;
         G4double posY = aStep->GetTrack()->GetPosition().y() / cm ;
         G4double posZ = aStep->GetTrack()->GetPosition().z() / cm ;
         G4double energy = secondaryParticleKineticEnergy * MeV;    // Hydrogen-1: A = 1, Z = 1
         HadrontherapyAnalysisManager::getInstance()->fillFragmentTuple(1, 1.0, energy, posX, posY, posZ);
+      }
       G4String secondaryParticleName =  def -> GetParticleName();
       //G4cout <<"Particle: " << secondaryParticleName << G4endl;
       //G4cout <<"Energy: " << secondaryParticleKineticEnergy << G4endl;
         HadrontherapyAnalysisManager* analysis =  HadrontherapyAnalysisManager::getInstance();
+#ifdef G4ANALYSIS_USE_ROOT
+        G4ParticleDefinition *def = aStep->GetTrack()->GetDefinition();
+        G4double secondaryParticleKineticEnergy =  aStep->GetTrack()->GetKineticEnergy();
+        G4String particleType = def->GetParticleType(); // particle type = nucleus for d, t, He3, alpha, and heavier nuclei
+        G4String particleName = def->GetParticleName(); // e.g. for alpha: the name = "alpha" and type = "nucleus"
+        if(particleType == "nucleus") {
+            G4int A = def->GetBaryonNumber();
+            G4double Z = def->GetPDGCharge();
+            G4double posX = aStep->GetTrack()->GetPosition().x() / cm;
+            G4double posY = aStep->GetTrack()->GetPosition().y() / cm;
+            G4double posZ = aStep->GetTrack()->GetPosition().z() / cm;
+            G4double energy = secondaryParticleKineticEnergy / A / MeV;
+            HadrontherapyAnalysisManager* analysisMgr =  HadrontherapyAnalysisManager::GetInstance();
+            analysisMgr->FillFragmentTuple(A, Z, energy, posX, posY, posZ);
+        } else if(particleName == "proton") {   // proton (hydrogen-1) is a special case
+            G4double posX = aStep->GetTrack()->GetPosition().x() / cm ;
+            G4double posY = aStep->GetTrack()->GetPosition().y() / cm ;
+            G4double posZ = aStep->GetTrack()->GetPosition().z() / cm ;
+            G4double energy = secondaryParticleKineticEnergy * MeV;    // Hydrogen-1: A = 1, Z = 1
+            HadrontherapyAnalysisManager::GetInstance()->FillFragmentTuple(1, 1.0, energy, posX, posY, posZ);
+        }
+        G4String secondaryParticleName =  def -> GetParticleName();
+        //G4cout <<"Particle: " << secondaryParticleName << G4endl;
+        //G4cout <<"Energy: " << secondaryParticleKineticEnergy << G4endl;
+        HadrontherapyAnalysisManager* analysis =  HadrontherapyAnalysisManager::GetInstance();
         //There is a bunch of stuff recorded with the energy 0, something should perhaps be done about this.
         if(secondaryParticleName == "proton") {
           analysis->hydrogenEnergy(secondaryParticleKineticEnergy / MeV);
+            analysis->hydrogenEnergy(secondaryParticleKineticEnergy / MeV);
+        }
         if(secondaryParticleName == "deuteron") {
           analysis->hydrogenEnergy((secondaryParticleKineticEnergy/2) / MeV);
+            analysis->hydrogenEnergy((secondaryParticleKineticEnergy/2) / MeV);
+        }
         if(secondaryParticleName == "triton") {
           analysis->hydrogenEnergy((secondaryParticleKineticEnergy/3) / MeV);
+            analysis->hydrogenEnergy((secondaryParticleKineticEnergy/3) / MeV);
+        }
         if(secondaryParticleName == "alpha") {
           analysis->heliumEnergy((secondaryParticleKineticEnergy/4) / MeV);
+            analysis->heliumEnergy((secondaryParticleKineticEnergy/4) / MeV);
+        }
         if(secondaryParticleName == "He3"){
           analysis->heliumEnergy((secondaryParticleKineticEnergy/3) / MeV);
+            analysis->heliumEnergy((secondaryParticleKineticEnergy/3) / MeV);
+        }
 #endif
 …
+    }
+  // Electromagnetic and hadronic processes of primary particles in the phantom
+  //setting phantomPhys correctly will break something here fixme
+  if ((aStep -> GetTrack() -> GetTrackID() == 1) &&
+    (aStep -> GetTrack() -> GetVolume() -> GetName() == "PhantomPhys") &&
+    (aStep -> GetPostStepPoint() -> GetProcessDefinedStep() != NULL))
+    // Electromagnetic and hadronic processes of primary particles in the phantom
+    //setting phantomPhys correctly will break something here fixme
+    if ((aStep -> GetTrack() -> GetTrackID() == 1) &&
+            (aStep -> GetTrack() -> GetVolume() -> GetName() == "PhantomPhys") &&
+            (aStep -> GetPostStepPoint() -> GetProcessDefinedStep() != NULL))
+    {
+        G4String process = aStep -> GetPostStepPoint() ->
+            GetProcessDefinedStep() -> GetProcessName();
+        if ((process == "Transportation") || (process == "StepLimiter")) {;}
+        else
+        {
+            if ((process == "msc") || (process == "hLowEIoni") || (process == "hIoni"))
+            {
+                runAction -> AddEMProcess();
+            }
+            else
+            {
+              G4String process = aStep -> GetPostStepPoint() ->
+                GetProcessDefinedStep() -> GetProcessName();
+              if ((process == "Transportation") || (process == "StepLimiter")) {;}
+              else {
+                if ((process == "msc") || (process == "hLowEIoni") || (process == "hIoni"))
+                  {
+                    runAction -> AddEMProcess();
+                  }
+                else
+                 {
+                   runAction -> AddHadronicProcess();
+                   if ( (process != "LElastic") && (process != "ProtonInelastic") && (process != "hElastic") )
+                     G4cout << "Warning! Unknown proton process: "<< process << G4endl;
+                 }
+              }
+                runAction -> AddHadronicProcess();
+                if ( (process != "LElastic") && (process != "ProtonInelastic") && (process != "hElastic") )
+                    G4cout << "Warning! Unknown proton process: "<< process << G4endl;
+            }
+ // Retrieve information about the secondaries originated in the phantom
+ G4SteppingManager*  steppingManager = fpSteppingManager;
+  // check if it is alive
+  //if(theTrack-> GetTrackStatus() == fAlive) { return; }
+  // Retrieve the secondary particles
+  G4TrackVector* fSecondary = steppingManager -> GetfSecondary();
+  for(size_t lp1=0;lp1<(*fSecondary).size(); lp1++)
+        }
+    }
+    // Retrieve information about the secondaries originated in the phantom
+    G4SteppingManager*  steppingManager = fpSteppingManager;
+    // check if it is alive
+    //if(theTrack-> GetTrackStatus() == fAlive) { return; }
+    // Retrieve the secondary particles
+    G4TrackVector* fSecondary = steppingManager -> GetfSecondary();
+    for(size_t lp1=0;lp1<(*fSecondary).size(); lp1++)
+    {
       G4String volumeName = (*fSecondary)[lp1] -> GetVolume() -> GetName();
       if (volumeName == "phantomPhys")
+        G4String volumeName = (*fSecondary)[lp1] -> GetVolume() -> GetName();
+        if (volumeName == "phantomPhys")
+        {
 #ifdef ANALYSIS_USE
           G4String secondaryParticleName =  (*fSecondary)[lp1]->GetDefinition() -> GetParticleName();
           G4double secondaryParticleKineticEnergy =  (*fSecondary)[lp1] -> GetKineticEnergy();
           HadrontherapyAnalysisManager* analysis =  HadrontherapyAnalysisManager::getInstance();
           if (secondaryParticleName == "e-")
             analysis -> electronEnergyDistribution(secondaryParticleKineticEnergy/MeV);
           if (secondaryParticleName == "gamma")
             analysis -> gammaEnergyDistribution(secondaryParticleKineticEnergy/MeV);
           if (secondaryParticleName == "deuteron")
             analysis -> deuteronEnergyDistribution(secondaryParticleKineticEnergy/MeV);
           if (secondaryParticleName == "triton")
             analysis -> tritonEnergyDistribution(secondaryParticleKineticEnergy/MeV);
           if (secondaryParticleName == "alpha")
             analysis -> alphaEnergyDistribution(secondaryParticleKineticEnergy/MeV);
           G4double z = (*fSecondary)[lp1]-> GetDynamicParticle() -> GetDefinition() -> GetPDGCharge();
           if (z > 0.)
+            {
               G4int a = (*fSecondary)[lp1]-> GetDynamicParticle() -> GetDefinition() -> GetBaryonNumber();
               G4int electronOccupancy = (*fSecondary)[lp1] ->  GetDynamicParticle() -> GetTotalOccupancy();
               // If a generic ion is originated in the detector, its baryonic number, PDG charge,
               // total number of electrons in the orbitals are stored in a ntuple
               analysis -> genericIonInformation(a, z, electronOccupancy, secondaryParticleKineticEnergy/MeV);
+#ifdef G4ANALYSIS_USE_ROOT
+            G4String secondaryParticleName =  (*fSecondary)[lp1]->GetDefinition() -> GetParticleName();
+            G4double secondaryParticleKineticEnergy =  (*fSecondary)[lp1] -> GetKineticEnergy();
+            HadrontherapyAnalysisManager* analysis =  HadrontherapyAnalysisManager::GetInstance();
+            if (secondaryParticleName == "e-")
+                analysis -> electronEnergyDistribution(secondaryParticleKineticEnergy/MeV);
+            if (secondaryParticleName == "gamma")
+                analysis -> gammaEnergyDistribution(secondaryParticleKineticEnergy/MeV);
+            if (secondaryParticleName == "deuteron")
+                analysis -> deuteronEnergyDistribution(secondaryParticleKineticEnergy/MeV);
+            if (secondaryParticleName == "triton")
+                analysis -> tritonEnergyDistribution(secondaryParticleKineticEnergy/MeV);
+            if (secondaryParticleName == "alpha")
+                analysis -> alphaEnergyDistribution(secondaryParticleKineticEnergy/MeV);
+            G4double z = (*fSecondary)[lp1]-> GetDynamicParticle() -> GetDefinition() -> GetPDGCharge();
+            if (z > 0.)
+            {
+                G4int a = (*fSecondary)[lp1]-> GetDynamicParticle() -> GetDefinition() -> GetBaryonNumber();
+                G4int electronOccupancy = (*fSecondary)[lp1] ->  GetDynamicParticle() -> GetTotalOccupancy();
+                // If a generic ion is originated in the detector, its baryonic number, PDG charge,
+                // total number of electrons in the orbitals are stored in a ntuple
+                analysis -> genericIonInformation(a, z, electronOccupancy, secondaryParticleKineticEnergy/MeV);
+            }
 #endif

Context Navigation

Legend:

trunk/examples/advanced/hadrontherapy/GNUmakefile

trunk/examples/advanced/hadrontherapy/Hadrontherapy.cc

trunk/examples/advanced/hadrontherapy/History

trunk/examples/advanced/hadrontherapy/README

trunk/examples/advanced/hadrontherapy/defaultMacro.mac

trunk/examples/advanced/hadrontherapy/include/CVS/Entries

trunk/examples/advanced/hadrontherapy/include/CVS/Tag

trunk/examples/advanced/hadrontherapy/include/HadrontherapyAnalysisFileMessenger.hh

trunk/examples/advanced/hadrontherapy/include/HadrontherapyAnalysisManager.hh

trunk/examples/advanced/hadrontherapy/include/HadrontherapyDetectorConstruction.hh

trunk/examples/advanced/hadrontherapy/include/HadrontherapyDetectorMessenger.hh

trunk/examples/advanced/hadrontherapy/include/HadrontherapyDetectorSD.hh

trunk/examples/advanced/hadrontherapy/include/HadrontherapyEventAction.hh

trunk/examples/advanced/hadrontherapy/include/HadrontherapyInteractionParameters.hh

trunk/examples/advanced/hadrontherapy/include/HadrontherapyMatrix.hh

trunk/examples/advanced/hadrontherapy/include/HadrontherapyPhysicsList.hh

trunk/examples/advanced/hadrontherapy/include/HadrontherapyPhysicsListMessenger.hh

trunk/examples/advanced/hadrontherapy/include/HadrontherapyPrimaryGeneratorAction.hh

trunk/examples/advanced/hadrontherapy/include/PassiveProtonBeamLine.hh

trunk/examples/advanced/hadrontherapy/include/PassiveProtonBeamLineMessenger.hh

trunk/examples/advanced/hadrontherapy/src/HadrontherapyAnalysisManager.cc

trunk/examples/advanced/hadrontherapy/src/HadrontherapyDetectorConstruction.cc

trunk/examples/advanced/hadrontherapy/src/HadrontherapyDetectorMessenger.cc

trunk/examples/advanced/hadrontherapy/src/HadrontherapyEventAction.cc

trunk/examples/advanced/hadrontherapy/src/HadrontherapyMatrix.cc

trunk/examples/advanced/hadrontherapy/src/HadrontherapyPhysicsList.cc

trunk/examples/advanced/hadrontherapy/src/HadrontherapyPhysicsListMessenger.cc

trunk/examples/advanced/hadrontherapy/src/HadrontherapyPrimaryGeneratorAction.cc

trunk/examples/advanced/hadrontherapy/src/HadrontherapyRunAction.cc

trunk/examples/advanced/hadrontherapy/src/HadrontherapySteppingAction.cc

Download in other formats: