LISTSERV 16.5 - LCD-CVS Archives

+package org.lcsim.detector; + +import static org.lcsim.units.clhep.SystemOfUnits.m; + +import org.lcsim.detector.material.IMaterial; +import org.lcsim.detector.material.MaterialElement; +import org.lcsim.detector.solids.Box; + +public class DetectorTestUtil +{ + private static IMaterial dummymat = new MaterialElement("dummymat",1,1,1.0); + + public static IPhysicalVolume createTestSetup() + { + IPhysicalVolume world = createWorld(); + createBox(world); + return world; + } + + public static final void createBox(IPhysicalVolume mom) + { + Box box = new Box("test_box",1.0*m,1.0*m,1.0*m); + LogicalVolume lvTest = new LogicalVolume("lvTest",box,dummymat); + new PhysicalVolume( + new CoordinateTransformation3D(), + "test", + lvTest, + mom.getLogicalVolume(), + 1); + } + + public static final IPhysicalVolume createWorld() + { + Box boxWorld = new Box( + "world_box", + 10.0*m, + 10.0*m, + 10.0*m); + + LogicalVolume lvWorld = + new LogicalVolume( + "world", + boxWorld, + dummymat); + + IPhysicalVolume pvTop = + new PhysicalVolume( + null, + "world", + lvWorld, + null, + 0); + + return pvTop; + } +}

+package org.lcsim.detector; + +public interface IPhysicalVolumeNavigator +{ + public void setTopPhysicalVolume(IPhysicalVolume physvol); + public IPhysicalVolumePath getPath(IPhysicalVolume physvol); + public IPhysicalVolumePath getPath(String path); + //public IPhysicalVolumePath getPath(String[] path); + //public IPhysicalVolumePath getPath(int[] copyNums); + //public ICoordinateTransformation3D getTransformation(IPhysicalVolume physvol); + //public ICoordinateTransformation3D getInverseTransformation(IPhysicalVolume physvol); + //public IPhysicalVolume getPhysicalVolume(Hep3Vector globalPoint); +}

+implements IPhysicalVolumeContainer +{ + public boolean uniqPhysVol=false; + public boolean uniqName=false; + public boolean uniqCopyNum=false; + + public PhysicalVolumeContainer() + {} + + public PhysicalVolumeContainer(boolean uniqPhysVol, boolean uniqName, boolean uniqCopyNum) + { + this.uniqPhysVol = uniqPhysVol; + this.uniqName = uniqName; + this.uniqCopyNum= uniqCopyNum; + } + + public boolean add(IPhysicalVolume physvol) + { + boolean added=false; + if (checkAdd(physvol)) + { + added = super.add(physvol); + } + else { + System.err.println("Unable to add the physvol <"+physvol.getName()+"> to the container."); + } + return added; + } + + private boolean checkAdd(IPhysicalVolume physvol) + { + boolean addit=true; + + if ( uniqPhysVol && contains(physvol)) + { + addit=false; + } + else if (uniqName && findByName(physvol.getName()).size() != 0) + { + addit=false; + } + else if (uniqCopyNum && findByCopyNum(physvol.getCopyNumber()).size() != 0) + { + addit=false; + } + return addit; + } + + public PhysicalVolumeContainer findByName(String name) + { + PhysicalVolumeContainer physvols = new PhysicalVolumeContainer(); + for (IPhysicalVolume physvol : this) + { + if (physvol.getName().equals(name)) + { + physvols.add(physvol); + } + } + return physvols; + } + + public PhysicalVolumeContainer findByCopyNum(int copyNum) + { + PhysicalVolumeContainer physvols = new PhysicalVolumeContainer(); + for (IPhysicalVolume physvol : this) + { + if (physvol.getCopyNumber() == copyNum) + { + physvols.add(physvol); + } + } + return physvols; + } +}

-package org.lcsim.detector; - -public interface IMaterial -{ - public String getName(); - public double getDensity(); - public double Z(); - public double A(); - public double getNuclearInteractionLength(); - public double getRadiationLength(); - public double getIonizationPotential(); - public State getState(); - enum State - { - Gas, - Liquid, - Solid, - Unknown - } -}

-package org.lcsim.detector; - -public class SolidStore -extends ObjectStore<ISolid> -implements ISolidStore -{ - private static ISolidStore store; - public static ISolidStore getInstance() - { - if ( store == null ) - { - store = new SolidStore(); - } - return store; - } -}

+package org.lcsim.detector.solids; + +import hep.physics.vec.Hep3Vector; + +import org.lcsim.detector.INamed; + +public interface ISolid +extends INamed +{ + public boolean isInside(Hep3Vector position); + public double getCubicVolume(); +}

+package org.lcsim.detector.solids; + +import org.lcsim.detector.ObjectStore; + +public class SolidStore +extends ObjectStore<ISolid> +implements ISolidStore +{ + private static ISolidStore store; + public static ISolidStore getInstance() + { + if ( store == null ) + { + store = new SolidStore(); + } + return store; + } +}

+ + public double computeRadiationLength(Material m) + { + double z = m.getZeff(); + double a = m.getAeff(); + double ionpot = + (4 * fine_structure_const * Avogadro * z * (z + 1) * + classic_electr_radius * classic_electr_radius * + log(183.0 * pow(z,-.33333334))) / a; + ionpot = 1.0/ionpot; + return ionpot; + } +

GeomConverter/src/org/lcsim/detector/material

BetheBlochCalculator.java added at 1.1

diff -N BetheBlochCalculator.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ BetheBlochCalculator.java	3 Mar 2007 01:19:39 -0000	1.1
@@ -0,0 +1,274 @@

+package org.lcsim.detector.material;
+
+import static java.lang.Math.PI;
+import static java.lang.Math.log;
+import static java.lang.Math.log10;
+import static java.lang.Math.pow;
+import static java.lang.Math.sqrt;
+import static org.lcsim.units.clhep.PhysicalConstants.Avogadro;
+import static org.lcsim.units.clhep.PhysicalConstants.classic_electr_radius;
+import static org.lcsim.units.clhep.PhysicalConstants.electron_mass_c2;
+import hep.physics.particle.Particle;
+import hep.physics.vec.Hep3Vector;
+
+import org.lcsim.detector.material.IMaterial.State;
+
+/**
+ *
+ * A calculator for dEdX using the Bethe-Bloch formula.
+ *
+ * @author Jeremy McCormick <[log in to unmask]>
+ * @author Caroline Milstene <[log in to unmask]>
+ */
+public class BetheBlochCalculator
+{	
+	/**
+	 * 
+	 * Calculate Bethe-Bloch using IMaterial and Particle.
+	 * 
+	 * @param material
+	 * @param particle
+	 * @param distance
+	 * @return
+	 */
+	public static double computeBetheBloch(
+			IMaterial material,
+			Particle particle,
+			double distance)
+	{
+		return computeBetheBloch(
+				material,
+				particle.getMomentum(),
+				particle.getMass(),
+				particle.getCharge(), 
+				distance);
+	}
+	
+	/**
+	 * Calculate Bethe-Bloch using an IMaterial.
+	 * 
+	 * @param m
+	 * @param p
+	 * @param mass
+	 * @param charge
+	 * @param distance
+	 * @return
+	 */
+	public static double computeBetheBloch(
+			IMaterial material,
+			Hep3Vector p,
+			double mass,
+			double charge,
+			double distance)
+	{
+		return computeBetheBloch(
+				material.getZ(),
+				material.getA(),
+				material.getDensity(),
+				material.getState(),
+				material.getPressure(),
+				material.getTemperature(),
+				p,
+				mass,
+				charge,
+				distance);
+	}
+	
+	/**
+	 * Calculate Bethe Bloch with numerical parameters.
+	 * 
+	 * The pressure and temperature arguments are only used if the
+	 * state is a gas.
+	 * 
+	 * @param Z
+	 * @param A
+	 * @param density
+	 * @param state
+	 * @param pressure
+	 * @param temperature
+	 * @param p
+	 * @param mass
+	 * @param charge
+	 * @param distance
+	 * @return
+	 */
+    public static double computeBetheBloch(
+    		double Z,
+    		double A,
+    		double density,
+    		State state,
+    		double pressure,
+    		double temperature,
+    		Hep3Vector p,
+            double mass, 
+            double charge, 
+            double distance)
+    {        
+        double zeff = Z;
+        double aeff = A;
+        double ZoverA = zeff / aeff;
+        
+        // K matches PDG, pg. 238 --> 0.307075 in MeV g-1 cm2 */
+        double K = ((4 * PI) * Avogadro * (classic_electr_radius * classic_electr_radius) * electron_mass_c2);
+        double z2 = charge * charge;
+        
+        // Convert p(GeV) to p(MeV).       
+        double[] pmev = new double[3];
+        for (int i = 0; i < 3; i++)
+        {
+            pmev[i] = p.v()[i] * 1e+3;
+        }
+        
+        double mag2 = pmev[0] * pmev[0] + pmev[1] * pmev[1] + pmev[2] * pmev[2];
+        
+        double beta2 = mag2 / (mass * mass + mag2);
+        //double beta = sqrt(beta2);
+        
+        double coeff1 = K * z2 * (ZoverA) * (1 / beta2);        
+        
+        double gamma = 1 / (sqrt(1 - beta2));
+        double gamma2 = gamma * gamma;
+                
+        // Compute T_max.
+        double numer_T_max = 2 * electron_mass_c2 * beta2 * gamma2;
+        double denom_T_max = 1 + (2 * gamma * electron_mass_c2 / mass) + pow((electron_mass_c2 / mass), 2);
+        double T_max = numer_T_max / denom_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;
+        
+        // Compute plasma E.
+        double eta = 1.0;
+        double rho_STP = density;
+        double rho = density;
+        
+        if (state == State.Gas)
+        {
+            eta = pressure * (273.15 / temperature);
+            rho_STP = rho / eta;
+        }
+        
+        double plasmaE = 28.816 * sqrt(rho_STP * ZoverA);
+        plasmaE *= 1e-6;
+        
+        // Cbar
+        double Cbar = 2.0 * log( I / plasmaE ) + 1.0;
+        
+        // Xa
+        double Xa = Cbar / 4.6052;
+        
+        double X1;
+        if (state == State.Gas)
+        {
+            if (Cbar < 12.25)
+            {
+                X1 = 4.0;
+            }
+            else
+            {
+                X1 = 5.0;
+            }
+        }
+        else
+        {
+            if ( I < 100.0)
+            {
+                X1 = 2.0;
+            }
+            else
+            {
+                X1 = 3.0;
+            }
+        }
+        
+        double m = 3.0;
+        
+        double X0 = 0;
+        if (state == State.Gas)
+        {
+            if (Cbar < 10.0)        X0 = 1.6;
+            else if (Cbar < 10.5)   X0 = 1.7;
+            else if (Cbar < 11.0)   X0 = 1.8;
+            else if (Cbar < 11.5)   X0 = 1.9;
+            else if (Cbar < 13.804) X0 = 2.0;
+            else                    X0 = 0.326 * Cbar - 2.5;
+        }
+        else
+        {
+            if (I < 100.0)
+            {
+                if (Cbar < 3.681) X0 = 0.2;
+                else              X0 = 0.326 * Cbar - 1.0;
+                X1 = 2.0;
+            }
+            else
+            {
+                if (Cbar < 5.215) X0 = 0.2;
+                else              X0 = 0.326 * Cbar - 1.5;
+                X1 = 3.0;
+            }
+        }
+        
+        double a = 4.6052 * (Xa - X0) / pow(X1 - X0, m);
+        
+        //double ASP = 0.1536 * ZoverA;        
+        //double BSP = log(511.0e9 / I2);
+        
+        if (state == State.Gas)
+        {
+            double eta_corr_1 = 0.5 * log10(eta);
+            double eta_corr_2 = 4.6052 * eta_corr_1;
+            
+            Cbar -= eta_corr_2;
+            X1   -= eta_corr_1;
+            X0   -= eta_corr_1;
+        }
+        
+        //double delta_estimate = log(plasmaE / I) + log((sqrt(beta2) * gamma)) - 0.5;
+        
+        double delta = 0;
+        double X = log10(sqrt(mag2) / (mass));
+        
+        if ( X0 < X && X < X1 )
+        {
+            delta = 4.6052 * X - Cbar + a * (X1 - X);
+        }
+        else if ( X > X1 )
+        {
+            delta = 4.6052 * X - Cbar;
+        }
+        else if ( X < X0 )
+        {
+            delta = 0;
+        }
+               
+        double coeff2 = 0.5 * (log((2 * electron_mass_c2 * beta2 * gamma2 * T_max) / I2));
+        
+        coeff2 -= beta2;             
+        coeff2 -= delta;     
+        
+        double result = coeff1 * coeff2;
+        
+        result = result * density;
+        result = result * distance;        
+        
+        return result;
+    }   
+}

\ No newline at end of file

+package org.lcsim.detector.material; + +public interface IMaterial +{ + public String getName(); + public double getDensity(); + public double getZ(); + public double getA(); + public double getMolecularWeight(); + public double getEffectiveNumberOfNucleons(); + public double getNuclearInteractionLength(); + public double getRadiationLength(); + public State getState(); + public double getMeltingPoint(); + public double getBoilingPoint(); + public double getTemperature(); + public double getPressure(); + enum State + { + Gas, + Liquid, + Solid, + Unknown + } +}

GeomConverter/src/org/lcsim/detector/material

MaterialElement.java added at 1.1

diff -N MaterialElement.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ MaterialElement.java	3 Mar 2007 01:19:39 -0000	1.1
@@ -0,0 +1,259 @@

+package org.lcsim.detector.material;
+
+import static java.lang.Math.log;
+import static java.lang.Math.pow;
+import static org.lcsim.units.clhep.PhysicalConstants.Avogadro;
+import static org.lcsim.units.clhep.PhysicalConstants.fine_structure_const;
+
+/**
+ * A refactored version of a material element class
+ * for the detector description package.
+ * 
+ * Uses CLHEP values of Avogadro, fine structure constant, and electron radius.
+ * 
+ * Computes radiation length using the Tsai formula ported from Lelaps.
+ * 
+ * Computes nuclear interaction length with a simple fit ported from Lelaps.
+ * 
+ * @author Jeremy McCormick <[log in to unmask]>
+ */
+public class MaterialElement 
+implements IMaterial
+{
+	// Defaults for non-mandatory parameters.
+	static private final double defaultTemperature=0.0;
+	static private final State defaultState=State.Unknown;
+	static private final double defaultBoilingPoint=0.0;
+	static private final double defaultMeltingPoint=0.0;
+	static private final double defaultPressure=0.0;
+	static private final double defaultIonizationPotential=0.0;
+	
+	// From constructor.
+	private String name;
+	private double Z;
+	private double A;
+	private double density;
+	
+	// Set by the detailed constructor.  Otherwise, use defaults.
+	private double ionizationPotential=defaultIonizationPotential;
+	private double temperature=defaultTemperature;
+	private double pressure=defaultPressure;
+	private double meltingPoint=defaultMeltingPoint;
+	private double boilingPoint=defaultBoilingPoint;
+	private State state=defaultState;
+	
+	// Derived quantities.
+	private double N;
+	private double molecularWeight;
+	private double nuclearInteractionLength;
+	private double radiationLength;
+	
+	public MaterialElement(
+			String name,
+			int Z,
+			int A,
+			double density)
+	{
+		// Basic parameters.
+		this.name = name;
+		this.Z = Z;
+		this.A = A;
+		this.density = density;		
+		
+		// Compute derived.
+		computeDerivedQuantities();
+	}	
+	
+	public MaterialElement(
+			String name, 
+			int Z, 
+			int A, 
+			double density, 
+			double ionizationPotential, 
+			State state,
+			double temperature,
+			double pressure,
+			double meltingPoint,
+			double boilingPoint)
+	{
+		// Basic parameters.
+		this.name = name;
+		this.Z = Z;
+		this.A = A;
+		this.density = density;
+		
+		// Detailed parameters.
+        this.density = density;
+        this.ionizationPotential = ionizationPotential;
+        this.state = state;
+		this.temperature = temperature;
+		this.pressure = pressure;
+		this.meltingPoint = meltingPoint;
+		this.boilingPoint = boilingPoint;
+                
+		// Compute derived.
+        computeDerivedQuantities();
+	}
+	
+	private void computeDerivedQuantities()
+	{		
+		computeEffectiveNumberOfNucleons();
+		computeMolecularWeigth();
+		computeRadiationLength();
+		computeNuclearInteractionLength();
+	}
+	
+	private void computeEffectiveNumberOfNucleons()
+	{
+		N = A / Avogadro;	
+	}
+	
+	// FIXME: Wrong for isotopes.
+	private void computeMolecularWeigth()
+	{
+		double natoms=1.0;
+		molecularWeight = A * natoms;
+	}
+	
+	private void computeRadiationLength()
+	{
+		radiationLength = computeRadiationLengthTsai(A,Z);
+	}
+	
+	private void computeNuclearInteractionLength()
+	{
+		nuclearInteractionLength = computeNuclearInteractionLength(A,Z);		
+	}
+
+	public double getA() 
+	{
+		return A;
+	}
+
+	public double getDensity() {
+		return density;
+	}
+
+	public String getName() {
+		return name;
+	}
+
+	public double getNuclearInteractionLength() 
+	{
+		return nuclearInteractionLength;
+	}
+
+	public double getRadiationLength() 
+	{
+		return radiationLength;
+	}
+
+	public State getState() 
+	{
+		return state;
+	}
+
+	public double getZ() 
+	{
+		return Z;
+	}
+
+	public double getBoilingPoint() 
+	{
+		return boilingPoint;
+	}
+
+	public double getMeltingPoint() 
+	{
+		return meltingPoint;
+	}
+
+	public double getPressure() 
+	{
+		return pressure;
+	}
+
+	public double getTemperature() 
+	{
+		return temperature;
+	}
+	
+	public double getMolecularWeight()
+	{
+		return molecularWeight;
+	}
+	
+	public double getIonizationPotential()
+	{
+		return ionizationPotential;
+	}
+	
+	public double getEffectiveNumberOfNucleons()
+	{
+		return N;
+	}
+	
+	public static double computeRadiationLengthTsai(double A, double Z)
+    {
+        double azsq = fine_structure_const * Z;
+        azsq *= azsq;
+        double f = azsq
+                * (1.0 / (1.0 + azsq) + 0.20206 - 0.0369 * azsq + 0.0083 * azsq
+                * azsq - 0.002 * azsq * azsq * azsq);
+        
+        double Lrad, LradP;
+        if (Z == 1)
+        {
+            Lrad = 5.31;
+            LradP = 6.144;
+        }
+        else if (Z == 2)
+        {
+            Lrad = 4.79;
+            LradP = 5.621;
+        }
+        else if (Z == 3)
+        {
+            Lrad = 4.74;
+            LradP = 5.805;
+        }
+        else if (Z == 4)
+        {
+            Lrad = 4.71;
+            LradP = 5.924;
+        }
+        else
+        {
+            Lrad = log(184.15 / pow(Z, 0.333333333));
+            LradP = log(1194.0 / pow(Z, 0.666666667));
+        }
+        double rlen = 716.408 * A / ((Z * Z * (Lrad - f) + Z * LradP));
+        return rlen;
+    }
+	
+    public static double computeNuclearInteractionLength(double A, double Z)
+    {
+        double NIL = 0;
+        if (Z == 1)
+        {
+            if (A < 1.5)
+            {
+                NIL = 50.8;
+            }
+            else
+            {
+                NIL = 54.7;
+            }
+        }
+        else if (Z == 2)
+        {
+            NIL = 65.2;
+        }
+        else
+        {
+            NIL = 40.8 * pow(A, 0.289);
+        }
+        
+        return NIL;
+    }	
+}

\ No newline at end of file

+package org.lcsim.detector.material; + +import java.util.ArrayList; + +public class MaterialStore +extends ArrayList<IMaterial> +implements IMaterialStore +{ + private static MaterialStore materialStore = null; + public static final IMaterialStore getInstance() + { + if ( materialStore == null ) + { + materialStore = new MaterialStore(); + } + return materialStore; + } +}

Commit in `GeomConverter/src/org/lcsim` on MAIN
`detector/DetectorTestUtil.java`	+56		added 1.1
`/IPhysicalVolumeNavigator.java`	+13		added 1.1
`/ILogicalVolume.java`	+3		1.3 -> 1.4
`/IPhysicalVolume.java`	+1	-1	1.2 -> 1.3
`/IPhysicalVolumeContainer.java`	+4	-1	1.2 -> 1.3
`/LogicalVolume.java`	+21	-4	1.1 -> 1.2
`/PhysicalVolume.java`	+18	-2	1.1 -> 1.2
`/PhysicalVolumeContainer.java`	+74	-1	1.1 -> 1.2
`/IMaterial.java`		-20	1.1 removed
`/ISolid.java`		-10	1.2 removed
`/ISolidStore.java`		-5	1.1 removed
`/SolidStore.java`		-16	1.1 removed
`detector/solids/ISolid.java`	+12		added 1.1
`/ISolidStore.java`	+7		added 1.1
`/SolidStore.java`	+18		added 1.1
`/AbstractSolid.java`		-2	1.1 -> 1.2
`/Tube.java`	+1	-2	1.1 -> 1.2
`material/MaterialCalculator.java`	+15		1.15 -> 1.16
`detector/material/BetheBlochCalculator.java`	+274		added 1.1
`/IMaterial.java`	+25		added 1.1
`/IMaterialStore.java`	+7		added 1.1
`/MaterialElement.java`	+259		added 1.1
`/MaterialStore.java`	+18		added 1.1
	+826	-64