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