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