File indexing completed on 2026-09-17 08:31:29
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
0047 #include "HistoManager.hh"
0048
0049 #include "DetectorConstruction.hh"
0050 #include "Histo.hh"
0051 #include "IonHIJINGPhysics.hh"
0052 #include "IonUrQMDPhysics.hh"
0053
0054 #include "G4Alpha.hh"
0055 #include "G4AntiProton.hh"
0056 #include "G4BuilderType.hh"
0057 #include "G4Deuteron.hh"
0058 #include "G4Electron.hh"
0059 #include "G4Gamma.hh"
0060 #include "G4He3.hh"
0061 #include "G4KaonMinus.hh"
0062 #include "G4KaonPlus.hh"
0063 #include "G4KaonZeroLong.hh"
0064 #include "G4KaonZeroShort.hh"
0065 #include "G4Material.hh"
0066 #include "G4MuonMinus.hh"
0067 #include "G4MuonPlus.hh"
0068 #include "G4Neutron.hh"
0069 #include "G4NistManager.hh"
0070 #include "G4PionMinus.hh"
0071 #include "G4PionPlus.hh"
0072 #include "G4PionZero.hh"
0073 #include "G4Positron.hh"
0074 #include "G4Proton.hh"
0075 #include "G4RunManager.hh"
0076 #include "G4SystemOfUnits.hh"
0077 #include "G4Triton.hh"
0078 #include "G4UnitsTable.hh"
0079 #include "G4VModularPhysicsList.hh"
0080 #include "G4VPhysicsConstructor.hh"
0081 #include "globals.hh"
0082
0083
0084
0085 HistoManager* HistoManager::fManager = 0;
0086
0087
0088
0089 HistoManager* HistoManager::GetPointer()
0090 {
0091 if (!fManager) {
0092 static HistoManager manager;
0093 fManager = &manager;
0094 }
0095 return fManager;
0096 }
0097
0098
0099
0100 HistoManager::HistoManager()
0101 {
0102 fVerbose = 0;
0103 fNSlices = 100;
0104 fNBinsE = 100;
0105 fNHisto = 20;
0106 fLength = 300. * mm;
0107 fEdepMax = 1.0 * GeV;
0108
0109 fPrimaryDef = 0;
0110 fPrimaryKineticEnergy = 0.0;
0111 fMaterial = 0;
0112 fBeamFlag = true;
0113
0114 fHisto = new Histo();
0115 fHisto->SetVerbose(fVerbose);
0116 fNeutron = G4Neutron::Neutron();
0117 fPhysList = 0;
0118 fIonPhysics = 0;
0119 Initialise();
0120 }
0121
0122
0123
0124 void HistoManager::Initialise()
0125 {
0126 fAbsZ0 = -0.5 * fLength;
0127 fNevt = 0;
0128 fNelec = 0;
0129 fNposit = 0;
0130 fNgam = 0;
0131 fNstep = 0;
0132 fNions = 0;
0133 fNdeut = 0;
0134 fNalpha = 0;
0135 fNkaons = 0;
0136 fNmuons = 0;
0137 fNcpions = 0;
0138 fNpi0 = 0;
0139 fNneutron = 0;
0140 fNproton = 0;
0141 fNaproton = 0;
0142
0143 fEdepEvt = 0.0;
0144 fEdepSum = 0.0;
0145 fEdepSum2 = 0.0;
0146 }
0147
0148
0149
0150 HistoManager::~HistoManager()
0151 {
0152 delete fHisto;
0153 }
0154
0155
0156
0157 void HistoManager::BookHisto()
0158 {
0159 fHisto->Add1D("0", "Energy deposition (MeV/mm/event) in the target", fNSlices, 0.0, fLength / mm,
0160 MeV / mm);
0161 fHisto->Add1D("1", "Log10 Energy (GeV) of gammas", fNBinsE, -5., 5., 1.0);
0162 fHisto->Add1D("2", "Log10 Energy (GeV) of electrons", fNBinsE, -5., 5., 1.0);
0163 fHisto->Add1D("3", "Log10 Energy (GeV) of positrons", fNBinsE, -5., 5., 1.0);
0164 fHisto->Add1D("4", "Log10 Energy (GeV) of protons", fNBinsE, -5., 5., 1.0);
0165 fHisto->Add1D("5", "Log10 Energy (GeV) of neutrons", fNBinsE, -5., 5., 1.0);
0166 fHisto->Add1D("6", "Log10 Energy (GeV) of charged pions", fNBinsE, -4., 6., 1.0);
0167 fHisto->Add1D("7", "Log10 Energy (GeV) of pi0", fNBinsE, -4., 6., 1.0);
0168 fHisto->Add1D("8", "Log10 Energy (GeV) of charged kaons", fNBinsE, -4., 6., 1.0);
0169 fHisto->Add1D("9", "Log10 Energy (GeV) of neutral kaons", fNBinsE, -4., 6., 1.0);
0170 fHisto->Add1D("10", "Log10 Energy (GeV) of deuterons and tritons", fNBinsE, -5., 5., 1.0);
0171 fHisto->Add1D("11", "Log10 Energy (GeV) of He3 and alpha", fNBinsE, -5., 5., 1.0);
0172 fHisto->Add1D("12", "Log10 Energy (GeV) of Generic Ions", fNBinsE, -5., 5., 1.0);
0173 fHisto->Add1D("13", "Log10 Energy (GeV) of muons", fNBinsE, -4., 6., 1.0);
0174 fHisto->Add1D("14", "Log10 Energy (GeV) of pi+", fNBinsE, -4., 6., 1.0);
0175 fHisto->Add1D("15", "Log10 Energy (GeV) of pi-", fNBinsE, -4., 6., 1.0);
0176 fHisto->Add1D("16", "X Section (mb) of Secondary Fragments Z with E>1 GeV (mb)", 25, 0.5, 25.5,
0177 1.0);
0178 fHisto->Add1D("17", "Secondary Fragment A E>1 GeV", 50, 0.5, 50.5, 1.0);
0179 fHisto->Add1D("18", "Secondary Fragment Z E<1 GeV", 25, 0.5, 25.5, 1.0);
0180 fHisto->Add1D("19", "Secondary Fragment A E<1 GeV", 50, 0.5, 50.5, 1.0);
0181 fHisto->Add1D("20", "X Section (mb) of Secondary Fragments Z (mb) ", 25, 0.5, 25.5, 1.0);
0182 fHisto->Add1D("21", "Secondary Fragment A ", 50, 0.5, 50.5, 1.0);
0183 }
0184
0185
0186
0187 void HistoManager::BeginOfRun()
0188 {
0189 Initialise();
0190 BookHisto();
0191 fHisto->Book();
0192
0193 if (fVerbose > 0) {
0194 G4cout << "HistoManager: Histograms are booked and run has been started" << G4endl
0195 << " BeginOfRun (After fHisto->book)" << G4endl;
0196 }
0197 }
0198
0199
0200
0201 void HistoManager::EndOfRun()
0202 {
0203 G4cout << "HistoManager: End of run actions are started" << G4endl;
0204
0205
0206 G4cout << "========================================================" << G4endl;
0207
0208 G4double x = (G4double)fNevt;
0209 if (fNevt > 0) {
0210 x = 1.0 / x;
0211 }
0212
0213 G4double xe = x * fNelec;
0214 G4double xg = x * fNgam;
0215 G4double xp = x * fNposit;
0216 G4double xs = x * fNstep;
0217 G4double xn = x * fNneutron;
0218 G4double xpn = x * fNproton;
0219 G4double xap = x * fNaproton;
0220 G4double xpc = x * fNcpions;
0221 G4double xp0 = x * fNpi0;
0222 G4double xpk = x * fNkaons;
0223 G4double xpm = x * fNmuons;
0224 G4double xid = x * fNdeut;
0225 G4double xia = x * fNalpha;
0226 G4double xio = x * fNions;
0227
0228 fEdepSum *= x;
0229 fEdepSum2 *= x;
0230 fEdepSum2 -= fEdepSum * fEdepSum;
0231 if (fEdepSum2 > 0.0)
0232 fEdepSum2 = std::sqrt(fEdepSum2);
0233 else
0234 fEdepSum2 = 0.0;
0235
0236 G4cout << "Beam particle " << fPrimaryDef->GetParticleName() << G4endl;
0237 G4cout << "Beam Energy(GeV) " << fPrimaryKineticEnergy / GeV << G4endl;
0238 G4cout << "Number of events " << fNevt << G4endl;
0239 G4cout << std::setprecision(4) << "Average energy deposit (GeV) " << fEdepSum / GeV
0240 << " RMS(GeV) " << fEdepSum2 / GeV << G4endl;
0241 G4cout << std::setprecision(4) << "Average number of steps " << xs << G4endl;
0242 G4cout << std::setprecision(4) << "Average number of gamma " << xg << G4endl;
0243 G4cout << std::setprecision(4) << "Average number of e- " << xe << G4endl;
0244 G4cout << std::setprecision(4) << "Average number of e+ " << xp << G4endl;
0245 G4cout << std::setprecision(4) << "Average number of neutrons " << xn << G4endl;
0246 G4cout << std::setprecision(4) << "Average number of protons " << xpn << G4endl;
0247 G4cout << std::setprecision(4) << "Average number of antiprotons " << xap << G4endl;
0248 G4cout << std::setprecision(4) << "Average number of pi+ & pi- " << xpc << G4endl;
0249 G4cout << std::setprecision(4) << "Average number of pi0 " << xp0 << G4endl;
0250 G4cout << std::setprecision(4) << "Average number of kaons " << xpk << G4endl;
0251 G4cout << std::setprecision(4) << "Average number of muons " << xpm << G4endl;
0252 G4cout << std::setprecision(4) << "Average number of deuterons+tritons " << xid << G4endl;
0253 G4cout << std::setprecision(4) << "Average number of He3+alpha " << xia << G4endl;
0254 G4cout << std::setprecision(4) << "Average number of ions " << xio << G4endl;
0255 G4cout << "========================================================" << G4endl;
0256 G4cout << G4endl;
0257
0258
0259 for (G4int i = 0; i < fNHisto; i++) {
0260 fHisto->ScaleH1(i, x);
0261 }
0262
0263
0264 G4double F = 1000 / (fLength * barn * fMaterial->GetTotNbOfAtomsPerVolume());
0265 if (F > 0.0) {
0266 fHisto->ScaleH1(16, F);
0267 fHisto->ScaleH1(20, F);
0268 }
0269
0270 fHisto->Save();
0271 }
0272
0273
0274
0275 void HistoManager::BeginOfEvent()
0276 {
0277 ++fNevt;
0278 fEdepEvt = 0.0;
0279 }
0280
0281
0282
0283 void HistoManager::EndOfEvent()
0284 {
0285 fEdepSum += fEdepEvt;
0286 fEdepSum2 += fEdepEvt * fEdepEvt;
0287 }
0288
0289
0290
0291 void HistoManager::ScoreNewTrack(const G4Track* track)
0292 {
0293 const G4ParticleDefinition* pd = track->GetDefinition();
0294 G4String name = pd->GetParticleName();
0295 G4double e = track->GetKineticEnergy();
0296
0297
0298 if (0 == track->GetParentID()) {
0299 fPrimaryKineticEnergy = e;
0300 fPrimaryDef = pd;
0301 G4ThreeVector dir = track->GetMomentumDirection();
0302 if (1 < fVerbose)
0303 G4cout << "### Primary " << name << " kinE(GeV)= " << e / GeV
0304 << "; m(GeV)= " << pd->GetPDGMass() / GeV << "; pos(mm)= " << track->GetPosition() / mm
0305 << "; dir= " << track->GetMomentumDirection() << G4endl;
0306
0307
0308 }
0309 else {
0310 if (1 < fVerbose) {
0311 G4cout << "=== Secondary " << name << " kinE(GeV)= " << e / GeV
0312 << "; m(GeV)= " << pd->GetPDGMass() / GeV << "; pos(mm)= " << track->GetPosition() / mm
0313 << "; dir= " << track->GetMomentumDirection() << G4endl;
0314 }
0315 e = std::log10(e / GeV);
0316 if (pd == G4Gamma::Gamma()) {
0317 fNgam++;
0318 fHisto->Fill(1, e, 1.0);
0319 }
0320 else if (pd == G4Electron::Electron()) {
0321 fNelec++;
0322 fHisto->Fill(2, e, 1.0);
0323 }
0324 else if (pd == G4Positron::Positron()) {
0325 fNposit++;
0326 fHisto->Fill(3, e, 1.0);
0327 }
0328 else if (pd == G4Proton::Proton()) {
0329 fNproton++;
0330 fHisto->Fill(4, e, 1.0);
0331 }
0332 else if (pd == fNeutron) {
0333 fNneutron++;
0334 fHisto->Fill(5, e, 1.0);
0335 }
0336 else if (pd == G4AntiProton::AntiProton()) {
0337 fNaproton++;
0338 }
0339 else if (pd == G4PionPlus::PionPlus()) {
0340 fNcpions++;
0341 fHisto->Fill(6, e, 1.0);
0342 fHisto->Fill(14, e, 1.0);
0343 }
0344 else if (pd == G4PionMinus::PionMinus()) {
0345 fNcpions++;
0346 fHisto->Fill(6, e, 1.0);
0347 fHisto->Fill(15, e, 1.0);
0348 }
0349 else if (pd == G4PionZero::PionZero()) {
0350 fNpi0++;
0351 fHisto->Fill(7, e, 1.0);
0352 }
0353 else if (pd == G4KaonPlus::KaonPlus() || pd == G4KaonMinus::KaonMinus()) {
0354 fNkaons++;
0355 fHisto->Fill(8, e, 1.0);
0356 }
0357 else if (pd == G4KaonZeroShort::KaonZeroShort() || pd == G4KaonZeroLong::KaonZeroLong()) {
0358 fNkaons++;
0359 fHisto->Fill(9, e, 1.0);
0360 }
0361 else if (pd == G4Deuteron::Deuteron() || pd == G4Triton::Triton()) {
0362 fNdeut++;
0363 fHisto->Fill(10, e, 1.0);
0364 }
0365 else if (pd == G4He3::He3() || pd == G4Alpha::Alpha()) {
0366 fNalpha++;
0367 fHisto->Fill(11, e, 1.0);
0368 }
0369 else if (pd->GetParticleType() == "nucleus") {
0370 fNions++;
0371 fHisto->Fill(12, e, 1.0);
0372 G4double Z = pd->GetPDGCharge() / eplus;
0373 G4double A = (G4double)pd->GetBaryonNumber();
0374 if (e > 0.0) {
0375 fHisto->Fill(16, Z, 1.0);
0376 fHisto->Fill(17, A, 1.0);
0377 }
0378 else {
0379 fHisto->Fill(18, Z, 1.0);
0380 fHisto->Fill(19, A, 1.0);
0381 }
0382 fHisto->Fill(20, Z, 1.0);
0383 fHisto->Fill(21, A, 1.0);
0384 }
0385 else if (pd == G4MuonPlus::MuonPlus() || pd == G4MuonMinus::MuonMinus()) {
0386 fNmuons++;
0387 fHisto->Fill(13, e, 1.0);
0388 }
0389 }
0390 }
0391
0392
0393
0394 void HistoManager::AddTargetStep(const G4Step* step)
0395 {
0396 ++fNstep;
0397 G4double edep = step->GetTotalEnergyDeposit();
0398
0399 if (edep > 0.0) {
0400 G4ThreeVector pos =
0401 (step->GetPreStepPoint()->GetPosition() + step->GetPostStepPoint()->GetPosition()) * 0.5;
0402
0403 G4double z = pos.z() - fAbsZ0;
0404
0405
0406 fEdepEvt += edep;
0407 fHisto->Fill(0, z, edep);
0408
0409 if (1 < fVerbose) {
0410 const G4Track* track = step->GetTrack();
0411 G4cout << "HistoManager::AddEnergy: e(keV)= " << edep / keV << "; z(mm)= " << z / mm
0412 << "; step(mm)= " << step->GetStepLength() / mm << " by "
0413 << track->GetDefinition()->GetParticleName()
0414 << " E(GeV)= " << track->GetKineticEnergy() / GeV << G4endl;
0415 }
0416 }
0417 }
0418
0419
0420
0421 void HistoManager::SetVerbose(G4int val)
0422 {
0423 fVerbose = val;
0424 fHisto->SetVerbose(val);
0425 }
0426
0427
0428
0429 void HistoManager::Fill(G4int id, G4double x, G4double w)
0430 {
0431 fHisto->Fill(id, x, w);
0432 }
0433
0434
0435
0436 void HistoManager::SetIonPhysics(const G4String& nam)
0437 {
0438 if (fIonPhysics) {
0439 G4cout << "### HistoManager WARNING: Ion Physics is already defined: <" << nam
0440 << "> is ignored!" << G4endl;
0441 }
0442 else if (nam == "HIJING") {
0443 #ifdef G4_USE_HIJING
0444 fIonPhysics = new IonHIJINGPhysics();
0445 fPhysList->ReplacePhysics(fIonPhysics);
0446 G4RunManager::GetRunManager()->PhysicsHasBeenModified();
0447 G4cout << "### SetIonPhysics: Ion Physics FTFP/Binary is added" << G4endl;
0448 #else
0449 G4cout << "Error: Ion Physics HIJING is requested but is not available" << G4endl;
0450 #endif
0451 }
0452 else if (nam == "UrQMD") {
0453 #ifdef G4_USE_URQMD
0454 fIonPhysics = new IonUrQMDPhysics(1);
0455 fPhysList->ReplacePhysics(fIonPhysics);
0456 G4RunManager::GetRunManager()->PhysicsHasBeenModified();
0457 G4cout << "### SetIonPhysics: Ion Physics UrQMD is added" << G4endl;
0458 #else
0459 G4cout << "Error: Ion Physics UrQMD is requested but is not available" << G4endl;
0460 #endif
0461 }
0462 else {
0463 G4cout << "### HistoManager WARNING: Ion Physics <" << nam << "> is unknown!" << G4endl;
0464 }
0465 }
0466
0467
