| 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
27 #include "G4AdjointeIonisationModel.hh"
28
29 #include "G4AdjointElectron.hh"
30 #include "G4Electron.hh"
31 #include "G4ParticleChange.hh"
32 #include "G4PhysicalConstants.hh"
33 #include "G4TrackStatus.hh"
34
35 //////////////////////////////////////////////
36 G4AdjointeIonisationModel::G4AdjointeIonisatio
37 : G4VEmAdjointModel("Inv_eIon_model")
38
39 {
40 fUseMatrix = true;
41 fUseMatrixPerElement = true;
42 fApplyCutInRange = true;
43 fOneMatrixForAllElements = true;
44
45 fAdjEquivDirectPrimPart = G4AdjointElectro
46 fAdjEquivDirectSecondPart = G4AdjointElectro
47 fDirectPrimaryPart = G4Electron::Elec
48 fSecondPartSameType = true;
49 }
50
51 //////////////////////////////////////////////
52 G4AdjointeIonisationModel::~G4AdjointeIonisati
53
54 //////////////////////////////////////////////
55 void G4AdjointeIonisationModel::SampleSecondar
56 const G4Track& aTrack, G4bool IsScatProjToPr
57 G4ParticleChange* fParticleChange)
58 {
59 const G4DynamicParticle* theAdjointPrimary =
60
61 // Elastic inverse scattering
62 G4double adjointPrimKinEnergy = theAdjointPr
63 G4double adjointPrimP = theAdjointPr
64
65 if(adjointPrimKinEnergy > GetHighEnergyLimit
66 {
67 return;
68 }
69
70 // Sample secondary energy
71 G4double projectileKinEnergy;
72 if(!fWithRapidSampling)
73 { // used by default
74 projectileKinEnergy =
75 SampleAdjSecEnergyFromCSMatrix(adjointPr
76
77 // Caution!!! this weight correction shoul
78 CorrectPostStepWeight(fParticleChange, aTr
79 adjointPrimKinEnergy
80 IsScatProjToProj);
81 }
82 else
83 { // only for testing
84 G4double Emin, Emax;
85 if(IsScatProjToProj)
86 {
87 Emin = GetSecondAdjEnergyMinForScatProjT
88
89 Emax = GetSecondAdjEnergyMaxForScatProjT
90 }
91 else
92 {
93 Emin = GetSecondAdjEnergyMinForProdToPro
94 Emax = GetSecondAdjEnergyMaxForProdToPro
95 }
96 projectileKinEnergy = Emin * std::pow(Emax
97
98 fLastCS = fLastAdjointCSForScatProjToProj;
99 if(!IsScatProjToProj)
100 fLastCS = fLastAdjointCSForProdToProj;
101
102 G4double new_weight = aTrack.GetWeight();
103 G4double used_diffCS =
104 fLastCS * std::log(Emax / Emin) / projec
105 G4double needed_diffCS = adjointPrimKinEne
106 if(!IsScatProjToProj)
107 needed_diffCS *= DiffCrossSectionPerVolu
108 fCurrentMaterial, projectileKinEnergy,
109 else
110 needed_diffCS *= DiffCrossSectionPerVolu
111 fCurrentMaterial, projectileKinEnergy,
112 new_weight *= needed_diffCS / used_diffCS;
113 fParticleChange->SetParentWeightByProcess(
114 fParticleChange->SetSecondaryWeightByProce
115 fParticleChange->ProposeParentWeight(new_w
116 }
117
118 // Kinematic:
119 // we consider a two body elastic scattering
120 // the projectile knock on an e- at rest and
121 G4double projectileM0 = fAdjEquivDi
122 G4double projectileTotalEnergy = projectileM
123 G4double projectileP2 =
124 projectileTotalEnergy * projectileTotalEne
125
126 // Companion
127 G4double companionM0 = fAdjEquivDirectPrimPa
128 if(IsScatProjToProj)
129 {
130 companionM0 = fAdjEquivDirectSecondPart->G
131 }
132 G4double companionTotalEnergy =
133 companionM0 + projectileKinEnergy - adjoin
134 G4double companionP2 =
135 companionTotalEnergy * companionTotalEnerg
136
137 // Projectile momentum
138 G4double P_parallel =
139 (adjointPrimP * adjointPrimP + projectileP
140 (2. * adjointPrimP);
141 G4double P_perp = std::sqrt(projectileP2 - P
142 G4ThreeVector dir_parallel = theAdjointPrima
143 G4double phi = G4UniformRand()
144 G4ThreeVector projectileMomentum =
145 G4ThreeVector(P_perp * std::cos(phi), P_pe
146 projectileMomentum.rotateUz(dir_parallel);
147
148 if(!IsScatProjToProj)
149 { // kill the primary and add a secondary
150 fParticleChange->ProposeTrackStatus(fStopA
151 fParticleChange->AddSecondary(
152 new G4DynamicParticle(fAdjEquivDirectPri
153 }
154 else
155 {
156 fParticleChange->ProposeEnergy(projectileK
157 fParticleChange->ProposeMomentumDirection(
158 }
159 }
160
161 //////////////////////////////////////////////
162 // The implementation here is correct for ener
163 // photoelectric and compton scattering the me
164 G4double G4AdjointeIonisationModel::DiffCrossS
165 G4double kinEnergyProj, G4double kinEnergyPr
166 {
167 G4double dSigmadEprod = 0.;
168 G4double Emax_proj = GetSecondAdjEnergyMa
169 G4double Emin_proj = GetSecondAdjEnergyMi
170
171 // the produced particle should have a kinet
172 // projectile
173 if(kinEnergyProj > Emin_proj && kinEnergyPro
174 {
175 dSigmadEprod = Z * DiffCrossSectionMoller(
176 }
177 return dSigmadEprod;
178 }
179
180 //////////////////////////////////////////////
181 G4double G4AdjointeIonisationModel::DiffCrossS
182 G4double kinEnergyProj, G4double kinEnergyPr
183 {
184 // G4double energy = kinEnergyProj + electro
185 G4double x = kinEnergyProd / kinEnergyP
186 G4double gam = (kinEnergyProj + electron_
187 G4double gamma2 = gam * gam;
188 G4double beta2 = 1.0 - 1.0 / gamma2;
189
190 G4double gg = (2.0 * gam - 1.0) / gamma2;
191 G4double y = 1.0 - x;
192 G4double fac = twopi_mc2_rcl2 / electron_mas
193 G4double dCS =
194 fac * (1. - gg + ((1.0 - gg * x) / (x * x)
195 (beta2 * (gam - 1.));
196 return dCS / kinEnergyProj;
197 }
198