Author: [log in to unmask]
Date: Fri Jan 9 14:44:03 2015
New Revision: 3486
Log:
Import fast mc module from CVS.
Added:
projects/lcsim/trunk/mc/
projects/lcsim/trunk/mc/pom.xml
projects/lcsim/trunk/mc/src/
projects/lcsim/trunk/mc/src/main/
projects/lcsim/trunk/mc/src/main/java/
projects/lcsim/trunk/mc/src/main/java/org/
projects/lcsim/trunk/mc/src/main/java/org/lcsim/
projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/
projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/
projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/MCFast.java
projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/Main.java
projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/cluster/
projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/cluster/ronan/
projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/cluster/ronan/ClusterResolutionTables.java
projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/cluster/ronan/MCFastRonan.java
projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/cluster/ronan/ReconCluster.java
projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/cluster/ronan/ReconEMCluster.java
projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/cluster/ronan/ReconHADCluster.java
projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/reconstructedparticle/
projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/reconstructedparticle/IDResolutionTables.java
projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/reconstructedparticle/MCFastParticleID.java
projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/reconstructedparticle/MCFastReconstructedParticle.java
projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/reconstructedparticle/MCFastReconstructedParticleDriver.java
projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/
projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/DocaTrackParameters.java
projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/LookupTable.java
projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/MCFastTracking.java
projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/ReconTrack.java
projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/ResolutionTable.java
projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/SimpleTables.java
projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/SmearTrack.java
projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/SmearTrackSimple.java
projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/TrackParameters.java
projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/TrackResolutionTables.java
projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/fix/
projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/fix/FastMCTrack.java
projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/fix/FastMCTrackDriver.java
projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/fix/FastMCTrackFactory.java
projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/util/
projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/util/CreateFinalStateMCParticleList.java
projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/util/MCParticleClassifier.java
projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/util/MonitorStdhep.java
Modified:
projects/lcsim/trunk/parent/pom.xml
projects/lcsim/trunk/pom.xml
Added: projects/lcsim/trunk/mc/pom.xml
=============================================================================
--- projects/lcsim/trunk/mc/pom.xml (added)
+++ projects/lcsim/trunk/mc/pom.xml Fri Jan 9 14:44:03 2015
@@ -0,0 +1,43 @@
+<project xmlns="http://maven.apache.org/POM/4.0.0" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://maven.apache.org/POM/4.0.0 http://maven.apache.org/maven-v4_0_0.xsd">
+ <modelVersion>4.0.0</modelVersion>
+ <artifactId>lcsim-mc</artifactId>
+ <name>mc</name>
+ <description>Monte Carlo packages including the MCFast event simulation</description>
+ <parent>
+ <groupId>org.lcsim</groupId>
+ <artifactId>lcsim-parent</artifactId>
+ <version>3.0.8-SNAPSHOT</version>
+ <relativePath>../parent/pom.xml</relativePath>
+ </parent>
+ <scm>
+ <url>http://java.freehep.org/svn/repos/lcdet/list/projects/lcsim/trunk/mc/</url>
+ <connection>scm:svn:svn://svn.freehep.org/lcdet/projects/lcsim/trunk/mc/</connection>
+ <developerConnection>scm:svn:svn://svn.freehep.org/lcdet/projects/lcsim/trunk/mc/</developerConnection>
+ </scm>
+ <dependencies>
+ <dependency>
+ <groupId>org.lcsim</groupId>
+ <artifactId>lcsim-tracking</artifactId>
+ </dependency>
+<!--
+ <dependency>
+ <groupId>org.lcsim</groupId>
+ <artifactId>lcsim-math</artifactId>
+ </dependency>
+ <dependency>
+ <groupId>org.lcsim</groupId>
+ <artifactId>lcsim-aida</artifactId>
+ </dependency>
+-->
+ <dependency>
+ <groupId>org.apache.commons</groupId>
+ <artifactId>commons-math</artifactId>
+ <version>2.2</version>
+ </dependency>
+ <dependency>
+ <groupId>commons-cli</groupId>
+ <artifactId>commons-cli</artifactId>
+ <version>1.0</version>
+ </dependency>
+ </dependencies>
+</project>
Added: projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/MCFast.java
=============================================================================
--- projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/MCFast.java (added)
+++ projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/MCFast.java Fri Jan 9 14:44:03 2015
@@ -0,0 +1,64 @@
+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 java.util.logging.Handler;
+import java.util.logging.Level;
+import java.util.logging.Logger;
+import org.lcsim.mc.fast.util.CreateFinalStateMCParticleList;
+import org.lcsim.util.Driver;
+
+/**
+ *
+ * @author Tony Johnson
+ */
+public class MCFast extends Driver {
+ /** Creates a new instance of MCFast */
+
+ public static Logger log;
+ String FSname = "GenFinalStateParticles";
+
+ public MCFast(boolean beamSpotConstraint, boolean simple, long seed, boolean printinfo, boolean refPoint000) {
+ this(beamSpotConstraint, simple, printinfo, refPoint000);
+ getRandom().setSeed(seed);
+ }
+
+ public MCFast(boolean beamSpotConstraint, boolean simple, long seed, boolean printinfo) {
+ this(beamSpotConstraint, simple, printinfo);
+ getRandom().setSeed(seed);
+ }
+
+ public MCFast(boolean beamSpotConstraint, boolean simple, boolean printinfo) {
+ this(beamSpotConstraint, simple, printinfo, false);
+ }
+
+ public MCFast(boolean beamSpotConstraint, boolean simple, boolean printinfo, boolean refPoint000) {
+ log = getLogger();
+ if (printinfo) {
+ log.setLevel(Level.INFO);
+ } else {
+ log.setLevel(Level.WARNING);
+ }
+ add(new CreateFinalStateMCParticleList("Gen"));
+ MCFastTracking mcft = new MCFastTracking(beamSpotConstraint, simple);
+ mcft.setFSList(FSname);
+ add(mcft);
+ MCFastRonan mcfr = new MCFastRonan();
+ mcfr.setFSList(FSname);
+ add(mcfr);
+ add(new MCFastReconstructedParticleDriver(refPoint000));
+ }
+
+ public MCFast(boolean beamSpotConstraint, boolean simple) {
+ this(beamSpotConstraint, simple, false);
+ }
+
+ public MCFast(boolean beamSpotConstraint, boolean simple, long seed) {
+ this(beamSpotConstraint, simple, seed, false);
+ }
+
+ public MCFast() {
+ this(false, false);
+ }
+}
Added: projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/Main.java
=============================================================================
--- projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/Main.java (added)
+++ projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/Main.java Fri Jan 9 14:44:03 2015
@@ -0,0 +1,196 @@
+package org.lcsim.mc.fast;
+
+import java.io.File;
+import java.io.IOException;
+import java.util.Random;
+
+import org.apache.commons.cli.CommandLine;
+import org.apache.commons.cli.HelpFormatter;
+import org.apache.commons.cli.Options;
+import org.apache.commons.cli.ParseException;
+import org.apache.commons.cli.PosixParser;
+import org.freehep.record.source.NoSuchRecordException;
+import org.lcsim.geometry.compact.Detector;
+import org.lcsim.geometry.util.DetectorLocator;
+import org.lcsim.mc.fast.MCFast;
+import org.lcsim.util.Driver;
+import org.lcsim.util.loop.LCIODriver;
+import org.lcsim.util.loop.LCSimLoop;
+import org.lcsim.util.loop.LCSimConditionsManagerImplementation;
+import org.lcsim.mc.fast.util.MonitorStdhep;
+
+/**
+ * This is a reworked class from an old main program provided by Tim and Norman. The jet and vertex reconstruction was removed, so this now only runs the Fast MC and writes the output LCIO file.
+ * Analysis and reconstruction would be performed in a subsequent lcsim XML job.
+ *
+ * @author Norman A. Graf
+ * @author Jeremy McCormick
+ */
+public class Main {
+ private static final String defaultDetector = "sidloi3";
+ private static final int defaultNumProcess = -1;
+ private static final int defaultNumSkip = -1;
+ private static final boolean defaultRefPoint000 = true;
+ private static final String defaultOutputFileName = "fastmc.slcio";
+ private static final boolean defaultBeamSpotConstraint = true;
+ private static final boolean defaultSimple = true;
+ private static final int defaultSeed = (new Random()).nextInt();
+ private static final boolean defaultDebug = false;
+
+ String inputFileName;
+ String outputFileName;
+ String detector;
+ int numToProcess;
+ int numToSkip;
+ int seed;
+ boolean refPoint000;
+ boolean beamSpotConstraint;
+ boolean simple;
+ boolean debug;
+
+ private static Options options = null;
+
+ private Main() {
+ createOptions();
+ }
+
+ Options getOptions() {
+ return options;
+ }
+
+ void parse(String[] args) {
+
+ // Create the parser and parse command line options.
+ PosixParser parser = new PosixParser();
+
+ CommandLine cmd = null;
+ try {
+ cmd = parser.parse(options, args);
+ } catch (ParseException x) {
+ System.out.println("Parsing failed: " + x.getMessage());
+ usage();
+ }
+
+ // Set the parameters from the parsed command line.
+ setParameters(cmd);
+
+ // Print out the parameters to System.out.
+ printParameters();
+ }
+
+ private void setParameters(CommandLine cmd) {
+ inputFileName = cmd.getOptionValue("i");
+ if (inputFileName == null)
+ usage();
+ outputFileName = cmd.getOptionValue("o", defaultOutputFileName);
+ detector = cmd.getOptionValue("d", defaultDetector);
+ numToProcess = Integer.valueOf(cmd.getOptionValue("r", String.valueOf(defaultNumProcess)));
+ numToSkip = Integer.valueOf(cmd.getOptionValue("s", String.valueOf(defaultNumSkip)));
+ seed = Integer.valueOf(cmd.getOptionValue("m", String.valueOf(defaultSeed)));
+ refPoint000 = Boolean.valueOf(cmd.getOptionValue("p", String.valueOf(defaultRefPoint000)));
+ beamSpotConstraint = Boolean.valueOf(cmd.getOptionValue("b", String.valueOf(defaultBeamSpotConstraint)));
+ simple = Boolean.valueOf(cmd.getOptionValue("S", String.valueOf(defaultSimple)));
+ debug = Boolean.valueOf(cmd.getOptionValue("v", String.valueOf(defaultDebug)));
+ }
+
+ private void printParameters() {
+ System.out.println("Received the following command line parameters:");
+ System.out.println('\t' + "inputFileName = " + inputFileName);
+ System.out.println('\t' + "outputFileName = " + outputFileName);
+ System.out.println('\t' + "detector = " + detector);
+ System.out.println('\t' + "events to process = " + numToProcess);
+ System.out.println('\t' + "events to skip = " + numToSkip);
+ System.out.println('\t' + "seed = " + seed);
+ System.out.println('\t' + "refPoint000 = " + refPoint000);
+ System.out.println('\t' + "beamSpotConstraint = " + beamSpotConstraint);
+ System.out.println('\t' + "simple mode = " + simple);
+ System.out.println('\t' + "debug mode = " + debug);
+ }
+
+ private void createOptions() {
+ options = new Options();
+ options.addOption("h", false, "print usage information");
+ options.addOption("i", true, "input file");
+ options.addOption("o", true, "output file name");
+ options.addOption("d", true, "detector name");
+ options.addOption("r", true, "number of events to process");
+ options.addOption("s", true, "number of events to skip");
+ options.addOption("m", true, "random seed");
+ options.addOption("p", false, "use default ref point");
+ options.addOption("b", false, "use beam spot constraint");
+ options.addOption("S", false, "use simple");
+ options.addOption("v", false, "print debug info");
+ }
+
+ private void usage() {
+ HelpFormatter formatter = new HelpFormatter();
+ formatter.printHelp(getClass().getCanonicalName(), options);
+ System.exit(1);
+ }
+
+ private void error(String message) {
+ System.out.println(message);
+ System.exit(1);
+ }
+
+ private void run() throws IOException, NoSuchRecordException {
+ // Check existence of detector.
+ Detector det = DetectorLocator.findDetector(detector);
+ if (det == null) {
+ error("Unknown detector: " + detector);
+ }
+
+ // Check existence of input file.
+ File input = new File(inputFileName);
+ if (!input.exists()) {
+ error("The input file " + input + " does not exist!");
+ }
+
+ // Setup the LCIO output driver.
+ LCIODriver writer = null;
+ File output = new File(outputFileName);
+ if (output.exists()) {
+ throw new RuntimeException("Output file already exists!");
+ }
+ writer = new LCIODriver(output);
+
+ // Initialize Fast MC driver.
+ Driver fast = new MCFast(beamSpotConstraint, simple, seed, debug, refPoint000);
+
+ // create the event loop
+ LCSimConditionsManagerImplementation.register();
+ LCSimLoop loop = new LCSimLoop();
+ if (input.getName().contains(".stdhep")) {
+ loop.setStdhepRecordSource(input, detector);
+ } else {
+ loop.setLCIORecordSource(input);
+ }
+
+ // Add drivers.
+ if (debug) {
+ MonitorStdhep analysis = new MonitorStdhep();
+ loop.add(analysis);
+ }
+ loop.add(fast);
+ loop.add(writer);
+
+ // Run the job.
+ if (numToSkip > 0) {
+ System.out.println("skipping " + numToSkip + " events");
+ loop.skip(numToSkip);
+ }
+ loop.loop(numToProcess);
+ loop.dispose();
+ System.out.println(getClass().getSimpleName() + ": processed " + loop.getTotalSupplied() + " events");
+ }
+
+ public static void main(String[] args) {
+ Main main = new Main();
+ main.parse(args);
+ try {
+ main.run();
+ } catch (Exception e) {
+ throw new RuntimeException(e);
+ }
+ }
+}
Added: projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/cluster/ronan/ClusterResolutionTables.java
=============================================================================
--- projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/cluster/ronan/ClusterResolutionTables.java (added)
+++ projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/cluster/ronan/ClusterResolutionTables.java Fri Jan 9 14:44:03 2015
@@ -0,0 +1,131 @@
+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) {
+ 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;
+ }
+}
Added: projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/cluster/ronan/MCFastRonan.java
=============================================================================
--- projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/cluster/ronan/MCFastRonan.java (added)
+++ projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/cluster/ronan/MCFastRonan.java Fri Jan 9 14:44:03 2015
@@ -0,0 +1,156 @@
+package org.lcsim.mc.fast.cluster.ronan;
+
+import java.util.ArrayList;
+import java.util.List;
+import java.util.Random;
+
+import org.lcsim.mc.fast.MCFast;
+import org.lcsim.conditions.ConditionsEvent;
+import org.lcsim.conditions.ConditionsListener;
+import org.lcsim.conditions.ConditionsSet;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.LCRelation;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.util.Driver;
+import org.lcsim.event.base.MyLCRelation;
+
+/**
+ * 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 boolean defaultMC = true;
+ private String fsname;
+
+ private ClusterResolutionTables clusterParm;
+
+ public void setFSList(String fslist) {
+ fsname = fslist;
+ defaultMC = false;
+ }
+
+ protected void process(EventHeader event) {
+ if (defaultMC) {
+ fsname = "MCParticle";
+ } else {
+ if (!event.hasCollection(MCParticle.class, fsname)) {
+ System.err.println("Collection " + fsname + " not found. Default Final State particles being used");
+ fsname = "MCParticle";
+ }
+ }
+ if (clusterParm == null) {
+ ConditionsSet conditions = getConditionsManager().getConditions("ClusterParameters");
+ conditions.addConditionsListener(this);
+ clusterParm = new ClusterResolutionTables(conditions);
+ }
+ List<Cluster> cl = new ArrayList<Cluster>();
+ List<LCRelation> lcrelationList = new ArrayList<LCRelation>();
+
+ boolean hist = getHistogramLevel() > 0;
+
+ List<MCParticle> particles = event.get(MCParticle.class, fsname);
+ for (MCParticle p : particles) {
+
+ // filter for FINALSTATE
+ if (defaultMC) {
+ 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();
+ if (Double.isNaN(E)) {
+ continue;
+ }
+
+ 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 (absPDGID == PhotonID || absPDGID == 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));
+ lcrelationList.add(new MyLCRelation(cl.get(cl.size() - 1), (MCParticle) p));
+
+ }
+
+ // Neutral hadrons
+ else if (absPDGID != 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));
+ lcrelationList.add(new MyLCRelation(cl.get(cl.size() - 1), (MCParticle) p));
+ }
+ }
+ 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());
+ }
+ MCFast.log.info(" 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);
+ event.put("ClustersToMCP", lcrelationList, LCRelation.class, 0);
+ }
+
+ public void conditionsChanged(ConditionsEvent event) {
+ ConditionsSet conditions = getConditionsManager().getConditions("ClusterParameters");
+ clusterParm = new ClusterResolutionTables(conditions);
+ }
+}
Added: projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/cluster/ronan/ReconCluster.java
=============================================================================
--- projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/cluster/ronan/ReconCluster.java (added)
+++ projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/cluster/ronan/ReconCluster.java Fri Jan 9 14:44:03 2015
@@ -0,0 +1,189 @@
+package org.lcsim.mc.fast.cluster.ronan;
+
+import org.lcsim.mc.fast.MCFast;
+import hep.physics.particle.Particle;
+import java.util.Collections;
+import java.util.List;
+import java.util.Random;
+import org.lcsim.event.Cluster;
+
+public 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 energyError;
+ protected double neg_energy;
+ protected double sigma;
+ protected double phi;
+ protected double radius;
+ protected double theta;
+ protected double transDist;
+
+ 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 getEnergyError() {
+ return energyError;
+ }
+
+ public void setEnergyError(double energyError) {
+ this.energyError = energyError;
+ }
+
+ public double getNegEnergy() {
+ return neg_energy;
+ }
+
+ public double getSigma() {
+ return sigma;
+ }
+
+ public void adjustEnergy(double neg_energy_total, double pos_energy_weight_total) {
+ MCFast.log.info(" 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;
+
+ MCFast.log.info(" 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 = Math.sqrt(Math.pow(a, 2.) * E + Math.pow(b * E, 2.));
+
+ 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 getParticleId() {
+ return 0;
+ }
+
+ public int getSize() {
+ return 0;
+ }
+}
Added: projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/cluster/ronan/ReconEMCluster.java
=============================================================================
--- projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/cluster/ronan/ReconEMCluster.java (added)
+++ projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/cluster/ronan/ReconEMCluster.java Fri Jan 9 14:44:03 2015
@@ -0,0 +1,26 @@
+package org.lcsim.mc.fast.cluster.ronan;
+
+import org.lcsim.util.aida.AIDA;
+import hep.physics.particle.Particle;
+import java.util.Random;
+
+public class ReconEMCluster extends ReconCluster {
+ ReconEMCluster(ClusterResolutionTables parm, Random rand, Particle mcp, boolean hist) {
+ super(parm, rand, mcp, hist);
+
+ a = parm.getEMResolution();
+ b = parm.getEMConstantTerm();
+ c = parm.getEMPositionError();
+ d = parm.getEMAlignmentError();
+
+ smear(rand, hist);
+ }
+
+ protected void smearPosition(Random rand, double E, boolean hist) {
+ double transSigma = c / Math.sqrt(E) + d;
+ transDist = transSigma * rand.nextGaussian();
+ if (hist)
+ AIDA.defaultInstance().cloud1D("EM: transDist").fill(transDist);
+ smearPosition(rand);
+ }
+}
Added: projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/cluster/ronan/ReconHADCluster.java
=============================================================================
--- projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/cluster/ronan/ReconHADCluster.java (added)
+++ projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/cluster/ronan/ReconHADCluster.java Fri Jan 9 14:44:03 2015
@@ -0,0 +1,26 @@
+package org.lcsim.mc.fast.cluster.ronan;
+
+import org.lcsim.util.aida.AIDA;
+import hep.physics.particle.Particle;
+import java.util.Random;
+
+public class ReconHADCluster extends ReconCluster {
+ ReconHADCluster(ClusterResolutionTables parm, Random rand, Particle mcp, boolean hist) {
+ super(parm, rand, mcp, hist);
+
+ a = parm.getHADResolution();
+ b = parm.getHADConstantTerm();
+ c = parm.getHADPositionError();
+ d = parm.getHADAlignmentError();
+
+ smear(rand, hist);
+ }
+
+ protected void smearPosition(Random rand, double E, boolean hist) {
+ double transSigma = c / Math.sqrt(E) + d;
+ transDist = transSigma * rand.nextGaussian();
+ if (hist)
+ AIDA.defaultInstance().cloud1D("HAD: transDist").fill(transDist);
+ smearPosition(rand);
+ }
+}
Added: projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/reconstructedparticle/IDResolutionTables.java
=============================================================================
--- projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/reconstructedparticle/IDResolutionTables.java (added)
+++ projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/reconstructedparticle/IDResolutionTables.java Fri Jan 9 14:44:03 2015
@@ -0,0 +1,59 @@
+/*
+ * IDResolutionTables.java
+ *
+ * Created on July 6, 2005, 7:13 PM
+ *
+ * @version "$Id: IDResolutionTables.java,v 1.2 2006/09/08 03:19:45 ngraf Exp $"
+ */
+
+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;
+ }
+}
Added: projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/reconstructedparticle/MCFastParticleID.java
=============================================================================
--- projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/reconstructedparticle/MCFastParticleID.java (added)
+++ projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/reconstructedparticle/MCFastParticleID.java Fri Jan 9 14:44:03 2015
@@ -0,0 +1,52 @@
+package org.lcsim.mc.fast.reconstructedparticle;
+
+import hep.physics.particle.properties.ParticleType;
+import org.lcsim.event.ParticleID;
+
+/**
+ * An implementation of ParticleID appropriate for the fast Monte Carlo
+ * @author ngraf
+ */
+public class MCFastParticleID implements ParticleID {
+ private ParticleType _type;
+
+ /** Creates a new instance of MCFastParticleID */
+ public MCFastParticleID(ParticleType type) {
+ _type = type;
+ }
+
+ /**
+ * Type - defined to be the pdgId.
+ */
+ public int getType() {
+ return _type.getPDGID();
+ }
+
+ /**
+ * The PDG code of this id - UnknownPDG ( 999999 ) if unknown.
+ */
+ public int getPDG() {
+ return _type.getPDGID();
+ }
+
+ /**
+ * The likelihood of this hypothesis - in a user defined normalization.
+ */
+ public double getLikelihood() {
+ return 1.;
+ }
+
+ /**
+ * Type of the algorithm/module that created this hypothesis. Check/set collection parameters PIDAlgorithmTypeName and PIDAlgorithmTypeID.
+ */
+ public int getAlgorithmType() {
+ return 0;
+ }
+
+ /**
+ * Parameters associated with this hypothesis. Check/set collection paramter PIDParameterNames for decoding the indices.
+ */
+ public double[] getParameters() {
+ return new double[1];
+ }
+}
Added: projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/reconstructedparticle/MCFastReconstructedParticle.java
=============================================================================
--- projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/reconstructedparticle/MCFastReconstructedParticle.java (added)
+++ projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/reconstructedparticle/MCFastReconstructedParticle.java Fri Jan 9 14:44:03 2015
@@ -0,0 +1,213 @@
+package org.lcsim.mc.fast.reconstructedparticle;
+
+import org.lcsim.event.Vertex;
+import org.lcsim.mc.fast.MCFast;
+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();
+
+ public MCFastReconstructedParticle(Track t, ParticleType type, Particle p, Cluster assoc_c, double wtcal, boolean refPoint000) {
+ MCFast.log.info(" PDGID= " + type.getPDGID() + " t.getPX,...= " + t.getPX() + " " + t.getPY() + " " + t.getPZ());
+ _mass = type.getMass();
+ addTrack(t);
+ _charge = t.getCharge();
+ // Use (0,0,0) for reference point if refPoint000=true
+ // Use true point of origin for reference point if refPoint000=false
+ if (refPoint000) {
+ _referencePoint = new BasicHep3Vector(0, 0, 0);
+ } else {
+ _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 (assoc_c != null) {
+ addCluster(assoc_c);
+ MCFast.log.info(" PDGID= " + type.getPDGID() + " e_track= " + e_track + " e_assoc_clus= " + assoc_c.getEnergy());
+ MCFast.log.info(" PDGID= " + type.getPDGID() + " _referencePoint= " + _referencePoint.x() + " " + _referencePoint.y() + " " + _referencePoint.z());
+ MCFast.log.info(" PDGID= " + type.getPDGID() + " p3_track= " + p3_track.x() + " " + p3_track.y() + " " + p3_track.z());
+ } else {
+ MCFast.log.info(" 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();
+ if (e < _mass) {
+ if (e > _mass - 1.e-5 && e > Double.MIN_VALUE) {
+ _mass = e - Double.MIN_VALUE;
+ } else {
+ MCFast.log.warning(" MCFastReconstructedParticle e < _mass e= " + e + " _mass= " + _mass);
+ MCFast.log.warning(" MCFastReconstructedParticle type = " + type);
+ MCFast.log.warning(" MCFastReconstructedParticle program will continue, but problem should be fixed ");
+ _mass = e - Double.MIN_VALUE;
+ // System.exit(0);
+ }
+
+ }
+
+ 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) {
+ _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);
+ MCFast.log.info(" PDGID= " + type.getPDGID() + " e_reco= " + e_reco + " mass= " + _mass);
+ MCFast.log.info(" PDGID= " + type.getPDGID() + " _referencePoint= " + _referencePoint.x() + " " + _referencePoint.y() + " " + _referencePoint.z());
+ MCFast.log.info(" 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();
+ }
+
+ public Vertex getStartVertex() {
+ return null;
+ }
+}
Added: projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/reconstructedparticle/MCFastReconstructedParticleDriver.java
=============================================================================
--- projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/reconstructedparticle/MCFastReconstructedParticleDriver.java (added)
+++ projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/reconstructedparticle/MCFastReconstructedParticleDriver.java Fri Jan 9 14:44:03 2015
@@ -0,0 +1,205 @@
+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 boolean refPoint000;
+ 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(true);
+ }
+
+ public MCFastReconstructedParticleDriver(boolean refPoint000) {
+ this(ParticlePropertyManager.getParticlePropertyProvider(), refPoint000);
+ }
+
+ public MCFastReconstructedParticleDriver(ParticlePropertyProvider ppp) {
+ this(ppp, true);
+ }
+
+ public MCFastReconstructedParticleDriver(ParticlePropertyProvider ppp, boolean refPoint000) {
+ this.refPoint000 = refPoint000;
+ 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(), refPoint000);
+ 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 (m_ct.get(c) != null)
+ 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();
+ int pdgid = p.getPDGID();
+ if ((abs(pdgid) == 2112) && (rand.nextDouble() < IDEff.getNeutronEff())) {
+ type = neutron;
+ } else {
+ type = klong;
+ }
+ if (m_ct.get(c) != null || c.getEnergy() < type.getMass())
+ continue;
+ if (hist) {
+ aida.histogram1D("hadronCLS-particle", 150, -3, 3).fill((emc.getEnergy() - emc.getMCParticle().getEnergy()) / (Math.sqrt(emc.getMCParticle().getEnergy())));
+ }
+
+ }
+ 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);
+ }
+
+}
Added: projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/DocaTrackParameters.java
=============================================================================
--- projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/DocaTrackParameters.java (added)
+++ projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/DocaTrackParameters.java Fri Jan 9 14:44:03 2015
@@ -0,0 +1,668 @@
+package org.lcsim.mc.fast.tracking;
+
+import Jama.util.Maths;
+import hep.physics.matrix.SymmetricMatrix;
+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;
+import org.lcsim.event.MCParticle;
+
+/**
+ * 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.9 2007/10/16 18:16:50 cassell Exp $
+ */
+public class DocaTrackParameters implements TrackParameters {
+ private Hep3Vector m_pdoca_ref = null;
+ private Hep3Vector m_xdoca_ref = null;
+ private Hep3Vector m_xref = null;
+ private SymmetricMatrix m_err = new SymmetricMatrix(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
+ //
+ // ====================================================
+ public DocaTrackParameters(double bField) {
+ reset();
+ m_Bz = bField;
+ }
+
+ public DocaTrackParameters(MCParticle p, double bField) {
+ this(bField, p.getMomentum(), p.getOrigin(), p.getCharge());
+ }
+
+ public DocaTrackParameters(double bField, Particle p) {
+ this(bField, p.getMomentum(), p.getOrigin(), p.getType().getCharge());
+ }
+
+ public DocaTrackParameters(double bField, Hep3Vector momentum, Hep3Vector x, double q) {
+ reset();
+ m_Bz = bField;
+ calculateDoca(momentum, x, q);
+ }
+
+ public DocaTrackParameters(double bField, double[] momentum, double[] x, double q) {
+ reset();
+ m_Bz = bField;
+ calculateDoca(momentum, x, q);
+ }
+
+ public DocaTrackParameters(double bField, double[] momentum, double[] x, double q, SymmetricMatrix errorMatrix) {
+ this(bField, momentum, x, q);
+ this.m_err = errorMatrix;
+ }
+
+ public DocaTrackParameters(double bField, double[] parameters) {
+ m_Bz = bField;
+ m_parm = parameters;
+ }
+
+ public DocaTrackParameters(double bField, double[] parameters, SymmetricMatrix errorMatrix) {
+ this(bField, parameters);
+ this.m_err = errorMatrix;
+ }
+
+ public DocaTrackParameters(double bField, double[] parameters, SymmetricMatrix errorMatrix, double chi2) {
+ this(bField, parameters);
+ this.m_err = errorMatrix;
+ setChi2(chi2);
+ }
+
+ public DocaTrackParameters(double bField, double[] parameters, SymmetricMatrix errorMatrix, double chi2, int ndf) {
+ this(bField, parameters);
+ this.m_err = errorMatrix;
+ setChi2(chi2);
+ setNDF(ndf);
+ }
+
+ public double getD0() {
+ return m_parm[0];
+ }
+
+ /**
+ * Get the error matrix as a 2-D array
+ * @see #getTrackParameter
+ */
+ public SymmetricMatrix getErrorMatrix() {
+ return m_err;
+ }
+
+ /**
+ * 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 + l * 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.
+ */
+ SymmetricMatrix 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 = Maths.toJamaMatrix(getErrorMatrix());
+ Matrix pErrorMatrix = tMatrix.transpose().times(errorMatrix.times(tMatrix));
+
+ return new SymmetricMatrix(Maths.fromJamaMatrix(pErrorMatrix));
+ }
+
+ /**
+ * 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 = Maths.toJamaMatrix(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();
+ }
+
+ // ====================================================
+ //
+ // 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 <= i; j++) {
+ m_err.setElement(i, j, 0);
+ }
+ }
+ m_l0 = 0.;
+ }
+
+ void setErrorMatrix(SymmetricMatrix error) {
+ m_err = error;
+ }
+}
Added: projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/LookupTable.java
=============================================================================
--- projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/LookupTable.java (added)
+++ projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/LookupTable.java Fri Jan 9 14:44:03 2015
@@ -0,0 +1,118 @@
+package org.lcsim.mc.fast.tracking;
+
+import java.io.BufferedReader;
+import java.io.IOException;
+
+import java.util.Arrays;
+import java.util.StringTokenizer;
+
+public 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");
+ }
+ }
+
+ public double interpolateVal(double val1, double val2) {
+ int index1 = binarySearch(m_key1, val1);
+ // implement cut-off
+ double t;
+ if (index1 < 0) {
+ t = m_key1[0];
+ index1 = 0;
+ } else if (index1 >= m_key1.length - 1) {
+ t = m_key1[m_key1.length - 1];
+ index1 = m_key1.length - 1;
+ } else
+ t = (val1 - m_key1[index1]) / (m_key1[index1 + 1] - m_key1[index1]);
+
+ double u;
+ int index2 = binarySearch(m_key2, val2);
+ if (index2 < 0) {
+ u = m_key2[0];
+ index2 = 0;
+ } else if (index2 >= m_key2.length - 1) {
+ u = m_key2[m_key2.length - 1];
+ index2 = m_key2.length - 1;
+ } else
+ 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)
+ {
+ if (value < key[0]) {
+ // throw new RuntimeException("Interpolation out of range: lower: "+value+" < "+key[0]);
+ return 0;
+ }
+
+ int pos = Arrays.binarySearch(key, value);
+ if (pos > 0) {
+ return Math.min(pos, key.length - 2);
+ } else {
+ 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]);
+ }
+}
Added: projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/MCFastTracking.java
=============================================================================
--- projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/MCFastTracking.java (added)
+++ projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/MCFastTracking.java Fri Jan 9 14:44:03 2015
@@ -0,0 +1,158 @@
+package org.lcsim.mc.fast.tracking;
+
+/**
+ * $Id: MCFastTracking.java,v 1.19 2009/05/29 22:49:55 timb Exp $
+ */
+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.event.LCRelation;
+import org.lcsim.event.MCParticle;
+import org.lcsim.util.Driver;
+import org.lcsim.event.base.MyLCRelation;
+
+/**
+ * Fast Monte Carlo tracking simulator
+ */
+public class MCFastTracking extends Driver implements ConditionsListener {
+ private TrackResolutionTables parm;
+ private SimpleTables SmTbl;
+ private boolean beamSpotConstraint;
+ private boolean simple;
+ private final static double[] IP = { 0, 0, 0 };
+ private boolean defaultMC = true;
+ private String fsname;
+ private String outputListName = null;
+
+ 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);
+ }
+ }
+
+ public void setFSList(String fslist) {
+ fsname = fslist;
+ defaultMC = false;
+ }
+
+ protected void process(EventHeader event) {
+ if (defaultMC) {
+ fsname = "MCParticle";
+ } else {
+ if (!event.hasCollection(MCParticle.class, fsname)) {
+ System.err.println("Collection " + fsname + " not found. Default Final State particles being used");
+ fsname = "MCParticle";
+ defaultMC = true;
+ }
+ }
+ if (parm == null) {
+ ConditionsSet conditions = getConditionsManager().getConditions("TrackParameters");
+ conditions.addConditionsListener(this);
+ parm = setTrackResolutionTables(conditions, beamSpotConstraint);
+ }
+
+ if (SmTbl == null) {
+ ConditionsSet conditions = getConditionsManager().getConditions("SimpleTrack");
+ conditions.addConditionsListener(this);
+ SmTbl = new SimpleTables(conditions);
+ }
+
+ double bField = event.getDetector().getFieldMap().getField(IP)[2];
+ boolean hist = getHistogramLevel() > 0;
+
+ List<Track> trackList = new ArrayList<Track>();
+ List<LCRelation> lcrelationList = new ArrayList<LCRelation>();
+ for (Iterator i = event.get(MCParticle.class, fsname).iterator(); i.hasNext();) {
+ Particle p = (Particle) i.next();
+
+ // filter for FINAL_STATE
+ if (defaultMC) {
+ if (p.getGeneratorStatus() != Particle.FINAL_STATE) {
+ continue;
+ }
+ }
+ double pCharge = p.getCharge();
+ if (pCharge == 0 || Double.isNaN(pCharge) || pCharge == Double.NEGATIVE_INFINITY || pCharge == Double.POSITIVE_INFINITY) {
+ continue;
+ }
+
+ double[] momentum = p.getMomentum().v();
+ if (Double.isNaN(momentum[0]) || Double.isNaN(momentum[1]) || Double.isNaN(momentum[2])) {
+ continue;
+ }
+
+ 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;
+ }
+
+ try {
+ ReconTrack newTrack = new ReconTrack(bField, parm, SmTbl, getRandom(), p, hist, simple);
+ trackList.add(newTrack);
+ lcrelationList.add(new MyLCRelation((Track) newTrack, (MCParticle) p));
+ } catch (hep.physics.particle.properties.UnknownParticleIDException x) {
+ System.out.println("WARNING: MCFastTracking ignored a particle of type " + p.getPDGID());
+ }
+ }
+ if (outputListName == null) {
+ outputListName = EventHeader.TRACKS;
+ }
+ event.put(outputListName, trackList, Track.class, 0);
+ event.put("TracksToMCP", lcrelationList, LCRelation.class, 0);
+ }
+
+ /** Specify the name under which to write out the list of tracks to the event. Default is EventHeader.TRACKS ("Tracks") */
+ public void setOutputList(String name) {
+ outputListName = name;
+ }
+
+ public void conditionsChanged(ConditionsEvent event) {
+ ConditionsSet conditions = getConditionsManager().getConditions("TrackParameters");
+ ConditionsSet simpleconditions = getConditionsManager().getConditions("SimpleTrack");
+ parm = setTrackResolutionTables(conditions, beamSpotConstraint);
+ SmTbl = new SimpleTables(simpleconditions);
+ }
+}
Added: projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/ReconTrack.java
=============================================================================
--- projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/ReconTrack.java (added)
+++ projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/ReconTrack.java Fri Jan 9 14:44:03 2015
@@ -0,0 +1,504 @@
+package org.lcsim.mc.fast.tracking;
+
+import Jama.*;
+import hep.physics.matrix.SymmetricMatrix;
+import org.lcsim.util.aida.AIDA;
+import hep.physics.particle.Particle;
+import hep.physics.vec.Hep3Vector;
+import org.lcsim.event.Track;
+import org.lcsim.event.TrackState;
+import org.lcsim.event.base.BaseTrackState;
+
+import java.io.*;
+import java.util.ArrayList;
+import java.util.Collections;
+import java.util.List;
+import java.util.Random;
+import org.lcsim.mc.fast.tracking.SimpleTables;
+
+/**
+ * Provides MC smeared track. <br>
+ * Handles to smeared and not-smeared set of DocaTrackParameters are provided. <br>
+ *
+ * @author Tony Johnson, Wolfgang Walkowiak
+ * @version $Id: ReconTrack.java,v 1.12 2012/07/26 16:46:14 grefe Exp $
+ */
+public class ReconTrack implements Track {
+ private static final String[][] matrixI = { { "(1,1):", "(1,2):", "(1,3):", "(1,4):", "(1,5):" }, { "(2,1):", "(2,2):", "(2,3):", "(2,4):", "(2,5):" }, { "(3,1):", "(3,2):", "(3,3):", "(3,4):", "(3,5):" }, { "(4,1):", "(4,2):", "(4,3):", "(4,4):", "(4,5):" }, { "(5,1):", "(5,2):", "(5,3):", "(5,4):", "(5,5):" } };
+
+ // =======================================================================
+ //
+ // private members
+ //
+ // =======================================================================
+ private DocaTrackParameters m_nosmear = null;
+ private DocaTrackParameters m_smear = null;
+ transient private Particle mc;
+ private int m_tcharge;
+ private double[] _refpoint = { 0, 0, 0 };
+ private List<TrackState> _trackStates;
+
+ ReconTrack(double bField, TrackResolutionTables parm, SimpleTables SmTbl, Random rand, Particle mc, boolean hist, boolean simple) {
+ this.mc = mc;
+
+ // get original momentum from MCParticle
+ // convert to helical parameters
+ m_nosmear = new DocaTrackParameters(bField, mc);
+
+ double pt = m_nosmear.getPt();
+
+ if (hist) {
+
+ double r = Math.abs(m_nosmear.getD0());
+
+ AIDA aida = AIDA.defaultInstance();
+ aida.cloud1D("ptsqr").fill(pt * pt);
+ aida.cloud1D("pt").fill(pt);
+ aida.cloud1D("phi").fill(m_nosmear.getPhi0());
+ aida.cloud1D("theta").fill(m_nosmear.getTheta());
+ aida.cloud1D("tanL").fill(m_nosmear.getTanL());
+ aida.cloud1D("r").fill(r);
+ aida.cloud1D("z").fill(mc.getOriginZ());
+ }
+ // get appropriate resolution table
+ double abscth = Math.abs(m_nosmear.getCosTheta());
+ double ptot = m_nosmear.getPtot();
+ ResolutionTable table = (abscth < parm.getPolarInner()) ? parm.getBarrelTable() : parm.getEndcapTable();
+
+ // get resolution values from interpolation and fill error matrix
+ m_nosmear.setErrorMatrix(getErrorMatrixFromTable(table, abscth, ptot));
+
+ // smear tracks according to error matrix
+ if (simple == true) {
+ m_smear = (DocaTrackParameters) SmearTrackSimple.smearTrackSimple(bField, m_nosmear, rand, SmTbl, pt, hist);
+ // double[] slice = {0, 0, 1, 0, 0};
+ // m_smear = (DocaTrackParameters) SmearTrackSimpleII.SmearTrackSimpleII(bField, m_nosmear, rand, SmTbl, slice, hist);
+ } else {
+ m_smear = (DocaTrackParameters) SmearTrack.smearTrack(bField, m_nosmear, rand);
+ }
+
+ if (hist) {
+ AIDA aida = AIDA.defaultInstance();
+ aida.cloud1D("ptNew").fill(m_smear.getPt());
+ aida.cloud1D("tanLNew").fill(m_smear.getTanL());
+ aida.cloud1D("rNew").fill(Math.abs(m_smear.getD0()));
+ aida.cloud1D("phiNew").fill(m_smear.getPhi0());
+ aida.cloud1D("zNew").fill(m_smear.getZ0());
+ }
+ m_tcharge = (int) (m_smear.getUnitCharge() * Math.abs(mc.getType().getCharge()));
+ _trackStates = new ArrayList<TrackState>();
+ _trackStates.add(new BaseTrackState(m_smear.getTrackParameters(), bField));
+ }
+
+ /**
+ * Get the full charge.
+ */
+ public int getCharge() {
+ return m_tcharge;
+ }
+
+ /**
+ * Get the chi2 from smearing.
+ */
+ public double getChi2() {
+ return m_smear.getChi2();
+ }
+
+ /**
+ * Get DOCA (2-dim) of smeared track. <br>
+ *
+ * Note: Use #getNotSmearedTrack().getDOCA() to access parameters of the not smeared track.
+ */
+ public double[] getDoca() {
+ return m_smear.getDoca();
+ }
+
+ public double[] getDocaMomentum(double[] refPoint) {
+ return m_smear.getDocaMomentum(refPoint);
+ }
+
+ /**
+ * Get momentum at DOCA (2-dim) of smeared track. <br>
+ *
+ * Note: Use #getNotSmearedTrack().getMomentum() to access parameters of the not smeared track.
+ */
+ public double[] getDocaMomentum() {
+ return m_smear.getMomentum();
+ }
+
+ public Hep3Vector getDocaMomentumVec(Hep3Vector refPoint) {
+ return m_smear.getDocaMomentumVec(refPoint);
+ }
+
+ /*
+ * Calculate and get Doca momentum on smeared track with respect to any space point.
+ */
+ public Hep3Vector getDocaMomentumVec(double[] refPoint) {
+ return m_smear.getDocaMomentumVec(refPoint);
+ }
+
+ /**
+ * Get x coordinate of momentum of the smeared track at DOCA.
+ */
+ public double getDocaMomentumX() {
+ return m_smear.getPX();
+ }
+
+ /**
+ * Get y coordinate of momentum of the smeared track at DOCA.
+ */
+ public double getDocaMomentumY() {
+ return m_smear.getPY();
+ }
+
+ /**
+ * Get z coordinate of momentum of the smeared track at DOCA.
+ */
+ public double getDocaMomentumZ() {
+ return m_smear.getPZ();
+ }
+
+ public double[] getDocaPosition(double[] refPoint) {
+ return m_smear.getDocaPosition(refPoint);
+ }
+
+ /**
+ * Calculate and get Doca position on the smeared track with respect to any space point.
+ */
+ public Hep3Vector getDocaPositionVec(Hep3Vector refPoint) {
+ return m_smear.getDocaPositionVec(refPoint);
+ }
+
+ public Hep3Vector getDocaPositionVec(double[] refPoint) {
+ return m_smear.getDocaPositionVec(refPoint);
+ }
+
+ /**
+ * Get transverse momentum of the smeared track at DOCA.
+ */
+ public double getDocaPt() {
+ return m_smear.getPt();
+ }
+
+ /**
+ * Calculate and get path length on the smeared 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.
+ */
+ public double getDocaTransversePathLength(Hep3Vector refPoint) {
+ return m_smear.getDocaTransversePathLength(refPoint);
+ }
+
+ public double getDocaTransversePathLength(double[] refPoint) {
+ return m_smear.getDocaTransversePathLength(refPoint);
+ }
+
+ /**
+ * Get x coordinate of DOCA of smeared track.
+ */
+ public double getDocaX() {
+ return m_smear.getDocaX();
+ }
+
+ /**
+ * Get y coordinate of DOCA of smeared track.
+ */
+ public double getDocaY() {
+ return m_smear.getDocaY();
+ }
+
+ /**
+ * Get z coordinate of DOCA of smeared track.
+ */
+ public double getDocaZ() {
+ return m_smear.getDocaZ();
+ }
+
+ /**
+ * Get the full error matrix.
+ * @see #getTrackParameter
+ */
+ public SymmetricMatrix getErrorMatrix() {
+ return m_smear.getErrorMatrix();
+ }
+
+ /**
+ * Get the MC particle for this track.
+ */
+ public Particle getMCParticle() {
+ return mc;
+ }
+
+ public double[] getMomentum(double l) {
+ return m_smear.getMomentum(l);
+ }
+
+ /**
+ * Get momentum of smeared track at original vertex point. <br>
+ *
+ * Note: Use #getNotSmearedTrack().getMomentum() to access parameters of the not smeared track.
+ */
+ public double[] getMomentum() {
+ return m_smear.getMomentum(m_smear.getL0());
+ }
+
+ /**
+ * Calculate and get momentum on track with respect to any path length l on track (l in xy plane).
+ */
+ public Hep3Vector getMomentumVec(double l) {
+ return m_smear.getMomentumVec(l);
+ }
+
+ public Hep3Vector getMomentumVec() {
+ return m_smear.getMomentumVec(m_smear.getL0());
+ }
+
+ /**
+ * Get x coordinate of momentum of the smeared track at original vertex.
+ */
+ public double getMomentumX() {
+ return m_smear.getMomentum(m_smear.getL0())[0];
+ }
+
+ /**
+ * Get y coordinate of momentum of the smeared track at original vertex.
+ */
+ public double getMomentumY() {
+ return m_smear.getMomentum(m_smear.getL0())[1];
+ }
+
+ /**
+ * Get z coordinate of momentum of the smeared track at original vertex.
+ */
+ public double getMomentumZ() {
+ return m_smear.getMomentum(m_smear.getL0())[2];
+ }
+
+ /**
+ * Get the number degrees of freedom.
+ */
+ public int getNDF() {
+ return m_smear.getNDF();
+ }
+
+ /**
+ * Get the complete parameter set for the not smeared track.
+ */
+ public DocaTrackParameters getNotSmearedTrack() {
+ return m_nosmear;
+ }
+
+ /**
+ * Get x coordinate of momentum of the smeared track at original vertex.
+ */
+ public double getPX() {
+ return getMomentumX();
+ }
+
+ /**
+ * Get y coordinate of momentum of the smeared track at original vertex.
+ */
+ public double getPY() {
+ return getMomentumY();
+ }
+
+ /**
+ * Get z coordinate of momentum of the smeared track at original vertex.
+ */
+ public double getPZ() {
+ return getMomentumZ();
+ }
+
+ public double[] getPosition(double l) {
+ return m_smear.getPosition(l);
+ }
+
+ /**
+ * Calculate and get position on track with respect to any path length l on track (l in xy plane).
+ */
+ public Hep3Vector getPositionVec(double l) {
+ return m_smear.getPositionVec(l);
+ }
+
+ /**
+ * Get transverse momentum of the smeared track at original vertex.
+ */
+ public double getPt() {
+ double[] p = getMomentum();
+
+ return Math.sqrt((p[0] * p[0]) + (p[1] * p[1]));
+ }
+
+ public double getRadiusOfInnermostHit() {
+ return 0; // FixMe:
+ }
+
+ /**
+ * Get the original vertex point of smeared MC track.
+ */
+
+ // Get the reference point used for track parameter calculations
+ //
+ public double[] getReferencePoint() {
+ return _refpoint;
+ }
+
+ public boolean isReferencePointPCA() {
+ return true;
+ }
+
+ /**
+ * Get x coordinate of the original vertex point of smeared MC track.
+ */
+ public double getReferencePointX() {
+ return getReferencePoint()[0];
+ }
+
+ /**
+ * Get y coordinate of the original vertex point of smeared MC track.
+ */
+ public double getReferencePointY() {
+ return getReferencePoint()[1];
+ }
+
+ /**
+ * Get z coordinate of the original vertex point of smeared MC track.
+ */
+ public double getReferencePointZ() {
+ return getReferencePoint()[2];
+ }
+
+ /**
+ * Get the complete parameter set for the smeared track.
+ */
+ public DocaTrackParameters getSmearedTrack() {
+ return m_smear;
+ }
+
+ public int[] getSubdetectorHitNumbers() {
+ return new int[0]; // FIXME
+ }
+
+ /**
+ * 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_smear.getTrackParameter(i);
+ }
+
+ /**
+ * Get the track parameters as an array
+ * @see #getTrackParameter
+ */
+ public double[] getTrackParameters() {
+ return m_smear.getTrackParameters();
+ }
+
+ public List getTrackerHits() {
+ return Collections.EMPTY_LIST;
+ }
+
+ public List getTracks() {
+ return Collections.EMPTY_LIST;
+ }
+
+ public int getType() {
+ return 0; // FIXME:
+ }
+
+ /**
+ * Calculate the error matrix for the momentum for a point on the smeared track specified by l.
+ */
+ public SymmetricMatrix calcMomentumErrorMatrix(double l) {
+ return m_smear.calcMomentumErrorMatrix(l);
+ }
+
+ /**
+ * Calculate the error matrix for the position coordinates for a point on the smeared track specified by l. Result is given as a 3x3 array for the matrix.
+ */
+ public double[][] calcPositionErrorMatrix(double l) {
+ return m_smear.calcPositionErrorMatrix(l);
+ }
+
+ public boolean fitSuccess() {
+ // FIXME Not implemented in fastMC
+ return false;
+ }
+
+ public double getdEdx() {
+ return 0; // FIXME
+ }
+
+ public double getdEdxError() {
+ return 0; // FIXME
+ }
+
+ public String toString() {
+ java.io.StringWriter buffer = new java.io.StringWriter();
+ java.io.PrintWriter out = new java.io.PrintWriter(buffer);
+
+ java.text.NumberFormat pf = java.text.NumberFormat.getInstance();
+ pf.setMinimumFractionDigits(12);
+ pf.setMaximumFractionDigits(12);
+ pf.setMinimumIntegerDigits(1);
+ pf.setMaximumIntegerDigits(3);
+ pf.setGroupingUsed(false);
+
+ Matrix mHlxPar = new Matrix(m_smear.getTrackParameters(), 1);
+
+ out.println("ReconTrack Parameters: d0 phi0 omega z0 tan(lambda)");
+ mHlxPar.print(out, pf, 16);
+ out.println("Error Matrix:");
+ out.println(m_smear.getErrorMatrix());
+
+ return (buffer.toString());
+ }
+
+ private SymmetricMatrix getErrorMatrixFromTable(ResolutionTable table, double abscth, double ptot) {
+ SymmetricMatrix errMatrix = new SymmetricMatrix(5);
+ for (int i = 0; i < 5; i++) {
+ for (int j = 0; j <= i; j++) {
+ errMatrix.setElement(i, j, table.findTable(matrixI[i][j]).interpolateVal(abscth, ptot));
+ }
+ }
+ return errMatrix;
+ }
+
+ public List<TrackState> getTrackStates() {
+ return _trackStates;
+ }
+}
Added: projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/ResolutionTable.java
=============================================================================
--- projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/ResolutionTable.java (added)
+++ projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/ResolutionTable.java Fri Jan 9 14:44:03 2015
@@ -0,0 +1,91 @@
+package org.lcsim.mc.fast.tracking;
+
+import java.io.BufferedReader;
+import java.io.IOException;
+import java.io.Reader;
+import java.util.HashMap;
+import java.util.Map;
+import java.util.StringTokenizer;
+
+public class ResolutionTable {
+ private Map hash = new HashMap();
+
+ public ResolutionTable(Reader in) throws IOException {
+ BufferedReader inn = new TrimReader(in);
+ setupTable(inn);
+ inn.close();
+ }
+
+ LookupTable findTable(String name) {
+ return (LookupTable) hash.get(name);
+ }
+
+ public LookupTable findTable(int i, int j) {
+ String name = String.format("(%s,%s):", i + 1, j + 1);
+ // System.out.printf("%s\n", name);
+ return findTable(name);
+ }
+
+ void setupTable(BufferedReader in) throws IOException {
+ // constant string to read Bruce's resolution table ini files
+ String tokenLine = "Cov matrix entry";
+
+ // Start by reading the parameter file
+ for (;;) {
+ String line = in.readLine();
+ if (line == null)
+ break;
+
+ //
+ // ww, 10/05/2000
+ //
+ // transform Bruce's cov matrix element index format into LCDs
+ // ie "Cov matrix entry n, n"
+ // to "(n,n):"
+ //
+ int elementI = 0;
+ int elementJ = 0;
+ if (line.indexOf(tokenLine) > -1) {
+ String nLine = line.substring(line.indexOf(tokenLine) + tokenLine.length());
+ StringBuffer sb = new StringBuffer("(");
+ StringTokenizer tokenizer = new StringTokenizer(nLine, " ,");
+ if (tokenizer.countTokens() == 2) {
+ String token = tokenizer.nextToken();
+ elementI = Integer.parseInt(token) - 1; // lcdtrk starts from 1
+ sb.append(token);
+ sb.append(",");
+ token = tokenizer.nextToken();
+ elementJ = Integer.parseInt(token) - 1; // lcdtrk starts from 1
+ sb.append(token);
+ sb.append("):");
+ } else {
+ sb.append("--none--)");
+ }
+ line = sb.toString();
+ }
+
+ // end of addition
+ if (!line.endsWith(":")) {
+ throw new IOException("Syntax error in ResolutionTable");
+ }
+ hash.put(line, new LookupTable(in, elementI, elementJ));
+ }
+ }
+
+ private class TrimReader extends BufferedReader {
+ TrimReader(Reader source) {
+ super(source);
+ }
+
+ public String readLine() throws IOException {
+ for (;;) {
+ String line = super.readLine();
+ if (line == null)
+ return null;
+ line = line.trim();
+ if (line.length() > 0)
+ return line;
+ }
+ }
+ }
+}
Added: projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/SimpleTables.java
=============================================================================
--- projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/SimpleTables.java (added)
+++ projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/SimpleTables.java Fri Jan 9 14:44:03 2015
@@ -0,0 +1,52 @@
+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 */
+ public 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;
+ }
+
+}
Added: projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/SmearTrack.java
=============================================================================
--- projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/SmearTrack.java (added)
+++ projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/SmearTrack.java Fri Jan 9 14:44:03 2015
@@ -0,0 +1,75 @@
+package org.lcsim.mc.fast.tracking;
+
+import Jama.EigenvalueDecomposition;
+import Jama.Matrix;
+import Jama.util.Maths;
+
+import java.util.Random;
+
+/**
+ * @author W.Walkowiak, 08/00
+ */
+class SmearTrack {
+ /**
+ * Smear track parameters according to the track's stored error matrix.
+ *
+ * @see TrackParameters
+ */
+ static DocaTrackParameters smearTrack(double bField, TrackParameters noSmear, Random rand) {
+
+ // get error matrix and do a sanity check
+ Matrix M = Maths.toJamaMatrix(noSmear.getErrorMatrix());
+ 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, noSmear.getErrorMatrix(), chi2);
+ if (noSmear instanceof DocaTrackParameters) {
+ smeared.setL0(((DocaTrackParameters) noSmear).getL0());
+ }
+
+ return smeared;
+ }
+}
Added: projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/SmearTrackSimple.java
=============================================================================
--- projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/SmearTrackSimple.java (added)
+++ projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/SmearTrackSimple.java Fri Jan 9 14:44:03 2015
@@ -0,0 +1,101 @@
+package org.lcsim.mc.fast.tracking;
+
+import Jama.EigenvalueDecomposition;
+import Jama.Matrix;
+import Jama.util.Maths;
+
+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 = Maths.toJamaMatrix(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().diagonal(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] = Math.abs(noSmear.getOmega()) * 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().e(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;
+ }
+}
Added: projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/TrackParameters.java
=============================================================================
--- projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/TrackParameters.java (added)
+++ projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/TrackParameters.java Fri Jan 9 14:44:03 2015
@@ -0,0 +1,77 @@
+package org.lcsim.mc.fast.tracking;
+
+import hep.physics.matrix.SymmetricMatrix;
+import hep.physics.vec.Hep3Vector;
+
+interface TrackParameters extends Hep3Vector {
+ double getD0();
+
+ /**
+ * Get the error matrix
+ * @see #getTrackParameter
+ */
+ SymmetricMatrix getErrorMatrix();
+
+ double[] getMomentum();
+
+ double getOmega();
+
+ double getPX();
+
+ double getPY();
+
+ double getPZ();
+
+ double getPhi0();
+
+ double getPtot();
+
+ double getTanL();
+
+ /**
+ * 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 q/abs(q) > 0)</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
+ */
+ double getTrackParameter(int i);
+
+ /**
+ * Get the track parameters as an array
+ * @see #getTrackParameter
+ */
+ double[] getTrackParameters();
+
+ int getUnitCharge();
+
+ double getZ0();
+}
Added: projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/TrackResolutionTables.java
=============================================================================
--- projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/TrackResolutionTables.java (added)
+++ projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/TrackResolutionTables.java Fri Jan 9 14:44:03 2015
@@ -0,0 +1,46 @@
+package org.lcsim.mc.fast.tracking;
+
+import org.lcsim.conditions.ConditionsSet;
+
+import java.io.*;
+
+public class TrackResolutionTables {
+ private ResolutionTable barrel;
+ private ResolutionTable endcap;
+ private double PolarInner;
+ private double PolarOuter;
+ private double PtMin;
+
+ public TrackResolutionTables(ConditionsSet set, boolean beamSpotConstraint) throws IOException {
+ PtMin = set.getDouble("PtMin");
+ PolarInner = set.getDouble("PolarInner");
+ PolarOuter = set.getDouble("PolarOuter");
+
+ String postfix = beamSpotConstraint ? "-bc" : "-nbc";
+ String barrelFile = set.getString("BarrelTableFile" + postfix);
+ String endcapFile = set.getString("EndcapTableFile" + postfix);
+
+ barrel = new ResolutionTable(set.getRawSubConditions(barrelFile).getReader());
+ endcap = new ResolutionTable(set.getRawSubConditions(endcapFile).getReader());
+ }
+
+ public ResolutionTable getBarrelTable() {
+ return barrel;
+ }
+
+ public ResolutionTable getEndcapTable() {
+ return endcap;
+ }
+
+ public double getPolarInner() {
+ return PolarInner;
+ }
+
+ public double getPolarOuter() {
+ return PolarOuter;
+ }
+
+ public double getPtMin() {
+ return PtMin;
+ }
+}
Added: projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/fix/FastMCTrack.java
=============================================================================
--- projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/fix/FastMCTrack.java (added)
+++ projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/fix/FastMCTrack.java Fri Jan 9 14:44:03 2015
@@ -0,0 +1,215 @@
+/**
+ * @version $Id: FastMCTrack.java,v 1.3 2012/07/26 16:46:14 grefe Exp $
+ */
+package org.lcsim.mc.fast.tracking.fix;
+
+import static java.lang.Math.sqrt;
+import hep.physics.matrix.SymmetricMatrix;
+import hep.physics.vec.Hep3Vector;
+
+import java.util.ArrayList;
+import java.util.List;
+
+import org.lcsim.event.LCIOParameters;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.Track;
+import org.lcsim.event.TrackState;
+import org.lcsim.event.TrackerHit;
+import org.lcsim.event.LCIOParameters.ParameterName;
+import org.lcsim.event.base.BaseTrackState;
+import org.lcsim.spacegeom.CartesianPoint;
+import org.lcsim.spacegeom.CartesianVector;
+import org.lcsim.spacegeom.SpacePoint;
+import org.lcsim.spacegeom.SpaceVector;
+
+/**
+ * @author jstrube The class to store the measurement information of the track of a charged particle in a magnetic field. This class represents the FastMC simulation. An invariant of this class is
+ * that the "origin" is the point of closest approach to the reference point.
+ */
+public class FastMCTrack implements Track {
+ protected LCIOParameters _parameters;
+ protected SymmetricMatrix _errorMatrix;
+ protected SpacePoint _referencePoint;
+ protected int _charge;
+ protected MCParticle _particle = null;
+ protected List<TrackState> _trackStates;
+
+ protected FastMCTrack(SpacePoint refPoint, LCIOParameters parameters, SymmetricMatrix errorMatrix, int charge) {
+ _referencePoint = refPoint;
+ _parameters = parameters;
+ _charge = charge;
+ _errorMatrix = errorMatrix;
+ _trackStates = new ArrayList<TrackState>();
+ _trackStates.add(new BaseTrackState(parameters.getValues(), errorMatrix.asPackedArray(true), refPoint.v(), 0));
+ }
+
+ protected FastMCTrack(SpacePoint refPoint, LCIOParameters parameters, SymmetricMatrix errorMatrix, int charge, MCParticle part) {
+ this(refPoint, parameters, errorMatrix, charge);
+ _particle = part;
+ }
+
+ public FastMCTrack(Track t) {
+ double[] p = t.getMomentum();
+ double pt = sqrt(p[0] * p[0] + p[1] * p[1]);
+ _parameters = new LCIOParameters(t.getTrackParameters(), pt);
+ _errorMatrix = t.getErrorMatrix();
+ _referencePoint = new CartesianPoint(t.getReferencePoint());
+ _charge = t.getCharge();
+ _trackStates = new ArrayList<TrackState>();
+ _trackStates.add(new BaseTrackState(_parameters.getValues(), _errorMatrix.asPackedArray(true), _referencePoint.v(), 0));
+ }
+
+ public boolean fitSuccess() {
+ return false;
+ }
+
+ /*
+ * (non-Javadoc)
+ *
+ * @see org.lcsim.contrib.JanStrube.tracking.Track#getCharge()
+ */
+ public int getCharge() {
+ return _charge;
+ }
+
+ public double getChi2() {
+ return -1;
+ }
+
+ public double getdEdx() {
+ return 0;
+ }
+
+ public double getdEdxError() {
+ return 0;
+ }
+
+ /*
+ * (non-Javadoc)
+ *
+ * @see org.lcsim.contrib.JanStrube.tracking.Track#getErrorMatrix()
+ */
+ public SymmetricMatrix getErrorMatrix() {
+ return _errorMatrix;
+ }
+
+ public MCParticle getMCParticle() {
+ return _particle;
+ }
+
+ public double[] getMomentum() {
+ return momentum().v();
+ }
+
+ public int getNDF() {
+ return 0;
+ }
+
+ /*
+ * (non-Javadoc)
+ *
+ * @see org.lcsim.contrib.JanStrube.tracking.Track#getParameter(org.lcsim.contrib.JanStrube.tracking.FastMCTrack.ParameterName)
+ */
+ public double getParameter(ParameterName name) {
+ return _parameters.get(name);
+ }
+
+ /*
+ * (non-Javadoc)
+ *
+ * @see org.lcsim.contrib.JanStrube.tracking.Track#getParameters()
+ */
+ public LCIOParameters getParameters() {
+ return _parameters;
+ }
+
+ /*
+ * (non-Javadoc)
+ *
+ * @see org.lcsim.contrib.JanStrube.tracking.Track#getPt()
+ */
+ public double getPt() {
+ return _parameters.getPt();
+ }
+
+ public double getPX() {
+ return momentum().x();
+ }
+
+ public double getPY() {
+ return momentum().y();
+ }
+
+ public double getPZ() {
+ return momentum().z();
+ }
+
+ public double getRadiusOfInnermostHit() {
+ return -1;
+ }
+
+ public double[] getReferencePoint() {
+ return _referencePoint.v();
+ }
+
+ public double getReferencePointX() {
+ return referencePoint().x();
+ }
+
+ public double getReferencePointY() {
+ return referencePoint().y();
+ }
+
+ public double getReferencePointZ() {
+ return referencePoint().z();
+ }
+
+ public int[] getSubdetectorHitNumbers() {
+ return null;
+ }
+
+ public List<TrackerHit> getTrackerHits() {
+ return null;
+ }
+
+ public double getTrackParameter(int i) {
+ return _parameters.getValues()[i];
+ }
+
+ public double[] getTrackParameters() {
+ return _parameters.getValues();
+ }
+
+ public List<Track> getTracks() {
+ return null;
+ }
+
+ public int getType() {
+ return 0;
+ }
+
+ public boolean isReferencePointPCA() {
+ return true;
+ }
+
+ public Hep3Vector momentum() {
+ return new CartesianVector(LCIOParameters.Parameters2Momentum(_parameters).v());
+ }
+
+ public SpacePoint position() {
+ return LCIOParameters.Parameters2Position(_parameters, _referencePoint);
+ }
+
+ /*
+ * (non-Javadoc)
+ *
+ * @see org.lcsim.contrib.JanStrube.tracking.Track#getReferencePoint()
+ */
+ public SpacePoint referencePoint() {
+ return _referencePoint;
+ }
+
+ public List<TrackState> getTrackStates() {
+ return _trackStates;
+ }
+}
Added: projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/fix/FastMCTrackDriver.java
=============================================================================
--- projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/fix/FastMCTrackDriver.java (added)
+++ projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/fix/FastMCTrackDriver.java Fri Jan 9 14:44:03 2015
@@ -0,0 +1,106 @@
+/**
+ * @version $Id: FastMCTrackDriver.java,v 1.1 2007/11/29 21:29:35 jstrube Exp $
+ */
+package org.lcsim.mc.fast.tracking.fix;
+
+import hep.physics.particle.Particle;
+
+import java.io.IOException;
+import java.util.ArrayList;
+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.MCParticle;
+import org.lcsim.event.Track;
+import org.lcsim.mc.fast.tracking.SimpleTables;
+import org.lcsim.mc.fast.tracking.TrackResolutionTables;
+import org.lcsim.util.Driver;
+
+/**
+ * A replacement for the current FastMC Tracks. The simple tables are currently not implemented
+ * @author jstrube
+ *
+ */
+public class FastMCTrackDriver extends Driver implements ConditionsListener {
+ private TrackResolutionTables parm;
+ private SimpleTables SmTbl;
+ private boolean beamSpotConstraint;
+ private boolean simple;
+ private final static double[] IP = { 0, 0, 0 };
+
+ public FastMCTrackDriver(boolean beamSpotConstraint) {
+ this.beamSpotConstraint = beamSpotConstraint;
+ }
+
+ public FastMCTrackDriver() {
+ this(false);
+ }
+
+ protected void process(EventHeader event) {
+ if (parm == null) {
+ ConditionsSet conditions = getConditionsManager().getConditions("TrackParameters");
+ conditions.addConditionsListener(this);
+ parm = setTrackResolutionTables(conditions, beamSpotConstraint);
+ }
+
+ if (SmTbl == null) {
+ ConditionsSet conditions = getConditionsManager().getConditions("SimpleTrack");
+ conditions.addConditionsListener(this);
+ SmTbl = new SimpleTables(conditions);
+ }
+
+ FastMCTrackFactory factory = new FastMCTrackFactory(event, beamSpotConstraint);
+
+ double bField = event.getDetector().getFieldMap().getField(IP)[2];
+ boolean hist = getHistogramLevel() > 0;
+
+ List<Track> trackList = new ArrayList<Track>();
+ for (MCParticle p : event.getMCParticles()) {
+ // filter for FINAL_STATE
+ if (p.getGeneratorStatus() != Particle.FINAL_STATE) {
+ continue;
+ }
+ double pCharge = p.getCharge();
+ if (pCharge == 0 || Double.isNaN(pCharge) || 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;
+ }
+
+ Track t = factory.getTrack(p);
+ trackList.add(t);
+ }
+ 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);
+ }
+
+ 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);
+ }
+ }
+
+}
Added: projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/fix/FastMCTrackFactory.java
=============================================================================
--- projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/fix/FastMCTrackFactory.java (added)
+++ projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/tracking/fix/FastMCTrackFactory.java Fri Jan 9 14:44:03 2015
@@ -0,0 +1,252 @@
+/**
+ * @version $Id: FastMCTrackFactory.java,v 1.3 2011/04/02 16:35:41 jstrube Exp $
+ */
+package org.lcsim.mc.fast.tracking.fix;
+
+import static java.lang.Math.PI;
+import static java.lang.Math.abs;
+import static java.lang.Math.sqrt;
+import static org.lcsim.event.LCIOParameters.ParameterName.omega;
+import static org.lcsim.event.LCIOParameters.ParameterName.phi0;
+import hep.physics.matrix.SymmetricMatrix;
+import hep.physics.vec.Hep3Vector;
+
+import java.io.IOException;
+import java.util.Random;
+
+import org.lcsim.conditions.ConditionsManager;
+import org.lcsim.conditions.ConditionsSet;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.LCIOParameters;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.Track;
+import org.lcsim.mc.fast.tracking.ResolutionTable;
+import org.lcsim.mc.fast.tracking.SimpleTables;
+import org.lcsim.mc.fast.tracking.TrackResolutionTables;
+import org.lcsim.spacegeom.SpacePoint;
+import org.lcsim.spacegeom.SpaceVector;
+import org.lcsim.util.swim.HelixSwimmer;
+
+import Jama.EigenvalueDecomposition;
+import Jama.Matrix;
+import Jama.util.Maths;
+
+/**
+ * @author jstrube This class creates a new FastMC Track. It is used as the interface between the track measurement and the detector. Since Track doesn't know anything about the magnetic field, and
+ * the material, it cannot transport its own parameters. Changing the reference point of a track requires swimming; it is therefore done in this class.
+ *
+ */
+public class FastMCTrackFactory {
+ private TrackResolutionTables _tables;
+ private SimpleTables _simpleTables;
+ private ConditionsManager _manager;
+ private double _Bz;
+ private HelixSwimmer _swimmer;
+ private static Random rDummy = new Random();
+ private static SpacePoint pDummy = new SpacePoint();
+
+ /**
+ * This constructor obtains the necessary information for construction like the field and the resolution tables from the event.
+ *
+ * @param event The current event
+ * @param beamConstraint A switch to obtain the resolution tables with or without beamconstraint
+ */
+ public FastMCTrackFactory(EventHeader event, boolean beamConstraint) {
+ this(event.getDetectorName(), event.getDetector().getFieldMap().getField(new double[3])[2], beamConstraint);
+ }
+
+ /**
+ * This constructor is only to be used by unit tests. It will instantiate the Factory with a detector name and a field.
+ *
+ */
+ public FastMCTrackFactory(String detectorName, double field, boolean beamConstraint) {
+ _Bz = field;
+ _manager = ConditionsManager.defaultInstance();
+ try {
+ // new detector, run 0
+ _manager.setDetector(detectorName, 0);
+ } catch (ConditionsManager.ConditionsNotFoundException e) {
+ System.err.print("Conditions for detector " + detectorName + " not found!");
+ }
+ ConditionsSet trackParameters = _manager.getConditions("TrackParameters");
+ ConditionsSet simpleTrack = _manager.getConditions("SimpleTrack");
+ try {
+ _tables = new TrackResolutionTables(trackParameters, beamConstraint);
+ _simpleTables = new SimpleTables(simpleTrack);
+ } catch (IOException e) {
+ }
+ _swimmer = new HelixSwimmer(field);
+ }
+
+ /**
+ * Creates a track from an MCParticle
+ * @param part The MCParticle that is transformed to a track
+ * @return A FastMCTrack instance that contains information about the particle that was used as input
+ */
+ public Track getTrack(MCParticle part) {
+ FastMCTrack t = (FastMCTrack) getTrack(part.getMomentum(), part.getOrigin(), (int) part.getCharge());
+ t._particle = part;
+ return t;
+ }
+
+ /**
+ * Creates a track from an MCParticle without smearing the parameters
+ * @param part The MCParticle that is transformed to a track
+ * @return A FastMCTrack instance that contains information about the particle that was used as input
+ */
+ public Track getUnsmearedTrack(MCParticle part) {
+ FastMCTrack t = (FastMCTrack) getTrack(new SpaceVector(part.getMomentum()), new SpacePoint(part.getOrigin()), pDummy, (int) part.getCharge(), rDummy, false);
+ t._particle = part;
+ return t;
+ }
+
+ /**
+ * Creates a new Track with the given parameters. See #{@link #getTrack(SpacePoint, SpacePoint, SpacePoint, int, Random)} for details.
+ *
+ * @param momentum The momentum at a given location
+ * @param location The location where the momentum is measured
+ * @param charge The charge of the Particle that created the Track
+ * @return A new NewTFastMCTrackect with the desired properties
+ */
+ public Track getTrack(SpaceVector momentum, SpacePoint location, int charge) {
+ return getTrack(momentum, location, pDummy, charge, rDummy);
+ }
+
+ /**
+ * Creates a new Track with the given parameters. See #{@link #getTrack(SpacePoint, SpacePoint, SpacePoint, int, Random)} for details.
+ *
+ * @param momentum The momentum at a given location
+ * @param location The location where the momentum is measured
+ * @param charge The charge of the Particle that created the Track
+ * @param random A random generator instance
+ * @return A new NewTrFastMCTrackct with the desired properties
+ */
+ public Track getTrack(SpaceVector momentum, SpacePoint location, int charge, Random random) {
+ return getTrack(momentum, location, pDummy, charge, random);
+ }
+
+ /**
+ * Creates a new Track with the given parameters. See #{@link #getTrack(SpacePoint, SpacePoint, SpacePoint, int, Random)} for details.
+ *
+ * @param momentum The momentum at a given location
+ * @param location The location where the momentum is measured
+ * @param referencePoint The point with respect to which the parameters are measured
+ * @param charge The charge of the Particle that created the Track
+ * @return A new NewTrFastMCTrackct with the desired properties
+ */
+ public Track getTrack(SpaceVector momentum, SpacePoint location, SpacePoint referencePoint, int charge) {
+ return getTrack(momentum, location, referencePoint, charge, rDummy);
+ }
+
+ /**
+ * This version is only to be used in unit tests.
+ *
+ * @param momentum The momentum at a given location
+ * @param location The location where the momentum is measured
+ * @param referencePoint The point with respect to which the parameters are measured
+ * @param charge The charge of the Particle that created the Track
+ * @param random A random generator instance
+ * @param shouldISmear This parameter switches smearing on/off. It should always be true except in Unit tests.
+ * @return A new NewTracFastMCTrack with the desired properties
+ */
+ public Track getTrack(SpaceVector momentum, SpacePoint location, SpacePoint referencePoint, int charge, Random random, boolean shouldISmear) {
+ _swimmer.setTrack(momentum, location, charge);
+ double alpha = _swimmer.getTrackLengthToPoint(referencePoint);
+ SpacePoint poca = _swimmer.getPointAtLength(alpha);
+ SpaceVector momentumAtPoca = _swimmer.getMomentumAtLength(alpha);
+ LCIOParameters parameters = LCIOParameters.SpaceMomentum2Parameters(poca, momentumAtPoca, referencePoint, charge, _Bz);
+ SymmetricMatrix errorMatrix = new SymmetricMatrix(5);
+ // this sets the measurement error
+ double cosTheta = abs(momentumAtPoca.cosTheta());
+ double p_mag = momentumAtPoca.magnitude();
+ ResolutionTable table = cosTheta < _tables.getPolarInner() ? _tables.getBarrelTable() : _tables.getEndcapTable();
+ for (int i = 0; i < 5; i++) {
+ for (int j = 0; j <= i; j++) {
+ double iVal = table.findTable(i, j).interpolateVal(cosTheta, p_mag);
+ errorMatrix.setElement(i, j, iVal);
+ }
+ }
+ // it's a bit inefficient to always have this condition here, although it's only used in tests.
+ LCIOParameters smearParams = shouldISmear ? smearParameters(parameters, errorMatrix, random) : parameters;
+
+ // System.out.println("Charge: " + charge);
+ // System.out.println("TrackFactory: POCA " + poca);
+ // System.out.println("TrackFactory: Momentum " + momentumAtPoca);
+ // System.out.println("TrackFactory: Parameters: " + smearParams);
+ // System.out.println("TrackFactory: POCA from parameters: " + LCIOTrackParameters.Parameters2Position(smearParams, referencePoint));
+ // System.out.println("TrackFactory: Parameters from POCA: " + LCIOTrackParameters.SpaceMomentum2Parameters(poca, momentumAtPoca, referencePoint, charge, _Bz));
+ return new FastMCTrack(referencePoint, smearParams, errorMatrix, charge);
+ }
+
+ /**
+ * Returns a new Track object initialized with the given values, and with its parameters smeared according to the Tables that are read from the detector. This method can take a random seed
+ *
+ * @param momentum The momentum at a given location
+ * @param location The location where the momentum is measured
+ * @param referencePoint The point with respect to which the parameters are measured
+ * @param charge The charge of the Particle that created the Track
+ * @param random A random generator instance
+ * @return A new NewTraFastMCTrackt with the desired properties
+ */
+ public Track getTrack(SpaceVector momentum, SpacePoint location, SpacePoint referencePoint, int charge, Random random) {
+ return getTrack(momentum, location, referencePoint, charge, random, true);
+ }
+
+ public Track getTrack(Hep3Vector momentum, Hep3Vector location, int charge) {
+ return getTrack(new SpaceVector(momentum), new SpacePoint(location), pDummy, charge, rDummy);
+ }
+
+ /**
+ * Swims the Track to a new reference point and calculates the parameters anew. It has to be done here, because it involves swimming, which has to be done outside the track
+ * @param track The track to be swum
+ * @param referencePoint The new reference point for the track to swim to
+ */
+ public void setNewReferencePoint(Track track, SpacePoint referencePoint) {
+ _swimmer.setTrack(track);
+ double alpha = _swimmer.getTrackLengthToPoint(referencePoint);
+
+ // TODO this involves transportation of the full covariance matrix.
+ // See Paul Avery's notes for details.
+ throw new RuntimeException("not yet implemented !");
+ }
+
+ /**
+ * Smears the measurement matrix with a Gaussian error
+ * @param oldParams The unsmeared Parameters
+ * @param errorMatrix The measurement error matrix
+ * @param random A random generator
+ * @return A new set of smeared parameters
+ */
+ private static LCIOParameters smearParameters(LCIOParameters oldParams, SymmetricMatrix sm, Random random) {
+ Matrix errorMatrix = Maths.toJamaMatrix(sm);
+ EigenvalueDecomposition eigen = errorMatrix.eig();
+ double[] realEigen = eigen.getRealEigenvalues();
+ double[] imaginaryEigen = eigen.getImagEigenvalues();
+ Matrix eigenValues = eigen.getV();
+ if (eigenValues.det() == 0) {
+ throw new RuntimeException("ErrorMatrix does not have orthogonal basis");
+ }
+ for (int i = 0; i < imaginaryEigen.length; i++) {
+ if (imaginaryEigen[i] != 0)
+ throw new RuntimeException("ErrorMatrix has imaginary eigenvalues");
+ }
+ Matrix x = new Matrix(5, 1);
+ for (int i = 0; i < 5; i++) {
+ if (realEigen[i] <= 0)
+ throw new RuntimeException("non-positive eigenvalue encountered");
+ x.set(i, 0, sqrt(realEigen[i]) * random.nextGaussian());
+ }
+ Matrix shift = eigenValues.times(x);
+ Matrix params = new Matrix(oldParams.getValues(), 5);
+ // calculate the new parameters
+ double[] parameters = params.plus(shift).getColumnPackedCopy();
+ double pt = oldParams.getPt() * oldParams.get(omega) / parameters[omega.ordinal()];
+ // adjust the new parameters if necessary
+ if (parameters[phi0.ordinal()] > PI) {
+ parameters[phi0.ordinal()] -= 2 * PI;
+ } else if (parameters[phi0.ordinal()] < -PI) {
+ parameters[phi0.ordinal()] += 2 * PI;
+ }
+ return new LCIOParameters(parameters, pt);
+ }
+}
Added: projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/util/CreateFinalStateMCParticleList.java
=============================================================================
--- projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/util/CreateFinalStateMCParticleList.java (added)
+++ projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/util/CreateFinalStateMCParticleList.java Fri Jan 9 14:44:03 2015
@@ -0,0 +1,365 @@
+package org.lcsim.mc.fast.util;
+
+import hep.physics.vec.Hep3Vector;
+import java.util.ArrayList;
+import java.util.List;
+import org.lcsim.event.MCParticle;
+import org.lcsim.util.Driver;
+import org.lcsim.event.EventHeader;
+import static org.lcsim.mc.fast.util.MCParticleClassifier.MCPClass;
+
+/**
+ * CreateFinalStateMCParticleList writes a list of Final state MCParticles to event. The choices are GeneratorFS(Gen) and SimulatorFS(Sim). The final state particles are those at the end of the
+ * decay/interaction chain for the choice(Gen or Sim). By default no particles created without destroying the parent are included. This can be overridden by the setKeepContinuousXXX methods. Decays
+ * and interactions beyond a radius or z position can be ignored by using the setRadiusCut and setZCut methods. The collection name of the final state particles can also be set.
+ * @author Ron Cassell
+ */
+public class CreateFinalStateMCParticleList extends Driver {
+ private MCPClass fs;
+ private String FStype;
+ private String CollectionName;
+ private double Rcut;
+ private double Zcut;
+ private boolean keepce;
+ private boolean keepcp;
+ private boolean keepch;
+ private final static String[] types = { "Gen", "Sim" };
+ private int itype;
+ private MCParticleClassifier cl;
+
+ /**
+ * Set up the defaults for creating the lists. Type parameter must be Gen or Sim.
+ *
+ * @param type
+ */
+ public CreateFinalStateMCParticleList(String type) {
+ fs = MCPClass.GEN_FINAL_STATE;
+ cl = new MCParticleClassifier();
+ FStype = type;
+ if (FStype.compareTo(types[0]) == 0) {
+ itype = 0;
+ } else if (FStype.compareTo(types[1]) == 0) {
+ itype = 1;
+ } else {
+ itype = 0;
+ System.out.println("CreateFinalStateMCParticleList created with invalid type " + FStype + ": Defaulting to " + types[0]);
+ FStype = types[0];
+ }
+ Rcut = 99999.;
+ Zcut = 99999.;
+ keepce = false;
+ keepcp = false;
+ keepch = false;
+ CollectionName = FStype + "FinalStateParticles";
+ }
+
+ /**
+ * Set the output collection name
+ *
+ * @param name
+ */
+ public void setCollectionName(String name) {
+ CollectionName = name;
+ }
+
+ /**
+ * Set a radius cut beyond which interactions and decays are ignored.
+ *
+ * @param rc
+ */
+ public void setRadiusCut(double rc) {
+ Rcut = rc;
+ }
+
+ /**
+ * Set a z cut beyond which interactions and decays are ignored.
+ *
+ * @param zc
+ */
+ public void setZCut(double zc) {
+ Zcut = zc;
+ }
+
+ /**
+ * Keep electrons created without destroying parent. (deltas, comptons)
+ */
+ public void setKeepContinuousElectrons() {
+ keepce = true;
+ }
+
+ /**
+ * Keep photons created without destroying parent. (brem, nuclear elastic and pseudo-elastic)
+ */
+ public void setKeepContinuousPhotons() {
+ keepcp = true;
+ }
+
+ /**
+ * Keep hadrons created without destroying parent. (nuclear elastic and pseudo-elastic)
+ */
+ public void setKeepContinuousHadrons() {
+ keepch = true;
+ }
+
+ /**
+ * Process and event. Create the final state particle list and store it in event.
+ *
+ * @param event
+ */
+ public void process(EventHeader event) {
+ //
+ // Create a list to hold the final state particles
+ //
+ List<MCParticle> fslist = new ArrayList<MCParticle>();
+ //
+ // Get the list of all MCParticles
+ //
+ List<MCParticle> all = (List<MCParticle>) event.get("MCParticle");
+ //
+ // If Generator final state particles requested, use the MCParticleClassifier
+ // to identify the final state particles, and add them to the list
+ //
+ if (itype == 0) {
+ boolean simulated = false;
+ for (MCParticle p : all) {
+ if (p.getSimulatorStatus().isDecayedInTracker() || p.getSimulatorStatus().isDecayedInCalorimeter() || p.getSimulatorStatus().hasLeftDetector() || p.getSimulatorStatus().isStopped()) {
+ simulated = true;
+ continue;
+ }
+ }
+ if (!simulated) {
+ for (MCParticle p : all) {
+ if (p.getGeneratorStatus() == MCParticle.FINAL_STATE)
+ fslist.add(p);
+ }
+ } else {
+ for (MCParticle p : all) {
+ if (cl.getClassification(p) == fs) {
+ fslist.add(p);
+ }
+ }
+ }
+ }
+ //
+ // Simulator final state particles requested
+ //
+ else {
+ //
+ // Check if we are keeping any continuous process created particles
+ //
+ boolean keepsome = keepce || keepcp || keepch;
+ //
+ // Loop over all MCParticles
+ //
+ for (MCParticle p : all) {
+ //
+ // Never keep backscatter particles
+ //
+ if (p.getSimulatorStatus().isBackscatter())
+ continue;
+ //
+ // Don't keep particles the Simulator never saw
+ //
+ boolean inSim = p.getSimulatorStatus().isDecayedInTracker() || p.getSimulatorStatus().isDecayedInCalorimeter() || p.getSimulatorStatus().hasLeftDetector() || p.getSimulatorStatus().isStopped();
+ if (!inSim)
+ continue;
+ //
+ // Find out if the particle has endpoint daughters
+ //
+ boolean hasepd = false;
+ for (MCParticle d : p.getDaughters()) {
+ if (d.getSimulatorStatus().isBackscatter())
+ continue;
+ if (d.getSimulatorStatus().vertexIsNotEndpointOfParent())
+ continue;
+ hasepd = true;
+ break;
+ }
+ //
+ // If the particle has endpoint daughters it is not final state
+ //
+ if (hasepd)
+ continue;
+ //
+ // This is a simulator final state particle. If it is a result of a
+ // continuous process, check to see if we keep it.
+ //
+ MCParticle pp = getFirstContinuousParticle(p);
+ if (pp != null) {
+ if (keepsome) {
+ if (keepThis(p, pp))
+ fslist.add(p);
+ }
+ }
+ //
+ // Add the final state particle to the list
+ //
+ else {
+ fslist.add(p);
+ }
+ }
+ }
+ //
+ // We now have the complete list of final state particles. Check for
+ // production vertex cuts
+ //
+ if ((Rcut < 9999.) || (Zcut < 9999.)) {
+ //
+ // Create a remove list and an add list
+ //
+ List<MCParticle> removelist = new ArrayList<MCParticle>();
+ List<MCParticle> addtolist = new ArrayList<MCParticle>();
+ //
+ // Loop over the final state list
+ //
+ for (MCParticle p : fslist) {
+ //
+ // Check vertex against cuts
+ //
+ Hep3Vector vtx = p.getOrigin();
+ if ((Math.abs(vtx.z()) > Zcut) || (Math.sqrt(vtx.x() * vtx.x() + vtx.y() * vtx.y()) > Rcut)) {
+ //
+ // Vertex failed cuts. Remove it.
+ //
+ removelist.add(p);
+ //
+ // Unless a continuous process particle, add the parent to the list of particles to add back in
+ //
+ if (!p.getSimulatorStatus().vertexIsNotEndpointOfParent()) {
+ if (!addtolist.contains(p.getParents().get(0)))
+ addtolist.add(p.getParents().get(0));
+ }
+ }
+ }
+ //
+ // Remove all the particles in the remove list from the final state list
+ //
+ for (MCParticle p : removelist) {
+ fslist.remove(p);
+ }
+ //
+ // Loop over the add list
+ //
+ for (MCParticle p : addtolist) {
+ Hep3Vector vtx = p.getOrigin();
+ MCParticle pp = p;
+ boolean hasparent = true;
+ //
+ // Trace parentage until cut is passed (or run out of parents or break chain with continuous particle)
+ //
+ while ((hasparent) && ((Math.abs(vtx.z()) > Zcut) || (Math.sqrt(vtx.x() * vtx.x() + vtx.y() * vtx.y()) > Rcut))) {
+ if (pp.getSimulatorStatus().vertexIsNotEndpointOfParent()) {
+ pp = null;
+ } else {
+ pp = pp.getParents().get(0);
+ }
+ if (pp == null) {
+ hasparent = false;
+ } else {
+ vtx = pp.getOrigin();
+ }
+ }
+ //
+ // Add the particle to the fs list if it exists
+ //
+ if (!(pp == null)) {
+ if (!fslist.contains(pp))
+ fslist.add(pp);
+ }
+ }
+ }
+ //
+ // Write the collection to event
+ //
+ event.put(CollectionName, fslist);
+ event.getMetaData(fslist).setSubset(true);
+ }
+
+ /**
+ * Decide if a particle resulting from continuous production should be kept as a final state particle.
+ *
+ * @param p - particle in question
+ * @param pp - first continuous production particle in parentage chain of p
+ */
+ public boolean keepThis(MCParticle p, MCParticle pp) {
+ boolean keepit = false;
+ //
+ // if first vneop particle is the particle in question, only need to check it for validity
+ //
+ if (p == pp) {
+ if (Math.abs(p.getPDGID()) == 11) {
+ if (keepce)
+ keepit = true;
+ } else if (p.getPDGID() == 22) {
+ if (keepcp)
+ keepit = true;
+ } else {
+ if (keepch)
+ keepit = true;
+ }
+ return keepit;
+ }
+ //
+ // otherwise, check the whole chain for validity
+ //
+ else {
+ //
+ // First check the particle in question
+ //
+ if (p.getSimulatorStatus().vertexIsNotEndpointOfParent()) {
+ if (Math.abs(p.getPDGID()) == 11) {
+ if (keepce)
+ keepit = true;
+ } else if (p.getPDGID() == 22) {
+ if (keepcp)
+ keepit = true;
+ } else {
+ if (keepch)
+ keepit = true;
+ }
+ if (!keepit)
+ return keepit;
+ }
+ //
+ // Particle itself is valid, check the chain
+ //
+ keepit = true;
+ MCParticle ppp = p;
+ while ((keepit) && (ppp != pp)) {
+ ppp = ppp.getParents().get(0);
+ if (ppp.getSimulatorStatus().vertexIsNotEndpointOfParent()) {
+ keepit = false;
+ if (Math.abs(ppp.getPDGID()) == 11) {
+ if (keepce)
+ keepit = true;
+ } else if (ppp.getPDGID() == 22) {
+ if (keepcp)
+ keepit = true;
+ } else {
+ if (keepch)
+ keepit = true;
+ }
+ }
+ }
+ return keepit;
+ }
+ }
+
+ /**
+ * Find and return the top VNEOP particle in parentage chain
+ *
+ * @param p - particle in question
+ */
+ public MCParticle getFirstContinuousParticle(MCParticle p) {
+ MCParticle rp = null;
+ if (p.getSimulatorStatus().vertexIsNotEndpointOfParent())
+ rp = p;
+ MCParticle pp = p;
+ while (pp.getParents().size() == 1) {
+ pp = pp.getParents().get(0);
+ if (pp.getSimulatorStatus().vertexIsNotEndpointOfParent())
+ rp = pp;
+ }
+ return rp;
+ }
+}
Added: projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/util/MCParticleClassifier.java
=============================================================================
--- projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/util/MCParticleClassifier.java (added)
+++ projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/util/MCParticleClassifier.java Fri Jan 9 14:44:03 2015
@@ -0,0 +1,91 @@
+package org.lcsim.mc.fast.util;
+
+import java.util.List;
+import org.lcsim.event.MCParticle;
+
+/**
+ * MCParticleClassifier has a single method, getClassification(MCParticle), that returns an enum classifying the particle.
+ *
+ * Possible return values are:
+ *
+ * GEN_INITIAL -> Generator particle the simulator never saw
+ *
+ * GEN_PREDECAY -> Generator particle passed to the simulator along with a predecay, and the simulator generated the decay products.
+ *
+ * GEN_FINAL_STATE -> Generator final state particle: either FinalState in the input file, or intermediate passed to the Simulator with a predecay that never happened.
+ *
+ * SIM_BACKSCATTER -> Simulator particle produced as a result of backscatter from a non-tracking region.
+ *
+ * SIM_VERTEX_NOT_PARENT_ENDPOINT -> Simulator particle produced without destroying parent.
+ *
+ * SIM_INTERACTED_OR_DECAYED -> Simulator particle produced as a result of an interaction or decay in a tracking region.
+ * @version $Id: MCParticleClassifier.java,v 1.5 2006/06/28 04:48:33 jstrube Exp $
+ */
+public class MCParticleClassifier {
+ public enum MCPClass {
+ GEN_INITIAL, GEN_PREDECAY, GEN_FINAL_STATE, SIM_BACKSCATTER, SIM_VERTEX_NOT_PARENT_ENDPOINT, SIM_INTERACTED_OR_DECAYED
+ };
+
+ public static MCPClass getClassification(MCParticle p) {
+ //
+ // Check if particle existed in Generator
+ //
+ if (p.getGeneratorStatus() > 0) {
+ //
+ // Remove the documentation particles in case the generator left them floating
+ //
+ if (p.getGeneratorStatus() == MCParticle.DOCUMENTATION)
+ return MCPClass.GEN_INITIAL;
+ //
+ // Check if the generator particle has daughters in the generator
+ //
+ boolean hasGeneratorDaughters = false;
+ List<MCParticle> daughters = p.getDaughters();
+ for (MCParticle d : daughters) {
+ if (d.getGeneratorStatus() > 0)
+ hasGeneratorDaughters = true;
+ }
+ //
+ // If no generator daughters return final state
+ //
+ if (!hasGeneratorDaughters) {
+ //
+ // Trap for strange generator problem: Hanging quark with INTERMEDIATE status
+ // Therefore must check that Simulator saw this particle, otherwise return
+ // GEN_INITIAL
+ //
+ boolean inSim = p.getSimulatorStatus().isDecayedInTracker() || p.getSimulatorStatus().isDecayedInCalorimeter() || p.getSimulatorStatus().hasLeftDetector() || p.getSimulatorStatus().isStopped();
+ if (inSim) {
+ return MCPClass.GEN_FINAL_STATE;
+ } else {
+ return MCPClass.GEN_INITIAL;
+ }
+ }
+ //
+ // A generator particle with gerenrator daughters:
+ // if the simulator saw it must be a predecay
+ //
+ if (p.getSimulatorStatus().isDecayedInTracker() | p.getSimulatorStatus().isDecayedInCalorimeter()) {
+ return MCPClass.GEN_PREDECAY;
+ }
+ //
+ // simulator never saw it: initial
+ //
+ return MCPClass.GEN_INITIAL;
+ } else {
+ //
+ // Particle didn't exist in generator: choices are backscatter, production without killing parent,
+ // or the result of an interaction or decay( killing parent)
+ //
+ if (p.getSimulatorStatus().isBackscatter()) {
+ return MCPClass.SIM_BACKSCATTER;
+ }
+
+ if (p.getSimulatorStatus().vertexIsNotEndpointOfParent()) {
+ return MCPClass.SIM_VERTEX_NOT_PARENT_ENDPOINT;
+ }
+
+ return MCPClass.SIM_INTERACTED_OR_DECAYED;
+ }
+ }
+}
Added: projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/util/MonitorStdhep.java
=============================================================================
--- projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/util/MonitorStdhep.java (added)
+++ projects/lcsim/trunk/mc/src/main/java/org/lcsim/mc/fast/util/MonitorStdhep.java Fri Jan 9 14:44:03 2015
@@ -0,0 +1,81 @@
+package org.lcsim.mc.fast.util;
+
+import static java.lang.Math.abs;
+
+import hep.physics.particle.properties.ParticlePropertyManager;
+import hep.physics.particle.properties.ParticlePropertyProvider;
+
+import java.io.IOException;
+import java.util.List;
+import java.util.Map;
+
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.util.Driver;
+
+public class MonitorStdhep extends Driver {
+
+ final int nEventSkip = 0;
+ // final int nEventSkip=95;
+ // final int nEventSkip=233;
+ // final int nEventSkip=684;
+ final int nprintMax = 20;
+ // final int nprintMax=Integer.MAX_VALUE;
+ final int nCloudMax = Integer.MAX_VALUE;
+ int ncnt = 0;
+ int ncnt_in = 0;
+ int ncnt_out = 0;
+ int eventNumber;
+ float eventWeight;
+ final ParticlePropertyProvider dPPP;
+ List<ReconstructedParticle> recoParticles;
+ int idrup;
+
+ int m_eventCount;
+
+ int ievt;
+ int nmax = 1000000;
+
+ public MonitorStdhep() throws IOException {
+
+ m_eventCount = 0;
+ dPPP = ParticlePropertyManager.getParticlePropertyProvider();
+ System.out.println(" MonitorStdhep constructor ");
+ }
+
+ protected void process(EventHeader event) {
+ // super.process(event);
+ ncnt_in++;
+ if (ncnt_in <= nEventSkip) {
+ // System.out.println(" event.getEventNumber= "+event.getEventNumber());
+ throw new Driver.NextEventException();
+ } else {
+ eventNumber = event.getEventNumber();
+ ncnt++;
+ if (ncnt <= nprintMax || ncnt % 100 == 0)
+ System.out.println(" ncnt_in= " + ncnt_in + " ncnt= " + ncnt + " eventNumber= " + eventNumber);
+ // super.process(event);
+ // idrup=event.getIntegerParameters().get("idrup")[0];
+
+ eventWeight = event.getWeight();
+ if (ncnt <= nprintMax) {
+ Map<String, float[]> headerFloatMap = event.getFloatParameters();
+ Map<String, int[]> headerIntMap = event.getIntegerParameters();
+ for (String headerFloatName : headerFloatMap.keySet())
+ System.out.println(" headerFloatName= " + headerFloatName + " value= " + headerFloatMap.get(headerFloatName)[0]);
+ for (String headerIntName : headerIntMap.keySet())
+ System.out.println(" headerIntName= " + headerIntName + " value= " + headerIntMap.get(headerIntName)[0]);
+ System.out.println(" idrup= " + idrup);
+ long eventTimeStamp = event.getTimeStamp();
+ System.out.println(" eventWeight= " + eventWeight);
+ }
+ List<MCParticle> particles = event.get(MCParticle.class, event.MC_PARTICLES);
+ for (MCParticle particle : particles) {
+ int iPdgId = abs(particle.getPDGID());
+ if (ncnt <= nprintMax)
+ System.out.println(" iPdgId= " + iPdgId);
+ }
+ }
+ }
+}
Modified: projects/lcsim/trunk/parent/pom.xml
=============================================================================
--- projects/lcsim/trunk/parent/pom.xml (original)
+++ projects/lcsim/trunk/parent/pom.xml Fri Jan 9 14:44:03 2015
@@ -174,6 +174,11 @@
<artifactId>lcsim-steering-files</artifactId>
<version>3.0.8-SNAPSHOT</version>
</dependency>
+ <dependency>
+ <groupId>org.lcsim</groupId>
+ <artifactId>lcsim-mc</artifactId>
+ <version>3.0.8-SNAPSHOT</version>
+ </dependency>
<!-- End org.lcsim dependencies. -->
<dependency>
<groupId>org.freehep</groupId>
Modified: projects/lcsim/trunk/pom.xml
=============================================================================
--- projects/lcsim/trunk/pom.xml (original)
+++ projects/lcsim/trunk/pom.xml Fri Jan 9 14:44:03 2015
@@ -139,6 +139,7 @@
<module>event-processing</module>
<module>job-manager</module>
<module>math</module>
+ <module>mc</module>
<module>plugin</module>
<module>recon-drivers</module>
<module>tracking</module>
@@ -163,6 +164,7 @@
<module>detector-framework</module>
<module>event-model</module>
<module>event-processing</module>
+ <module>mc</module>
<module>job-manager</module>
<module>tracking</module>
<module>vertexing</module>
########################################################################
Use REPLY-ALL to reply to list
To unsubscribe from the LCDET-SVN list, click the following link:
https://listserv.slac.stanford.edu/cgi-bin/wa?SUBED1=LCDET-SVN&A=1
|
