| Geant4 Cross Reference |
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 //
26 /// \file electromagnetic/TestEm17/src/RunAction.cc
27 /// \brief Implementation of the RunAction class
28 //
29 //
30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32
33 #include "RunAction.hh"
34
35 #include "DetectorConstruction.hh"
36 #include "HistoManager.hh"
37 #include "MuCrossSections.hh"
38 #include "PrimaryGeneratorAction.hh"
39
40 #include "G4EmCalculator.hh"
41 #include "G4PhysicalConstants.hh"
42 #include "G4ProductionCutsTable.hh"
43 #include "G4Run.hh"
44 #include "G4RunManager.hh"
45 #include "G4SystemOfUnits.hh"
46 #include "G4UnitsTable.hh"
47 #include "Randomize.hh"
48
49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
50
51 RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* prim, HistoManager* HistM)
52 : G4UserRunAction(), fDetector(det), fPrimary(prim), fProcCounter(0), fHistoManager(HistM)
53 {
54 fMucs = new MuCrossSections();
55 }
56
57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
58
59 RunAction::~RunAction()
60 {
61 delete fMucs;
62 }
63
64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
65
66 void RunAction::BeginOfRunAction(const G4Run* aRun)
67 {
68 G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
69
70 // save Rndm status
71 CLHEP::HepRandom::showEngineStatus();
72
73 fProcCounter = new ProcessesCount();
74 fHistoManager->Book();
75 }
76
77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
78
79 void RunAction::CountProcesses(const G4String& procName)
80 {
81 // does the process already encounted ?
82 size_t n = fProcCounter->size();
83 for (size_t i = 0; i < n; ++i) {
84 if ((*fProcCounter)[i]->GetName() == procName) {
85 (*fProcCounter)[i]->Count();
86 return;
87 }
88 }
89 OneProcessCount* count = new OneProcessCount(procName);
90 count->Count();
91 fProcCounter->push_back(count);
92 }
93
94 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
95
96 void RunAction::EndOfRunAction(const G4Run* aRun)
97 {
98 G4int NbOfEvents = aRun->GetNumberOfEvent();
99 if (NbOfEvents == 0) return;
100
101 // std::ios::fmtflags mode = G4cout.flags();
102 G4int prec = G4cout.precision(2);
103
104 const G4Material* material = fDetector->GetMaterial();
105 G4double length = fDetector->GetSize();
106 G4double density = material->GetDensity();
107
108 G4String particle = fPrimary->GetParticleGun()->GetParticleDefinition()->GetParticleName();
109 G4double energy = fPrimary->GetParticleGun()->GetParticleEnergy();
110
111 G4cout << "\n The run consists of " << NbOfEvents << " " << particle << " of "
112 << G4BestUnit(energy, "Energy") << " through " << G4BestUnit(length, "Length") << " of "
113 << material->GetName() << " (density: " << G4BestUnit(density, "Volumic Mass") << ")"
114 << G4endl;
115
116 // total number of process calls
117 G4double countTot = 0.;
118 G4cout << "\n Number of process calls --->";
119 for (size_t i = 0; i < fProcCounter->size(); ++i) {
120 G4String procName = (*fProcCounter)[i]->GetName();
121 if (procName != "Transportation") {
122 G4int count = (*fProcCounter)[i]->GetCounter();
123 G4cout << "\t" << procName << " : " << count;
124 countTot += count;
125 }
126 }
127
128 // compute totalCrossSection, meanFreePath and massicCrossSection
129 //
130 G4double totalCrossSection = countTot / (NbOfEvents * length);
131 G4double MeanFreePath = 1. / totalCrossSection;
132 G4double massCrossSection = totalCrossSection / density;
133
134 G4cout.precision(5);
135 G4cout << "\n Simulation: "
136 << "total CrossSection = " << totalCrossSection * cm << " /cm"
137 << "\t MeanFreePath = " << G4BestUnit(MeanFreePath, "Length")
138 << "\t massicCrossSection = " << massCrossSection * g / cm2 << " cm2/g" << G4endl;
139
140 // compute theoretical predictions
141 //
142 if (particle == "mu+" || particle == "mu-") {
143 totalCrossSection = 0.;
144 for (size_t i = 0; i < fProcCounter->size(); ++i) {
145 G4String procName = (*fProcCounter)[i]->GetName();
146 if (procName != "Transportation") {
147 totalCrossSection += ComputeTheory(procName, NbOfEvents);
148 FillCrossSectionHisto(procName, NbOfEvents);
149 }
150 }
151
152 MeanFreePath = 1. / totalCrossSection;
153 massCrossSection = totalCrossSection / density;
154
155 G4cout << " Theory: "
156 << "total CrossSection = " << totalCrossSection * cm << " /cm"
157 << "\t MeanFreePath = " << G4BestUnit(MeanFreePath, "Length")
158 << "\t massicCrossSection = " << massCrossSection * g / cm2 << " cm2/g" << G4endl;
159 }
160
161 // G4cout.setf(mode,std::ios::floatfield);
162 G4cout.precision(prec);
163
164 // delete and remove all contents in fProcCounter
165 size_t n = fProcCounter->size();
166 for (size_t i = 0; i < n; ++i) {
167 delete (*fProcCounter)[i];
168 }
169 delete fProcCounter;
170
171 fHistoManager->Save();
172
173 // show Rndm status
174 // CLHEP::HepRandom::showEngineStatus();
175 }
176
177 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
178
179 G4double RunAction::ComputeTheory(const G4String& process, G4int NbOfMu)
180 {
181 const G4Material* material = fDetector->GetMaterial();
182 G4double ekin = fPrimary->GetParticleGun()->GetParticleEnergy();
183 G4double particleMass = fPrimary->GetParticleGun()->GetParticleDefinition()->GetPDGMass();
184
185 G4int id = 0;
186 G4double cut = 1.e-10 * ekin;
187 if (process == "muIoni") {
188 id = 11;
189 cut = GetEnergyCut(material, 1);
190 }
191 else if (process == "muPairProd") {
192 id = 12;
193 cut = 2 * (GetEnergyCut(material, 1) + electron_mass_c2);
194 }
195 else if (process == "muBrems") {
196 id = 13;
197 cut = GetEnergyCut(material, 0);
198 }
199 else if (process == "muonNuclear") {
200 id = 14;
201 cut = 100 * MeV;
202 }
203 else if (process == "muToMuonPairProd") {
204 id = 18;
205 cut = 2 * particleMass;
206 }
207 if (id == 0) {
208 return 0.;
209 }
210
211 G4int nbOfBins = 100;
212 // G4double binMin = -10.;
213 G4double binMin = std::log10(cut / ekin);
214 G4double binMax = 0.;
215 G4double binWidth = (binMax - binMin) / G4double(nbOfBins);
216
217 // create histo for theoretical crossSections, with same bining as simulation
218 //
219 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
220
221 G4H1* histoTh = 0;
222 if (fHistoManager->HistoExist(id)) {
223 histoTh = analysisManager->GetH1(fHistoManager->GetHistoID(id));
224 nbOfBins = fHistoManager->GetNbins(id);
225 binMin = fHistoManager->GetVmin(id);
226 binMax = fHistoManager->GetVmax(id);
227 binWidth = fHistoManager->GetBinWidth(id);
228 }
229
230 // compute and plot differential crossSection, as function of energy transfert.
231 // compute and return integrated crossSection for a given process.
232 //(note: to compare with simulation, the integrated crossSection is function
233 // of the energy cut.)
234 //
235 G4double lgeps, etransf, sigmaE, dsigma;
236 G4double sigmaTot = 0.;
237 const G4double ln10 = std::log(10.);
238 G4double length = fDetector->GetSize();
239
240 // G4cout << "MU: " << process << " E= " << ekin
241 // <<" binMin= " << binMin << " binW= " << binWidth << G4endl;
242
243 for (G4int ibin = 0; ibin < nbOfBins; ibin++) {
244 lgeps = binMin + (ibin + 0.5) * binWidth;
245 etransf = ekin * std::pow(10., lgeps);
246 sigmaE = fMucs->CR_Macroscopic(process, material, ekin, etransf);
247 dsigma = sigmaE * etransf * binWidth * ln10;
248 if (etransf > cut) sigmaTot += dsigma;
249 if (histoTh) {
250 G4double NbProcess = NbOfMu * length * dsigma;
251 histoTh->fill(lgeps, NbProcess);
252 }
253 }
254
255 // return integrated crossSection
256 //
257 return sigmaTot;
258 }
259
260 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
261
262 void RunAction::FillCrossSectionHisto(const G4String& process, G4int)
263 {
264 const G4Material* material = fDetector->GetMaterial();
265 G4double ekin = fPrimary->GetParticleGun()->GetParticleEnergy();
266 G4ParticleDefinition* particle = fPrimary->GetParticleGun()->GetParticleDefinition();
267 G4double particleMass = particle->GetPDGMass();
268
269 G4EmCalculator emCal;
270
271 G4int id = 0;
272 G4double cut = 1.e-10 * ekin;
273 if (process == "muIoni") {
274 id = 21;
275 cut = GetEnergyCut(material, 1);
276 }
277 else if (process == "muPairProd") {
278 id = 22;
279 cut = 2 * (GetEnergyCut(material, 1) + electron_mass_c2);
280 }
281 else if (process == "muBrems") {
282 id = 23;
283 cut = GetEnergyCut(material, 0);
284 }
285 else if (process == "muonNuclear") {
286 id = 24;
287 cut = 100 * MeV;
288 }
289 else if (process == "muToMuonPairProd") {
290 id = 28;
291 cut = 2 * particleMass;
292 }
293 if (id == 0) {
294 return;
295 }
296
297 G4int nbOfBins = 100;
298 G4double binMin = cut;
299 G4double binMax = ekin;
300 G4double binWidth = (binMax - binMin) / G4double(nbOfBins);
301
302 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
303
304 G4H1* histoTh = 0;
305 if (fHistoManager->HistoExist(id)) {
306 histoTh = analysisManager->GetH1(fHistoManager->GetHistoID(id));
307 nbOfBins = fHistoManager->GetNbins(id);
308 binMin = fHistoManager->GetVmin(id);
309 binMax = fHistoManager->GetVmax(id);
310 binWidth = fHistoManager->GetBinWidth(id);
311 }
312
313 G4double sigma, primaryEnergy;
314
315 for (G4int ibin = 0; ibin < nbOfBins; ibin++) {
316 primaryEnergy = binMin + (ibin + 0.5) * binWidth;
317 sigma = emCal.GetCrossSectionPerVolume(primaryEnergy, particle, process, material);
318 if (histoTh) {
319 histoTh->fill(primaryEnergy, sigma);
320 }
321 }
322 }
323
324 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
325
326 G4double RunAction::GetEnergyCut(const G4Material* material, G4int idParticle)
327 {
328 G4ProductionCutsTable* table = G4ProductionCutsTable::GetProductionCutsTable();
329
330 size_t index = 0;
331 while ((table->GetMaterialCutsCouple(index)->GetMaterial() != material)
332 && (index < table->GetTableSize()))
333 index++;
334
335 return (*(table->GetEnergyCutsVector(idParticle)))[index];
336 }
337
338 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
339