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