| 1 |
dlange |
1.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);
|
