lcsim/src/org/lcsim/recon/cluster/mipfinder
diff -N MipQualityDecision.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ MipQualityDecision.java 13 Jul 2008 23:31:31 -0000 1.1
@@ -0,0 +1,129 @@
+package org.lcsim.recon.cluster.mipfinder;
+
+import java.util.*;
+
+import org.lcsim.util.*;
+import org.lcsim.util.decision.*;
+import org.lcsim.geometry.subdetector.CylindricalCalorimeter;
+import org.lcsim.geometry.*;
+import org.lcsim.event.*;
+
+/**
+ * A class to look at whether a calorimeter cluster is
+ * consistent with being a track/MIP. Criteria are based on
+ *
+ * a) number of layers in cluster
+ * b) number of layers with >1 hit
+ * c) chi^2 of helix fit to cluster
+ *
+ * @author [log in to unmask]
+ * @version $Id: MipQualityDecision.java,v 1.1 2008/07/13 23:31:31 mcharles Exp $
+ */
+
+
+public class MipQualityDecision implements DecisionMakerSingle<Cluster>
+{
+ public MipQualityDecision() {
+ // No user-settable options for now
+ }
+
+ public boolean valid(Cluster mip) {
+ // How many layers in this MIP?
+ // And how many instances of >1 hit in a layer?
+ int countDuplicates = 0;
+ Set<Integer> layersSeen = new HashSet<Integer>();
+ for (CalorimeterHit hit : mip.getCalorimeterHits()) {
+ int layer = getLayer(hit);
+ Integer layerInt = new Integer(layer);
+ if (!layersSeen.contains(layerInt)) {
+ layersSeen.add(layerInt);
+ } else {
+ countDuplicates++;
+ }
+ }
+ int countHits = mip.getCalorimeterHits().size();
+ int countLayers = layersSeen.size();
+ double fractionDuplicates = ((double)(countDuplicates)) / ((double)(countLayers));
+
+ // Now check chisq. Start by getting the hits:
+ List<org.lcsim.fit.helicaltrack.HelicalTrackHit> hitsToFit = new Vector<org.lcsim.fit.helicaltrack.HelicalTrackHit>();
+ for (CalorimeterHit calHit : mip.getCalorimeterHits()) {
+ double pos[] = calHit.getPosition();
+ String subdetName = calHit.getSubdetector().getName();
+ double err_rphi = 0.0;
+ double err_z = 0.0;
+ if (subdetName.compareTo("EMBarrel")==0) {
+ // Barrel -- r fixed but phi uncertain
+ // rphi error = cell size / sqrt(12)
+ // z error = cell size / sqrt(12)
+ err_rphi = 3.5 / 3.464;
+ err_z = 3.5 / 3.464;
+ } else if (subdetName.compareTo("EMEndcap")==0) {
+ // Endcap -- z fixed but r and phi uncertain
+ // rphi error = cell size / sqrt(12) ??
+ // z error = silicon width / sqrt(12)
+ err_rphi = 3.5 / 3.464;
+ err_z = 0.32 / 3.464; // 320 micron silicon
+ } else if (subdetName.compareTo("HADBarrel")==0) {
+ err_rphi = 10.0 / 3.464;
+ err_z = 10.0 / 3.464;
+ } else if (subdetName.compareTo("HADEndcap")==0) {
+ err_rphi = 10.0 / 3.464;
+ err_z = 3.0 / 3.464; // 1.2 mm gas OR 5.0 mm scint... pick 3mm as a compromise
+ } else {
+ throw new AssertionError("Subdet not recognized: "+subdetName);
+ }
+ org.lcsim.fit.helicaltrack.HelicalTrack3DHit fitHit = new org.lcsim.fit.helicaltrack.HelicalTrack3DHit(pos[0], pos[1], pos[2], err_rphi, err_z);
+ hitsToFit.add(fitHit);
+ }
+
+ // Now do the fit:
+ org.lcsim.fit.helicaltrack.HelicalTrackFitter fitter = new org.lcsim.fit.helicaltrack.HelicalTrackFitter();
+ org.lcsim.fit.helicaltrack.HelicalTrackFitter.FitStatus fitOK = fitter.fit(hitsToFit);
+
+ // Check fit quality:
+ if (fitOK == org.lcsim.fit.helicaltrack.HelicalTrackFitter.FitStatus.Success) {
+ org.lcsim.fit.helicaltrack.HelicalTrackFit fitResult = fitter.getFit();
+ double chisqArray[] = fitResult.chisq();
+ int ndfArray[] = fitResult.ndf();
+ if (chisqArray.length != ndfArray.length) {
+ throw new AssertionError("Chisq length = "+chisqArray.length+" but ndf length = "+ndfArray.length);
+ }
+ if (ndfArray.length >= 1) {
+ double chisq = chisqArray[0];
+ double ndf = ndfArray[0];
+ if (ndf > 0) {
+ double prob = 1.0;
+ if (chisq > 0.0) {
+ prob = org.lcsim.math.chisq.ChisqProb.gammq(ndf,chisq);
+ }
+ double quality = prob * (1.0 - fractionDuplicates);
+ double cut = 0.05; // require 5% prob by default
+ if (countLayers < 6) {
+ cut = 0.1; // More stringent for short MIPs which are more likely to be combinatorics
+ } else if (countLayers > 15) {
+ cut = 0.01; // Less stringent for long MIPs which are likely to be sound but may be multiple-scattered
+ }
+ boolean qualityPass = (quality >= cut);
+ return qualityPass;
+ } else {
+ // Not enough points for chisq check
+ return false;
+ }
+ } else {
+ // No fit output
+ return false;
+ }
+ } else {
+ // Fit failed
+ return false;
+ }
+ }
+
+ protected int getLayer(CalorimeterHit hit) {
+ org.lcsim.geometry.IDDecoder id = hit.getIDDecoder();
+ id.setID(hit.getCellID());
+ int layer = id.getLayer();
+ return layer;
+ }
+}
lcsim/src/org/lcsim/recon/cluster/mipfinder
diff -N ShowerPointFinder.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ShowerPointFinder.java 13 Jul 2008 23:31:31 -0000 1.1
@@ -0,0 +1,268 @@
+package org.lcsim.recon.cluster.mipfinder;
+
+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; //to test neighbour
+ protected boolean debug = false;
+
+ //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
+
+ 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>();
+ BasicCluster mipclus = new 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) {//No any mip found
+ MapTrkToMIP.put(tr,mipclus);
+ continue; //go to next track
+ }
+
+ //new cluster for finding shower point
+ Hep3Vector interceptPoint = m_findCluster.performExtrapolation(tr); //not need except debuging
+ //Hep3Vector tangent = m_findCluster.getTangent();
+ //Hep3Vector tangentUnit = VecOp.unit(tangent);
+
+ if(debug){
+ Hep3Vector v = new BasicHep3Vector(tr.getMomentum());
+ System.out.println("Given track");
+ 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= " + v.magnitude());
+ 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);
+ Hep3Vector first = hitPosition(firstHit);
+ double distance = first.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);
+ Hep3Vector v = hitPosition(hit);
+ double pmag = v.magnitude();
+ sortedMIPbyPos.put(pmag, hit);
+ }
+
+ //test if the mip is enough closet to track otherwise take it as no mip.
+ CalorimeterHit firstmip = sortedMIPbyPos.get(sortedMIPbyPos.firstKey());
+ IDDecoder firstmipid = firstmip.getIDDecoder();
+ firstmipid.setID(firstmip.getCellID());
+ if(firstmipid.getLayer() > 3) { // if the mip is not in first 3 layers, put 0 in mipclus
+ MapTrkToMIP.put(tr,mipclus);
+ continue; //go to next track
+ }
+
+ if(debug) System.out.println("Debug: Starting mip");
+ Iterator itermip = sortedMIPbyPos.entrySet().iterator();
+ for(int i=0 ; i < sortedMIPbyPos.size() ;i++){
+ Map.Entry entry = (Map.Entry) itermip.next();
+ Double key = (Double)entry.getKey();
+ CalorimeterHit hit = (CalorimeterHit)entry.getValue();
+ examPosition(hit);
+ mipclus.addHit(hit);
+ }
+ if(debug) System.out.println("Debug: The end of the mip" );
+
+ int lastIterationWithFoundHit = -1;
+ int maxSkippedLayers = 3;
+ for( int iIteration = 0 ; iIteration < 70 ; iIteration++){ // 70? until finding shower starting point
+ //To get the information from last hit
+ CalorimeterHit hit = mipclus.getCalorimeterHits().get( mipclus.getCalorimeterHits().size()-1);
+ 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);
+ double lastlayer = id.getLayer();
+
+ // Check for neighbouring hits in a 5x5x(2n+1) grid
+ // where n is the number of iterations since we last saw a hit.
+ // For example, if we saw a hit in the previous layer/iteration
+ // it'll be a 5x5x3 grid.
+ Set<Long> nearIDs = new HashSet<Long>();
+ int countIterationsSinceLastHit = (iIteration - lastIterationWithFoundHit);
+ if (countIterationsSinceLastHit <= 0) { throw new AssertionError("Book-keeping error"); }
+ if (countIterationsSinceLastHit > maxSkippedLayers+1) { throw new AssertionError("Book-keeping error"); }
+ long[] nearbyHitArray = id.getNeighbourIDs(countIterationsSinceLastHit, 2, 2);
+ for (long lastnearID : nearbyHitArray ) { nearIDs.add(lastnearID); }
+
+ Hep3Vector last0 = hitPosition(hit);
+ CalorimeterHit hit1 = mipclus.getCalorimeterHits().get( mipclus.getCalorimeterHits().size()-2);
+ Hep3Vector last1 = hitPosition(hit1);
+
+ //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);
+ int hitLayer = idall.getLayer();
+ //use geometric information and layer number to find next hit
+ boolean EcalToHcal= lastlayer == 30 && hitmag > p && idall.getLayer() == 0;
+ boolean InSide = (hitmag > p || hitLayer > lastlayer) && nearIDs.contains(hitall.getCellID());
+ boolean nearInnerEdgeBarrelEcal = (r < ecalRs && r > ecalRs-4);
+ boolean nearOuterEdgeEndcapEcal = (Math.abs(z) < ecalZe && Math.abs(z) > ecalZs);
+ boolean nearInnerEdgeBarrelHcal = (r < hcalRs && r > hcalRs-4);
+ boolean nearOuterEdgeEndcapHcal = (Math.abs(z) < hcalZe && Math.abs(z) > hcalZs);
+ boolean EndToBarrelEcal = nearInnerEdgeBarrelEcal && nearOuterEdgeEndcapEcal && (hitmag > p) && (idall.getLayer()==0);
+ boolean EndToBarrelHcal = nearInnerEdgeBarrelHcal && nearOuterEdgeEndcapHcal && (hitmag > p) && (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) {
+ // No hit found in this layer.
+ if (countIterationsSinceLastHit <= maxSkippedLayers) {
+ continue; // go to next layer
+ } else {
+ break; // stop -- ran out of hits and haven't seen any for several iterations
+ }
+ } else {
+ // Found one or more hits
+ lastIterationWithFoundHit = iIteration;
+ }
+ CalorimeterHit hitAdd = sortedHitbyPos.get(sortedHitbyPos.firstKey()); //add the hit with a smallest angle
+ examPosition(hitAdd);
+ 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++) {
+ CalorimeterHit toRemove = mipclus.getCalorimeterHits().get(mipclus.getCalorimeterHits().size()-1);
+ mipclus.removeHit(toRemove);
+ }
+ break;
+ }
+ } //layer loop after mip
+
+ 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);
+ Hep3Vector v = hitPosition(hit);
+ System.out.println("Debug: last "+i+" pos= "+v.magnitude()+" (" +v.x() + " " + v.y() + " " + v.z()+ ")" );
+ }
+ System.out.println("MIP hit size= " + mipclus.getCalorimeterHits().size());
+ }
+
+ MapTrkToMIP.put(tr,mipclus);
+ } //track loop
+
+ return MapTrkToMIP;
+ }
+
+ //examine neighbour hits
+ protected void examPosition(CalorimeterHit hit){
+ Set<Long> allHitIDs = new HashSet<Long>();
+ for (CalorimeterHit allhit : m_allhits) allHitIDs.add(allhit.getCellID());
+
+ IDDecoder id = hit.getIDDecoder();
+ id.setID(hit.getCellID());
+ 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){
+ Hep3Vector v = hitPosition(hit);
+ int p = (int) (v.magnitude());
+ 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 " + id.getLayer() + " layer of " + hit.getSubdetector().getName());
+ }
+ }
+
+ //give the vector of hit
+ protected Hep3Vector hitPosition(CalorimeterHit hit){
+ 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);
+ return v;
+ }
+}