diff -N HighPElectronFinder.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ HighPElectronFinder.java	4 Dec 2008 20:11:58 -0000	1.1
@@ -0,0 +1,170 @@

+/*
+ * HighPElectronFinder.java
+ *
+ * Created on December 4, 2008, 8:15 AM
+ *
+ * Find High momentum electrons using PFA output and modify
+ * ReconstructedParticles
+ */
+
+package org.lcsim.recon.postrecon.leptonID.electron;
+import hep.physics.particle.properties.ParticlePropertyManager;
+import hep.physics.particle.properties.ParticlePropertyProvider;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.BasicHepLorentzVector;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+import java.util.Map;
+import java.util.HashMap;
+import java.io.IOException;
+import java.util.List;
+import java.util.ArrayList;
+import org.lcsim.geometry.Detector;
+import org.lcsim.geometry.compact.Subdetector;
+import org.lcsim.geometry.subdetector.CylindricalCalorimeter;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.Track;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.base.BaseReconstructedParticle;
+import org.lcsim.event.base.BaseParticleID;
+import org.lcsim.util.Driver;
+import org.lcsim.util.loop.LCIODriver;
+import org.lcsim.util.swim.HelixSwimmer;
+import static java.lang.Math.abs;
+import static java.lang.Math.sqrt;
+import static java.lang.Math.acos;
+import static java.lang.Math.pow;
+
+/**
+ *
+ * @author cassell
+ * algorithm from Tim Barklow
+ */
+public class HighPElectronFinder extends Driver {
+    
+    final double trackMomMagMin=15.;
+    final double cosangMaxMin=0.99998;
+    final double[] shapeCut = {300.,2000.,2000.};
+    final double epDiffCut = 0.3;
+    final double mElectron;
+    boolean modifyRP = true;
+    int ncnt=0;
+    int ncnt_in=0;
+    int eventNumber;
+    Hep3Vector nucr=new BasicHep3Vector(0,0,0);
+    double m_ECAL_barrel_r;
+    double m_ECAL_endcap_z;
+    final ParticlePropertyProvider dPPP;
+    Map<Track,Cluster> electronTrackClusterMap;
+    /** Creates a new instance of HighPElectronFinder */
+    public HighPElectronFinder() {
+        dPPP = ParticlePropertyManager.getParticlePropertyProvider();
+        mElectron = dPPP.get(11).getMass();
+        System.out.println(" HighPElectronFinder constructor   mElectron= "+mElectron);
+    }
+    protected void process(EventHeader event) {
+        eventNumber=event.getEventNumber();
+        ncnt++;
+//	    if(ncnt_in%100 == 1) System.out.println(" ncnt_in= "+ncnt_in+" ncnt= "+ncnt+" eventNumber= "+eventNumber);
+        List<Track> tracks = new ArrayList<Track>(event.get(Track.class,"Tracks"));
+        List<Cluster> clusters = new ArrayList<Cluster>(event.get(Cluster.class,"Clusters"));
+        Detector thisDetector = event.getDetector();
+        double thisBfield=thisDetector.getFieldMap().getField(new BasicHep3Vector(0,0,0)).magnitude();
+        CylindricalCalorimeter emb = ((CylindricalCalorimeter) thisDetector.getSubdetectors().get("EMBarrel"));
+        CylindricalCalorimeter eme = ((CylindricalCalorimeter) thisDetector.getSubdetectors().get("EMEndcap"));
+        m_ECAL_barrel_r = emb.getLayering().getDistanceToLayerSensorMid(0);
+        m_ECAL_endcap_z = eme.getLayering().getDistanceToLayerSensorMid(0);
+        HelixSwimmer helixswim=new HelixSwimmer(thisBfield);
+        List<ReconstructedParticle> recoParticles = event.get(ReconstructedParticle.class,"ReconstructedParticles");
+        electronTrackClusterMap = new HashMap<Track,Cluster>();
+        for(ReconstructedParticle p:recoParticles)
+        {
+            if(Math.abs(p.getMass()-mElectron) < .00001)
+            {
+                tracks.remove(p.getTracks().get(0));
+                if(p.getClusters().size() > 0)clusters.remove(p.getClusters().get(0));
+            }
+        }
+        for(Track track : tracks) {
+            Hep3Vector trackMom = new BasicHep3Vector(track.getMomentum());
+            double[] trackPar=track.getTrackParameters();
+            helixswim.setTrack(track);
+            double distecalbarr=helixswim.getDistanceToRadius(m_ECAL_barrel_r);
+            double distecalendcap=helixswim.getDistanceToZ(m_ECAL_endcap_z);
+            if(trackMom.magnitude() > trackMomMagMin && ((!Double.isNaN(distecalbarr)) || (!Double.isNaN(distecalendcap)))) {
+                int nHitsElectron = 0;
+                if((!Double.isNaN(distecalbarr)) && (Double.isNaN(distecalendcap))) {
+                    nucr=helixswim.getPointAtLength(distecalbarr);
+                } else if((Double.isNaN(distecalbarr)) && (!Double.isNaN(distecalendcap))) {
+                    nucr=helixswim.getPointAtLength(distecalendcap);
+                } else if(helixswim.getPointAtLength(distecalbarr).magnitude() < helixswim.getPointAtLength(distecalendcap).magnitude()) {
+                    nucr=helixswim.getPointAtLength(distecalbarr);
+                } else {
+                    nucr=helixswim.getPointAtLength(distecalendcap);
+                }
+                Hep3Vector normTrackMom=VecOp.mult(1./trackMom.magnitude(),trackMom);
+                Hep3Vector normEcal=VecOp.mult(1./nucr.magnitude(),nucr);
+                double cosangMax=-1000.;
+                int iclus=-1;
+                Hep3Vector vecClus=new BasicHep3Vector(0,0,0);
+                double emClusEOpt=0.;
+                for (Cluster cluster : clusters) {
+                    double[] subE = cluster.getSubdetectorEnergies();
+                    double eSubTot = 0.;
+                    for(int i=0;i<subE.length;i++) eSubTot += subE[i];
+                    double emClusE = subE[2]+subE[6];
+                    Hep3Vector vec1=new BasicHep3Vector(cluster.getPosition());
+                    double cosang = VecOp.dot(normEcal,vec1)/vec1.magnitude();
+                    double cosangTrackMom = VecOp.dot(normTrackMom,vec1)/vec1.magnitude();
+                    if(cosang > cosangMax) {
+                        cosangMax=cosang;
+                        iclus=clusters.indexOf(cluster);
+                        vecClus=vec1;
+                        emClusEOpt=emClusE;
+                    }
+                }
+                if(iclus < 0)continue;
+                Cluster cluster = clusters.get(iclus);
+                double[] clusShape = cluster.getShape();
+                double[] subE = cluster.getSubdetectorEnergies();
+                if(cosangMax > cosangMaxMin) {
+                    if(abs((emClusEOpt-trackMom.magnitude())/trackMom.magnitude()) < epDiffCut) {
+                        if(clusShape[0]<shapeCut[0] && clusShape[1]<shapeCut[1] && clusShape[2]<shapeCut[2] ) {
+                            
+                            electronTrackClusterMap.put(track,cluster);
+                            clusters.remove(cluster);
+                        }
+                    }
+                }
+            }
+        }
+        
+        if(modifyRP) {
+            for(Track track : electronTrackClusterMap.keySet()) {
+                Hep3Vector trackMom = new BasicHep3Vector(track.getMomentum());
+                Cluster cluster = electronTrackClusterMap.get(track);
+                List<ReconstructedParticle> removeRP = new ArrayList<ReconstructedParticle>();
+                for(ReconstructedParticle recoParticle : recoParticles) {
+                    if(recoParticle.getTracks().indexOf(track) >= 0) {
+                        removeRP.add(recoParticle);
+                    }
+                    if(recoParticle.getClusters().indexOf(cluster) >= 0) {
+                        removeRP.add(recoParticle);
+                    }
+                }
+                for(ReconstructedParticle recoParticle : removeRP) recoParticles.remove(recoParticle);
+                BaseReconstructedParticle electronRP = new BaseReconstructedParticle(mElectron);
+                electronRP.set4Vector(new BasicHepLorentzVector(sqrt(pow(mElectron,2)+trackMom.magnitudeSquared()),trackMom));
+                electronRP.setCharge(track.getCharge());
+                electronRP.setParticleIdUsed(new BaseParticleID(dPPP.get(-11*track.getCharge())));
+                electronRP.setReferencePoint(new BasicHep3Vector(track.getReferencePoint()));
+                electronRP.addTrack(track);
+                electronRP.addCluster(cluster);
+                recoParticles.add(electronRP);
+            }
+        }
+        
+    }
+    
+}