Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
examples/extended/electromagnetic/TestEm17/src/RunAction.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: RunAction.cc 99671 2016-09-30 08:05:17Z gcosmo $
30 //
31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
33 
34 #include "RunAction.hh"
35 
36 #include "DetectorConstruction.hh"
37 #include "PrimaryGeneratorAction.hh"
38 #include "HistoManager.hh"
39 #include "MuCrossSections.hh"
40 
41 #include "G4Run.hh"
42 #include "G4RunManager.hh"
43 #include "G4UnitsTable.hh"
44 
45 #include "G4PhysicalConstants.hh"
46 #include "G4SystemOfUnits.hh"
47 #include "Randomize.hh"
48 
49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
50 
52  HistoManager* HistM)
53  : G4UserRunAction(),
54  fDetector(det), fPrimary(prim), fProcCounter(0), fHistoManager(HistM)
55 {}
56 
57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
58 
60 {}
61 
62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
63 
64 void RunAction::BeginOfRunAction(const G4Run* aRun)
65 {
66  G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
67 
68  // save Rndm status
71 
74 }
75 
76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
77 
79 {
80  //does the process already encounted ?
81  size_t nbProc = fProcCounter->size();
82  size_t i = 0;
83  while ((i<nbProc)&&((*fProcCounter)[i]->GetName()!=procName)) i++;
84  if (i == nbProc) fProcCounter->push_back( new OneProcessCount(procName));
85 
86  (*fProcCounter)[i]->Count();
87 }
88 
89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
90 
91 void RunAction::EndOfRunAction(const G4Run* aRun)
92 {
93  G4int NbOfEvents = aRun->GetNumberOfEvent();
94  if (NbOfEvents == 0) return;
95 
96  std::ios::fmtflags mode = G4cout.flags();
97  G4int prec = G4cout.precision(2);
98 
99  G4Material* material = fDetector->GetMaterial();
100  G4double length = fDetector->GetSize();
101  G4double density = material->GetDensity();
102 
104  ->GetParticleName();
106 
107  G4cout << "\n The run consists of " << NbOfEvents << " "<< particle << " of "
108  << G4BestUnit(energy,"Energy") << " through "
109  << G4BestUnit(length,"Length") << " of "
110  << material->GetName() << " (density: "
111  << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
112 
113  //total number of process calls
114  G4double countTot = 0.;
115  G4cout << "\n Number of process calls --->";
116  for (size_t i=0; i< fProcCounter->size();i++) {
117  G4String procName = (*fProcCounter)[i]->GetName();
118  if (procName != "Transportation") {
119  G4int count = (*fProcCounter)[i]->GetCounter();
120  G4cout << "\t" << procName << " : " << count;
121  countTot += count;
122  }
123  }
124  G4cout << G4endl;
125 
126  //compute totalCrossSection, meanFreePath and massicCrossSection
127  //
128  G4double totalCrossSection = countTot/(NbOfEvents*length);
129  G4double MeanFreePath = 1./totalCrossSection;
130  G4double massCrossSection =totalCrossSection/density;
131 
132  G4cout.precision(5);
133  G4cout << "\n Simulation: "
134  << "total CrossSection = " << totalCrossSection*cm << " /cm"
135  << "\t MeanFreePath = " << G4BestUnit(MeanFreePath,"Length")
136  << "\t massicCrossSection = " << massCrossSection*g/cm2 << " cm2/g"
137  << G4endl;
138 
139  //compute theoritical predictions
140  //
141  if(particle == "mu+" || particle == "mu-") {
142  totalCrossSection = 0.;
143  for (size_t i=0; i< fProcCounter->size();i++) {
144  G4String procName = (*fProcCounter)[i]->GetName();
145  if (procName != "Transportation") {
146  totalCrossSection += ComputeTheory(procName, NbOfEvents);
147  }
148  }
149 
150  MeanFreePath = 1./totalCrossSection;
151  massCrossSection = totalCrossSection/density;
152 
153  G4cout << " Theory: "
154  << "total CrossSection = " << totalCrossSection*cm << " /cm"
155  << "\t MeanFreePath = " << G4BestUnit(MeanFreePath,"Length")
156  << "\t massicCrossSection = " << massCrossSection*g/cm2 << " cm2/g"
157  << G4endl;
158  }
159 
160  G4cout.setf(mode,std::ios::floatfield);
161  G4cout.precision(prec);
162 
163  // delete and remove all contents in fProcCounter
164  while (fProcCounter->size()>0){
165  OneProcessCount* aProcCount=fProcCounter->back();
166  fProcCounter->pop_back();
167  delete aProcCount;
168  }
169  delete fProcCounter;
170 
171  fHistoManager->Save();
172 
173  // show Rndm status
175 }
176 
177 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
178 
180 {
181  G4Material* material = fDetector->GetMaterial();
183  MuCrossSections crossSections;
184 
185  G4int id = 0; G4double cut = 0.;
186  if (process == "muIoni") {id = 11; cut = GetEnergyCut(material,1);}
187  else if (process == "muPairProd") {id = 12; cut = 2*(GetEnergyCut(material,1)
188  + electron_mass_c2); }
189  else if (process == "muBrems") {id = 13; cut = GetEnergyCut(material,0);}
190  else if (process == "muonNuclear"){id = 14; }
191  if (id == 0) { return 0.; }
192 
193  G4int nbOfBins = 100;
194  G4double binMin = -10.;
195  G4double binMax = 0.;
196  G4double binWidth = (binMax-binMin)/G4double(nbOfBins);
197 
198  //create histo for theoritical crossSections, with same bining as simulation
199  //
200  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
201 
202  G4AnaH1* histoTh = 0;
203  // G4AnaH1* histoMC = 0;
204  if (fHistoManager->HistoExist(id)) {
205  histoTh = analysisManager->GetH1(fHistoManager->GetHistoID(id));
206  //histoMC = analysisManager->GetH1(fHistoManager->GetHistoID(id-10));
207  nbOfBins = fHistoManager->GetNbins(id);
208  binMin = fHistoManager->GetVmin (id);
209  binMax = fHistoManager->GetVmax (id);
210  binWidth = fHistoManager->GetBinWidth(id);
211  }
212 
213  //compute and plot differential crossSection, as function of energy transfert.
214  //compute and return integrated crossSection for a given process.
215  //(note: to compare with simulation, the integrated crossSection is function
216  // of the energy cut.)
217  //
218  G4double lgeps, etransf, sigmaE, dsigma;
219  G4double sigmaTot = 0.;
220  const G4double ln10 = std::log(10.);
221  G4double length = fDetector->GetSize();
222 
223  for (G4int ibin=0; ibin<nbOfBins; ibin++) {
224  lgeps = binMin + (ibin+0.5)*binWidth;
225  etransf = ekin*std::pow(10.,lgeps);
226  sigmaE = crossSections.CR_Macroscopic(process,material,ekin,etransf);
227  dsigma = sigmaE*etransf*binWidth*ln10;
228  if (etransf > cut) sigmaTot += dsigma;
229  if (histoTh) {
230  G4double NbProcess = NbOfMu*length*dsigma;
231  histoTh->fill(lgeps, NbProcess);
232  }
233  }
234 
235  //return integrated crossSection
236  //
237  return sigmaTot;
238 }
239 
240 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
241 
242 #include "G4ProductionCutsTable.hh"
243 
245 {
247 
248  size_t index = 0;
249  while ( (table->GetMaterialCutsCouple(index)->GetMaterial() != material) &&
250  (index < table->GetTableSize())) index++;
251 
252  return (*(table->GetEnergyCutsVector(idParticle)))[index];
253 }
254 
255 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
256 
static G4RunManager * GetRunManager()
Definition: G4RunManager.cc:80
static const double prec
Definition: RanecuEngine.cc:58
std::vector< OneProcessCount * > ProcessesCount
#define G4endl
Definition: G4ios.hh:61
tools::histo::h1d G4AnaH1
Definition: g4csv_defs.hh:45
const G4String & GetParticleName() const
const std::vector< G4double > * GetEnergyCutsVector(size_t pcIdx) const
static const G4double ln10
G4int GetNumberOfEvent() const
Definition: G4Run.hh:79
const G4String & GetName() const
Definition: G4Material.hh:179
static constexpr double g
Definition: G4SIunits.hh:183
double G4double
Definition: G4Types.hh:76
static constexpr double cm2
Definition: G4SIunits.hh:120
static constexpr double electron_mass_c2
G4int GetRunID() const
Definition: G4Run.hh:76
double energy
Definition: plottest35.C:25
The primary generator action class with particle gun.
G4CsvAnalysisManager G4AnalysisManager
Definition: g4csv_defs.hh:77
Definition: G4Run.hh:46
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
static G4ProductionCutsTable * GetProductionCutsTable()
void SetRandomNumberStore(G4bool flag)
int G4int
Definition: G4Types.hh:78
G4double GetParticleEnergy() const
static constexpr double cm
Definition: G4SIunits.hh:119
G4GLOB_DLL std::ostream G4cout
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
const G4Material * GetMaterial() const
Definition of the MuCrossSections class.
G4ParticleDefinition * GetParticleDefinition() const
static void showEngineStatus()
Definition: Random.cc:303
Simple detector construction with a box volume placed in a world.
G4double CR_Macroscopic(const G4String &, G4Material *, G4double, G4double)
G4double GetDensity() const
Definition: G4Material.hh:181