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