lcsim/src/org/lcsim/contrib/uiowa
diff -u -r1.3 -r1.4
--- MIPReassignmentAlgorithm.java 8 Jul 2008 17:26:58 -0000 1.3
+++ MIPReassignmentAlgorithm.java 10 Jul 2008 20:35:54 -0000 1.4
@@ -32,14 +32,11 @@
import hep.aida.IHistogram1D;
import hep.aida.ICloud1D;
import hep.aida.ICloud2D;
+import org.lcsim.event.EventHeader;
public class MIPReassignmentAlgorithm implements ReassignClustersAlgorithm {
protected double m_limit;
protected LocalHelixExtrapolator m_findCluster;
- protected Map<Track, Set<Cluster>> m_newMapTrackToShowerComponents;
- protected List<Cluster> m_mips;
- protected Set<CalorimeterHit> m_allhits;
- protected int exam = 0;
protected boolean nomip = false;
protected boolean debug = true;
protected double m_distance = 0;
@@ -52,23 +49,13 @@
protected ICloud2D m_histo_degree_distance;
protected ICloud2D m_histo_costheta_distance;
protected ICloud2D m_histo_trans_longi;
- //sid01
- protected double ecalRs=1270; //ECAL inner R
- protected double ecalRe=1410; //ECAL outer R
- protected double hcalRs=1410; //HCAL inner R
- protected double hcalRe=2362; //HCAL outer R
- protected double ecalZs=1680; //ECAL inner Z
- protected double ecalZe=1820; //ECAL outer Z
- protected double hcalZs=1820; //HCAL inner Z
- protected double hcalZe=2772; //HCAL inner Z
+ protected Map<Track, BasicCluster> m_newMapMIP;
- public MIPReassignmentAlgorithm(double limit, LocalHelixExtrapolator findCluster, Map<Track, Set<Cluster>> newMapTrackToShowerComponents, List<Cluster> mips, Set<CalorimeterHit> allHits, ICloud2D histo1, ICloud2D histo2, ICloud2D histo3, ICloud2D histo4) {
+ public MIPReassignmentAlgorithm(double limit, LocalHelixExtrapolator findCluster, Map<Track, BasicCluster> MapTrkToMIP, ICloud2D histo1, ICloud2D histo2, ICloud2D histo3, ICloud2D histo4) {
m_limit = limit;
m_findCluster = findCluster;
- m_newMapTrackToShowerComponents = newMapTrackToShowerComponents;
- m_mips = mips;
- m_allhits = allHits;
+ m_newMapMIP = MapTrkToMIP;
m_histo_radian_distance = histo1;
m_histo_degree_distance = histo2;
m_histo_costheta_distance = histo3;
@@ -84,9 +71,22 @@
Hep3Vector displacement = new BasicHep3Vector(); //displacement between cluster and last point
Hep3Vector displacementUnit = new BasicHep3Vector();
Hep3Vector interceptPoint = m_findCluster.performExtrapolation(tr);
- List<Hep3Vector> lastTwoPositions = computePointsOfShower(tr);
- if(lastTwoPositions == null){ //These is no mip close to track.
+ BasicCluster newmip = m_newMapMIP.get(tr);
+ List<Hep3Vector> lastTwoPositions = new ArrayList<Hep3Vector>();
+ for(int i=0; i<2 ; i++){
+ CalorimeterHit hit = newmip.getCalorimeterHits().get(newmip.getCalorimeterHits().size()-1-i);
+ IDDecoder id = hit.getIDDecoder();
+ id.setID(hit.getCellID());
+ double x = id.getPosition()[0];
+ double y = id.getPosition()[1];
+ double z = id.getPosition()[2];
+ Hep3Vector v = new BasicHep3Vector(x,y,z);
+ lastTwoPositions.add(v);
+ }
+
+
+ if(newmip.getCalorimeterHits().size() < 5 || newmip == null){ //These is no mip close to track.
if (interceptPoint != null) {
tangent = m_findCluster.getTangent(); //tangent obtained using extrapolated track.
tangentUnit = VecOp.unit(tangent);
@@ -97,8 +97,6 @@
}else{ //There is a mip attached to track.
Hep3Vector last0pos = lastTwoPositions.get(0);
Hep3Vector last1pos = lastTwoPositions.get(1);
- if(debug) System.out.println("Debug: last0pos= " + last0pos.magnitude() );
- if(debug) System.out.println("Debug: last1pos= " + last1pos.magnitude() );
tangent = VecOp.sub(last0pos, last1pos); //tangent obtained using last two points.
tangentUnit = VecOp.unit(tangent);
@@ -106,9 +104,9 @@
displacement = VecOp.sub(clusterPosition, last1pos);
displacementUnit = VecOp.unit(displacement);
if(debug) {
- System.out.println("Debug: From last two points");
- System.out.println("Debug: 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.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()+")");
}
}
@@ -117,14 +115,15 @@
m_degree=m_angle*180/3.14;
m_trans = displacement.magnitude()*Math.sin(m_angle);
m_costheta = VecOp.dot(tangentUnit, displacementUnit);
+
m_histo_radian_distance.fill(m_angle,m_distance);
m_histo_degree_distance.fill(m_degree,m_distance);
m_histo_costheta_distance.fill(m_costheta,m_distance);
m_histo_trans_longi.fill(m_trans, m_distance);
if(debug){
- System.out.println("Debug: Angle= " + m_angle);
- System.out.println("Debug: Distance= " + m_distance);
+ System.out.print(" Angle= " + m_angle);
+ System.out.println(" Distance= " + m_distance);
}
if (m_angle < m_limit) return new Double(m_angle); // return the angle if it is below cone angle.
@@ -132,212 +131,4 @@
return null;
}
-
- //find last two points just before track started to shower.
- protected List<Hep3Vector> computePointsOfShower(Track tr) {
- Set<Cluster> clusoftrk = m_newMapTrackToShowerComponents.get(tr);
- Set<Cluster> mipsoftrk = new HashSet<Cluster>();
- for(Cluster cl : clusoftrk)
- if(m_mips.contains(cl))
- mipsoftrk.add(cl);
- if(mipsoftrk.size() == 0) return null; //no any mip found
-
- if(debug) System.out.println("Given track");
-
- Hep3Vector interceptPoint = m_findCluster.performExtrapolation(tr);
- Hep3Vector tangent = m_findCluster.getTangent();
- Hep3Vector tangentUnit = VecOp.unit(tangent);
-
- Hep3Vector first = new BasicHep3Vector();
- Hep3Vector trackThreeMomentum = new BasicHep3Vector(tr.getMomentum());
- double trackMomentum = trackThreeMomentum.magnitude();
-
- if(debug){
- System.out.println("Debug: Size of Cluster of track= " + clusoftrk.size());
- System.out.println("Debug: Size of Cluster of Mip of track= " + mipsoftrk.size());
- System.out.println("Debug: trk momentum= " + trackMomentum);
- if(interceptPoint != null) System.out.println("Debug: extra trk: pos= " + interceptPoint.magnitude() + " ("+ interceptPoint.x() + " " + interceptPoint.y() + " " + interceptPoint.z() + ")" );
- }
-
- //find closest mip to track
- SortedMap<Double, Cluster> mipsSortedByDistance = new TreeMap();
- for(Cluster cl : mipsoftrk){
- CalorimeterHit firstHit = cl.getCalorimeterHits().get(0);
- IDDecoder id = firstHit.getIDDecoder();
- id.setID(firstHit.getCellID());
- double x = id.getPosition()[0];
- double y = id.getPosition()[1];
- double z = id.getPosition()[2];
- first = new BasicHep3Vector(x, y, z);
- Hep3Vector displacement = new BasicHep3Vector();
- if (interceptPoint != null) {
- displacement = VecOp.sub(first, interceptPoint);
- double distance = displacement.magnitude();
- mipsSortedByDistance.put(distance, cl);
- }
- }
-
- //arrange by the length of the hit position
- SortedMap<Double, CalorimeterHit> sortedMIPbyPos = new TreeMap();
- Cluster cl = mipsSortedByDistance.get(mipsSortedByDistance.firstKey());
- for( int i = 0 ; i < cl.getCalorimeterHits().size() ; i++){
- CalorimeterHit hit=cl.getCalorimeterHits().get(i);
- IDDecoder idMIP = hit.getIDDecoder();
- idMIP.setID(hit.getCellID());
- double x = idMIP.getPosition()[0];
- double y = idMIP.getPosition()[1];
- double z = idMIP.getPosition()[2];
- double pmag = Math.sqrt(x*x+y*y+z*z);
- sortedMIPbyPos.put(pmag, hit);
- }
-
- List<Hep3Vector> miphit = new ArrayList<Hep3Vector>();//vector list of each hit
-
- Set<Long> allHitIDs = new HashSet<Long>();
- for (CalorimeterHit hit : m_allhits) allHitIDs.add(hit.getCellID());
-
- int lastlayer = 0;
- long[] lastlayer5x5x2 = new long[26];
- String lastDet = new String();
- int sizeofmip = sortedMIPbyPos.size();
-
- if(debug) System.out.println("Debug: Starting mip");
- Iterator itermip = sortedMIPbyPos.entrySet().iterator();
- for(int i=0 ; i < sizeofmip ;i++){
- Map.Entry entry = (Map.Entry) itermip.next();
- Double key = (Double)entry.getKey();
- CalorimeterHit hit = (CalorimeterHit)entry.getValue();
- IDDecoder id = hit.getIDDecoder();
- id.setID(hit.getCellID());
- if(i==0){
- if(id.getLayer()>3){
- System.out.println("No mip found near interceptpoint! It's secondary mip.");
- return null;
- }
- }
- Hep3Vector pos = hitPosition(hit,allHitIDs);
- miphit.add(pos);
- }
- if(debug) System.out.println("Debug: The end of the mip" );
-
- for( int i = 0 ; i < 100 ; i++){ // 100? until finding shower starting point
- //To get the information from last hit
- CalorimeterHit hit = sortedMIPbyPos.get(sortedMIPbyPos.lastKey());
- IDDecoder id = hit.getIDDecoder();
- id.setID(hit.getCellID());
- double x = id.getPosition()[0];
- double y = id.getPosition()[1];
- double z = id.getPosition()[2];
- double p = Math.sqrt(x*x+y*y+z*z);
- double r = Math.sqrt(x*x+y*y);
- lastlayer = id.getLayer();
- lastDet = hit.getSubdetector().getName();
- lastlayer5x5x2 = id.getNeighbourIDs(1, 2, 2);
-
- //Searh next hit in next layer and arrange in angular order
- SortedMap<Double, CalorimeterHit> sortedHitbyPos= new TreeMap();
- for (CalorimeterHit hitall : m_allhits){
- IDDecoder idall = hitall.getIDDecoder();
- idall.setID(hitall.getCellID());
- double xall = idall.getPosition()[0];
- double yall = idall.getPosition()[1];
- double zall = idall.getPosition()[2];
- double hitmag = Math.sqrt(xall*xall+yall*yall+zall*zall);
- Set<Long> nearIDs = new HashSet<Long>();
- for (long lastnearID : lastlayer5x5x2 ) nearIDs.add(lastnearID);
-
- //use geometric information and layer number to find next hit
- boolean EcalToHcal= lastlayer == 30 && hitmag > p && idall.getLayer() == 0;
- boolean InSide = hitmag > p && nearIDs.contains(hitall.getCellID());
- boolean EndToBarrelEcal = ((r < ecalRs && r > ecalRs-4) && hitmag > p && ((ecalZs < z) || (z < ecalZe)) && (idall.getLayer()==0));
- boolean EndToBarrelHcal = ((r < ecalRs && r > ecalRs-4) && hitmag > p && ((hcalZs < z) || (z < hcalZe)) && (idall.getLayer()==0));
-
-
- if(debug){
- if(EndToBarrelEcal || EndToBarrelHcal ) System.out.println("This is a boundary condition");
- }
-
- if(EcalToHcal || InSide || EndToBarrelEcal || EndToBarrelHcal){
- Hep3Vector hitpos = new BasicHep3Vector(xall,yall,zall);
- Hep3Vector test = VecOp.sub(hitpos, miphit.get(miphit.size()-1));
- Hep3Vector test_unit = VecOp.unit(test);
- Hep3Vector lastmip = VecOp.sub(miphit.get(miphit.size()-1), miphit.get(miphit.size()-2));
- Hep3Vector lastmip_unit = VecOp.unit(lastmip);
- double a = Math.acos(VecOp.dot(lastmip_unit, test_unit));
- sortedHitbyPos.put(a, hitall);
- }
-
- }
-
- if(sortedHitbyPos.size() == 0) continue;
- CalorimeterHit hitAdd = sortedHitbyPos.get(sortedHitbyPos.firstKey()); //add the hit with a smallest angle
- IDDecoder idAdd = hitAdd.getIDDecoder();
- idAdd.setID(hitAdd.getCellID());
- double xAdd = idAdd.getPosition()[0];
- double yAdd = idAdd.getPosition()[1];
- double zAdd = idAdd.getPosition()[2];
- double Addmag=Math.sqrt(xAdd*xAdd+yAdd*yAdd+zAdd*zAdd);
- sortedMIPbyPos.put(Addmag,hitAdd);
-
- Hep3Vector tempPos = hitPosition(hitAdd,allHitIDs);
- miphit.add(tempPos);
-
- if (exam > 2) {
- if(debug) System.out.println("Debug: showering is already started " + exam + " layers ago" );
- if(miphit.size() < 5) break; // if size of miphit is fewer than 5, not need to go back to three layers before
- else for(int k=0; k < 3 ; k++) miphit.remove(miphit.size()-1);
- break;
- }
- }
-
- if(debug){
- int k=0;
- if(miphit.size() < 4) k = miphit.size();
- else k = 4;
- for(int i =0; i < k ;i++){
- System.out.println("Debug: last " + i + " pos= " + miphit.get(miphit.size()-1-i).magnitude() + " (" +miphit.get(miphit.size()-1-i).x() + " " + miphit.get(miphit.size()-1-i).y() + " " + miphit.get(miphit.size()-1-i).z()+ ")" );
- }
- }
-
- List<Hep3Vector> lastTwoPositions = new ArrayList<Hep3Vector>();
- if(miphit.size() < 5){ //since there are only 4 hits
- lastTwoPositions.add(miphit.get(1));
- lastTwoPositions.add(miphit.get(0));
- } else {
- lastTwoPositions.add(miphit.get(miphit.size()-1));
- lastTwoPositions.add(miphit.get(miphit.size()-2));
- }
-
- return lastTwoPositions;
- }
-
- //give the hit position vector and examine neighbour hits
- protected Hep3Vector hitPosition(CalorimeterHit hit, Set<Long> allHitIDs){
-
- IDDecoder id = hit.getIDDecoder();
- id.setID(hit.getCellID());
- int lay = id.getLayer();
- double x = id.getPosition()[0];
- double y = id.getPosition()[1];
- double z = id.getPosition()[2];
- Hep3Vector v = new BasicHep3Vector(x,y,z);
- int p = (int) (v.magnitude());
-
- int count[] = new int[4];
- for (int i =0 ; i < 4 ; i++){
- long[] neighbours = id.getNeighbourIDs(0,i+1,i+1);
- for(long neighbourID : neighbours)
- if (allHitIDs.contains(neighbourID)) count[i]++;
- }
-
- if ( count[3] >=2 ) exam++ ; // if number of neighbour hits in 9x9 is larger than 1, count 1
- else exam = 0; // if number of neighbour is none, set to zero and start again from 0
-
- if(debug){
- System.out.println(" hit: pos=" + p + " , " + count[0] + " hits in 3x3, " + count[1] + " hits in 5x5 " + count[2] + " hits in 7x7 " + count[3] + " hits in 9x9 " + lay + " layer of " + hit.getSubdetector().getName());
- }
-
- return v;
- }
-
}
