File indexing completed on 2025-12-16 09:30:22
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 "EmAcceptance.hh"
0036 #include "PrimaryGeneratorAction.hh"
0037
0038 #include "G4Run.hh"
0039 #include "G4SystemOfUnits.hh"
0040 #include "G4UnitsTable.hh"
0041
0042 #include <iomanip>
0043
0044
0045
0046 Run::Run(DetectorConstruction* det, PrimaryGeneratorAction* kin) : fDet(det), fKin(kin)
0047 {
0048 Reset();
0049 }
0050
0051
0052
0053 void Run::Reset()
0054 {
0055 f_nLbin = fDet->GetnLtot();
0056 f_dEdL.resize(f_nLbin);
0057 fSumELongit.resize(f_nLbin);
0058 fSumELongitCumul.resize(f_nLbin);
0059 fSumE2Longit.resize(f_nLbin);
0060 fSumE2LongitCumul.resize(f_nLbin);
0061
0062 f_nRbin = fDet->GetnRtot();
0063 f_dEdR.resize(f_nRbin);
0064 fSumERadial.resize(f_nRbin);
0065 fSumERadialCumul.resize(f_nRbin);
0066 fSumE2Radial.resize(f_nRbin);
0067 fSumE2RadialCumul.resize(f_nRbin);
0068
0069 fChargedStep = 0.0;
0070 fNeutralStep = 0.0;
0071
0072 fVerbose = 0;
0073
0074
0075
0076 for (G4int i = 0; i < f_nLbin; ++i) {
0077 fSumELongit[i] = fSumE2Longit[i] = fSumELongitCumul[i] = fSumE2LongitCumul[i] = 0.;
0078 }
0079 for (G4int j = 0; j < f_nRbin; ++j) {
0080 fSumERadial[j] = fSumE2Radial[j] = fSumERadialCumul[j] = fSumE2RadialCumul[j] = 0.;
0081 }
0082
0083 fSumChargTrLength = fSum2ChargTrLength = fSumNeutrTrLength = fSum2NeutrTrLength = 0.;
0084 }
0085
0086
0087
0088 void Run::InitializePerEvent()
0089 {
0090
0091 for (G4int i = 0; i < f_nLbin; ++i) {
0092 f_dEdL[i] = 0.;
0093 }
0094
0095 for (G4int j = 0; j < f_nRbin; ++j) {
0096 f_dEdR[j] = 0.;
0097 }
0098
0099
0100 fChargTrLength = fNeutrTrLength = 0.;
0101 }
0102
0103
0104
0105 void Run::FillPerEvent()
0106 {
0107
0108
0109 G4double dLCumul = 0.;
0110 for (G4int i = 0; i < f_nLbin; ++i) {
0111 fSumELongit[i] += f_dEdL[i];
0112 fSumE2Longit[i] += f_dEdL[i] * f_dEdL[i];
0113 dLCumul += f_dEdL[i];
0114 fSumELongitCumul[i] += dLCumul;
0115 fSumE2LongitCumul[i] += dLCumul * dLCumul;
0116 }
0117
0118 G4double dRCumul = 0.;
0119 for (G4int j = 0; j < f_nRbin; j++) {
0120 fSumERadial[j] += f_dEdR[j];
0121 fSumE2Radial[j] += f_dEdR[j] * f_dEdR[j];
0122 dRCumul += f_dEdR[j];
0123 fSumERadialCumul[j] += dRCumul;
0124 fSumE2RadialCumul[j] += dRCumul * dRCumul;
0125 }
0126
0127 fSumChargTrLength += fChargTrLength;
0128 fSum2ChargTrLength += fChargTrLength * fChargTrLength;
0129 fSumNeutrTrLength += fNeutrTrLength;
0130 fSum2NeutrTrLength += fNeutrTrLength * fNeutrTrLength;
0131
0132
0133
0134
0135 G4double Ekin = fKin->GetParticleGun()->GetParticleEnergy();
0136 G4double mass = fKin->GetParticleGun()->GetParticleDefinition()->GetPDGMass();
0137 G4double radl = fDet->GetMaterial()->GetRadlen();
0138
0139 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0140 analysisManager->FillH1(1, 100. * dLCumul / (Ekin + mass));
0141 analysisManager->FillH1(2, fChargTrLength / radl);
0142 analysisManager->FillH1(3, fNeutrTrLength / radl);
0143
0144
0145 G4double norm = 100. / (Ekin + mass);
0146 G4double dLradl = fDet->GetdLradl();
0147 for (G4int i = 0; i < f_nLbin; ++i) {
0148 G4double bin = (i + 0.5) * dLradl;
0149 analysisManager->FillP1(0, bin, norm * f_dEdL[i] / dLradl);
0150 }
0151 G4double dRradl = fDet->GetdRradl();
0152 for (G4int j = 0; j < f_nRbin; ++j) {
0153 G4double bin = (j + 0.5) * dRradl;
0154 analysisManager->FillP1(1, bin, norm * f_dEdR[j] / dRradl);
0155 }
0156 }
0157
0158
0159
0160 void Run::Merge(const G4Run* run)
0161 {
0162 const Run* localRun = static_cast<const Run*>(run);
0163
0164 fChargedStep += localRun->fChargedStep;
0165 fNeutralStep += localRun->fNeutralStep;
0166
0167 for (G4int i = 0; i < f_nLbin; ++i) {
0168 fSumELongit[i] += localRun->fSumELongit[i];
0169 fSumE2Longit[i] += localRun->fSumE2Longit[i];
0170 fSumELongitCumul[i] += localRun->fSumELongitCumul[i];
0171 fSumE2LongitCumul[i] += localRun->fSumE2LongitCumul[i];
0172 }
0173
0174 for (G4int j = 0; j < f_nRbin; ++j) {
0175 fSumERadial[j] += localRun->fSumERadial[j];
0176 fSumE2Radial[j] += localRun->fSumE2Radial[j];
0177 fSumERadialCumul[j] += localRun->fSumERadialCumul[j];
0178 fSumE2RadialCumul[j] += localRun->fSumE2RadialCumul[j];
0179 }
0180
0181 fSumChargTrLength += localRun->fSumChargTrLength;
0182 fSum2ChargTrLength += localRun->fSum2ChargTrLength;
0183 fSumNeutrTrLength += localRun->fSumNeutrTrLength;
0184 fSum2NeutrTrLength += localRun->fSum2NeutrTrLength;
0185
0186 G4Run::Merge(run);
0187 }
0188
0189
0190
0191 void Run::EndOfRun(G4double edep, G4double rms, G4double& limit)
0192 {
0193 G4int NbOfEvents = GetNumberOfEvent();
0194
0195 G4double kinEnergy = fKin->GetParticleGun()->GetParticleEnergy();
0196 assert(NbOfEvents * kinEnergy > 0);
0197
0198 fChargedStep /= G4double(NbOfEvents);
0199 fNeutralStep /= G4double(NbOfEvents);
0200
0201 G4double mass = fKin->GetParticleGun()->GetParticleDefinition()->GetPDGMass();
0202 G4double norme = 100. / (NbOfEvents * (kinEnergy + mass));
0203
0204
0205
0206 G4double dLradl = fDet->GetdLradl();
0207
0208 MyVector MeanELongit(f_nLbin), rmsELongit(f_nLbin);
0209 MyVector MeanELongitCumul(f_nLbin), rmsELongitCumul(f_nLbin);
0210
0211 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0212
0213 G4int i;
0214 for (i = 0; i < f_nLbin; ++i) {
0215 MeanELongit[i] = norme * fSumELongit[i];
0216 rmsELongit[i] =
0217 norme * std::sqrt(std::abs(NbOfEvents * fSumE2Longit[i] - fSumELongit[i] * fSumELongit[i]));
0218
0219 MeanELongitCumul[i] = norme * fSumELongitCumul[i];
0220 rmsELongitCumul[i] = norme
0221 * std::sqrt(std::abs(NbOfEvents * fSumE2LongitCumul[i]
0222 - fSumELongitCumul[i] * fSumELongitCumul[i]));
0223 G4double bin = (i + 0.5) * dLradl;
0224 analysisManager->FillH1(4, bin, MeanELongit[i] / dLradl);
0225 analysisManager->FillH1(5, bin, rmsELongit[i] / dLradl);
0226 bin = (i + 1) * dLradl;
0227 analysisManager->FillH1(6, bin, MeanELongitCumul[i]);
0228 analysisManager->FillH1(7, bin, rmsELongitCumul[i]);
0229 }
0230
0231
0232
0233 G4double dRradl = fDet->GetdRradl();
0234
0235 MyVector MeanERadial(f_nRbin), rmsERadial(f_nRbin);
0236 MyVector MeanERadialCumul(f_nRbin), rmsERadialCumul(f_nRbin);
0237
0238 for (i = 0; i < f_nRbin; ++i) {
0239 MeanERadial[i] = norme * fSumERadial[i];
0240 rmsERadial[i] =
0241 norme * std::sqrt(std::abs(NbOfEvents * fSumE2Radial[i] - fSumERadial[i] * fSumERadial[i]));
0242
0243 MeanERadialCumul[i] = norme * fSumERadialCumul[i];
0244 rmsERadialCumul[i] = norme
0245 * std::sqrt(std::abs(NbOfEvents * fSumE2RadialCumul[i]
0246 - fSumERadialCumul[i] * fSumERadialCumul[i]));
0247
0248 G4double bin = (i + 0.5) * dRradl;
0249 analysisManager->FillH1(8, bin, MeanERadial[i] / dRradl);
0250 analysisManager->FillH1(9, bin, rmsERadial[i] / dRradl);
0251 bin = (i + 1) * dRradl;
0252 analysisManager->FillH1(10, bin, MeanERadialCumul[i]);
0253 analysisManager->FillH1(11, bin, rmsERadialCumul[i]);
0254 }
0255
0256
0257
0258 const G4double EMoliere = 90.;
0259 G4double iMoliere = 0.;
0260 if ((MeanERadialCumul[0] <= EMoliere) && (MeanERadialCumul[f_nRbin - 1] >= EMoliere)) {
0261 G4int imin = 0;
0262 while ((imin < f_nRbin - 1) && (MeanERadialCumul[imin] < EMoliere)) {
0263 ++imin;
0264 }
0265
0266 G4double del = MeanERadialCumul[imin + 1] - MeanERadialCumul[imin];
0267 G4double ratio = (del > 0.0) ? (EMoliere - MeanERadialCumul[imin]) / del : 0.0;
0268 iMoliere = 1. + imin + ratio;
0269 }
0270
0271
0272
0273 norme = 1. / (NbOfEvents * (fDet->GetMaterial()->GetRadlen()));
0274 G4double MeanChargTrLength = norme * fSumChargTrLength;
0275 G4double rmsChargTrLength =
0276 norme
0277 * std::sqrt(std::abs(NbOfEvents * fSum2ChargTrLength - fSumChargTrLength * fSumChargTrLength));
0278
0279 G4double MeanNeutrTrLength = norme * fSumNeutrTrLength;
0280 G4double rmsNeutrTrLength =
0281 norme
0282 * std::sqrt(std::abs(NbOfEvents * fSum2NeutrTrLength - fSumNeutrTrLength * fSumNeutrTrLength));
0283
0284
0285 std::ios::fmtflags mode = G4cout.flags();
0286 G4cout.setf(std::ios::fixed, std::ios::floatfield);
0287 G4int prec = G4cout.precision(2);
0288
0289 if (fVerbose) {
0290 G4cout << " LOGITUDINAL PROFILE "
0291 << " CUMULATIVE LOGITUDINAL PROFILE" << G4endl << G4endl;
0292
0293 G4cout << " bin "
0294 << " Mean rms "
0295 << " bin "
0296 << " Mean rms \n"
0297 << G4endl;
0298
0299 for (i = 0; i < f_nLbin; ++i) {
0300 G4double inf = i * dLradl, sup = inf + dLradl;
0301
0302 G4cout << std::setw(8) << inf << "->" << std::setw(5) << sup << " radl: " << std::setw(7)
0303 << MeanELongit[i] << "% " << std::setw(9) << rmsELongit[i] << "% "
0304 << " 0->" << std::setw(5) << sup << " radl: " << std::setw(7)
0305 << MeanELongitCumul[i] << "% " << std::setw(7) << rmsELongitCumul[i] << "% "
0306 << G4endl;
0307 }
0308
0309 G4cout << G4endl << G4endl << G4endl;
0310
0311 G4cout << " RADIAL PROFILE "
0312 << " CUMULATIVE RADIAL PROFILE" << G4endl << G4endl;
0313
0314 G4cout << " bin "
0315 << " Mean rms "
0316 << " bin "
0317 << " Mean rms \n"
0318 << G4endl;
0319
0320 for (i = 0; i < f_nRbin; ++i) {
0321 G4double inf = i * dRradl, sup = inf + dRradl;
0322
0323 G4cout << std::setw(8) << inf << "->" << std::setw(5) << sup << " radl: " << std::setw(7)
0324 << MeanERadial[i] << "% " << std::setw(9) << rmsERadial[i] << "% "
0325 << " 0->" << std::setw(5) << sup << " radl: " << std::setw(7)
0326 << MeanERadialCumul[i] << "% " << std::setw(7) << rmsERadialCumul[i] << "% "
0327 << G4endl;
0328 }
0329 }
0330
0331 G4cout << "\n ===== SUMMARY ===== \n" << G4endl;
0332
0333 G4cout << " Total number of events: " << NbOfEvents << "\n"
0334 << " Mean number of charged steps: " << fChargedStep << G4endl;
0335 G4cout << " Mean number of neutral steps: " << fNeutralStep << "\n" << G4endl;
0336
0337 G4cout << " energy deposit : " << std::setw(7) << MeanELongitCumul[f_nLbin - 1] << " % E0 +- "
0338 << std::setw(7) << rmsELongitCumul[f_nLbin - 1] << " % E0" << G4endl;
0339 G4cout << " charged traklen: " << std::setw(7) << MeanChargTrLength << " radl +- " << std::setw(7)
0340 << rmsChargTrLength << " radl" << G4endl;
0341 G4cout << " neutral traklen: " << std::setw(7) << MeanNeutrTrLength << " radl +- " << std::setw(7)
0342 << rmsNeutrTrLength << " radl" << G4endl;
0343
0344 if (iMoliere > 0.) {
0345 G4double RMoliere1 = iMoliere * fDet->GetdRradl();
0346 G4double RMoliere2 = iMoliere * fDet->GetdRlength();
0347 G4cout << "\n " << EMoliere << " % confinement: radius = " << RMoliere1 << " radl ("
0348 << G4BestUnit(RMoliere2, "Length") << ")"
0349 << "\n"
0350 << G4endl;
0351 }
0352
0353 G4cout.setf(mode, std::ios::floatfield);
0354 G4cout.precision(prec);
0355
0356
0357
0358 G4int nLbin = fDet->GetnLtot();
0359 if (limit < DBL_MAX) {
0360 EmAcceptance acc;
0361 acc.BeginOfAcceptance("Total Energy in Absorber", NbOfEvents);
0362 G4double e = MeanELongitCumul[nLbin - 1] / 100.;
0363 G4double r = rmsELongitCumul[nLbin - 1] / 100.;
0364 acc.EmAcceptanceGauss("Edep", NbOfEvents, e, edep, rms, limit);
0365 acc.EmAcceptanceGauss("Erms", NbOfEvents, r, rms, rms, 2.0 * limit);
0366 acc.EndOfAcceptance();
0367 }
0368 limit = DBL_MAX;
0369 }
0370
0371
