mfp/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 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     fTotalCount(0),
0055     fSumTrack(0.),
0056     fSumTrack2(0.),
0057     fEnTransfer(0.)
0058 {}
0059
0060 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0061
0062 Run::~Run() {}
0063
0064 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0065
0066 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
0067 {
0068   fParticle = particle;
0069   fEkin = energy;
0070 }
0071
0072 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0073
0074 void Run::CountProcesses(G4String procName)
0075 {
0076   std::map<G4String, G4int>::iterator it = fProcCounter.find(procName);
0077   if (it == fProcCounter.end()) {
0078     fProcCounter[procName] = 1;
0079   }
0080   else {
0081     fProcCounter[procName]++;
0082   }
0083 }
0084
0085 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0086
0087 void Run::SumTrack(G4double track)
0088 {
0089   fTotalCount++;
0090   fSumTrack += track;
0091   fSumTrack2 += track * track;
0092 }
0093
0094 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0095
0096 void Run::SumeTransf(G4double energy)
0097 {
0098   fEnTransfer += energy;
0099 }
0100
0101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0102
0103 void Run::Merge(const G4Run* run)
0104 {
0105   const Run* localRun = static_cast<const Run*>(run);
0106
0107   // Pass information about primary particle
0108   fParticle = localRun->fParticle;
0109   fEkin = localRun->fEkin;
0110
0111   // map: processes count
0112   std::map<G4String, G4int>::const_iterator it;
0113   for (it = localRun->fProcCounter.begin(); it != localRun->fProcCounter.end(); ++it) {
0114     G4String procName = it->first;
0115     G4int localCount = it->second;
0116
0117     if (fProcCounter.find(procName) == fProcCounter.end()) {
0118       fProcCounter[procName] = localCount;
0119     }
0120     else {
0121       fProcCounter[procName] += localCount;
0122     }
0123   }
0124
0125   fTotalCount += localRun->fTotalCount;
0126   fSumTrack += localRun->fSumTrack;
0127   fSumTrack2 += localRun->fSumTrack2;
0128   fEnTransfer += localRun->fEnTransfer;
0129
0130   G4Run::Merge(run);
0131 }
0132
0133 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0134
0135 void Run::EndOfRun()
0136 {
0137   std::ios::fmtflags mode = G4cout.flags();
0138   G4cout.setf(std::ios::fixed, std::ios::floatfield);
0139   G4int prec = G4cout.precision(2);
0140
0141   // Run conditions
0142   G4Material* material = fDetector->GetAbsorMaterial();
0143   G4double density = material->GetDensity();
0144   G4String partName = fParticle->GetParticleName();
0145
0146   G4cout << "\n ======================== run summary =====================\n";
0147   G4cout << "\n The run is " << numberOfEvent << " " << partName << " of "
0148          << G4BestUnit(fEkin, "Energy") << " through a sphere of radius "
0149          << G4BestUnit(fDetector->GetAbsorRadius(), "Length") << "of " << material->GetName()
0150          << " (density: " << G4BestUnit(density, "Volumic Mass") << ")" << G4endl;
0151
0152   if (numberOfEvent == 0) {
0153     G4cout.setf(mode, std::ios::floatfield);
0154     G4cout.precision(prec);
0155     return;
0156   }
0157
0158   // Frequency of processes
0159   G4int survive = 0;
0160   G4cout << "\n Process calls frequency --->";
0161   std::map<G4String, G4int>::iterator it;
0162   for (it = fProcCounter.begin(); it != fProcCounter.end(); it++) {
0163     G4String procName = it->first;
0164     G4int count = it->second;
0165     G4cout << "\t" << procName << " = " << count;
0166     if (procName == "Transportation") survive = count;
0167   }
0168
0169   if (survive > 0) {
0170     G4cout << "\n\n Nb of incident particles surviving after "
0171            << "a radius of " << G4BestUnit(fDetector->GetAbsorRadius(), "Length") << " of "
0172            << material->GetName() << " : " << survive << G4endl;
0173   }
0174
0175   if (fTotalCount == 0) fTotalCount = 1;  // force printing anyway
0176
0177   // Compute mean free path and related quantities
0178   G4double MeanFreePath = fSumTrack / fTotalCount;
0179   G4double MeanTrack2 = fSumTrack2 / fTotalCount;
0180   G4double rmsBis = std::sqrt(std::fabs(MeanTrack2 - MeanFreePath * MeanFreePath));
0181   G4double CrossSection = 1. / MeanFreePath;
0182   G4double massicMFP = MeanFreePath * density;
0183   G4double massicCS = 1. / massicMFP;
0184
0185   G4cout << "\n\n MeanFreePath:\t" << G4BestUnit(MeanFreePath, "Length") << " +- "
0186          << G4BestUnit(rmsBis, "Length")
0187          << "\t\t\tmassic: " << G4BestUnit(massicMFP, "Mass/Surface") << "\n CrossSection:\t"
0188          << CrossSection * cm << " cm^-1 "
0189          << "\t\t\tmassic: " << G4BestUnit(massicCS, "Surface/Mass") << G4endl;
0190
0191   // Compute energy transfer coefficient
0192   G4double MeanTransfer = fEnTransfer / fTotalCount;
0193   G4double massTransfCoef = massicCS * MeanTransfer / fEkin;
0194
0195   G4cout << "\n mean energy of charged secondaries: " << G4BestUnit(MeanTransfer, "Energy")
0196          << "\tmass_energy_transfer coef: " << G4BestUnit(massTransfCoef, "Surface/Mass") << G4endl;
0197
0198   // Output file
0199   FILE* myFile;
0200   myFile = fopen("mfp.txt", "a");
0201   fprintf(myFile, "%e %e %e \n", fEkin / eV, MeanFreePath / nm, rmsBis / nm);
0202   fclose(myFile);
0203
0204   // Remove all contents in fProcCounter
0205   fProcCounter.clear();
0206
0207   // Reset default formats
0208   G4cout.setf(mode, std::ios::floatfield);
0209   G4cout.precision(prec);
0210 }