Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 /// \file electromagnetic/TestEm11/src/Run.cc 27 /// \brief Implementation of the Run class 28 // 29 // 30 //....oooOO0OOooo........oooOO0OOooo........oo 31 //....oooOO0OOooo........oooOO0OOooo........oo 32 33 #include "Run.hh" 34 35 #include "DetectorConstruction.hh" 36 #include "EventAction.hh" 37 #include "HistoManager.hh" 38 #include "PrimaryGeneratorAction.hh" 39 40 #include "G4Event.hh" 41 #include "G4Material.hh" 42 #include "G4SystemOfUnits.hh" 43 #include "G4UnitsTable.hh" 44 45 #include <iomanip> 46 47 //....oooOO0OOooo........oooOO0OOooo........oo 48 49 Run::Run(DetectorConstruction* detector) : fDe 50 { 51 for (G4int i = 0; i < 3; ++i) { 52 fStatus[i] = 0; 53 fTotEdep[i] = fEleak[i] = fEtotal[i] = 0.; 54 } 55 fTotEdep[1] = fEleak[1] = fEtotal[1] = joule 56 57 for (G4int i = 0; i < kMaxAbsor; ++i) { 58 fEdeposit[i] = 0.; 59 fEmin[i] = joule; 60 fEmax[i] = 0.; 61 fCsdaRange[i] = 0.; 62 fXfrontNorm[i] = 0.; 63 } 64 } 65 66 //....oooOO0OOooo........oooOO0OOooo........oo 67 68 void Run::SetPrimary(G4ParticleDefinition* par 69 { 70 fParticle = particle; 71 fEkin = energy; 72 } 73 74 //....oooOO0OOooo........oooOO0OOooo........oo 75 76 void Run::AddEdep(G4int i, G4double e) 77 { 78 if (e > 0.) { 79 fEdeposit[i] += e; 80 if (e < fEmin[i]) fEmin[i] = e; 81 if (e > fEmax[i]) fEmax[i] = e; 82 } 83 } 84 85 //....oooOO0OOooo........oooOO0OOooo........oo 86 87 void Run::AddTotEdep(G4double e) 88 { 89 if (e > 0.) { 90 fTotEdep[0] += e; 91 if (e < fTotEdep[1]) fTotEdep[1] = e; 92 if (e > fTotEdep[2]) fTotEdep[2] = e; 93 } 94 } 95 96 //....oooOO0OOooo........oooOO0OOooo........oo 97 98 void Run::AddEleak(G4double e) 99 { 100 if (e > 0.) { 101 fEleak[0] += e; 102 if (e < fEleak[1]) fEleak[1] = e; 103 if (e > fEleak[2]) fEleak[2] = e; 104 } 105 } 106 107 //....oooOO0OOooo........oooOO0OOooo........oo 108 109 void Run::AddEtotal(G4double e) 110 { 111 if (e > 0.) { 112 fEtotal[0] += e; 113 if (e < fEtotal[1]) fEtotal[1] = e; 114 if (e > fEtotal[2]) fEtotal[2] = e; 115 } 116 } 117 118 //....oooOO0OOooo........oooOO0OOooo........oo 119 120 void Run::AddTrackLength(G4double t) 121 { 122 fTrackLen += t; 123 fTrackLen2 += t * t; 124 } 125 126 //....oooOO0OOooo........oooOO0OOooo........oo 127 128 void Run::AddProjRange(G4double x) 129 { 130 fProjRange += x; 131 fProjRange2 += x * x; 132 } 133 134 //....oooOO0OOooo........oooOO0OOooo........oo 135 136 void Run::AddStepSize(G4int nb, G4double st) 137 { 138 fNbOfSteps += nb; 139 fNbOfSteps2 += nb * nb; 140 fStepSize += st; 141 fStepSize2 += st * st; 142 } 143 144 //....oooOO0OOooo........oooOO0OOooo........oo 145 146 void Run::AddTrackStatus(G4int i) 147 { 148 fStatus[i]++; 149 } 150 151 //....oooOO0OOooo........oooOO0OOooo........oo 152 153 void Run::SetCsdaRange(G4int i, G4double value 154 { 155 fCsdaRange[i] = value; 156 } 157 158 //....oooOO0OOooo........oooOO0OOooo........oo 159 160 void Run::SetXfrontNorm(G4int i, G4double valu 161 { 162 fXfrontNorm[i] = value; 163 } 164 165 //....oooOO0OOooo........oooOO0OOooo........oo 166 167 G4double Run::GetCsdaRange(G4int i) 168 { 169 return fCsdaRange[i]; 170 } 171 172 //....oooOO0OOooo........oooOO0OOooo........oo 173 174 G4double Run::GetXfrontNorm(G4int i) 175 { 176 return fXfrontNorm[i]; 177 } 178 179 //....oooOO0OOooo........oooOO0OOooo........oo 180 181 void Run::Merge(const G4Run* run) 182 { 183 const Run* localRun = static_cast<const Run* 184 185 // pass information about primary particle 186 fParticle = localRun->fParticle; 187 fEkin = localRun->fEkin; 188 189 // accumulate sums 190 fTrackLen += localRun->fTrackLen; 191 fTrackLen2 += localRun->fTrackLen2; 192 fProjRange += localRun->fProjRange; 193 fProjRange2 += localRun->fProjRange2; 194 fNbOfSteps += localRun->fNbOfSteps; 195 fNbOfSteps2 += localRun->fNbOfSteps2; 196 fStepSize += localRun->fStepSize; 197 fStepSize2 += localRun->fStepSize2; 198 199 G4int nbOfAbsor = fDetector->GetNbOfAbsor(); 200 for (G4int i = 1; i <= nbOfAbsor; ++i) { 201 fEdeposit[i] += localRun->fEdeposit[i]; 202 fCsdaRange[i] = localRun->fCsdaRange[i]; 203 fXfrontNorm[i] = localRun->fXfrontNorm[i]; 204 // min, max 205 G4double min, max; 206 min = localRun->fEmin[i]; 207 max = localRun->fEmax[i]; 208 if (fEmin[i] > min) fEmin[i] = min; 209 if (fEmax[i] < max) fEmax[i] = max; 210 } 211 212 for (G4int i = 0; i < 3; ++i) 213 fStatus[i] += localRun->fStatus[i]; 214 215 // total Edep 216 fTotEdep[0] += localRun->fTotEdep[0]; 217 G4double min, max; 218 min = localRun->fTotEdep[1]; 219 max = localRun->fTotEdep[2]; 220 if (fTotEdep[1] > min) fTotEdep[1] = min; 221 if (fTotEdep[2] < max) fTotEdep[2] = max; 222 223 // Eleak 224 fEleak[0] += localRun->fEleak[0]; 225 min = localRun->fEleak[1]; 226 max = localRun->fEleak[2]; 227 if (fEleak[1] > min) fEleak[1] = min; 228 if (fEleak[2] < max) fEleak[2] = max; 229 230 // Etotal 231 fEtotal[0] += localRun->fEtotal[0]; 232 min = localRun->fEtotal[1]; 233 max = localRun->fEtotal[2]; 234 if (fEtotal[1] > min) fEtotal[1] = min; 235 if (fEtotal[2] < max) fEtotal[2] = max; 236 237 G4Run::Merge(run); 238 } 239 240 //....oooOO0OOooo........oooOO0OOooo........oo 241 242 void Run::EndOfRun() 243 { 244 std::ios::fmtflags mode = G4cout.flags(); 245 G4cout.setf(std::ios::fixed, std::ios::float 246 G4int prec = G4cout.precision(2); 247 248 // run conditions 249 // 250 G4String partName = fParticle->GetParticleNa 251 G4int nbOfAbsor = fDetector->GetNbOfAbsor(); 252 253 G4cout << "\n ======================== run s 254 G4cout << "\n The run is " << numberOfEvent 255 << G4BestUnit(fEkin, "Energy") << " t 256 for (G4int i = 1; i <= nbOfAbsor; i++) { 257 G4Material* material = fDetector->GetAbsor 258 G4double thickness = fDetector->GetAbsorTh 259 G4double density = material->GetDensity(); 260 G4cout << std::setw(5) << i << std::setw(1 261 << material->GetName() << " (densit 262 << G4endl; 263 } 264 265 if (numberOfEvent == 0) { 266 G4cout.setf(mode, std::ios::floatfield); 267 G4cout.precision(prec); 268 return; 269 } 270 271 G4cout.precision(3); 272 G4double rms(0); 273 274 // Edep in absorbers 275 // 276 for (G4int i = 1; i <= nbOfAbsor; i++) { 277 fEdeposit[i] /= numberOfEvent; 278 279 G4cout << "\n Edep in absorber " << i << " 280 << G4BestUnit(fEmin[i], "Energy") < 281 } 282 G4cout << G4endl; 283 284 if (nbOfAbsor > 1) { 285 fTotEdep[0] /= numberOfEvent; 286 G4cout << "\n Edep in all absorbers = " << 287 << G4BestUnit(fTotEdep[1], "Energy" 288 << G4endl; 289 } 290 291 // Eleak 292 // 293 fEleak[0] /= numberOfEvent; 294 G4cout << " Energy leakage = " << G4Best 295 << G4BestUnit(fEleak[1], "Energy") << 296 << G4endl; 297 298 // Etotal 299 // 300 fEtotal[0] /= numberOfEvent; 301 G4cout << " Energy total = " << G4Best 302 << G4BestUnit(fEtotal[1], "Energy") < 303 << G4endl; 304 305 // compute track length of primary track 306 // 307 fTrackLen /= numberOfEvent; 308 fTrackLen2 /= numberOfEvent; 309 rms = fTrackLen2 - fTrackLen * fTrackLen; 310 if (rms > 0.) 311 rms = std::sqrt(rms); 312 else 313 rms = 0.; 314 315 G4cout.precision(3); 316 G4cout << "\n Track length of primary track 317 << G4BestUnit(rms, "Length"); 318 319 // compare with csda range 320 // 321 G4int NbOfAbsor = fDetector->GetNbOfAbsor(); 322 if (NbOfAbsor == 1) { 323 G4cout << "\n Range from EmCalculator = " 324 << " (from full dE/dx)" << G4endl; 325 } 326 327 // compute projected range of primary track 328 // 329 fProjRange /= numberOfEvent; 330 fProjRange2 /= numberOfEvent; 331 rms = fProjRange2 - fProjRange * fProjRange; 332 if (rms > 0.) 333 rms = std::sqrt(rms); 334 else 335 rms = 0.; 336 337 G4cout << "\n Projected range 338 << G4BestUnit(rms, "Length") << G4end 339 340 // nb of steps and step size of primary trac 341 // 342 G4double dNofEvents = double(numberOfEvent); 343 G4double fNbSteps = fNbOfSteps / dNofEvents, 344 rms = fNbSteps2 - fNbSteps * fNbSteps; 345 if (rms > 0.) 346 rms = std::sqrt(rms); 347 else 348 rms = 0.; 349 350 G4cout.precision(2); 351 G4cout << "\n Nb of steps of primary track 352 353 fStepSize /= numberOfEvent; 354 fStepSize2 /= numberOfEvent; 355 rms = fStepSize2 - fStepSize * fStepSize; 356 if (rms > 0.) 357 rms = std::sqrt(rms); 358 else 359 rms = 0.; 360 361 G4cout.precision(3); 362 G4cout << "\t Step size= " << G4BestUnit(fSt 363 << G4BestUnit(rms, "Length") << G4end 364 365 // transmission coefficients 366 // 367 G4double absorbed = 100. * fStatus[0] / dNof 368 G4double transmit = 100. * fStatus[1] / dNof 369 G4double reflected = 100. * fStatus[2] / dNo 370 371 G4cout.precision(2); 372 G4cout << "\n absorbed = " << absorbed << " 373 << " transmit = " << transmit << " 374 << " reflected = " << reflected << 375 << G4endl; 376 377 // normalize histograms of longitudinal ener 378 // 379 G4AnalysisManager* analysisManager = G4Analy 380 G4int ih = 1; 381 G4double binWidth = analysisManager->GetH1Wi 382 G4double fac = (1. / (numberOfEvent * binWid 383 analysisManager->ScaleH1(ih, fac); 384 385 ih = 8; 386 binWidth = analysisManager->GetH1Width(ih); 387 fac = (1. / (numberOfEvent * binWidth)) * (g 388 analysisManager->ScaleH1(ih, fac); 389 390 // reset default formats 391 G4cout.setf(mode, std::ios::floatfield); 392 G4cout.precision(prec); 393 } 394 395 //....oooOO0OOooo........oooOO0OOooo........oo 396