35 #include "G4ProcessManager.hh" 36 #include "G4UnitsTable.hh" 37 #include "G4EmCalculator.hh" 38 #include "G4Electron.hh" 39 #include "G4SystemOfUnits.hh" 46 :detector(det), primary(kin)
59 G4int prec = G4cout.precision(6);
63 ->GetParticleDefinition();
64 G4String partName = particle->GetParticleName();
65 G4double charge = particle->GetPDGCharge();
70 G4String matName = material->GetName();
71 G4double density = material->GetDensity();
72 G4double radl = material->GetRadlen();
74 G4cout <<
"\n " << partName <<
" (" 75 << G4BestUnit(energy,
"Energy") <<
") in " 76 << material->GetName() <<
" (density: " 77 << G4BestUnit(density,
"Volumic Mass") <<
"; radiation length: " 78 << G4BestUnit(radl,
"Length") <<
")" << G4endl;
83 G4cout <<
"\n Range cuts : \t gamma " 84 << std::setw(8) << G4BestUnit(
rangeCut[0],
"Length")
85 <<
"\t e- " << std::setw(8) << G4BestUnit(
rangeCut[1],
"Length");
86 G4cout <<
"\n Energy cuts : \t gamma " 87 << std::setw(8) << G4BestUnit(
energyCut[0],
"Energy")
88 <<
"\t e- " << std::setw(8) << G4BestUnit(
energyCut[1],
"Energy")
93 G4ProcessVector* plist = particle->GetProcessManager()->GetProcessList();
96 std::vector<G4String> emName;
97 std::vector<G4double> enerCut;
98 size_t length = plist->size();
99 for (
size_t j=0; j<length; j++) {
100 procName = (*plist)[j]->GetProcessName();
102 if ((procName ==
"eBrem")||(procName ==
"muBrems")) cut =
energyCut[0];
103 if (((*plist)[j]->GetProcessType() == fElectromagnetic) &&
104 (procName !=
"msc")) {
105 emName.push_back(procName);
106 enerCut.push_back(cut);
111 G4cout <<
"\n processes : ";
112 for (
size_t j=0; j<emName.size();j++)
113 G4cout <<
"\t" << std::setw(13) << emName[j] <<
"\t";
114 G4cout <<
"\t" << std::setw(13) <<
"total";
117 G4EmCalculator emCal;
121 if (material->GetNumberOfElements() == 1) {
122 G4double
Z = material->GetZ();
123 G4double A = material->GetA();
125 std::vector<G4double> sigma0;
126 G4double sig, sigtot = 0.;
128 for (
size_t j=0; j<emName.size();j++) {
129 sig = emCal.ComputeCrossSectionPerAtom
130 (energy,particle,emName[j],Z,A,enerCut[j]);
132 sigma0.push_back(sig);
134 sigma0.push_back(sigtot);
136 G4cout <<
"\n \n cross section per atom : ";
137 for (
size_t j=0; j<sigma0.size();j++) {
138 G4cout <<
"\t" << std::setw(13) << G4BestUnit(sigma0[j],
"Surface");
144 std::vector<G4double> sigma1;
145 std::vector<G4double> sigma2;
146 G4double Sig, Sigtot = 0.;
148 for (
size_t j=0; j<emName.size();j++) {
149 Sig = emCal.GetCrossSectionPerVolume(energy,particle,emName[j],material);
150 if (Sig == 0.) Sig = emCal.ComputeCrossSectionPerVolume
151 (energy,particle,emName[j],material,enerCut[j]);
153 sigma1.push_back(Sig);
154 sigma2.push_back(Sig/density);
156 sigma1.push_back(Sigtot);
157 sigma2.push_back(Sigtot/density);
160 G4cout <<
"\n \n cross section per volume : ";
161 for (
size_t j=0; j<sigma1.size();j++) {
162 G4cout <<
"\t" << std::setw(13) << sigma1[j]*cm <<
" cm^-1";
165 G4cout <<
"\n cross section per mass : ";
166 for (
size_t j=0; j<sigma2.size();j++) {
167 G4cout <<
"\t" << std::setw(13) << G4BestUnit(sigma2[j],
"Surface/Mass");
174 G4cout <<
"\n \n mean free path : ";
175 for (
size_t j=0; j<sigma1.size();j++) {
177 if (sigma1[j] > 0.) lambda = 1/sigma1[j];
178 G4cout <<
"\t" << std::setw(13) << G4BestUnit( lambda,
"Length");
182 G4cout <<
"\n (g/cm2) : ";
183 for (
size_t j=0; j<sigma2.size();j++) {
185 if (sigma2[j] > 0.) lambda = 1/sigma2[j];
186 G4cout <<
"\t" << std::setw(13) << G4BestUnit( lambda,
"Mass/Surface");
191 G4cout.precision(prec);
192 G4cout <<
"\n-------------------------------------------------------------\n" 198 std::vector<G4double> dedx1;
199 std::vector<G4double> dedx2;
200 G4double dedx, dedxtot = 0.;
202 for (
size_t j=0; j<emName.size();j++) {
203 dedx = emCal.ComputeDEDX(energy,particle,emName[j],material,enerCut[j]);
204 dedx1.push_back(dedx);
205 dedx2.push_back(dedx/density);
207 dedxtot = emCal.GetDEDX(energy,particle,material);
208 dedx1.push_back(dedxtot);
209 dedx2.push_back(dedxtot/density);
212 G4cout <<
"\n \n restricted dE/dx : ";
213 for (
size_t j=0; j<sigma1.size();j++) {
214 G4cout <<
"\t" << std::setw(13) << G4BestUnit(dedx1[j],
"Energy/Length");
217 G4cout <<
"\n (MeV/g/cm2) : ";
218 for (
size_t j=0; j<sigma2.size();j++) {
219 G4cout <<
"\t" << std::setw(13) << G4BestUnit(dedx2[j],
"Energy*Surface/Mass");
223 G4double range1 = emCal.GetRangeFromRestricteDEDX(energy,particle,material);
224 G4double range2 = range1*density;
227 G4double Range1 = emCal.GetCSDARange(energy,particle,material);
228 G4double Range2 = Range1*density;
231 G4cout <<
"\n \n range from restrict dE/dx: " 232 <<
"\t" << std::setw(8) << G4BestUnit(range1,
"Length")
233 <<
" (" << std::setw(8) << G4BestUnit(range2,
"Mass/Surface") <<
")";
235 G4cout <<
"\n range from full dE/dx : " 236 <<
"\t" << std::setw(8) << G4BestUnit(Range1,
"Length")
237 <<
" (" << std::setw(8) << G4BestUnit(Range2,
"Mass/Surface") <<
")";
240 G4double MSmfp1 = emCal.GetMeanFreePath(energy,particle,
"msc",material);
241 G4double MSmfp2 = MSmfp1*density;
244 G4cout <<
"\n \n transport mean free path : " 245 <<
"\t" << std::setw(8) << G4BestUnit(MSmfp1,
"Length")
246 <<
" (" << std::setw(8) << G4BestUnit(MSmfp2,
"Mass/Surface") <<
")";
250 G4cout <<
"\n-------------------------------------------------------------\n";
254 G4cout.precision(prec);
264 #include "G4ProductionCutsTable.hh" 270 G4ProductionCutsTable* theCoupleTable =
271 G4ProductionCutsTable::GetProductionCutsTable();
273 size_t numOfCouples = theCoupleTable->GetTableSize();
274 const G4MaterialCutsCouple* couple = 0;
276 for (
size_t i=0; i<numOfCouples; i++) {
277 couple = theCoupleTable->GetMaterialCutsCouple(i);
282 (*(theCoupleTable->GetRangeCutsVector(idxG4GammaCut)))[index];
284 (*(theCoupleTable->GetRangeCutsVector(idxG4ElectronCut)))[index];
286 (*(theCoupleTable->GetRangeCutsVector(idxG4PositronCut)))[index];
289 (*(theCoupleTable->GetEnergyCutsVector(idxG4GammaCut)))[index];
291 (*(theCoupleTable->GetEnergyCutsVector(idxG4ElectronCut)))[index];
293 (*(theCoupleTable->GetEnergyCutsVector(idxG4PositronCut)))[index];
304 G4EmCalculator emCal;
307 const G4double radl = material->GetRadlen();
308 G4double ekin = 5*MeV;
310 G4double err = 1., errmax = 0.001;
311 G4int iter = 0 , itermax = 10;
312 while (err > errmax && iter < itermax) {
315 emCal.ComputeDEDX(ekin,G4Electron::Electron(),
"eIoni",material);
319 G4cout <<
"\n \n critical energy (Rossi) : " 320 <<
"\t" << std::setw(8) << G4BestUnit(ekin,
"Energy");
323 G4double pdga[2] = { 610*MeV, 710*MeV };
324 G4double pdgb[2] = { 1.24, 0.92 };
327 if (material->GetNumberOfElements() == 1) {
329 if (material->GetState() == kStateGas) istat = 1;
330 G4double Zeff = material->GetZ() + pdgb[istat];
331 EcPdg = pdga[istat]/Zeff;
332 G4cout <<
"\t\t\t (from Pdg formula : " 333 << std::setw(8) << G4BestUnit(EcPdg,
"Energy") <<
")";
336 const G4double Es = 21.2052*MeV;
337 G4double rMolier1 = Es/ekin, rMolier2 = rMolier1*radl;
338 G4cout <<
"\n Moliere radius : " 339 <<
"\t" << std::setw(8) << rMolier1 <<
" X0 " 340 <<
"= " << std::setw(8) << G4BestUnit(rMolier2,
"Length");
342 if (material->GetNumberOfElements() == 1) {
343 G4double rMPdg = radl*Es/EcPdg;
344 G4cout <<
"\t (from Pdg formula : " 345 << std::setw(8) << G4BestUnit(rMPdg,
"Length") <<
")";
void EndOfRunAction(const G4Run *)
DetectorConstruction * detector
constexpr auto abs(T v)
Returns the absolute value of the argument.
G4ParticleGun * GetParticleGun()
void BeginOfRunAction(const G4Run *)
PrimaryGeneratorAction * primary
G4Material * GetMaterial()
RunAction(DetectorConstruction *, PrimaryGeneratorAction *)