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