ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
HadrontherapyInteractionParameters.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file HadrontherapyInteractionParameters.cc
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 
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 
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 
65 {
66  if (pMessenger) delete pMessenger;
67  delete nistMat;
68  delete nistEle;
69 }
70 
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 }
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);
98  massDedx.push_back ( (dedxtot / density)/(MeV*cm2/g) );
99  }
100  }
101  else // one point only
102  {
103  energy.push_back(kinEmin/MeV);
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 == "") ? "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
131 {
133  if (Pmaterial) density = Pmaterial -> GetDensity();
134  return Pmaterial;
135 }
136 // Parse arguments line
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
192 {
193  // first check if RunManager is above G4State_Idle
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;
201  beamFlag = true;
202  }
203 
204 }
205 // print a list of Nist elements and materials
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