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

Timestamp:

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

Author:

garnier

Message:

update processes

File:

• trunk/source/processes/hadronic/models/coherent_elastic/src/G4HadronElastic.cc (modified) (13 diffs)

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