lcsim/src/org/lcsim/contrib/Barklow
diff -N HighPElectronID.java
--- HighPElectronID.java 26 Nov 2008 12:53:46 -0000 1.2
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,171 +0,0 @@
-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);
- }
- }
-
- }
-
-}