File indexing completed on 2026-05-26 07:40:20
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 "EmAcceptance.hh"
0033 #include "HistoManager.hh"
0034 #include "PrimaryGeneratorAction.hh"
0035
0036 #include "G4Electron.hh"
0037 #include "G4Gamma.hh"
0038 #include "G4ParticleDefinition.hh"
0039 #include "G4ParticleTable.hh"
0040 #include "G4Positron.hh"
0041 #include "G4SystemOfUnits.hh"
0042 #include "G4Track.hh"
0043 #include "G4UnitsTable.hh"
0044
0045 #include <iomanip>
0046
0047
0048
0049 Run::Run(DetectorConstruction* det) : fDetector(det)
0050 {
0051
0052
0053 for (G4int k = 0; k < kMaxAbsor; k++) {
0054 fSumEAbs[k] = fSum2EAbs[k] = fSumLAbs[k] = fSum2LAbs[k] = 0.;
0055 fEnergyDeposit[k].clear();
0056 fEdeptrue[k] = fRmstrue[k] = 1.;
0057 fLimittrue[k] = DBL_MAX;
0058 }
0059
0060 fEdepTot = fEdepTot2 = 0.;
0061 fEleakTot = fEleakTot2 = 0.;
0062 fEtotal = fEtotal2 = 0.;
0063
0064
0065
0066 G4int nbPlanes = (fDetector->GetNbOfLayers()) * (fDetector->GetNbOfAbsor()) + 2;
0067 fEnergyFlow.resize(nbPlanes);
0068 fLateralEleak.resize(nbPlanes);
0069 for (G4int k = 0; k < nbPlanes; k++) {
0070 fEnergyFlow[k] = fLateralEleak[k] = 0.;
0071 }
0072 }
0073
0074
0075
0076 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
0077 {
0078 fParticle = particle;
0079 fEkin = energy;
0080 }
0081
0082
0083
0084 void Run::FillPerEvent(G4int kAbs, G4double EAbs, G4double LAbs)
0085 {
0086
0087
0088 if (fApplyLimit) fEnergyDeposit[kAbs].push_back(EAbs);
0089 fSumEAbs[kAbs] += EAbs;
0090 fSum2EAbs[kAbs] += EAbs * EAbs;
0091 fSumLAbs[kAbs] += LAbs;
0092 fSum2LAbs[kAbs] += LAbs * LAbs;
0093 }
0094
0095
0096
0097 void Run::SumEnergies(G4double edeptot, G4double eleak)
0098 {
0099 fEdepTot += edeptot;
0100 fEdepTot2 += edeptot * edeptot;
0101 fEleakTot += eleak;
0102 fEleakTot2 += eleak * eleak;
0103
0104 G4double etotal = edeptot + eleak;
0105 fEtotal += etotal;
0106 fEtotal2 += etotal * etotal;
0107 }
0108
0109
0110
0111 void Run::SumEnergyFlow(G4int plane, G4double Eflow)
0112 {
0113 fEnergyFlow[plane] += Eflow;
0114 }
0115
0116
0117
0118 void Run::SumLateralEleak(G4int cell, G4double Eflow)
0119 {
0120 fLateralEleak[cell] += Eflow;
0121 }
0122
0123
0124
0125 void Run::AddChargedStep()
0126 {
0127 fChargedStep += 1.0;
0128 }
0129
0130
0131
0132 void Run::AddNeutralStep()
0133 {
0134 fNeutralStep += 1.0;
0135 }
0136
0137
0138
0139 void Run::AddSecondaryTrack(const G4Track* track)
0140 {
0141 const G4ParticleDefinition* d = track->GetDefinition();
0142 if (d == G4Gamma::Gamma()) {
0143 ++fN_gamma;
0144 }
0145 else if (d == G4Electron::Electron()) {
0146 ++fN_elec;
0147 }
0148 else if (d == G4Positron::Positron()) {
0149 ++fN_pos;
0150 }
0151 }
0152
0153
0154
0155 void Run::Merge(const G4Run* run)
0156 {
0157 const Run* localRun = static_cast<const Run*>(run);
0158
0159
0160 fParticle = localRun->fParticle;
0161 fEkin = localRun->fEkin;
0162
0163
0164
0165 for (G4int k = 0; k < kMaxAbsor; k++) {
0166 fSumEAbs[k] += localRun->fSumEAbs[k];
0167 fSum2EAbs[k] += localRun->fSum2EAbs[k];
0168 fSumLAbs[k] += localRun->fSumLAbs[k];
0169 fSum2LAbs[k] += localRun->fSum2LAbs[k];
0170 }
0171
0172 fEdepTot += localRun->fEdepTot;
0173 fEdepTot2 += localRun->fEdepTot2;
0174
0175 fEleakTot += localRun->fEleakTot;
0176 fEleakTot2 += localRun->fEleakTot2;
0177
0178 fEtotal += localRun->fEtotal;
0179 fEtotal2 += localRun->fEtotal2;
0180
0181 G4int nbPlanes = (fDetector->GetNbOfLayers()) * (fDetector->GetNbOfAbsor()) + 2;
0182 for (G4int k = 0; k < nbPlanes; k++) {
0183 fEnergyFlow[k] += localRun->fEnergyFlow[k];
0184 fLateralEleak[k] += localRun->fLateralEleak[k];
0185 }
0186
0187 for (G4int k=0; k<kMaxAbsor; k++) {
0188 fEnergyDeposit[k].insert(fEnergyDeposit[k].end(),
0189 localRun->fEnergyDeposit[k].begin(), localRun->fEnergyDeposit[k].end());
0190 }
0191
0192 fChargedStep += localRun->fChargedStep;
0193 fNeutralStep += localRun->fNeutralStep;
0194
0195 fN_gamma += localRun->fN_gamma;
0196 fN_elec += localRun->fN_elec;
0197 fN_pos += localRun->fN_pos;
0198
0199 fApplyLimit = localRun->fApplyLimit;
0200
0201 for (G4int k = 0; k < kMaxAbsor; k++) {
0202 fEdeptrue[k] = localRun->fEdeptrue[k];
0203 fRmstrue[k] = localRun->fRmstrue[k];
0204 fLimittrue[k] = localRun->fLimittrue[k];
0205 }
0206
0207 G4Run::Merge(run);
0208 }
0209
0210
0211
0212 void Run::EndOfRun()
0213 {
0214 G4int nEvt = numberOfEvent;
0215 G4double norm = G4double(nEvt);
0216 if (norm > 0) norm = 1. / norm;
0217 G4double qnorm = std::sqrt(norm);
0218
0219 fChargedStep *= norm;
0220 fNeutralStep *= norm;
0221
0222
0223
0224 G4double beamEnergy = fEkin;
0225 G4double sqbeam = std::sqrt(beamEnergy / GeV);
0226
0227 G4double MeanEAbs, MeanEAbs2, rmsEAbs, resolution, rmsres;
0228 G4double MeanLAbs, MeanLAbs2, rmsLAbs;
0229
0230 std::ios::fmtflags mode = G4cout.flags();
0231 G4int prec = G4cout.precision(2);
0232 G4cout << "\n------------------------------------------------------------\n";
0233 G4cout << std::setw(14) << "material" << std::setw(17) << "Edep RMS" << std::setw(33)
0234 << "sqrt(E0(GeV))*rmsE/Emean" << std::setw(23) << "total tracklen \n \n";
0235
0236 for (G4int k = 1; k <= fDetector->GetNbOfAbsor(); k++) {
0237 MeanEAbs = fSumEAbs[k] * norm;
0238 MeanEAbs2 = fSum2EAbs[k] * norm;
0239 rmsEAbs = std::sqrt(std::abs(MeanEAbs2 - MeanEAbs * MeanEAbs));
0240
0241
0242 if (fApplyLimit) {
0243 G4int nn = 0;
0244 G4double sume = 0.0;
0245 G4double sume2 = 0.0;
0246
0247 G4double lim = rmsEAbs * 2.5;
0248 for (G4int i = 0; i < nEvt; i++) {
0249 G4double e = (fEnergyDeposit[k])[i];
0250 if (std::abs(e - MeanEAbs) < lim) {
0251 sume += e;
0252 sume2 += e * e;
0253 nn++;
0254 }
0255 }
0256 G4double norm1 = G4double(nn);
0257 if (norm1 > 0.0) norm1 = 1.0 / norm1;
0258 MeanEAbs = sume * norm1;
0259 MeanEAbs2 = sume2 * norm1;
0260 rmsEAbs = std::sqrt(std::abs(MeanEAbs2 - MeanEAbs * MeanEAbs));
0261 }
0262
0263 resolution = (MeanEAbs > 0.) ? 100. * sqbeam * rmsEAbs / MeanEAbs : 0.0;
0264 rmsres = resolution * qnorm;
0265
0266
0267 fSumEAbs[k] = MeanEAbs;
0268 fSum2EAbs[k] = rmsEAbs;
0269
0270 MeanLAbs = fSumLAbs[k] * norm;
0271 MeanLAbs2 = fSum2LAbs[k] * norm;
0272 rmsLAbs = std::sqrt(std::abs(MeanLAbs2 - MeanLAbs * MeanLAbs));
0273
0274
0275
0276 G4cout << std::setw(14) << fDetector->GetAbsorMaterial(k)->GetName() << ": "
0277 << std::setprecision(5) << std::setw(6) << G4BestUnit(MeanEAbs, "Energy") << " : "
0278 << std::setprecision(4) << std::setw(5) << G4BestUnit(rmsEAbs, "Energy") << std::setw(10)
0279 << resolution << " +- " << std::setw(5) << rmsres << " %" << std::setprecision(3)
0280 << std::setw(10) << G4BestUnit(MeanLAbs, "Length") << " +- " << std::setw(4)
0281 << G4BestUnit(rmsLAbs, "Length") << G4endl;
0282 }
0283
0284
0285
0286 fEdepTot *= norm;
0287 fEdepTot2 *= norm;
0288 G4double rmsEdep = std::sqrt(std::abs(fEdepTot2 - fEdepTot * fEdepTot));
0289
0290 G4cout << "\n Total energy deposited = " << std::setprecision(4) << G4BestUnit(fEdepTot, "Energy")
0291 << " +- " << G4BestUnit(rmsEdep, "Energy") << G4endl;
0292
0293
0294
0295 fEleakTot *= norm;
0296 fEleakTot2 *= norm;
0297 G4double rmsEleak = std::sqrt(std::abs(fEleakTot2 - fEleakTot * fEleakTot));
0298
0299 G4cout << " Energy leakage = " << G4BestUnit(fEleakTot, "Energy") << " +- "
0300 << G4BestUnit(rmsEleak, "Energy") << G4endl;
0301
0302
0303
0304 fEtotal *= norm;
0305 fEtotal2 *= norm;
0306 G4double rmsEtotal = std::sqrt(std::abs(fEtotal2 - fEtotal * fEtotal));
0307
0308 G4cout << " Total energy : Edep + Eleak = " << G4BestUnit(fEtotal, "Energy") << " +- "
0309 << G4BestUnit(rmsEtotal, "Energy") << G4endl;
0310
0311 G4cout << "------------------------------------------------------------\n";
0312
0313 G4cout << " Beam particle " << fParticle->GetParticleName()
0314 << " E = " << G4BestUnit(beamEnergy, "Energy") << G4endl;
0315 G4cout << " Mean number of gamma " << (G4double)fN_gamma * norm << G4endl;
0316 G4cout << " Mean number of e- " << (G4double)fN_elec * norm << G4endl;
0317 G4cout << " Mean number of e+ " << (G4double)fN_pos * norm << G4endl;
0318 G4cout << std::setprecision(6) << " Mean number of charged steps " << fChargedStep << G4endl;
0319 G4cout << " Mean number of neutral steps " << fNeutralStep << G4endl;
0320 G4cout << "------------------------------------------------------------\n";
0321
0322
0323
0324 G4AnalysisManager* analysis = G4AnalysisManager::Instance();
0325 G4int Idmax = (fDetector->GetNbOfLayers()) * (fDetector->GetNbOfAbsor());
0326 for (G4int Id = 1; Id <= Idmax + 1; Id++) {
0327 analysis->FillH1(2 * kMaxAbsor + 1, (G4double)Id, fEnergyFlow[Id]);
0328 analysis->FillH1(2 * kMaxAbsor + 2, (G4double)Id, fLateralEleak[Id]);
0329 }
0330
0331
0332
0333 G4double EdepTot[kMaxAbsor];
0334 for (G4int k = 0; k < kMaxAbsor; k++)
0335 EdepTot[k] = 0.;
0336
0337 G4int nbOfAbsor = fDetector->GetNbOfAbsor();
0338 for (G4int Id = 1; Id <= Idmax; Id++) {
0339 G4int iAbsor = Id % nbOfAbsor;
0340 if (iAbsor == 0) iAbsor = nbOfAbsor;
0341 EdepTot[iAbsor] += (fEnergyFlow[Id] - fEnergyFlow[Id + 1] - fLateralEleak[Id]);
0342 }
0343
0344 G4cout << std::setprecision(3) << "\n Energy deposition from Energy flow balance : \n"
0345 << std::setw(10) << " material \t Total Edep \n \n";
0346 G4cout.precision(6);
0347
0348 for (G4int k = 1; k <= nbOfAbsor; k++) {
0349 EdepTot[k] *= norm;
0350 G4cout << std::setw(10) << fDetector->GetAbsorMaterial(k)->GetName() << ":"
0351 << "\t " << G4BestUnit(EdepTot[k], "Energy") << "\n";
0352 }
0353
0354 G4cout << "\n------------------------------------------------------------\n" << G4endl;
0355
0356
0357 EmAcceptance acc;
0358 G4bool isStarted = false;
0359 for (G4int j = 1; j <= fDetector->GetNbOfAbsor(); j++) {
0360 if (fLimittrue[j] < DBL_MAX) {
0361 if (!isStarted) {
0362 acc.BeginOfAcceptance("Sampling Calorimeter", nEvt);
0363 isStarted = true;
0364 }
0365 MeanEAbs = fSumEAbs[j];
0366 rmsEAbs = fSum2EAbs[j];
0367 G4String mat = fDetector->GetAbsorMaterial(j)->GetName();
0368 acc.EmAcceptanceGauss("Edep" + mat, nEvt, MeanEAbs, fEdeptrue[j], fRmstrue[j], fLimittrue[j]);
0369 acc.EmAcceptanceGauss("Erms" + mat, nEvt, rmsEAbs, fRmstrue[j], fRmstrue[j],
0370 2.0 * fLimittrue[j]);
0371 }
0372 }
0373 if (isStarted) acc.EndOfAcceptance();
0374
0375
0376
0377 for (G4int ih = kMaxAbsor + 1; ih < kMaxHisto - 2; ih++) {
0378 analysis->ScaleH1(ih, norm / MeV);
0379 }
0380
0381 G4cout.setf(mode, std::ios::floatfield);
0382 G4cout.precision(prec);
0383 }
0384
0385
0386
0387 void Run::SetEdepAndRMS(G4int i, G4double edep, G4double rms, G4double lim)
0388 {
0389 if (i >= 0 && i < kMaxAbsor) {
0390 fEdeptrue[i] = edep;
0391 fRmstrue[i] = rms;
0392 fLimittrue[i] = lim;
0393 }
0394 }
0395
0396
0397
0398 void Run::SetApplyLimit(G4bool val)
0399 {
0400 fApplyLimit = val;
0401 }
0402
0403
