Changeset 1230 for trunk/examples/advanced/hadrontherapy/src

Timestamp:

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

Author:

garnier

Message:

update to geant4.9.3

Location:

trunk/examples/advanced/hadrontherapy/src

Files:

• HadrontherapyAnalysisManager.cc (modified) (7 diffs)
• HadrontherapyDetectorConstruction.cc (modified) (2 diffs)
• HadrontherapyDetectorMessenger.cc (modified) (2 diffs)
• HadrontherapyEventAction.cc (modified) (3 diffs)
• HadrontherapyMatrix.cc (modified) (5 diffs)
• HadrontherapyModulator.cc (modified) (6 diffs)
• HadrontherapyParticles.cc (modified) (1 diff)
• HadrontherapyPhysicsList.cc (modified) (1 diff)
• HadrontherapyPhysicsListMessenger.cc (modified) (1 diff)
• HadrontherapyPrimaryGeneratorAction.cc (modified) (7 diffs)
• HadrontherapyPrimaryGeneratorMessenger.cc (modified) (1 diff)
• HadrontherapyRunAction.cc (modified) (1 diff)
• HadrontherapySteppingAction.cc (modified) (5 diffs)

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
 }