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/TestEm15/src/RunActi 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 "HistoManager.hh" 37 #include "PrimaryGeneratorAction.hh" 38 39 #include "G4EmCalculator.hh" 40 #include "G4Run.hh" 41 #include "G4RunManager.hh" 42 #include "G4SystemOfUnits.hh" 43 #include "G4UnitsTable.hh" 44 #include "Randomize.hh" 45 46 #include <iomanip> 47 48 //....oooOO0OOooo........oooOO0OOooo........oo 49 50 RunAction::RunAction(DetectorConstruction* det 51 : G4UserRunAction(), fDetector(det), fPrimar 52 { 53 fHistoManager = new HistoManager(); 54 } 55 56 //....oooOO0OOooo........oooOO0OOooo........oo 57 58 RunAction::~RunAction() 59 { 60 delete fHistoManager; 61 } 62 63 //....oooOO0OOooo........oooOO0OOooo........oo 64 65 void RunAction::BeginOfRunAction(const G4Run*) 66 { 67 // save Rndm status 68 ////G4RunManager::GetRunManager()->SetRandom 69 CLHEP::HepRandom::showEngineStatus(); 70 71 fProcCounter = new ProcessesCount; 72 fTotalCount = 0; 73 74 fTruePL = fTruePL2 = fGeomPL = fGeomPL2 = 0. 75 fLDispl = fLDispl2 = fPsiSpa = fPsiSpa2 = 0. 76 fTetPrj = fTetPrj2 = 0.; 77 fPhiCor = fPhiCor2 = 0.; 78 79 // histograms 80 // 81 G4AnalysisManager* analysisManager = G4Analy 82 if (analysisManager->IsActive()) { 83 analysisManager->OpenFile(); 84 } 85 } 86 87 //....oooOO0OOooo........oooOO0OOooo........oo 88 89 void RunAction::CountProcesses(G4String procNa 90 { 91 // does the process already encounted ? 92 size_t nbProc = fProcCounter->size(); 93 size_t i = 0; 94 while ((i < nbProc) && ((*fProcCounter)[i]-> 95 i++; 96 if (i == nbProc) fProcCounter->push_back(new 97 98 (*fProcCounter)[i]->Count(); 99 } 100 101 //....oooOO0OOooo........oooOO0OOooo........oo 102 103 void RunAction::EndOfRunAction(const G4Run* aR 104 { 105 G4int NbOfEvents = aRun->GetNumberOfEvent(); 106 if (NbOfEvents == 0) return; 107 108 G4int prec = G4cout.precision(5); 109 110 G4Material* material = fDetector->GetMateria 111 G4double density = material->GetDensity(); 112 113 G4ParticleDefinition* particle = fPrimary->G 114 G4String Particle = particle->GetParticleNam 115 G4double energy = fPrimary->GetParticleGun() 116 G4cout << "\n The run consists of " << NbOfE 117 << G4BestUnit(energy, "Energy") << " 118 << G4BestUnit(fDetector->GetBoxSize() 119 << " (density: " << G4BestUnit(densit 120 121 // frequency of processes 122 G4cout << "\n Process calls frequency --->"; 123 for (size_t i = 0; i < fProcCounter->size(); 124 G4String procName = (*fProcCounter)[i]->Ge 125 G4int count = (*fProcCounter)[i]->GetCount 126 G4cout << "\t" << procName << " = " << cou 127 } 128 129 if (fTotalCount > 0) { 130 // compute path length and related quantit 131 // 132 G4double MeanTPL = fTruePL / fTotalCount; 133 G4double MeanTPL2 = fTruePL2 / fTotalCount 134 G4double rmsTPL = std::sqrt(std::fabs(Mean 135 136 G4double MeanGPL = fGeomPL / fTotalCount; 137 G4double MeanGPL2 = fGeomPL2 / fTotalCount 138 G4double rmsGPL = std::sqrt(std::fabs(Mean 139 140 G4double MeanLaD = fLDispl / fTotalCount; 141 G4double MeanLaD2 = fLDispl2 / fTotalCount 142 G4double rmsLaD = std::sqrt(std::fabs(Mean 143 144 G4double MeanPsi = fPsiSpa / (fTotalCount) 145 G4double MeanPsi2 = fPsiSpa2 / (fTotalCoun 146 G4double rmsPsi = std::sqrt(std::fabs(Mean 147 148 G4double MeanTeta = fTetPrj / (2 * fTotalC 149 G4double MeanTeta2 = fTetPrj2 / (2 * fTota 150 G4double rmsTeta = std::sqrt(std::fabs(Mea 151 152 G4double MeanCorrel = fPhiCor / (fTotalCou 153 G4double MeanCorrel2 = fPhiCor2 / (fTotalC 154 G4double rmsCorrel = std::sqrt(std::fabs(M 155 156 G4cout << "\n\n truePathLength :\t" << G4B 157 << G4BestUnit(rmsTPL, "Length") << 158 << G4BestUnit(MeanGPL, "Length") << 159 << "\n lateralDisplac :\t" << G4Bes 160 << G4BestUnit(rmsLaD, "Length") << 161 << " +- " << rmsPsi / mrad << " mra 162 << " (" << MeanPsi / deg << " deg" 163 << " +- " << rmsPsi / deg << " deg) 164 165 G4cout << "\n Theta_plane :\t" << rmsTe 166 << " (" << rmsTeta / deg << " deg) 167 << "\n phi correlation:\t" << MeanC 168 << " (std::cos(phi_pos - phi_dir)) 169 170 // cross check from G4EmCalculator 171 // 172 G4cout << "\n Verification from G4EmCalcul 173 174 G4EmCalculator emCal; 175 176 // get transport mean free path (for multi 177 G4double MSmfp = emCal.GetMeanFreePath(ene 178 179 // get range from restricted dedx 180 G4double range = emCal.GetRangeFromRestric 181 182 // effective facRange 183 G4double efFacrange = MeanTPL / std::max(M 184 if (MeanTPL / range >= 0.99) efFacrange = 185 186 G4cout << "\n transport mean free path :\t 187 << "\n range from restrict dE/dx:\t 188 << "\n ---> effective facRange :\t 189 190 G4cout << "\n compute theta0 from Highland 191 << " (" << ComputeMscHighland(Mean 192 } 193 else 194 G4cout << G4endl; 195 196 // restore default format 197 G4cout.precision(prec); 198 199 // delete and remove all contents in fProcCo 200 while (fProcCounter->size() > 0) { 201 OneProcessCount* aProcCount = fProcCounter 202 fProcCounter->pop_back(); 203 delete aProcCount; 204 } 205 delete fProcCounter; 206 207 // save histograms 208 G4AnalysisManager* analysisManager = G4Analy 209 if (analysisManager->IsActive()) { 210 analysisManager->Write(); 211 analysisManager->CloseFile(); 212 } 213 214 // show Rndm status 215 CLHEP::HepRandom::showEngineStatus(); 216 } 217 218 //....oooOO0OOooo........oooOO0OOooo........oo 219 220 G4double RunAction::ComputeMscHighland(G4doubl 221 { 222 // compute the width of the Gaussian central 223 // projected angular distribution. 224 // Eur. Phys. Jour. C15 (2000) page 166, for 225 226 G4double t = pathLength / (fDetector->GetMat 227 if (t < DBL_MIN) return 0.; 228 229 G4ParticleGun* particle = fPrimary->GetParti 230 G4double T = particle->GetParticleEnergy(); 231 G4double M = particle->GetParticleDefinition 232 G4double z = std::abs(particle->GetParticleD 233 234 G4double bpc = T * (T + 2 * M) / (T + M); 235 G4double teta0 = 13.6 * MeV * z * std::sqrt( 236 return teta0; 237 } 238 239 //....oooOO0OOooo........oooOO0OOooo........oo 240