File indexing completed on 2026-03-31 07:49:58
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 #include "HistoManager.hh"
0043
0044 #include "EmAcceptance.hh"
0045 #include "Histo.hh"
0046
0047 #include "G4Electron.hh"
0048 #include "G4EmProcessSubType.hh"
0049 #include "G4Gamma.hh"
0050 #include "G4GammaGeneralProcess.hh"
0051 #include "G4MaterialCutsCouple.hh"
0052 #include "G4Positron.hh"
0053 #include "G4SystemOfUnits.hh"
0054 #include "G4UnitsTable.hh"
0055 #include "G4VEmProcess.hh"
0056 #include "G4VEnergyLossProcess.hh"
0057 #include "G4VProcess.hh"
0058
0059
0060
0061 HistoManager* HistoManager::fManager = nullptr;
0062
0063
0064
0065 HistoManager* HistoManager::GetPointer()
0066 {
0067 if (nullptr == fManager) {
0068 fManager = new HistoManager();
0069 }
0070 return fManager;
0071 }
0072
0073
0074
0075 HistoManager::HistoManager()
0076 : fGamma(G4Gamma::Gamma()),
0077 fElectron(G4Electron::Electron()),
0078 fPositron(G4Positron::Positron()),
0079 fHisto(new Histo())
0080 {
0081 fVerbose = 1;
0082 fEvt1 = -1;
0083 fEvt2 = -1;
0084 fNmax = 3;
0085 fMaxEnergy = 50.0 * MeV;
0086 fBeamEnergy = 1. * GeV;
0087 fMaxEnergyAbs = 10. * MeV;
0088 fBinsE = 100;
0089 fBinsEA = 40;
0090 fBinsED = 100;
0091 fNHisto = 13;
0092
0093 BookHisto();
0094 }
0095
0096
0097
0098 HistoManager::~HistoManager()
0099 {
0100 delete fHisto;
0101 }
0102
0103
0104
0105 void HistoManager::BookHisto()
0106 {
0107 fHisto->Add1D("10", "Evis/E0 in central crystal", fBinsED, 0.0, 1, 1.0);
0108 fHisto->Add1D("11", "Evis/E0 in 3x3", fBinsED, 0.0, 1.0, 1.0);
0109 fHisto->Add1D("12", "Evis/E0 in 5x5", fBinsED, 0.0, 1.0, 1.0);
0110 fHisto->Add1D("13", "Energy (MeV) of delta-electrons", fBinsE, 0.0, fMaxEnergy, MeV);
0111 fHisto->Add1D("14", "Energy (MeV) of gammas", fBinsE, 0.0, fMaxEnergy, MeV);
0112 fHisto->Add1D("15", "Energy (MeV) in abs1", fBinsEA, 0.0, fMaxEnergyAbs, MeV);
0113 fHisto->Add1D("16", "Energy (MeV) in abs2", fBinsEA, 0.0, fMaxEnergyAbs, MeV);
0114 fHisto->Add1D("17", "Energy (MeV) in abs3", fBinsEA, 0.0, fMaxEnergyAbs, MeV);
0115 fHisto->Add1D("18", "Energy (MeV) in abs4", fBinsEA, 0.0, fMaxEnergyAbs, MeV);
0116 fHisto->Add1D("19", "Number of vertex hits", 20, -0.5, 19.5, 1.0);
0117 fHisto->Add1D("20", "E1/E9 Ratio", fBinsED, 0.0, 1, 1.0);
0118 fHisto->Add1D("21", "E1/E25 Ratio", fBinsED, 0.0, 1.0, 1.0);
0119 fHisto->Add1D("22", "E9/E25 Ratio", fBinsED, 0.0, 1.0, 1.0);
0120 }
0121
0122
0123
0124 void HistoManager::BeginOfRun()
0125 {
0126
0127 fEvt = 0;
0128 fElec = 0;
0129 fPosit = 0;
0130 fGam = 0;
0131 fStep = 0;
0132 fLowe = 0;
0133
0134 for (G4int i = 0; i < 6; i++) {
0135 fStat[i] = 0;
0136 fEdep[i] = 0.0;
0137 fErms[i] = 0.0;
0138 if (i < 3) {
0139 fEdeptr[i] = 0.0;
0140 fErmstr[i] = 0.0;
0141 }
0142 }
0143
0144
0145 fBrem.resize(93, 0.0);
0146 fPhot.resize(93, 0.0);
0147 fComp.resize(93, 0.0);
0148 fConv.resize(93, 0.0);
0149
0150
0151 for (G4int i = 0; i < fNmax; i++) {
0152 fEdeptrue[i] = 1.0;
0153 fRmstrue[i] = 1.0;
0154 fLimittrue[i] = 10.;
0155 }
0156
0157 if (fHisto->IsActive()) {
0158 for (G4int i = 0; i < fNHisto; ++i) {
0159 fHisto->Activate(i, true);
0160 }
0161 fHisto->Book();
0162
0163 if (fVerbose > 0) {
0164 G4cout << "HistoManager: Histograms are booked and run has been started" << G4endl;
0165 }
0166 }
0167 }
0168
0169
0170
0171 void HistoManager::EndOfRun(G4int runID)
0172 {
0173 G4cout << "HistoManager: End of run actions are started RunID# " << runID << G4endl;
0174 G4String nam[6] = {"1x1", "3x3", "5x5", "E1/E9 ", "E1/E25", "E9/E25"};
0175
0176
0177
0178 G4cout << "=================================================================" << G4endl;
0179 G4double x = (G4double)fEvt;
0180 if (fEvt > 0) x = 1.0 / x;
0181 G4int j;
0182 for (j = 0; j < fNmax; j++) {
0183
0184 fEdep[j] *= x;
0185 G4double y = fErms[j] * x - fEdep[j] * fEdep[j];
0186 if (y < 0.0) y = 0.0;
0187 fErms[j] = std::sqrt(y);
0188
0189
0190 G4double xx = G4double(fStat[j]);
0191 if (xx > 0.0) xx = 1.0 / xx;
0192 fEdeptr[j] *= xx;
0193 y = fErmstr[j] * xx - fEdeptr[j] * fEdeptr[j];
0194 if (y < 0.0) y = 0.0;
0195 fErmstr[j] = std::sqrt(y);
0196 }
0197 G4double xe = x * (G4double)fElec;
0198 G4double xg = x * (G4double)fGam;
0199 G4double xp = x * (G4double)fPosit;
0200 G4double xs = x * fStep;
0201
0202 G4double f = 100. * std::sqrt(fBeamEnergy / GeV);
0203
0204 G4cout << "Number of events " << fEvt << G4endl;
0205 G4cout << std::setprecision(4) << "Average number of e- " << xe << G4endl;
0206 G4cout << std::setprecision(4) << "Average number of gamma " << xg << G4endl;
0207 G4cout << std::setprecision(4) << "Average number of e+ " << xp << G4endl;
0208 G4cout << std::setprecision(4) << "Average number of steps " << xs << G4endl;
0209
0210 for (j = 0; j < 3; ++j) {
0211 G4double ex = fEdeptr[j];
0212 G4double sx = fErmstr[j];
0213 G4double xx = G4double(fStat[j]);
0214 if (xx > 0.0) xx = 1.0 / xx;
0215 G4double r = sx * std::sqrt(xx);
0216 G4cout << std::setprecision(4) << "Edep " << nam[j] << " = " << ex << " +- "
0217 << r;
0218 if (ex > 0.1) G4cout << " res= " << f * sx / ex << " % " << fStat[j];
0219 G4cout << G4endl;
0220 }
0221 if (fLimittrue[0] < 10. || fLimittrue[1] < 10. || fLimittrue[2] < 10.) {
0222 G4cout << "=========== Mean values without trancating =====================" << G4endl;
0223 for (j = 0; j < fNmax; j++) {
0224 G4double ex = fEdep[j];
0225 G4double sx = fErms[j];
0226 G4double rx = sx * std::sqrt(x);
0227 G4cout << std::setprecision(4) << "Edep " << nam[j] << " = " << ex << " +- "
0228 << rx;
0229 if (ex > 0.0) G4cout << " res= " << f * sx / ex << " %";
0230 G4cout << G4endl;
0231 }
0232 }
0233 G4cout << "=========== Ratios without trancating ===========================" << G4endl;
0234 for (j = 3; j < 6; ++j) {
0235 G4double e = fEdep[j];
0236 G4double xx = G4double(fStat[j]);
0237 if (xx > 0.0) xx = 1.0 / xx;
0238 e *= xx;
0239 G4double y = fErms[j] * xx - e * e;
0240 G4double r = 0.0;
0241 if (y > 0.0) r = std::sqrt(y * xx);
0242 G4cout << " " << nam[j] << " = " << e << " +- " << r;
0243 G4cout << G4endl;
0244 }
0245 G4cout << std::setprecision(4) << "Beam Energy " << fBeamEnergy / GeV << " GeV"
0246 << G4endl;
0247 if (fLowe > 0) G4cout << "Number of events E/E0<0.8 " << fLowe << G4endl;
0248 G4cout << "==================================================================" << G4endl;
0249 G4cout << G4endl;
0250
0251
0252 if (fHisto->IsActive()) {
0253 for (G4int i = 0; i < fNHisto; ++i) {
0254 fHisto->ScaleH1(i, x);
0255 }
0256 fHisto->Save();
0257 }
0258 if (0 < runID) {
0259 return;
0260 }
0261
0262
0263 EmAcceptance acc;
0264 G4bool isStarted = false;
0265 for (j = 0; j < fNmax; j++) {
0266 G4double ltrue = fLimittrue[j];
0267 if (ltrue < DBL_MAX) {
0268 if (!isStarted) {
0269 acc.BeginOfAcceptance("Crystal Calorimeter", fEvt);
0270 isStarted = true;
0271 }
0272 G4double etrue = fEdeptrue[j];
0273 G4double rtrue = fRmstrue[j];
0274 acc.EmAcceptanceGauss("Edep" + nam[j], fEvt, fEdeptr[j], etrue, rtrue, ltrue);
0275 acc.EmAcceptanceGauss("Erms" + nam[j], fEvt, fErmstr[j], rtrue, rtrue, 2.0 * ltrue);
0276 }
0277 }
0278 if (isStarted) acc.EndOfAcceptance();
0279
0280
0281 G4cout << " Z bremsstrahlung photoeffect compton conversion" << G4endl;
0282 for (j = 1; j < 93; ++j) {
0283 G4int n1 = G4int(fBrem[j] * x);
0284 G4int n2 = G4int(fPhot[j] * x);
0285 G4int n3 = G4int(fComp[j] * x);
0286 G4int n4 = G4int(fConv[j] * x);
0287 if (n1 + n2 + n3 + n4 > 0) {
0288 G4cout << std::setw(4) << j << std::setw(12) << n1 << std::setw(12) << n2 << std::setw(12)
0289 << n3 << std::setw(12) << n4 << G4endl;
0290 }
0291 }
0292 }
0293
0294
0295
0296 void HistoManager::BeginOfEvent()
0297 {
0298 ++fEvt;
0299
0300 fEabs1 = 0.0;
0301 fEabs2 = 0.0;
0302 fEabs3 = 0.0;
0303 fEabs4 = 0.0;
0304 fEvertex.clear();
0305 fNvertex.clear();
0306 for (G4int i = 0; i < 25; i++) {
0307 fE[i] = 0.0;
0308 }
0309 }
0310
0311
0312
0313 void HistoManager::EndOfEvent()
0314 {
0315 G4double e9 = 0.0;
0316 G4double e25 = 0.0;
0317 for (G4int i = 0; i < 25; i++) {
0318 fE[i] /= fBeamEnergy;
0319 e25 += fE[i];
0320 if ((6 <= i && 8 >= i) || (11 <= i && 13 >= i) || (16 <= i && 18 >= i)) e9 += fE[i];
0321 }
0322
0323 if (1 < fVerbose && e25 < 0.8) {
0324 ++fLowe;
0325 G4cout << "### in the event# " << fEvt << " E25= " << e25 << G4endl;
0326 }
0327
0328
0329 G4double e0 = fE[12];
0330 G4double e19 = 0.0;
0331 G4double e125 = 0.0;
0332 G4double e925 = 0.0;
0333 if (e9 > 0.0) {
0334 e19 = e0 / e9;
0335 e125 = e0 / e25;
0336 e925 = e9 / e25;
0337 fEdep[3] += e19;
0338 fErms[3] += e19 * e19;
0339 fEdep[4] += e125;
0340 fErms[4] += e125 * e125;
0341 fEdep[5] += e925;
0342 fErms[5] += e925 * e925;
0343 fStat[3] += 1;
0344 fStat[4] += 1;
0345 fStat[5] += 1;
0346 }
0347
0348
0349 fHisto->Fill(0, e0, 1.0);
0350 fHisto->Fill(1, e9, 1.0);
0351 fHisto->Fill(2, e25, 1.0);
0352 fHisto->Fill(5, fEabs1, 1.0);
0353 fHisto->Fill(6, fEabs2, 1.0);
0354 fHisto->Fill(7, fEabs3, 1.0);
0355 fHisto->Fill(8, fEabs4, 1.0);
0356 fHisto->Fill(9, G4double(fNvertex.size()), 1.0);
0357 fHisto->Fill(10, e19, 1.0);
0358 fHisto->Fill(11, e125, 1.0);
0359 fHisto->Fill(12, e925, 1.0);
0360
0361
0362 fEdep[0] += e0;
0363 fErms[0] += e0 * e0;
0364 fEdep[1] += e9;
0365 fErms[1] += e9 * e9;
0366 fEdep[2] += e25;
0367 fErms[2] += e25 * e25;
0368
0369
0370 if (std::abs(e0 - fEdeptrue[0]) < fRmstrue[0] * fLimittrue[0]) {
0371 fStat[0] += 1;
0372 fEdeptr[0] += e0;
0373 fErmstr[0] += e0 * e0;
0374 }
0375 if (std::abs(e9 - fEdeptrue[1]) < fRmstrue[1] * fLimittrue[1]) {
0376 fStat[1] += 1;
0377 fEdeptr[1] += e9;
0378 fErmstr[1] += e9 * e9;
0379 }
0380 if (std::abs(e25 - fEdeptrue[2]) < fRmstrue[2] * fLimittrue[2]) {
0381 fStat[2] += 1;
0382 fEdeptr[2] += e25;
0383 fErmstr[2] += e25 * e25;
0384 }
0385 }
0386
0387
0388
0389 void HistoManager::ScoreNewTrack(const G4Track* aTrack)
0390 {
0391
0392 ResetTrackLength();
0393 const G4ParticleDefinition* particle = aTrack->GetDefinition();
0394 const G4DynamicParticle* dynParticle = aTrack->GetDynamicParticle();
0395
0396 G4int pid = aTrack->GetParentID();
0397 G4double kinE = dynParticle->GetKineticEnergy();
0398 G4ThreeVector pos = aTrack->GetVertexPosition();
0399
0400
0401 if (0 == pid) {
0402 G4double mass = 0.0;
0403 if (particle) {
0404 mass = particle->GetPDGMass();
0405 }
0406
0407 G4ThreeVector dir = dynParticle->GetMomentumDirection();
0408 if (1 < fVerbose) {
0409 G4cout << "TrackingAction: Primary kinE(MeV)= " << kinE / MeV << "; m(MeV)= " << mass / MeV
0410 << "; pos= " << pos << "; dir= " << dir << G4endl;
0411 }
0412
0413
0414 }
0415 else {
0416 const G4VProcess* proc = aTrack->GetCreatorProcess();
0417 G4int type = proc->GetProcessSubType();
0418
0419 if (type == fBremsstrahlung) {
0420 auto elm = static_cast<const G4VEnergyLossProcess*>(proc)->GetCurrentElement();
0421 if (nullptr != elm) {
0422 G4int Z = elm->GetZasInt();
0423 if (Z > 0 && Z < 93) {
0424 fBrem[Z] += 1.0;
0425 }
0426 }
0427 }
0428 else if (type == fPhotoElectricEffect) {
0429 auto elm = static_cast<const G4VEmProcess*>(proc)->GetCurrentElement();
0430 if (nullptr != elm) {
0431 G4int Z = elm->GetZasInt();
0432 if (Z > 0 && Z < 93) {
0433 fPhot[Z] += 1.0;
0434 }
0435 }
0436 }
0437 else if (type == fGammaConversion) {
0438 auto elm = static_cast<const G4VEmProcess*>(proc)->GetCurrentElement();
0439 if (nullptr != elm) {
0440 G4int Z = elm->GetZasInt();
0441 if (Z > 0 && Z < 93) {
0442 fConv[Z] += 1.0;
0443 }
0444 }
0445 }
0446 else if (type == fComptonScattering) {
0447 auto elm = static_cast<const G4VEmProcess*>(proc)->GetCurrentElement();
0448 if (nullptr != elm) {
0449 G4int Z = elm->GetZasInt();
0450 if (Z > 0 && Z < 93) {
0451 fComp[Z] += 1.0;
0452 }
0453 }
0454 }
0455
0456
0457 if (particle == fElectron) {
0458 if (1 < fVerbose) {
0459 G4cout << "TrackingAction: Secondary electron " << G4endl;
0460 }
0461 AddDeltaElectron(dynParticle);
0462 }
0463 else if (particle == fPositron) {
0464 if (1 < fVerbose) {
0465 G4cout << "TrackingAction: Secondary positron " << G4endl;
0466 }
0467 AddPositron(dynParticle);
0468 }
0469 else if (particle == fGamma) {
0470 if (1 < fVerbose) {
0471 G4cout << "TrackingAction: Secondary gamma; parentID= " << pid << " E= " << kinE << G4endl;
0472 }
0473 AddPhoton(dynParticle);
0474 }
0475 }
0476 }
0477
0478
0479
0480 void HistoManager::AddEnergy(G4double edep, G4int volIndex, G4int copyNo)
0481 {
0482 if (1 < fVerbose) {
0483 G4cout << "HistoManager::AddEnergy: e(keV)= " << edep / keV << "; volIdx= " << volIndex
0484 << "; copyNo= " << copyNo << G4endl;
0485 }
0486 if (0 == volIndex) {
0487 fE[copyNo] += edep;
0488 }
0489 else if (1 == volIndex) {
0490 fEabs1 += edep;
0491 }
0492 else if (2 == volIndex) {
0493 fEabs2 += edep;
0494 }
0495 else if (3 == volIndex) {
0496 fEabs3 += edep;
0497 }
0498 else if (4 == volIndex) {
0499 fEabs4 += edep;
0500 }
0501 else if (5 == volIndex) {
0502 G4int n = fNvertex.size();
0503 G4bool newPad = true;
0504 if (n > 0) {
0505 for (G4int i = 0; i < n; i++) {
0506 if (copyNo == fNvertex[i]) {
0507 newPad = false;
0508 fEvertex[i] += edep;
0509 break;
0510 }
0511 }
0512 }
0513 if (newPad) {
0514 fNvertex.push_back(copyNo);
0515 fEvertex.push_back(edep);
0516 }
0517 }
0518 }
0519
0520
0521
0522 void HistoManager::AddDeltaElectron(const G4DynamicParticle* elec)
0523 {
0524 G4double e = elec->GetKineticEnergy() / MeV;
0525 if (e > 0.0) {
0526 ++fElec;
0527 fHisto->Fill(3, e, 1.0);
0528 }
0529 }
0530
0531
0532
0533 void HistoManager::AddPhoton(const G4DynamicParticle* ph)
0534 {
0535 G4double e = ph->GetKineticEnergy() / MeV;
0536 if (e > 0.0) {
0537 ++fGam;
0538 fHisto->Fill(4, e, 1.0);
0539 }
0540 }
0541
0542
0543
0544 void HistoManager::SetEdepAndRMS(G4int i, const G4ThreeVector& val)
0545 {
0546 if (i < fNmax && i >= 0) {
0547 if (val[0] > 0.0) fEdeptrue[i] = val[0];
0548 if (val[1] > 0.0) fRmstrue[i] = val[1];
0549 if (val[2] > 0.0) fLimittrue[i] = val[2];
0550 }
0551 }
0552
0553
