| 1 |
--- source/processes/hadronic/models/coherent_elastic/src/G4HadronElastic.cc 2007-12-11 12:39:35.000000000 +0100
|
| 2 |
+++ source/processes/hadronic/models/coherent_elastic/src/G4HadronElastic.cc 2008-04-17 15:48:24.000000000 +0200
|
| 3 |
@@ -23,8 +23,8 @@
|
| 4 |
// * acceptance of all terms of the Geant4 Software license. *
|
| 5 |
// ********************************************************************
|
| 6 |
//
|
| 7 |
-// $Id: G4HadronElastic.cc,v 1.56 2007/12/11 11:39:35 vnivanch Exp $
|
| 8 |
-// GEANT4 tag $Name: geant4-09-01-patch-01 $
|
| 9 |
+// $Id: G4HadronElastic.cc,v 1.58 2008/04/10 14:57:23 vnivanch Exp $
|
| 10 |
+// GEANT4 tag $Name: $
|
| 11 |
//
|
| 12 |
//
|
| 13 |
// Physics model class G4HadronElastic (derived from G4LElastic)
|
| 14 |
@@ -87,7 +87,7 @@
|
| 15 |
lowEnergyRecoilLimit = 100.*keV;
|
| 16 |
lowEnergyLimitQ = 0.0*GeV;
|
| 17 |
lowEnergyLimitHE = 1.0*GeV;
|
| 18 |
- lowestEnergyLimit= 0.0*keV;
|
| 19 |
+ lowestEnergyLimit= 1.e-6*eV;
|
| 20 |
plabLowLimit = 20.0*MeV;
|
| 21 |
|
| 22 |
qCManager = G4QElasticCrossSection::GetPointer();
|
| 23 |
@@ -133,12 +133,12 @@
|
| 24 |
G4double zTarget = targetNucleus.GetZ();
|
| 25 |
|
| 26 |
G4double plab = aParticle->GetTotalMomentum();
|
| 27 |
- if (verboseLevel >1)
|
| 28 |
+ if (verboseLevel >1) {
|
| 29 |
G4cout << "G4HadronElastic::DoIt: Incident particle plab="
|
| 30 |
<< plab/GeV << " GeV/c "
|
| 31 |
<< " ekin(MeV) = " << ekin/MeV << " "
|
| 32 |
<< aParticle->GetDefinition()->GetParticleName() << G4endl;
|
| 33 |
-
|
| 34 |
+ }
|
| 35 |
// Scattered particle referred to axis of incident particle
|
| 36 |
const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
|
| 37 |
G4double m1 = theParticle->GetPDGMass();
|
| 38 |
@@ -195,10 +195,11 @@
|
| 39 |
// Sample t
|
| 40 |
//
|
| 41 |
if(gtype == fQElastic) {
|
| 42 |
- if (verboseLevel >1)
|
| 43 |
+ if (verboseLevel >1) {
|
| 44 |
G4cout << "G4HadronElastic: Z= " << Z << " N= "
|
| 45 |
<< N << " pdg= " << projPDG
|
| 46 |
<< " mom(GeV)= " << plab/GeV << " " << qCManager << G4endl;
|
| 47 |
+ }
|
| 48 |
if(Z == 1 && N == 2) N = 1;
|
| 49 |
else if(Z == 2 && N == 1) N = 2;
|
| 50 |
G4double cs = qCManager->GetCrossSection(false,plab,Z,N,projPDG);
|
| 51 |
@@ -245,16 +246,11 @@
|
| 52 |
G4double cost = 1. - 2.0*t/tmax;
|
| 53 |
G4double sint;
|
| 54 |
|
| 55 |
- if( cost >= 1.0 )
|
| 56 |
+ if( cost >= 1.0 || cost < -1 )
|
| 57 |
{
|
| 58 |
cost = 1.0;
|
| 59 |
sint = 0.0;
|
| 60 |
}
|
| 61 |
- else if( cost <= -1.0)
|
| 62 |
- {
|
| 63 |
- cost = -1.0;
|
| 64 |
- sint = 0.0;
|
| 65 |
- }
|
| 66 |
else
|
| 67 |
{
|
| 68 |
sint = std::sqrt((1.0-cost)*(1.0+cost));
|
| 69 |
@@ -269,11 +265,13 @@
|
| 70 |
nlv1.boost(bst);
|
| 71 |
|
| 72 |
G4double eFinal = nlv1.e() - m1;
|
| 73 |
- if (verboseLevel > 1)
|
| 74 |
+ if (verboseLevel > 1) {
|
| 75 |
G4cout << "Scattered: "
|
| 76 |
<< nlv1<<" m= " << m1 << " ekin(MeV)= " << eFinal
|
| 77 |
<< " Proj: 4-mom " << lv1
|
| 78 |
<<G4endl;
|
| 79 |
+ }
|
| 80 |
+ if(eFinal <= lowestEnergyLimit) {
|
| 81 |
if(eFinal < 0.0) {
|
| 82 |
G4cout << "G4HadronElastic WARNING ekin= " << eFinal
|
| 83 |
<< " after scattering of "
|
| 84 |
@@ -281,20 +279,22 @@
|
| 85 |
<< " p(GeV/c)= " << plab
|
| 86 |
<< " on " << theDef->GetParticleName()
|
| 87 |
<< G4endl;
|
| 88 |
- eFinal = 0.0;
|
| 89 |
- nlv1.setE(m1);
|
| 90 |
}
|
| 91 |
+ theParticleChange.SetEnergyChange(0.0);
|
| 92 |
+ nlv1 = G4LorentzVector(0.0,0.0,0.0,m1);
|
| 93 |
|
| 94 |
+ } else {
|
| 95 |
theParticleChange.SetMomentumChange(nlv1.vect().unit());
|
| 96 |
theParticleChange.SetEnergyChange(eFinal);
|
| 97 |
+ }
|
| 98 |
|
| 99 |
G4LorentzVector nlv0 = lv - nlv1;
|
| 100 |
G4double erec = nlv0.e() - m2;
|
| 101 |
- if (verboseLevel > 1)
|
| 102 |
+ if (verboseLevel > 1) {
|
| 103 |
G4cout << "Recoil: "
|
| 104 |
<< nlv0<<" m= " << m2 << " ekin(MeV)= " << erec
|
| 105 |
<<G4endl;
|
| 106 |
-
|
| 107 |
+ }
|
| 108 |
if(erec > lowEnergyRecoilLimit) {
|
| 109 |
G4DynamicParticle * aSec = new G4DynamicParticle(theDef, nlv0);
|
| 110 |
theParticleChange.AddSecondary(aSec);
|
