lcsim/src/org/lcsim/recon/pfa/identifier
diff -u -r1.2 -r1.3
--- LocalHelixExtrapolator.java 13 Jun 2008 01:03:04 -0000 1.2
+++ LocalHelixExtrapolator.java 20 Jun 2008 18:53:10 -0000 1.3
@@ -37,13 +37,21 @@
protected boolean m_init = false;
protected double m_ECAL_barrel_zmin;
+ protected double m_HCAL_barrel_zmin;
protected double m_ECAL_barrel_zmax;
+ protected double m_HCAL_barrel_zmax;
protected double m_ECAL_barrel_r;
+ protected double m_HCAL_barrel_r;
protected double m_ECAL_endcap_z;
+ protected double m_HCAL_endcap_z;
protected Vector<Double> m_ECAL_barrel_layering_r;
+ protected Vector<Double> m_HCAL_barrel_layering_r;
protected Vector<Double> m_ECAL_endcap_layering_z;
+ protected Vector<Double> m_HCAL_endcap_layering_z;
protected double m_ECAL_endcap_rmin;
+ protected double m_HCAL_endcap_rmin;
protected double m_ECAL_endcap_rmax;
+ protected double m_HCAL_endcap_rmax;
protected double[] m_fieldStrength;
boolean m_barrelValid = false;
boolean m_endcapValid = false;
@@ -70,12 +78,20 @@
Detector det = event.getDetector();
CylindricalCalorimeter emb = ((CylindricalCalorimeter) det.getSubdetectors().get("EMBarrel"));
CylindricalCalorimeter eme = ((CylindricalCalorimeter) det.getSubdetectors().get("EMEndcap"));
+ CylindricalCalorimeter hdb = ((CylindricalCalorimeter) det.getSubdetectors().get("HADBarrel"));
+ CylindricalCalorimeter hde = ((CylindricalCalorimeter) det.getSubdetectors().get("HADEndcap"));
m_ECAL_barrel_zmin = emb.getZMin();
+ m_HCAL_barrel_zmin = hdb.getZMin();
m_ECAL_barrel_zmax = emb.getZMax();
+ m_HCAL_barrel_zmax = hdb.getZMax();
m_ECAL_barrel_r = emb.getLayering().getDistanceToLayerSensorMid(0);
+ m_HCAL_barrel_r = hdb.getLayering().getDistanceToLayerSensorMid(0);
m_ECAL_endcap_z = eme.getLayering().getDistanceToLayerSensorMid(0);
+ m_HCAL_endcap_z = hde.getLayering().getDistanceToLayerSensorMid(0);
m_ECAL_endcap_rmin = eme.getInnerRadius();
+ m_HCAL_endcap_rmin = hde.getInnerRadius();
m_ECAL_endcap_rmax = eme.getOuterRadius();
+ m_HCAL_endcap_rmax = hde.getOuterRadius();
double[] zero = {0, 0, 0};
m_fieldStrength = det.getFieldMap().getField(zero);
m_ECAL_barrel_layering_r = new Vector<Double>();
@@ -88,6 +104,16 @@
double z = eme.getLayering().getDistanceToLayerSensorMid(iLayer);
m_ECAL_endcap_layering_z.add(new Double(z));
}
+ m_HCAL_barrel_layering_r = new Vector<Double>();
+ for (int iLayer=0; iLayer<hdb.getLayering().getLayerCount(); iLayer++) {
+ double r = hdb.getLayering().getDistanceToLayerSensorMid(iLayer);
+ m_HCAL_barrel_layering_r.add(new Double(r));
+ }
+ m_HCAL_endcap_layering_z = new Vector<Double>();
+ for (int iLayer=0; iLayer<hde.getLayering().getLayerCount(); iLayer++) {
+ double z = hde.getLayering().getDistanceToLayerSensorMid(iLayer);
+ m_HCAL_endcap_layering_z.add(new Double(z));
+ }
m_init = true;
}
@@ -332,8 +358,8 @@
return tangent;
}
- /** Assumes extrapolation has already been done. */
- public Hep3Vector extendToECALLayer(int layer) {
+ private Hep3Vector extendToBarrelLayer(int layer, Vector<Double> barrel_layering_r, double barrel_zmin, double barrel_zmax )
+ {
double dphi = (m_track_phi0 - m_track_phi1);
while (dphi < -Math.PI) { dphi += 2.0*Math.PI; }
while (dphi > Math.PI) { dphi -= 2.0*Math.PI; }
@@ -342,15 +368,14 @@
double radius = m_trackParam_radius;
double dz_by_dphi = m_trackParam_dz_by_dphi;
- if (m_barrelValid) {
- double a = radius;
- double b = Math.sqrt(x_c*x_c + y_c*y_c);
- double r = m_ECAL_barrel_layering_r.get(layer);
- double cos_phi_minus_theta = (r*r - a*a - b*b) / (2.0*a*b); // Obviously, this blows up if a or b is zero
- double theta = Math.atan2(y_c, x_c); // in the range (-pi, +pi)
- double dphi_to_barrel = 0.0;
- if (cos_phi_minus_theta < -1.0) {
- // No solution
+ double a = radius;
+ double b = Math.sqrt(x_c*x_c + y_c*y_c);
+ double r = barrel_layering_r.get(layer);
+ double cos_phi_minus_theta = (r*r - a*a - b*b) / (2.0*a*b); // Obviously, this blows up if a or b is zero
+ double theta = Math.atan2(y_c, x_c); // in the range (-pi, +pi)
+ double dphi_to_barrel = 0.0;
+ if (cos_phi_minus_theta < -1.0) {
+ // No solution
return null;
} else if (cos_phi_minus_theta == -1.0) {
// Unique solution: phi = theta + pi
@@ -396,29 +421,60 @@
double found_barrel_y = y_c + radius * Math.sin(found_barrel_phi);
double found_barrel_polar_r = Math.sqrt(found_barrel_x*found_barrel_x + found_barrel_y*found_barrel_y);
double found_barrel_polar_phi = Math.atan2(found_barrel_y, found_barrel_x);
- boolean validSolution =(found_barrel_z >= m_ECAL_barrel_zmin-m_cutSeparation && found_barrel_z <= m_ECAL_barrel_zmax+m_cutSeparation);
+ boolean validSolution =(found_barrel_z >= barrel_zmin-m_cutSeparation && found_barrel_z <= barrel_zmax+m_cutSeparation);
if (validSolution) {
return new BasicHep3Vector(found_barrel_x, found_barrel_y, found_barrel_z);
} else {
return null;
}
+ }
+
+ private Hep3Vector extendToEndcapLayer(int layer, Vector<Double> endcap_layering_z, double endcap_rmin, double endcap_rmax )
+ {
+ double dphi = (m_track_phi0 - m_track_phi1);
+ while (dphi < -Math.PI) { dphi += 2.0*Math.PI; }
+ while (dphi > Math.PI) { dphi -= 2.0*Math.PI; }
+ double x_c = m_trackParam_xc;
+ double y_c = m_trackParam_yc;
+ double radius = m_trackParam_radius;
+ double dz_by_dphi = m_trackParam_dz_by_dphi;
+
+ double previous_z = m_trackPoint_z;
+ double new_z = Math.abs(endcap_layering_z.get(layer));
+ if (previous_z < 0) { new_z = -new_z; }
+ double dz = new_z - previous_z;
+ double deltaPhi = dz / m_trackParam_dz_by_dphi;
+ double phi = m_trackPoint_phi + deltaPhi;
+ double found_endcap_x = x_c + radius * Math.cos(phi);
+ double found_endcap_y = y_c + radius * Math.sin(phi);
+ double found_endcap_z = new_z;
+ double found_endcap_polar_r = Math.sqrt(found_endcap_x*found_endcap_x + found_endcap_y*found_endcap_y);
+ boolean validSolution = (found_endcap_polar_r >= endcap_rmin-m_cutSeparation && found_endcap_polar_r <= endcap_rmax+m_cutSeparation);
+ if (validSolution) {
+ return new BasicHep3Vector(found_endcap_x, found_endcap_y, found_endcap_z);
+ } else {
+ return null;
+ }
+ }
+
+ /** Assumes extrapolation has already been done. */
+ public Hep3Vector extendToECALLayer(int layer) {
+ if (m_barrelValid) {
+ return extendToBarrelLayer(layer, m_ECAL_barrel_layering_r, m_ECAL_barrel_zmin, m_ECAL_barrel_zmax);
} else if (m_endcapValid) {
- double previous_z = m_trackPoint_z;
- double new_z = Math.abs(m_ECAL_endcap_layering_z.get(layer));
- if (previous_z < 0) { new_z = -new_z; }
- double dz = new_z - previous_z;
- double deltaPhi = dz / m_trackParam_dz_by_dphi;
- double phi = m_trackPoint_phi + deltaPhi;
- double found_endcap_x = x_c + radius * Math.cos(phi);
- double found_endcap_y = y_c + radius * Math.sin(phi);
- double found_endcap_z = new_z;
- double found_endcap_polar_r = Math.sqrt(found_endcap_x*found_endcap_x + found_endcap_y*found_endcap_y);
- boolean validSolution = (found_endcap_polar_r >= m_ECAL_endcap_rmin-m_cutSeparation && found_endcap_polar_r <= m_ECAL_endcap_rmax+m_cutSeparation);
- if (validSolution) {
- return new BasicHep3Vector(found_endcap_x, found_endcap_y, found_endcap_z);
- } else {
- return null;
- }
+ return extendToEndcapLayer(layer, m_ECAL_endcap_layering_z, m_ECAL_endcap_rmin, m_ECAL_endcap_rmax);
+ } else {
+ // No solution
+ return null;
+ }
+ }
+
+ /** Assumes extrapolation has already been done. */
+ public Hep3Vector extendToHCALLayer(int layer) {
+ if (m_barrelValid) {
+ return extendToBarrelLayer(layer, m_HCAL_barrel_layering_r, m_HCAL_barrel_zmin, m_HCAL_barrel_zmax);
+ } else if (m_endcapValid) {
+ return extendToEndcapLayer(layer, m_HCAL_endcap_layering_z, m_HCAL_endcap_rmin, m_HCAL_endcap_rmax);
} else {
// No solution
return null;