12 removed files
lcsim/src/org/lcsim/contrib/timb/mc/fast
diff -N MCFast.java
--- MCFast.java 26 May 2006 07:21:47 -0000 1.1
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,32 +0,0 @@
-package org.lcsim.mc.fast;
-
-import org.lcsim.mc.fast.cluster.ronan.MCFastRonan;
-import org.lcsim.mc.fast.reconstructedparticle.MCFastReconstructedParticleDriver;
-import org.lcsim.mc.fast.tracking.MCFastTracking;
-import org.lcsim.util.Driver;
-
-/**
- *
- * @author Tony Johnson
- */
-public class MCFast extends Driver
-{
- /** Creates a new instance of MCFast */
- public MCFast(boolean beamSpotConstraint, boolean simple)
- {
- add(new MCFastTracking(beamSpotConstraint, simple));
- add(new MCFastRonan());
- add(new MCFastReconstructedParticleDriver());
- }
-
- public MCFast(boolean beamSpotConstraint, boolean simple, long seed)
- {
- this(beamSpotConstraint, simple);
- getRandom().setSeed(seed);
- }
-
- public MCFast()
- {
- this(false,false);
- }
-}
lcsim/src/org/lcsim/contrib/timb/mc/fast/cluster/ronan
diff -N ClusterResolutionTables.java
--- ClusterResolutionTables.java 26 May 2006 07:22:04 -0000 1.1
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,153 +0,0 @@
-package org.lcsim.mc.fast.cluster.ronan;
-
-import org.lcsim.conditions.ConditionsSet;
-import static java.lang.Math.sqrt;
-import static java.lang.Math.pow;
-
-public class ClusterResolutionTables
-{
- private boolean JETParameterization;
- private double JETResolution;
- private double JETHadDegradeFraction;
- private double JETEMEnergyFraction;
- private double JETHadEnergyFraction;
- private double Lambda_j;
-
- private double EMAlignmentError;
- private double EMConstantTerm;
- private double EMPositionError;
- private double EMResolution;
- private double EMOnset;
- private double EMSharpness;
-
- private double HADAlignmentError;
- private double HADConstantTerm;
- private double HADPositionError;
- private double HADResolution;
- private double HADOnset;
- private double HADSharpness;
-
- private double PolarEMInner;
- private double PolarEMOuter;
- private double PolarHADInner;
- private double PolarHADOuter;
-
- ClusterResolutionTables(ConditionsSet set)
- {
- JETParameterization = Boolean.parseBoolean(set.getString("JETParameterization"));
- JETResolution = set.getDouble("JETResolution");
- JETHadDegradeFraction = set.getDouble("JETHadDegradeFraction");
- JETEMEnergyFraction = set.getDouble("JETEMEnergyFraction");
- JETHadEnergyFraction = set.getDouble("JETHadEnergyFraction");
-
- EMOnset = set.getDouble("EMOnset");
- EMSharpness = set.getDouble("EMSharpness");
- PolarEMInner = set.getDouble("PolarEMInner");
- PolarEMOuter = set.getDouble("PolarEMOuter");
-
- EMResolution = set.getDouble("EMResolution");
- EMConstantTerm = set.getDouble("EMConstantTerm");
- EMPositionError = set.getDouble("EMPositionError");
- EMAlignmentError = set.getDouble("EMAlignmentError");
-
- HADOnset = set.getDouble("HADOnset");
- HADSharpness = set.getDouble("HADSharpness");
- PolarHADInner = set.getDouble("PolarHADInner");
- PolarHADOuter = set.getDouble("PolarHADOuter");
-
- HADResolution = set.getDouble("HADResolution");
- HADConstantTerm = set.getDouble("HADConstantTerm");
- HADPositionError = set.getDouble("HADPositionError");
- HADAlignmentError = set.getDouble("HADAlignmentError");
- if (JETParameterization)
- {
- EMConstantTerm=0.;
- HADConstantTerm=0.;
- Lambda_j=(pow(JETResolution,2)-JETEMEnergyFraction*pow(EMResolution,2)-JETHadEnergyFraction*pow(HADResolution,2))
- /((1.-JETHadDegradeFraction)*JETEMEnergyFraction*pow(EMResolution,2)+JETHadDegradeFraction*JETHadEnergyFraction*pow(HADResolution,2));
- EMResolution*=sqrt(1.+Lambda_j*(1.-JETHadDegradeFraction));
- HADResolution*=sqrt(1.+Lambda_j*JETHadDegradeFraction);
- System.out.println(" JETParameterization settings Lamda_j= "+Lambda_j+" EMResolution= "+EMResolution+" HADResolution= "+HADResolution);
- }
- }
-
- public double getEMAlignmentError()
- {
- return EMAlignmentError;
- }
-
- public double getEMConstantTerm()
- {
- return EMConstantTerm;
- }
-
- public double getEMPositionError()
- {
- return EMPositionError;
- }
-
- public double getEMResolution()
- {
- return EMResolution;
- }
-
- public double getEMOnset()
- {
- return EMOnset;
- }
-
- public double getEMSharpness()
- {
- return EMSharpness;
- }
-
- public double getHADAlignmentError()
- {
- return HADAlignmentError;
- }
-
- public double getHADConstantTerm()
- {
- return HADConstantTerm;
- }
-
- public double getHADPositionError()
- {
- return HADPositionError;
- }
-
- public double getHADResolution()
- {
- return HADResolution;
- }
-
- public double getHADOnset()
- {
- return HADOnset;
- }
-
- public double getHADSharpness()
- {
- return HADSharpness;
- }
-
- public double getPolarEMInner()
- {
- return PolarEMInner;
- }
-
- public double getPolarEMOuter()
- {
- return PolarEMOuter;
- }
-
- public double getPolarHADInner()
- {
- return PolarHADInner;
- }
-
- public double getPolarHADOuter()
- {
- return PolarHADOuter;
- }
-}
lcsim/src/org/lcsim/contrib/timb/mc/fast/cluster/ronan
diff -N MCFastRonan.java
--- MCFastRonan.java 26 May 2006 07:22:06 -0000 1.1
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,149 +0,0 @@
-package org.lcsim.mc.fast.cluster.ronan;
-
-import java.util.ArrayList;
-import java.util.List;
-import java.util.Random;
-
-import org.lcsim.conditions.ConditionsEvent;
-import org.lcsim.conditions.ConditionsListener;
-import org.lcsim.conditions.ConditionsSet;
-import org.lcsim.event.Cluster;
-import org.lcsim.event.EventHeader;
-import org.lcsim.event.MCParticle;
-import org.lcsim.util.Driver;
-
-
-
-
-/**
- * Fast Monte Carlo cluster simulator
- * @author M.Ronan Oct 2000 - Added "refined" cluster simulation
- * @version
- */
-public class MCFastRonan extends Driver implements ConditionsListener
-{
- private final static int ElecID = 11;
- private final static int NuEID = 12;
- private final static int MuID = 13;
- private final static int NuMuID = 14;
- private final static int NuTauID = 16;
- private final static int PhotonID = 22;
- private final static int Neutralino1 = 1000022;
- private final static int Neutralino2 = 1000023;
- private final static int Neutralino3 = 1000025;
- private final static int Neutralino4 = 1000035;
- private static long iprint=0;
-
-
- private ClusterResolutionTables clusterParm;
-
- protected void process(EventHeader event)
- {
- if (clusterParm == null)
- {
- ConditionsSet conditions = getConditionsManager().getConditions("ClusterParameters");
- conditions.addConditionsListener(this);
- clusterParm = new ClusterResolutionTables(conditions);
- }
- List<Cluster> cl = new ArrayList<Cluster>();
-
- boolean hist = getHistogramLevel() > 0;
-
- List<MCParticle> particles = event.getMCParticles();
- for (MCParticle p : particles)
- {
-
- // filter for FINALSTATE
- if (p.getGeneratorStatus() != p.FINAL_STATE)
- {
- continue;
- }
-
- int PDGID = p.getPDGID();
- int absPDGID = Math.abs(PDGID);
- double charge = p.getCharge();
-
-
- // filter neutrinos
- boolean neutrino = absPDGID == NuEID || absPDGID == NuMuID || absPDGID == NuTauID
- || absPDGID == Neutralino1 || absPDGID == Neutralino2 || absPDGID == Neutralino3 || absPDGID == Neutralino4;
- if (neutrino)
- {
- continue;
- }
-
- double E = p.getEnergy();
- double pt2 = p.getMomentum().magnitudeSquared()-p.getPZ()*p.getPZ();
- double pt = Math.sqrt(pt2);
- double ptot = p.getMomentum().magnitude();
- double cosTheta = p.getPZ() / p.getMomentum().magnitude();
-
- Random rand = getRandom();
-
- // Photons
- if (PDGID == PhotonID || PDGID == ElecID)
- {
- // within acceptance
- // double thing = (1 - 1 / ( 1 + Math.exp( (E-clusterParm.getEMOnset())*clusterParm.getEMSharpness() ) ));
- //if (rand.nextDouble() > thing)
- //{
- // continue;
- //}
- if (E < clusterParm.getEMOnset())
- {
- continue;
- }
- if (Math.abs(cosTheta) > clusterParm.getPolarEMOuter())
- {
- continue;
- }
-
- cl.add(new ReconEMCluster(clusterParm, rand, p, hist));
-
- }
-
- // Neutral hadrons
- else if (PDGID != MuID)
- {
- // within acceptance
-
- //double thing = (1 - 1 / ( 1 + Math.exp( (E-clusterParm.getHADOnset())*clusterParm.getHADSharpness() ) ));
- //if (rand.nextDouble() > thing)
- //{
- //continue;
- //}
- if (E < clusterParm.getHADOnset())
- {
- continue;
- }
- if (Math.abs(cosTheta) > clusterParm.getPolarHADOuter())
- {
- continue;
- }
-
- cl.add(new ReconHADCluster(clusterParm, rand, p, hist));
- }
- }
- double neg_energy_total = 0.;
- double pos_energy_weight_total = 0.;
- for (Cluster rcl : cl )
- {
- if( Math.abs(((ReconCluster)rcl).getMCParticle().getCharge()) > Double.MIN_VALUE ) continue ;
- neg_energy_total+=((ReconCluster)rcl).getNegEnergy();
- pos_energy_weight_total+=((ReconCluster)rcl).getNegEnergy() < 0. ? 0. : Math.min(((ReconCluster)rcl).getSigma(),((ReconCluster)rcl).getEnergy());
- }
- iprint++;
- if(iprint < -20) System.out.println(" MCFast neg_energy_total= "+neg_energy_total+" pos_energy_weight_total= "+pos_energy_weight_total);
- if( neg_energy_total < -Double.MIN_VALUE ) for (Cluster rcl : cl ) if( Math.abs(((ReconCluster)rcl).getMCParticle().getCharge()) < Double.MIN_VALUE && ((ReconCluster)rcl).getNegEnergy() >= 0. )
- ((ReconCluster)rcl).adjustEnergy(neg_energy_total,pos_energy_weight_total);
- event.put(EventHeader.CLUSTERS, cl, Cluster.class, 0);
- }
-
-
- public void conditionsChanged(ConditionsEvent event)
- {
- ConditionsSet conditions = getConditionsManager().getConditions("ClusterParameters");
- clusterParm = new ClusterResolutionTables(conditions);
- }
-}
-
lcsim/src/org/lcsim/contrib/timb/mc/fast/cluster/ronan
diff -N ReconCluster.java
--- ReconCluster.java 26 May 2006 07:22:08 -0000 1.1
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,204 +0,0 @@
-package org.lcsim.mc.fast.cluster.ronan;
-
-import hep.physics.particle.Particle;
-import java.util.Collections;
-import java.util.List;
-import java.util.Random;
-import org.lcsim.event.Cluster;
-
-abstract class ReconCluster implements Cluster
-{
- protected ClusterResolutionTables parm;
- protected Particle mcp;
- protected double a = 0;
- protected double b = 0;
- protected double c = 0;
- protected double d = 0;
- protected double energy;
- protected double neg_energy;
- protected double sigma;
- protected double phi;
- protected double radius;
- protected double theta;
- protected double transDist;
- private static long iprint=0;
-
- ReconCluster(ClusterResolutionTables parm, Random rand, Particle mcp, boolean hist)
- {
- this.parm = parm;
- this.mcp = mcp;
- }
-
- /** Best estimate for total energy of cluster */
- public double getEnergy()
- {
- return energy;
- }
-
- public double getNegEnergy()
- {
- return neg_energy;
- }
-
- public double getSigma()
- {
- return sigma;
- }
-
- public void adjustEnergy(double neg_energy_total, double pos_energy_weight_total)
- {
- iprint++;
- if(iprint < -200) System.out.println(" min(sigma,energy)="+Math.min(sigma,energy)+" ratio= "+(Math.min(sigma,energy)/pos_energy_weight_total)+" before adjust energy= "+energy);
- energy+=neg_energy_total*Math.min(sigma,energy)/pos_energy_weight_total;
-
- if (energy <= mcp.getMass()) energy=mcp.getMass()+Double.MIN_VALUE;
-
- if(iprint < -200) System.out.println(" neg_energy_total= "+neg_energy_total+" after adjust energy= "+energy);
- }
-
- protected void smear(Random rand, boolean hist)
- {
- // Get true energy from MCParticle
- double E = mcp.getEnergy();
-
- // Smear reconstructed energy
-
- smearEnergy(rand, E, hist);
-
- // Smear reconstructed position
- smearPosition(rand, E, hist);
- }
-
- protected void smearEnergy(Random rand, double E, boolean hist)
- {
- sigma = ((a / Math.sqrt(E)) + b) * E;
-
- energy = E + (sigma * rand.nextGaussian());
- if (energy <= mcp.getMass())
- {
- neg_energy=energy-mcp.getMass();
- energy=mcp.getMass()+Double.MIN_VALUE;
- }
- else
- {
- neg_energy=0.;
- }
- }
-
- protected void smearPosition(Random rand)
- {
- // Get true direction from MCParticle
- double Px = mcp.getPX();
- double Py = mcp.getPY();
- double Pz = mcp.getPZ();
-
- double P = Math.sqrt((Px * Px) + (Py * Py) + (Pz * Pz));
- double Phi = Math.atan2(Py, Px);
- if (Phi < 0)
- {
- Phi += (2 * Math.PI);
- }
-
- double Theta = Math.acos(Pz / P);
-
- // Simulate position smearing on a sphere of radius 2 meters
- radius = 2000.0;
-
- double x = (radius * Px) / P;
- double y = (radius * Py) / P;
- double z = (radius * Pz) / P;
-
- //these vectors vt and vs (orthonorm.) span a plane perpendicular to the momentum vector,
- //so smearing with a transdist will involve a lin. comb. of these
- double[] vt = {-Math.cos(Theta)*Math.cos(Phi),-Math.cos(Theta)*Math.sin(Phi),Math.sin(Theta)};
- double[] vs = {Math.sin(Phi),-Math.cos(Phi),0};
-
- //restricted to [0,PI] since transdist can be negative
- double alpha = rand.nextDouble() * Math.PI;
- x = x + transDist * ( Math.cos(alpha) * vt[0] + Math.sin(alpha) * vs[0] );
- y = y + transDist * ( Math.cos(alpha) * vt[1] + Math.sin(alpha) * vs[1] );
- z = z + transDist * ( Math.cos(alpha) * vt[2] + Math.sin(alpha) * vs[2] );
-
- phi = Math.atan2(y, x);
- if (phi < 0)
- {
- phi += (2 * Math.PI);
- }
- theta = Math.acos(z / radius);
- }
-
- public Particle getMCParticle()
- {
- return mcp;
- }
-
- abstract void smearPosition(Random rand, double E, boolean hist);
-
- public double[] getHitContributions()
- {
- return null;
- }
-
- public List getClusters()
- {
- return Collections.EMPTY_LIST;
- }
-
- public double[] getSubdetectorEnergies()
- {
- return null;
- }
-
- public double[] getPositionError()
- {
- return null; // fixme:
- }
-
- public int getType()
- {
- return 0; // Fixme:
- }
-
- public double getITheta()
- {
- return 0; // Fixme:
- }
-
- public double getIPhi()
- {
- return 0; // Fixme:
- }
-
- public double[] getDirectionError()
- {
- return null; // Fixme:
- }
-
- public List getCalorimeterHits()
- {
- return Collections.EMPTY_LIST;
- }
-
- public double[] getShape()
- {
- return null;
- }
-
- public double[] getPosition()
- {
- double x = radius * Math.sin(theta) * Math.cos(phi);
- double y = radius * Math.sin(theta) * Math.sin(phi);
- double z = radius * Math.cos(theta);
- return new double[] { x,y,z };
- }
-
- public double[] getParticleType()
- {
- return null; // Fixme:
- }
-
- public int getSize()
- {
- return 0;
- }
-}
lcsim/src/org/lcsim/contrib/timb/mc/fast/reconstructedparticle
diff -N IDResolutionTables.java
--- IDResolutionTables.java 26 May 2006 07:21:49 -0000 1.1
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,68 +0,0 @@
-/*
- * IDResolutionTables.java
- *
- * Created on July 6, 2005, 7:13 PM
- *
- * To change this template, choose Tools | Options and locate the template under
- * the Source Creation and Management node. Right-click the template and choose
- * Open. You can then make changes to the template in the Source Editor.
- */
-
-package org.lcsim.mc.fast.reconstructedparticle;
-
-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/timb/mc/fast/reconstructedparticle
diff -N MCFastReconstructedParticle.java
--- MCFastReconstructedParticle.java 26 May 2006 07:21:50 -0000 1.1
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,225 +0,0 @@
-package org.lcsim.mc.fast.reconstructedparticle;
-
-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 org.lcsim.mc.fast.tracking.ReconTrack;
-
-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 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();
- private static long iprint=0;
-
- public MCFastReconstructedParticle(Track t, ParticleType type, Particle p, Cluster assoc_c, double wtcal)
- {
- iprint++;
- if(iprint < -200) System.out.println(" PDGID= "+type.getPDGID()+" t.getPX,...= "+t.getPX()+" "+t.getPY()+" "+t.getPZ());
- _mass = type.getMass();
- addTrack(t);
- _charge = t.getCharge();
- // Use true point of origin for reference point for now.
- _referencePoint = p.getOrigin();
- e_track = sqrt(((ReconTrack)t).getDocaMomentumVec(_referencePoint).magnitudeSquared()+ _mass*_mass);
- p_track.setV3(e_track, ((ReconTrack)t).getDocaMomentumVec(_referencePoint));
- p3_track = p_track.v3();
- if (iprint < -200)
- {
- if (assoc_c != null)
- {
- System.out.println(" PDGID= "+type.getPDGID()+" e_track= "+e_track+" e_assoc_clus= "+assoc_c.getEnergy());
- System.out.println(" PDGID= "+type.getPDGID()+" _referencePoint= "+_referencePoint.x()+" "+_referencePoint.y()+" "+_referencePoint.z());
- System.out.println(" PDGID= "+type.getPDGID()+" p3_track= "+p3_track.x()+" "+p3_track.y()+" "+p3_track.z());
- }
- else
- {
- System.out.println(" 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));
- }
-
- public MCFastReconstructedParticle(Cluster c, ParticleType type, Particle p)
- {
- _mass = type.getMass();
- addCluster(c);
- double e = c.getEnergy();
- 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)
- {
- iprint++;
- _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);
- if (iprint < -200) {
- System.out.println(" PDGID= "+type.getPDGID()+" e_reco= "+e_reco+" mass= "+_mass);
- System.out.println(" PDGID= "+type.getPDGID()+" _referencePoint= "+_referencePoint.x()+" "+_referencePoint.y()+" "+_referencePoint.z());
- System.out.println(" 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 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();
- }
-}
lcsim/src/org/lcsim/contrib/timb/mc/fast/reconstructedparticle
diff -N MCFastReconstructedParticleDriver.java
--- MCFastReconstructedParticleDriver.java 26 May 2006 07:21:53 -0000 1.1
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,219 +0,0 @@
-package org.lcsim.mc.fast.reconstructedparticle;
-
-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.mc.fast.tracking.ReconTrack;
-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.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(((ReconTrack)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(((ReconTrack)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 ReconTrack)
- {
- ReconTrack rt = (ReconTrack) 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(abs(p.getCharge()) > 0.) 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();
- if(abs(p.getCharge()) > 0.) continue;
- int pdgid = p.getPDGID();
- if (hist)
- {
- aida.histogram1D("hadronCLS-particle", 150, -3, 3).fill((emc.getEnergy()-emc.getMCParticle().getEnergy())/(Math.sqrt(emc.getMCParticle().getEnergy())));
- }
-
- if ((abs(pdgid)==2112) && (rand.nextDouble() < IDEff.getNeutronEff()))
- {
- type = neutron;
- }
- else
- {
- type = klong;
- }
- }
- 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/timb/mc/fast/tracking
diff -N DocaTrackParameters.java
--- DocaTrackParameters.java 26 May 2006 07:21:55 -0000 1.1
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,759 +0,0 @@
-package org.lcsim.mc.fast.tracking;
-
-import Jama.Matrix;
-
-import hep.physics.particle.Particle;
-
-import hep.physics.vec.BasicHep3Vector;
-import hep.physics.vec.Hep3Vector;
-import hep.physics.vec.VecOp;
-import org.lcsim.constants.Constants;
-
-
-/**
- * Holds DOCA parameters and error matrix of track. Can be initialized
- * with a MC truth particle. <br>
- *
- * @author Tony Johnson, Wolfgang Walkowiak
- * @version $Id: DocaTrackParameters.java,v 1.1 2006/05/26 07:21:55 timb Exp $
- */
-class DocaTrackParameters implements TrackParameters
-{
- private Hep3Vector m_pdoca_ref = null;
- private Hep3Vector m_xdoca_ref = null;
- private Hep3Vector m_xref = null;
- private double[][] m_err = new double[5][5];
- private double[] m_parm = new double[5];
- private double m_Bz = 0.;
- private double m_chi2 = -1.;
- private double m_l0 = 0.;
- private double m_l_ref = 0.;
- private int m_ndf = 5;
-
- //====================================================
- //
- // Constructors
- //
- //====================================================
- DocaTrackParameters(double bField)
- {
- reset();
- m_Bz = bField;
- }
-
- DocaTrackParameters(double bField, Particle p)
- {
- this(bField, p.getMomentum(), p.getOrigin(), p.getType().getCharge());
- }
-
- DocaTrackParameters(double bField, Hep3Vector momentum, Hep3Vector x, double q)
- {
- reset();
- m_Bz = bField;
- calculateDoca(momentum, x, q);
- }
-
- DocaTrackParameters(double bField, double[] momentum, double[] x, double q)
- {
- reset();
- m_Bz = bField;
- calculateDoca(momentum, x, q);
- }
-
- DocaTrackParameters(double bField, double[] momentum, double[] x, double q, double[][] errorMatrix)
- {
- this(bField, momentum, x, q);
- fillErrorMatrix(errorMatrix);
- }
-
- DocaTrackParameters(double bField, double[] parameters)
- {
- m_Bz = bField;
- m_parm = parameters;
- }
-
- DocaTrackParameters(double bField, double[] parameters, double[][] errorMatrix)
- {
- this(bField, parameters);
- fillErrorMatrix(errorMatrix);
- }
-
- DocaTrackParameters(double bField, double[] parameters, double[][] errorMatrix, double chi2)
- {
- this(bField, parameters);
- fillErrorMatrix(errorMatrix);
- setChi2(chi2);
- }
-
- DocaTrackParameters(double bField, double[] parameters, double[][] errorMatrix, double chi2, int ndf)
- {
- this(bField, parameters);
- fillErrorMatrix(errorMatrix);
- setChi2(chi2);
- setNDF(ndf);
- }
-
- public double getD0()
- {
- return m_parm[0];
- }
-
- /**
- * Get the error matrix as a 2-D array
- * @see #getTrackParameter
- */
- public double[][] getErrorMatrix()
- {
- return m_err;
- }
-
- /**
- * Get an individual element of the error matrix
- * @see #getTrackParameter
- */
- public double getErrorMatrixElement(int i, int j)
- {
- return (i < j) ? m_err[j][i] : m_err[i][j];
- }
-
- /**
- * Get momentum at DOCA.
- */
- public double[] getMomentum()
- {
- return new double[] { getPX(), getPY(), getPZ() };
- }
-
- public double getOmega()
- {
- return m_parm[2];
- }
-
- public double getPX()
- {
- return getPt() * Math.cos(m_parm[1]);
- }
-
- public double getPY()
- {
- return getPt() * Math.sin(m_parm[1]);
- }
-
- public double getPZ()
- {
- return getPt() * m_parm[4];
- }
-
- public double getPhi0()
- {
- return m_parm[1];
- }
-
- /**
- * Get total momentum at DOCA.
- */
- public double getPtot()
- {
- return Math.sqrt((getPX() * getPX()) + (getPY() * getPY()) + (getPZ() * getPZ()));
- }
-
- public double getTanL()
- {
- return m_parm[4];
- }
-
- /**
- * Get an individual track parameter. <br>
- *
- * The track parameters for LCD are defined as follows
- * <table>
- * <tr><th>Index</th><th>Meaning</th></tr>
- * <tr><td> 0 </td><td> d0 = XY impact parameter </td><tr>
- * <tr><td> 1 </td><td> phi0 </td><tr> </td><tr>
- * <tr><td> 2 </td><td> omega = 1/curv.radius (negative for negative tracks) </td><tr>
- * <tr><td> 3 </td><td> z0 = z of track (z impact parameter) </td><tr>
- * <tr><td> 4 </td><td> s = tan lambda </td><tr>
- * </table>
- * @param i The index of the track parameter
- * @return The track parameter with the specified index
- *
- * All parameters are given at the DOCA.
- */
- public double getTrackParameter(int i)
- {
- return m_parm[i];
- }
-
- /**
- * Get the track parameters as an array
- * @see #getTrackParameter
- */
- public double[] getTrackParameters()
- {
- return m_parm;
- }
-
- /**
- * Get the unit charge, ie +1, 0, -1.
- */
- public int getUnitCharge()
- {
- if (m_Bz != 0)
- {
- return (int) Math.signum(m_parm[2]);
- }
- else
- {
- return 0;
- }
- }
-
- public double getZ0()
- {
- return m_parm[3];
- }
-
- /**
- * Get cos(Theta) as calculated from the
- * momentum vector at the DOCA. <br>
- *
- * Note: This is the same as getCosTheta()
- */
-
- public double magnitude()
- {
- return Math.sqrt(VecOp.dot(getDocaVec(), getDocaVec()));
- }
- public double magnitudeSquared()
- {
- return VecOp.dot(getDocaVec(), getDocaVec());
- }
-
- public double[] v()
- {
- return getDoca();
- }
-
- // IHep3Vector methods
- public double x()
- {
- return getDocaX();
- }
-
- public double y()
- {
- return getDocaY();
- }
-
- public double z()
- {
- return getDocaZ();
- }
-
- /**
- * Store the chi2.
- */
- void setChi2(double chi2)
- {
- m_chi2 = chi2;
- }
-
- /**
- * Return the chi2 from smearing. <br>
- *
- * Note: The chi2 is to be calculated and
- * stored by the smearing routine
- * using setChi2().
- */
- double getChi2()
- {
- return m_chi2;
- }
-
- /**
- * get cos(theta) at DOCA.
- */
- double getCosTheta()
- {
- return getPZ() / getPtot();
- }
-
- /**
- * Get coordinates of DOCA.
- */
- double[] getDoca()
- {
- return new double[] { getDocaX(), getDocaY(), getDocaZ() };
- }
-
- double[] getDocaMomentum(double[] refPoint)
- {
- return getDocaMomentumVec(refPoint).v();
- }
-
- /**
- * Calculate and get Doca momentum on track with respect to any space point.
- */
- Hep3Vector getDocaMomentumVec(Hep3Vector refPoint)
- {
- if ((refPoint.x() != 0.) || (refPoint.y() != 0.))
- {
- checkCalcDoca(refPoint);
-
- return m_pdoca_ref;
- }
- else
- {
- return getMomentumVec();
- }
- }
-
- Hep3Vector getDocaMomentumVec(double[] refPoint)
- {
- return getDocaMomentumVec(new BasicHep3Vector(refPoint[0], refPoint[1], refPoint[2]));
- }
-
- double[] getDocaPosition(double[] refPoint)
- {
- return getDocaPositionVec(refPoint).v();
- }
-
- //====================================================
- //
- // methods
- //
- //====================================================
-
- /**
- * Calculate and get Doca position on track with respect to
- * any space point.
- */
- Hep3Vector getDocaPositionVec(Hep3Vector refPoint)
- {
- if ((refPoint.x() != 0.) || (refPoint.y() != 0.))
- {
- checkCalcDoca(refPoint);
-
- return m_xdoca_ref;
- }
- else
- {
- return getDocaVec();
- }
- }
-
- Hep3Vector getDocaPositionVec(double[] refPoint)
- {
- return getDocaPositionVec(new BasicHep3Vector(refPoint[0], refPoint[1], refPoint[2]));
- }
-
- /**
- * Calculate and get path length on track for a doca to any space point
- * in respect to the track defining doca (with respect to the origin).
- * The length l is given in the transverse plane. <br>
- * Use L = l*tan(lambda) to convert.
- */
- double getDocaTransversePathLength(Hep3Vector refPoint)
- {
- if ((refPoint.x() != 0.) || (refPoint.y() != 0))
- {
- checkCalcDoca(refPoint);
-
- return m_l_ref;
- }
- else
- {
- return 0.;
- }
- }
-
- double getDocaTransversePathLength(double[] refPoint)
- {
- return getDocaTransversePathLength(new BasicHep3Vector(refPoint[0], refPoint[1], refPoint[2]));
- }
-
- /**
- * Get coordinates of DOCA.
- */
- Hep3Vector getDocaVec()
- {
- return new BasicHep3Vector(getDocaX(), getDocaY(), getDocaZ());
- }
-
- double getDocaX()
- {
- return (-m_parm[0] * Math.sin(m_parm[1]));
- }
-
- double getDocaY()
- {
- return (m_parm[0] * Math.cos(m_parm[1]));
- }
-
- double getDocaZ()
- {
- return (m_parm[3]);
- }
-
- /**
- * Set the (transverse) path length l0 to original track vertex.
- */
- void setL0(double l0)
- {
- m_l0 = l0;
- }
-
- /**
- * Get the (transverse) path length l0 to original track vertex.
- */
- double getL0()
- {
- return m_l0;
- }
-
- double[] getMomentum(double l)
- {
- return getMomentumVec(l).v();
- }
-
- /**
- * Calculate and get momentum on track with respect to
- * any path length l on track (l in xy plane).
- */
- Hep3Vector getMomentumVec(double l)
- {
- double phi0 = m_parm[1];
- double omega = m_parm[2];
- double tanl = m_parm[4];
-
- int iq = getUnitCharge();
-
- double phi = phi0 + (omega * l);
- double pt = Constants.fieldConversion * iq * m_Bz / omega;
-
- double px = pt * Math.cos(phi);
- double py = pt * Math.sin(phi);
- double pz = pt * tanl;
-
- // System.out.println("l: "+l+" p: ("+px+", "+py+", "+pz+")");
- return new BasicHep3Vector(px, py, pz);
- }
-
- Hep3Vector getMomentumVec()
- {
- return new BasicHep3Vector(getPX(), getPY(), getPZ());
- }
-
- /**
- * Change the number degrees of freedom.
- */
- void setNDF(int ndf)
- {
- m_ndf = ndf;
- }
-
- /**
- * Get the number degrees of freedom.
- *
- * Default is 5 unless changed with setNDF().
- */
- int getNDF()
- {
- return m_ndf;
- }
-
- double[] getPosition(double l)
- {
- return getPositionVec(l).v();
- }
-
- /**
- * Calculate and get position on track with respect to
- * any path length l on track (l in xy plane).
- */
- Hep3Vector getPositionVec(double l)
- {
- double d0 = m_parm[0];
- double phi0 = m_parm[1];
- double omega = m_parm[2];
- double z0 = m_parm[3];
- double tanl = m_parm[4];
-
- double phi = phi0 + omega;
- double rho = 1 / omega;
-
- double x = (rho * Math.sin(phi)) - ((rho + d0) * Math.sin(phi0));
- double y = (-rho * Math.cos(phi)) + ((rho + d0) * Math.cos(phi0));
- double z = z0 + (l * tanl);
-
- return new BasicHep3Vector(x, y, z);
- }
-
- /**
- * Get transverse momentum at DOCA.
- */
- double getPt()
- {
- if (m_parm[2] != 0.)
- {
- return Math.abs(Constants.fieldConversion * m_Bz / m_parm[2]);
- }
- else
- {
- return 0.;
- }
- }
-
- /**
- * Get theta angle at DOCA.
- */
- double getTheta()
- {
- return Math.atan2(getPt(), getPZ());
- }
-
- /**
- * Calculate the error matrix for the momentum
- * for a point on the track specified by l.
- * Result is given as a 3x3 array for the matrix.
- */
- double[][] calcMomentumErrorMatrix(double l)
- {
- double rho = 1. / getOmega();
- double phi = getPhi0() + (getOmega() * l);
- double tanl = getTanL();
- double c = Constants.fieldConversion * Math.abs(m_Bz);
- double sphi = Math.sin(phi);
- double cphi = Math.cos(phi);
-
- Matrix tMatrix = new Matrix(5, 3, 0.);
- tMatrix.set(1, 0, -c * rho * sphi);
- tMatrix.set(1, 1, c * rho * cphi);
- tMatrix.set(2, 0, (-c * rho * rho * cphi) - (c * rho * l * sphi));
- tMatrix.set(2, 1, (-c * rho * rho * sphi) + (c * rho * l * cphi));
- tMatrix.set(2, 2, -c * rho * rho * tanl);
- tMatrix.set(4, 2, c * rho);
-
- Matrix errorMatrix = new Matrix(getErrorMatrix());
- Matrix pErrorMatrix = tMatrix.transpose().times(errorMatrix.times(tMatrix));
-
- return pErrorMatrix.getArrayCopy();
- }
-
- /**
- * Calculate the error matrix for the position coordinates
- * for a point on the track specified by l.
- * Result is given as a 3x3 array for the matrix.
- */
- double[][] calcPositionErrorMatrix(double l)
- {
- double d0 = getD0();
- double rho = 1. / getOmega();
- double phi = getPhi0() + (getOmega() * l);
- double sphi0 = Math.sin(getPhi0());
- double cphi0 = Math.cos(getPhi0());
- double sphi = Math.sin(phi);
- double cphi = Math.cos(phi);
-
- Matrix tMatrix = new Matrix(5, 3, 0.);
- tMatrix.set(0, 0, -sphi0);
- tMatrix.set(0, 1, cphi0);
- tMatrix.set(1, 0, (rho * cphi) - ((rho + d0) * cphi0));
- tMatrix.set(1, 1, (rho * sphi) - ((rho + d0) * sphi0));
- tMatrix.set(2, 0, (rho * l * cphi) - (rho * rho * (sphi - sphi0)));
- tMatrix.set(2, 1, (rho * l * sphi) + (rho * rho * (cphi - cphi0)));
- tMatrix.set(3, 2, 1.);
- tMatrix.set(4, 2, l);
-
- Matrix errorMatrix = new Matrix(getErrorMatrix());
-
- // MyContext.println(MyContext.getHeader());
- // MyContext.printMatrix("Error matrix:",errorMatrix,10,15);
- // MyContext.printMatrix("Transf matrix:",tMatrix,10,15);
- Matrix xErrorMatrix = tMatrix.transpose().times(errorMatrix.times(tMatrix));
-
- return xErrorMatrix.getArrayCopy();
- }
-
- void fillErrorMatrix(double[][] errorMatrix)
- {
- m_err = errorMatrix;
-
- // fill upper triangle from lower triangle
- for (int i = 0; i < 5; i++)
- for (int j = 0; j < i; j++)
- m_err[j][i] = m_err[i][j];
- }
-
- //====================================================
- //
- // private methods
- //
- //====================================================
-
- /*
- * Calculate the DOCA for a set of parameters with respect to the origin.
- */
- private void calculateDoca(double[] momentum, double[] trackPoint, double q)
- {
- Hep3Vector p = new BasicHep3Vector(momentum[0], momentum[1], momentum[2]);
- Hep3Vector x = new BasicHep3Vector(trackPoint[0], trackPoint[1], trackPoint[2]);
- calculateDoca(p, x, q);
- }
-
- /*
- * Calculate the DOCA for a set of parameters in vectors with respect to the origin.
- */
- private void calculateDoca(Hep3Vector momentum, Hep3Vector trackPoint, double charge)
- {
- reset();
-
- Hep3Vector xOrigin = new BasicHep3Vector(0., 0., 0.);
-
- Hep3Vector[] result = calculateDoca(momentum, trackPoint, charge, xOrigin);
- Hep3Vector xdoca = result[0];
- Hep3Vector pdoca = result[1];
- Hep3Vector dphdl = result[2];
-
- int iq = (int) (charge / Math.abs(charge));
- double pt = Math.sqrt((pdoca.x() * pdoca.x()) + (pdoca.y() * pdoca.y()));
-
- // now calculate track parameters
- double d0 = Math.sqrt((xdoca.x() * xdoca.x()) + (xdoca.y() * xdoca.y()));
- if (VecOp.cross(xdoca, pdoca).z() > 0)
- {
- d0 = -d0;
- }
-
- double phi0 = Math.atan2(pdoca.y(), pdoca.x());
- if (phi0 > Math.PI)
- {
- phi0 -= (2 * Math.PI);
- }
- if (phi0 < -Math.PI)
- {
- phi0 += (2 * Math.PI);
- }
-
- double omega = Constants.fieldConversion * iq * m_Bz / pt;
- double tanl = pdoca.z() / pt;
- double z0 = xdoca.z();
-
- // now fill trackparameters
- m_parm[0] = d0;
- m_parm[1] = phi0;
- m_parm[2] = omega;
- m_parm[3] = z0;
- m_parm[4] = tanl;
-
- // save the distance to orignial track vertex
- m_l0 = -dphdl.y();
-
- // System.out.println("DocaTrackParameters: xdoca = ("+
- // xdoca.x()+", "+xdoca.y()+", "+xdoca.z()+")");
- // System.out.println("DocaTrackParameters: pdoca = ("+
- // pdoca.x()+", "+pdoca.y()+", "+pdoca.z()+")");
- // System.out.println("DocaTrackParameters: d0: "+m_parm[0]+
- // " phi0: "+m_parm[1]+" omega: "+m_parm[2]+
- // " z0: "+m_parm[3]+" tanl: "+m_parm[4]);
- // System.out.println("DocaTrackParameters: m_l0 = "+m_l0);
- }
-
- /*
- * Calculate DOCA position and momentum vectors with respect to any
- * space point.
- */
- private Hep3Vector[] calculateDoca(Hep3Vector momentum, Hep3Vector trackPoint, double charge, Hep3Vector refPoint)
- {
- // subtract refPoint
- Hep3Vector xp = VecOp.sub(trackPoint, refPoint);
-
- int iq = (int) (charge / Math.abs(charge));
- double pt = Math.sqrt((momentum.x() * momentum.x()) + (momentum.y() * momentum.y()));
- double tanl = momentum.z() / pt;
- double rho = pt / (m_Bz * Constants.fieldConversion);
-
- BasicHep3Vector xdoca;
- Hep3Vector pdoca;
- Hep3Vector dphdl;
-
- // System.out.println("calculateDoca: m_flip: "+m_flip+" iq: "+iq);
- if (xp.magnitude() > 0.) // no need for calculation if xp = (0,0,0) !
- {
- // calculate position and momentum at doca
- Hep3Vector nzv = new BasicHep3Vector(0., 0., iq);
- Hep3Vector xxc = new BasicHep3Vector(VecOp.cross(momentum, nzv).x(), VecOp.cross(momentum, nzv).y(), 0.);
- Hep3Vector nxxc = VecOp.unit(xxc);
- BasicHep3Vector xc = (BasicHep3Vector) VecOp.add(xp, VecOp.mult(rho, nxxc));
- xc.setV(xc.x(), xc.y(), 0.);
-
- Hep3Vector nxc = VecOp.unit(xc);
-
- BasicHep3Vector catMC = (BasicHep3Vector) VecOp.cross(nzv, nxc);
- catMC.setV(catMC.x(), catMC.y(), 0.);
-
- Hep3Vector ncatMC = VecOp.unit(catMC);
-
- pdoca = new BasicHep3Vector(pt * ncatMC.x(), pt * ncatMC.y(), momentum.z());
-
- double dphi = Math.asin(VecOp.cross(nxxc, nxc).z());
- double dl = -dphi * rho * iq;
-
- xdoca = (BasicHep3Vector) VecOp.add(xc, VecOp.mult(-rho, nxc));
- xdoca.setV(xdoca.x(), xdoca.y(), xp.z() + (dl * tanl));
-
- // save dphi and dl
- dphdl = new BasicHep3Vector(dphi, dl, 0.);
- }
- else
- {
- xdoca = (BasicHep3Vector) xp;
- pdoca = momentum;
- dphdl = new BasicHep3Vector();
- }
-
- // add refPoint back in again
- xdoca = (BasicHep3Vector) VecOp.add(xdoca, refPoint);
-
- return new Hep3Vector[] { xdoca, pdoca, dphdl };
- }
-
- /*
- * Calculate Doca for this track with respect to any space point,
- * if this has not been done before.
- */
- private void checkCalcDoca(Hep3Vector refPoint)
- {
- // hassle with calculation only, if not done before yet!
- if ((m_xref == null) || (refPoint.x() != m_xref.x()) || (refPoint.y() != m_xref.y()) || (refPoint.z() != m_xref.z()))
- {
- m_xref = refPoint;
-
- // find doca vectors
- Hep3Vector xdoca = getDocaVec();
- Hep3Vector pdoca = getMomentumVec();
- int iq = getUnitCharge();
-
- // calculate doca to refPoint
- Hep3Vector[] result = calculateDoca(pdoca, xdoca, (double) iq, refPoint);
- m_xdoca_ref = result[0];
- m_pdoca_ref = result[1];
-
- // distance along track to original point
- m_l_ref = result[2].y();
- }
- }
-
- private void reset()
- {
- for (int i = 0; i < 5; i++)
- {
- m_parm[i] = 0;
- for (int j = 0; j < 5; j++)
- {
- m_err[i][j] = 0;
- }
- }
- m_l0 = 0.;
- }
-}
lcsim/src/org/lcsim/contrib/timb/mc/fast/tracking
diff -N LookupTable.java
--- LookupTable.java 26 May 2006 07:21:58 -0000 1.1
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,126 +0,0 @@
-package org.lcsim.mc.fast.tracking;
-
-import java.io.BufferedReader;
-import java.io.IOException;
-
-import java.util.Arrays;
-import java.util.StringTokenizer;
-
-
-class LookupTable
-{
- // cosine theta
- private double[] m_key1;
- // momentum
- private double[] m_key2;
- private double[][] m_matrix;
- private int m_numBins1;
- private int m_numBins2;
-
- // vector units
- // 1 dr [ cm ] 10.
- // 2 dphi [ ] 1.
- // 3 domega [ cm-1 ] .1
- // 4 dz [ cm ] 10.
- // 5 dlambda [ ] 1.
- //
- double[] conversionFromCmToMm = {10.0, 1.0, 0.1, 10.0, 1.0};
-
- double[][] conversionFromCmToMmMatrix = { {100.00, 10.00, 1.00, 100.00, 10.00},
- { 10.00, 1.00, 0.10, 10.00, 1.00},
- { 1.00, 0.10, 0.01, 1.00, 0.10},
- {100.00, 10.00, 1.00, 100.00, 10.00},
- { 10.00, 1.00, 0.10, 10.00, 1.00} };
-
- LookupTable(BufferedReader in, int iTerm, int jTerm) throws IOException
- {
- // read in the number of cosine theta points
- int m_numBins1 = Integer.parseInt(in.readLine());
-
- // read in the number of momentum points
- int m_numBins2 = Integer.parseInt(in.readLine());
-
- m_matrix = new double[m_numBins1][m_numBins2];
- m_key1 = new double[m_numBins1];
- m_key2 = new double[m_numBins2];
-
- for (int i = 0; i < m_numBins1; i++) // i is # of cosine theta bin
- {
- m_key1[i] = Double.valueOf(in.readLine()).doubleValue(); // cosine theta
- for (int j = 0; j < m_numBins2; j++) // j is # of momentum bin
- {
- StringTokenizer t = new StringTokenizer(in.readLine());
- m_key2[j] = Double.valueOf(t.nextToken()).doubleValue();
- m_matrix[i][j] = Double.valueOf(t.nextToken()).doubleValue()*conversionFromCmToMmMatrix[iTerm][jTerm]; // momentum
- }
- }
- if (!in.readLine().equals("end"))
- {
- throw new IOException("Missing end in lookup table");
- }
- }
-
- double interpolateVal(double val1, double val2)
- {
- int index1 = binarySearch(m_key1, val1);
- int index2 = binarySearch(m_key2, val2);
- // System.out.println(" index1= "+index1+" index2= "+index2);
-
- double t = (val1 - m_key1[index1]) / (m_key1[index1 + 1] - m_key1[index1]);
-
- double u = (val2 - m_key2[index2]) / (m_key2[index2 + 1] - m_key2[index2]);
-
- double y1 = m_matrix[index1][index2];
- double y2 = m_matrix[index1 + 1][index2];
- double y3 = m_matrix[index1 + 1][index2 + 1];
- double y4 = m_matrix[index1][index2 + 1];
-
- return ((1 - t) * (1 - u) * y1) + (t * (1 - u) * y2) + (t * u * y3) + ((1 - t) * u * y4);
- }
- private int binarySearch(double[] key, double value)
-// {
-// int result = binarySearchX(key,value);
-// System.out.print("Looking for "+value+" in [");
-// for (int i=0; i<key.length; i++) System.out.print(key[i]+",");
-// System.out.print("] ");
-// System.out.println("result="+result);
-// return result;
-// }
-//
-// private int binarySearchX(double[] key, double value)
- {
- // System.out.print("Looking for "+value+" in [");
- // for (int i=0; i<key.length; i++) System.out.print(key[i]+",");
- // System.out.print("] ");
- // System.out.println(" key.length= "+key.length);
- if (value < key[0])
- {
- // System.out.println("Interpolation out of range: lower: "+value+" < "+key[0]);
- //throw new RuntimeException("Interpolation out of range: lower: "+value+" < "+key[0]);
- return 0;
- }
- // else if(value >= key[key.length-1]){
- // System.out.println("Interpolation out of range: upper: "+value+" >= "+key[key.length-1]);
- // }
-
- int pos = Arrays.binarySearch(key, value);
- if (pos > 0)
- {
- // System.out.println("pos= "+pos);
- return Math.min(pos,key.length-2);
- }
- else
- {
- // System.out.println("-pos-2= "+(-pos-2));
- return Math.min(-pos-2, key.length-2);
- }
-
- // // Ok, this isn't really a binary search, probably doesn't matter
- // for (int i=1; i<key.length; i++) if (value<key[i]) return i-1;
- //
- //
- // throw new LCDException("Interpolation out of range: upper: "
- // +value+" >= "+key[key.length-1]);
- // System.out.println("Interpolation out of range: upper: "+value+" >= "+key[key.length-1]);
- }
-}
lcsim/src/org/lcsim/contrib/timb/mc/fast/tracking
diff -N MCFastTracking.java
--- MCFastTracking.java 26 May 2006 07:21:59 -0000 1.1
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,178 +0,0 @@
-package org.lcsim.mc.fast.tracking;
-
-import hep.physics.particle.Particle;
-import java.io.IOException;
-import java.util.ArrayList;
-import java.util.Iterator;
-import java.util.List;
-
-import org.lcsim.conditions.ConditionsEvent;
-import org.lcsim.conditions.ConditionsListener;
-import org.lcsim.conditions.ConditionsSet;
-import org.lcsim.event.EventHeader;
-import org.lcsim.event.Track;
-import org.lcsim.util.Driver;
-import org.lcsim.particle.ParticleUtilities;
-
-
-/**
- * Fast Monte Carlo tracking simulator
- */
-public class MCFastTracking extends Driver implements ConditionsListener
-{
- private List<Particle> bquark_list = new ArrayList<Particle>();
- private List<Particle> anti_bquark_list = new ArrayList<Particle>();
- private Particle cquark;
- private Particle anti_cquark;
- private TrackResolutionTables parm;
- private SimpleTables SmTbl;
- private boolean beamSpotConstraint;
- private boolean simple;
- private final static double[] IP = { 0, 0, 0 };
- private static int nprint = 0 ;
- private static final int nprint_max = -2 ;
-
- public MCFastTracking()
- {
- this(false);
- }
-
- public MCFastTracking(boolean beamSpotConstraint)
- {
- this.beamSpotConstraint = beamSpotConstraint;
- }
-
- public MCFastTracking(boolean beamSpotConstraint, boolean simple)
- {
- this.beamSpotConstraint = beamSpotConstraint;
- this.simple = simple;
- }
-
- public void setBeamSpotConstraint(boolean beamSpotConstraint)
- {
- this.beamSpotConstraint = beamSpotConstraint;
- if (parm != null)
- {
- ConditionsSet conditions = getConditionsManager().getConditions("TrackParameters");
- parm = setTrackResolutionTables(conditions,beamSpotConstraint);
- }
- }
-
- public boolean isBeamSpotConstraint()
- {
- return this.beamSpotConstraint;
- }
- private TrackResolutionTables setTrackResolutionTables(ConditionsSet conditions,boolean beamSpotConstraint)
- {
- try
- {
- return new TrackResolutionTables(conditions, beamSpotConstraint);
- }
- catch (IOException x)
- {
- throw new RuntimeException("Error reading track resolution tables",x);
- }
- }
-
- protected void process(EventHeader event)
- {
- if (parm == null)
- {
- ConditionsSet conditions = getConditionsManager().getConditions("TrackParameters");
- conditions.addConditionsListener(this);
- parm = setTrackResolutionTables(conditions, beamSpotConstraint);
- }
-
- if (SmTbl == null)
- {
- ConditionsSet simpleconditions = getConditionsManager().getConditions("SimpleTrack");
- simpleconditions.addConditionsListener(this);
- SmTbl = new SimpleTables(simpleconditions);
- }
-
- double bField = event.getDetector().getFieldMap().getField(IP)[2];
- boolean hist = getHistogramLevel() > 0;
-
- List<Track> trackList = new ArrayList<Track>();
- bquark_list.clear();
- anti_bquark_list.clear();
- cquark=null;
- anti_cquark=null;
- for (Iterator i = event.getMCParticles().iterator(); i.hasNext();)
- {
- Particle p = (Particle) i.next();
- if(Math.abs(p.getPDGID()) == 5) {
- boolean addp;
-
- if(p.getParents().size() == 0) addp = false;
- else addp=Math.abs(((Particle)p.getParents().get(0)).getPDGID()) != 5;
-
- if(p.getPDGID() == 5 && addp) bquark_list.add(p);
- if(p.getPDGID() == -5 && addp) anti_bquark_list.add(p);
- }
- if(p.getPDGID() == 4) cquark=p;
- if(p.getPDGID() == -4) anti_cquark=p;
-
- // filter for FINAL_STATE
- if (p.getGeneratorStatus() != Particle.FINAL_STATE)
- {
- continue;
- }
- double pCharge = p.getCharge();
- if (pCharge == 0 || pCharge == Double.NaN || pCharge == Double.NEGATIVE_INFINITY || pCharge == Double.POSITIVE_INFINITY)
- {
- continue;
- }
-
- double[] momentum = p.getMomentum().v();
- double pt2 = (momentum[0] * momentum[0]) + (momentum[1] * momentum[1]);
- double pt = Math.sqrt(pt2);
- double ptot = Math.sqrt(pt2 + (momentum[2] * momentum[2]));
- double cosTheta = momentum[2] / ptot;
-
- // within acceptance
- if (pt < parm.getPtMin())
- {
- continue;
- }
- if (Math.abs(cosTheta) > parm.getPolarOuter())
- {
- continue;
- }
-
- trackList.add(new ReconTrack(bField, parm, SmTbl, getRandom(), p, hist, simple));
- }
- nprint++;
- if(nprint <= nprint_max) {
-
- for(Particle bquark : bquark_list) {
- System.out.println(" MCFastTracking nprint= "+nprint+" bquark.getPDGID= "+bquark.getPDGID()+" bquark.getGeneratorStatus= "+bquark.getGeneratorStatus()
- +" ((Particle)bquark.getParents().get(0)).getPDGID= "+((Particle)bquark.getParents().get(0)).getPDGID());
- ParticleUtilities.dumpParticleHierarchy(bquark);
- }
- for(Particle anti_bquark : anti_bquark_list) {
- System.out.println(" MCFastTracking nprint= "+nprint+" anti_bquark.getPDGID= "+anti_bquark.getPDGID()+" anti_bquark.getGeneratorStatus= "+anti_bquark.getGeneratorStatus()
- +" ((Particle)anti_bquark.getParents().get(0)).getPDGID= "+((Particle)anti_bquark.getParents().get(0)).getPDGID());
- ParticleUtilities.dumpParticleHierarchy(anti_bquark);
- }
- if(cquark != null && bquark_list.size() == 0) {
- System.out.println(" MCFastTracking nprint= "+nprint+" cquark.getPDGID= "+cquark.getPDGID()+" cquark.getGeneratorStatus= "+cquark.getGeneratorStatus());
- ParticleUtilities.dumpParticleHierarchy(cquark);
- }
- if(anti_cquark != null && anti_bquark_list.size() == 0) {
- System.out.println(" MCFastTracking nprint= "+nprint+" anti_cquark.getPDGID= "+anti_cquark.getPDGID()+" anti_cquark.getGeneratorStatus= "+anti_cquark.getGeneratorStatus());
- ParticleUtilities.dumpParticleHierarchy(anti_cquark);
- }
- }
-
- event.put(EventHeader.TRACKS, trackList, Track.class, 0);
- }
-
- public void conditionsChanged(ConditionsEvent event)
- {
- ConditionsSet conditions = getConditionsManager().getConditions("TrackParameters");
- ConditionsSet simpleconditions = getConditionsManager().getConditions("SimpleTrack");
- parm = setTrackResolutionTables(conditions, beamSpotConstraint);
- SmTbl = new SimpleTables(simpleconditions);
- }
-}
lcsim/src/org/lcsim/contrib/timb/mc/fast/tracking
diff -N SimpleTables.java
--- SimpleTables.java 26 May 2006 07:22:01 -0000 1.1
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,61 +0,0 @@
-
-package org.lcsim.mc.fast.tracking;
-
-
-import org.lcsim.conditions.ConditionsSet;
-
-/**
- *
- * @author Daniel
- */
-public class SimpleTables {
-
- private double ConstantTerm;
- private double ThetaTerm;
- private double TanLambdaErrorScale;
- private double PhiErrorScale;
- private double D0ErrorScale;
- private double Z0ErrorScale;
-
- /** Creates a new instance of SimpleTables */
- SimpleTables(ConditionsSet set)
- {
- ConstantTerm = set.getDouble("ConstantTerm");
- ThetaTerm = set.getDouble("ThetaTerm");
- TanLambdaErrorScale = set.getDouble("TanLambdaErrorScale");
- PhiErrorScale = set.getDouble("PhiErrorScale");
- D0ErrorScale = set.getDouble("D0ErrorScale");
- Z0ErrorScale = set.getDouble("Z0ErrorScale");
- }
-
- public double getConstantTerm()
- {
- return ConstantTerm;
- }
-
- public double getThetaTerm()
- {
- return ThetaTerm;
- }
- public double getTanLambdaErrorScale()
- {
- return TanLambdaErrorScale;
- }
-
- public double getPhiErrorScale()
- {
- return PhiErrorScale;
- }
-
- public double getD0ErrorScale()
- {
- return D0ErrorScale;
- }
-
- public double getZ0ErrorScale()
- {
- return Z0ErrorScale;
- }
-
-
-}
lcsim/src/org/lcsim/contrib/timb/mc/fast/tracking
diff -N SmearTrackSimple.java
--- SmearTrackSimple.java 26 May 2006 07:22:03 -0000 1.1
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,112 +0,0 @@
-package org.lcsim.mc.fast.tracking;
-
-import Jama.EigenvalueDecomposition;
-import Jama.Matrix;
-
-import org.lcsim.mc.fast.tracking.SimpleTables;
-import java.util.Random;
-
-
-/**
- * @author T. Barklow
- */
-class SmearTrackSimple
-{
- /**
- * Smear track parameters according to modified version of track's stored error matrix.
- *
- * @see TrackParameters
- */
- static DocaTrackParameters smearTrackSimple(double bField, TrackParameters noSmear, Random rand, SimpleTables SmTbl, double pt, boolean hist)
- {
-
- final double errScale=0.0001 ;
- final double eMScale=1.e14 ;
- // get copy of error matrix and prepare for modification
- Matrix eM = Matrix.constructWithCopy(noSmear.getErrorMatrix());
- double[] errscale = { SmTbl.getD0ErrorScale(), SmTbl.getPhiErrorScale(), 1., SmTbl.getZ0ErrorScale(), SmTbl.getTanLambdaErrorScale() };
- double[] oldDiagErr = new double[5];
- double[] newDiagErr = new double[5];
- for(int i=0; i<5 ; i++) {
- oldDiagErr[i]= Math.sqrt(noSmear.getErrorMatrix()[i][i]);
- if(i == 2) {
- double th = Math.atan(1/(noSmear.getTanL()));
- double a = SmTbl.getConstantTerm();
- double b = SmTbl.getThetaTerm()/(pt*Math.sin(th));
- newDiagErr[i]= pt * Math.sqrt(a*a + b*b);
- }
- else {
- newDiagErr[i]= errscale[i]*oldDiagErr[i];
- }
- eM.set(i,i,Math.pow(newDiagErr[i],2.));
- for(int j=0; j<i; j++) {
- eM.set(i,j,noSmear.getErrorMatrix()[i][j]*newDiagErr[i]*newDiagErr[j]/oldDiagErr[i]/oldDiagErr[j]);
- eM.set(j,i,eM.get(i,j));
- }
- }
- Matrix M = eM.copy();
- if (M.det() <= 0.)
- {
- throw new RuntimeException("Error matrix not positive definite!");
- }
-
- // run Eigenvalue decomposition and get matrices and vectors
- EigenvalueDecomposition eig = M.eig();
-
- Matrix T = eig.getV();
- double[] er = eig.getRealEigenvalues();
- double[] ei = eig.getImagEigenvalues();
-
- // sanity check: det(T) != 0
- if (T.det() == 0.)
- {
- throw new RuntimeException("Non orthogonal basis!");
- }
-
- // sanity check: no imaginary eigenvalues
- for (int i = 0; i < ei.length; i++)
- if (ei[i] != 0.)
- {
- throw new RuntimeException("Imaginary Eigenvalues seen!");
- }
-
- // now do the real smearing
- double[] dev = new double[5];
- for (int i = 0; i < er.length; i++)
- {
- if (er[i] <= 0)
- {
- throw new RuntimeException("Non-positive Eigenvalue seen!");
- }
- dev[i] = Math.sqrt(er[i]) * rand.nextGaussian();
- }
-
- Matrix shift = T.times(new Matrix(dev, 5));
- Matrix val = new Matrix(noSmear.getTrackParameters(), 5);
- double[] newval = (val.plus(shift)).getColumnPackedCopy();
-
- // adjust new phi value to [-pi,pi] if necessary
- if (newval[1] > Math.PI)
- {
- newval[1] -= (2. * Math.PI);
- }
- if (newval[1] < -Math.PI)
- {
- newval[1] += (2. * Math.PI);
- }
-
- // Chi2 calculation
- double chi2 = ((shift.transpose()).times((M.inverse()).times(shift))).get(0, 0);
-
- // DocaTrackParameters smeared = new DocaTrackParameters(bField, newval, (eM.timesEquals(eMScale)).getArray(), chi2);
- // DocaTrackParameters smeared = new DocaTrackParameters(bField, newval, eM.getArray(), chi2);
- // DocaTrackParameters smeared = new DocaTrackParameters(bField, newval, (Matrix.constructWithCopy(noSmear.getErrorMatrix()).timesEquals(errScale)).getArray(), chi2);
- DocaTrackParameters smeared = new DocaTrackParameters(bField, newval, noSmear.getErrorMatrix(), chi2);
- if (noSmear instanceof DocaTrackParameters)
- {
- smeared.setL0(((DocaTrackParameters) noSmear).getL0());
- }
-
- return smeared;
- }
-}
CVSspam 0.2.8
