Changeset 1230 for trunk/examples/advanced/hadrontherapy

Timestamp:

Jan 8, 2010, 3:02:48 PM (14 years ago)

Author:

garnier

Message:

update to geant4.9.3

Location:

trunk/examples/advanced/hadrontherapy

Files:

• GNUmakefile (modified) (2 diffs)
• Hadrontherapy.cc (modified) (3 diffs)
• History (modified) (2 diffs)
• defaultMacro.mac (modified) (1 diff)
• include/CVS (added)
• include/CVS/Entries (added)
• include/CVS/Repository (added)
• include/CVS/Root (added)
• include/CVS/Tag (added)
• include/Decay.hh (modified) (1 diff)
• include/HadrontherapyAnalysisFileMessenger.hh (added)
• include/HadrontherapyAnalysisManager.hh (modified) (2 diffs)
• include/HadrontherapyDetectorConstruction.hh (modified) (1 diff)
• include/HadrontherapyDetectorHit.hh (added)
• include/HadrontherapyDetectorMessenger.hh (modified) (3 diffs)
• include/HadrontherapyDetectorROGeometry.hh (added)
• include/HadrontherapyDetectorSD.hh (added)
• include/HadrontherapyDummySD.hh (modified) (1 diff)
• include/HadrontherapyEventAction.hh (modified) (3 diffs)
• include/HadrontherapyEventActionMessenger.hh (added)
• include/HadrontherapyGeometryController.hh (added)
• include/HadrontherapyGeometryMessenger.hh (added)
• include/HadrontherapyInteractionParameters.hh (added)
• include/HadrontherapyMatrix.hh (modified) (4 diffs)
• include/HadrontherapyParameterMessenger.hh (added)
• include/HadrontherapyParticles.hh (modified) (1 diff)
• include/HadrontherapyPhysicsList.hh (modified) (2 diffs)
• include/HadrontherapyPhysicsListMessenger.hh (modified) (2 diffs)
• include/HadrontherapyPrimaryGeneratorAction.hh (modified) (2 diffs)
• include/HadrontherapyPrimaryGeneratorMessenger.hh (modified) (1 diff)
• include/HadrontherapyRunAction.hh (modified) (2 diffs)
• include/HadrontherapyStepMax.hh (added)
• include/HadrontherapyStepMaxMessenger.hh (added)
• include/HadrontherapySteppingAction.hh (modified) (1 diff)
• include/IAEADetectorConstruction.hh (added)
• include/IAEADetectorMessenger.hh (added)
• include/IAEAScoreWriter.hh (added)
• include/LocalINCLIonIonInelasticPhysic.hh (added)
• include/LocalIonIonInelasticPhysic.hh (added)
• include/PassiveProtonBeamLine.hh (added)
• include/PassiveProtonBeamLineMessenger.hh (added)
• src/HadrontherapyAnalysisManager.cc (modified) (7 diffs)
• src/HadrontherapyDetectorConstruction.cc (modified) (2 diffs)
• src/HadrontherapyDetectorMessenger.cc (modified) (2 diffs)
• src/HadrontherapyEventAction.cc (modified) (3 diffs)
• src/HadrontherapyMatrix.cc (modified) (5 diffs)
• src/HadrontherapyModulator.cc (modified) (6 diffs)
• src/HadrontherapyParticles.cc (modified) (1 diff)
• src/HadrontherapyPhysicsList.cc (modified) (1 diff)
• src/HadrontherapyPhysicsListMessenger.cc (modified) (1 diff)
• src/HadrontherapyPrimaryGeneratorAction.cc (modified) (7 diffs)
• src/HadrontherapyPrimaryGeneratorMessenger.cc (modified) (1 diff)
• src/HadrontherapyRunAction.cc (modified) (1 diff)
• src/HadrontherapySteppingAction.cc (modified) (5 diffs)

trunk/examples/advanced/hadrontherapy/GNUmakefile

@@ -1 @@
-# $Id: GNUmakefile,v 1.5 2004/11/30 09:06:18 guatelli Exp $
+# $Id: GNUmakefile,v 1.12 2009/08/13 20:48:04 kaitanie Exp $
 # --------------------------------------------------------------
 # GNUmakefile for examples module.  Gabriele Cosmo, 06/04/98.
@@ -15 +15 @@
 all: lib bin
 
+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
+CPPFLAGS += $(shell root-config --cflags)
+LDFLAGS  += $(shell root-config --glibs)
+endif
+endif
+
 include $(G4INSTALL)/config/binmake.gmk
 
-ifdef G4ANALYSIS_USE  
- CPPFLAGS += `aida-config --include`
- LDFLAGS  += `aida-config --lib`
- LOADLIBS += `aida-config --lib` 
+ifdef G4ANALYSIS_USE 
 endif

trunk/examples/advanced/hadrontherapy/Hadrontherapy.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: Hadrontherapy.cc Main of the Hadrontherapy example; Version 4.0 May 2005
+// Hadrontherapy.cc 
+//
+// Main of the Hadrontherapy example; 
+// Released with the Geant4 9.3 version (December 2009)
+//
+// Last modified: G.A.P.Cirrone
+// 
+// See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
+//
 // ----------------------------------------------------------------------------
 //                 GEANT 4 - Hadrontherapy example
@@ -30 +38 @@
 // Code developed by:
 //
-// G.A.P. Cirrone(a)*, F. Di Rosa(a), S. Guatelli(b), G. Russo(a)
+// 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)
+//
+// 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, Genova, Italy
+//
+// (b) INFN Section of Genova, Italy
 // 
-// * cirrone@lns.infn.it
+// (c) University of Wallongong, Australia
+//
+// (d) Helsinki Institute of Physics, Helsinki, Finland
+//
+// (e) CERN, (CH)
+//
+//  *Corresponding author, email to cirrone@lns.infn.it
 // ----------------------------------------------------------------------------
+
 #include "G4RunManager.hh"
 #include "G4UImanager.hh"
 #include "G4UIterminal.hh"
 #include "G4UItcsh.hh"
-#ifdef G4UI_USE_XM
-#include "G4UIXm.hh"
-#endif
-#ifdef G4VIS_USE
-#include "G4VisExecutive.hh"
-#endif
 #include "HadrontherapyEventAction.hh"
-#include "HadrontherapyDetectorConstruction.hh"
 #include "HadrontherapyPhysicsList.hh"
-#include "HadrontherapyPhantomSD.hh"
+#include "HadrontherapyDetectorSD.hh"
 #include "HadrontherapyPrimaryGeneratorAction.hh"
 #include "HadrontherapyRunAction.hh"
@@ -61 +76 @@
 #include "globals.hh"
 #include "HadrontherapySteppingAction.hh"
-#ifdef  G4ANALYSIS_USE
 #include "HadrontherapyAnalysisManager.hh"
-#endif
-
+#include "HadrontherapyGeometryController.hh"
+#include "HadrontherapyGeometryMessenger.hh"
+#include "HadrontherapyInteractionParameters.hh"
+#include "G4ScoringManager.hh"
+#include "IAEAScoreWriter.hh"
+
+#if defined(G4UI_USE_TCSH)
+#include "G4UIterminal.hh"
+#include "G4UItcsh.hh"
+#endif
+
+#ifdef G4UI_USE_XM
+#include "G4UIXm.hh"
+#endif
+
+#ifdef G4VIS_USE
+#include "G4VisExecutive.hh"
+#endif
+
+#ifdef G4UI_USE_QT
+#include "G4UIQt.hh"
+#include "G4Qt.hh"
+#endif
+
+//////////////////////////////////////////////////////////////////////////////////////////////
 int main(int argc ,char ** argv)
 {
-
   // Set the Random engine
-
   CLHEP::HepRandom::setTheEngine(new CLHEP::RanecuEngine());
 
   G4RunManager* runManager = new G4RunManager;
 
-  // Initialize the geometry
-  runManager -> SetUserInitialization(new HadrontherapyDetectorConstruction());
-  
+  //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);
   
-  // Initialize the primary particles  
-  runManager -> SetUserAction(new HadrontherapyPrimaryGeneratorAction());
-
-  // Initialize matrix 
-  HadrontherapyMatrix* matrix = new HadrontherapyMatrix();
-  matrix -> Initialize();
-
   // Optional UserActions: run, event, stepping
   HadrontherapyRunAction* pRunAction = new HadrontherapyRunAction();
   runManager -> SetUserAction(pRunAction);
 
-  HadrontherapyEventAction* pEventAction = new HadrontherapyEventAction(matrix);
+  HadrontherapyEventAction* pEventAction = new HadrontherapyEventAction();
   runManager -> SetUserAction(pEventAction);
-
 
   HadrontherapySteppingAction* steppingAction = new HadrontherapySteppingAction(pRunAction); 
   runManager -> SetUserAction(steppingAction);    
 
-
-#ifdef G4ANALYSIS_USE
-  HadrontherapyAnalysisManager* analysis = 
-    HadrontherapyAnalysisManager::getInstance();
-  analysis -> book();
-#endif
-  
+  // Interaction data: stopping powers
+  HadrontherapyInteractionParameters* pInteraction = new HadrontherapyInteractionParameters();
+
 #ifdef G4VIS_USE
   // Visualization manager
   G4VisManager* visManager = new G4VisExecutive;
   visManager -> Initialize();
-#endif
+#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
+#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
+     // 
+#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
+#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;
+#endif
   
-  G4UIsession* session = 0;
-  if (argc == 1)   // Define UI session for interactive mode.
-    {
-      session = new G4UIterminal();
-    } 
-
-  // Get the pointer to the User Interface manager 
-  G4UImanager* UI = G4UImanager::GetUIpointer();  
-  if (session)   // Define UI session for interactive mode.
-    { 
-      G4cout<<" UI session starts ..."<< G4endl;
-      UI -> ApplyCommand("/control/execute defaultMacro.mac");    
-      session -> SessionStart();
-      delete session;
-    }
-  else           // Batch mode
-    { 
-      G4String command = "/control/execute ";
-      G4String fileName = argv[1];
-      UI -> ApplyCommand(command + fileName);
-    }  
-
-  matrix -> TotalEnergyDeposit();
-
-#ifdef G4ANALYSIS_USE
-  analysis -> finish();
-#endif
-  
-  // Job termination
-#ifdef G4VIS_USE
-  delete visManager;
-#endif
-
+  delete geometryMessenger;
+  delete geometryController;
+  delete pInteraction; 
   delete runManager;
-
   return 0;
 }

trunk/examples/advanced/hadrontherapy/History

@@ -1 @@
------------------------------------------------------------
-$Id: History, v 1.6 2004/02/27  G.A.P. Cirrone
------------------------------------------------------------
-
-     ====================================================
-                     Geant4 - Hadrontherapy
-     ====================================================
-
-                      Category History file
+ -------------------------------------------------------------------------------
+History File, 2004/02/27 G.A.P. Cirrone, Created
+cirrone@lns.infn.it
+http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
+-------------------------------------------------------------------------------
+
+         ====================================================
+             History file of the Hadrontherapy application
+         ====================================================
+
+25.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. 
+
+22.11.2009 S.E.Mazzaglia; Tag: Hadrontherapy-V09-02-39
+           - Correction in the initialization of the passiveProtonBeamLine class.
+
+20.11.2009 P.Kaitaniemi; Tag: hadrontherapy-V09-02-38
+           - Fixes and updates to the analysis scripts in RootScripts/iaeaBenchmark
+           - Updated RootScripts/README
+
+19.11.2009 G.A.P.Cirrone; Tag: hadrontherapy-V09-02-37
+           - Minor revisions;
+
+18.11.2009 G.A.P.Cirrone; Tag: hadrontherapy-V09-02-36
+           - Correction for a missing function in the HadrontherapyPhysicsListMessenger.cc class file
+
+18.11.2009 G.A.P.Cirrone; Tag: hadrontherapy-V09-02-35
+           - Updated the README file and general code revision for the 
+             Geant4 9.3 release
+
+17.11.2009 S.E.Mazzaglia; Tag: hadrontherapy-V09-02-34
+           - Added some functionalities in order to change, via messengers, the geometry, the voxelization 
+             of the detector, and the disposition in the space of the detector/phantom.
+           - Added the possibility to calculate the stopping powers for ions too.
+           - Modified the HadrontherapyDetectorROGeometry class constructor.
+           - Various additions and fixes to the matrix class.
+
+10.11.2009 G.A.P.Cirrone; Tag: hadrontherapy-V09-02-33
+           - Added the possibility to make a graphical user interface (GUI) using the QT libraries. 
+             To start a GUI the correct enviroment variables must be configured (see the Geant4 installation 
+             manual) and a QT version must be installed  
+
+05.10.2009 P.Kaitaniemi; Tag: hadrontherapy-V09-02-32
+           - Fixed a compilation error with GCC 4.4
+
+28.09.2009 S.E.Mazzaglia; Tag: hadrontherapy-V09-02-31
+           - Now the HadrontherapyDetectorConstruction class implements only phantom and detector [RO]geometry.
+             World volume and the rest of the geometry is inside another class whose messenger allows 
+             modification by users with the same old syntax
+           - Removed HadrontherapyInteractionParameters from the HadrontherapyGeometryController class 
+                          
+20.09.2009 P.Kaitaniemi; Tag: hadrontherapy-V09-02-30
+           - Added ability to use command based scoring
+           - IAEA geometry: produce Bragg peak using command based scoring
+           - Various additions and fixes to the IAEA ROOT scripts
+           - Additional data extracted from E. Haettner's thesis
+
+20.09.2009 P.Kaitaniemi; Tag: hadrontherapy-V09-02-29
+           - Moved HadrontherapyInteractionParameters initialization to the HadrontherapyGeometryController class
+
+11.09.2009 G.A.P.Cirrone; Tag: hadrontherapy-V09-02-28
+           - Added messengers to control the event number and to draw only particular tracks.
+             The new command are accessible via the command /event/drawTracks and /event/PrintEventNumber  
+
+08.09.2009 S.E.Mazzaglia; Tag: hadrontherapy-V09-02-27
+           - Added a method to retrieve stopping power values for protons, alphas and electrons.
+             This method is implemented in the new class HadrontherapyInteractionParameters
+
+13.08.2009 P.Kaitaniemi; Tag: hadrontherapy-V09-02-26
+           - Fixed compilation errors when AIDA analysis is used
+
+03.08.2009 P.Kaitaniemi; Tag: hadrontherapy-V09-02-25
+           - Added ability to select the geometry using G4 macro commands
+           - Improved plotting scripts and improved normalization for the
+             fragment energy distribution
+
+27.07.2009 P.Kaitaniemi; Tag: hadrontherapy-V09-02-24
+           - IAEA geometry: added ability to remove the phantom by setting its thickness to zero
+           - Collect simulation metadata: number of events, distance of the detector (IAEA geometry),
+             depth of the phantom (IAEA geometry), beam energy, energy error
+           - Added ability to produce angular distribution plots
+
+17.07.2009 P.Kaitaniemi; Tag: hadrontherapy-V09-02-23
+           - Tuned geometry of the E. Haettner experimet (IAEADetectorConstruction)
+           - Adopted G4ANALYSIS_USE_ROOT flag to activate ROOT analysis
+           - Improved plotting scripts
+
+13.07.2009 P.Kaitaniemi; Tag: hadrontherapy-V09-02-22
+           - Added the first version of the IAEA benchmark geometry based on
+              E. Haettner's thesis
+           - Collect fragment energy distributions
+           - Added fragment energy distribution data
+           - ROOT script preparing an IAEA benchmark figure with data
+
+08.07.2009 G.A.P.Cirrone; Tag: hadrontherapy-V09-02-21
+           - Removed the README file in ASCII format 
+
+27.06.2009 G.A.P.Cirrone; Tag: hadrontherapy-V09-02-20
+           - Eliminated not necessary dependences in the SteppingAction class
+           - Added folders containing experimental data (its name is 'experimentalData') and 
+             ROOT scripts ('RootScripts') where Root scripts are stored to
+             perform a fast comparison with experimental data.
+             A folder where simulation results are stored is also created. Its
+             name is 'simulationResults'. 
+
+27.06.2009 P. Kaitaniemi; Tag hadrontherapy-V09-02-19
+           - Added ability to change the name of the output file between runs
+
+26.06.2009 P. Kaitaniemi; Tag hadrontherapy-V09-02-18
+           - Fixed a bug in the physics list. Local ion-ion hadronic physics was not loaded due
+             to an uninitialized variable (locIonIonInelasticIsRegistered)
+           - Ability to use /analysis/setAnalysisFile <filename> to set the name of the output file
+           - Added Doxygen documentation tags to the source code and Doxyfile for
+             documentation settings
+           - Support for direct use of ROOT for analysis in addition to the default AIDA one
+           - Local INCL/ABLA physics list for deuterons, tritons and alphas
+
+26.06.2009 G.A.P.Cirrone; Tag hadrontherapy-V09-02-17
+           - Corrected the definition of total inelastic cross section for light ions in the 
+             LocalIonIonInelasticPhysic.cc file
+
+26.06.2009 G.A.P.Cirrone; Tag: hadrontherapy-V09-02-16
+           - Momentarely removed the class for LET calculation 
+             for a conflict with the general structure of Hadrontherapy
+
+10.06.2009 G.A.P.Cirrone; Tag: hadrontherapy-V09-02-15
+           - Corrected a bug in the detector construction
+
+05.06.2009 G.A.P.Cirrone and S.Mazzaglia; Tag: hadrontherapy-V09-02-14
+           - Added a preliminary version of classes for LET calculation.
+
+30.05.2009 G.A.P.Cirrone; Tag hadrontherapy-V09-02-13
+           - README_Hadrontherapy.pdf file updated and improved
+
+29.05.2009 G.A.P.Cirrone; Tag hadrontherapy-V09-02-12
+           - Implemented the new Low energy models (Livermore and Penelope)
+             now migrated to the new interface (common to the Standard models)
+             Livermore and Penelope models can be implemented:
+             Activating the buit-in physics lists (G4EmLivermorePhysics and G4EmPenelopePhysics)
+             Activation can be done via macro commands in the usual way
+
+19.05.2009 F.Romano; Tag hadrontherapy-V09-02-11
+           - Corrected the stepMax value in each macro in order to avoid 
+             a wrong dose deposition in the first slice.
+           - Modified and revised the README and macro files.
+
+15.05.2009 G.A.P.Cirrone; Tag hadrontherapy-V09-02-10
+           - Corrected a but in the call of a physic list
+           - Corrected a bug in the proton_therapy.mac file
+
+14.05.2009 G.A.P.Cirrone; Tag hadrontherapy-V09-02-09
+           - Definitively added the StepMax class to change the max step lenght
+
+14.05.2009  G.A.P.Cirrone; Tag: hadrontherapy-V09-02-08
+           - README file improved.
+
+14.05.2009  G.A.P.Cirrone; Tag: hadrontherapy-V09-02-07
+           - Physic implementation completely changed. Now Hadrontherapy can be launched 
+             with physics lists, packages and built-in physic models;
+             In the README we give some suggestion in the physic models to use.
+             All models can be activated via macro command.
+           - Improved HadrontherapyModulator.cc file;
+
+29.03.2009  G.A.P.Cirrone; Tag: hadrontehrapy-V09-02-06
+           - Extended limits of Binary Cascade in HIProtonNeutronBinary.cc
+           - Corrected and improved "default" macro file.
+           - Improved HadrontherapyDetectorConstruction.cc file;
+           - Comments on HIProtonneutronPrecompound.cc file;
+           - Improved physicsHadronicPrecompound.mac file;
+
+18.03.2009  G.A.P.Cirrone; Tag: hadrontherapy-V09-02-05
+           - Corrected macro file for the use of the QGSP_BIC package
+
+18.03.2009: G.A.P.Cirrone; Tag: hadrontherapy-V09-02-04
+           - Added commands to Detector Messenger to give the possibility to choose beetween
+             different beam lines
+           - Added comments to HadrontherapyMatrix;
+           _ Improved code and added comments to HIProtonNeutronPrecompound and
+             HIProtonNeutronBinary;
+           - Corrected macro file using the Precompound inelastic model;
+           - Removed the class file HadrontherapyMaterial and improved all the geometry files
+
+05.03.2009: G.A.P.Cirrone; Tag: hadrontherapy-V09-02-03
+           - Updated README
+
+02.03.2009: G.A.P.Cirrone; Tag: hadrontherapy-V09-02-02
+           - Changed name of HadrontherapyBeamLine file to PassiveProtonBeamLine
+
+02.03.2009: G.A.P.Cirrone; Tag:  hadrontherapy-V09-02-01
+           - Added generation of ASCII file with dose deposited in the phantom voxels
+
+22.02.2009: G.Folger; Tag: hadrontherapy-V09-02-00
+           - Fix a compilation warning on used ionShenCrossSection in
+             HIIonLEP.cc
+
+20.11.2008: G.A.P.Cirrone and M.Russo; Tag: hadrontherapy-V09-01-11
+           - Fixed path of macro files
+
+20.11.2008: G.A.P.Cirrone; Tag: hadrontherapy-V09-01-10
+           - Updated the History file
+           - Corrected cross sections definitions for ions
+           - Revised the definition and use of the electromagnetic options
+             for the use with the Standard models
+             
+20.11.2008: G.A.P.Cirrone and M.Russo; Tag: hadrontherapy-V09-01-09
+           - Updated readme and improved the comments. 
+
+20.11.2008: G.A.P.Cirrone and M.Russo; Tag: hadrontherapy-V09-01-08
+           - Add new approach for the choice of the physic models.
+             Now packaged physic lists can be used alternatively 
+             to the the physic models implemented in the class files
+             EM, HE and HI. 
+           - Improved the electromagnetic models for the generic ions
+
+22.09.2008 G.A.P.Cirrone; Tag: hadrontherapy-V09-01-07
+           - Corrected the G4eBremsstrahlung() process in the file 
+             EMElectronStandard.cc;
+           - Updated the head of the History file;
+
+17.09.2008 A.Lechner; Tag: hadrontherapy-V09-01-06
+           - Corrections in the Low Energy Electromagnetic physic lists.
+
+15.06.2008 G.A.P.Cirrone; Tag: hadrontherapy-V09-01-05
+           - Removed AIDA call from GNUmakefile 
+
+19.05.2008 G.A.P.Cirrone tag hadrontherapy-V09-01-04
+           - Added in the beam line the MOPI detector. MOPI is a microstrip
+             detector that, in the real case, is able to check during
+             the treatment, the beam simmetry of the therapy beam.
+             Its physical structure is here exactly simulated so that 
+             the its contribute to the energy loss can be take into account;
+             A detailed description if the detector can be found in 
+             NIM A 572 (2007) 1094-1101 and its references.
+           - Corrected the position of the Phantom and Detector;
+           - Added variables to the HadrontherapyBeamLine.cc file;
+           - Added comments to the HadrontherapyBeamLine.cc file
+             to improve the clearness.
+           - Updated the README file.
+           - Changed the default dimensions of histogram bins 
+             (from 200 um to 100 um).
+
+09.03.2008 G.A.P.Cirrone tag hadrontherapy-V09-01-03
+           - Completed the update of the new beam line
+
+09.03.2008 G.A.P.Cirrone tag hadrontherapy-V09-01-02
+           - Added comments to the PhysicsList class;
+           - Eliminated not used production cuts in PhysicsList;
+           - Added NIST definition materials in Material class;
+           - Code review of the DetectorConstruction class;
+           - Changed name of the volume where the energy deposited is collected
+             from "phantom" to "detector". "Detector" is a more appropiate 
+             name.
+           - Changed name of the volume where the detector is inserted from 
+             "patient" to the more appropriate "Water Phantom";
+
+03.03.2008 G.A.P.Cirrone tag hadrontherapy-V09-01-01
+           - Added the generation of .root file;
+           - Removed a segmentation due to an uncorect pointer
+             in the EMHadronIonStandard class;
+           - Added options for an accurate use of Standard electromagnetic models
+             in the EMHadronIonStandard, EMElectronStandard, 
+             EMPositronStandard, EMPhotonStandard  and EMMuonStandard classes;
+           - Added a macro file (physicsElectromagneticStandard.mac) 
+             for the use of Hadrontherapy with the Standard Electromagnetic models;
+           - Corrected in the defaultMacro.mac, a wrong command for the 
+             activation of the Standard Electromagnetic models;
+         
+29.02.2007 G.A.P.Cirrone tag hadrontherapy-V09-01-00
+           - Updated README
 
 16.11.2007 Anton Lechner tag hadrontherapy-V09-00-00
@@ -31 +292 @@
 
 07.05.2007 G.A.P. Cirrone (hadrontherapy-V08-02-02)
-         - Geometry upgrade(hadrontherapyBeamLine class) according to the experimental CATANA 
-            proton therapy beam line; 
+         - Geometry upgrade(hadrontherapyBeamLine class) according 
+           to the experimental CATANA proton therapy beam line; 
 
 23.04.2007 S. Guatelli (hadrontherapy-V08-02-01)

trunk/examples/advanced/hadrontherapy/defaultMacro.mac

@@ -1 @@
-#----------------------------------------------------------------------------
-# DEFAULT MACRO FOR THE  
-# HADRONTHERAPY EXAMPLE
+# G.A.P.Cirrone
 #
+# Default macro file. It is called if no argument is provided at run
+# 
+# i.e. simply typing $G4WORKDIR/bin/Linux-++/Hadrontherapy <no argument here!> 
 #
-# THIS MACRO SIMPLY PERMIT TO  RUN A SIMULATION 
-# WITHOUT THE VISUALISATION  
+# This macro can be used for a proton beam in water. Both electrmagnetic and
+# hadronic models are swiched on
+
+#########################
+# Set of the verboses 
 #
-# THE RANGE SHIFTER MATERIAL AND THICKNESS CAN BE SPECIFIED 
-#
-# NOTE THAT THE MODULATOR MATERIAL IS POLTMETHYLMETHACRILATE 
-# (PMMA) FOR DEFAULT. IF ONE WANT CARRY OUT A SIMULATION WITHOUT 
-# THE MODULATOR HE/SHE MUST SET "Air" the <<ModMater>> in the
-# <<GetMater>> function of the HadrontherapyModulator.cc class
-#
-# USERS SHOULD GIVE A LOOK TO THE HELP OF THE IDLE TO KNOW
-# THE ACTIVATED MESSSENGERS FOR THE GEOMETRY
-#
-# ADDITIONAL INFORMATIONS ON THE MESSENGER AVAILABLE CAN BE FOUND
-# INSIDE THE HADRONTHERAPY DOCUMENTATION (http://www.ge.infn.it/geant4/examples/).
-# 
-# ANYWAY SEND ME AN E-MAIL FOR ANY QUESTION: cirrone@lns.infn.it.   
-# --------------------------------------------------------------------------------  
-
-
 /control/verbose 1
 /tracking/verbose 0
-/run/verbose 0
+/run/verbose 1
 /event/verbose 0
 
-# SETTING FOR THE PHYSICS MODELS
-/physics/addPhysics Decay
-/physics/addPhysics EM-Photon-EPDL
-/physics/addPhysics EM-Electron-EEDL
-/physics/addPhysics EM-Positron-Standard
-/physics/addPhysics EM-HadronIon-LowE
-/physics/addPhysics EM-Muon-Standard
-/physics/addPhysics HadronicEl-HadronIon-LElastic
-/physics/addPhysics HadronicInel-ProtonNeutron-LEP
-/physics/addPhysics HadronicInel-Ion-LEP
-/physics/addPhysics HadronicInel-Pion-LEP
-/physics/addPhysics HadronicAtRest-MuonMinus-Capture
+##########################
+# Set of the physic models
+#
+/physic/addPhysics emstandard_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)
 
-# FIX THE FOLLOWIG PARAMETERS
-# TO SET THE RANGE SHIFTER
-#/beamLine/RangeShifter/thickness 4 cm
-#/beamLine/RangeShifter/RSMat Water
-#/tracking/verbose 1
-# SET OF THE VISUALISATION CHARACTERISTICS
+##########################
+# Initialisation procedure
+#
+/run/initialize
+
+##########################
+# Visualisation 
+#
 /vis/scene/create
+#/vis/open OGLIQt # only if QT library are installed
 /vis/open OGLIX
 /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
+/vis/scene/endOfEventAction accumulate -1
 /vis/viewer/update
-/run/beamOn 100
 
+##########################
+# Set here the cut and the step max for the tracking. 
+# Suggested values of cut and step:
+#
+/physic/setCuts 0.01 mm
+/Step/waterPhantomStepMax 0.01 mm 
 
+#########################
+# Set the primary particle type,
+# energy and position along the X direction
+#
+/gun/particle proton
+/beam/energy/meanEnergy 62 MeV
+/beam/energy/sigmaEnergy 400 keV
+/beam/position/Xposition -2700 mm
 
+#########################
+# Display the event number
+# during the run
+# 
+/event/printEventNumber 10
 
-
-
-
-
-
-
-
+#########################
+# Start of the run
+#
+/run/beamOn 10

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: HadrontherapyDetectorMessenger.hh; May 2005
-// ----------------------------------------------------------------------------
-//                 GEANT 4 - Hadrontherapy example
-// ----------------------------------------------------------------------------
-// Code developed by:
-//
-// G.A.P. Cirrone(a)*, F. Di Rosa(a), S. Guatelli(b), G. Russo(a)
-// 
-// (a) Laboratori Nazionali del Sud 
-//     of the INFN, Catania, Italy
-// (b) INFN Section of Genova, Genova, Italy
-// 
-// * cirrone@lns.infn.it
-// ----------------------------------------------------------------------------
+// $Id: HadrontherapyDetectorMessenger.hh;
+// See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
 
 // ==============================

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// ----------------------------------------------------------------------------
-// $Id: HadrontherapyAnalysisManager.hh; May 2005
-// ----------------------------------------------------------------------------
-//                 GEANT 4 - Hadrontherapy example
-// ----------------------------------------------------------------------------
-// Code developed by:
-//
-// G.A.P. Cirrone(a)*, F. Di Rosa(a), S. Guatelli(b), G. Russo(a)
-// 
-// (a) Laboratori Nazionali del Sud 
-//     of the INFN, Catania, Italy
-// (b) INFN Section of Genova, Genova, Italy
-// 
-// * cirrone@lns.infn.it
-// ----------------------------------------------------------------------------
-
-#ifdef G4ANALYSIS_USE
+// HadrontherapyAnalysisManager.hh; May 2005
+// See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
+
 #ifndef HADRONTHERAPYANALYSISMANAGER_HH
 #define HADRONTHERAPYANALYSISMANAGER_HH 1
 
 #include "globals.hh"
-# include <AIDA/AIDA.h>
+
+#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 ITree;
   class IAnalysisFactory;
   class ITreeFactory;
 }
-
+#endif
+
+#ifdef G4ANALYSIS_USE_ROOT ///< If analysis is done directly with ROOT
+#include "TROOT.h"
+#include "TFile.h"
+#include "TNtuple.h"
+#include "TH1F.h"
+#endif
+
+/**
+ * Messenger class for analysis-settings for HadronTherapyAnalysisManager 
+ */
+class HadrontherapyAnalysisFileMessenger;
+
+/**
+ * A class for connecting the simulation to an analysis package.
+ */
 class HadrontherapyAnalysisManager
 {
 private:
+        /**
+         * 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();
+         *
+         * @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);
   
-  static HadrontherapyAnalysisManager* getInstance();
-  
-  void book();
-  // Book the histograms and ntuples in a .hbk file
-  
-  void FillEnergyDeposit(G4int voxelXId, G4int voxelYId, G4int voxelZId, 
+  /**
+  * Fill the ntuple with the energy deposit in the phantom
+  */
+  void FillEnergyDeposit(G4int voxelXId, G4int voxelYId, G4int voxelZId,
                          G4double energyDeposit);
-  // Fill the ntuple with the energy deposit in the phantom 
-
-  void BraggPeak(G4int, G4double);
-  // Fill 1D histogram with the Bragg peak in the phantom
+
+  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
+  ///< 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
+  ///< 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
+  ///< 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
+  ///< 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
+  ///< 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
+  ///< 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
+  ///< 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
+  ///< Fill 1D histogram with the energy deposit of secondary pions
 
   void electronEnergyDistribution(G4double secondaryParticleKineticEnergy);
-  // Energy distribution of secondary electrons originated in the phantom
+  ///< Energy distribution of secondary electrons originated in the phantom
 
   void gammaEnergyDistribution(G4double secondaryParticleKineticEnergy);
-  // Energy distribution of secondary gamma originated in the phantom
+  ///< Energy distribution of secondary gamma originated in the phantom
 
   void deuteronEnergyDistribution(G4double secondaryParticleKineticEnergy);
-  // Energy distribution of secondary deuterons originated in the phantom
+  ///< Energy distribution of secondary deuterons originated in the phantom
 
   void tritonEnergyDistribution(G4double secondaryParticleKineticEnergy);
-  // Energy distribution of secondary tritons originated in the phantom
+  ///< Energy distribution of secondary tritons originated in the phantom
 
   void alphaEnergyDistribution(G4double secondaryParticleKineticEnergy);
-  // Energy distribution of secondary alpha originated in the phantom
+  ///< 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
+  ///< Close the .hbk file with the histograms and the ntuples
+
+void flush();
+
+#ifdef G4ANALYSIS_USE_ROOT
+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
 
 private:
   static HadrontherapyAnalysisManager* instance;
+  HadrontherapyAnalysisFileMessenger* fMess;
+  G4String analysisFileName;
+#ifdef G4ANALYSIS_USE
   AIDA::IAnalysisFactory* aFact;
-  AIDA::ITree* theTree; 
+  AIDA::ITree* theTree;
   AIDA::IHistogramFactory *histFact;
   AIDA::ITupleFactory *tupFact;
@@ -129 +185 @@
   AIDA::IHistogram1D *h5;
   AIDA::IHistogram1D *h6;
-  AIDA::IHistogram1D *h7; 
-  AIDA::IHistogram1D *h8; 
+  AIDA::IHistogram1D *h7;
+  AIDA::IHistogram1D *h8;
   AIDA::IHistogram1D *h9;
   AIDA::IHistogram1D *h10;
   AIDA::IHistogram1D *h11;
-  AIDA::IHistogram1D *h12; 
-  AIDA::IHistogram1D *h13; 
+  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
-#endif
-
+
+#endif
+

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: HadrontherapyDetectorConstruction.hh; Version 4.0 May 2005
-// ----------------------------------------------------------------------------
-//                 GEANT 4 - Hadrontherapy example
-// ----------------------------------------------------------------------------
-// Code developed by:
-//
-// G.A.P. Cirrone(a)*, F. Di Rosa(a), S. Guatelli(b), G. Russo(a)
-// 
-// (a) Laboratori Nazionali del Sud 
-//     of the INFN, Catania, Italy
-// (b) INFN Section of Genova, Genova, Italy
-// 
-// * cirrone@lns.infn.it
-// ----------------------------------------------------------------------------
+// HadrontherapyDetectorConstruction.hh; 
+// See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy//
 
 #ifndef HadrontherapyDetectorConstruction_H
 #define HadrontherapyDetectorConstruction_H 1
 
+#include "G4Box.hh"
 #include "globals.hh"
-#include "G4VUserDetectorConstruction.hh"
+#include "G4VisAttributes.hh"
+#include "G4LogicalVolume.hh"
+#include "G4UnitsTable.hh"
 
 class G4VPhysicalVolume;
 class G4LogicalVolume;
-class HadrontherapyPhantomROGeometry;
-class HadrontherapyBeamLine;
+class HadrontherapyDetectorROGeometry;
 class HadrontherapyDetectorMessenger;
-class HadrontherapyModulator;
-class HadrontherapyPhantomSD;
-class HadrontherapyMaterial;
+class HadrontherapyDetectorSD;
+class HadrontherapyMatrix;
 
-class HadrontherapyDetectorConstruction : public G4VUserDetectorConstruction
+class HadrontherapyDetectorConstruction 
 {
 public:
 
-  HadrontherapyDetectorConstruction();
+  HadrontherapyDetectorConstruction(G4VPhysicalVolume*);
 
   ~HadrontherapyDetectorConstruction();
 
-  G4VPhysicalVolume* Construct();  
 
 private: 
 
-  void ConstructBeamLine();
-  // This method allows to define the beam line geometry in the
-  // experimental set-up
+  void ConstructPhantom();
+  void ConstructDetector();
+  void ConstructSensitiveDetector(G4ThreeVector position_respect_to_WORLD);
+  
+public: 
+// Get detector position relative to WORLD
+inline G4ThreeVector GetDetectorToWorldPosition()
+  {
+    return phantomPosition + detectorPosition;
+  }
+/////////////////////////////////////////////////////////////////////////////
+// Get displacement between phantom and detector by detector position, phantom and detector sizes
+inline G4ThreeVector GetDetectorToPhantomPosition()
+{
+    return G4ThreeVector(phantomSizeX - detectorSizeX + detectorPosition.getX(),
+                         phantomSizeY - detectorSizeY + detectorPosition.getY(),
+                         phantomSizeZ - detectorSizeZ + detectorPosition.getZ()
+                          );
+}
 
- void ConstructPhantom(); 
- // This method allows to define the phantom geometry in the
- // experimental set-up
- 
- void ConstructSensitiveDetector();
-  // The sensitive detector is associated to the phantom volume
+/////////////////////////////////////////////////////////////////////////////
+// Calculate (and set) detector position by displacement, phantom and detector sizes
+inline void SetDetectorPosition()
+  {
+          // Adjust detector position
+          detectorPosition.setX(detectorToPhantomPosition.getX() - phantomSizeX + detectorSizeX);
+          detectorPosition.setY(detectorToPhantomPosition.getY() - phantomSizeY + detectorSizeY);
+          detectorPosition.setZ(detectorToPhantomPosition.getZ() - phantomSizeZ + detectorSizeZ);
+     
+      if (detectorPhysicalVolume) detectorPhysicalVolume -> SetTranslation(detectorPosition); 
+  }
+/////////////////////////////////////////////////////////////////////////////
+// Check whether detector is inside phantom
+inline bool IsInside(G4double detectorHalfX,
+                     G4double detectorHalfY,
+                     G4double detectorHalfZ,
+                     G4double phantomHalfX,
+                     G4double phantomHalfY,
+                     G4double phantomHalfZ,
+                     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) 
+                 {
+                    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()) 
+                 {
+                    G4cout << "Error: X dimension doesn't fit with detector to phantom relative position" << G4endl;
+                    return false;
+                 }
+        }
 
-public: 
+        if (detectorHalfY > 0. && phantomHalfY > 0.&& detectorToPhantomPosition.getY() >=0.)
+        {
+            if (detectorHalfY > phantomHalfY) 
+                 {
+                    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()) 
+             {
+                   G4cout << "Error: Y dimension doesn't fit with detector to phantom relative position" << G4endl;
+                   return false;
+             }
+        }                        
 
-  void SetModulatorAngle(G4double angle);
-  // This method allows moving the modulator through UI commands
+        if (detectorHalfZ > 0. && phantomHalfZ > 0.&& detectorToPhantomPosition.getZ() >=0.)
+        {
+            if (detectorHalfZ > phantomHalfZ) 
+                 {
+                    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()) 
+             {
+                   G4cout << "Error: Z dimension doesn't fit with detector to phantom relative position" << G4endl;
+                   return false;
+             }
+        }
+/*
+    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;
+}
+/////////////////////////////////////////////////////////////////////////////
 
-  void SetRangeShifterXPosition(G4double translation);
-  // This method allows to move the Range Shifter along
-  // the X axis through UI commands
+  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);
+  G4LogicalVolume* GetDetectorLogicalVolume(){ return detectorLogicalVolume;}
 
-  void SetRangeShifterXSize(G4double halfSize);
-  // This method allows to change the size of the range shifter along
-  // the X axis through UI command.
 
-  void SetFirstScatteringFoilSize(G4double halfSize);
-  // This method allows to change the size of the first scattering foil
-  // along the X axis through UI command.
+private:
 
-  void SetSecondScatteringFoilSize (G4double halfSize); 
-  // This method allows to change the size of the second scattering foil
-  // along the X axis through UI command.
+  HadrontherapyDetectorMessenger* detectorMessenger; 
 
-  void SetOuterRadiusStopper (G4double value); 
-  // This method allows to change the size of the outer radius of the stopper
-  // through UI command.
+  G4VisAttributes* skyBlue;
+  G4VisAttributes* red;
 
-  void SetInnerRadiusFinalCollimator (G4double value);
-  // This method allows to change the size of the inner radius of the 
-  // final collimator through UI command.
+  G4VPhysicalVolume* motherPhys;
 
-  void SetRSMaterial(G4String material);
-  // This method allows to change the material 
-  // of the range shifter through UI command.
+  HadrontherapyDetectorSD*         detectorSD; // Pointer to sensitive detector
+  HadrontherapyDetectorROGeometry* detectorROGeometry; // Pointer to ROGeometry 
+  HadrontherapyMatrix*             matrix;
 
-  G4double ComputeVoxelSize() {return phantomSizeX/numberOfVoxelsAlongX;};
-  // Returns the size of the voxel along the X axis
- 
-private:
+  G4VPhysicalVolume* phantomPhysicalVolume;
+  G4LogicalVolume*   phantomLogicalVolume; 
+  G4LogicalVolume*   detectorLogicalVolume;
+  G4VPhysicalVolume* detectorPhysicalVolume;
   
-  HadrontherapyPhantomSD* phantomSD; // Pointer to sensitive detector
-
-  HadrontherapyPhantomROGeometry* phantomROGeometry; // Pointer to ROGeometry 
-
-  HadrontherapyBeamLine* beamLine; // Pointer to the beam line 
-                                   // geometry component
-
-  HadrontherapyModulator* modulator; // Pointer to the modulator 
-                                     // geometry component
-
-  G4VPhysicalVolume* physicalTreatmentRoom;
-  G4VPhysicalVolume* patientPhysicalVolume;
-  G4LogicalVolume* phantomLogicalVolume;
-  G4VPhysicalVolume* phantomPhysicalVolume;
-  
-  HadrontherapyDetectorMessenger* detectorMessenger; 
-  HadrontherapyMaterial* material;
-
   G4double phantomSizeX; 
   G4double phantomSizeY; 
   G4double phantomSizeZ;
-   
+
+  G4double detectorSizeX; 
+  G4double detectorSizeY; 
+  G4double detectorSizeZ;
+
+  G4ThreeVector phantomPosition, detectorPosition, detectorToPhantomPosition; //  phantom center, detector center, detector to phantom relative position
+
+  G4double sizeOfVoxelAlongX; 
+  G4double sizeOfVoxelAlongY; 
+  G4double sizeOfVoxelAlongZ; 
+
   G4int numberOfVoxelsAlongX; 
   G4int numberOfVoxelsAlongY;
   G4int numberOfVoxelsAlongZ;  
+
+  G4Box* phantom;
+  G4Box* detector;
+
 };
 #endif

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: HadrontherapyDetectorMessenger.hh; May 2005
-// ----------------------------------------------------------------------------
-//                 GEANT 4 - Hadrontherapy example
-// ----------------------------------------------------------------------------
-// Code developed by:
-//
-// G.A.P. Cirrone(a)*, F. Di Rosa(a), S. Guatelli(b), G. Russo(a)
-// 
-// (a) Laboratori Nazionali del Sud 
-//     of the INFN, Catania, Italy
-// (b) INFN Section of Genova, Genova, Italy
-// 
-// * cirrone@lns.infn.it
-// ----------------------------------------------------------------------------
+// HadrontherapyDetectorMessenger.hh;
+// See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy//
+
 #ifndef HadrontherapyDetectorMessenger_h
 #define HadrontherapyDetectorMessenger_h 1
@@ -48 +37 @@
 class G4UIcmdWithADoubleAndUnit;
 class G4UIcmdWithAString;
+class G4UIcmdWith3VectorAndUnit;
 
 class HadrontherapyDetectorMessenger: public G4UImessenger
@@ -55 +45 @@
   ~HadrontherapyDetectorMessenger();
     
-    void SetNewValue(G4UIcommand*, G4String);
+  void SetNewValue(G4UIcommand*, G4String);
     
 private:
 
-  // Pointer to the detector component
+  // Pointer to the phantom/detector 
   HadrontherapyDetectorConstruction* hadrontherapyDetector;
-  
-  G4UIdirectory* modulatorDir; // Control of the modulator 
-  G4UIdirectory* beamLineDir;  // Control of the beam line
- 
-  G4UIdirectory* rangeShifterDir; 
-  // Control of the range shifter component of the beam line
 
-  G4UIdirectory* firstScatteringFoilDir;
-  // Control of the first scattering foil component of the beam line
-  
-  G4UIdirectory* secondScatteringFoilDir;
-  // Control of the first scattering foil component of the beam line
-  
-  G4UIdirectory* rangeStopperDir;
-  // Control of the range stopper component of the beam line
-  
-  G4UIdirectory* finalCollimatorDir;
-  // Control of the final collimator component of the beam line
-  
-  G4UIcmdWithADoubleAndUnit* modulatorAngleCmd;
-  // UI command to rotate the modulator wheel
+  G4UIdirectory *changeThePhantomDir,  *changeTheDetectorDir; 
 
-  G4UIcmdWithAString*   rangeShifterMatCmd;
-  // UI command to set the material of the rangeShifter component of 
-  // the beam line 
-
-  G4UIcmdWithADoubleAndUnit* rangeShifterXSizeCmd;
-  // UI command to set half of the X size of the rangeShifter component of 
-  // the beam line 
-
-  G4UIcmdWithADoubleAndUnit* rangeShifterXPositionCmd;
-  // UI command to change the X position of the rangeShifter component of 
-  // the beam line 
-
-  G4UIcmdWithADoubleAndUnit* firstScatteringFoilXSizeCmd;
-  // UI command to set half X size of the first scattering foil of 
-  // the beam line 
-
-  G4UIcmdWithADoubleAndUnit* secondScatteringFoilXSizeCmd;
-  // UI command to set half X size of the second scattering foil 
-  // the beam line 
-
-  G4UIcmdWithADoubleAndUnit* outerRadiusStopperCmd;
-  // UI command to set the outer radius of the range stopper component of 
-  // the beam line 
-
-  G4UIcmdWithADoubleAndUnit* innerRadiusFinalCollimatorCmd;
-  // UI command to set the inner radius of the final collimator component of 
-  // the beam line 
+  G4UIcmdWith3VectorAndUnit *changeThePhantomSizeCmd,
+                            *changeThePhantomPositionCmd, 
+                            *changeTheDetectorSizeCmd, 
+                            *changeTheDetectorToPhantomPositionCmd,
+                            *changeTheDetectorVoxelCmd;
 };
 #endif

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// ----------------------------------------------------------------------------
-//                 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
-// --------------------------------------------------------------
+// HadrontherapyDummySD.hh
+// See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
 
 #ifndef HadrontherapyDummySD_h

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: HadrontherapyEventAction.hh; May 2005
-// ----------------------------------------------------------------------------
-//                 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
-// --------------------------------------------------------------
+// HadrontherapyEventAction.hh; 
+// See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
 
 #ifndef HadrontherapyEventAction_h
@@ -46 +34 @@
 
 class HadrontherapyMatrix;
+class HadrontherapyEventActionMessenger;
 
 class HadrontherapyEventAction : public G4UserEventAction
 {
 public:
-  HadrontherapyEventAction(HadrontherapyMatrix*);
+  HadrontherapyEventAction();
   ~HadrontherapyEventAction();
 
@@ -56 +45 @@
   void BeginOfEventAction(const G4Event*);
   void EndOfEventAction(const G4Event*);
-    
+
+  void SetPrintModulo(G4int val)
+  {
+    printModulo = val;
+  };
+
+  void SetDrawFlag(G4String val)
+  {
+    drawFlag = val;
+  };
+
 private: 
   G4String drawFlag; //Visualisation flag
   G4int hitsCollectionID;
   HadrontherapyMatrix *matrix; 
+  G4int printModulo;  
+  HadrontherapyEventActionMessenger* pointerEventMessenger;
 };
 

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: HadrontherapyMatrix.hh; May 2005
-// ----------------------------------------------------------------------------
-//                 GEANT 4 - Hadrontherapy example
-// ----------------------------------------------------------------------------
-// Code developed by:
-//
-// G.A.P. Cirrone(a)*, 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
-// ----------------------------------------------------------------------------
+// HadrontherapyMatrix.hh;
+// See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
 
 #ifndef HadrontherapyMatrix_H
@@ -43 +31 @@
 
 #include "globals.hh"
+#include <vector>
 
 // The information: energy deposit and position in the phantom
@@ -49 +38 @@
 class HadrontherapyMatrix 
 {
+private:
+  HadrontherapyMatrix(G4int voxelX, G4int voxelY, G4int voxelZ); //<--- this is supposed to be a singleton
+
 public:
-  HadrontherapyMatrix();
+
   ~HadrontherapyMatrix();
+// Get object instance only
+  static HadrontherapyMatrix* getInstance();
+// Make & Get instance
+  static HadrontherapyMatrix* getInstance(G4int nX, G4int nY, G4int nZ);
+ 
+  void flush();
  
   void Initialize(); 
   // All the elements of the matrix are initialised to zero
- 
+   
   void Fill(G4int i, G4int j, G4int k, G4double energyDeposit);
   // The matrix is filled with the energy deposit 
@@ -64 +62 @@
   // 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
 
 private:
+
+  static HadrontherapyMatrix* instance;
   G4int numberVoxelX;
   G4int numberVoxelY;

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: HadrontherapyParticles.hh; May 2005
-// ----------------------------------------------------------------------------
-//                 GEANT 4 - Hadrontherapy example
-// ----------------------------------------------------------------------------
-// Code developed by:
-//
-// G.A.P. Cirrone(a)*, 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
-// ----------------------------------------------------------------------------
+// HadrontherapyParticles.hh;
+// See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
 
 #ifndef HADRONTHERAPYPARTICLES_HH

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: HadrontherapyPhysicsList.hh,v 1.0
-// ----------------------------------------------------------------------------
-//                 GEANT 4 - Hadrontherapy example
-// ----------------------------------------------------------------------------
-// Code developed by:
-//
-// G.A.P. Cirrone(a)*, 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
-// ----------------------------------------------------------------------------
-#ifndef HADRONTHERAPYPHYSICSLIST_H
-#define HADRONTHERAPYPHYSICSLIST_H 1
+// HadrontherapyPhysicsList.hh
+// See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
+
+#ifndef HadrontherapyPhysicsList_h
+#define HadrontherapyPhysicsList_h 1
 
 #include "G4VModularPhysicsList.hh"
+#include "G4EmConfigurator.hh"
 #include "globals.hh"
 
+class G4VPhysicsConstructor;
+class HadrontherapyStepMax;
 class HadrontherapyPhysicsListMessenger;
 
@@ -49 +41 @@
 {
 public:
+
   HadrontherapyPhysicsList();
   virtual ~HadrontherapyPhysicsList();
 
-  virtual void SetCuts();
-  void AddPhysicsList(const G4String& name);  
-  
+  void ConstructParticle();
+
+  void SetCuts();
+  void SetCutForGamma(G4double);
+  void SetCutForElectron(G4double);
+  void SetCutForPositron(G4double);
+
+  void AddPhysicsList(const G4String& name);
+  void ConstructProcess();
+
+  void AddStepMax();
+  HadrontherapyStepMax* GetStepMaxProcess() {return stepMaxProcess;};
+  void AddPackage(const G4String& name);
+
 private:
-  G4bool decayIsRegistered;
-  G4bool emElectronIsRegistered;
-  G4bool emPositronIsRegistered;
-  G4bool emPhotonIsRegistered;
-  G4bool emIonIsRegistered;
-  G4bool emMuonIsRegistered;
-  G4bool hadrElasticHadronIonIsRegistered;
-  G4bool hadrInelasticPionIsRegistered;
-  G4bool hadrInelasticIonIsRegistered;
-  G4bool hadrInelasticProtonNeutronIsRegistered;
-  G4bool hadrAtRestMuonIsRegistered;
 
-  HadrontherapyPhysicsListMessenger* messenger;
+  G4EmConfigurator em_config;
+
+  G4double cutForGamma;
+  G4double cutForElectron;
+  G4double cutForPositron;
+
+  G4bool helIsRegisted;
+  G4bool bicIsRegisted;
+  G4bool biciIsRegisted;
+  G4bool locIonIonInelasticIsRegistered;
+  G4bool radioactiveDecayIsRegisted;
+
+  G4String                             emName;
+  G4VPhysicsConstructor*               emPhysicsList;
+  G4VPhysicsConstructor*               decPhysicsList;
+  std::vector<G4VPhysicsConstructor*>  hadronPhys;
+
+  HadrontherapyStepMax* stepMaxProcess;
+
+  HadrontherapyPhysicsListMessenger* pMessenger;
 };
 
 #endif
-
-
-

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: HadrontherapyPhisicsListMessenger.hh; May 2005
-// ----------------------------------------------------------------------------
-//                 GEANT 4 - Hadrontherapy example
-// ----------------------------------------------------------------------------
-// Code developed by:
-//
-// G.A.P. Cirrone(a)*, 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
-// ----------------------------------------------------------------------------
+// HadrontherapyPhysicsListsMessenger.hh
+// See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
 
-#ifndef HADRONTHERAPYPHYSICSLISTMESSENGER_HH
-#define HADRONTHERAPYPHYSICSLISTMESSENGER_HH 1
+#ifndef HadrontherapyPhysicsListMessenger_h
+#define HadrontherapyPhysicsListMessenger_h 1
 
 #include "globals.hh"
@@ -47 +35 @@
 class HadrontherapyPhysicsList;
 class G4UIdirectory;
-class G4UIcmdWithoutParameter;
-class G4UIcmdWithADouble;
 class G4UIcmdWithADoubleAndUnit;
-class G4UIcmdWithABool;
 class G4UIcmdWithAString;
 
-class HadrontherapyPhysicsListMessenger: public G4UImessenger {
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
-public:
+class HadrontherapyPhysicsListMessenger: public G4UImessenger
+{
+  public:
   
-  HadrontherapyPhysicsListMessenger(HadrontherapyPhysicsList* physList);
+    HadrontherapyPhysicsListMessenger(HadrontherapyPhysicsList* );
+   ~HadrontherapyPhysicsListMessenger();
+    
+    void SetNewValue(G4UIcommand*, G4String);
+    
+  private:
   
-  ~HadrontherapyPhysicsListMessenger();
-  
-  void SetNewValue(G4UIcommand*, G4String);
-  
-private:
-  
-  HadrontherapyPhysicsList* physicsList;   
-  G4UIdirectory* listDir;
-  G4UIcmdWithAString* physicsListCmd;
+  HadrontherapyPhysicsList* pPhysicsList;
+    
+  G4UIdirectory*             physDir;        
+  G4UIcmdWithADoubleAndUnit* gammaCutCmd;
+  G4UIcmdWithADoubleAndUnit* electCutCmd;
+  G4UIcmdWithADoubleAndUnit* protoCutCmd;    
+  G4UIcmdWithADoubleAndUnit* allCutCmd;    
+  G4UIcmdWithAString*        pListCmd;
+  G4UIcmdWithAString* packageListCmd;    
 };
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
 #endif
 
-
-
-
-
-
-
-

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: HadrontherapyPrimaryGeneratorAction.hh; May 2005
-// ----------------------------------------------------------------------------
-//                 GEANT 4 - Hadrontherapy example
-// ----------------------------------------------------------------------------
-// Code developed by:
-//
-// G.A.P. Cirrone(a)*, 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
-// ----------------------------------------------------------------------------
+// HadrontherapyPrimaryGeneratorAction.hh; May 2005
+// See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
 
 #ifndef HadrontherapyPrimaryGeneratorAction_h
@@ -68 +56 @@
   void SetsigmaMomentumY(G4double);
   void SetsigmaMomentumZ(G4double);
+  G4double GetmeanKineticEnergy(void);
     
 private:

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: HadrontherapyPrimaryGeneratorMessenger.hh; May 2005
-// ----------------------------------------------------------------------------
-//                 GEANT 4 - Hadrontherapy example
-// ----------------------------------------------------------------------------
-// Code developed by:
-//
-// G.A.P. Cirrone(a)*, 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
-// ----------------------------------------------------------------------------
+// HadrontherapyPrimaryGeneratorMessenger.hh;
+// See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
 
 #ifndef HadrontherapyPrimaryGeneratorMessenger_h

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: HadrontherapyRunAction.hh,v 3.0, September 2004;
-// ----------------------------------------------------------------------------
-//                 GEANT 4 - Hadrontherapy example
-// ----------------------------------------------------------------------------
-// Code developed by:
-//
-// G.A.P. Cirrone(a)*, F. Di Rosa(a), S. Guatelli(b), G. Russo(a)
-// 
-// (a) Laboratori Nazionali del Sud 
-//     of the INFN, Catania, Italy
-// (b) INFN Section of Genova, Genova, Italy
-// 
-// * cirrone@lns.infn.it
-// ----------------------------------------------------------------------------
+// HadrontherapyRunAction.hh
+// See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
 
 #ifndef HadrontherapyRunAction_h
@@ -51 +39 @@
 class HadrontherapyRunMessenger;
 class HadrontherapyFactory;
-class HadrontherapyFactoryIr;
-class HadrontherapyFactoryI;
-
-
 
 class HadrontherapyRunAction : public G4UserRunAction

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: HadrontherapyProtonSteppingAction.hh; May 2005
-// ----------------------------------------------------------------------------
-//                 GEANT 4 - Hadrontherapy example
-// ----------------------------------------------------------------------------
-// Code developed by:
-//
-// G.A.P. Cirrone(a)*, F. Di Rosa(a), S. Guatelli(b), G. Russo(a)
-// 
-// (a) Laboratori Nazionali del Sud 
-//     of the INFN, Catania, Italy
-// (b) INFN Section of Genova, Genova, Italy
-// 
-// * cirrone@lns.infn.it
-// ----------------------------------------------------------------------------
+// HadrontherapyProtonSteppingAction.hh;
+// See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
 
 #ifndef HadrontherapySteppingAction_h

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: HadrontherapyAnalisysManager.cc;  May 2005
-// ----------------------------------------------------------------------------
-//                 GEANT 4 - Hadrontherapy example
-// ----------------------------------------------------------------------------
-// Code developed by:
-//
-// G.A.P. Cirrone(a)*, F. Di Rosa(a), S. Guatelli(b), G. Russo(a)
-// 
-// (a) Laboratori Nazionali del Sud 
-//     of the INFN, Catania, Italy
-// (b) INFN Section of Genova, Genova, Italy
-// 
-// * cirrone@lns.infn.it
-// ----------------------------------------------------------------------------
-
-#ifdef  G4ANALYSIS_USE
+// $Id: HadrontherapyAnalisysManager.cc;
+// See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
+
 #include "HadrontherapyAnalysisManager.hh"
-
+#include "HadrontherapyMatrix.hh"
+#include "HadrontherapyAnalysisFileMessenger.hh"
+#include <time.h>
+#ifdef ANALYSIS_USE
 HadrontherapyAnalysisManager* HadrontherapyAnalysisManager::instance = 0;
 
-HadrontherapyAnalysisManager::HadrontherapyAnalysisManager() : 
-  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), ntuple(0),
-  ionTuple(0)
-{  
-}
-
-HadrontherapyAnalysisManager::~HadrontherapyAnalysisManager() 
-{ 
+#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);
+}
+#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;
@@ -100 +125 @@
   delete h1;
   h1 = 0;
-  
+
   delete tupFact;
   tupFact = 0;
@@ -112 +137 @@
   delete aFact;
   aFact = 0;
-}
-
+#endif
+#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;
+#endif
+}
+/////////////////////////////////////////////////////////////////////////////
 HadrontherapyAnalysisManager* HadrontherapyAnalysisManager::getInstance()
 {
@@ -120 +207 @@
 }
 
-void HadrontherapyAnalysisManager::book() 
-{
+/////////////////////////////////////////////////////////////////////////////
+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 file
-  G4String fileName = "hadrontherapy.hbk";
+  // 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;
 
@@ -136 +238 @@
 
   // Create the histograms with the enrgy deposit along the X axis
-  h1 = histFact -> createHistogram1D("10","slice, energy", 200, 0., 200. );
-
-  h2 = histFact -> createHistogram1D("20","Secondary protons - slice, energy", 200, 0., 200. );
- 
-  h3 = histFact -> createHistogram1D("30","Secondary neutrons - slice, energy", 200, 0., 200. );
-
-  h4 = histFact -> createHistogram1D("40","Secondary alpha - slice, energy", 200, 0., 200. );
-
-  h5 = histFact -> createHistogram1D("50","Secondary gamma - slice, energy", 200, 0., 200. );
-
-  h6 = histFact -> createHistogram1D("60","Secondary electron - slice, energy", 200, 0., 200. );
-
-  h7 = histFact -> createHistogram1D("70","Secondary triton - slice, energy", 200, 0., 200. );
-
-  h8 = histFact -> createHistogram1D("80","Secondary deuteron - slice, energy", 200, 0., 200. );
-
-  h9 = histFact -> createHistogram1D("90","Secondary pion - slice, energy", 200, 0., 200. );
- 
+  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
@@ -173 +279 @@
   G4String options2 = "";
   if (tupFact) ionTuple = tupFact -> create("2","2", columnNames2, options2);
-}
-
-void HadrontherapyAnalysisManager::FillEnergyDeposit(G4int i, 
-                                                     G4int j, 
+
+  // 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",
+                           70, 0., 500.);
+  histo16 = createHistogram1D("hydrogenEnergyAfterPhantom","Energy distribution of secondary helium fragments after the phantom",
+                           70, 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");
@@ -185 +328 @@
        G4int kSlice = ntuple -> findColumn("k");
        G4int iEnergy = ntuple -> findColumn("energy");
-      
+
        ntuple -> fill(iSlice,i);
-       ntuple -> fill(jSlice,j); 
+       ntuple -> fill(jSlice,j);
        ntuple -> fill(kSlice,k);
        ntuple -> fill(iEnergy, energy);     }
 
-  ntuple -> addRow(); 
-}
-
+  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);
-}
-
-void HadrontherapyAnalysisManager::genericIonInformation(G4int a, 
-                                                         G4double z, 
+#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) 
-{
+                                                         G4double energy)
+{
+#ifdef G4ANALYSIS_USE
   if (ionTuple)    {
        G4int aIndex = ionTuple -> findColumn("a");
        G4int zIndex = ionTuple -> findColumn("z");
-       G4int electronIndex = ionTuple -> findColumn("occupancy");  
+       G4int electronIndex = ionTuple -> findColumn("occupancy");
        G4int energyIndex = ionTuple -> findColumn("energy");
-      
+
        ionTuple -> fill(aIndex,a);
-      ionTuple -> fill(zIndex,z);  
-      ionTuple -> fill(electronIndex, electronOccupancy); 
+      ionTuple -> fill(zIndex,z);
+      ionTuple -> fill(electronIndex, electronOccupancy);
        ionTuple -> fill(energyIndex, energy);
      }
-   ionTuple -> addRow(); 
-}
-
-void HadrontherapyAnalysisManager::finish() 
-{  
+   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
-
-
-
-
-
-
-
-
-
-
-
+#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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: HadrontherapyDetectorConstruction.cc; Version 4.0 May 2005
-// ----------------------------------------------------------------------------
-//                 GEANT 4 - Hadrontherapy example
-// ----------------------------------------------------------------------------
-// Code developed by:
+// HadrontherapyDetectorConstruction.cc
 //
-// G.A.P. Cirrone(a)*, F. Di Rosa(a), S. Guatelli(b), G. Russo(a)
-// 
-// (a) Laboratori Nazionali del Sud 
-//     of the INFN, Catania, Italy
-// (b) INFN Section of Genova, Genova, Italy
-// 
-// * cirrone@lns.infn.it
-// ----------------------------------------------------------------------------
+// See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
 
 #include "G4SDManager.hh"
@@ -50 +39 @@
 #include "G4Colour.hh"
 #include "G4UserLimits.hh"
+#include "G4UnitsTable.hh"
 #include "G4VisAttributes.hh"
-#include "HadrontherapyPhantomROGeometry.hh"
+#include "G4NistManager.hh"
+#include "HadrontherapyDetectorROGeometry.hh"
 #include "HadrontherapyDetectorMessenger.hh"
-#include "HadrontherapyPhantomSD.hh"
+#include "HadrontherapyDetectorSD.hh"
 #include "HadrontherapyDetectorConstruction.hh"
-#include "HadrontherapyMaterial.hh"
-#include "HadrontherapyBeamLine.hh"
-#include "HadrontherapyModulator.hh"
-
-HadrontherapyDetectorConstruction::HadrontherapyDetectorConstruction()
-  : phantomSD(0), phantomROGeometry(0), beamLine(0), modulator(0),
-    physicalTreatmentRoom(0),
-    patientPhysicalVolume(0), 
-    phantomLogicalVolume(0), 
-    phantomPhysicalVolume(0)
-{
-  // Messenger to change parameters of the geometry
+#include "HadrontherapyMatrix.hh"
+
+/////////////////////////////////////////////////////////////////////////////
+HadrontherapyDetectorConstruction::HadrontherapyDetectorConstruction(G4VPhysicalVolume* physicalTreatmentRoom)
+  : motherPhys(physicalTreatmentRoom),
+    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
+{
+  // NOTE! that the HadrontherapyDetectorConstruction class
+  // does NOT inherit from G4VUserDetectorConstruction G4 class
+  // So the Construct() mandatory virtual method is inside another geometric class
+  // (like the passiveProtonBeamLIne, ...)
+
+  // Messenger to change parameters of the phantom/detector geometry
   detectorMessenger = new HadrontherapyDetectorMessenger(this);
 
-  material = new HadrontherapyMaterial();
-
-  // Phantom sizes
-  phantomSizeX = 20.*mm;
-  phantomSizeY = 20.*mm;
-  phantomSizeZ = 20.*mm;
-
-  // Number of the phantom voxels  
-  numberOfVoxelsAlongX = 200;
-  numberOfVoxelsAlongY = 200;
-  numberOfVoxelsAlongZ = 200;
-}
-
+  // 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);
+}
+
+/////////////////////////////////////////////////////////////////////////////
 HadrontherapyDetectorConstruction::~HadrontherapyDetectorConstruction()
 { 
-  delete material;
-  if (phantomROGeometry) delete phantomROGeometry;  
-  delete detectorMessenger;
-}
-
-G4VPhysicalVolume* HadrontherapyDetectorConstruction::Construct()
+    delete detectorROGeometry;// This should be safe in C++ even if the argument is a NULL pointer  
+    delete matrix;  
+    delete detectorMessenger;
+}
+
+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);
+    phantomLogicalVolume = new G4LogicalVolume(phantom, 
+                                             waterNist, 
+                                             "phantomLog", 0, 0, 0);
+  
+    phantomPhysicalVolume = new G4PVPlacement(0,
+                                            phantomPosition,
+                                            "phantomPhys",
+                                            phantomLogicalVolume,
+                                            motherPhys,
+                                            false,
+                                            0);
+
+// Visualisation attributes of the phantom
+    red = new G4VisAttributes(G4Colour(255/255., 0/255. ,0/255.));
+    red -> SetVisibility(true);
+    red -> SetForceSolid(true);
+//red -> SetForceWireframe(true);
+    phantomLogicalVolume -> SetVisAttributes(red);
+}
+
+/////////////////////////////////////////////////////////////////////////////
+void HadrontherapyDetectorConstruction::ConstructDetector()
+{
+  //-----------
+  // Detector
+  //-----------
+    G4bool isotopes =  false; 
+    G4Material* waterNist = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER", isotopes);
+    detector = new G4Box("Detector",detectorSizeX,detectorSizeY,detectorSizeZ);
+    detectorLogicalVolume = new G4LogicalVolume(detector,
+                                                waterNist,
+                                                "DetectorLog",
+                                                0,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);
+  
+// Visualisation attributes of the detector 
+    skyBlue = new G4VisAttributes( G4Colour(135/255. , 206/255. ,  235/255. ));
+    skyBlue -> SetVisibility(true);
+    skyBlue -> SetForceSolid(true);
+//skyBlue -> SetForceWireframe(true);
+    detectorLogicalVolume -> SetVisAttributes(skyBlue);
+  
+}
+/////////////////////////////////////////////////////////////////////////////
+void  HadrontherapyDetectorConstruction::ConstructSensitiveDetector(G4ThreeVector detectorToWorldPosition)
 {  
-  // Define the materials of the experimental set-up
-  material -> DefineMaterials();
-  
-  // Define the geometry components
-  ConstructBeamLine();
-  ConstructPhantom();
-  
-  // Set the sensitive detector where the energy deposit is collected
-  ConstructSensitiveDetector();
+    // Install new Sensitive Detector and ROGeometry 
+    delete detectorROGeometry; // this should be safe in C++ also if we have a NULL pointer
+
+    // Sensitive Detector and ReadOut geometry definition
+    G4SDManager* sensitiveDetectorManager = G4SDManager::GetSDMpointer();
+
+    G4String sensitiveDetectorName = "Detector"; 
+    if (!detectorSD)
+        {
+            // The sensitive detector is instantiated
+            detectorSD = new HadrontherapyDetectorSD(sensitiveDetectorName);
+        }
+    // The Read Out Geometry is instantiated
+    G4String ROGeometryName = "DetectorROGeometry";
+    detectorROGeometry = new HadrontherapyDetectorROGeometry(ROGeometryName,
+                                                            detectorToWorldPosition,
+                                                            detectorSizeX,
+                                                            detectorSizeY,
+                                                            detectorSizeZ,
+                                                            numberOfVoxelsAlongX,
+                                                            numberOfVoxelsAlongY,
+                                                            numberOfVoxelsAlongZ);
+
+    G4cout << "Instantiating new Read Out Geometry \"" << ROGeometryName << "\""<< G4endl;
+    // This will invoke Build() HadrontherapyDetectorROGeometry virtual method 
+    detectorROGeometry -> BuildROGeometry();
+    // Attach ROGeometry to SDetector
+    detectorSD -> SetROgeometry(detectorROGeometry);
+    //sensitiveDetectorManager -> Activate(sensitiveDetectorName, true);
+    if (!sensitiveDetectorManager -> FindSensitiveDetector(sensitiveDetectorName, false))
+        {
+            G4cout << "Registering new DetectorSD \"" << sensitiveDetectorName << "\""<< G4endl;
+            // Register user SD
+            sensitiveDetectorManager -> AddNewDetector(detectorSD);
+            // Attach SD to detector logical volume
+            detectorLogicalVolume -> SetSensitiveDetector(detectorSD);
+        }
+}
+
+/////////////////
+// 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)
+    {
+        if (sizeX > 2*detectorSizeX)
+        {
+            G4cout << "WARNING: Voxel X size must be smaller or equal than that of detector X" << G4endl;
+            return false;
+        }
+        // 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)
+    {
+        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)
+    {
+        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;
+
+    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);
+                   }
  
-  return physicalTreatmentRoom;
-}
-
-void HadrontherapyDetectorConstruction::ConstructBeamLine()
-{ 
-  G4Material* air = material -> GetMat("Air") ;
-  G4Material* water = material -> GetMat("Water");
-
-  // ---------------------
-  // Treatment room - World volume
-  //---------------------
-
-  // Treatment room sizes
-  const G4double worldX = 400.0 *cm;
-  const G4double worldY = 400.0 *cm;
-  const G4double worldZ = 400.0 *cm;
-
-  G4Box* treatmentRoom = new G4Box("TreatmentRoom",worldX,worldY,worldZ);
-
-  G4LogicalVolume* logicTreatmentRoom = new G4LogicalVolume(treatmentRoom, 
-                                                            air, 
-                                                            "logicTreatmentRoom", 
-                                                            0,0,0);
-
-
-
-  physicalTreatmentRoom = new G4PVPlacement(0,
-                                            G4ThreeVector(),
-                                            "physicalTreatmentRoom", 
-                                            logicTreatmentRoom, 
-                                            0,false,0);
-
-  G4double maxStepTreatmentRoom = 0.1 *mm;
-  logicTreatmentRoom -> SetUserLimits(new G4UserLimits(maxStepTreatmentRoom));
-
-  // The treatment room is invisible in the Visualisation
-  logicTreatmentRoom -> SetVisAttributes (G4VisAttributes::Invisible);
- 
-  beamLine = new HadrontherapyBeamLine(physicalTreatmentRoom);
-  beamLine -> HadrontherapyBeamLineSupport();
-  beamLine -> HadrontherapyBeamScatteringFoils();
-  beamLine -> HadrontherapyBeamCollimators();
-  beamLine -> HadrontherapyBeamMonitoring();
-  beamLine -> HadrontherapyBeamNozzle();
-  beamLine -> HadrontherapyBeamFinalCollimator();
-
-  modulator = new HadrontherapyModulator();
-  modulator -> BuildModulator(physicalTreatmentRoom);
-
-  // Patient - Mother volume of the phantom
-  G4Box* patient = new G4Box("patient",20 *cm, 20 *cm, 20 *cm);
-
-  G4LogicalVolume* patientLogicalVolume = new G4LogicalVolume(patient,
-                                                              water, 
-                                                              "patientLog", 0, 0, 0);
-
-  patientPhysicalVolume = new G4PVPlacement(0,G4ThreeVector(0., 0., 0.),
-                                            "patientPhys",
-                                            patientLogicalVolume,
-                                            physicalTreatmentRoom,
-                                            false,0);
- 
-  // Visualisation attributes of the patient 
-  G4VisAttributes * redWire = new G4VisAttributes(G4Colour(1. ,0. ,0.));
-  redWire -> SetVisibility(true);
-  redWire -> SetForceWireframe(true);
-  patientLogicalVolume -> SetVisAttributes(redWire); 
-}
-
-void HadrontherapyDetectorConstruction::ConstructPhantom()
-{
-  G4Colour  lightBlue   (0.0, 0.0, .75);
-
-  G4Material* water = material -> GetMat("Water");
-
-  //ComputeVoxelSize();
-
-  //----------------------
-  // Water phantom  
-  //----------------------
-  G4Box* phantom = new G4Box("Phantom",phantomSizeX,phantomSizeY,phantomSizeZ);
-
-  phantomLogicalVolume = new G4LogicalVolume(phantom,
-                                             water,
-                                             "PhantomLog",
-                                             0,0,0);
-
-  // Fixing the max step allowed in the phantom
-  G4double maxStep = 0.01 *mm;
-  phantomLogicalVolume -> SetUserLimits(new G4UserLimits(maxStep));
-
-  G4double phantomXtranslation = -180.*mm;
-  phantomPhysicalVolume = new G4PVPlacement(0,
-                                            G4ThreeVector(phantomXtranslation, 0.0 *mm, 0.0 *mm),
-                                            "PhantomPhys",
-                                            phantomLogicalVolume,
-                                            patientPhysicalVolume,
-                                            false,0);
- 
-  // Visualisation attributes of the phantom
-  G4VisAttributes* simpleBoxVisAttributes = new G4VisAttributes(lightBlue);
-  simpleBoxVisAttributes -> SetVisibility(true);
-  simpleBoxVisAttributes -> SetForceSolid(true);
-  phantomLogicalVolume -> SetVisAttributes(simpleBoxVisAttributes);
-
-  // **************
-  // Cut per Region
-  // **************
-  
-  // A smaller cut is fixed in the phantom to calculate the energy deposit with the
-  // required accuracy 
-  G4Region* aRegion = new G4Region("PhantomLog");
-  phantomLogicalVolume -> SetRegion(aRegion);
-  aRegion -> AddRootLogicalVolume(phantomLogicalVolume);
-}
-
-void  HadrontherapyDetectorConstruction::ConstructSensitiveDetector()
-{  
-  // Sensitive Detector and ReadOut geometry definition
-  G4SDManager* sensitiveDetectorManager = G4SDManager::GetSDMpointer();
-
-  G4String sensitiveDetectorName = "Phantom"; 
-
-  if(!phantomSD)
-    {
-      // The sensitive detector is instantiated
-      phantomSD = new HadrontherapyPhantomSD(sensitiveDetectorName);
-      
-      // The Read Out Geometry is instantiated
-      G4String ROGeometryName = "PhantomROGeometry";
-      phantomROGeometry = new HadrontherapyPhantomROGeometry(ROGeometryName,
-                                                             phantomSizeX,
-                                                             phantomSizeY,
-                                                             phantomSizeZ,
-                                                             numberOfVoxelsAlongX,
-                                                             numberOfVoxelsAlongY,
-                                                             numberOfVoxelsAlongZ);
-      phantomROGeometry -> BuildROGeometry();
-      phantomSD -> SetROgeometry(phantomROGeometry);
-      sensitiveDetectorManager -> AddNewDetector(phantomSD);
-      phantomLogicalVolume -> SetSensitiveDetector(phantomSD);
-    }
-}
-
-void HadrontherapyDetectorConstruction::SetModulatorAngle(G4double value)
-{  
-  modulator -> SetModulatorAngle(value);
-  G4RunManager::GetRunManager() -> GeometryHasBeenModified();
-}
-
-void HadrontherapyDetectorConstruction::SetRangeShifterXPosition(G4double value)
-{
-  beamLine -> SetRangeShifterXPosition(value);
-  G4RunManager::GetRunManager() -> GeometryHasBeenModified();
-}
-
-void HadrontherapyDetectorConstruction::SetRangeShifterXSize(G4double value)
-{
-  beamLine -> SetRangeShifterXSize(value);
-  G4RunManager::GetRunManager() -> GeometryHasBeenModified();
-}
-
-void HadrontherapyDetectorConstruction::SetFirstScatteringFoilSize(G4double value)
-{
-  beamLine -> SetFirstScatteringFoilXSize(value);  
-  G4RunManager::GetRunManager() -> GeometryHasBeenModified();
-}
-
-void HadrontherapyDetectorConstruction::SetSecondScatteringFoilSize(G4double value)
-{
-  beamLine -> SetSecondScatteringFoilXSize(value);
-  G4RunManager::GetRunManager() -> GeometryHasBeenModified();
-}
-
-void HadrontherapyDetectorConstruction::SetOuterRadiusStopper(G4double value)
-{
-  beamLine -> SetOuterRadiusStopper(value); 
-  G4RunManager::GetRunManager() -> GeometryHasBeenModified(); 
-}
-
-void HadrontherapyDetectorConstruction::SetInnerRadiusFinalCollimator(G4double value)
-{
-  beamLine -> SetInnerRadiusFinalCollimator(value);
-  G4RunManager::GetRunManager() -> GeometryHasBeenModified();
-}
-
-void HadrontherapyDetectorConstruction::SetRSMaterial(G4String materialChoice)
-{
-  beamLine -> SetRSMaterial(materialChoice);
-}
-
-
-
+
+    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;
+}

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: HadrontherapyDetectorMessenger.cc; May 2005
-// ----------------------------------------------------------------------------
-//                 GEANT 4 - Hadrontherapy example
-// ----------------------------------------------------------------------------
-// Code developed by:
-//
-// G.A.P. Cirrone(a)*, F. Di Rosa(a), S. Guatelli(b), G. Russo(a)
-// 
-// (a) Laboratori Nazionali del Sud 
-//     of the INFN, Catania, Italy
-// (b) INFN Section of Genova, Genova, Italy
-// 
-// * cirrone@lns.infn.it
-// ----------------------------------------------------------------------------
+// $Id: HadrontherapyDetectorMessenger.cc; 
+// See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
+
+
 #include "HadrontherapyDetectorMessenger.hh"
 #include "HadrontherapyDetectorConstruction.hh"
@@ -43 +33 @@
 #include "G4UIcmdWithADoubleAndUnit.hh"
 #include "G4UIcmdWithAString.hh"
+#include "G4UIcmdWith3VectorAndUnit.hh"
 
-HadrontherapyDetectorMessenger::HadrontherapyDetectorMessenger(
-                                                               HadrontherapyDetectorConstruction* detector)
+/////////////////////////////////////////////////////////////////////////////
+HadrontherapyDetectorMessenger::HadrontherapyDetectorMessenger(HadrontherapyDetectorConstruction* detector)
   :hadrontherapyDetector(detector)
-{ 
-  modulatorDir = new G4UIdirectory("/modulator/");
-  modulatorDir -> SetGuidance("Command to rotate the modulator wheel");
-  
-  beamLineDir = new G4UIdirectory("/beamLine/");
-  beamLineDir -> SetGuidance("set specification of range shifter");  
+{
+    // Change Phantom size
+    changeThePhantomDir = new G4UIdirectory("/changePhantom/");
+    changeThePhantomDir -> SetGuidance("Command to change the Phantom Size/position");
+    changeThePhantomSizeCmd = new G4UIcmdWith3VectorAndUnit("/changePhantom/size", this);
+    changeThePhantomSizeCmd -> SetGuidance("Insert sizes X Y and Z"
+                                           "\n   0 or negative values mean <<Don't change it!>>");
+    changeThePhantomSizeCmd -> SetParameterName("PhantomSizeAlongX", 
+                                                "PhantomSizeAlongY", 
+                                                "PhantomSizeAlongZ", false);
+    changeThePhantomSizeCmd -> SetDefaultUnit("mm");
+    changeThePhantomSizeCmd -> SetUnitCandidates("um mm cm"); 
+    changeThePhantomSizeCmd -> AvailableForStates(G4State_Idle);
 
-  rangeShifterDir = new G4UIdirectory("/beamLine/RangeShifter/");
-  rangeShifterDir -> SetGuidance("set specification of range shifter");  
 
-  firstScatteringFoilDir = new G4UIdirectory("/beamLine/ScatteringFoil1/");
-  firstScatteringFoilDir -> SetGuidance("set specification of first scattering foil");  
- 
-  secondScatteringFoilDir = new G4UIdirectory("/beamLine/ScatteringFoil2/");
-  secondScatteringFoilDir -> SetGuidance("set specification of second scattering foil");  
- 
-  rangeStopperDir = new G4UIdirectory("/beamLine/Stopper/");
-  rangeStopperDir -> SetGuidance("set specification of stopper");  
+    // Change Phantom position
+    changeThePhantomPositionCmd = new G4UIcmdWith3VectorAndUnit("/changePhantom/position", this);
+    changeThePhantomPositionCmd -> SetGuidance("Insert X Y and Z dimensions for the position of the center of the Phantom"
+                                                " respect to that of the \"World\""); 
+    changeThePhantomPositionCmd -> SetParameterName("PositionAlongX", 
+                                                    "PositionAlongY", 
+                                                    "PositionAlongZ", false);
+    changeThePhantomPositionCmd -> SetDefaultUnit("mm");
+    changeThePhantomPositionCmd -> SetUnitCandidates("mm cm m"); 
+    changeThePhantomPositionCmd -> AvailableForStates(G4State_Idle);
 
-  finalCollimatorDir = new G4UIdirectory("/beamLine/FinalCollimator/");
-  finalCollimatorDir -> SetGuidance("set specification of final collimator");  
 
-  modulatorAngleCmd = new G4UIcmdWithADoubleAndUnit("/modulator/angle",this);
-  modulatorAngleCmd -> SetGuidance("Set Modulator Angle");
-  modulatorAngleCmd -> SetParameterName("Size",false);
-  modulatorAngleCmd -> SetRange("Size>=0.");
-  modulatorAngleCmd -> SetUnitCategory("Angle");  
-  modulatorAngleCmd -> AvailableForStates(G4State_Idle);
-  
-  rangeShifterMatCmd = new G4UIcmdWithAString("/beamLine/RangeShifter/RSMat",this);
-  rangeShifterMatCmd -> SetGuidance("Set material of range shifter");
-  rangeShifterMatCmd -> SetParameterName("choice",false);
-  rangeShifterMatCmd -> AvailableForStates(G4State_Idle);
-  
-  rangeShifterXSizeCmd = new G4UIcmdWithADoubleAndUnit("/beamLine/RangeShifter/thickness",this);
-  rangeShifterXSizeCmd -> SetGuidance("Set half of the thickness of range shifter along X axis");
-  rangeShifterXSizeCmd -> SetParameterName("Size",false);
-  rangeShifterXSizeCmd -> SetDefaultUnit("mm");  
-  rangeShifterXSizeCmd -> SetUnitCandidates("mm cm m");  
-  rangeShifterXSizeCmd -> AvailableForStates(G4State_Idle);
-  
-  rangeShifterXPositionCmd = new G4UIcmdWithADoubleAndUnit("/beamLine/RangeShifter/position",this);
-  rangeShifterXPositionCmd -> SetGuidance("Set position of range shifter");
-  rangeShifterXPositionCmd -> SetParameterName("Size",false);
-  rangeShifterXPositionCmd -> SetDefaultUnit("mm");  
-  rangeShifterXPositionCmd -> SetUnitCandidates("mm cm m");  
-  rangeShifterXPositionCmd -> AvailableForStates(G4State_Idle);
-  
-  firstScatteringFoilXSizeCmd = new G4UIcmdWithADoubleAndUnit("/beamLine/ScatteringFoil1/thickness",this);
-  firstScatteringFoilXSizeCmd -> SetGuidance("Set hlaf thickness of first scattering foil");
-  firstScatteringFoilXSizeCmd -> SetParameterName("Size",false);
-  firstScatteringFoilXSizeCmd -> SetDefaultUnit("mm");  
-  firstScatteringFoilXSizeCmd -> SetUnitCandidates("mm cm m");  
-  firstScatteringFoilXSizeCmd -> AvailableForStates(G4State_Idle);
-  
-  secondScatteringFoilXSizeCmd = new G4UIcmdWithADoubleAndUnit("/beamLine/ScatteringFoil2/thickness",this);
-  secondScatteringFoilXSizeCmd -> SetGuidance("Set half thickness of second scattering foil");
-  secondScatteringFoilXSizeCmd -> SetParameterName("Size",false);
-  secondScatteringFoilXSizeCmd -> SetDefaultUnit("mm");  
-  secondScatteringFoilXSizeCmd -> SetUnitCandidates("mm cm m");  
-  secondScatteringFoilXSizeCmd -> AvailableForStates(G4State_Idle);
-  
-  outerRadiusStopperCmd = new G4UIcmdWithADoubleAndUnit("/beamLine/Stopper/outRadius",this);
-  outerRadiusStopperCmd -> SetGuidance("Set size of outer radius");
-  outerRadiusStopperCmd -> SetParameterName("Size",false);
-  outerRadiusStopperCmd -> SetDefaultUnit("mm");  
-  outerRadiusStopperCmd -> SetUnitCandidates("mm cm m");  
-  outerRadiusStopperCmd -> AvailableForStates(G4State_Idle);
-  
-  innerRadiusFinalCollimatorCmd = new G4UIcmdWithADoubleAndUnit("/beamLine/FinalCollimator/halfInnerRad",this);
-  innerRadiusFinalCollimatorCmd -> SetGuidance("Set size of inner radius ( max 21.5 mm)");
-  innerRadiusFinalCollimatorCmd -> SetParameterName("Size",false);
-  innerRadiusFinalCollimatorCmd -> SetDefaultUnit("mm");  
-  innerRadiusFinalCollimatorCmd -> SetUnitCandidates("mm cm m");  
-  innerRadiusFinalCollimatorCmd -> AvailableForStates(G4State_Idle);
+    //  Change detector size
+    changeTheDetectorDir = new G4UIdirectory("/changeDetector/");
+    changeTheDetectorDir -> SetGuidance("Command to change the Detector's Size/position/Voxels");
+    
+    changeTheDetectorSizeCmd = new G4UIcmdWith3VectorAndUnit("/changeDetector/size",this);
+    changeTheDetectorSizeCmd -> SetGuidance("Insert sizes for X Y and Z dimensions of the Detector"
+                                            "\n   0 or negative values mean <<Don't change it>>");
+    changeTheDetectorSizeCmd -> SetParameterName("DetectorSizeAlongX", "DetectorSizeAlongY", "DetectorSizeAlongZ", false);
+    changeTheDetectorSizeCmd -> SetDefaultUnit("mm");
+    changeTheDetectorSizeCmd -> SetUnitCandidates("um mm cm"); 
+    changeTheDetectorSizeCmd -> AvailableForStates(G4State_Idle);
+
+    //  Change the detector to phantom displacement
+    changeTheDetectorToPhantomPositionCmd = new G4UIcmdWith3VectorAndUnit("/changeDetector/displacement",this);
+    changeTheDetectorToPhantomPositionCmd -> SetGuidance("Insert X Y and Z displacements between Detector and Phantom"
+                                                         "\nNegative values mean <<Don't change it!>>"); 
+    changeTheDetectorToPhantomPositionCmd -> SetParameterName("DisplacementAlongX",
+                                                              "DisplacementAlongY", 
+                                                              "DisplacementAlongZ", false);
+    changeTheDetectorToPhantomPositionCmd -> SetDefaultUnit("mm");
+    changeTheDetectorToPhantomPositionCmd -> SetUnitCandidates("um mm cm"); 
+    changeTheDetectorToPhantomPositionCmd -> AvailableForStates(G4State_Idle);
+    
+    // Change voxels by its size
+    changeTheDetectorVoxelCmd = new G4UIcmdWith3VectorAndUnit("/changeDetector/voxelSize",this);
+    changeTheDetectorVoxelCmd -> SetGuidance("Insert Voxel sizes for X Y and Z dimensions"
+                                             "\n   0 or negative values mean <<Don't change it!>>");
+    changeTheDetectorVoxelCmd -> SetParameterName("VoxelSizeAlongX", "VoxelSizeAlongY", "VoxelSizeAlongZ", false);
+    changeTheDetectorVoxelCmd -> SetDefaultUnit("mm");
+    changeTheDetectorVoxelCmd -> SetUnitCandidates("um mm cm");
+    changeTheDetectorVoxelCmd -> AvailableForStates(G4State_Idle);
+   }
+
+/////////////////////////////////////////////////////////////////////////////
+HadrontherapyDetectorMessenger::~HadrontherapyDetectorMessenger()
+{
+    delete changeThePhantomDir; 
+    delete changeThePhantomSizeCmd; 
+    delete changeThePhantomPositionCmd; 
+    delete changeTheDetectorDir; 
+    delete changeTheDetectorSizeCmd; 
+    delete changeTheDetectorToPhantomPositionCmd; 
+    delete changeTheDetectorVoxelCmd; 
 }
 
-HadrontherapyDetectorMessenger::~HadrontherapyDetectorMessenger()
-{ 
-  delete innerRadiusFinalCollimatorCmd;  
-  delete outerRadiusStopperCmd;  
-  delete secondScatteringFoilXSizeCmd; 
-  delete firstScatteringFoilXSizeCmd; 
-  delete rangeShifterXPositionCmd;
-  delete rangeShifterXSizeCmd;
-  delete rangeShifterMatCmd;
-  delete modulatorAngleCmd;
-  delete finalCollimatorDir; 
-  delete rangeStopperDir;
-  delete secondScatteringFoilDir;
-  delete firstScatteringFoilDir; 
-  delete rangeShifterDir;  
-  delete beamLineDir; 
-  delete modulatorDir;   
+/////////////////////////////////////////////////////////////////////////////
+void HadrontherapyDetectorMessenger::SetNewValue(G4UIcommand* command,G4String newValue)
+{
+        
+  if( command == changeThePhantomSizeCmd)
+    {
+        G4ThreeVector size = changeThePhantomSizeCmd -> GetNew3VectorValue(newValue);
+        hadrontherapyDetector -> SetPhantomSize(size.getX(),size.getY(),size.getZ());
+    }
+  else if (command == changeThePhantomPositionCmd )
+  {
+         G4ThreeVector size = changeThePhantomPositionCmd -> GetNew3VectorValue(newValue);
+         hadrontherapyDetector -> SetPhantomPosition(size);
+  }
+  else if (command == changeTheDetectorSizeCmd)
+  {
+        G4ThreeVector size = changeTheDetectorSizeCmd  -> GetNew3VectorValue(newValue);
+        hadrontherapyDetector -> SetDetectorSize(size.getX(),size.getY(),size.getZ());
+  }
+  else if (command == changeTheDetectorToPhantomPositionCmd)
+  {
+        G4ThreeVector size = changeTheDetectorToPhantomPositionCmd-> GetNew3VectorValue(newValue);
+        hadrontherapyDetector -> SetDetectorToPhantomPosition(size);
+  }
+  else if (command == changeTheDetectorVoxelCmd)
+  {
+        G4ThreeVector size = changeTheDetectorVoxelCmd  -> GetNew3VectorValue(newValue);
+        hadrontherapyDetector -> SetNumberOfVoxelBySize(size.getX(),size.getY(),size.getZ());
+  }
 }
-
-void HadrontherapyDetectorMessenger::SetNewValue(G4UIcommand* command,G4String newValue)
-{ 
-  if( command == modulatorAngleCmd )
-    { hadrontherapyDetector -> SetModulatorAngle
-           (modulatorAngleCmd -> GetNewDoubleValue(newValue));}
-
-  if( command == rangeShifterMatCmd )
-    { hadrontherapyDetector -> SetRSMaterial(newValue);}
-
-  if( command == rangeShifterXSizeCmd )
-    { hadrontherapyDetector -> SetRangeShifterXSize
-            (rangeShifterXSizeCmd -> GetNewDoubleValue(newValue));}
-
-  if( command == rangeShifterXPositionCmd )
-    { hadrontherapyDetector -> SetRangeShifterXPosition
-                  (rangeShifterXPositionCmd -> GetNewDoubleValue(newValue));}
-
-  if( command == firstScatteringFoilXSizeCmd )
-    { hadrontherapyDetector -> SetFirstScatteringFoilSize
-                  (firstScatteringFoilXSizeCmd -> GetNewDoubleValue(newValue));}
-
-  if( command == secondScatteringFoilXSizeCmd )
-    { hadrontherapyDetector -> SetSecondScatteringFoilSize
-                  (secondScatteringFoilXSizeCmd -> GetNewDoubleValue(newValue));}
-
-  if( command == outerRadiusStopperCmd )
-    { hadrontherapyDetector -> SetOuterRadiusStopper(
-                    outerRadiusStopperCmd -> GetNewDoubleValue(newValue));}
-
-  if( command == innerRadiusFinalCollimatorCmd )
-    { hadrontherapyDetector -> SetInnerRadiusFinalCollimator
-                  (innerRadiusFinalCollimatorCmd -> GetNewDoubleValue(newValue));}
-}
-

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: HadrontherapyEventAction.cc; May 2005
-// ----------------------------------------------------------------------------
-//                 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
-// --------------------------------------------------------------
+// $Id: HadrontherapyEventAction.cc; 
+// See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
+
 #include "G4Event.hh"
 #include "G4EventManager.hh"
@@ -48 +37 @@
 #include "G4VVisManager.hh"
 #include "HadrontherapyEventAction.hh"
-#include "HadrontherapyPhantomHit.hh"
-#include "HadrontherapyPhantomSD.hh"
+#include "HadrontherapyDetectorHit.hh"
+#include "HadrontherapyDetectorSD.hh"
 #include "HadrontherapyDetectorConstruction.hh"
 #include "HadrontherapyMatrix.hh"
+#include "HadrontherapyEventActionMessenger.hh"
 
-HadrontherapyEventAction::HadrontherapyEventAction(HadrontherapyMatrix* matrixPointer) :
-  drawFlag("all" )
+/////////////////////////////////////////////////////////////////////////////
+HadrontherapyEventAction::HadrontherapyEventAction() :
+  drawFlag("all" ),printModulo(1000), pointerEventMessenger(0)
 { 
   hitsCollectionID = -1;
-  matrix = matrixPointer; 
+  pointerEventMessenger = new HadrontherapyEventActionMessenger(this);
 }
 
+/////////////////////////////////////////////////////////////////////////////
 HadrontherapyEventAction::~HadrontherapyEventAction()
 {
+ delete pointerEventMessenger;
 }
 
-void HadrontherapyEventAction::BeginOfEventAction(const G4Event* )
-{
- G4SDManager* pSDManager = G4SDManager::GetSDMpointer();
- if(hitsCollectionID == -1)
-        hitsCollectionID = pSDManager -> GetCollectionID("HadrontherapyPhantomHitsCollection");
+/////////////////////////////////////////////////////////////////////////////
+void HadrontherapyEventAction::BeginOfEventAction(const G4Event* evt)
+{ 
+  G4int evtNb = evt->GetEventID();
+  
+  //printing survey
+  if (evtNb%printModulo == 0)
+    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;
-
+  return;
   G4HCofThisEvent* HCE = evt -> GetHCofThisEvent();
-  HadrontherapyPhantomHitsCollection* CHC = NULL; 
  
   if(HCE)
-    CHC = (HadrontherapyPhantomHitsCollection*)(HCE -> GetHC(hitsCollectionID));
-  
-  if(CHC)
-    {
-      if(matrix)
-        {
-          // Fill the matrix with the information: voxel and associated energy deposit 
-          // in the phantom at the end of the event
+  {
+    HadrontherapyDetectorHitsCollection* CHC = (HadrontherapyDetectorHitsCollection*)(HCE -> GetHC(hitsCollectionID));
+    if(CHC)
+     {
+           matrix = HadrontherapyMatrix::getInstance();  
+       if(matrix)
+          { 
+              // Fill the matrix with the information: voxel and associated energy deposit 
+          // in the detector at the end of the event
 
           G4int HitCount = CHC -> entries();
@@ -98 +98 @@
               matrix -> Fill(i, j, k, energyDeposit/MeV);              
             }
-        }
+          }
     }
-
+  }
   // Extract the trajectories and draw them in the visualisation
 

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: HadrontherapyPhantomSD.cc; May 2005
-// ----------------------------------------------------------------------------
-//                 GEANT 4 - Hadrontherapy example
-// ----------------------------------------------------------------------------
-// Code developed by:
-//
-// G.A.P. Cirrone(a)*, 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
-// ----------------------------------------------------------------------------
+// $Id: HadrontherapyMatrix.cc;
+// See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy//
 
 #include "HadrontherapyMatrix.hh"
@@ -44 +32 @@
 #include <fstream>
 
-HadrontherapyMatrix::HadrontherapyMatrix()
+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)
 {  
-  // Number of the voxels of the phantom
-  numberVoxelX = 200;
-  numberVoxelY = 200;
-  numberVoxelZ = 200; 
- 
-  // Create the matrix
+// 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!");
 }
 
 HadrontherapyMatrix::~HadrontherapyMatrix()
 {
-  delete[] matrix;
+    delete[] matrix;
+}
+
+void HadrontherapyMatrix::flush()
+{
+        if(matrix)
+        for(int i=0; i<numberVoxelX*numberVoxelY*numberVoxelZ; i++)
+        {
+          matrix[i] = 0;
+        }
 }
 void HadrontherapyMatrix::Initialize()
 { 
-  // Initialise the elemnts of the matrix to zero
-
+// 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++)
+           {
+              for(G4int k = 0; k < numberVoxelZ; k++)
 
-            matrix[(i*numberVoxelY+j)*numberVoxelZ+k] = 0.;
-        }
+              matrix[Index(i,j,k)] = 0.;
+           } 
     }
 }
@@ -78 +95 @@
 {
   if (matrix)
-    matrix[(i*numberVoxelY+j)*numberVoxelZ+k] += energyDeposit;
+    matrix[Index(i,j,k)] += energyDeposit;
   
-  // Store the energy deposit in the matrix elemnt corresponding 
-  // to the phantom voxel  
+  // Store the energy deposit in the matrix element corresponding 
+  // to the phantom voxel  i, j, k
 }
 
@@ -95 +112 @@
   if (matrix)
     {  
-        for(G4int l = 0; l < numberVoxelZ; l++) 
-          {
-            k = l;
-            
-            for(G4int m = 0; m < numberVoxelY; m++) 
-              { 
-                j = m * numberVoxelZ + k; 
+      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++)
@@ -108 +128 @@
                     if(matrix[i] != 0)
                       { 
-                                        
-#ifdef G4ANALYSIS_USE   
+                        ofs << n << '\t' << m << '\t' <<
+                          k << '\t' << matrix[i] << G4endl;
+                        
+#ifdef ANALYSIS_USE
                         HadrontherapyAnalysisManager* analysis = 
-                          HadrontherapyAnalysisManager::getInstance();
+                        HadrontherapyAnalysisManager::getInstance();
                         analysis -> FillEnergyDeposit(n, m, k, matrix[i]);
                         analysis -> BraggPeak(n, matrix[i]);
 #endif
                       }
-                  }       
+}       
               }
           }

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: HadrontherapyModulator.cc; May 2005
-// ----------------------------------------------------------------------------
-//                 GEANT 4 - Hadrontherapy example
-// ----------------------------------------------------------------------------
-// Code developed by:
-//
-// G.A.P. Cirrone(a)*, 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
-// ----------------------------------------------------------------------------
+// $Id: HadrontherapyModulator.cc; 
+// See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
 
 #include "G4Material.hh"
@@ -51 +39 @@
 #include "G4VisAttributes.hh"
 #include "G4Colour.hh"
-#include "HadrontherapyMaterial.hh"
 #include "HadrontherapyModulator.hh"
-#include "HadrontherapyMaterial.hh"
 #include "G4Transform3D.hh"
 #include "G4ios.hh"
 #include <fstream>
 #include "G4RunManager.hh"
+#include "G4NistManager.hh"
 
 HadrontherapyModulator::HadrontherapyModulator():physiMotherMod(0),
@@ -141 +128 @@
   rm -> rotateY(phi); 
 }
-
+/////////////////////////////////////////////////////////////////////////////
 HadrontherapyModulator::~HadrontherapyModulator() 
 {
   delete rm;
 }
-
+/////////////////////////////////////////////////////////////////////////////
 void HadrontherapyModulator::BuildModulator(G4VPhysicalVolume* motherVolume)
 {
-
-
-  //Materials used for the modulator wheel
-  HadrontherapyMaterial* material = new HadrontherapyMaterial();
-
-  G4Material* Mod0Mater = material -> GetMat("Air");
-  G4Material* ModMater = material -> GetMat("Air");
-  delete material;
-
+  G4bool isotopes = false;
+  G4Material* airNist =  G4NistManager::Instance()->FindOrBuildMaterial("G4_AIR", isotopes);
+
+  G4Material* Mod0Mater = airNist;
+  G4Material* ModMater = airNist;
+ 
   G4double innerRadiusOfTheTube = 2.5 *cm;
   G4double outerRadiusOfTheTube = 9.5 *cm;
@@ -163 +147 @@
 
   // Mother of the modulator wheel  
-
-  G4ThreeVector positionMotherMod = G4ThreeVector(-2260.50 *mm, 30 *mm, 50 *mm);
+  G4ThreeVector positionMotherMod = G4ThreeVector(-1960.50 *mm, 30 *mm, 50 *mm);
  
   G4Box* solidMotherMod = new G4Box("MotherMod", 12 *cm, 12 *cm, 12 *cm);
  
   G4LogicalVolume * logicMotherMod = new G4LogicalVolume(solidMotherMod, Mod0Mater,"MotherMod",0,0,0);
-
- 
 
   physiMotherMod = new G4PVPlacement(rm,positionMotherMod,  "MotherMod", 
@@ -200 +181 @@
   logicMod0 = new G4LogicalVolume(solidMod0, Mod0Mater, "Mod0",0,0,0);
   
-  
   physiMod0 = new G4PVPlacement(G4Transform3D(rm2, positionMod0), 
                                 logicMod0,    
@@ -212 +192 @@
   // First modulator sclice
   //----------------------------------------------------------
- 
  
   G4double startAngleOfTheTube1 = 54.267*deg;

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: HadrontherapyParticles.cc; May 2005
-// ----------------------------------------------------------------------------
-//                 GEANT 4 - Hadrontherapy example
-// ----------------------------------------------------------------------------
-// Code developed by:
-//
-// G.A.P. Cirrone(a)*, 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
-// ----------------------------------------------------------------------------
+// $Id: HadrontherapyParticles.cc; 
+// See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
+
 #include "HadrontherapyParticles.hh"
 #include "G4ParticleDefinition.hh"

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: HadrontherapyPhysicsList.cc,v 1.0
-// ----------------------------------------------------------------------------
-//                 GEANT 4 - Hadrontherapy example
-// ----------------------------------------------------------------------------
-// Code developed by:
-//
-// G.A.P. Cirrone(a)*, 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 "globals.hh"
-#include "G4ProcessManager.hh"
-#include "G4Region.hh"
-#include "G4RegionStore.hh"
-#include "G4ParticleDefinition.hh"
-#include "G4ParticleTypes.hh"
-#include "G4ParticleTable.hh"
+// HadrontherapyPhysicsList.cc
+// See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
+
+// This class provides all the physic models that can be activated inside Hadrontherapy;
+// Each model can be setted via macro commands;
+// Inside Hadrontherapy the models can be activate with three different complementar methods:
+//
+// 1. Use of the *Packages*.
+//    Packages (that are contained inside the
+//    Geant4 distribution at $G4INSTALL/source/physics_lists/lists) provide a full set
+//    of models (both electromagnetic and hadronic).
+//    The User can use this method simply add the line /physic/addPackage <nameOfPackage>
+//    in his/her macro file. No other action is required.
+//    For Hadrontherapy applications we suggest the use of the QGSP_BIC package
+//    for proton beams. The same can be used
+//    also for ligth ion beam.
+//    Example of use of package can be found in the packageQGSP_BIC.mac file.
+//
+// 2. Use of the *Physic Lists*.
+//    Physic lists are also already ready to use inside the Geant4 distribution
+//    ($G4INSTALL/source/physics_lists/builders). To use them the simple
+//    /physic/addPhysics <nameOfPhysicList> command must be used in the macro.
+//    In Hadrontherapy we provide physics list to activate Electromagnetic,
+//    Hadronic elastic and Hadronic inelastic models.
+//
+//    For Hadrontherapy we suggest the use of:
+//
+//    /physic/addPhysic/emstandard_option3 (electromagnetic model)
+//    /physic/addPhysic/QElastic (hadronic elastic model)
+//    /physic/addPhysic/binary (hadronic inelastic models for proton and neutrons)
+//    /physic/addPhysic/binary_ion (hadronic inelastic models for ions)
+//
+//    Example of the use of physics lists can be found in the macro files included in the
+//    'macro' folder .
+//
+// 3. Use of a *local* physics. In this case the models are implemented in local files
+//    contained in the Hadrontherapy folder. The use of local physic is recommended
+//    to more expert Users.
+//    We provide as local, only the LocalStandardICRU73EmPhysic.cc (an Elecromagnetic
+//    implementation containing the new ICRU73 data table for ions stopping powers)
+//    and the LocalIonIonInelasticPhysic.cc (physic list to use for the ion-ion interaction
+//    case)
+//    The *local* physics can be activated with the same /physic/addPhysic <nameOfPhysic> command;
+//
+//    While Packages approch must be used exclusively, Physics List and Local physics can
+//    be activated, if necessary, contemporaneously in the same simulation run.
+//
+//    AT MOMENT, IF ACCURATE RESULTS ARE NEDED, WE STRONGLY RECOMMEND THE USE OF THE MACROS:
+//    proton_therapy.mac: use of the built-in Geant4 physics list for proton beams)
+//    ion_therapy.mac   : use of mixed combination of native Geant4 physic lists
+//                        and local physic for ion-ion enelastic processes)
+
 #include "HadrontherapyPhysicsList.hh"
 #include "HadrontherapyPhysicsListMessenger.hh"
-#include "HadrontherapyParticles.hh"
-#include "Decay.hh"
-#include "EMPhotonStandard.hh"
-#include "EMPhotonEPDL.hh"
-#include "EMPhotonPenelope.hh"
-#include "EMElectronStandard.hh"
-#include "EMElectronEEDL.hh"
-#include "EMElectronPenelope.hh"
-#include "EMPositronStandard.hh"
-#include "EMPositronPenelope.hh"
-#include "EMMuonStandard.hh"
-#include "EMHadronIonLowEICRU49.hh"
-#include "EMHadronIonLowEZiegler1977.hh"
-#include "EMHadronIonLowEZiegler1985.hh"
-#include "EMHadronIonStandard.hh"
-#include "HIProtonNeutronPrecompound.hh"
-#include "HIProtonNeutronBertini.hh"
-#include "HIProtonNeutronBinary.hh"
-#include "HIProtonNeutronLEP.hh"
-#include "HIProtonNeutronPrecompoundGEM.hh"
-#include "HIProtonNeutronPrecompoundFermi.hh"
-#include "HIProtonNeutronPrecompoundGEMFermi.hh"
-#include "HIPionBertini.hh"
-#include "HIPionLEP.hh"
-#include "HIIonLEP.hh"
-#include "HEHadronIonLElastic.hh"
-#include "HEHadronIonBertiniElastic.hh"
-#include "HEHadronIonQElastic.hh"
-#include "HEHadronIonUElastic.hh"
-#include "HRMuonMinusCapture.hh"
-
-
-HadrontherapyPhysicsList::HadrontherapyPhysicsList(): G4VModularPhysicsList(),
-                                                      decayIsRegistered(false),
-                                                      emElectronIsRegistered(false), 
-                                                      emPositronIsRegistered(false), 
-                                                      emPhotonIsRegistered(false), 
-                                                      emIonIsRegistered(false),
-                                                      emMuonIsRegistered(false),
-                                                      hadrElasticHadronIonIsRegistered(false),
-                                                      hadrInelasticPionIsRegistered(false),
-                                                      hadrInelasticIonIsRegistered(false),
-                                                      hadrInelasticProtonNeutronIsRegistered(false),
-                                                      hadrAtRestMuonIsRegistered(false)
-{
-  // The secondary production threshold is set to 10. mm
-  // for all the particles in all the experimental set-up
-  // The phantom is defined as a Geant4 Region. Here the cut is fixed to 0.001 mm
-  defaultCutValue = 0.01 * mm;
-
-  // Messenger: it is possible to activate physics processes and models interactively 
-  messenger = new HadrontherapyPhysicsListMessenger(this);
+#include "HadrontherapyStepMax.hh"
+#include "G4PhysListFactory.hh"
+#include "G4VPhysicsConstructor.hh"
+
+// Local physic directly implemented in the Hadronthrapy directory
+#include "LocalIonIonInelasticPhysic.hh"             // Physic dedicated to the ion-ion inelastic processes
+#include "LocalINCLIonIonInelasticPhysic.hh"         // Physic dedicated to the ion-ion inelastic processes using INCL/ABLA
+
+// Physic lists (contained inside the Geant4 distribution)
+#include "G4EmStandardPhysics_option3.hh"
+#include "G4EmLivermorePhysics.hh"
+#include "G4EmPenelopePhysics.hh"
+#include "G4DecayPhysics.hh"
+#include "G4HadronElasticPhysics.hh"
+#include "G4HadronDElasticPhysics.hh"
+#include "G4HadronHElasticPhysics.hh"
+#include "G4HadronQElasticPhysics.hh"
+#include "G4HadronInelasticQBBC.hh"
+#include "G4IonBinaryCascadePhysics.hh"
+#include "G4Decay.hh"
+
+#include "G4LossTableManager.hh"
+#include "G4UnitsTable.hh"
+#include "G4ProcessManager.hh"
+
+#include "G4IonFluctuations.hh"
+#include "G4IonParametrisedLossModel.hh"
+#include "G4EmProcessOptions.hh"
+
+#include "G4RadioactiveDecayPhysics.hh"
+
+/////////////////////////////////////////////////////////////////////////////
+HadrontherapyPhysicsList::HadrontherapyPhysicsList() : G4VModularPhysicsList()
+{
+  G4LossTableManager::Instance();
+  defaultCutValue = 1.*mm;
+  cutForGamma     = defaultCutValue;
+  cutForElectron  = defaultCutValue;
+  cutForPositron  = defaultCutValue;
+
+  helIsRegisted  = false;
+  bicIsRegisted  = false;
+  biciIsRegisted = false;
+  locIonIonInelasticIsRegistered = false;
+  radioactiveDecayIsRegisted = false;
+
+  stepMaxProcess  = 0;
+
+  pMessenger = new HadrontherapyPhysicsListMessenger(this);
 
   SetVerboseLevel(1);
 
-  // Register all the particles involved in the experimental set-up
-  RegisterPhysics( new HadrontherapyParticles("particles") );
-}
-
+  // EM physics
+  emPhysicsList = new G4EmStandardPhysics_option3(1);
+  emName = G4String("emstandard_opt3");
+
+  // Deacy physics and all particles
+  decPhysicsList = new G4DecayPhysics();
+}
+
+/////////////////////////////////////////////////////////////////////////////
 HadrontherapyPhysicsList::~HadrontherapyPhysicsList()
-{ 
-  delete messenger;
-}
-
+{
+  delete pMessenger;
+  delete emPhysicsList;
+  delete decPhysicsList;
+  for(size_t i=0; i<hadronPhys.size(); i++) {delete hadronPhys[i];}
+}
+
+/////////////////////////////////////////////////////////////////////////////
+void HadrontherapyPhysicsList::AddPackage(const G4String& name)
+{
+  G4PhysListFactory factory;
+  G4VModularPhysicsList* phys =factory.GetReferencePhysList(name);
+  G4int i=0;
+  const G4VPhysicsConstructor* elem= phys->GetPhysics(i);
+  G4VPhysicsConstructor* tmp = const_cast<G4VPhysicsConstructor*> (elem);
+  while (elem !=0)
+        {
+          RegisterPhysics(tmp);
+          elem= phys->GetPhysics(++i) ;
+          tmp = const_cast<G4VPhysicsConstructor*> (elem);
+        }
+}
+
+/////////////////////////////////////////////////////////////////////////////
+void HadrontherapyPhysicsList::ConstructParticle()
+{
+  decPhysicsList->ConstructParticle();
+}
+
+/////////////////////////////////////////////////////////////////////////////
+void HadrontherapyPhysicsList::ConstructProcess()
+{
+  // transportation
+  //
+  AddTransportation();
+
+  // electromagnetic physics list
+  //
+  emPhysicsList->ConstructProcess();
+  em_config.AddModels();
+
+  // decay physics list
+  //
+  decPhysicsList->ConstructProcess();
+
+  // hadronic physics lists
+  for(size_t i=0; i<hadronPhys.size(); i++) {
+    hadronPhys[i]->ConstructProcess();
+  }
+
+  // step limitation (as a full process)
+  //
+  AddStepMax();
+}
+
+/////////////////////////////////////////////////////////////////////////////
 void HadrontherapyPhysicsList::AddPhysicsList(const G4String& name)
 {
-  G4cout << "Adding PhysicsList component " << name << G4endl;
-  
-
-  // ****************
-  // *** A. DECAY ***
-  // ****************
-
-
-  if (name == "Decay") 
-    {
-      if (decayIsRegistered) 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name  
-                 << " cannot be registered ---- decay List already existing" 
-                 << G4endl;
-        } 
-      else 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name 
-                 << " is registered" << G4endl;
-          RegisterPhysics( new Decay(name) );
-          decayIsRegistered = true;
-        }
-    }
-
-
-  // ***************************************
-  // *** B. ELECTROMAGNETIC INTERACTIONS ***
-  // ***************************************
-
-
-  // ***************
-  // *** Photons ***
-  // ***************
-  
-  // *** Option I: Standard 
-
-  if (name == "EM-Photon-Standard") 
-    {
-      if (emPhotonIsRegistered) 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name  
-                 << " cannot be registered ---- photon List already existing" 
-                 << G4endl;
-        } 
-      else 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name 
-                 << " is registered" << G4endl;
-          RegisterPhysics( new EMPhotonStandard(name) );
-          emPhotonIsRegistered = true;
-        }
-    }
-
-
-  // *** Option II: Low Energy based on the Livermore libraries
-
-  if (name == "EM-Photon-EPDL") 
-    {
-      if (emPhotonIsRegistered) 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name  
-                 << " cannot be registered ---- photon List already existing"
-                 << G4endl;
-        } 
-      else 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name 
-                 << " is registered" << G4endl;
-          RegisterPhysics( new EMPhotonEPDL(name) );
-          emPhotonIsRegistered = true;
-        }
-    } 
-
-
-  // *** Option III: Low Energy Penelope
-
-  if (name == "EM-Photon-Penelope")
-    {
-      if (emPhotonIsRegistered) 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name  
-                 << " cannot be registered ---- photon List already existing" 
-                 << G4endl;
-        } 
-      else 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name 
-                 << " is registered" << G4endl;
-          RegisterPhysics( new EMPhotonPenelope(name) );
-          emPhotonIsRegistered = true;
-        }
-    }
-
-
-  // *****************
-  // *** Electrons ***
-  // *****************
-  
-  // *** Option I: Standard
-
-  if (name == "EM-Electron-Standard") 
-    {
-      if (emElectronIsRegistered) 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name  
-                 << " cannot be registered ---- electron List already existing" 
-                 << G4endl;
-        } 
-      else 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name 
-                 << " is registered" << G4endl;
-          RegisterPhysics( new EMElectronStandard(name) );        
-          emElectronIsRegistered = true;
-        }
-    }
-
-
-  // *** Option II: Low Energy based on the Livermore libraries
-
-  if (name == "EM-Electron-EEDL") 
-    {
-      if (emElectronIsRegistered) 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name  
-                 << " cannot be registered ---- electron List already existing"                  
-                 << G4endl;
-        } 
-      else 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name 
-                 << " is registered" << G4endl;
-          RegisterPhysics( new EMElectronEEDL(name) );
-          emElectronIsRegistered = true;
-        }
-    } 
-
-
-  // *** Option III: Low Energy Penelope
-
-  if (name == "EM-Electron-Penelope")
-    {
-      if (emElectronIsRegistered) 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name 
-                 << " cannot be registered ---- electron List already existing"                  
-                 << G4endl;
-        } 
-      else 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name 
-                 << " is registered" << G4endl;
-          RegisterPhysics( new EMElectronPenelope(name) );
-          emElectronIsRegistered = true;
-        }
-    }
-
-
-  // *****************
-  // *** Positrons ***
-  // *****************
-
-  // *** Option I: Standard
-  if (name == "EM-Positron-Standard") 
-    {
-      if (emPositronIsRegistered) 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name  
-                 << " cannot be registered ---- positron List already existing"                  
-                 << G4endl;
-        } 
-      else 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name 
-                 << " is registered" << G4endl;
-          RegisterPhysics( new EMPositronStandard(name) );
-          emPositronIsRegistered = true;
-        }
-    }
-
-
-  // *** Option II: Low Energy Penelope
-
-  if (name == "EM-Positron-Penelope") 
-    {
-      if (emPositronIsRegistered) 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name  
-                 << " cannot be registered ---- positron List already existing"   
-                 << G4endl;
-        } 
-      else 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name 
-                 << " is registered" << G4endl;
-          RegisterPhysics( new EMPositronPenelope(name) );
-          emPositronIsRegistered = true;
-        }
-    }
- 
-
-  // ************************
-  // *** Hadrons and Ions ***
-  // ************************
-
-  // *** Option I: Low Energy with ICRU49 stopping power parametrisation
-  
-  if (name == "EM-HadronIon-LowE") 
-    {
-      if (emIonIsRegistered) 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name  
-                 << " cannot be registered ---- proton List already existing" 
-                 << G4endl;
-        } 
-      else 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name 
-                 << " is registered" << G4endl;
-          RegisterPhysics( new EMHadronIonLowEICRU49(name) );
-          emIonIsRegistered = true;
-        }
-    }
-
-
-  // *** Option II: Low Energy with Ziegler 1977 stopping power parametrisation
-
-  if (name == "EM-HadronIon-LowEZiegler1977") 
-    {
-      if (emIonIsRegistered) 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name  
-                 << " cannot be registered ---- proton List already existing" 
-                 << G4endl;
-        } 
-      else 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name 
-                 << " is registered" << G4endl;
-          RegisterPhysics( new EMHadronIonLowEZiegler1977(name) );
-          emIonIsRegistered = true;
-        }
-    }
-
-
-  // *** Option III: Low Energy with Ziegler 1985 stopping power parametrisation
-
-  if (name == "EM-HadronIon-LowEZiegler1985") 
-    {
-      if (emIonIsRegistered) 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name  
-                 << " cannot be registered ---- proton List already existing" 
-                 << G4endl;
-        } 
-      else 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name 
-                 << " is registered" << G4endl;
-          RegisterPhysics( new EMHadronIonLowEZiegler1985(name) );
-          emIonIsRegistered = true;
-        }
-    }
-
-
-  // *** Option IV: Standard
-
-  if (name == "EM-HadronIon-Standard") 
-    {
-      if (emIonIsRegistered) 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name  
-                 << " cannot be registered ---- ion List already existing" 
-                 << G4endl;
-        } 
-      else 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name 
-                 << " is registered" << G4endl;
-          RegisterPhysics( new EMHadronIonStandard(name) );
-          emIonIsRegistered = true;
-        }
-    }
-
-
-  // *************
-  // *** Muons ***
-  // *************
-
-  // *** Option I: Standard 
-
-  if (name == "EM-Muon-Standard") 
-    {
-      if (emMuonIsRegistered) 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name  
-                 << " cannot be registered ---- decay List already existing" 
-                 << G4endl;
-        } 
-      else 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name 
-                 << " is registered" << G4endl;
-          RegisterPhysics( new EMMuonStandard(name) );
-          emMuonIsRegistered = true;
-        }
-    }
-
-
-  // ********************************
-  // *** C. HADRONIC INTERACTIONS ***
-  // ********************************
-
-
-  // ******************************
-  // *** C.1. ELASTIC PROCESSES ***
-  // ******************************
-
-
-  // ************************
-  // *** Hadrons and Ions ***
-  // ************************
-
-  // *** Option I: GHEISHA like LEP model
-
-  if (name == "HadronicEl-HadronIon-LElastic") 
-    {
-      if (hadrElasticHadronIonIsRegistered) 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name  
-                 << " cannot be registered ---- hadronic Elastic Scattering List already existing" 
-                 << G4endl;
-        } 
-      else 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name 
-                 << " is registered" << G4endl;
-          RegisterPhysics( new HEHadronIonLElastic(name) );
-          hadrElasticHadronIonIsRegistered = true;
-        }
-    }
-  
-
-  // *** Option II: Bertini model
-
-  if (name == "HadronicEl-HadronIon-Bert") 
-    {
-      if (hadrElasticHadronIonIsRegistered) 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name  
-                 << " cannot be registered ---- hadronic Elastic Scattering List already existing" 
-                 << G4endl;
-        } 
-      else 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name 
-                 << " is registered" << G4endl;
-          RegisterPhysics( new HEHadronIonBertiniElastic(name) );
-          hadrElasticHadronIonIsRegistered = true;
-        }
-    }
-
-
-  // *** Option III: Process G4QElastic
-
-  if (name == "HadronicEl-HadronIon-QElastic") 
-    {
-      if (hadrElasticHadronIonIsRegistered) 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name  
-                 << " cannot be registered ---- hadronic Elastic Scattering List already existing" 
-                 << G4endl;
-        } 
-      else 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name 
-                 << " is registered" << G4endl;
-          RegisterPhysics( new HEHadronIonQElastic(name) );
-          hadrElasticHadronIonIsRegistered = true;
-        }
-    }
-
-
-  // *** Option III: Process G4UHadronElasticProcess
-
-  if (name == "HadronicEl-HadronIon-UElastic") 
-    {
-     if (hadrElasticHadronIonIsRegistered) 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name  
-                 << " cannot be registered ---- hadronic Elastic Scattering List already existing" 
-                 << G4endl;
-        } 
-     else 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name 
-                 << " is registered" << G4endl;
-          RegisterPhysics( new HEHadronIonUElastic(name) );
-          hadrElasticHadronIonIsRegistered = true;
-        }
-    }
-   
-   
-  // ********************************
-  // *** C.2. INELASTIC PROCESSES ***
-  // ********************************
-
-
-  // *************
-  // *** Pions ***
-  // *************
-
-  // *** Option I: Bertini model
-
-  if (name == "HadronicInel-Pion-Bertini") 
-    {
-      if (hadrInelasticPionIsRegistered) 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name  
-                 << " cannot be registered ---- decay List already existing" 
-                 << G4endl;
-        } 
-      else 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name 
-                 << " is registered" << G4endl;
-          RegisterPhysics( new HIPionBertini(name) );
-          hadrInelasticPionIsRegistered = true;
-        }
-    }
-
-
-  // *** Option II: GHEISHA like LEP model
-
-  if (name == "HadronicInel-Pion-LEP") 
-    {
-      if (hadrInelasticPionIsRegistered) 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name  
-                 << " cannot be registered ---- decay List already existing" 
-                 << G4endl;
-        } 
-      else 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name 
-                 << " is registered" << G4endl;
-          RegisterPhysics( new HIPionLEP(name) );
-          hadrInelasticPionIsRegistered = true;
-        }
-    }
-
-
-  // ************
-  // *** Ions ***
-  // ************
-
-  // *** Option I: GHEISHA like LEP model
-
-  if (name == "HadronicInel-Ion-LEP") 
-    {
-      if (hadrInelasticIonIsRegistered) 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name  
-                 << " cannot be registered ---- decay List already existing" 
-                 << G4endl;
-        } 
-      else 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name 
-                 << " is registered" << G4endl;
-          RegisterPhysics( new HIIonLEP(name) );
-          hadrInelasticIonIsRegistered = true;
-        }
-    }
-
-
-  // *************************
-  // *** Protons, Neutrons ***
-  // *************************
-
-  // *** Option I: GHEISHA like LEP model
-
-  if (name == "HadronicInel-ProtonNeutron-LEP") 
-    {
-      if (hadrInelasticProtonNeutronIsRegistered) 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name  
-                 << " cannot be registered ---- decay List already existing" 
-                 << G4endl;
-        } 
-      else 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name 
-                 << " is registered" << G4endl;
-          RegisterPhysics( new HIProtonNeutronLEP(name) );
-          hadrInelasticProtonNeutronIsRegistered = true;
-        }
-    }
-
-
-  // *** Option II: Bertini Cascade Model
-
-  if (name == "HadronicInel-ProtonNeutron-Bert") 
-    {
-      if (hadrInelasticProtonNeutronIsRegistered) 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name  
-                 << " cannot be registered ---- decay List already existing" 
-                 << G4endl;
-        } 
-      else 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name 
-                 << " is registered" << G4endl;
-          RegisterPhysics( new HIProtonNeutronBertini(name) );
-          hadrInelasticProtonNeutronIsRegistered = true;
-        }
-    }
-
-
-  // *** Option III: Binary Cascade Model
-
-  if (name == "HadronicInel-ProtonNeutron-Bin") 
-    {
-      if (hadrInelasticProtonNeutronIsRegistered) 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name  
-                 << " cannot be registered ---- decay List already existing" 
-                 << G4endl;
-        } 
-      else 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name 
-                 << " is registered" << G4endl;
-          RegisterPhysics( new HIProtonNeutronBinary(name) );
-          hadrInelasticProtonNeutronIsRegistered = true;
-        }
-    }
-
-
-  // *** Option IV: Precompound Model combined with Default Evaporation
-
-  if (name == "HadronicInel-ProtonNeutron-Prec") 
-    {
-      if (hadrInelasticProtonNeutronIsRegistered) 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name  
-                 << " cannot be registered ---- decay List already existing" 
-                 << G4endl;
-        } 
-      else 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name 
-                 << " is registered" << G4endl;
-          RegisterPhysics( new HIProtonNeutronPrecompound(name) );
-          hadrInelasticProtonNeutronIsRegistered = true;
-        }
-    }
-
-
-  // *** Option V: Precompound Model combined with GEM Evaporation
-
-  if (name == "HadronicInel-ProtonNeutron-PrecGEM") 
-    {
-      if (hadrInelasticProtonNeutronIsRegistered) 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name  
-                 << " cannot be registered ---- decay List already existing" 
-                 << G4endl;
-        } 
-      else 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name 
-                 << " is registered" << G4endl;
-          RegisterPhysics( new HIProtonNeutronPrecompoundGEM(name) );
-          hadrInelasticProtonNeutronIsRegistered = true;
-        }
-    }
-
-
-  // *** Option VI: Precompound Model combined with default Evaporation 
-  //                and Fermi Break-up model
-
-  if (name == "HadronicInel-ProtonNeutron-PrecFermi") 
-    {
-      if (hadrInelasticProtonNeutronIsRegistered) 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name  
-                 << " cannot be registered ---- decay List already existing" 
-                 << G4endl;
-        } 
-      else 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name 
-                 << " is registered" << G4endl;
-          RegisterPhysics( new HIProtonNeutronPrecompoundFermi(name) );
-          hadrInelasticProtonNeutronIsRegistered = true;
-        }
-    }
-
-
-  // *** Option VII: Precompound Model combined with GEM Evaporation 
-  //                 and Fermi Break-up model
-
-  if (name == "HadronicInel-ProtonNeutron-PrecGEMFermi") 
-    {
-      if (hadrInelasticProtonNeutronIsRegistered) 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name  
-                 << " cannot be registered ---- decay List already existing" 
-                 << G4endl;
-        } 
-      else 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name 
-                 << " is registered" << G4endl;
-          RegisterPhysics( new HIProtonNeutronPrecompoundGEMFermi(name) );
-          hadrInelasticProtonNeutronIsRegistered = true;
-        }
-    }   
-
-
-  // ******************************
-  // *** C.3. AT-REST PROCESSES ***
-  // ******************************
-
-
-  // **************
-  // *** Muons- ***
-  // **************
-
-  // *** Option I: Muon Minus Capture at Rest
-
-  if (name == "HadronicAtRest-MuonMinus-Capture") 
-    {
-      if (hadrAtRestMuonIsRegistered) 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name  
-                 << " cannot be registered ---- decay List already existing" 
-                 << G4endl;
-        } 
-      else 
-        {
-          G4cout << "HadrontherapyPhysicsList::AddPhysicsList: " << name 
-                 << " is registered" << G4endl;
-          RegisterPhysics( new HRMuonMinusCapture(name) );
-          hadrAtRestMuonIsRegistered = true;
-        }
-    }
-
-}
-
+
+  if (verboseLevel>1) {
+    G4cout << "PhysicsList::AddPhysicsList: <" << name << ">" << G4endl;
+  }
+  if (name == emName) return;
+
+  /////////////////////////////////////////////////////////////////////////////
+  //   ELECTROMAGNETIC MODELS
+  /////////////////////////////////////////////////////////////////////////////
+
+  if (name == "standard_opt3") {
+    emName = name;
+    delete emPhysicsList;
+    emPhysicsList = new G4EmStandardPhysics_option3();
+    G4cout << "THE FOLLOWING ELECTROMAGNETIC PHYSICS LIST HAS BEEN ACTIVATED: G4EmStandardPhysics_option3" << G4endl;
+
+  } else if (name == "LowE_Livermore") {
+    emName = name;
+    delete emPhysicsList;
+    emPhysicsList = new G4EmLivermorePhysics();
+    G4cout << "THE FOLLOWING ELECTROMAGNETIC PHYSICS LIST HAS BEEN ACTIVATED: G4EmLivermorePhysics" << G4endl;
+
+  } else if (name == "LowE_Penelope") {
+    emName = name;
+    delete emPhysicsList;
+    emPhysicsList = new G4EmPenelopePhysics();
+    G4cout << "THE FOLLOWING ELECTROMAGNETIC PHYSICS LIST HAS BEEN ACTIVATED: G4EmLivermorePhysics" << G4endl;
+
+    /////////////////////////////////////////////////////////////////////////////
+    //   HADRONIC MODELS
+    /////////////////////////////////////////////////////////////////////////////
+  } else if (name == "elastic" && !helIsRegisted) {
+    G4cout << "THE FOLLOWING HADRONIC ELASTIC PHYSICS LIST HAS BEEN ACTIVATED: G4HadronElasticPhysics()" << G4endl;
+    hadronPhys.push_back( new G4HadronElasticPhysics());
+    helIsRegisted = true;
+
+  } else if (name == "DElastic" && !helIsRegisted) {
+    hadronPhys.push_back( new G4HadronDElasticPhysics());
+    helIsRegisted = true;
+
+  } else if (name == "HElastic" && !helIsRegisted) {
+    hadronPhys.push_back( new G4HadronHElasticPhysics());
+    helIsRegisted = true;
+
+  } else if (name == "QElastic" && !helIsRegisted) {
+    hadronPhys.push_back( new G4HadronQElasticPhysics());
+    helIsRegisted = true;
+
+  } else if (name == "binary" && !bicIsRegisted) {
+    hadronPhys.push_back(new G4HadronInelasticQBBC());
+    bicIsRegisted = true;
+    G4cout << "THE FOLLOWING HADRONIC INELASTIC PHYSICS LIST HAS BEEN ACTIVATED: G4HadronInelasticQBBC()" << G4endl;
+
+  } else if (name == "binary_ion" && !biciIsRegisted) {
+    hadronPhys.push_back(new G4IonBinaryCascadePhysics());
+    biciIsRegisted = true;
+
+  } else if (name == "local_ion_ion_inelastic" && !locIonIonInelasticIsRegistered) {
+    hadronPhys.push_back(new LocalIonIonInelasticPhysic());
+    locIonIonInelasticIsRegistered = true;
+
+  } else if (name == "local_incl_ion_ion_inelastic" && !locIonIonInelasticIsRegistered) {
+    hadronPhys.push_back(new LocalINCLIonIonInelasticPhysic());
+    locIonIonInelasticIsRegistered = true;
+
+  } else if (name == "radioactive_decay" && !radioactiveDecayIsRegisted ) {
+    hadronPhys.push_back(new G4RadioactiveDecayPhysics());
+    radioactiveDecayIsRegisted = true;
+
+  } else {
+
+    G4cout << "PhysicsList::AddPhysicsList: <" << name << ">"
+           << " is not defined"
+           << G4endl;
+  }
+}
+
+/////////////////////////////////////////////////////////////////////////////
+void HadrontherapyPhysicsList::AddStepMax()
+{
+  // Step limitation seen as a process
+  stepMaxProcess = new HadrontherapyStepMax();
+
+  theParticleIterator->reset();
+  while ((*theParticleIterator)()){
+    G4ParticleDefinition* particle = theParticleIterator->value();
+    G4ProcessManager* pmanager = particle->GetProcessManager();
+
+    if (stepMaxProcess->IsApplicable(*particle) && pmanager)
+      {
+        pmanager ->AddDiscreteProcess(stepMaxProcess);
+      }
+  }
+}
+
+/////////////////////////////////////////////////////////////////////////////
 void HadrontherapyPhysicsList::SetCuts()
-{  
-  // Set the threshold of production equal to the defaultCutValue
-  // in the experimental set-up
-  G4VUserPhysicsList::SetCutsWithDefault();
-     
-  G4double lowlimit=250*eV;
-  G4ProductionCutsTable::GetProductionCutsTable() ->SetEnergyRange(lowlimit, 100.*GeV);
-  // Definition of a smaller threshold of production in the phantom region
-  // where high accuracy is required in the energy deposit calculation
-
-  G4String regionName = "PhantomLog";
-  G4Region* region = G4RegionStore::GetInstance()->GetRegion(regionName);
-  G4ProductionCuts* cuts = new G4ProductionCuts ;
-  G4double regionCut = 0.01*mm;
-  cuts -> SetProductionCut(regionCut,G4ProductionCuts::GetIndex("gamma"));
-  cuts -> SetProductionCut(regionCut,G4ProductionCuts::GetIndex("e-"));
-  cuts -> SetProductionCut(regionCut,G4ProductionCuts::GetIndex("e+"));
-  cuts -> SetProductionCut(regionCut,G4ProductionCuts::GetIndex("proton"));
-  cuts -> SetProductionCut(regionCut,G4ProductionCuts::GetIndex("genericIons"));
-  region -> SetProductionCuts(cuts);
+{
+
+  if (verboseLevel >0){
+    G4cout << "PhysicsList::SetCuts:";
+    G4cout << "CutLength : " << G4BestUnit(defaultCutValue,"Length") << G4endl;
+  }
+
+  // set cut values for gamma at first and for e- second and next for e+,
+  // because some processes for e+/e- need cut values for gamma
+  SetCutValue(cutForGamma, "gamma");
+  SetCutValue(cutForElectron, "e-");
+  SetCutValue(cutForPositron, "e+");
 
   if (verboseLevel>0) DumpCutValuesTable();
 }
 
-
+/////////////////////////////////////////////////////////////////////////////
+void HadrontherapyPhysicsList::SetCutForGamma(G4double cut)
+{
+  cutForGamma = cut;
+  SetParticleCuts(cutForGamma, G4Gamma::Gamma());
+}
+
+/////////////////////////////////////////////////////////////////////////////
+void HadrontherapyPhysicsList::SetCutForElectron(G4double cut)
+{
+  cutForElectron = cut;
+  SetParticleCuts(cutForElectron, G4Electron::Electron());
+}
+
+/////////////////////////////////////////////////////////////////////////////
+void HadrontherapyPhysicsList::SetCutForPositron(G4double cut)
+{
+  cutForPositron = cut;
+  SetParticleCuts(cutForPositron, G4Positron::Positron());
+}

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: HadrontherapyPhisicsListMessenger.cc; May 2005
-// ----------------------------------------------------------------------------
-//                 GEANT 4 - Hadrontherapy example
-// ----------------------------------------------------------------------------
-// Code developed by:
-//
-// G.A.P. Cirrone(a)*, 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
-// ----------------------------------------------------------------------------
+// HadrontherapyPhysicsListMessenger.cc
+// See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
+
 #include "HadrontherapyPhysicsListMessenger.hh"
+
 #include "HadrontherapyPhysicsList.hh"
 #include "G4UIdirectory.hh"
-#include "G4UIcmdWithoutParameter.hh"
-#include "G4UIcmdWithADouble.hh"
 #include "G4UIcmdWithADoubleAndUnit.hh"
-#include "G4UIcmdWithABool.hh"
 #include "G4UIcmdWithAString.hh"
 
-HadrontherapyPhysicsListMessenger::HadrontherapyPhysicsListMessenger(HadrontherapyPhysicsList * physList)
-:physicsList(physList)
-{  
- listDir = new G4UIdirectory("/physics/");
-  // Building modular PhysicsList
+/////////////////////////////////////////////////////////////////////////////
+HadrontherapyPhysicsListMessenger::HadrontherapyPhysicsListMessenger(HadrontherapyPhysicsList* pPhys)
+:pPhysicsList(pPhys)
+{
+  physDir = new G4UIdirectory("/physic/");
+  physDir->SetGuidance("Commands to activate physics models and set cuts");
+   
+  gammaCutCmd = new G4UIcmdWithADoubleAndUnit("/physic/setGCut",this);  
+  gammaCutCmd->SetGuidance("Set gamma cut.");
+  gammaCutCmd->SetParameterName("Gcut",false);
+  gammaCutCmd->SetUnitCategory("Length");
+  gammaCutCmd->SetRange("Gcut>0.0");
+  gammaCutCmd->AvailableForStates(G4State_PreInit,G4State_Idle);
 
- physicsListCmd = new G4UIcmdWithAString("/physics/addPhysics",this);  
- physicsListCmd->SetGuidance("Add chunks of PhysicsList.");
- physicsListCmd->SetParameterName("physList",false);
- physicsListCmd->AvailableForStates(G4State_PreInit);
+  electCutCmd = new G4UIcmdWithADoubleAndUnit("/physic/setECut",this);  
+  electCutCmd->SetGuidance("Set electron cut.");
+  electCutCmd->SetParameterName("Ecut",false);
+  electCutCmd->SetUnitCategory("Length");
+  electCutCmd->SetRange("Ecut>0.0");
+  electCutCmd->AvailableForStates(G4State_PreInit,G4State_Idle);
+  
+  protoCutCmd = new G4UIcmdWithADoubleAndUnit("/physic/setPCut",this);  
+  protoCutCmd->SetGuidance("Set positron cut.");
+  protoCutCmd->SetParameterName("Pcut",false);
+  protoCutCmd->SetUnitCategory("Length");
+  protoCutCmd->SetRange("Pcut>0.0");
+  protoCutCmd->AvailableForStates(G4State_PreInit,G4State_Idle);  
+
+  allCutCmd = new G4UIcmdWithADoubleAndUnit("/physic/setCuts",this);  
+  allCutCmd->SetGuidance("Set cut for all.");
+  allCutCmd->SetParameterName("cut",false);
+  allCutCmd->SetUnitCategory("Length");
+  allCutCmd->SetRange("cut>0.0");
+  allCutCmd->AvailableForStates(G4State_PreInit,G4State_Idle);  
+
+  pListCmd = new G4UIcmdWithAString("/physic/addPhysics",this);  
+  pListCmd->SetGuidance("Add physics list.");
+  pListCmd->SetParameterName("PList",false);
+  pListCmd->AvailableForStates(G4State_PreInit);  
+
+  packageListCmd = new G4UIcmdWithAString("/physic/addPackage",this);
+  packageListCmd->SetGuidance("Add physics package.");
+  packageListCmd->SetParameterName("package",false);
+  packageListCmd->AvailableForStates(G4State_PreInit);
 }
 
+/////////////////////////////////////////////////////////////////////////////
 HadrontherapyPhysicsListMessenger::~HadrontherapyPhysicsListMessenger()
 {
-  delete physicsListCmd;
-  delete listDir;
+  delete gammaCutCmd;
+  delete electCutCmd;
+  delete protoCutCmd;
+  delete allCutCmd;
+  delete pListCmd;
+  delete physDir;    
+  delete packageListCmd;
 }
 
-void HadrontherapyPhysicsListMessenger::SetNewValue(G4UIcommand* command,G4String newValue)
-{
- if (command == physicsListCmd)
-    { physicsList->AddPhysicsList(newValue);} 
-}
+/////////////////////////////////////////////////////////////////////////////
+void HadrontherapyPhysicsListMessenger::SetNewValue(G4UIcommand* command,
+                                          G4String newValue)
+{       
+  if( command == gammaCutCmd )
+   { pPhysicsList->SetCutForGamma(gammaCutCmd->GetNewDoubleValue(newValue));}
+     
+  if( command == electCutCmd )
+   { pPhysicsList->SetCutForElectron(electCutCmd->GetNewDoubleValue(newValue));}
+     
+  if( command == protoCutCmd )
+   { pPhysicsList->SetCutForPositron(protoCutCmd->GetNewDoubleValue(newValue));}
+
+  if( command == allCutCmd )
+    {
+      G4double cut = allCutCmd->GetNewDoubleValue(newValue);
+      pPhysicsList->SetCutForGamma(cut);
+      pPhysicsList->SetCutForElectron(cut);
+      pPhysicsList->SetCutForPositron(cut);
+    } 
+
+  if( command == pListCmd )
+   { pPhysicsList->AddPhysicsList(newValue);}
 
 
+  if( command == packageListCmd )
+   { pPhysicsList->AddPackage(newValue);}
 
 
-
-
+}

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: HadrontherapyPositronPrimaryGeneratorAction.cc; May 2005
+// HadrontherapyPrimarygeneratorAction.cc;
+// See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
 // ----------------------------------------------------------------------------
 //                 GEANT 4 - Hadrontherapy example
@@ -30 +31 @@
 // Code developed by:
 //
-// G.A.P. Cirrone(a)*, F. Di Rosa(a), S. Guatelli(b), G. Russo(a)
+// G.A.P. Cirrone(a)*, F.Romano(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
+//     of the INFN, Catania, Italy
 // 
 // * cirrone@lns.infn.it
-// ----------------------------------------------------------------------------
+//
+// ------------------------------------------------------------------------------
+
 #include "HadrontherapyPrimaryGeneratorAction.hh"
 #include "HadrontherapyPrimaryGeneratorMessenger.hh"
@@ -45 +47 @@
 #include "G4ParticleDefinition.hh"
 #include "Randomize.hh"
+#include "HadrontherapyAnalysisManager.hh"
 
 HadrontherapyPrimaryGeneratorAction::HadrontherapyPrimaryGeneratorAction()
@@ -75 +78 @@
 
   // Define the energy of primary particles:
-  // gaussian distribution with mean energy = 64.55 *MeV
-  // and sigma = 300.0 *keV
-  G4double defaultMeanKineticEnergy = 63.50 *MeV;
+  // gaussian distribution with mean energy = 62.0 *MeV
+  // and sigma = 400.0 *keV
+  G4double defaultMeanKineticEnergy = 62.0 *MeV;
   meanKineticEnergy = defaultMeanKineticEnergy;
 
-  G4double defaultsigmaEnergy = 300.0 *keV;
+  G4double defaultsigmaEnergy = 400.0 *keV;
   sigmaEnergy = defaultsigmaEnergy;
+  
+#ifdef ANALYSIS_USE
+  // Write these values into the analysis if needed. Have to be written separately on change.
+  HadrontherapyAnalysisManager::getInstance()->setBeamMetaData(meanKineticEnergy, sigmaEnergy);
+#endif
 
   // Define the parameters of the initial position: 
   // the y, z coordinates have a gaussian distribution
-  G4double defaultX0 = -3248.59 *mm;  
+  
+  G4double defaultX0 = -2700.0 *mm;
   X0 = defaultX0;
 
@@ -111 +120 @@
 void HadrontherapyPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
 {
+#ifdef ANALYSIS_USE
+  // Increment the event counter
+  HadrontherapyAnalysisManager::getInstance()->startNewEvent();
+#endif
+
   // ****************************************
   // Set the beam angular apread 
@@ -162 +176 @@
 
 void HadrontherapyPrimaryGeneratorAction::SetmeanKineticEnergy (G4double val )  
-{ meanKineticEnergy = val;} 
+{
+        meanKineticEnergy = val;
+#ifdef ANALYSIS_USE
+  // Update the beam-data in the analysis manager
+  HadrontherapyAnalysisManager::getInstance()->setBeamMetaData(meanKineticEnergy, sigmaEnergy);
+#endif
+
+} 
 
 void HadrontherapyPrimaryGeneratorAction::SetsigmaEnergy (G4double val )  
-{ sigmaEnergy = val;}
+{ 
+        sigmaEnergy = val;
+#ifdef ANALYSIS_USE
+  // Update the sigmaenergy in the metadata.
+  HadrontherapyAnalysisManager::getInstance()->setBeamMetaData(meanKineticEnergy, sigmaEnergy);
+#endif
+}
 
 void HadrontherapyPrimaryGeneratorAction::SetXposition (G4double val )  
@@ -187 +214 @@
 void HadrontherapyPrimaryGeneratorAction::SetsigmaMomentumZ (G4double val )  
 { sigmaMomentumZ = val;}
+
+G4double HadrontherapyPrimaryGeneratorAction::GetmeanKineticEnergy(void)
+{ return meanKineticEnergy;}

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: HadrontherapyPrimaryGeneratorMessenger.cc; May 2005
-// ----------------------------------------------------------------------------
-//                 GEANT 4 - Hadrontherapy example
-// ----------------------------------------------------------------------------
-// Code developed by:
-//
-// G.A.P. Cirrone(a)*, 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
-// ----------------------------------------------------------------------------
+// HadrontherapyPrimaryGeneratorMessenger.cc;
+// See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
 
 #include "HadrontherapyPrimaryGeneratorMessenger.hh"

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: HadrontherapyRunAction.cc,v 3.0, September 2004;
-// ----------------------------------------------------------------------------
-//                 GEANT 4 - Hadrontherapy example
-// ----------------------------------------------------------------------------
-// Code developed by:
-//
-// G.A.P. Cirrone(a)*, F. Di Rosa(a), S. Guatelli(b), G. Russo(a)
-// 
-// (a) Laboratori Nazionali del Sud 
-//     of the INFN, Catania, Italy
-// (b) INFN Section of Genova, Genova, Italy
-// 
-// * cirrone@lns.infn.it
-// ----------------------------------------------------------------------------
+// $Id: HadrontherapyRunAction.cc
+// See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
 
 #include "HadrontherapyRunAction.hh"

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: HadrontherapyProtonSteppingAction.cc; May 2005
-// ----------------------------------------------------------------------------
-//                 GEANT 4 - Hadrontherapy example
-// ----------------------------------------------------------------------------
-// Code developed by:
-//
-// G.A.P. Cirrone(a)*, F. Di Rosa(a), S. Guatelli(b), G. Russo(a)
-// 
-// (a) Laboratori Nazionali del Sud 
-//     of the INFN, Catania, Italy
-// (b) INFN Section of Genova, Genova, Italy
-// 
-// * cirrone@lns.infn.it
-// ----------------------------------------------------------------------------
+// HadrontherapyProtonSteppingAction.cc; 
+// See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
 
 #include "G4SteppingManager.hh"
@@ -52 +40 @@
 #include "G4ParticleTypes.hh"
 
-#ifdef G4ANALYSIS_USE   
 #include "HadrontherapyAnalysisManager.hh"
-#endif
 
 #include "HadrontherapyRunAction.hh"
 
-HadrontherapySteppingAction::HadrontherapySteppingAction( HadrontherapyRunAction* run)
+/////////////////////////////////////////////////////////////////////////////
+HadrontherapySteppingAction::HadrontherapySteppingAction( HadrontherapyRunAction *run)
 {
   runAction = run;
 }
 
+/////////////////////////////////////////////////////////////////////////////
 HadrontherapySteppingAction::~HadrontherapySteppingAction()
 {
 }
 
+/////////////////////////////////////////////////////////////////////////////
 void HadrontherapySteppingAction::UserSteppingAction(const G4Step* aStep)
 { 
+    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();   
+        //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);
+        }
+        if(secondaryParticleName == "deuteron") {
+          analysis->hydrogenEnergy((secondaryParticleKineticEnergy/2) / MeV);
+        }
+        if(secondaryParticleName == "triton") {
+          analysis->hydrogenEnergy((secondaryParticleKineticEnergy/3) / MeV);
+        }
+        if(secondaryParticleName == "alpha") {
+          analysis->heliumEnergy((secondaryParticleKineticEnergy/4) / MeV);
+        }
+        if(secondaryParticleName == "He3"){
+          analysis->heliumEnergy((secondaryParticleKineticEnergy/3) / MeV);             
+        }
+#endif
+
+        aStep->GetTrack()->SetTrackStatus(fKillTrackAndSecondaries);
+    }
+
   // Electromagnetic and hadronic processes of primary particles in the phantom
-  
- if ((aStep -> GetTrack() -> GetTrackID() == 1) &&
+  //setting phantomPhys correctly will break something here fixme
+  if ((aStep -> GetTrack() -> GetTrackID() == 1) &&
     (aStep -> GetTrack() -> GetVolume() -> GetName() == "PhantomPhys") &&
     (aStep -> GetPostStepPoint() -> GetProcessDefinedStep() != NULL))
@@ -96 +134 @@
  // Retrieve information about the secondaries originated in the phantom
 
-#ifdef G4ANALYSIS_USE   
-  G4SteppingManager*  steppingManager = fpSteppingManager;
-  G4Track* theTrack = aStep -> GetTrack();
-
+ G4SteppingManager*  steppingManager = fpSteppingManager;
+  
   // check if it is alive
-  if(theTrack-> GetTrackStatus() == fAlive) { return; }
+  //if(theTrack-> GetTrackStatus() == fAlive) { return; }
 
   // Retrieve the secondary particles
@@ -110 +146 @@
       G4String volumeName = (*fSecondary)[lp1] -> GetVolume() -> GetName(); 
  
-      if (volumeName == "PhantomPhys")
+      if (volumeName == "phantomPhys")
         {
+#ifdef ANALYSIS_USE   
           G4String secondaryParticleName =  (*fSecondary)[lp1]->GetDefinition() -> GetParticleName();  
           G4double secondaryParticleKineticEnergy =  (*fSecondary)[lp1] -> GetKineticEnergy();     
-   
+
           HadrontherapyAnalysisManager* analysis =  HadrontherapyAnalysisManager::getInstance();   
         
@@ -137 +174 @@
               G4int a = (*fSecondary)[lp1]-> GetDynamicParticle() -> GetDefinition() -> GetBaryonNumber();
               G4int electronOccupancy = (*fSecondary)[lp1] ->  GetDynamicParticle() -> GetTotalOccupancy(); 
-              // If a generic ion is originated in the phantom, its baryonic number, PDG charge, 
+              
+              // 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);                   
+              analysis -> genericIonInformation(a, z, electronOccupancy, secondaryParticleKineticEnergy/MeV);
             }
+#endif
         }
     }
-#endif
 }