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 Run.cc 26 /// \file Run.cc 27 /// \brief Implementation of the Run class 27 /// \brief Implementation of the Run class 28 // 28 // 29 // << 29 // $Id: Run.cc 71376 2013-06-14 07:44:50Z maire $ >> 30 // 30 //....oooOO0OOooo........oooOO0OOooo........oo 31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 31 //....oooOO0OOooo........oooOO0OOooo........oo 32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 32 33 33 #include "Run.hh" 34 #include "Run.hh" 34 << 35 #include "DetectorConstruction.hh" 35 #include "DetectorConstruction.hh" 36 #include "HistoManager.hh" << 37 #include "PrimaryGeneratorAction.hh" 36 #include "PrimaryGeneratorAction.hh" >> 37 #include "HistoManager.hh" 38 38 39 #include "G4HadronicProcess.hh" << 40 #include "G4HadronicProcessStore.hh" << 41 #include "G4Neutron.hh" << 42 #include "G4ProcessTable.hh" 39 #include "G4ProcessTable.hh" 43 #include "G4SystemOfUnits.hh" << 40 #include "G4HadronicProcessStore.hh" 44 #include "G4UnitsTable.hh" 41 #include "G4UnitsTable.hh" >> 42 #include "G4SystemOfUnits.hh" 45 43 46 //....oooOO0OOooo........oooOO0OOooo........oo 44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 47 45 48 Run::Run(DetectorConstruction* det) : fDetecto << 46 Run::Run(DetectorConstruction* det) >> 47 : G4Run(), >> 48 fDetector(det), fParticle(0), fEkin(0.), >> 49 fTotalCount(0), fGammaCount(0), >> 50 fSumTrack(0.), fSumTrack2(0.), >> 51 fTargetXXX(false) 49 { 52 { 50 for (G4int i = 0; i < 3; i++) { << 53 for (G4int i=0; i<3; i++) { fPbalance[i] = 0. ; } 51 fPbalance[i] = 0.; << 54 for (G4int i=0; i<3; i++) { fNbGamma[i] = 0 ; } 52 } << 53 for (G4int i = 0; i < 3; i++) { << 54 fNbGamma[i] = 0; << 55 } << 56 fPbalance[1] = DBL_MAX; << 57 fNbGamma[1] = 10000; << 58 } 55 } >> 56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 57 >> 58 Run::~Run() >> 59 { } 59 60 60 //....oooOO0OOooo........oooOO0OOooo........oo 61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 61 62 62 void Run::SetPrimary(G4ParticleDefinition* par 63 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy) 63 { << 64 { 64 fParticle = particle; 65 fParticle = particle; 65 fEkin = energy; 66 fEkin = energy; 66 } << 67 } 67 68 68 //....oooOO0OOooo........oooOO0OOooo........oo 69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 69 70 70 void Run::SetTargetXXX(G4bool flag) 71 void Run::SetTargetXXX(G4bool flag) 71 { << 72 { 72 fTargetXXX = flag; 73 fTargetXXX = flag; 73 } 74 } 74 << 75 75 //....oooOO0OOooo........oooOO0OOooo........oo 76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 76 77 77 void Run::CountProcesses(G4VProcess* process) << 78 void Run::CountProcesses(G4VProcess* process) 78 { 79 { 79 if (process == nullptr) return; << 80 G4String procName = process->GetProcessName( 80 G4String procName = process->GetProcessName(); 81 std::map<G4String, G4int>::iterator it = fPr << 81 std::map<G4String,G4int>::iterator it = fProcCounter.find(procName); 82 if (it == fProcCounter.end()) { << 82 if ( it == fProcCounter.end()) { 83 fProcCounter[procName] = 1; 83 fProcCounter[procName] = 1; 84 } 84 } 85 else { 85 else { 86 fProcCounter[procName]++; << 86 fProcCounter[procName]++; 87 } 87 } 88 } << 88 } 89 << 90 //....oooOO0OOooo........oooOO0OOooo........oo 89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 91 90 92 void Run::SumTrack(G4double trackl) 91 void Run::SumTrack(G4double trackl) 93 { 92 { 94 fTotalCount++; 93 fTotalCount++; 95 fSumTrack += trackl; << 94 fSumTrack += trackl; fSumTrack2 += trackl*trackl; 96 fSumTrack2 += trackl * trackl; << 97 } 95 } 98 96 99 //....oooOO0OOooo........oooOO0OOooo........oo 97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 100 98 101 void Run::CountNuclearChannel(G4String name, G 99 void Run::CountNuclearChannel(G4String name, G4double Q) 102 { 100 { 103 std::map<G4String, NuclChannel>::iterator it 101 std::map<G4String, NuclChannel>::iterator it = fNuclChannelMap.find(name); 104 if (it == fNuclChannelMap.end()) { << 102 if ( it == fNuclChannelMap.end()) { 105 fNuclChannelMap[name] = NuclChannel(1, Q); 103 fNuclChannelMap[name] = NuclChannel(1, Q); 106 } 104 } 107 else { 105 else { 108 NuclChannel& data = it->second; 106 NuclChannel& data = it->second; 109 data.fCount++; 107 data.fCount++; 110 data.fQ += Q; 108 data.fQ += Q; 111 } << 109 } 112 } 110 } 113 << 111 114 //....oooOO0OOooo........oooOO0OOooo........oo 112 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 115 113 116 void Run::ParticleCount(G4String name, G4doubl 114 void Run::ParticleCount(G4String name, G4double Ekin) 117 { 115 { 118 std::map<G4String, ParticleData>::iterator i 116 std::map<G4String, ParticleData>::iterator it = fParticleDataMap.find(name); 119 if (it == fParticleDataMap.end()) { << 117 if ( it == fParticleDataMap.end()) { 120 fParticleDataMap[name] = ParticleData(1, E 118 fParticleDataMap[name] = ParticleData(1, Ekin, Ekin, Ekin); 121 } 119 } 122 else { 120 else { 123 ParticleData& data = it->second; 121 ParticleData& data = it->second; 124 data.fCount++; 122 data.fCount++; 125 data.fEmean += Ekin; 123 data.fEmean += Ekin; 126 // update min max << 124 //update min max 127 G4double emin = data.fEmin; 125 G4double emin = data.fEmin; 128 if (Ekin < emin) data.fEmin = Ekin; 126 if (Ekin < emin) data.fEmin = Ekin; 129 G4double emax = data.fEmax; 127 G4double emax = data.fEmax; 130 if (Ekin > emax) data.fEmax = Ekin; << 128 if (Ekin > emax) data.fEmax = Ekin; 131 } << 129 } 132 } 130 } 133 //....oooOO0OOooo........oooOO0OOooo........oo 131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 134 132 135 void Run::Balance(G4double Pbal) 133 void Run::Balance(G4double Pbal) 136 { << 134 { 137 fPbalance[0] += Pbal; 135 fPbalance[0] += Pbal; 138 // update min max << 136 //update min max 139 if (fTotalCount == 1) fPbalance[1] = fPbalan << 137 if (fTotalCount == 1) fPbalance[1] = fPbalance[2] = Pbal; 140 if (Pbal < fPbalance[1]) fPbalance[1] = Pbal 138 if (Pbal < fPbalance[1]) fPbalance[1] = Pbal; 141 if (Pbal > fPbalance[2]) fPbalance[2] = Pbal << 139 if (Pbal > fPbalance[2]) fPbalance[2] = Pbal; 142 } 140 } 143 141 144 //....oooOO0OOooo........oooOO0OOooo........oo 142 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 145 143 146 void Run::CountGamma(G4int nGamma) 144 void Run::CountGamma(G4int nGamma) 147 { << 145 { 148 fGammaCount++; 146 fGammaCount++; 149 fNbGamma[0] += nGamma; 147 fNbGamma[0] += nGamma; 150 // update min max << 148 //update min max 151 if (fGammaCount == 1) fNbGamma[1] = fNbGamma << 149 if (fGammaCount == 1) fNbGamma[1] = fNbGamma[2] = nGamma; 152 if (nGamma < fNbGamma[1]) fNbGamma[1] = nGam 150 if (nGamma < fNbGamma[1]) fNbGamma[1] = nGamma; 153 if (nGamma > fNbGamma[2]) fNbGamma[2] = nGam << 151 if (nGamma > fNbGamma[2]) fNbGamma[2] = nGamma; 154 } 152 } 155 153 156 //....oooOO0OOooo........oooOO0OOooo........oo 154 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 157 155 158 void Run::Merge(const G4Run* run) 156 void Run::Merge(const G4Run* run) 159 { 157 { 160 const Run* localRun = static_cast<const Run* 158 const Run* localRun = static_cast<const Run*>(run); 161 << 159 162 // primary particle info << 160 //primary particle info 163 // 161 // 164 fParticle = localRun->fParticle; 162 fParticle = localRun->fParticle; 165 fEkin = localRun->fEkin; << 163 fEkin = localRun->fEkin; 166 << 164 167 // accumulate sums 165 // accumulate sums 168 // 166 // 169 fTotalCount += localRun->fTotalCount; << 167 fTotalCount += localRun->fTotalCount; 170 fGammaCount += localRun->fGammaCount; << 168 fGammaCount += localRun->fGammaCount; 171 fSumTrack += localRun->fSumTrack; 169 fSumTrack += localRun->fSumTrack; 172 fSumTrack2 += localRun->fSumTrack2; << 170 fSumTrack2 += localRun->fSumTrack2; >> 171 >> 172 >> 173 //map: processes count >> 174 std::map<G4String,G4int>::const_iterator itp; >> 175 for ( itp = localRun->fProcCounter.begin(); >> 176 itp != localRun->fProcCounter.end(); ++itp ) { 173 177 174 fPbalance[0] += localRun->fPbalance[0]; << 175 G4double min, max; << 176 min = localRun->fPbalance[1]; << 177 max = localRun->fPbalance[2]; << 178 if (fPbalance[1] > min) fPbalance[1] = min; << 179 if (fPbalance[2] < max) fPbalance[2] = max; << 180 << 181 fNbGamma[0] += localRun->fNbGamma[0]; << 182 G4int nbmin, nbmax; << 183 nbmin = localRun->fNbGamma[1]; << 184 nbmax = localRun->fNbGamma[2]; << 185 if (fNbGamma[1] > nbmin) fNbGamma[1] = nbmin << 186 if (fNbGamma[2] < nbmax) fNbGamma[2] = nbmax << 187 << 188 // map: processes count << 189 std::map<G4String, G4int>::const_iterator it << 190 for (itp = localRun->fProcCounter.begin(); i << 191 G4String procName = itp->first; 178 G4String procName = itp->first; 192 G4int localCount = itp->second; 179 G4int localCount = itp->second; 193 if (fProcCounter.find(procName) == fProcCo << 180 if ( fProcCounter.find(procName) == fProcCounter.end()) { 194 fProcCounter[procName] = localCount; 181 fProcCounter[procName] = localCount; 195 } 182 } 196 else { 183 else { 197 fProcCounter[procName] += localCount; 184 fProcCounter[procName] += localCount; 198 } << 185 } 199 } 186 } 200 << 187 201 // map: nuclear channels << 188 //map: nuclear channels 202 std::map<G4String, NuclChannel>::const_itera << 189 std::map<G4String,NuclChannel>::const_iterator itc; 203 for (itc = localRun->fNuclChannelMap.begin() << 190 for (itc = localRun->fNuclChannelMap.begin(); >> 191 itc != localRun->fNuclChannelMap.end(); ++itc) { >> 192 204 G4String name = itc->first; 193 G4String name = itc->first; 205 const NuclChannel& localData = itc->second << 194 const NuclChannel& localData = itc->second; 206 if (fNuclChannelMap.find(name) == fNuclCha << 195 if ( fNuclChannelMap.find(name) == fNuclChannelMap.end()) { 207 fNuclChannelMap[name] = NuclChannel(loca << 196 fNuclChannelMap[name] >> 197 = NuclChannel(localData.fCount, localData.fQ); 208 } 198 } 209 else { 199 else { 210 NuclChannel& data = fNuclChannelMap[name << 200 NuclChannel& data = fNuclChannelMap[name]; 211 data.fCount += localData.fCount; 201 data.fCount += localData.fCount; 212 data.fQ += localData.fQ; << 202 data.fQ += localData.fQ; 213 } << 203 } 214 } << 204 } 215 << 205 216 // map: particles count << 206 //map: particles count 217 std::map<G4String, ParticleData>::const_iter << 207 std::map<G4String,ParticleData>::const_iterator itn; 218 for (itn = localRun->fParticleDataMap.begin( << 208 for (itn = localRun->fParticleDataMap.begin(); >> 209 itn != localRun->fParticleDataMap.end(); ++itn) { >> 210 219 G4String name = itn->first; 211 G4String name = itn->first; 220 const ParticleData& localData = itn->secon << 212 const ParticleData& localData = itn->second; 221 if (fParticleDataMap.find(name) == fPartic << 213 if ( fParticleDataMap.find(name) == fParticleDataMap.end()) { 222 fParticleDataMap[name] = << 214 fParticleDataMap[name] 223 ParticleData(localData.fCount, localDa << 215 = ParticleData(localData.fCount, >> 216 localData.fEmean, >> 217 localData.fEmin, >> 218 localData.fEmax); 224 } 219 } 225 else { 220 else { 226 ParticleData& data = fParticleDataMap[na << 221 ParticleData& data = fParticleDataMap[name]; 227 data.fCount += localData.fCount; 222 data.fCount += localData.fCount; 228 data.fEmean += localData.fEmean; 223 data.fEmean += localData.fEmean; 229 G4double emin = localData.fEmin; 224 G4double emin = localData.fEmin; 230 if (emin < data.fEmin) data.fEmin = emin 225 if (emin < data.fEmin) data.fEmin = emin; 231 G4double emax = localData.fEmax; 226 G4double emax = localData.fEmax; 232 if (emax > data.fEmax) data.fEmax = emax << 227 if (emax > data.fEmax) data.fEmax = emax; 233 } << 228 } 234 } 229 } 235 << 230 236 G4Run::Merge(run); << 231 G4Run::Merge(run); 237 } << 232 } 238 233 239 //....oooOO0OOooo........oooOO0OOooo........oo 234 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 240 235 241 void Run::EndOfRun(G4bool print) << 236 void Run::EndOfRun(G4bool print) 242 { 237 { 243 G4int prec = 5, wid = prec + 2; << 238 G4int prec = 5, wid = prec + 2; 244 G4int dfprec = G4cout.precision(prec); 239 G4int dfprec = G4cout.precision(prec); 245 << 240 246 // run condition << 241 //run condition 247 // 242 // 248 const G4Material* material = fDetector->GetM << 243 G4Material* material = fDetector->GetMaterial(); 249 G4double density = material->GetDensity(); 244 G4double density = material->GetDensity(); 250 << 245 251 G4String Particle = fParticle->GetParticleNa << 246 G4String Particle = fParticle->GetParticleName(); 252 G4cout << "\n The run is " << numberOfEvent << 247 G4cout << "\n The run is " << numberOfEvent << " "<< Particle << " of " 253 << G4BestUnit(fEkin, "Energy") << " t << 248 << G4BestUnit(fEkin,"Energy") << " through " 254 << " of " << material->GetName() << " << 249 << G4BestUnit(fDetector->GetSize(),"Length") << " of " 255 << ")" << G4endl; << 250 << material->GetName() << " (density: " 256 << 251 << G4BestUnit(density,"Volumic Mass") << ")" << G4endl; 257 if (numberOfEvent == 0) { << 252 258 G4cout.precision(dfprec); << 253 if (numberOfEvent == 0) { G4cout.precision(dfprec); return;} 259 return; << 254 260 } << 255 //frequency of processes 261 << 262 // frequency of processes << 263 // 256 // 264 G4cout << "\n Process calls frequency:" << G << 257 G4cout << "\n Process calls frequency:" << G4endl; 265 G4int survive = 0; 258 G4int survive = 0; 266 std::map<G4String, G4int>::iterator it; << 259 std::map<G4String,G4int>::iterator it; 267 for (it = fProcCounter.begin(); it != fProcC 260 for (it = fProcCounter.begin(); it != fProcCounter.end(); it++) { 268 G4String procName = it->first; << 261 G4String procName = it->first; 269 G4int count = it->second; << 262 G4int count = it->second; 270 G4cout << "\t" << procName << "= " << coun << 263 G4cout << "\t" << procName << "= " << count; 271 if (procName == "Transportation") survive << 264 if (procName == "Transportation") survive = count; 272 } 265 } 273 G4cout << G4endl; 266 G4cout << G4endl; 274 << 267 275 if (survive > 0) { 268 if (survive > 0) { 276 G4cout << "\n Nb of incident particles sur 269 G4cout << "\n Nb of incident particles surviving after " 277 << G4BestUnit(fDetector->GetSize(), << 270 << G4BestUnit(fDetector->GetSize(),"Length") << " of " 278 << survive << G4endl; << 271 << material->GetName() << " : " << survive << G4endl; 279 } 272 } 280 << 273 281 if (fTotalCount == 0) fTotalCount = 1; // f << 274 if (fTotalCount == 0) fTotalCount = 1; //force printing anyway 282 << 275 283 // compute mean free path and related quanti << 276 //compute mean free path and related quantities 284 // << 277 // 285 G4double MeanFreePath = fSumTrack / fTotalCo << 278 G4double MeanFreePath = fSumTrack /fTotalCount; 286 G4double MeanTrack2 = fSumTrack2 / fTotalCou << 279 G4double MeanTrack2 = fSumTrack2/fTotalCount; 287 G4double rms = std::sqrt(std::fabs(MeanTrack << 280 G4double rms = std::sqrt(std::fabs(MeanTrack2 - MeanFreePath*MeanFreePath)); 288 G4double CrossSection = 0.0; 281 G4double CrossSection = 0.0; 289 if (MeanFreePath > 0.0) { << 282 if(MeanFreePath > 0.0) { CrossSection = 1./MeanFreePath; } 290 CrossSection = 1. / MeanFreePath; << 283 G4double massicMFP = MeanFreePath*density; 291 } << 284 G4double massicCS = 0.0; 292 G4double massicMFP = MeanFreePath * density; << 285 if(massicMFP > 0.0) { massicCS = 1./massicMFP; } 293 G4double massicCS = 0.0; << 286 294 if (massicMFP > 0.0) { << 287 G4cout << "\n\n MeanFreePath:\t" << G4BestUnit(MeanFreePath,"Length") 295 massicCS = 1. / massicMFP; << 288 << " +- " << G4BestUnit( rms,"Length") 296 } << 289 << "\tmassic: " << G4BestUnit(massicMFP, "Mass/Surface") 297 << 290 << "\n CrossSection:\t" << CrossSection*cm << " cm^-1 " 298 G4cout << "\n\n MeanFreePath:\t" << G4BestUn << 291 << "\t\tmassic: " << G4BestUnit(massicCS, "Surface/Mass") 299 << G4BestUnit(rms, "Length") << "\tma << 292 << G4endl; 300 << "\n CrossSection:\t" << CrossSecti << 293 301 << "\t\tmassic: " << G4BestUnit(massi << 294 //cross section per atom (only for single material) 302 << 303 // cross section per atom (only for single m << 304 // 295 // 305 if (material->GetNumberOfElements() == 1) { 296 if (material->GetNumberOfElements() == 1) { 306 G4double nbAtoms = material->GetTotNbOfAto 297 G4double nbAtoms = material->GetTotNbOfAtomsPerVolume(); 307 G4double crossSection = CrossSection / nbA << 298 G4double crossSection = CrossSection/nbAtoms; 308 G4cout << " crossSection per atom:\t" << G << 299 G4cout << " crossSection per atom:\t" 309 } << 300 << G4BestUnit(crossSection,"Surface") << G4endl; 310 // check cross section from G4HadronicProces << 301 } >> 302 //check cross section from G4HadronicProcessStore 311 // 303 // 312 G4cout << "\n Verification: " 304 G4cout << "\n Verification: " 313 << "crossSections from G4HadronicProc << 305 << "crossSections from G4HadronicProcessStore:"; 314 << 306 315 G4ProcessTable* processTable = G4ProcessTabl << 307 G4ProcessTable* processTable = G4ProcessTable::GetProcessTable(); 316 G4HadronicProcessStore* store = G4HadronicPr 308 G4HadronicProcessStore* store = G4HadronicProcessStore::Instance(); 317 G4double sumc1 = 0.0, sumc2 = 0.0; << 309 G4double sumc1 = 0.0, sumc2 = 0.0; 318 const G4Element* element = << 310 if (material->GetNumberOfElements() == 1) { 319 (material->GetNumberOfElements() == 1) ? m << 311 const G4Element* element = material->GetElement(0); 320 for (it = fProcCounter.begin(); it != fProcC << 312 for (it = fProcCounter.begin(); it != fProcCounter.end(); it++) { 321 G4String procName = it->first; << 313 G4String procName = it->first; 322 const G4VProcess* process = processTable-> << 314 G4VProcess* process = processTable->FindProcess(procName, fParticle); 323 PrintXS(process, material, element, store, << 315 G4double xs1 = 324 } << 316 store->GetCrossSectionPerVolume(fParticle,fEkin,process,material); 325 if (sumc1 > 0.0) { << 317 G4double massSigma = xs1/density; 326 G4cout << "\n" << 318 sumc1 += massSigma; 327 << std::setw(20) << "total" << 319 G4double xs2 = 328 << " = " << G4BestUnit(sumc1, "Surf << 320 store->GetCrossSectionPerAtom(fParticle,fEkin,process,element,material); 329 if (sumc2 > 0.0) { << 321 sumc2 += xs2; 330 G4cout << G4BestUnit(sumc2, "Surface"); << 322 G4cout << "\n" << std::setw(20) << procName << "= " 331 } << 323 << G4BestUnit(massSigma, "Surface/Mass") << "\t" 332 G4cout << G4endl; << 324 << G4BestUnit(xs2, "Surface"); 333 } << 325 334 else { << 326 } 335 G4cout << " not available" << G4endl; << 327 G4cout << "\n" << std::setw(20) << "total" << "= " 336 } << 328 << G4BestUnit(sumc1, "Surface/Mass") << "\t" 337 << 329 << G4BestUnit(sumc2, "Surface") << G4endl; 338 // nuclear channel count << 330 } else { 339 // << 331 for (it = fProcCounter.begin(); it != fProcCounter.end(); it++) { 340 G4cout << "\n List of nuclear reactions: \n" << 332 G4String procName = it->first; 341 std::map<G4String, NuclChannel>::iterator ic << 333 G4VProcess* process = processTable->FindProcess(procName, fParticle); 342 for (ic = fNuclChannelMap.begin(); ic != fNu << 334 G4double xs = 343 G4String name = ic->first; << 335 store->GetCrossSectionPerVolume(fParticle,fEkin,process,material); >> 336 G4double massSigma = xs/density; >> 337 sumc1 += massSigma; >> 338 G4cout << "\n" << std::setw(20) << procName << "= " >> 339 << G4BestUnit(massSigma, "Surface/Mass"); >> 340 } >> 341 G4cout << "\n" << std::setw(20) << "total" << "= " >> 342 << G4BestUnit(sumc1, "Surface/Mass") << G4endl; >> 343 } >> 344 >> 345 //nuclear channel count >> 346 // >> 347 G4cout << "\n List of nuclear reactions: \n" << G4endl; >> 348 std::map<G4String,NuclChannel>::iterator ic; >> 349 for (ic = fNuclChannelMap.begin(); ic != fNuclChannelMap.end(); ic++) { >> 350 G4String name = ic->first; 344 NuclChannel data = ic->second; 351 NuclChannel data = ic->second; 345 G4int count = data.fCount; 352 G4int count = data.fCount; 346 G4double Q = data.fQ / count; << 353 G4double Q = data.fQ/count; 347 if (print) << 354 if (print) 348 G4cout << " " << std::setw(60) << name << 355 G4cout << " " << std::setw(50) << name << ": " << std::setw(7) << count 349 << " Q = " << std::setw(wid) << << 356 << " Q = " << std::setw(wid) << G4BestUnit(Q, "Energy") 350 } << 357 << G4endl; 351 << 358 } 352 // Gamma count << 359 353 // << 360 //Gamma count 354 if (print && (fGammaCount > 0)) { << 361 // 355 G4cout << "\n" << 362 if (print && (fGammaCount > 0)) { 356 << std::setw(58) << "number of gamm << 363 G4cout << "\n" << std::setw(58) << "Number of gamma: N = " 357 << fNbGamma[2] << G4endl; << 364 << fNbGamma[1] << " --> " << fNbGamma[2] << G4endl; 358 } << 365 } 359 << 366 360 if (print && fTargetXXX) { << 367 if (print && fTargetXXX) { 361 G4cout << "\n --> NOTE: XXXX because neu << 368 G4cout 362 } << 369 << "\n --> NOTE: XXXX because neutronHP is unable to return target nucleus" 363 << 370 << G4endl; 364 // particles count << 371 } 365 // << 372 366 G4cout << "\n List of generated particles:" << 373 //particles count 367 << 374 // 368 std::map<G4String, ParticleData>::iterator i << 375 G4cout << "\n List of generated particles:" << G4endl; 369 for (itn = fParticleDataMap.begin(); itn != << 376 >> 377 std::map<G4String,ParticleData>::iterator itn; >> 378 for (itn = fParticleDataMap.begin(); itn != fParticleDataMap.end(); itn++) { 370 G4String name = itn->first; 379 G4String name = itn->first; 371 ParticleData data = itn->second; 380 ParticleData data = itn->second; 372 G4int count = data.fCount; 381 G4int count = data.fCount; 373 G4double eMean = data.fEmean / count; << 382 G4double eMean = data.fEmean/count; 374 G4double eMin = data.fEmin; 383 G4double eMin = data.fEmin; 375 G4double eMax = data.fEmax; << 384 G4double eMax = data.fEmax; 376 if (print) << 385 if (print) 377 G4cout << " " << std::setw(13) << name << 386 G4cout << " " << std::setw(13) << name << ": " << std::setw(7) << count 378 << " Emean = " << std::setw(wid) << 387 << " Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy") 379 << G4BestUnit(eMin, "Energy") << << 388 << "\t( " << G4BestUnit(eMin, "Energy") 380 << G4endl; << 389 << " --> " << G4BestUnit(eMax, "Energy") 381 } << 390 << ")" << G4endl; 382 << 391 } 383 // energy momentum balance << 392 384 // << 393 //energy momentum balance 385 if (fTotalCount > 1) { << 394 // 386 G4double Pbmean = fPbalance[0] / fTotalCou << 395 if (fTotalCount > 1) { 387 G4cout << "\n Momentum balance: Pmean = << 396 G4double Pbmean = fPbalance[0]/fTotalCount; 388 << "\t( " << G4BestUnit(fPbalance[1 << 397 G4cout << "\n Momentum balance: Pmean = " 389 << G4BestUnit(fPbalance[2], "Energy << 398 << std::setw(wid) << G4BestUnit(Pbmean, "Energy") 390 << G4endl; << 399 << "\t( " << G4BestUnit(fPbalance[1], "Energy") 391 } << 400 << " --> " << G4BestUnit(fPbalance[2], "Energy") 392 << 401 << ") \n" << G4endl; 393 // normalize histograms << 402 } >> 403 >> 404 //normalize histograms 394 ////G4AnalysisManager* analysisManager = G4A 405 ////G4AnalysisManager* analysisManager = G4AnalysisManager::Instance(); 395 ////G4double factor = 1./numberOfEvent; 406 ////G4double factor = 1./numberOfEvent; 396 ////analysisManager->ScaleH1(3,factor); 407 ////analysisManager->ScaleH1(3,factor); 397 << 408 398 // remove all contents in fProcCounter, fCou << 409 //remove all contents in fProcCounter, fCount 399 fProcCounter.clear(); 410 fProcCounter.clear(); 400 fNuclChannelMap.clear(); << 411 fNuclChannelMap.clear(); 401 fParticleDataMap.clear(); 412 fParticleDataMap.clear(); 402 << 413 403 // restore default format << 414 //restore default format 404 G4cout.precision(dfprec); << 415 G4cout.precision(dfprec); 405 } << 406 << 407 //....oooOO0OOooo........oooOO0OOooo........oo << 408 << 409 void Run::PrintXS(const G4VProcess* proc, cons << 410 G4HadronicProcessStore* stor << 411 { << 412 if (nullptr == proc) { << 413 return; << 414 } << 415 G4double xs1 = store->GetCrossSectionPerVolu << 416 G4double massSigma = xs1 / density; << 417 sum1 += massSigma; << 418 if (nullptr != elm) { << 419 G4double xs2 = store->GetCrossSectionPerAt << 420 sum2 += xs2; << 421 G4cout << "\n" << 422 << std::setw(20) << proc->GetProces << 423 << G4BestUnit(massSigma, "Surface/M << 424 } << 425 else { << 426 G4cout << "\n" << 427 << std::setw(20) << proc->GetProces << 428 << G4BestUnit(massSigma, "Surface/M << 429 } << 430 } 416 } 431 417 432 //....oooOO0OOooo........oooOO0OOooo........oo 418 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 433 419