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/TestEm0/src/RunActio 26 /// \file electromagnetic/TestEm0/src/RunAction.cc 27 /// \brief Implementation of the RunAction cla 27 /// \brief Implementation of the RunAction 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 "RunAction.hh" 33 #include "RunAction.hh" 34 << 35 #include "DetectorConstruction.hh" 34 #include "DetectorConstruction.hh" 36 #include "PrimaryGeneratorAction.hh" 35 #include "PrimaryGeneratorAction.hh" 37 36 38 #include "G4Electron.hh" << 37 #include "G4Run.hh" 39 #include "G4EmCalculator.hh" << 38 #include "G4ProcessManager.hh" 40 #include "G4LossTableManager.hh" 39 #include "G4LossTableManager.hh" >> 40 #include "G4UnitsTable.hh" >> 41 #include "G4EmCalculator.hh" 41 #include "G4PhysicalConstants.hh" 42 #include "G4PhysicalConstants.hh" 42 #include "G4Positron.hh" << 43 #include "G4ProcessManager.hh" << 44 #include "G4Run.hh" << 45 #include "G4SystemOfUnits.hh" 43 #include "G4SystemOfUnits.hh" 46 #include "G4UnitsTable.hh" << 44 #include "G4Electron.hh" >> 45 #include "G4Positron.hh" 47 46 48 #include <vector> 47 #include <vector> 49 48 50 //....oooOO0OOooo........oooOO0OOooo........oo 49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 51 50 52 RunAction::RunAction(DetectorConstruction* det 51 RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* kin) 53 : fDetector(det), fPrimary(kin) << 52 : G4UserRunAction(), fDetector(det), fPrimary(kin) 54 {} << 53 { } >> 54 >> 55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 56 >> 57 RunAction::~RunAction() >> 58 { } >> 59 55 //....oooOO0OOooo........oooOO0OOooo........oo 60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 56 61 57 void RunAction::BeginOfRunAction(const G4Run*) 62 void RunAction::BeginOfRunAction(const G4Run*) 58 { 63 { 59 // set precision for printing << 64 //set precision for printing 60 G4int prec = G4cout.precision(6); 65 G4int prec = G4cout.precision(6); 61 << 66 62 // instanciate EmCalculator << 67 //instanciate EmCalculator 63 G4EmCalculator emCal; 68 G4EmCalculator emCal; 64 // emCal.SetVerbose(2); 69 // emCal.SetVerbose(2); 65 << 70 66 // get particle << 71 // get particle 67 G4ParticleDefinition* particle = fPrimary->G << 72 G4ParticleDefinition* particle = fPrimary->GetParticleGun() >> 73 ->GetParticleDefinition(); 68 G4String partName = particle->GetParticleNam 74 G4String partName = particle->GetParticleName(); 69 G4double charge = particle->GetPDGCharge(); << 75 G4double charge = particle->GetPDGCharge(); 70 G4double energy = fPrimary->GetParticleGun() << 76 G4double energy = fPrimary->GetParticleGun()->GetParticleEnergy(); 71 << 77 72 // get material 78 // get material 73 const G4Material* material = fDetector->GetM 79 const G4Material* material = fDetector->GetMaterial(); 74 G4String matName = material->GetName(); << 80 G4String matName = material->GetName(); 75 G4double density = material->GetDensity(); << 81 G4double density = material->GetDensity(); 76 G4double radl = material->GetRadlen(); << 82 G4double radl = material->GetRadlen(); 77 << 83 78 G4cout << "\n " << partName << " (" << G4Bes << 84 G4cout << "\n " << partName << " (" 79 << material->GetName() << " (density: << 85 << G4BestUnit(energy,"Energy") << ") in " 80 << "; radiation length: " << G4Best << 86 << material->GetName() << " (density: " >> 87 << G4BestUnit(density,"Volumic Mass") << "; radiation length: " >> 88 << G4BestUnit(radl, "Length") << ")" << G4endl; 81 89 82 // get cuts << 90 // get cuts 83 GetCuts(); 91 GetCuts(); 84 if (charge != 0.) { 92 if (charge != 0.) { 85 G4cout << "\n Range cuts: \t gamma " << s << 93 G4cout << "\n Range cuts: \t gamma " 86 << "\t e- " << std::setw(12) << G4B << 94 << std::setw(12) << G4BestUnit(fRangeCut[0],"Length") 87 G4cout << "\n Energy cuts: \t gamma " << s << 95 << "\t e- " << std::setw(12) << G4BestUnit(fRangeCut[1],"Length"); 88 << "\t e- " << std::setw(12) << G4B << 96 G4cout << "\n Energy cuts: \t gamma " 89 } << 97 << std::setw(12) << G4BestUnit(fEnergyCut[0],"Energy") 90 << 98 << "\t e- " << std::setw(12) << G4BestUnit(fEnergyCut[1],"Energy") >> 99 << G4endl; >> 100 } >> 101 91 // max energy transfert 102 // max energy transfert 92 if (charge != 0.) { 103 if (charge != 0.) { 93 G4double Mass_c2 = particle->GetPDGMass(); << 104 G4double Mass_c2 = particle->GetPDGMass(); 94 G4double moverM = electron_mass_c2 / Mass_ << 105 G4double moverM = electron_mass_c2/Mass_c2; 95 G4double gamM1 = energy / Mass_c2, gam = g << 106 G4double gamM1 = energy/Mass_c2, gam = gamM1 + 1., gamP1 = gam + 1.; 96 G4double Tmax = energy; << 107 G4double Tmax = energy; 97 if (particle == G4Electron::Electron()) { << 108 if(particle == G4Electron::Electron()) { 98 Tmax *= 0.5; << 109 Tmax *= 0.5; 99 } << 110 } else if(particle != G4Positron::Positron()) { 100 else if (particle != G4Positron::Positron( << 111 Tmax = (2*electron_mass_c2*gamM1*gamP1)/(1.+2*gam*moverM+moverM*moverM); 101 Tmax = (2 * electron_mass_c2 * gamM1 * g << 112 } 102 } << 113 G4double range = emCal.GetCSDARange(Tmax,G4Electron::Electron(),material); 103 G4double range = emCal.GetCSDARange(Tmax, << 114 104 << 115 G4cout << "\n Max_energy _transferable : " << G4BestUnit(Tmax,"Energy") 105 G4cout << "\n Max_energy _transferable : << 116 << " (" << G4BestUnit(range,"Length") << ")" << G4endl; 106 << G4BestUnit(range, "Length") << " << 107 } 117 } 108 << 118 109 // get processList and extract EM processes 119 // get processList and extract EM processes (but not MultipleScattering) 110 G4ProcessVector* plist = particle->GetProces 120 G4ProcessVector* plist = particle->GetProcessManager()->GetProcessList(); 111 G4String procName; 121 G4String procName; 112 G4double cut; 122 G4double cut; 113 std::vector<G4String> emName; 123 std::vector<G4String> emName; 114 std::vector<G4double> enerCut; 124 std::vector<G4double> enerCut; 115 size_t length = plist->size(); 125 size_t length = plist->size(); 116 for (size_t j = 0; j < length; j++) { << 126 for (size_t j=0; j<length; j++) { 117 procName = (*plist)[j]->GetProcessName(); << 127 procName = (*plist)[j]->GetProcessName(); 118 cut = fEnergyCut[1]; << 128 cut = fEnergyCut[1]; 119 if ((procName == "eBrem") || (procName == << 129 if ((procName == "eBrem")||(procName == "muBrems")) cut = fEnergyCut[0]; 120 if (((*plist)[j]->GetProcessType() == fEle << 130 if (((*plist)[j]->GetProcessType() == fElectromagnetic) && 121 emName.push_back(procName); << 131 (procName != "msc")) { 122 enerCut.push_back(cut); << 132 emName.push_back(procName); 123 } << 133 enerCut.push_back(cut); >> 134 } 124 } 135 } 125 << 136 126 // write html documentation, if requested 137 // write html documentation, if requested 127 char* htmlDocName = std::getenv("G4PhysListN << 138 char* htmlDocName = std::getenv("G4PhysListName"); // file name 128 char* htmlDocDir = std::getenv("G4PhysListDo << 139 char* htmlDocDir = std::getenv("G4PhysListDocDir"); // directory 129 if (htmlDocName && htmlDocDir) { 140 if (htmlDocName && htmlDocDir) { 130 G4LossTableManager::Instance()->DumpHtml() 141 G4LossTableManager::Instance()->DumpHtml(); 131 } 142 } 132 << 143 133 // print list of processes 144 // print list of processes 134 G4cout << "\n processes : "; 145 G4cout << "\n processes : "; 135 for (size_t j = 0; j < emName.size(); ++j) { << 146 for (size_t j=0; j<emName.size(); ++j) { 136 G4cout << "\t" << std::setw(14) << emName[ 147 G4cout << "\t" << std::setw(14) << emName[j] << "\t"; 137 } 148 } 138 G4cout << "\t" << std::setw(14) << "total"; << 149 G4cout << "\t" << std::setw(14) <<"total"; 139 << 150 140 // compute cross section per atom (only for << 151 //compute cross section per atom (only for single material) 141 if (material->GetNumberOfElements() == 1) { 152 if (material->GetNumberOfElements() == 1) { 142 G4double Z = material->GetZ(); 153 G4double Z = material->GetZ(); 143 G4double A = material->GetA(); 154 G4double A = material->GetA(); 144 << 155 145 std::vector<G4double> sigma0; 156 std::vector<G4double> sigma0; 146 G4double sig, sigtot = 0.; 157 G4double sig, sigtot = 0.; 147 158 148 for (size_t j = 0; j < emName.size(); j++) << 159 for (size_t j=0; j<emName.size();j++) { 149 sig = emCal.ComputeCrossSectionPerAtom(e << 160 sig = emCal.ComputeCrossSectionPerAtom 150 sigtot += sig; << 161 (energy,particle,emName[j],Z,A,enerCut[j]); 151 sigma0.push_back(sig); << 162 sigtot += sig; >> 163 sigma0.push_back(sig); 152 } 164 } 153 sigma0.push_back(sigtot); 165 sigma0.push_back(sigtot); 154 166 155 G4cout << "\n \n cross section per atom 167 G4cout << "\n \n cross section per atom : "; 156 for (size_t j = 0; j < sigma0.size(); ++j) << 168 for (size_t j=0; j<sigma0.size(); ++j) { 157 G4cout << "\t" << std::setw(9) << G4Best 169 G4cout << "\t" << std::setw(9) << G4BestUnit(sigma0[j], "Surface"); 158 } 170 } 159 G4cout << G4endl; 171 G4cout << G4endl; 160 } 172 } 161 << 173 162 // get cross section per volume << 174 //get cross section per volume 163 std::vector<G4double> sigma0; 175 std::vector<G4double> sigma0; 164 std::vector<G4double> sigma1; 176 std::vector<G4double> sigma1; 165 std::vector<G4double> sigma2; 177 std::vector<G4double> sigma2; 166 G4double Sig, SigtotComp = 0., Sigtot = 0.; 178 G4double Sig, SigtotComp = 0., Sigtot = 0.; 167 179 168 for (size_t j = 0; j < emName.size(); ++j) { << 180 for (size_t j=0; j<emName.size(); ++j) { 169 Sig = emCal.ComputeCrossSectionPerVolume(e << 181 Sig = emCal.ComputeCrossSectionPerVolume 170 SigtotComp += Sig; << 182 (energy,particle,emName[j],material,enerCut[j]); >> 183 SigtotComp += Sig; 171 sigma0.push_back(Sig); 184 sigma0.push_back(Sig); 172 Sig = emCal.GetCrossSectionPerVolume(energ << 185 Sig = emCal.GetCrossSectionPerVolume(energy,particle,emName[j],material); 173 Sigtot += Sig; << 186 Sigtot += Sig; 174 sigma1.push_back(Sig); 187 sigma1.push_back(Sig); 175 sigma2.push_back(Sig / density); << 188 sigma2.push_back(Sig/density); 176 } 189 } 177 sigma0.push_back(SigtotComp); 190 sigma0.push_back(SigtotComp); 178 sigma1.push_back(Sigtot); 191 sigma1.push_back(Sigtot); 179 sigma2.push_back(Sigtot / density); << 192 sigma2.push_back(Sigtot/density); 180 << 193 181 // print cross sections << 194 //print cross sections 182 G4cout << "\n compCrossSectionPerVolume: "; << 195 G4cout << "\n \n compCrossSectionPerVolume: "; 183 for (size_t j = 0; j < sigma0.size(); ++j) { << 196 for (size_t j=0; j<sigma0.size(); ++j) { 184 G4cout << "\t" << std::setw(9) << sigma0[j << 197 G4cout << "\t" << std::setw(9) << sigma0[j]*cm << std::setw(6) << " cm^-1"; 185 } 198 } 186 G4cout << "\n cross section per volume : "; 199 G4cout << "\n cross section per volume : "; 187 for (size_t j = 0; j < sigma1.size(); ++j) { << 200 for (size_t j=0; j<sigma1.size(); ++j) { 188 G4cout << "\t" << std::setw(9) << sigma1[j << 201 G4cout << "\t" << std::setw(9) << sigma1[j]*cm << std::setw(6) << " cm^-1"; 189 } 202 } 190 << 203 191 G4cout << "\n cross section per mass : "; 204 G4cout << "\n cross section per mass : "; 192 for (size_t j = 0; j < sigma2.size(); ++j) { << 205 for (size_t j=0; j<sigma2.size(); ++j) { 193 G4cout << "\t" << std::setw(9) << G4BestUn << 206 G4cout << "\t" << std::setw(9) 194 } << 207 << G4BestUnit(sigma2[j], "Surface/Mass"); 195 << 208 } 196 // print mean free path << 209 197 << 210 //print mean free path >> 211 198 G4double lambda; 212 G4double lambda; 199 << 213 200 G4cout << "\n \n mean free path : 214 G4cout << "\n \n mean free path : "; 201 for (size_t j = 0; j < sigma1.size(); ++j) { << 215 for (size_t j=0; j<sigma1.size(); ++j) { 202 lambda = DBL_MAX; << 216 lambda = DBL_MAX; 203 if (sigma1[j] > 0.) lambda = 1 / sigma1[j] << 217 if (sigma1[j] > 0.) lambda = 1/sigma1[j]; 204 G4cout << "\t" << std::setw(9) << G4BestUn << 218 G4cout << "\t" << std::setw(9) << G4BestUnit( lambda, "Length"); 205 } << 219 } 206 << 220 207 // mean free path (g/cm2) << 221 //mean free path (g/cm2) 208 G4cout << "\n (g/cm2) : "; << 222 G4cout << "\n (g/cm2) : "; 209 for (size_t j = 0; j < sigma2.size(); ++j) { << 223 for (size_t j=0; j<sigma2.size(); ++j) { 210 lambda = DBL_MAX; << 224 lambda = DBL_MAX; 211 if (sigma2[j] > 0.) lambda = 1 / sigma2[j] << 225 if (sigma2[j] > 0.) lambda = 1/sigma2[j]; 212 G4cout << "\t" << std::setw(9) << G4BestUn << 226 G4cout << "\t" << std::setw(9) << G4BestUnit( lambda, "Mass/Surface"); 213 } 227 } 214 G4cout << G4endl; 228 G4cout << G4endl; 215 << 229 216 if (charge == 0.) { 230 if (charge == 0.) { 217 G4cout.precision(prec); 231 G4cout.precision(prec); 218 G4cout << "\n----------------------------- << 232 G4cout << "\n-----------------------------------------------------------\n" >> 233 << G4endl; 219 return; 234 return; 220 } 235 } 221 << 236 222 // get stopping power << 237 //get stopping power 223 std::vector<G4double> dedx1; 238 std::vector<G4double> dedx1; 224 std::vector<G4double> dedx2; << 239 std::vector<G4double> dedx2; 225 G4double dedx, dedxtot = 0.; 240 G4double dedx, dedxtot = 0.; 226 size_t nproc = emName.size(); 241 size_t nproc = emName.size(); 227 242 228 for (size_t j = 0; j < nproc; ++j) { << 243 for (size_t j=0; j<nproc; ++j) { 229 dedx = emCal.ComputeDEDX(energy, particle, << 244 dedx = emCal.ComputeDEDX(energy,particle,emName[j],material,enerCut[j]); 230 dedxtot += dedx; 245 dedxtot += dedx; 231 dedx1.push_back(dedx); 246 dedx1.push_back(dedx); 232 dedx2.push_back(dedx / density); << 247 dedx2.push_back(dedx/density); 233 } 248 } 234 dedx1.push_back(dedxtot); 249 dedx1.push_back(dedxtot); 235 dedx2.push_back(dedxtot / density); << 250 dedx2.push_back(dedxtot/density); 236 << 251 237 // print stopping power << 252 //print stopping power 238 G4cout << "\n \n restricted dE/dx : 253 G4cout << "\n \n restricted dE/dx : "; 239 for (size_t j = 0; j <= nproc; ++j) { << 254 for (size_t j=0; j<=nproc; ++j) { 240 G4cout << "\t" << std::setw(9) << G4BestUn << 255 G4cout << "\t" << std::setw(14) >> 256 << G4BestUnit(dedx1[j],"Energy/Length"); 241 } 257 } 242 << 258 243 G4cout << "\n (MeV/g/cm2) : "; 259 G4cout << "\n (MeV/g/cm2) : "; 244 for (size_t j = 0; j <= nproc; ++j) { << 260 for (size_t j=0; j<=nproc; ++j) { 245 G4cout << "\t" << std::setw(9) << G4BestUn << 261 G4cout << "\t" << std::setw(14) >> 262 << G4BestUnit(dedx2[j],"Energy*Surface/Mass"); 246 } 263 } 247 dedxtot = 0.; 264 dedxtot = 0.; 248 265 249 for (size_t j = 0; j < nproc; ++j) { << 266 for (size_t j=0; j<nproc; ++j) { 250 dedx = emCal.ComputeDEDX(energy, particle, << 267 dedx = emCal.ComputeDEDX(energy,particle,emName[j],material,energy); 251 dedxtot += dedx; 268 dedxtot += dedx; 252 dedx1[j] = dedx; 269 dedx1[j] = dedx; 253 dedx2[j] = dedx / density; << 270 dedx2[j] = dedx/density; 254 } 271 } 255 dedx1[nproc] = dedxtot; 272 dedx1[nproc] = dedxtot; 256 dedx2[nproc] = dedxtot / density; << 273 dedx2[nproc] = dedxtot/density; 257 << 274 258 // print stopping power << 275 //print stopping power 259 G4cout << "\n \n unrestricted dE/dx : 276 G4cout << "\n \n unrestricted dE/dx : "; 260 for (size_t j = 0; j <= nproc; ++j) { << 277 for (size_t j=0; j<=nproc; ++j) { 261 G4cout << "\t" << std::setw(9) << G4BestUn << 278 G4cout << "\t" << std::setw(14) << G4BestUnit(dedx1[j],"Energy/Length"); 262 } 279 } 263 << 280 264 G4cout << "\n (MeV/g/cm2) : "; 281 G4cout << "\n (MeV/g/cm2) : "; 265 for (size_t j = 0; j <= nproc; ++j) { << 282 for (size_t j=0; j<=nproc; ++j) { 266 G4cout << "\t" << std::setw(9) << G4BestUn << 283 G4cout << "\t" << std::setw(14) 267 } << 284 << G4BestUnit(dedx2[j],"Energy*Surface/Mass"); 268 << 285 } 269 // get range from restricted dedx << 286 270 G4double range1 = emCal.GetRangeFromRestrict << 287 //get range from restricted dedx 271 G4double range2 = range1 * density; << 288 G4double range1 = emCal.GetRangeFromRestricteDEDX(energy,particle,material); 272 << 289 G4double range2 = range1*density; 273 // print range << 290 274 G4cout << "\n \n range from restrict dE/dx: << 291 //print range 275 << "\t" << std::setw(9) << G4BestUnit << 292 G4cout << "\n \n range from restrict dE/dx: " 276 << G4BestUnit(range2, "Mass/Surface") << 293 << "\t" << std::setw(9) << G4BestUnit(range1,"Length") 277 << 294 << " (" << std::setw(9) << G4BestUnit(range2,"Mass/Surface") << ")"; 278 // get range from full dedx << 295 >> 296 //get range from full dedx 279 G4double EmaxTable = G4EmParameters::Instanc 297 G4double EmaxTable = G4EmParameters::Instance()->MaxEnergyForCSDARange(); 280 if (energy < EmaxTable) { << 298 if(energy < EmaxTable) { 281 G4double Range1 = emCal.GetCSDARange(energ << 299 G4double Range1 = emCal.GetCSDARange(energy,particle,material); 282 G4double Range2 = Range1 * density; << 300 G4double Range2 = Range1*density; 283 << 301 284 G4cout << "\n range from full dE/dx : << 302 G4cout << "\n range from full dE/dx : " 285 << "\t" << std::setw(9) << G4BestUn << 303 << "\t" << std::setw(9) << G4BestUnit(Range1,"Length") 286 << G4BestUnit(Range2, "Mass/Surface << 304 << " (" << std::setw(9) << G4BestUnit(Range2,"Mass/Surface") << ")"; 287 } << 305 } 288 << 306 289 // get transport mean free path (for multipl << 307 //get transport mean free path (for multiple scattering) 290 G4double MSmfp1 = emCal.GetMeanFreePath(ener << 308 G4double MSmfp1 = emCal.GetMeanFreePath(energy,particle,"msc",material); 291 G4double MSmfp2 = MSmfp1 * density; << 309 G4double MSmfp2 = MSmfp1*density; 292 << 310 293 // print transport mean free path << 311 //print transport mean free path 294 G4cout << "\n \n transport mean free path : << 312 G4cout << "\n \n transport mean free path : " 295 << "\t" << std::setw(9) << G4BestUnit << 313 << "\t" << std::setw(9) << G4BestUnit(MSmfp1,"Length") 296 << G4BestUnit(MSmfp2, "Mass/Surface") << 314 << " (" << std::setw(9) << G4BestUnit(MSmfp2,"Mass/Surface") << ")"; 297 315 298 if (particle == G4Electron::Electron()) Crit 316 if (particle == G4Electron::Electron()) CriticalEnergy(); 299 << 317 300 G4cout << "\n------------------------------- 318 G4cout << "\n-------------------------------------------------------------\n"; 301 G4cout << G4endl; 319 G4cout << G4endl; 302 << 320 303 // reset default precision << 321 // reset default precision 304 G4cout.precision(prec); << 322 G4cout.precision(prec); 305 } 323 } 306 324 307 //....oooOO0OOooo........oooOO0OOooo........oo 325 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 308 326 309 void RunAction::EndOfRunAction(const G4Run*) { << 327 void RunAction::EndOfRunAction(const G4Run* ) >> 328 { } 310 329 311 //....oooOO0OOooo........oooOO0OOooo........oo 330 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 312 331 313 #include "G4ProductionCutsTable.hh" 332 #include "G4ProductionCutsTable.hh" 314 333 315 //....oooOO0OOooo........oooOO0OOooo........oo 334 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 316 335 317 void RunAction::GetCuts() 336 void RunAction::GetCuts() 318 { << 337 { 319 G4ProductionCutsTable* theCoupleTable = G4Pr << 338 G4ProductionCutsTable* theCoupleTable = 320 << 339 G4ProductionCutsTable::GetProductionCutsTable(); >> 340 321 size_t numOfCouples = theCoupleTable->GetTab 341 size_t numOfCouples = theCoupleTable->GetTableSize(); 322 const G4MaterialCutsCouple* couple = 0; 342 const G4MaterialCutsCouple* couple = 0; 323 G4int index = 0; 343 G4int index = 0; 324 for (size_t i = 0; i < numOfCouples; i++) { << 344 for (size_t i=0; i<numOfCouples; i++) { 325 couple = theCoupleTable->GetMaterialCutsCo << 345 couple = theCoupleTable->GetMaterialCutsCouple(i); 326 if (couple->GetMaterial() == fDetector->Ge << 346 if (couple->GetMaterial() == fDetector->GetMaterial()) {index = i; break;} 327 index = i; << 347 } 328 break; << 348 329 } << 349 fRangeCut[0] = 330 } << 350 (*(theCoupleTable->GetRangeCutsVector(idxG4GammaCut)))[index]; >> 351 fRangeCut[1] = >> 352 (*(theCoupleTable->GetRangeCutsVector(idxG4ElectronCut)))[index]; >> 353 fRangeCut[2] = >> 354 (*(theCoupleTable->GetRangeCutsVector(idxG4PositronCut)))[index]; >> 355 >> 356 fEnergyCut[0] = >> 357 (*(theCoupleTable->GetEnergyCutsVector(idxG4GammaCut)))[index]; >> 358 fEnergyCut[1] = >> 359 (*(theCoupleTable->GetEnergyCutsVector(idxG4ElectronCut)))[index]; >> 360 fEnergyCut[2] = >> 361 (*(theCoupleTable->GetEnergyCutsVector(idxG4PositronCut)))[index]; 331 362 332 fRangeCut[0] = (*(theCoupleTable->GetRangeCu << 333 fRangeCut[1] = (*(theCoupleTable->GetRangeCu << 334 fRangeCut[2] = (*(theCoupleTable->GetRangeCu << 335 << 336 fEnergyCut[0] = (*(theCoupleTable->GetEnergy << 337 fEnergyCut[1] = (*(theCoupleTable->GetEnergy << 338 fEnergyCut[2] = (*(theCoupleTable->GetEnergy << 339 } 363 } 340 364 341 //....oooOO0OOooo........oooOO0OOooo........oo 365 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 342 366 343 void RunAction::CriticalEnergy() 367 void RunAction::CriticalEnergy() 344 { 368 { 345 // compute e- critical energy (Rossi definit 369 // compute e- critical energy (Rossi definition) and Moliere radius. 346 // Review of Particle Physics - Eur. Phys. J 370 // Review of Particle Physics - Eur. Phys. J. C3 (1998) page 147 347 // 371 // 348 G4EmCalculator emCal; 372 G4EmCalculator emCal; 349 << 373 350 const G4Material* material = fDetector->GetM 374 const G4Material* material = fDetector->GetMaterial(); 351 const G4double radl = material->GetRadlen(); 375 const G4double radl = material->GetRadlen(); 352 G4double ekin = 5 * MeV; << 376 G4double ekin = 5*MeV; 353 G4double deioni; 377 G4double deioni; 354 G4double err = 1., errmax = 0.001; << 378 G4double err = 1., errmax = 0.001; 355 G4int iter = 0, itermax = 10; << 379 G4int iter = 0 , itermax = 10; 356 while (err > errmax && iter < itermax) { 380 while (err > errmax && iter < itermax) { 357 iter++; << 381 iter++; 358 deioni = radl * emCal.ComputeDEDX(ekin, G4 << 382 deioni = radl* 359 err = std::abs(deioni - ekin) / ekin; << 383 emCal.ComputeDEDX(ekin,G4Electron::Electron(),"eIoni",material); >> 384 err = std::abs(deioni - ekin)/ekin; 360 ekin = deioni; 385 ekin = deioni; 361 } 386 } 362 G4cout << "\n \n critical energy (Rossi) : << 387 G4cout << "\n \n critical energy (Rossi) : " 363 << "\t" << std::setw(8) << G4BestUnit << 388 << "\t" << std::setw(8) << G4BestUnit(ekin,"Energy"); 364 << 389 365 // Pdg formula (only for single material) << 390 //Pdg formula (only for single material) 366 G4double pdga[2] = {610 * MeV, 710 * MeV}; << 391 G4double pdga[2] = { 610*MeV, 710*MeV }; 367 G4double pdgb[2] = {1.24, 0.92}; << 392 G4double pdgb[2] = { 1.24, 0.92 }; 368 G4double EcPdg = 0.; 393 G4double EcPdg = 0.; 369 << 394 370 if (material->GetNumberOfElements() == 1) { 395 if (material->GetNumberOfElements() == 1) { 371 G4int istat = 0; 396 G4int istat = 0; 372 if (material->GetState() == kStateGas) ist << 397 if (material->GetState() == kStateGas) istat = 1; 373 G4double Zeff = material->GetZ() + pdgb[is 398 G4double Zeff = material->GetZ() + pdgb[istat]; 374 EcPdg = pdga[istat] / Zeff; << 399 EcPdg = pdga[istat]/Zeff; 375 G4cout << "\t\t\t (from Pdg formula : " << << 400 G4cout << "\t\t\t (from Pdg formula : " 376 } << 401 << std::setw(8) << G4BestUnit(EcPdg,"Energy") << ")"; 377 << 402 } 378 const G4double Es = 21.2052 * MeV; << 403 379 G4double rMolier1 = Es / ekin, rMolier2 = rM << 404 const G4double Es = 21.2052*MeV; 380 G4cout << "\n Moliere radius : " << 405 G4double rMolier1 = Es/ekin, rMolier2 = rMolier1*radl; 381 << "\t" << std::setw(8) << rMolier1 < << 406 G4cout << "\n Moliere radius : " 382 << "= " << std::setw(8) << G4BestUnit << 407 << "\t" << std::setw(8) << rMolier1 << " X0 " 383 << 408 << "= " << std::setw(8) << G4BestUnit(rMolier2,"Length"); 384 if (material->GetNumberOfElements() == 1) { << 409 385 G4double rMPdg = radl * Es / EcPdg; << 410 if (material->GetNumberOfElements() == 1) { 386 G4cout << "\t (from Pdg formula : " << std << 411 G4double rMPdg = radl*Es/EcPdg; 387 } << 412 G4cout << "\t (from Pdg formula : " >> 413 << std::setw(8) << G4BestUnit(rMPdg,"Length") << ")"; >> 414 } 388 } 415 } 389 416 390 //....oooOO0OOooo........oooOO0OOooo........oo 417 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 391 418