101{
103 G4double Mproj = particle->GetPDGMass();
106
108
110
112
113 if (
Z == 1 &&
A == 1) theTargetDef = theProton;
114 else if (
Z == 1 &&
A == 2) theTargetDef = theDeuteron;
117 else if (
Z == 2 &&
A == 4) theTargetDef = theAlpha;
118
119
121
122
123
125 lv += Pproj;
127
129 Pproj.boost(-bst);
130
133
134 fbst = bst;
135 fptot= ptot;
136 fTmax = 4.0*ptot*ptot;
137
138 if(Plab < (std::abs(particle->GetBaryonNumber())*100)*MeV)
140
141
146 else { PrNucleonMass = NucleonMass; }
147 G4double energyPerN = std::sqrt(
sqr(PlabPerN) +
sqr(PrNucleonMass));
148 energyPerN -= PrNucleonMass;
149
150
151 G4double Z1 = particle->GetPDGCharge();
153
157 fWaveVector = ptot;
158
160 const G4double mevToBarn = 0.38938e+6;
161 G4double XsCoulomb = mevToBarn*
sqr(n/fWaveVector)*
pi*(1+ctet1)/(1.+Am)/(1.+2.*Am-ctet1);
162
165
166 XsElastHadronic/=millibarn; XsTotalHadronic/=millibarn;
167
168 G4double CoulombProb = XsCoulomb/(XsCoulomb+XsElastHadronic);
169
171 {
172
175
176 G4double par1 = 2.*(1.+Am)/(1.+ctet1);
177
178
179
180 G4double cosThetaCMS = (par1*ctet1- Ksi*(1.+2.*Am))/(par1-Ksi);
181
183 Fproj.setPz(PtZ);
184 G4double PtProjCMS = ptot*std::sqrt(1.0 - cosThetaCMS*cosThetaCMS);
185 G4double PtX= PtProjCMS * std::cos(phi);
186 G4double PtY= PtProjCMS * std::sin(phi);
187 Fproj.setPx(PtX);
188 Fproj.setPy(PtY);
189 Fproj.setE(std::sqrt(PtX*PtX+PtY*PtY+PtZ*PtZ+Mproj*Mproj));
190 T = -(Pproj-Fproj).mag2();
191 }
192 else
193 {
194
195
197
200
202
207
208 if((
A>=12.) && (
A<27) ) fRa=fRa*0.85;
209 if((
A>=27.) && (
A<48) ) fRa=fRa*0.90;
210 if((
A>=48.) && (
A<65) ) fRa=fRa*0.95;
211
215
216 if ((theParticle == theAProton) || (theParticle == theANeutron))
217 {
218 if(theTargetDef == theProton)
219 {
220
221 if(Plab < 610.)
222 { rho = 1.3347-10.342*Plab/1000.+22.277*Plab/1000.*Plab/1000.-
223 13.634*Plab/1000.*Plab/1000.*Plab/1000. ;}
224 if((Plab < 5500.)&&(Plab >= 610.) )
225 { rho = 0.22; }
226 if((Plab >= 5500.)&&(Plab < 12300.) )
227 { rho = -0.32; }
228 if( Plab >= 12300.)
229 { rho = 0.135-2.26/(std::sqrt(
S)) ;}
230 Ref2 = 0.35 + 0.9/std::sqrt(std::sqrt(
S-4.*0.88))+0.04*
G4Log(
S) ;
231 ceff2 = 0.375 - 2./
S + 0.44/(
sqr(
S-4.)+1.5) ;
232 Ref2 =Ref2*Ref2;
233 ceff2 = ceff2*ceff2;
234 }
235
237 {
238 Ref2 = fRa*fRa - 0.28 + 0.019 * sig_pbarp + 2.06e-6 * sig_pbarp*sig_pbarp;
239 ceff2 = 0.297 + 7.853e-04*sig_pbarp + 0.2899*
G4Exp(-0.03*sig_pbarp);
240 }
242 {
243 Ref2 = fRa*fRa - 1.36 + 0.025 * sig_pbarp - 3.69e-7 * sig_pbarp*sig_pbarp;
244 ceff2 = 0.149 + 7.091e-04*sig_pbarp + 0.3743*
G4Exp(-0.03*sig_pbarp);
245 }
247 {
248 Ref2 = fRa*fRa - 1.36 + 0.025 * sig_pbarp - 3.69e-7 * sig_pbarp*sig_pbarp;
249 ceff2 = 0.149 + 7.091e-04*sig_pbarp + 0.3743*
G4Exp(-0.03*sig_pbarp);
250 }
252 {
253 Ref2 = fRa*fRa -0.46 +0.03*sig_pbarp - 2.98e-6*sig_pbarp*sig_pbarp;
254 ceff2= 0.078 + 6.657e-4*sig_pbarp + 0.3359*
G4Exp(-0.03*sig_pbarp);
255 }
257 {
258 Ref2 = fRa*fRa +2.48*0.01*sig_pbarp*fRa - 2.23e-6*sig_pbarp*sig_pbarp*fRa*fRa;
259 ceff2 = 0.16+3.3e-4*sig_pbarp+0.35*
G4Exp(-0.03*sig_pbarp);
260 }
261 }
262
263 if (theParticle == theADeuteron)
264 {
265 if(theTargetDef == theProton)
266 {
267 ceff2 = 0.297 + 7.853e-04*sig_pbarp + 0.2899*
G4Exp(-0.03*sig_pbarp);
268 }
269 if(theTargetDef == theDeuteron)
270 {
271 ceff2 = 0.65 + 3.0e-4*sig_pbarp + 0.55 *
G4Exp(-0.03*sig_pbarp);
272 }
274 {
275 ceff2 = 0.57 + 2.5e-4*sig_pbarp + 0.65 *
G4Exp(-0.02*sig_pbarp);
276 }
277 if(theTargetDef == theAlpha)
278 {
279 ceff2 = 0.40 + 3.5e-4 *sig_pbarp + 0.45 *
G4Exp(-0.02*sig_pbarp);
280 }
282 {
283 ceff2 = 0.38 + 2.0e-4 *sig_pbarp + 0.5 *
G4Exp(-0.03*sig_pbarp);
284 }
285 }
286
287 if( (theParticle ==theAHe3) || (theParticle ==theATriton) )
288 {
289 if(theTargetDef == theProton)
290 {
291 ceff2 = 0.149 + 7.091e-04*sig_pbarp + 0.3743*
G4Exp(-0.03*sig_pbarp);
292 }
293 if(theTargetDef == theDeuteron)
294 {
295 ceff2 = 0.57 + 2.5e-4*sig_pbarp + 0.65 *
G4Exp(-0.02*sig_pbarp);
296 }
298 {
299 ceff2 = 0.39 + 2.7e-4*sig_pbarp + 0.7 *
G4Exp(-0.02*sig_pbarp);
300 }
301 if(theTargetDef == theAlpha)
302 {
303 ceff2 = 0.24 + 3.5e-4*sig_pbarp + 0.75 *
G4Exp(-0.03*sig_pbarp);
304 }
306 {
307 ceff2 = 0.26 + 2.2e-4*sig_pbarp + 0.33*
G4Exp(-0.03*sig_pbarp);
308 }
309 }
310
311 if ( (theParticle == theAAlpha) || (ceff2 == 0.0) )
312 {
313 if(theTargetDef == theProton)
314 {
315 ceff2= 0.078 + 6.657e-4*sig_pbarp + 0.3359*
G4Exp(-0.03*sig_pbarp);
316 }
317 if(theTargetDef == theDeuteron)
318 {
319 ceff2 = 0.40 + 3.5e-4 *sig_pbarp + 0.45 *
G4Exp(-0.02*sig_pbarp);
320 }
322 {
323 ceff2 = 0.24 + 3.5e-4*sig_pbarp + 0.75 *
G4Exp(-0.03*sig_pbarp);
324 }
325 if(theTargetDef == theAlpha)
326 {
327 ceff2 = 0.17 + 3.5e-4*sig_pbarp + 0.45 *
G4Exp(-0.03*sig_pbarp);
328 }
330 {
331 ceff2 = 0.22 + 2.0e-4*sig_pbarp + 0.2 *
G4Exp(-0.03*sig_pbarp);
332 }
333 }
334
335 fRef=std::sqrt(Ref2);
336 fceff = std::sqrt(ceff2);
337
340
341 const G4int maxNumberOfLoops = 10000;
342 G4int loopCounter = 0;
343 do
344 {
349 BracFunct = BracFunct * Q;
350 }
352 ++loopCounter < maxNumberOfLoops );
353 if ( loopCounter >= maxNumberOfLoops ) {
354 fTetaCMS = 0.0;
355 return 0.0;
356 }
357
359 T*=3.893913e+4;
360
361 }
362
363 return T;
364}
G4double S(G4double temp)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
G4double G4Log(G4double x)
G4double BesselOneByArg(G4double z)
G4double CalculateParticleBeta(const G4ParticleDefinition *particle, G4double momentum)
G4double CalculateZommerfeld(G4double beta, G4double Z1, G4double Z2)
G4double DampFactor(G4double z)
G4double GetcosTeta1(G4double plab, G4int A)
G4double BesselJzero(G4double z)
G4double CalculateAm(G4double momentum, G4double n, G4double Z)
virtual G4double GetTotalElementCrossSection(const G4ParticleDefinition *aParticle, G4double kinEnergy, G4int Z, G4double A)
virtual G4double GetElasticElementCrossSection(const G4ParticleDefinition *aParticle, G4double kinEnergy, G4int Z, G4double A)
G4double GetAntiHadronNucleonTotCrSc(const G4ParticleDefinition *aParticle, G4double kinEnergy)
static G4double GetNuclearMass(const G4double A, const G4double Z)
G4int GetBaryonNumber() const
G4double Z13(G4int Z) const
static G4Triton * Triton()
G4double energy(const ThreeVector &p, const G4double m)