| 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/TestEm1/src/Run.cc 26 /// \file electromagnetic/TestEm1/src/Run.cc
27 /// \brief Implementation of the Run class 27 /// \brief Implementation of the Run class
28 // 28 //
29 // << 29 // $Id: Run.cc 71376 2013-06-14 07:44:50Z maire $
>> 30 //
30 //....oooOO0OOooo........oooOO0OOooo........oo 31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31 //....oooOO0OOooo........oooOO0OOooo........oo 32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 33
33 #include "Run.hh" 34 #include "Run.hh"
34 <<
35 #include "DetectorConstruction.hh" 35 #include "DetectorConstruction.hh"
36 #include "PrimaryGeneratorAction.hh" 36 #include "PrimaryGeneratorAction.hh"
37 37
38 #include "G4EmCalculator.hh" <<
39 #include "G4SystemOfUnits.hh" <<
40 #include "G4UnitsTable.hh" 38 #include "G4UnitsTable.hh"
>> 39 #include "G4SystemOfUnits.hh"
>> 40 #include "G4EmCalculator.hh"
41 41
42 #include <iomanip> 42 #include <iomanip>
43 43
44 //....oooOO0OOooo........oooOO0OOooo........oo 44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
45 45
46 Run::Run(const DetectorConstruction* det) : fD << 46 Run::Run(DetectorConstruction* det)
>> 47 : G4Run(),
>> 48 fDetector(det),
>> 49 fParticle(0), fEkin(0.),
>> 50 fNbOfTraks0(0), fNbOfTraks1(0),
>> 51 fNbOfSteps0(0), fNbOfSteps1(0),
>> 52 fEdep(0),
>> 53 fTrueRange(0.), fTrueRange2(0.),
>> 54 fProjRange(0.), fProjRange2(0.),
>> 55 fTransvDev(0.), fTransvDev2(0.)
>> 56 { }
47 57
48 //....oooOO0OOooo........oooOO0OOooo........oo 58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
49 59
50 void Run::SetPrimary(const G4ParticleDefinitio << 60 Run::~Run()
51 { << 61 { }
>> 62
>> 63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 64
>> 65 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
>> 66 {
52 fParticle = particle; 67 fParticle = particle;
53 fEkin = energy; 68 fEkin = energy;
54 } 69 }
55 //....oooOO0OOooo........oooOO0OOooo........oo 70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56 71
57 void Run::CountProcesses(const G4String& procN << 72 void Run::CountProcesses(G4String procName)
58 { 73 {
59 std::map<G4String, G4int>::iterator it = fPr << 74 std::map<G4String,G4int>::iterator it = fProcCounter.find(procName);
60 if (it == fProcCounter.end()) { << 75 if ( it == fProcCounter.end()) {
61 fProcCounter[procName] = 1; 76 fProcCounter[procName] = 1;
62 } 77 }
63 else { 78 else {
64 fProcCounter[procName]++; << 79 fProcCounter[procName]++;
65 } 80 }
66 } 81 }
67 << 82
68 //....oooOO0OOooo........oooOO0OOooo........oo 83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
69 84
70 void Run::Merge(const G4Run* run) 85 void Run::Merge(const G4Run* run)
71 { 86 {
72 const Run* localRun = static_cast<const Run* 87 const Run* localRun = static_cast<const Run*>(run);
73 88
74 // pass information about primary particle 89 // pass information about primary particle
75 fParticle = localRun->fParticle; 90 fParticle = localRun->fParticle;
76 fEkin = localRun->fEkin; << 91 fEkin = localRun->fEkin;
77 92
78 // accumulate sums 93 // accumulate sums
79 // 94 //
80 fNbOfTraks0 += localRun->fNbOfTraks0; << 95 fNbOfTraks0 += localRun->fNbOfTraks0;
81 fNbOfTraks1 += localRun->fNbOfTraks1; << 96 fNbOfTraks1 += localRun->fNbOfTraks1;
82 fNbOfSteps0 += localRun->fNbOfSteps0; 97 fNbOfSteps0 += localRun->fNbOfSteps0;
83 fNbOfSteps1 += localRun->fNbOfSteps1; << 98 fNbOfSteps1 += localRun->fNbOfSteps1;
84 fEdep += localRun->fEdep; << 99 fEdep += localRun->fEdep;
85 fEleak += localRun->fEleak; << 100 fTrueRange += localRun->fTrueRange;
86 fNIEL += localRun->fNIEL; <<
87 fTrueRange += localRun->fTrueRange; <<
88 fTrueRange2 += localRun->fTrueRange2; 101 fTrueRange2 += localRun->fTrueRange2;
89 fProjRange += localRun->fProjRange; << 102 fProjRange += localRun->fProjRange;
90 fProjRange2 += localRun->fProjRange2; 103 fProjRange2 += localRun->fProjRange2;
91 fTransvDev += localRun->fTransvDev; << 104 fTransvDev += localRun->fTransvDev;
92 fTransvDev2 += localRun->fTransvDev2; << 105 fTransvDev2 += localRun->fTransvDev2;
93 << 106
94 // map: processes count << 107 //map: processes count
95 std::map<G4String, G4int>::const_iterator it << 108 std::map<G4String,G4int>::const_iterator it;
96 for (it = localRun->fProcCounter.begin(); it << 109 for (it = localRun->fProcCounter.begin();
>> 110 it !=localRun->fProcCounter.end(); ++it) {
>> 111
97 G4String procName = it->first; 112 G4String procName = it->first;
98 G4int localCount = it->second; << 113 G4int localCount = it->second;
99 if (fProcCounter.find(procName) == fProcCo << 114 if ( fProcCounter.find(procName) == fProcCounter.end()) {
100 fProcCounter[procName] = localCount; 115 fProcCounter[procName] = localCount;
101 } 116 }
102 else { 117 else {
103 fProcCounter[procName] += localCount; 118 fProcCounter[procName] += localCount;
104 } << 119 }
105 } 120 }
106 << 121
107 G4Run::Merge(run); << 122 G4Run::Merge(run);
108 } << 123 }
109 124
110 //....oooOO0OOooo........oooOO0OOooo........oo 125 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
111 126
112 void Run::EndOfRun() 127 void Run::EndOfRun()
113 { 128 {
114 G4int prec = 5, wid = prec + 2; << 129 G4int prec = 5, wid = prec + 2;
115 G4int dfprec = G4cout.precision(prec); 130 G4int dfprec = G4cout.precision(prec);
116 << 131
117 // run condition << 132 //run condition
118 // << 133 //
119 G4String partName = fParticle->GetParticleNa << 134 G4String partName = fParticle->GetParticleName();
120 const G4Material* material = fDetector->GetM << 135 G4Material* material = fDetector->GetMaterial();
121 G4double density = material->GetDensity(); << 136 G4double density = material->GetDensity();
122 << 137
123 G4cout << "\n ======================== run s 138 G4cout << "\n ======================== run summary ======================\n";
124 G4cout << "\n The run is: " << numberOfEvent 139 G4cout << "\n The run is: " << numberOfEvent << " " << partName << " of "
125 << G4BestUnit(fEkin, "Energy") << " t << 140 << G4BestUnit(fEkin,"Energy") << " through "
126 << " of " << material->GetName() << " << 141 << G4BestUnit(fDetector->GetSize(),"Length") << " of "
127 << ")" << G4endl; << 142 << material->GetName() << " (density: "
128 << 143 << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
129 if (numberOfEvent == 0) { << 144
130 G4cout.precision(dfprec); << 145 if (numberOfEvent == 0) { G4cout.precision(dfprec); return;}
131 return; << 146
132 } << 147 G4double dNbOfEvents = double(numberOfEvent);
133 << 148 G4cout << "\n total energy deposit: "
134 G4double dNbOfEvents = (G4double)numberOfEve << 149 << G4BestUnit(fEdep/dNbOfEvents, "Energy") << G4endl;
135 G4cout << "\n Energy deposit: " << 150
136 << G4BestUnit(fEdep/dNbOfEvents, "En << 151 //nb of tracks and steps per event
137 G4cout << " Energy leakage: " << 152 //
138 << G4BestUnit(fEleak/dNbOfEvents, "En << 153 G4cout << "\n nb tracks/event"
139 G4cout << " Edep + Eleak: " << 154 << " neutral: " << std::setw(wid) << fNbOfTraks0/dNbOfEvents
140 << G4BestUnit((fEdep+fEleak)/dNbOfEve << 155 << " charged: " << std::setw(wid) << fNbOfTraks1/dNbOfEvents
141 G4cout << " \n NIEL energy calculated: " << 156 << "\n nb steps/event"
142 << G4BestUnit(fNIEL/dNbOfEvents, "En << 157 << " neutral: " << std::setw(wid) << fNbOfSteps0/dNbOfEvents
143 << 158 << " charged: " << std::setw(wid) << fNbOfSteps1/dNbOfEvents
144 // nb of tracks and steps per event << 159 << G4endl;
145 // << 160
146 G4cout << "\n Nb tracks/event" << 161 //frequency of processes
147 << " neutral: " << std::setw(wid) < <<
148 << " charged: " << std::setw(wid) < <<
149 << " neutral: " << std::setw(wid) < <<
150 << " charged: " << std::setw(wid) < <<
151 <<
152 // frequency of processes <<
153 // 162 //
154 G4cout << "\n Process calls frequency :" << 163 G4cout << "\n Process calls frequency :" << G4endl;
155 G4int index = 0; << 164 G4int index = 0;
156 std::map<G4String, G4int>::iterator it; << 165 std::map<G4String,G4int>::iterator it;
157 for (it = fProcCounter.begin(); it != fProcC 166 for (it = fProcCounter.begin(); it != fProcCounter.end(); it++) {
158 G4String procName = it->first; << 167 G4String procName = it->first;
159 G4int count = it->second; << 168 G4int count = it->second;
160 G4String space = " "; << 169 G4String space = " "; if (++index%3 == 0) space = "\n";
161 if (++index % 3 == 0) space = "\n"; << 170 G4cout << " " << std::setw(20) << procName << "="<< std::setw(7) << count
162 G4cout << " " << std::setw(20) << procName << 171 << space;
163 } << 172 }
164 G4cout << G4endl; 173 G4cout << G4endl;
165 << 174
166 // compute true and projected ranges, and tr << 175 //compute true and projected ranges, and transverse dispersion
167 // <<
168 fTrueRange /= numberOfEvent; <<
169 fTrueRange2 /= numberOfEvent; <<
170 G4double trueRms = fTrueRange2 - fTrueRange <<
171 if (trueRms > 0.) <<
172 trueRms = std::sqrt(trueRms); <<
173 else <<
174 trueRms = 0.; <<
175 <<
176 fProjRange /= numberOfEvent; <<
177 fProjRange2 /= numberOfEvent; <<
178 G4double projRms = fProjRange2 - fProjRange <<
179 if (projRms > 0.) <<
180 projRms = std::sqrt(projRms); <<
181 else <<
182 projRms = 0.; <<
183 <<
184 fTransvDev /= 2 * numberOfEvent; <<
185 fTransvDev2 /= 2 * numberOfEvent; <<
186 G4double trvsRms = fTransvDev2 - fTransvDev <<
187 if (trvsRms > 0.) <<
188 trvsRms = std::sqrt(trvsRms); <<
189 else <<
190 trvsRms = 0.; <<
191 <<
192 // compare true range with csda range from P <<
193 // 176 //
>> 177 fTrueRange /= numberOfEvent; fTrueRange2 /= numberOfEvent;
>> 178 G4double trueRms = fTrueRange2 - fTrueRange*fTrueRange;
>> 179 if (trueRms>0.) trueRms = std::sqrt(trueRms); else trueRms = 0.;
>> 180
>> 181 fProjRange /= numberOfEvent; fProjRange2 /= numberOfEvent;
>> 182 G4double projRms = fProjRange2 - fProjRange*fProjRange;
>> 183 if (projRms>0.) projRms = std::sqrt(projRms); else projRms = 0.;
>> 184
>> 185 fTransvDev /= 2*numberOfEvent; fTransvDev2 /= 2*numberOfEvent;
>> 186 G4double trvsRms = fTransvDev2 - fTransvDev*fTransvDev;
>> 187 if (trvsRms>0.) trvsRms = std::sqrt(trvsRms); else trvsRms = 0.;
>> 188
>> 189 //compare true range with csda range from PhysicsTables
>> 190 //
194 G4EmCalculator emCalculator; 191 G4EmCalculator emCalculator;
195 G4double rangeTable = 0.; 192 G4double rangeTable = 0.;
196 if (fParticle->GetPDGCharge() != 0.) 193 if (fParticle->GetPDGCharge() != 0.)
197 rangeTable = emCalculator.GetCSDARange(fEk << 194 rangeTable = emCalculator.GetCSDARange(fEkin,fParticle,material);
198 << 195
199 G4cout << "\n------------------------------- 196 G4cout << "\n---------------------------------------------------------\n";
200 G4cout << " Primary particle : "; << 197 G4cout << " Primary particle : " ;
201 G4cout << "\n true Range = " << G4BestUnit(f << 198 G4cout << "\n true Range = " << G4BestUnit(fTrueRange,"Length")
202 << " rms = " << G4BestUnit(trueRms, << 199 << " rms = " << G4BestUnit(trueRms, "Length");
203 << 200
204 G4cout << "\n proj Range = " << G4BestUnit(f << 201 G4cout << "\n proj Range = " << G4BestUnit(fProjRange,"Length")
205 << " rms = " << G4BestUnit(projRms, << 202 << " rms = " << G4BestUnit(projRms, "Length");
206 << 203
207 G4cout << "\n proj/true = " << fProjRange / << 204 G4cout << "\n proj/true = " << fProjRange/fTrueRange;
208 << 205
209 G4cout << "\n transverse dispersion at end = << 206 G4cout << "\n transverse dispersion at end = "
210 << 207 << G4BestUnit(trvsRms,"Length");
211 G4cout << "\n mass true Range from simu << 208
212 << G4BestUnit(fTrueRange * density, " << 209 G4cout << "\n mass true Range from simulation = "
213 << "\n from PhysicsTable (csda << 210 << G4BestUnit(fTrueRange*density, "Mass/Surface")
214 << G4BestUnit(rangeTable * density, " << 211 << "\n from PhysicsTable (csda range) = "
>> 212 << G4BestUnit(rangeTable*density, "Mass/Surface");
215 G4cout << "\n------------------------------- 213 G4cout << "\n---------------------------------------------------------\n";
216 G4cout << G4endl; 214 G4cout << G4endl;
217 << 215
218 // remove all contents in fProcCounter << 216 // remove all contents in fProcCounter
219 fProcCounter.clear(); 217 fProcCounter.clear();
220 << 218
221 // restore default format << 219 //restore default format
222 G4cout.precision(dfprec); << 220 G4cout.precision(dfprec);
223 } 221 }
224 222
225 //....oooOO0OOooo........oooOO0OOooo........oo 223 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
226 224