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