wvalue/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/wvalue/src/Run.cc
0037 /// \brief Implementation of the Run class
0038
0039 #include "Run.hh"
0040
0041 #include "DetectorConstruction.hh"
0042 #include "HistoManager.hh"
0043 #include "PrimaryGeneratorAction.hh"
0044
0045 #include "G4Material.hh"
0046 #include "G4SystemOfUnits.hh"
0047 #include "G4UnitsTable.hh"
0048
0049 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0050
0051 Run::Run(const DetectorConstruction* detector)
0052   : G4Run(),
0053     fDetector(detector),
0054     fParticle(0),
0055     fEkin(0.),
0056     fNbInelastic(0),
0057     fNbInelastic2(0),
0058     fEdeposit(0.),
0059     fEdeposit2(0.),
0060     fTrackLen(0.),
0061     fTrackLen2(0.),
0062     fProjRange(0.),
0063     fProjRange2(0.),
0064     fNbOfSteps(0),
0065     fNbOfSteps2(0),
0066     fStepSize(0.),
0067     fStepSize2(0.)
0068 {}
0069
0070 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0071
0072 Run::~Run() {}
0073
0074 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0075
0076 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
0077 {
0078   fParticle = particle;
0079   fEkin = energy;
0080 }
0081
0082 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0083
0084 void Run::AddInelastic(G4int nb)
0085 {
0086   fNbInelastic += nb;
0087   fNbInelastic2 += nb * nb;
0088 }
0089
0090 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0091
0092 void Run::AddEdep(G4double e)
0093 {
0094   fEdeposit += e;
0095   fEdeposit2 += e * e;
0096 }
0097
0098 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0099
0100 void Run::AddTrackLength(G4double t)
0101 {
0102   fTrackLen += t;
0103   fTrackLen2 += t * t;
0104 }
0105
0106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0107
0108 void Run::AddProjRange(G4double x)
0109 {
0110   fProjRange += x;
0111   fProjRange2 += x * x;
0112 }
0113
0114 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0115
0116 void Run::AddStepSize(G4int nb, G4double st)
0117 {
0118   fNbOfSteps += nb;
0119   fNbOfSteps2 += nb * nb;
0120   fStepSize += st;
0121   fStepSize2 += st * st;
0122 }
0123
0124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0125
0126 void Run::Merge(const G4Run* run)
0127 {
0128   const Run* localRun = static_cast<const Run*>(run);
0129
0130   // Pass information about primary particle
0131
0132   fParticle = localRun->fParticle;
0133   fEkin = localRun->fEkin;
0134
0135   // Accumulate sums
0136
0137   fNbInelastic += localRun->fNbInelastic;
0138   fNbInelastic2 += localRun->fNbInelastic2;
0139   fEdeposit += localRun->fEdeposit;
0140   fEdeposit2 += localRun->fEdeposit2;
0141   fTrackLen += localRun->fTrackLen;
0142   fTrackLen2 += localRun->fTrackLen2;
0143   fProjRange += localRun->fProjRange;
0144   fProjRange2 += localRun->fProjRange2;
0145   fNbOfSteps += localRun->fNbOfSteps;
0146   fNbOfSteps2 += localRun->fNbOfSteps2;
0147   fStepSize += localRun->fStepSize;
0148   fStepSize2 += localRun->fStepSize2;
0149
0150   G4Run::Merge(run);
0151 }
0152
0153 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0154
0155 void Run::EndOfRun()
0156 {
0157   std::ios::fmtflags mode = G4cout.flags();
0158   G4cout.setf(std::ios::fixed, std::ios::floatfield);
0159   G4int prec = G4cout.precision(2);
0160
0161   // Run conditions
0162
0163   G4Material* material = fDetector->GetAbsorMaterial();
0164   G4double density = material->GetDensity();
0165   G4String partName = fParticle->GetParticleName();
0166
0167   G4cout << "\n ======================== run summary =====================\n";
0168   G4cout << "\n The run is " << numberOfEvent << " " << partName << " of "
0169          << G4BestUnit(fEkin, "Energy") << " through a sphere of radius "
0170          << G4BestUnit(fDetector->GetAbsorRadius(), "Length") << "of " << material->GetName()
0171          << " (density: " << G4BestUnit(density, "Volumic Mass") << ")" << G4endl;
0172
0173   if (numberOfEvent == 0) {
0174     G4cout.setf(mode, std::ios::floatfield);
0175     G4cout.precision(prec);
0176     return;
0177   }
0178
0179   fNbInelastic /= numberOfEvent;
0180   fNbInelastic2 /= numberOfEvent;
0181
0182   G4double rms = fNbInelastic2 - fNbInelastic * fNbInelastic;
0183   if (rms > 0.)
0184     rms = std::sqrt(rms);
0185   else
0186     rms = 0.;
0187
0188   G4cout.precision(3);
0189   G4cout << "\n Nb of ionisations = " << fNbInelastic << " +- " << rms << G4endl;
0190
0191   G4cout.precision(3);
0192   G4cout << "\n w = " << G4BestUnit((fEkin) / fNbInelastic, "Energy") << " +- "
0193          << G4BestUnit((fEkin)*rms / (fNbInelastic * fNbInelastic), "Energy") << G4endl;
0194
0195   // Output file
0196
0197   if (fNbInelastic > 0.) {
0198     FILE* myFile;
0199     myFile = fopen("wvalue.txt", "a");
0200     fprintf(myFile, "%e %e %e %e %e \n", fEkin / eV, fNbInelastic, rms, fEkin / eV / fNbInelastic,
0201             (fEkin / eV) * rms / (fNbInelastic * fNbInelastic));
0202     fclose(myFile);
0203   }
0204   //
0205
0206   fEdeposit /= numberOfEvent;
0207   fEdeposit2 /= numberOfEvent;
0208   rms = fEdeposit2 - fEdeposit * fEdeposit;
0209   if (rms > 0.)
0210     rms = std::sqrt(rms);
0211   else
0212     rms = 0.;
0213
0214   G4cout.precision(3);
0215   G4cout << "\n Total Energy deposited        = " << G4BestUnit(fEdeposit, "Energy") << " +- "
0216          << G4BestUnit(rms, "Energy") << G4endl;
0217
0218   // Compute track length of primary track
0219
0220   fTrackLen /= numberOfEvent;
0221   fTrackLen2 /= numberOfEvent;
0222   rms = fTrackLen2 - fTrackLen * fTrackLen;
0223   if (rms > 0.)
0224     rms = std::sqrt(rms);
0225   else
0226     rms = 0.;
0227
0228   G4cout.precision(3);
0229   G4cout << "\n Track length of primary track = " << G4BestUnit(fTrackLen, "Length") << " +- "
0230          << G4BestUnit(rms, "Length");
0231
0232   // Compute projected range of primary track
0233
0234   fProjRange /= numberOfEvent;
0235   fProjRange2 /= numberOfEvent;
0236   rms = fProjRange2 - fProjRange * fProjRange;
0237   if (rms > 0.)
0238     rms = std::sqrt(rms);
0239   else
0240     rms = 0.;
0241
0242   G4cout << "\n Projected range               = " << G4BestUnit(fProjRange, "Length") << " +- "
0243          << G4BestUnit(rms, "Length") << G4endl;
0244
0245   // Nb of steps and step size of primary track
0246
0247   G4double dNofEvents = double(numberOfEvent);
0248   G4double fNbSteps = fNbOfSteps / dNofEvents, fNbSteps2 = fNbOfSteps2 / dNofEvents;
0249   rms = fNbSteps2 - fNbSteps * fNbSteps;
0250   if (rms > 0.)
0251     rms = std::sqrt(rms);
0252   else
0253     rms = 0.;
0254
0255   G4cout.precision(2);
0256   G4cout << "\n Nb of steps of primary track  = " << fNbSteps << " +- " << rms << G4endl;
0257
0258   fStepSize /= numberOfEvent;
0259   fStepSize2 /= numberOfEvent;
0260   rms = fStepSize2 - fStepSize * fStepSize;
0261   if (rms > 0.)
0262     rms = std::sqrt(rms);
0263   else
0264     rms = 0.;
0265
0266   G4cout.precision(3);
0267   G4cout << "\n Step size                     = " << G4BestUnit(fStepSize, "Length") << " +- "
0268          << G4BestUnit(rms, "Length") << G4endl;
0269
0270   // Normalize histograms of longitudinal energy profile
0271
0272   G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0273   G4int ih = 1;
0274   G4double binWidth = analysisManager->GetH1Width(ih);
0275   G4double fac = (1. / (numberOfEvent * binWidth)) * (mm / MeV);
0276   analysisManager->ScaleH1(ih, fac);
0277
0278   // Reset default formats
0279
0280   G4cout.setf(mode, std::ios::floatfield);
0281   G4cout.precision(prec);
0282 }