Print

Print


Commit in lcsim/src/org/lcsim/recon/pfa/structural on MAIN
ChargedHadronClusterEnergyCalculator.java+64-441.2 -> 1.3
MJC: Make dE/dx estimation for MIP hits more general -- should work now for [some] geometries other than sid01

lcsim/src/org/lcsim/recon/pfa/structural
ChargedHadronClusterEnergyCalculator.java 1.2 -> 1.3
diff -u -r1.2 -r1.3
--- ChargedHadronClusterEnergyCalculator.java	27 Mar 2008 18:57:25 -0000	1.2
+++ ChargedHadronClusterEnergyCalculator.java	27 Mar 2008 23:34:07 -0000	1.3
@@ -10,6 +10,17 @@
 import org.lcsim.event.*;
 import org.lcsim.event.util.*;
 
+import org.lcsim.detector.IDetectorElement;
+import org.lcsim.detector.IDetectorElementContainer;
+import org.lcsim.detector.IParameters;
+import org.lcsim.detector.IGeometryInfo;
+import org.lcsim.detector.ILogicalVolume;
+import org.lcsim.detector.IPhysicalVolume;
+import org.lcsim.detector.IPhysicalVolumeContainer;
+import org.lcsim.detector.material.IMaterial;
+import org.lcsim.detector.material.IMaterialMixture;
+import org.lcsim.detector.solids.ISolid;
+
 /**
   * ECAL/HCAL calibration designed for charged showers. The hits are divided into two classes: those
   * which appear to come from MIPs, and all others. The MIP hits are assigned energy based on naive
@@ -17,13 +28,11 @@
   * Since the number of non-MIP hits could be small or zero, the neutral calibration should not have
   * a large zero offset.
   *
-  * The dE/dx energy calculation is currently very crude and hard-coded to the sid01 geometry.
-  * Some improvements needed:
-  *   1) Read in geometry and material information from detector at run-time;
-  *   2) Increase energy loss for tracks that pass through the layer at an angle (currently
+  * Still to fix:
+  *   1) Increase energy loss for tracks that pass through the layer at an angle (currently
   *      it assumes all tracks are perpendicular to the material planes).
   *
-  * @version $Id: ChargedHadronClusterEnergyCalculator.java,v 1.2 2008/03/27 18:57:25 mcharles Exp $
+  * @version $Id: ChargedHadronClusterEnergyCalculator.java,v 1.3 2008/03/27 23:34:07 mcharles Exp $
   */
 
 public class ChargedHadronClusterEnergyCalculator extends Driver implements ClusterEnergyCalculator
@@ -69,7 +78,7 @@
 	}
 	return eSum;
     }
-    // This is ugly for now!
+
     private double getMIPEnergy(CalorimeterHit hit) {
 	// First, grab hit info (ID, layer)
 	org.lcsim.geometry.IDDecoder id = hit.getIDDecoder();
@@ -78,49 +87,60 @@
 	// Next, figure out which subdetector we're in.
 	org.lcsim.geometry.Subdetector subdet = hit.getSubdetector();
 	if ( ! subdet.isCalorimeter() ) { throw new AssertionError("Cluster hit outside calorimeter"); }
-	String name = subdet.getName();
 	// Set up dE/dx constants:
 	org.lcsim.util.step.DeDx materialDatabase = org.lcsim.util.step.DeDx.instance();
-	double dedx_Silicon = materialDatabase.getDeDx("Silicon");
-	double dedx_Copper = materialDatabase.getDeDx("Copper");
-	double dedx_Kapton = materialDatabase.getDeDx("Kapton");
-	double dedx_Air = materialDatabase.getDeDx("Air");
-	double dedx_Tungsten = materialDatabase.getDeDx("Tungsten");
-	double dedx_Steel235 = materialDatabase.getDeDx("Steel235");
-	double dedx_G10 = materialDatabase.getDeDx("G10");
-	double dedx_RPC = materialDatabase.getDeDx("RPCGasDefault");
-	double dedx_PyrexGlass = materialDatabase.getDeDx("PyrexGlass");
-	// Compute energy loss:
-	double energyLoss = 0.0;
-	if (name.compareTo("EMBarrel")==0 || name.compareTo("EMEndcap")==0) {
-	    if (layer == 0) {
-		energyLoss += 0.032 * dedx_Silicon;
-		energyLoss += 0.005 * dedx_Copper;
-		energyLoss += 0.030 * dedx_Kapton;
-		energyLoss += 0.033 * dedx_Air;
-	    } else if (layer > 0 && layer <= 20) {
-		energyLoss += 0.250 * dedx_Tungsten;
-		energyLoss += 0.032 * dedx_Silicon;
-		energyLoss += 0.005 * dedx_Copper;
-		energyLoss += 0.030 * dedx_Kapton;
-		energyLoss += 0.033 * dedx_Air;
-	    } else if (layer > 20 && layer <= 30) {
-		energyLoss += 0.500 * dedx_Tungsten;
-		energyLoss += 0.032 * dedx_Silicon;
-		energyLoss += 0.005 * dedx_Copper;
-		energyLoss += 0.030 * dedx_Kapton;
-		energyLoss += 0.033 * dedx_Air;
+
+	// 1) There is one SUBDETECTOR -- subdetectorElement
+	//      * If the SUBDETECTOR is a barrel, it has one corresponding SUBDETECTOR LOGICAL VOLUME
+	//      * If the SUBDETECTOR is an endcap, it has two daughter detector elements (for the +z
+	//           and -z endcaps), each of which has one corresponding SUBDETECTOR LOGICAL VOLUME.
+	//           Since the endcaps are symmetric, we pick the first one.
+	// 2) The SUBDETECTOR LOGICAL VOLUME has multiple LAYER PHYSICAL VOLUME daughters (e.g. 31 for ecal, 34 for hcal)
+	// 3) Each LAYER PHYSICAL VOLUME has a corresponding LAYER LOGICAL VOLUME
+	// 4) Each LAYER LOGICAL VOLUME has multiple SUBLAYER PHYSICAL VOLUME daughters (e.g. TungstenDens24, Silicon, Copper, Kapton, Air)
+	//
+	// Because of the way the ordering works, the total material since the last active element
+	//    is the same as the sum of the physical volume daughters in the current layer *except*
+	//    for the very first layer in the subdetector, for which we ignore anything that comes
+	//    after the active element. Neglect the difference for that layer for now.
+
+	double estimatedEnergyLoss = 0.0;
+
+	IDetectorElement subDetectorElement = subdet.getDetectorElement();	
+	IGeometryInfo subDetectorGeometry = null;
+	if (!subdet.isEndcap()) {
+	    // Barrel -- just one geometry:
+	    subDetectorGeometry = subDetectorElement.getGeometry();
+	} else {
+	    // Endcap -- extra level of indirection
+	    IDetectorElement oneEndcap = subDetectorElement.getChildren().get(0);
+	    subDetectorGeometry = oneEndcap.getGeometry();
+	}
+	ILogicalVolume subDetectorLogicalVolume = subDetectorGeometry.getLogicalVolume();
+	IPhysicalVolume layerPhysicalVolume = subDetectorLogicalVolume.getDaughter(layer);
+	ILogicalVolume layerLogicalVolume = layerPhysicalVolume.getLogicalVolume();
+	for (IPhysicalVolume subLayerPhysicalVolume : layerLogicalVolume.getDaughters()) {
+	    ILogicalVolume subLayerLogicalVolume = subLayerPhysicalVolume.getLogicalVolume();
+	    IMaterial subLayerMaterial = subLayerLogicalVolume.getMaterial();
+	    double estimated_dedx = materialDatabase.getDeDx(subLayerMaterial.getName());
+	    ISolid subLayerSolid = subLayerLogicalVolume.getSolid();
+	    if (subLayerSolid instanceof org.lcsim.detector.solids.Tube) {
+		org.lcsim.detector.solids.Tube subLayerTube = (org.lcsim.detector.solids.Tube)(subLayerSolid);
+		double layerThickness;
+		if (subdet.isEndcap()) {
+		    layerThickness = 2.0 * subLayerTube.getZHalfLength();
+		} else {
+		    layerThickness = subLayerTube.getOuterRadius() - subLayerTube.getInnerRadius();
+		}
+
+		// Watch out for factor of 10 -- this comes from the dE/dx being in units of MeV/cm and distance being in units of mm
+		estimatedEnergyLoss += (layerThickness * estimated_dedx * 0.1);
+	    } else {
+		throw new AssertionError("Don't know how to get layer thickness for geometry type "+subLayerSolid.getClass().getName());
 	    }
-	} else if (name.compareTo("HADBarrel")==0 || name.compareTo("HADEndcap")==0) {
-	    energyLoss += 2.00 * dedx_Steel235;
-	    energyLoss += 0.30 * dedx_G10;
-	    energyLoss += 0.11 * dedx_PyrexGlass;
-	    energyLoss += 0.12 * dedx_RPC;
-	    energyLoss += 0.11 * dedx_PyrexGlass;
-	    energyLoss += 0.16 * dedx_Air;
 	}
 
-	return energyLoss;
+	return estimatedEnergyLoss;
     }
 
 }
CVSspam 0.2.8