GeomConverter/src/org/lcsim/material
diff -u -r1.5 -r1.6
--- MaterialCalculator.java 22 Jul 2005 00:11:32 -0000 1.5
+++ MaterialCalculator.java 11 Jan 2006 00:18:44 -0000 1.6
@@ -1,20 +1,17 @@
-/*
- * CoulombCorrectionCalculator.java
- *
- * Created on July 1, 2005, 11:18 AM
- */
-
package org.lcsim.material;
-import static java.lang.Math.log;
+import static java.lang.Math.log10;
import static java.lang.Math.pow;
+import static java.lang.Math.PI;
+import static java.lang.Math.sqrt;
+import static java.lang.Math.log;
/**
*
* @author jeremym
*/
abstract public class MaterialCalculator
-{
+{
/* Fine structure constant from PDG July 2004 */
public static final double FINE_STRUCTURE_CONSTANT = 1 / 137.03599911;
@@ -99,8 +96,155 @@
return NIL;
}
+
+ /* electron mass x c^2 in MeV */
+ public static final double M_e = 0.510998918;
+
+ /* Avogadro's number */
+ public static final double N_A = 6.0221415e23;
+
+ /* classical electron radius in centimeters */
+ public static final double r_e = 2.817940325e-13;
+
+ //public static final double clight = 2.99792458e08; // m/s
+
+ // FIXME: add path length as in Lelaps
+ // FIXME: add delta term
+ // FIXME: debug and check output (Matprop)
+ // FIXME: take out debug prints when completed
+ public static double computeBetheBlock(Material material, double[] particleMomentum, double particleMass, double particleCharge)
+ {
+ assert( particleMomentum.length == 3);
+
+ double zeff = material.getZeff();
+ double aeff = material.getAeff();
+
+ /* K matches PDG, pg. 238 --> 0.307075 in MeV g-1 cm2 */
+ double K = ( (4*PI) * N_A * (r_e * r_e) * M_e ) / aeff;
+ double ZoverA = zeff / aeff;
+
+ double z2 = particleCharge * particleCharge;
+
+ double mag2 =
+ particleMomentum[0] * particleMomentum[0] +
+ particleMomentum[1] * particleMomentum[1] +
+ particleMomentum[2] * particleMomentum[2];
+
+ double beta2 =
+ mag2 / (particleMass * particleMass + mag2);
+
+ System.out.println("beta2 = " + beta2);
+
+ double coeff1 = K * z2 * ( ZoverA ) * ( 1 / beta2 );
+ System.out.println("coeff1 = " + coeff1);
+ // end first coefficient calc
+
+ double gamma = 1 / ( sqrt(1 + beta2) );
+ double gamma2 = gamma * gamma;
+
+ System.out.println("gamma2 = " + gamma2);
+
+ /* compute T_max */
+ double numer_T_max = 2 * M_e * beta2 * gamma2;
+
+ System.out.println("numer_T_max = " + numer_T_max);
+ double denom_T_max =
+ 1 + ( 2 * gamma * M_e / particleMass ) + pow( ( M_e / particleMass), 2);
+ System.out.println("denom_T_max = " + denom_T_max);
+
+ double T_max = numer_T_max / denom_T_max;
+ System.out.println("T_max = " + T_max);
+ /* end compute T_max */
+
+ /* compute I using lelaps/CEPack/cematerial.cc*/
+ double I = 0.0;
+ if ( zeff > 12 )
+ {
+ I = zeff * ( 9.76 + 58.8 * pow(zeff, -1.19));
+ }
+ else
+ {
+ if (zeff == 1)
+ {
+ I = 18.7;
+ }
+ else
+ {
+ I = 13.0 * zeff;
+ }
+ }
+ I *= 1e-6; // convert I to MeV
+ double I2 = I * I;
+
+ System.out.println("I = " + I);
+ System.out.println("I2 = " + I2);
+
+ // plasma E
+ double eta = 1.0;
+ double rho_STP = material.getDensity();
+ double rho = material.getDensity();
+ if ( material.getState() == org.lcsim.material.MaterialState.GAS )
+ {
+ eta = material.getPressure() * (273.15 / material.getTemperature() );
+ rho_STP = rho / eta;
+ }
+
+ double plasmaE = 28.816 * sqrt(rho_STP * ZoverA);
+
+ System.out.println("plasmaE = " + plasmaE);
+ plasmaE *= 1e-6;
+ //
+
+ double delta = log(plasmaE / I) + log((sqrt(beta2) * gamma)) - 1/2;
+ System.out.println("delta = " + delta);
+
+ double coeff2 = 0.5 * ( log ( ( 2 * M_e * beta2 * gamma2 * T_max ) / I2 ));
+
+ double _num = 2 * M_e * beta2 * gamma2 * T_max;
+ System.out.println("2me c2*beta2 ... = " + _num);
+ _num = _num / I2;
+ System.out.println("numerator / i2 = " + _num);
+ System.out.println("coeff2 = " + coeff2);
+ coeff2 -= beta2;
+// coeff2 -= ( delta / 2 );
+ System.out.println("coeff2 (final) = " + coeff2);
+
+ double result = coeff1 * coeff2;
+ System.out.println("final dEdX = " + result);
+ System.out.println("----");
+ return result;
+ }
}
-
+/*
+ beta2 = 0.5
+coeff1 = 0.6141498919582697
+gamma2 = 0.6666666666666669
+numer_T_max = 0.3406659453333334
+denom_T_max = 2.0955776330067057
+T_max = 0.16256422094205625
+I = 1.8699999999999997E-5
+I2 = 3.496899999999999E-10
+plasmaE = 28.816
+delta = -0.11689887976099544
+coeff2 = 9.440219502206741
+coeff2 (final) = 8.998668942087239
+final results = 5.526531558551115
+result = 5.526531558551115
+ */
+
+/*
+
+excitation potential from lelaps/CEPack/ceelement.c
+
+ if (e.Z > 12) {
+ I = Z * (9.76 + 58.8 * pow(Z, -1.19));
+ }
+ else {
+ if (e.Z == 1) I = 18.7; // ********should be 19.0
+ else I = 13.0 * e.Z;
+ }
+ */
+
/*
public static double computeCoulombCorrection(double Z)
{
GeomConverter/test/org/lcsim/material
diff -N BetheBlockTest.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ BetheBlockTest.java 11 Jan 2006 00:18:44 -0000 1.1
@@ -0,0 +1,49 @@
+/*
+ * MaterialManagerTest.java
+ *
+ * Created on July 1, 2005, 2:25 PM
+ */
+
+package org.lcsim.material;
+
+import junit.framework.Test;
+import junit.framework.TestCase;
+import junit.framework.TestSuite;
+
+import org.lcsim.material.Material;
+import org.lcsim.material.MaterialElement;
+
+/**
+ *
+ * @author jeremym
+ */
+public class BetheBlockTest extends TestCase
+{
+ /** Creates a new instance of MaterialManagerTest */
+ public BetheBlockTest()
+ {}
+
+ public BetheBlockTest(String testName)
+ {
+ super(testName);
+ }
+
+ public static TestSuite suite()
+ {
+ return new TestSuite(BetheBlockTest.class);
+ }
+
+ public void test1()
+ {
+ MaterialElement element = new MaterialElement("DummyElement", 1.0, 1.0);
+ Material material = new Material("DummyMaterial", 1, 1.0, org.lcsim.material.MaterialState.SOLID, 1.0, 1.0);
+ material.addElement(element, 1.0);
+
+ double pT[] = {1.0, 0, 0};
+ double mass = 1.0;
+ //double mass = 1.0;
+
+ double dEdx = MaterialCalculator.computeBetheBlock(material, pT, mass, 1.0);
+ System.out.println("result = " + dEdx);
+ }
+}
\ No newline at end of file
