84{
86 G4int idx = fShellIdx;
87
88
89 if(idx < 0 || idx >= nShells) {
90 if(nShells> nprob) {
91 nprob = nShells;
92 prob.resize(nprob,0.0);
93 }
95 for(idx=0; idx<nShells; ++idx) {
98 prob[idx] = sum;
99 }
101 for(idx=0; idx<nShells; ++idx) {
102 if(sum <= prob[idx]) { break; }
103 }
104 }
107
108
109
110
111
114 static const G4int nmax = 100;
115 do {
117
121 G4double e = kinEnergyFinal + ePotEnergy + electron_mass_c2;
122 G4double p = sqrt((e + electron_mass_c2)*(e - electron_mass_c2));
123
127 totEnergy += ePotEnergy;
128 totMomentum = sqrt((totEnergy + electron_mass_c2)
129 *(totEnergy - electron_mass_c2));
130 }
131
132 G4double eTotEnergy = eKinEnergy + electron_mass_c2;
133 G4double eTotMomentum = sqrt(eKinEnergy*(eTotEnergy + electron_mass_c2));
135 G4double sintet = sqrt((1 - costet)*(1 + costet));
136
137 cost = 1.0;
138 if(n >= nmax) {
139
140
141
142
143
144
145
146
147
148
149
150 if(0.0 == bindingEnergy) { isOK = true; }
152 }
153
154 G4double x0 = p*(totMomentum + eTotMomentum*costet);
155
156
157
158
159
160
161 if(x0 > 0.0) {
162 G4double x1 = p*eTotMomentum*sintet;
163 G4double x2 = totEnergy*(eTotEnergy - e) - e*eTotEnergy
164 - totMomentum*eTotMomentum*costet + electron_mass_c2*electron_mass_c2;
166 if(std::abs(y) <= 1.0) {
167 cost = -(x2 + x1*sqrt(1. - y*y))/x0;
168 if(std::abs(cost) <= 1.0) { isOK = true; }
169 else { cost = 1.0; }
170 }
171
172
173
174
175
176
177
178
179
180 }
181
182
183 } while(!isOK);
184
185 G4double sint = sqrt((1 - cost)*(1 + cost));
187
190
192}
G4double G4Log(G4double x)
void set(double x, double y, double z)
Hep3Vector & rotateUz(const Hep3Vector &)
static G4int GetNumberOfElectrons(G4int Z, G4int SubshellNb)
static G4double GetBindingEnergy(G4int Z, G4int SubshellNb)
static G4int GetNumberOfShells(G4int Z)
const G4ThreeVector & GetMomentumDirection() const
const G4ParticleDefinition * GetParticleDefinition() const
G4double GetTotalEnergy() const
G4double GetTotalMomentum() const
G4ThreeVector fLocalDirection
G4double bindingEnergy(G4int A, G4int Z)