| 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 // Geant4 Header : G4AntiNuclElastic
28 //
29 //
30
31 #include "G4AntiNuclElastic.hh"
32
33 #include "G4PhysicalConstants.hh"
34 #include "G4SystemOfUnits.hh"
35 #include "G4ParticleTable.hh"
36 #include "G4ParticleDefinition.hh"
37 #include "G4IonTable.hh"
38 #include "Randomize.hh"
39 #include "G4AntiProton.hh"
40 #include "G4AntiNeutron.hh"
41 #include "G4AntiDeuteron.hh"
42 #include "G4AntiAlpha.hh"
43 #include "G4AntiTriton.hh"
44 #include "G4AntiHe3.hh"
45 #include "G4Proton.hh"
46 #include "G4Neutron.hh"
47 #include "G4Deuteron.hh"
48 #include "G4Alpha.hh"
49 #include "G4Pow.hh"
50 #include "G4Exp.hh"
51 #include "G4Log.hh"
52
53 #include "G4NucleiProperties.hh"
54 #include "G4CrossSectionDataSetRegistry.hh"
55
56 G4AntiNuclElastic::G4AntiNuclElastic()
57 : G4HadronElastic("AntiAElastic")
58 {
59 //V.Ivanchenko commented out
60 //SetMinEnergy( 0.1*GeV );
61 //SetMaxEnergy( 10.*TeV );
62
63 theAProton = G4AntiProton::AntiProton(
64 theANeutron = G4AntiNeutron::AntiNeutro
65 theADeuteron = G4AntiDeuteron::AntiDeute
66 theATriton = G4AntiTriton::AntiTriton(
67 theAAlpha = G4AntiAlpha::AntiAlpha();
68 theAHe3 = G4AntiHe3::AntiHe3();
69
70 theProton = G4Proton::Proton();
71 theNeutron = G4Neutron::Neutron();
72 theDeuteron = G4Deuteron::Deuteron();
73 theAlpha = G4Alpha::Alpha();
74
75 G4CrossSectionDataSetRegistry* reg = G4Cross
76 cs = static_cast<G4ComponentAntiNuclNuclearX
77 if(!cs) { cs = new G4ComponentAntiNuclNuclea
78
79 fParticle = 0;
80 fWaveVector = 0.;
81 fBeta = 0.;
82 fZommerfeld = 0.;
83 fAm = 0.;
84 fTetaCMS = 0.;
85 fRa = 0.;
86 fRef = 0.;
87 fceff = 0.;
88 fptot = 0.;
89 fTmax = 0.;
90 fThetaLab = 0.;
91 }
92
93 //////////////////////////////////////////////
94 G4AntiNuclElastic::~G4AntiNuclElastic()
95 {}
96
97 //////////////////////////////////////////////
98 // sample momentum transfer in the CMS system
99 G4double G4AntiNuclElastic::SampleInvariantT(c
100 G4double Plab, G4int Z, G4int
101 {
102 G4double T;
103 G4double Mproj = particle->GetPDGMass();
104 G4LorentzVector Pproj(0.,0.,Plab,std::sqrt(P
105 G4double ctet1 = GetcosTeta1(Plab, A);
106
107 G4double energy=Pproj.e()-Mproj;
108
109 const G4ParticleDefinition* theParticle = pa
110
111 G4ParticleDefinition * theTargetDef = 0;
112
113 if (Z == 1 && A == 1) theTargetDef = th
114 else if (Z == 1 && A == 2) theTargetDef = th
115 else if (Z == 1 && A == 3) theTargetDef = G4
116 else if (Z == 2 && A == 3) theTargetDef = G4
117 else if (Z == 2 && A == 4) theTargetDef = th
118
119
120 G4double TargMass =G4NucleiProperties::GetNu
121
122 //transform to CMS
123
124 G4LorentzVector lv(0.0,0.0,0.0,TargMass);
125 lv += Pproj;
126 G4double S = lv.mag2()/(GeV*GeV);
127
128 G4ThreeVector bst = lv.boostVector();
129 Pproj.boost(-bst);
130
131 G4ThreeVector p1 = Pproj.vect();
132 G4double ptot = p1.mag();
133
134 fbst = bst;
135 fptot= ptot;
136 fTmax = 4.0*ptot*ptot; // In (MeV/c)^2
137
138 if(Plab < (std::abs(particle->GetBaryonNumbe
139 {return fTmax*G4UniformRand();}
140
141 // Calculation of NN collision properties
142 G4double PlabPerN = Plab/std::abs(theParticl
143 G4double NucleonMass = 0.5*( theProton->GetP
144 G4double PrNucleonMass(0.); // Projectile a
145 if( std::abs(theParticle->GetBaryonNumber())
146 else
147 G4double energyPerN = std::sqrt( sqr(PlabPer
148 energyPerN -= PrNucleonMass;
149 //---
150
151 G4double Z1 = particle->GetPDGCharge();
152 G4double Z2 = Z;
153
154 G4double beta = CalculateParticleBeta(partic
155 G4double n = CalculateZommerfeld( beta, Z1
156 G4double Am = CalculateAm( ptot, n, Z2 );
157 fWaveVector = ptot; // /hbarc;
158
159 G4LorentzVector Fproj(0.,0.,0.,0.);
160 const G4double mevToBarn = 0.38938e+6;
161 G4double XsCoulomb = mevToBarn*sqr(n/fWaveVe
162
163 G4double XsElastHadronic =cs->GetElasticElem
164 G4double XsTotalHadronic =cs->GetTotalElemen
165
166 XsElastHadronic/=millibarn; XsTotalHadronic/
167
168 G4double CoulombProb = XsCoulomb/(XsCoulomb
169
170 if(G4UniformRand() < CoulombProb)
171 { // Simulation of Coulomb scattering
172
173 G4double phi = twopi * G4UniformRand();
174 G4double Ksi = G4UniformRand();
175
176 G4double par1 = 2.*(1.+Am)/(1.+ctet1);
177
178 // ////sample ThetaCMS in Coulomb part
179
180 G4double cosThetaCMS = (par1*ctet1- Ksi*(1
181
182 G4double PtZ=ptot*cosThetaCMS;
183 Fproj.setPz(PtZ);
184 G4double PtProjCMS = ptot*std::sqrt(1.0 -
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*P
190 T = -(Pproj-Fproj).mag2();
191 }
192 else
193 {
194 // Simulation of strong interaction scatte
195
196 G4double Qmax = 2.*ptot/197.33; // in fm
197
198 G4double Amag = 1.0; // A1 in Majora
199 G4double SlopeMag = 0.5; // A2 in Majora
200
201 G4double sig_pbarp = cs->GetAntiHadronNucl
202
203 fRa = 1.113*G4Pow::GetInstance()->Z13(A) -
204 0.227/G4Pow::GetInstance()->Z13(A);
205 if(A == 3) fRa=1.81;
206 if(A == 4) fRa=1.37;
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
212 G4double Ref2 = XsTotalHadronic/10./2./pi;
213 G4double ceff2 = 0.0;
214 G4double rho = 0.0;
215
216 if ((theParticle == theAProton) || (thePa
217 {
218 if(theTargetDef == theProton)
219 {
220 // Determination of the real part of P
221 if(Plab < 610.)
222 { rho = 1.3347-10.342*Plab/1000.+22.27
223 13.634*Plab/1000.*Plab/1000.*P
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
231 ceff2 = 0.375 - 2./S + 0.44/(sqr(S-4.)
232 Ref2 =Ref2*Ref2;
233 ceff2 = ceff2*ceff2;
234 }
235
236 if( (Z==1)&&(A==2) )
237 {
238 Ref2 = fRa*fRa - 0.28 + 0.019 * sig_pb
239 ceff2 = 0.297 + 7.853e-04*sig_pbarp +
240 }
241 if( (Z==1)&&(A==3) )
242 {
243 Ref2 = fRa*fRa - 1.36 + 0.025 * sig_pb
244 ceff2 = 0.149 + 7.091e-04*sig_pbarp +
245 }
246 if( (Z==2)&&(A==3) )
247 {
248 Ref2 = fRa*fRa - 1.36 + 0.025 * sig_pb
249 ceff2 = 0.149 + 7.091e-04*sig_pbarp +
250 }
251 if( (Z==2)&&(A==4) )
252 {
253 Ref2 = fRa*fRa -0.46 +0.03*sig_pbarp -
254 ceff2= 0.078 + 6.657e-4*sig_pbarp + 0.
255 }
256 if(Z>2)
257 {
258 Ref2 = fRa*fRa +2.48*0.01*sig_pbarp*fR
259 ceff2 = 0.16+3.3e-4*sig_pbarp+0.35*G4E
260 }
261 } // End of if ((theParticle == theAProto
262
263 if (theParticle == theADeuteron)
264 {
265 if(theTargetDef == theProton)
266 {
267 ceff2 = 0.297 + 7.853e-04*sig_pbarp +
268 }
269 if(theTargetDef == theDeuteron)
270 {
271 ceff2 = 0.65 + 3.0e-4*sig_pbarp + 0.55
272 }
273 if( (theTargetDef == G4Triton::Triton())
274 {
275 ceff2 = 0.57 + 2.5e-4*sig_pbarp + 0.65
276 }
277 if(theTargetDef == theAlpha)
278 {
279 ceff2 = 0.40 + 3.5e-4 *sig_pbarp + 0.4
280 }
281 if(Z>2)
282 {
283 ceff2 = 0.38 + 2.0e-4 *sig_pbarp + 0.5
284 }
285 }
286
287 if( (theParticle ==theAHe3) || (theParticl
288 {
289 if(theTargetDef == theProton)
290 {
291 ceff2 = 0.149 + 7.091e-04*sig_pbarp +
292 }
293 if(theTargetDef == theDeuteron)
294 {
295 ceff2 = 0.57 + 2.5e-4*sig_pbarp + 0.65
296 }
297 if( (theTargetDef == G4Triton::Triton())
298 {
299 ceff2 = 0.39 + 2.7e-4*sig_pbarp + 0.7
300 }
301 if(theTargetDef == theAlpha)
302 {
303 ceff2 = 0.24 + 3.5e-4*sig_pbarp + 0.75
304 }
305 if(Z>2)
306 {
307 ceff2 = 0.26 + 2.2e-4*sig_pbarp + 0.33
308 }
309 }
310
311 if ( (theParticle == theAAlpha) || (ceff2
312 {
313 if(theTargetDef == theProton)
314 {
315 ceff2= 0.078 + 6.657e-4*sig_pbarp + 0.
316 }
317 if(theTargetDef == theDeuteron)
318 {
319 ceff2 = 0.40 + 3.5e-4 *sig_pbarp + 0.4
320 }
321 if( (theTargetDef == G4Triton::Triton())
322 {
323 ceff2 = 0.24 + 3.5e-4*sig_pbarp + 0.75
324 }
325 if(theTargetDef == theAlpha)
326 {
327 ceff2 = 0.17 + 3.5e-4*sig_pbarp + 0.45
328 }
329 if(Z>2)
330 {
331 ceff2 = 0.22 + 2.0e-4*sig_pbarp + 0.2
332 }
333 }
334
335 fRef=std::sqrt(Ref2);
336 fceff = std::sqrt(ceff2);
337
338 G4double Q = 0.0 ;
339 G4double BracFunct;
340
341 const G4int maxNumberOfLoops = 10000;
342 G4int loopCounter = 0;
343 do
344 {
345 Q = -G4Log(1.-(1.- G4Exp(-SlopeMag * Qma
346 G4double x = fRef * Q;
347 BracFunct = ( ( sqr(BesselOneByArg(x))+s
348 * sqr(DampFactor(pi*fceff*Q))) /(Ama
349 BracFunct = BracFunct * Q;
350 }
351 while ( (G4UniformRand()>BracFunct) &&
352 ++loopCounter < maxNumberOfLoops )
353 if ( loopCounter >= maxNumberOfLoops ) {
354 fTetaCMS = 0.0;
355 return 0.0;
356 }
357
358 T= sqr(Q);
359 T*=3.893913e+4; // fm^(-2) -> MeV^2
360
361 } // End of simulation of strong interactio
362
363 return T;
364 }
365
366 //////////////////////////////////////////////
367 // Sample of Theta in CMS
368 G4double G4AntiNuclElastic::SampleThetaCMS(co
369
370 {
371 G4double T;
372 T = SampleInvariantT( p, plab, Z, A);
373
374 // NaN finder
375 if(!(T < 0.0 || T >= 0.0))
376 {
377 if (verboseLevel > 0)
378 {
379 G4cout << "G4DiffuseElastic:WARNING: A =
380 << " mom(GeV)= " << plab/GeV
381 << " S-wave will be sampled"
382 << G4endl;
383 }
384 T = G4UniformRand()*fTmax;
385
386 }
387
388 if(fptot > 0.)
389 {
390 G4double cosTet=1.0-T/(2.*fptot*fptot);
391 if(cosTet > 1.0 ) cosTet= 1.;
392 if(cosTet < -1.0 ) cosTet=-1.;
393 fTetaCMS=std::acos(cosTet);
394 return fTetaCMS;
395 } else
396 {
397 return 2.*G4UniformRand()-1.;
398 }
399 }
400
401
402 //////////////////////////////////////////////
403 // Sample of Theta in Lab System
404 G4double G4AntiNuclElastic::SampleThetaLab(co
405
406 {
407 G4double T;
408 T = SampleInvariantT( p, plab, Z, A);
409
410 // NaN finder
411 if(!(T < 0.0 || T >= 0.0))
412 {
413 if (verboseLevel > 0)
414 {
415 G4cout << "G4DiffuseElastic:WARNING: A =
416 << " mom(GeV)= " << plab/GeV
417 << " S-wave will be sampled"
418 << G4endl;
419 }
420 T = G4UniformRand()*fTmax;
421 }
422
423 G4double phi = G4UniformRand()*twopi;
424
425 G4double cost(1.);
426 if(fTmax > 0.) {cost = 1. - 2.0*T/fTmax;}
427
428 G4double sint;
429 if( cost >= 1.0 )
430 {
431 cost = 1.0;
432 sint = 0.0;
433 }
434 else if( cost <= -1.0)
435 {
436 cost = -1.0;
437 sint = 0.0;
438 }
439 else
440 {
441 sint = std::sqrt((1.0-cost)*(1.0+cost));
442 }
443
444 G4double m1 = p->GetPDGMass();
445 G4ThreeVector v(sint*std::cos(phi),sint*std:
446 v *= fptot;
447 G4LorentzVector nlv(v.x(),v.y(),v.z(),std::s
448
449 nlv.boost(fbst);
450
451 G4ThreeVector np = nlv.vect();
452 G4double theta = np.theta();
453 fThetaLab = theta;
454
455 return theta;
456 }
457
458 //////////////////////////////////////////////
459 // Calculation of Damp factor
460 G4double G4AntiNuclElastic::DampFactor(G4doub
461 {
462 G4double df;
463 G4double f3 = 6.; // first factorials
464
465 if( std::fabs(x) < 0.01 )
466 {
467 df=1./(1.+x*x/f3);
468 }
469 else
470 {
471 df = x/std::sinh(x);
472 }
473 return df;
474 }
475
476
477 //////////////////////////////////////////////
478 // Calculation of particle velocity Beta
479
480 G4double G4AntiNuclElastic::CalculateParticle
481 G4double mom
482 {
483 G4double mass = particle->GetPDGMass();
484 G4double a = momentum/mass;
485 fBeta = a/std::sqrt(1+a*a);
486
487 return fBeta;
488 }
489
490
491 //////////////////////////////////////////////
492 // Calculation of parameter Zommerfeld
493
494 G4double G4AntiNuclElastic::CalculateZommerfe
495 {
496 fZommerfeld = fine_structure_const*Z1*Z2/bet
497
498 return fZommerfeld;
499 }
500
501 //////////////////////////////////////////////
502 //
503 G4double G4AntiNuclElastic::CalculateAm( G4dou
504 {
505 G4double k = momentum/hbarc;
506 G4double ch = 1.13 + 3.76*n*n;
507 G4double zn = 1.77*k/G4Pow::GetInstance()->
508 G4double zn2 = zn*zn;
509 fAm = ch/zn2;
510
511 return fAm;
512 }
513
514 //////////////////////////////////////////////
515 //
516 // Bessel J0 function based on rational approx
517 // J.F. Hart, Computer Approximations, New Yor
518
519 G4double G4AntiNuclElastic::BesselJzero(G4doub
520 {
521 G4double modvalue, value2, fact1, fact2, arg
522
523 modvalue = std::fabs(value);
524
525 if ( value < 8.0 && value > -8.0 )
526 {
527 value2 = value*value;
528
529 fact1 = 57568490574.0 + value2*(-13362590
530 + value2*( 65161964
531 + value2*(-11214424
532 + value2*( 77392.33
533 + value2*(-184.9052
534
535 fact2 = 57568490411.0 + value2*( 10295329
536 + value2*( 9494680.
537 + value2*(59272.648
538 + value2*(267.85327
539 + value2*1.0
540
541 bessel = fact1/fact2;
542 }
543 else
544 {
545 arg = 8.0/modvalue;
546
547 value2 = arg*arg;
548
549 shift = modvalue-0.785398164;
550
551 fact1 = 1.0 + value2*(-0.1098628627e-2
552 + value2*(0.2734510407e-4
553 + value2*(-0.2073370639e-5
554 + value2*0.2093887211e-6 )
555 fact2 = -0.1562499995e-1 + value2*(0.143048
556 + value2*(-0.691
557 + value2*(0.7621
558 - value2*0.93494
559
560 bessel = std::sqrt(0.636619772/modvalue)*(
561 }
562 return bessel;
563 }
564
565
566 //////////////////////////////////////////////
567 // Bessel J1 function based on rational approx
568 // J.F. Hart, Computer Approximations, New Yor
569
570 G4double G4AntiNuclElastic::BesselJone(G4doub
571 {
572 G4double modvalue, value2, fact1, fact2, arg
573
574 modvalue = std::fabs(value);
575
576 if ( modvalue < 8.0 )
577 {
578 value2 = value*value;
579 fact1 = value*(72362614232.0 + value2*(-7
580 + value2*( 2
581 + value2*(-2
582 + value2*( 1
583 + value2*(-3
584
585 fact2 = 144725228442.0 + value2*(23005351
586 + value2*(18583304
587 + value2*(99447.43
588 + value2*(376.9991
589 + value2*1.0
590 bessel = fact1/fact2;
591 }
592 else
593 {
594 arg = 8.0/modvalue;
595 value2 = arg*arg;
596
597 shift = modvalue - 2.356194491;
598
599 fact1 = 1.0 + value2*( 0.183105e-2
600 + value2*(-0.3516396496e-4
601 + value2*(0.2457520174e-5
602 + value2*(-0.240337019e-6
603
604 fact2 = 0.04687499995 + value2*(-0.2002690
605 + value2*( 0.8449199
606 + value2*(-0.8822898
607 + value2*0.105787412
608
609 bessel = std::sqrt( 0.636619772/modvalue)*
610 if (value < 0.0) bessel = -bessel;
611 }
612 return bessel;
613 }
614
615 //////////////////////////////////////////////
616 // return J1(x)/x with special case for small
617 G4double G4AntiNuclElastic::BesselOneByArg(G4d
618 {
619 G4double x2, result;
620
621 if( std::fabs(x) < 0.01 )
622 {
623 x *= 0.5;
624 x2 = x*x;
625 result = (2.- x2 + x2*x2/6.)/4.;
626 }
627 else
628 {
629 result = BesselJone(x)/x;
630 }
631 return result;
632 }
633
634 //////////////////////////////////////////////
635 // return angle from which Coulomb scattering
636 G4double G4AntiNuclElastic::GetcosTeta1(G4doub
637 {
638
639 // G4double p0 =G4LossTableManager::Instance()
640 G4double p0 = 1.*hbarc/fermi;
641 //G4double cteta1 = 1.0 - p0*p0/2.0 * pow(A,2.
642 G4double cteta1 = 1.0 - p0*p0/2.0 * G4Pow::G
643 //////////////////
644 if(cteta1 < -1.) cteta1 = -1.0;
645 return cteta1;
646 }
647
648
649
650
651
652
653
654