| 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 //
28 // Physics model class G4NuclNuclDiffuseElasti
29 //
30 //
31 // G4 Model: optical diffuse elastic scatterin
32 //
33 // 24-May-07 V. Grichine
34 //
35
36 #include "G4NuclNuclDiffuseElastic.hh"
37 #include "G4ParticleTable.hh"
38 #include "G4ParticleDefinition.hh"
39 #include "G4IonTable.hh"
40 #include "G4NucleiProperties.hh"
41
42 #include "Randomize.hh"
43 #include "G4Integrator.hh"
44 #include "globals.hh"
45 #include "G4PhysicalConstants.hh"
46 #include "G4SystemOfUnits.hh"
47
48 #include "G4Proton.hh"
49 #include "G4Neutron.hh"
50 #include "G4Deuteron.hh"
51 #include "G4Alpha.hh"
52 #include "G4PionPlus.hh"
53 #include "G4PionMinus.hh"
54
55 #include "G4Element.hh"
56 #include "G4ElementTable.hh"
57 #include "G4NistManager.hh"
58 #include "G4PhysicsTable.hh"
59 #include "G4PhysicsLogVector.hh"
60 #include "G4PhysicsFreeVector.hh"
61 #include "G4HadronicParameters.hh"
62
63 //////////////////////////////////////////////
64 //
65 // Test Constructor. Just to check xsc
66
67
68 G4NuclNuclDiffuseElastic::G4NuclNuclDiffuseEla
69 : G4HadronElastic("NNDiffuseElastic"), fPart
70 {
71 SetMinEnergy( 50*MeV );
72 SetMaxEnergy( G4HadronicParameters::Instance
73 verboseLevel = 0;
74 lowEnergyRecoilLimit = 100.*keV;
75 lowEnergyLimitQ = 0.0*GeV;
76 lowEnergyLimitHE = 0.0*GeV;
77 lowestEnergyLimit= 0.0*keV;
78 plabLowLimit = 20.0*MeV;
79
80 theProton = G4Proton::Proton();
81 theNeutron = G4Neutron::Neutron();
82 theDeuteron = G4Deuteron::Deuteron();
83 theAlpha = G4Alpha::Alpha();
84 thePionPlus = G4PionPlus::PionPlus();
85 thePionMinus= G4PionMinus::PionMinus();
86
87 fEnergyBin = 300; // Increased from the ori
88 fAngleBin = 200;
89
90 fEnergyVector = new G4PhysicsLogVector( the
91 fAngleTable = 0;
92
93 fParticle = 0;
94 fWaveVector = 0.;
95 fAtomicWeight = 0.;
96 fAtomicNumber = 0.;
97 fNuclearRadius = 0.;
98 fBeta = 0.;
99 fZommerfeld = 0.;
100 fAm = 0.;
101 fAddCoulomb = false;
102 // Ranges of angle table relative to current
103
104 // Empirical parameters
105
106 fCofAlphaMax = 1.5;
107 fCofAlphaCoulomb = 0.5;
108
109 fProfileDelta = 1.;
110 fProfileAlpha = 0.5;
111
112 fCofLambda = 1.0;
113 fCofDelta = 0.04;
114 fCofAlpha = 0.095;
115
116 fNuclearRadius1 = fNuclearRadius2 = fNuclear
117 = fRutherfordRatio = fCoulombPhase0 = fHal
118 = fRutherfordTheta = fProfileLambda = fCof
119 = fSumSigma = fEtaRatio = fReZ = 0.0;
120 fMaxL = 0;
121
122 fNuclearRadiusCof = 1.0;
123 fCoulombMuC = 0.0;
124 }
125
126 //////////////////////////////////////////////
127 //
128 // Destructor
129
130 G4NuclNuclDiffuseElastic::~G4NuclNuclDiffuseEl
131 {
132 if ( fEnergyVector ) {
133 delete fEnergyVector;
134 fEnergyVector = 0;
135 }
136
137 for ( std::vector<G4PhysicsTable*>::iterator
138 it != fAngleBank.end(); ++it ) {
139 if ( (*it) ) (*it)->clearAndDestroy();
140 delete *it;
141 *it = 0;
142 }
143 fAngleTable = 0;
144 }
145
146 //////////////////////////////////////////////
147 //
148 // Initialisation for given particle using ele
149
150 void G4NuclNuclDiffuseElastic::Initialise()
151 {
152
153 // fEnergyVector = new G4PhysicsLogVector( t
154
155 const G4ElementTable* theElementTable = G4El
156 std::size_t jEl, numOfEl = G4Element::GetNum
157
158 // projectile radius
159
160 G4double A1 = G4double( fParticle->GetBaryon
161 G4double R1 = CalculateNuclearRad(A1);
162
163 for(jEl = 0 ; jEl < numOfEl; ++jEl) // appli
164 {
165 fAtomicNumber = (*theElementTable)[jEl]->G
166 fAtomicWeight = G4NistManager::Instance()-
167
168 fNuclearRadius = CalculateNuclearRad(fAtom
169 fNuclearRadius += R1;
170
171 if(verboseLevel > 0)
172 {
173 G4cout<<"G4NuclNuclDiffuseElastic::Initi
174 <<(*theElementTable)[jEl]->GetName()<<G4
175 }
176 fElementNumberVector.push_back(fAtomicNumb
177 fElementNameVector.push_back((*theElementT
178
179 BuildAngleTable();
180 fAngleBank.push_back(fAngleTable);
181 }
182 }
183
184
185 //////////////////////////////////////////////
186 //
187 // return differential elastic cross section d
188
189 G4double
190 G4NuclNuclDiffuseElastic::GetDiffuseElasticXsc
191 G4doub
192 G4double momentum,
193 G4doub
194 {
195 fParticle = particle;
196 fWaveVector = momentum/hbarc;
197 fAtomicWeight = A;
198 fAddCoulomb = false;
199 fNuclearRadius = CalculateNuclearRad(A);
200
201 G4double sigma = fNuclearRadius*fNuclearRadi
202
203 return sigma;
204 }
205
206 //////////////////////////////////////////////
207 //
208 // return invariant differential elastic cross
209
210 G4double
211 G4NuclNuclDiffuseElastic::GetInvElasticXsc( co
212 G4doub
213 G4double plab,
214 G4doub
215 {
216 G4double m1 = particle->GetPDGMass();
217 G4LorentzVector lv1(0.,0.,plab,std::sqrt(pla
218
219 G4int iZ = static_cast<G4int>(Z+0.5);
220 G4int iA = static_cast<G4int>(A+0.5);
221 G4ParticleDefinition * theDef = 0;
222
223 if (iZ == 1 && iA == 1) theDef = thePro
224 else if (iZ == 1 && iA == 2) theDef = theDeu
225 else if (iZ == 1 && iA == 3) theDef = G4Trit
226 else if (iZ == 2 && iA == 3) theDef = G4He3:
227 else if (iZ == 2 && iA == 4) theDef = theAlp
228 else theDef = G4ParticleTable::GetParticleTa
229
230 G4double tmass = theDef->GetPDGMass();
231
232 G4LorentzVector lv(0.0,0.0,0.0,tmass);
233 lv += lv1;
234
235 G4ThreeVector bst = lv.boostVector();
236 lv1.boost(-bst);
237
238 G4ThreeVector p1 = lv1.vect();
239 G4double ptot = p1.mag();
240 G4double ptot2 = ptot*ptot;
241 G4double cost = 1 - 0.5*std::fabs(tMand)/pto
242
243 if( cost >= 1.0 ) cost = 1.0;
244 else if( cost <= -1.0) cost = -1.0;
245
246 G4double thetaCMS = std::acos(cost);
247
248 G4double sigma = GetDiffuseElasticXsc( parti
249
250 sigma *= pi/ptot2;
251
252 return sigma;
253 }
254
255 //////////////////////////////////////////////
256 //
257 // return differential elastic cross section d
258 // correction
259
260 G4double
261 G4NuclNuclDiffuseElastic::GetDiffuseElasticSum
262 G4doub
263 G4double momentum,
264 G4doub
265 {
266 fParticle = particle;
267 fWaveVector = momentum/hbarc;
268 fAtomicWeight = A;
269 fAtomicNumber = Z;
270 fNuclearRadius = CalculateNuclearRad(A);
271 fAddCoulomb = false;
272
273 G4double z = particle->GetPDGCharge();
274
275 G4double kRt = fWaveVector*fNuclearRadius*
276 G4double kRtC = 1.9;
277
278 if( z && (kRt > kRtC) )
279 {
280 fAddCoulomb = true;
281 fBeta = CalculateParticleBeta( parti
282 fZommerfeld = CalculateZommerfeld( fBeta,
283 fAm = CalculateAm( momentum, fZomm
284 }
285 G4double sigma = fNuclearRadius*fNuclearRadi
286
287 return sigma;
288 }
289
290 //////////////////////////////////////////////
291 //
292 // return invariant differential elastic cross
293 // correction
294
295 G4double
296 G4NuclNuclDiffuseElastic::GetInvElasticSumXsc(
297 G4doub
298 G4double plab,
299 G4doub
300 {
301 G4double m1 = particle->GetPDGMass();
302
303 G4LorentzVector lv1(0.,0.,plab,std::sqrt(pla
304
305 G4int iZ = static_cast<G4int>(Z+0.5);
306 G4int iA = static_cast<G4int>(A+0.5);
307
308 G4ParticleDefinition* theDef = 0;
309
310 if (iZ == 1 && iA == 1) theDef = thePro
311 else if (iZ == 1 && iA == 2) theDef = theDeu
312 else if (iZ == 1 && iA == 3) theDef = G4Trit
313 else if (iZ == 2 && iA == 3) theDef = G4He3:
314 else if (iZ == 2 && iA == 4) theDef = theAlp
315 else theDef = G4ParticleTable::GetParticleTa
316
317 G4double tmass = theDef->GetPDGMass();
318
319 G4LorentzVector lv(0.0,0.0,0.0,tmass);
320 lv += lv1;
321
322 G4ThreeVector bst = lv.boostVector();
323 lv1.boost(-bst);
324
325 G4ThreeVector p1 = lv1.vect();
326 G4double ptot = p1.mag();
327 G4double ptot2 = ptot*ptot;
328 G4double cost = 1 - 0.5*std::fabs(tMand)/
329
330 if( cost >= 1.0 ) cost = 1.0;
331 else if( cost <= -1.0) cost = -1.0;
332
333 G4double thetaCMS = std::acos(cost);
334
335 G4double sigma = GetDiffuseElasticSumXsc( pa
336
337 sigma *= pi/ptot2;
338
339 return sigma;
340 }
341
342 //////////////////////////////////////////////
343 //
344 // return invariant differential elastic cross
345 // correction
346
347 G4double
348 G4NuclNuclDiffuseElastic::GetInvCoulombElastic
349 G4doub
350 G4double plab,
351 G4doub
352 {
353 G4double m1 = particle->GetPDGMass();
354 G4LorentzVector lv1(0.,0.,plab,std::sqrt(pla
355
356 G4int iZ = static_cast<G4int>(Z+0.5);
357 G4int iA = static_cast<G4int>(A+0.5);
358 G4ParticleDefinition * theDef = 0;
359
360 if (iZ == 1 && iA == 1) theDef = thePro
361 else if (iZ == 1 && iA == 2) theDef = theDeu
362 else if (iZ == 1 && iA == 3) theDef = G4Trit
363 else if (iZ == 2 && iA == 3) theDef = G4He3:
364 else if (iZ == 2 && iA == 4) theDef = theAlp
365 else theDef = G4ParticleTable::GetParticleTa
366
367 G4double tmass = theDef->GetPDGMass();
368
369 G4LorentzVector lv(0.0,0.0,0.0,tmass);
370 lv += lv1;
371
372 G4ThreeVector bst = lv.boostVector();
373 lv1.boost(-bst);
374
375 G4ThreeVector p1 = lv1.vect();
376 G4double ptot = p1.mag();
377 G4double ptot2 = ptot*ptot;
378 G4double cost = 1 - 0.5*std::fabs(tMand)/pto
379
380 if( cost >= 1.0 ) cost = 1.0;
381 else if( cost <= -1.0) cost = -1.0;
382
383 G4double thetaCMS = std::acos(cost);
384
385 G4double sigma = GetCoulombElasticXsc( parti
386
387 sigma *= pi/ptot2;
388
389 return sigma;
390 }
391
392 //////////////////////////////////////////////
393 //
394 // return differential elastic probability d(p
395
396 G4double
397 G4NuclNuclDiffuseElastic::GetDiffElasticProb(
398 G4doub
399 // G4double momentum,
400 // G4double A
401 )
402 {
403 G4double sigma, bzero, bzero2, bonebyarg, bo
404 G4double delta, diffuse, gamma;
405 G4double e1, e2, bone, bone2;
406
407 // G4double wavek = momentum/hbarc; // wave
408 // G4double r0 = 1.08*fermi;
409 // G4double rad = r0*G4Pow::GetInstance()-
410
411 G4double kr = fWaveVector*fNuclearRadius;
412 G4double kr2 = kr*kr;
413 G4double krt = kr*theta;
414
415 bzero = BesselJzero(krt);
416 bzero2 = bzero*bzero;
417 bone = BesselJone(krt);
418 bone2 = bone*bone;
419 bonebyarg = BesselOneByArg(krt);
420 bonebyarg2 = bonebyarg*bonebyarg;
421
422 // VI - Coverity complains
423 /*
424 if (fParticle == theProton)
425 {
426 diffuse = 0.63*fermi;
427 gamma = 0.3*fermi;
428 delta = 0.1*fermi*fermi;
429 e1 = 0.3*fermi;
430 e2 = 0.35*fermi;
431 }
432 else // as proton, if were not defined
433 {
434 */
435 diffuse = 0.63*fermi;
436 gamma = 0.3*fermi;
437 delta = 0.1*fermi*fermi;
438 e1 = 0.3*fermi;
439 e2 = 0.35*fermi;
440 //}
441 G4double lambda = 15.; // 15 ok
442
443 // G4double kgamma = fWaveVector*gamma;
444
445 G4double kgamma = lambda*(1.-G4Exp(-fWave
446 G4double kgamma2 = kgamma*kgamma;
447
448 // G4double dk2t = delta*fWaveVector*fWaveV
449 // G4double dk2t2 = dk2t*dk2t;
450 // G4double pikdt = pi*fWaveVector*diffuse*t
451
452 G4double pikdt = lambda*(1.-G4Exp(-pi*fWa
453
454 damp = DampFactor(pikdt);
455 damp2 = damp*damp;
456
457 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector
458 G4double e2dk3t = -2.*e2*delta*fWaveVector*
459
460
461 sigma = kgamma2;
462 // sigma += dk2t2;
463 sigma *= bzero2;
464 sigma += mode2k2*bone2 + e2dk3t*bzero*bone;
465 sigma += kr2*bonebyarg2;
466 sigma *= damp2; // *rad*rad;
467
468 return sigma;
469 }
470
471 //////////////////////////////////////////////
472 //
473 // return differential elastic probability d(p
474 // Coulomb correction
475
476 G4double
477 G4NuclNuclDiffuseElastic::GetDiffElasticSumPro
478 G4doub
479 // G4double momentum,
480 // G4double A
481 )
482 {
483 G4double sigma, bzero, bzero2, bonebyarg, bo
484 G4double delta, diffuse, gamma;
485 G4double e1, e2, bone, bone2;
486
487 // G4double wavek = momentum/hbarc; // wave
488 // G4double r0 = 1.08*fermi;
489 // G4double rad = r0*G4Pow::GetInstance()-
490
491 G4double kr = fWaveVector*fNuclearRadius;
492 G4double kr2 = kr*kr;
493 G4double krt = kr*theta;
494
495 bzero = BesselJzero(krt);
496 bzero2 = bzero*bzero;
497 bone = BesselJone(krt);
498 bone2 = bone*bone;
499 bonebyarg = BesselOneByArg(krt);
500 bonebyarg2 = bonebyarg*bonebyarg;
501
502 if (fParticle == theProton)
503 {
504 diffuse = 0.63*fermi;
505 // diffuse = 0.6*fermi;
506 gamma = 0.3*fermi;
507 delta = 0.1*fermi*fermi;
508 e1 = 0.3*fermi;
509 e2 = 0.35*fermi;
510 }
511 else // as proton, if were not defined
512 {
513 diffuse = 0.63*fermi;
514 gamma = 0.3*fermi;
515 delta = 0.1*fermi*fermi;
516 e1 = 0.3*fermi;
517 e2 = 0.35*fermi;
518 }
519 G4double lambda = 15.; // 15 ok
520 // G4double kgamma = fWaveVector*gamma;
521 G4double kgamma = lambda*(1.-G4Exp(-fWave
522
523 // G4cout<<"kgamma = "<<kgamma<<G4endl;
524
525 if(fAddCoulomb) // add Coulomb correction
526 {
527 G4double sinHalfTheta = std::sin(0.5*thet
528 G4double sinHalfTheta2 = sinHalfTheta*sinH
529
530 kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta
531 // kgamma += 0.65*fZommerfeld/kr/(sinHalfThe
532 }
533
534 G4double kgamma2 = kgamma*kgamma;
535
536 // G4double dk2t = delta*fWaveVector*fWaveV
537 // G4cout<<"dk2t = "<<dk2t<<G4endl;
538 // G4double dk2t2 = dk2t*dk2t;
539 // G4double pikdt = pi*fWaveVector*diffuse*t
540
541 G4double pikdt = lambda*(1.-G4Exp(-pi*fWa
542
543 // G4cout<<"pikdt = "<<pikdt<<G4endl;
544
545 damp = DampFactor(pikdt);
546 damp2 = damp*damp;
547
548 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector
549 G4double e2dk3t = -2.*e2*delta*fWaveVector*
550
551 sigma = kgamma2;
552 // sigma += dk2t2;
553 sigma *= bzero2;
554 sigma += mode2k2*bone2;
555 sigma += e2dk3t*bzero*bone;
556
557 // sigma += kr2*(1 + 8.*fZommerfeld*fZommerf
558 sigma += kr2*bonebyarg2; // correction at J
559
560 sigma *= damp2; // *rad*rad;
561
562 return sigma;
563 }
564
565
566 //////////////////////////////////////////////
567 //
568 // return differential elastic probability d(p
569 // Coulomb correction
570
571 G4double
572 G4NuclNuclDiffuseElastic::GetDiffElasticSumPro
573 {
574 G4double theta;
575
576 theta = std::sqrt(alpha);
577
578 // theta = std::acos( 1 - alpha/2. );
579
580 G4double sigma, bzero, bzero2, bonebyarg, bo
581 G4double delta, diffuse, gamma;
582 G4double e1, e2, bone, bone2;
583
584 // G4double wavek = momentum/hbarc; // wave
585 // G4double r0 = 1.08*fermi;
586 // G4double rad = r0*G4Pow::GetInstance()-
587
588 G4double kr = fWaveVector*fNuclearRadius;
589 G4double kr2 = kr*kr;
590 G4double krt = kr*theta;
591
592 bzero = BesselJzero(krt);
593 bzero2 = bzero*bzero;
594 bone = BesselJone(krt);
595 bone2 = bone*bone;
596 bonebyarg = BesselOneByArg(krt);
597 bonebyarg2 = bonebyarg*bonebyarg;
598
599 if (fParticle == theProton)
600 {
601 diffuse = 0.63*fermi;
602 // diffuse = 0.6*fermi;
603 gamma = 0.3*fermi;
604 delta = 0.1*fermi*fermi;
605 e1 = 0.3*fermi;
606 e2 = 0.35*fermi;
607 }
608 else // as proton, if were not defined
609 {
610 diffuse = 0.63*fermi;
611 gamma = 0.3*fermi;
612 delta = 0.1*fermi*fermi;
613 e1 = 0.3*fermi;
614 e2 = 0.35*fermi;
615 }
616 G4double lambda = 15.; // 15 ok
617 // G4double kgamma = fWaveVector*gamma;
618 G4double kgamma = lambda*(1.-G4Exp(-fWave
619
620 // G4cout<<"kgamma = "<<kgamma<<G4endl;
621
622 if(fAddCoulomb) // add Coulomb correction
623 {
624 G4double sinHalfTheta = theta*0.5; // std
625 G4double sinHalfTheta2 = sinHalfTheta*sinH
626
627 kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta
628 // kgamma += 0.65*fZommerfeld/kr/(sinHalfThe
629 }
630
631 G4double kgamma2 = kgamma*kgamma;
632
633 // G4double dk2t = delta*fWaveVector*fWaveV
634 // G4cout<<"dk2t = "<<dk2t<<G4endl;
635 // G4double dk2t2 = dk2t*dk2t;
636 // G4double pikdt = pi*fWaveVector*diffuse*t
637
638 G4double pikdt = lambda*(1.-G4Exp(-pi*fWa
639
640 // G4cout<<"pikdt = "<<pikdt<<G4endl;
641
642 damp = DampFactor(pikdt);
643 damp2 = damp*damp;
644
645 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector
646 G4double e2dk3t = -2.*e2*delta*fWaveVector*
647
648 sigma = kgamma2;
649 // sigma += dk2t2;
650 sigma *= bzero2;
651 sigma += mode2k2*bone2;
652 sigma += e2dk3t*bzero*bone;
653
654 // sigma += kr2*(1 + 8.*fZommerfeld*fZommerf
655 sigma += kr2*bonebyarg2; // correction at J
656
657 sigma *= damp2; // *rad*rad;
658
659 return sigma;
660 }
661
662
663 //////////////////////////////////////////////
664 //
665 // return differential elastic probability 2*p
666
667 G4double
668 G4NuclNuclDiffuseElastic::GetIntegrandFunction
669 {
670 G4double result;
671
672 result = GetDiffElasticSumProbA(alpha);
673
674 // result *= 2*pi*std::sin(theta);
675
676 return result;
677 }
678
679 //////////////////////////////////////////////
680 //
681 // return integral elastic cross section d(sig
682
683 G4double
684 G4NuclNuclDiffuseElastic::IntegralElasticProb(
685 G4doub
686 G4double momentum,
687 G4doub
688 {
689 G4double result;
690 fParticle = particle;
691 fWaveVector = momentum/hbarc;
692 fAtomicWeight = A;
693
694 fNuclearRadius = CalculateNuclearRad(A);
695
696
697 G4Integrator<G4NuclNuclDiffuseElastic,G4doub
698
699 // result = integral.Legendre10(this,&G4Nucl
700 result = integral.Legendre96(this,&G4NuclNuc
701
702 return result;
703 }
704
705 //////////////////////////////////////////////
706 //
707 // Return inv momentum transfer -t > 0
708
709 G4double G4NuclNuclDiffuseElastic::SampleT( co
710 G4double p, G4double A)
711 {
712 G4double theta = SampleThetaCMS( aParticle,
713 G4double t = 2*p*p*( 1 - std::cos(theta)
714 return t;
715 }
716
717 //////////////////////////////////////////////
718 //
719 // Return scattering angle sampled in cms
720
721
722 G4double
723 G4NuclNuclDiffuseElastic::SampleThetaCMS(const
724 G4doubl
725 {
726 G4int i, iMax = 100;
727 G4double norm, theta1, theta2, thetaMax;
728 G4double result = 0., sum = 0.;
729
730 fParticle = particle;
731 fWaveVector = momentum/hbarc;
732 fAtomicWeight = A;
733
734 fNuclearRadius = CalculateNuclearRad(A);
735
736 thetaMax = 10.174/fWaveVector/fNuclearRadius
737
738 if (thetaMax > pi) thetaMax = pi;
739
740 G4Integrator<G4NuclNuclDiffuseElastic,G4doub
741
742 // result = integral.Legendre10(this,&G4Nucl
743 norm = integral.Legendre96(this,&G4NuclNuclD
744
745 norm *= G4UniformRand();
746
747 for(i = 1; i <= iMax; i++)
748 {
749 theta1 = (i-1)*thetaMax/iMax;
750 theta2 = i*thetaMax/iMax;
751 sum += integral.Legendre10(this,&G4NuclN
752
753 if ( sum >= norm )
754 {
755 result = 0.5*(theta1 + theta2);
756 break;
757 }
758 }
759 if (i > iMax ) result = 0.5*(theta1 + theta2
760
761 G4double sigma = pi*thetaMax/iMax;
762
763 result += G4RandGauss::shoot(0.,sigma);
764
765 if(result < 0.) result = 0.;
766 if(result > thetaMax) result = thetaMax;
767
768 return result;
769 }
770
771 //////////////////////////////////////////////
772 ///////////////////// Table preparation and r
773
774 //////////////////////////////////////////////
775 //
776 // Interface function. Return inv momentum tra
777
778 G4double G4NuclNuclDiffuseElastic::SampleInvar
779
780 {
781 fParticle = aParticle;
782 fAtomicNumber = G4double(Z);
783 fAtomicWeight = G4double(A);
784
785 G4double t(0.), m1 = fParticle->GetPDGMass()
786 G4double totElab = std::sqrt(m1*m1+p*p);
787 G4double mass2 = G4NucleiProperties::GetNucl
788 G4LorentzVector lv1(p,0.0,0.0,totElab);
789 G4LorentzVector lv(0.0,0.0,0.0,mass2);
790 lv += lv1;
791
792 G4ThreeVector bst = lv.boostVector();
793 lv1.boost(-bst);
794
795 G4ThreeVector p1 = lv1.vect();
796 G4double momentumCMS = p1.mag();
797
798 // t = SampleTableT( aParticle, momentumCMS
799
800 t = SampleCoulombMuCMS( aParticle, momentum
801
802
803 return t;
804 }
805
806 //////////////////////////////////////////////
807 //
808 // Return inv momentum transfer -t > 0 as Coul
809
810 G4double G4NuclNuclDiffuseElastic::SampleCoulo
811 {
812 G4double t(0.), rand(0.), mu(0.);
813
814 G4double A1 = G4double( aParticle->GetBaryon
815 G4double R1 = CalculateNuclearRad(A1);
816
817 fNuclearRadius = CalculateNuclearRad(fAtomi
818 fNuclearRadius += R1;
819
820 InitDynParameters(fParticle, p);
821
822 fCoulombMuC = fHalfRutThetaTg2/(1.+fHalfRutT
823
824 rand = G4UniformRand();
825
826 // sample (1-cosTheta) in 0 < Theta < ThetaC
827
828 mu = fCoulombMuC*rand*fAm;
829 mu /= fAm + fCoulombMuC*( 1. - rand );
830
831 t = 4.*p*p*mu;
832
833 return t;
834 }
835
836
837 //////////////////////////////////////////////
838 //
839 // Return inv momentum transfer -t > 0 from in
840
841 G4double G4NuclNuclDiffuseElastic::SampleTable
842
843 {
844 G4double alpha = SampleTableThetaCMS( aParti
845 // G4double t = 2*p*p*( 1 - std::cos(std
846 G4double t = p*p*alpha; // -
847 return t;
848 }
849
850 //////////////////////////////////////////////
851 //
852 // Return scattering angle2 sampled in cms acc
853
854
855 G4double
856 G4NuclNuclDiffuseElastic::SampleTableThetaCMS(
857 G4doubl
858 {
859 std::size_t iElement;
860 G4int iMomentum, iAngle;
861 G4double randAngle, position, theta1, theta2
862 G4double m1 = particle->GetPDGMass();
863
864 for(iElement = 0; iElement < fElementNumberV
865 {
866 if( std::fabs(Z - fElementNumberVector[iEl
867 }
868 if ( iElement == fElementNumberVector.size()
869 {
870 InitialiseOnFly(Z,A); // table preparation
871
872 // iElement--;
873
874 // G4cout << "G4NuclNuclDiffuseElastic: El
875 // << " is not found, return zero angle" <
876 // return 0.; // no table for this element
877 }
878 // G4cout<<"iElement = "<<iElement<<G4endl;
879
880 fAngleTable = fAngleBank[iElement];
881
882 G4double kinE = std::sqrt(momentum*momentum
883
884 for( iMomentum = 0; iMomentum < fEnergyBin;
885 {
886 // G4cout<<iMomentum<<", kinE = "<<kinE/Me
887
888 if( kinE < fEnergyVector->GetLowEdgeEnergy
889 }
890 // G4cout<<"iMomentum = "<<iMomentum<<", fEn
891
892
893 if ( iMomentum >= fEnergyBin ) iMomentum = f
894 if ( iMomentum < 0 ) iMomentum = 0
895
896
897 if (iMomentum == fEnergyBin -1 || iMomentum
898 {
899 position = (*(*fAngleTable)(iMomentum))(fA
900
901 // G4cout<<"position = "<<position<<G4endl
902
903 for(iAngle = 0; iAngle < fAngleBin-1; iAng
904 {
905 if( position < (*(*fAngleTable)(iMomentu
906 }
907 if (iAngle >= fAngleBin-1) iAngle = fAngle
908
909 // G4cout<<"iAngle = "<<iAngle<<G4endl;
910
911 randAngle = GetScatteringAngle(iMomentum,
912
913 // G4cout<<"randAngle = "<<randAngle<<G4en
914 }
915 else // kinE inside between energy table ed
916 {
917 // position = (*(*fAngleTable)(iMomentum))
918 position = (*(*fAngleTable)(iMomentum))(0)
919
920 // G4cout<<"position = "<<position<<G4endl
921
922 for(iAngle = 0; iAngle < fAngleBin-1; iAng
923 {
924 // if( position < (*(*fAngleTable)(iMome
925 if( position > (*(*fAngleTable)(iMomentu
926 }
927 if (iAngle >= fAngleBin-1) iAngle = fAngle
928
929 // G4cout<<"iAngle = "<<iAngle<<G4endl;
930
931 theta2 = GetScatteringAngle(iMomentum, iA
932
933 // G4cout<<"theta2 = "<<theta2<<G4endl;
934
935 E2 = fEnergyVector->GetLowEdgeEnergy(iMome
936
937 // G4cout<<"E2 = "<<E2<<G4endl;
938
939 iMomentum--;
940
941 // position = (*(*fAngleTable)(iMomentum))
942
943 // G4cout<<"position = "<<position<<G4endl
944
945 for(iAngle = 0; iAngle < fAngleBin-1; iAng
946 {
947 // if( position < (*(*fAngleTable)(iMome
948 if( position > (*(*fAngleTable)(iMomentu
949 }
950 if (iAngle >= fAngleBin-1) iAngle = fAngle
951
952 theta1 = GetScatteringAngle(iMomentum, iA
953
954 // G4cout<<"theta1 = "<<theta1<<G4endl;
955
956 E1 = fEnergyVector->GetLowEdgeEnergy(iMome
957
958 // G4cout<<"E1 = "<<E1<<G4endl;
959
960 W = 1.0/(E2 - E1);
961 W1 = (E2 - kinE)*W;
962 W2 = (kinE - E1)*W;
963
964 randAngle = W1*theta1 + W2*theta2;
965
966 // randAngle = theta2;
967 }
968 // G4double angle = randAngle;
969 // if (randAngle > 0.) randAngle /= 2*pi*std
970 // G4cout<<"randAngle = "<<randAngle/degree<
971
972 return randAngle;
973 }
974
975 //////////////////////////////////////////////
976 //
977 // Initialisation for given particle on fly us
978
979 void G4NuclNuclDiffuseElastic::InitialiseOnFly
980 {
981 fAtomicNumber = Z; // atomic number
982 fAtomicWeight = G4NistManager::Instance()->
983
984 G4double A1 = G4double( fParticle->GetBaryon
985 G4double R1 = CalculateNuclearRad(A1);
986
987 fNuclearRadius = CalculateNuclearRad(fAtomic
988
989 if( verboseLevel > 0 )
990 {
991 G4cout<<"G4NuclNuclDiffuseElastic::Initial
992 <<Z<<"; and A = "<<A<<G4endl;
993 }
994 fElementNumberVector.push_back(fAtomicNumber
995
996 BuildAngleTable();
997
998 fAngleBank.push_back(fAngleTable);
999
1000 return;
1001 }
1002
1003 /////////////////////////////////////////////
1004 //
1005 // Build for given particle and element table
1006 // For the moment in lab system.
1007
1008 void G4NuclNuclDiffuseElastic::BuildAngleTabl
1009 {
1010 G4int i, j;
1011 G4double partMom, kinE, m1 = fParticle->Get
1012 G4double alpha1, alpha2, alphaMax, alphaCou
1013
1014 // G4cout<<"particle z = "<<z<<"; particle
1015
1016 G4Integrator<G4NuclNuclDiffuseElastic,G4dou
1017
1018 fAngleTable = new G4PhysicsTable(fEnergyBin
1019
1020 for( i = 0; i < fEnergyBin; i++)
1021 {
1022 kinE = fEnergyVector->GetLowEdgeEn
1023
1024 // G4cout<<G4endl;
1025 // G4cout<<"kinE = "<<kinE/MeV<<" MeV"<<G
1026
1027 partMom = std::sqrt( kinE*(kinE + 2*m
1028
1029 InitDynParameters(fParticle, partMom);
1030
1031 alphaMax = fRutherfordTheta*fCofAlphaMax;
1032
1033 if(alphaMax > pi) alphaMax = pi;
1034
1035 // VI: Coverity complain
1036 //alphaMax = pi2;
1037
1038 alphaCoulomb = fRutherfordTheta*fCofAlpha
1039
1040 // G4cout<<"alphaCoulomb = "<<alphaCoulom
1041
1042 G4PhysicsFreeVector* angleVector = new G4
1043
1044 // G4PhysicsLogVector* angleBins = new G
1045
1046 G4double delth = (alphaMax-alphaCoulomb)/
1047
1048 sum = 0.;
1049
1050 // fAddCoulomb = false;
1051 fAddCoulomb = true;
1052
1053 // for(j = 1; j < fAngleBin; j++)
1054 for(j = fAngleBin-1; j >= 1; j--)
1055 {
1056 // alpha1 = angleBins->GetLowEdgeEnergy
1057 // alpha2 = angleBins->GetLowEdgeEnergy
1058
1059 // alpha1 = alphaMax*(j-1)/fAngleBin;
1060 // alpha2 = alphaMax*( j )/fAngleBin;
1061
1062 alpha1 = alphaCoulomb + delth*(j-1);
1063 // if(alpha1 < kRlim2) alpha1 = kRlim2;
1064 alpha2 = alpha1 + delth;
1065
1066 delta = integral.Legendre10(this, &G4Nu
1067 // delta = integral.Legendre96(this, &G
1068
1069 sum += delta;
1070
1071 angleVector->PutValue( j-1 , alpha1, su
1072 // G4cout<<"j-1 = "<<j-1<<"; alpha2 = "
1073 }
1074 fAngleTable->insertAt(i,angleVector);
1075
1076 // delete[] angleVector;
1077 // delete[] angleBins;
1078 }
1079 return;
1080 }
1081
1082 /////////////////////////////////////////////
1083 //
1084 //
1085
1086 G4double
1087 G4NuclNuclDiffuseElastic:: GetScatteringAngle
1088 {
1089 G4double x1, x2, y1, y2, randAngle;
1090
1091 if( iAngle == 0 )
1092 {
1093 randAngle = (*fAngleTable)(iMomentum)->Ge
1094 // iAngle++;
1095 }
1096 else
1097 {
1098 if ( iAngle >= G4int((*fAngleTable)(iMome
1099 {
1100 iAngle = G4int((*fAngleTable)(iMomentum
1101 }
1102 y1 = (*(*fAngleTable)(iMomentum))(iAngle-
1103 y2 = (*(*fAngleTable)(iMomentum))(iAngle)
1104
1105 x1 = (*fAngleTable)(iMomentum)->GetLowEdg
1106 x2 = (*fAngleTable)(iMomentum)->GetLowEdg
1107
1108 if ( x1 == x2 ) randAngle = x2;
1109 else
1110 {
1111 if ( y1 == y2 ) randAngle = x1 + ( x2 -
1112 else
1113 {
1114 randAngle = x1 + ( position - y1 )*(
1115 }
1116 }
1117 }
1118 return randAngle;
1119 }
1120
1121
1122
1123 /////////////////////////////////////////////
1124 //
1125 // Return scattering angle sampled in lab sys
1126
1127
1128
1129 G4double
1130 G4NuclNuclDiffuseElastic::SampleThetaLab( con
1131 G4dou
1132 {
1133 const G4ParticleDefinition* theParticle = a
1134 G4double m1 = theParticle->GetPDGMass();
1135 G4double plab = aParticle->GetTotalMomentum
1136 G4LorentzVector lv1 = aParticle->Get4Moment
1137 G4LorentzVector lv(0.0,0.0,0.0,tmass);
1138 lv += lv1;
1139
1140 G4ThreeVector bst = lv.boostVector();
1141 lv1.boost(-bst);
1142
1143 G4ThreeVector p1 = lv1.vect();
1144 G4double ptot = p1.mag();
1145 G4double tmax = 4.0*ptot*ptot;
1146 G4double t = 0.0;
1147
1148
1149 //
1150 // Sample t
1151 //
1152
1153 t = SampleT( theParticle, ptot, A);
1154
1155 // NaN finder
1156 if(!(t < 0.0 || t >= 0.0))
1157 {
1158 if (verboseLevel > 0)
1159 {
1160 G4cout << "G4NuclNuclDiffuseElastic:WAR
1161 << " mom(GeV)= " << plab/GeV
1162 << " S-wave will be sampled"
1163 << G4endl;
1164 }
1165 t = G4UniformRand()*tmax;
1166 }
1167 if(verboseLevel>1)
1168 {
1169 G4cout <<" t= " << t << " tmax= " << tmax
1170 << " ptot= " << ptot << G4endl;
1171 }
1172 // Sampling of angles in CM system
1173
1174 G4double phi = G4UniformRand()*twopi;
1175 G4double cost = 1. - 2.0*t/tmax;
1176 G4double sint;
1177
1178 if( cost >= 1.0 )
1179 {
1180 cost = 1.0;
1181 sint = 0.0;
1182 }
1183 else if( cost <= -1.0)
1184 {
1185 cost = -1.0;
1186 sint = 0.0;
1187 }
1188 else
1189 {
1190 sint = std::sqrt((1.0-cost)*(1.0+cost));
1191 }
1192 if (verboseLevel>1)
1193 {
1194 G4cout << "cos(t)=" << cost << " std::sin
1195 }
1196 G4ThreeVector v1(sint*std::cos(phi),sint*st
1197 v1 *= ptot;
1198 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),s
1199
1200 nlv1.boost(bst);
1201
1202 G4ThreeVector np1 = nlv1.vect();
1203
1204 // G4double theta = std::acos( np1.z()/np
1205
1206 G4double theta = np1.theta();
1207
1208 return theta;
1209 }
1210
1211 /////////////////////////////////////////////
1212 //
1213 // Return scattering angle in lab system (tar
1214
1215
1216
1217 G4double
1218 G4NuclNuclDiffuseElastic::ThetaCMStoThetaLab(
1219 G4dou
1220 {
1221 const G4ParticleDefinition* theParticle = a
1222 G4double m1 = theParticle->GetPDGMass();
1223 // G4double plab = aParticle->GetTotalMomen
1224 G4LorentzVector lv1 = aParticle->Get4Moment
1225 G4LorentzVector lv(0.0,0.0,0.0,tmass);
1226
1227 lv += lv1;
1228
1229 G4ThreeVector bst = lv.boostVector();
1230
1231 lv1.boost(-bst);
1232
1233 G4ThreeVector p1 = lv1.vect();
1234 G4double ptot = p1.mag();
1235
1236 G4double phi = G4UniformRand()*twopi;
1237 G4double cost = std::cos(thetaCMS);
1238 G4double sint;
1239
1240 if( cost >= 1.0 )
1241 {
1242 cost = 1.0;
1243 sint = 0.0;
1244 }
1245 else if( cost <= -1.0)
1246 {
1247 cost = -1.0;
1248 sint = 0.0;
1249 }
1250 else
1251 {
1252 sint = std::sqrt((1.0-cost)*(1.0+cost));
1253 }
1254 if (verboseLevel>1)
1255 {
1256 G4cout << "cos(tcms)=" << cost << " std::
1257 }
1258 G4ThreeVector v1(sint*std::cos(phi),sint*st
1259 v1 *= ptot;
1260 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),s
1261
1262 nlv1.boost(bst);
1263
1264 G4ThreeVector np1 = nlv1.vect();
1265
1266
1267 G4double thetaLab = np1.theta();
1268
1269 return thetaLab;
1270 }
1271
1272 /////////////////////////////////////////////
1273 //
1274 // Return scattering angle in CMS system (tar
1275
1276
1277
1278 G4double
1279 G4NuclNuclDiffuseElastic::ThetaLabToThetaCMS(
1280 G4dou
1281 {
1282 const G4ParticleDefinition* theParticle = a
1283 G4double m1 = theParticle->GetPDGMass();
1284 G4double plab = aParticle->GetTotalMomentum
1285 G4LorentzVector lv1 = aParticle->Get4Moment
1286 G4LorentzVector lv(0.0,0.0,0.0,tmass);
1287
1288 lv += lv1;
1289
1290 G4ThreeVector bst = lv.boostVector();
1291
1292 // lv1.boost(-bst);
1293
1294 // G4ThreeVector p1 = lv1.vect();
1295 // G4double ptot = p1.mag();
1296
1297 G4double phi = G4UniformRand()*twopi;
1298 G4double cost = std::cos(thetaLab);
1299 G4double sint;
1300
1301 if( cost >= 1.0 )
1302 {
1303 cost = 1.0;
1304 sint = 0.0;
1305 }
1306 else if( cost <= -1.0)
1307 {
1308 cost = -1.0;
1309 sint = 0.0;
1310 }
1311 else
1312 {
1313 sint = std::sqrt((1.0-cost)*(1.0+cost));
1314 }
1315 if (verboseLevel>1)
1316 {
1317 G4cout << "cos(tlab)=" << cost << " std::
1318 }
1319 G4ThreeVector v1(sint*std::cos(phi),sint*st
1320 v1 *= plab;
1321 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),s
1322
1323 nlv1.boost(-bst);
1324
1325 G4ThreeVector np1 = nlv1.vect();
1326
1327
1328 G4double thetaCMS = np1.theta();
1329
1330 return thetaCMS;
1331 }
1332
1333 /////////////////////////////////////////////
1334 //
1335 // Test for given particle and element table
1336 // For the moment in lab system.
1337
1338 void G4NuclNuclDiffuseElastic::TestAngleTable
1339 G
1340 {
1341 fAtomicNumber = Z; // atomic number
1342 fAtomicWeight = A; // number of nucleo
1343 fNuclearRadius = CalculateNuclearRad(fAtomi
1344
1345
1346
1347 G4cout<<"G4NuclNuclDiffuseElastic::TestAngl
1348 <<Z<<"; and A = "<<A<<G4endl;
1349
1350 fElementNumberVector.push_back(fAtomicNumbe
1351
1352
1353
1354
1355 G4int i=0, j;
1356 G4double a = 0., z = theParticle->GetPDGCha
1357 G4double alpha1=0., alpha2=0., alphaMax=0.,
1358 G4double deltaL10 = 0., deltaL96 = 0., delt
1359 G4double sumL10 = 0.,sumL96 = 0.,sumAG = 0.
1360 G4double epsilon = 0.001;
1361
1362 G4Integrator<G4NuclNuclDiffuseElastic,G4dou
1363
1364 fAngleTable = new G4PhysicsTable(fEnergyBin
1365
1366 fWaveVector = partMom/hbarc;
1367
1368 G4double kR = fWaveVector*fNuclearRadiu
1369 G4double kR2 = kR*kR;
1370 G4double kRmax = 10.6; // 10.6, 18, 10.174
1371 G4double kRcoul = 1.2; // 1.4, 2.5; // on t
1372
1373 alphaMax = kRmax*kRmax/kR2;
1374
1375 if (alphaMax > 4.) alphaMax = 4.; // vmg05
1376
1377 alphaCoulomb = kRcoul*kRcoul/kR2;
1378
1379 if( z )
1380 {
1381 a = partMom/m1; // beta*gamma
1382 fBeta = a/std::sqrt(1+a*a);
1383 fZommerfeld = CalculateZommerfeld( fBet
1384 fAm = CalculateAm( partMom, fZo
1385 }
1386 G4PhysicsFreeVector* angleVector = new G4Ph
1387
1388 // G4PhysicsLogVector* angleBins = new G4P
1389
1390
1391 fAddCoulomb = false;
1392
1393 for(j = 1; j < fAngleBin; j++)
1394 {
1395 // alpha1 = angleBins->GetLowEdgeEnergy
1396 // alpha2 = angleBins->GetLowEdgeEnergy
1397
1398 alpha1 = alphaMax*(j-1)/fAngleBin;
1399 alpha2 = alphaMax*( j )/fAngleBin;
1400
1401 if( ( alpha2 > alphaCoulomb ) && z ) fAdd
1402
1403 deltaL10 = integral.Legendre10(this, &G4N
1404 deltaL96 = integral.Legendre96(this, &G4N
1405 deltaAG = integral.AdaptiveGauss(this, &
1406 alpha1
1407
1408 // G4cout<<alpha1<<"\t"<<std::sqrt(alph
1409 // <<deltaL10<<"\t"<<deltaL96<<"\t"
1410
1411 sumL10 += deltaL10;
1412 sumL96 += deltaL96;
1413 sumAG += deltaAG;
1414
1415 G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/d
1416 <<sumL10<<"\t"<<sumL96<<"\t"<<sum
1417
1418 angleVector->PutValue( j-1 , alpha1, sumL
1419 }
1420 fAngleTable->insertAt(i,angleVector);
1421 fAngleBank.push_back(fAngleTable);
1422
1423 /*
1424 // Integral over all angle range - Bad accu
1425
1426 sumL10 = integral.Legendre10(this, &G4NuclN
1427 sumL96 = integral.Legendre96(this, &G4NuclN
1428 sumAG = integral.AdaptiveGauss(this, &G4Nu
1429 0., al
1430 G4cout<<G4endl;
1431 G4cout<<alpha2<<"\t"<<std::sqrt(alpha2)/deg
1432 <<sumL10<<"\t"<<sumL96<<"\t"<<sum
1433 */
1434 return;
1435 }
1436
1437 /////////////////////////////////////////////
1438 //
1439 //
1440
1441 G4double G4NuclNuclDiffuseElastic::GetLegend
1442 {
1443 G4double legPol, epsilon = 1.e-6;
1444 G4double x = std::cos(theta);
1445
1446 if ( n < 0 ) legPol = 0.;
1447 else if( n == 0 ) legPol = 1.;
1448 else if( n == 1 ) legPol = x;
1449 else if( n == 2 ) legPol = (3.*x*x-1.)/2.;
1450 else if( n == 3 ) legPol = (5.*x*x*x-3.*x)/
1451 else if( n == 4 ) legPol = (35.*x*x*x*x-30.
1452 else if( n == 5 ) legPol = (63.*x*x*x*x*x-7
1453 else if( n == 6 ) legPol = (231.*x*x*x*x*x*
1454 else
1455 {
1456 // legPol = ( (2*n-1)*x*GetLegendrePol(n-
1457
1458 legPol = std::sqrt( 2./(n*CLHEP::pi*std::
1459 }
1460 return legPol;
1461 }
1462
1463 /////////////////////////////////////////////
1464 //
1465 //
1466
1467 G4complex G4NuclNuclDiffuseElastic::GetErfCom
1468 {
1469 G4int n;
1470 G4double n2, cofn, shny, chny, fn, gn;
1471
1472 G4double x = z.real();
1473 G4double y = z.imag();
1474
1475 G4double outRe = 0., outIm = 0.;
1476
1477 G4double twox = 2.*x;
1478 G4double twoxy = twox*y;
1479 G4double twox2 = twox*twox;
1480
1481 G4double cof1 = G4Exp(-x*x)/CLHEP::pi;
1482
1483 G4double cos2xy = std::cos(twoxy);
1484 G4double sin2xy = std::sin(twoxy);
1485
1486 G4double twoxcos2xy = twox*cos2xy;
1487 G4double twoxsin2xy = twox*sin2xy;
1488
1489 for( n = 1; n <= nMax; n++)
1490 {
1491 n2 = n*n;
1492
1493 cofn = G4Exp(-0.5*n2)/(n2+twox2); // /(n
1494
1495 chny = std::cosh(n*y);
1496 shny = std::sinh(n*y);
1497
1498 fn = twox - twoxcos2xy*chny + n*sin2xy*s
1499 gn = twoxsin2xy*chny + n*cos2xy*s
1500
1501 fn *= cofn;
1502 gn *= cofn;
1503
1504 outRe += fn;
1505 outIm += gn;
1506 }
1507 outRe *= 2*cof1;
1508 outIm *= 2*cof1;
1509
1510 if(std::abs(x) < 0.0001)
1511 {
1512 outRe += GetErf(x);
1513 outIm += cof1*y;
1514 }
1515 else
1516 {
1517 outRe += GetErf(x) + cof1*(1-cos2xy)/twox
1518 outIm += cof1*sin2xy/twox;
1519 }
1520 return G4complex(outRe, outIm);
1521 }
1522
1523
1524 /////////////////////////////////////////////
1525 //
1526 //
1527
1528 G4complex G4NuclNuclDiffuseElastic::GetErfInt
1529 {
1530 G4double outRe, outIm;
1531
1532 G4double x = z.real();
1533 G4double y = z.imag();
1534 fReZ = x;
1535
1536 G4Integrator<G4NuclNuclDiffuseElastic,G4dou
1537
1538 outRe = integral.Legendre96(this,&G4NuclNuc
1539 outIm = integral.Legendre96(this,&G4NuclNuc
1540
1541 outRe *= 2./std::sqrt(CLHEP::pi);
1542 outIm *= 2./std::sqrt(CLHEP::pi);
1543
1544 outRe += GetErf(x);
1545
1546 return G4complex(outRe, outIm);
1547 }
1548
1549 /////////////////////////////////////////////
1550 //
1551 //
1552
1553
1554 G4complex G4NuclNuclDiffuseElastic::GammaLess
1555 {
1556 G4double sinThetaR = 2.*fHalfRutThetaT
1557 G4double cosHalfThetaR2 = 1./(1. + fHalfRut
1558
1559 G4double u = std::sqrt(0.5*fPr
1560 G4double kappa = u/std::sqrt(CLHEP
1561 G4double dTheta = theta - fRutherfo
1562 u *= dTheta;
1563 G4double u2 = u*u;
1564 G4double u2m2p3 = u2*2./3.;
1565
1566 G4complex im = G4complex(0.,1.);
1567 G4complex order = G4complex(u,u);
1568 order /= std::sqrt(2.);
1569
1570 G4complex gamma = CLHEP::pi*kappa*G
1571 G4complex a0 = 0.5*(1. + 4.*(1.+
1572 G4complex a1 = 0.5*(1. + 2.*(1.+
1573 G4complex out = gamma*(1. - a1*dT
1574
1575 return out;
1576 }
1577
1578 /////////////////////////////////////////////
1579 //
1580 //
1581
1582 G4complex G4NuclNuclDiffuseElastic::GammaMore
1583 {
1584 G4double sinThetaR = 2.*fHalfRutThetaT
1585 G4double cosHalfThetaR2 = 1./(1. + fHalfRut
1586
1587 G4double u = std::sqrt(0.5*fPr
1588 G4double kappa = u/std::sqrt(CLHEP
1589 G4double dTheta = theta - fRutherfo
1590 u *= dTheta;
1591 G4double u2 = u*u;
1592 G4double u2m2p3 = u2*2./3.;
1593
1594 G4complex im = G4complex(0.,1.);
1595 G4complex order = G4complex(u,u);
1596 order /= std::sqrt(2.);
1597 G4complex gamma = CLHEP::pi*kappa*G
1598 G4complex a0 = 0.5*(1. + 4.*(1.+
1599 G4complex a1 = 0.5*(1. + 2.*(1.+
1600 G4complex out = -gamma*(1. - a1*d
1601
1602 return out;
1603 }
1604
1605 /////////////////////////////////////////////
1606 //
1607 //
1608
1609 G4complex G4NuclNuclDiffuseElastic::Amplitud
1610 {
1611 G4double kappa = std::sqrt(0.5*fProfileLamb
1612 G4complex out = G4complex(kappa/fWaveVector
1613
1614 out *= PhaseNear(theta);
1615
1616 if( theta <= fRutherfordTheta )
1617 {
1618 out *= GammaLess(theta) + ProfileNear(the
1619 // out *= GammaMore(theta) + ProfileNear(
1620 out += CoulombAmplitude(theta);
1621 }
1622 else
1623 {
1624 out *= GammaMore(theta) + ProfileNear(the
1625 // out *= GammaLess(theta) + ProfileNear(
1626 }
1627 return out;
1628 }
1629
1630 /////////////////////////////////////////////
1631 //
1632 //
1633
1634 G4complex G4NuclNuclDiffuseElastic::Amplitude
1635 {
1636 G4double sinThetaR = 2.*fHalfRutThetaTg/(1
1637 G4double dTheta = 0.5*(theta - fRutherf
1638 G4double sindTheta = std::sin(dTheta);
1639 G4double persqrt2 = std::sqrt(0.5);
1640
1641 G4complex order = G4complex(persqrt2,pe
1642 order *= std::sqrt(0.5*fProfil
1643 // order *= std::sqrt(0.5*fPro
1644
1645 G4complex out;
1646
1647 if ( theta <= fRutherfordTheta )
1648 {
1649 out = 1. - 0.5*GetErfcInt(-order)*Profile
1650 }
1651 else
1652 {
1653 out = 0.5*GetErfcInt(order)*ProfileNear(t
1654 }
1655
1656 out *= CoulombAmplitude(theta);
1657
1658 return out;
1659 }
1660
1661 /////////////////////////////////////////////
1662 //
1663 //
1664
1665 G4complex G4NuclNuclDiffuseElastic::Amplitud
1666 {
1667 G4int n;
1668 G4double T12b, b, b2; // cosTheta = std::co
1669 G4complex out = G4complex(0.,0.), shiftC, s
1670 G4complex im = G4complex(0.,1.);
1671
1672 for( n = 0; n < fMaxL; n++)
1673 {
1674 shiftC = std::exp( im*2.*CalculateCoulomb
1675 // b = ( fZommerfeld + std::sqrt( fZommer
1676 b = ( std::sqrt( G4double(n*(n+1)) ) )/fW
1677 b2 = b*b;
1678 T12b = fSumSigma*G4Exp(-b2/fNuclearRadius
1679 shiftN = std::exp( -0.5*(1.-im*fEtaRatio)
1680 out += (2.*n+1.)*shiftC*shiftN*GetLegend
1681 }
1682 out /= 2.*im*fWaveVector;
1683 out += CoulombAmplitude(theta);
1684 return out;
1685 }
1686
1687
1688 /////////////////////////////////////////////
1689 //
1690 //
1691
1692 G4complex G4NuclNuclDiffuseElastic::Amplitud
1693 {
1694 G4int n;
1695 G4double T12b, a, aTemp, b2, sinThetaH = st
1696 G4double sinThetaH2 = sinThetaH*sinThetaH;
1697 G4complex out = G4complex(0.,0.);
1698 G4complex im = G4complex(0.,1.);
1699
1700 a = -fSumSigma/CLHEP::twopi/fNuclearRadius
1701 b2 = fWaveVector*fWaveVector*fNuclearRadius
1702
1703 aTemp = a;
1704
1705 for( n = 1; n < fMaxL; n++)
1706 {
1707 T12b = aTemp*G4Exp(-b2/n)/n;
1708 aTemp *= a;
1709 out += T12b;
1710 G4cout<<"out = "<<out<<G4endl;
1711 }
1712 out *= -4.*im*fWaveVector/CLHEP::pi;
1713 out += CoulombAmplitude(theta);
1714 return out;
1715 }
1716
1717
1718 /////////////////////////////////////////////
1719 //
1720 // Test for given particle and element table
1721 // For the partMom in CMS.
1722
1723 void G4NuclNuclDiffuseElastic::InitParameters
1724 G4d
1725 {
1726 fAtomicNumber = Z; // atomic number
1727 fAtomicWeight = A; // number of nucleo
1728
1729 fNuclearRadius2 = CalculateNuclearRad(fAtom
1730 G4double A1 = G4double( theParticle->Ge
1731 fNuclearRadius1 = CalculateNuclearRad(A1);
1732 // fNuclearRadius = std::sqrt(fNuclearRadiu
1733 fNuclearRadius = fNuclearRadius1 + fNuclear
1734
1735 G4double a = 0.;
1736 G4double z = theParticle->GetPDGCharge();
1737 G4double m1 = theParticle->GetPDGMass();
1738
1739 fWaveVector = partMom/CLHEP::hbarc;
1740
1741 G4double lambda = fCofLambda*fWaveVector*fN
1742 G4cout<<"kR = "<<lambda<<G4endl;
1743
1744 if( z )
1745 {
1746 a = partMom/m1; // beta*gamma f
1747 fBeta = a/std::sqrt(1+a*a);
1748 fZommerfeld = CalculateZommerfeld( fBeta,
1749 fRutherfordRatio = fZommerfeld/fWaveVecto
1750 fAm = CalculateAm( partMom, fZomm
1751 }
1752 G4cout<<"fZommerfeld = "<<fZommerfeld<<G4en
1753 fProfileLambda = lambda; // *std::sqrt(1.-2
1754 G4cout<<"fProfileLambda = "<<fProfileLambda
1755 fProfileDelta = fCofDelta*fProfileLambda;
1756 fProfileAlpha = fCofAlpha*fProfileLambda;
1757
1758 CalculateCoulombPhaseZero();
1759 CalculateRutherfordAnglePar();
1760
1761 return;
1762 }
1763 /////////////////////////////////////////////
1764 //
1765 // Test for given particle and element table
1766 // For the partMom in CMS.
1767
1768 void G4NuclNuclDiffuseElastic::InitDynParamet
1769 G4d
1770 {
1771 G4double a = 0.;
1772 G4double z = theParticle->GetPDGCharge();
1773 G4double m1 = theParticle->GetPDGMass();
1774
1775 fWaveVector = partMom/CLHEP::hbarc;
1776
1777 G4double lambda = fCofLambda*fWaveVector*fN
1778
1779 if( z )
1780 {
1781 a = partMom/m1; // beta*gamma f
1782 fBeta = a/std::sqrt(1+a*a);
1783 fZommerfeld = CalculateZommerfeld( fBeta,
1784 fRutherfordRatio = fZommerfeld/fWaveVecto
1785 fAm = CalculateAm( partMom, fZomm
1786 }
1787 fProfileLambda = lambda; // *std::sqrt(1.-2
1788 fProfileDelta = fCofDelta*fProfileLambda;
1789 fProfileAlpha = fCofAlpha*fProfileLambda;
1790
1791 CalculateCoulombPhaseZero();
1792 CalculateRutherfordAnglePar();
1793
1794 return;
1795 }
1796
1797
1798 /////////////////////////////////////////////
1799 //
1800 // Test for given particle and element table
1801 // For the partMom in CMS.
1802
1803 void G4NuclNuclDiffuseElastic::InitParameters
1804 G4d
1805 {
1806 fAtomicNumber = Z; // target atomic nu
1807 fAtomicWeight = A; // target number of
1808
1809 fNuclearRadius2 = CalculateNuclearRad(fAtom
1810 G4double A1 = G4double( aParticle->GetD
1811 fNuclearRadius1 = CalculateNuclearRad(A1);
1812 fNuclearRadiusSquare = fNuclearRadius1*fNuc
1813
1814
1815 G4double a = 0., kR12;
1816 G4double z = aParticle->GetDefinition()->G
1817 G4double m1 = aParticle->GetDefinition()->G
1818
1819 fWaveVector = partMom/CLHEP::hbarc;
1820
1821 G4double pN = A1 - z;
1822 if( pN < 0. ) pN = 0.;
1823
1824 G4double tN = A - Z;
1825 if( tN < 0. ) tN = 0.;
1826
1827 G4double pTkin = aParticle->GetKineticEnerg
1828 pTkin /= A1;
1829
1830
1831 fSumSigma = (Z*z+pN*tN)*GetHadronNucleonXsc
1832 (z*tN+pN*Z)*GetHadronNucleonXsc
1833
1834 G4cout<<"fSumSigma = "<<fSumSigma/CLHEP::mi
1835 G4cout<<"pi*R2 = "<<CLHEP::pi*fNuclearRadiu
1836 kR12 = fWaveVector*std::sqrt(fNuclearRadius
1837 G4cout<<"k*sqrt(R2) = "<<kR12<<" "<<G4endl;
1838 fMaxL = (G4int(kR12)+1)*4;
1839 G4cout<<"fMaxL = "<<fMaxL<<" "<<G4endl;
1840
1841 if( z )
1842 {
1843 a = partMom/m1; // beta*gamma
1844 fBeta = a/std::sqrt(1+a*a);
1845 fZommerfeld = CalculateZommerfeld( fBet
1846 fAm = CalculateAm( partMom, fZo
1847 }
1848
1849 CalculateCoulombPhaseZero();
1850
1851
1852 return;
1853 }
1854
1855
1856 /////////////////////////////////////////////
1857 //
1858 // Returns nucleon-nucleon cross-section base
1859 // data from mainly http://wwwppds.ihep.su:80
1860 // projectile nucleon is pParticle with pTkin
1861
1862 G4double
1863 G4NuclNuclDiffuseElastic::GetHadronNucleonXsc
1864
1865
1866 {
1867 G4double xsection(0), /*Delta,*/ A0, B0;
1868 G4double hpXsc(0);
1869 G4double hnXsc(0);
1870
1871
1872 G4double targ_mass = tParticle->GetPDGM
1873 G4double proj_mass = pParticle->GetPDGM
1874
1875 G4double proj_energy = proj_mass + pTkin;
1876 G4double proj_momentum = std::sqrt(pTkin*(p
1877
1878 G4double sMand = CalcMandelstamS ( proj_mas
1879
1880 sMand /= CLHEP::GeV*CLHEP::GeV; //
1881 proj_momentum /= CLHEP::GeV;
1882 proj_energy /= CLHEP::GeV;
1883 proj_mass /= CLHEP::GeV;
1884 G4double logS = G4Log(sMand);
1885
1886 // General PDG fit constants
1887
1888
1889 // fEtaRatio=Re[f(0)]/Im[f(0)]
1890
1891 if( proj_momentum >= 1.2 )
1892 {
1893 fEtaRatio = 0.13*(logS - 5.8579332)*G4Po
1894 }
1895 else if( proj_momentum >= 0.6 )
1896 {
1897 fEtaRatio = -75.5*(G4Pow::GetInstance()->
1898 (G4Pow::GetInstance()->powA(3*proj_moment
1899 }
1900 else
1901 {
1902 fEtaRatio = 15.5*proj_momentum/(27*proj_m
1903 }
1904 G4cout<<"fEtaRatio = "<<fEtaRatio<<G4endl;
1905
1906 // xsc
1907
1908 if( proj_momentum >= 10. ) // high energy:
1909 // if( proj_momentum >= 2.)
1910 {
1911 //Delta = 1.;
1912
1913 //if( proj_energy < 40. ) Delta = 0.916+0
1914
1915 //AR-12Aug2016 if( proj_momentum >= 10.)
1916 {
1917 B0 = 7.5;
1918 A0 = 100. - B0*G4Log(3.0e7);
1919
1920 xsection = A0 + B0*G4Log(proj_energy)
1921 + 103*G4Pow::GetInstance()-
1922 0.93827*0.93827,-0.165);
1923 }
1924 }
1925 else // low energy pp = nn != np
1926 {
1927 if(pParticle == tParticle) // pp or nn
1928 {
1929 if( proj_momentum < 0.73 )
1930 {
1931 hnXsc = 23 + 50*( G4Pow::GetInstanc
1932 }
1933 else if( proj_momentum < 1.05 )
1934 {
1935 hnXsc = 23 + 40*(G4Log(proj_momentu
1936 (G4Log(proj_momentum
1937 }
1938 else // if( proj_momentum < 10. )
1939 {
1940 hnXsc = 39.0 +
1941 75*(proj_momentum - 1.2)/(G4Pow
1942 }
1943 xsection = hnXsc;
1944 }
1945 else // pn to be np
1946 {
1947 if( proj_momentum < 0.8 )
1948 {
1949 hpXsc = 33+30*G4Pow::GetInstance()-
1950 }
1951 else if( proj_momentum < 1.4 )
1952 {
1953 hpXsc = 33+30*G4Pow::GetInstance()-
1954 }
1955 else // if( proj_momentum < 10. )
1956 {
1957 hpXsc = 33.3+
1958 20.8*(G4Pow::GetInstance()->pow
1959 (G4Pow::GetInstance()->powA(
1960 }
1961 xsection = hpXsc;
1962 }
1963 }
1964 xsection *= CLHEP::millibarn; // parametris
1965 G4cout<<"xsection = "<<xsection/CLHEP::mill
1966 return xsection;
1967 }
1968
1969 /////////////////////////////////////////////
1970 //
1971 // The ratio el/ruth for Fresnel smooth nucle
1972
1973 G4double G4NuclNuclDiffuseElastic::GetRatioGe
1974 {
1975 G4double sinThetaR = 2.*fHalfRutThetaTg/(1
1976 G4double dTheta = 0.5*(theta - fRutherf
1977 G4double sindTheta = std::sin(dTheta);
1978
1979 G4double prof = Profile(theta);
1980 G4double prof2 = prof*prof;
1981 // G4double profmod = std::abs(prof);
1982 G4double order = std::sqrt(fProfileLam
1983
1984 order = std::abs(order); // since sin chang
1985 // G4cout<<"order = "<<order<<G4endl;
1986
1987 G4double cosFresnel = GetCint(order);
1988 G4double sinFresnel = GetSint(order);
1989
1990 G4double out;
1991
1992 if ( theta <= fRutherfordTheta )
1993 {
1994 out = 1. + 0.5*( (0.5-cosFresnel)*(0.5-c
1995 out += ( cosFresnel + sinFresnel - 1. )*p
1996 }
1997 else
1998 {
1999 out = 0.5*( (0.5-cosFresnel)*(0.5-cosFres
2000 }
2001
2002 return out;
2003 }
2004
2005 /////////////////////////////////////////////
2006 //
2007 // For the calculation of arg Gamma(z) one ne
2008 // of ln(Gamma(z))
2009
2010 G4complex G4NuclNuclDiffuseElastic::GammaLoga
2011 {
2012 const G4double cof[6] = { 76.18009172947146
2013 24.0140982408309
2014 0.1208650973866
2015 G4int j;
2016 G4complex z = zz - 1.0;
2017 G4complex tmp = z + 5.5;
2018 tmp -= (z + 0.5) * std::log(tmp);
2019 G4complex ser = G4complex(1.000000000190015
2020
2021 for ( j = 0; j <= 5; j++ )
2022 {
2023 z += 1.0;
2024 ser += cof[j]/z;
2025 }
2026 return -tmp + std::log(2.5066282746310005*s
2027 }
2028
2029 /////////////////////////////////////////////
2030 //
2031 // Bessel J0 function based on rational appro
2032 // J.F. Hart, Computer Approximations, New Yo
2033
2034 G4double G4NuclNuclDiffuseElastic::BesselJzer
2035 {
2036 G4double modvalue, value2, fact1, fact2, ar
2037
2038 modvalue = std::fabs(value);
2039
2040 if ( value < 8.0 && value > -8.0 )
2041 {
2042 value2 = value*value;
2043
2044 fact1 = 57568490574.0 + value2*(-1336259
2045 + value2*( 6516196
2046 + value2*(-1121442
2047 + value2*( 77392.3
2048 + value2*(-184.905
2049
2050 fact2 = 57568490411.0 + value2*( 1029532
2051 + value2*( 9494680
2052 + value2*(59272.64
2053 + value2*(267.8532
2054 + value2*1.0
2055
2056 bessel = fact1/fact2;
2057 }
2058 else
2059 {
2060 arg = 8.0/modvalue;
2061
2062 value2 = arg*arg;
2063
2064 shift = modvalue-0.785398164;
2065
2066 fact1 = 1.0 + value2*(-0.1098628627e-2
2067 + value2*(0.2734510407e-4
2068 + value2*(-0.2073370639e-5
2069 + value2*0.2093887211e-6
2070
2071 fact2 = -0.1562499995e-1 + value2*(0.143
2072 + value2*(-0.69
2073 + value2*(0.762
2074 - value2*0.9349
2075
2076 bessel = std::sqrt(0.636619772/modvalue)*
2077 }
2078 return bessel;
2079 }
2080
2081 /////////////////////////////////////////////
2082 //
2083 // Bessel J1 function based on rational appro
2084 // J.F. Hart, Computer Approximations, New Yo
2085
2086 G4double G4NuclNuclDiffuseElastic::BesselJone
2087 {
2088 G4double modvalue, value2, fact1, fact2, ar
2089
2090 modvalue = std::fabs(value);
2091
2092 if ( modvalue < 8.0 )
2093 {
2094 value2 = value*value;
2095
2096 fact1 = value*(72362614232.0 + value2*(-
2097 + value2*(
2098 + value2*(-
2099 + value2*(
2100 + value2*(-
2101
2102 fact2 = 144725228442.0 + value2*(2300535
2103 + value2*(1858330
2104 + value2*(99447.4
2105 + value2*(376.999
2106 + value2*1.0
2107 bessel = fact1/fact2;
2108 }
2109 else
2110 {
2111 arg = 8.0/modvalue;
2112
2113 value2 = arg*arg;
2114
2115 shift = modvalue - 2.356194491;
2116
2117 fact1 = 1.0 + value2*( 0.183105e-2
2118 + value2*(-0.3516396496e-4
2119 + value2*(0.2457520174e-5
2120 + value2*(-0.240337019e-6
2121
2122 fact2 = 0.04687499995 + value2*(-0.200269
2123 + value2*( 0.844919
2124 + value2*(-0.882289
2125 + value2*0.10578741
2126
2127 bessel = std::sqrt( 0.636619772/modvalue)
2128
2129 if (value < 0.0) bessel = -bessel;
2130 }
2131 return bessel;
2132 }
2133
2134 //
2135 //
2136 /////////////////////////////////////////////
2137