lcsim/src/org/lcsim/contrib/uiowa
diff -N ShowerPointFinder.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ShowerPointFinder.java 10 Jul 2008 20:34:52 -0000 1.1
@@ -0,0 +1,290 @@
+package org.lcsim.contrib.uiowa;
+
+import java.util.*;
+import org.lcsim.util.*;
+import org.lcsim.event.util.*;
+import org.lcsim.event.*;
+import hep.physics.vec.*;
+
+import java.io.IOException;
+import hep.physics.particle.properties.*;
+import org.lcsim.util.*;
+import org.lcsim.util.hitmap.*;
+import org.lcsim.recon.cluster.mipfinder.*;
+import org.lcsim.recon.cluster.util.*;
+import org.lcsim.recon.cluster.clumpfinder.*;
+import org.lcsim.recon.pfa.identifier.*;
+import org.lcsim.mc.fast.tracking.ReconTrack;
+import org.lcsim.event.base.*;
+import hep.physics.particle.Particle;
+import org.lcsim.recon.cluster.structural.likelihood.*;
+import org.lcsim.recon.cluster.structural.*;
+import org.lcsim.recon.cluster.mipfinder.*;
+import org.lcsim.recon.cluster.clumpfinder.*;
+import org.lcsim.util.swim.Line;
+import org.lcsim.geometry.Detector;
+import org.lcsim.geometry.subdetector.CylindricalCalorimeter;
+import org.lcsim.geometry.*;
+import org.lcsim.util.decision.*;
+
+ public class ShowerPointFinder{
+ 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;
+ protected double m_angle = 0;
+ protected double m_degree = 0;
+ protected double m_costheta = 0;
+ protected double m_trans = 0;
+
+ //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 EventHeader m_event;
+
+ public ShowerPointFinder(LocalHelixExtrapolator findCluster, Map<Track, Set<Cluster>> newMapTrackToShowerComponents, List<Cluster> mips, Set<CalorimeterHit> allHits) {
+ m_findCluster = findCluster;
+ m_newMapTrackToShowerComponents = newMapTrackToShowerComponents;
+ m_mips = mips;
+ m_allhits = allHits;
+ }
+
+
+ public Map<Track, BasicCluster> findMips(List<Track> tracksSortedByMomentum) {
+ Map<Track, BasicCluster> MapTrkToMIP = new HashMap<Track, BasicCluster>();
+ for(Track tr : tracksSortedByMomentum){
+ 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) continue; //no any mip found
+
+ if(debug) System.out.println("Given track");
+
+ BasicCluster mipclus = new BasicCluster();
+
+ 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());
+
+ 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();
+ Hep3Vector pos = hitPosition(hit,allHitIDs);
+ miphit.add(pos);
+ mipclus.addHit(hit);
+ }
+ 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
+ double p=0;
+ double r=0;
+ int lastlayer = 0;
+ double x=0;
+ double y=0;
+ double z=0;
+ Hep3Vector last0 = new BasicHep3Vector();
+ Hep3Vector last1 = new BasicHep3Vector();
+
+ for(int s=0; s<2 ; s++){
+ CalorimeterHit hit = mipclus.getCalorimeterHits().get( mipclus.getCalorimeterHits().size()-1-s);
+ IDDecoder id = hit.getIDDecoder();
+ id.setID(hit.getCellID());
+ x = id.getPosition()[0];
+ y = id.getPosition()[1];
+ z = id.getPosition()[2];
+ if(s==0){
+ p = Math.sqrt(x*x+y*y+z*z);
+ r = Math.sqrt(x*x+y*y);
+ lastlayer = id.getLayer();
+ lastDet = hit.getSubdetector().getName();
+ lastlayer5x5x2 = id.getNeighbourIDs(1, 2, 2);
+ last0 = new BasicHep3Vector(x,y,z);
+ } else {last1 = new BasicHep3Vector(x,y,z);}
+ }
+
+
+ //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);
+ Hep3Vector hitpos = new BasicHep3Vector(xall,yall,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 test = VecOp.sub(hitpos, last0);
+ Hep3Vector test_unit = VecOp.unit(test);
+ Hep3Vector lastmip = VecOp.sub(last0, last1);
+ 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);
+ mipclus.addHit(hitAdd);
+
+ if (exam > 2) {
+ if(debug) System.out.println("Debug: showering is already started " + exam + " layers ago" );
+ if(mipclus.getCalorimeterHits().size() < 5) break; // if hit size is fewer than 5, don't remove
+ else for(int k=0; k < 3 ; k++) {
+ miphit.remove(miphit.size()-1);
+ CalorimeterHit toRemove = mipclus.getCalorimeterHits().get(mipclus.getCalorimeterHits().size()-1);
+ mipclus.removeHit(toRemove);
+ }
+ break;
+ }
+ }
+
+ if(debug){
+ int k=0;
+ if(mipclus.getCalorimeterHits().size() < 4) k = mipclus.getCalorimeterHits().size();
+ else k = 4;
+ for(int i =0; i < k ;i++){
+ CalorimeterHit hit = mipclus.getCalorimeterHits().get( mipclus.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);
+ System.out.println("Debug: last "+i+" pos= "+v.magnitude()+" (" +v.x() + " " + v.y() + " " + v.z()+ ")" );
+ }
+ }
+
+ if(debug) System.out.println("MIP hit size= " + mipclus.getCalorimeterHits().size());
+
+ CalorimeterHit hit = mipclus.getCalorimeterHits().get(0);//first fit
+ IDDecoder id = hit.getIDDecoder();
+ id.setID(hit.getCellID());
+ if(id.getLayer() <= 3) MapTrkToMIP.put(tr,mipclus); //remove possible secondary mip.
+ }
+ return MapTrkToMIP;
+ }
+
+ //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;
+ }
+
+}