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