Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 /// \file electromagnetic/TestEm15/src/RunActi 26 /// \file electromagnetic/TestEm15/src/RunAction.cc 27 /// \brief Implementation of the RunAction cla 27 /// \brief Implementation of the RunAction class 28 // 28 // 29 // << 29 // 30 //....oooOO0OOooo........oooOO0OOooo........oo 30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 31 //....oooOO0OOooo........oooOO0OOooo........oo 31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 32 32 33 #include "RunAction.hh" 33 #include "RunAction.hh" 34 34 35 #include "DetectorConstruction.hh" 35 #include "DetectorConstruction.hh" 36 #include "HistoManager.hh" << 37 #include "PrimaryGeneratorAction.hh" 36 #include "PrimaryGeneratorAction.hh" >> 37 #include "HistoManager.hh" 38 38 39 #include "G4EmCalculator.hh" << 40 #include "G4Run.hh" 39 #include "G4Run.hh" 41 #include "G4RunManager.hh" 40 #include "G4RunManager.hh" 42 #include "G4SystemOfUnits.hh" << 43 #include "G4UnitsTable.hh" 41 #include "G4UnitsTable.hh" 44 #include "Randomize.hh" << 42 #include "G4EmCalculator.hh" 45 43 >> 44 #include "Randomize.hh" >> 45 #include "G4SystemOfUnits.hh" 46 #include <iomanip> 46 #include <iomanip> 47 47 48 //....oooOO0OOooo........oooOO0OOooo........oo 48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 49 49 50 RunAction::RunAction(DetectorConstruction* det 50 RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* prim) 51 : G4UserRunAction(), fDetector(det), fPrimar << 51 : G4UserRunAction(),fDetector(det), fPrimary(prim), fProcCounter(0), >> 52 fHistoManager(0) 52 { 53 { 53 fHistoManager = new HistoManager(); << 54 fHistoManager = new HistoManager(); 54 } 55 } 55 56 56 //....oooOO0OOooo........oooOO0OOooo........oo 57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 57 58 58 RunAction::~RunAction() 59 RunAction::~RunAction() 59 { 60 { 60 delete fHistoManager; << 61 delete fHistoManager; 61 } 62 } 62 63 63 //....oooOO0OOooo........oooOO0OOooo........oo 64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 64 65 65 void RunAction::BeginOfRunAction(const G4Run*) 66 void RunAction::BeginOfRunAction(const G4Run*) 66 { << 67 { 67 // save Rndm status 68 // save Rndm status 68 ////G4RunManager::GetRunManager()->SetRandom 69 ////G4RunManager::GetRunManager()->SetRandomNumberStore(true); 69 CLHEP::HepRandom::showEngineStatus(); 70 CLHEP::HepRandom::showEngineStatus(); 70 71 71 fProcCounter = new ProcessesCount; 72 fProcCounter = new ProcessesCount; 72 fTotalCount = 0; 73 fTotalCount = 0; 73 << 74 74 fTruePL = fTruePL2 = fGeomPL = fGeomPL2 = 0. 75 fTruePL = fTruePL2 = fGeomPL = fGeomPL2 = 0.; 75 fLDispl = fLDispl2 = fPsiSpa = fPsiSpa2 = 0. 76 fLDispl = fLDispl2 = fPsiSpa = fPsiSpa2 = 0.; 76 fTetPrj = fTetPrj2 = 0.; 77 fTetPrj = fTetPrj2 = 0.; 77 fPhiCor = fPhiCor2 = 0.; 78 fPhiCor = fPhiCor2 = 0.; 78 << 79 79 // histograms << 80 //histograms 80 // 81 // 81 G4AnalysisManager* analysisManager = G4Analy 82 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance(); 82 if (analysisManager->IsActive()) { << 83 if ( analysisManager->IsActive() ) { 83 analysisManager->OpenFile(); 84 analysisManager->OpenFile(); 84 } << 85 } 85 } 86 } 86 87 87 //....oooOO0OOooo........oooOO0OOooo........oo 88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 88 89 89 void RunAction::CountProcesses(G4String procNa 90 void RunAction::CountProcesses(G4String procName) 90 { 91 { 91 // does the process already encounted ? << 92 //does the process already encounted ? 92 size_t nbProc = fProcCounter->size(); << 93 size_t nbProc = fProcCounter->size(); 93 size_t i = 0; << 94 size_t i = 0; 94 while ((i < nbProc) && ((*fProcCounter)[i]-> << 95 while ((i<nbProc)&&((*fProcCounter)[i]->GetName()!=procName)) i++; 95 i++; << 96 if (i == nbProc) fProcCounter->push_back( new OneProcessCount(procName)); 96 if (i == nbProc) fProcCounter->push_back(new << 97 97 98 (*fProcCounter)[i]->Count(); << 98 (*fProcCounter)[i]->Count(); 99 } 99 } 100 100 101 //....oooOO0OOooo........oooOO0OOooo........oo 101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 102 102 103 void RunAction::EndOfRunAction(const G4Run* aR 103 void RunAction::EndOfRunAction(const G4Run* aRun) 104 { 104 { 105 G4int NbOfEvents = aRun->GetNumberOfEvent(); 105 G4int NbOfEvents = aRun->GetNumberOfEvent(); 106 if (NbOfEvents == 0) return; 106 if (NbOfEvents == 0) return; 107 << 107 108 G4int prec = G4cout.precision(5); << 108 G4int prec = G4cout.precision(5); 109 << 109 110 G4Material* material = fDetector->GetMateria 110 G4Material* material = fDetector->GetMaterial(); 111 G4double density = material->GetDensity(); 111 G4double density = material->GetDensity(); 112 << 112 113 G4ParticleDefinition* particle = fPrimary->G << 113 G4ParticleDefinition* particle = 114 G4String Particle = particle->GetParticleNam << 114 fPrimary->GetParticleGun()->GetParticleDefinition(); >> 115 G4String Particle = particle->GetParticleName(); 115 G4double energy = fPrimary->GetParticleGun() 116 G4double energy = fPrimary->GetParticleGun()->GetParticleEnergy(); 116 G4cout << "\n The run consists of " << NbOfE << 117 G4cout << "\n The run consists of " << NbOfEvents << " "<< Particle << " of " 117 << G4BestUnit(energy, "Energy") << " << 118 << G4BestUnit(energy,"Energy") << " through " 118 << G4BestUnit(fDetector->GetBoxSize() << 119 << G4BestUnit(fDetector->GetBoxSize(),"Length") << " of " 119 << " (density: " << G4BestUnit(densit << 120 << material->GetName() << " (density: " 120 << 121 << G4BestUnit(density,"Volumic Mass") << ")" << G4endl; 121 // frequency of processes << 122 >> 123 //frequency of processes 122 G4cout << "\n Process calls frequency --->"; 124 G4cout << "\n Process calls frequency --->"; 123 for (size_t i = 0; i < fProcCounter->size(); << 125 for (size_t i=0; i< fProcCounter->size();i++) { 124 G4String procName = (*fProcCounter)[i]->Ge << 126 G4String procName = (*fProcCounter)[i]->GetName(); 125 G4int count = (*fProcCounter)[i]->GetCount << 127 G4int count = (*fProcCounter)[i]->GetCounter(); 126 G4cout << "\t" << procName << " = " << cou << 128 G4cout << "\t" << procName << " = " << count; 127 } 129 } 128 << 130 129 if (fTotalCount > 0) { 131 if (fTotalCount > 0) { 130 // compute path length and related quantit << 132 >> 133 //compute path length and related quantities 131 // 134 // 132 G4double MeanTPL = fTruePL / fTotalCount; << 135 G4double MeanTPL = fTruePL /fTotalCount; 133 G4double MeanTPL2 = fTruePL2 / fTotalCount << 136 G4double MeanTPL2 = fTruePL2/fTotalCount; 134 G4double rmsTPL = std::sqrt(std::fabs(Mean << 137 G4double rmsTPL = std::sqrt(std::fabs(MeanTPL2 - MeanTPL*MeanTPL)); 135 << 138 136 G4double MeanGPL = fGeomPL / fTotalCount; << 139 G4double MeanGPL = fGeomPL /fTotalCount; 137 G4double MeanGPL2 = fGeomPL2 / fTotalCount << 140 G4double MeanGPL2 = fGeomPL2/fTotalCount; 138 G4double rmsGPL = std::sqrt(std::fabs(Mean << 141 G4double rmsGPL = std::sqrt(std::fabs(MeanGPL2 - MeanGPL*MeanGPL)); 139 << 142 140 G4double MeanLaD = fLDispl / fTotalCount; << 143 G4double MeanLaD = fLDispl /fTotalCount; 141 G4double MeanLaD2 = fLDispl2 / fTotalCount << 144 G4double MeanLaD2 = fLDispl2/fTotalCount; 142 G4double rmsLaD = std::sqrt(std::fabs(Mean << 145 G4double rmsLaD = std::sqrt(std::fabs(MeanLaD2 - MeanLaD*MeanLaD)); 143 << 146 144 G4double MeanPsi = fPsiSpa / (fTotalCount) << 147 G4double MeanPsi = fPsiSpa /(fTotalCount); 145 G4double MeanPsi2 = fPsiSpa2 / (fTotalCoun << 148 G4double MeanPsi2 = fPsiSpa2/(fTotalCount); 146 G4double rmsPsi = std::sqrt(std::fabs(Mean << 149 G4double rmsPsi = std::sqrt(std::fabs(MeanPsi2 - MeanPsi*MeanPsi)); 147 << 150 148 G4double MeanTeta = fTetPrj / (2 * fTotalC << 151 G4double MeanTeta = fTetPrj /(2*fTotalCount); 149 G4double MeanTeta2 = fTetPrj2 / (2 * fTota << 152 G4double MeanTeta2 = fTetPrj2/(2*fTotalCount); 150 G4double rmsTeta = std::sqrt(std::fabs(Mea << 153 G4double rmsTeta = std::sqrt(std::fabs(MeanTeta2 - MeanTeta*MeanTeta)); 151 << 154 152 G4double MeanCorrel = fPhiCor / (fTotalCou << 155 G4double MeanCorrel = fPhiCor /(fTotalCount); 153 G4double MeanCorrel2 = fPhiCor2 / (fTotalC << 156 G4double MeanCorrel2 = fPhiCor2/(fTotalCount); 154 G4double rmsCorrel = std::sqrt(std::fabs(M << 157 G4double rmsCorrel = 155 << 158 std::sqrt(std::fabs(MeanCorrel2-MeanCorrel*MeanCorrel)); 156 G4cout << "\n\n truePathLength :\t" << G4B << 159 157 << G4BestUnit(rmsTPL, "Length") << << 160 G4cout << "\n\n truePathLength :\t" << G4BestUnit(MeanTPL,"Length") 158 << G4BestUnit(MeanGPL, "Length") << << 161 << " +- " << G4BestUnit( rmsTPL,"Length") 159 << "\n lateralDisplac :\t" << G4Bes << 162 << "\n geomPathLength :\t" << G4BestUnit(MeanGPL,"Length") 160 << G4BestUnit(rmsLaD, "Length") << << 163 << " +- " << G4BestUnit( rmsGPL,"Length") 161 << " +- " << rmsPsi / mrad << " mra << 164 << "\n lateralDisplac :\t" << G4BestUnit(MeanLaD,"Length") 162 << " (" << MeanPsi / deg << " deg" << 165 << " +- " << G4BestUnit( rmsLaD,"Length") 163 << " +- " << rmsPsi / deg << " deg) << 166 << "\n Psi :\t" << MeanPsi/mrad << " mrad" 164 << 167 << " +- " << rmsPsi /mrad << " mrad" 165 G4cout << "\n Theta_plane :\t" << rmsTe << 168 << " (" << MeanPsi/deg << " deg" 166 << " (" << rmsTeta / deg << " deg) << 169 << " +- " << rmsPsi /deg << " deg)" 167 << "\n phi correlation:\t" << MeanC << 170 << G4endl; 168 << " (std::cos(phi_pos - phi_dir)) << 171 169 << 172 G4cout << "\n Theta_plane :\t" << rmsTeta/mrad << " mrad" 170 // cross check from G4EmCalculator << 173 << " (" << rmsTeta/deg << " deg)" >> 174 << "\n phi correlation:\t" << MeanCorrel >> 175 << " +- " << rmsCorrel >> 176 << " (std::cos(phi_pos - phi_dir))" >> 177 << G4endl; >> 178 >> 179 >> 180 //cross check from G4EmCalculator 171 // 181 // 172 G4cout << "\n Verification from G4EmCalcul 182 G4cout << "\n Verification from G4EmCalculator. \n"; 173 << 183 174 G4EmCalculator emCal; 184 G4EmCalculator emCal; 175 << 185 176 // get transport mean free path (for multi << 186 //get transport mean free path (for multiple scattering) 177 G4double MSmfp = emCal.GetMeanFreePath(ene << 187 G4double MSmfp = emCal.GetMeanFreePath(energy,particle,"msc",material); 178 << 188 179 // get range from restricted dedx << 189 //get range from restricted dedx 180 G4double range = emCal.GetRangeFromRestric << 190 G4double range = emCal.GetRangeFromRestricteDEDX(energy,particle,material); 181 << 191 182 // effective facRange << 192 //effective facRange 183 G4double efFacrange = MeanTPL / std::max(M << 193 G4double efFacrange = MeanTPL/std::max(MSmfp, range); 184 if (MeanTPL / range >= 0.99) efFacrange = << 194 if (MeanTPL/range >= 0.99) efFacrange = 1.; 185 << 195 186 G4cout << "\n transport mean free path :\t << 196 G4cout << "\n transport mean free path :\t" << G4BestUnit(MSmfp,"Length") 187 << "\n range from restrict dE/dx:\t << 197 << "\n range from restrict dE/dx:\t" << G4BestUnit(range,"Length") 188 << "\n ---> effective facRange :\t << 198 << "\n ---> effective facRange :\t" << efFacrange 189 << 199 << G4endl; 190 G4cout << "\n compute theta0 from Highland << 200 191 << " (" << ComputeMscHighland(Mean << 201 G4cout << "\n compute theta0 from Highland :\t" 192 } << 202 << ComputeMscHighland(MeanTPL)/mrad << " mrad" 193 else << 203 << " (" << ComputeMscHighland(MeanTPL)/deg << " deg)" 194 G4cout << G4endl; << 204 << G4endl; 195 << 205 196 // restore default format << 206 } else 197 G4cout.precision(prec); << 207 G4cout<< G4endl; 198 << 208 199 // delete and remove all contents in fProcCo << 209 //restore default format 200 while (fProcCounter->size() > 0) { << 210 G4cout.precision(prec); 201 OneProcessCount* aProcCount = fProcCounter << 211 >> 212 // delete and remove all contents in fProcCounter >> 213 while (fProcCounter->size()>0){ >> 214 OneProcessCount* aProcCount=fProcCounter->back(); 202 fProcCounter->pop_back(); 215 fProcCounter->pop_back(); 203 delete aProcCount; 216 delete aProcCount; 204 } 217 } 205 delete fProcCounter; 218 delete fProcCounter; 206 << 219 207 // save histograms << 220 //save histograms 208 G4AnalysisManager* analysisManager = G4Analy << 221 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance(); 209 if (analysisManager->IsActive()) { << 222 if ( analysisManager->IsActive() ) { 210 analysisManager->Write(); << 223 analysisManager->Write(); 211 analysisManager->CloseFile(); << 224 analysisManager->CloseFile(); 212 } << 225 } 213 226 214 // show Rndm status 227 // show Rndm status 215 CLHEP::HepRandom::showEngineStatus(); 228 CLHEP::HepRandom::showEngineStatus(); 216 } 229 } 217 230 218 //....oooOO0OOooo........oooOO0OOooo........oo 231 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 219 232 220 G4double RunAction::ComputeMscHighland(G4doubl 233 G4double RunAction::ComputeMscHighland(G4double pathLength) 221 { 234 { 222 // compute the width of the Gaussian central << 235 //compute the width of the Gaussian central part of the MultipleScattering 223 // projected angular distribution. << 236 //projected angular distribution. 224 // Eur. Phys. Jour. C15 (2000) page 166, for << 237 //Eur. Phys. Jour. C15 (2000) page 166, formule 23.9 225 << 238 226 G4double t = pathLength / (fDetector->GetMat << 239 G4double t = pathLength/(fDetector->GetMaterial()->GetRadlen()); 227 if (t < DBL_MIN) return 0.; << 240 if (t < DBL_MIN) return 0.; 228 << 241 229 G4ParticleGun* particle = fPrimary->GetParti << 242 G4ParticleGun* particle = fPrimary->GetParticleGun(); 230 G4double T = particle->GetParticleEnergy(); << 243 G4double T = particle->GetParticleEnergy(); 231 G4double M = particle->GetParticleDefinition << 244 G4double M = particle->GetParticleDefinition()->GetPDGMass(); 232 G4double z = std::abs(particle->GetParticleD << 245 G4double z = std::abs(particle->GetParticleDefinition()->GetPDGCharge()/eplus); 233 << 246 234 G4double bpc = T * (T + 2 * M) / (T + M); << 247 G4double bpc = T*(T+2*M)/(T+M); 235 G4double teta0 = 13.6 * MeV * z * std::sqrt( << 248 G4double teta0 = 13.6*MeV*z*std::sqrt(t)*(1.+0.038*std::log(t))/bpc; 236 return teta0; << 249 return teta0; 237 } 250 } 238 251 239 //....oooOO0OOooo........oooOO0OOooo........oo 252 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 240 253