File indexing completed on 2026-09-20 08:29:46
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 "HistoManager.hh"
0033 #include "PrimaryGeneratorAction.hh"
0034
0035 #include "G4EmCalculator.hh"
0036 #include "G4SystemOfUnits.hh"
0037 #include "G4UnitsTable.hh"
0038
0039 #include <iomanip>
0040
0041
0042
0043 Run::Run(DetectorConstruction* det) : fDetector(det) {}
0044
0045
0046
0047 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
0048 {
0049 fParticle = particle;
0050 fEkin = energy;
0051
0052
0053 fMscThetaCentral = 3 * ComputeMscHighland();
0054 }
0055
0056
0057
0058 void Run::Merge(const G4Run* run)
0059 {
0060 const Run* localRun = static_cast<const Run*>(run);
0061
0062
0063 fParticle = localRun->fParticle;
0064 fEkin = localRun->fEkin;
0065
0066 fMscThetaCentral = localRun->fMscThetaCentral;
0067
0068
0069
0070 fEnergyDeposit += localRun->fEnergyDeposit;
0071 fEnergyDeposit2 += localRun->fEnergyDeposit2;
0072 fTrakLenCharged += localRun->fTrakLenCharged;
0073 fTrakLenCharged2 += localRun->fTrakLenCharged2;
0074 fTrakLenNeutral += localRun->fTrakLenNeutral;
0075 fTrakLenNeutral2 += localRun->fTrakLenNeutral2;
0076 fNbStepsCharged += localRun->fNbStepsCharged;
0077 fNbStepsCharged2 += localRun->fNbStepsCharged2;
0078 fNbStepsNeutral += localRun->fNbStepsNeutral;
0079 fNbStepsNeutral2 += localRun->fNbStepsNeutral2;
0080 fMscProjecTheta += localRun->fMscProjecTheta;
0081 fMscProjecTheta2 += localRun->fMscProjecTheta2;
0082
0083 fTypes[0] += localRun->fTypes[0];
0084 fTypes[1] += localRun->fTypes[1];
0085 fTypes[2] += localRun->fTypes[2];
0086 fTypes[3] += localRun->fTypes[3];
0087
0088 fNbGamma += localRun->fNbGamma;
0089 fNbElect += localRun->fNbElect;
0090 fNbPosit += localRun->fNbPosit;
0091
0092 fTransmit[0] += localRun->fTransmit[0];
0093 fTransmit[1] += localRun->fTransmit[1];
0094 fReflect[0] += localRun->fReflect[0];
0095 fReflect[1] += localRun->fReflect[1];
0096
0097 fMscEntryCentral += localRun->fMscEntryCentral;
0098
0099 fEnergyLeak[0] += localRun->fEnergyLeak[0];
0100 fEnergyLeak[1] += localRun->fEnergyLeak[1];
0101 fEnergyLeak2[0] += localRun->fEnergyLeak2[0];
0102 fEnergyLeak2[1] += localRun->fEnergyLeak2[1];
0103
0104 G4Run::Merge(run);
0105 }
0106
0107
0108
0109 void Run::EndOfRun()
0110 {
0111
0112
0113 G4int TotNbofEvents = numberOfEvent;
0114 if (TotNbofEvents == 0) return;
0115
0116 G4double EnergyBalance = fEnergyDeposit + fEnergyLeak[0] + fEnergyLeak[1];
0117 EnergyBalance /= TotNbofEvents;
0118
0119 fEnergyDeposit /= TotNbofEvents;
0120 fEnergyDeposit2 /= TotNbofEvents;
0121 G4double rmsEdep = fEnergyDeposit2 - fEnergyDeposit * fEnergyDeposit;
0122 if (rmsEdep > 0.)
0123 rmsEdep = std::sqrt(rmsEdep / TotNbofEvents);
0124 else
0125 rmsEdep = 0.;
0126
0127 fTrakLenCharged /= TotNbofEvents;
0128 fTrakLenCharged2 /= TotNbofEvents;
0129 G4double rmsTLCh = fTrakLenCharged2 - fTrakLenCharged * fTrakLenCharged;
0130 if (rmsTLCh > 0.)
0131 rmsTLCh = std::sqrt(rmsTLCh / TotNbofEvents);
0132 else
0133 rmsTLCh = 0.;
0134
0135 fTrakLenNeutral /= TotNbofEvents;
0136 fTrakLenNeutral2 /= TotNbofEvents;
0137 G4double rmsTLNe = fTrakLenNeutral2 - fTrakLenNeutral * fTrakLenNeutral;
0138 if (rmsTLNe > 0.)
0139 rmsTLNe = std::sqrt(rmsTLNe / TotNbofEvents);
0140 else
0141 rmsTLNe = 0.;
0142
0143 fNbStepsCharged /= TotNbofEvents;
0144 fNbStepsCharged2 /= TotNbofEvents;
0145 G4double rmsStCh = fNbStepsCharged2 - fNbStepsCharged * fNbStepsCharged;
0146 if (rmsStCh > 0.)
0147 rmsStCh = std::sqrt(rmsStCh / TotNbofEvents);
0148 else
0149 rmsStCh = 0.;
0150
0151 fNbStepsNeutral /= TotNbofEvents;
0152 fNbStepsNeutral2 /= TotNbofEvents;
0153 G4double rmsStNe = fNbStepsNeutral2 - fNbStepsNeutral * fNbStepsNeutral;
0154 if (rmsStNe > 0.)
0155 rmsStNe = std::sqrt(rmsStNe / TotNbofEvents);
0156 else
0157 rmsStNe = 0.;
0158
0159 G4double Gamma = (G4double)fNbGamma / TotNbofEvents;
0160 G4double Elect = (G4double)fNbElect / TotNbofEvents;
0161 G4double Posit = (G4double)fNbPosit / TotNbofEvents;
0162
0163 G4double transmit[2];
0164 transmit[0] = 100. * fTransmit[0] / TotNbofEvents;
0165 transmit[1] = 100. * fTransmit[1] / TotNbofEvents;
0166
0167 G4double reflect[2];
0168 reflect[0] = 100. * fReflect[0] / TotNbofEvents;
0169 reflect[1] = 100. * fReflect[1] / TotNbofEvents;
0170
0171 G4double rmsMsc = 0., tailMsc = 0.;
0172 if (fMscEntryCentral > 0) {
0173 fMscProjecTheta /= fMscEntryCentral;
0174 fMscProjecTheta2 /= fMscEntryCentral;
0175 rmsMsc = fMscProjecTheta2 - fMscProjecTheta * fMscProjecTheta;
0176 if (rmsMsc > 0.) {
0177 rmsMsc = std::sqrt(rmsMsc);
0178 }
0179 if (fTransmit[1] > 0.0) {
0180 tailMsc = 100. - (100. * fMscEntryCentral) / (2 * fTransmit[1]);
0181 }
0182 }
0183
0184 fEnergyLeak[0] /= TotNbofEvents;
0185 fEnergyLeak2[0] /= TotNbofEvents;
0186 G4double rmsEl0 = fEnergyLeak2[0] - fEnergyLeak[0] * fEnergyLeak[0];
0187 if (rmsEl0 > 0.)
0188 rmsEl0 = std::sqrt(rmsEl0 / TotNbofEvents);
0189 else
0190 rmsEl0 = 0.;
0191
0192 fEnergyLeak[1] /= TotNbofEvents;
0193 fEnergyLeak2[1] /= TotNbofEvents;
0194 G4double rmsEl1 = fEnergyLeak2[1] - fEnergyLeak[1] * fEnergyLeak[1];
0195 if (rmsEl1 > 0.)
0196 rmsEl1 = std::sqrt(rmsEl1 / TotNbofEvents);
0197 else
0198 rmsEl1 = 0.;
0199
0200
0201
0202 const G4Material* material = fDetector->GetAbsorberMaterial();
0203 G4double length = fDetector->GetAbsorberThickness();
0204 G4double density = material->GetDensity();
0205 G4String partName = fParticle->GetParticleName();
0206
0207 G4EmCalculator emCalculator;
0208 G4double dEdxTable = 0., dEdxFull = 0.;
0209 if (fParticle->GetPDGCharge() != 0.) {
0210 dEdxTable = emCalculator.GetDEDX(fEkin, fParticle, material);
0211 dEdxFull = emCalculator.ComputeTotalDEDX(fEkin, fParticle, material);
0212 }
0213 G4double stopTable = dEdxTable / density;
0214 G4double stopFull = dEdxFull / density;
0215
0216
0217
0218 G4double meandEdx = fEnergyDeposit / length;
0219 G4double stopPower = meandEdx / density;
0220
0221 G4cout << "\n ======================== run summary ======================\n";
0222
0223 G4int prec = G4cout.precision(3);
0224
0225 G4cout << "\n The run was " << TotNbofEvents << " " << partName << " of "
0226 << G4BestUnit(fEkin, "Energy") << " through " << G4BestUnit(length, "Length") << " of "
0227 << material->GetName() << " (density: " << G4BestUnit(density, "Volumic Mass") << ")"
0228 << G4endl;
0229
0230 G4cout.precision(4);
0231
0232 G4cout << "\n Total energy deposit in absorber per event = "
0233 << G4BestUnit(fEnergyDeposit, "Energy") << " +- " << G4BestUnit(rmsEdep, "Energy")
0234 << G4endl;
0235
0236 G4cout << "\n -----> Mean dE/dx = " << meandEdx / (MeV / cm) << " MeV/cm"
0237 << "\t(" << stopPower / (MeV * cm2 / g) << " MeV*cm2/g)" << G4endl;
0238
0239 G4cout << "\n From formulas :" << G4endl;
0240 G4cout << " restricted dEdx = " << dEdxTable / (MeV / cm) << " MeV/cm"
0241 << "\t(" << stopTable / (MeV * cm2 / g) << " MeV*cm2/g)" << G4endl;
0242
0243 G4cout << " full dEdx = " << dEdxFull / (MeV / cm) << " MeV/cm"
0244 << "\t(" << stopFull / (MeV * cm2 / g) << " MeV*cm2/g)" << G4endl;
0245
0246 G4cout << "\n Leakage : primary = " << G4BestUnit(fEnergyLeak[0], "Energy") << " +- "
0247 << G4BestUnit(rmsEl0, "Energy")
0248 << " secondaries = " << G4BestUnit(fEnergyLeak[1], "Energy") << " +- "
0249 << G4BestUnit(rmsEl1, "Energy") << G4endl;
0250
0251 G4cout << " Energy balance : edep + eleak = " << G4BestUnit(EnergyBalance, "Energy") << G4endl;
0252
0253 G4cout << "\n Total track length (charged) in absorber per event = "
0254 << G4BestUnit(fTrakLenCharged, "Length") << " +- " << G4BestUnit(rmsTLCh, "Length")
0255 << G4endl;
0256
0257 G4cout << " Total track length (neutral) in absorber per event = "
0258 << G4BestUnit(fTrakLenNeutral, "Length") << " +- " << G4BestUnit(rmsTLNe, "Length")
0259 << G4endl;
0260
0261 G4cout << "\n Number of steps (charged) in absorber per event = " << fNbStepsCharged << " +- "
0262 << rmsStCh << G4endl;
0263
0264 G4cout << " Number of steps (neutral) in absorber per event = " << fNbStepsNeutral << " +- "
0265 << rmsStNe << G4endl;
0266
0267 G4cout << "\n Number of secondaries per event : Gammas = " << Gamma << "; electrons = " << Elect
0268 << "; positrons = " << Posit << G4endl;
0269
0270 G4cout << "\n Number of events with the primary particle transmitted = " << transmit[1] << " %"
0271 << G4endl;
0272
0273 G4cout << " Number of events with at least 1 particle transmitted "
0274 << "(same charge as primary) = " << transmit[0] << " %" << G4endl;
0275
0276 G4cout << "\n Number of events with the primary particle reflected = " << reflect[1] << " %"
0277 << G4endl;
0278
0279 G4cout << " Number of events with at least 1 particle reflected "
0280 << "(same charge as primary) = " << reflect[0] << " %" << G4endl;
0281
0282
0283
0284 G4cout << "\n MultipleScattering:"
0285 << "\n rms proj angle of transmit primary particle = " << rmsMsc / mrad
0286 << " mrad (central part only)" << G4endl;
0287
0288 G4cout << " computed theta0 (Highland formula) = " << ComputeMscHighland() / mrad
0289 << " mrad" << G4endl;
0290
0291 G4cout << " central part defined as +- " << fMscThetaCentral / mrad << " mrad; "
0292 << " Tail ratio = " << tailMsc << " %" << G4endl;
0293
0294
0295
0296 G4cout << "\n Gamma process counts:" << G4endl;
0297 G4cout << " Photoeffect " << fTypes[0] << G4endl;
0298 G4cout << " Compton " << fTypes[1] << G4endl;
0299 G4cout << " Conversion " << fTypes[2] << G4endl;
0300 G4cout << " Rayleigh " << fTypes[3] << G4endl;
0301
0302
0303
0304 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0305
0306 G4int ih = 1;
0307 G4double binWidth = analysisManager->GetH1Width(ih);
0308 G4double fac = 1. / (TotNbofEvents * binWidth);
0309 analysisManager->ScaleH1(ih, fac);
0310
0311 ih = 10;
0312 binWidth = analysisManager->GetH1Width(ih);
0313 fac = 1. / (TotNbofEvents * binWidth);
0314 analysisManager->ScaleH1(ih, fac);
0315
0316 ih = 12;
0317 analysisManager->ScaleH1(ih, 1. / TotNbofEvents);
0318
0319
0320 G4cout.precision(prec);
0321 }
0322
0323
0324
0325 G4double Run::ComputeMscHighland()
0326 {
0327
0328
0329
0330
0331 G4double t =
0332 (fDetector->GetAbsorberThickness()) / (fDetector->GetAbsorberMaterial()->GetRadlen());
0333 if (t < DBL_MIN) return 0.;
0334
0335 G4double T = fEkin;
0336 G4double M = fParticle->GetPDGMass();
0337 G4double z = std::abs(fParticle->GetPDGCharge() / eplus);
0338
0339 G4double bpc = T * (T + 2 * M) / (T + M);
0340 G4double teta0 = 13.6 * MeV * z * std::sqrt(t) * (1. + 0.038 * std::log(t)) / bpc;
0341 return teta0;
0342 }
0343
0344
