Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 /// \file electromagnetic/TestEm9/src/HistoMan 27 /// \brief Implementation of the HistoManager 28 // 29 // 30 //-------------------------------------------- 31 // 32 // ClassName: HistoManager 33 // 34 // 35 // Author: V.Ivanchenko 30/01/01 36 // 37 //-------------------------------------------- 38 // 39 40 //....oooOO0OOooo........oooOO0OOooo........oo 41 //....oooOO0OOooo........oooOO0OOooo........oo 42 43 #include "HistoManager.hh" 44 45 #include "EmAcceptance.hh" 46 #include "Histo.hh" 47 48 #include "G4Electron.hh" 49 #include "G4EmProcessSubType.hh" 50 #include "G4Gamma.hh" 51 #include "G4GammaGeneralProcess.hh" 52 #include "G4MaterialCutsCouple.hh" 53 #include "G4Positron.hh" 54 #include "G4SystemOfUnits.hh" 55 #include "G4UnitsTable.hh" 56 #include "G4VEmProcess.hh" 57 #include "G4VEnergyLossProcess.hh" 58 #include "G4VProcess.hh" 59 60 //....oooOO0OOooo........oooOO0OOooo........oo 61 62 HistoManager* HistoManager::fManager = nullptr 63 64 //....oooOO0OOooo........oooOO0OOooo........oo 65 66 HistoManager* HistoManager::GetPointer() 67 { 68 if (nullptr == fManager) { 69 fManager = new HistoManager(); 70 } 71 return fManager; 72 } 73 74 //....oooOO0OOooo........oooOO0OOooo........oo 75 76 HistoManager::HistoManager() 77 : fGamma(G4Gamma::Gamma()), 78 fElectron(G4Electron::Electron()), 79 fPositron(G4Positron::Positron()), 80 fHisto(new Histo()) 81 { 82 fVerbose = 1; 83 fEvt1 = -1; 84 fEvt2 = -1; 85 fNmax = 3; 86 fMaxEnergy = 50.0 * MeV; 87 fBeamEnergy = 1. * GeV; 88 fMaxEnergyAbs = 10. * MeV; 89 fBinsE = 100; 90 fBinsEA = 40; 91 fBinsED = 100; 92 fNHisto = 13; 93 94 BookHisto(); 95 } 96 97 //....oooOO0OOooo........oooOO0OOooo........oo 98 99 HistoManager::~HistoManager() 100 { 101 delete fHisto; 102 } 103 104 //....oooOO0OOooo........oooOO0OOooo........oo 105 106 void HistoManager::BookHisto() 107 { 108 fHisto->Add1D("10", "Evis/E0 in central crys 109 fHisto->Add1D("11", "Evis/E0 in 3x3", fBinsE 110 fHisto->Add1D("12", "Evis/E0 in 5x5", fBinsE 111 fHisto->Add1D("13", "Energy (MeV) of delta-e 112 fHisto->Add1D("14", "Energy (MeV) of gammas" 113 fHisto->Add1D("15", "Energy (MeV) in abs1", 114 fHisto->Add1D("16", "Energy (MeV) in abs2", 115 fHisto->Add1D("17", "Energy (MeV) in abs3", 116 fHisto->Add1D("18", "Energy (MeV) in abs4", 117 fHisto->Add1D("19", "Number of vertex hits", 118 fHisto->Add1D("20", "E1/E9 Ratio", fBinsED, 119 fHisto->Add1D("21", "E1/E25 Ratio", fBinsED, 120 fHisto->Add1D("22", "E9/E25 Ratio", fBinsED, 121 } 122 123 //....oooOO0OOooo........oooOO0OOooo........oo 124 125 void HistoManager::BeginOfRun() 126 { 127 // initilise scoring 128 fEvt = 0; 129 fElec = 0; 130 fPosit = 0; 131 fGam = 0; 132 fStep = 0; 133 fLowe = 0; 134 135 for (G4int i = 0; i < 6; i++) { 136 fStat[i] = 0; 137 fEdep[i] = 0.0; 138 fErms[i] = 0.0; 139 if (i < 3) { 140 fEdeptr[i] = 0.0; 141 fErmstr[i] = 0.0; 142 } 143 } 144 145 // initialise counters 146 fBrem.resize(93, 0.0); 147 fPhot.resize(93, 0.0); 148 fComp.resize(93, 0.0); 149 fConv.resize(93, 0.0); 150 151 // initialise acceptance - by default is not 152 for (G4int i = 0; i < fNmax; i++) { 153 fEdeptrue[i] = 1.0; 154 fRmstrue[i] = 1.0; 155 fLimittrue[i] = 10.; 156 } 157 158 if (fHisto->IsActive()) { 159 for (G4int i = 0; i < fNHisto; ++i) { 160 fHisto->Activate(i, true); 161 } 162 fHisto->Book(); 163 164 if (fVerbose > 0) { 165 G4cout << "HistoManager: Histograms are 166 } 167 } 168 } 169 170 //....oooOO0OOooo........oooOO0OOooo........oo 171 172 void HistoManager::EndOfRun(G4int runID) 173 { 174 G4cout << "HistoManager: End of run actions 175 G4String nam[6] = {"1x1", "3x3", "5x5", "E1/ 176 177 // average 178 179 G4cout << "================================= 180 G4double x = (G4double)fEvt; 181 if (fEvt > 0) x = 1.0 / x; 182 G4int j; 183 for (j = 0; j < fNmax; j++) { 184 // total mean 185 fEdep[j] *= x; 186 G4double y = fErms[j] * x - fEdep[j] * fEd 187 if (y < 0.0) y = 0.0; 188 fErms[j] = std::sqrt(y); 189 190 // trancated mean 191 G4double xx = G4double(fStat[j]); 192 if (xx > 0.0) xx = 1.0 / xx; 193 fEdeptr[j] *= xx; 194 y = fErmstr[j] * xx - fEdeptr[j] * fEdeptr 195 if (y < 0.0) y = 0.0; 196 fErmstr[j] = std::sqrt(y); 197 } 198 G4double xe = x * (G4double)fElec; 199 G4double xg = x * (G4double)fGam; 200 G4double xp = x * (G4double)fPosit; 201 G4double xs = x * fStep; 202 203 G4double f = 100. * std::sqrt(fBeamEnergy / 204 205 G4cout << "Number of events " << 206 G4cout << std::setprecision(4) << "Average n 207 G4cout << std::setprecision(4) << "Average n 208 G4cout << std::setprecision(4) << "Average n 209 G4cout << std::setprecision(4) << "Average n 210 211 for (j = 0; j < 3; ++j) { 212 G4double ex = fEdeptr[j]; 213 G4double sx = fErmstr[j]; 214 G4double xx = G4double(fStat[j]); 215 if (xx > 0.0) xx = 1.0 / xx; 216 G4double r = sx * std::sqrt(xx); 217 G4cout << std::setprecision(4) << "Edep " 218 << r; 219 if (ex > 0.1) G4cout << " res= " << f * 220 G4cout << G4endl; 221 } 222 if (fLimittrue[0] < 10. || fLimittrue[1] < 1 223 G4cout << "=========== Mean values withou 224 for (j = 0; j < fNmax; j++) { 225 G4double ex = fEdep[j]; 226 G4double sx = fErms[j]; 227 G4double rx = sx * std::sqrt(x); 228 G4cout << std::setprecision(4) << "Edep 229 << rx; 230 if (ex > 0.0) G4cout << " res= " << f 231 G4cout << G4endl; 232 } 233 } 234 G4cout << "=========== Ratios without tranc 235 for (j = 3; j < 6; ++j) { 236 G4double e = fEdep[j]; 237 G4double xx = G4double(fStat[j]); 238 if (xx > 0.0) xx = 1.0 / xx; 239 e *= xx; 240 G4double y = fErms[j] * xx - e * e; 241 G4double r = 0.0; 242 if (y > 0.0) r = std::sqrt(y * xx); 243 G4cout << " " << nam[j] << " = 244 G4cout << G4endl; 245 } 246 G4cout << std::setprecision(4) << "Beam Ener 247 << G4endl; 248 if (fLowe > 0) G4cout << "Number of events E 249 G4cout << "================================= 250 G4cout << G4endl; 251 252 // normalise histograms 253 if (fHisto->IsActive()) { 254 for (G4int i = 0; i < fNHisto; ++i) { 255 fHisto->ScaleH1(i, x); 256 } 257 fHisto->Save(); 258 } 259 if (0 < runID) { 260 return; 261 } 262 263 // Acceptance only for the first run 264 EmAcceptance acc; 265 G4bool isStarted = false; 266 for (j = 0; j < fNmax; j++) { 267 G4double ltrue = fLimittrue[j]; 268 if (ltrue < DBL_MAX) { 269 if (!isStarted) { 270 acc.BeginOfAcceptance("Crystal Calorim 271 isStarted = true; 272 } 273 G4double etrue = fEdeptrue[j]; 274 G4double rtrue = fRmstrue[j]; 275 acc.EmAcceptanceGauss("Edep" + nam[j], f 276 acc.EmAcceptanceGauss("Erms" + nam[j], f 277 } 278 } 279 if (isStarted) acc.EndOfAcceptance(); 280 281 // atom frequency 282 G4cout << " Z bremsstrahlung photoeffect 283 for (j = 1; j < 93; ++j) { 284 G4int n1 = G4int(fBrem[j] * x); 285 G4int n2 = G4int(fPhot[j] * x); 286 G4int n3 = G4int(fComp[j] * x); 287 G4int n4 = G4int(fConv[j] * x); 288 if (n1 + n2 + n3 + n4 > 0) { 289 G4cout << std::setw(4) << j << std::setw 290 << n3 << std::setw(12) << n4 << G 291 } 292 } 293 } 294 295 //....oooOO0OOooo........oooOO0OOooo........oo 296 297 void HistoManager::BeginOfEvent() 298 { 299 ++fEvt; 300 301 fEabs1 = 0.0; 302 fEabs2 = 0.0; 303 fEabs3 = 0.0; 304 fEabs4 = 0.0; 305 fEvertex.clear(); 306 fNvertex.clear(); 307 for (G4int i = 0; i < 25; i++) { 308 fE[i] = 0.0; 309 } 310 } 311 312 //....oooOO0OOooo........oooOO0OOooo........oo 313 314 void HistoManager::EndOfEvent() 315 { 316 G4double e9 = 0.0; 317 G4double e25 = 0.0; 318 for (G4int i = 0; i < 25; i++) { 319 fE[i] /= fBeamEnergy; 320 e25 += fE[i]; 321 if ((6 <= i && 8 >= i) || (11 <= i && 13 > 322 } 323 324 if (1 < fVerbose && e25 < 0.8) { 325 ++fLowe; 326 G4cout << "### in the event# " << fEvt << 327 } 328 329 // compute ratios 330 G4double e0 = fE[12]; 331 G4double e19 = 0.0; 332 G4double e125 = 0.0; 333 G4double e925 = 0.0; 334 if (e9 > 0.0) { 335 e19 = e0 / e9; 336 e125 = e0 / e25; 337 e925 = e9 / e25; 338 fEdep[3] += e19; 339 fErms[3] += e19 * e19; 340 fEdep[4] += e125; 341 fErms[4] += e125 * e125; 342 fEdep[5] += e925; 343 fErms[5] += e925 * e925; 344 fStat[3] += 1; 345 fStat[4] += 1; 346 fStat[5] += 1; 347 } 348 349 // Fill histo 350 fHisto->Fill(0, e0, 1.0); 351 fHisto->Fill(1, e9, 1.0); 352 fHisto->Fill(2, e25, 1.0); 353 fHisto->Fill(5, fEabs1, 1.0); 354 fHisto->Fill(6, fEabs2, 1.0); 355 fHisto->Fill(7, fEabs3, 1.0); 356 fHisto->Fill(8, fEabs4, 1.0); 357 fHisto->Fill(9, G4double(fNvertex.size()), 1 358 fHisto->Fill(10, e19, 1.0); 359 fHisto->Fill(11, e125, 1.0); 360 fHisto->Fill(12, e925, 1.0); 361 362 // compute sums 363 fEdep[0] += e0; 364 fErms[0] += e0 * e0; 365 fEdep[1] += e9; 366 fErms[1] += e9 * e9; 367 fEdep[2] += e25; 368 fErms[2] += e25 * e25; 369 370 // trancated mean 371 if (std::abs(e0 - fEdeptrue[0]) < fRmstrue[0 372 fStat[0] += 1; 373 fEdeptr[0] += e0; 374 fErmstr[0] += e0 * e0; 375 } 376 if (std::abs(e9 - fEdeptrue[1]) < fRmstrue[1 377 fStat[1] += 1; 378 fEdeptr[1] += e9; 379 fErmstr[1] += e9 * e9; 380 } 381 if (std::abs(e25 - fEdeptrue[2]) < fRmstrue[ 382 fStat[2] += 1; 383 fEdeptr[2] += e25; 384 fErmstr[2] += e25 * e25; 385 } 386 } 387 388 //....oooOO0OOooo........oooOO0OOooo........oo 389 390 void HistoManager::ScoreNewTrack(const G4Track 391 { 392 // Save primary parameters 393 ResetTrackLength(); 394 const G4ParticleDefinition* particle = aTrac 395 const G4DynamicParticle* dynParticle = aTrac 396 397 G4int pid = aTrack->GetParentID(); 398 G4double kinE = dynParticle->GetKineticEnerg 399 G4ThreeVector pos = aTrack->GetVertexPositio 400 401 // primary 402 if (0 == pid) { 403 G4double mass = 0.0; 404 if (particle) { 405 mass = particle->GetPDGMass(); 406 } 407 408 G4ThreeVector dir = dynParticle->GetMoment 409 if (1 < fVerbose) { 410 G4cout << "TrackingAction: Primary kinE( 411 << "; pos= " << pos << "; dir= " 412 } 413 414 // secondary 415 } 416 else { 417 const G4VProcess* proc = aTrack->GetCreato 418 G4int type = proc->GetProcessSubType(); 419 420 if (type == fBremsstrahlung) { 421 auto elm = static_cast<const G4VEnergyLo 422 if (nullptr != elm) { 423 G4int Z = elm->GetZasInt(); 424 if (Z > 0 && Z < 93) { 425 fBrem[Z] += 1.0; 426 } 427 } 428 } 429 else if (type == fPhotoElectricEffect) { 430 auto elm = static_cast<const G4VEmProces 431 if (nullptr != elm) { 432 G4int Z = elm->GetZasInt(); 433 if (Z > 0 && Z < 93) { 434 fPhot[Z] += 1.0; 435 } 436 } 437 } 438 else if (type == fGammaConversion) { 439 auto elm = static_cast<const G4VEmProces 440 if (nullptr != elm) { 441 G4int Z = elm->GetZasInt(); 442 if (Z > 0 && Z < 93) { 443 fConv[Z] += 1.0; 444 } 445 } 446 } 447 else if (type == fComptonScattering) { 448 auto elm = static_cast<const G4VEmProces 449 if (nullptr != elm) { 450 G4int Z = elm->GetZasInt(); 451 if (Z > 0 && Z < 93) { 452 fComp[Z] += 1.0; 453 } 454 } 455 } 456 457 // delta-electron 458 if (particle == fElectron) { 459 if (1 < fVerbose) { 460 G4cout << "TrackingAction: Secondary e 461 } 462 AddDeltaElectron(dynParticle); 463 } 464 else if (particle == fPositron) { 465 if (1 < fVerbose) { 466 G4cout << "TrackingAction: Secondary p 467 } 468 AddPositron(dynParticle); 469 } 470 else if (particle == fGamma) { 471 if (1 < fVerbose) { 472 G4cout << "TrackingAction: Secondary g 473 } 474 AddPhoton(dynParticle); 475 } 476 } 477 } 478 479 //....oooOO0OOooo........oooOO0OOooo........oo 480 481 void HistoManager::AddEnergy(G4double edep, G4 482 { 483 if (1 < fVerbose) { 484 G4cout << "HistoManager::AddEnergy: e(keV) 485 << "; copyNo= " << copyNo << G4endl 486 } 487 if (0 == volIndex) { 488 fE[copyNo] += edep; 489 } 490 else if (1 == volIndex) { 491 fEabs1 += edep; 492 } 493 else if (2 == volIndex) { 494 fEabs2 += edep; 495 } 496 else if (3 == volIndex) { 497 fEabs3 += edep; 498 } 499 else if (4 == volIndex) { 500 fEabs4 += edep; 501 } 502 else if (5 == volIndex) { 503 G4int n = fNvertex.size(); 504 G4bool newPad = true; 505 if (n > 0) { 506 for (G4int i = 0; i < n; i++) { 507 if (copyNo == fNvertex[i]) { 508 newPad = false; 509 fEvertex[i] += edep; 510 break; 511 } 512 } 513 } 514 if (newPad) { 515 fNvertex.push_back(copyNo); 516 fEvertex.push_back(edep); 517 } 518 } 519 } 520 521 //....oooOO0OOooo........oooOO0OOooo........oo 522 523 void HistoManager::AddDeltaElectron(const G4Dy 524 { 525 G4double e = elec->GetKineticEnergy() / MeV; 526 if (e > 0.0) { 527 ++fElec; 528 fHisto->Fill(3, e, 1.0); 529 } 530 } 531 532 //....oooOO0OOooo........oooOO0OOooo........oo 533 534 void HistoManager::AddPhoton(const G4DynamicPa 535 { 536 G4double e = ph->GetKineticEnergy() / MeV; 537 if (e > 0.0) { 538 ++fGam; 539 fHisto->Fill(4, e, 1.0); 540 } 541 } 542 543 //....oooOO0OOooo........oooOO0OOooo........oo 544 545 void HistoManager::SetEdepAndRMS(G4int i, cons 546 { 547 if (i < fNmax && i >= 0) { 548 if (val[0] > 0.0) fEdeptrue[i] = val[0]; 549 if (val[1] > 0.0) fRmstrue[i] = val[1]; 550 if (val[2] > 0.0) fLimittrue[i] = val[2]; 551 } 552 } 553 554 //....oooOO0OOooo........oooOO0OOooo........oo 555