File indexing completed on 2025-04-04 08:05: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
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046 #include "HistoManager.hh"
0047
0048 #include "Histo.hh"
0049
0050 #include "G4Alpha.hh"
0051 #include "G4AntiProton.hh"
0052 #include "G4Deuteron.hh"
0053 #include "G4Electron.hh"
0054 #include "G4Gamma.hh"
0055 #include "G4He3.hh"
0056 #include "G4KaonMinus.hh"
0057 #include "G4KaonPlus.hh"
0058 #include "G4KaonZeroLong.hh"
0059 #include "G4KaonZeroShort.hh"
0060 #include "G4MuonMinus.hh"
0061 #include "G4MuonPlus.hh"
0062 #include "G4Neutron.hh"
0063 #include "G4PhysicalConstants.hh"
0064 #include "G4PionMinus.hh"
0065 #include "G4PionPlus.hh"
0066 #include "G4PionZero.hh"
0067 #include "G4Positron.hh"
0068 #include "G4Proton.hh"
0069 #include "G4Step.hh"
0070 #include "G4SystemOfUnits.hh"
0071 #include "G4Track.hh"
0072 #include "G4Triton.hh"
0073 #include "G4UnitsTable.hh"
0074 #include "globals.hh"
0075
0076
0077
0078 HistoManager* HistoManager::fManager = nullptr;
0079
0080
0081
0082 HistoManager* HistoManager::GetPointer()
0083 {
0084 if (nullptr == fManager) {
0085 static HistoManager manager;
0086 fManager = &manager;
0087 }
0088 return fManager;
0089 }
0090
0091
0092
0093 HistoManager::HistoManager()
0094 : fPrimaryDef(nullptr),
0095 fRadius(10 * CLHEP::cm),
0096 fLength(300. * CLHEP::mm),
0097 fEdepMax(1.0 * CLHEP::GeV)
0098 {
0099 fHisto = new Histo();
0100 fHisto->SetVerbose(fVerbose);
0101 fNeutron = G4Neutron::Neutron();
0102 }
0103
0104
0105
0106 HistoManager::~HistoManager()
0107 {
0108 delete fHisto;
0109 }
0110
0111
0112
0113 void HistoManager::BookHisto()
0114 {
0115 fHistoBooked = true;
0116 fHisto->Add1D("1", "Energy deposition (MeV/mm/event) in the target", fNSlices, 0.0, fLength / mm,
0117 MeV / mm);
0118 fHisto->Add1D("2", "Log10 Energy (MeV) of gammas", fNBinsE, -5., 5., 1.0);
0119 fHisto->Add1D("3", "Log10 Energy (MeV) of electrons", fNBinsE, -5., 5., 1.0);
0120 fHisto->Add1D("4", "Log10 Energy (MeV) of positrons", fNBinsE, -5., 5., 1.0);
0121 fHisto->Add1D("5", "Log10 Energy (MeV) of protons", fNBinsE, -5., 5., 1.0);
0122 fHisto->Add1D("6", "Log10 Energy (MeV) of neutrons", fNBinsE, -5., 5., 1.0);
0123 fHisto->Add1D("7", "Log10 Energy (MeV) of charged pions", fNBinsE, -4., 6., 1.0);
0124 fHisto->Add1D("8", "Log10 Energy (MeV) of pi0", fNBinsE, -4., 6., 1.0);
0125 fHisto->Add1D("9", "Log10 Energy (MeV) of charged kaons", fNBinsE, -4., 6., 1.0);
0126 fHisto->Add1D("10", "Log10 Energy (MeV) of neutral kaons", fNBinsE, -4., 6., 1.0);
0127 fHisto->Add1D("11", "Log10 Energy (MeV) of deuterons and tritons", fNBinsE, -5., 5., 1.0);
0128 fHisto->Add1D("12", "Log10 Energy (MeV) of He3 and alpha", fNBinsE, -5., 5., 1.0);
0129 fHisto->Add1D("13", "Log10 Energy (MeV) of Generic Ions", fNBinsE, -5., 5., 1.0);
0130 fHisto->Add1D("14", "Log10 Energy (MeV) of muons", fNBinsE, -4., 6., 1.0);
0131 fHisto->Add1D("15", "log10 Energy (MeV) of side-leaked neutrons", fNBinsE, -5., 5., 1.0);
0132 fHisto->Add1D("16", "log10 Energy (MeV) of forward-leaked neutrons", fNBinsE, -5., 5., 1.0);
0133 fHisto->Add1D("17", "log10 Energy (MeV) of backward-leaked neutrons", fNBinsE, -5., 5., 1.0);
0134 fHisto->Add1D("18", "log10 Energy (MeV) of leaking protons", fNBinsE, -4., 6., 1.0);
0135 fHisto->Add1D("19", "log10 Energy (MeV) of leaking charged pions", fNBinsE, -4., 6., 1.0);
0136 fHisto->Add1D("20", "Log10 Energy (MeV) of pi+", fNBinsE, -4., 6., 1.0);
0137 fHisto->Add1D("21", "Log10 Energy (MeV) of pi-", fNBinsE, -4., 6., 1.0);
0138 fHisto->Add1D("22", "Energy deposition in the target normalized to beam energy", 110, 0.0, 1.1,
0139 1.0);
0140 fHisto->Add1D("23", "EM energy deposition in the target normalized to beam energy", 110, 0.0, 1.1,
0141 1.0);
0142 fHisto->Add1D("24", "Pion energy deposition in the target normalized to beam energy", 110, 0.0,
0143 1.1, 1.0);
0144 fHisto->Add1D("25", "Proton energy deposition in the target normalized to beam energy", 110, 0.0,
0145 1.1, 1.0);
0146 fHisto->Add1D("26", "Energy (MeV) of pi+", fNBinsE, 0., 500., 1.0);
0147 fHisto->Add1D("27", "Energy (MeV) of pi-", fNBinsE, 0., 500., 1.0);
0148 fHisto->Add1D("28", "Energy (MeV) of pi0", fNBinsE, 0., 500., 1.0);
0149 }
0150
0151
0152
0153 void HistoManager::BeginOfRun()
0154 {
0155 fR2 = fRadius * fRadius;
0156 fAbsZ0 = 0.5 * fLength;
0157 fNevt = 0;
0158 fNelec = 0;
0159 fNposit = 0;
0160 fNgam = 0;
0161 fNstep = 0;
0162 fNprot_leak = 0;
0163 fNpiofNleak = 0;
0164 fNions = 0;
0165 fNdeut = 0;
0166 fNalpha = 0;
0167 fNkaons = 0;
0168 fNmuons = 0;
0169 fNcpions = 0;
0170 fNpi0 = 0;
0171 fNneutron = 0;
0172 fNproton = 0;
0173 fNaproton = 0;
0174 fNneu_forw = 0;
0175 fNneu_leak = 0;
0176 fNneu_back = 0;
0177
0178 fEdepSum = 0.0;
0179 fEdepSum2 = 0.0;
0180
0181 if (!fHistoBooked) {
0182 BookHisto();
0183 }
0184 fHisto->Book();
0185
0186 if (fVerbose > 0) {
0187 G4cout << "HistoManager: Histograms are booked and run has been started" << G4endl;
0188 }
0189 }
0190
0191
0192
0193 void HistoManager::EndOfRun()
0194 {
0195 G4cout << "HistoManager: End of run actions are started" << G4endl;
0196
0197
0198 G4cout << "========================================================" << G4endl;
0199
0200 G4double x = (G4double)fNevt;
0201 if (fNevt > 0) {
0202 x = 1.0 / x;
0203 }
0204
0205 G4double xe = x * (G4double)fNelec;
0206 G4double xg = x * (G4double)fNgam;
0207 G4double xp = x * (G4double)fNposit;
0208 G4double xs = x * (G4double)fNstep;
0209 G4double xn = x * (G4double)fNneutron;
0210 G4double xpn = x * (G4double)fNproton;
0211 G4double xap = x * (G4double)fNaproton;
0212 G4double xnf = x * (G4double)fNneu_forw;
0213 G4double xnb = x * (G4double)fNneu_leak;
0214 G4double xnbw = x * (G4double)fNneu_back;
0215 G4double xpl = x * (G4double)fNprot_leak;
0216 G4double xal = x * (G4double)fNpiofNleak;
0217 G4double xpc = x * (G4double)fNcpions;
0218 G4double xp0 = x * (G4double)fNpi0;
0219 G4double xpk = x * (G4double)fNkaons;
0220 G4double xpm = x * (G4double)fNmuons;
0221 G4double xid = x * (G4double)fNdeut;
0222 G4double xia = x * (G4double)fNalpha;
0223 G4double xio = x * (G4double)fNions;
0224
0225 fEdepSum *= x;
0226 fEdepSum2 *= x;
0227 fEdepSum2 -= fEdepSum * fEdepSum;
0228 if (fEdepSum2 > 0.0) {
0229 fEdepSum2 = std::sqrt(fEdepSum2);
0230 }
0231 else {
0232 fEdepSum2 = 0.0;
0233 }
0234
0235 G4cout << "Beam particle " << fPrimaryDef->GetParticleName() << G4endl;
0236 G4cout << "Beam Energy(MeV) " << fPrimaryKineticEnergy / MeV << G4endl;
0237 G4cout << "Number of events " << fNevt << G4endl;
0238 G4cout << std::setprecision(4) << "Average energy deposit (MeV) " << fEdepSum / MeV
0239 << " RMS(MeV) " << fEdepSum2 / MeV << G4endl;
0240 G4cout << std::setprecision(4) << "Average number of steps " << xs << G4endl;
0241 G4cout << std::setprecision(4) << "Average number of gamma " << xg << G4endl;
0242 G4cout << std::setprecision(4) << "Average number of e- " << xe << G4endl;
0243 G4cout << std::setprecision(4) << "Average number of e+ " << xp << G4endl;
0244 G4cout << std::setprecision(4) << "Average number of neutrons " << xn << G4endl;
0245 G4cout << std::setprecision(4) << "Average number of protons " << xpn << G4endl;
0246 G4cout << std::setprecision(4) << "Average number of antiprotons " << xap << G4endl;
0247 G4cout << std::setprecision(4) << "Average number of pi+ & pi- " << xpc << G4endl;
0248 G4cout << std::setprecision(4) << "Average number of pi0 " << xp0 << G4endl;
0249 G4cout << std::setprecision(4) << "Average number of kaons " << xpk << G4endl;
0250 G4cout << std::setprecision(4) << "Average number of muons " << xpm << G4endl;
0251 G4cout << std::setprecision(4) << "Average number of deuterons+tritons " << xid << G4endl;
0252 G4cout << std::setprecision(4) << "Average number of He3+alpha " << xia << G4endl;
0253 G4cout << std::setprecision(4) << "Average number of ions " << xio << G4endl;
0254 G4cout << std::setprecision(4) << "Average number of forward neutrons " << xnf << G4endl;
0255 G4cout << std::setprecision(4) << "Average number of reflected neutrons " << xnb << G4endl;
0256 G4cout << std::setprecision(4) << "Average number of leaked neutrons " << xnbw << G4endl;
0257 G4cout << std::setprecision(4) << "Average number of proton leak " << xpl << G4endl;
0258 G4cout << std::setprecision(4) << "Average number of pion leak " << xal << G4endl;
0259 G4cout << "========================================================" << G4endl;
0260 G4cout << G4endl;
0261
0262
0263 for (G4int i = 0; i < fNHisto; i++) {
0264 fHisto->ScaleH1(i, x);
0265 }
0266
0267 fHisto->Save();
0268 }
0269
0270
0271
0272 void HistoManager::BeginOfEvent()
0273 {
0274 fEdepEvt = 0.0;
0275 fEdepEM = 0.0;
0276 fEdepPI = 0.0;
0277 fEdepP = 0.0;
0278 }
0279
0280
0281
0282 void HistoManager::EndOfEvent()
0283 {
0284 fEdepSum += fEdepEvt;
0285 fEdepSum2 += fEdepEvt * fEdepEvt;
0286 fHisto->Fill(21, fEdepEvt / fPrimaryKineticEnergy, 1.0);
0287 fHisto->Fill(22, fEdepEM / fPrimaryKineticEnergy, 1.0);
0288 fHisto->Fill(23, fEdepPI / fPrimaryKineticEnergy, 1.0);
0289 fHisto->Fill(24, fEdepP / fPrimaryKineticEnergy, 1.0);
0290 }
0291
0292
0293
0294 void HistoManager::ScoreNewTrack(const G4Track* track)
0295 {
0296 const G4ParticleDefinition* pd = track->GetDefinition();
0297 const G4String name = pd->GetParticleName();
0298 const G4double ee = track->GetKineticEnergy();
0299 G4double e = ee;
0300
0301
0302 if (0 == track->GetParentID()) {
0303 fNevt++;
0304 fPrimaryKineticEnergy = e;
0305 fPrimaryDef = pd;
0306 G4ThreeVector dir = track->GetMomentumDirection();
0307 if (1 < fVerbose)
0308 G4cout << "### Primary " << name << " kinE(MeV)= " << e / MeV
0309 << "; m(MeV)= " << pd->GetPDGMass() / MeV << "; pos(mm)= " << track->GetPosition() / mm
0310 << "; dir= " << track->GetMomentumDirection() << G4endl;
0311
0312
0313 }
0314 else {
0315 if (1 < fVerbose)
0316 G4cout << "=== Secondary " << name << " kinE(MeV)= " << e / MeV
0317 << "; m(MeV)= " << pd->GetPDGMass() / MeV << "; pos(mm)= " << track->GetPosition() / mm
0318 << "; dir= " << track->GetMomentumDirection() << G4endl;
0319 e = (e > 0.0) ? std::log10(e / MeV) : -100.;
0320 if (pd == G4Gamma::Gamma()) {
0321 fNgam++;
0322 fHisto->Fill(1, e, 1.0);
0323 }
0324 else if (pd == G4Electron::Electron()) {
0325 fNelec++;
0326 fHisto->Fill(2, e, 1.0);
0327 }
0328 else if (pd == G4Positron::Positron()) {
0329 fNposit++;
0330 fHisto->Fill(3, e, 1.0);
0331 }
0332 else if (pd == G4Proton::Proton()) {
0333 fNproton++;
0334 fHisto->Fill(4, e, 1.0);
0335 }
0336 else if (pd == fNeutron) {
0337 fNneutron++;
0338 fHisto->Fill(5, e, 1.0);
0339 }
0340 else if (pd == G4AntiProton::AntiProton()) {
0341 fNaproton++;
0342 }
0343 else if (pd == G4PionPlus::PionPlus()) {
0344 fNcpions++;
0345 fHisto->Fill(6, e, 1.0);
0346 fHisto->Fill(19, e, 1.0);
0347 fHisto->Fill(25, ee, 1.0);
0348 }
0349 else if (pd == G4PionMinus::PionMinus()) {
0350 fNcpions++;
0351 fHisto->Fill(6, e, 1.0);
0352 fHisto->Fill(20, e, 1.0);
0353 fHisto->Fill(26, ee, 1.0);
0354 }
0355 else if (pd == G4PionZero::PionZero()) {
0356 fNpi0++;
0357 fHisto->Fill(7, e, 1.0);
0358 fHisto->Fill(27, ee, 1.0);
0359 }
0360 else if (pd == G4KaonPlus::KaonPlus() || pd == G4KaonMinus::KaonMinus()) {
0361 fNkaons++;
0362 fHisto->Fill(8, e, 1.0);
0363 }
0364 else if (pd == G4KaonZeroShort::KaonZeroShort() || pd == G4KaonZeroLong::KaonZeroLong()) {
0365 fNkaons++;
0366 fHisto->Fill(9, e, 1.0);
0367 }
0368 else if (pd == G4Deuteron::Deuteron() || pd == G4Triton::Triton()) {
0369 fNdeut++;
0370 fHisto->Fill(10, e, 1.0);
0371 }
0372 else if (pd == G4He3::He3() || pd == G4Alpha::Alpha()) {
0373 fNalpha++;
0374 fHisto->Fill(11, e, 1.0);
0375 }
0376 else if (pd->GetParticleType() == "nucleus") {
0377 fNions++;
0378 fHisto->Fill(12, e, 1.0);
0379 }
0380 else if (pd == G4MuonPlus::MuonPlus() || pd == G4MuonMinus::MuonMinus()) {
0381 fNmuons++;
0382 fHisto->Fill(13, e, 1.0);
0383 }
0384 }
0385 }
0386
0387
0388
0389 void HistoManager::AddTargetStep(const G4Step* step)
0390 {
0391 ++fNstep;
0392 const G4double fEdep = step->GetTotalEnergyDeposit();
0393 const G4Track* track = step->GetTrack();
0394 const G4ParticleDefinition* pd = track->GetDefinition();
0395 if (1 < fVerbose) {
0396 G4cout << "TargetSD::ProcessHits: " << pd->GetParticleName()
0397 << " E(MeV)=" << track->GetKineticEnergy() / CLHEP::MeV
0398 << " beta1= " << step->GetPreStepPoint()->GetVelocity() / CLHEP::c_light
0399 << " beta2= " << step->GetPostStepPoint()->GetVelocity() / CLHEP::c_light
0400 << " weight= " << step->GetTrack()->GetWeight()
0401 << " t(ns)=" << track->GetGlobalTime() / CLHEP::ns << G4endl;
0402 }
0403 if (fEdep > 0.0) {
0404 G4ThreeVector pos =
0405 (step->GetPreStepPoint()->GetPosition() + step->GetPostStepPoint()->GetPosition()) * 0.5;
0406
0407 G4double z = pos.z() + fAbsZ0;
0408
0409
0410 fEdepEvt += fEdep;
0411 fHisto->Fill(0, z, fEdep);
0412
0413 if (pd == G4Gamma::Gamma() || pd == G4Electron::Electron() || pd == G4Positron::Positron()) {
0414 fEdepEM += fEdep;
0415 }
0416 else if (pd == G4PionPlus::PionPlus() || pd == G4PionMinus::PionMinus()) {
0417 fEdepPI += fEdep;
0418 }
0419 else if (pd == G4Proton::Proton() || pd == G4AntiProton::AntiProton()) {
0420 fEdepP += fEdep;
0421 }
0422
0423 if (1 < fVerbose) {
0424 G4cout << "HistoManager::AddEnergy: e(keV)= " << fEdep / keV << "; z(mm)= " << z / mm
0425 << "; step(mm)= " << step->GetStepLength() / mm << " by " << pd->GetParticleName()
0426 << " E(MeV)= " << track->GetKineticEnergy() / MeV << G4endl;
0427 }
0428 }
0429 }
0430
0431
0432
0433 void HistoManager::AddLeakingParticle(const G4Track* track)
0434 {
0435 const G4ParticleDefinition* pd = track->GetDefinition();
0436 const G4StepPoint* sp = track->GetStep()->GetPreStepPoint();
0437 const G4double ekin = sp->GetKineticEnergy();
0438 G4double e = -100.;
0439 if (ekin > 0.0) {
0440 e = std::log10(ekin / CLHEP::MeV);
0441 }
0442 else {
0443 G4cout << "### HistoManager::AddLeakingParticle " << pd->GetParticleName()
0444 << " Eprestep(MeV)=" << ekin << " trackID=" << track->GetTrackID() << G4endl;
0445 return;
0446 }
0447
0448 const G4ThreeVector& pos = sp->GetPosition();
0449 const G4ThreeVector& dir = sp->GetMomentumDirection();
0450 G4double x = pos.x();
0451 G4double y = pos.y();
0452 G4double z = pos.z();
0453 G4double vx = dir.x();
0454 G4double vy = dir.y();
0455 G4double vz = dir.z();
0456
0457 G4bool isLeaking = false;
0458
0459
0460 const G4double del = 0.001 * CLHEP::mm;
0461 if (std::abs(z - fAbsZ0) < del && vz > 0.0) {
0462 isLeaking = true;
0463 if (pd == fNeutron) {
0464 ++fNneu_forw;
0465 fHisto->Fill(15, e, 1.0);
0466 }
0467 else
0468 isLeaking = true;
0469
0470
0471 }
0472 else if (std::abs(z + fAbsZ0) < del && vz < 0.0) {
0473 isLeaking = true;
0474 if (pd == fNeutron) {
0475 ++fNneu_back;
0476 fHisto->Fill(16, e, 1.0);
0477 }
0478 else
0479 isLeaking = true;
0480
0481
0482 }
0483 else if (std::abs(z) <= fAbsZ0 + del && x * vx + y * vy > 0.0
0484 && std::abs(x * x + y * y - fR2) < del * fRadius)
0485 {
0486 isLeaking = true;
0487 if (pd == fNeutron) {
0488 ++fNneu_leak;
0489 fHisto->Fill(14, e, 1.0);
0490 }
0491 else
0492 isLeaking = true;
0493 }
0494
0495
0496 if (isLeaking) {
0497 if (pd == G4Proton::Proton()) {
0498 fHisto->Fill(17, e, 1.0);
0499 ++fNprot_leak;
0500 }
0501 else if (pd == G4PionPlus::PionPlus() || pd == G4PionMinus::PionMinus()) {
0502 fHisto->Fill(18, e, 1.0);
0503 ++fNpiofNleak;
0504 }
0505 }
0506 }
0507
0508
0509
0510 void HistoManager::SetVerbose(G4int val)
0511 {
0512 fVerbose = val;
0513 fHisto->SetVerbose(val);
0514 }
0515
0516
0517
0518 void HistoManager::Fill(G4int id, G4double x, G4double w)
0519 {
0520 fHisto->Fill(id, x, w);
0521 }
0522
0523
