Changeset 962 for trunk/source/processes/hadronic/models/coherent_elastic/src

Timestamp:

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

Author:

garnier

Message:

update processes

Location:

trunk/source/processes/hadronic/models/coherent_elastic/src

Files:

• G4ChargeExchange.cc (modified) (10 diffs)
• G4ChargeExchangeProcess.cc (modified) (11 diffs)
• G4DiffuseElastic.cc (modified) (1 diff)
• G4ElasticHadrNucleusHE.cc (modified) (1 diff)
• G4HadronElastic.cc (modified) (13 diffs)
• G4UHadronElasticProcess.cc (modified) (4 diffs)

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

@@ -25 +25 @@
 //
 //
-// $Id: G4ChargeExchange.cc,v 1.11 2007/05/25 17:46:52 dennis Exp $
-// GEANT4 tag $Name:  $
+// $Id: G4ChargeExchange.cc,v 1.14 2008/12/18 13:01:48 gunter Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 //
@@ -41 +41 @@
 #include "G4ParticleDefinition.hh"
 #include "G4IonTable.hh"
-#include "G4QElasticCrossSection.hh"
-#include "G4VQCrossSection.hh"
-#include "G4ElasticHadrNucleusHE.hh"
 #include "Randomize.hh"
-#include "G4HadronElastic.hh"
-
-
-G4ChargeExchange::G4ChargeExchange(G4HadronElastic* hel, G4double elim,
-                                   G4double ehigh)
-: G4HadronicInteraction("G4ChargeExchange"),
-  fElastic(hel),
-  native(false),
-  ekinlim(elim),
-  ekinhigh(ehigh)
+#include "G4NucleiProperties.hh"
+
+G4ChargeExchange::G4ChargeExchange() : G4HadronicInteraction("Charge Exchange")
 {
   SetMinEnergy( 0.0*GeV );
   SetMaxEnergy( 100.*TeV );
-  ekinlow = 19.0*MeV;
-
-  verboseLevel= 0;
-  if(!fElastic) {
-    native = true;
-    fElastic = new G4HadronElastic();
-  }
-  qCManager   = fElastic->GetCS();
-  hElastic    = fElastic->GetHElastic();
+
+  lowEnergyRecoilLimit = 100.*keV;  
+  lowestEnergyLimit    = 1.*MeV;  
 
   theProton   = G4Proton::Proton();
@@ -100 +84 @@
 
 G4ChargeExchange::~G4ChargeExchange()
-{
-  if(native) delete fElastic;
-}
+{}
 
 G4HadFinalState* G4ChargeExchange::ApplyYourself(
@@ -117 +99 @@
   G4int A = static_cast<G4int>(aTarget+0.5);
 
-  if(ekin == 0.0 || A < 3) {
+  if(ekin <= lowestEnergyLimit || A < 3) {
     theParticleChange.SetEnergyChange(ekin);
     theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
@@ -133 +115 @@
   // Scattered particle referred to axis of incident particle
   const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
-  G4double m1 = theParticle->GetPDGMass();
 
   G4int N = A - Z;
@@ -145 +126 @@
   G4ParticleDefinition * theDef = 0;
 
-  if (Z == 1 && A == 3) theDef = theT;
-  else if (Z == 2 && A == 3) theDef = theHe3;
-  else if (Z == 2 && A == 4) theDef = theA;
-  else theDef = G4ParticleTable::GetParticleTable()->FindIon(Z,A,0,Z);
-
-  G4double m2 = theDef->GetPDGMass();
+  G4double m2 = G4NucleiProperties::GetNuclearMass((G4double)A, (G4double)Z);
   G4LorentzVector lv1 = aParticle->Get4Momentum();
   G4LorentzVector lv0(0.0,0.0,0.0,m2);
@@ -259 +235 @@
 
   G4double etot = lv0.e() + lv1.e();
-  if(etot < m11 + m21) return &theParticleChange;
+
+  // kinematiacally impossible
+  if(etot < m11 + m21) {
+    theParticleChange.SetEnergyChange(ekin);
+    theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
+    return &theParticleChange;
+  }
 
   G4ThreeVector p1 = lv1.vect();
-  G4double e1 = 0.5*etot*(1.0 + (m21*m21 - m11*m11)/(etot*etot));
+  G4double e1 = 0.5*etot*(1.0 - (m21*m21 - m11*m11)/(etot*etot));
   //  G4double e2 = etot - e1;
   G4double ptot = std::sqrt(e1*e1 - m11*m11);
 
   G4double tmax = 4.0*ptot*ptot;
-  G4double t = 0.0;
-
-  // Choose generator
-  G4ElasticGenerator gtype = fLElastic;
-  if ((theParticle == theProton || theParticle == theNeutron) && 
-      Z <= 2 && ekin >= ekinlow) {
-    gtype = fQElastic;
-  } else {
-    if(ekin >= ekinlow)       gtype = fSWave;
-    else if(ekin >= ekinhigh) gtype = fHElastic;
-  }
-
-  // Sample t
-  if(gtype == fQElastic) {
-    if (verboseLevel > 1)
-      G4cout << "G4ChargeExchange: Z= " << Z << " N= "
-             << N << " pdg= " <<  projPDG
-             << " mom(GeV)= " << plab/GeV << "  " << qCManager << G4endl;
-    if(Z == 1 && N == 2) N = 1;
-    else if (Z == 2 && N == 1) N = 2;
-    G4double cs = qCManager->GetCrossSection(false,plab,Z,N,projPDG);
-    if(cs > 0.0) t = qCManager->GetExchangeT(Z,N,projPDG);
-    else gtype = fSWave;
-  }
-
-  if(gtype == fHElastic) {
-    t = hElastic->SampleT(theParticle,plab,Z,A);
-    if(t > tmax) gtype = fSWave;
-  }
-
-  if(gtype == fLElastic) {
-    t = GeV*GeV*fElastic->SampleT(ptot,m1,m2,aTarget);
-    if(t > tmax) gtype = fSWave;
-  }
-
-  // NaN finder
-  if(!(t < 0.0 || t >= 0.0)) {
-    if (verboseLevel > -1) {
-      G4cout << "G4ChargeExchange:WARNING: Z= " << Z << " N= " 
-             << N << " pdg= " <<  projPDG
-             << " mom(GeV)= " << plab/GeV 
-             << " the model type " << gtype;
-      if(gtype ==  fQElastic) G4cout << " CHIPS ";
-      else if(gtype ==  fLElastic) G4cout << " LElastic ";
-      else if(gtype ==  fHElastic) G4cout << " HElastic ";
-      G4cout << " t= " << t
-             << " S-wave will be sampled" 
-             << G4endl; 
-    }
-    gtype = fSWave;
-  }
-
-  if(gtype == fSWave) t = G4UniformRand()*tmax;
+  G4double g2 = GeV*GeV; 
+
+  G4double t = g2*SampleT(tmax/g2,aTarget);
 
   if(verboseLevel>1)
-    G4cout <<"type= " << gtype <<" t= " << t << " tmax= " << tmax
+    G4cout <<"## G4ChargeExchange t= " << t << " tmax= " << tmax
            << " ptot= " << ptot << G4endl;
 
@@ -328 +260 @@
   G4double phi  = G4UniformRand()*twopi;
   G4double cost = 1. - 2.0*t/tmax;
-  if(std::abs(cost) > 1.0) cost = -1.0 + 2.0*G4UniformRand();
+  if(std::abs(cost) > 1.0) cost = 1.0;
   G4double sint = std::sqrt((1.0-cost)*(1.0+cost));
 
-  if (verboseLevel > 1)
-    G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl;
+  //if (verboseLevel > 1)
+  //  G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl;
 
   G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
@@ -343 +275 @@
 
   theParticleChange.SetStatusChange(stopAndKill);
+  theParticleChange.SetEnergyChange(0.0);
   G4DynamicParticle * aSec = new G4DynamicParticle(theSecondary, nlv1);
   theParticleChange.AddSecondary(aSec);
 
   G4double erec =  nlv0.e() - m21;
+
+  //G4cout << "erec= " <<erec << " Esec= " << aSec->GetKineticEnergy() << G4endl;  
+
   if(theHyperon) {
     theParticleChange.SetLocalEnergyDeposit(erec);
@@ -352 +288 @@
     aSec->SetDefinition(theRecoil);
     aSec->SetKineticEnergy(0.0);
-  } else if(erec > ekinlim) {
+  } else if(erec > lowEnergyRecoilLimit) {
     aSec = new G4DynamicParticle(theRecoil, nlv0);
     theParticleChange.AddSecondary(aSec);
   } else {
+    if(erec < 0.0) erec = 0.0;
     theParticleChange.SetLocalEnergyDeposit(erec);
   }
   return &theParticleChange;
 }
+
+G4double G4ChargeExchange::SampleT(G4double tmax, G4double A)
+{
+  G4double aa, bb, cc, dd;
+  if (A <= 62.) {
+    aa = std::pow(A, 1.63);
+    bb = 14.5*std::pow(A, 0.66);
+    cc = 1.4*std::pow(A, 0.33);
+    dd = 10.;
+  } else {
+    aa = std::pow(A, 1.33);
+    bb = 60.*std::pow(A, 0.33);
+    cc = 0.4*std::pow(A, 0.40);
+    dd = 10.;
+  }
+  G4double x1 = (1.0 - std::exp(-tmax*bb))*aa/bb;
+  G4double x2 = (1.0 - std::exp(-tmax*dd))*cc/dd;
+  
+  G4double t;
+  G4double y = bb;
+  if(G4UniformRand()*(x1 + x2) < x2) y = dd;
+
+  do {t = -std::log(G4UniformRand())/y;} while (t > tmax);
+
+  return t;
+}
+

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);
 }
 

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

@@ -25 +25 @@
 //
 // $Id: G4DiffuseElastic.cc,v 1.20 2008/01/14 10:39:13 vnivanch Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 //

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

@@ -26 +26 @@
 //
 // $Id: G4ElasticHadrNucleusHE.cc,v 1.79 2008/01/14 10:39:13 vnivanch Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 //

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4HadronElastic.cc,v 1.56.2.1 2008/04/23 14:14:55 gcosmo Exp $
-// GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: G4HadronElastic.cc,v 1.61 2008/08/05 07:37:39 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 //
@@ -100 +100 @@
   thePionPlus = G4PionPlus::PionPlus();
   thePionMinus= G4PionMinus::PionMinus();
+
+  nnans = 0;
+  npos  = 0;
+  nneg  = 0;
+  neneg = 0;
 }
 
@@ -105 +110 @@
 {
   delete hElastic;
+  if( (nnans + npos + nneg + neneg) > 0 ) {
+    G4cout << "### G4HadronElastic destructor Warnings: ";
+    if(nnans > 0) G4cout << "###          N(nans)    = " << nnans;
+    if(npos > 0)  G4cout << "###          N(cost > 1)= " << npos;
+    if(nneg > 0)  G4cout << "###          N(cost <-1)= " << nneg;
+    if(neneg > 0) G4cout << "###          N(E < 0)=    " << neneg;
+    G4cout << "###" << G4endl;
+  }
 }
 
@@ -148 +161 @@
   G4int N = A - Z;
   G4int projPDG = theParticle->GetPDGEncoding();
-  if (verboseLevel>1) 
+  if (verboseLevel>1) {
     G4cout << "G4HadronElastic for " << theParticle->GetParticleName()
            << " PDGcode= " << projPDG << " on nucleus Z= " << Z 
            << " A= " << A << " N= " << N 
            << G4endl;
-
+  }
   G4ParticleDefinition * theDef = 0;
 
@@ -180 +193 @@
 
   // Q-elastic for p,n scattering on H and He
-  if (theParticle == theProton || theParticle == theNeutron)
+  if (theParticle == theProton || theParticle == theNeutron) {
     //     && Z <= 2 && ekin >= lowEnergyLimitQ)  
     gtype = fQElastic;
 
-  else {
+  } else {
     // S-wave for very low energy
     if(plab < plabLowLimit) gtype = fSWave;
     // HE-elastic for energetic projectile mesons
     else if(ekin >= lowEnergyLimitHE && theParticle->GetBaryonNumber() == 0) 
-      gtype = fHElastic;
+      { gtype = fHElastic; }
   }
 
@@ -212 +225 @@
 
   if(gtype == fLElastic) {
-    t = GeV*GeV*SampleT(ptot,m1,m2,aTarget);
+    G4double g2 = GeV*GeV; 
+    t = g2*SampleT(tmax/g2,m1,m2,aTarget);
   }
 
@@ -234 +248 @@
     }
     t = 0.0;
+    nnans++;
   }
 
   if(gtype == fSWave) t = G4UniformRand()*tmax;
 
-  if(verboseLevel>1)
+  if(verboseLevel>1) {
     G4cout <<"type= " << gtype <<" t= " << t << " tmax= " << tmax 
            << " ptot= " << ptot << G4endl;
-
+  }
   // Sampling in CM system
   G4double phi  = G4UniformRand()*twopi;
@@ -247 +262 @@
   G4double sint;
 
-  if( cost >= 1.0 || cost < -1 ) 
-  {
+  // problem in sampling
+  if(cost >= 1.0) {
     cost = 1.0;
     sint = 0.0;
-  }
-  else  
-  {
+    npos++;
+  } else if(cost < -1 ) {
+    /*
+    G4cout << "G4HadronElastic:WARNING: Z= " << Z << " N= " 
+           << N << " " << aParticle->GetDefinition()->GetParticleName()
+           << " mom(GeV)= " << plab/GeV 
+           << " the model type " << gtype;
+    if(gtype ==  fQElastic) G4cout << " CHIPS ";
+    else if(gtype ==  fLElastic) G4cout << " LElastic ";
+    else if(gtype ==  fHElastic) G4cout << " HElastic ";
+    G4cout << " cost= " << cost 
+           << G4endl; 
+    */
+    cost = 1.0;
+    sint = 0.0;
+    nneg++;
+
+    // normal situation
+  } else  {
     sint = std::sqrt((1.0-cost)*(1.0+cost));
   }    
-  if (verboseLevel>1) 
+  if (verboseLevel>1) {
     G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl;
-
+  }
   G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
   v1 *= ptot;
@@ -273 +304 @@
   }
   if(eFinal <= lowestEnergyLimit) {
-    if(eFinal < 0.0) {
+    if(eFinal < 0.0 && verboseLevel > 0) {
+      neneg++;
       G4cout << "G4HadronElastic WARNING ekin= " << eFinal
              << " after scattering of " 
@@ -308 +340 @@
 
 G4double 
-G4HadronElastic::SampleT(G4double, G4double, G4double, G4double atno2)
+G4HadronElastic::SampleT(G4double tmax, G4double, G4double, G4double atno2)
 {
   // G4cout << "Entering elastic scattering 2"<<G4endl;
@@ -315 +347 @@
   // parameterized functions and a Monte Carlo method to invert the CDF.
 
-  G4double ran = G4UniformRand();
+  //  G4double ran = G4UniformRand();
   G4double aa, bb, cc, dd, rr;
   if (atno2 <= 62.) {
@@ -330 +362 @@
   aa = aa/bb;
   cc = cc/dd;
+  G4double ran, t1, t2;
+  do {
+    ran = G4UniformRand();
+    t1 = -std::log(ran)/bb;
+    t2 = -std::log(ran)/dd;
+  } while(t1 > tmax || t2 > tmax);
+
   rr = (aa + cc)*ran;
+
   if (verboseLevel > 1) {
     G4cout << "DoIt: aa,bb,cc,dd,rr" << G4endl;
     G4cout << aa << " " << bb << " " << cc << " " << dd << " " << rr << G4endl;
-  }
-  G4double t1 = -std::log(ran)/bb;
-  G4double t2 = -std::log(ran)/dd;
-  if (verboseLevel > 1) {
     G4cout << "t1,Fctcos " << t1 << " " << Fctcos(t1, aa, bb, cc, dd, rr) << G4endl;
     G4cout << "t2,Fctcos " << t2 << " " << Fctcos(t2, aa, bb, cc, dd, rr) << G4endl;
@@ -348 +384 @@
               aa, bb, cc, dd, rr);
   if (verboseLevel > 1) {
-    G4cout << "From Rtmi, ier1=" << ier1 << G4endl;
+    G4cout << "From Rtmi, ier1=" << ier1 << " t= " << t << G4endl;
     G4cout << "t, Fctcos " << t << " " << Fctcos(t, aa, bb, cc, dd, rr) << G4endl;
   }

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4UHadronElasticProcess.cc,v 1.35.2.1 2008/04/23 14:14:55 gcosmo Exp $
-// GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: G4UHadronElasticProcess.cc,v 1.39 2008/10/22 08:16:40 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // Geant4 Hadron Elastic Scattering Process -- header file
@@ -58 +58 @@
   : G4HadronicProcess(pName), lowestEnergy(0.0), first(true)
 {
+  SetProcessSubType(fHadronElastic);
   AddDataSet(new G4HadronElasticDataSet);
   theProton   = G4Proton::Proton();
@@ -89 +90 @@
     if(theParticle->GetPDGCharge() != 0.0) lowestEnergy = eV;
      
-    if(verboseLevel>1 || 
-       (verboseLevel==1 && theParticle == theNeutron)) {
-      G4cout << G4endl;
+    //    if(verboseLevel>1 || 
+    //   (verboseLevel==1 && theParticle == theNeutron)) {
+    if(verboseLevel>1 && theParticle == theNeutron) {
+    //      G4cout << G4endl;
       G4cout << "G4UHadronElasticProcess for " 
              << theParticle->GetParticleName()
@@ -100 +102 @@
     } 
   }
-  store->BuildPhysicsTable(aParticleType);
+  G4HadronicProcess::BuildPhysicsTable(aParticleType);
+  //store->BuildPhysicsTable(aParticleType);
 }