lcsim/src/org/lcsim/contrib/Barklow
diff -u -r1.1 -r1.2
--- HighPElectronID.java 26 Nov 2008 12:35:40 -0000 1.1
+++ HighPElectronID.java 26 Nov 2008 12:53:46 -0000 1.2
@@ -1,170 +1,171 @@
-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 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;
-/*
- * Post-PFA identification of high momentum electrons
- * Track and cluster for electron candidates zre stored in Map<Track,Cluster> electronTrackClusterMap
- * (one track and one cluster per electron)
- * Rebuilds ReconstructedParticle collection if modifyRP = true.
- *
- */
-public class HighPElectronID extends Driver
-{
- final int nEventSkip=0;
- final double trackMomMagMin=15.;
- final double cosangMaxMin=0.99998;
- 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;
- public HighPElectronID() throws IOException
- {
- if(modifyRP) {
- LCIODriver writer = new LCIODriver("after_electronID.slcio");
- add(writer);
- }
- dPPP = ParticlePropertyManager.getParticlePropertyProvider();
- mElectron = dPPP.get(11).getMass();
- System.out.println(" HighPElectronID constructor mElectron= "+mElectron);
- }
- protected void process(EventHeader event)
- {
- ncnt_in++;
- if (ncnt_in <= nEventSkip) {
- throw new Driver.NextEventException();
- }
- else {
- eventNumber=event.getEventNumber();
- ncnt++;
- if(ncnt_in%100 == 1) System.out.println(" ncnt_in= "+ncnt_in+" ncnt= "+ncnt+" eventNumber= "+eventNumber);
- List<Track> tracks = 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(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=0;
- 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;
- }
- }
- Cluster cluster = clusters.get(iclus);
- double[] clusShape = cluster.getShape();
- double[] subE = cluster.getSubdetectorEnergies();
- if(cosangMax > cosangMaxMin) {
- if(abs((emClusEOpt-trackMom.magnitude())/trackMom.magnitude()) < 0.3) {
- if(clusShape[0]<300. && clusShape[1]<2000. && clusShape[2]<2000. ) {
-
- 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);
- }
- super.process(event);
- }
- }
-
- }
-
-}
+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;
+/*
+ * Post-PFA identification of high momentum electrons
+ * Track and cluster for electron candidates zre stored in Map<Track,Cluster> electronTrackClusterMap
+ * (one track and one cluster per electron)
+ * Rebuilds ReconstructedParticle collection if modifyRP = true.
+ *
+ */
+public class HighPElectronID extends Driver
+{
+ final int nEventSkip=0;
+ final double trackMomMagMin=15.;
+ final double cosangMaxMin=0.99998;
+ 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;
+ public HighPElectronID() throws IOException
+ {
+ if(modifyRP) {
+ LCIODriver writer = new LCIODriver("after_electronID.slcio");
+ add(writer);
+ }
+ dPPP = ParticlePropertyManager.getParticlePropertyProvider();
+ mElectron = dPPP.get(11).getMass();
+ System.out.println(" HighPElectronID constructor mElectron= "+mElectron);
+ }
+ protected void process(EventHeader event)
+ {
+ ncnt_in++;
+ if (ncnt_in <= nEventSkip) {
+ throw new Driver.NextEventException();
+ }
+ else {
+ eventNumber=event.getEventNumber();
+ ncnt++;
+ if(ncnt_in%100 == 1) System.out.println(" ncnt_in= "+ncnt_in+" ncnt= "+ncnt+" eventNumber= "+eventNumber);
+ List<Track> tracks = 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(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=0;
+ 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;
+ }
+ }
+ Cluster cluster = clusters.get(iclus);
+ double[] clusShape = cluster.getShape();
+ double[] subE = cluster.getSubdetectorEnergies();
+ if(cosangMax > cosangMaxMin) {
+ if(abs((emClusEOpt-trackMom.magnitude())/trackMom.magnitude()) < 0.3) {
+ if(clusShape[0]<300. && clusShape[1]<2000. && clusShape[2]<2000. ) {
+
+ 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);
+ }
+ super.process(event);
+ }
+ }
+
+ }
+
+}
