lcsim/src/org/lcsim/contrib/uiowa
diff -u -r1.9 -r1.10
--- MIPReassignmentAlgorithm.java 22 Jul 2008 18:52:02 -0000 1.9
+++ MIPReassignmentAlgorithm.java 25 Jul 2008 22:44:34 -0000 1.10
@@ -32,40 +32,40 @@
import hep.aida.IHistogram1D;
import hep.aida.ICloud1D;
import hep.aida.ICloud2D;
+import org.lcsim.util.aida.AIDA;
import org.lcsim.event.EventHeader;
-public class MIPReassignmentAlgorithm implements ReassignClustersAlgorithm {
+
+ public class MIPReassignmentAlgorithm implements ReassignClustersAlgorithm {
+ private AIDA aida = AIDA.defaultInstance();
protected double m_limit;
- protected LocalHelixExtrapolator m_findCluster;
- protected boolean nomip = false;
+ protected LocalHelixExtrapolator m_extrap;
protected boolean debug = false;
protected double m_distance = 0;
protected double m_angle = 0;
protected double m_degree = 0;
protected double m_costheta = 0;
protected double m_trans = 0;
-
- protected ICloud2D m_histo_radian_distance;
- protected ICloud2D m_histo_degree_distance;
- protected ICloud2D m_histo_costheta_distance;
- protected ICloud2D m_histo_trans_longi;
-
+ protected double m_height = 0;
+ protected double m_cos_two = 0;
+
+ protected int m_shape;
+
protected Map<Track, BasicCluster> m_newMapMIP;
- public MIPReassignmentAlgorithm(double limit, LocalHelixExtrapolator findCluster, Map<Track, BasicCluster> MapTrkToMIP, ICloud2D histo1, ICloud2D histo2, ICloud2D histo3, ICloud2D histo4) {
+ public MIPReassignmentAlgorithm(double limit, LocalHelixExtrapolator findCluster, Map<Track, BasicCluster> MapTrkToMIP, int shape) {
+
m_limit = limit;
- m_findCluster = findCluster;
+ m_extrap = findCluster;
m_newMapMIP = MapTrkToMIP;
- m_histo_radian_distance = histo1;
- m_histo_degree_distance = histo2;
- m_histo_costheta_distance = histo3;
- m_histo_trans_longi = histo4;
+ m_shape = shape;
}
public MIPReassignmentAlgorithm(double limit, LocalHelixExtrapolator findCluster, Map<Track, BasicCluster> MapTrkToMIP) {
m_limit = limit;
- m_findCluster = findCluster;
+ m_extrap = findCluster;
m_newMapMIP = MapTrkToMIP;
+ m_shape = 7;
}
@@ -98,14 +98,15 @@
Hep3Vector clusterPosition = new BasicHep3Vector(); //cluster position
Hep3Vector displacement = new BasicHep3Vector(); //displacement between cluster and last point
Hep3Vector displacementUnit = new BasicHep3Vector();
- Hep3Vector interceptPoint = m_findCluster.performExtrapolation(tr);
-
+ Hep3Vector interceptPoint = m_extrap.performExtrapolation(tr);
+ Hep3Vector P = new BasicHep3Vector(tr.getMomentum());
+
BasicCluster newmip = m_newMapMIP.get(tr);
List<Hep3Vector> lastTwoPositions = new ArrayList<Hep3Vector>();
if(newmip.getCalorimeterHits().size() < 5){ //These is no mip close to track.
if (interceptPoint != null) {
- tangent = m_findCluster.getTangent(); //tangent obtained using extrapolated track.
+ tangent = m_extrap.getTangent(); //tangent obtained using extrapolated track.
tangentUnit = VecOp.unit(tangent);
clusterPosition = new BasicHep3Vector(clus.getPosition());
displacement = VecOp.sub(clusterPosition, interceptPoint);
@@ -124,17 +125,31 @@
}
Hep3Vector last0pos = lastTwoPositions.get(0);
Hep3Vector last1pos = lastTwoPositions.get(1);
-
- tangent = VecOp.sub(last0pos, last1pos); //tangent obtained using last two points.
+
+ //Chose one
+ Hep3Vector extrapTangent = m_extrap.getTangent(last0pos, tr);
+ Hep3Vector extrapTangentUnit = VecOp.unit(extrapTangent);
+ Hep3Vector twoPointsTangent = VecOp.sub(last0pos, last1pos);
+ Hep3Vector twoPointsTangentUnit = VecOp.unit(twoPointsTangent);
+
+ boolean usingExtrap = true;
+ if(usingExtrap) tangent = extrapTangent; //tangent obtained using extrapoltated track
+ else tangent = twoPointsTangent;
+
tangentUnit = VecOp.unit(tangent);
clusterPosition = new BasicHep3Vector(clus.getPosition());
displacement = VecOp.sub(clusterPosition, last1pos);
displacementUnit = VecOp.unit(displacement);
+ //For debug
+ m_cos_two = VecOp.dot(twoPointsTangentUnit, extrapTangentUnit);
+
if(debug) {
System.out.print("Debug: From last two points (" + last0pos.magnitude() + ", " +last1pos.magnitude()+ " )");
System.out.print(" Displacement from (r="+Math.sqrt(last1pos.x()*last1pos.x()+last1pos.y()*last1pos.y())+", z="+last1pos.z()+")");
- System.out.print(" to (r="+Math.sqrt(clusterPosition.x()*clusterPosition.x()+clusterPosition.y()*clusterPosition.y())+", z="+clusterPosition.z()+")");
- }
+ System.out.println(" to (r="+Math.sqrt(clusterPosition.x()*clusterPosition.x()+clusterPosition.y()*clusterPosition.y())+", z="+clusterPosition.z()+")");
+ System.out.println("Costheta between two tangents= " + m_cos_two + " at " + P.magnitude() + " GeV");
+
+ }
}
m_distance=displacement.magnitude()*VecOp.dot(tangentUnit, displacementUnit); //longitudinal distance
@@ -143,19 +158,49 @@
m_trans = displacement.magnitude()*Math.sin(m_angle);
m_costheta = VecOp.dot(tangentUnit, displacementUnit);
- if (m_histo_radian_distance != null) { m_histo_radian_distance.fill(m_angle,m_distance); }
- if (m_histo_degree_distance != null) { m_histo_degree_distance.fill(m_degree,m_distance); }
- if (m_histo_costheta_distance != null) { m_histo_costheta_distance.fill(m_costheta,m_distance); }
- if (m_histo_trans_longi != null) { m_histo_trans_longi.fill(m_trans, m_distance); }
-
- if(debug){
- System.out.print(" Angle= " + m_angle);
- System.out.println(" Distance= " + m_distance);
+ aida.histogram1D("Track P distribution", 200, 0., 200.).fill(P.magnitude());
+ aida.cloud2D("cos between two vectors", 10000).fill(P.magnitude(),m_cos_two);
+ aida.cloud2D("m_histo_radian_distance",10000).fill(m_angle,m_distance);
+ aida.cloud2D("m_histo_degree_distance",10000).fill(m_degree,m_distance);
+ aida.cloud2D("m_histo_costheta_distance",10000).fill(m_costheta,m_distance);
+ aida.cloud2D("m_histo_trans_longi",10000).fill(m_trans, m_distance);
+ aida.cloud2D("m_histo_degree_trans",10000).fill(m_degree,m_trans);
+ aida.cloud2D("m_histo_p_ratio",10000).fill(P.magnitude(), m_distance/m_trans);
+ aida.cloud2D("m_histo_p_distance",10000).fill(P.magnitude(), m_distance);
+ aida.cloud2D("m_histo_p_trans",10000).fill(P.magnitude(), m_trans);
+
+ if ( P.magnitude() < 5 ) m_height = 200;
+ if ( P.magnitude() >= 5 && P.magnitude() < 30 ) m_height = 600;
+ if ( P.magnitude() >= 30) m_height = 1200;
+
+ switch (m_shape) {
+ case 1: //CylinderWithHeight_Angul
+ if (m_trans < m_limit && m_angle < 1.57 && m_distance < m_height) return new Double(m_angle);
+ break;
+ case 2://CylinderWithHeight_Trans
+ if (m_trans < m_limit && m_angle < 1.57 && m_distance < m_height) return new Double(m_trans);
+ break;
+ case 3://Cylinder_Angul
+ if (m_trans < m_limit && m_angle < 1.57 ) return new Double(m_angle);
+ break;
+ case 4://Cylinder_Trans
+ if (m_trans < m_limit && m_angle < 1.57 ) return new Double(m_trans);
+ break;
+ case 5://ConeWithHeight_Angul
+ if (m_angle < m_limit && m_distance < m_height) return new Double(m_angle);
+ break;
+ case 6://ConeWithHeight_Trans
+ if (m_angle < m_limit && m_distance < m_height) return new Double(m_trans);
+ break;
+ case 7://Cone_Angul
+ if (m_angle < m_limit) return new Double(m_angle);
+ break;
+ case 8://Cone_Trans
+ if (m_angle < m_limit) return new Double(m_trans);
+ break;
}
- if (m_angle < m_limit) return new Double(m_angle); // return the angle if it is below cone angle.
-
- return null;
+ return null;
}
}
