File indexing completed on 2025-10-31 08:24:00
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, bool isMaster)
0052 : G4Run(),
0053 fDetector(det),
0054 fKinematic(kin),
0055 fProcCounter(0),
0056 fEdepCavity(0.),
0057 fEdepCavity2(0.),
0058 fTrkSegmCavity(0.),
0059 fNbEventCavity(0),
0060 fStepWall(0.),
0061 fStepWall2(0.),
0062 fStepCavity(0.),
0063 fStepCavity2(0.),
0064 fNbStepWall(0),
0065 fNbStepCavity(0),
0066 fEnergyGun(0.),
0067 fMassWall(0.),
0068 fMassCavity(0.),
0069 fIsMaster(isMaster)
0070 {
0071
0072
0073 G4ParticleDefinition* particleGun = fKinematic->GetParticleGun()->GetParticleDefinition();
0074 G4String partName = particleGun->GetParticleName();
0075 fEnergyGun = fKinematic->GetParticleGun()->GetParticleEnergy();
0076
0077
0078
0079 G4double cavityThickness = fDetector->GetCavityThickness();
0080 G4Material* mateCavity = fDetector->GetCavityMaterial();
0081 G4double densityCavity = mateCavity->GetDensity();
0082 fMassCavity = cavityThickness * densityCavity;
0083
0084 G4double wallThickness = fDetector->GetWallThickness();
0085 G4Material* mateWall = fDetector->GetWallMaterial();
0086 G4double densityWall = mateWall->GetDensity();
0087
0088 G4EmCalculator emCal;
0089 G4double RangeWall = emCal.GetCSDARange(fEnergyGun, particleGun, mateWall);
0090 G4double factor = 1.2;
0091 G4double effWallThick = factor * RangeWall;
0092 if ((effWallThick > wallThickness) || (effWallThick <= 0.)) effWallThick = wallThickness;
0093 fMassWall = 2 * effWallThick * densityWall;
0094
0095 G4double massTotal = fMassWall + fMassCavity;
0096 G4double fMassWallRatio = fMassWall / massTotal;
0097 fKinematic->RunInitialisation(effWallThick, fMassWallRatio);
0098
0099 G4double massRatio = fMassCavity / fMassWall;
0100
0101
0102
0103 G4double worldRadius = fDetector->GetWallRadius();
0104 G4double RangeCavity = emCal.GetCSDARange(fEnergyGun, particleGun, mateCavity);
0105
0106
0107
0108 std::ios::fmtflags mode = G4cout.flags();
0109 G4cout.setf(std::ios::fixed, std::ios::floatfield);
0110 G4int prec = G4cout.precision(3);
0111
0112 G4cout << "\n ===================== run conditions =====================\n";
0113
0114 G4cout << "\n The run will be " << numberOfEvent << " " << partName << " of "
0115 << G4BestUnit(fEnergyGun, "Energy") << " through 2*" << G4BestUnit(effWallThick, "Length")
0116 << " of " << mateWall->GetName()
0117 << " (density: " << G4BestUnit(densityWall, "Volumic Mass")
0118 << "); Mass/cm2 = " << G4BestUnit(fMassWall * cm2, "Mass")
0119 << "\n csdaRange: " << G4BestUnit(RangeWall, "Length") << G4endl;
0120
0121 G4cout << "\n the cavity is " << G4BestUnit(cavityThickness, "Length") << " of "
0122 << mateCavity->GetName() << " (density: " << G4BestUnit(densityCavity, "Volumic Mass")
0123 << "); Mass/cm2 = " << G4BestUnit(fMassCavity * cm2, "Mass")
0124 << " --> massRatio = " << std::setprecision(6) << massRatio << G4endl;
0125
0126 G4cout.precision(3);
0127 G4cout << " Wall radius: " << G4BestUnit(worldRadius, "Length")
0128 << "; range in cavity: " << G4BestUnit(RangeCavity, "Length") << G4endl;
0129
0130 G4cout << "\n ==========================================================\n";
0131
0132
0133
0134 G4double dedxWall = emCal.GetDEDX(fEnergyGun, G4Electron::Electron(), mateWall);
0135 dedxWall /= densityWall;
0136 G4double dedxCavity = emCal.GetDEDX(fEnergyGun, G4Electron::Electron(), mateCavity);
0137 dedxCavity /= densityCavity;
0138
0139 G4cout << std::setprecision(4)
0140 << "\n StoppingPower in wall = " << G4BestUnit(dedxWall, "Energy*Surface/Mass")
0141 << "\n in cavity = " << G4BestUnit(dedxCavity, "Energy*Surface/Mass")
0142 << G4endl;
0143
0144
0145
0146 fProcCounter = new ProcessesCount;
0147
0148
0149
0150 fPartFlowCavity[0] = fPartFlowCavity[1] = 0;
0151 fEnerFlowCavity[0] = fEnerFlowCavity[1] = 0.;
0152
0153
0154
0155 fEdepCavity = fEdepCavity2 = fTrkSegmCavity = 0.;
0156 fNbEventCavity = 0;
0157
0158
0159
0160 fStepWall = fStepWall2 = fStepCavity = fStepCavity2 = 0.;
0161 fNbStepWall = fNbStepCavity = 0;
0162
0163
0164 G4cout.setf(mode, std::ios::floatfield);
0165 G4cout.precision(prec);
0166
0167
0168
0169 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0170 if (analysisManager->IsActive()) {
0171 analysisManager->OpenFile();
0172 }
0173 }
0174
0175
0176
0177 Run::~Run() {}
0178
0179
0180
0181 void Run::CountProcesses(G4String procName)
0182 {
0183
0184 size_t nbProc = fProcCounter->size();
0185 size_t i = 0;
0186 while ((i < nbProc) && ((*fProcCounter)[i]->GetName() != procName))
0187 i++;
0188 if (i == nbProc) fProcCounter->push_back(new OneProcessCount(procName));
0189
0190 (*fProcCounter)[i]->Count();
0191 }
0192
0193
0194
0195 void Run::SurveyConvergence(G4int NbofEvents)
0196 {
0197 if (NbofEvents == 0) return;
0198
0199
0200
0201 G4int Nwall = fKinematic->GetWallCount();
0202 G4int Ncavity = fKinematic->GetCavityCount();
0203 G4double Iwall = Nwall / fMassWall;
0204 G4double Icavity = Ncavity / fMassCavity;
0205 G4double Iratio = Icavity / Iwall;
0206 G4double Itot = NbofEvents / (fMassWall + fMassCavity);
0207 G4double energyFluence = fEnergyGun * Itot;
0208
0209
0210
0211 G4double doseCavity = fEdepCavity / fMassCavity;
0212 G4double ratio = doseCavity / energyFluence;
0213 G4double err = 100 * (ratio - 1.);
0214
0215 std::ios::fmtflags mode = G4cout.flags();
0216 G4cout.setf(std::ios::fixed, std::ios::floatfield);
0217 G4int prec = G4cout.precision(5);
0218
0219 G4cout << "--->evntNb= " << NbofEvents << " Nwall= " << Nwall << " Ncav= " << Ncavity
0220 << " Ic/Iw= " << Iratio << " Ne-_cav= " << fPartFlowCavity[0]
0221 << " doseCavity/Ebeam= " << ratio << " (100*(ratio-1) = " << err << " %) \n"
0222 << G4endl;
0223
0224
0225 G4cout.setf(mode, std::ios::floatfield);
0226 G4cout.precision(prec);
0227 }
0228
0229
0230
0231 void Run::EndOfRun()
0232 {
0233
0234 std::ios::fmtflags mode = G4cout.flags();
0235 G4cout.setf(std::ios::fixed, std::ios::floatfield);
0236 G4int prec = G4cout.precision(3);
0237
0238 if (numberOfEvent == 0) return;
0239
0240
0241
0242 G4cout << "\n Process calls frequency --->";
0243 for (size_t i = 0; i < fProcCounter->size(); i++) {
0244 G4String procName = (*fProcCounter)[i]->GetName();
0245 G4int count = (*fProcCounter)[i]->GetCounter();
0246 G4cout << " " << procName << "= " << count;
0247 }
0248 G4cout << G4endl;
0249
0250
0251
0252 G4cout << "\n Charged particle flow in cavity :"
0253 << "\n Enter --> nbParticles = " << fPartFlowCavity[0]
0254 << "\t Energy = " << G4BestUnit(fEnerFlowCavity[0], "Energy")
0255 << "\n Exit --> nbParticles = " << fPartFlowCavity[1]
0256 << "\t Energy = " << G4BestUnit(fEnerFlowCavity[1], "Energy") << G4endl;
0257
0258 if (fPartFlowCavity[0] == 0) return;
0259
0260 G4int Nwall = fKinematic->GetWallCount();
0261 G4int Ncavity = fKinematic->GetCavityCount();
0262
0263 G4double Iwall = Nwall / fMassWall;
0264 G4double Icavity = Ncavity / fMassCavity;
0265 G4double Iratio = Icavity / Iwall;
0266 G4double Itot = numberOfEvent / (fMassWall + fMassCavity);
0267 G4double energyFluence = fEnergyGun * Itot;
0268
0269 G4cout.precision(5);
0270 G4cout << "\n beamFluence in wall = " << Nwall << "\t in cavity = " << Ncavity
0271 << "\t Icav/Iwall = " << Iratio
0272 << "\t energyFluence = " << energyFluence / (MeV * cm2 / mg) << " MeV*cm2/mg" << G4endl;
0273
0274
0275
0276 if (fNbEventCavity == 0) return;
0277 G4double meanEdep = fEdepCavity / fNbEventCavity;
0278 G4double meanEdep2 = fEdepCavity2 / fNbEventCavity;
0279 G4double varianceEdep = meanEdep2 - meanEdep * meanEdep;
0280 G4double dEoverE = 0.;
0281 if (varianceEdep > 0.) dEoverE = std::sqrt(varianceEdep / fNbEventCavity) / meanEdep;
0282
0283
0284
0285 G4double doseCavity = fEdepCavity / fMassCavity;
0286
0287 G4double ratio = doseCavity / energyFluence, error = ratio * dEoverE;
0288
0289 G4cout << "\n Total edep in cavity = " << G4BestUnit(fEdepCavity, "Energy") << " +- "
0290 << 100 * dEoverE << " %"
0291 << "\n Total dose in cavity = " << doseCavity / (MeV * cm2 / mg) << " MeV*cm2/mg"
0292 << " +- " << 100 * dEoverE << " %"
0293 << "\n\n DoseCavity/EnergyFluence = " << ratio << " +- " << error << G4endl;
0294
0295
0296 G4double meantrack = fTrkSegmCavity / fPartFlowCavity[0];
0297
0298 G4cout.precision(4);
0299 G4cout << "\n Total charged trackLength in cavity = " << G4BestUnit(fTrkSegmCavity, "Length")
0300 << " (mean value = " << G4BestUnit(meantrack, "Length") << ")" << G4endl;
0301
0302
0303
0304 fStepWall /= fNbStepWall;
0305 fStepWall2 /= fNbStepWall;
0306 G4double rms = fStepWall2 - fStepWall * fStepWall;
0307 if (rms > 0.)
0308 rms = std::sqrt(rms);
0309 else
0310 rms = 0.;
0311 G4double nbTrackWall = fKinematic->GetWallCount();
0312
0313 G4cout << "\n StepSize of ch. tracks in wall = " << G4BestUnit(fStepWall, "Length") << " +- "
0314 << G4BestUnit(rms, "Length") << "\t (nbSteps/track = " << double(fNbStepWall) / nbTrackWall
0315 << ")";
0316
0317 fStepCavity /= fNbStepCavity;
0318 fStepCavity2 /= fNbStepCavity;
0319 rms = fStepCavity2 - fStepCavity * fStepCavity;
0320 if (rms > 0.)
0321 rms = std::sqrt(rms);
0322 else
0323 rms = 0.;
0324
0325 G4cout << "\n StepSize of ch. tracks in cavity = " << G4BestUnit(fStepCavity, "Length") << " +- "
0326 << G4BestUnit(rms, "Length")
0327 << "\t (nbSteps/track = " << double(fNbStepCavity) / fPartFlowCavity[0] << ")";
0328
0329 G4cout << G4endl;
0330
0331
0332 G4cout.setf(mode, std::ios::floatfield);
0333 G4cout.precision(prec);
0334
0335
0336 while (fProcCounter->size() > 0) {
0337 OneProcessCount* aProcCount = fProcCounter->back();
0338 fProcCounter->pop_back();
0339 delete aProcCount;
0340 }
0341 delete fProcCounter;
0342
0343
0344 CLHEP::HepRandom::showEngineStatus();
0345 }
0346
0347
0348
0349 void Run::Merge(const G4Run* run)
0350 {
0351 const Run* localRun = static_cast<const Run*>(run);
0352
0353
0354 fPartFlowCavity[0] += localRun->fPartFlowCavity[0];
0355 fPartFlowCavity[1] += localRun->fPartFlowCavity[1];
0356 fEnerFlowCavity[0] += localRun->fEnerFlowCavity[0];
0357 fEnerFlowCavity[1] += localRun->fEnerFlowCavity[1];
0358 fEdepCavity += localRun->fEdepCavity;
0359 fEdepCavity2 += localRun->fEdepCavity2;
0360 fTrkSegmCavity += localRun->fTrkSegmCavity;
0361 fNbEventCavity += localRun->fNbEventCavity;
0362 fStepWall += localRun->fStepWall;
0363 fStepWall2 += localRun->fStepWall2;
0364 fStepCavity += localRun->fStepCavity;
0365 fStepCavity2 += localRun->fStepCavity2;
0366 fNbStepWall += localRun->fNbStepWall;
0367 fNbStepCavity += localRun->fNbStepCavity;
0368
0369
0370 fKinematic->AddWallCount(localRun->fKinematic->GetWallCount());
0371 fKinematic->AddCavityCount(localRun->fKinematic->GetCavityCount());
0372
0373
0374 std::vector<OneProcessCount*>::iterator it;
0375 for (it = localRun->fProcCounter->begin(); it != localRun->fProcCounter->end(); it++) {
0376 OneProcessCount* process = *it;
0377 for (G4int i = 0; i < process->GetCounter(); i++)
0378 this->CountProcesses(process->GetName());
0379 }
0380
0381 G4Run::Merge(run);
0382 }
0383
0384
