lcsim/src/org/lcsim/recon/cheater
diff -u -r1.7 -r1.8
--- ReconCheater.java 1 Nov 2005 16:44:37 -0000 1.7
+++ ReconCheater.java 2 Nov 2005 23:15:45 -0000 1.8
@@ -71,13 +71,13 @@
double ECalResolution, HCalResolution;
double ECalSampling, HCalSampling, HCalDigital; // Hits/GeV
double pTrackMin, EClusterMin;
- double ECalEnergyMin, HCalEnergyMin, NDigitalMin;
+ double ECalEnergyMin, HCalEnergyMin, HitEnergyMin = 0.00005, NDigitalMin;
boolean usePerfectEnergyFlow;
double Distance2XCluster, Distance4XCluster;
- static boolean hist = false;
+ boolean hist = false, calibrate = false;
/** */
@@ -102,6 +102,7 @@
}
public void setHist(boolean hist) { this.hist = hist; initHist(); }
+ public void setCalibrate(boolean calibrate) { this.calibrate = calibrate; }
void init()
@@ -134,6 +135,7 @@
EClusterMin = cheatingTable.getEClusterMin();
ECalEnergyMin = cheatingTable.getECalEnergyMin();
HCalEnergyMin = cheatingTable.getHCalEnergyMin();
+ HitEnergyMin = cheatingTable.getHitEnergyMin();
NDigitalMin = cheatingTable.getNDigitalMin();
usePerfectEnergyFlow = cheatingTable.usePerfectEnergyFlow();
@@ -163,7 +165,7 @@
protected int nEvt = 0;
boolean first = true, firstEvents = true;
- boolean debug = true;
+ boolean debug = false;
int imcp = -1;
/** Reconstruct an event by using the Monte Carlo truth to cheat. */
@@ -693,7 +695,7 @@
if (Klong || neutron) {
Hep3Vector endPoint = mcp.getEndPoint();
double energy = getMeasuredEnergy(cluster);
- if (energy<EClusterMin) continue;
+ if (energy<HCalEnergyMin) continue;
// Check for decay or interaction products.
MCParticle primaryParent = findPrimaryParent(mcp); primaryMCParents.put(mcp,primaryParent);
@@ -712,7 +714,7 @@
MCParticle mcp2 = cluster2.getMCParticle(); int PDGID2 = mcp2.getPDGID();
Hep3Vector endPoint2 = mcp2.getEndPoint();
double energy2 = getMeasuredEnergy(cluster2);
- if (energy2<EClusterMin) continue;
+ if (energy2<HCalEnergyMin) continue;
double dx = endPoint2.x()-endPoint.x(), dy = endPoint2.y()-endPoint.y(), dz = endPoint2.z()-endPoint.z();
double d = Math.sqrt(dx*dx + dy*dy + dz*dz);
@@ -728,6 +730,9 @@
}
+ int NClusterHist = 0;
+ int MaxClusterHist = 5000;
+
void analyzeReconstructedClusters(List<CheatCluster> clusterList)
{
AIDA aida = AIDA.defaultInstance();
@@ -743,7 +748,7 @@
if (!Klong && !neutron) continue;
Hep3Vector endPoint = mcp.getEndPoint();
double energy = getMeasuredEnergy(cluster);
- if (energy<EClusterMin) continue;
+ if (energy<HCalEnergyMin) continue;
// Check for decay or interaction products.
MCParticle primaryParent = findPrimaryParent(mcp); primaryMCParents.put(mcp,primaryParent);
@@ -756,6 +761,14 @@
continue;
nNeutrals++;
+ // Neutral clusters
+ if (calibrate && NClusterHist<MaxClusterHist) {
+ for (CalorimeterHit hit : cluster.getCalorimeterHits()) {
+ NClusterHist++;
+ aida.cloud1D("hit: energy").fill(hit.getRawEnergy());
+ }
+ }
+
// Plot distance to nearest large cluster.
double distance = 1000., distanceLarge2XCluster = 1000., distanceLarge4XCluster = 1000.;
for (CheatCluster cluster2 : clusterList) {
@@ -763,7 +776,7 @@
MCParticle mcp2 = cluster2.getMCParticle(); int PDGID2 = mcp2.getPDGID();
Hep3Vector endPoint2 = mcp2.getEndPoint();
double energy2 = getMeasuredEnergy(cluster2);
- if (energy2<EClusterMin) continue;
+ if (energy2<HCalEnergyMin) continue;
double dx = endPoint2.x()-endPoint.x(), dy = endPoint2.y()-endPoint.y(), dz = endPoint2.z()-endPoint.z();
double d = Math.sqrt(dx*dx + dy*dy + dz*dz);
@@ -797,7 +810,18 @@
aida.cloud1D("# particles").fill(particles.size());
aida.cloud1D("Final Energy").fill(finalEnergy);
aida.cloud1D("Energy diff").fill(finalEnergy-totalEventEnergy);
-
+/*
+ // Reconstructed charged particles.
+ for (MCParticle mcp : charged.keySet()) {
+ Hep3Vector pV = mcp.getMomentum();
+ double p = Math.sqrt(pV.x()*pV.x() + pV.y()*pV.y() + pV.z()*pV.z());
+ //int PDGID = mcp.getPDGID();
+ CheatTrack track = charged.get(mcp);
+ double[] P = track.getMomentum();
+ double momentum = Math.sqrt(P[0]*P[0] + P[1]*P[1] + P[2]*P[2]);
+ aida.cloud2D("charged Mom. diff. vs. Mom.").fill(momentum-p,p);
+ }
+*/
// Reconstructed neutrals.
for (MCParticle mcp : neutrals.keySet()) {
double e = mcp.getEnergy();
@@ -805,27 +829,46 @@
CheatCluster cluster = neutrals.get(mcp);
double eRatio = cluster.getRawEnergy()/mcp.getEnergy();
int nhits = 0;
- for (CalorimeterHit hit : cluster.getCalorimeterHits())
+ for (CalorimeterHit hit : cluster.getCalorimeterHits()) {
+ if (hit.getRawEnergy()<HitEnergyMin) continue;
if (hit.getTime()<100.) nhits++;
+ }
double ratio = nhits/e;
+ double energy = getMeasuredEnergy(cluster);
+ double diff = Math.sqrt(e) * (energy - e) / e;
if (PDGID==22) {
+ if (calibrate) {
aida.cloud1D("photon ratio").fill(eRatio);
aida.cloud2D("photon Energy vs. ratio").fill(eRatio,e);
+ }
+ if (Math.abs(diff)>1.0) continue;
+ if (e>1. && e<10.) aida.cloud1D("photon Energy diff.").fill(diff);
+ aida.cloud2D("photon Energy diff. vs. Energy").fill(diff,e);
}
else if (PDGID==130) {
+ if (calibrate) {
aida.cloud1D("Klong E ratio").fill(eRatio);
aida.cloud2D("Klong Energy vs. E ratio").fill(eRatio,e);
aida.cloud1D("Klong Hit ratio").fill(ratio);
aida.cloud2D("Klong Energy vs. Hit ratio").fill(ratio,e);
if (e<1.0) aida.cloud1D("Klong low Energy").fill(e);
+ }
+ if (Math.abs(diff)>2.0) continue;
+ if (e>1. && e<10.) aida.cloud1D("Klong Energy diff.").fill(diff);
+ aida.cloud2D("Klong Energy diff. vs. Energy").fill(diff,e);
}
else if (Math.abs(PDGID)==2112) {
+ if (calibrate) {
aida.cloud1D("neutron E ratio").fill(eRatio);
aida.cloud2D("neutron Energy vs. E ratio").fill(eRatio,e);
aida.cloud1D("neutron Hit ratio").fill(ratio);
aida.cloud2D("neutron Energy vs. Hit ratio").fill(ratio,e);
if (e<1.0) aida.cloud1D("neutron low Energy").fill(e);
+ }
+ if (Math.abs(diff)>3.0) continue;
+ if (e>1. && e<10.) aida.cloud1D("neutron Energy diff.").fill(diff);
+ aida.cloud2D("neutron Energy diff. vs. Energy").fill(diff,e);
}
}
@@ -1241,8 +1284,10 @@
}
else {
int nhits = 0;
- for (CalorimeterHit hit : cluster.getCalorimeterHits())
+ for (CalorimeterHit hit : cluster.getCalorimeterHits()) {
+ if (hit.getRawEnergy()<HitEnergyMin) continue;
if (hit.getTime()<100.) nhits++;
+ }
energy = nhits/HCalDigital;
if (nhits<NDigitalMin || energy<HCalEnergyMin) energy = 0.;
}