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