Changeset 961 for trunk/source/processes/electromagnetic/standard/src/G4CoulombScatteringModel.cc

Timestamp:

Apr 6, 2009, 12:21:12 PM (15 years ago)

Author:

garnier

Message:

update processes

File:

• trunk/source/processes/electromagnetic/standard/src/G4CoulombScatteringModel.cc (modified) (6 diffs)

trunk/source/processes/electromagnetic/standard/src/G4CoulombScatteringModel.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4CoulombScatteringModel.cc,v 1.29 2007/11/09 11:45:45 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: G4CoulombScatteringModel.cc,v 1.37 2008/07/31 13:11:34 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // -------------------------------------------------------------------
@@ -44 +44 @@
 // 19.10.06 V.Ivanchenko use inheritance from G4eCoulombScatteringModel
 // 09.10.07 V.Ivanchenko reorganized methods, add cut dependence in scattering off e- 
+// 09.06.08 V.Ivanchenko SelectIsotope is moved to the base class
 //
 // Class Description:
@@ -55 +56 @@
 #include "Randomize.hh"
 #include "G4ParticleChangeForGamma.hh"
-#include "G4NistManager.hh"
 #include "G4ParticleTable.hh"
 #include "G4IonTable.hh"
@@ -64 +64 @@
 using namespace std;
 
-G4CoulombScatteringModel::G4CoulombScatteringModel(
-  G4double thetaMin, G4double thetaMax, G4bool build, 
-  G4double tlim, const G4String& nam)
-  : G4eCoulombScatteringModel(thetaMin,thetaMax,build,tlim,nam)
-{
-  theMatManager    = G4NistManager::Instance();
-  theParticleTable = G4ParticleTable::GetParticleTable();
-}
+G4CoulombScatteringModel::G4CoulombScatteringModel(const G4String& nam)
+  : G4eCoulombScatteringModel(nam)
+{}
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
@@ -84 +79 @@
                                 G4double kinEnergy, 
                                 G4double Z, 
-                                G4double A, 
+                                G4double, 
                                 G4double cutEnergy,
                                 G4double)
 {
-  if(p == particle && kinEnergy == tkin && Z == targetZ &&
-     A == targetA && cutEnergy == ecut) return nucXSection;
-
-  // Lab system
-  G4double ekin = std::max(keV, kinEnergy);
-  nucXSection = ComputeElectronXSectionPerAtom(p,ekin,Z,A,cutEnergy);
+  SetupParticle(p);
+  G4double ekin = std::max(lowEnergyLimit, kinEnergy);
+  SetupKinematic(ekin, cutEnergy);
+
+  // save lab system kinematics
+  G4double xtkin = tkin;
+  G4double xmom2 = mom2;
+  G4double xinvb = invbeta2;
 
   // CM system
   G4int iz      = G4int(Z);
-  G4double m1   = theMatManager->GetAtomicMassAmu(iz)*amu_c2;
+  G4double m2   = fNistManager->GetAtomicMassAmu(iz)*amu_c2;
   G4double etot = tkin + mass;
   G4double ptot = sqrt(mom2);
-  G4double bet  = ptot/(etot + m1);
-  G4double gam  = 1.0/sqrt((1.0 - bet)*(1.0 + bet));
-  G4double momCM= gam*(ptot - bet*etot);
-
-  //  G4cout << "ptot= " << ptot << " etot= " << etot << " beta= " 
-  //     << bet << " gam= " << gam << " Z= " << Z << " A= " << A << G4endl;
-  // G4cout << " CM. mom= " << momCM  << " m= " << m 
-  // << " m1= " << m1 << " iz= " << iz <<G4endl;
-
-  G4double momCM2 = momCM*momCM;
-  cosTetMaxNuc = std::max(cosThetaMax, 1.0 - 0.5*q2Limit/momCM2);
-  if(1.5 > targetA && p == theProton && cosTetMaxNuc < 0.0) cosTetMaxNuc = 0.0;
-  //G4cout << " ctmax= " << cosTetMaxNuc 
-  //<< " ctmin= " << cosThetaMin << G4endl;  
-
-  // Cross section in CM system 
-  if(cosTetMaxNuc < cosThetaMin) {
-    G4double effmass = mass*m1/(mass + m1);
-    G4double x1 = 1.0 - cosThetaMin;
-    G4double x2 = 1.0 - cosTetMaxNuc;
-    G4double z1 = x1 + screenZ;
-    G4double z2 = x2 + screenZ;
-    G4double d  = 1.0/formfactA;
-    G4double zn1= x1 + d;
-    G4double zn2= x2 + d;
-    nucXSection += coeff*Z*Z*chargeSquare*(1.0 +  effmass*effmass/momCM2)
-      *(1./z1 - 1./z2 + 1./zn1 - 1./zn2 + 
-        2.0*formfactA*std::log(z1*zn2/(z2*zn1)))/momCM2;
-    //G4cout << "XS: x1= " << x1 << " x2= " << x2 
-    //<< " cross= " << cross << G4endl;
-    //G4cout << "momCM2= " << momCM2 << " invbeta2= " << invbeta2 
-    //       << " coeff= " << coeff << G4endl;
-  }
-  if(nucXSection < 0.0) nucXSection = 0.0;
-  return nucXSection;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
-G4double G4CoulombScatteringModel::SelectIsotope(const G4Element* elm)
-{
-  G4double N = elm->GetN(); 
-  G4int ni   = elm->GetNumberOfIsotopes();
-  if(ni > 0) {
-    G4double* ab = elm->GetRelativeAbundanceVector();
-    G4double x = G4UniformRand();
-    G4int idx;
-    for(idx=0; idx<ni; idx++) { 
-      x -= ab[idx];
-      if (x <= 0.0) break;
-    }
-    if(idx >= ni) {
-      G4cout << "G4CoulombScatteringModel::SelectIsotope WARNING: "
-             << "abandance vector for"
-             << elm->GetName() << " is not normalised to unit" << G4endl;
-    } else {
-      N = G4double(elm->GetIsotope(idx)->GetN());
-    }
-  }
-  return N;
+
+  G4double m12  = mass*mass;
+  G4double momCM= ptot*m2/sqrt(m12 + m2*m2 + 2.0*etot*m2);
+
+  mom2 = momCM*momCM;
+  tkin = sqrt(mom2 + m12) - mass;
+  invbeta2 = 1.0 +  m12/mom2;
+
+  SetupTarget(Z, tkin);
+
+  G4double xsec = CrossSectionPerAtom();
+
+  // restore Lab system kinematics
+  tkin = xtkin;
+  mom2 = xmom2;
+  invbeta2 = xinvb;
+
+  return xsec;
 }
 
@@ -169 +124 @@
                                const G4DynamicParticle* dp,
                                G4double cutEnergy, 
-                               G4double maxEnergy)
+                               G4double)
 {
-  const G4Material* aMaterial = couple->GetMaterial();
-  const G4ParticleDefinition* p = dp->GetDefinition();
   G4double kinEnergy = dp->GetKineticEnergy();
-
-  // Select isotope and setup
-  SetupParticle(p);
-  const G4Element* elm = 
-    SelectRandomAtom(aMaterial,p,kinEnergy,cutEnergy,maxEnergy);
-  G4double Z  = elm->GetZ();
-  G4double A  = SelectIsotope(elm);
+  if(kinEnergy <= DBL_MIN) return;
+  DefineMaterial(couple);
+  SetupParticle(dp->GetDefinition());
+  G4double ekin = std::max(lowEnergyLimit, kinEnergy);
+  SetupKinematic(ekin, cutEnergy);
+
+  // Choose nucleus
+  currentElement = SelectRandomAtom(couple,particle,ekin,ecut,tkin);
+
+  G4double Z  = currentElement->GetZ();
   G4int iz    = G4int(Z);
-  G4int ia    = G4int(A + 0.5);
-
-  G4double cross = 
-    ComputeCrossSectionPerAtom(p,kinEnergy,Z,A,cutEnergy,maxEnergy);
-
-  G4double costm = cosTetMaxNuc;
-  G4double formf = formfactA;
-  if(G4UniformRand()*cross < elecXSection) {
-    costm = cosTetMaxElec;
-    formf = 0.0;
-  }
-
-  //  G4cout << "SampleSec: Ekin= " << kinEnergy << " m1= " << m1 
-  // << " Z= "<< Z << " A= " <<A<< G4endl; 
-
-  if(costm >= cosThetaMin) return; 
-
-  // kinematics in CM system
-  G4double m1   = theParticleTable->GetIonTable()->GetNucleusMass(iz, ia);
-  G4double etot = kinEnergy + mass;
+  G4int ia    = SelectIsotopeNumber(currentElement);
+  G4double m2 = theParticleTable->GetIonTable()->GetNucleusMass(iz, ia);
+
+  // CM system
+  G4double etot = tkin + mass;
   G4double ptot = sqrt(mom2);
-  G4double bet  = ptot/(etot + m1);
+
+  G4double momCM= ptot*m2/sqrt(mass*mass + m2*m2 + 2.0*etot*m2);
+  mom2 = momCM*momCM;
+  G4double m12 = mass*mass;
+  G4double eCM = sqrt(mom2 + m12);
+
+  // a correction for heavy projectile
+  G4double fm = m2/(mass + m2);
+  invbeta2 = 1.0 +  m12*fm*fm/mom2;
+
+  // sample scattering angle in CM system
+  SetupTarget(Z, eCM - mass);
+
+  G4double cost = SampleCosineTheta();
+  G4double z1   = 1.0 - cost;
+  if(z1 < 0.0) return;
+
+  G4double sint = sqrt(z1*(1.0 + cost));
+  G4double phi  = twopi * G4UniformRand();
+
+  // kinematics in the Lab system
+  G4double bet  = ptot/(etot + m2);
   G4double gam  = 1.0/sqrt((1.0 - bet)*(1.0 + bet));
-  G4double pCM  = gam*(ptot - bet*etot);
-  G4double eCM  = gam*(etot - bet*ptot);
-
-  G4double x1 = 1. - cosThetaMin + screenZ;
-  G4double x2 = 1. - costm;
-  G4double x3 = cosThetaMin - costm;
-
-  G4double grej,  z, z1; 
-  do {
-    z  = G4UniformRand()*x3;
-    z1 = (x1*x2 - screenZ*z)/(x1 + z);
-    if(z1 < 0.0) z1 = 0.0;
-    else if(z1 > 2.0) z1 = 2.0;
-    grej = 1.0/(1.0 + formf*z1);
-  } while ( G4UniformRand() > grej*grej );  
-  
-  G4double cost = 1.0 - z1;
-  G4double sint= sqrt(z1*(2.0 - z1));
-
-  G4double phi = twopi * G4UniformRand();
-
-  // projectile after scattering
-  G4double pzCM = pCM*cost;
-  G4ThreeVector v1(pCM*cos(phi)*sint,pCM*sin(phi)*sint,gam*(pzCM + bet*eCM));
+  G4double pzCM = momCM*cost;
+
+  G4ThreeVector v1(momCM*cos(phi)*sint,momCM*sin(phi)*sint,gam*(pzCM + bet*eCM));
   G4ThreeVector dir = dp->GetMomentumDirection(); 
   G4ThreeVector newDirection = v1.unit();
   newDirection.rotateUz(dir);   
   fParticleChange->ProposeMomentumDirection(newDirection);   
+
   G4double elab = gam*(eCM + bet*pzCM);
-  G4double ekin = elab - mass;
+  ekin = elab - mass;
   if(ekin < 0.0) ekin = 0.0;
-  G4double plab = sqrt(ekin*(ekin + 2.0*mass));
   fParticleChange->SetProposedKineticEnergy(ekin);
 
   // recoil
   G4double erec = kinEnergy - ekin;
-  if(erec > Z*aMaterial->GetIonisation()->GetMeanExcitationEnergy()) {
+  G4double th = 
+    std::min(recoilThreshold,
+             Z*currentElement->GetIonisation()->GetMeanExcitationEnergy());
+
+  if(erec > th) {
     G4ParticleDefinition* ion = theParticleTable->FindIon(iz, ia, 0, iz);
+    G4double plab = sqrt(ekin*(ekin + 2.0*mass));
     G4ThreeVector p2 = (ptot*dir - plab*newDirection).unit();
     G4DynamicParticle* newdp  = new G4DynamicParticle(ion, p2, erec);