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