10 modified files
lcsim/src/org/lcsim/mc/fast/cluster/ronan
diff -u -r1.5 -r1.6
--- ReconCluster.java 25 Aug 2005 23:36:53 -0000 1.5
+++ ReconCluster.java 8 Sep 2006 03:19:43 -0000 1.6
@@ -15,10 +15,13 @@
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)
{
@@ -32,6 +35,27 @@
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
@@ -47,13 +71,18 @@
protected void smearEnergy(Random rand, double E, boolean hist)
{
- double sigma = ((a / Math.sqrt(E)) + b) * E;
+ sigma = ((a / Math.sqrt(E)) + b) * E;
- energy = 0;
- while (energy <= mcp.getMass())
- {
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)
lcsim/src/org/lcsim/mc/fast/cluster/ronan
diff -u -r1.9 -r1.10
--- MCFastRonan.java 8 Nov 2005 17:04:32 -0000 1.9
+++ MCFastRonan.java 8 Sep 2006 03:19:43 -0000 1.10
@@ -12,14 +12,19 @@
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 $Id: MCFastRonan.java,v 1.9 2005/11/08 17:04:32 tonyj Exp $
+ * @version $Id: MCFastRonan.java,v 1.10 2006/09/08 03:19:43 ngraf Exp $
*/
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;
@@ -27,6 +32,7 @@
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;
@@ -39,7 +45,7 @@
conditions.addConditionsListener(this);
clusterParm = new ClusterResolutionTables(conditions);
}
- List cl = new ArrayList();
+ List<Cluster> cl = new ArrayList<Cluster>();
boolean hist = getHistogramLevel() > 0;
@@ -75,11 +81,15 @@
Random rand = getRandom();
// Photons
- if (PDGID == PhotonID)
+ if (PDGID == PhotonID || PDGID == ElecID)
{
// within acceptance
- double thing = (1 - 1 / ( 1 + Math.exp( (E-clusterParm.getEMOnset())*clusterParm.getEMSharpness() ) ));
- if (rand.nextDouble() > thing)
+ // double thing = (1 - 1 / ( 1 + Math.exp( (E-clusterParm.getEMOnset())*clusterParm.getEMSharpness() ) ));
+ //if (rand.nextDouble() > thing)
+ //{
+ // continue;
+ //}
+ if (E < clusterParm.getEMOnset())
{
continue;
}
@@ -93,23 +103,39 @@
}
// Neutral hadrons
- else if (charge == 0)
+ else if (PDGID != MuID)
{
// within acceptance
- double thing = (1 - 1 / ( 1 + Math.exp( (E-clusterParm.getHADOnset())*clusterParm.getHADSharpness() ) ));
- if (rand.nextDouble() > thing)
- {
+ //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;
+ 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);
}
lcsim/src/org/lcsim/mc/fast/cluster/ronan
diff -u -r1.3 -r1.4
--- ClusterResolutionTables.java 9 Aug 2005 18:34:45 -0000 1.3
+++ ClusterResolutionTables.java 8 Sep 2006 03:19:43 -0000 1.4
@@ -1,9 +1,18 @@
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;
@@ -25,6 +34,12 @@
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");
@@ -44,6 +59,16 @@
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()
lcsim/src/org/lcsim/mc/fast/tracking
diff -u -r1.1 -r1.2
--- SmearTrackSimple.java 9 Aug 2005 18:34:03 -0000 1.1
+++ SmearTrackSimple.java 8 Sep 2006 03:19:44 -0000 1.2
@@ -1,87 +1,112 @@
package org.lcsim.mc.fast.tracking;
-import java.util.Random;
+import Jama.EigenvalueDecomposition;
+import Jama.Matrix;
+
import org.lcsim.mc.fast.tracking.SimpleTables;
-import org.lcsim.util.aida.AIDA;
+import java.util.Random;
+
/**
- *
- * @author Daniel
+ * @author T. Barklow
*/
-public class SmearTrackSimple
- {
-
- /** Creates a new instance of SmearTrackSimple */
+class SmearTrackSimple
+{
/**
- * Smear track parameters according to the track's stored error matrix.
- *
- * @see TrackParameters
- */
- private static AIDA aida = AIDA.defaultInstance();
-
+ * 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)
{
- // get error matrix and do a sanity check
-
- double[] dev = noSmear.getTrackParameters();
- for (int i = 0; i < 5; i++)
- {
- if (noSmear.getErrorMatrix()[i][i] <= 0)
- {
- throw new RuntimeException("Non-positive Eigenvalue seen!");
- }
-
-
- if (i == 2)
- {
- double th = Math.atan(1/(noSmear.getTanL()));
-
- double a = SmTbl.getConstantTerm();
- double b = SmTbl.getThetaTerm()/(pt*Math.sin(th));
-
- double smmat = (Math.sqrt(noSmear.getErrorMatrix()[i][i])) / (Math.abs(pt * dev[i]));
- double p = Math.sqrt( noSmear.getPX()*noSmear.getPX() + noSmear.getPY()*noSmear.getPY() + noSmear.getPZ()*noSmear.getPZ() );
-
- if ((p < 102) && (p > 98) &&(hist))
- {
- aida.histogram2D("theta dependence of matrix element", 400, 0, 2, 400, .000005, .0005).fill(-Math.log10(1-Math.cos(th)),smmat);
- aida.histogram2D("theta dependence of parameterization", 400, 0, 2, 400, .000005, .0005).fill(-Math.log10(1-Math.cos(th)),Math.sqrt(a*a + b*b));
- }
-
- dev[i] = dev[i] + dev[i] * pt * Math.sqrt(a*a + b*b) * rand.nextGaussian();
-
- //dev[i] = dev[i] + Math.sqrt(noSmear.getErrorMatrix()[i][i]) * rand.nextGaussian();
-
- } else {
-
- dev[i] = dev[i] + Math.sqrt(noSmear.getErrorMatrix()[i][i]) * rand.nextGaussian();
-
- }
-
-
- }
-
-
- // adjust new phi value to [-pi,pi] if necessary
- if (dev[1] > Math.PI)
- {
- dev[1] -= (2. * Math.PI);
- }
- if (dev[1] < -Math.PI)
- {
- dev[1] += (2. * Math.PI);
- }
-
- // Chi2 calculation
- double chi2 = dev[0]*dev[0]*noSmear.getErrorMatrix()[0][0];
-
- DocaTrackParameters smeared = new DocaTrackParameters(bField, dev, noSmear.getErrorMatrix(), chi2);
- if (noSmear instanceof DocaTrackParameters)
- {
- smeared.setL0(((DocaTrackParameters) noSmear).getL0());
- }
-
- return smeared;
- }
+ 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;
+ }
}
lcsim/src/org/lcsim/mc/fast/tracking
diff -u -r1.6 -r1.7
--- DocaTrackParameters.java 8 Jul 2006 10:21:55 -0000 1.6
+++ DocaTrackParameters.java 8 Sep 2006 03:19:44 -0000 1.7
@@ -15,7 +15,7 @@
* with a MC truth particle. <br>
*
* @author Tony Johnson, Wolfgang Walkowiak
- * @version $Id: DocaTrackParameters.java,v 1.6 2006/07/08 10:21:55 jstrube Exp $
+ * @version $Id: DocaTrackParameters.java,v 1.7 2006/09/08 03:19:44 ngraf Exp $
*/
class DocaTrackParameters implements TrackParameters
{
@@ -296,7 +296,7 @@
*/
Hep3Vector getDocaMomentumVec(Hep3Vector refPoint)
{
- if ((refPoint.x() != 0.) && (refPoint.y() != 0.))
+ if ((refPoint.x() != 0.) || (refPoint.y() != 0.))
{
checkCalcDoca(refPoint);
@@ -330,7 +330,7 @@
*/
Hep3Vector getDocaPositionVec(Hep3Vector refPoint)
{
- if ((refPoint.x() != 0.) && (refPoint.y() != 0.))
+ if ((refPoint.x() != 0.) || (refPoint.y() != 0.))
{
checkCalcDoca(refPoint);
@@ -355,7 +355,7 @@
*/
double getDocaTransversePathLength(Hep3Vector refPoint)
{
- if ((refPoint.x() != 0.) && (refPoint.y() != 0))
+ if ((refPoint.x() != 0.) || (refPoint.y() != 0))
{
checkCalcDoca(refPoint);
lcsim/src/org/lcsim/mc/fast/tracking
diff -u -r1.3 -r1.4
--- LookupTable.java 15 Jul 2006 09:35:59 -0000 1.3
+++ LookupTable.java 8 Sep 2006 03:19:44 -0000 1.4
@@ -97,11 +97,11 @@
int pos = Arrays.binarySearch(key, value);
if (pos > 0)
{
- return pos;
+ return Math.min(pos,key.length-2);
}
else
{
- return -pos - 2;
+ return Math.min(-pos-2, key.length-2);
}
// // Ok, this isn't really a binary search, probably doesn't matter
lcsim/src/org/lcsim/mc/fast/tracking
diff -u -r1.2 -r1.3
--- SimpleTables.java 15 Jul 2006 09:35:59 -0000 1.2
+++ SimpleTables.java 8 Sep 2006 03:19:44 -0000 1.3
@@ -8,16 +8,25 @@
*
* @author Daniel
*/
-public class SimpleTables {
-
- private double ConstantTerm;
- private double ThetaTerm;
+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 */
- public SimpleTables(ConditionsSet set)
+ public SimpleTables(ConditionsSet set)
{
ConstantTerm = set.getDouble("ConstantTerm");
- ThetaTerm = set.getDouble("ThetaTerm");
+ ThetaTerm = set.getDouble("ThetaTerm");
+ TanLambdaErrorScale = set.getDouble("TanLambdaErrorScale");
+ PhiErrorScale = set.getDouble("PhiErrorScale");
+ D0ErrorScale = set.getDouble("D0ErrorScale");
+ Z0ErrorScale = set.getDouble("Z0ErrorScale");
}
public double getConstantTerm()
@@ -29,5 +38,24 @@
{
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/mc/fast/reconstructedparticle
diff -u -r1.10 -r1.11
--- MCFastReconstructedParticleDriver.java 24 Aug 2005 00:06:29 -0000 1.10
+++ MCFastReconstructedParticleDriver.java 8 Sep 2006 03:19:45 -0000 1.11
@@ -5,6 +5,8 @@
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;
@@ -40,6 +42,10 @@
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()
@@ -61,6 +67,10 @@
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)
@@ -70,16 +80,29 @@
if (IDEff == null)
{
- ConditionsSet conditions = getConditionsManager().getConditions("IDEfficiency");
- conditions.addConditionsListener(this);
- IDEff = new IDResolutionTables(conditions);
+ 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...
- List<Track> tracks = event.getTracks();
for(Track t : tracks)
{
ParticleType type = null;
@@ -90,7 +113,7 @@
int pdgid = p.getPDGID();
- // electrons and muons are special
+ // charged track id
if((abs(pdgid)== 11) && (rand.nextDouble() < IDEff.getElectronEff()))
{
type = rt.getCharge() > 0 ? eplus : eminus;
@@ -99,9 +122,17 @@
{
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 ? pizero : rt.getCharge() > 0 ? piplus : piminus;
+ type = rt.getCharge() > 0 ? piplus : piminus;
}
if ((p.getEnergy() > 10) && (hist))
@@ -117,7 +148,7 @@
}
// assume pion for remaining charged tracks
- MCFastReconstructedParticle rp = new MCFastReconstructedParticle(t, type, p);
+ MCFastReconstructedParticle rp = new MCFastReconstructedParticle(t, type, p, m_tc.get(t), IDEff.getWtChgTrkCal());
rpList.add(rp);
if (hist)
@@ -128,9 +159,9 @@
}
// loop over clusters...
- List<Cluster> clusters = event.getClusters();
for(Cluster c : clusters)
{
+ //if(m_ct.get(c) != null) continue;
Particle p = null;
ParticleType type = null;
// photons for EM
@@ -138,6 +169,7 @@
{
ReconEMCluster emc = (ReconEMCluster) c;
p = emc.getMCParticle();
+ if(abs(p.getCharge()) > 0.) continue;
type = photon;
if (hist)
{
@@ -149,6 +181,7 @@
{
ReconHADCluster emc = (ReconHADCluster) c;
p = emc.getMCParticle();
+ if(abs(p.getCharge()) > 0.) continue;
int pdgid = p.getPDGID();
if (hist)
{
@@ -179,8 +212,8 @@
public void conditionsChanged(ConditionsEvent event)
{
- ConditionsSet conditions = getConditionsManager().getConditions("IDEfficiency");
- IDEff = new IDResolutionTables(conditions);
+ ConditionsSet idconditions = getConditionsManager().getConditions("IDEfficiency");
+ IDEff = new IDResolutionTables(idconditions);
}
}
lcsim/src/org/lcsim/mc/fast/reconstructedparticle
diff -u -r1.7 -r1.8
--- MCFastReconstructedParticle.java 20 Aug 2005 21:54:44 -0000 1.7
+++ MCFastReconstructedParticle.java 8 Sep 2006 03:19:45 -0000 1.8
@@ -4,6 +4,7 @@
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;
@@ -12,9 +13,11 @@
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;
/**
@@ -27,26 +30,56 @@
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 _fourVec = new BasicHepLorentzVector();
+ 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)
+ 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);
- double e = sqrt(t.getPX()*t.getPX()+t.getPY()*t.getPY()+t.getPZ()*t.getPZ() + _mass*_mass);
- _fourVec.setV3(e, t.getPX(), t.getPY(), t.getPZ());
_charge = t.getCharge();
-
- // Fixme: Probably wrong, this is not a point on the track, nor should we assume that the track
- // reference point is at the POCA.
- _referencePoint = new BasicHep3Vector(t.getReferencePointX(), t.getReferencePointY(), t.getReferencePointZ());
+ // 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));
}
@@ -66,12 +99,30 @@
double px = (pm/len)*(point[0]);
double py = (pm/len)*(point[1]);
double pz = (pm/len)*(point[2]);
- _fourVec.setV3(e, px, py, pz);
+ 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()
@@ -81,12 +132,12 @@
public Hep3Vector getMomentum()
{
- return _fourVec.v3();
+ return p_reco.v3();
}
public double getEnergy()
{
- return _fourVec.t();
+ return p_reco.t();
}
public double[] getCovMatrix()
@@ -162,7 +213,7 @@
public HepLorentzVector asFourVector()
{
- return _fourVec;
+ return p_reco;
}
public String toString()
lcsim/src/org/lcsim/mc/fast/reconstructedparticle
diff -u -r1.1 -r1.2
--- IDResolutionTables.java 15 Jul 2005 22:09:29 -0000 1.1
+++ IDResolutionTables.java 8 Sep 2006 03:19:45 -0000 1.2
@@ -3,9 +3,7 @@
*
* 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.
+ * @version "$Id: IDResolutionTables.java,v 1.2 2006/09/08 03:19:45 ngraf Exp $"
*/
package org.lcsim.mc.fast.reconstructedparticle;
@@ -18,16 +16,22 @@
*/
public class IDResolutionTables {
- private double ElectronEff;
- private double MuonEff;
- private double NeutronEff;
+ 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()
@@ -37,12 +41,26 @@
public double getMuonEff()
{
- return ElectronEff;
+ return MuonEff;
+ }
+
+ public double getProtonEff()
+ {
+ return ProtonEff;
+ }
+
+ public double getKaonEff()
+ {
+ return KaonEff;
}
public double getNeutronEff()
{
- return ElectronEff;
+ return NeutronEff;
}
+ public double getWtChgTrkCal()
+ {
+ return WtChgTrkCal;
+ }
}
CVSspam 0.2.8
