lcsim/src/org/lcsim/recon/tracking/seedtracker
diff -u -r1.2 -r1.3
--- MaterialManager.java 13 Oct 2008 01:05:27 -0000 1.2
+++ MaterialManager.java 7 Jan 2009 02:41:39 -0000 1.3
@@ -6,7 +6,6 @@
* To change this template, choose Tools | Template Manager
* and open the template in the editor.
*/
-
package org.lcsim.recon.tracking.seedtracker;
import hep.physics.vec.Hep3Vector;
@@ -39,241 +38,209 @@
* @version 1.0
*/
public class MaterialManager {
-
+
private static final boolean DEBUG = false; //enable debug output
private static final boolean TUBE_ONLY = true; //only use Tube elements for calculating volume.
-
private List<MaterialPolyconeSegment> _matpc = new ArrayList<MaterialPolyconeSegment>();
private List<MaterialCylinder> _matcyl = new ArrayList<MaterialCylinder>();
private List<MaterialDisk> _matdsk = new ArrayList<MaterialDisk>();
- private HashMap<ISolid,Double> solid_vol_map = new HashMap<ISolid,Double>(400);
+ private HashMap<ISolid, Double> solid_vol_map = new HashMap<ISolid, Double>(400);
private double _rmax;
private double _zmax = 1800.;
-
+
/** Creates a new instance of MaterialManager */
public MaterialManager() {
}
-
+
public void BuildModel(Detector det) {
-
+
// Build the model of tracker material
// Each volume defined in the compact.xml file is modelled as either
// a thin cylinder or disk with a thickness in radiation lengths
// that gives the correct total amount of material.
//
// First find the logical volume associated with the tracker
-
+
IPhysicalVolumeNavigator nav = det.getNavigator();
ILogicalVolume ltrkr = det.getTrackingVolume().getLogicalVolume();
// Loop over the volumes defined at the compact.xml level
for (IPhysicalVolume pvtree : ltrkr.getDaughters()) {
// Flatten the geometry tree to get all daughters with material
- //List<IPhysicalVolume> pvflat = Flatten(pvtree);
- List<UniquePV> pvflat = Flatten(pvtree,nav);
-
-// System.out.println(pvflat.size());
+ List<UniquePV> pvflat = Flatten(pvtree, nav);
+
// Calculate the total volume of material, skip this object if 0
VolumeGroupInfo vgi = performVolumeGroupCalculations(pvflat);
-
+
double vtot;
- if (TUBE_ONLY) vtot = vgi.vtot_tube_only;
- else vtot = vgi.vtot;
-
-
- if (pvtree.getLogicalVolume().getSolid() instanceof Polycone){
+ if (TUBE_ONLY) {
+ vtot = vgi.vtot_tube_only;
+ } else {
+ vtot = vgi.vtot;
+ }
+
+
+ if (pvtree.getLogicalVolume().getSolid() instanceof Polycone) {
handlePolycone(pvtree);
continue;
}
-
+
if (vtot > 0.) {
-
+
// Calculate the average radiation length for this volume
-
-
+
+
// Determine if this volume should be modeled as barrel or disk
- if (isCylinder(vgi.rmin,vgi.rmax,vgi.zmin,vgi.zmax)) {
+ if (isCylinder(vgi.rmin, vgi.rmax, vgi.zmin, vgi.zmax)) {
// Calculate the weighted radius of the elements
double zlen = vgi.zmax - vgi.zmin;
double thickness = vtot / (2. * Math.PI * vgi.weighted_r * zlen * vgi.X0);
-
- if (DEBUG){
+
+ if (DEBUG) {
System.out.println(pvtree.getName());
- System.out.println("x0: "+vgi.X0 + "| zmin: "+vgi.zmin + "| zmax: "+vgi.zmax + "| vtot: "+vtot + "| thickness: "+thickness +
- "| rmin: "+vgi.rmin + "| rmax: "+vgi.rmax);
+ System.out.println("x0: " + vgi.X0 + "| zmin: " + vgi.zmin + "| zmax: " + vgi.zmax + "| vtot: " + vtot + "| thickness: " + thickness +
+ "| rmin: " + vgi.rmin + "| rmax: " + vgi.rmax);
System.out.println();
}
-
+
_matcyl.add(new MaterialCylinder(pvtree, vgi.weighted_r, vgi.zmin, vgi.zmax, thickness));
} else {
-
- double thickness = vtot / (Math.PI * (vgi.rmax*vgi.rmax - vgi.rmin*vgi.rmin) * vgi.X0);
-
- if (DEBUG){
+
+ double thickness = vtot / (Math.PI * (vgi.rmax * vgi.rmax - vgi.rmin * vgi.rmin) * vgi.X0);
+
+ if (DEBUG) {
System.out.println(pvtree.getName());
- System.out.println("x0: "+vgi.X0 + "| zmin: "+vgi.zmin + "| zmax: "+vgi.zmax + "| vtot: "+vtot + "| thickness: "+thickness +
- "| rmin: "+vgi.rmin + "| rmax: "+vgi.rmax);
+ System.out.println("x0: " + vgi.X0 + "| zmin: " + vgi.zmin + "| zmax: " + vgi.zmax + "| vtot: " + vtot + "| thickness: " + thickness +
+ "| rmin: " + vgi.rmin + "| rmax: " + vgi.rmax);
System.out.println();
}
-
+
_matdsk.add(new MaterialDisk(pvtree, vgi.rmin, vgi.rmax, vgi.weighted_z, thickness));
}
}
}
-
+
solid_vol_map.clear();
-
+
// Find the envelope of the tracking volume
- // First loop over all the subdetector elements
- IDetectorElement de = det.getDetectorElement();
- for (IDetectorElement child : de.getChildren()) {
-
- // Get the identifier for this subdetector
- IIdentifierHelper helper = child.getIdentifierHelper();
- if (helper instanceof DetectorIdentifierHelper) {
- DetectorIdentifierHelper dehelper = (DetectorIdentifierHelper) helper;
- IIdentifier id = child.getIdentifier();
-
- // Check if this is the ECalBarrel
- if (dehelper.isEcal(id)) {
-
- // Check for geometry info to make sure we don't have an ECal endcap
- if (child.hasGeometryInfo()) {
-
- // Get the solid and make sure its a tube
- ISolid solid = child.getGeometry().getLogicalVolume().getSolid();
- if (solid instanceof Tube) {
- Tube tube = (Tube) solid;
-
- // Take the inner ECal radius was the outer tracker radius
- _rmax = tube.getInnerRadius();
- System.out.println("Ecal radius = "+_rmax);
- }
- } else {
- // Have the ECal Endcap volume - loop over the endcaps
- for (IDetectorElement grandchild : child.getChildren()) {
- if (dehelper.isEcalEndcapPositive(grandchild.getIdentifier())) {
- if (grandchild.hasGeometryInfo()) {
- ISolid solid = grandchild.getGeometry().getLogicalVolume().getSolid();
- if (solid instanceof Tube) {
- Tube tube = (Tube) solid;
- _zmax = grandchild.getGeometry().getPosition().z() - tube.getZHalfLength();
- System.out.println("ECal inner Z = "+_zmax);
- }
- }
- }
- }
- }
- }
- }
+
+ ISolid trkvol = det.getTrackingVolume().getLogicalVolume().getSolid();
+ if (trkvol instanceof Tube) {
+ Tube trktube = (Tube) trkvol;
+ _rmax = trktube.getOuterRadius();
+ _zmax = trktube.getZHalfLength();
+ System.out.println("Ecal radius = "+_rmax);
+ System.out.println("ECal inner Z = "+_zmax);
}
+
return;
}
-
+
public List<MaterialCylinder> getMaterialCylinders() {
return _matcyl;
}
-
+
public List<MaterialDisk> getMaterialDisks() {
return _matdsk;
}
-
+
public List<MaterialPolyconeSegment> getMaterialPolyconeSegments() {
return _matpc;
}
-
+
public double getRMax() {
return _rmax;
}
-
+
public double getZMax() {
return _zmax;
}
-
- private List<UniquePV> Flatten(IPhysicalVolume vol, IPhysicalVolumeNavigator nav){
-
+
+ private List<UniquePV> Flatten(IPhysicalVolume vol, IPhysicalVolumeNavigator nav) {
+
LinkedList<UniquePV> pvtree = new LinkedList<UniquePV>();
List<UniquePV> pvflat = new ArrayList<UniquePV>();
- pvtree.add(new UniquePV(vol,nav));
-
- while(pvtree.size()>0) {
-
+ pvtree.add(new UniquePV(vol, nav));
+
+ while (pvtree.size() > 0) {
+
UniquePV upv = pvtree.poll();
IPhysicalVolume pv = upv.getPV();
-
- if (pv.getLogicalVolume().getNumberOfDaughters()==0) {
+
+ if (pv.getLogicalVolume().getNumberOfDaughters() == 0) {
pvflat.add(upv);
- }
-
- else {
- for(IPhysicalVolume p : pv.getLogicalVolume().getDaughters()) {
+ } else {
+ for (IPhysicalVolume p : pv.getLogicalVolume().getDaughters()) {
pvtree.add(upv.createDaughterUniquePV(p));
}
+
}
}
-
+
return pvflat;
}
-
-
+
private boolean isCylinder(double rmin, double rmax, double zmin, double zmax) {
return (rmax - rmin) * Math.abs(zmax + zmin) < (zmax - zmin) * (rmax + rmin);
}
-
- // special handling for Polycone...
- private void handlePolycone(IPhysicalVolume pv){
+
+// special handling for Polycone...
+ private void handlePolycone(IPhysicalVolume pv) {
Polycone pc = (Polycone) pv.getLogicalVolume().getSolid();
IMaterial mat = pv.getLogicalVolume().getMaterial();
-
+
//Loop through each segment
- for (int i = 0; i < pc.getNumberOfZPlanes()-1; i++){
+ for (int i = 0; i <
+ pc.getNumberOfZPlanes() - 1; i++) {
ZPlane zp1 = pc.getZPlane(i);
- ZPlane zp2 = pc.getZPlane(i+1);
-
+ ZPlane zp2 = pc.getZPlane(i + 1);
+
double z1 = zp1.getZ();
double z2 = zp2.getZ();
double vol = Polycone.getSegmentVolume(zp1, zp2);
- double zlen = Math.abs(z2-z1);
- double ravg = 0.25*(zp1.getRMax() + zp1.getRMin() + zp2.getRMax() + zp2.getRMin());
- double ang = Math.atan2( 0.5 * (zp1.getRMax() + zp1.getRMin() - zp2.getRMax() - zp2.getRMin())
- , zlen );
- double X0 = 10*mat.getRadiationLength()/mat.getDensity();
+ double zlen = Math.abs(z2 - z1);
+ double ravg = 0.25 * (zp1.getRMax() + zp1.getRMin() + zp2.getRMax() + zp2.getRMin());
+ double ang = Math.atan2(0.5 * (zp1.getRMax() + zp1.getRMin() - zp2.getRMax() - zp2.getRMin()), zlen);
+ double X0 = 10 * mat.getRadiationLength() / mat.getDensity();
double thickness = Math.cos(ang) * vol / (2 * Math.PI * ravg * zlen * X0);
-
+
//This is a cylinder
- if(zp1.getRMax()==zp2.getRMax() && zp1.getRMin()==zp2.getRMin()) {
- _matcyl.add(new MaterialCylinder(pv, ravg, Math.min(z1,z2), Math.max(z1,z2), thickness));
- if(DEBUG){
- System.out.println("Cylindrical segment of "+pv.getName());
- System.out.println("zmin = "+ z1 + "| zmax = "+z2 + "| ravg = "+ravg + "| thickness = "+thickness);
+ if (zp1.getRMax() == zp2.getRMax() && zp1.getRMin() == zp2.getRMin()) {
+ _matcyl.add(new MaterialCylinder(pv, ravg, Math.min(z1, z2), Math.max(z1, z2), thickness));
+ if (DEBUG) {
+ System.out.println("Cylindrical segment of " + pv.getName());
+ System.out.println("zmin = " + z1 + "| zmax = " + z2 + "| ravg = " + ravg + "| thickness = " + thickness);
}
- }
-
- //Otherwise this is a non-cylindrical polycone segment
+
+ } //Otherwise this is a non-cylindrical polycone segment
else {
_matpc.add(new MaterialPolyconeSegment(pv, zp1, zp2, thickness, ang));
- if(DEBUG){
- System.out.println("Non-Cylindrical segment of "+pv.getName());
- System.out.println("ZPlane 1: "+zp1.toString() + "| ZPlane 2: "+zp2.toString() + "| thickness = "+thickness);
+ if (DEBUG) {
+ System.out.println("Non-Cylindrical segment of " + pv.getName());
+ System.out.println("ZPlane 1: " + zp1.toString() + "| ZPlane 2: " + zp2.toString() + "| thickness = " + thickness);
}
+
+
+
}
}
}
-
+
/**
* A "struct" holding geometry information about a single physical volume
*/
class VolumeInfo {
+
double rmax = 0.0;
double rmin = 1.e10;
double zmin = 1.e10;
double zmax = -1.e10;
}
-
-
+
/**
* A "struct" holding geometry information about lists of physical volumes
*/
- class VolumeGroupInfo{
+ class VolumeGroupInfo {
+
double rmax = 0.0;
double rmin = 1.e10;
double zmin = 1.e10;
@@ -284,73 +251,82 @@
double vtot_tube_only = 0.;
double vtot = 0.0;
}
-
-
- //This function performs all the calculations on lists of physical volumes
- private VolumeGroupInfo performVolumeGroupCalculations(List<UniquePV> volgroup) {
-
+
+//This function performs all the calculations on lists of physical volumes
+ private VolumeGroupInfo performVolumeGroupCalculations(
+ List<UniquePV> volgroup) {
+
VolumeGroupInfo vgi = new VolumeGroupInfo();
-
+
//If we have a top-level polycone, don't bother doing anything, because it'll be handled specially
- if (volgroup.size()==1 && volgroup.get(0).getSolid() instanceof Polycone){
+ if (volgroup.size() == 1 && volgroup.get(0).getSolid() instanceof Polycone) {
return vgi;
}
-
- //The normal case
+
+//The normal case
double totwgt = 0.0;
- if (DEBUG && volgroup.isEmpty()) System.out.println("Empty volume group...");
- for(UniquePV pv : volgroup) {
-
+ if (DEBUG && volgroup.isEmpty()) {
+ System.out.println("Empty volume group...");
+ }
+ for (UniquePV pv : volgroup) {
+
//increment total volume
- double vol = this.getVolumeOfSolid(pv.getSolid());
- if (pv.getSolid() instanceof Tube) vgi.vtot_tube_only += vol;
- vgi.vtot+=vol;
+ double vol = this.getVolumeOfSolid(pv.getSolid());
+ if (pv.getSolid() instanceof Tube) {
+ vgi.vtot_tube_only += vol;
+ }
+ vgi.vtot += vol;
//calculate weighted R / Z / Radiation Length
VolumeInfo vi = performVolumeCalculations(pv);
IMaterial mat = pv.getPV().getLogicalVolume().getMaterial();
double matX0 = 10.0 * mat.getRadiationLength() / mat.getDensity();
double wgt = vol / matX0;
double z0 = pv.getLtoGTransform().getTranslation().z();
- vgi.weighted_r+= 0.5 * (vi.rmin + vi.rmax) * wgt;
- vgi.weighted_z+= z0 * wgt;
- totwgt += wgt;
-
+ vgi.weighted_r += 0.5 * (vi.rmin + vi.rmax) * wgt;
+ vgi.weighted_z += z0 * wgt;
+ totwgt +=
+ wgt;
+
//grab (z/r)(mins/maxes)
- vgi.zmin = Math.min(vi.zmin,vgi.zmin);
- vgi.zmax = Math.max(vi.zmax,vgi.zmax);
- vgi.rmin = Math.min(vi.rmin,vgi.rmin);
- vgi.rmax = Math.max(vi.rmax,vgi.rmax);
-
+ vgi.zmin = Math.min(vi.zmin, vgi.zmin);
+ vgi.zmax = Math.max(vi.zmax, vgi.zmax);
+ vgi.rmin = Math.min(vi.rmin, vgi.rmin);
+ vgi.rmax = Math.max(vi.rmax, vgi.rmax);
+
}
-
- //finish weighted R/Z calculations + perform X0 calculation
+
+//finish weighted R/Z calculations + perform X0 calculation
if (totwgt > 0.) {
vgi.weighted_r /= totwgt;
vgi.weighted_z /= totwgt;
vgi.X0 = vgi.vtot / totwgt;
}
-
+
return vgi;
}
-
- private double getVolumeOfSolid(ISolid solid){
- if (solid_vol_map.containsKey(solid)){
+
+ private double getVolumeOfSolid(ISolid solid) {
+ if (solid_vol_map.containsKey(solid)) {
return solid_vol_map.get(solid).doubleValue();
- }
-
- else {
+ } else {
double vol;
- try{ vol = solid.getCubicVolume(); } catch(Exception e) { vol = 0.0; }
- solid_vol_map.put(solid,vol);
+ try {
+ vol = solid.getCubicVolume();
+ } catch (Exception e) {
+ vol = 0.0;
+ }
+
+ solid_vol_map.put(solid, vol);
return vol;
}
+
}
-
- private VolumeInfo performVolumeCalculations(UniquePV pv){
-
+
+ private VolumeInfo performVolumeCalculations(UniquePV pv) {
+
VolumeInfo vi = new VolumeInfo();
ISolid solid = pv.getSolid();
-
+
//ASSUMPTION: tube is along z-axis and has center at r = 0
if (solid instanceof Tube) {
Tube tube = (Tube) solid;
@@ -359,11 +335,9 @@
vi.zmin = z0 - tube.getZHalfLength();
vi.rmin = tube.getInnerRadius();
vi.rmax = tube.getOuterRadius();
- }
-
- else if (solid instanceof Box) {
+ } else if (solid instanceof Box) {
Box box = (Box) solid;
- for (Point3D p : box.getVertices()){
+ for (Point3D p : box.getVertices()) {
Hep3Vector transformed = pv.localToGlobal(p.getHep3Vector());
vi.zmin = Math.min(transformed.z(), vi.zmin);
vi.zmax = Math.max(transformed.z(), vi.zmax);
@@ -371,9 +345,8 @@
vi.rmin = Math.min(vi.rmin, r);
vi.rmax = Math.max(vi.rmax, r);
}
- }
-
- //Note: this information will NOT be used most of the time...
+
+ } //Note: this information will NOT be used most of the time...
// Polycones that are top-level elements (e.g. the beampipe) are
// handled specially (since the radiation length is a function of z).
// The information here will only be used in case a top-level element
@@ -382,101 +355,99 @@
else if (solid instanceof Polycone) {
Polycone pc = (Polycone) solid;
List<Polycone.ZPlane> zplanes = pc.getZPlanes();
-
+
//For now, just take the minimum rmin and rmax of the polycone
- for (Polycone.ZPlane z : zplanes){
- if (z.getRMax()>0 && z.getRMin()>0) {
- vi.rmin = Math.min(vi.rmin,z.getRMin());
- vi.rmax = vi.rmax > 0. ? Math.min(vi.rmax,z.getRMax()) : z.getRMax();
+ for (Polycone.ZPlane z : zplanes) {
+ if (z.getRMax() > 0 && z.getRMin() > 0) {
+ vi.rmin = Math.min(vi.rmin, z.getRMin());
+ vi.rmax = vi.rmax > 0. ? Math.min(vi.rmax, z.getRMax()) : z.getRMax();
}
+
}
-
+
vi.zmin = pc.getZPlanes().get(0).getZ();
- vi.zmax = pc.getZPlanes().get(pc.getZPlanes().size()-1).getZ();
-
+ vi.zmax = pc.getZPlanes().get(pc.getZPlanes().size() - 1).getZ();
+
//check for wrong order
if (vi.zmin > vi.zmax) {
double temp = vi.zmin;
vi.zmin = vi.zmax;
vi.zmax = temp;
}
+
}
-
+
return vi;
}
-
-
-
-
+
/**
* A UniquePV is a wrapper around IPhysicalVolumePath which provides
* some convenience methods and caches transformations.
*/
- class UniquePV{
-
+ class UniquePV {
+
IPhysicalVolumePath path;
IPhysicalVolumeNavigator nav;
ITransform3D transform = null;
-
+
/**
* Generates a top-level UniquePV.
* @param root The top-level IPhysicalVolume
* @param navigator The IPhysicalVolumeNavigator associated with the detector
*/
- public UniquePV(IPhysicalVolume root, IPhysicalVolumeNavigator navigator){
+ public UniquePV(IPhysicalVolume root, IPhysicalVolumeNavigator navigator) {
path = new PhysicalVolumePath();
nav = navigator;
path.add(root);
}
-
+
/**
* Generates a UniquePV from a path. (Shallow copy of path)
* @param path
* @param navigator
*/
- public UniquePV(IPhysicalVolumePath path, IPhysicalVolumeNavigator navigator){
+ public UniquePV(IPhysicalVolumePath path, IPhysicalVolumeNavigator navigator) {
this.path = path;
nav = navigator;
}
-
+
/**
* Returns the IPhysicalVolume (the last element of the path)
*/
- public IPhysicalVolume getPV(){
+ public IPhysicalVolume getPV() {
return path.getLeafVolume();
}
-
+
/**
* Creates a UniquePV that is a daughter of the current UniquePV (deep copy made)
* @param daughter
* @return
*/
- public UniquePV createDaughterUniquePV(IPhysicalVolume daughter){
+ public UniquePV createDaughterUniquePV(IPhysicalVolume daughter) {
IPhysicalVolumePath np = new PhysicalVolumePath();
np.addAll(path);
np.add(daughter);
- return new UniquePV(np,nav);
+ return new UniquePV(np, nav);
}
-
-
+
/**
* Transforms the given vector from local to global coords.
* @param v the untransformed local Hep3Vector
* @return the transformed global Hep3Vector
*/
public Hep3Vector localToGlobal(Hep3Vector v) {
-
+
return getLtoGTransform().transformed(v);
}
-
+
/**
* Returns the solid associated with the physical volume.
* @return
*/
- public ISolid getSolid(){
+ public ISolid getSolid() {
return this.getPV().getLogicalVolume().getSolid();
}
-
+
/**
* Returns the local-to-global transform
* @return an ITransform3D from local coordinates to global coordinates.
@@ -487,13 +458,10 @@
}
return transform;
}
-
+
@Override
- public String toString(){
+ public String toString() {
return path.toString();
}
}
-
-
-
}