range/src/Run.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 // This example is provided by the Geant4-DNA collaboration
0027 // Any report or published results obtained using the Geant4-DNA software
0028 // shall cite the following Geant4-DNA collaboration publications:
0029 // Med. Phys. 45 (2018) e722-e739
0030 // Phys. Med. 31 (2015) 861-874
0031 // Med. Phys. 37 (2010) 4692-4708
0032 // Int. J. Model. Simul. Sci. Comput. 1 (2010) 157–178
0033 //
0034 // The Geant4-DNA web site is available at http://geant4-dna.org
0035 //
0036 /// \file medical/dna/range/src/Run.cc
0037 /// \brief Implementation of the Run class
0038
0039 #include "Run.hh"
0040
0041 #include "PrimaryGeneratorAction.hh"
0042
0043 #include "G4Material.hh"
0044 #include "G4SystemOfUnits.hh"
0045 #include "G4UnitsTable.hh"
0046
0047 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0048
0049 Run::Run(const DetectorConstruction* detector)
0050   : G4Run(),
0051     fDetector(detector),
0052     fParticle(0),
0053     fEkin(0.),
0054     fEdeposit(0.),
0055     fEdeposit2(0.),
0056     fTrackLen(0.),
0057     fTrackLen2(0.),
0058     fProjRange(0.),
0059     fProjRange2(0.),
0060     fPenetration(0.),
0061     fPenetration2(0.),
0062     fNbOfSteps(0),
0063     fNbOfSteps2(0),
0064     fStepSize(0.),
0065     fStepSize2(0.)
0066 {}
0067
0068 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0069
0070 Run::~Run() {}
0071
0072 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0073
0074 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
0075 {
0076   fParticle = particle;
0077   fEkin = energy;
0078 }
0079
0080 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0081
0082 void Run::AddEdep(G4double e)
0083 {
0084   fEdeposit += e;
0085   fEdeposit2 += e * e;
0086 }
0087
0088 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0089
0090 void Run::AddTrackLength(G4double t)
0091 {
0092   fTrackLen += t;
0093   fTrackLen2 += t * t;
0094 }
0095
0096 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0097
0098 void Run::AddProjRange(G4double x)
0099 {
0100   fProjRange += x;
0101   fProjRange2 += x * x;
0102 }
0103
0104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0105
0106 void Run::AddPenetration(G4double x)
0107 {
0108   fPenetration += x;
0109   fPenetration2 += x * x;
0110 }
0111
0112 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0113
0114 void Run::AddStepSize(G4int nb, G4double st)
0115 {
0116   fNbOfSteps += nb;
0117   fNbOfSteps2 += nb * nb;
0118   fStepSize += st;
0119   fStepSize2 += st * st;
0120 }
0121
0122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0123
0124 void Run::Merge(const G4Run* run)
0125 {
0126   const Run* localRun = static_cast<const Run*>(run);
0127
0128   // Pass information about primary particle
0129   fParticle = localRun->fParticle;
0130   fEkin = localRun->fEkin;
0131
0132   // Accumulate sums
0133   fEdeposit += localRun->fEdeposit;
0134   fEdeposit2 += localRun->fEdeposit2;
0135   fTrackLen += localRun->fTrackLen;
0136   fTrackLen2 += localRun->fTrackLen2;
0137   fProjRange += localRun->fProjRange;
0138   fProjRange2 += localRun->fProjRange2;
0139   fPenetration += localRun->fPenetration;
0140   fPenetration2 += localRun->fPenetration2;
0141   fNbOfSteps += localRun->fNbOfSteps;
0142   fNbOfSteps2 += localRun->fNbOfSteps2;
0143   fStepSize += localRun->fStepSize;
0144   fStepSize2 += localRun->fStepSize2;
0145
0146   G4Run::Merge(run);
0147 }
0148
0149 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0150
0151 void Run::EndOfRun()
0152 {
0153   std::ios::fmtflags mode = G4cout.flags();
0154   G4cout.setf(std::ios::fixed, std::ios::floatfield);
0155   G4int prec = G4cout.precision(2);
0156
0157   // Run conditions
0158   G4Material* material = fDetector->GetAbsorMaterial();
0159   G4double density = material->GetDensity();
0160   G4String partName = fParticle->GetParticleName();
0161
0162   G4cout << "\n ======================= run summary ====================\n";
0163   G4cout << "\n The run is " << numberOfEvent << " " << partName << " of "
0164          << G4BestUnit(fEkin, "Energy") << " through a sphere of radius "
0165          << G4BestUnit(fDetector->GetAbsorRadius(), "Length") << "of " << material->GetName()
0166          << " (density: " << G4BestUnit(density, "Volumic Mass") << ")" << G4endl;
0167
0168   if (numberOfEvent == 0) {
0169     G4cout.setf(mode, std::ios::floatfield);
0170     G4cout.precision(prec);
0171     return;
0172   }
0173
0174   // Compute track length of primary track
0175   fTrackLen /= numberOfEvent;
0176   fTrackLen2 /= numberOfEvent;
0177   G4double rmsTrack = fTrackLen2 - fTrackLen * fTrackLen;
0178
0179   if (rmsTrack > 0.)
0180     rmsTrack = std::sqrt(rmsTrack);
0181   else
0182     rmsTrack = 0.;
0183
0184   G4cout.precision(3);
0185   G4cout << "\n Track length of primary track = " << G4BestUnit(fTrackLen, "Length") << " +- "
0186          << G4BestUnit(rmsTrack, "Length");
0187
0188   // Compute projected range of primary track
0189   fProjRange /= numberOfEvent;
0190   fProjRange2 /= numberOfEvent;
0191   G4double rmsProj = fProjRange2 - fProjRange * fProjRange;
0192   if (rmsProj > 0.)
0193     rmsProj = std::sqrt(rmsProj);
0194   else
0195     rmsProj = 0.;
0196
0197   G4cout << "\n Projected range               = " << G4BestUnit(fProjRange, "Length") << " +- "
0198          << G4BestUnit(rmsProj, "Length");
0199
0200   // Compute penetration of primary track
0201   fPenetration /= numberOfEvent;
0202   fPenetration2 /= numberOfEvent;
0203   G4double rmsPene = fPenetration2 - fPenetration * fPenetration;
0204   if (rmsPene > 0.)
0205     rmsPene = std::sqrt(rmsPene);
0206   else
0207     rmsPene = 0.;
0208
0209   G4cout << "\n Penetration                   = " << G4BestUnit(fPenetration, "Length") << " +- "
0210          << G4BestUnit(rmsPene, "Length") << G4endl;
0211
0212   //
0213
0214   // Output file
0215   FILE* myFile;
0216   myFile = fopen("range.txt", "a");
0217   fprintf(myFile, "%e %e %e %e %e %e %e\n", fEkin / eV, fTrackLen / nm, rmsTrack / nm,
0218           fProjRange / nm, rmsProj / nm, fPenetration / nm, rmsPene / nm);
0219   fclose(myFile);
0220
0221   // Reset default formats
0222   G4cout.setf(mode, std::ios::floatfield);
0223   G4cout.precision(prec);
0224 }