Changeset 1315 for trunk/source/processes/hadronic/models/cascade/cascade/src/paraMaker.cc

Timestamp:

Jun 18, 2010, 11:42:07 AM (14 years ago)

Author:

garnier

Message:

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08

File:

• trunk/source/processes/hadronic/models/cascade/cascade/src/paraMaker.cc (modified) (3 diffs)

trunk/source/processes/hadronic/models/cascade/cascade/src/paraMaker.cc

@@ -23 +23 @@
 // * acceptance of all terms of the Geant4 Software license.          *
 // ********************************************************************
+// $Id: paraMaker.cc,v 1.17 2010/06/01 17:20:31 mkelsey Exp $
+// Geant4 tag: $Name: geant4-09-04-beta-cand-01 $
 //
+// 20100412  M. Kelsey -- Modify paraMaker[Truncated] to take buffer as argument
+// 20100517  M. Kelsey -- BUG FIX:  Must check for array boundary "if (Z>=70)"
+// 20100517  M. Kelsey -- Use G4CascadeInterpolator, which handles boundaries
+// 20100601  M. Kelsey -- Bug fix from Gunter Folger; resize(6,0.), not clear()
+
 #include "G4InuclSpecialFunctions.hh"
+#include "G4CascadeInterpolator.hh"
 
-std::pair<std::vector<G4double>, std::vector<G4double> > G4InuclSpecialFunctions::paraMaker(G4double Z) {
-  G4int verboseLevel = 1;
+void 
+G4InuclSpecialFunctions::paraMaker(G4double Z,
+           std::pair<std::vector<G4double>, std::vector<G4double> >& parms) {
+  G4int verboseLevel(0);
 
   if (verboseLevel > 3) {
@@ -36 +46 @@
   // coulumb barier, c.s. etc needed for evaporators
 
-  const G4double Z1[5] = {10.0, 20.0, 30.0, 50.0, 70.0};
-  const G4double AP[5] = {0.42, 0.58, 0.68, 0.77, 0.80};
-  const G4double CP[5] = {0.50, 0.28, 0.20, 0.15, 0.10};
-  const G4double AA[5] = {0.68, 0.82, 0.91, 0.97, 0.98};
-  const G4double CA[5] = {0.10, 0.10, 0.10, 0.08, 0.06};
-  std::vector<G4double> AK(6);
-  std::vector<G4double> CPA(6);
+  static const G4double Z1[5] = {10.0, 20.0, 30.0, 50.0, 70.0};
+  static const G4double AP[5] = {0.42, 0.58, 0.68, 0.77, 0.80};
+  static const G4double CP[5] = {0.50, 0.28, 0.20, 0.15, 0.10};
+  static const G4double AA[5] = {0.68, 0.82, 0.91, 0.97, 0.98};
+  static const G4double CA[5] = {0.10, 0.10, 0.10, 0.08, 0.06};
+
+  // Set up input buffer for results
+  std::vector<G4double>& AK = parms.first; 
+  AK.resize(6,0.);
+
+  std::vector<G4double>& CPA = parms.second;
+  CPA.resize(6,0.);
+
   AK[0] = 0.0;
   CPA[0] = 0.0;
-  G4double AK2 = 0.0;
-  G4double CP2 = 0.0;
-  G4double AK6 = 0.0; 
-  G4double CP6 = 0.0;
 
-  if (Z < 10.0) {
-    AK2=0.42;
-    CP2=0.5;
-    AK6=0.68;
-    CP6=0.1;
+  static G4CascadeInterpolator<5> interp(Z1, false);    // Do not extrapolate
+  AK[1]  = interp.interpolate(Z, AP);
+  AK[5]  = interp.interpolate(Z, AA);
+  CPA[1] = interp.interpolate(Z, CP);
+  CPA[5] = interp.interpolate(Z, CA);
+  
+  AK[2] = AK[1] + 0.06;
+  AK[3] = AK[1] + 0.12;
+  AK[4] = AK[5] - 0.06;
 
-  } else if(Z > 70.0) {
-    AK6=0.98;  // normal
-    CP6=0.06;
-    AK2=0.8;
-    CP2=0.1;
-    // AK6=1.1; // modified
-    // CP6=0.0;
+  CPA[2] = CPA[1] * 0.5;
+  CPA[3] = CPA[1] / 3.0;  
+  CPA[4] = 4.0 * CPA[5] / 3.0;
 
-  } else {
-
-    for (G4int i = 1; i < 5; i++) {
-
-      if (Z <= Z1[i]) {
-        G4double Z2 = 1.0 / (Z1[i] - Z1[i - 1]);
-        AK2 = ((AP[i] - AP[i - 1]) * Z + AP[i - 1] * Z1[i] - AP[i] * Z1[i - 1]) * Z2;
-        CP2 = ((CP[i] - CP[i - 1]) * Z + CP[i - 1] * Z1[i] - CP[i] * Z1[i - 1]) * Z2;
-        AK6 = ((AA[i] - AA[i - 1]) * Z + AA[i - 1] * Z1[i] - AA[i] * Z1[i - 1]) * Z2;
-        CP6 = ((CA[i] - CA[i - 1]) * Z + CA[i - 1] * Z1[i] - CA[i] * Z1[i - 1]) * Z2;
-
-        break;
-      };
-    };
-  };
-
-  AK[1] = AK2;
-  AK[5] = AK6;
-  CPA[1] = CP2;
-  CPA[5] = CP6;
-  AK[2] = AK2 + 0.06;
-  CPA[2] = CP2 * 0.5;
-  AK[3] = AK2 + 0.12;
-  CPA[3] = CP2 / 3.0;  
-  AK[4] = AK6 - 0.06;
-  CPA[4] = 4.0 * CP6 / 3.0;
-
-  return std::pair<std::vector<G4double>, std::vector<G4double> >(AK, CPA);
+  return;       // Buffer filled
 }
 
-std::pair<G4double, G4double> G4InuclSpecialFunctions::paraMakerTruncated(G4double Z) {
-  G4int verboseLevel = 1;
+void 
+G4InuclSpecialFunctions::paraMakerTruncated(G4double Z,
+                                    std::pair<G4double,G4double>& parms) {
+  G4int verboseLevel(0);
 
   if (verboseLevel > 3) {
@@ -102 +89 @@
 
   // truncated version of the previous one
-  const G4double Z1[5] = {10.0, 20.0, 30.0, 50.0, 70.0};
-  const G4double AP[5] = {0.42, 0.58, 0.68, 0.77, 0.8};
-  const G4double CP[5] = {0.5, 0.28, 0.2, 0.15, 0.1};
-  G4double AK2=0.;
-  G4double CP2=0.;
+  static const G4double Z1[5] = {10.0, 20.0, 30.0, 50.0, 70.0};
+  static const G4double AP[5] = {0.42, 0.58, 0.68, 0.77, 0.8};
+  static const G4double CP[5] = {0.5, 0.28, 0.2, 0.15, 0.1};
 
-  if (Z < 10.0) {
-    AK2=0.42;
-    CP2=0.5;
+  // Set up buffers for output
+  G4double& AK2=parms.first;
+  G4double& CP2=parms.second;
 
-  } else if (Z > 70.0) {
-    AK2=0.8;
-    CP2=0.1;
+  static G4CascadeInterpolator<5> interp(Z1, false);            // Do not extrapolate
+  AK2 = interp.interpolate(Z, AP);
+  CP2 = interp.interpolate(Z, CP);
 
-  } else {
-
-    for (G4int i = 1; i < 5; i++) {
-
-      if (Z < Z1[i]) {
-        G4double Z2 = 1.0 / (Z1[i] - Z1[i - 1]);
-        AK2 = ((AP[i] - AP[i - 1]) * Z + AP[i - 1] * Z1[i] - AP[i] * Z1[i - 1]) * Z2;
-        CP2 = ((CP[i] - CP[i - 1]) * Z + CP[i - 1] * Z1[i] - CP[i] * Z1[i - 1]) * Z2;
-
-        break;
-      };
-    };
-  };
-
-  return std::pair<G4double, G4double>(AK2, CP2);
+  return;       // Buffer filled
 }