File indexing completed on 2025-02-23 09:20:59
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 fChargedStep += localRun->fChargedStep;
0192 fNeutralStep += localRun->fNeutralStep;
0193
0194 fN_gamma += localRun->fN_gamma;
0195 fN_elec += localRun->fN_elec;
0196 fN_pos += localRun->fN_pos;
0197
0198 fApplyLimit = localRun->fApplyLimit;
0199
0200 for (G4int k = 0; k < kMaxAbsor; k++) {
0201 fEdeptrue[k] = localRun->fEdeptrue[k];
0202 fRmstrue[k] = localRun->fRmstrue[k];
0203 fLimittrue[k] = localRun->fLimittrue[k];
0204 }
0205
0206 G4Run::Merge(run);
0207 }
0208
0209
0210
0211 void Run::EndOfRun()
0212 {
0213 G4int nEvt = numberOfEvent;
0214 G4double norm = G4double(nEvt);
0215 if (norm > 0) norm = 1. / norm;
0216 G4double qnorm = std::sqrt(norm);
0217
0218 fChargedStep *= norm;
0219 fNeutralStep *= norm;
0220
0221
0222
0223 G4double beamEnergy = fEkin;
0224 G4double sqbeam = std::sqrt(beamEnergy / GeV);
0225
0226 G4double MeanEAbs, MeanEAbs2, rmsEAbs, resolution, rmsres;
0227 G4double MeanLAbs, MeanLAbs2, rmsLAbs;
0228
0229 std::ios::fmtflags mode = G4cout.flags();
0230 G4int prec = G4cout.precision(2);
0231 G4cout << "\n------------------------------------------------------------\n";
0232 G4cout << std::setw(14) << "material" << std::setw(17) << "Edep RMS" << std::setw(33)
0233 << "sqrt(E0(GeV))*rmsE/Emean" << std::setw(23) << "total tracklen \n \n";
0234
0235 for (G4int k = 1; k <= fDetector->GetNbOfAbsor(); k++) {
0236 MeanEAbs = fSumEAbs[k] * norm;
0237 MeanEAbs2 = fSum2EAbs[k] * norm;
0238 rmsEAbs = std::sqrt(std::abs(MeanEAbs2 - MeanEAbs * MeanEAbs));
0239
0240
0241 if (fApplyLimit) {
0242 G4int nn = 0;
0243 G4double sume = 0.0;
0244 G4double sume2 = 0.0;
0245
0246 G4double lim = rmsEAbs * 2.5;
0247 for (G4int i = 0; i < nEvt; i++) {
0248 G4double e = (fEnergyDeposit[k])[i];
0249 if (std::abs(e - MeanEAbs) < lim) {
0250 sume += e;
0251 sume2 += e * e;
0252 nn++;
0253 }
0254 }
0255 G4double norm1 = G4double(nn);
0256 if (norm1 > 0.0) norm1 = 1.0 / norm1;
0257 MeanEAbs = sume * norm1;
0258 MeanEAbs2 = sume2 * norm1;
0259 rmsEAbs = std::sqrt(std::abs(MeanEAbs2 - MeanEAbs * MeanEAbs));
0260 }
0261
0262 resolution = (MeanEAbs > 0.) ? 100. * sqbeam * rmsEAbs / MeanEAbs : 0.0;
0263 rmsres = resolution * qnorm;
0264
0265
0266 fSumEAbs[k] = MeanEAbs;
0267 fSum2EAbs[k] = rmsEAbs;
0268
0269 MeanLAbs = fSumLAbs[k] * norm;
0270 MeanLAbs2 = fSum2LAbs[k] * norm;
0271 rmsLAbs = std::sqrt(std::abs(MeanLAbs2 - MeanLAbs * MeanLAbs));
0272
0273
0274
0275 G4cout << std::setw(14) << fDetector->GetAbsorMaterial(k)->GetName() << ": "
0276 << std::setprecision(5) << std::setw(6) << G4BestUnit(MeanEAbs, "Energy") << " : "
0277 << std::setprecision(4) << std::setw(5) << G4BestUnit(rmsEAbs, "Energy") << std::setw(10)
0278 << resolution << " +- " << std::setw(5) << rmsres << " %" << std::setprecision(3)
0279 << std::setw(10) << G4BestUnit(MeanLAbs, "Length") << " +- " << std::setw(4)
0280 << G4BestUnit(rmsLAbs, "Length") << G4endl;
0281 }
0282
0283
0284
0285 fEdepTot *= norm;
0286 fEdepTot2 *= norm;
0287 G4double rmsEdep = std::sqrt(std::abs(fEdepTot2 - fEdepTot * fEdepTot));
0288
0289 G4cout << "\n Total energy deposited = " << std::setprecision(4) << G4BestUnit(fEdepTot, "Energy")
0290 << " +- " << G4BestUnit(rmsEdep, "Energy") << G4endl;
0291
0292
0293
0294 fEleakTot *= norm;
0295 fEleakTot2 *= norm;
0296 G4double rmsEleak = std::sqrt(std::abs(fEleakTot2 - fEleakTot * fEleakTot));
0297
0298 G4cout << " Energy leakage = " << G4BestUnit(fEleakTot, "Energy") << " +- "
0299 << G4BestUnit(rmsEleak, "Energy") << G4endl;
0300
0301
0302
0303 fEtotal *= norm;
0304 fEtotal2 *= norm;
0305 G4double rmsEtotal = std::sqrt(std::abs(fEtotal2 - fEtotal * fEtotal));
0306
0307 G4cout << " Total energy : Edep + Eleak = " << G4BestUnit(fEtotal, "Energy") << " +- "
0308 << G4BestUnit(rmsEtotal, "Energy") << G4endl;
0309
0310 G4cout << "------------------------------------------------------------\n";
0311
0312 G4cout << " Beam particle " << fParticle->GetParticleName()
0313 << " E = " << G4BestUnit(beamEnergy, "Energy") << G4endl;
0314 G4cout << " Mean number of gamma " << (G4double)fN_gamma * norm << G4endl;
0315 G4cout << " Mean number of e- " << (G4double)fN_elec * norm << G4endl;
0316 G4cout << " Mean number of e+ " << (G4double)fN_pos * norm << G4endl;
0317 G4cout << std::setprecision(6) << " Mean number of charged steps " << fChargedStep << G4endl;
0318 G4cout << " Mean number of neutral steps " << fNeutralStep << G4endl;
0319 G4cout << "------------------------------------------------------------\n";
0320
0321
0322
0323 G4AnalysisManager* analysis = G4AnalysisManager::Instance();
0324 G4int Idmax = (fDetector->GetNbOfLayers()) * (fDetector->GetNbOfAbsor());
0325 for (G4int Id = 1; Id <= Idmax + 1; Id++) {
0326 analysis->FillH1(2 * kMaxAbsor + 1, (G4double)Id, fEnergyFlow[Id]);
0327 analysis->FillH1(2 * kMaxAbsor + 2, (G4double)Id, fLateralEleak[Id]);
0328 }
0329
0330
0331
0332 G4double EdepTot[kMaxAbsor];
0333 for (G4int k = 0; k < kMaxAbsor; k++)
0334 EdepTot[k] = 0.;
0335
0336 G4int nbOfAbsor = fDetector->GetNbOfAbsor();
0337 for (G4int Id = 1; Id <= Idmax; Id++) {
0338 G4int iAbsor = Id % nbOfAbsor;
0339 if (iAbsor == 0) iAbsor = nbOfAbsor;
0340 EdepTot[iAbsor] += (fEnergyFlow[Id] - fEnergyFlow[Id + 1] - fLateralEleak[Id]);
0341 }
0342
0343 G4cout << std::setprecision(3) << "\n Energy deposition from Energy flow balance : \n"
0344 << std::setw(10) << " material \t Total Edep \n \n";
0345 G4cout.precision(6);
0346
0347 for (G4int k = 1; k <= nbOfAbsor; k++) {
0348 EdepTot[k] *= norm;
0349 G4cout << std::setw(10) << fDetector->GetAbsorMaterial(k)->GetName() << ":"
0350 << "\t " << G4BestUnit(EdepTot[k], "Energy") << "\n";
0351 }
0352
0353 G4cout << "\n------------------------------------------------------------\n" << G4endl;
0354
0355
0356 EmAcceptance acc;
0357 G4bool isStarted = false;
0358 for (G4int j = 1; j <= fDetector->GetNbOfAbsor(); j++) {
0359 if (fLimittrue[j] < DBL_MAX) {
0360 if (!isStarted) {
0361 acc.BeginOfAcceptance("Sampling Calorimeter", nEvt);
0362 isStarted = true;
0363 }
0364 MeanEAbs = fSumEAbs[j];
0365 rmsEAbs = fSum2EAbs[j];
0366 G4String mat = fDetector->GetAbsorMaterial(j)->GetName();
0367 acc.EmAcceptanceGauss("Edep" + mat, nEvt, MeanEAbs, fEdeptrue[j], fRmstrue[j], fLimittrue[j]);
0368 acc.EmAcceptanceGauss("Erms" + mat, nEvt, rmsEAbs, fRmstrue[j], fRmstrue[j],
0369 2.0 * fLimittrue[j]);
0370 }
0371 }
0372 if (isStarted) acc.EndOfAcceptance();
0373
0374
0375
0376 for (G4int ih = kMaxAbsor + 1; ih < kMaxHisto - 2; ih++) {
0377 analysis->ScaleH1(ih, norm / MeV);
0378 }
0379
0380 G4cout.setf(mode, std::ios::floatfield);
0381 G4cout.precision(prec);
0382 }
0383
0384
0385
0386 void Run::SetEdepAndRMS(G4int i, G4double edep, G4double rms, G4double lim)
0387 {
0388 if (i >= 0 && i < kMaxAbsor) {
0389 fEdeptrue[i] = edep;
0390 fRmstrue[i] = rms;
0391 fLimittrue[i] = lim;
0392 }
0393 }
0394
0395
0396
0397 void Run::SetApplyLimit(G4bool val)
0398 {
0399 fApplyLimit = val;
0400 }
0401
0402
