File indexing completed on 2026-09-15 08:28:19
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029 #include "RunAction.hh"
0030
0031 #include "DetectorConstruction.hh"
0032 #include "HistoManager.hh"
0033 #include "PrimaryGeneratorAction.hh"
0034
0035 #include "G4EmCalculator.hh"
0036 #include "G4Run.hh"
0037 #include "G4UnitsTable.hh"
0038 #include "Randomize.hh"
0039
0040 #include <iomanip>
0041
0042
0043
0044 RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* kin)
0045 : fDetector(det), fPrimary(kin)
0046 {
0047 fHistoManager = new HistoManager();
0048 }
0049
0050
0051
0052 RunAction::~RunAction()
0053 {
0054 delete fHistoManager;
0055 }
0056
0057
0058
0059 void RunAction::BeginOfRunAction(const G4Run*)
0060 {
0061
0062
0063 fNbSteps = 0;
0064 fTrackLength = 0.;
0065 fStepMin = DBL_MAX;
0066 fStepMax = 0.;
0067
0068 fEdepPrimary = fEdepSecondary = fEdepTotal = 0.;
0069 fEdepPrimMin = fEdepSecMin = fEdepTotMin = DBL_MAX;
0070 fEdepPrimMax = fEdepSecMax = fEdepTotMax = 0.;
0071
0072 fEnergyTransfered = 0.;
0073 fEtransfMin = DBL_MAX;
0074 fEtransfMax = 0.;
0075
0076 fEnergyLost = 0.;
0077 fElostMin = DBL_MAX;
0078 fElostMax = 0.;
0079
0080 fEnergyBalance = 0.;
0081 fEbalMin = DBL_MAX;
0082 fEbalMax = 0.;
0083
0084
0085
0086 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0087 if (analysisManager->IsActive()) {
0088 analysisManager->OpenFile();
0089 }
0090
0091
0092 CLHEP::HepRandom::showEngineStatus();
0093 }
0094
0095
0096
0097 void RunAction::CountProcesses(G4String procName)
0098 {
0099 std::map<G4String, G4int>::iterator it = fProcCounter.find(procName);
0100 if (it == fProcCounter.end()) {
0101 fProcCounter[procName] = 1;
0102 }
0103 else {
0104 fProcCounter[procName]++;
0105 }
0106 }
0107
0108
0109
0110 void RunAction::TrackLength(G4double step)
0111 {
0112 fTrackLength += step;
0113 fNbSteps++;
0114 if (step < fStepMin) fStepMin = step;
0115 if (step > fStepMax) fStepMax = step;
0116 }
0117
0118
0119
0120 void RunAction::EnergyDeposited(G4double edepPrim, G4double edepSecond)
0121 {
0122 fEdepPrimary += edepPrim;
0123 if (edepPrim < fEdepPrimMin) fEdepPrimMin = edepPrim;
0124 if (edepPrim > fEdepPrimMax) fEdepPrimMax = edepPrim;
0125
0126 fEdepSecondary += edepSecond;
0127 if (edepSecond < fEdepSecMin) fEdepSecMin = edepSecond;
0128 if (edepSecond > fEdepSecMax) fEdepSecMax = edepSecond;
0129 }
0130
0131
0132
0133 void RunAction::EnergyTransferedByProcess(G4String process, G4double energy)
0134 {
0135 std::map<G4String, MinMaxData>::iterator it = fEtransfByProcess.find(process);
0136 if (it == fEtransfByProcess.end()) {
0137 fEtransfByProcess[process] = MinMaxData(1, energy, energy, energy);
0138 }
0139 else {
0140 MinMaxData& data = it->second;
0141 data.fCount++;
0142 data.fVsum += energy;
0143
0144 G4double emin = data.fVmin;
0145 if (energy < emin) data.fVmin = energy;
0146 G4double emax = data.fVmax;
0147 if (energy > emax) data.fVmax = energy;
0148 }
0149 }
0150
0151
0152
0153 void RunAction::EnergyTransfered(G4double energy)
0154 {
0155 fEnergyTransfered += energy;
0156 if (energy < fEtransfMin) fEtransfMin = energy;
0157 if (energy > fEtransfMax) fEtransfMax = energy;
0158 }
0159
0160
0161
0162 void RunAction::TotalEnergyLost(G4double energy)
0163 {
0164 fEnergyLost += energy;
0165 if (energy < fElostMin) fElostMin = energy;
0166 if (energy > fElostMax) fElostMax = energy;
0167 }
0168
0169
0170
0171 void RunAction::EnergyBalance(G4double energy)
0172 {
0173 fEnergyBalance += energy;
0174 if (energy < fEbalMin) fEbalMin = energy;
0175 if (energy > fEbalMax) fEbalMax = energy;
0176 }
0177
0178
0179
0180 void RunAction::TotalEnergyDeposit(G4double energy)
0181 {
0182 fEdepTotal += energy;
0183 if (energy < fEdepTotMin) fEdepTotMin = energy;
0184 if (energy > fEdepTotMax) fEdepTotMax = energy;
0185 }
0186
0187
0188
0189 void RunAction::EnergySpectrumOfSecondaries(G4String particle, G4double energy)
0190 {
0191 std::map<G4String, MinMaxData>::iterator it = fEkinOfSecondaries.find(particle);
0192 if (it == fEkinOfSecondaries.end()) {
0193 fEkinOfSecondaries[particle] = MinMaxData(1, energy, energy, energy);
0194 }
0195 else {
0196 MinMaxData& data = it->second;
0197 data.fCount++;
0198 data.fVsum += energy;
0199
0200 G4double emin = data.fVmin;
0201 if (energy < emin) data.fVmin = energy;
0202 G4double emax = data.fVmax;
0203 if (energy > emax) data.fVmax = energy;
0204 }
0205 }
0206
0207
0208
0209 void RunAction::EndOfRunAction(const G4Run* aRun)
0210 {
0211 G4int nbEvents = aRun->GetNumberOfEvent();
0212 if (nbEvents == 0) return;
0213
0214 G4Material* material = fDetector->GetMaterial();
0215 G4double length = fDetector->GetSize();
0216 G4double density = material->GetDensity();
0217
0218 G4ParticleDefinition* particle = fPrimary->GetParticleGun()->GetParticleDefinition();
0219 G4String partName = particle->GetParticleName();
0220 G4double ePrimary = fPrimary->GetParticleGun()->GetParticleEnergy();
0221
0222 G4int prec = G4cout.precision(3);
0223 G4cout << "\n ======================== run summary ======================\n";
0224 G4cout << "\n The run was " << nbEvents << " " << partName << " of "
0225 << G4BestUnit(ePrimary, "Energy") << " through " << G4BestUnit(length, "Length") << " of "
0226 << material->GetName() << " (density: " << G4BestUnit(density, "Volumic Mass") << ")";
0227 G4cout << G4endl;
0228
0229 if (particle->GetPDGCharge() == 0.) return;
0230
0231 G4cout.precision(4);
0232
0233
0234
0235 G4cout << "\n Process defining step :" << G4endl;
0236 G4int index = 0;
0237 for (const auto& procCounter : fProcCounter) {
0238 G4String procName = procCounter.first;
0239 G4int count = procCounter.second;
0240 G4String space = " ";
0241 if (++index % 4 == 0) space = "\n";
0242 G4cout << " " << std::setw(15) << procName << "=" << std::setw(7) << count << space;
0243 }
0244 G4cout << G4endl;
0245
0246
0247
0248 G4double trackLPerEvent = fTrackLength / nbEvents;
0249 G4double nbStepPerEvent = double(fNbSteps) / nbEvents;
0250 G4double stepSize = fTrackLength / fNbSteps;
0251
0252 G4cout << "\n TrackLength = " << G4BestUnit(trackLPerEvent, "Length")
0253 << " nb of steps = " << nbStepPerEvent
0254 << " stepSize = " << G4BestUnit(stepSize, "Length") << " ("
0255 << G4BestUnit(fStepMin, "Length") << "--> " << G4BestUnit(fStepMax, "Length") << ")"
0256 << G4endl;
0257
0258
0259
0260 G4double energyPerEvent = fEdepPrimary / nbEvents;
0261
0262 G4cout << "\n Energy continuously deposited along primary track"
0263 << " (restricted dE/dx) dE1 = " << G4BestUnit(energyPerEvent, "Energy") << " ("
0264 << G4BestUnit(fEdepPrimMin, "Energy") << " --> " << G4BestUnit(fEdepPrimMax, "Energy")
0265 << ")" << G4endl;
0266
0267
0268
0269 G4EmCalculator emCal;
0270
0271 G4double r0 = emCal.GetRangeFromRestricteDEDX(ePrimary, particle, material);
0272 G4double r1 = r0 - trackLPerEvent;
0273 G4double etry = ePrimary - energyPerEvent;
0274 G4double efinal = 0.;
0275 if (r1 > 0.) efinal = GetEnergyFromRestrictedRange(r1, particle, material, etry);
0276 G4double dEtable = ePrimary - efinal;
0277 G4double ratio = 0.;
0278 if (dEtable > 0.) ratio = energyPerEvent / dEtable;
0279
0280 G4cout << "\n Evaluation of dE1 from reading restricted Range table : dE1_table = "
0281 << G4BestUnit(dEtable, "Energy") << " ---> dE1/dE1_table = " << ratio << G4endl;
0282
0283
0284
0285 G4cout << "\n Energy transfered to secondary particles :" << G4endl;
0286 std::map<G4String, MinMaxData>::iterator it1;
0287 for (it1 = fEtransfByProcess.begin(); it1 != fEtransfByProcess.end(); it1++) {
0288 G4String name = it1->first;
0289 MinMaxData data = it1->second;
0290 energyPerEvent = data.fVsum / nbEvents;
0291 G4double eMin = data.fVmin;
0292 G4double eMax = data.fVmax;
0293
0294 G4cout << " " << std::setw(17) << "due to " + name << ": dE2 = " << std::setw(6)
0295 << G4BestUnit(energyPerEvent, "Energy") << " (" << G4BestUnit(eMin, "Energy") << " --> "
0296 << G4BestUnit(eMax, "Energy") << ")" << G4endl;
0297 }
0298
0299
0300
0301 energyPerEvent = fEnergyTransfered / nbEvents;
0302
0303 G4cout << "\n Total energy transfered to secondaries : dE3 = sum of dE2 = "
0304 << G4BestUnit(energyPerEvent, "Energy") << " (" << G4BestUnit(fEtransfMin, "Energy")
0305 << " --> " << G4BestUnit(fEtransfMax, "Energy") << ")" << G4endl;
0306
0307
0308
0309 energyPerEvent = fEnergyLost / nbEvents;
0310
0311 G4cout << "\n Total energy lost by incident particle : dE4 = dE1 + dE3 = "
0312 << G4BestUnit(energyPerEvent, "Energy") << " (" << G4BestUnit(fElostMin, "Energy")
0313 << " --> " << G4BestUnit(fElostMax, "Energy") << ")" << G4endl;
0314
0315
0316
0317 energyPerEvent = fEnergyBalance / nbEvents;
0318
0319 G4cout << "\n calcul of dE4 from energy balance : dE4_bal = E_in - E_out = "
0320 << G4BestUnit(energyPerEvent, "Energy") << " (" << G4BestUnit(fEbalMin, "Energy")
0321 << " --> " << G4BestUnit(fEbalMax, "Energy") << ")" << G4endl;
0322
0323
0324
0325 r0 = emCal.GetCSDARange(ePrimary, particle, material);
0326 r1 = r0 - trackLPerEvent;
0327 etry = ePrimary - energyPerEvent;
0328 efinal = 0.;
0329 if (r1 > 0.) efinal = GetEnergyFromCSDARange(r1, particle, material, etry);
0330 dEtable = ePrimary - efinal;
0331 ratio = 0.;
0332 if (dEtable > 0.) ratio = energyPerEvent / dEtable;
0333
0334 G4cout << "\n Evaluation of dE4 from reading full Range table : dE4_table = "
0335 << G4BestUnit(dEtable, "Energy") << " ---> dE4/dE4_table = " << ratio << G4endl;
0336
0337
0338
0339 G4cout << "\n Energy spectrum of secondary particles :" << G4endl;
0340 std::map<G4String, MinMaxData>::iterator it2;
0341 for (it2 = fEkinOfSecondaries.begin(); it2 != fEkinOfSecondaries.end(); it2++) {
0342 G4String name = it2->first;
0343 MinMaxData data = it2->second;
0344 G4int count = data.fCount;
0345 G4double eMean = data.fVsum / count;
0346 G4double eMin = data.fVmin;
0347 G4double eMax = data.fVmax;
0348
0349 G4cout << " " << std::setw(13) << name << ": " << std::setw(7) << count
0350 << " Emean = " << std::setw(6) << G4BestUnit(eMean, "Energy") << " ("
0351 << G4BestUnit(eMin, "Energy") << " --> " << G4BestUnit(eMax, "Energy") << ")" << G4endl;
0352 }
0353 G4cout << G4endl;
0354
0355
0356
0357
0358 if (fEdepSecondary > 0.) {
0359 energyPerEvent = fEdepSecondary / nbEvents;
0360
0361 G4cout << "\n Energy continuously deposited along secondary tracks"
0362 << " (restricted dE/dx) dE5 = " << G4BestUnit(energyPerEvent, "Energy") << " ("
0363 << G4BestUnit(fEdepSecMin, "Energy") << " --> " << G4BestUnit(fEdepSecMax, "Energy")
0364 << ")" << G4endl;
0365
0366
0367
0368 energyPerEvent = fEdepTotal / nbEvents;
0369
0370 G4cout << "\n Total energy deposited : dE6 = dE1 + dE5 = "
0371 << G4BestUnit(energyPerEvent, "Energy") << " (" << G4BestUnit(fEdepTotMin, "Energy")
0372 << " --> " << G4BestUnit(fEdepTotMax, "Energy") << ") \n"
0373 << G4endl;
0374 }
0375
0376 G4cout.precision(prec);
0377
0378
0379
0380 fProcCounter.clear();
0381 fEtransfByProcess.clear();
0382 fEkinOfSecondaries.clear();
0383
0384
0385 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0386 if (analysisManager->IsActive()) {
0387 analysisManager->Write();
0388 analysisManager->CloseFile();
0389 }
0390
0391
0392 CLHEP::HepRandom::showEngineStatus();
0393 }
0394
0395
0396
0397 G4double RunAction::GetEnergyFromRestrictedRange(G4double range, G4ParticleDefinition* particle,
0398 G4Material* material, G4double Etry)
0399 {
0400 G4EmCalculator emCal;
0401
0402 G4double Energy = Etry, dE = 0., dEdx;
0403 G4double r, dr;
0404 G4double err = 1., errmax = 0.00001;
0405 G4int iter = 0, itermax = 10;
0406 while (err > errmax && iter < itermax) {
0407 iter++;
0408 Energy -= dE;
0409 r = emCal.GetRangeFromRestricteDEDX(Energy, particle, material);
0410 dr = r - range;
0411 dEdx = emCal.GetDEDX(Energy, particle, material);
0412 dE = dEdx * dr;
0413 err = std::abs(dE) / Energy;
0414 }
0415 if (iter == itermax) {
0416 G4cout << "\n ---> warning: RunAction::GetEnergyFromRestRange() did not converge"
0417 << " Etry = " << G4BestUnit(Etry, "Energy")
0418 << " Energy = " << G4BestUnit(Energy, "Energy") << " err = " << err
0419 << " iter = " << iter << G4endl;
0420 }
0421
0422 return Energy;
0423 }
0424
0425
0426
0427 G4double RunAction::GetEnergyFromCSDARange(G4double range, G4ParticleDefinition* particle,
0428 G4Material* material, G4double Etry)
0429 {
0430 G4EmCalculator emCal;
0431
0432 G4double Energy = Etry, dE = 0., dEdx;
0433 G4double r, dr;
0434 G4double err = 1., errmax = 0.00001;
0435 G4int iter = 0, itermax = 10;
0436 while (err > errmax && iter < itermax) {
0437 iter++;
0438 Energy -= dE;
0439 r = emCal.GetCSDARange(Energy, particle, material);
0440 dr = r - range;
0441 dEdx = emCal.ComputeTotalDEDX(Energy, particle, material);
0442 dE = dEdx * dr;
0443 err = std::abs(dE) / Energy;
0444 }
0445 if (iter == itermax) {
0446 G4cout << "\n ---> warning: RunAction::GetEnergyFromCSDARange() did not converge"
0447 << " Etry = " << G4BestUnit(Etry, "Energy")
0448 << " Energy = " << G4BestUnit(Energy, "Energy") << " err = " << err
0449 << " iter = " << iter << G4endl;
0450 }
0451
0452 return Energy;
0453 }
0454
0455
