Changeset 1228 for trunk/source/processes/hadronic/models/abrasion

Timestamp:

Jan 8, 2010, 11:56:51 AM (15 years ago)

Author:

garnier

Message:

update geant4.9.3 tag

Location:

trunk/source/processes/hadronic/models/abrasion

Files:

• include/G4WilsonAbrasionModel.hh (modified) (3 diffs)
• src/G4WilsonAbrasionModel.cc (modified) (7 diffs)

trunk/source/processes/hadronic/models/abrasion/include/G4WilsonAbrasionModel.hh

@@ -40 +40 @@
 // MODULE:              G4WilsonAbrasionModel.hh
 //
-// Version:             B.1
-// Date:                15/04/04
+// Version:             1.0
+// Date:                08/12/2009
 // Author:              P R Truscott
 // Organisation:        QinetiQ Ltd, UK
@@ -57 +57 @@
 // 15 March 2004, P R Truscott, QinetiQ Ltd, UK
 // Beta release
+//
+// 08 December 2009, P R Truscott, QinetiQ Ltd, Ltd
+// ver 1.0
+// Variable fradius defined. See .cc file for more details.
 //
 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -119 +123 @@
     G4double               B;
     G4double               third;
+    G4double               fradius;
 };
 ////////////////////////////////////////////////////////////////////////////////

trunk/source/processes/hadronic/models/abrasion/src/G4WilsonAbrasionModel.cc

@@ -38 +38 @@
 // MODULE:              G4WilsonAbrasionModel.cc
 //
-// Version:             B.2
-// Date:                18/01/05
+// Version:             1.0
+// Date:                08/12/2009
 // Author:              P R Truscott
 // Organisation:        QinetiQ Ltd, UK
@@ -60 +60 @@
 // handler deleted to prevent memory leaks.  Also particle change of projectile
 // fragment previously not properly defined.
+//
+// 08 December 2009, P R Truscott, QinetiQ Ltd, Ltd
+// ver 1.0
+// There was originally a possibility of the minimum impact parameter AFTER
+// considering Couloumb repulsion, rm, being too large.  Now if:
+//     rm >= fradius * (rP + rT)
+// where fradius is currently 0.99, then it is assumed the primary track is
+// unchanged.  Additional conditions to escape from while-loop over impact
+// parameter: if the loop counter evtcnt reaches 1000, the primary track is
+// continued.
+// Additional clauses have been included in
+//    G4WilsonAbrasionModel::GetNucleonInducedExcitation
+// Previously it was possible to get sqrt of negative number as Wilson algorithm
+// not properly defined if either:
+//    rT > rP && rsq < rTsq - rPsq) or (rP > rT && rsq < rPsq - rTsq)
 //
 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -151 +166 @@
   B                = 10.0 * MeV;
   third            = 1.0 / 3.0;
+  fradius          = 0.99;
   conserveEnergy   = false;
   conserveMomentum = true;
@@ -197 +213 @@
   B                = 10.0 * MeV;
   third            = 1.0 / 3.0;
+  fradius          = 0.99;
   conserveEnergy   = false;
   conserveMomentum = true;
@@ -208 +225 @@
 // The destructor doesn't have to do a great deal!
 //
-  delete theExcitationHandler;
-  delete theExcitationHandlerx;
+  if (theExcitationHandler)  delete theExcitationHandler;
+  if (theExcitationHandlerx) delete theExcitationHandlerx;
+  if (useAblation)           delete theAblation;
+//  delete theExcitationHandler;
+//  delete theExcitationHandlerx;
 }
 ////////////////////////////////////////////////////////////////////////////////
@@ -304 +324 @@
     G4double rPT   = rP + rT;
     G4double rPTsq = rPT * rPT;
+//
+//
+// This is a "catch" to make sure we don't go into an infinite loop because the
+// energy is too low to overcome nuclear repulsion.  PRT 20091023.  If the
+// value of rm < fradius * rPT then we're unlikely to sample a small enough
+// impact parameter (energy of incident particle is too low).  Return primary
+// and don't complete nuclear interaction analysis. 
+//
+    if (rm >= fradius * rPT) {
+      theParticleChange.SetStatusChange(isAlive);
+      theParticleChange.SetEnergyChange(theTrack.GetKineticEnergy());
+      theParticleChange.SetMomentumChange(theTrack.Get4Momentum().vect().unit());
+      if (verboseLevel >= 2) {
+        G4cout <<"Particle energy too low to overcome repulsion." <<G4endl;
+        G4cout <<"Event rejected and original track maintained" <<G4endl;
+        G4cout <<"########################################"
+               <<"########################################"
+               <<G4endl;
+      }
+      return &theParticleChange;
+    }
+//
+//
+// Now sample impact parameter until the criterion is met that projectile
+// and target overlap, but repulsion is taken into consideration.
+//
+    G4int evtcnt   = 0;
     r              = 1.1 * rPT;
-    while (r > rPT)
+    while (r > rPT && ++evtcnt < 1000)
     {
       G4double bsq = rPTsq * G4UniformRand();
       r            = (rm + std::sqrt(rm*rm + 4.0*bsq)) / 2.0;
     }
+//
+//
+// We've tried to sample this 1000 times, but failed.  Assume nuclei do not
+// collide.
+//
+    if (evtcnt >= 1000) {
+      theParticleChange.SetStatusChange(isAlive);
+      theParticleChange.SetEnergyChange(theTrack.GetKineticEnergy());
+      theParticleChange.SetMomentumChange(theTrack.Get4Momentum().vect().unit());
+      if (verboseLevel >= 2) {
+        G4cout <<"Particle energy too low to overcome repulsion." <<G4endl;
+        G4cout <<"Event rejected and original track maintained" <<G4endl;
+        G4cout <<"########################################"
+               <<"########################################"
+               <<G4endl;
+      }
+      return &theParticleChange;
+    }
+
+
     rsq = r * r;
 //
@@ -764 +831 @@
   if (r > rT) Cl = 2.0*std::sqrt(rPsq + 2.0*r*rT - rsq - rTsq);
   else        Cl = 2.0*rP;
-  
-  G4double bP = (rPsq+rsq-rTsq)/2.0/r;
-  G4double Ct = 2.0*std::sqrt(rPsq - bP*bP);
+//
+//
+// The next lines have been changed to include a "catch" to make sure if the
+// projectile and target are too close, Ct is set to twice rP or twice rT.
+// Otherwise the standard Wilson algorithm should work fine.
+// PRT 20091023.
+//
+  G4double Ct = 0.0;
+  if      (rT > rP && rsq < rTsq - rPsq) Ct = 2.0 * rP;
+  else if (rP > rT && rsq < rPsq - rTsq) Ct = 2.0 * rT;
+  else {
+    G4double bP = (rPsq+rsq-rTsq)/2.0/r;
+    G4double x  = rPsq - bP*bP;
+    if (x < 0.0) {
+      G4cerr <<"########################################"
+             <<"########################################"
+             <<G4endl;
+      G4cerr <<"ERROR IN G4WilsonAbrasionModel::GetNucleonInducedExcitation"
+             <<G4endl;
+      G4cerr <<"rPsq - bP*bP < 0.0 and cannot be square-rooted" <<G4endl;
+      G4cerr <<"Set to zero instead" <<G4endl;
+      G4cerr <<"########################################"
+             <<"########################################"
+             <<G4endl;
+    }
+    Ct = 2.0*std::sqrt(x);
+  }
   
   G4double Ex = 13.0 * Cl / fermi;