Warning, file /geant4/examples/extended/electromagnetic/TestEm3/src/Run.cc was not indexed
or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).
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
