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 medical/fanoCavity/src/RunAction.cc << 26 // $Id: RunAction.cc,v 1.4 2009-01-22 18:34:06 maire Exp $ 27 /// \brief Implementation of the RunAction cla << 27 // GEANT4 tag $Name: geant4-09-04-patch-01 $ 28 // << 28 // 29 // << 30 //....oooOO0OOooo........oooOO0OOooo........oo 29 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 31 //....oooOO0OOooo........oooOO0OOooo........oo 30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 32 31 33 #include "RunAction.hh" 32 #include "RunAction.hh" 34 << 35 #include "DetectorConstruction.hh" 33 #include "DetectorConstruction.hh" 36 #include "HistoManager.hh" << 37 #include "PrimaryGeneratorAction.hh" 34 #include "PrimaryGeneratorAction.hh" 38 #include "Run.hh" << 35 #include "HistoManager.hh" 39 36 40 #include "G4Electron.hh" << 41 #include "G4EmCalculator.hh" << 42 #include "G4PhysicalConstants.hh" << 43 #include "G4Run.hh" 37 #include "G4Run.hh" 44 #include "G4RunManager.hh" 38 #include "G4RunManager.hh" 45 #include "G4SystemOfUnits.hh" << 46 #include "G4UnitsTable.hh" 39 #include "G4UnitsTable.hh" 47 #include "Randomize.hh" << 40 #include "G4EmCalculator.hh" >> 41 #include "G4Electron.hh" 48 42 >> 43 #include "Randomize.hh" 49 #include <iomanip> 44 #include <iomanip> 50 45 51 //....oooOO0OOooo........oooOO0OOooo........oo 46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 52 47 53 RunAction::RunAction(DetectorConstruction* det << 48 RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* kin, 54 : fDetector(det), fKinematic(kin), fRun(null << 49 HistoManager* histo) 55 { << 50 :detector(det),kinematic(kin),ProcCounter(0),histoManager(histo) 56 fHistoManager = new HistoManager(); << 51 { } 57 } << 58 52 59 //....oooOO0OOooo........oooOO0OOooo........oo 53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 60 54 61 RunAction::~RunAction() 55 RunAction::~RunAction() 62 { << 56 { } 63 delete fHistoManager; << 64 } << 65 57 66 //....oooOO0OOooo........oooOO0OOooo........oo 58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 67 59 68 G4Run* RunAction::GenerateRun() << 60 void RunAction::BeginOfRunAction(const G4Run* aRun) 69 { << 61 { 70 fRun = new Run(fDetector, fKinematic); << 62 G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl; 71 return fRun; << 63 >> 64 // save Rndm status >> 65 G4RunManager::GetRunManager()->SetRandomNumberStore(false); >> 66 CLHEP::HepRandom::showEngineStatus(); >> 67 >> 68 //geometry >> 69 // >> 70 wallThickness = detector->GetWallThickness(); >> 71 wallRadius = detector->GetWallRadius(); >> 72 mateWall = detector->GetWallMaterial(); >> 73 densityWall = mateWall->GetDensity(); >> 74 >> 75 cavityThickness = detector->GetCavityThickness(); >> 76 cavityRadius = detector->GetCavityRadius(); >> 77 surfaceCavity = pi*cavityRadius*cavityRadius; >> 78 volumeCavity = surfaceCavity*cavityThickness; >> 79 mateCavity = detector->GetCavityMaterial(); >> 80 densityCavity = mateCavity->GetDensity(); >> 81 massCavity = volumeCavity*densityCavity; >> 82 >> 83 //process counter >> 84 // >> 85 ProcCounter = new ProcessesCount; >> 86 >> 87 //kinetic energy of charged secondary a creation >> 88 // >> 89 Esecondary = Esecondary2 = 0.; >> 90 nbSec = 0; >> 91 >> 92 //charged particles and energy flow in cavity >> 93 // >> 94 PartFlowCavity[0] = PartFlowCavity[1] = 0; >> 95 EnerFlowCavity[0] = EnerFlowCavity[1] = 0.; >> 96 >> 97 //total energy deposit and charged track segment in cavity >> 98 // >> 99 EdepCavity = EdepCavity2 = trkSegmCavity = 0.; >> 100 nbEventCavity = 0; >> 101 >> 102 //stepLenth of charged particles >> 103 // >> 104 stepWall = stepWall2 = stepCavity = stepCavity2 =0.; >> 105 nbStepWall = nbStepCavity = 0; >> 106 >> 107 //histograms >> 108 // >> 109 histoManager->book(); >> 110 72 } 111 } >> 112 73 //....oooOO0OOooo........oooOO0OOooo........oo 113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 74 114 75 void RunAction::BeginOfRunAction(const G4Run* << 115 void RunAction::CountProcesses(G4String procName) 76 { 116 { 77 G4cout << "### Run " << aRun->GetRunID() << << 117 //does the process already encounted ? >> 118 size_t nbProc = ProcCounter->size(); >> 119 size_t i = 0; >> 120 while ((i<nbProc)&&((*ProcCounter)[i]->GetName()!=procName)) i++; >> 121 if (i == nbProc) ProcCounter->push_back( new OneProcessCount(procName)); 78 122 79 // do not save Rndm status << 123 (*ProcCounter)[i]->Count(); 80 G4RunManager::GetRunManager()->SetRandomNumb << 81 // if (!isMaster) CLHEP::HepRandom::showEngi << 82 } 124 } 83 125 84 //....oooOO0OOooo........oooOO0OOooo........oo 126 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 85 127 86 void RunAction::EndOfRunAction(const G4Run*) << 128 void RunAction::SurveyConvergence(G4int NbofEvents) 87 { << 129 { 88 // compute and print statistic << 130 if (NbofEvents == 0) return; 89 if (isMaster) fRun->EndOfRun(); << 131 >> 132 //mean kinetic energy of secondary electrons >> 133 // >> 134 G4double meanEsecond = 0.; >> 135 if (nbSec > 0) meanEsecond = Esecondary/nbSec; >> 136 G4double rateEmean = 0.; >> 137 // compute variation rate (%), iteration to iteration >> 138 if (oldEmean > 0.) rateEmean = 100*(meanEsecond/oldEmean - 1.); >> 139 oldEmean = meanEsecond; >> 140 >> 141 //beam energy fluence >> 142 // >> 143 G4double rBeam = wallRadius*(kinematic->GetBeamRadius()); >> 144 G4double surfaceBeam = pi*rBeam*rBeam; >> 145 G4double beamEnergy = kinematic->GetParticleGun()->GetParticleEnergy(); >> 146 >> 147 //total dose in cavity >> 148 // >> 149 G4double doseCavity = EdepCavity/massCavity; >> 150 G4double doseOverBeam = doseCavity*surfaceBeam/(NbofEvents*beamEnergy); >> 151 G4double rateDose = 0.; >> 152 // compute variation rate (%), iteration to iteration >> 153 if (oldDose > 0.) rateDose = 100*(doseOverBeam/oldDose - 1.); >> 154 oldDose = doseOverBeam; >> 155 >> 156 std::ios::fmtflags mode = G4cout.flags(); >> 157 G4cout.setf(std::ios::fixed,std::ios::floatfield); >> 158 G4int prec = G4cout.precision(3); >> 159 >> 160 G4cout << "\n ---> NbofEvents= " << NbofEvents >> 161 << " NbOfelectr= " << nbSec >> 162 << " Tkin= " << G4BestUnit(meanEsecond,"Energy") >> 163 << " (" << rateEmean << " %)" >> 164 << " NbOfelec in cav= " << PartFlowCavity[0] >> 165 << " Dose/EnFluence= " << G4BestUnit(doseOverBeam,"Surface/Mass") >> 166 << " (" << rateDose << " %)" >> 167 << G4endl; >> 168 >> 169 // reset default formats >> 170 G4cout.setf(mode,std::ios::floatfield); >> 171 G4cout.precision(prec); >> 172 } 90 173 91 // save histograms << 174 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 92 G4AnalysisManager* analysisManager = G4Analy << 93 if (analysisManager->IsActive()) { << 94 analysisManager->Write(); << 95 analysisManager->CloseFile(); << 96 } << 97 175 >> 176 void RunAction::EndOfRunAction(const G4Run* aRun) >> 177 { >> 178 std::ios::fmtflags mode = G4cout.flags(); >> 179 G4cout.setf(std::ios::fixed,std::ios::floatfield); >> 180 >> 181 G4int NbofEvents = aRun->GetNumberOfEvent(); >> 182 if (NbofEvents == 0) return; >> 183 >> 184 //run conditions >> 185 // >> 186 G4ParticleDefinition* particle = kinematic->GetParticleGun() >> 187 ->GetParticleDefinition(); >> 188 G4String partName = particle->GetParticleName(); >> 189 G4double energy = kinematic->GetParticleGun()->GetParticleEnergy(); >> 190 >> 191 G4cout << "\n ======================== run summary ======================\n"; >> 192 >> 193 G4int prec = G4cout.precision(3); >> 194 >> 195 G4cout << "\n The run consists of " << NbofEvents << " "<< partName << " of " >> 196 << G4BestUnit(energy,"Energy") << " through 2*" >> 197 << G4BestUnit(wallThickness,"Length") << " of " >> 198 << mateWall->GetName() << " (density: " >> 199 << G4BestUnit(densityWall,"Volumic Mass") << ")" << G4endl; >> 200 >> 201 G4cout << "\n the cavity is " >> 202 << G4BestUnit(cavityThickness,"Length") << " of " >> 203 << mateCavity->GetName() << " (density: " >> 204 << G4BestUnit(densityCavity,"Volumic Mass") << "); Mass = " >> 205 << G4BestUnit(massCavity,"Mass") << G4endl; >> 206 >> 207 G4cout << "\n ============================================================\n"; >> 208 >> 209 //frequency of processes >> 210 // >> 211 G4cout << "\n Process calls frequency --->"; >> 212 for (size_t i=0; i< ProcCounter->size();i++) { >> 213 G4String procName = (*ProcCounter)[i]->GetName(); >> 214 G4int count = (*ProcCounter)[i]->GetCounter(); >> 215 G4cout << " " << procName << "= " << count; >> 216 } >> 217 G4cout << G4endl; >> 218 >> 219 //extract cross sections with G4EmCalculator >> 220 // >> 221 G4EmCalculator emCalculator; >> 222 G4cout << "\n Gamma crossSections in wall material :"; >> 223 G4double sumc = 0.0; >> 224 for (size_t i=0; i< ProcCounter->size();i++) { >> 225 G4String procName = (*ProcCounter)[i]->GetName(); >> 226 G4double massSigma = >> 227 emCalculator.ComputeCrossSectionPerVolume(energy,particle, >> 228 procName,mateWall)/densityWall; >> 229 if (massSigma > 0.) { >> 230 sumc += massSigma; >> 231 G4cout << " " << procName << "= " >> 232 << G4BestUnit(massSigma, "Surface/Mass"); >> 233 } >> 234 } >> 235 G4cout << " --> total= " >> 236 << G4BestUnit(sumc, "Surface/Mass") << G4endl; >> 237 >> 238 //mean kinetic energy of secondary electrons >> 239 // >> 240 if (nbSec == 0) return; >> 241 G4double meanEsecond = Esecondary/nbSec, meanEsecond2 = Esecondary2/nbSec; >> 242 G4double varianceEsec = meanEsecond2 - meanEsecond*meanEsecond; >> 243 G4double dToverT = 0.; >> 244 if (varianceEsec>0.) dToverT = std::sqrt(varianceEsec/nbSec)/meanEsecond; >> 245 G4double csdaRange = >> 246 emCalculator.GetCSDARange(meanEsecond,G4Electron::Electron(),mateWall); >> 247 >> 248 G4cout.precision(4); >> 249 G4cout >> 250 << "\n Mean energy of secondary e- = " << G4BestUnit(meanEsecond,"Energy") >> 251 << " +- " << 100*dToverT << " %" >> 252 << " (--> range in wall material = " << G4BestUnit(csdaRange,"Length") >> 253 << ")" >> 254 << G4endl; >> 255 >> 256 //compute mass energy transfer coefficient >> 257 // >> 258 G4double massTransfCoef = sumc*meanEsecond/energy; >> 259 >> 260 G4cout << " Mass_energy_transfer coef: " >> 261 << G4BestUnit(massTransfCoef, "Surface/Mass") >> 262 << G4endl; >> 263 >> 264 //stopping power from EmCalculator >> 265 // >> 266 G4double dedxWall = >> 267 emCalculator.GetDEDX(meanEsecond,G4Electron::Electron(),mateWall); >> 268 dedxWall /= densityWall; >> 269 G4double dedxCavity = >> 270 emCalculator.GetDEDX(meanEsecond,G4Electron::Electron(),mateCavity); >> 271 dedxCavity /= densityCavity; >> 272 >> 273 G4cout >> 274 << "\n StoppingPower in wall = " >> 275 << G4BestUnit(dedxWall,"Energy*Surface/Mass") >> 276 << "\n in cavity = " >> 277 << G4BestUnit(dedxCavity,"Energy*Surface/Mass") >> 278 << G4endl; >> 279 >> 280 //charged particle flow in cavity >> 281 // >> 282 G4cout >> 283 << "\n Charged particle flow in cavity :" >> 284 << "\n Enter --> nbParticles = " << PartFlowCavity[0] >> 285 << "\t Energy = " << G4BestUnit (EnerFlowCavity[0], "Energy") >> 286 << "\n Exit --> nbParticles = " << PartFlowCavity[1] >> 287 << "\t Energy = " << G4BestUnit (EnerFlowCavity[1], "Energy") >> 288 << G4endl; >> 289 >> 290 if (PartFlowCavity[0] == 0) return; >> 291 >> 292 //beam energy fluence >> 293 // >> 294 G4double rBeam = wallRadius*(kinematic->GetBeamRadius()); >> 295 G4double surfaceBeam = pi*rBeam*rBeam; >> 296 >> 297 //error on Edep in cavity >> 298 // >> 299 if (nbEventCavity == 0) return; >> 300 G4double meanEdep = EdepCavity/nbEventCavity; >> 301 G4double meanEdep2 = EdepCavity2/nbEventCavity; >> 302 G4double varianceEdep = meanEdep2 - meanEdep*meanEdep; >> 303 G4double dEoverE = 0.; >> 304 if(varianceEdep>0.) dEoverE = std::sqrt(varianceEdep/nbEventCavity)/meanEdep; >> 305 >> 306 //total dose in cavity >> 307 // >> 308 G4double doseCavity = EdepCavity/massCavity; >> 309 G4double doseOverBeam = doseCavity*surfaceBeam/(NbofEvents*energy); >> 310 >> 311 //track length in cavity >> 312 G4double meantrack = trkSegmCavity/PartFlowCavity[0]; >> 313 >> 314 G4cout.precision(4); >> 315 G4cout >> 316 << "\n Total edep in cavity = " << G4BestUnit(EdepCavity,"Energy") >> 317 << " +- " << 100*dEoverE << " %" >> 318 << "\t Total charged trackLength = " << G4BestUnit(trkSegmCavity,"Length") >> 319 << " (mean value = " << G4BestUnit(meantrack,"Length") << ")" >> 320 << "\n Total dose in cavity = " << doseCavity/(MeV/mg) << " MeV/mg" >> 321 << "\n Dose/EnergyFluence = " << G4BestUnit(doseOverBeam,"Surface/Mass") >> 322 << G4endl; >> 323 >> 324 //ratio simulation/theory >> 325 // >> 326 G4double ratio = doseOverBeam/massTransfCoef; >> 327 G4double error = ratio*std::sqrt(dEoverE*dEoverE + dToverT*dToverT); >> 328 >> 329 G4cout.precision(5); >> 330 G4cout >> 331 << "\n (Dose/EnergyFluence)/Mass_energy_transfer = " << ratio >> 332 << " +- " << error << G4endl; >> 333 >> 334 //compute mean step size of charged particles >> 335 // >> 336 stepWall /= nbStepWall; stepWall2 /= nbStepWall; >> 337 G4double rms = stepWall2 - stepWall*stepWall; >> 338 if (rms>0.) rms = std::sqrt(rms); else rms = 0.; >> 339 >> 340 G4cout.precision(4); >> 341 G4cout >> 342 << "\n StepSize of ch. tracks in wall = " >> 343 << G4BestUnit(stepWall,"Length") << " +- " << G4BestUnit( rms,"Length") >> 344 << "\t (nbSteps/track = " << double(nbStepWall)/nbSec << ")"; >> 345 >> 346 stepCavity /= nbStepCavity; stepCavity2 /= nbStepCavity; >> 347 rms = stepCavity2 - stepCavity*stepCavity; >> 348 if (rms>0.) rms = std::sqrt(rms); else rms = 0.; >> 349 >> 350 G4cout >> 351 << "\n StepSize of ch. tracks in cavity = " >> 352 << G4BestUnit(stepCavity,"Length") << " +- " << G4BestUnit( rms,"Length") >> 353 << "\t (nbSteps/track = " << double(nbStepCavity)/PartFlowCavity[0] << ")"; >> 354 >> 355 G4cout << G4endl; >> 356 >> 357 // reset default formats >> 358 G4cout.setf(mode,std::ios::floatfield); >> 359 G4cout.precision(prec); >> 360 >> 361 // delete and remove all contents in ProcCounter >> 362 while (ProcCounter->size()>0){ >> 363 OneProcessCount* aProcCount=ProcCounter->back(); >> 364 ProcCounter->pop_back(); >> 365 delete aProcCount; >> 366 } >> 367 delete ProcCounter; >> 368 >> 369 // save histograms >> 370 histoManager->save(); >> 371 98 // show Rndm status 372 // show Rndm status 99 // if(!isMaster )CLHEP::HepRandom::showEngin << 373 CLHEP::HepRandom::showEngineStatus(); 100 } 374 } 101 375 102 //....oooOO0OOooo........oooOO0OOooo........oo 376 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 103 377