Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/advanced/hadrontherapy/src/HadrontherapyInteractionParameters.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  1 //
  2 // ********************************************************************
  3 // * License and Disclaimer                                           *
  4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.                             *
 10 // *                                                                  *
 11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                                                  *
 18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // ********************************************************************
 25 //
 26 // Hadrontherapy advanced example for Geant4
 27 // See more at: https://twiki.cern.ch/twiki/bin/view/Geant4/AdvancedExamplesHadrontherapy
 28 
 29 #include <fstream>
 30 #include <iostream>
 31 #include <sstream>
 32 #include <cmath>
 33 #include <vector>
 34 
 35 #include "HadrontherapyInteractionParameters.hh"
 36 #include "HadrontherapyParameterMessenger.hh"
 37 #include "HadrontherapyDetectorConstruction.hh"
 38 
 39 #include "globals.hh"
 40 #include "G4SystemOfUnits.hh"
 41 #include "G4UnitsTable.hh"
 42 #include "G4UImanager.hh"
 43 #include "G4RunManager.hh"
 44 #include "G4LossTableManager.hh"
 45 #include "G4Material.hh"
 46 #include "G4MaterialCutsCouple.hh"
 47 #include "G4ParticleDefinition.hh"
 48 #include "G4ParticleTable.hh"
 49 #include "G4NistManager.hh"
 50 #include "G4Element.hh"
 51 #include "G4StateManager.hh"
 52 
 53 HadrontherapyInteractionParameters::HadrontherapyInteractionParameters(G4bool wantMessenger): 
 54     nistEle(new G4NistElementBuilder(0)),                     
 55     nistMat(new G4NistMaterialBuilder(nistEle, 0)),                   
 56     data(G4cout.rdbuf()), 
 57     pMessenger(0),
 58     beamFlag(false)
 59 
 60 {
 61     if (wantMessenger) pMessenger = new HadrontherapyParameterMessenger(this); 
 62 }
 63 
 64 HadrontherapyInteractionParameters::~HadrontherapyInteractionParameters()
 65 {
 66     if (pMessenger) delete pMessenger; 
 67     delete nistMat; 
 68     delete nistEle; 
 69 }
 70 
 71 G4double HadrontherapyInteractionParameters::GetStopping (G4double ene, 
 72                                                     const G4ParticleDefinition* pDef, 
 73                       const G4Material* pMat,
 74                 G4double dens)
 75 {
 76     if (dens) return ComputeTotalDEDX(ene, pDef, pMat)/dens;
 77     return ComputeTotalDEDX(ene, pDef, pMat);
 78 }
 79 bool HadrontherapyInteractionParameters::GetStoppingTable(const G4String& vararg)
 80 {
 81   // Check arguments
 82   if ( !ParseArg(vararg)) return false;
 83   // Clear previous energy & mass sp vectors
 84         energy.clear(); 
 85         massDedx.clear();
 86   // log scale 
 87   if (kinEmin != kinEmax && npoints >1)
 88   {
 89       G4double logmin = std::log10(kinEmin);
 90       G4double logmax = std::log10(kinEmax); 
 91       G4double en;
 92       // uniform log space
 93             for (G4double c = 0.; c < npoints; c++)
 94          {
 95         en = std::pow(10., logmin + ( c*(logmax-logmin)  / (npoints - 1.)) );  
 96         energy.push_back(en/MeV);
 97         dedxtot =  ComputeTotalDEDX (en, particle, material);
 98               massDedx.push_back ( (dedxtot / density)/(MeV*cm2/g) );
 99          }
100   }
101   else // one point only
102   {
103       energy.push_back(kinEmin/MeV);
104       dedxtot =  ComputeTotalDEDX (kinEmin, particle, material);
105       massDedx.push_back ( (dedxtot / density)/(MeV*cm2/g) );
106   }
107 
108     G4cout.precision(6);  
109     data <<  "MeV             " << "MeV*cm2/g      " << particle << " (into " << 
110       material << ", density = " << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
111     data << G4endl;
112     data << std::left << std::setfill(' ');
113     for (size_t i=0; i<energy.size(); i++){
114     data << std::setw(16) << energy[i] << massDedx[i] << G4endl;
115   }
116     outfile.close();
117     // This will plot
118 
119 // Info to user
120     G4String ofName = (filename == "") ? G4String("User terminal"): filename;
121     G4cout << "User choice:\n";
122     G4cout << "Kinetic energy lower limit= "<< G4BestUnit(kinEmin,"Energy") << 
123         ", Kinetic energy upper limit= " << G4BestUnit(kinEmax,"Energy") << 
124            ", npoints= "<< npoints << ", particle= \"" << particle << 
125      "\", material= \"" << material << "\", filename= \""<< 
126      ofName << "\"" << G4endl;
127     return true;
128 }
129 // Search for user material choice inside G4NistManager database
130 G4Material* HadrontherapyInteractionParameters::GetNistMaterial(G4String mat)
131 {
132     Pmaterial = G4NistManager::Instance()->FindOrBuildMaterial(mat);
133     if (Pmaterial) density = Pmaterial -> GetDensity(); 
134     return Pmaterial;
135 }
136 // Parse arguments line
137 bool HadrontherapyInteractionParameters::ParseArg(const G4String& vararg)
138 {
139   kinEmin = kinEmax = npoints = 0.;
140   particle = material = filename = "";
141   // set internal variables
142   std::istringstream strParam(vararg);
143   // TODO here check for number and parameters consistency 
144   strParam >> std::skipws >> material >> kinEmin >> kinEmax >> npoints >> particle >> filename;
145   // npoints must be an integer!
146   npoints = std::floor(npoints); 
147 
148 // Check that  kinEmax >= kinEmin > 0 &&  npoints >= 1 
149 // TODO NIST points and linear scale
150    if (kinEmax == 0. && kinEmin > 0. ) kinEmax = kinEmin;
151    if (kinEmax == 0. && kinEmin == 0. ) kinEmax = kinEmin = 1.*MeV;
152    if (kinEmax < kinEmin) 
153    {
154      G4cout << "WARNING: kinEmin must not exceed kinEmax!" << G4endl;
155      G4cout << "Usage: /parameter/command  material EkinMin EKinMax nPoints [particle] [output filename]" << G4endl;  
156      return false;
157    }
158   if (npoints < 1) npoints = 1;
159 
160     // check if element/material is into database
161     if (!GetNistMaterial(material) )
162      {
163        G4cout << "WARNING: material \"" << material << "\" doesn't exist in NIST elements/materials"
164            " table [$G4INSTALL/source/materials/src/G4NistMaterialBuilder.cc]" << G4endl; 
165        G4cout << "Use command \"/parameter/nist\" to see full materials list" << G4endl; 
166        return false;
167      }
168     // Check for particle
169     if (particle == "") particle = "proton"; // default to "proton"
170     else if ( !FindParticle(particle) )
171      {
172     G4cout << "WARNING: Particle \"" << particle << "\" isn't supported." << G4endl;
173     G4cout << "Try the command \"/particle/list\" to get full supported particles list." << G4endl;
174     G4cout << "If you are interested in an ion that isn't in this list you must give it to the particle gun."
175               "\nTry the commands:\n/gun/particle ion"
176         "\n/gun/ion <atomic number> <mass number> <[charge]>" << G4endl << G4endl;
177     return false;
178      }
179     // start physics by forcing a G4RunManager::BeamOn(): 
180     BeamOn();
181     // Set output file
182     if( filename != "" ) 
183        {
184           outfile.open(filename,std::ios_base::trunc); // overwrite existing file
185           data.rdbuf(outfile.rdbuf());
186        }
187     else data.rdbuf(G4cout.rdbuf());  // output is G4cout!                
188     return true;
189 }
190 // Force physics tables build
191 void HadrontherapyInteractionParameters::BeamOn()
192 {
193     // first check if RunManager is above G4State_Idle 
194     G4StateManager* mState = G4StateManager::GetStateManager();
195     G4ApplicationState  aState = mState -> GetCurrentState(); 
196     if ( aState <= G4State_Idle && beamFlag == false)
197    {
198       G4cout << "Issuing a G4RunManager::beamOn()... "; 
199       G4cout << "Current Run State is " << mState -> GetStateString( aState ) << G4endl; 
200       G4RunManager::GetRunManager() -> BeamOn(0);
201       beamFlag = true;
202    }
203 
204 }
205 // print a list of Nist elements and materials
206 void HadrontherapyInteractionParameters::ListOfNistMaterials(const G4String& vararg)
207 {
208 /*
209  $G4INSTALL/source/materials/src/G4NistElementBuilder.cc
210  You can also construct a new material by the ConstructNewMaterial method:
211  see $G4INSTALL/source/materials/src/G4NistMaterialBuilder.cc
212 */
213     // Get simplest full list
214      if (vararg =="list")
215     {
216       const std::vector<G4String>& vec =  nistMat -> GetMaterialNames(); 
217       for (size_t i=0; i<vec.size(); i++)
218        {
219           G4cout << std::setw(12) << std::left << i+1 << vec[i] << G4endl;
220        }
221       G4cout << G4endl;
222           }
223     else if (vararg =="all" || vararg =="simple" || vararg =="compound" || vararg =="hep" )
224   {
225     nistMat -> ListMaterials(vararg);
226   }
227 }
228 
229 
230