lcsim/src/org/lcsim/recon/pfa/structural
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;
}
}