Changeset 962 for trunk/source/processes/hadronic/models/coherent_elastic/src/G4ChargeExchangeProcess.cc

Timestamp:

Apr 6, 2009, 12:30:29 PM (16 years ago)

Author:

garnier

Message:

update processes

File:

• trunk/source/processes/hadronic/models/coherent_elastic/src/G4ChargeExchangeProcess.cc (modified) (11 diffs)

trunk/source/processes/hadronic/models/coherent_elastic/src/G4ChargeExchangeProcess.cc

@@ -25 +25 @@
 //
 //
-// $Id: G4ChargeExchangeProcess.cc,v 1.9 2007/01/30 10:23:26 vnivanch Exp $
-// GEANT4 tag $Name:  $
+// $Id: G4ChargeExchangeProcess.cc,v 1.15 2008/11/27 16:43:00 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 //
-// Geant4 Hadron Elastic Scattering Process -- header file
+// Geant4 Hadron Charge Exchange Process -- source file
 //
 // Created 21 April 2006 V.Ivanchenko
@@ -45 +45 @@
 #include "G4CrossSectionDataStore.hh"
 #include "G4HadronElasticDataSet.hh"
-#include "G4VQCrossSection.hh"
-#include "G4QElasticCrossSection.hh"
-#include "G4QCHIPSWorld.hh"
 #include "G4Element.hh"
 #include "G4ElementVector.hh"
@@ -53 +50 @@
 #include "G4Neutron.hh"
 #include "G4Proton.hh"
-#include "G4HadronElastic.hh"
 #include "G4PhysicsLinearVector.hh"
 
@@ -59 +55 @@
   : G4HadronicProcess(procName), first(true)
 {
-  thEnergy = 19.*MeV;
+  SetProcessSubType(fChargeExchange);
+  thEnergy = 20.*MeV;
   verboseLevel= 1;
-  qCManager = 0;
   AddDataSet(new G4HadronElasticDataSet);
   theProton   = G4Proton::Proton();
@@ -99 +95 @@
 }
 
-void G4ChargeExchangeProcess::SetQElasticCrossSection(G4VQCrossSection* p)
-{
-  qCManager = p;
-}
-
 void G4ChargeExchangeProcess::
 BuildPhysicsTable(const G4ParticleDefinition& aParticleType)
@@ -120 +111 @@
 
       G4double F[n] = {0.33,0.27,0.29,0.31,0.27,0.18,0.13,0.1,0.09,0.07};
-      factors = new G4PhysicsLinearVector(0.0,1.8*GeV,n);
+      factors = new G4PhysicsLinearVector(0.0,2.0*GeV,n);
       for(size_t i=0; i<n; i++) {factors->PutValue(i,F[i]);}
 
@@ -126 +117 @@
 
       G4double F[n] = {0.50,0.45,0.40,0.35,0.30,0.25,0.06,0.04,0.005,0.0};
-      factors = new G4PhysicsLinearVector(0.0,3.6*GeV,n);
+      factors = new G4PhysicsLinearVector(0.0,4.0*GeV,n);
       for(size_t i=0; i<n; i++) {factors->PutValue(i,F[i]);}
     }
+    //factors->SetSpline(true);
 
     if(verboseLevel>1)
@@ -135 +127 @@
              << G4endl;
   }
-  store->BuildPhysicsTable(aParticleType);
-}
-
-G4double G4ChargeExchangeProcess::GetMeanFreePath(const G4Track& track,
-                                                  G4double,
-                                                  G4ForceCondition* cond)
-{
-  *cond = NotForced;
-  const G4DynamicParticle* dp = track.GetDynamicParticle();
-  const G4Material* material = track.GetMaterial();
-  cross = 0.0;
-  G4double x = DBL_MAX;
-
-  // The process is effective only above the threshold
-  if(dp->GetKineticEnergy() < thEnergy) return x;
-
-  // Compute cross sesctions
-  const G4ElementVector* theElementVector = material->GetElementVector();
-  const G4double* theAtomNumDensityVector = material->GetVecNbOfAtomsPerVolume();
-  G4double temp = material->GetTemperature();
-  G4int nelm   = material->GetNumberOfElements();
-  if(verboseLevel>1)
-    G4cout << "G4ChargeExchangeProcess get mfp for "
-           << theParticle->GetParticleName()
-           << "  p(GeV)= " << dp->GetTotalMomentum()/GeV
-           << " in " << material->GetName()
-           << G4endl;
-  for (G4int i=0; i<nelm; i++) {
-    const G4Element* elm = (*theElementVector)[i];
-    G4double x = GetMicroscopicCrossSection(dp, elm, temp);
-    cross += theAtomNumDensityVector[i]*x;
-    xsec[i] = cross;
-  }
-  if(verboseLevel>1)
-    G4cout << "G4ChargeExchangeProcess cross(1/mm)= " << cross
-           << "  E(MeV)= " << dp->GetKineticEnergy()
-           << "  " << theParticle->GetParticleName()
-           << "  in " << material->GetName()
-           << G4endl;
-  if(cross > DBL_MIN) x = 1./cross;
-
-  return x;
+  G4HadronicProcess::BuildPhysicsTable(aParticleType);
 }
 
@@ -188 +139 @@
   G4int iz = G4int(Z);
   G4double x = 0.0;
-  if(iz == 1) return x;
+
+  // The process is effective only above the threshold
+  if(iz == 1 || dp->GetKineticEnergy() < thEnergy) return x;
 
   if(verboseLevel>1)
@@ -195 +148 @@
            << G4endl;
   x = store->GetCrossSection(dp, elm, temp);
-
-  // NaN finder
-  if(!(x < 0.0 || x >= 0.0)) {
-    if (verboseLevel > -1) {
-      G4cout << "G4ChargeExchangeProcess WARNING: Z= " << iz  
-             << " pdg= " <<  pPDG
-             << " mom(GeV)= " << dp->GetTotalMomentum()/GeV 
-             << " cross= " << x 
-             << " set to zero"
-             << G4endl; 
-    }
-    x = 0.0;
-  }
 
   if(verboseLevel>1)
@@ -217 +157 @@
   G4bool b;
   G4double A = elm->GetN();
-  x *= factors->GetValue(dp->GetTotalMomentum(), b)/std::pow(A, 0.42);
+  G4double ptot = dp->GetTotalMomentum();
+  x *= factors->GetValue(ptot, b)/std::pow(A, 0.42);
   if(theParticle == thePiPlus || theParticle == theProton ||
      theParticle == theKPlus  || theParticle == theANeutron)
-    x *= (1.0 - Z/A);
+    { x *= (1.0 - Z/A); }
 
   else if(theParticle == thePiMinus || theParticle == theNeutron ||
           theParticle == theKMinus  || theParticle == theAProton)
-    x *= Z/A;
+    { x *= Z/A; }
+
+  if(theParticle->GetPDGMass() < GeV) {
+    if(ptot > 2.*GeV) x *= 4.0*GeV*GeV/(ptot*ptot);
+  }
+
+  if(verboseLevel>1) 
+    G4cout << "Corrected cross(mb)= " << x/millibarn << G4endl;
 
   return x;
-}
-
-G4VParticleChange* G4ChargeExchangeProcess::PostStepDoIt(
-                                  const G4Track& track,
-                                  const G4Step& step)
-{
-  G4ForceCondition* cn = 0;
-  aParticleChange.Initialize(track);
-  G4double mfp = GetMeanFreePath(track, 0.0, cn);
-  if(mfp == DBL_MAX) return G4VDiscreteProcess::PostStepDoIt(track,step);
-
-  G4double kineticEnergy = track.GetKineticEnergy();
-  G4Material* material = track.GetMaterial();
-
-  // Select element
-  const G4ElementVector* theElementVector = material->GetElementVector();
-  G4Element* elm = (*theElementVector)[0];
-  G4int nelm = material->GetNumberOfElements() - 1;
-  if (nelm > 0) {
-    G4double x = G4UniformRand()*cross;
-    G4int i = -1;
-    do {i++;} while (x > xsec[i] && i < nelm);
-    elm = (*theElementVector)[i];
-  }
-  G4double Z = elm->GetZ();
-  G4double A = G4double(G4int(elm->GetN()+0.5));
-
-  // Select isotope
-  G4IsotopeVector* isv = elm->GetIsotopeVector(); 
-  G4int ni = 0;
-  if(isv) ni = isv->size();
-
-  if(ni == 1) { 
-    A = G4double((*isv)[0]->GetN());
-  } else if(ni > 1) {
-
-    G4double* ab = elm->GetRelativeAbundanceVector();
-    G4int j = -1;
-    ni--;
-    G4double y = G4UniformRand();
-    do {
-      j++;
-      y -= ab[j];
-    } while (y > 0.0 && j < ni);
-    A = G4double((*isv)[j]->GetN());
-  } 
-  G4HadronicInteraction* hadi =
-    ChooseHadronicInteraction( kineticEnergy, material, elm);
-
-  // Initialize the hadronic projectile from the track
-  G4HadProjectile thePro(track);
-  if(verboseLevel>1)
-    G4cout << "G4ChargeExchangeProcess::PostStepDoIt for "
-           << theParticle->GetParticleName()
-           << " Target Z= " << Z
-           << " A= " << A << G4endl;
-  targetNucleus.SetParameters(A, Z);
-
-  aParticleChange.Initialize(track);
-  G4HadFinalState* result = hadi->ApplyYourself(thePro, targetNucleus);
-  G4ThreeVector indir = track.GetMomentumDirection();
-  G4int nsec = result->GetNumberOfSecondaries();
-
-  if(verboseLevel>1)
-    G4cout << "Efin= " << result->GetEnergyChange()
-           << " de= " << result->GetLocalEnergyDeposit()
-           << " nsec= " << nsec
-           << G4endl;
-
-
-  if(nsec > 0) {
-    aParticleChange.ProposeEnergy(0.0);
-    aParticleChange.ProposeTrackStatus(fStopAndKill);
-    aParticleChange.ProposeLocalEnergyDeposit(result->GetLocalEnergyDeposit());
-    aParticleChange.SetNumberOfSecondaries(nsec);
-    for(G4int j=0; j<nsec; j++) {
-      G4DynamicParticle* p = result->GetSecondary(j)->GetParticle();
-      G4ThreeVector pdir = p->GetMomentumDirection();
-      // G4cout << "recoil " << pdir << G4endl;
-      pdir = pdir.rotateUz(indir);
-      // G4cout << "recoil rotated " << pdir << G4endl;
-      p->SetMomentumDirection(pdir);
-      aParticleChange.AddSecondary(p);
-    }
-  }
-  result->Clear();
-
-  return G4VDiscreteProcess::PostStepDoIt(track,step);
 }