Geant4 Cross Reference |
1 // 2 // ******************************************************************** 3 // * License and Disclaimer * 4 // * * 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. * 10 // * * 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitation of liability. * 17 // * * 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************************************** 25 // 26 /// \file electromagnetic/TestEm1/src/Run.cc 27 /// \brief Implementation of the Run class 28 // 29 // 30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 32 33 #include "Run.hh" 34 35 #include "DetectorConstruction.hh" 36 #include "PrimaryGeneratorAction.hh" 37 38 #include "G4EmCalculator.hh" 39 #include "G4SystemOfUnits.hh" 40 #include "G4UnitsTable.hh" 41 42 #include <iomanip> 43 44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 45 46 Run::Run(const DetectorConstruction* det) : fDetector(det) {} 47 48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 49 50 void Run::SetPrimary(const G4ParticleDefinition* particle, G4double energy) 51 { 52 fParticle = particle; 53 fEkin = energy; 54 } 55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 56 57 void Run::CountProcesses(const G4String& procName) 58 { 59 std::map<G4String, G4int>::iterator it = fProcCounter.find(procName); 60 if (it == fProcCounter.end()) { 61 fProcCounter[procName] = 1; 62 } 63 else { 64 fProcCounter[procName]++; 65 } 66 } 67 68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 69 70 void Run::Merge(const G4Run* run) 71 { 72 const Run* localRun = static_cast<const Run*>(run); 73 74 // pass information about primary particle 75 fParticle = localRun->fParticle; 76 fEkin = localRun->fEkin; 77 78 // accumulate sums 79 // 80 fNbOfTraks0 += localRun->fNbOfTraks0; 81 fNbOfTraks1 += localRun->fNbOfTraks1; 82 fNbOfSteps0 += localRun->fNbOfSteps0; 83 fNbOfSteps1 += localRun->fNbOfSteps1; 84 fEdep += localRun->fEdep; 85 fEleak += localRun->fEleak; 86 fNIEL += localRun->fNIEL; 87 fTrueRange += localRun->fTrueRange; 88 fTrueRange2 += localRun->fTrueRange2; 89 fProjRange += localRun->fProjRange; 90 fProjRange2 += localRun->fProjRange2; 91 fTransvDev += localRun->fTransvDev; 92 fTransvDev2 += localRun->fTransvDev2; 93 94 // map: processes count 95 std::map<G4String, G4int>::const_iterator it; 96 for (it = localRun->fProcCounter.begin(); it != localRun->fProcCounter.end(); ++it) { 97 G4String procName = it->first; 98 G4int localCount = it->second; 99 if (fProcCounter.find(procName) == fProcCounter.end()) { 100 fProcCounter[procName] = localCount; 101 } 102 else { 103 fProcCounter[procName] += localCount; 104 } 105 } 106 107 G4Run::Merge(run); 108 } 109 110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 111 112 void Run::EndOfRun() 113 { 114 G4int prec = 5, wid = prec + 2; 115 G4int dfprec = G4cout.precision(prec); 116 117 // run condition 118 // 119 G4String partName = fParticle->GetParticleName(); 120 const G4Material* material = fDetector->GetMaterial(); 121 G4double density = material->GetDensity(); 122 123 G4cout << "\n ======================== run summary ======================\n"; 124 G4cout << "\n The run is: " << numberOfEvent << " " << partName << " of " 125 << G4BestUnit(fEkin, "Energy") << " through " << G4BestUnit(fDetector->GetSize(), "Length") 126 << " of " << material->GetName() << " (density: " << G4BestUnit(density, "Volumic Mass") 127 << ")" << G4endl; 128 129 if (numberOfEvent == 0) { 130 G4cout.precision(dfprec); 131 return; 132 } 133 134 G4double dNbOfEvents = (G4double)numberOfEvent; 135 G4cout << "\n Energy deposit: " 136 << G4BestUnit(fEdep/dNbOfEvents, "Energy") << G4endl; 137 G4cout << " Energy leakage: " 138 << G4BestUnit(fEleak/dNbOfEvents, "Energy") << G4endl; 139 G4cout << " Edep + Eleak: " 140 << G4BestUnit((fEdep+fEleak)/dNbOfEvents, "Energy") << G4endl; 141 G4cout << " \n NIEL energy calculated: " 142 << G4BestUnit(fNIEL/dNbOfEvents, "Energy") << G4endl; 143 144 // nb of tracks and steps per event 145 // 146 G4cout << "\n Nb tracks/event" 147 << " neutral: " << std::setw(wid) << fNbOfTraks0 / dNbOfEvents 148 << " charged: " << std::setw(wid) << fNbOfTraks1 / dNbOfEvents << "\n Nb steps/event" 149 << " neutral: " << std::setw(wid) << fNbOfSteps0 / dNbOfEvents 150 << " charged: " << std::setw(wid) << fNbOfSteps1 / dNbOfEvents << G4endl; 151 152 // frequency of processes 153 // 154 G4cout << "\n Process calls frequency :" << G4endl; 155 G4int index = 0; 156 std::map<G4String, G4int>::iterator it; 157 for (it = fProcCounter.begin(); it != fProcCounter.end(); it++) { 158 G4String procName = it->first; 159 G4int count = it->second; 160 G4String space = " "; 161 if (++index % 3 == 0) space = "\n"; 162 G4cout << " " << std::setw(20) << procName << "=" << std::setw(7) << count << space; 163 } 164 G4cout << G4endl; 165 166 // compute true and projected ranges, and transverse dispersion 167 // 168 fTrueRange /= numberOfEvent; 169 fTrueRange2 /= numberOfEvent; 170 G4double trueRms = fTrueRange2 - fTrueRange * fTrueRange; 171 if (trueRms > 0.) 172 trueRms = std::sqrt(trueRms); 173 else 174 trueRms = 0.; 175 176 fProjRange /= numberOfEvent; 177 fProjRange2 /= numberOfEvent; 178 G4double projRms = fProjRange2 - fProjRange * fProjRange; 179 if (projRms > 0.) 180 projRms = std::sqrt(projRms); 181 else 182 projRms = 0.; 183 184 fTransvDev /= 2 * numberOfEvent; 185 fTransvDev2 /= 2 * numberOfEvent; 186 G4double trvsRms = fTransvDev2 - fTransvDev * fTransvDev; 187 if (trvsRms > 0.) 188 trvsRms = std::sqrt(trvsRms); 189 else 190 trvsRms = 0.; 191 192 // compare true range with csda range from PhysicsTables 193 // 194 G4EmCalculator emCalculator; 195 G4double rangeTable = 0.; 196 if (fParticle->GetPDGCharge() != 0.) 197 rangeTable = emCalculator.GetCSDARange(fEkin, fParticle, material); 198 199 G4cout << "\n---------------------------------------------------------\n"; 200 G4cout << " Primary particle : "; 201 G4cout << "\n true Range = " << G4BestUnit(fTrueRange, "Length") 202 << " rms = " << G4BestUnit(trueRms, "Length"); 203 204 G4cout << "\n proj Range = " << G4BestUnit(fProjRange, "Length") 205 << " rms = " << G4BestUnit(projRms, "Length"); 206 207 G4cout << "\n proj/true = " << fProjRange / fTrueRange; 208 209 G4cout << "\n transverse dispersion at end = " << G4BestUnit(trvsRms, "Length"); 210 211 G4cout << "\n mass true Range from simulation = " 212 << G4BestUnit(fTrueRange * density, "Mass/Surface") 213 << "\n from PhysicsTable (csda range) = " 214 << G4BestUnit(rangeTable * density, "Mass/Surface"); 215 G4cout << "\n---------------------------------------------------------\n"; 216 G4cout << G4endl; 217 218 // remove all contents in fProcCounter 219 fProcCounter.clear(); 220 221 // restore default format 222 G4cout.precision(dfprec); 223 } 224 225 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 226