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/TestEm5/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 "HistoManager.hh" 37 #include "PrimaryGeneratorAction.hh" 38 39 #include "G4EmCalculator.hh" 40 #include "G4SystemOfUnits.hh" 41 #include "G4UnitsTable.hh" 42 43 #include <iomanip> 44 45 //....oooOO0OOooo........oooOO0OOooo........oo 46 47 Run::Run(DetectorConstruction* det) : fDetecto 48 49 //....oooOO0OOooo........oooOO0OOooo........oo 50 51 void Run::SetPrimary(G4ParticleDefinition* par 52 { 53 fParticle = particle; 54 fEkin = energy; 55 56 // compute theta0 57 fMscThetaCentral = 3 * ComputeMscHighland(); 58 } 59 60 //....oooOO0OOooo........oooOO0OOooo........oo 61 62 void Run::Merge(const G4Run* run) 63 { 64 const Run* localRun = static_cast<const Run* 65 66 // pass information about primary particle 67 fParticle = localRun->fParticle; 68 fEkin = localRun->fEkin; 69 70 fMscThetaCentral = localRun->fMscThetaCentra 71 72 // accumulate sums 73 // 74 fEnergyDeposit += localRun->fEnergyDeposit; 75 fEnergyDeposit2 += localRun->fEnergyDeposit2 76 fTrakLenCharged += localRun->fTrakLenCharged 77 fTrakLenCharged2 += localRun->fTrakLenCharge 78 fTrakLenNeutral += localRun->fTrakLenNeutral 79 fTrakLenNeutral2 += localRun->fTrakLenNeutra 80 fNbStepsCharged += localRun->fNbStepsCharged 81 fNbStepsCharged2 += localRun->fNbStepsCharge 82 fNbStepsNeutral += localRun->fNbStepsNeutral 83 fNbStepsNeutral2 += localRun->fNbStepsNeutra 84 fMscProjecTheta += localRun->fMscProjecTheta 85 fMscProjecTheta2 += localRun->fMscProjecThet 86 87 fTypes[0] += localRun->fTypes[0]; 88 fTypes[1] += localRun->fTypes[1]; 89 fTypes[2] += localRun->fTypes[2]; 90 fTypes[3] += localRun->fTypes[3]; 91 92 fNbGamma += localRun->fNbGamma; 93 fNbElect += localRun->fNbElect; 94 fNbPosit += localRun->fNbPosit; 95 96 fTransmit[0] += localRun->fTransmit[0]; 97 fTransmit[1] += localRun->fTransmit[1]; 98 fReflect[0] += localRun->fReflect[0]; 99 fReflect[1] += localRun->fReflect[1]; 100 101 fMscEntryCentral += localRun->fMscEntryCentr 102 103 fEnergyLeak[0] += localRun->fEnergyLeak[0]; 104 fEnergyLeak[1] += localRun->fEnergyLeak[1]; 105 fEnergyLeak2[0] += localRun->fEnergyLeak2[0] 106 fEnergyLeak2[1] += localRun->fEnergyLeak2[1] 107 108 G4Run::Merge(run); 109 } 110 111 //....oooOO0OOooo........oooOO0OOooo........oo 112 113 void Run::EndOfRun() 114 { 115 // compute mean and rms 116 // 117 G4int TotNbofEvents = numberOfEvent; 118 if (TotNbofEvents == 0) return; 119 120 G4double EnergyBalance = fEnergyDeposit + fE 121 EnergyBalance /= TotNbofEvents; 122 123 fEnergyDeposit /= TotNbofEvents; 124 fEnergyDeposit2 /= TotNbofEvents; 125 G4double rmsEdep = fEnergyDeposit2 - fEnergy 126 if (rmsEdep > 0.) 127 rmsEdep = std::sqrt(rmsEdep / TotNbofEvent 128 else 129 rmsEdep = 0.; 130 131 fTrakLenCharged /= TotNbofEvents; 132 fTrakLenCharged2 /= TotNbofEvents; 133 G4double rmsTLCh = fTrakLenCharged2 - fTrakL 134 if (rmsTLCh > 0.) 135 rmsTLCh = std::sqrt(rmsTLCh / TotNbofEvent 136 else 137 rmsTLCh = 0.; 138 139 fTrakLenNeutral /= TotNbofEvents; 140 fTrakLenNeutral2 /= TotNbofEvents; 141 G4double rmsTLNe = fTrakLenNeutral2 - fTrakL 142 if (rmsTLNe > 0.) 143 rmsTLNe = std::sqrt(rmsTLNe / TotNbofEvent 144 else 145 rmsTLNe = 0.; 146 147 fNbStepsCharged /= TotNbofEvents; 148 fNbStepsCharged2 /= TotNbofEvents; 149 G4double rmsStCh = fNbStepsCharged2 - fNbSte 150 if (rmsStCh > 0.) 151 rmsStCh = std::sqrt(rmsStCh / TotNbofEvent 152 else 153 rmsStCh = 0.; 154 155 fNbStepsNeutral /= TotNbofEvents; 156 fNbStepsNeutral2 /= TotNbofEvents; 157 G4double rmsStNe = fNbStepsNeutral2 - fNbSte 158 if (rmsStNe > 0.) 159 rmsStNe = std::sqrt(rmsStNe / TotNbofEvent 160 else 161 rmsStNe = 0.; 162 163 G4double Gamma = (G4double)fNbGamma / TotNbo 164 G4double Elect = (G4double)fNbElect / TotNbo 165 G4double Posit = (G4double)fNbPosit / TotNbo 166 167 G4double transmit[2]; 168 transmit[0] = 100. * fTransmit[0] / TotNbofE 169 transmit[1] = 100. * fTransmit[1] / TotNbofE 170 171 G4double reflect[2]; 172 reflect[0] = 100. * fReflect[0] / TotNbofEve 173 reflect[1] = 100. * fReflect[1] / TotNbofEve 174 175 G4double rmsMsc = 0., tailMsc = 0.; 176 if (fMscEntryCentral > 0) { 177 fMscProjecTheta /= fMscEntryCentral; 178 fMscProjecTheta2 /= fMscEntryCentral; 179 rmsMsc = fMscProjecTheta2 - fMscProjecThet 180 if (rmsMsc > 0.) { 181 rmsMsc = std::sqrt(rmsMsc); 182 } 183 if (fTransmit[1] > 0.0) { 184 tailMsc = 100. - (100. * fMscEntryCentra 185 } 186 } 187 188 fEnergyLeak[0] /= TotNbofEvents; 189 fEnergyLeak2[0] /= TotNbofEvents; 190 G4double rmsEl0 = fEnergyLeak2[0] - fEnergyL 191 if (rmsEl0 > 0.) 192 rmsEl0 = std::sqrt(rmsEl0 / TotNbofEvents) 193 else 194 rmsEl0 = 0.; 195 196 fEnergyLeak[1] /= TotNbofEvents; 197 fEnergyLeak2[1] /= TotNbofEvents; 198 G4double rmsEl1 = fEnergyLeak2[1] - fEnergyL 199 if (rmsEl1 > 0.) 200 rmsEl1 = std::sqrt(rmsEl1 / TotNbofEvents) 201 else 202 rmsEl1 = 0.; 203 204 // Stopping Power from input Table. 205 // 206 const G4Material* material = fDetector->GetA 207 G4double length = fDetector->GetAbsorberThic 208 G4double density = material->GetDensity(); 209 G4String partName = fParticle->GetParticleNa 210 211 G4EmCalculator emCalculator; 212 G4double dEdxTable = 0., dEdxFull = 0.; 213 if (fParticle->GetPDGCharge() != 0.) { 214 dEdxTable = emCalculator.GetDEDX(fEkin, fP 215 dEdxFull = emCalculator.ComputeTotalDEDX(f 216 } 217 G4double stopTable = dEdxTable / density; 218 G4double stopFull = dEdxFull / density; 219 220 // Stopping Power from simulation. 221 // 222 G4double meandEdx = fEnergyDeposit / length; 223 G4double stopPower = meandEdx / density; 224 225 G4cout << "\n ======================== run s 226 227 G4int prec = G4cout.precision(3); 228 229 G4cout << "\n The run was " << TotNbofEvents 230 << G4BestUnit(fEkin, "Energy") << " t 231 << material->GetName() << " (density: 232 << G4endl; 233 234 G4cout.precision(4); 235 236 G4cout << "\n Total energy deposit in absorb 237 << G4BestUnit(fEnergyDeposit, "Energy 238 << G4endl; 239 240 G4cout << "\n -----> Mean dE/dx = " << meand 241 << "\t(" << stopPower / (MeV * cm2 / 242 243 G4cout << "\n From formulas :" << G4endl; 244 G4cout << " restricted dEdx = " << dEdxTab 245 << "\t(" << stopTable / (MeV * cm2 / 246 247 G4cout << " full dEdx = " << dEdxFul 248 << "\t(" << stopFull / (MeV * cm2 / g 249 250 G4cout << "\n Leakage : primary = " << G4Be 251 << G4BestUnit(rmsEl0, "Energy") 252 << " secondaries = " << G4BestUnit( 253 << G4BestUnit(rmsEl1, "Energy") << G4 254 255 G4cout << " Energy balance : edep + eleak = 256 257 G4cout << "\n Total track length (charged) i 258 << G4BestUnit(fTrakLenCharged, "Lengt 259 << G4endl; 260 261 G4cout << " Total track length (neutral) in 262 << G4BestUnit(fTrakLenNeutral, "Lengt 263 << G4endl; 264 265 G4cout << "\n Number of steps (charged) in a 266 << rmsStCh << G4endl; 267 268 G4cout << " Number of steps (neutral) in abs 269 << rmsStNe << G4endl; 270 271 G4cout << "\n Number of secondaries per even 272 << "; positrons = " << Posit << G4e 273 274 G4cout << "\n Number of events with the prim 275 << G4endl; 276 277 G4cout << " Number of events with at least 278 << "(same charge as primary) = " << t 279 280 G4cout << "\n Number of events with the prim 281 << G4endl; 282 283 G4cout << " Number of events with at least 284 << "(same charge as primary) = " << r 285 286 // compute width of the Gaussian central par 287 // 288 G4cout << "\n MultipleScattering:" 289 << "\n rms proj angle of transmit pr 290 << " mrad (central part only)" << G4e 291 292 G4cout << " computed theta0 (Highland formu 293 << " mrad" << G4endl; 294 295 G4cout << " central part defined as +- " << 296 << " Tail ratio = " << tailMsc << " 297 298 // gamma process counts 299 // 300 G4cout << "\n Gamma process counts:" << G4en 301 G4cout << " Photoeffect " << fTypes[0] << 302 G4cout << " Compton " << fTypes[1] << 303 G4cout << " Conversion " << fTypes[2] << 304 G4cout << " Rayleigh " << fTypes[3] << 305 306 // normalize histograms 307 // 308 G4AnalysisManager* analysisManager = G4Analy 309 310 G4int ih = 1; 311 G4double binWidth = analysisManager->GetH1Wi 312 G4double fac = 1. / (TotNbofEvents * binWidt 313 analysisManager->ScaleH1(ih, fac); 314 315 ih = 10; 316 binWidth = analysisManager->GetH1Width(ih); 317 fac = 1. / (TotNbofEvents * binWidth); 318 analysisManager->ScaleH1(ih, fac); 319 320 ih = 12; 321 analysisManager->ScaleH1(ih, 1. / TotNbofEve 322 323 // reset default precision 324 G4cout.precision(prec); 325 } 326 327 //....oooOO0OOooo........oooOO0OOooo........oo 328 329 G4double Run::ComputeMscHighland() 330 { 331 // compute the width of the Gaussian central 332 // projected angular distribution. 333 // Eur. Phys. Jour. C15 (2000) page 166, for 334 335 G4double t = 336 (fDetector->GetAbsorberThickness()) / (fDe 337 if (t < DBL_MIN) return 0.; 338 339 G4double T = fEkin; 340 G4double M = fParticle->GetPDGMass(); 341 G4double z = std::abs(fParticle->GetPDGCharg 342 343 G4double bpc = T * (T + 2 * M) / (T + M); 344 G4double teta0 = 13.6 * MeV * z * std::sqrt( 345 return teta0; 346 } 347 348 //....oooOO0OOooo........oooOO0OOooo........oo 349