lcsim-contrib/src/main/java/org/lcsim/contrib/jeremym/pfa/cheat
diff -u -r1.1 -r1.2
--- PerfectClusteringDriver.java 2 Oct 2009 01:10:03 -0000 1.1
+++ PerfectClusteringDriver.java 3 Oct 2009 00:06:36 -0000 1.2
@@ -8,7 +8,7 @@
import org.lcsim.event.EventHeader;
import org.lcsim.event.MCParticle;
import org.lcsim.event.SimCalorimeterHit;
-import org.lcsim.recon.cluster.util.BasicCluster;
+import org.lcsim.recon.cluster.cheat.CheatCluster;
import org.lcsim.util.Driver;
/**
@@ -169,7 +169,7 @@
if (contribHits.size() >= this.minHits)
{
// Create a new cluster.
- BasicCluster perfectCluster = new BasicCluster();
+ CheatCluster perfectCluster = new CheatCluster(particle);
// Add contributing hit list to cluster.
for (CalorimeterHit hit : contribHits)
lcsim-contrib/src/main/java/org/lcsim/contrib/jeremym/pfa/cheat
diff -u -r1.2 -r1.3
--- SimpleAnalysis.java 2 Oct 2009 01:44:14 -0000 1.2
+++ SimpleAnalysis.java 3 Oct 2009 00:06:36 -0000 1.3
@@ -2,10 +2,15 @@
import java.io.IOException;
+import org.lcsim.event.CalorimeterHit;
import org.lcsim.event.Cluster;
import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
import org.lcsim.event.ReconstructedParticle;
import org.lcsim.event.Track;
+import org.lcsim.event.base.BaseTrackMC;
+import org.lcsim.event.util.ParticleTypeClassifier;
+import org.lcsim.recon.cluster.cheat.CheatCluster;
import org.lcsim.util.Driver;
import org.lcsim.util.aida.AIDA;
@@ -16,6 +21,14 @@
AIDA aida = AIDA.defaultInstance();
int nevents = 0;
double avgE = 0.;
+ double mcpETotal = 0;
+ double allMcpETotal = 0;
+ double eventEnergy = 500; // 500 GeV default event E
+
+ public void setEventEnergy(double eventEnergy)
+ {
+ this.eventEnergy = eventEnergy;
+ }
public SimpleAnalysis()
{}
@@ -32,14 +45,72 @@
public void process(EventHeader event)
{
+ double allMcpESum = 0;
+ for (MCParticle mcp : event.get(MCParticle.class, EventHeader.MC_PARTICLES))
+ {
+ if (mcp.getGeneratorStatus() == mcp.FINAL_STATE)
+ {
+ allMcpESum += mcp.getEnergy();
+ allMcpETotal += mcp.getEnergy();
+ }
+ }
double energySum = 0.;
int nparticles = 0;
int nclusters = 0;
int ntracks = 0;
int ntrackHits = 0;
int nhits = 0;
+ double photonESum = 0;
+ int nphoton = 0;
+ double chargedHadronESum = 0;
+ int nchargedHadron = 0;
+ double neutralHadronESum = 0;
+ int nneutralHadron = 0;
+ double calHitESum = 0;
+ int nTracksLowHits = 0;
+ int nTracksLowHits2 =0;
+ int nTracksLowPt = 0;
+ int nTracksLowPt2 = 0;
+ double muonESum = 0;
+ int nmuon = 0;
+ double electronESum = 0;
+ int nelectron = 0;
+ double mcparticleESum = 0;
+ double mcparticleMassSum = 0;
+ int nclustersLowHits = 0;
+ int nclustersLowHits2 = 0;
+ int nneutrino = 0;
+ int neutrinoESum = 0;
+ double positronESum = 0;
+ int npositron = 0;
for (ReconstructedParticle p : event.get(ReconstructedParticle.class, particleCollection))
{
+ int pdgid = 0;
+ MCParticle mcp = null;
+ // MCParticle from BaseTrackMC.
+ if (p.getTracks().size() != 0)
+ {
+ mcp = ((BaseTrackMC)p.getTracks().get(0)).getMCParticle();
+ }
+ // MCParticle from CheatCluster.
+ else if (p.getClusters().size() != 0)
+ {
+ mcp = ((CheatCluster)p.getClusters().get(0)).getMCParticle();
+ }
+ // No tracks or clusters. This should never happen!
+ else
+ {
+ System.out.println("Skipping particle without tracks or clusters!");
+ continue;
+ }
+ // Set PDG id.
+ pdgid = mcp.getPDGID();
+
+ mcparticleESum += mcp.getEnergy();
+ mcpETotal += mcp.getEnergy();
+ mcparticleMassSum += mcp.getMass();
+
+ double charge = p.getCharge();
energySum += p.getEnergy();
++nparticles;
nclusters += p.getClusters().size();
@@ -48,26 +119,143 @@
for (Cluster c : p.getClusters())
{
nhits += c.getCalorimeterHits().size();
+ if (c.getCalorimeterHits().size() < 3)
+ {
+ ++nclustersLowHits;
+ }
+ else if (c.getCalorimeterHits().size() < 10)
+ {
+ ++nclustersLowHits2;
+ }
nhitsP += c.getCalorimeterHits().size();
+ for (CalorimeterHit calhit : c.getCalorimeterHits())
+ {
+ calHitESum += calhit.getCorrectedEnergy();
+ }
}
int ntrackHitsP = 0;
for (Track t : p.getTracks())
{
ntrackHits += t.getTrackerHits().size();
ntrackHitsP += t.getTrackerHits().size();
+ double mom[] = t.getMomentum();
+ double pt = Math.sqrt(mom[0]*mom[0]+mom[1]*mom[1]);
+ if (pt < 0.2)
+ {
+ ++nTracksLowPt;
+ }
+ else if (pt < 1.0)
+ {
+ ++nTracksLowPt2;
+ }
+ }
+ if (ntrackHitsP < 7)
+ {
+ ++nTracksLowHits;
+ }
+ else if (ntrackHitsP < 3)
+ {
+ ++nTracksLowHits2;
+ }
+ aida.cloud1D("Tracks per Particle").fill(p.getTracks().size());
+ aida.cloud1D("Clusters per Particle").fill(p.getClusters().size());
+ aida.cloud1D("CalHits per Particle").fill(nhitsP);
+ aida.cloud1D("TrackerHits per Particle").fill(ntrackHitsP);
+ aida.cloud1D("Energy per Particle").fill(p.getEnergy());
+ aida.cloud1D("Mass per Particle").fill(p.getMass());
+ if (ParticleTypeClassifier.isPhoton(pdgid))
+ {
+ photonESum += p.getEnergy();
+ ++nphoton;
+ }
+ // FIXME: ParticleTypeClassifier doesn't know K0L or K0S are hadrons.
+ else if (ParticleTypeClassifier.isHadron(pdgid) || pdgid == 130 || pdgid == 310)
+ {
+ if (charge == 0)
+ {
+ neutralHadronESum += p.getEnergy();
+ ++nneutralHadron;
+ }
+ else
+ {
+ chargedHadronESum += p.getEnergy();
+ ++nchargedHadron;
+ }
+ }
+ else if (ParticleTypeClassifier.isMuon(pdgid))
+ {
+ muonESum += p.getEnergy();
+ ++nmuon;
+ }
+ else if (ParticleTypeClassifier.isElectron(pdgid))
+ {
+ electronESum += p.getEnergy();
+ ++nelectron;
+ }
+ else if (ParticleTypeClassifier.isNeutrino(pdgid))
+ {
+ neutrinoESum += p.getEnergy();
+ ++nneutrino;
+ }
+ else if (ParticleTypeClassifier.isPositron(pdgid))
+ {
+ positronESum += p.getEnergy();
+ ++npositron;
+ }
+ else
+ {
+ System.out.println("WARNING - Unhandled PDGid " + pdgid);
}
- aida.cloud1D("Number of Tracks for Particle").fill(p.getTracks().size());
- aida.cloud1D("Number of Clusters for Particle").fill(p.getClusters().size());
- aida.cloud1D("Number of CalorimeterHits for Particle").fill(nhitsP);
- aida.cloud1D("Number of TrackerHits for Particle").fill(ntrackHitsP);
- aida.cloud1D("Particle Energy").fill(p.getEnergy());
}
- aida.cloud1D("Total Energy per Event").fill(energySum);
- aida.cloud1D("Total Number of Particles per Event").fill(nparticles);
- aida.cloud1D("Total Number of Clusters per Event").fill(nclusters);
- aida.cloud1D("Total Number of Tracks per Event").fill(ntracks);
- aida.cloud1D("Total Number of Cal Hits per Event").fill(nhits);
- aida.cloud1D("Total Number of Tracker Hits per Event").fill(ntrackHits);
+ aida.cloud1D("Energy per Event").fill(energySum);
+ aida.cloud1D("Particles per Event").fill(nparticles);
+ aida.cloud1D("Clusters per Event").fill(nclusters);
+ aida.cloud1D("Clusters with < 3 CalHits per Event").fill(nclustersLowHits);
+ aida.cloud1D("Clusters with < 10 CalHits per Event").fill(nclustersLowHits2);
+ aida.cloud1D("Tracks per Event").fill(ntracks);
+ aida.cloud1D("CalHits per Event").fill(nhits);
+ aida.cloud1D("CalHit Energy per Event").fill(calHitESum);
+ aida.cloud1D("Tracker Hits per Event").fill(ntrackHits);
+ aida.cloud1D("Tracks with < 7 TrackerHits per Event").fill(nTracksLowHits);
+ aida.cloud1D("Tracks with < 3 TrackerHits per Event").fill(nTracksLowHits2);
+ aida.cloud1D("Tracks with pT < 200 MeV per Event").fill(nTracksLowPt);
+ aida.cloud1D("Tracks with pT < 1 GeV per Event").fill(nTracksLowPt2);
+ aida.cloud1D("Photon Energy per Event").fill(photonESum);
+ aida.cloud1D("Photons per Event").fill(nphoton);
+ aida.cloud1D("Charged Hadron Energy per Event").fill(chargedHadronESum);
+ aida.cloud1D("Charged Hadrons per Event").fill(nchargedHadron);
+ aida.cloud1D("Neutral Hadron Energy per Event").fill(neutralHadronESum);
+ aida.cloud1D("Neutral Hadrons per Event").fill(nneutralHadron);
+ aida.cloud1D("Photon Energy Fraction").fill(photonESum/energySum);
+ aida.cloud1D("Neutral Hadron Energy Fraction").fill(neutralHadronESum/energySum);
+ aida.cloud1D("Charged Hadron Energy Fraction").fill(chargedHadronESum/energySum);
+ aida.cloud1D("Muon Energy per Event").fill(muonESum);
+ aida.cloud1D("Muons per Event").fill(nmuon);
+ aida.cloud1D("Muon Energy Fraction per Event").fill(muonESum/energySum);
+ aida.cloud1D("Electron Energy per Event").fill(electronESum);
+ aida.cloud1D("Electrons per Event").fill(nelectron);
+ aida.cloud1D("Electron Energy Fraction per Event").fill(electronESum/energySum);
+ aida.cloud1D("Positron Energy per Event").fill(positronESum);
+ aida.cloud1D("Positrons per Event").fill(npositron);
+ aida.cloud1D("Positron Energy Fraction per Event").fill(positronESum/energySum);
+ aida.cloud1D("Reconstructed MCParticle Energy per Event").fill(mcparticleESum);
+ aida.cloud1D("Reconstructed MCParticle Mass Sum per Event").fill(mcparticleMassSum);
+ aida.cloud1D("All FS MCParticle Energy per Event").fill(allMcpESum);
+ aida.cloud1D("All MCParticle Energy Minus Reconstructed MCParticle Energy").fill(allMcpESum - mcparticleESum);
+ double visE =
+ photonESum +
+ neutralHadronESum +
+ chargedHadronESum +
+ muonESum +
+ positronESum +
+ electronESum;
+ aida.cloud1D("Visible Particle Energy per Event").fill(visE);
+ aida.cloud1D("Visible Particle Energy Fraction per Event").fill(visE/energySum);
+ aida.cloud1D("Neutrino Energy Sum per Event").fill(neutrinoESum);
+ aida.cloud1D("Neutrinos per Event").fill(nneutrino);
+ aida.cloud1D("Neutrino Energy Fraction per Event").fill(neutrinoESum/energySum);
+ aida.cloud1D("Visible plus Neutrino Energy per Event").fill(visE + neutrinoESum);
+ aida.cloud1D("Visible plus Neutrino Energy Fraction per Event").fill((visE + neutrinoESum)/energySum);
avgE += energySum;
++nevents;
}
@@ -75,7 +263,14 @@
public void endOfData()
{
avgE = avgE / nevents;
+ aida.cloud1D("Total Events Processed").fill(nevents);
aida.cloud1D("Average Energy for " + nevents + " Events").fill(avgE);
+ double resolution = aida.cloud1D("Energy per Event").rms() / aida.cloud1D("Energy per Event").mean();
+ aida.cloud1D("Simple Energy Resolution for " + nevents + " Events").fill(resolution);
+ aida.cloud1D("Energy Error from FS MCParticle Energy for " + nevents + " Events").fill(avgE - (mcpETotal / nevents));
+ aida.cloud1D("Average Reconstructed MCParticle Energy for " + nevents + " Events").fill(this.mcpETotal / nevents);
+ aida.cloud1D("Average FS MCParticle Energy for " + nevents + " Events").fill(this.allMcpETotal / nevents);
+ aida.cloud1D("Energy Error from Reconstructed MCParticle Energy for " + nevents + " Events").fill(avgE - (mcpETotal / nevents));
try
{
aida.saveAs(filename);