Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
extended/medical/fanoCavity2/src/Run.cc
이 파일의 문서화 페이지로 가기
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
28 //
29 // $Id: Run.cc 71035 2013-06-10 09:17:35Z gcosmo $
30 //
31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
33 
34 #include "Run.hh"
35 #include "DetectorConstruction.hh"
36 #include "PrimaryGeneratorAction.hh"
37 #include "HistoManager.hh"
38 
39 #include "G4Run.hh"
40 #include "G4RunManager.hh"
41 #include "G4UnitsTable.hh"
42 #include "G4EmCalculator.hh"
43 #include "G4Electron.hh"
44 
45 #include "G4SystemOfUnits.hh"
46 #include "Randomize.hh"
47 #include <iomanip>
48 
49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
50 
52 :G4Run(),fDetector(det), fKinematic(kin), fProcCounter(0),
53  fEdepCavity(0.), fEdepCavity2(0.),
54  fTrkSegmCavity(0.), fNbEventCavity(0),
55  fStepWall(0.), fStepWall2(0.),
56  fStepCavity(0.), fStepCavity2(0.),
57  fNbStepWall(0), fNbStepCavity(0),
58  fEnergyGun(0.), fMassWall(0.),
59  fMassCavity(0.),fIsMaster(isMaster)
60 {
61 
62  //run conditions
63  //
64  G4ParticleDefinition* particleGun
66  G4String partName = particleGun->GetParticleName();
68 
69  //geometry : effective wall volume
70  //
71  G4double cavityThickness = fDetector->GetCavityThickness();
72  G4Material* mateCavity = fDetector->GetCavityMaterial();
73  G4double densityCavity = mateCavity->GetDensity();
74  fMassCavity = cavityThickness*densityCavity;
75 
76  G4double wallThickness = fDetector->GetWallThickness();
77  G4Material* mateWall = fDetector->GetWallMaterial();
78  G4double densityWall = mateWall->GetDensity();
79 
80  G4EmCalculator emCal;
81  G4double RangeWall = emCal.GetCSDARange(fEnergyGun,particleGun,mateWall);
82  G4double factor = 1.2;
83  G4double effWallThick = factor*RangeWall;
84  if ((effWallThick > wallThickness)||(effWallThick <= 0.))
85  effWallThick = wallThickness;
86  fMassWall = 2*effWallThick*densityWall;
87 
88  G4double massTotal = fMassWall + fMassCavity;
89  G4double fMassWallRatio = fMassWall/massTotal;
90  fKinematic->RunInitialisation(effWallThick, fMassWallRatio );
91 
92  G4double massRatio = fMassCavity/fMassWall;
93 
94  //check radius
95  //
96  G4double worldRadius =fDetector->GetWorldRadius();
97  G4double RangeCavity =emCal.GetCSDARange(fEnergyGun,particleGun,mateCavity);
98 
99  //G4String partName = particleGun->GetParticleName();
100 
101 
102  std::ios::fmtflags mode = G4cout.flags();
103  G4cout.setf(std::ios::fixed,std::ios::floatfield);
104  G4int prec = G4cout.precision(3);
105 
106  G4cout << "\n ===================== run conditions =====================\n";
107 
108  G4cout << "\n The run will be " << numberOfEvent << " "<< partName << " of "
109  << G4BestUnit(fEnergyGun,"Energy") << " through 2*"
110  << G4BestUnit(effWallThick,"Length") << " of "
111  << mateWall->GetName() << " (density: "
112  << G4BestUnit(densityWall,"Volumic Mass") << "); Mass/cm2 = "
113  << G4BestUnit(fMassWall*cm2,"Mass")
114  << "\n csdaRange: " << G4BestUnit(RangeWall,"Length") << G4endl;
115 
116  G4cout << "\n the cavity is "
117  << G4BestUnit(cavityThickness,"Length") << " of "
118  << mateCavity->GetName() << " (density: "
119  << G4BestUnit(densityCavity,"Volumic Mass") << "); Mass/cm2 = "
120  << G4BestUnit(fMassCavity*cm2,"Mass")
121  << " --> massRatio = "<< std::setprecision(6) << massRatio << G4endl;
122 
123  G4cout.precision(3);
124  G4cout << " World radius: " << G4BestUnit(worldRadius,"Length")
125  << "; range in cavity: " << G4BestUnit(RangeCavity,"Length")
126  << G4endl;
127 
128  G4cout << "\n ==========================================================\n";
129 
130  //stopping power from EmCalculator
131  //
132  G4double dedxWall =
133  emCal.GetDEDX(fEnergyGun,G4Electron::Electron(),mateWall);
134  dedxWall /= densityWall;
135  G4double dedxCavity =
136  emCal.GetDEDX(fEnergyGun,G4Electron::Electron(),mateCavity);
137  dedxCavity /= densityCavity;
138 
139  G4cout << std::setprecision(4)
140  << "\n StoppingPower in wall = "
141  << G4BestUnit(dedxWall,"Energy*Surface/Mass")
142  << "\n in cavity = "
143  << G4BestUnit(dedxCavity,"Energy*Surface/Mass")
144  << G4endl;
145 
146  //process counter
147  //
149 
150  //charged particles and energy flow in cavity
151  //
152  fPartFlowCavity[0] = fPartFlowCavity[1] = 0;
153  fEnerFlowCavity[0] = fEnerFlowCavity[1] = 0.;
154 
155  //total energy deposit and charged track segment in cavity
156  //
158  fNbEventCavity = 0;
159 
160  //stepLenth of charged particles
161  //
164 
165 
166  // reset default formats
167  G4cout.setf(mode,std::ios::floatfield);
168  G4cout.precision(prec);
169 
170  //histograms
171  //
172  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
173  if ( analysisManager->IsActive() ) {
174  analysisManager->OpenFile();
175  }
176 
177 }
178 
179 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
180 
181 Run::~Run()
182 {
183 }
184 
185 
186 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
187 
188 void Run::CountProcesses(G4String procName)
189 {
190  //does the process already encounted ?
191  size_t nbProc = fProcCounter->size();
192  size_t i = 0;
193  while ((i<nbProc)&&((*fProcCounter)[i]->GetName()!=procName)) i++;
194  if (i == nbProc) fProcCounter->push_back( new OneProcessCount(procName));
195 
196  (*fProcCounter)[i]->Count();
197 }
198 
199 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
200 
201 void Run::SurveyConvergence(G4int NbofEvents)
202 {
203  if (NbofEvents == 0) return;
204 
205 
206 
207 
208  //beam fluence
209  //
210  G4int Nwall = fKinematic->GetWallCount();
211  G4int Ncavity = fKinematic->GetCavityCount();
212  G4double Iwall = Nwall/fMassWall;
213  G4double Icavity = Ncavity/fMassCavity;
214  G4double Iratio = Icavity/Iwall;
215  G4double Itot = NbofEvents/(fMassWall+fMassCavity);
216  G4double energyFluence = fEnergyGun*Itot;
217 
218  //total dose in cavity
219  //
220  G4double doseCavity = fEdepCavity/fMassCavity;
221  G4double ratio = doseCavity/energyFluence;
222  G4double err = 100*(ratio-1.);
223 
224  std::ios::fmtflags mode = G4cout.flags();
225  G4cout.setf(std::ios::fixed,std::ios::floatfield);
226  G4int prec = G4cout.precision(5);
227 
228  G4cout << "\n--->evntNb= " << NbofEvents
229  << " Nwall= " << Nwall
230  << " Ncav= " << Ncavity
231  << " Ic/Iw= " << Iratio
232  << " Ne-_cav= " << fPartFlowCavity[0]
233  << " doseCavity/Ebeam= " << ratio
234  << " (100*(ratio-1) = " << err << " %)"
235  << G4endl;
236 
237  // reset default formats
238  G4cout.setf(mode,std::ios::floatfield);
239  G4cout.precision(prec);
240 }
241 
242 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
243 
244 void Run::EndOfRun()
245 { // Only call by Master thread
246 
247 
248  std::ios::fmtflags mode = G4cout.flags();
249  G4cout.setf(std::ios::fixed,std::ios::floatfield);
250  G4int prec = G4cout.precision(3);
251 
252 
253  if (numberOfEvent == 0) return;
254 
255  //frequency of processes
256  //
257  G4cout << "\n Process calls frequency --->";
258  for (size_t i=0; i< fProcCounter->size();i++) {
259  G4String procName = (*fProcCounter)[i]->GetName();
260  G4int count = (*fProcCounter)[i]->GetCounter();
261  G4cout << " " << procName << "= " << count;
262  }
263  G4cout << G4endl;
264 
265  //charged particle flow in cavity
266  //
267  G4cout
268  << "\n Charged particle flow in cavity :"
269  << "\n Enter --> nbParticles = " << fPartFlowCavity[0]
270  << "\t Energy = " << G4BestUnit (fEnerFlowCavity[0], "Energy")
271  << "\n Exit --> nbParticles = " << fPartFlowCavity[1]
272  << "\t Energy = " << G4BestUnit (fEnerFlowCavity[1], "Energy")
273  << G4endl;
274 
275  if (fPartFlowCavity[0] == 0) return;
276 
277  G4int Nwall = fKinematic->GetWallCount();
278  G4int Ncavity = fKinematic->GetCavityCount();
279 
280 
281  G4double Iwall = Nwall/fMassWall;
282  G4double Icavity = Ncavity/fMassCavity;
283  G4double Iratio = Icavity/Iwall;
285  G4double energyFluence = fEnergyGun*Itot;
286 
287  G4cout.precision(5);
288  G4cout
289  << "\n beamFluence in wall = " << Nwall
290  << "\t in cavity = " << Ncavity
291  << "\t Icav/Iwall = " << Iratio
292  << "\t energyFluence = " << energyFluence/(MeV*cm2/mg) << " MeV*cm2/mg"
293  << G4endl;
294 
295  //error on Edep in cavity
296  //
297  if (fNbEventCavity == 0) return;
300  G4double varianceEdep = meanEdep2 - meanEdep*meanEdep;
301  G4double dEoverE = 0.;
302  if(varianceEdep>0.) dEoverE = std::sqrt(varianceEdep/fNbEventCavity)/meanEdep;
303 
304  //total dose in cavity
305  //
306  G4double doseCavity = fEdepCavity/fMassCavity;
307 
308  G4double ratio = doseCavity/energyFluence, error = ratio*dEoverE;
309 
310  G4cout
311  << "\n Total edep in cavity = " << G4BestUnit(fEdepCavity,"Energy")
312  << " +- " << 100*dEoverE << " %"
313  << "\n Total dose in cavity = " << doseCavity/(MeV*cm2/mg) << " MeV*cm2/mg"
314  << " +- " << 100*dEoverE << " %"
315  << "\n\n DoseCavity/EnergyFluence = " << ratio
316  << " +- " << error << G4endl;
317 
318 
319  //track length in cavity
320  G4double meantrack = fTrkSegmCavity/fPartFlowCavity[0];
321 
322  G4cout.precision(4);
323  G4cout
324  << "\n Total charged trackLength in cavity = "
325  << G4BestUnit(fTrkSegmCavity,"Length")
326  << " (mean value = " << G4BestUnit(meantrack,"Length") << ")"
327  << G4endl;
328 
329  //compute mean step size of charged particles
330  //
333  if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
334  G4double nbTrackWall = fKinematic->GetWallCount();
335 
336  G4cout
337  << "\n StepSize of ch. tracks in wall = "
338  << G4BestUnit(fStepWall,"Length") << " +- " << G4BestUnit( rms,"Length")
339  << "\t (nbSteps/track = " << double(fNbStepWall)/nbTrackWall << ")";
340 
343  if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
344 
345  G4cout
346  << "\n StepSize of ch. tracks in cavity = "
347  << G4BestUnit(fStepCavity,"Length") << " +- " << G4BestUnit( rms,"Length")
348  << "\t (nbSteps/track = " <<double(fNbStepCavity)/fPartFlowCavity[0] << ")";
349 
350  G4cout << G4endl;
351 
352  // reset default formats
353  G4cout.setf(mode,std::ios::floatfield);
354  G4cout.precision(prec);
355 
356  // delete and remove all contents in fProcCounter
357  while (fProcCounter->size()>0){
358  OneProcessCount* aProcCount=fProcCounter->back();
359  fProcCounter->pop_back();
360  delete aProcCount;
361  }
362  delete fProcCounter;
363 
364  // show Rndm status
366 }
367 
368 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
369 
370 void Run::Merge(const G4Run* run) {
371 
372  const Run* localRun = static_cast<const Run*>(run);
373 
374  // Merge Run variables
375  fPartFlowCavity[0]+= localRun->fPartFlowCavity[0];
376  fPartFlowCavity[1]+= localRun->fPartFlowCavity[1];
377  fEnerFlowCavity[0]+= localRun->fEnerFlowCavity[0];
378  fEnerFlowCavity[1]+= localRun->fEnerFlowCavity[1];
379  fEdepCavity += localRun->fEdepCavity;
380  fEdepCavity2 += localRun->fEdepCavity2;
381  fTrkSegmCavity += localRun->fTrkSegmCavity;
382  fNbEventCavity += localRun->fNbEventCavity;
383  fStepWall += localRun->fStepWall;
384  fStepWall2 += localRun->fStepWall2;
385  fStepCavity += localRun->fStepCavity;
386  fStepCavity2 += localRun->fStepCavity2;
387  fNbStepWall += localRun->fNbStepWall;
388  fNbStepCavity += localRun->fNbStepCavity;
389 
390  // Merge PrimaryGenerator variables
393 
394  // Merge ProcessCount varaibles
395  std::vector<OneProcessCount*>::iterator it;
396  for (it = localRun->fProcCounter->begin();
397  it !=localRun->fProcCounter->end(); it++ )
398  {
399  OneProcessCount* process = *it;
400  for ( G4int i = 0; i < process->GetCounter() ; i++)
401  this->CountProcesses(process->GetName());
402  }
403 
404  G4Run::Merge(run);
405 
406 }
407 
408 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
static const double prec
Definition: RanecuEngine.cc:58
DetectorConstruction * fDetector
static constexpr double MeV
Definition: G4SIunits.hh:214
G4int numberOfEvent
Definition: G4Run.hh:59
std::vector< OneProcessCount * > ProcessesCount
#define G4endl
Definition: G4ios.hh:61
G4double GetDEDX(G4double kinEnergy, const G4ParticleDefinition *, const G4Material *, const G4Region *r=nullptr)
const G4String & GetParticleName() const
const G4String & GetName() const
Definition: G4Material.hh:179
double G4double
Definition: G4Types.hh:76
static constexpr double cm2
Definition: G4SIunits.hh:120
static constexpr double mg
Definition: G4SIunits.hh:184
virtual void Merge(const G4Run *)
The primary generator action class with particle gun.
G4CsvAnalysisManager G4AnalysisManager
Definition: g4csv_defs.hh:77
static G4Electron * Electron()
Definition: G4Electron.cc:94
std::map< G4String, G4int > fProcCounter
Definition: G4Run.hh:46
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
int G4int
Definition: G4Types.hh:78
void CountProcesses(G4String procName)
G4double GetParticleEnergy() const
G4double GetCSDARange(G4double kinEnergy, const G4ParticleDefinition *, const G4Material *, const G4Region *r=nullptr)
G4GLOB_DLL std::ostream G4cout
static PROLOG_HANDLER error
Definition: xmlrole.cc:112
virtual void Merge(const G4Run *)
Definition: G4Run.cc:54
G4ParticleDefinition * GetParticleDefinition() const
static void showEngineStatus()
Definition: Random.cc:303
PrimaryGeneratorAction * fKinematic
Simple detector construction with a box volume placed in a world.
G4double GetDensity() const
Definition: G4Material.hh:181