lcsim/src/org/lcsim/contrib/uiowa
diff -u -r1.10 -r1.11
--- MIPReassignmentAlgorithm.java 25 Jul 2008 22:44:34 -0000 1.10
+++ MIPReassignmentAlgorithm.java 12 Aug 2008 22:45:06 -0000 1.11
@@ -34,24 +34,24 @@
import hep.aida.ICloud2D;
import org.lcsim.util.aida.AIDA;
import org.lcsim.event.EventHeader;
+import org.lcsim.contrib.uiowa.ReclusterDriver;
+ public class MIPReassignmentAlgorithm extends ReclusterDriver implements ReassignClustersAlgorithm {
+ private double m_limit;
+ private LocalHelixExtrapolator m_extrap;
+ private boolean debug = false;
+ private double m_dist = 0;
+ private double m_longi = 0;
+ private double m_angle = 0;
+ private double m_degree = 0;
+ private double m_costheta = 0;
+ private double m_trans = 0;
+ private double m_height = 0;
+ private double m_cos_two = 0;
- public class MIPReassignmentAlgorithm implements ReassignClustersAlgorithm {
- private AIDA aida = AIDA.defaultInstance();
- protected double m_limit;
- 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 double m_height = 0;
- protected double m_cos_two = 0;
+ private int m_shape;
- protected int m_shape;
-
- protected Map<Track, BasicCluster> m_newMapMIP;
+ private Map<Track, BasicCluster> m_newMapMIP;
public MIPReassignmentAlgorithm(double limit, LocalHelixExtrapolator findCluster, Map<Track, BasicCluster> MapTrkToMIP, int shape) {
@@ -74,7 +74,6 @@
//Check extrapolation for all the tracks and
// take the best score among them.
if ( tr.getTracks().size() != 0) {
- System.out.println("multiple");
Double bestAmongDaughters = null;
for (Track daughterTrack : tr.getTracks()) {
Double daughterFigureOfMerit = this.computeFigureOfMerit(daughterTrack, clus);
@@ -92,26 +91,20 @@
// All done
return bestAmongDaughters;
}
-
+
+ Hep3Vector last0pos = new BasicHep3Vector();
+ Hep3Vector last1pos = new BasicHep3Vector();
Hep3Vector tangent = new BasicHep3Vector(); //extrapolated direction
- Hep3Vector tangentUnit = new BasicHep3Vector();
- Hep3Vector clusterPosition = new BasicHep3Vector(); //cluster position
- Hep3Vector displacement = new BasicHep3Vector(); //displacement between cluster and last point
- Hep3Vector displacementUnit = new BasicHep3Vector();
- 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_extrap.getTangent(); //tangent obtained using extrapolated track.
- tangentUnit = VecOp.unit(tangent);
- clusterPosition = new BasicHep3Vector(clus.getPosition());
- displacement = VecOp.sub(clusterPosition, interceptPoint);
- displacementUnit = VecOp.unit(displacement);
- }
+ if(newmip.getCalorimeterHits().size() < 3){ //These is no mip close to track. There should at least two hits.
+ last0pos = m_extrap.performExtrapolation(tr);
+ if (last0pos != null) {
+ tangent = m_extrap.getTangent(last0pos,tr); //tangent obtained using extrapolated track.
+ }
}else{ //There is a mip attached to track.
for(int i=0; i<2 ; i++){
CalorimeterHit hit = newmip.getCalorimeterHits().get(newmip.getCalorimeterHits().size()-1-i);
@@ -123,74 +116,51 @@
Hep3Vector v = new BasicHep3Vector(x,y,z);
lastTwoPositions.add(v);
}
- Hep3Vector last0pos = lastTwoPositions.get(0);
- Hep3Vector last1pos = lastTwoPositions.get(1);
-
+ last0pos = lastTwoPositions.get(0);
+ last1pos = lastTwoPositions.get(1);
//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.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");
-
- }
+ if(usingExtrap) tangent = m_extrap.getTangent(last0pos, tr); //tangent obtained using extrapoltated track
+ else tangent = VecOp.sub(last0pos, last1pos);
}
-
- m_distance=displacement.magnitude()*VecOp.dot(tangentUnit, displacementUnit); //longitudinal distance
+
+ Hep3Vector tangentUnit = VecOp.unit(tangent);
+ Hep3Vector temp = VecOp.mult(0,tangentUnit);
+ //going back 500 mm or 800 mm back to get the last1pos
+ Hep3Vector imgpos = VecOp.sub(last0pos,temp);
+ Hep3Vector clusterPosition = new BasicHep3Vector(clus.getPosition());
+ Hep3Vector displacement = VecOp.sub(clusterPosition, imgpos);
+ Hep3Vector displacementUnit = VecOp.unit(displacement);
+
+ m_dist = displacement.magnitude();
+ m_longi=displacement.magnitude()*VecOp.dot(tangentUnit, displacementUnit); //longitudinal distance
m_angle=Math.acos(VecOp.dot(tangentUnit, displacementUnit));
- m_degree=m_angle*180/3.14;
+ m_degree=m_angle*180/Math.PI;
m_trans = displacement.magnitude()*Math.sin(m_angle);
m_costheta = VecOp.dot(tangentUnit, displacementUnit);
-
- 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);
+ if (m_trans < m_limit && m_longi > 0 && m_longi < 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);
+ if (m_trans < m_limit && m_longi < 0 && m_longi < 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);
+ if (m_trans < m_limit && m_longi > 0 ) return new Double(m_angle);
break;
case 4://Cylinder_Trans
- if (m_trans < m_limit && m_angle < 1.57 ) return new Double(m_trans);
+ if (m_trans < m_limit && m_longi > 0 ) return new Double(m_trans);
break;
case 5://ConeWithHeight_Angul
- if (m_angle < m_limit && m_distance < m_height) return new Double(m_angle);
+ if (m_angle < m_limit && m_longi < 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);
+ if (m_angle < m_limit && m_longi < m_height) return new Double(m_trans);
break;
case 7://Cone_Angul
if (m_angle < m_limit) return new Double(m_angle);
@@ -198,6 +168,14 @@
case 8://Cone_Trans
if (m_angle < m_limit) return new Double(m_trans);
break;
+ case 9://Cone_dist
+ if (m_angle < m_limit) return new Double(m_dist);
+ case 10://Hemisphere
+ if (m_dist < m_limit) return new Double(m_dist);
+ case 11://Complex
+ if (m_angle < m_limit && m_longi < m_height && m_trans < 200) return new Double(m_angle);
+ case 12://Complex_backward
+ if ((m_angle < m_limit && m_longi < m_height && m_trans < 200) || (m_trans < 100 && m_longi > -200 && m_longi < 0)) return new Double(m_angle);
}
return null;