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/TestEm12/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 "G4Material.hh" 40 #include "G4SystemOfUnits.hh" 41 #include "G4UnitsTable.hh" 42 43 //....oooOO0OOooo........oooOO0OOooo........oo 44 45 Run::Run(DetectorConstruction* detector) : fDe 46 47 //....oooOO0OOooo........oooOO0OOooo........oo 48 49 void Run::SetPrimary(G4ParticleDefinition* par 50 { 51 fParticle = particle; 52 fEkin = energy; 53 } 54 55 //....oooOO0OOooo........oooOO0OOooo........oo 56 57 void Run::AddEdep(G4double e) 58 { 59 fEdeposit += e; 60 fEdeposit2 += e * e; 61 } 62 63 //....oooOO0OOooo........oooOO0OOooo........oo 64 65 void Run::AddTrackLength(G4double t) 66 { 67 fTrackLen += t; 68 fTrackLen2 += t * t; 69 } 70 71 //....oooOO0OOooo........oooOO0OOooo........oo 72 73 void Run::AddProjRange(G4double x) 74 { 75 fProjRange += x; 76 fProjRange2 += x * x; 77 } 78 79 //....oooOO0OOooo........oooOO0OOooo........oo 80 81 void Run::AddStepSize(G4int nb, G4double st) 82 { 83 fNbOfSteps += nb; 84 fNbOfSteps2 += nb * nb; 85 fStepSize += st; 86 fStepSize2 += st * st; 87 } 88 89 //....oooOO0OOooo........oooOO0OOooo........oo 90 91 void Run::SetCsdaRange(G4double value) 92 { 93 fCsdaRange = value; 94 } 95 96 //....oooOO0OOooo........oooOO0OOooo........oo 97 98 G4double Run::GetCsdaRange() 99 { 100 return fCsdaRange; 101 } 102 103 //....oooOO0OOooo........oooOO0OOooo........oo 104 105 void Run::Merge(const G4Run* run) 106 { 107 const Run* localRun = static_cast<const Run* 108 109 // pass information about primary particle 110 fParticle = localRun->fParticle; 111 fEkin = localRun->fEkin; 112 113 // accumulate sums 114 fEdeposit += localRun->fEdeposit; 115 fEdeposit2 += localRun->fEdeposit2; 116 fTrackLen += localRun->fTrackLen; 117 fTrackLen2 += localRun->fTrackLen2; 118 fProjRange += localRun->fProjRange; 119 fProjRange2 += localRun->fProjRange2; 120 fNbOfSteps += localRun->fNbOfSteps; 121 fNbOfSteps2 += localRun->fNbOfSteps2; 122 fStepSize += localRun->fStepSize; 123 fStepSize2 += localRun->fStepSize2; 124 125 fCsdaRange = localRun->fCsdaRange; 126 127 G4Run::Merge(run); 128 } 129 130 //....oooOO0OOooo........oooOO0OOooo........oo 131 132 void Run::EndOfRun() 133 { 134 std::ios::fmtflags mode = G4cout.flags(); 135 G4cout.setf(std::ios::fixed, std::ios::float 136 G4int prec = G4cout.precision(2); 137 138 // run conditions 139 // 140 G4Material* material = fDetector->GetAbsorMa 141 G4double density = material->GetDensity(); 142 G4String partName = fParticle->GetParticleNa 143 144 G4cout << "\n ======================== run s 145 G4cout << "\n The run is " << numberOfEvent 146 << G4BestUnit(fEkin, "Energy") << " t 147 << G4BestUnit(fDetector->GetAbsorRadi 148 << " (density: " << G4BestUnit(densit 149 150 if (numberOfEvent == 0) { 151 G4cout.setf(mode, std::ios::floatfield); 152 G4cout.precision(prec); 153 return; 154 } 155 156 fEdeposit /= numberOfEvent; 157 fEdeposit2 /= numberOfEvent; 158 G4double rms = fEdeposit2 - fEdeposit * fEde 159 if (rms > 0.) 160 rms = std::sqrt(rms); 161 else 162 rms = 0.; 163 164 G4cout.precision(3); 165 G4cout << "\n Total Energy deposited 166 << G4BestUnit(rms, "Energy") << G4end 167 168 // compute track length of primary track 169 // 170 fTrackLen /= numberOfEvent; 171 fTrackLen2 /= numberOfEvent; 172 rms = fTrackLen2 - fTrackLen * fTrackLen; 173 if (rms > 0.) 174 rms = std::sqrt(rms); 175 else 176 rms = 0.; 177 178 G4cout.precision(3); 179 G4cout << "\n Track length of primary track 180 << G4BestUnit(rms, "Length"); 181 182 // compare with csda range 183 // 184 G4cout << "\n Range from EmCalculator = " << 185 << " (from full dE/dx)" << G4endl; 186 187 // compute projected range of primary track 188 // 189 fProjRange /= numberOfEvent; 190 fProjRange2 /= numberOfEvent; 191 rms = fProjRange2 - fProjRange * fProjRange; 192 if (rms > 0.) 193 rms = std::sqrt(rms); 194 else 195 rms = 0.; 196 197 G4cout << "\n Projected range 198 << G4BestUnit(rms, "Length") << G4end 199 200 // nb of steps and step size of primary trac 201 // 202 G4double dNofEvents = double(numberOfEvent); 203 G4double fNbSteps = fNbOfSteps / dNofEvents, 204 rms = fNbSteps2 - fNbSteps * fNbSteps; 205 if (rms > 0.) 206 rms = std::sqrt(rms); 207 else 208 rms = 0.; 209 210 G4cout.precision(2); 211 G4cout << "\n Nb of steps of primary track 212 213 fStepSize /= numberOfEvent; 214 fStepSize2 /= numberOfEvent; 215 rms = fStepSize2 - fStepSize * fStepSize; 216 if (rms > 0.) 217 rms = std::sqrt(rms); 218 else 219 rms = 0.; 220 221 G4cout.precision(3); 222 G4cout << "\t Step size= " << G4BestUnit(fSt 223 << G4BestUnit(rms, "Length") << G4end 224 225 // normalize histograms of longitudinal ener 226 // 227 G4AnalysisManager* analysisManager = G4Analy 228 G4int ih = 1; 229 G4double binWidth = analysisManager->GetH1Wi 230 G4double fac = (1. / (numberOfEvent * binWid 231 analysisManager->ScaleH1(ih, fac); 232 233 // normalize histogram d(E/E0)/d(r/r0) 234 // 235 ih = 8; 236 binWidth = analysisManager->GetH1Width(ih); 237 fac = 1. / (numberOfEvent * binWidth * fEkin 238 analysisManager->ScaleH1(ih, fac); 239 240 // reset default formats 241 G4cout.setf(mode, std::ios::floatfield); 242 G4cout.precision(prec); 243 } 244 245 //....oooOO0OOooo........oooOO0OOooo........oo 246