TestEm1/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 electromagnetic/TestEm1/src/Run.cc
0027 /// \brief Implementation of the Run class
0028 //
0029 //
0030 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0031 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0032
0033 #include "Run.hh"
0034
0035 #include "DetectorConstruction.hh"
0036 #include "PrimaryGeneratorAction.hh"
0037
0038 #include "G4EmCalculator.hh"
0039 #include "G4SystemOfUnits.hh"
0040 #include "G4UnitsTable.hh"
0041
0042 #include <iomanip>
0043
0044 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0045
0046 Run::Run(const DetectorConstruction* det) : fDetector(det) {}
0047
0048 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0049
0050 void Run::SetPrimary(const G4ParticleDefinition* particle, G4double energy)
0051 {
0052   fParticle = particle;
0053   fEkin = energy;
0054 }
0055 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0056
0057 void Run::CountProcesses(const G4String& procName)
0058 {
0059   std::map<G4String, G4int>::iterator it = fProcCounter.find(procName);
0060   if (it == fProcCounter.end()) {
0061     fProcCounter[procName] = 1;
0062   }
0063   else {
0064     fProcCounter[procName]++;
0065   }
0066 }
0067
0068 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0069
0070 void Run::Merge(const G4Run* run)
0071 {
0072   const Run* localRun = static_cast<const Run*>(run);
0073
0074   // pass information about primary particle
0075   fParticle = localRun->fParticle;
0076   fEkin = localRun->fEkin;
0077
0078   // accumulate sums
0079   //
0080   fNbOfTraks0 += localRun->fNbOfTraks0;
0081   fNbOfTraks1 += localRun->fNbOfTraks1;
0082   fNbOfSteps0 += localRun->fNbOfSteps0;
0083   fNbOfSteps1 += localRun->fNbOfSteps1;
0084   fEdep += localRun->fEdep;
0085   fEleak += localRun->fEleak;
0086   fNIEL += localRun->fNIEL;
0087   fTrueRange += localRun->fTrueRange;
0088   fTrueRange2 += localRun->fTrueRange2;
0089   fProjRange += localRun->fProjRange;
0090   fProjRange2 += localRun->fProjRange2;
0091   fTransvDev += localRun->fTransvDev;
0092   fTransvDev2 += localRun->fTransvDev2;
0093
0094   // map: processes count
0095   std::map<G4String, G4int>::const_iterator it;
0096   for (it = localRun->fProcCounter.begin(); it != localRun->fProcCounter.end(); ++it) {
0097     G4String procName = it->first;
0098     G4int localCount = it->second;
0099     if (fProcCounter.find(procName) == fProcCounter.end()) {
0100       fProcCounter[procName] = localCount;
0101     }
0102     else {
0103       fProcCounter[procName] += localCount;
0104     }
0105   }
0106
0107   G4Run::Merge(run);
0108 }
0109
0110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0111
0112 void Run::EndOfRun()
0113 {
0114   G4int prec = 5, wid = prec + 2;
0115   G4int dfprec = G4cout.precision(prec);
0116
0117   // run condition
0118   //
0119   G4String partName = fParticle->GetParticleName();
0120   const G4Material* material = fDetector->GetMaterial();
0121   G4double density = material->GetDensity();
0122
0123   G4cout << "\n ======================== run summary ======================\n";
0124   G4cout << "\n The run is: " << numberOfEvent << " " << partName << " of "
0125          << G4BestUnit(fEkin, "Energy") << " through " << G4BestUnit(fDetector->GetSize(), "Length")
0126          << " of " << material->GetName() << " (density: " << G4BestUnit(density, "Volumic Mass")
0127          << ")" << G4endl;
0128
0129   if (numberOfEvent == 0) {
0130     G4cout.precision(dfprec);
0131     return;
0132   }
0133
0134   G4double dNbOfEvents = (G4double)numberOfEvent;
0135   G4cout << "\n Energy deposit:   "
0136          << G4BestUnit(fEdep/dNbOfEvents,  "Energy") << G4endl;
0137   G4cout << " Energy leakage:   "
0138          << G4BestUnit(fEleak/dNbOfEvents, "Energy") << G4endl;
0139   G4cout << " Edep + Eleak:     "
0140          << G4BestUnit((fEdep+fEleak)/dNbOfEvents, "Energy") << G4endl;
0141   G4cout << " \n NIEL energy calculated: "
0142          << G4BestUnit(fNIEL/dNbOfEvents,  "Energy") << G4endl;
0143
0144   // nb of tracks and steps per event
0145   //
0146   G4cout << "\n Nb tracks/event"
0147          << "   neutral: " << std::setw(wid) << fNbOfTraks0 / dNbOfEvents
0148          << "   charged: " << std::setw(wid) << fNbOfTraks1 / dNbOfEvents << "\n Nb  steps/event"
0149          << "   neutral: " << std::setw(wid) << fNbOfSteps0 / dNbOfEvents
0150          << "   charged: " << std::setw(wid) << fNbOfSteps1 / dNbOfEvents << G4endl;
0151
0152   // frequency of processes
0153   //
0154   G4cout << "\n Process calls frequency :" << G4endl;
0155   G4int index = 0;
0156   std::map<G4String, G4int>::iterator it;
0157   for (it = fProcCounter.begin(); it != fProcCounter.end(); it++) {
0158     G4String procName = it->first;
0159     G4int count = it->second;
0160     G4String space = " ";
0161     if (++index % 3 == 0) space = "\n";
0162     G4cout << " " << std::setw(20) << procName << "=" << std::setw(7) << count << space;
0163   }
0164   G4cout << G4endl;
0165
0166   // compute true and projected ranges, and transverse dispersion
0167   //
0168   fTrueRange /= numberOfEvent;
0169   fTrueRange2 /= numberOfEvent;
0170   G4double trueRms = fTrueRange2 - fTrueRange * fTrueRange;
0171   if (trueRms > 0.)
0172     trueRms = std::sqrt(trueRms);
0173   else
0174     trueRms = 0.;
0175
0176   fProjRange /= numberOfEvent;
0177   fProjRange2 /= numberOfEvent;
0178   G4double projRms = fProjRange2 - fProjRange * fProjRange;
0179   if (projRms > 0.)
0180     projRms = std::sqrt(projRms);
0181   else
0182     projRms = 0.;
0183
0184   fTransvDev /= 2 * numberOfEvent;
0185   fTransvDev2 /= 2 * numberOfEvent;
0186   G4double trvsRms = fTransvDev2 - fTransvDev * fTransvDev;
0187   if (trvsRms > 0.)
0188     trvsRms = std::sqrt(trvsRms);
0189   else
0190     trvsRms = 0.;
0191
0192   // compare true range with csda range from PhysicsTables
0193   //
0194   G4EmCalculator emCalculator;
0195   G4double rangeTable = 0.;
0196   if (fParticle->GetPDGCharge() != 0.)
0197     rangeTable = emCalculator.GetCSDARange(fEkin, fParticle, material);
0198
0199   G4cout << "\n---------------------------------------------------------\n";
0200   G4cout << " Primary particle : ";
0201   G4cout << "\n true Range = " << G4BestUnit(fTrueRange, "Length")
0202          << "   rms = " << G4BestUnit(trueRms, "Length");
0203
0204   G4cout << "\n proj Range = " << G4BestUnit(fProjRange, "Length")
0205          << "   rms = " << G4BestUnit(projRms, "Length");
0206
0207   G4cout << "\n proj/true  = " << fProjRange / fTrueRange;
0208
0209   G4cout << "\n transverse dispersion at end = " << G4BestUnit(trvsRms, "Length");
0210
0211   G4cout << "\n      mass true Range from simulation = "
0212          << G4BestUnit(fTrueRange * density, "Mass/Surface")
0213          << "\n       from PhysicsTable (csda range) = "
0214          << G4BestUnit(rangeTable * density, "Mass/Surface");
0215   G4cout << "\n---------------------------------------------------------\n";
0216   G4cout << G4endl;
0217
0218   // remove all contents in fProcCounter
0219   fProcCounter.clear();
0220
0221   // restore default format
0222   G4cout.precision(dfprec);
0223 }
0224
0225 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......