File indexing completed on 2025-02-23 09:22:15
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 "HistoManager.hh"
0037 #include "PrimaryGeneratorAction.hh"
0038
0039 #include "G4Electron.hh"
0040 #include "G4EmCalculator.hh"
0041 #include "G4Run.hh"
0042 #include "G4RunManager.hh"
0043 #include "G4SystemOfUnits.hh"
0044 #include "G4UnitsTable.hh"
0045 #include "Randomize.hh"
0046
0047 #include <iomanip>
0048
0049
0050
0051 Run::Run(DetectorConstruction* det, PrimaryGeneratorAction* kin)
0052 : fDetector(det), fKinematic(kin), fProcCounter(0), fMateWall(0), fMateCavity(0)
0053
0054 {
0055
0056
0057 fWallThickness = fDetector->GetWallThickness();
0058 fWallRadius = fDetector->GetWallRadius();
0059 fMateWall = fDetector->GetWallMaterial();
0060 fDensityWall = fMateWall->GetDensity();
0061
0062 fCavityThickness = fDetector->GetCavityThickness();
0063 fCavityRadius = fDetector->GetCavityRadius();
0064 fSurfaceCavity = CLHEP::pi * fCavityRadius * fCavityRadius;
0065 fVolumeCavity = fSurfaceCavity * fCavityThickness;
0066 fMateCavity = fDetector->GetCavityMaterial();
0067 fDensityCavity = fMateCavity->GetDensity();
0068 fMassCavity = fVolumeCavity * fDensityCavity;
0069
0070
0071
0072 fProcCounter = new ProcessesCount;
0073
0074
0075
0076 fEsecondary = fEsecondary2 = 0.;
0077 fNbSec = 0;
0078
0079
0080
0081 fPartFlowCavity[0] = fPartFlowCavity[1] = 0;
0082 fEnerFlowCavity[0] = fEnerFlowCavity[1] = 0.;
0083
0084
0085
0086 fEdepCavity = fEdepCavity2 = fTrkSegmCavity = 0.;
0087 fNbEventCavity = 0;
0088
0089
0090
0091 fOldEmean = fOldDose = 0.;
0092
0093
0094
0095 fStepWall = fStepWall2 = fStepCavity = fStepCavity2 = 0.;
0096 fNbStepWall = fNbStepCavity = 0;
0097
0098
0099
0100 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0101 if (analysisManager->IsActive()) {
0102 analysisManager->OpenFile();
0103 }
0104 }
0105
0106
0107
0108 Run::~Run() {}
0109
0110
0111
0112 void Run::EndOfRun()
0113 {
0114 std::ios::fmtflags mode = G4cout.flags();
0115 G4cout.setf(std::ios::fixed, std::ios::floatfield);
0116
0117 if (numberOfEvent == 0) return;
0118
0119
0120
0121 G4ParticleDefinition* particle = fKinematic->GetParticleGun()->GetParticleDefinition();
0122 G4String partName = particle->GetParticleName();
0123 G4double energy = fKinematic->GetParticleGun()->GetParticleEnergy();
0124
0125 G4cout << "\n ======================== run summary ======================\n";
0126
0127 G4int prec = G4cout.precision(3);
0128
0129 G4cout << "\n The run consists of " << numberOfEvent << " " << partName << " of "
0130 << G4BestUnit(energy, "Energy") << " through 2*" << G4BestUnit(fWallThickness, "Length")
0131 << " of " << fMateWall->GetName()
0132 << " (density: " << G4BestUnit(fDensityWall, "Volumic Mass") << ")" << G4endl;
0133
0134 G4cout << "\n the cavity is " << G4BestUnit(fCavityThickness, "Length") << " of "
0135 << fMateCavity->GetName() << " (density: " << G4BestUnit(fDensityCavity, "Volumic Mass")
0136 << "); Mass = " << G4BestUnit(fMassCavity, "Mass") << G4endl;
0137
0138 G4cout << "\n ============================================================\n";
0139
0140
0141
0142 G4cout << "\n Process calls frequency --->";
0143 for (size_t i = 0; i < fProcCounter->size(); i++) {
0144 G4String procName = (*fProcCounter)[i]->GetName();
0145 G4int count = (*fProcCounter)[i]->GetCounter();
0146 G4cout << " " << procName << "= " << count;
0147 }
0148 G4cout << G4endl;
0149
0150
0151
0152 G4EmCalculator emCalculator;
0153 G4cout << "\n Gamma crossSections in wall material :";
0154 G4double sumc = 0.0;
0155 for (size_t i = 0; i < fProcCounter->size(); i++) {
0156 G4String procName = (*fProcCounter)[i]->GetName();
0157 G4double massSigma =
0158 emCalculator.ComputeCrossSectionPerVolume(energy, particle, procName, fMateWall)
0159 / fDensityWall;
0160 if (massSigma > 0.) {
0161 sumc += massSigma;
0162 G4cout << " " << procName << "= " << G4BestUnit(massSigma, "Surface/Mass");
0163 }
0164 }
0165 G4cout << " --> total= " << G4BestUnit(sumc, "Surface/Mass") << G4endl;
0166
0167
0168
0169 if (fNbSec == 0) return;
0170 G4double meanEsecond = fEsecondary / fNbSec, meanEsecond2 = fEsecondary2 / fNbSec;
0171 G4double varianceEsec = meanEsecond2 - meanEsecond * meanEsecond;
0172 G4double dToverT = 0.;
0173 if (varianceEsec > 0.) dToverT = std::sqrt(varianceEsec / fNbSec) / meanEsecond;
0174 G4double csdaRange = emCalculator.GetCSDARange(meanEsecond, G4Electron::Electron(), fMateWall);
0175
0176 G4cout.precision(4);
0177 G4cout << "\n Mean energy of secondary e- = " << G4BestUnit(meanEsecond, "Energy") << " +- "
0178 << 100 * dToverT << " %"
0179 << " (--> range in wall material = " << G4BestUnit(csdaRange, "Length") << ")" << G4endl;
0180
0181
0182
0183 G4double massTransfCoef = sumc * meanEsecond / energy;
0184
0185 G4cout << " Mass_energy_transfer coef: " << G4BestUnit(massTransfCoef, "Surface/Mass") << G4endl;
0186
0187
0188
0189 G4double dedxWall = emCalculator.GetDEDX(meanEsecond, G4Electron::Electron(), fMateWall);
0190 dedxWall /= fDensityWall;
0191 G4double dedxCavity = emCalculator.GetDEDX(meanEsecond, G4Electron::Electron(), fMateCavity);
0192 dedxCavity /= fDensityCavity;
0193
0194 G4cout << "\n StoppingPower in wall = " << G4BestUnit(dedxWall, "Energy*Surface/Mass")
0195 << "\n in cavity = " << G4BestUnit(dedxCavity, "Energy*Surface/Mass")
0196 << G4endl;
0197
0198
0199
0200 G4cout << "\n Charged particle flow in cavity :"
0201 << "\n Enter --> nbParticles = " << fPartFlowCavity[0]
0202 << "\t Energy = " << G4BestUnit(fEnerFlowCavity[0], "Energy")
0203 << "\n Exit --> nbParticles = " << fPartFlowCavity[1]
0204 << "\t Energy = " << G4BestUnit(fEnerFlowCavity[1], "Energy") << G4endl;
0205
0206 if (fPartFlowCavity[0] == 0) return;
0207
0208
0209
0210 G4double rBeam = fWallRadius * (fKinematic->GetBeamRadius());
0211 G4double surfaceBeam = CLHEP::pi * rBeam * rBeam;
0212
0213
0214
0215 if (fNbEventCavity == 0) return;
0216 G4double meanEdep = fEdepCavity / fNbEventCavity;
0217 G4double meanEdep2 = fEdepCavity2 / fNbEventCavity;
0218 G4double varianceEdep = meanEdep2 - meanEdep * meanEdep;
0219 G4double dEoverE = 0.;
0220 if (varianceEdep > 0.) dEoverE = std::sqrt(varianceEdep / fNbEventCavity) / meanEdep;
0221
0222
0223
0224 G4double doseCavity = fEdepCavity / fMassCavity;
0225 G4double doseOverBeam = doseCavity * surfaceBeam / (numberOfEvent * energy);
0226
0227
0228 G4double meantrack = fTrkSegmCavity / fPartFlowCavity[0];
0229
0230 G4cout.precision(4);
0231 G4cout << "\n Total edep in cavity = " << G4BestUnit(fEdepCavity, "Energy") << " +- "
0232 << 100 * dEoverE << " %"
0233 << "\t Total charged trackLength = " << G4BestUnit(fTrkSegmCavity, "Length")
0234 << " (mean value = " << G4BestUnit(meantrack, "Length") << ")"
0235 << "\n Total dose in cavity = " << doseCavity / (MeV / mg) << " MeV/mg"
0236 << "\n Dose/EnergyFluence = " << G4BestUnit(doseOverBeam, "Surface/Mass") << G4endl;
0237
0238
0239
0240 G4double ratio = doseOverBeam / massTransfCoef;
0241 G4double error = ratio * std::sqrt(dEoverE * dEoverE + dToverT * dToverT);
0242
0243 G4cout.precision(5);
0244 G4cout << "\n (Dose/EnergyFluence)/Mass_energy_transfer = " << ratio << " +- " << error << G4endl;
0245
0246
0247
0248 fStepWall /= fNbStepWall;
0249 fStepWall2 /= fNbStepWall;
0250 G4double rms = fStepWall2 - fStepWall * fStepWall;
0251 if (rms > 0.)
0252 rms = std::sqrt(rms);
0253 else
0254 rms = 0.;
0255
0256 G4cout.precision(4);
0257 G4cout << "\n StepSize of ch. tracks in wall = " << G4BestUnit(fStepWall, "Length") << " +- "
0258 << G4BestUnit(rms, "Length") << "\t (nbSteps/track = " << double(fNbStepWall) / fNbSec
0259 << ")";
0260
0261 fStepCavity /= fNbStepCavity;
0262 fStepCavity2 /= fNbStepCavity;
0263 rms = fStepCavity2 - fStepCavity * fStepCavity;
0264 if (rms > 0.)
0265 rms = std::sqrt(rms);
0266 else
0267 rms = 0.;
0268
0269 G4cout << "\n StepSize of ch. tracks in cavity = " << G4BestUnit(fStepCavity, "Length") << " +- "
0270 << G4BestUnit(rms, "Length")
0271 << "\t (nbSteps/track = " << double(fNbStepCavity) / fPartFlowCavity[0] << ")";
0272
0273 G4cout << G4endl;
0274
0275
0276 G4cout.setf(mode, std::ios::floatfield);
0277 G4cout.precision(prec);
0278
0279
0280 while (fProcCounter->size() > 0) {
0281 OneProcessCount* aProcCount = fProcCounter->back();
0282 fProcCounter->pop_back();
0283 delete aProcCount;
0284 }
0285 delete fProcCounter;
0286 }
0287
0288
0289
0290 void Run::SurveyConvergence(G4int NbofEvents)
0291 {
0292 if (NbofEvents == 0) return;
0293
0294
0295
0296 G4double meanEsecond = 0.;
0297 if (fNbSec > 0) meanEsecond = fEsecondary / fNbSec;
0298 G4double rateEmean = 0.;
0299
0300 if (fOldEmean > 0.) rateEmean = 100 * (meanEsecond / fOldEmean - 1.);
0301 fOldEmean = meanEsecond;
0302
0303
0304
0305 G4double rBeam = fWallRadius * (fKinematic->GetBeamRadius());
0306 G4double surfaceBeam = CLHEP::pi * rBeam * rBeam;
0307 G4double beamEnergy = fKinematic->GetParticleGun()->GetParticleEnergy();
0308
0309
0310
0311 G4double doseCavity = fEdepCavity / fMassCavity;
0312 G4double doseOverBeam = doseCavity * surfaceBeam / (NbofEvents * beamEnergy);
0313 G4double rateDose = 0.;
0314
0315 if (fOldDose > 0.) rateDose = 100 * (doseOverBeam / fOldDose - 1.);
0316 fOldDose = doseOverBeam;
0317
0318 std::ios::fmtflags mode = G4cout.flags();
0319 G4cout.setf(std::ios::fixed, std::ios::floatfield);
0320 G4int prec = G4cout.precision(3);
0321
0322 G4cout << " ---> NbofEvents= " << NbofEvents << " NbOfelectr= " << fNbSec
0323 << " Tkin= " << G4BestUnit(meanEsecond, "Energy") << " (" << rateEmean << " %)"
0324 << " NbOfelec in cav= " << fPartFlowCavity[0]
0325 << " Dose/EnFluence= " << G4BestUnit(doseOverBeam, "Surface/Mass") << " (" << rateDose
0326 << " %) \n"
0327 << G4endl;
0328
0329
0330 G4cout.setf(mode, std::ios::floatfield);
0331 G4cout.precision(prec);
0332 }
0333
0334
0335
0336 void Run::Merge(const G4Run* run)
0337 {
0338 const Run* localRun = static_cast<const Run*>(run);
0339
0340 fPartFlowCavity[0] += localRun->fPartFlowCavity[0];
0341 fPartFlowCavity[1] += localRun->fPartFlowCavity[1];
0342 fEnerFlowCavity[0] += localRun->fEnerFlowCavity[0];
0343 fEnerFlowCavity[1] += localRun->fEnerFlowCavity[1];
0344 fEdepCavity += localRun->fEdepCavity;
0345 fEdepCavity2 += localRun->fEdepCavity2;
0346 fTrkSegmCavity += localRun->fTrkSegmCavity;
0347 fNbEventCavity += localRun->fNbEventCavity;
0348
0349 fStepWall += localRun->fStepWall;
0350 fStepWall2 += localRun->fStepWall2;
0351 fStepCavity += localRun->fStepCavity;
0352 fStepCavity2 += localRun->fStepCavity2;
0353 fNbStepWall += localRun->fNbStepWall;
0354 fNbStepCavity += localRun->fNbStepCavity;
0355
0356 fEsecondary += localRun->fEsecondary;
0357 fEsecondary2 += localRun->fEsecondary2;
0358
0359 fNbSec += localRun->fNbSec;
0360
0361
0362
0363
0364
0365 std::vector<OneProcessCount*>::iterator it;
0366 for (it = localRun->fProcCounter->begin(); it != localRun->fProcCounter->end(); it++) {
0367 OneProcessCount* process = *it;
0368 for (G4int i = 0; i < process->GetCounter(); i++)
0369 this->CountProcesses(process->GetName());
0370 }
0371
0372 G4Run::Merge(run);
0373 }
0374
0375
0376
0377 void Run::CountProcesses(G4String procName)
0378 {
0379
0380 size_t nbProc = fProcCounter->size();
0381 size_t i = 0;
0382 while ((i < nbProc) && ((*fProcCounter)[i]->GetName() != procName))
0383 i++;
0384 if (i == nbProc) fProcCounter->push_back(new OneProcessCount(procName));
0385
0386 (*fProcCounter)[i]->Count();
0387 }
0388
0389
