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