| 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.3 2007/10/29 12:36:26 maire Exp $
27 /// \brief Implementation of the RunAction cla << 27 // GEANT4 tag $Name: geant4-09-01-patch-02 $
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*(dEoverE + 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