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/TestEm18/src/RunActi 26 /// \file electromagnetic/TestEm18/src/RunAction.cc 27 /// \brief Implementation of the RunAction cla 27 /// \brief Implementation of the RunAction 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 "RunAction.hh" 33 #include "RunAction.hh" 34 << 35 #include "DetectorConstruction.hh" 34 #include "DetectorConstruction.hh" 36 #include "HistoManager.hh" << 37 #include "PrimaryGeneratorAction.hh" 35 #include "PrimaryGeneratorAction.hh" >> 36 #include "HistoManager.hh" 38 37 39 #include "G4EmCalculator.hh" << 40 #include "G4Run.hh" 38 #include "G4Run.hh" 41 #include "G4UnitsTable.hh" 39 #include "G4UnitsTable.hh" 42 #include "Randomize.hh" << 40 #include "G4EmCalculator.hh" 43 41 >> 42 #include "Randomize.hh" 44 #include <iomanip> 43 #include <iomanip> 45 44 46 //....oooOO0OOooo........oooOO0OOooo........oo 45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 47 46 48 RunAction::RunAction(DetectorConstruction* det 47 RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* kin) 49 : fDetector(det), fPrimary(kin) << 48 :G4UserRunAction(),fDetector(det), fPrimary(kin), fHistoManager(0) 50 { << 49 { 51 fHistoManager = new HistoManager(); << 50 fHistoManager = new HistoManager(); 52 } 51 } 53 52 54 //....oooOO0OOooo........oooOO0OOooo........oo 53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 55 54 56 RunAction::~RunAction() 55 RunAction::~RunAction() 57 { << 56 { 58 delete fHistoManager; << 57 delete fHistoManager; 59 } 58 } 60 59 61 //....oooOO0OOooo........oooOO0OOooo........oo 60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 62 61 63 void RunAction::BeginOfRunAction(const G4Run*) 62 void RunAction::BeginOfRunAction(const G4Run*) 64 { 63 { 65 // initialisation << 64 //initialisation 66 // 65 // 67 fNbSteps = 0; 66 fNbSteps = 0; 68 fTrackLength = 0.; 67 fTrackLength = 0.; 69 fStepMin = DBL_MAX; 68 fStepMin = DBL_MAX; 70 fStepMax = 0.; 69 fStepMax = 0.; 71 70 72 fEdepPrimary = fEdepSecondary = fEdepTotal = 71 fEdepPrimary = fEdepSecondary = fEdepTotal = 0.; 73 fEdepPrimMin = fEdepSecMin = fEdepTotMin = D 72 fEdepPrimMin = fEdepSecMin = fEdepTotMin = DBL_MAX; 74 fEdepPrimMax = fEdepSecMax = fEdepTotMax = 0 73 fEdepPrimMax = fEdepSecMax = fEdepTotMax = 0.; 75 74 76 fEnergyTransfered = 0.; 75 fEnergyTransfered = 0.; 77 fEtransfMin = DBL_MAX; 76 fEtransfMin = DBL_MAX; 78 fEtransfMax = 0.; 77 fEtransfMax = 0.; 79 78 80 fEnergyLost = 0.; 79 fEnergyLost = 0.; 81 fElostMin = DBL_MAX; 80 fElostMin = DBL_MAX; 82 fElostMax = 0.; 81 fElostMax = 0.; 83 82 84 fEnergyBalance = 0.; 83 fEnergyBalance = 0.; 85 fEbalMin = DBL_MAX; 84 fEbalMin = DBL_MAX; 86 fEbalMax = 0.; 85 fEbalMax = 0.; 87 86 88 // histograms << 87 //histograms 89 // 88 // 90 G4AnalysisManager* analysisManager = G4Analy 89 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance(); 91 if (analysisManager->IsActive()) { << 90 if ( analysisManager->IsActive() ) { 92 analysisManager->OpenFile(); 91 analysisManager->OpenFile(); 93 } 92 } 94 93 95 // show Rndm status 94 // show Rndm status 96 CLHEP::HepRandom::showEngineStatus(); 95 CLHEP::HepRandom::showEngineStatus(); 97 } 96 } 98 97 99 //....oooOO0OOooo........oooOO0OOooo........oo 98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 100 99 101 void RunAction::CountProcesses(G4String procNa << 100 void RunAction::CountProcesses(G4String procName) 102 { 101 { 103 std::map<G4String, G4int>::iterator it = fPr << 102 std::map<G4String,G4int>::iterator it = fProcCounter.find(procName); 104 if (it == fProcCounter.end()) { << 103 if ( it == fProcCounter.end()) { 105 fProcCounter[procName] = 1; 104 fProcCounter[procName] = 1; 106 } 105 } 107 else { 106 else { 108 fProcCounter[procName]++; << 107 fProcCounter[procName]++; 109 } 108 } 110 } 109 } 111 110 112 //....oooOO0OOooo........oooOO0OOooo........oo 111 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 113 112 114 void RunAction::TrackLength(G4double step) << 113 void RunAction::TrackLength (G4double step) 115 { 114 { 116 fTrackLength += step; << 115 fTrackLength += step; fNbSteps++; 117 fNbSteps++; << 116 if (step<fStepMin) fStepMin = step; 118 if (step < fStepMin) fStepMin = step; << 117 if (step>fStepMax) fStepMax = step; 119 if (step > fStepMax) fStepMax = step; << 120 } 118 } 121 119 122 //....oooOO0OOooo........oooOO0OOooo........oo 120 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 123 121 124 void RunAction::EnergyDeposited(G4double edepP << 122 void RunAction::EnergyDeposited (G4double edepPrim, G4double edepSecond) 125 { 123 { 126 fEdepPrimary += edepPrim; 124 fEdepPrimary += edepPrim; 127 if (edepPrim < fEdepPrimMin) fEdepPrimMin = << 125 if (edepPrim<fEdepPrimMin) fEdepPrimMin = edepPrim; 128 if (edepPrim > fEdepPrimMax) fEdepPrimMax = << 126 if (edepPrim>fEdepPrimMax) fEdepPrimMax = edepPrim; 129 << 127 130 fEdepSecondary += edepSecond; 128 fEdepSecondary += edepSecond; 131 if (edepSecond < fEdepSecMin) fEdepSecMin = << 129 if (edepSecond<fEdepSecMin) fEdepSecMin = edepSecond; 132 if (edepSecond > fEdepSecMax) fEdepSecMax = << 130 if (edepSecond>fEdepSecMax) fEdepSecMax = edepSecond; 133 } 131 } 134 132 135 //....oooOO0OOooo........oooOO0OOooo........oo 133 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 136 134 137 void RunAction::EnergyTransferedByProcess(G4St 135 void RunAction::EnergyTransferedByProcess(G4String process, G4double energy) 138 { 136 { 139 std::map<G4String, MinMaxData>::iterator it 137 std::map<G4String, MinMaxData>::iterator it = fEtransfByProcess.find(process); 140 if (it == fEtransfByProcess.end()) { << 138 if ( it == fEtransfByProcess.end()) { 141 fEtransfByProcess[process] = MinMaxData(1, 139 fEtransfByProcess[process] = MinMaxData(1, energy, energy, energy); 142 } 140 } 143 else { 141 else { 144 MinMaxData& data = it->second; 142 MinMaxData& data = it->second; 145 data.fCount++; 143 data.fCount++; 146 data.fVsum += energy; 144 data.fVsum += energy; 147 // update min max << 145 //update min max 148 G4double emin = data.fVmin; 146 G4double emin = data.fVmin; 149 if (energy < emin) data.fVmin = energy; 147 if (energy < emin) data.fVmin = energy; 150 G4double emax = data.fVmax; 148 G4double emax = data.fVmax; 151 if (energy > emax) data.fVmax = energy; << 149 if (energy > emax) data.fVmax = energy; 152 } 150 } 153 } 151 } 154 152 155 //....oooOO0OOooo........oooOO0OOooo........oo 153 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 156 154 157 void RunAction::EnergyTransfered(G4double ener << 155 void RunAction::EnergyTransfered (G4double energy) 158 { 156 { 159 fEnergyTransfered += energy; 157 fEnergyTransfered += energy; 160 if (energy < fEtransfMin) fEtransfMin = ener << 158 if (energy<fEtransfMin) fEtransfMin = energy; 161 if (energy > fEtransfMax) fEtransfMax = ener << 159 if (energy>fEtransfMax) fEtransfMax = energy; 162 } 160 } 163 161 164 //....oooOO0OOooo........oooOO0OOooo........oo 162 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 165 163 166 void RunAction::TotalEnergyLost(G4double energ << 164 void RunAction::TotalEnergyLost (G4double energy) 167 { 165 { 168 fEnergyLost += energy; 166 fEnergyLost += energy; 169 if (energy < fElostMin) fElostMin = energy; << 167 if (energy<fElostMin) fElostMin = energy; 170 if (energy > fElostMax) fElostMax = energy; << 168 if (energy>fElostMax) fElostMax = energy; 171 } 169 } 172 170 173 //....oooOO0OOooo........oooOO0OOooo........oo 171 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 174 172 175 void RunAction::EnergyBalance(G4double energy) << 173 void RunAction::EnergyBalance (G4double energy) 176 { 174 { 177 fEnergyBalance += energy; 175 fEnergyBalance += energy; 178 if (energy < fEbalMin) fEbalMin = energy; << 176 if (energy<fEbalMin) fEbalMin = energy; 179 if (energy > fEbalMax) fEbalMax = energy; << 177 if (energy>fEbalMax) fEbalMax = energy; 180 } 178 } 181 179 182 //....oooOO0OOooo........oooOO0OOooo........oo 180 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 183 181 184 void RunAction::TotalEnergyDeposit(G4double en << 182 void RunAction::TotalEnergyDeposit (G4double energy) 185 { 183 { 186 fEdepTotal += energy; 184 fEdepTotal += energy; 187 if (energy < fEdepTotMin) fEdepTotMin = ener << 185 if (energy<fEdepTotMin) fEdepTotMin = energy; 188 if (energy > fEdepTotMax) fEdepTotMax = ener << 186 if (energy>fEdepTotMax) fEdepTotMax = energy; 189 } 187 } 190 188 191 //....oooOO0OOooo........oooOO0OOooo........oo 189 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 192 190 193 void RunAction::EnergySpectrumOfSecondaries(G4 191 void RunAction::EnergySpectrumOfSecondaries(G4String particle, G4double energy) 194 { 192 { 195 std::map<G4String, MinMaxData>::iterator it << 193 std::map<G4String,MinMaxData>::iterator it = fEkinOfSecondaries.find(particle); 196 if (it == fEkinOfSecondaries.end()) { << 194 if ( it == fEkinOfSecondaries.end()) { 197 fEkinOfSecondaries[particle] = MinMaxData( 195 fEkinOfSecondaries[particle] = MinMaxData(1, energy, energy, energy); 198 } 196 } 199 else { 197 else { 200 MinMaxData& data = it->second; 198 MinMaxData& data = it->second; 201 data.fCount++; 199 data.fCount++; 202 data.fVsum += energy; 200 data.fVsum += energy; 203 // update min max << 201 //update min max 204 G4double emin = data.fVmin; 202 G4double emin = data.fVmin; 205 if (energy < emin) data.fVmin = energy; 203 if (energy < emin) data.fVmin = energy; 206 G4double emax = data.fVmax; 204 G4double emax = data.fVmax; 207 if (energy > emax) data.fVmax = energy; << 205 if (energy > emax) data.fVmax = energy; 208 } 206 } 209 } 207 } 210 208 211 //....oooOO0OOooo........oooOO0OOooo........oo 209 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 212 210 213 void RunAction::EndOfRunAction(const G4Run* aR 211 void RunAction::EndOfRunAction(const G4Run* aRun) 214 { 212 { 215 G4int nbEvents = aRun->GetNumberOfEvent(); 213 G4int nbEvents = aRun->GetNumberOfEvent(); 216 if (nbEvents == 0) return; 214 if (nbEvents == 0) return; 217 << 215 218 G4Material* material = fDetector->GetMateria 216 G4Material* material = fDetector->GetMaterial(); 219 G4double length = fDetector->GetSize(); << 217 G4double length = fDetector->GetSize(); 220 G4double density = material->GetDensity(); 218 G4double density = material->GetDensity(); 221 << 219 222 G4ParticleDefinition* particle = fPrimary->G << 220 G4ParticleDefinition* particle = fPrimary->GetParticleGun() >> 221 ->GetParticleDefinition(); 223 G4String partName = particle->GetParticleNam 222 G4String partName = particle->GetParticleName(); 224 G4double ePrimary = fPrimary->GetParticleGun 223 G4double ePrimary = fPrimary->GetParticleGun()->GetParticleEnergy(); 225 << 224 226 G4int prec = G4cout.precision(3); 225 G4int prec = G4cout.precision(3); 227 G4cout << "\n ======================== run s 226 G4cout << "\n ======================== run summary ======================\n"; 228 G4cout << "\n The run was " << nbEvents << " 227 G4cout << "\n The run was " << nbEvents << " " << partName << " of " 229 << G4BestUnit(ePrimary, "Energy") << << 228 << G4BestUnit(ePrimary,"Energy") << " through " 230 << material->GetName() << " (density: << 229 << G4BestUnit(length,"Length") << " of " >> 230 << material->GetName() << " (density: " >> 231 << G4BestUnit(density,"Volumic Mass") << ")"; 231 G4cout << G4endl; 232 G4cout << G4endl; 232 233 233 if (particle->GetPDGCharge() == 0.) return; 234 if (particle->GetPDGCharge() == 0.) return; 234 235 235 G4cout.precision(4); 236 G4cout.precision(4); 236 237 237 // frequency of processes << 238 //frequency of processes 238 // 239 // 239 G4cout << "\n Process defining step :" << G4 240 G4cout << "\n Process defining step :" << G4endl; 240 G4int index = 0; 241 G4int index = 0; 241 for (const auto& procCounter : fProcCounter) << 242 for ( const auto& procCounter : fProcCounter ) { 242 G4String procName = procCounter.first; << 243 G4String procName = procCounter.first; 243 G4int count = procCounter.second; << 244 G4int count = procCounter.second; 244 G4String space = " "; << 245 G4String space = " "; if (++index%4 == 0) space = "\n"; 245 if (++index % 4 == 0) space = "\n"; << 246 G4cout << " " << std::setw(15) << procName << "="<< std::setw(7) << count 246 G4cout << " " << std::setw(15) << procName << 247 << space; 247 } 248 } 248 G4cout << G4endl; 249 G4cout << G4endl; 249 250 250 // track length << 251 //track length 251 // << 252 G4double trackLPerEvent = fTrackLength / nbE << 253 G4double nbStepPerEvent = double(fNbSteps) / << 254 G4double stepSize = fTrackLength / fNbSteps; << 255 << 256 G4cout << "\n TrackLength = " << G4BestUnit( << 257 << " nb of steps = " << nbStepPerEve << 258 << " stepSize = " << G4BestUnit(step << 259 << G4BestUnit(fStepMin, "Length") << << 260 << G4endl; << 261 << 262 // continuous energy deposited by primary tr << 263 // 252 // 264 G4double energyPerEvent = fEdepPrimary / nbE << 253 G4double trackLPerEvent = fTrackLength/nbEvents; 265 << 254 G4double nbStepPerEvent = double(fNbSteps)/nbEvents; 266 G4cout << "\n Energy continuously deposited << 255 G4double stepSize = fTrackLength/fNbSteps; 267 << " (restricted dE/dx) dE1 = " << G << 256 268 << G4BestUnit(fEdepPrimMin, "Energy") << 257 G4cout 269 << ")" << G4endl; << 258 << "\n TrackLength = " 270 << 259 << G4BestUnit(trackLPerEvent, "Length") 271 // eveluation of dE1 from reading restricted << 260 << " nb of steps = " << nbStepPerEvent >> 261 << " stepSize = " << G4BestUnit(stepSize, "Length") >> 262 << " (" << G4BestUnit(fStepMin, "Length") >> 263 << "--> " << G4BestUnit(fStepMax, "Length") << ")" >> 264 << G4endl; >> 265 >> 266 //continuous energy deposited by primary track dE1 >> 267 // >> 268 G4double energyPerEvent = fEdepPrimary/nbEvents; >> 269 >> 270 G4cout >> 271 << "\n Energy continuously deposited along primary track" >> 272 << " (restricted dE/dx) dE1 = " >> 273 << G4BestUnit(energyPerEvent, "Energy") >> 274 << " (" << G4BestUnit(fEdepPrimMin, "Energy") >> 275 << " --> " << G4BestUnit(fEdepPrimMax, "Energy") << ")" >> 276 << G4endl; >> 277 >> 278 //eveluation of dE1 from reading restricted Range table 272 // 279 // 273 G4EmCalculator emCal; 280 G4EmCalculator emCal; 274 << 281 275 G4double r0 = emCal.GetRangeFromRestricteDED << 282 G4double r0 = emCal.GetRangeFromRestricteDEDX(ePrimary,particle,material); 276 G4double r1 = r0 - trackLPerEvent; 283 G4double r1 = r0 - trackLPerEvent; 277 G4double etry = ePrimary - energyPerEvent; << 284 G4double etry = ePrimary - energyPerEvent; 278 G4double efinal = 0.; 285 G4double efinal = 0.; 279 if (r1 > 0.) efinal = GetEnergyFromRestricte << 286 if (r1 > 0.) efinal = GetEnergyFromRestrictedRange(r1,particle,material,etry); 280 G4double dEtable = ePrimary - efinal; 287 G4double dEtable = ePrimary - efinal; 281 G4double ratio = 0.; 288 G4double ratio = 0.; 282 if (dEtable > 0.) ratio = energyPerEvent / d << 289 if (dEtable > 0.) ratio = energyPerEvent/dEtable; 283 << 290 284 G4cout << "\n Evaluation of dE1 from reading << 291 G4cout 285 << G4BestUnit(dEtable, "Energy") << " << 292 << "\n Evaluation of dE1 from reading restricted Range table : dE1_table = " 286 << 293 << G4BestUnit(dEtable, "Energy") >> 294 << " ---> dE1/dE1_table = " << ratio >> 295 << G4endl; >> 296 287 // energy transfered to secondary particles 297 // energy transfered to secondary particles by process : dE2 288 // 298 // 289 G4cout << "\n Energy transfered to secondary 299 G4cout << "\n Energy transfered to secondary particles :" << G4endl; 290 std::map<G4String, MinMaxData>::iterator it1 << 300 std::map<G4String,MinMaxData>::iterator it1; 291 for (it1 = fEtransfByProcess.begin(); it1 != 301 for (it1 = fEtransfByProcess.begin(); it1 != fEtransfByProcess.end(); it1++) { 292 G4String name = it1->first; << 302 G4String name = it1->first; 293 MinMaxData data = it1->second; << 303 MinMaxData data = it1->second; 294 energyPerEvent = data.fVsum / nbEvents; << 304 energyPerEvent = data.fVsum/nbEvents; 295 G4double eMin = data.fVmin; << 305 G4double eMin = data.fVmin; 296 G4double eMax = data.fVmax; << 306 G4double eMax = data.fVmax; 297 << 307 298 G4cout << " " << std::setw(17) << "due to << 308 G4cout << " " << std::setw(17) << "due to " + name << ": dE2 = " 299 << G4BestUnit(energyPerEvent, "Ener << 309 << std::setw(6) << G4BestUnit(energyPerEvent, "Energy") 300 << G4BestUnit(eMax, "Energy") << ") << 310 << " (" << G4BestUnit(eMin, "Energy") >> 311 << " --> " << G4BestUnit(eMax, "Energy") >> 312 << ")" << G4endl; 301 } 313 } 302 314 303 // total energy tranfered : dE3 = sum of dE2 315 // total energy tranfered : dE3 = sum of dE2 304 // 316 // 305 energyPerEvent = fEnergyTransfered / nbEvent << 317 energyPerEvent = fEnergyTransfered/nbEvents; 306 << 318 307 G4cout << "\n Total energy transfered to sec << 319 G4cout 308 << G4BestUnit(energyPerEvent, "Energy << 320 << "\n Total energy transfered to secondaries : dE3 = sum of dE2 = " 309 << " --> " << G4BestUnit(fEtransfMax, << 321 << G4BestUnit(energyPerEvent, "Energy") >> 322 << " (" << G4BestUnit(fEtransfMin, "Energy") >> 323 << " --> " << G4BestUnit(fEtransfMax, "Energy") << ")" >> 324 << G4endl; 310 325 311 // total energy lost by incident particle : 326 // total energy lost by incident particle : dE4 = dE1 + dE3 312 // 327 // 313 energyPerEvent = fEnergyLost / nbEvents; << 328 energyPerEvent = fEnergyLost/nbEvents; 314 << 329 315 G4cout << "\n Total energy lost by incident << 330 G4cout 316 << G4BestUnit(energyPerEvent, "Energy << 331 << "\n Total energy lost by incident particle : dE4 = dE1 + dE3 = " 317 << " --> " << G4BestUnit(fElostMax, " << 332 << G4BestUnit(energyPerEvent, "Energy") 318 << 333 << " (" << G4BestUnit(fElostMin, "Energy") >> 334 << " --> " << G4BestUnit(fElostMax, "Energy") << ")" >> 335 << G4endl; >> 336 319 // calcul of energy lost from energy balance 337 // calcul of energy lost from energy balance : dE4_bal = E_in - E_out 320 // 338 // 321 energyPerEvent = fEnergyBalance / nbEvents; << 339 energyPerEvent = fEnergyBalance/nbEvents; 322 << 340 323 G4cout << "\n calcul of dE4 from energy bala << 341 G4cout 324 << G4BestUnit(energyPerEvent, "Energy << 342 << "\n calcul of dE4 from energy balance : dE4_bal = E_in - E_out = " 325 << " --> " << G4BestUnit(fEbalMax, "E << 343 << G4BestUnit(energyPerEvent, "Energy") 326 << 344 << " (" << G4BestUnit(fEbalMin, "Energy") 327 // eveluation of dE4 from reading full Range << 345 << " --> " << G4BestUnit(fEbalMax, "Energy") << ")" >> 346 << G4endl; >> 347 >> 348 //eveluation of dE4 from reading full Range table 328 // 349 // 329 r0 = emCal.GetCSDARange(ePrimary, particle, << 350 r0 = emCal.GetCSDARange(ePrimary,particle,material); 330 r1 = r0 - trackLPerEvent; 351 r1 = r0 - trackLPerEvent; 331 etry = ePrimary - energyPerEvent; << 352 etry = ePrimary - energyPerEvent; 332 efinal = 0.; 353 efinal = 0.; 333 if (r1 > 0.) efinal = GetEnergyFromCSDARange << 354 if (r1 > 0.) efinal = GetEnergyFromCSDARange(r1,particle,material,etry); 334 dEtable = ePrimary - efinal; 355 dEtable = ePrimary - efinal; 335 ratio = 0.; 356 ratio = 0.; 336 if (dEtable > 0.) ratio = energyPerEvent / d << 357 if (dEtable > 0.) ratio = energyPerEvent/dEtable; 337 << 358 338 G4cout << "\n Evaluation of dE4 from reading << 359 G4cout 339 << G4BestUnit(dEtable, "Energy") << " << 360 << "\n Evaluation of dE4 from reading full Range table : dE4_table = " 340 << 361 << G4BestUnit(dEtable, "Energy") 341 // energy spectrum of secondary particles << 362 << " ---> dE4/dE4_table = " << ratio >> 363 << G4endl; >> 364 >> 365 //energy spectrum of secondary particles 342 // 366 // 343 G4cout << "\n Energy spectrum of secondary p 367 G4cout << "\n Energy spectrum of secondary particles :" << G4endl; 344 std::map<G4String, MinMaxData>::iterator it2 << 368 std::map<G4String,MinMaxData>::iterator it2; 345 for (it2 = fEkinOfSecondaries.begin(); it2 ! << 369 for (it2 = fEkinOfSecondaries.begin();it2 != fEkinOfSecondaries.end(); it2++){ 346 G4String name = it2->first; << 370 G4String name = it2->first; 347 MinMaxData data = it2->second; << 371 MinMaxData data = it2->second; 348 G4int count = data.fCount; << 372 G4int count = data.fCount; 349 G4double eMean = data.fVsum / count; << 373 G4double eMean = data.fVsum/count; 350 G4double eMin = data.fVmin; << 374 G4double eMin = data.fVmin; 351 G4double eMax = data.fVmax; << 375 G4double eMax = data.fVmax; 352 << 376 353 G4cout << " " << std::setw(13) << name << << 377 G4cout << " " << std::setw(13) << name << ": " << std::setw(7) << count 354 << " Emean = " << std::setw(6) << << 378 << " Emean = " << std::setw(6) << G4BestUnit(eMean, "Energy") 355 << G4BestUnit(eMin, "Energy") << " << 379 << " (" << G4BestUnit(eMin, "Energy") >> 380 << " --> " << G4BestUnit(eMax, "Energy") >> 381 << ")" << G4endl; 356 } 382 } 357 G4cout << G4endl; 383 G4cout << G4endl; 358 << 384 359 // continuous energy deposited by secondary << 385 //continuous energy deposited by secondary tracks dE5 360 // (only if secondary particles are tracked << 386 // (only if secondary particles are tracked) 361 // 387 // 362 if (fEdepSecondary > 0.) { 388 if (fEdepSecondary > 0.) { 363 energyPerEvent = fEdepSecondary / nbEvents << 389 energyPerEvent = fEdepSecondary/nbEvents; 364 390 365 G4cout << "\n Energy continuously deposite << 391 G4cout 366 << " (restricted dE/dx) dE5 = " << << 392 << "\n Energy continuously deposited along secondary tracks" 367 << G4BestUnit(fEdepSecMin, "Energy" << 393 << " (restricted dE/dx) dE5 = " 368 << ")" << G4endl; << 394 << G4BestUnit(energyPerEvent, "Energy") >> 395 << " (" << G4BestUnit(fEdepSecMin, "Energy") >> 396 << " --> " << G4BestUnit(fEdepSecMax, "Energy") << ")" >> 397 << G4endl; 369 398 370 // total energy deposited : dE6 = dE1 + dE 399 // total energy deposited : dE6 = dE1 + dE5 371 // 400 // 372 energyPerEvent = fEdepTotal / nbEvents; << 401 energyPerEvent = fEdepTotal/nbEvents; 373 << 374 G4cout << "\n Total energy deposited : dE6 << 375 << G4BestUnit(energyPerEvent, "Ener << 376 << " --> " << G4BestUnit(fEdepTotMa << 377 << G4endl; << 378 } << 379 402 >> 403 G4cout >> 404 << "\n Total energy deposited : dE6 = dE1 + dE5 = " >> 405 << G4BestUnit(energyPerEvent, "Energy") >> 406 << " (" << G4BestUnit(fEdepTotMin, "Energy") >> 407 << " --> " << G4BestUnit(fEdepTotMax, "Energy") << ") \n" >> 408 << G4endl; >> 409 } >> 410 380 G4cout.precision(prec); 411 G4cout.precision(prec); 381 << 412 382 // clear maps << 413 //clear maps 383 // 414 // 384 fProcCounter.clear(); 415 fProcCounter.clear(); 385 fEtransfByProcess.clear(); 416 fEtransfByProcess.clear(); 386 fEkinOfSecondaries.clear(); 417 fEkinOfSecondaries.clear(); 387 418 388 // save histograms << 419 //save histograms 389 G4AnalysisManager* analysisManager = G4Analy << 420 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance(); 390 if (analysisManager->IsActive()) { << 421 if ( analysisManager->IsActive() ) { 391 analysisManager->Write(); 422 analysisManager->Write(); 392 analysisManager->CloseFile(); 423 analysisManager->CloseFile(); 393 } 424 } 394 425 395 // show Rndm status 426 // show Rndm status 396 CLHEP::HepRandom::showEngineStatus(); 427 CLHEP::HepRandom::showEngineStatus(); 397 } 428 } 398 429 399 //....oooOO0OOooo........oooOO0OOooo........oo 430 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 400 431 401 G4double RunAction::GetEnergyFromRestrictedRan << 432 G4double RunAction::GetEnergyFromRestrictedRange(G4double range, 402 << 433 G4ParticleDefinition* particle, G4Material* material, G4double Etry) 403 { 434 { 404 G4EmCalculator emCal; 435 G4EmCalculator emCal; 405 << 436 406 G4double Energy = Etry, dE = 0., dEdx; 437 G4double Energy = Etry, dE = 0., dEdx; 407 G4double r, dr; 438 G4double r, dr; 408 G4double err = 1., errmax = 0.00001; << 439 G4double err = 1., errmax = 0.00001; 409 G4int iter = 0, itermax = 10; << 440 G4int iter = 0 , itermax = 10; 410 while (err > errmax && iter < itermax) { 441 while (err > errmax && iter < itermax) { 411 iter++; 442 iter++; 412 Energy -= dE; 443 Energy -= dE; 413 r = emCal.GetRangeFromRestricteDEDX(Energy << 444 r = emCal.GetRangeFromRestricteDEDX(Energy,particle,material); 414 dr = r - range; << 445 dr = r - range; 415 dEdx = emCal.GetDEDX(Energy, particle, mat << 446 dEdx = emCal.GetDEDX(Energy,particle,material); 416 dE = dEdx * dr; << 447 dE = dEdx*dr; 417 err = std::abs(dE) / Energy; << 448 err = std::abs(dE)/Energy; 418 } 449 } 419 if (iter == itermax) { 450 if (iter == itermax) { 420 G4cout << "\n ---> warning: RunAction::Ge << 451 G4cout 421 << " Etry = " << G4BestUnit(Etry, << 452 << "\n ---> warning: RunAction::GetEnergyFromRestRange() did not converge" 422 << " Energy = " << G4BestUnit(Ene << 453 << " Etry = " << G4BestUnit(Etry,"Energy") 423 << " iter = " << iter << G4endl; << 454 << " Energy = " << G4BestUnit(Energy,"Energy") 424 } << 455 << " err = " << err 425 << 456 << " iter = " << iter << G4endl; >> 457 } >> 458 426 return Energy; 459 return Energy; 427 } 460 } 428 461 429 //....oooOO0OOooo........oooOO0OOooo........oo 462 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 430 463 431 G4double RunAction::GetEnergyFromCSDARange(G4d << 464 G4double RunAction::GetEnergyFromCSDARange(G4double range, 432 G4M << 465 G4ParticleDefinition* particle, G4Material* material, G4double Etry) 433 { 466 { 434 G4EmCalculator emCal; 467 G4EmCalculator emCal; 435 << 468 436 G4double Energy = Etry, dE = 0., dEdx; 469 G4double Energy = Etry, dE = 0., dEdx; 437 G4double r, dr; 470 G4double r, dr; 438 G4double err = 1., errmax = 0.00001; << 471 G4double err = 1., errmax = 0.00001; 439 G4int iter = 0, itermax = 10; << 472 G4int iter = 0 , itermax = 10; 440 while (err > errmax && iter < itermax) { 473 while (err > errmax && iter < itermax) { 441 iter++; 474 iter++; 442 Energy -= dE; 475 Energy -= dE; 443 r = emCal.GetCSDARange(Energy, particle, m << 476 r = emCal.GetCSDARange(Energy,particle,material); 444 dr = r - range; << 477 dr = r - range; 445 dEdx = emCal.ComputeTotalDEDX(Energy, part << 478 dEdx = emCal.ComputeTotalDEDX(Energy,particle,material); 446 dE = dEdx * dr; << 479 dE = dEdx*dr; 447 err = std::abs(dE) / Energy; << 480 err = std::abs(dE)/Energy; 448 } 481 } 449 if (iter == itermax) { 482 if (iter == itermax) { 450 G4cout << "\n ---> warning: RunAction::Ge << 483 G4cout 451 << " Etry = " << G4BestUnit(Etry, << 484 << "\n ---> warning: RunAction::GetEnergyFromCSDARange() did not converge" 452 << " Energy = " << G4BestUnit(Ene << 485 << " Etry = " << G4BestUnit(Etry,"Energy") 453 << " iter = " << iter << G4endl; << 486 << " Energy = " << G4BestUnit(Energy,"Energy") 454 } << 487 << " err = " << err 455 << 488 << " iter = " << iter << G4endl; >> 489 } >> 490 456 return Energy; 491 return Energy; 457 } 492 } 458 493 459 //....oooOO0OOooo........oooOO0OOooo........oo 494 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 460 495