File indexing completed on 2026-06-14 07:53:03
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 "Run.hh"
0030
0031 #include "DetectorConstruction.hh"
0032 #include "EventAction.hh"
0033 #include "HistoManager.hh"
0034 #include "PrimaryGeneratorAction.hh"
0035
0036 #include "G4Event.hh"
0037 #include "G4Material.hh"
0038 #include "G4SystemOfUnits.hh"
0039 #include "G4UnitsTable.hh"
0040
0041 #include <iomanip>
0042
0043
0044
0045 Run::Run(DetectorConstruction* detector) : fDetector(detector)
0046 {
0047 for (G4int i = 0; i < 3; ++i) {
0048 fStatus[i] = 0;
0049 fTotEdep[i] = fEleak[i] = fEtotal[i] = 0.;
0050 }
0051 fTotEdep[1] = fEleak[1] = fEtotal[1] = joule;
0052
0053 for (G4int i = 0; i < kMaxAbsor; ++i) {
0054 fEdeposit[i] = 0.;
0055 fEmin[i] = joule;
0056 fEmax[i] = 0.;
0057 fCsdaRange[i] = 0.;
0058 fXfrontNorm[i] = 0.;
0059 }
0060 }
0061
0062
0063
0064 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
0065 {
0066 fParticle = particle;
0067 fEkin = energy;
0068 }
0069
0070
0071
0072 void Run::AddEdep(G4int i, G4double e)
0073 {
0074 if (e > 0.) {
0075 fEdeposit[i] += e;
0076 if (e < fEmin[i]) fEmin[i] = e;
0077 if (e > fEmax[i]) fEmax[i] = e;
0078 }
0079 }
0080
0081
0082
0083 void Run::AddTotEdep(G4double e)
0084 {
0085 if (e > 0.) {
0086 fTotEdep[0] += e;
0087 if (e < fTotEdep[1]) fTotEdep[1] = e;
0088 if (e > fTotEdep[2]) fTotEdep[2] = e;
0089 }
0090 }
0091
0092
0093
0094 void Run::AddEleak(G4double e)
0095 {
0096 if (e > 0.) {
0097 fEleak[0] += e;
0098 if (e < fEleak[1]) fEleak[1] = e;
0099 if (e > fEleak[2]) fEleak[2] = e;
0100 }
0101 }
0102
0103
0104
0105 void Run::AddEtotal(G4double e)
0106 {
0107 if (e > 0.) {
0108 fEtotal[0] += e;
0109 if (e < fEtotal[1]) fEtotal[1] = e;
0110 if (e > fEtotal[2]) fEtotal[2] = e;
0111 }
0112 }
0113
0114
0115
0116 void Run::AddTrackLength(G4double t)
0117 {
0118 fTrackLen += t;
0119 fTrackLen2 += t * t;
0120 }
0121
0122
0123
0124 void Run::AddProjRange(G4double x)
0125 {
0126 fProjRange += x;
0127 fProjRange2 += x * x;
0128 }
0129
0130
0131
0132 void Run::AddStepSize(G4int nb, G4double st)
0133 {
0134 fNbOfSteps += nb;
0135 fNbOfSteps2 += nb * nb;
0136 fStepSize += st;
0137 fStepSize2 += st * st;
0138 }
0139
0140
0141
0142 void Run::AddTrackStatus(G4int i)
0143 {
0144 fStatus[i]++;
0145 }
0146
0147
0148
0149 void Run::SetCsdaRange(G4int i, G4double value)
0150 {
0151 fCsdaRange[i] = value;
0152 }
0153
0154
0155
0156 void Run::SetXfrontNorm(G4int i, G4double value)
0157 {
0158 fXfrontNorm[i] = value;
0159 }
0160
0161
0162
0163 G4double Run::GetCsdaRange(G4int i)
0164 {
0165 return fCsdaRange[i];
0166 }
0167
0168
0169
0170 G4double Run::GetXfrontNorm(G4int i)
0171 {
0172 return fXfrontNorm[i];
0173 }
0174
0175
0176
0177 void Run::Merge(const G4Run* run)
0178 {
0179 const Run* localRun = static_cast<const Run*>(run);
0180
0181
0182 fParticle = localRun->fParticle;
0183 fEkin = localRun->fEkin;
0184
0185
0186 fTrackLen += localRun->fTrackLen;
0187 fTrackLen2 += localRun->fTrackLen2;
0188 fProjRange += localRun->fProjRange;
0189 fProjRange2 += localRun->fProjRange2;
0190 fNbOfSteps += localRun->fNbOfSteps;
0191 fNbOfSteps2 += localRun->fNbOfSteps2;
0192 fStepSize += localRun->fStepSize;
0193 fStepSize2 += localRun->fStepSize2;
0194
0195 G4int nbOfAbsor = fDetector->GetNbOfAbsor();
0196 for (G4int i = 1; i <= nbOfAbsor; ++i) {
0197 fEdeposit[i] += localRun->fEdeposit[i];
0198 fCsdaRange[i] = localRun->fCsdaRange[i];
0199 fXfrontNorm[i] = localRun->fXfrontNorm[i];
0200
0201 G4double min, max;
0202 min = localRun->fEmin[i];
0203 max = localRun->fEmax[i];
0204 if (fEmin[i] > min) fEmin[i] = min;
0205 if (fEmax[i] < max) fEmax[i] = max;
0206 }
0207
0208 for (G4int i = 0; i < 3; ++i)
0209 fStatus[i] += localRun->fStatus[i];
0210
0211
0212 fTotEdep[0] += localRun->fTotEdep[0];
0213 G4double min, max;
0214 min = localRun->fTotEdep[1];
0215 max = localRun->fTotEdep[2];
0216 if (fTotEdep[1] > min) fTotEdep[1] = min;
0217 if (fTotEdep[2] < max) fTotEdep[2] = max;
0218
0219
0220 fEleak[0] += localRun->fEleak[0];
0221 min = localRun->fEleak[1];
0222 max = localRun->fEleak[2];
0223 if (fEleak[1] > min) fEleak[1] = min;
0224 if (fEleak[2] < max) fEleak[2] = max;
0225
0226
0227 fEtotal[0] += localRun->fEtotal[0];
0228 min = localRun->fEtotal[1];
0229 max = localRun->fEtotal[2];
0230 if (fEtotal[1] > min) fEtotal[1] = min;
0231 if (fEtotal[2] < max) fEtotal[2] = max;
0232
0233 G4Run::Merge(run);
0234 }
0235
0236
0237
0238 void Run::EndOfRun()
0239 {
0240 std::ios::fmtflags mode = G4cout.flags();
0241 G4cout.setf(std::ios::fixed, std::ios::floatfield);
0242 G4int prec = G4cout.precision(2);
0243
0244
0245
0246 G4String partName = fParticle->GetParticleName();
0247 G4int nbOfAbsor = fDetector->GetNbOfAbsor();
0248
0249 G4cout << "\n ======================== run summary =====================\n";
0250 G4cout << "\n The run is " << numberOfEvent << " " << partName << " of "
0251 << G4BestUnit(fEkin, "Energy") << " through " << nbOfAbsor << " absorbers: \n";
0252 for (G4int i = 1; i <= nbOfAbsor; i++) {
0253 G4Material* material = fDetector->GetAbsorMaterial(i);
0254 G4double thickness = fDetector->GetAbsorThickness(i);
0255 G4double density = material->GetDensity();
0256 G4cout << std::setw(5) << i << std::setw(10) << G4BestUnit(thickness, "Length") << " of "
0257 << material->GetName() << " (density: " << G4BestUnit(density, "Volumic Mass") << ")"
0258 << G4endl;
0259 }
0260
0261 if (numberOfEvent == 0) {
0262 G4cout.setf(mode, std::ios::floatfield);
0263 G4cout.precision(prec);
0264 return;
0265 }
0266
0267 G4cout.precision(3);
0268 G4double rms(0);
0269
0270
0271
0272 for (G4int i = 1; i <= nbOfAbsor; i++) {
0273 fEdeposit[i] /= numberOfEvent;
0274
0275 G4cout << "\n Edep in absorber " << i << " = " << G4BestUnit(fEdeposit[i], "Energy") << "\t("
0276 << G4BestUnit(fEmin[i], "Energy") << "-->" << G4BestUnit(fEmax[i], "Energy") << ")";
0277 }
0278 G4cout << G4endl;
0279
0280 if (nbOfAbsor > 1) {
0281 fTotEdep[0] /= numberOfEvent;
0282 G4cout << "\n Edep in all absorbers = " << G4BestUnit(fTotEdep[0], "Energy") << "\t("
0283 << G4BestUnit(fTotEdep[1], "Energy") << "-->" << G4BestUnit(fTotEdep[2], "Energy") << ")"
0284 << G4endl;
0285 }
0286
0287
0288
0289 fEleak[0] /= numberOfEvent;
0290 G4cout << " Energy leakage = " << G4BestUnit(fEleak[0], "Energy") << "\t("
0291 << G4BestUnit(fEleak[1], "Energy") << "-->" << G4BestUnit(fEleak[2], "Energy") << ")"
0292 << G4endl;
0293
0294
0295
0296 fEtotal[0] /= numberOfEvent;
0297 G4cout << " Energy total = " << G4BestUnit(fEtotal[0], "Energy") << "\t("
0298 << G4BestUnit(fEtotal[1], "Energy") << "-->" << G4BestUnit(fEtotal[2], "Energy") << ")"
0299 << G4endl;
0300
0301
0302
0303 fTrackLen /= numberOfEvent;
0304 fTrackLen2 /= numberOfEvent;
0305 rms = fTrackLen2 - fTrackLen * fTrackLen;
0306 if (rms > 0.)
0307 rms = std::sqrt(rms);
0308 else
0309 rms = 0.;
0310
0311 G4cout.precision(3);
0312 G4cout << "\n Track length of primary track = " << G4BestUnit(fTrackLen, "Length") << " +- "
0313 << G4BestUnit(rms, "Length");
0314
0315
0316
0317 G4int NbOfAbsor = fDetector->GetNbOfAbsor();
0318 if (NbOfAbsor == 1) {
0319 G4cout << "\n Range from EmCalculator = " << G4BestUnit(fCsdaRange[1], "Length")
0320 << " (from full dE/dx)" << G4endl;
0321 }
0322
0323
0324
0325 fProjRange /= numberOfEvent;
0326 fProjRange2 /= numberOfEvent;
0327 rms = fProjRange2 - fProjRange * fProjRange;
0328 if (rms > 0.)
0329 rms = std::sqrt(rms);
0330 else
0331 rms = 0.;
0332
0333 G4cout << "\n Projected range = " << G4BestUnit(fProjRange, "Length") << " +- "
0334 << G4BestUnit(rms, "Length") << G4endl;
0335
0336
0337
0338 G4double dNofEvents = double(numberOfEvent);
0339 G4double fNbSteps = fNbOfSteps / dNofEvents, fNbSteps2 = fNbOfSteps2 / dNofEvents;
0340 rms = fNbSteps2 - fNbSteps * fNbSteps;
0341 if (rms > 0.)
0342 rms = std::sqrt(rms);
0343 else
0344 rms = 0.;
0345
0346 G4cout.precision(2);
0347 G4cout << "\n Nb of steps of primary track = " << fNbSteps << " +- " << rms;
0348
0349 fStepSize /= numberOfEvent;
0350 fStepSize2 /= numberOfEvent;
0351 rms = fStepSize2 - fStepSize * fStepSize;
0352 if (rms > 0.)
0353 rms = std::sqrt(rms);
0354 else
0355 rms = 0.;
0356
0357 G4cout.precision(3);
0358 G4cout << "\t Step size= " << G4BestUnit(fStepSize, "Length") << " +- "
0359 << G4BestUnit(rms, "Length") << G4endl;
0360
0361
0362
0363 G4double absorbed = 100. * fStatus[0] / dNofEvents;
0364 G4double transmit = 100. * fStatus[1] / dNofEvents;
0365 G4double reflected = 100. * fStatus[2] / dNofEvents;
0366
0367 G4cout.precision(2);
0368 G4cout << "\n absorbed = " << absorbed << " %"
0369 << " transmit = " << transmit << " %"
0370 << " reflected = " << reflected << " % \n"
0371 << G4endl;
0372
0373
0374
0375 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0376 G4int ih = 1;
0377 G4double binWidth = analysisManager->GetH1Width(ih) * analysisManager->GetH1Unit(ih);
0378 G4double fac = (1. / (numberOfEvent * binWidth)) * (mm / MeV);
0379 analysisManager->ScaleH1(ih, fac);
0380
0381 ih = 8;
0382 binWidth = analysisManager->GetH1Width(ih);
0383 fac = (1. / (numberOfEvent * binWidth)) * (g / (MeV * cm2));
0384 analysisManager->ScaleH1(ih, fac);
0385
0386
0387 G4cout.setf(mode, std::ios::floatfield);
0388 G4cout.precision(prec);
0389 }
0390
0391
