RE07/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 #include "Run.hh"
0030
0031 #include "DetectorConstruction.hh"
0032 #include "PrimaryGeneratorAction.hh"
0033
0034 #include "G4Electron.hh"
0035 #include "G4Gamma.hh"
0036 #include "G4ParticleDefinition.hh"
0037 #include "G4ParticleTable.hh"
0038 #include "G4Positron.hh"
0039 #include "G4SystemOfUnits.hh"
0040 #include "G4Track.hh"
0041 #include "G4UnitsTable.hh"
0042
0043 #include <iomanip>
0044
0045 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0046
0047 Run::Run(DetectorConstruction* det)
0048   : G4Run(),
0049     fDetector(det),
0050     fParticle(nullptr),
0051     fEkin(0.),
0052     fChargedStep(0),
0053     fNeutralStep(0),
0054     fN_gamma(0),
0055     fN_elec(0),
0056     fN_pos(0)
0057 {
0058   // initialize cumulative quantities
0059   //
0060   for (G4int k = 0; k < kMaxAbsor; k++) {
0061     fSumEAbs[k] = fSum2EAbs[k] = fSumLAbs[k] = fSum2LAbs[k] = 0.;
0062     fEnergyDeposit[k].clear();
0063   }
0064 }
0065
0066 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0067
0068 Run::~Run() {}
0069
0070 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0071
0072 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
0073 {
0074   fParticle = particle;
0075   fEkin = energy;
0076 }
0077
0078 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0079
0080 void Run::FillPerEvent(G4int kAbs, G4double EAbs, G4double LAbs)
0081 {
0082   // accumulate statistic with restriction
0083   //
0084   fEnergyDeposit[kAbs].push_back(EAbs);
0085   fSumEAbs[kAbs] += EAbs;
0086   fSum2EAbs[kAbs] += EAbs * EAbs;
0087   fSumLAbs[kAbs] += LAbs;
0088   fSum2LAbs[kAbs] += LAbs * LAbs;
0089 }
0090
0091 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0092
0093 void Run::AddChargedStep()
0094 {
0095   fChargedStep += 1.0;
0096 }
0097
0098 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0099
0100 void Run::AddNeutralStep()
0101 {
0102   fNeutralStep += 1.0;
0103 }
0104
0105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0106
0107 void Run::AddSecondaryTrack(const G4Track* track)
0108 {
0109   const G4ParticleDefinition* d = track->GetDefinition();
0110   if (d == G4Gamma::Gamma()) {
0111     ++fN_gamma;
0112   }
0113   else if (d == G4Electron::Electron()) {
0114     ++fN_elec;
0115   }
0116   else if (d == G4Positron::Positron()) {
0117     ++fN_pos;
0118   }
0119 }
0120
0121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0122
0123 void Run::Merge(const G4Run* run)
0124 {
0125   const Run* localRun = static_cast<const Run*>(run);
0126
0127   // pass information about primary particle
0128   fParticle = localRun->fParticle;
0129   fEkin = localRun->fEkin;
0130
0131   // accumulate sums
0132   //
0133   for (G4int k = 0; k < kMaxAbsor; k++) {
0134     fSumEAbs[k] += localRun->fSumEAbs[k];
0135     fSum2EAbs[k] += localRun->fSum2EAbs[k];
0136     fSumLAbs[k] += localRun->fSumLAbs[k];
0137     fSum2LAbs[k] += localRun->fSum2LAbs[k];
0138   }
0139
0140   fChargedStep += localRun->fChargedStep;
0141   fNeutralStep += localRun->fNeutralStep;
0142
0143   fN_gamma += localRun->fN_gamma;
0144   fN_elec += localRun->fN_elec;
0145   fN_pos += localRun->fN_pos;
0146
0147   G4Run::Merge(run);
0148 }
0149
0150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0151
0152 void Run::EndOfRun()
0153 {
0154   G4int nEvt = numberOfEvent;
0155   G4double norm = G4double(nEvt);
0156   if (norm > 0) norm = 1. / norm;
0157   G4double qnorm = std::sqrt(norm);
0158
0159   fChargedStep *= norm;
0160   fNeutralStep *= norm;
0161
0162   // compute and print statistic
0163   //
0164   G4double beamEnergy = fEkin;
0165   G4double sqbeam = std::sqrt(beamEnergy / GeV);
0166
0167   G4double MeanEAbs, MeanEAbs2, rmsEAbs, resolution, rmsres;
0168   G4double MeanLAbs, MeanLAbs2, rmsLAbs;
0169
0170   std::ios::fmtflags mode = G4cout.flags();
0171   G4int prec = G4cout.precision(2);
0172   G4cout << "\n------------------------------------------------------------\n";
0173   G4cout << std::setw(14) << "material" << std::setw(17) << "Edep       RMS" << std::setw(33)
0174          << "sqrt(E0(GeV))*rmsE/Emean" << std::setw(23) << "total tracklen \n \n";
0175
0176   for (G4int k = 1; k <= fDetector->GetNbOfAbsor(); k++) {
0177     MeanEAbs = fSumEAbs[k] * norm;
0178     MeanEAbs2 = fSum2EAbs[k] * norm;
0179     rmsEAbs = std::sqrt(std::abs(MeanEAbs2 - MeanEAbs * MeanEAbs));
0180
0181     resolution = 100. * sqbeam * rmsEAbs / MeanEAbs;
0182     rmsres = resolution * qnorm;
0183
0184     // Save mean and RMS
0185     fSumEAbs[k] = MeanEAbs;
0186     fSum2EAbs[k] = rmsEAbs;
0187
0188     MeanLAbs = fSumLAbs[k] * norm;
0189     MeanLAbs2 = fSum2LAbs[k] * norm;
0190     rmsLAbs = std::sqrt(std::abs(MeanLAbs2 - MeanLAbs * MeanLAbs));
0191
0192     // print
0193     //
0194     G4cout << std::setw(14) << fDetector->GetAbsorMaterial(k)->GetName() << ": "
0195            << std::setprecision(5) << std::setw(6) << G4BestUnit(MeanEAbs, "Energy") << " :  "
0196            << std::setprecision(4) << std::setw(5) << G4BestUnit(rmsEAbs, "Energy") << std::setw(10)
0197            << resolution << " +- " << std::setw(5) << rmsres << " %" << std::setprecision(3)
0198            << std::setw(10) << G4BestUnit(MeanLAbs, "Length") << " +- " << std::setw(4)
0199            << G4BestUnit(rmsLAbs, "Length") << G4endl;
0200   }
0201   G4cout << "\n------------------------------------------------------------\n";
0202
0203   G4cout << " Beam particle " << fParticle->GetParticleName()
0204          << "  E = " << G4BestUnit(beamEnergy, "Energy") << G4endl;
0205   G4cout << " Mean number of gamma       " << (G4double)fN_gamma * norm << G4endl;
0206   G4cout << " Mean number of e-          " << (G4double)fN_elec * norm << G4endl;
0207   G4cout << " Mean number of e+          " << (G4double)fN_pos * norm << G4endl;
0208   G4cout << std::setprecision(6) << " Mean number of charged steps  " << fChargedStep << G4endl;
0209   G4cout << " Mean number of neutral steps  " << fNeutralStep << G4endl;
0210   G4cout << "------------------------------------------------------------\n" << G4endl;
0211
0212   G4cout.setf(mode, std::ios::floatfield);
0213   G4cout.precision(prec);
0214 }
0215
0216 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......