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