124{
125 static const G4int nE=336;
126 static const G4int mL=nE-1;
129 static const G4double lEMi=std::log(EMi);
130 static const G4double lEMa=std::log(EMa);
131 static const G4double dlnE=(lEMa-lEMi)/mL;
132 static const G4double alop=1./137.036/3.14159265;
133 static const G4double mel=0.5109989;
134 static const G4double lmel=std::log(mel);
135
136 static std::vector <G4int> colN;
137 static std::vector <G4int> colZ;
138 static std::vector <G4int> colF;
139 static std::vector <G4double> colTH;
140 static std::vector <G4double> colH;
141
142
143 const G4double Energy = aPart->GetKineticEnergy()/MeV;
144 const G4int targetAtomicNumber = AA;
145 const G4int targZ = ZZ;
146 const G4int targN = targetAtomicNumber-targZ;
147 if (Energy<=EMi) return 0.;
148
149 G4int PDG=aPart->GetDefinition()->GetPDGEncoding();
150 if (PDG == 11 || PDG == -11)
151 {
153 if(targN!=lastN || targZ!=lastZ)
154 {
155 lastE = 0.;
156 lastG = 0.;
157 lastN = targN;
158 lastZ = targZ;
161 if(n)
for(
G4int i=0; i<
n; i++)
if(colN[i]==targN && colZ[i]==targZ)
162 {
163 in=true;
164 lastTH =colTH[i];
165 lastF =colF[i];
166 lastH =colH[i];
167 lastJ1 =J1[i];
168 lastJ2 =J2[i];
169 lastJ3 =J3[i];
170 }
171 if(!in)
172 {
176 lastF = GetFunctions(A,lastJ1,lastJ2,lastJ3);
177 lastH = alop*A*(1.-.072*std::log(A));
178 lastTH = ThresholdEnergy(targZ, targN);
179#ifdef pdebug
180 G4cout<<
"G4ElNucCS::GetCrossSection: lastH="<<lastH<<
",A="<<A<<
G4endl;
181#endif
182 colN.push_back(targN);
183 colZ.push_back(targZ);
184 colF.push_back(lastF);
185 J1.push_back(lastJ1);
186 J2.push_back(lastJ2);
187 J3.push_back(lastJ3);
188 colH.push_back(lastH);
189 colTH.push_back(lastTH);
190 }
191 }
192
193
194 else if(lastE == Energy) return lastSig*millibarn;
195
196 lastE=Energy;
197 if (Energy<=lastTH)
198 {
199 lastSig=0.;
200 return 0.;
201 }
203 lastG=lE-lmel;
206 if(lE<lEMa)
207 {
209 G4int blast=
static_cast<int>(shift);
210 if(blast<0) blast=0;
211 if(blast>=mL) blast=mL-1;
212 shift-=blast;
213 lastL=blast+1;
214 G4double YNi=dlg1*lastJ1[blast]-lgoe*(lastJ2[blast]+lastJ2[blast]-lastJ3[blast]/lastE);
215 G4double YNj=dlg1*lastJ1[lastL]-lgoe*(lastJ2[lastL]+lastJ2[lastL]-lastJ3[lastL]/lastE);
216 lastSig= YNi+shift*(YNj-YNi);
217 if(lastSig>YNj)lastSig=YNj;
218#ifdef pdebug
219 G4cout<<
"G4ElNucCS::GetCS:S="<<lastSig<<
",E="<<lE<<
",Yi="<<YNi<<
",Yj="<<YNj<<
",M="<<lEMa<<
G4endl;
220 G4cout<<
"G4EN::GCS:s="<<shift<<
",Jb="<<lastJ1[blast]<<
",J="<<lastJ1[lastL]<<
",b="<<blast<<
G4endl;
221#endif
222 }
223 else
224 {
225 lastL=mL;
226 G4double term1=lastJ1[mL]+lastH*HighEnergyJ1(lE);
227 G4double term2=lastJ2[mL]+lastH*HighEnergyJ2(lE);
228 G4double term3=lastJ3[mL]+lastH*HighEnergyJ3(lE);
229 lastSig=dlg1*term1-lgoe*(term2+term2-term3/lastE);
230#ifdef pdebug
231 G4cout<<
"G4ElNucCS::GetCrossSec:S="<<lastSig<<
",lE="<<lE<<
",J1="<<lastH*HighEnergyJ1(lE)<<
",Pm="
232 <<lastJ1[mL]<<
",Fm="<<lastJ2[mL]<<
",Fh="<<lastH*HighEnergyJ2(lE)<<
",EM="<<lEMa<<
G4endl;
233#endif
234 }
235 }
236 else return 0.;
237 if(lastSig<0.) lastSig = 0.;
238 lastE=Energy;
239 return lastSig*millibarn;
240}