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 polarisation/Pol01/src/RunAction.cc 27 /// \brief Implementation of the RunAction cla 28 // 29 // 30 //....oooOO0OOooo........oooOO0OOooo........oo 31 //....oooOO0OOooo........oooOO0OOooo........oo 32 33 #include "RunAction.hh" 34 35 #include "DetectorConstruction.hh" 36 #include "PrimaryGeneratorAction.hh" 37 #include "ProcessesCount.hh" 38 39 #include "G4AccumulableManager.hh" 40 #include "G4Electron.hh" 41 #include "G4EmCalculator.hh" 42 #include "G4Gamma.hh" 43 #include "G4ParticleDefinition.hh" 44 #include "G4PhysicalConstants.hh" 45 #include "G4Positron.hh" 46 #include "G4Run.hh" 47 #include "G4RunManager.hh" 48 #include "G4SystemOfUnits.hh" 49 #include "G4UnitsTable.hh" 50 51 #include <iomanip> 52 53 //....oooOO0OOooo........oooOO0OOooo........oo 54 55 RunAction::RunAction(DetectorConstruction* det 56 : fDetector(det), fPrimary(prim), fAnalysisM 57 { 58 fGamma = G4Gamma::Gamma(); 59 fElectron = G4Electron::Electron(); 60 fPositron = G4Positron::Positron(); 61 62 auto accumulableManager = G4AccumulableManag 63 auto fPhotonStats = new ParticleStatistics(" 64 auto fElectronStats = new ParticleStatistics 65 auto fPositronStats = new ParticleStatistics 66 auto fProcCounter = new ProcessesCount("Proc 67 68 accumulableManager->Register(fPhotonStats); 69 accumulableManager->Register(fElectronStats) 70 accumulableManager->Register(fPositronStats) 71 72 accumulableManager->Register(fProcCounter); 73 74 BookHisto(); 75 } 76 77 //....oooOO0OOooo........oooOO0OOooo........oo 78 79 RunAction::~RunAction() {} 80 81 //....oooOO0OOooo........oooOO0OOooo........oo 82 83 void RunAction::BeginOfRunAction(const G4Run* 84 { 85 G4cout << "### Run " << aRun->GetRunID() << 86 87 auto accumulableManager = G4AccumulableManag 88 accumulableManager->Reset(); 89 90 // save Rndm status 91 // G4RunManager::GetRunManager()->SetRandom 92 // CLHEP::HepRandom::showEngineStatus(); 93 94 fTotalEventCount = 0; 95 96 // Open file histogram file 97 fAnalysisManager->OpenFile(); 98 } 99 100 //....oooOO0OOooo........oooOO0OOooo........oo 101 102 void RunAction::FillData(const G4ParticleDefin 103 G4double costheta, G4 104 { 105 auto accManager = G4AccumulableManager::Inst 106 G4int id = -1; 107 if (particle == fGamma) { 108 dynamic_cast<ParticleStatistics*>(accManag 109 ->FillData(kinEnergy, costheta, longitud 110 if (fAnalysisManager) { 111 id = 1; 112 } 113 } 114 else if (particle == fElectron) { 115 dynamic_cast<ParticleStatistics*>(accManag 116 ->FillData(kinEnergy, costheta, longitud 117 if (fAnalysisManager) { 118 id = 5; 119 } 120 } 121 else if (particle == fPositron) { 122 dynamic_cast<ParticleStatistics*>(accManag 123 ->FillData(kinEnergy, costheta, longitud 124 if (fAnalysisManager) { 125 id = 9; 126 } 127 } 128 if (id > 0) { 129 fAnalysisManager->FillH1(id, kinEnergy, 1. 130 fAnalysisManager->FillH1(id + 1, costheta, 131 fAnalysisManager->FillH1(id + 2, phi, 1.0) 132 fAnalysisManager->FillH1(id + 3, longitudi 133 } 134 } 135 136 //....oooOO0OOooo........oooOO0OOooo........oo 137 138 void RunAction::BookHisto() 139 { 140 // Always creating analysis manager 141 fAnalysisManager = G4AnalysisManager::Instan 142 fAnalysisManager->SetDefaultFileType("root") 143 fAnalysisManager->SetActivation(true); 144 fAnalysisManager->SetVerboseLevel(1); 145 146 // Open file histogram file 147 fAnalysisManager->SetFileName("pol01"); 148 149 fAnalysisManager->SetFirstHistoId(1); 150 151 // Creating an 1-dimensional histograms in t 152 const G4String id[] = {"h1", "h2", "h3", "h4 153 const G4String title[] = { 154 "Gamma Energy distribution", // 1 155 "Gamma Cos(Theta) distribution", // 2 156 "Gamma Phi angular distribution", // 3 157 "Gamma longitudinal Polarization", // 4 158 "Electron Energy distribution", // 5 159 "Electron Cos(Theta) distribution", // 6 160 "Electron Phi angular distribution", // 7 161 "Electron longitudinal Polarization", // 162 "Positron Energy distribution", // 9 163 "Positron Cos(Theta) distribution", // 10 164 "Positron Phi angular distribution", // 1 165 "Positron longitudinal Polarization" // 1 166 }; 167 G4double vmin, vmax; 168 G4int nbins = 120; 169 for (int i = 0; i < 12; ++i) { 170 G4int j = i - i / 4 * 4; 171 if (0 == j) { 172 vmin = 0.; 173 vmax = 12. * MeV; 174 } 175 else if (1 == j) { 176 vmin = -1.; 177 vmax = 1.; 178 } 179 else if (2 == j) { 180 vmin = 0.; 181 vmax = pi; 182 } 183 else { 184 vmin = -1.5; 185 vmax = 1.5; 186 } 187 G4int ih = fAnalysisManager->CreateH1(id[i 188 fAnalysisManager->SetH1Activation(ih, fals 189 } 190 } 191 192 //....oooOO0OOooo........oooOO0OOooo........oo 193 194 void RunAction::SaveHisto(G4int nevents) 195 { 196 if (fAnalysisManager) { 197 if (IsMaster()) { 198 G4double norm = 1.0 / G4double(nevents); 199 for (int i = 0; i < 12; ++i) { 200 fAnalysisManager->ScaleH1(i, norm); 201 } 202 } 203 fAnalysisManager->Write(); 204 fAnalysisManager->CloseFile(); 205 } 206 } 207 208 //....oooOO0OOooo........oooOO0OOooo........oo 209 210 void RunAction::CountProcesses(G4String& procN 211 { 212 auto accManager = G4AccumulableManager::Inst 213 dynamic_cast<ProcessesCount*>(accManager->Ge 214 } 215 216 //....oooOO0OOooo........oooOO0OOooo........oo 217 218 void RunAction::EndOfRunAction(const G4Run* aR 219 { 220 // Total number of events in run (all threa 221 G4int NbOfEvents = aRun->GetNumberOfEventToB 222 // G4int NbOfEvents = aRun->GetNumberOfEven 223 224 if (NbOfEvents == 0) return; 225 226 G4int prec = G4cout.precision(5); 227 228 G4Material* material = fDetector->GetMateria 229 G4double density = material->GetDensity(); 230 231 if (fPrimary != nullptr) { 232 G4ParticleDefinition* particle = fPrimary- 233 G4String Particle = particle->GetParticleN 234 G4double energy = fPrimary->GetParticleGun 235 G4cout << "\n The run consists of " << fTo 236 << G4BestUnit(energy, "Energy") << 237 << G4BestUnit(fDetector->GetBoxSize 238 << " (density: " << G4BestUnit(dens 239 } 240 // cross check from G4EmCalculator 241 // G4cout << "\n Verification from G4EmCal 242 // G4EmCalculator emCal; 243 244 auto accManager = G4AccumulableManager::Inst 245 accManager->Merge(); 246 247 if (IsMaster()) { 248 // frequency of processes 249 G4cout << "\n Process calls frequency ---> 250 dynamic_cast<ProcessesCount*>(accManager-> 251 G4cout << " Gamma: \n"; 252 dynamic_cast<ParticleStatistics*>(accManag 253 ->PrintResults(NbOfEvents); 254 G4cout << " Electron: \n"; 255 dynamic_cast<ParticleStatistics*>(accManag 256 ->PrintResults(NbOfEvents); 257 G4cout << " Positron: \n"; 258 dynamic_cast<ParticleStatistics*>(accManag 259 ->PrintResults(NbOfEvents); 260 G4cout << G4endl; 261 } 262 263 // restore default format 264 G4cout.precision(prec); 265 266 // write out histograms 267 SaveHisto(NbOfEvents); 268 269 if (IsMaster()) { 270 // show Rndm status 271 CLHEP::HepRandom::showEngineStatus(); 272 } 273 } 274 275 //....oooOO0OOooo........oooOO0OOooo........oo 276 277 void RunAction::EventFinished() 278 { 279 auto accManager = G4AccumulableManager::Inst 280 281 ++fTotalEventCount; 282 283 dynamic_cast<ParticleStatistics*>(accManager 284 dynamic_cast<ParticleStatistics*>(accManager 285 dynamic_cast<ParticleStatistics*>(accManager 286 } 287 288 //....oooOO0OOooo........oooOO0OOooo........oo 289 290 RunAction::ParticleStatistics::ParticleStatist 291 : G4VAccumulable(name), 292 fCurrentNumber(0), 293 fTotalNumber(0), 294 fTotalNumber2(0), 295 fSumEnergy(0), 296 fSumEnergy2(0), 297 fSumPolarization(0), 298 fSumPolarization2(0), 299 fSumCosTheta(0), 300 fSumCosTheta2(0) 301 {} 302 303 //....oooOO0OOooo........oooOO0OOooo........oo 304 305 RunAction::ParticleStatistics::~ParticleStatis 306 307 //....oooOO0OOooo........oooOO0OOooo........oo 308 309 void RunAction::ParticleStatistics::EventFinis 310 { 311 fTotalNumber += fCurrentNumber; 312 fTotalNumber2 += fCurrentNumber * fCurrentNu 313 fCurrentNumber = 0; 314 } 315 //....oooOO0OOooo........oooOO0OOooo........oo 316 317 void RunAction::ParticleStatistics::FillData(G 318 G 319 { 320 ++fCurrentNumber; 321 fSumEnergy += kinEnergy; 322 fSumEnergy2 += kinEnergy * kinEnergy; 323 fSumPolarization += longitudinalPolarization 324 fSumPolarization2 += longitudinalPolarizatio 325 fSumCosTheta += costheta; 326 fSumCosTheta2 += costheta * costheta; 327 } 328 329 //....oooOO0OOooo........oooOO0OOooo........oo 330 331 void RunAction::ParticleStatistics::PrintResul 332 { 333 G4cout << "Mean Number per Event :" << G4dou 334 << "\n"; 335 if (fTotalNumber == 0) fTotalNumber = 1; 336 G4double energyMean = fSumEnergy / fTotalNum 337 G4double energyRms = std::sqrt(fSumEnergy2 / 338 G4cout << "Mean Energy :" << G4BestUnit(ener 339 << G4BestUnit(energyRms, "Energy") << 340 G4double polarizationMean = fSumPolarization 341 G4double polarizationRms = 342 std::sqrt(fSumPolarization2 / fTotalNumber 343 G4cout << "Mean Polarization :" << polarizat 344 } 345 346 //....oooOO0OOooo........oooOO0OOooo........oo 347 348 void RunAction::ParticleStatistics::Reset() 349 { 350 fCurrentNumber = 0; 351 fTotalNumber = fTotalNumber2 = 0; 352 fSumEnergy = fSumEnergy2 = 0; 353 fSumPolarization = fSumPolarization2 = 0; 354 fSumCosTheta = fSumCosTheta2 = 0; 355 } 356 357 //....oooOO0OOooo........oooOO0OOooo........oo 358 359 void RunAction::ParticleStatistics::Merge(cons 360 { 361 auto rstat = dynamic_cast<const RunAction::P 362 363 fCurrentNumber += rstat.fCurrentNumber; 364 fTotalNumber += rstat.fTotalNumber; 365 fTotalNumber2 += rstat.fTotalNumber2; 366 fSumEnergy += rstat.fSumEnergy; 367 fSumEnergy2 += rstat.fSumEnergy2; 368 fSumPolarization += rstat.fSumPolarization; 369 fSumPolarization2 += rstat.fSumPolarization2 370 fSumCosTheta += rstat.fSumCosTheta; 371 fSumCosTheta2 += rstat.fSumCosTheta2; 372 } 373 374 //....oooOO0OOooo........oooOO0OOooo........oo 375