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