| 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 electromagnetic/TestEm17/src/RunActi << 26 // $Id: RunAction.cc,v 1.4 2008/03/31 10:22:59 vnivanch Exp $
27 /// \brief Implementation of the RunAction cla << 27 // GEANT4 tag $Name: geant4-09-03-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 33
35 #include "DetectorConstruction.hh" 34 #include "DetectorConstruction.hh"
>> 35 #include "PrimaryGeneratorAction.hh"
36 #include "HistoManager.hh" 36 #include "HistoManager.hh"
37 #include "MuCrossSections.hh" 37 #include "MuCrossSections.hh"
38 #include "PrimaryGeneratorAction.hh" <<
39 38
40 #include "G4EmCalculator.hh" <<
41 #include "G4PhysicalConstants.hh" <<
42 #include "G4ProductionCutsTable.hh" <<
43 #include "G4Run.hh" 39 #include "G4Run.hh"
44 #include "G4RunManager.hh" 40 #include "G4RunManager.hh"
45 #include "G4SystemOfUnits.hh" <<
46 #include "G4UnitsTable.hh" 41 #include "G4UnitsTable.hh"
>> 42
47 #include "Randomize.hh" 43 #include "Randomize.hh"
48 44
>> 45 #ifdef G4ANALYSIS_USE
>> 46 #include "AIDA/AIDA.h"
>> 47 #endif
>> 48
49 //....oooOO0OOooo........oooOO0OOooo........oo 49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
50 50
51 RunAction::RunAction(DetectorConstruction* det << 51 RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* prim,
52 : G4UserRunAction(), fDetector(det), fPrimar << 52 HistoManager* HistM)
53 { << 53 : detector(det), primary(prim), ProcCounter(0), histoManager(HistM)
54 fMucs = new MuCrossSections(); << 54 {}
55 } <<
56 55
57 //....oooOO0OOooo........oooOO0OOooo........oo 56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
58 57
59 RunAction::~RunAction() 58 RunAction::~RunAction()
60 { << 59 {}
61 delete fMucs; <<
62 } <<
63 60
64 //....oooOO0OOooo........oooOO0OOooo........oo 61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
65 62
66 void RunAction::BeginOfRunAction(const G4Run* 63 void RunAction::BeginOfRunAction(const G4Run* aRun)
67 { << 64 {
68 G4cout << "### Run " << aRun->GetRunID() << 65 G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
69 << 66
70 // save Rndm status 67 // save Rndm status
>> 68 G4RunManager::GetRunManager()->SetRandomNumberStore(true);
71 CLHEP::HepRandom::showEngineStatus(); 69 CLHEP::HepRandom::showEngineStatus();
72 70
73 fProcCounter = new ProcessesCount(); << 71 ProcCounter = new ProcessesCount;
74 fHistoManager->Book(); << 72
>> 73 histoManager->book();
75 } 74 }
76 75
77 //....oooOO0OOooo........oooOO0OOooo........oo 76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
78 77
79 void RunAction::CountProcesses(const G4String& << 78 void RunAction::CountProcesses(G4String procName)
80 { 79 {
81 // does the process already encounted ? << 80 //does the process already encounted ?
82 size_t n = fProcCounter->size(); << 81 size_t nbProc = ProcCounter->size();
83 for (size_t i = 0; i < n; ++i) { << 82 size_t i = 0;
84 if ((*fProcCounter)[i]->GetName() == procN << 83 while ((i<nbProc)&&((*ProcCounter)[i]->GetName()!=procName)) i++;
85 (*fProcCounter)[i]->Count(); << 84 if (i == nbProc) ProcCounter->push_back( new OneProcessCount(procName));
86 return; << 85
87 } << 86 (*ProcCounter)[i]->Count();
88 } <<
89 OneProcessCount* count = new OneProcessCount <<
90 count->Count(); <<
91 fProcCounter->push_back(count); <<
92 } 87 }
93 88
94 //....oooOO0OOooo........oooOO0OOooo........oo 89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
95 90
96 void RunAction::EndOfRunAction(const G4Run* aR 91 void RunAction::EndOfRunAction(const G4Run* aRun)
97 { 92 {
98 G4int NbOfEvents = aRun->GetNumberOfEvent(); 93 G4int NbOfEvents = aRun->GetNumberOfEvent();
99 if (NbOfEvents == 0) return; 94 if (NbOfEvents == 0) return;
100 << 95
101 // std::ios::fmtflags mode = G4cout.flags() << 96 std::ios::fmtflags mode = G4cout.flags();
102 G4int prec = G4cout.precision(2); << 97 G4int prec = G4cout.precision(2);
103 << 98
104 const G4Material* material = fDetector->GetM << 99 G4Material* material = detector->GetMaterial();
105 G4double length = fDetector->GetSize(); << 100 G4double length = detector->GetSize();
106 G4double density = material->GetDensity(); 101 G4double density = material->GetDensity();
107 << 102
108 G4String particle = fPrimary->GetParticleGun << 103 G4String particle = primary->GetParticleGun()->GetParticleDefinition()
109 G4double energy = fPrimary->GetParticleGun() << 104 ->GetParticleName();
110 << 105 G4double energy = primary->GetParticleGun()->GetParticleEnergy();
111 G4cout << "\n The run consists of " << NbOfE << 106
112 << G4BestUnit(energy, "Energy") << " << 107 G4cout << "\n The run consists of " << NbOfEvents << " "<< particle << " of "
113 << material->GetName() << " (density: << 108 << G4BestUnit(energy,"Energy") << " through "
114 << G4endl; << 109 << G4BestUnit(length,"Length") << " of "
115 << 110 << material->GetName() << " (density: "
116 // total number of process calls << 111 << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
>> 112
>> 113 //total number of process calls
117 G4double countTot = 0.; 114 G4double countTot = 0.;
118 G4cout << "\n Number of process calls --->"; 115 G4cout << "\n Number of process calls --->";
119 for (size_t i = 0; i < fProcCounter->size(); << 116 for (size_t i=0; i< ProcCounter->size();i++) {
120 G4String procName = (*fProcCounter)[i]->Ge << 117 G4String procName = (*ProcCounter)[i]->GetName();
121 if (procName != "Transportation") { << 118 if (procName != "Transportation") {
122 G4int count = (*fProcCounter)[i]->GetCou << 119 G4int count = (*ProcCounter)[i]->GetCounter();
123 G4cout << "\t" << procName << " : " << c << 120 G4cout << "\t" << procName << " : " << count;
124 countTot += count; << 121 countTot += count;
125 } << 122 }
126 } << 123 }
127 << 124 G4cout << G4endl;
128 // compute totalCrossSection, meanFreePath a << 125
>> 126 //compute totalCrossSection, meanFreePath and massicCrossSection
129 // 127 //
130 G4double totalCrossSection = countTot / (NbO << 128 G4double totalCrossSection = countTot/(NbOfEvents*length);
131 G4double MeanFreePath = 1. / totalCrossSecti << 129 G4double MeanFreePath = 1./totalCrossSection;
132 G4double massCrossSection = totalCrossSectio << 130 G4double massCrossSection =totalCrossSection/density;
133 << 131
134 G4cout.precision(5); 132 G4cout.precision(5);
135 G4cout << "\n Simulation: " 133 G4cout << "\n Simulation: "
136 << "total CrossSection = " << totalCr << 134 << "total CrossSection = " << totalCrossSection*cm << " /cm"
137 << "\t MeanFreePath = " << G4BestUnit << 135 << "\t MeanFreePath = " << G4BestUnit(MeanFreePath,"Length")
138 << "\t massicCrossSection = " << mass << 136 << "\t massicCrossSection = " << massCrossSection*g/cm2 << " cm2/g"
139 << 137 << G4endl;
140 // compute theoretical predictions << 138
>> 139 //compute theoritical predictions
141 // 140 //
142 if (particle == "mu+" || particle == "mu-") << 141 if(particle == "mu+" || particle == "mu-") {
143 totalCrossSection = 0.; 142 totalCrossSection = 0.;
144 for (size_t i = 0; i < fProcCounter->size( << 143 for (size_t i=0; i< ProcCounter->size();i++) {
145 G4String procName = (*fProcCounter)[i]-> << 144 G4String procName = (*ProcCounter)[i]->GetName();
146 if (procName != "Transportation") { << 145 if (procName != "Transportation")
147 totalCrossSection += ComputeTheory(pro << 146 totalCrossSection += ComputeTheory(procName, NbOfEvents);
148 FillCrossSectionHisto(procName, NbOfEv <<
149 } <<
150 } 147 }
151 << 148
152 MeanFreePath = 1. / totalCrossSection; << 149 MeanFreePath = 1./totalCrossSection;
153 massCrossSection = totalCrossSection / den << 150 massCrossSection = totalCrossSection/density;
154 << 151
155 G4cout << " Theory: " 152 G4cout << " Theory: "
156 << "total CrossSection = " << total << 153 << "total CrossSection = " << totalCrossSection*cm << " /cm"
157 << "\t MeanFreePath = " << G4BestUn << 154 << "\t MeanFreePath = " << G4BestUnit(MeanFreePath,"Length")
158 << "\t massicCrossSection = " << ma << 155 << "\t massicCrossSection = " << massCrossSection*g/cm2 << " cm2/g"
159 } << 156 << G4endl;
160 << 157 }
161 // G4cout.setf(mode,std::ios::floatfield); << 158
162 G4cout.precision(prec); << 159 G4cout.setf(mode,std::ios::floatfield);
163 << 160 G4cout.precision(prec);
164 // delete and remove all contents in fProcCo << 161
165 size_t n = fProcCounter->size(); << 162 // delete and remove all contents in ProcCounter
166 for (size_t i = 0; i < n; ++i) { << 163 while (ProcCounter->size()>0){
167 delete (*fProcCounter)[i]; << 164 OneProcessCount* aProcCount=ProcCounter->back();
168 } << 165 ProcCounter->pop_back();
169 delete fProcCounter; << 166 delete aProcCount;
170 << 167 }
171 fHistoManager->Save(); << 168 delete ProcCounter;
172 << 169
>> 170 histoManager->save();
>> 171
173 // show Rndm status 172 // show Rndm status
174 // CLHEP::HepRandom::showEngineStatus(); << 173 CLHEP::HepRandom::showEngineStatus();
175 } 174 }
176 175
177 //....oooOO0OOooo........oooOO0OOooo........oo 176 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
178 177
179 G4double RunAction::ComputeTheory(const G4Stri << 178 G4double RunAction::ComputeTheory(G4String process, G4int NbOfMu)
180 { << 179 {
181 const G4Material* material = fDetector->GetM << 180 G4Material* material = detector->GetMaterial();
182 G4double ekin = fPrimary->GetParticleGun()-> << 181 G4double length = detector->GetSize();
183 G4double particleMass = fPrimary->GetParticl << 182 G4double ekin = primary->GetParticleGun()->GetParticleEnergy();
184 << 183 MuCrossSections crossSections;
185 G4int id = 0; << 184
186 G4double cut = 1.e-10 * ekin; << 185 G4int id = 0; G4double cut = 0.;
187 if (process == "muIoni") { << 186 if (process == "muIoni") {id = 1; cut = GetEnergyCut(material,1);}
188 id = 11; << 187 else if (process == "muPairProd") {id = 2; cut = 2*(GetEnergyCut(material,1)
189 cut = GetEnergyCut(material, 1); << 188 + electron_mass_c2); }
190 } << 189 else if (process == "muBrems") {id = 3; cut = GetEnergyCut(material,0);}
191 else if (process == "muPairProd") { << 190 else if (process == "muNucl") id = 4;
192 id = 12; << 191 else if (process == "hIoni") {id = 5; cut = GetEnergyCut(material,1);}
193 cut = 2 * (GetEnergyCut(material, 1) + ele << 192 else if (process == "hPairProd") {id = 6; cut = 2*(GetEnergyCut(material,1)
194 } << 193 + electron_mass_c2); }
195 else if (process == "muBrems") { << 194 else if (process == "hBrems") {id = 7; cut = GetEnergyCut(material,0);}
196 id = 13; << 195 if (id == 0) return 0.;
197 cut = GetEnergyCut(material, 0); << 196
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; 197 G4int nbOfBins = 100;
212 // G4double binMin = -10.; << 198 G4double binMin = -10., binMax = 0., binWidth = (binMax-binMin)/nbOfBins;
213 G4double binMin = std::log10(cut / ekin); << 199
214 G4double binMax = 0.; << 200 #ifdef G4ANALYSIS_USE
215 G4double binWidth = (binMax - binMin) / G4do << 201 //create histo for theoritical crossSections, with same bining as simulation
216 <<
217 // create histo for theoretical crossSection <<
218 // 202 //
219 G4AnalysisManager* analysisManager = G4Analy << 203 const G4String label[] = { "0", "1", "2", "3", "4", "5", "6", "7", "8", "9",
220 << 204 "10", "11", "12", "13", "14", "15", "16", "17", "18", "19"};
221 G4H1* histoTh = 0; << 205
222 if (fHistoManager->HistoExist(id)) { << 206 AIDA::IHistogram1D* histoMC = 0; AIDA::IHistogram1D* histoTh = 0;
223 histoTh = analysisManager->GetH1(fHistoMan << 207 if (histoManager->HistoExist(id)) {
224 nbOfBins = fHistoManager->GetNbins(id); << 208 histoMC = histoManager->GetHisto(id);
225 binMin = fHistoManager->GetVmin(id); << 209 nbOfBins = histoManager->GetNbins(id);
226 binMax = fHistoManager->GetVmax(id); << 210 binMin = histoManager->GetVmin (id);
227 binWidth = fHistoManager->GetBinWidth(id); << 211 binMax = histoManager->GetVmax (id);
228 } << 212 binWidth = histoManager->GetBinWidth(id);
229 << 213
230 // compute and plot differential crossSectio << 214 G4String labelTh = label[MaxHisto + id];
231 // compute and return integrated crossSectio << 215 G4String titleTh = histoManager->GetTitle(id) + " (Th)";
>> 216 histoTh = histoManager->GetHistogramFactory()
>> 217 ->createHistogram1D(labelTh,titleTh,nbOfBins,binMin,binMax);
>> 218 }
>> 219 #endif
>> 220
>> 221 //compute and plot differential crossSection, as function of energy transfert.
>> 222 //compute and return integrated crossSection for a given process.
232 //(note: to compare with simulation, the int 223 //(note: to compare with simulation, the integrated crossSection is function
233 // of the energy cut.) << 224 // of the energy cut.)
234 // << 225 //
235 G4double lgeps, etransf, sigmaE, dsigma; << 226 G4double lgeps, etransf, sigmaE, dsigma, NbProcess;
236 G4double sigmaTot = 0.; 227 G4double sigmaTot = 0.;
237 const G4double ln10 = std::log(10.); 228 const G4double ln10 = std::log(10.);
238 G4double length = fDetector->GetSize(); << 229
239 << 230 for (G4int ibin=0; ibin<nbOfBins; ibin++) {
240 // G4cout << "MU: " << process << " E= " << << 231 lgeps = binMin + (ibin+0.5)*binWidth;
241 // <<" binMin= " << binMin << " binW << 232 etransf = ekin*std::pow(10.,lgeps);
242 << 233 sigmaE = crossSections.CR_Macroscopic(process,material,ekin,etransf);
243 for (G4int ibin = 0; ibin < nbOfBins; ibin++ << 234 dsigma = sigmaE*etransf*binWidth*ln10;
244 lgeps = binMin + (ibin + 0.5) * binWidth; << 235 if (etransf > cut) sigmaTot += dsigma;
245 etransf = ekin * std::pow(10., lgeps); << 236 NbProcess = NbOfMu*length*dsigma;
246 sigmaE = fMucs->CR_Macroscopic(process, ma << 237 #ifdef G4ANALYSIS_USE
247 dsigma = sigmaE * etransf * binWidth * ln1 << 238 if (histoTh) histoTh->fill(lgeps,NbProcess);
248 if (etransf > cut) sigmaTot += dsigma; << 239 #endif
249 if (histoTh) { << 240 }
250 G4double NbProcess = NbOfMu * length * d << 241
251 histoTh->fill(lgeps, NbProcess); << 242 #ifdef G4ANALYSIS_USE
252 } << 243 //compare simulation and theory
253 } <<
254 <<
255 // return integrated crossSection <<
256 // 244 //
257 return sigmaTot; << 245 if (histoMC && histoTh) histoManager->GetHistogramFactory()
258 } << 246 ->divide(label[2*MaxHisto+id], *histoMC, *histoTh);
259 << 247 #endif
260 //....oooOO0OOooo........oooOO0OOooo........oo << 248
261 << 249 //return integrated crossSection
262 void RunAction::FillCrossSectionHisto(const G4 << 250 //
263 { << 251 return sigmaTot;
264 const G4Material* material = fDetector->GetM <<
265 G4double ekin = fPrimary->GetParticleGun()-> <<
266 G4ParticleDefinition* particle = fPrimary->G <<
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) + ele <<
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) / G4do <<
301 <<
302 G4AnalysisManager* analysisManager = G4Analy <<
303 <<
304 G4H1* histoTh = 0; <<
305 if (fHistoManager->HistoExist(id)) { <<
306 histoTh = analysisManager->GetH1(fHistoMan <<
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) * bi <<
317 sigma = emCal.GetCrossSectionPerVolume(pri <<
318 if (histoTh) { <<
319 histoTh->fill(primaryEnergy, sigma); <<
320 } <<
321 } <<
322 } 252 }
323 253
324 //....oooOO0OOooo........oooOO0OOooo........oo 254 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
325 255
326 G4double RunAction::GetEnergyCut(const G4Mater << 256 #include "G4ProductionCutsTable.hh"
327 { <<
328 G4ProductionCutsTable* table = G4ProductionC <<
329 257
330 size_t index = 0; << 258 G4double RunAction::GetEnergyCut(G4Material* material, G4int idParticle)
331 while ((table->GetMaterialCutsCouple(index)- << 259 {
332 && (index < table->GetTableSize())) << 260 G4ProductionCutsTable* table = G4ProductionCutsTable::GetProductionCutsTable();
333 index++; << 261
>> 262 size_t index = 0;
>> 263 while ( (table->GetMaterialCutsCouple(index)->GetMaterial() != material) &&
>> 264 (index < table->GetTableSize())) index++;
334 265
335 return (*(table->GetEnergyCutsVector(idParti << 266 return (*(table->GetEnergyCutsVector(idParticle)))[index];
336 } << 267 }
337 268
338 //....oooOO0OOooo........oooOO0OOooo........oo 269 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 270
339 271