hadrontherapy/src/HadrontherapyInteractionParameters.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // Hadrontherapy advanced example for Geant4
0027 // See more at: https://twiki.cern.ch/twiki/bin/view/Geant4/AdvancedExamplesHadrontherapy
0028
0029 #include <fstream>
0030 #include <iostream>
0031 #include <sstream>
0032 #include <cmath>
0033 #include <vector>
0034
0035 #include "HadrontherapyInteractionParameters.hh"
0036 #include "HadrontherapyParameterMessenger.hh"
0037 #include "HadrontherapyDetectorConstruction.hh"
0038
0039 #include "globals.hh"
0040 #include "G4SystemOfUnits.hh"
0041 #include "G4UnitsTable.hh"
0042 #include "G4UImanager.hh"
0043 #include "G4RunManager.hh"
0044 #include "G4LossTableManager.hh"
0045 #include "G4Material.hh"
0046 #include "G4MaterialCutsCouple.hh"
0047 #include "G4ParticleDefinition.hh"
0048 #include "G4ParticleTable.hh"
0049 #include "G4NistManager.hh"
0050 #include "G4Element.hh"
0051 #include "G4StateManager.hh"
0052
0053 HadrontherapyInteractionParameters::HadrontherapyInteractionParameters(G4bool wantMessenger):
0054     nistEle(new G4NistElementBuilder(0)),
0055     nistMat(new G4NistMaterialBuilder(nistEle, 0)),
0056     data(G4cout.rdbuf()),
0057     pMessenger(0),
0058     beamFlag(false)
0059
0060 {
0061     if (wantMessenger) pMessenger = new HadrontherapyParameterMessenger(this);
0062 }
0063
0064 HadrontherapyInteractionParameters::~HadrontherapyInteractionParameters()
0065 {
0066     if (pMessenger) delete pMessenger;
0067     delete nistMat;
0068     delete nistEle;
0069 }
0070
0071 G4double HadrontherapyInteractionParameters::GetStopping (G4double ene,
0072                                                       const G4ParticleDefinition* pDef,
0073                                   const G4Material* pMat,
0074                               G4double dens)
0075 {
0076     if (dens) return ComputeTotalDEDX(ene, pDef, pMat)/dens;
0077     return ComputeTotalDEDX(ene, pDef, pMat);
0078 }
0079 bool HadrontherapyInteractionParameters::GetStoppingTable(const G4String& vararg)
0080 {
0081     // Check arguments
0082     if ( !ParseArg(vararg)) return false;
0083     // Clear previous energy & mass sp vectors
0084         energy.clear();
0085         massDedx.clear();
0086     // log scale
0087     if (kinEmin != kinEmax && npoints >1)
0088     {
0089         G4double logmin = std::log10(kinEmin);
0090         G4double logmax = std::log10(kinEmax);
0091         G4double en;
0092         // uniform log space
0093             for (G4double c = 0.; c < npoints; c++)
0094            {
0095             en = std::pow(10., logmin + ( c*(logmax-logmin)  / (npoints - 1.)) );
0096             energy.push_back(en/MeV);
0097             dedxtot =  ComputeTotalDEDX (en, particle, material);
0098                 massDedx.push_back ( (dedxtot / density)/(MeV*cm2/g) );
0099            }
0100     }
0101     else // one point only
0102     {
0103         energy.push_back(kinEmin/MeV);
0104         dedxtot =  ComputeTotalDEDX (kinEmin, particle, material);
0105         massDedx.push_back ( (dedxtot / density)/(MeV*cm2/g) );
0106     }
0107
0108     G4cout.precision(6);
0109     data <<  "MeV             " << "MeV*cm2/g      " << particle << " (into " <<
0110         material << ", density = " << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
0111     data << G4endl;
0112     data << std::left << std::setfill(' ');
0113     for (size_t i=0; i<energy.size(); i++){
0114         data << std::setw(16) << energy[i] << massDedx[i] << G4endl;
0115     }
0116     outfile.close();
0117     // This will plot
0118
0119 // Info to user
0120     G4String ofName = (filename == "") ? G4String("User terminal"): filename;
0121     G4cout << "User choice:\n";
0122     G4cout << "Kinetic energy lower limit= "<< G4BestUnit(kinEmin,"Energy") <<
0123           ", Kinetic energy upper limit= " << G4BestUnit(kinEmax,"Energy") <<
0124              ", npoints= "<< npoints << ", particle= \"" << particle <<
0125          "\", material= \"" << material << "\", filename= \""<<
0126          ofName << "\"" << G4endl;
0127     return true;
0128 }
0129 // Search for user material choice inside G4NistManager database
0130 G4Material* HadrontherapyInteractionParameters::GetNistMaterial(G4String mat)
0131 {
0132     Pmaterial = G4NistManager::Instance()->FindOrBuildMaterial(mat);
0133     if (Pmaterial) density = Pmaterial -> GetDensity();
0134     return Pmaterial;
0135 }
0136 // Parse arguments line
0137 bool HadrontherapyInteractionParameters::ParseArg(const G4String& vararg)
0138 {
0139   kinEmin = kinEmax = npoints = 0.;
0140   particle = material = filename = "";
0141   // set internal variables
0142   std::istringstream strParam(vararg);
0143   // TODO here check for number and parameters consistency
0144   strParam >> std::skipws >> material >> kinEmin >> kinEmax >> npoints >> particle >> filename;
0145   // npoints must be an integer!
0146   npoints = std::floor(npoints);
0147
0148 // Check that  kinEmax >= kinEmin > 0 &&  npoints >= 1
0149 // TODO NIST points and linear scale
0150    if (kinEmax == 0. && kinEmin > 0. ) kinEmax = kinEmin;
0151    if (kinEmax == 0. && kinEmin == 0. ) kinEmax = kinEmin = 1.*MeV;
0152    if (kinEmax < kinEmin)
0153    {
0154      G4cout << "WARNING: kinEmin must not exceed kinEmax!" << G4endl;
0155      G4cout << "Usage: /parameter/command  material EkinMin EKinMax nPoints [particle] [output filename]" << G4endl;
0156      return false;
0157    }
0158   if (npoints < 1) npoints = 1;
0159
0160     // check if element/material is into database
0161     if (!GetNistMaterial(material) )
0162        {
0163          G4cout << "WARNING: material \"" << material << "\" doesn't exist in NIST elements/materials"
0164                " table [$G4INSTALL/source/materials/src/G4NistMaterialBuilder.cc]" << G4endl;
0165          G4cout << "Use command \"/parameter/nist\" to see full materials list" << G4endl;
0166          return false;
0167        }
0168     // Check for particle
0169     if (particle == "") particle = "proton"; // default to "proton"
0170     else if ( !FindParticle(particle) )
0171        {
0172         G4cout << "WARNING: Particle \"" << particle << "\" isn't supported." << G4endl;
0173         G4cout << "Try the command \"/particle/list\" to get full supported particles list." << G4endl;
0174         G4cout << "If you are interested in an ion that isn't in this list you must give it to the particle gun."
0175                   "\nTry the commands:\n/gun/particle ion"
0176               "\n/gun/ion <atomic number> <mass number> <[charge]>" << G4endl << G4endl;
0177         return false;
0178        }
0179     // start physics by forcing a G4RunManager::BeamOn():
0180     BeamOn();
0181     // Set output file
0182     if( filename != "" )
0183        {
0184           outfile.open(filename,std::ios_base::trunc); // overwrite existing file
0185           data.rdbuf(outfile.rdbuf());
0186        }
0187     else data.rdbuf(G4cout.rdbuf());    // output is G4cout!
0188     return true;
0189 }
0190 // Force physics tables build
0191 void HadrontherapyInteractionParameters::BeamOn()
0192 {
0193     // first check if RunManager is above G4State_Idle
0194     G4StateManager* mState = G4StateManager::GetStateManager();
0195     G4ApplicationState  aState = mState -> GetCurrentState();
0196     if ( aState <= G4State_Idle && beamFlag == false)
0197      {
0198         G4cout << "Issuing a G4RunManager::beamOn()... ";
0199         G4cout << "Current Run State is " << mState -> GetStateString( aState ) << G4endl;
0200         G4RunManager::GetRunManager() -> BeamOn(0);
0201         beamFlag = true;
0202      }
0203
0204 }
0205 // print a list of Nist elements and materials
0206 void HadrontherapyInteractionParameters::ListOfNistMaterials(const G4String& vararg)
0207 {
0208 /*
0209  $G4INSTALL/source/materials/src/G4NistElementBuilder.cc
0210  You can also construct a new material by the ConstructNewMaterial method:
0211  see $G4INSTALL/source/materials/src/G4NistMaterialBuilder.cc
0212 */
0213     // Get simplest full list
0214      if (vararg =="list")
0215       {
0216         const std::vector<G4String>& vec =  nistMat -> GetMaterialNames();
0217         for (size_t i=0; i<vec.size(); i++)
0218          {
0219             G4cout << std::setw(12) << std::left << i+1 << vec[i] << G4endl;
0220          }
0221         G4cout << G4endl;
0222           }
0223     else if (vararg =="all" || vararg =="simple" || vararg =="compound" || vararg =="hep" )
0224     {
0225         nistMat -> ListMaterials(vararg);
0226     }
0227 }
0228
0229