lcsim/src/org/lcsim/contrib/uiowa/devel
diff -N PFA.java
--- PFA.java 12 Jan 2006 22:21:43 -0000 1.2
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,282 +0,0 @@
-package devel;
-
-import java.util.List;
-import java.util.Vector;
-
-import hep.physics.vec.*;
-
-import org.lcsim.event.EventHeader;
-import org.lcsim.util.Driver;
-import org.lcsim.event.Cluster;
-import org.lcsim.event.MCParticle;
-import org.lcsim.recon.cluster.mst.*;
-
-import structural.*;
-import structural.likelihood.*;
-
-/**
- * An example PFA using the structual algorithm.
- *
- * @version $Id: PFA.java,v 1.2 2006/01/12 22:21:43 mcharles Exp $
- */
-
-public class PFA extends Driver
-{
- // Special treatment for likelihood quantities that need per-event info
- // (e.g. geometry)
- List<StructuralLikelihoodQuantityWithEventInfo> m_perEventQuantities = null;
-
- TestFragmentIdentifier m_wrappedID = null;
- TestFragmentIdentifier m_wrappedID2 = null;
-
- public void suspend()
- {
- if (m_wrappedID != null) {
- m_wrappedID.commit();
- }
- if (m_wrappedID2 != null) {
- m_wrappedID2.commit();
- }
- super.suspend();
- }
-
- public PFA()
- {
- boolean writeLikelihood = false;
-
- // Step 1 (here):
- // For each MC particle, make a helix swimmer.
- // Requires B field.
- // Output: Helix swimmer and MC particle (linked)
- // Step 2:
- // For each helix, extrapolate it to the calorimeter
- // surfaces (barrel and endcap).
- // Requires geometry information.
- // Output: Position and momentum for each particle/swimmer
- // Step 3:
- // For each helix that reaches an ECAL, look for
- // a MIP stub or cluster nearby.
- // Output: Bi-directional maps from MCParticle/swimmer to cluster
- // Step 4:
- // (a) Handle MCParticles which weren't mapped properly
- // (b) Look at E/P for MCParticles which were mapped.
- //
- String nameOfHelixList = new String("HelixList");
- String nameOfExtrapolationInfoList = new String("HelixExtrapolatedToECALList");
- String nameOfHelixToClusterMap = new String("HelixToClusterMap");
- String nameOfClusterToHelixMap = new String("ClusterToHelixMap");
- add(new MakeHelixSwimmersFromTruth(nameOfHelixList)); // step 1 and 2
- add(new SwimToECAL(nameOfHelixList, nameOfExtrapolationInfoList)); // step 2
-
- // Begin with a big-scale cluster set, made with the MST:
- Metrics geomDist = new GeometricalDistance();
- Metrics hitHitDist = new MinimumHitToHitDistance();
-
- MSTClusterDriver mstDriverEcal = new MSTClusterDriver("EMCal");
- MSTClusterDriver mstDriverHcal = new MSTClusterDriver("HCal");
- mstDriverEcal.registerMetrics(geomDist);
- mstDriverHcal.registerMetrics(geomDist);
- mstDriverEcal.setThreshold(30.0); // 3cm
- mstDriverHcal.setThreshold(100.0); // 10cm
- mstDriverEcal.setClusterName("MSTCluster EMCal");
- mstDriverHcal.setClusterName("MSTCluster HCal");
- add(mstDriverEcal);
- add(mstDriverHcal);
-
- // Find MIPs within clusters. Have to do this before merging ECAL and HCAL
- // because the MIP-finder is layer-based -- we don't want to confuse layer n
- // of the ECAL with layer n of the HCAL.
- TrackSegmentFinder findTracksEcal = new TrackSegmentFinder("MSTCluster EMCal", "Track segments EMCal");
- TrackSegmentFinder findTracksHcal = new TrackSegmentFinder("MSTCluster HCal", "Track segments HCal");
- add(findTracksEcal);
- add(findTracksHcal);
-
- // Try to link the tracks from the tracking system to MIP segments and/or clusters:
- add(new MatchHelixToCluster(nameOfExtrapolationInfoList, "Track segments EMCal", "MSTCluster EMCal", nameOfHelixToClusterMap, nameOfClusterToHelixMap));
-
- // Now, link them together across the ECAL-HCAL boundary:
- MSTClusterDriver mstDriverLink = new MSTClusterDriver("User");
- mstDriverLink.registerMetrics(hitHitDist);
- mstDriverLink.setThreshold(50.0); // 5cm
- mstDriverLink.addUserInputList("MSTCluster EMCal");
- mstDriverLink.addUserInputList("MSTCluster HCal");
- mstDriverLink.setClusterName("MSTCluster linked");
- add(mstDriverLink);
-
- // Special thing: Map MIPs to clusters correctly after this cluster linking...
- Remapper<Cluster> remapMIPs = new Remapper<Cluster>("MSTCluster linked", "Track segments linked");
- remapMIPs.addInputClusters("MSTCluster EMCal", "Track segments EMCal");
- remapMIPs.addInputClusters("MSTCluster HCal", "Track segments HCal");
- remapMIPs.setDebug(true);
- add(remapMIPs);
- String nameOfClusterToHelixMapLinked = new String("ClusterToHelixMapLinked");
- Remapper<TrackExtrapolationInfo> remapTracks = new Remapper<TrackExtrapolationInfo>("MSTCluster linked", nameOfClusterToHelixMapLinked);
- remapTracks.addInputClusters("MSTCluster EMCal", nameOfClusterToHelixMap);
- remapTracks.addInputClusters("MSTCluster HCal", nameOfClusterToHelixMap);
- remapTracks.setDebug(true);
- add(remapTracks);
-
- // Find clumps within clusters
- ClumpFinder findClumps = new ClumpFinder("MSTCluster linked", "Track segments linked", "Clumps");
- add(findClumps);
-
- // Run likelihood structural analysis
- if (true) {
- ClusterAssociator assoc = new ClusterEnergyAssociator();
- if (writeLikelihood) {
- // Obtain and write likelihood histograms
- System.out.println("ExamplePFA: I will obtain and write out likelihood histograms.");
- LikelihoodEvaluator eval = new LikelihoodEvaluator();
- eval.addLikelihoodQuantityTrackToTrack(new TrackToTrackDOCA(), 50, 0.0, 100.0, false, true);
- eval.addLikelihoodQuantityTrackToTrack(new TrackToTrackPOCAInCalorimeter(), 2, -0.5, 1.5, false, false);
- eval.addLikelihoodQuantityTrackToTrack(new TrackToTrackSmallestDistanceToPOCA(), 25, 0.0, 250.0, false, true);
- //eval.addLikelihoodQuantityTrackToTrack(new TrackToTrackIntermediateHitsCount(), 10, -0.5, 9.5, false, true);
- //eval.addLikelihoodQuantityTrackToTrack(new TrackToTrackIntermediateHitsFraction(), 11, -0.05, 1.05, false, false);
- eval.addLikelihoodQuantityTrackToClump(new TrackToClumpDOCA(), 50, 0.0, 300.0, false, true);
- eval.addLikelihoodQuantityTrackToClump(new ClusterToClusterMinDistance(), 25, 0.0, 250.0, false, true);
- eval.addLikelihoodQuantityClumpToClump(new ClumpToClumpDOCA(), 20, 0.0, 200.0, false, true);
- eval.addLikelihoodQuantityClumpToClump(new ClusterToClusterMinDistance(), 20, 0.0, 200.0, false, true);
- // Handle things that have per-event info:
- makeEventInfoList(eval);
-
- LikelihoodFindingStructuralDriver likelihoodWriter = new LikelihoodFindingStructuralDriver(eval, assoc, "MSTCluster linked", "Track segments linked", "Clumps");
- likelihoodWriter.setIgnoreClusterDecision(new ClusterSizeDecision(10));
- add(likelihoodWriter);
- Driver checkpoint = new LikelihoodEvaluatorCheckpointDriver(eval, 10);
- add(checkpoint);
- } else {
- // Use pre-existing likelihood histograms to check clusters.
- // Output is, for each large cluster, a list of skeletons clusters inside it.
- System.out.println("ExamplePFA: I will read in likelihood histograms and use them to make clusters.");
- LikelihoodEvaluator eval = LikelihoodEvaluator.readFromFile("likelihood.bin");
- LikelihoodLinkPlotterDriver likelihoodPlotter = new LikelihoodLinkPlotterDriver(eval, 0.5, 0.6, 0.8, assoc, "MSTCluster linked", "Track segments linked", "Clumps", "MapClustersToSkeletons");
- likelihoodPlotter.setIgnoreClusterDecision(new ClusterSizeDecision(10));
- likelihoodPlotter.initPlots("likelihoodPerformance.aida");
- add(likelihoodPlotter);
- // Some likelihood quantities need per-event info:
- makeEventInfoList(eval);
- // Assign hits that are within the local cluster but not part of a clump/MIP.
- add (new HaloAssigner("MapClustersToSkeletons"));
- // Output is a list of clusters where:
- // A large cluster with multiple skeletons has a separate entry for each skeleton (+ associated halo)
- // A large cluster with one skeleton will have one entry
- // A cluster with no skeletons will have one entry
- // The total number of hits is the same as the total number of hits in "MSTCluster linked"
- add (new MakeSeparatedClusters("MSTCluster linked", "MapClustersToSkeletons", "MSTCluster separated"));
-
- // Or, cheat like a bandit
- add (new StructuralCheater("MSTCluster linked", "MSTCluster separated cheated"));
-
- // Handle fragments:
- FragmentIdentifier nonCheatID = new SimpleFragmentIdentifier(nameOfHelixToClusterMap);
- FragmentIdentifier cheatID = new CheatFragmentIdentifier("MSTCluster separated");
- FragmentIdentifier cheatID2 = new CheatFragmentIdentifier("MSTCluster separated cheated");
- m_wrappedID = new TestFragmentIdentifier(nonCheatID, cheatID);
- //m_wrappedID2 = new TestFragmentIdentifier(nonCheatID, cheatID2);
- add (new FragmentMerger("MSTCluster separated", "MSTCluster fragments merged", m_wrappedID)); // or nonCheatID
- //add (new FragmentMerger("MSTCluster separated", "MSTCluster fragments merged", nonCheatID));
- add (new FragmentMerger("MSTCluster separated", "MSTCluster fragments merged-cheat", cheatID));
- //add (new TestFragmentMerger("MSTCluster separated", "MSTCluster fragments merged-cheat", cheatID, new CheatFragmentMerger("MSTCluster separated", "MSTCluster fragments cheated xxx", cheatID))); // test!
- add (new FragmentRemover("MSTCluster separated", "MSTCluster fragments removed", nonCheatID));
- add (new FragmentRemover("MSTCluster separated", "MSTCluster fragments removed-cheat", cheatID));
- add (new CheatFragmentMerger("MSTCluster separated", "MSTCluster fragments cheated", cheatID));
- add (new FragmentMerger("MSTCluster separated cheated", "MSTCluster fragments merged (cheated on separation)", nonCheatID));
- //add (new FragmentMerger("MSTCluster separated cheated", "MSTCluster fragments merged (cheated on separation)", m_wrappedID2));
-
- // When done, check the total energy in the event
- /*
- add (new EventEnergySum("MSTCluster fragments merged", nameOfClusterToHelixMapLinked, nameOfHelixToClusterMap, "EnergySumHistos-merged.aida"));
- add (new EventEnergySum("MSTCluster fragments merged-cheat", nameOfClusterToHelixMapLinked, nameOfHelixToClusterMap, "EnergySumHistos-merged-cheat.aida"));
- add (new EventEnergySum("MSTCluster fragments removed", nameOfClusterToHelixMapLinked, nameOfHelixToClusterMap, "EnergySumHistos-removed.aida"));
- add (new EventEnergySum("MSTCluster fragments removed-cheat", nameOfClusterToHelixMapLinked, nameOfHelixToClusterMap, "EnergySumHistos-removed-cheat.aida"));
- add (new EventEnergySum("MSTCluster fragments cheated", nameOfClusterToHelixMapLinked, nameOfHelixToClusterMap, "EnergySumHistos-cheated.aida"));
- add (new EventEnergySum("MSTCluster linked", nameOfClusterToHelixMapLinked, nameOfHelixToClusterMap, "EnergySumHistos-all.aida"));
- //add (new test.Test("MSTCluster linked")); // TEST!
- add (new EventEnergySum("MSTCluster fragments merged (cheated on separation)", nameOfClusterToHelixMapLinked, nameOfHelixToClusterMap, "EnergySumHistos-merged-cheated-on-separation.aida"));
- */
-
- List<String> knownClusterLists = new Vector<String>();
- knownClusterLists.add("MSTCluster EMCal");
- knownClusterLists.add("MSTCluster HCal");
- knownClusterLists.add("MSTCluster linked");
- knownClusterLists.add("MSTCluster separated");
-
-// add (new CheckStatusOfHitList("EcalBarrHits"));
-// add (new CheckStatusOfHitList("EcalEndcapHits"));
-// add (new CheckStatusOfHitList("HcalBarrHits"));
-// add (new CheckStatusOfHitList("HcalEndcapHits"));
-// add (new CheckStatusOfClusterList("MSTCluster EMCal", knownClusterLists));
-// add (new CheckStatusOfClusterList("MSTCluster HCal", knownClusterLists));
-// add (new CheckStatusOfClusterList("MSTCluster linked", knownClusterLists));
-// add (new CheckStatusOfClusterList("MSTCluster separated", knownClusterLists));
- }
- }
- }
-
- int m_count = 0;
- public void process(EventHeader event) {
- System.out.println("DEBUG: "+this.getClass().getName()+": Event "+m_count+":"); m_count++;
- // Special handling of things that need per-event info:
- if (m_perEventQuantities != null) {
- for (StructuralLikelihoodQuantityWithEventInfo quant : m_perEventQuantities) {
- quant.setEventInfo(event);
- }
- }
- // Here we can do a veto
- if (checkEnergyBarrel(event, 0.9)) {
- if (checkEnergyNeutrinos(event, 100.0)) {
- super.process(event);
- }
- }
- }
-
- protected boolean checkEnergyNeutrinos(EventHeader event, double maxNeutrinoEnergy)
- {
- double truthNeutrinoEnergySum = 0.0;
- List<MCParticle> eventMCParticles = event.getMCParticles();
- for (MCParticle p : eventMCParticles) {
- int pdg = p.getPDGID();
- if (pdg==12 || pdg==14 || pdg==16 || pdg==18 || pdg==-12 || pdg==-14 || pdg==-16 || pdg==-18) {
- truthNeutrinoEnergySum += p.getEnergy();
- }
- }
- return (truthNeutrinoEnergySum <= maxNeutrinoEnergy); // <= in case they put 0 for the maximum
- }
-
- protected boolean checkEnergyBarrel(EventHeader event, double threshold)
- {
- // Require that threshold (e.g. 90%) of the final-state particles are within the barrel ( cos(theta) < 0.8 )
- double energySumBarrel = 0.0;
- double energySumNonBarrel = 0.0;
- List<MCParticle> mcps = event.getMCParticles();
- for (MCParticle mcp : mcps) {
- if (mcp.getGeneratorStatus() == mcp.FINAL_STATE) {
- Hep3Vector momentum = mcp.getMomentum();
- double cosTheta = momentum.z() / momentum.magnitude();
- if (Math.abs(cosTheta) < 0.8) {
- // barrel
- energySumBarrel += mcp.getEnergy();
- } else {
- // non-barrel
- energySumNonBarrel += mcp.getEnergy();
- }
- }
- }
- double energySumTotal = energySumBarrel + energySumNonBarrel;
- return (energySumBarrel / energySumTotal > threshold);
- }
-
- protected void makeEventInfoList(LikelihoodEvaluator eval)
- {
- // Handle things that have per-event info:
- m_perEventQuantities = new Vector<StructuralLikelihoodQuantityWithEventInfo>();
- List<StructuralLikelihoodQuantity> quantities = eval.getLikelihoodQuantities();
- for (StructuralLikelihoodQuantity quant : quantities) {
- if (quant instanceof StructuralLikelihoodQuantityWithEventInfo) {
- StructuralLikelihoodQuantityWithEventInfo quantTmp = (StructuralLikelihoodQuantityWithEventInfo) (quant);
- m_perEventQuantities.add(quantTmp);
- }
- }
- }
-
-}
