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