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/TestEm17/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 "MuCrossSections.hh" 38 #include "PrimaryGeneratorAction.hh" 39 40 #include "G4EmCalculator.hh" 41 #include "G4PhysicalConstants.hh" 42 #include "G4ProductionCutsTable.hh" 43 #include "G4Run.hh" 44 #include "G4RunManager.hh" 45 #include "G4SystemOfUnits.hh" 46 #include "G4UnitsTable.hh" 47 #include "Randomize.hh" 48 49 //....oooOO0OOooo........oooOO0OOooo........oo 50 51 RunAction::RunAction(DetectorConstruction* det 52 : G4UserRunAction(), fDetector(det), fPrimar 53 { 54 fMucs = new MuCrossSections(); 55 } 56 57 //....oooOO0OOooo........oooOO0OOooo........oo 58 59 RunAction::~RunAction() 60 { 61 delete fMucs; 62 } 63 64 //....oooOO0OOooo........oooOO0OOooo........oo 65 66 void RunAction::BeginOfRunAction(const G4Run* 67 { 68 G4cout << "### Run " << aRun->GetRunID() << 69 70 // save Rndm status 71 CLHEP::HepRandom::showEngineStatus(); 72 73 fProcCounter = new ProcessesCount(); 74 fHistoManager->Book(); 75 } 76 77 //....oooOO0OOooo........oooOO0OOooo........oo 78 79 void RunAction::CountProcesses(const G4String& 80 { 81 // does the process already encounted ? 82 size_t n = fProcCounter->size(); 83 for (size_t i = 0; i < n; ++i) { 84 if ((*fProcCounter)[i]->GetName() == procN 85 (*fProcCounter)[i]->Count(); 86 return; 87 } 88 } 89 OneProcessCount* count = new OneProcessCount 90 count->Count(); 91 fProcCounter->push_back(count); 92 } 93 94 //....oooOO0OOooo........oooOO0OOooo........oo 95 96 void RunAction::EndOfRunAction(const G4Run* aR 97 { 98 G4int NbOfEvents = aRun->GetNumberOfEvent(); 99 if (NbOfEvents == 0) return; 100 101 // std::ios::fmtflags mode = G4cout.flags() 102 G4int prec = G4cout.precision(2); 103 104 const G4Material* material = fDetector->GetM 105 G4double length = fDetector->GetSize(); 106 G4double density = material->GetDensity(); 107 108 G4String particle = fPrimary->GetParticleGun 109 G4double energy = fPrimary->GetParticleGun() 110 111 G4cout << "\n The run consists of " << NbOfE 112 << G4BestUnit(energy, "Energy") << " 113 << material->GetName() << " (density: 114 << G4endl; 115 116 // total number of process calls 117 G4double countTot = 0.; 118 G4cout << "\n Number of process calls --->"; 119 for (size_t i = 0; i < fProcCounter->size(); 120 G4String procName = (*fProcCounter)[i]->Ge 121 if (procName != "Transportation") { 122 G4int count = (*fProcCounter)[i]->GetCou 123 G4cout << "\t" << procName << " : " << c 124 countTot += count; 125 } 126 } 127 128 // compute totalCrossSection, meanFreePath a 129 // 130 G4double totalCrossSection = countTot / (NbO 131 G4double MeanFreePath = 1. / totalCrossSecti 132 G4double massCrossSection = totalCrossSectio 133 134 G4cout.precision(5); 135 G4cout << "\n Simulation: " 136 << "total CrossSection = " << totalCr 137 << "\t MeanFreePath = " << G4BestUnit 138 << "\t massicCrossSection = " << mass 139 140 // compute theoretical predictions 141 // 142 if (particle == "mu+" || particle == "mu-") 143 totalCrossSection = 0.; 144 for (size_t i = 0; i < fProcCounter->size( 145 G4String procName = (*fProcCounter)[i]-> 146 if (procName != "Transportation") { 147 totalCrossSection += ComputeTheory(pro 148 FillCrossSectionHisto(procName, NbOfEv 149 } 150 } 151 152 MeanFreePath = 1. / totalCrossSection; 153 massCrossSection = totalCrossSection / den 154 155 G4cout << " Theory: " 156 << "total CrossSection = " << total 157 << "\t MeanFreePath = " << G4BestUn 158 << "\t massicCrossSection = " << ma 159 } 160 161 // G4cout.setf(mode,std::ios::floatfield); 162 G4cout.precision(prec); 163 164 // delete and remove all contents in fProcCo 165 size_t n = fProcCounter->size(); 166 for (size_t i = 0; i < n; ++i) { 167 delete (*fProcCounter)[i]; 168 } 169 delete fProcCounter; 170 171 fHistoManager->Save(); 172 173 // show Rndm status 174 // CLHEP::HepRandom::showEngineStatus(); 175 } 176 177 //....oooOO0OOooo........oooOO0OOooo........oo 178 179 G4double RunAction::ComputeTheory(const G4Stri 180 { 181 const G4Material* material = fDetector->GetM 182 G4double ekin = fPrimary->GetParticleGun()-> 183 G4double particleMass = fPrimary->GetParticl 184 185 G4int id = 0; 186 G4double cut = 1.e-10 * ekin; 187 if (process == "muIoni") { 188 id = 11; 189 cut = GetEnergyCut(material, 1); 190 } 191 else if (process == "muPairProd") { 192 id = 12; 193 cut = 2 * (GetEnergyCut(material, 1) + ele 194 } 195 else if (process == "muBrems") { 196 id = 13; 197 cut = GetEnergyCut(material, 0); 198 } 199 else if (process == "muonNuclear") { 200 id = 14; 201 cut = 100 * MeV; 202 } 203 else if (process == "muToMuonPairProd") { 204 id = 18; 205 cut = 2 * particleMass; 206 } 207 if (id == 0) { 208 return 0.; 209 } 210 211 G4int nbOfBins = 100; 212 // G4double binMin = -10.; 213 G4double binMin = std::log10(cut / ekin); 214 G4double binMax = 0.; 215 G4double binWidth = (binMax - binMin) / G4do 216 217 // create histo for theoretical crossSection 218 // 219 G4AnalysisManager* analysisManager = G4Analy 220 221 G4H1* histoTh = 0; 222 if (fHistoManager->HistoExist(id)) { 223 histoTh = analysisManager->GetH1(fHistoMan 224 nbOfBins = fHistoManager->GetNbins(id); 225 binMin = fHistoManager->GetVmin(id); 226 binMax = fHistoManager->GetVmax(id); 227 binWidth = fHistoManager->GetBinWidth(id); 228 } 229 230 // compute and plot differential crossSectio 231 // compute and return integrated crossSectio 232 //(note: to compare with simulation, the int 233 // of the energy cut.) 234 // 235 G4double lgeps, etransf, sigmaE, dsigma; 236 G4double sigmaTot = 0.; 237 const G4double ln10 = std::log(10.); 238 G4double length = fDetector->GetSize(); 239 240 // G4cout << "MU: " << process << " E= " << 241 // <<" binMin= " << binMin << " binW 242 243 for (G4int ibin = 0; ibin < nbOfBins; ibin++ 244 lgeps = binMin + (ibin + 0.5) * binWidth; 245 etransf = ekin * std::pow(10., lgeps); 246 sigmaE = fMucs->CR_Macroscopic(process, ma 247 dsigma = sigmaE * etransf * binWidth * ln1 248 if (etransf > cut) sigmaTot += dsigma; 249 if (histoTh) { 250 G4double NbProcess = NbOfMu * length * d 251 histoTh->fill(lgeps, NbProcess); 252 } 253 } 254 255 // return integrated crossSection 256 // 257 return sigmaTot; 258 } 259 260 //....oooOO0OOooo........oooOO0OOooo........oo 261 262 void RunAction::FillCrossSectionHisto(const G4 263 { 264 const G4Material* material = fDetector->GetM 265 G4double ekin = fPrimary->GetParticleGun()-> 266 G4ParticleDefinition* particle = fPrimary->G 267 G4double particleMass = particle->GetPDGMass 268 269 G4EmCalculator emCal; 270 271 G4int id = 0; 272 G4double cut = 1.e-10 * ekin; 273 if (process == "muIoni") { 274 id = 21; 275 cut = GetEnergyCut(material, 1); 276 } 277 else if (process == "muPairProd") { 278 id = 22; 279 cut = 2 * (GetEnergyCut(material, 1) + ele 280 } 281 else if (process == "muBrems") { 282 id = 23; 283 cut = GetEnergyCut(material, 0); 284 } 285 else if (process == "muonNuclear") { 286 id = 24; 287 cut = 100 * MeV; 288 } 289 else if (process == "muToMuonPairProd") { 290 id = 28; 291 cut = 2 * particleMass; 292 } 293 if (id == 0) { 294 return; 295 } 296 297 G4int nbOfBins = 100; 298 G4double binMin = cut; 299 G4double binMax = ekin; 300 G4double binWidth = (binMax - binMin) / G4do 301 302 G4AnalysisManager* analysisManager = G4Analy 303 304 G4H1* histoTh = 0; 305 if (fHistoManager->HistoExist(id)) { 306 histoTh = analysisManager->GetH1(fHistoMan 307 nbOfBins = fHistoManager->GetNbins(id); 308 binMin = fHistoManager->GetVmin(id); 309 binMax = fHistoManager->GetVmax(id); 310 binWidth = fHistoManager->GetBinWidth(id); 311 } 312 313 G4double sigma, primaryEnergy; 314 315 for (G4int ibin = 0; ibin < nbOfBins; ibin++ 316 primaryEnergy = binMin + (ibin + 0.5) * bi 317 sigma = emCal.GetCrossSectionPerVolume(pri 318 if (histoTh) { 319 histoTh->fill(primaryEnergy, sigma); 320 } 321 } 322 } 323 324 //....oooOO0OOooo........oooOO0OOooo........oo 325 326 G4double RunAction::GetEnergyCut(const G4Mater 327 { 328 G4ProductionCutsTable* table = G4ProductionC 329 330 size_t index = 0; 331 while ((table->GetMaterialCutsCouple(index)- 332 && (index < table->GetTableSize())) 333 index++; 334 335 return (*(table->GetEnergyCutsVector(idParti 336 } 337 338 //....oooOO0OOooo........oooOO0OOooo........oo 339