15 removed files
lcsim/src/org/lcsim/contrib/JanStrube/fastMC
diff -N IDResolutionTables.java
--- IDResolutionTables.java 29 Nov 2007 02:25:56 -0000 1.1
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,66 +0,0 @@
-/*
- * IDResolutionTables.java
- *
- * Created on July 6, 2005, 7:13 PM
- *
- * @version "$Id: IDResolutionTables.java,v 1.1 2007/11/29 02:25:56 jstrube Exp $"
- */
-
-package org.lcsim.contrib.JanStrube.fastMC;
-
-import org.lcsim.conditions.ConditionsSet;
-
-/**
- *
- * @author Daniel
- */
-public class IDResolutionTables {
-
- private double ElectronEff;
- private double MuonEff;
- private double ProtonEff;
- private double KaonEff;
- private double NeutronEff;
- private double WtChgTrkCal;
-
- /** Creates a new instance of IDResolutionTables */
- IDResolutionTables(ConditionsSet set)
- {
- ElectronEff = set.getDouble("Electron");
- MuonEff = set.getDouble("Muon");
- ProtonEff = set.getDouble("Proton");
- KaonEff = set.getDouble("Kaon");
- NeutronEff = set.getDouble("Neutron");
- WtChgTrkCal = set.getDouble("wt_charged_track_calorimeter_energy");
- }
-
- public double getElectronEff()
- {
- return ElectronEff;
- }
-
- public double getMuonEff()
- {
- return MuonEff;
- }
-
- public double getProtonEff()
- {
- return ProtonEff;
- }
-
- public double getKaonEff()
- {
- return KaonEff;
- }
-
- public double getNeutronEff()
- {
- return NeutronEff;
- }
-
- public double getWtChgTrkCal()
- {
- return WtChgTrkCal;
- }
-}
lcsim/src/org/lcsim/contrib/JanStrube/fastMC
diff -N MCFast.java
--- MCFast.java 29 Nov 2007 21:29:36 -0000 1.2
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,58 +0,0 @@
-package org.lcsim.contrib.JanStrube.fastMC;
-
-import org.lcsim.mc.fast.cluster.ronan.MCFastRonan;
-import org.lcsim.mc.fast.tracking.fix.FastMCTrackDriver;
-
-import java.util.logging.Handler;
-import java.util.logging.Level;
-import java.util.logging.Logger;
-import org.lcsim.util.Driver;
-
-/**
- *
- * @author Tony Johnson
- */
-public class MCFast extends Driver
-{
- /** Creates a new instance of MCFast */
-
- public static Logger log;
-
-
- public MCFast(boolean beamSpotConstraint, boolean simple, long seed, boolean printinfo)
- {
- this(beamSpotConstraint, simple, printinfo);
- getRandom().setSeed(seed);
- }
-
- public MCFast(boolean beamSpotConstraint, boolean simple, boolean printinfo)
- {
- log=getLogger();
- if(printinfo)
- {
- log.setLevel(Level.INFO);
- }
- else
- {
- log.setLevel(Level.WARNING);
- }
- add(new FastMCTrackDriver(beamSpotConstraint));
- add(new MCFastRonan());
- add(new MCFastReconstructedParticleDriver());
- }
-
- public MCFast(boolean beamSpotConstraint, boolean simple)
- {
- this(beamSpotConstraint, simple, false);
- }
-
- public MCFast(boolean beamSpotConstraint, boolean simple, long seed)
- {
- this(beamSpotConstraint, simple, seed, false);
- }
-
- public MCFast()
- {
- this(false,false,false);
- }
-}
lcsim/src/org/lcsim/contrib/JanStrube/fastMC
diff -N MCFastParticleID.java
--- MCFastParticleID.java 29 Nov 2007 02:25:56 -0000 1.1
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,55 +0,0 @@
-package org.lcsim.contrib.JanStrube.fastMC;
-
-import hep.physics.particle.properties.ParticleType;
-import org.lcsim.event.ParticleID;
-
-/**
- * An implementation of ParticleID appropriate for the fast Monte Carlo
- * @author ngraf
- */
-public class MCFastParticleID implements ParticleID
-{
- private ParticleType _type;
- /** Creates a new instance of MCFastParticleID */
- public MCFastParticleID(ParticleType type)
- {
- _type = type;
- }
-
- /** Type - defined to be the pdgId.
- */
- public int getType()
- {
- return _type.getPDGID();
- }
-
- /** The PDG code of this id - UnknownPDG ( 999999 ) if unknown.
- */
- public int getPDG()
- {
- return _type.getPDGID();
- }
-
- /**The likelihood of this hypothesis - in a user defined normalization.
- */
- public double getLikelihood()
- {
- return 1.;
- }
-
- /** Type of the algorithm/module that created this hypothesis.
- * Check/set collection parameters PIDAlgorithmTypeName and PIDAlgorithmTypeID.
- */
- public int getAlgorithmType()
- {
- return 0;
- }
-
- /** Parameters associated with this hypothesis.
- * Check/set collection paramter PIDParameterNames for decoding the indices.
- */
- public double[] getParameters()
- {
- return new double[1];
- }
-}
lcsim/src/org/lcsim/contrib/JanStrube/fastMC
diff -N MCFastReconstructedParticle.java
--- MCFastReconstructedParticle.java 29 Nov 2007 02:25:56 -0000 1.1
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,251 +0,0 @@
-package org.lcsim.contrib.JanStrube.fastMC;
-
-import org.lcsim.event.Vertex;
-//import org.lcsim.mc.fast.MCFast;
-import org.lcsim.spacegeom.CartesianPoint;
-import org.lcsim.spacegeom.SpacePoint;
-
-import hep.physics.vec.BasicHep3Vector;
-import hep.physics.vec.BasicHepLorentzVector;
-import hep.physics.vec.Hep3Vector;
-import hep.physics.vec.HepLorentzVector;
-import hep.physics.vec.VecOp;
-import hep.physics.particle.Particle;
-import hep.physics.particle.properties.ParticleType;
-import java.util.ArrayList;
-import java.util.List;
-import org.lcsim.event.Cluster;
-import org.lcsim.event.ParticleID;
-import org.lcsim.event.ReconstructedParticle;
-import org.lcsim.event.Track;
-
-import static java.lang.Math.sqrt;
-import static java.lang.Math.pow;
-import static java.lang.Math.abs;
-
-
-/**
- *
- * @author ngraf
- */
-public class MCFastReconstructedParticle implements ReconstructedParticle
-{
- // ReconstructedParticle attributes
- private double[] _covMatrix = new double[10];
- private double _mass;
- private double _charge;
- private double e_track;
- private double e_reco;
- private Particle _particle;
- private Hep3Vector _referencePoint;
- private Hep3Vector p3_track;
- private List<ParticleID> _particleIds = new ArrayList<ParticleID>();
- private ParticleID _particleIdUsed;
- private double _goodnessOfPid;
- private List<ReconstructedParticle> _particles = new ArrayList<ReconstructedParticle>();
- private List<Cluster> _clusters = new ArrayList<Cluster>();
- private List<Track> _tracks = new ArrayList<Track>();
- private BasicHepLorentzVector p_reco = new BasicHepLorentzVector();
- private BasicHepLorentzVector p_track = new BasicHepLorentzVector();
-
- public MCFastReconstructedParticle(Track t, ParticleType type, Particle p, Cluster assoc_c, double wtcal)
- {
-// MCFast.log.info(" PDGID= "+type.getPDGID()+" t.getPX,...= "+t.getPX()+" "+t.getPY()+" "+t.getPZ());
- _mass = type.getMass();
- addTrack(t);
- _charge = t.getCharge();
- // Use origin for reference point for now.
- _referencePoint = new SpacePoint();
- double[] trackP = t.getMomentum();
- e_track = sqrt(trackP[0]*trackP[0]+trackP[1]*trackP[1]+trackP[2]*trackP[2] + _mass*_mass);
- p_track.setV3(e_track, new CartesianPoint(t.getMomentum()));
- p3_track = p_track.v3();
- if (assoc_c != null)
- {
- addCluster(assoc_c);
- MCFast.log.info(" PDGID= "+type.getPDGID()+" e_track= "+e_track+" e_assoc_clus= "+assoc_c.getEnergy());
- MCFast.log.info(" PDGID= "+type.getPDGID()+" _referencePoint= "+_referencePoint.x()+" "+_referencePoint.y()+" "+_referencePoint.z());
- MCFast.log.info(" PDGID= "+type.getPDGID()+" p3_track= "+p3_track.x()+" "+p3_track.y()+" "+p3_track.z());
- }
- else
- {
- MCFast.log.info(" assoc_c = null PDGID= "+type.getPDGID());
- }
- if ( assoc_c != null && wtcal > 0. && abs(type.getPDGID()) != 11 )
- {
- e_reco = (1.-wtcal)*e_track + wtcal*assoc_c.getEnergy();
- if(e_reco < _mass) e_reco = _mass;
- p_reco.setV3(e_reco, VecOp.mult(sqrt(e_reco*e_reco - _mass*_mass),VecOp.unit(p3_track)));
- }
- else
- {
- e_reco = e_track;
- p_reco.setV3(e_reco, p3_track);
- }
-
- addParticleID(new MCFastParticleID(type));
- _particle = p;
- }
-
- public MCFastReconstructedParticle(Cluster c, ParticleType type, Particle p)
- {
- _mass = type.getMass();
- addCluster(c);
- double e = c.getEnergy();
- if (e < _mass)
- {
- if(e > _mass-1.e-5 && e > Double.MIN_VALUE)
- {
- _mass = e-Double.MIN_VALUE;
- }
- else
- {
- MCFast.log.warning(" MCFastReconstructedParticle e < _mass e= "+e+" _mass= "+_mass);
- MCFast.log.warning(" MCFastReconstructedParticle type = "+type) ;
- MCFast.log.warning(" MCFastReconstructedParticle program will continue, but problem should be fixed ");
- _mass = e-Double.MIN_VALUE;
- // System.exit(0);
- }
-
- }
-
- double pm = sqrt(e*e-_mass*_mass);
- // get direction from position of cluster and assume it comes from the origin
- double[] point = c.getPosition();
- double len = sqrt(point[0]*point[0] + point[1]*point[1] + point[2]*point[2]);
-
- _referencePoint = new BasicHep3Vector(0, 0, 0);
-
- double px = (pm/len)*(point[0]);
- double py = (pm/len)*(point[1]);
- double pz = (pm/len)*(point[2]);
- p_reco.setV3(e, px, py, pz);
- _charge = 0.;
-
- addParticleID(new MCFastParticleID(type));
- }
-
- public MCFastReconstructedParticle(double[] vxd, double[] mom, double mass, double charge, ParticleType type)
- {
- _mass = mass;
- _charge = charge;
- _referencePoint = new BasicHep3Vector(vxd);
- p3_track = new BasicHep3Vector(mom);
- e_reco = sqrt(pow(_mass,2)+p3_track.magnitudeSquared());
- p_reco.setV3(e_reco,p3_track);
- MCFast.log.info(" PDGID= "+type.getPDGID()+" e_reco= "+e_reco+" mass= "+_mass);
- MCFast.log.info(" PDGID= "+type.getPDGID()+" _referencePoint= "+_referencePoint.x()+" "+_referencePoint.y()+" "+_referencePoint.z());
- MCFast.log.info(" PDGID= "+type.getPDGID()+" p3_track= "+p3_track.x()+" "+p3_track.y()+" "+p3_track.z());
-
- addParticleID(new MCFastParticleID(type));
- }
-
- // ReconstructedParticle interface
-
- public int getType()
- {
- return _particleIdUsed.getType();
- }
-
- public Hep3Vector getMomentum()
- {
- return p_reco.v3();
- }
-
- public double getEnergy()
- {
- return p_reco.t();
- }
-
- public double[] getCovMatrix()
- {
- return _covMatrix;
- }
-
- public double getMass()
- {
- return _mass;
- }
-
- public double getCharge()
- {
- return _charge;
- }
-
- public Particle getParticle() {
- return _particle;
- }
-
- public Hep3Vector getReferencePoint()
- {
- return _referencePoint;
- }
-
- public List<ParticleID> getParticleIDs()
- {
- return _particleIds;
- }
-
- public ParticleID getParticleIDUsed()
- {
- return _particleIdUsed;
- }
-
- public double getGoodnessOfPID()
- {
- return _goodnessOfPid;
- }
-
- public List<ReconstructedParticle> getParticles()
- {
- return _particles;
- }
-
- public List<Cluster> getClusters()
- {
- return _clusters;
- }
-
- public List<Track> getTracks()
- {
- return _tracks;
- }
-
- public void addParticleID(ParticleID pid)
- {
- _particleIds.add(pid);
- _particleIdUsed = pid;
- }
-
- public void addParticle(ReconstructedParticle particle)
- {
- _particles.add(particle);
- }
-
- public void addCluster(Cluster cluster)
- {
- _clusters.add(cluster);
- }
-
- public void addTrack(Track track)
- {
- _tracks.add(track);
- }
-
- public HepLorentzVector asFourVector()
- {
- return p_reco;
- }
-
- public String toString()
- {
- StringBuffer sb = new StringBuffer("MCFastReconstructedParticle: \n");
- sb.append("E: "+getEnergy());
- return sb.toString();
- }
-
- public Vertex getStartVertex()
- {
- return null;
- }
-}
lcsim/src/org/lcsim/contrib/JanStrube/fastMC
diff -N MCFastReconstructedParticleDriver.java
--- MCFastReconstructedParticleDriver.java 29 Nov 2007 21:29:36 -0000 1.2
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,219 +0,0 @@
-package org.lcsim.contrib.JanStrube.fastMC;
-
-import hep.physics.particle.Particle;
-import hep.physics.particle.properties.ParticlePropertyManager;
-import hep.physics.particle.properties.ParticlePropertyProvider;
-import hep.physics.particle.properties.ParticleType;
-import java.util.ArrayList;
-import java.util.Map;
-import java.util.HashMap;
-import java.util.List;
-import java.util.Random;
-import org.lcsim.event.Cluster;
-import org.lcsim.event.EventHeader;
-import org.lcsim.event.ReconstructedParticle;
-import org.lcsim.event.Track;
-import org.lcsim.mc.fast.cluster.ronan.ReconHADCluster;
-import org.lcsim.util.Driver;
-import org.lcsim.conditions.ConditionsEvent;
-import org.lcsim.conditions.ConditionsListener;
-import org.lcsim.conditions.ConditionsSet;
-
-import static java.lang.Math.abs;
-import org.lcsim.mc.fast.cluster.ronan.ReconEMCluster;
-import org.lcsim.mc.fast.tracking.fix.FastMCTrack;
-import org.lcsim.util.aida.AIDA;
-/**
- *
- * @author ngraf
- */
-public class MCFastReconstructedParticleDriver extends Driver implements ConditionsListener
-{
- private ParticlePropertyProvider ppp;
- private IDResolutionTables IDEff;
- private AIDA aida = AIDA.defaultInstance();
- private ParticleType eminus;
- private ParticleType eplus;
- private ParticleType klong;
- private ParticleType muminus;
- private ParticleType muplus;
- private ParticleType neutron;
- private ParticleType photon;
- private ParticleType pizero;
- private ParticleType piplus;
- private ParticleType piminus;
- private ParticleType pplus;
- private ParticleType pminus;
- private ParticleType kplus;
- private ParticleType kminus;
-
- /** Creates a new instance of MCFastReconstructedParticleDriver */
- public MCFastReconstructedParticleDriver()
- {
- this(ParticlePropertyManager.getParticlePropertyProvider());
- }
-
- public MCFastReconstructedParticleDriver(ParticlePropertyProvider ppp)
- {
- this.ppp = ppp;
- //
- eminus = ppp.get(11);
- eplus = ppp.get(-11);
- klong = ppp.get(130);
- muminus = ppp.get(13);
- muplus = ppp.get(-13);
- neutron = ppp.get(2112);
- photon = ppp.get(22);
- pizero = ppp.get(111);
- piplus = ppp.get(211);
- piminus = ppp.get(-211);
- pplus = ppp.get(2212);
- pminus = ppp.get(-2212);
- kplus = ppp.get(321);
- kminus = ppp.get(-321);
- }
-
- protected void process(EventHeader event)
- {
-
- boolean hist = getHistogramLevel() > 0;
-
- if (IDEff == null)
- {
- ConditionsSet idconditions = getConditionsManager().getConditions("IDEfficiency");
- idconditions.addConditionsListener(this);
- IDEff = new IDResolutionTables(idconditions);
- }
-
- Random rand = getRandom();
-
- List<Track> tracks = event.getTracks();
- List<Cluster> clusters = event.getClusters();
-
- // Set up Track-Cluster association; for now cheat using MCParticle
- Map<Particle, Track> m_pt = new HashMap<Particle, Track>();
- Map<Particle, Cluster> m_pc = new HashMap<Particle, Cluster>();
- Map<Cluster, Track> m_ct = new HashMap<Cluster, Track>();
- Map<Track, Cluster> m_tc = new HashMap<Track, Cluster>();
-
- for(Track t : tracks) m_pt.put(((FastMCTrack)t).getMCParticle(),t);
- for(Cluster c : clusters) m_pc.put((c instanceof ReconEMCluster ? ((ReconEMCluster)c).getMCParticle() : ((ReconHADCluster)c).getMCParticle() ),c);
- for(Track t : tracks) m_tc.put(t,m_pc.get(((FastMCTrack)t).getMCParticle()));
- for(Cluster c : clusters) m_ct.put(c,m_pt.get(c instanceof ReconEMCluster ? ((ReconEMCluster)c).getMCParticle() : ((ReconHADCluster)c).getMCParticle() ));
-
- List<ReconstructedParticle> rpList = new ArrayList<ReconstructedParticle>();
- // start with the smeared tracks...
- for(Track t : tracks)
- {
- ParticleType type = null;
- if(t instanceof FastMCTrack)
- {
- FastMCTrack rt = (FastMCTrack) t;
- Particle p = rt.getMCParticle();
- int pdgid = p.getPDGID();
-
-
- // charged track id
- if((abs(pdgid)== 11) && (rand.nextDouble() < IDEff.getElectronEff()))
- {
- type = rt.getCharge() > 0 ? eplus : eminus;
- }
- else if((abs(pdgid)== 13) && (rand.nextDouble() < IDEff.getMuonEff()))
- {
- type = rt.getCharge() > 0 ? muplus : muminus;
- }
- else if((abs(pdgid)== 2212) && (rand.nextDouble() < IDEff.getProtonEff()))
- {
- type = rt.getCharge() > 0 ? pplus : pminus;
- }
- else if((abs(pdgid)== 321) && (rand.nextDouble() < IDEff.getKaonEff()))
- {
- type = rt.getCharge() > 0 ? kplus : kminus;
- }
- else
- {
- type = rt.getCharge() > 0 ? piplus : piminus;
- }
-
- if ((p.getEnergy() > 10) && (hist))
- {
- if (Math.abs(t.getTrackParameter(4)) < 1)
- {
- aida.histogram1D("track-particle", 150, -0.0003, 0.0003).fill( ( Math.sqrt(t.getPX()*t.getPX() + t.getPY()*t.getPY() + t.getPZ()*t.getPZ()) - Math.sqrt(p.getPX()*p.getPX() + p.getPY()*p.getPY() + p.getPZ()*p.getPZ()) ) / ( p.getPX()*p.getPX() + p.getPY()*p.getPY() + p.getPZ()*p.getPZ()) );
- }
- if ((Math.abs(t.getTrackParameter(4)) < 2.16) && (Math.abs(t.getTrackParameter(4)) > 1.96) && (p.getMomentum().magnitude() > 80) && (p.getMomentum().magnitude() < 120))
- {
- aida.histogram1D("track-particle-cut", 150, -0.01, 0.01).fill( ( Math.sqrt(t.getPX()*t.getPX() + t.getPY()*t.getPY() + t.getPZ()*t.getPZ()) - Math.sqrt(p.getPX()*p.getPX() + p.getPY()*p.getPY() + p.getPZ()*p.getPZ()) ) / (Math.sqrt( p.getPX()*p.getPX() + p.getPY()*p.getPY() + p.getPZ()*p.getPZ())) );
- }
- }
-
- // assume pion for remaining charged tracks
- MCFastReconstructedParticle rp = new MCFastReconstructedParticle(t, type, p, m_tc.get(t), IDEff.getWtChgTrkCal());
- rpList.add(rp);
-
- if (hist)
- {
- aida.histogram1D("recon-particle", 150, -5, 5).fill((rp.getEnergy()-p.getEnergy())/(Math.sqrt(p.getEnergy())));
- }
- }
- }
-
- // loop over clusters...
- for(Cluster c : clusters)
- {
- //if(m_ct.get(c) != null) continue;
- Particle p = null;
- ParticleType type = null;
- // photons for EM
- if( c instanceof ReconEMCluster)
- {
- ReconEMCluster emc = (ReconEMCluster) c;
- p = emc.getMCParticle();
- if(m_ct.get(c) != null) continue;
- type = photon;
- if (hist)
- {
- aida.histogram1D("photonCLS-particle", 150, -3, 3).fill((emc.getEnergy()-emc.getMCParticle().getEnergy())/(Math.sqrt(emc.getMCParticle().getEnergy())));
- }
- }
- // assume a KZeroLong here for had cluster
- else if(c instanceof ReconHADCluster)
- {
- ReconHADCluster emc = (ReconHADCluster) c;
- p = emc.getMCParticle();
- int pdgid = p.getPDGID();
- if ((abs(pdgid)==2112) && (rand.nextDouble() < IDEff.getNeutronEff()))
- {
- type = neutron;
- }
- else
- {
- type = klong;
- }
- if(m_ct.get(c) != null || c.getEnergy() < type.getMass()) continue;
- if (hist)
- {
- aida.histogram1D("hadronCLS-particle", 150, -3, 3).fill((emc.getEnergy()-emc.getMCParticle().getEnergy())/(Math.sqrt(emc.getMCParticle().getEnergy())));
- }
-
- }
- MCFastReconstructedParticle rp = new MCFastReconstructedParticle(c, type, p);
- rpList.add(rp);
- if (hist)
- {
- aida.histogram1D("recon-particle", 150, -10, 10).fill(rp.getEnergy()-p.getEnergy());
- }
-
- }
- // add the reconstructedparticles to the event
- event.put(event.MCFASTRECONSTRUCTEDPARTICLES, rpList, ReconstructedParticle.class, 0);
-
- }
-
- public void conditionsChanged(ConditionsEvent event)
- {
- ConditionsSet idconditions = getConditionsManager().getConditions("IDEfficiency");
- IDEff = new IDResolutionTables(idconditions);
- }
-
-}
lcsim/src/org/lcsim/contrib/JanStrube/tracking
diff -N TrackTester.java
--- TrackTester.java 29 Nov 2007 02:25:57 -0000 1.2
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,23 +0,0 @@
-/**
- * @version $Id: TrackTester.java,v 1.2 2007/11/29 02:25:57 jstrube Exp $
- */
-package org.lcsim.contrib.JanStrube.tracking;
-
-import java.util.Collection;
-
-import org.lcsim.event.Track;
-import org.lcsim.util.swim.HelixSwimmer;
-
-/**
- * @author jstrube
- *
- */
-public class TrackTester {
- public static void testTrackList(Collection<? extends Track> trackList) {
- for (Track t : trackList) {
- HelixSwimmer swimmer = new HelixSwimmer(5);
- swimmer.setTrack(t);
- t.getMomentum();
- }
- }
-}
lcsim/src/org/lcsim/contrib/JanStrube/vtxFitter
diff -N Fitter.java
--- Fitter.java 29 Nov 2007 21:29:39 -0000 1.18
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,365 +0,0 @@
-package org.lcsim.contrib.JanStrube.vtxFitter;
-
-/**
- * @version $Id: Fitter.java,v 1.18 2007/11/29 21:29:39 jstrube Exp $
- */
-
-import Jama.util.Maths;
-import static java.lang.Math.abs;
-import static java.lang.Math.atan2;
-import static java.lang.Math.sqrt;
-import static org.lcsim.constants.Constants.fieldConversion;
-import static org.lcsim.event.LCIOParameters.SpaceMomentum2Parameters;
-
-
-import java.util.Collection;
-
-import org.lcsim.event.EventHeader;
-import org.lcsim.event.Track;
-import org.lcsim.mc.fast.tracking.fix.FastMCTrackFactory;
-import org.lcsim.spacegeom.CartesianPoint;
-import org.lcsim.spacegeom.CartesianVector;
-import org.lcsim.spacegeom.SpacePoint;
-import org.lcsim.spacegeom.SpaceVector;
-import org.lcsim.util.swim.HelixSwimmer;
-
-import Jama.Matrix;
-
-/**
- * This is an implementation of a vertex fitter following the paper by Grab and
- * Luchsinger. The parametrization of the Tracks follows the L3 convention. The
- * fitter takes a collection of tracks as input, then performs the sequence of
- * filtering and smoothing until convergence of the vertex in chi<sup>2</sup>.
- * The expensive computation of the full covariance matrix is only done once,
- * after convergence.
- *
- * @author jstrube
- *
- */
-public class Fitter {
- // The Vertex that is being fitted
- // TODO: maybe the filter and smoothe functions should be part of Vertex...
- Vertex fitVertex = new Vertex();
- FastMCTrackFactory _trackFactory;
- double B_z;
- // convergence criterion
- double delta_chi2 = 1e-7;
- int max_iter = 30;
- // position
- Matrix x_prev = new Matrix(new double[] {0,0,0}, 3);
- // covariance for x
- Matrix C_prev = new Matrix(new double[][] { { 10000, 0, 0 }, { 0, 10000, 0 },
- { 0, 0, 10000 } });
-
-
- HelixSwimmer _swimmer;
-
- // chi2
- double chi2_prev;
-
- /**
- * Definitions:<br />
- * x_k = estimate of the vertex position after using the
- * information of k tracks<br />
- * xt = true vertex position<br />
- * C_k = cov(x_k) = cov(xk - xt)<br />
- * qk = estimate of the momentum of particle k at xt<br />
- * qkt = true momentum of particle k at xt<br />
- * D_k = cov(qk) = cov(qk - qkt)<br />
- * Ek = cov(xk, qk) = cov((xk-xt)(qk-qkt))<br />
- * mk = measurement vector<br />
- * vk = measurement noise<br />
- * Vk = cov(vk)<br />
- * Gk = Vk^-1 weightmatrix of particle k
- *
- * @param p
- */
-
- public Fitter(EventHeader event) {
- B_z = event.getDetector().getFieldMap().getField(
- new double[] { 0, 0, 0 })[2];
- _swimmer = new HelixSwimmer(B_z);
- _trackFactory = new FastMCTrackFactory(event, false);
- }
-
- // This is for unit tests only
- // NEVER USE THIS CONSTRUCTOR
- Fitter(double bz) {
- B_z = bz;
- _trackFactory = new FastMCTrackFactory("sid01", bz, false);
- _swimmer = new HelixSwimmer(bz);
- }
-
- public Vertex fit(Collection<? extends Track> tracks, SpacePoint initialPosition) {
- x_prev = new Matrix(initialPosition.getCartesianArray(), 3);
- // The Vertex object
- Vertex particle = new Vertex(tracks, initialPosition);
- // set the conditions to enter the loop
- chi2_prev = 1000;
- double chi2_this = 0;
- int n_iter = 0;
- while (abs(chi2_prev - chi2_this) > delta_chi2 && n_iter < max_iter) {
- // System.out.printf("iteration: %d\tdelta chi2: %.3f\n", n_iter,
- // chi2_prev - chi2_this);
- n_iter++;
- // the Vertex that is being worked on in this iteration
- fitVertex = new Vertex();
- fitVertex._location = particle._location;
- // System.out.printf("before: %f\t after: %f\n", chi2_this,
- // chi2_prev);
- chi2_this = chi2_prev;
- chi2_prev = 0;
- for (VtxTrack t : particle.getVtxTracks()) {
- fitVertex.addTrack(t, filter(t));
- }
- // FIXME: This only changes the parameters of the tracks, not the error matrix.
- for (VtxTrack t : fitVertex.getVtxTracks()) {
- smoothe(t, false);
- }
- // the particle is now the currently best approximation
- particle = fitVertex;
- }
- particle._spatialCovarianceMatrix = C_prev;
- particle._chi2 = chi2_this;
- // final iteration, this time compute the full covariance matrix
-// for (VtxTrack t : particle.getVtxTracks()) {
-// smoothe(t, true);
-// }
- return particle;
- }
-
- // Turn everything into a Matrix
- double filter(VtxTrack t) {
- // First, get the Matrices at the point of linearization x0, p0
- // a good choice for x0 is the previous best estimate and
- // p0 is the momentum of the track at the POCA to x0
- _swimmer.setTrack(t);
- SpacePoint vtxPosition = fitVertex.location();
- // at the moment, particle.origin is the best estimate so far
- double alpha = _swimmer.getTrackLengthToPoint(vtxPosition);
-// System.out.printf("alpha: %f\n", alpha);
- SpaceVector p = _swimmer.getMomentumAtLength(alpha);
- // x is the POCA of the current track to the best estimate
-
- SpacePoint x = _swimmer.getPointAtLength(alpha);
-// System.out.printf("x: %s\n", x);
- Matrix A_k = getSpatialDerivativeMatrix(vtxPosition, p, t.referencePoint(), t.getCharge(), B_z);
- Matrix B_k = getMomentumDerivativeMatrix(vtxPosition, p, t.referencePoint(), t.getCharge(), B_z);
- Matrix q0 = new Matrix(p.getCartesianArray(), 3);
-// System.out.printf("q0: %s\n", q0);
- // turn the measurement into a Matrix object
- Matrix h = new Matrix(SpaceMomentum2Parameters(vtxPosition,
- p, t.referencePoint(), t.getCharge(), B_z).getValues(), 5);
- // the measurement of the track is taken at the POCA to the best
- // estimate
-// System.out.printf("h: %s\n", h);
- Matrix m_k = new Matrix(SpaceMomentum2Parameters(x, p, t.referencePoint(), t.getCharge(), B_z).getValues(), 5);
- Matrix c_k0 = h.minus(A_k.times(x_prev).plus(B_k.times(q0)));
-
- // Now for the real filtering equations
- Matrix G_k = Maths.toJamaMatrix(t.getErrorMatrix()).inverse();
- Matrix W_k = B_k.transpose().times(G_k.times(B_k)).inverse();
- Matrix G_kB = G_k.minus((G_k.times(B_k.times(W_k.times(B_k.transpose()
- .times(G_k))))));
-
- // cov(x)
-// System.out.printf("C_k before: %s\n", C_prev);
- Matrix C_k = C_prev.inverse().plus(
- A_k.transpose().times(G_kB.times(A_k))).inverse();
-// System.out.printf("C_k after: %s\n", C_k);
-// System.out.printf("G_kb: %s\n", G_kB);
-// System.out.printf("m_k: %s\n", m_k);
-// System.out.printf("c_k0: %s \n", c_k0);
-// System.out.printf("x before: %s\n", x_prev);
- Matrix x_k = C_k.times(C_prev.inverse().times(x_prev).plus(
- A_k.transpose().times(G_kB.times(m_k.minus(c_k0)))));
-// System.out.printf("x after: %s\n", x_k);
- Matrix q_k = W_k.times(B_k.transpose().times(
- G_k.times(m_k.minus(c_k0.plus(A_k.times(x_k))))));
-
- SpaceVector updatedMomentum = new CartesianVector(q_k
- .getRowPackedCopy());
- // cov(q)
- // now the chi2
- Matrix r_k = m_k.minus(c_k0.plus(A_k.times(x_k).plus(B_k.times(q_k))));
- // System.out.printf("chi2: %s + %s",
- // r_k.transpose().times(G_k.times(r_k)),
- // x_k.minus(x_prev).transpose().times(C_prev.inverse().times(x_k.minus(x_prev))));
- double chi2_kf = r_k.transpose().times(G_k.times(r_k)).get(0, 0)
- + x_k.minus(x_prev).transpose().times(
- C_prev.inverse().times(x_k.minus(x_prev))).get(0, 0);
-
- chi2_prev += chi2_kf;
-// System.out.printf("before: %s", x_prev);
- x_prev = x_k;
-// System.out.printf("after: %s", x_prev);
- C_prev = C_k;
-// System.out.println(C_prev);
- fitVertex._location = new CartesianPoint(x_k.getRowPackedCopy());
- // FIXME: The error matrix is not updated
- //t.updateMomentum(updatedMomentum, fitVertex._location, B_z);
- // FIXME: This is the easiest way to get an updated error matrix, but it should be done the Right Way (TM)
- //t = new VtxTrack(_trackFactory.getTrack(updatedMomentum, vtxPosition, t.getCharge()));
- // System.out.println(particle._origin);
- // if (fullFit) {
- // Matrix D_k = W_k.plus(W_k.times(B_k.transpose().times(
- // G_k.times(A_k.times(C_k.times(A_k.transpose().times(G_k.times(B_k.times(W_k)))))))));
- // // System.out.printf("D_k: %g\t%g\t%g\n", D_k.get(0, 0), D_k.get(1,
- // 1), D_k.get(2, 2));
- // // cov(x, q)
- // Matrix E_k =
- // C_k.times(A_k.transpose().times(G_k.times(B_k.times(W_k)))).uminus();
- // }
- return chi2_kf;
- }
-
- // smoothing is performed after every loop over the set of tracks.
- void smoothe(VtxTrack t, boolean fullFit) {
- _swimmer.setTrack(t);
- SpacePoint vtxPosition = fitVertex.location();
- // at the moment, particle.origin is the best estimate so far
- double alpha = _swimmer.getTrackLengthToPoint(vtxPosition);
- SpaceVector p = _swimmer.getMomentumAtLength(alpha);
- // x is the POCA of the current track to the best estimate
- Matrix q0 = new Matrix(p.getCartesianArray(), 3);
- SpacePoint x = _swimmer.getPointAtLength(alpha);
-
- Matrix A_k = getSpatialDerivativeMatrix(vtxPosition, p, t.referencePoint(), t.getCharge(), B_z);
- Matrix B_k = getMomentumDerivativeMatrix(vtxPosition, p, t.referencePoint(), t.getCharge(), B_z);
-
- Matrix G_k = Maths.toJamaMatrix(t.getErrorMatrix()).inverse();
- Matrix W_k = B_k.transpose().times(G_k.times(B_k)).inverse();
- // the measurement of the track is taken at the POCA to the best
- // estimate
- Matrix m_k = new Matrix(SpaceMomentum2Parameters(x, p, t.referencePoint(), t.getCharge(), B_z).getValues(), 5);
- Matrix h = new Matrix(SpaceMomentum2Parameters(vtxPosition, p, t.referencePoint(), t.getCharge(), B_z).getValues(), 5);
-
- Matrix c_k0 = h.minus(A_k.times(x_prev).plus(B_k.times(q0)));
- Matrix q_kN = W_k.times(B_k.transpose().times(G_k)).times(
- m_k.minus(c_k0.plus(A_k.times(x_prev))));
- //SpaceVector updatedMomentum = new CartesianVector(q_kN
-// .getRowPackedCopy());
-// System.out.printf("old momentum: %s\n", q0);
-// System.out.printf("new momentum: %s\n", q_kN);
-
- //t.updateMomentum((FastMCTrack)_trackFactory.getTrack(updatedMomentum, vtxPosition, vtxPosition, t.getCharge()));
- if (fullFit) {
- Matrix D_kN = W_k.plus(W_k.times(B_k.transpose().times(
- G_k.times(A_k.times(C_prev.times(A_k.transpose().times(
- G_k.times(B_k.times(W_k)))))))));
- Matrix E_kN = C_prev.times(
- A_k.transpose().times(G_k.times(B_k.times(W_k)))).uminus();
- }
- }
-
- // // returns y = A.x
- // static double[] multiply(Matrix A, double[] x) {
- // if (A.getColumnDimension() != x.length)
- // throw new IllegalArgumentException("dimensions do not match");
- // double[] result = new double[A.getRowDimension()];
- //
- // for (int i=0; i<A.getRowDimension(); i++) {
- // for (int j=0; j<A.getColumnDimension(); j++) {
- // result[i] += A.get(i, j) * x[j];
- // }
- // }
- // return result;
- // }
- //
- // // for notational convenience
- // static Matrix multiply(Matrix a, Matrix b) {
- // return a.times(b);
- // }
-
- public static Matrix getMomentumDerivativeMatrix(SpacePoint x,
- SpaceVector p, SpacePoint referencePoint, int charge, double Bz) {
- double field = Bz * fieldConversion;
- double pt = p.rxy();
- double px0 = p.x() - charge * field * (x.y()-referencePoint.y());
- double py0 = p.y() + charge * field * (x.x()-referencePoint.x());
- double pt0 = sqrt(px0 * px0 + py0 * py0);
- double phi = atan2(p.y(), p.x());
- double phi0 = atan2(py0, px0);
- double l = (phi - phi0) * pt / (charge * field);
-
- double[] dTanLambda_dp = new double[] {
- -p.z() * p.x() / (pt * pt * pt),
- -p.z() * p.y() / (pt * pt * pt), 1 / pt };
- double[] dPhi_dp = new double[] { -py0 / (pt0 * pt0),
- px0 / (pt0 * pt0), 0 };
- double[] dOmega_dp = new double[] {
- -charge * field * p.x() / (pt * pt * pt),
- -charge * field * p.y() / (pt * pt * pt), 0 };
- double[] dD0_dp = new double[] { -py0 / pt0, px0 / pt0, 0 };
- double[] dZ_dp = new double[] {
- -p.z() / (field * charge)
- * (p.y() / (pt * pt) - py0 / (pt0 * pt0)),
- -p.z() / (field * charge)
- * (p.x() / (pt * pt) - px0 / (pt0 * pt0)), -l / pt };
- return new Matrix(new double[][] { dD0_dp, dPhi_dp, dOmega_dp, dZ_dp,
- dTanLambda_dp });
- }
-
- /**
- * Calculates the derivatives of the measurement vector with respect to
- * cartesian spatial coordinates.
- *
- * @param x
- * @param p
- * @param charge
- * @return the derivatives Matrix. Assumes an ordering of (d0, phi0, omega,
- * z, tanLambda) in the measurement vector
- */
- public static Matrix getSpatialDerivativeMatrix(SpacePoint x,
- SpaceVector p, SpacePoint referencePoint, int charge, double Bz) {
- double field = Bz * fieldConversion;
-
- double px0 = p.x() - charge * field * (x.y() - referencePoint.y());
- double py0 = p.y() + charge * field * (x.x() - referencePoint.x());
- double pt0 = sqrt(px0 * px0 + py0 * py0);
-
- // be very explicit about the derivatives
- // dr = (dx, dy, dz)^T
- double[] dTanLambda_dr = new double[] { 0, 0, 0 };
- double[] dPhi_dr = new double[] { px0 * field * charge / (pt0 * pt0),
- py0 * field * charge / (pt0 * pt0), 0 };
- double[] dOmega_dr = new double[] { 0, 0, 0 };
- double[] dD0_dr = new double[] { -py0 / pt0, px0 / pt0, 0 };
- double[] dZ_dr = new double[] { -p.z() * px0 / (pt0 * pt0),
- -p.z() * py0 / (pt0 * pt0), 1 };
-
- return new Matrix(new double[][] { dD0_dr, dPhi_dr, dOmega_dr, dZ_dr,
- dTanLambda_dr });
- }
-
- // for notational convenience
- static Matrix add(Matrix a, Matrix b) {
- return a.plus(b);
- }
-
- // matrix A_c in the Grab and Luchsinger paper
- static Matrix chargedSpatialLinearity(Track t) {
-
- return new Matrix(0, 0);
- }
-
- static Matrix chargedMomentumLinearity(Track t) {
- return new Matrix(0, 0);
- }
-
- static Matrix calcSpatialLinearityMatrix(Track t) {
- return new Matrix(5, 5);
- }
-
- static Matrix calcMomentumLinearityMatrix(Track t) {
- return new Matrix(5, 5);
- }
-
- public double getDelta_chi2() {
- return delta_chi2;
- }
-
- public void setDelta_chi2(double d) {
- delta_chi2 = d;
- }
-}
lcsim/src/org/lcsim/contrib/JanStrube/vtxFitter
diff -N FitterTestDriver.java
--- FitterTestDriver.java 29 Nov 2007 21:29:38 -0000 1.3
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,173 +0,0 @@
-package org.lcsim.contrib.JanStrube.vtxFitter;
-
-import static java.lang.Math.abs;
-import static java.lang.Math.sqrt;
-
-import hep.physics.matrix.SymmetricMatrix;
-
-import java.util.ArrayList;
-import java.util.List;
-import java.util.Random;
-
-import org.lcsim.contrib.JanStrube.vtxFitter.Fitter;
-import org.lcsim.contrib.JanStrube.vtxFitter.Vertex;
-import org.lcsim.event.EventHeader;
-import org.lcsim.event.MCParticle;
-import org.lcsim.event.LCIOParameters;
-import org.lcsim.event.LCIOParameters.ParameterName;
-import org.lcsim.event.Track;
-import org.lcsim.mc.fast.tracking.fix.FastMCTrack;
-import org.lcsim.mc.fast.tracking.fix.FastMCTrackFactory;
-import org.lcsim.spacegeom.CartesianPoint;
-import org.lcsim.spacegeom.CartesianVector;
-import org.lcsim.spacegeom.SpacePoint;
-import org.lcsim.spacegeom.SpaceVector;
-import org.lcsim.util.Driver;
-import org.lcsim.util.aida.AIDA;
-
-import Jama.EigenvalueDecomposition;
-import Jama.Matrix;
-
-public class FitterTestDriver extends Driver {
- private AIDA aida = AIDA.defaultInstance();
- private SpacePoint referencePoint;
- private SpacePoint origin;
- private Random rand;
- private int count = 0;
-
- public FitterTestDriver() {
- referencePoint = new CartesianPoint(0, 0, 0);
- origin = new CartesianPoint(0, 0, 0);
- rand = new Random(123);
- }
-
- public void process(EventHeader event) {
- count++;
- if (count % 500 == 0)
- System.out.printf("Processing event %d\n", count);
- // beamconstraint when smearing the tracks ?
- FastMCTrackFactory tf = new FastMCTrackFactory(event, false);
- for (MCParticle part : event.getMCParticles()) {
- // K0: 310, J/Psi: 443, B+: 521, B0: 511
- if (abs(part.getPDGID()) != 521 && abs(part.getPDGID()) != 511)
- continue;
- List<Track> unsmearedInput = new ArrayList<Track>();
- List<Track> smearedInput = new ArrayList<Track>();
- origin = new SpacePoint(part.getEndPoint());
-
- for (MCParticle daughter : part.getDaughters()) {
- if (daughter.getCharge() == 0)
- continue;
- FastMCTrack unsmearedTrack = (FastMCTrack) tf.getTrack(new SpaceVector(daughter
- .getMomentum()), new SpacePoint(daughter.getOrigin()),
- new SpacePoint(), (int) daughter
- .getCharge(), rand, false);
- // System.err.println("Particle: charge " +
- // daughter.getCharge());
- // System.err.println("Position:" + new
- // SpacePoint(daughter.getOrigin()));
- // System.err.println("Momentum:" + new
- // SpacePoint(daughter.getMomentum()));
- SpacePoint unsmearedPosition = LCIOParameters
- .Parameters2Position(unsmearedTrack.getParameters(),
- referencePoint);
- SpacePoint unsmearedMomentum = new CartesianVector(LCIOParameters
- .Parameters2Momentum(unsmearedTrack.getParameters()).v());
- // System.err.println("unsmeared Position: " +
- // unsmearedPosition);
- // System.err.println("unsmeared Momentum: " +
- // unsmearedMomentum);
- unsmearedInput.add(unsmearedTrack);
- aida.cloud1D("track_unsmeared_pt").fill(unsmearedTrack.getPt());
- aida.cloud1D("track_unsmeared_diffx").fill(
- unsmearedPosition.x() - daughter.getOriginX());
- aida.cloud1D("track_unsmeared_diffy").fill(
- unsmearedPosition.y() - daughter.getOriginY());
- aida.cloud1D("track_unsmeared_diffz").fill(
- unsmearedPosition.z() - daughter.getOriginZ());
- aida.cloud1D("track_unsmeared_diffpx").fill(
- unsmearedMomentum.x() - daughter.getPX());
- aida.cloud1D("track_unsmeared_diffpy").fill(
- unsmearedMomentum.y() - daughter.getPY());
- aida.cloud1D("track_unsmeared_diffpz").fill(
- unsmearedMomentum.z() - daughter.getPZ());
-
- FastMCTrack smearedTrack = (FastMCTrack) tf.getTrack(new SpaceVector(daughter
- .getMomentum()), new SpacePoint(daughter.getOrigin()),
- new SpacePoint(), (int) daughter
- .getCharge(), rand, true);
- SpacePoint smearedPosition = LCIOParameters
- .Parameters2Position(smearedTrack.getParameters(),
- referencePoint);
- SpacePoint smearedMomentum = new CartesianVector(LCIOParameters
- .Parameters2Momentum(smearedTrack.getParameters()).v());
- // System.err.println("smeared Position: " + smearedPosition);
- // System.err.println("smeared Momentum: " + smearedMomentum);
- smearedInput.add(smearedTrack);
- SymmetricMatrix error = smearedTrack.getErrorMatrix();
- //error.invert();
- aida.cloud1D("track_pull_d0").fill((unsmearedTrack.getParameter(ParameterName.d0)-smearedTrack.getParameter(ParameterName.d0))/sqrt(error.diagonal(ParameterName.d0.ordinal())));
- aida.cloud1D("track_pull_phi0").fill((unsmearedTrack.getParameter(ParameterName.phi0)-smearedTrack.getParameter(ParameterName.phi0))/sqrt(error.diagonal(ParameterName.phi0.ordinal())));
- aida.cloud1D("track_pull_omega").fill((unsmearedTrack.getParameter(ParameterName.omega)-smearedTrack.getParameter(ParameterName.omega))/sqrt(error.diagonal(ParameterName.omega.ordinal())));
- aida.cloud1D("track_pull_z").fill((unsmearedTrack.getParameter(ParameterName.z0)-smearedTrack.getParameter(ParameterName.z0))/sqrt(error.diagonal(ParameterName.z0.ordinal())));
- aida.cloud1D("track_pull_tanLambda").fill((unsmearedTrack.getParameter(ParameterName.tanLambda)-smearedTrack.getParameter(ParameterName.tanLambda))/sqrt(error.diagonal(ParameterName.tanLambda.ordinal())));
-
- aida.cloud1D("track_smeared_pt").fill(smearedTrack.getPt());
- aida.cloud1D("track_smeared_diffx").fill(
- smearedPosition.x() - daughter.getOriginX());
- aida.cloud1D("track_smeared_diffy").fill(
- smearedPosition.y() - daughter.getOriginY());
- aida.cloud1D("track_smeared_diffz").fill(
- smearedPosition.z() - daughter.getOriginZ());
- aida.cloud1D("track_smeared_diffpx").fill(
- smearedMomentum.x() - daughter.getPX());
- aida.cloud1D("track_smeared_diffpy").fill(
- smearedMomentum.y() - daughter.getPY());
- aida.cloud1D("track_smeared_diffpz").fill(
- smearedMomentum.z() - daughter.getPZ());
- aida.cloud1D("track_diff_d0").fill(
- unsmearedTrack.getParameters().getValues()[0]
- - smearedTrack.getParameters().getValues()[0]);
- aida.cloud1D("track_diff_phi0").fill(
- unsmearedTrack.getParameters().getValues()[1]
- - smearedTrack.getParameters().getValues()[1]);
- aida.cloud1D("track_diff_omega").fill(
- unsmearedTrack.getParameters().getValues()[2]
- - smearedTrack.getParameters().getValues()[2]);
- aida.cloud1D("track_diff_z").fill(
- unsmearedTrack.getParameters().getValues()[3]
- - smearedTrack.getParameters().getValues()[3]);
- aida.cloud1D("track_diff_tanLambda").fill(
- unsmearedTrack.getParameters().getValues()[4]
- - smearedTrack.getParameters().getValues()[4]);
- }
- //Fitter fitter = new Fitter(event);
- //Vertex newVtx1 = fitter.fit(unsmearedInput, new SpacePoint());
- Fitter fitter2 = new Fitter(event);
- Vertex newVtx2 = fitter2.fit(smearedInput, new SpacePoint());
- double sigmaX = sqrt(newVtx2.getSpatialCovarianceMatrix().get(0, 0));
- double sigmaY = sqrt(newVtx2.getSpatialCovarianceMatrix().get(1, 1));
- double sigmaZ = sqrt(newVtx2.getSpatialCovarianceMatrix().get(2, 2));
- double deltaX = newVtx2.location().x() - part.getEndPoint().x();
- double deltaY = newVtx2.location().y() - part.getEndPoint().y();
- double deltaZ = newVtx2.location().z() - part.getEndPoint().z();
- aida.cloud1D("decayLength").fill(
- new SpacePoint(part.getEndPoint()).magnitude());
- aida.cloud1D("vtx_smeared_deltaX").fill(deltaX);
- aida.cloud1D("vtx_smeared_deltaY").fill(deltaY);
- aida.cloud1D("vtx_smeared_deltaZ").fill(deltaZ);
- aida.cloud1D("vtx_sigma_x").fill(sigmaX);
- aida.cloud1D("vtx_sigma_y").fill(sigmaY);
- aida.cloud1D("vtx_sigma_z").fill(sigmaZ);
- aida.cloud1D("vtx_pull_x").fill(deltaX / sigmaX);
- aida.cloud1D("vtx_pull_y").fill(deltaY / sigmaY);
- aida.cloud1D("vtx_pull_z").fill(deltaZ / sigmaZ);
-// aida.cloud1D("vtx_unsmeared_deltaX").fill(
-// newVtx1.location().x() - part.getEndPoint().x());
-// aida.cloud1D("vtx_unsmeared_deltaY").fill(
-// newVtx1.location().y() - part.getEndPoint().y());
-// aida.cloud1D("vtx_unsmeared_deltaZ").fill(
-// newVtx1.location().z() - part.getEndPoint().z());
- }
- }
-}
lcsim/src/org/lcsim/contrib/JanStrube/vtxFitter
diff -N Vertex.java
--- Vertex.java 29 Nov 2007 02:25:59 -0000 1.8
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,128 +0,0 @@
-package org.lcsim.contrib.JanStrube.vtxFitter;
-
-import java.util.Collection;
-import java.util.HashMap;
-import java.util.HashSet;
-import java.util.Map;
-import java.util.Set;
-
-import org.lcsim.event.Track;
-import org.lcsim.spacegeom.SpacePoint;
-
-import Jama.Matrix;
-
-public class Vertex {
- SpacePoint _location;
- SpacePoint _locationError;
- Map<VtxTrack, Double> _tracks;
-
- double _chi2;
- double _ndf;
- Matrix _fullCovarianceMatrix;
- Matrix _spatialCovarianceMatrix;
- Matrix _momentumCovarianceMatrix;
-
- public Vertex() {
- _location = new SpacePoint();
- _locationError = new SpacePoint();
- _tracks = new HashMap<VtxTrack, Double>();
- }
-
- public Vertex(Collection<? extends Track> tracks, SpacePoint initialPos) {
- this();
- setTracks(tracks);
- _location = initialPos;
- }
-
- public void setLocation(SpacePoint x) {
- _location = x;
- }
-
- public SpacePoint location() {
- return _location;
- }
-
- void addTrack(VtxTrack t, double chi2) {
- _tracks.put(t, chi2);
- }
-
- void setTracks(Map<VtxTrack, Double> tracks) {
- _tracks = tracks;
- }
-
- void setTracks(Collection<? extends Track> tracks) {
- for (Track t : tracks) {
- addTrack(new VtxTrack(t), 1000);
- }
- }
-
- // public void setTracks(List<Track> tracks) {
- // _tracks.clear();
- // for (Track t : tracks) {
- // _tracks.put(t, 0.);
- // }
- // }
-
- public double getChi2() {
- return _chi2;
- }
-
- public void setMaxChi2(double chi2) {
- _chi2 = chi2;
- for (Track t : _tracks.keySet()) {
- if (_tracks.get(t) > _chi2) {
- _tracks.remove(t);
- }
- }
- }
-
- public double getNdf() {
- return _ndf;
- }
-
- public void setNdf(double ndf) {
- _ndf = ndf;
- }
-
- public Set<VtxTrack> getVtxTracks() {
- return _tracks.keySet();
- }
-
- public Set<Track> getTracks() {
- Set<Track> set = new HashSet<Track>();
- for (Track t : _tracks.keySet()) {
- set.add(t);
- }
- return set;
- }
-
- /**
- * @return The full covarianceMatrix. The covariance matrices for position
- * and momentum are submatrices of this Matrix.
- */
- public Matrix getCovarianceMatrix() {
- throw new UnsupportedOperationException();
- }
-
- /**
- * @return The spatial components of the full covariance matrix
- */
- public Matrix getSpatialCovarianceMatrix() {
- return _spatialCovarianceMatrix;
- }
-
- /**
- * @return The momentum components of the full covariance matrix
- */
- public Matrix getMomentumCovarianceMatrix() {
- throw new UnsupportedOperationException();
- }
-
- public String toString() {
- StringBuilder output = new StringBuilder();
- output.append(String.format("Covariance Matrix:\n%s\n", _spatialCovarianceMatrix));
- output.append(String.format("Chi2: %.4f\n", _chi2));
- output.append(String.format("nDOF: %.0f\n", _ndf));
- return output.toString();
- }
-}
lcsim/src/org/lcsim/contrib/JanStrube/vtxFitter
diff -N VtxTrack.java
--- VtxTrack.java 29 Nov 2007 21:29:39 -0000 1.6
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,29 +0,0 @@
-package org.lcsim.contrib.JanStrube.vtxFitter;
-
-import hep.physics.matrix.SymmetricMatrix;
-
-import org.lcsim.event.Track;
-import org.lcsim.mc.fast.tracking.fix.FastMCTrack;
-import org.lcsim.spacegeom.SpacePoint;
-import org.lcsim.spacegeom.SpaceVector;
-
-import Jama.Matrix;
-import Jama.util.Maths;
-
-public class VtxTrack extends FastMCTrack {
-
- /*
- * Copy constructor. The members of Track have to be cloned,
- * because they will be modified in an irreversible way
- */
- VtxTrack(Track t) {
- super(t);
- }
-
- void updateMomentum(FastMCTrack t) {
- _parameters = t.getParameters();
- _errorMatrix = t.getErrorMatrix();
- _referencePoint = t.referencePoint();
- _charge = t.getCharge();
- }
-}
lcsim/src/org/lcsim/contrib/JanStrube/zvtop
diff -N ZvMaximum.java
--- ZvMaximum.java 29 Nov 2007 02:25:56 -0000 1.2
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,67 +0,0 @@
-package org.lcsim.contrib.JanStrube.zvtop;
-
-import org.lcsim.event.Track;
-import org.lcsim.spacegeom.SpacePoint;
-
-/**
- * Maximum in the overlap of two tracks
- * @author jstrube
- * @version $Id: ZvMaximum.java,v 1.2 2007/11/29 02:25:56 jstrube Exp $
- */
-class ZvMaximum implements Comparable {
- private Track _track_i;
- private Track _track_j;
- private SpacePoint _location;
- private double _value;
-
- public ZvMaximum(Track track_i, Track track_j, SpacePoint loc, double overlap) {
- _track_i = track_i;
- _track_j = track_j;
- _location = loc;
- // TODO for debugging only
- _value = overlap;
- }
-
- /**
- * Sets the new location of the ZvMaximum (used by clusterCandidates)
- * @param loc new location of the maximum. It is assumed that the
- * overlap doesn't change enough between the two points for the order to
- * change in the ZvMaximumMatrix
- * Not intended for public consumption
- */
- void setLocation(SpacePoint loc) {
- _location = loc;
- }
-
- public SpacePoint getLocation() {
- return _location;
- }
-
- public double getValue() {
- return _value;
- }
-
- public Track getTrack_i() {
- return _track_i;
- }
-
- public Track getTrack_j() {
- return _track_j;
- }
-
- public int compareTo(Object o) {
- if (o.getClass() != ZvMaximum.class)
- throw new IllegalArgumentException("Can only compare ZvMaximum Objects");
- ZvMaximum max2 = (ZvMaximum) o;
- if (_value > max2._value)
- return 1;
- if (_value < max2._value)
- return -1;
- return 0;
- }
-
- public String toString() {
- String s = String.format("value: %f\nlocation: %s\n", _value, _location);
- return s;
- }
-}
lcsim/src/org/lcsim/contrib/JanStrube/zvtop
diff -N ZvMaximumMatrix.java
--- ZvMaximumMatrix.java 9 Oct 2006 17:25:27 -0000 1.1
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,72 +0,0 @@
-package org.lcsim.contrib.JanStrube.zvtop;
-
-import java.util.Comparator;
-import java.util.Iterator;
-import java.util.Set;
-import java.util.SortedMap;
-import java.util.TreeMap;
-import java.util.Map.Entry;
-
-/**
- * Container of ZvMaxima
- * This enables for the ZvMaxima to be sorted by V(r)
- * while they still can be accessed by index (i, j)
- * @author jstrube
- * @version $Id: ZvMaximumMatrix.java,v 1.1 2006/10/09 17:25:27 jstrube Exp $
- *
- */
-class ZvMaximumMatrix implements Iterable {
-
- private class MatrixIndex {
- MatrixIndex(int i, int j) {
- _i = i;
- _j = j;
- }
- int _i;
- int _j;
- }
- SortedMap<ZvMaximum, MatrixIndex> _table;
-
- /**
- *
- */
- ZvMaximumMatrix(Comparator comp) {
- _table = new TreeMap<ZvMaximum, MatrixIndex>(comp);
- }
-
- void add(ZvMaximum max, int i, int j) {
- _table.put(max, new MatrixIndex(i, j));
- }
-
- /* (non-Javadoc)
- * @see java.lang.Iterable#iterator()
- */
- public Iterator iterator() {
- return _table.keySet().iterator();
- }
-
- /**
- * Accessor of ZvMaximum by indices of the constituent tracks in the list
- * @param i index of the first Track
- * @param j index of the second Track
- * @return the ZvMaximum
- */
- ZvMaximum get(int i, int j) throws IllegalAccessException {
- // Slow Access - OK, because primary access is by order in V(r)
- Set<Entry<ZvMaximum, MatrixIndex>> entries = _table.entrySet();
- for (Entry<ZvMaximum, MatrixIndex> entry : entries) {
- MatrixIndex index = entry.getValue();
- // Symmetry in indices
- if ((index._i == i && index._j == j)
- || (index._i == j && index._j == i) ) {
- return entry.getKey();
- }
- }
- throw new IllegalAccessException("ZvMaximumMatrix.get: index not found");
- }
-
- Set<ZvMaximum> getMaxima() {
- return _table.keySet();
- }
-
-}
lcsim/src/org/lcsim/contrib/JanStrube/zvtop
diff -N ZvTop.java
--- ZvTop.java 29 Nov 2007 02:25:56 -0000 1.3
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,501 +0,0 @@
-package org.lcsim.contrib.JanStrube.zvtop;
-
-import static org.lcsim.recon.vertexing.zvtop4.VectorArithmetic.add;
-import static org.lcsim.recon.vertexing.zvtop4.VectorArithmetic.dot;
-import static org.lcsim.recon.vertexing.zvtop4.VectorArithmetic.multiply;
-import static org.lcsim.recon.vertexing.zvtop4.VectorArithmetic.subtract;
-import static org.lcsim.recon.vertexing.zvtop4.VectorArithmetic.unit;
-import static org.lcsim.recon.vertexing.zvtop4.VectorArithmetic.angle;
-import static org.lcsim.contrib.JanStrube.zvtop.ZvUtil.twoTrackMax;
-
-import java.util.ArrayList;
-import java.util.Comparator;
-import java.util.HashMap;
-import java.util.HashSet;
-import java.util.List;
-import java.util.Map;
-import java.util.Set;
-import java.util.SortedSet;
-import java.util.TreeSet;
-
-import org.lcsim.contrib.JanStrube.vtxFitter.Fitter;
-import org.lcsim.contrib.JanStrube.vtxFitter.Vertex;
-import org.lcsim.event.Track;
-import org.lcsim.spacegeom.CartesianPoint;
-import org.lcsim.spacegeom.SpacePoint;
-import org.lcsim.spacegeom.SpacePointVector;
-import org.lcsim.spacegeom.SpaceVector;
-
-// TODO think about implementation in terms of static members
-/**
- * Main class to handle the vertexing Vertices are found by the function
- * @see #findVertices(List)
- * @author jstrube
- * @version $Id: ZvTop.java,v 1.3 2007/11/29 02:25:56 jstrube Exp $
- */
-public class ZvTop {
- /**
- * @author jstrube
- */
- class MaxComp implements Comparator {
-
- /**
- *
- */
- MaxComp() {
- super();
- // TODO Auto-generated constructor stub
- }
-
- /*
- * (non-Javadoc)
- * @see java.util.Comparator#compare(T, T)
- */
- public int compare(Object arg1, Object arg2) {
- if (arg1.getClass() != ZvMaximum.class
- && arg2.getClass() != ZvMaximum.class)
- throw new IllegalArgumentException(
- "Can only compare ZvMaximum Objects");
- ZvMaximum max1 = (ZvMaximum) arg1;
- ZvMaximum max2 = (ZvMaximum) arg2;
- double overlap1 = max1.getValue();
- double overlap2 = max2.getValue();
- if (overlap1 > overlap2)
- return 1;
- if (overlap1 < overlap2)
- return -1;
- return 0;
- }
-
- }
-
- /*
- * ZVTOP groups overlaps that cannot be resolved.
- * This class holds for each of those clusters the set of tracks associated with them
- * as well as the position of the overlap with the highest value
- */
- private class TrackCluster {
- public Set<Track> tracks;
- public SpacePoint position;
- public TrackCluster(Set<Track> t, SpacePoint p) {
- tracks = t;
- position = p;
- }
- }
-
- // if a two-track overlap is less than this cut, do not add max to the list
- private double _globalOverlapCut = 1E-10;
- private SpaceVector _jetAxis;
- private double _maxDistFromIP = 25;
- // List of vertex candidates.
- // The two-track maxima are transformed to the nearest location of a maximum
- // in the overlap function.
- // The highest of those maxima is taken as the location of the vertex and
- // all tracks belonging to ZvMaxima that are unresolved from this highest
- // maximum are added to the vertex.
- // At this stage tracks can be assigned to more than one vertex.
- // The fitVertices function takes care of assigning tracks uniquely to a
- // vertex.
- // TODO need to distinguish between pruning and fitting
- private List<TrackCluster> _vertexCandidateList;
- // Special weight for the angle of the spatial point (argument of overlap)
- // wrt. the jet axis
- private double angularWeight = 5;
- private double chiSquareCut = 15;
- // These data members must be final, because changing their values
- // has an effect on many of the other classes
- // Special weight for the IP in the overlap function
- private double ipWeight = 1;
- // resolution function is expensive in CPU
- private Map<ZvMaximum, Set<ZvMaximum>> isResolvedFromMap;
-
- // Number of points to check along the line from point 1 to point 2
- int iterationMax = 3;
-
- // stores the f_i(r)*f_j(r), which are later translated into the
- // V(r_ij)
- ZvMaximumMatrix maximumMatrix;
- private double resolvedCut = 0.7;
-
- private double sigmaL;
-
- private double sigmaT;
-
- List<Track> trackList;
-
- private Map<ZvMaximum, Set<ZvMaximum>> unresolvedMaximaMap;
-
- private Fitter vtxFitter;
-
- /**
- * Constructor Creates a list of ZvTracks from the argument list. Calculates
- * the overlap and resolution functions for the tracks
- * @param jetAxisValue
- */
- public ZvTop(SpaceVector jetAxisValue, Fitter f) {
- // FIXME Using arbitrary values for the measurement errors
- sigmaT = 0.01;
- sigmaL = 0.1;
- _jetAxis = jetAxisValue;
- maximumMatrix = new ZvMaximumMatrix(new MaxComp());
- isResolvedFromMap = new HashMap<ZvMaximum, Set<ZvMaximum>>();
- unresolvedMaximaMap = new HashMap<ZvMaximum, Set<ZvMaximum>>();
- trackList = new ArrayList<Track>();
- _vertexCandidateList = new ArrayList<TrackCluster>();
- vtxFitter = f;
- return;
- }
-
- /**
- * Calculates the resolution for two points in space. Two points in space
- * are resolved, if the maximum overlap along the line that connects them
- * divided by the greater of the overlap values at the two points is less
- * than a cutoff value R0
- * @param sp1 point 1
- * @param sp2 point 2
- * @return true, if the resolution function is less than the cutoff, false
- * otherwise
- */
- public boolean areResolved(SpacePoint sp1, SpacePoint sp2) {
- double overlap1 = overlap(sp1);
- double overlap2 = overlap(sp2);
- double minOverlap = 200;
- SpaceVector vec = subtract(sp2, sp1);
- for (int i = 1; i < iterationMax; ++i) {
- SpacePoint rVec = add(sp1, multiply(vec, 1.0*i / iterationMax));
- double overlap = overlap(rVec);
- if (overlap < minOverlap)
- minOverlap = overlap;
- }
- if (overlap1 < minOverlap)
- minOverlap = overlap1;
- if (overlap2 < minOverlap)
- minOverlap = overlap2;
-// System.out.println("Value: " + minOverlap
-// / ((overlap2 < overlap1) ? overlap2 : overlap1));
- return minOverlap / ((overlap2 < overlap1) ? overlap2 : overlap1) < resolvedCut;
- }
-
- /**
- * Calculates for each ZvMaximum whether of not it is resolved from other
- * ZvMaxima. Fills two maps to cache those values. Each ZvMaximum is mapped
- * to a set of resolved/unresolved ZvMaxima
- */
- private void assignResolution() {
-// System.out.println("assignResolution:");
- final Set<ZvMaximum> maxSet = maximumMatrix.getMaxima();
- for (ZvMaximum iMax : maxSet) {
- Set<ZvMaximum> isResolvedFromImax = new HashSet<ZvMaximum>();
- Set<ZvMaximum> isNotResolvedFromImax = new HashSet<ZvMaximum>();
- for (ZvMaximum jMax : maxSet) {
- if (iMax == jMax)
- continue;
-// System.out.printf("%s and %s ", iMax, jMax);
- if (areResolved(iMax.getLocation(), jMax.getLocation())) {
- isResolvedFromImax.add(jMax);
-// System.out.println("are resolved");
- } else {
- isNotResolvedFromImax.add(jMax);
-// System.out.println("are not resolved");
- }
- }
- isResolvedFromMap.put(iMax, isResolvedFromImax);
- unresolvedMaximaMap.put(iMax, isNotResolvedFromImax);
- }
-// System.out.printf("isResolvedFromMap: %d\tunresolvedMaximaMap: %d\n",
-// isResolvedFromMap.size(), unresolvedMaximaMap.size());
-// System.out.println("ZvMaxima");
-// for (ZvMaximum iMax : unresolvedMaximaMap.keySet()) {
-// System.out.println(iMax);
-// System.out.printf("Is resolved from %d other Maxima\n",
-// isResolvedFromMap.get(iMax).size());
-// System.out.printf("Is unresolved from %d other Maxima\n",
-// unresolvedMaximaMap.get(iMax).size());
-// }
-// System.out.println("done assignResolution");
- }
-
- /**
- * Creates a list of Vertex candidates. Each candidate consists of a number
- * of track-IP overlaps that are not spatially resolved. The assignment is
- * unique in the sense that each overlap is assigned to exactly one
- * candidate.
- */
- private void clusterCandidates() {
-// System.err.println("entering clusterCandidates");
- SortedSet<ZvMaximum> isAvailable = new TreeSet<ZvMaximum>(maximumMatrix
- .getMaxima());
- // each cluster consists of several ZvMaxima
- // add the maximum and the set of unresolved maxima to the list of
- // candidates
- // __is__ already sorted by overlap value
- while (!isAvailable.isEmpty()) {
- ZvMaximum highestRemaining = isAvailable.first();
- Set<ZvMaximum> unresolvedSet = unresolvedMaximaMap
- .get(highestRemaining);
- unresolvedSet.add(highestRemaining);
- Set<Track> tracks = new HashSet<Track>();
- for (ZvMaximum imax : unresolvedSet) {
- tracks.add(imax.getTrack_i());
- tracks.add(imax.getTrack_j());
- }
- _vertexCandidateList.add(new TrackCluster(tracks, highestRemaining.getLocation()));
- isAvailable.removeAll(unresolvedSet);
- }
-// System.err.printf("done clusterCandidates: %d Candidates",
-// _vertexCandidateList.size());
- return;
- }
-
- /**
- * Transforms each of the maxima in the two-track overlap to the nearest
- * global maximum in V(r). Finds the global maxima first.
- */
- // TODO maybe the maxima don't even have to be transformed
- // TODO the sequence of for-loops is a little clumsy
- private void findGlobalMaxima() {
-// System.out.println("Entering findGlobalMaxima");
- // search in micrometer steps in xy
- final double xyStep = 0.001;
- // search in 10 micrometer steps in z
- final double zStep = 0.01;
- final int maxSteps = 500;
- for (ZvMaximum iMax : maximumMatrix.getMaxima()) {
- SpacePoint searchPoint = iMax.getLocation();
- SpacePoint maxPoint = new SpacePoint(searchPoint);
- double xPos = searchPoint.x();
- double yPos = searchPoint.y();
- double zPos = searchPoint.z();
- double globalMax = iMax.getValue();
- // perform the iteration __consecutively__, not simultaneously !
- // for every direction: go in positive direction until maximum is
- // found,
- // reset to initial position and go in negative direction
- // +x
- for (int x = 1; x <= maxSteps; ++x) {
- xPos += xyStep;
- // TODO this is maybe where a double[] is acceptable
- searchPoint = new CartesianPoint(xPos, yPos, zPos);
- double overlap = overlap(searchPoint);
- if (overlap > globalMax) {
- globalMax = overlap;
- maxPoint = new SpacePoint(searchPoint);
- } else {
- xPos -= xyStep;
- break;
- }
- }
- // reset to original pos
- xPos = iMax.getLocation().x();
- // -x
- for (int x = 1; x <= maxSteps; ++x) {
- xPos -= xyStep;
- searchPoint = new CartesianPoint(xPos, yPos, zPos);
- double overlap = overlap(searchPoint);
- if (overlap > globalMax) {
- globalMax = overlap;
- maxPoint = new SpacePoint(searchPoint);
- } else {
- xPos += xyStep;
- break;
- }
- }
- // search in x is done. set to found maximum
- xPos = maxPoint.x();
- // +y
- for (int y = 1; y <= maxSteps; ++y) {
- yPos += xyStep;
- searchPoint = new CartesianPoint(xPos, yPos, zPos);
- double overlap = overlap(searchPoint);
- if (overlap > globalMax) {
- globalMax = overlap;
- maxPoint = new SpacePoint(searchPoint);
- } else {
- yPos -= xyStep;
- break;
- }
- }
- // reset to original pos
- yPos = iMax.getLocation().y();
- // -y
- for (int y = 1; y <= maxSteps; ++y) {
- yPos -= xyStep;
- searchPoint = new CartesianPoint(xPos, yPos, zPos);
- double overlap = overlap(searchPoint);
- if (overlap > globalMax) {
- globalMax = overlap;
- maxPoint = new SpacePoint(searchPoint);
- } else {
- yPos += xyStep;
- break;
- }
- }
- // search in y is done. set to found maximum
- yPos = maxPoint.y();
- // +z
- for (int z = 1; z <= maxSteps; ++z) {
- zPos += zStep;
- searchPoint = new CartesianPoint(xPos, yPos, zPos);
- double overlap = overlap(searchPoint);
- if (overlap > globalMax) {
- globalMax = overlap;
- maxPoint = new SpacePoint(searchPoint);
- } else {
- zPos -= zStep;
- break;
- }
- }
- // reset to original pos
- zPos = iMax.getLocation().z();
- // -z
- for (int z = 1; z <= maxSteps; ++z) {
- zPos -= zStep;
- searchPoint = new CartesianPoint(xPos, yPos, zPos);
- double overlap = overlap(searchPoint);
- if (overlap > globalMax) {
- globalMax = overlap;
- maxPoint = new SpacePoint(searchPoint);
- } else {
- zPos += zStep;
- break;
- }
- }
- // search in x,y,z is done. set to found maximum
- iMax.setLocation(maxPoint);
- }
-// System.out.println("Global Maxima done");
- }
-
- /**
- * Finds the Maxima
- * @param list List of tracks
- * @return The list of vertices
- */
- public List<Vertex> findVertices(List<Track> list) {
- trackList.clear();
-// System.err
-// .printf("entering findVertives with %d Tracks\n", list.size());
- // FIXME Put in IP as constraint. At the moment, no vertex contains IP
- for (Track iTrack : list) {
- trackList.add(iTrack);
- }
- int size = trackList.size();
- for (int i = 0; i < size; ++i) {
- Track iTrack = trackList.get(i);
- for (int j = i + 1; j < size; ++j) {
- Track jTrack = trackList.get(j);
- SpacePoint location = twoTrackMax(iTrack, jTrack);
- double overlap = overlap(location);
- // consider only significant overlap
- if (overlap < _globalOverlapCut)
- continue;
- ZvMaximum max = new ZvMaximum(iTrack, jTrack, location, overlap);
- maximumMatrix.add(max, i, j);
- }
- }
-// System.out.println("maximumMatrix: " + maximumMatrix._table.size());
- findGlobalMaxima();
- assignResolution();
- clusterCandidates();
-// System.err.println("done findVertices\n");
- return fitVertices();
- }
-
- /**
- * fits the tracks to the vertices. A ZvVertex can only have tracks with a
- * chi2 contribution below a certain cut value. They are assigned from a
- * sorted set of maxima, so that the tracks are always assigned to the
- * vertex with the highest value of V(r), if possible
- * @return A List of Vertices
- */
- private List<Vertex> fitVertices() {
- List<Vertex> l = new ArrayList<Vertex>();
- Set<Track> unavailableTracks = new HashSet<Track>();
-// System.err.printf("found %d Candidates\n", _vertexCandidateList.size());
-// for (Set<Track> iCand : _vertexCandidateList) {
-// System.err.printf(iCand.toString());
-// }
- // The Maxima are clustered according to the resolution function
- // A Cluster is a Set of ZvMaxima
- // Add all tracks from a cluster to the vertex
- // The vertex decides based on the chiSquared value
- // which of those tracks it wants to keep
- for (TrackCluster iVtx : _vertexCandidateList) {
-// System.out.println("using chi2 cut of " + chiSquareCut);
- for (Track t : unavailableTracks) {
- iVtx.tracks.remove(t);
- }
- // FIXME using a vertex as the container isn't so good
- // because a) tracks can drop out
- // so that a newly created vertex should be returned...
- // that makes for messy bookkeeping
- Vertex v = vtxFitter.fit(iVtx.tracks, iVtx.position);
- l.add(v);
- unavailableTracks.addAll(v.getTracks());
- }
- return l;
- }
-
- public double getGlobalOverlapCut() {
- return _globalOverlapCut;
- }
-
- /**
- * Calculates the global overlap at a given point in space The global
- * overlap is the sum S of all track tubes minus the sum of the squares of
- * the track tubes divided by S
- * @param location The point in space at which to evaluate the function
- * @return The value at that point
- */
- public double overlap(SpacePoint location) {
- double numerator = 0;
- double denominator = 0;
- // calculate the angular weight
- // the angle is taken not from (000),
- // but from a point -0.1 mm along the jet axis
- double[] origin = multiply(unit(_jetAxis), -0.1).v();
- SpacePointVector vec2Location = new SpacePointVector(
- new CartesianPoint(origin[0], origin[1], origin[2]), location);
- double angle = angle(_jetAxis, vec2Location);
- double angularFactor = Math.exp(-angularWeight * angle * angle);
- if (dot(unit(_jetAxis), vec2Location.getDirection()) > _maxDistFromIP) {
- return 0;
- }
- // FIXME dealing with IP
- double ipFactor = 0; // ipWeight *
- // trackList.get(0).getTubeValue(location);
- denominator += ipFactor;
- numerator += 0; // ipFactor * trackList.get(0).getTubeValue(location);
- for (int i = 0; i < trackList.size(); ++i) {
- Track iTrack = trackList.get(i);
- double iTubeVal = ZvTrack.getProbability(iTrack, location);
- numerator += iTubeVal * iTubeVal;
- denominator += iTubeVal;
- }
- return angularFactor * (denominator - numerator / denominator);
- }
-
- public void setAngularWeight(double w) {
- angularWeight = w;
- }
-
- public void setChiSquareCut(double c) {
- chiSquareCut = c;
- }
-
- public void setGlobalOverlapCut(double globalOverlapCut) {
- _globalOverlapCut = globalOverlapCut;
- }
-
- public void setIpWeight(double w) {
- ipWeight = w;
- }
-
- public void setMaxDistFromIP(double d) {
- _maxDistFromIP = d;
- }
-
- public void setResolvedCut(double c) {
- resolvedCut = c;
- }
-}
lcsim/src/org/lcsim/contrib/JanStrube/zvtop
diff -N ZvTrack.java
--- ZvTrack.java 29 Nov 2007 21:29:35 -0000 1.4
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,24 +0,0 @@
-/**
- * @version $Id: ZvTrack.java,v 1.4 2007/11/29 21:29:35 jstrube Exp $
- */
-package org.lcsim.contrib.JanStrube.zvtop;
-
-import org.lcsim.event.Track;
-import org.lcsim.mc.fast.tracking.fix.FastMCTrack;
-import org.lcsim.spacegeom.SpacePoint;
-
-/**
- * @author jstrube
- *
- */
-public class ZvTrack {
- /**
- * returns the probability that a track goes through a certain point
- * Assuming gaussian errors.
- * After Lyons, p. 47
- */
- public static double getProbability(Track t, SpacePoint l) {
- throw new RuntimeException("Not implemented");
-// return 0;
- }
-}
lcsim/src/org/lcsim/contrib/JanStrube/zvtop
diff -N ZvUtil.java
--- ZvUtil.java 29 Nov 2007 02:25:56 -0000 1.3
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,119 +0,0 @@
-package org.lcsim.contrib.JanStrube.zvtop;
-
-import static org.lcsim.recon.vertexing.zvtop4.VectorArithmetic.unit;
-
-import java.util.Collection;
-
-import org.lcsim.event.LCIOParameters.ParameterName;
-import org.lcsim.event.Track;
-import org.lcsim.spacegeom.CartesianVector;
-import org.lcsim.spacegeom.SpacePoint;
-import org.lcsim.spacegeom.SpaceVector;
-import org.lcsim.util.swim.HelixSwimmer;
-
-/**
- * Repository of utility functions
- * @author jstrube
- * @version $Id: ZvUtil.java,v 1.3 2007/11/29 02:25:56 jstrube Exp $
- */
-final class ZvUtil {
-
- /**
- * Private constructor ensures the class is never instantiated
- */
- private ZvUtil() {
- }
-
- /**
- * Calculates the location where the product of the two "track probabilities"
- * is at a maximum
- * @param iTrack
- * @param jTrack
- * @return the location of the maximum
- */
- public static SpacePoint twoTrackMax(Track iTrack, Track jTrack) {
- return swimmerTwoTrackMax(iTrack, jTrack);
- }
-
- /*
- * calculates the maximum of the two track overlap using the swimmer
- * The error is computed using FIXME
- */
- private static SpacePoint swimmerTwoTrackMax(Track iTrack, Track jTrack) {
- HelixSwimmer swimmer = new HelixSwimmer(5);
- swimmer.setTrack(iTrack);
-// swimmer.getTrackLengthToPoint(jTrack.position());
- throw new RuntimeException("Implement me !");
-// return new SpacePoint();
- }
-
-
- /**
- * Sums the momenta of all tracks in the list
- * @param list
- * @return the vector sum of the momenta
- */
- public static SpaceVector sumTrackMomenta(Collection<Track> list) {
- double x = 0.;
- double y = 0.;
- double z = 0.;
- for (Track iTrack: list) {
- double[] p3 = iTrack.getMomentum();
- x += p3[0];
- y += p3[1];
- z += p3[2];
- }
- return new CartesianVector(x, y, z);
- }
-
- /**
- * Sums the charge of all tracks in the list
- * @param list
- * @return the sum of track charges
- */
- public static double sumTrackCharge(Collection<Track> list) {
- double result = 0;
- for (Track iTrack : list) {
- result += iTrack.getCharge();
- }
- return result;
- }
-
-
- /**
- * Projection of error matrix along specific direction <br>
- * (Does normalization of direction vector.)
- * @param A error matrix to be projected (3x3)
- * @param direction direction to project along (3-dim)
- * @return square of sigma along direction
- */
- /*
- * TODO Legacy code
- */
- public static double projectSigma2(double[][]A, SpaceVector direction)
- {
- SpaceVector dir = unit(direction);
- double[] d = new double[] {dir.x(), dir.y(), dir.z()};
- double sig2 = 0;
- for (int i=0; i<3; i++) {
- for (int j=0; j<3; j++)
- sig2 += d[i]*A[i][j]*d[j];
- }
- return sig2;
- }
-
- /**
- * calculate p/tangent unit vector
- */
- public static SpaceVector getUnitTangent(double[] hlxPar)
- {
- int tanLambda = ParameterName.tanLambda.ordinal();
- int phi = ParameterName.phi0.ordinal();
- double norm = Math.sqrt(1.+hlxPar[tanLambda]*hlxPar[tanLambda]);
- double x = Math.cos(hlxPar[phi])/norm;
- double y = Math.sin(hlxPar[phi])/norm;
- double z = hlxPar[tanLambda]/norm;
- return new CartesianVector(x, y, z);
- }
-
-}
CVSspam 0.2.8