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