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 // This example is provided by the Geant4-DNA << 27 // Any report or published results obtained us << 28 // shall cite the following Geant4-DNA collabo << 29 // Med. Phys. 45 (2018) e722-e739 << 30 // Phys. Med. 31 (2015) 861-874 << 31 // Med. Phys. 37 (2010) 4692-4708 << 32 // Int. J. Model. Simul. Sci. Comput. 1 (2010) << 33 // << 34 // The Geant4-DNA web site is available at htt << 35 // << 36 /// \file medical/dna/range/src/Run.cc 26 /// \file medical/dna/range/src/Run.cc 37 /// \brief Implementation of the Run class 27 /// \brief Implementation of the Run class >> 28 // >> 29 // $Id: Run.cc 71376 2013-06-14 07:44:50Z maire $ >> 30 // >> 31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 38 33 39 #include "Run.hh" 34 #include "Run.hh" >> 35 #include "DetectorConstruction.hh" 40 36 >> 37 #include "HistoManager.hh" 41 #include "PrimaryGeneratorAction.hh" 38 #include "PrimaryGeneratorAction.hh" 42 39 43 #include "G4Material.hh" 40 #include "G4Material.hh" 44 #include "G4SystemOfUnits.hh" 41 #include "G4SystemOfUnits.hh" 45 #include "G4UnitsTable.hh" 42 #include "G4UnitsTable.hh" 46 43 47 //....oooOO0OOooo........oooOO0OOooo........oo 44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 48 45 49 Run::Run(const DetectorConstruction* detector) 46 Run::Run(const DetectorConstruction* detector) 50 : G4Run(), << 47 : G4Run(), 51 fDetector(detector), << 48 fDetector(detector), 52 fParticle(0), << 49 fParticle(0), fEkin(0.), 53 fEkin(0.), << 50 fEdeposit(0.), fEdeposit2(0.), 54 fEdeposit(0.), << 51 fTrackLen(0.), fTrackLen2(0.), 55 fEdeposit2(0.), << 52 fProjRange(0.), fProjRange2(0.), 56 fTrackLen(0.), << 53 fPenetration(0.), fPenetration2(0.), 57 fTrackLen2(0.), << 54 fNbOfSteps(0), fNbOfSteps2(0), 58 fProjRange(0.), << 55 fStepSize(0.), fStepSize2(0.) 59 fProjRange2(0.), << 56 { } 60 fPenetration(0.), << 61 fPenetration2(0.), << 62 fNbOfSteps(0), << 63 fNbOfSteps2(0), << 64 fStepSize(0.), << 65 fStepSize2(0.) << 66 {} << 67 57 68 //....oooOO0OOooo........oooOO0OOooo........oo 58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 69 59 70 Run::~Run() {} << 60 Run::~Run() >> 61 { } 71 62 72 //....oooOO0OOooo........oooOO0OOooo........oo 63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 73 64 74 void Run::SetPrimary(G4ParticleDefinition* par << 65 void Run::SetPrimary (G4ParticleDefinition* particle, G4double energy) 75 { << 66 { 76 fParticle = particle; 67 fParticle = particle; 77 fEkin = energy; << 68 fEkin = energy; 78 } 69 } 79 70 80 //....oooOO0OOooo........oooOO0OOooo........oo 71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 81 72 82 void Run::AddEdep(G4double e) << 73 void Run::AddEdep (G4double e) 83 { 74 { 84 fEdeposit += e; << 75 fEdeposit += e; 85 fEdeposit2 += e * e; << 76 fEdeposit2 += e*e; 86 } 77 } 87 78 88 //....oooOO0OOooo........oooOO0OOooo........oo 79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 89 << 80 90 void Run::AddTrackLength(G4double t) << 81 void Run::AddTrackLength (G4double t) 91 { 82 { 92 fTrackLen += t; << 83 fTrackLen += t; 93 fTrackLen2 += t * t; << 84 fTrackLen2 += t*t; 94 } 85 } 95 86 96 //....oooOO0OOooo........oooOO0OOooo........oo 87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 97 << 88 98 void Run::AddProjRange(G4double x) << 89 void Run::AddProjRange (G4double x) 99 { 90 { 100 fProjRange += x; << 91 fProjRange += x; 101 fProjRange2 += x * x; << 92 fProjRange2 += x*x; 102 } 93 } 103 94 104 //....oooOO0OOooo........oooOO0OOooo........oo 95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 105 << 96 106 void Run::AddPenetration(G4double x) << 97 void Run::AddPenetration (G4double x) 107 { 98 { 108 fPenetration += x; << 99 fPenetration += x; 109 fPenetration2 += x * x; << 100 fPenetration2 += x*x; 110 } 101 } 111 102 112 //....oooOO0OOooo........oooOO0OOooo........oo 103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 113 << 104 114 void Run::AddStepSize(G4int nb, G4double st) << 105 void Run::AddStepSize (G4int nb, G4double st) 115 { 106 { 116 fNbOfSteps += nb; << 107 fNbOfSteps += nb; 117 fNbOfSteps2 += nb * nb; << 108 fNbOfSteps2 += nb*nb; 118 fStepSize += st; << 109 fStepSize += st ; 119 fStepSize2 += st * st; << 110 fStepSize2 += st*st; 120 } 111 } 121 112 122 //....oooOO0OOooo........oooOO0OOooo........oo 113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 123 114 124 void Run::Merge(const G4Run* run) 115 void Run::Merge(const G4Run* run) 125 { 116 { 126 const Run* localRun = static_cast<const Run* 117 const Run* localRun = static_cast<const Run*>(run); 127 << 118 128 // Pass information about primary particle << 119 // pass information about primary particle 129 fParticle = localRun->fParticle; 120 fParticle = localRun->fParticle; 130 fEkin = localRun->fEkin; << 121 fEkin = localRun->fEkin; 131 122 132 // Accumulate sums << 123 // accumulate sums 133 fEdeposit += localRun->fEdeposit; << 124 fEdeposit += localRun->fEdeposit; 134 fEdeposit2 += localRun->fEdeposit2; << 125 fEdeposit2 += localRun->fEdeposit2; 135 fTrackLen += localRun->fTrackLen; << 126 fTrackLen += localRun->fTrackLen; 136 fTrackLen2 += localRun->fTrackLen2; << 127 fTrackLen2 += localRun->fTrackLen2; 137 fProjRange += localRun->fProjRange; << 128 fProjRange += localRun->fProjRange; 138 fProjRange2 += localRun->fProjRange2; 129 fProjRange2 += localRun->fProjRange2; 139 fPenetration += localRun->fPenetration; << 130 fPenetration += localRun->fPenetration; 140 fPenetration2 += localRun->fPenetration2; 131 fPenetration2 += localRun->fPenetration2; 141 fNbOfSteps += localRun->fNbOfSteps; << 132 fNbOfSteps += localRun->fNbOfSteps ; 142 fNbOfSteps2 += localRun->fNbOfSteps2; 133 fNbOfSteps2 += localRun->fNbOfSteps2; 143 fStepSize += localRun->fStepSize; << 134 fStepSize += localRun->fStepSize; 144 fStepSize2 += localRun->fStepSize2; << 135 fStepSize2 += localRun->fStepSize2; 145 136 146 G4Run::Merge(run); << 137 G4Run::Merge(run); 147 } << 138 } 148 139 149 //....oooOO0OOooo........oooOO0OOooo........oo 140 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 150 141 151 void Run::EndOfRun() << 142 void Run::EndOfRun() 152 { 143 { 153 std::ios::fmtflags mode = G4cout.flags(); 144 std::ios::fmtflags mode = G4cout.flags(); 154 G4cout.setf(std::ios::fixed, std::ios::float << 145 G4cout.setf(std::ios::fixed,std::ios::floatfield); 155 G4int prec = G4cout.precision(2); 146 G4int prec = G4cout.precision(2); 156 << 147 157 // Run conditions << 148 //run conditions >> 149 // 158 G4Material* material = fDetector->GetAbsorMa 150 G4Material* material = fDetector->GetAbsorMaterial(); 159 G4double density = material->GetDensity(); << 151 G4double density = material->GetDensity(); 160 G4String partName = fParticle->GetParticleNa 152 G4String partName = fParticle->GetParticleName(); 161 << 153 162 G4cout << "\n ======================= run su << 154 G4cout << "\n ======================= run summary ====================\n"; 163 G4cout << "\n The run is " << numberOfEvent << 155 G4cout 164 << G4BestUnit(fEkin, "Energy") << " t << 156 << "\n The run is " << numberOfEvent << " "<< partName << " of " 165 << G4BestUnit(fDetector->GetAbsorRadi << 157 << G4BestUnit(fEkin,"Energy") << " through a sphere of radius " 166 << " (density: " << G4BestUnit(densit << 158 << G4BestUnit(fDetector->GetAbsorRadius(),"Length") << "of " >> 159 << material->GetName() << " (density: " >> 160 << G4BestUnit(density,"Volumic Mass") << ")" << G4endl; 167 161 168 if (numberOfEvent == 0) { 162 if (numberOfEvent == 0) { 169 G4cout.setf(mode, std::ios::floatfield); << 163 G4cout.setf(mode,std::ios::floatfield); 170 G4cout.precision(prec); << 164 G4cout.precision(prec); 171 return; 165 return; 172 } 166 } 173 << 167 174 // Compute track length of primary track << 168 //compute track length of primary track 175 fTrackLen /= numberOfEvent; << 169 // 176 fTrackLen2 /= numberOfEvent; << 170 fTrackLen /= numberOfEvent; fTrackLen2 /= numberOfEvent; 177 G4double rmsTrack = fTrackLen2 - fTrackLen * << 171 G4double rmsTrack = fTrackLen2 - fTrackLen*fTrackLen; 178 << 172 179 if (rmsTrack > 0.) << 173 if (rmsTrack>0.) rmsTrack = std::sqrt(rmsTrack); else rmsTrack = 0.; 180 rmsTrack = std::sqrt(rmsTrack); << 174 181 else << 175 G4cout.precision(3); 182 rmsTrack = 0.; << 176 G4cout 183 << 177 << "\n Track length of primary track = " << G4BestUnit(fTrackLen,"Length") 184 G4cout.precision(3); << 178 << " +- " 185 G4cout << "\n Track length of primary track << 179 << G4BestUnit( rmsTrack,"Length"); 186 << G4BestUnit(rmsTrack, "Length"); << 180 187 << 181 //compute projected range of primary track 188 // Compute projected range of primary track << 182 // 189 fProjRange /= numberOfEvent; << 183 fProjRange /= numberOfEvent; fProjRange2 /= numberOfEvent; 190 fProjRange2 /= numberOfEvent; << 184 G4double rmsProj = fProjRange2 - fProjRange*fProjRange; 191 G4double rmsProj = fProjRange2 - fProjRange << 185 if (rmsProj>0.) rmsProj = std::sqrt(rmsProj); else rmsProj = 0.; 192 if (rmsProj > 0.) << 186 193 rmsProj = std::sqrt(rmsProj); << 187 G4cout 194 else << 188 << "\n Projected range = " 195 rmsProj = 0.; << 189 << G4BestUnit(fProjRange,"Length") 196 << 190 << " +- " << G4BestUnit( rmsProj,"Length"); 197 G4cout << "\n Projected range << 191 198 << G4BestUnit(rmsProj, "Length"); << 192 //compute penetration of primary track 199 << 193 // 200 // Compute penetration of primary track << 194 fPenetration /= numberOfEvent; fPenetration2 /= numberOfEvent; 201 fPenetration /= numberOfEvent; << 195 G4double rmsPene = fPenetration2 - fPenetration*fPenetration; 202 fPenetration2 /= numberOfEvent; << 196 if (rmsPene>0.) rmsPene = std::sqrt(rmsPene); else rmsPene = 0.; 203 G4double rmsPene = fPenetration2 - fPenetrat << 197 204 if (rmsPene > 0.) << 198 G4cout 205 rmsPene = std::sqrt(rmsPene); << 199 << "\n Penetration = " 206 else << 200 << G4BestUnit(fPenetration,"Length") 207 rmsPene = 0.; << 201 << " +- " << G4BestUnit( rmsPene,"Length") 208 << 202 << G4endl; 209 G4cout << "\n Penetration << 203 210 << G4BestUnit(rmsPene, "Length") << G << 211 << 212 // 204 // 213 205 214 // Output file << 206 //output file 215 FILE* myFile; << 207 FILE *myFile; 216 myFile = fopen("range.txt", "a"); << 208 myFile = fopen ("range.txt","a"); 217 fprintf(myFile, "%e %e %e %e %e %e %e\n", fE << 209 fprintf (myFile, "%e %e %e %e %e %e %e\n", 218 fProjRange / nm, rmsProj / nm, fPene << 210 fEkin/eV, 219 fclose(myFile); << 211 fTrackLen/nm, >> 212 rmsTrack/nm, >> 213 fProjRange/nm, >> 214 rmsProj/nm, >> 215 fPenetration/nm, >> 216 rmsPene/nm >> 217 ); >> 218 fclose (myFile); 220 219 221 // Reset default formats << 220 // reset default formats 222 G4cout.setf(mode, std::ios::floatfield); << 221 G4cout.setf(mode,std::ios::floatfield); 223 G4cout.precision(prec); 222 G4cout.precision(prec); >> 223 224 } 224 } >> 225 >> 226 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 225 227