Geant4  10.01.p02
HadrontherapyInteractionParameters.cc
Go to the documentation of this file.
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 // Visit the Hadrontherapy web site (http://www.lns.infn.it/link/Hadrontherapy) to request
26 // the *COMPLETE* version of this program, together with its documentation;
27 // Hadrontherapy (both basic and full version) are supported by the Italian INFN
28 // Institute in the framework of the MC-INFN Group
29 //
30 
31 #include <fstream>
32 #include <iostream>
33 #include <sstream>
34 #include <cmath>
35 #include <vector>
36 
40 
41 #include "globals.hh"
42 #include "G4SystemOfUnits.hh"
43 #include "G4UnitsTable.hh"
44 #include "G4UImanager.hh"
45 #include "G4RunManager.hh"
46 #include "G4LossTableManager.hh"
47 #include "G4Material.hh"
48 #include "G4MaterialCutsCouple.hh"
49 #include "G4ParticleDefinition.hh"
50 #include "G4ParticleTable.hh"
51 #include "G4NistManager.hh"
52 #include "G4Element.hh"
53 #include "G4StateManager.hh"
54 
56  nistEle(new G4NistElementBuilder(0)),
57  nistMat(new G4NistMaterialBuilder(nistEle, 0)),
58  data(G4cout.rdbuf()),
59  pMessenger(0),
60  beamFlag(false)
61 #ifdef G4ANALYSIS_USE_ROOT
62  ,theRootCanvas(0),
63  theRootGraph(0)
64 #endif
65 {
66  if (wantMessenger) pMessenger = new HadrontherapyParameterMessenger(this);
67 }
68 
70 {
71  if (pMessenger) delete pMessenger;
72  delete nistMat;
73  delete nistEle;
74 }
75 
77  const G4ParticleDefinition* pDef,
78  const G4Material* pMat,
79  G4double dens)
80 {
81  if (dens) return ComputeTotalDEDX(ene, pDef, pMat)/dens;
82  return ComputeTotalDEDX(ene, pDef, pMat);
83 }
85 {
86  // Check arguments
87  if ( !ParseArg(vararg)) return false;
88  // Clear previous energy & mass sp vectors
89  energy.clear();
90  massDedx.clear();
91  // log scale
92  if (kinEmin != kinEmax && npoints >1)
93  {
94  G4double logmin = std::log10(kinEmin);
95  G4double logmax = std::log10(kinEmax);
96  G4double en;
97  // uniform log space
98  for (G4double c = 0.; c < npoints; c++)
99  {
100  en = std::pow(10., logmin + ( c*(logmax-logmin) / (npoints - 1.)) );
101  energy.push_back(en/MeV);
103  massDedx.push_back ( (dedxtot / density)/(MeV*cm2/g) );
104  }
105  }
106  else // one point only
107  {
108  energy.push_back(kinEmin/MeV);
110  massDedx.push_back ( (dedxtot / density)/(MeV*cm2/g) );
111  }
112 
113  G4cout.precision(6);
114  data << "MeV " << "MeV*cm2/g " << particle << " (into " <<
115  material << ", density = " << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
116  data << G4endl;
117  data << std::left << std::setfill(' ');
118  for (size_t i=0; i<energy.size(); i++){
119  data << std::setw(16) << energy[i] << massDedx[i] << G4endl;
120  }
121  outfile.close();
122  // This will plot
123 #ifdef G4ANALYSIS_USE_ROOT
124  PlotStopping("pdf");
125 #endif
126 
127 // Info to user
128  G4String ofName = (filename == "") ? "User terminal": filename;
129  G4cout << "User choice:\n";
130  G4cout << "Kinetic energy lower limit= "<< G4BestUnit(kinEmin,"Energy") <<
131  ", Kinetic energy upper limit= " << G4BestUnit(kinEmax,"Energy") <<
132  ", npoints= "<< npoints << ", particle= \"" << particle <<
133  "\", material= \"" << material << "\", filename= \""<<
134  ofName << "\"" << G4endl;
135  return true;
136 }
138 // Save Plot
139 #ifdef G4ANALYSIS_USE_ROOT
140 void HadrontherapyInteractionParameters::PlotStopping(const G4String& filetype)
141 {
142  if (!theRootCanvas)
143  {
144  gROOT->Reset();
145  gROOT->SetStyle("Plain");
146  theRootCanvas = new TCanvas("theRootCanvas","Interaction Parameters",200, 10, 600,400);
147  theRootCanvas -> SetFillColor(20);
148  theRootCanvas -> SetBorderMode(1);
149  theRootCanvas -> SetBorderSize(1);
150  theRootCanvas -> SetFrameBorderMode(0);
151  theRootCanvas -> SetGrid();
152  // Use global pad: root manual pgg 109,...
153  }
154 
155  if (theRootGraph) delete theRootGraph;
156  theRootGraph = new TGraph(energy.size(), &energy[0], &massDedx[0]);
157  //theRootGraph = new TGraph();
158  axisX = theRootGraph -> GetXaxis(),
159  axisY = theRootGraph -> GetYaxis();
160  axisX -> SetTitle("MeV");
161  axisY -> SetTitle("Stopping Power (MeV cm2/g)");
162  //axisX -> SetNdivisions(500,kTRUE);
163  //axisX -> SetTickLength(0.03);
164  //axisX -> SetLabelOffset(2.005);
165  axisX -> SetAxisColor(2);
166  axisY -> SetAxisColor(2);
167  gPad -> SetLogx(1);
168  gPad -> SetLogy(1);
169  theRootGraph -> SetMarkerColor(4);
170  theRootGraph -> SetMarkerStyle(20);// circle
171  theRootGraph -> SetMarkerSize(.5);
172 
173  G4String gName = particle.substr(0, particle.find("[") ); // cut excitation energy
174  gName = gName + "_" + material;
175  G4String fName = "./referenceData/interaction/" + gName + "." + filetype;
176  theRootGraph -> SetTitle(gName);
177  theRootGraph -> Draw("AP");
178  //theRootCanvas -> Update();
179  //theRootCanvas -> Draw();
180  theRootCanvas -> SaveAs(fName);
181 }
182 #endif
183 
184 // Search for user material choice inside G4NistManager database
186 {
188  if (Pmaterial) density = Pmaterial -> GetDensity();
189  return Pmaterial;
190 }
191 // Parse arguments line
193 {
194  kinEmin = kinEmax = npoints = 0.;
195  particle = material = filename = "";
196  // set internal variables
197  std::istringstream strParam(vararg);
198  // TODO here check for number and parameters consistency
199  strParam >> std::skipws >> material >> kinEmin >> kinEmax >> npoints >> particle >> filename;
200  // npoints must be an integer!
201  npoints = std::floor(npoints);
202 
203 // Check that kinEmax >= kinEmin > 0 && npoints >= 1
204 // TODO NIST points and linear scale
205  if (kinEmax == 0. && kinEmin > 0. ) kinEmax = kinEmin;
206  if (kinEmax == 0. && kinEmin == 0. ) kinEmax = kinEmin = 1.*MeV;
207  if (kinEmax < kinEmin)
208  {
209  G4cout << "WARNING: kinEmin must not exceed kinEmax!" << G4endl;
210  G4cout << "Usage: /parameter/command material kinetic Emin kinetic Emax nPoints [particle] [output filename]" << G4endl;
211  return false;
212  }
213  if (npoints < 1) npoints = 1;
214 
215  // check if element/material is into database
216  if (!GetNistMaterial(material) )
217  {
218  G4cout << "WARNING: material \"" << material << "\" doesn't exist in NIST elements/materials"
219  " table [$G4INSTALL/source/materials/src/G4NistMaterialBuilder.cc]" << G4endl;
220  G4cout << "Use command \"/parameter/nist\" to see full materials list" << G4endl;
221  return false;
222  }
223  // Check for particle
224  if (particle == "") particle = "proton"; // default to "proton"
225  else if ( !FindParticle(particle) )
226  {
227  G4cout << "WARNING: Particle \"" << particle << "\" isn't supported." << G4endl;
228  G4cout << "Try the command \"/particle/list\" to get full supported particles list." << G4endl;
229  G4cout << "If you are interested in an ion that isn't in this list you must give it to the particle gun."
230  "\nTry the commands:\n/gun/particle ion"
231  "\n/gun/ion <atomic number> <mass number> <[charge]>" << G4endl << G4endl;
232  return false;
233  }
234  // start physics by forcing a G4RunManager::BeamOn():
235  BeamOn();
236  // Set output file
237  if( filename != "" )
238  {
239  outfile.open(filename,std::ios_base::trunc); // overwrite existing file
240  data.rdbuf(outfile.rdbuf());
241  }
242  else data.rdbuf(G4cout.rdbuf()); // output is G4cout!
243  return true;
244 }
245 // Force physics tables build
247 {
248  // first check if RunManager is above G4State_Idle
250  G4ApplicationState aState = mState -> GetCurrentState();
251  if ( aState <= G4State_Idle && beamFlag == false)
252  {
253  G4cout << "Issuing a G4RunManager::beamOn()... ";
254  G4cout << "Current Run State is " << mState -> GetStateString( aState ) << G4endl;
256  beamFlag = true;
257  }
258 
259 }
260 // print a list of Nist elements and materials
262 {
263 /*
264  $G4INSTALL/source/materials/src/G4NistElementBuilder.cc
265  You can also construct a new material by the ConstructNewMaterial method:
266  see $G4INSTALL/source/materials/src/G4NistMaterialBuilder.cc
267 */
268  // Get simplest full list
269  if (vararg =="list")
270  {
271  const std::vector<G4String>& vec = nistMat -> GetMaterialNames();
272  for (size_t i=0; i<vec.size(); i++)
273  {
274  G4cout << std::setw(12) << std::left << i+1 << vec[i] << G4endl;
275  }
276  G4cout << G4endl;
277  }
278  else if (vararg =="all" || vararg =="simple" || vararg =="compound" || vararg =="hep" )
279  {
280  nistMat -> ListMaterials(vararg);
281  }
282 }
283 
284 
G4Material * FindOrBuildMaterial(const G4String &name, G4bool isotopes=true, G4bool warning=false)
static const double MeV
Definition: G4SIunits.hh:193
static const double cm2
Definition: G4SIunits.hh:107
G4String fName
Definition: G4AttUtils.hh:55
G4double ComputeTotalDEDX(G4double kinEnergy, const G4ParticleDefinition *, const G4Material *, G4double cut=DBL_MAX)
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
static G4NistManager * Instance()
G4double GetStopping(G4double energy, const G4ParticleDefinition *, const G4Material *, G4double density=0.)
static G4StateManager * GetStateManager()
G4GLOB_DLL std::ostream G4cout
bool G4bool
Definition: G4Types.hh:79
static G4RunManager * GetRunManager()
Definition: G4RunManager.cc:79
const G4ParticleDefinition * FindParticle(const G4String &)
static const double g
Definition: G4SIunits.hh:162
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4ApplicationState