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   fNIEL += localRun->fNIEL;
0086   fTrueRange += localRun->fTrueRange;
0087   fTrueRange2 += localRun->fTrueRange2;
0088   fProjRange += localRun->fProjRange;
0089   fProjRange2 += localRun->fProjRange2;
0090   fTransvDev += localRun->fTransvDev;
0091   fTransvDev2 += localRun->fTransvDev2;
0092
0093   // map: processes count
0094   std::map<G4String, G4int>::const_iterator it;
0095   for (it = localRun->fProcCounter.begin(); it != localRun->fProcCounter.end(); ++it) {
0096     G4String procName = it->first;
0097     G4int localCount = it->second;
0098     if (fProcCounter.find(procName) == fProcCounter.end()) {
0099       fProcCounter[procName] = localCount;
0100     }
0101     else {
0102       fProcCounter[procName] += localCount;
0103     }
0104   }
0105
0106   G4Run::Merge(run);
0107 }
0108
0109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0110
0111 void Run::EndOfRun()
0112 {
0113   G4int prec = 5, wid = prec + 2;
0114   G4int dfprec = G4cout.precision(prec);
0115
0116   // run condition
0117   //
0118   G4String partName = fParticle->GetParticleName();
0119   const G4Material* material = fDetector->GetMaterial();
0120   G4double density = material->GetDensity();
0121
0122   G4cout << "\n ======================== run summary ======================\n";
0123   G4cout << "\n The run is: " << numberOfEvent << " " << partName << " of "
0124          << G4BestUnit(fEkin, "Energy") << " through " << G4BestUnit(fDetector->GetSize(), "Length")
0125          << " of " << material->GetName() << " (density: " << G4BestUnit(density, "Volumic Mass")
0126          << ")" << G4endl;
0127
0128   if (numberOfEvent == 0) {
0129     G4cout.precision(dfprec);
0130     return;
0131   }
0132
0133   G4double dNbOfEvents = (G4double)numberOfEvent;
0134   G4cout << "\n Total energy deposit:   " << G4BestUnit(fEdep / dNbOfEvents, "Energy") << G4endl;
0135   G4cout << " NIEL energy calculated: " << G4BestUnit(fNIEL / dNbOfEvents, "Energy") << G4endl;
0136
0137   // nb of tracks and steps per event
0138   //
0139   G4cout << "\n Nb tracks/event"
0140          << "   neutral: " << std::setw(wid) << fNbOfTraks0 / dNbOfEvents
0141          << "   charged: " << std::setw(wid) << fNbOfTraks1 / dNbOfEvents << "\n Nb  steps/event"
0142          << "   neutral: " << std::setw(wid) << fNbOfSteps0 / dNbOfEvents
0143          << "   charged: " << std::setw(wid) << fNbOfSteps1 / dNbOfEvents << G4endl;
0144
0145   // frequency of processes
0146   //
0147   G4cout << "\n Process calls frequency :" << G4endl;
0148   G4int index = 0;
0149   std::map<G4String, G4int>::iterator it;
0150   for (it = fProcCounter.begin(); it != fProcCounter.end(); it++) {
0151     G4String procName = it->first;
0152     G4int count = it->second;
0153     G4String space = " ";
0154     if (++index % 3 == 0) space = "\n";
0155     G4cout << " " << std::setw(20) << procName << "=" << std::setw(7) << count << space;
0156   }
0157   G4cout << G4endl;
0158
0159   // compute true and projected ranges, and transverse dispersion
0160   //
0161   fTrueRange /= numberOfEvent;
0162   fTrueRange2 /= numberOfEvent;
0163   G4double trueRms = fTrueRange2 - fTrueRange * fTrueRange;
0164   if (trueRms > 0.)
0165     trueRms = std::sqrt(trueRms);
0166   else
0167     trueRms = 0.;
0168
0169   fProjRange /= numberOfEvent;
0170   fProjRange2 /= numberOfEvent;
0171   G4double projRms = fProjRange2 - fProjRange * fProjRange;
0172   if (projRms > 0.)
0173     projRms = std::sqrt(projRms);
0174   else
0175     projRms = 0.;
0176
0177   fTransvDev /= 2 * numberOfEvent;
0178   fTransvDev2 /= 2 * numberOfEvent;
0179   G4double trvsRms = fTransvDev2 - fTransvDev * fTransvDev;
0180   if (trvsRms > 0.)
0181     trvsRms = std::sqrt(trvsRms);
0182   else
0183     trvsRms = 0.;
0184
0185   // compare true range with csda range from PhysicsTables
0186   //
0187   G4EmCalculator emCalculator;
0188   G4double rangeTable = 0.;
0189   if (fParticle->GetPDGCharge() != 0.)
0190     rangeTable = emCalculator.GetCSDARange(fEkin, fParticle, material);
0191
0192   G4cout << "\n---------------------------------------------------------\n";
0193   G4cout << " Primary particle : ";
0194   G4cout << "\n true Range = " << G4BestUnit(fTrueRange, "Length")
0195          << "   rms = " << G4BestUnit(trueRms, "Length");
0196
0197   G4cout << "\n proj Range = " << G4BestUnit(fProjRange, "Length")
0198          << "   rms = " << G4BestUnit(projRms, "Length");
0199
0200   G4cout << "\n proj/true  = " << fProjRange / fTrueRange;
0201
0202   G4cout << "\n transverse dispersion at end = " << G4BestUnit(trvsRms, "Length");
0203
0204   G4cout << "\n      mass true Range from simulation = "
0205          << G4BestUnit(fTrueRange * density, "Mass/Surface")
0206          << "\n       from PhysicsTable (csda range) = "
0207          << G4BestUnit(rangeTable * density, "Mass/Surface");
0208   G4cout << "\n---------------------------------------------------------\n";
0209   G4cout << G4endl;
0210
0211   // remove all contents in fProcCounter
0212   fProcCounter.clear();
0213
0214   // restore default format
0215   G4cout.precision(dfprec);
0216 }
0217
0218 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......