Author: [log in to unmask]
Date: Wed Jan 6 08:21:16 2016
New Revision: 4091
Log:
actually add the stuff in my sandbox
Added:
java/sandbox/DatabaseConditionsManager.java
java/sandbox/HPSTrack.java
java/sandbox/MultipleScattering.java
java/sandbox/SvtBiasConstant.java
java/sandbox/SvtOldHeaderAnalysisDriver.java
java/sandbox/SvtOldHeaderDataInfo.java
java/sandbox/TrackUtils.java
Added: java/sandbox/DatabaseConditionsManager.java
=============================================================================
--- java/sandbox/DatabaseConditionsManager.java (added)
+++ java/sandbox/DatabaseConditionsManager.java Wed Jan 6 08:21:16 2016
@@ -0,0 +1,1134 @@
+package org.hps.conditions.database;
+
+import java.io.File;
+import java.io.FileInputStream;
+import java.io.FileNotFoundException;
+import java.io.IOException;
+import java.io.InputStream;
+import java.sql.Connection;
+import java.sql.DriverManager;
+import java.sql.PreparedStatement;
+import java.sql.ResultSet;
+import java.sql.SQLException;
+import java.sql.Statement;
+import java.util.ArrayList;
+import java.util.LinkedHashSet;
+import java.util.List;
+import java.util.Set;
+import java.util.logging.Level;
+import java.util.logging.Logger;
+
+import org.hps.conditions.api.AbstractConditionsObjectConverter;
+import org.hps.conditions.api.ConditionsObject;
+import org.hps.conditions.api.ConditionsObjectCollection;
+import org.hps.conditions.api.ConditionsObjectException;
+import org.hps.conditions.api.ConditionsRecord;
+import org.hps.conditions.api.ConditionsRecord.ConditionsRecordCollection;
+import org.hps.conditions.api.ConditionsSeries;
+import org.hps.conditions.api.TableMetaData;
+import org.hps.conditions.api.TableRegistry;
+import org.hps.conditions.ecal.EcalConditions;
+import org.hps.conditions.ecal.EcalConditionsConverter;
+import org.hps.conditions.ecal.TestRunEcalConditionsConverter;
+import org.hps.conditions.svt.SvtConditions;
+import org.hps.conditions.svt.SvtConditionsConverter;
+import org.hps.conditions.svt.SvtDetectorSetup;
+import org.hps.conditions.svt.TestRunSvtConditionsConverter;
+import org.jdom.Document;
+import org.jdom.Element;
+import org.jdom.JDOMException;
+import org.jdom.input.SAXBuilder;
+import org.lcsim.conditions.ConditionsConverter;
+import org.lcsim.conditions.ConditionsManager;
+import org.lcsim.conditions.ConditionsManagerImplementation;
+import org.lcsim.geometry.Detector;
+import org.lcsim.geometry.Subdetector;
+import org.lcsim.util.log.DefaultLogFormatter;
+import org.lcsim.util.log.LogUtil;
+import org.lcsim.util.loop.DetectorConditionsConverter;
+
+/**
+ * <p>
+ * This class provides the top-level API for accessing database conditions, as well as configuring the database
+ * connection, initializing all required components, and loading required converters and table meta data. It is
+ * registered as the global <code>ConditionsManager</code> in the constructor.
+ * <p>
+ * Differences between Test Run and Engineering Run configurations are handled automatically.
+ *
+ * @see org.lcsim.conditions.ConditionsManager
+ * @author Jeremy McCormick, SLAC
+ */
+@SuppressWarnings("rawtypes")
+public final class DatabaseConditionsManager extends ConditionsManagerImplementation {
+
+ /**
+ * Name of system property that can be used to specify custom database connection parameters.
+ */
+ private static final String CONNECTION_PROPERTY = "org.hps.conditions.connection.file";
+
+ /**
+ * The default XML config.
+ */
+ private static final String DEFAULT_CONFIG = "/org/hps/conditions/config/conditions_database_prod.xml";
+
+ /**
+ * The connection properties resource for connecting to the default JLAB database.
+ */
+ private static final String DEFAULT_CONNECTION_PROPERTIES_RESOURCE = "/org/hps/conditions/config/jlab_dev_connection.prop";
+ //private static final String DEFAULT_CONNECTION_PROPERTIES_RESOURCE = "/org/hps/conditions/config/jlab_connection.prop";
+
+ /**
+ * The Eng Run XML config.
+ */
+ private static final String ENGRUN_CONFIG = "/org/hps/conditions/config/conditions_database_engrun.xml";
+
+ /**
+ * Initialize the logger.
+ */
+ private static Logger logger = LogUtil.create(DatabaseConditionsManager.class.getName(), new DefaultLogFormatter(),
+ Level.FINE);
+
+ /**
+ * The Test Run XML config.
+ */
+ private static final String TEST_RUN_CONFIG = "/org/hps/conditions/config/conditions_database_testrun_2012.xml";
+
+ /**
+ * The max value for a run to be considered Test Run.
+ */
+ private static final int TEST_RUN_MAX_RUN = 1365;
+
+ static {
+ // Set default login timeout of 5 seconds.
+ DriverManager.setLoginTimeout(5);
+ }
+
+ /**
+ * Get the static instance of this class.
+ *
+ * @return the static instance of the manager
+ */
+ public static synchronized DatabaseConditionsManager getInstance() {
+
+ // Is there no manager installed yet?
+ if (!ConditionsManager.isSetup() || !(ConditionsManager.defaultInstance() instanceof DatabaseConditionsManager)) {
+ logger.finest("creating new DatabaseConditionsManager instance");
+ // Create a new instance if necessary, which will install it globally as the default.
+ new DatabaseConditionsManager();
+ }
+
+ // Get the instance back from the default conditions system and check that the type is correct now.
+ final ConditionsManager manager = ConditionsManager.defaultInstance();
+ if (!(manager instanceof DatabaseConditionsManager)) {
+ logger.severe("default conditions manager has wrong type: " + manager.getClass());
+ throw new RuntimeException("Default conditions manager has the wrong type: "
+ + ConditionsManager.defaultInstance().getClass().getName());
+ }
+
+ return (DatabaseConditionsManager) manager;
+ }
+
+ /**
+ * Get the Logger for this class, which can be used by related sub-classes if they do not have their own logger.
+ *
+ * @return the Logger for this class
+ */
+ public static Logger getLogger() {
+ return logger;
+ }
+
+ /**
+ * Utility method to determine if a run number is from the 2012 Test Run.
+ *
+ * @param runNumber the run number
+ * @return <code>true</code> if run number is from the Test Run
+ */
+ public static boolean isTestRun(final int runNumber) {
+ return runNumber > 0 && runNumber <= TEST_RUN_MAX_RUN;
+ }
+
+ /**
+ * Reset the global static instance of the conditions manager to a new object.
+ */
+ public static synchronized void resetInstance() {
+ logger.finest("DatabaseConditionsManager instance is being reset");
+ new DatabaseConditionsManager();
+ }
+
+ /**
+ * True to cache all known conditions sets (from keys) during initialization.
+ */
+ private boolean cacheAllConditions = false;
+
+ /**
+ * True to close the connection after initialization.
+ */
+ private boolean closeConnectionAfterInitialize = true;
+
+ /**
+ * The current database connection.
+ */
+ private Connection connection;
+
+ /**
+ * The current connection parameters.
+ */
+ private ConnectionParameters connectionParameters;
+
+ /**
+ * The connection properties file, if one is being used from the command line.
+ */
+ private File connectionPropertiesFile;
+
+ /**
+ * Create the global registry of conditions object converters.
+ */
+ private final ConverterRegistry converters = ConverterRegistry.create();
+
+ /**
+ * The converter for creating the combined ECAL conditions object.
+ */
+ private ConditionsConverter ecalConverter;
+
+ /**
+ * The default ECAL detector name in the detector geometry.
+ */
+ private String ecalName = "Ecal";
+
+ /**
+ * True to freeze the system after initialization.
+ */
+ private boolean freezeAfterInitialize = false;
+
+ /**
+ * True if the conditions manager was configured from an XML configuration resource or file.
+ */
+ private boolean isConfigured = false;
+
+ /**
+ * True if manager is connected to the database.
+ */
+ private boolean isConnected = false;
+
+ /**
+ * True if the conditions system has been frozen and will ignore updates after it is initialized.
+ */
+ private boolean isFrozen = false;
+
+ /**
+ * True if the manager has been initialized, e.g. the {@link #setDetector(String, int)} method was called.
+ */
+ private boolean isInitialized = false;
+
+ /**
+ * True if current run number is from Test Run.
+ */
+ private boolean isTestRun = false;
+
+ /**
+ * Flag used to print connection parameters one time.
+ */
+ private boolean loggedConnectionParameters = false;
+
+ /**
+ * True to setup the SVT detector model with conditions.
+ */
+ private boolean setupSvtDetector = true;
+
+ /**
+ * The converter for creating the combined SVT conditions object.
+ */
+ private ConditionsConverter svtConverter;
+
+ /**
+ * The default SVT name in the detector geometry.
+ */
+ private String svtName = "Tracker";
+
+ /**
+ * The helper for setting up the SVT detector with its conditions information.
+ */
+ private final SvtDetectorSetup svtSetup = new SvtDetectorSetup(this.svtName);
+
+ /**
+ * Create the global registry of table meta data.
+ */
+ private final TableRegistry tableRegistry = TableRegistry.getTableRegistry();
+
+ /**
+ * The currently active conditions tag.
+ */
+ private String tag = null;
+
+ /**
+ * Class constructor. Calling this will automatically register this manager as the global default.
+ */
+ private DatabaseConditionsManager() {
+ this.registerConditionsConverter(new DetectorConditionsConverter());
+ this.setupConnectionFromSystemProperty();
+ ConditionsManager.setDefaultConditionsManager(this);
+ this.setRun(-1);
+ for (final AbstractConditionsObjectConverter converter : this.converters.values()) {
+ // logger.fine("registering converter for " + converter.getType());
+ this.registerConditionsConverter(converter);
+ }
+ this.addConditionsListener(this.svtSetup);
+ }
+
+ /**
+ * Cache conditions sets for all known tables.
+ */
+ private void cacheConditionsSets() {
+ for (final TableMetaData meta : this.tableRegistry.values()) {
+ try {
+ logger.fine("caching conditions " + meta.getKey() + " with type "
+ + meta.getCollectionClass().getCanonicalName());
+ this.getCachedConditions(meta.getCollectionClass(), meta.getKey());
+ } catch (final Exception e) {
+ logger.warning("could not cache conditions " + meta.getKey());
+ }
+ }
+ }
+
+ /**
+ * Close the database connection.
+ */
+ public synchronized void closeConnection() {
+ logger.fine("closing connection");
+ if (this.connection != null) {
+ try {
+ if (!this.connection.isClosed()) {
+ this.connection.close();
+ }
+ } catch (final SQLException e) {
+ throw new RuntimeException(e);
+ }
+ }
+ this.connection = null;
+ this.isConnected = false;
+ logger.fine("connection closed");
+ }
+
+ /**
+ * Close the database connection but only if there was a connection opened based on the flag. Otherwise, it should
+ * be left open. Used in conjunction with return value of {@link #openConnection()}.
+ *
+ * @param connectionOpened <code>true</code> to close the connection; <code>false</code> to leave it open
+ */
+ public synchronized void closeConnection(final boolean connectionOpened) {
+ if (connectionOpened) {
+ this.closeConnection();
+ }
+ }
+
+ /**
+ * This method will return <code>true</code> if the given collection ID already exists in the table.
+ *
+ * @param tableName the name of the table
+ * @param collectionID the collection ID value
+ * @return <code>true</code> if collection exists
+ */
+ public boolean collectionExists(final String tableName, final int collectionID) {
+ final String sql = "SELECT * FROM " + tableName + " where collection_id = " + collectionID;
+ final ResultSet resultSet = this.selectQuery(sql);
+ try {
+ resultSet.last();
+ } catch (final SQLException e) {
+ e.printStackTrace();
+ }
+ int rowCount = 0;
+ try {
+ rowCount = resultSet.getRow();
+ } catch (final SQLException e) {
+ e.printStackTrace();
+ }
+ return rowCount != 0;
+ }
+
+ /**
+ * Configure this class from an <code>InputStream</code> which should point to an XML document.
+ *
+ * @param in the InputStream which should be in XML format
+ */
+ private void configure(final InputStream in) {
+ if (!this.isConfigured) {
+ final SAXBuilder builder = new SAXBuilder();
+ Document config = null;
+ try {
+ config = builder.build(in);
+ } catch (JDOMException | IOException e) {
+ throw new RuntimeException(e);
+ }
+ this.loadConfiguration(config);
+ try {
+ in.close();
+ } catch (final IOException e) {
+ logger.warning(e.getMessage());
+ }
+ this.isConfigured = true;
+ } else {
+ logger.warning("System is already configured, so call to configure is ignored!");
+ }
+ }
+
+ /**
+ * Find a collection of conditions validity records by key name. The key name is distinct from the table name, but
+ * they are usually set to the same value.
+ *
+ * @param name the conditions key name
+ * @return the set of matching conditions records
+ */
+ public ConditionsRecordCollection findConditionsRecords(final String name) {
+ final ConditionsRecordCollection runConditionsRecords = this.getCachedConditions(
+ ConditionsRecordCollection.class, "conditions").getCachedData();
+ logger.fine("searching for conditions with name " + name + " in " + runConditionsRecords.size() + " records");
+ final ConditionsRecordCollection foundConditionsRecords = new ConditionsRecordCollection();
+ for (final ConditionsRecord record : runConditionsRecords) {
+ if (record.getName().equals(name)) {
+ if (this.matchesTag(record)) {
+ try {
+ foundConditionsRecords.add(record);
+ } catch (final ConditionsObjectException e) {
+ throw new RuntimeException(e);
+ }
+ logger.finer("found matching conditions record " + record.getRowId());
+ } else {
+ logger.finer("conditions record " + record.getRowId() + " rejected from non-matching tag "
+ + record.getTag());
+ }
+ }
+ }
+ logger.fine("found " + foundConditionsRecords.size() + " conditions records matching tag " + this.tag);
+ return foundConditionsRecords;
+ }
+
+ /**
+ * Find table information from the collection type.
+ *
+ * @param type the collection type
+ * @return the table information or <code>null</code> if does not exist
+ */
+ public List<TableMetaData> findTableMetaData(final Class<?> type) {
+ return this.tableRegistry.findByCollectionType(type);
+ }
+
+ /**
+ * Find table information from the name.
+ *
+ * @param name the name of the table
+ * @return the table information or <code>null</code> if does not exist
+ */
+ public TableMetaData findTableMetaData(final String name) {
+ return this.tableRegistry.findByTableName(name);
+ }
+
+ /**
+ * This method can be called to "freeze" the conditions system so that any subsequent updates to run number or
+ * detector name will be ignored.
+ */
+ public synchronized void freeze() {
+ if (this.getDetector() != null && this.getRun() != -1) {
+ this.isFrozen = true;
+ logger.config("conditions system is frozen");
+ } else {
+ logger.warning("conditions system cannot be frozen because it is not initialized yet");
+ }
+ }
+
+ /**
+ * Add a row for a new collection and return the new collection ID assigned to it.
+ *
+ * @param tableName the name of the table
+ * @param comment an optional comment about this new collection
+ * @return the collection's ID
+ * @throws SQLException
+ */
+ public synchronized int getCollectionId(final ConditionsObjectCollection<?> collection, final String description)
+ throws SQLException {
+
+ final String caller = Thread.currentThread().getStackTrace()[2].getClassName();
+ final String log = "created by " + System.getProperty("user.name") + " using "
+ + caller.substring(caller.lastIndexOf('.') + 1);
+ final boolean opened = this.openConnection();
+ PreparedStatement statement = null;
+ ResultSet resultSet = null;
+ int collectionId = -1;
+ try {
+ statement = this.connection.prepareStatement(
+ "INSERT INTO collections (table_name, log, description, created) VALUES (?, ?, ?, NOW())",
+ Statement.RETURN_GENERATED_KEYS);
+ statement.setString(1, collection.getTableMetaData().getTableName());
+ statement.setString(2, log);
+ if (description == null) {
+ statement.setNull(3, java.sql.Types.VARCHAR);
+ } else {
+ statement.setString(3, description);
+ }
+ statement.execute();
+ resultSet = statement.getGeneratedKeys();
+ resultSet.next();
+ collectionId = resultSet.getInt(1);
+ } finally {
+ if (resultSet != null) {
+ resultSet.close();
+ }
+ if (statement != null) {
+ statement.close();
+ }
+ this.closeConnection(opened);
+ }
+ collection.setCollectionId(collectionId);
+ return collectionId;
+ }
+
+ /**
+ * Get a list of all the {@link ConditionsRecord} objects.
+ *
+ * @return the list of all the {@link ConditionsRecord} objects
+ */
+ // FIXME: This should use a cache that is created during initialization, rather than look these up every time.
+ public ConditionsRecordCollection getConditionsRecords() {
+ logger.finer("getting conditions records ...");
+ final ConditionsRecordCollection conditionsRecords = new ConditionsRecordCollection();
+ for (final TableMetaData tableMetaData : this.tableRegistry.values()) {
+ try {
+ final ConditionsRecordCollection foundConditionsRecords = this.findConditionsRecords(tableMetaData
+ .getKey());
+ logger.finer("found " + foundConditionsRecords.size() + " collections with name "
+ + tableMetaData.getKey());
+ conditionsRecords.addAll(foundConditionsRecords);
+ } catch (final Exception e) {
+ e.printStackTrace();
+ logger.warning(e.getMessage());
+ }
+ }
+ logger.finer("found " + conditionsRecords + " conditions records");
+ logger.getHandlers()[0].flush();
+ return conditionsRecords;
+ }
+
+ /**
+ * Get a conditions series with one or more collections.
+ *
+ * @param collectionType the type of the collection
+ * @param tableName the name of the data table
+ * @param <ObjectType> the type of the conditions object
+ * @param <CollectionType> the type of the conditions collection
+ * @return the conditions series
+ */
+ @SuppressWarnings("unchecked")
+ public <ObjectType extends ConditionsObject, CollectionType extends ConditionsObjectCollection<ObjectType>> ConditionsSeries<ObjectType, CollectionType> getConditionsSeries(
+ final Class<CollectionType> collectionType, final String tableName) {
+
+ final TableMetaData metaData = this.tableRegistry.get(tableName);
+ if (metaData == null) {
+ throw new IllegalArgumentException("No table metadata found for type " + collectionType.getName());
+ }
+ if (!metaData.getCollectionClass().equals(collectionType)) {
+ throw new IllegalArgumentException("The type " + collectionType.getName() + " does not match the class "
+ + metaData.getCollectionClass().getName() + " from the meta data");
+ }
+ final Class<? extends ConditionsObject> objectType = metaData.getObjectClass();
+ final ConditionsSeriesConverter<ObjectType, CollectionType> converter = new ConditionsSeriesConverter(
+ objectType, collectionType);
+ return converter.createSeries(tableName);
+ }
+
+ /**
+ * Get the JDBC connection.
+ *
+ * @return the JDBC connection
+ */
+ public Connection getConnection() {
+ if (!this.isConnected()) {
+ this.openConnection();
+ }
+ return this.connection;
+ }
+
+ /**
+ * Get the current LCSim compact <code>Detector</code> object with the geometry and detector model.
+ *
+ * @return the detector object
+ */
+ public Detector getDetectorObject() {
+ return this.getCachedConditions(Detector.class, "compact.xml").getCachedData();
+ }
+
+ /**
+ * Get the combined ECAL conditions for this run.
+ *
+ * @return the combined ECAL conditions
+ */
+ public EcalConditions getEcalConditions() {
+ return this.getCachedConditions(EcalConditions.class, "ecal_conditions").getCachedData();
+ }
+
+ /**
+ * Get the name of the ECAL in the detector geometry.
+ *
+ * @return the name of the ECAL
+ */
+ public String getEcalName() {
+ return this.ecalName;
+ }
+
+ /**
+ * Get the subdetector object of the ECAL.
+ *
+ * @return the ECAL subdetector
+ */
+ public Subdetector getEcalSubdetector() {
+ return this.getDetectorObject().getSubdetector(this.ecalName);
+ }
+
+ /**
+ * Get the combined SVT conditions for this run.
+ *
+ * @return the combined SVT conditions
+ */
+ public SvtConditions getSvtConditions() {
+ return this.getCachedConditions(SvtConditions.class, "svt_conditions").getCachedData();
+ }
+
+ /**
+ * Get the set of available conditions tags from the conditions table
+ *
+ * @return the set of available conditions tags
+ */
+ public Set<String> getTags() {
+ logger.fine("getting list of available conditions tags");
+ final boolean openedConnection = this.openConnection();
+ final Set<String> tags = new LinkedHashSet<String>();
+ final ResultSet rs = this
+ .selectQuery("select distinct(tag) from conditions where tag is not null order by tag");
+ try {
+ while (rs.next()) {
+ tags.add(rs.getString(1));
+ }
+ } catch (final SQLException e) {
+ throw new RuntimeException(e);
+ }
+ try {
+ rs.close();
+ } catch (final SQLException e) {
+ logger.log(Level.WARNING, "error closing ResultSet", e);
+ }
+ final StringBuffer sb = new StringBuffer();
+ sb.append("found unique conditions tags: ");
+ for (final String tag : tags) {
+ sb.append(tag + " ");
+ }
+ sb.setLength(sb.length() - 1);
+ logger.fine(sb.toString());
+ this.closeConnection(openedConnection);
+ return tags;
+ }
+
+ /**
+ * True if there is a conditions record with the given name.
+ *
+ * @param name the conditions record name (usually will match to table name)
+ * @return <code>true</code> if a conditions record exists with the given name
+ */
+ public boolean hasConditionsRecord(final String name) {
+ return this.findConditionsRecords(name).size() != 0;
+ }
+
+ /**
+ * Perform all necessary initialization, including setup of the XML configuration and loading of conditions onto the
+ * Detector. This is called from the {@link #setDetector(String, int)} method to setup the manager for a new run or
+ * detector.
+ *
+ * @param detectorName the name of the detector model
+ * @param runNumber the run number
+ * @throws ConditionsNotFoundException if there is a conditions system error
+ */
+ private void initialize(final String detectorName, final int runNumber) throws ConditionsNotFoundException {
+
+ logger.config("initializing with detector " + detectorName + " and run " + runNumber);
+
+ // Is not configured yet?
+ if (!this.isConfigured) {
+ if (isTestRun(runNumber)) {
+ // This looks like the Test Run so use the custom configuration for it.
+ this.setXmlConfig(DatabaseConditionsManager.TEST_RUN_CONFIG);
+ } else if (runNumber > TEST_RUN_MAX_RUN) {
+ // Run numbers greater than max of Test Run assumed to be Eng Run (for now!).
+ this.setXmlConfig(DatabaseConditionsManager.ENGRUN_CONFIG);
+ } else if (runNumber == 0) {
+ // Use the default configuration because the run number is basically meaningless.
+ this.setXmlConfig(DatabaseConditionsManager.DEFAULT_CONFIG);
+ }
+ }
+
+ // Register the converters for this initialization.
+ logger.fine("registering converters");
+ this.registerConverters();
+
+ // Enable or disable the setup of the SVT detector.
+ logger.fine("enabling SVT setup: " + this.setupSvtDetector);
+ this.svtSetup.setEnabled(this.setupSvtDetector);
+
+ // Open the database connection.
+ this.openConnection();
+
+ // Call the super class's setDetector method to construct the detector object and activate conditions listeners.
+ logger.fine("activating default conditions manager");
+ super.setDetector(detectorName, runNumber);
+
+ // Should all conditions sets be cached?
+ if (this.cacheAllConditions) {
+ // Cache the conditions sets of all registered converters.
+ logger.fine("caching all conditions sets ...");
+ this.cacheConditionsSets();
+ }
+
+ if (this.closeConnectionAfterInitialize) {
+ logger.fine("closing connection after initialization");
+ // Close the connection.
+ this.closeConnection();
+ }
+
+ // Should the conditions system be frozen now?
+ if (this.freezeAfterInitialize) {
+ // Freeze the conditions system so subsequent updates will be ignored.
+ this.freeze();
+ logger.config("system was frozen after initialization");
+ }
+
+ this.isInitialized = true;
+
+ logger.info("conditions system initialized successfully");
+
+ // Flush logger after initialization.
+ logger.getHandlers()[0].flush();
+ }
+
+ /**
+ * Check if connected to the database.
+ *
+ * @return <code>true</code> if connected
+ */
+ public boolean isConnected() {
+ return this.isConnected;
+ }
+
+ /**
+ * True if conditions system is frozen
+ *
+ * @return <code>true</code> if conditions system is currently frozen
+ */
+ public boolean isFrozen() {
+ return this.isFrozen;
+ }
+
+ /**
+ * True if conditions manager is properly initialized.
+ *
+ * @return <code>true</code> if the manager is initialized
+ */
+ public boolean isInitialized() {
+ return this.isInitialized;
+ }
+
+ /**
+ * Return <code>true</code> if Test Run configuration is active
+ *
+ * @return <code>true</code> if Test Run configuration is active
+ */
+ public boolean isTestRun() {
+ return this.isTestRun;
+ }
+
+ /**
+ * Load configuration information from an XML document.
+ *
+ * @param document the XML document
+ */
+ private void loadConfiguration(final Document document) {
+
+ final Element node = document.getRootElement().getChild("configuration");
+
+ if (node == null) {
+ return;
+ }
+
+ Element element = node.getChild("setupSvtDetector");
+ if (element != null) {
+ this.setupSvtDetector = Boolean.parseBoolean(element.getText());
+ logger.config("setupSvtDetector = " + this.setupSvtDetector);
+ }
+
+ element = node.getChild("ecalName");
+ if (element != null) {
+ this.setEcalName(element.getText());
+ }
+
+ element = node.getChild("svtName");
+ if (element != null) {
+ this.setSvtName(element.getText());
+ }
+
+ element = node.getChild("freezeAfterInitialize");
+ if (element != null) {
+ this.freezeAfterInitialize = Boolean.parseBoolean(element.getText());
+ logger.config("freezeAfterInitialize = " + this.freezeAfterInitialize);
+ }
+
+ element = node.getChild("cacheAllCondition");
+ if (element != null) {
+ this.cacheAllConditions = Boolean.parseBoolean(element.getText());
+ logger.config("cacheAllConditions = " + this.cacheAllConditions);
+ }
+
+ element = node.getChild("isTestRun");
+ if (element != null) {
+ this.isTestRun = Boolean.parseBoolean(element.getText());
+ logger.config("isTestRun = " + this.isTestRun);
+ }
+
+ element = node.getChild("logLevel");
+ if (element != null) {
+ this.setLogLevel(Level.parse(element.getText()));
+ }
+
+ element = node.getChild("closeConnectionAfterInitialize");
+ if (element != null) {
+ this.closeConnectionAfterInitialize = Boolean.parseBoolean(element.getText());
+ logger.config("closeConnectionAfterInitialize = " + this.closeConnectionAfterInitialize);
+ }
+
+ element = node.getChild("loginTimeout");
+ if (element != null) {
+ final Integer timeout = Integer.parseInt(element.getText());
+ DriverManager.setLoginTimeout(timeout);
+ logger.config("loginTimeout = " + timeout);
+ }
+ }
+
+ /**
+ * Return <code>true</code> if the conditions record matches the current tag
+ *
+ * @param record the conditions record
+ * @return <code>true</code> if conditions record matches the currently used tag
+ */
+ private boolean matchesTag(final ConditionsRecord record) {
+ if (this.tag == null) {
+ // If there is no tag set then all records pass.
+ return true;
+ }
+ final String recordTag = record.getTag();
+ if (recordTag == null) {
+ // If there is a tag set but the record has no tag, it is rejected.
+ return false;
+ }
+ return this.tag.equals(recordTag);
+ }
+
+ public <CollectionType extends ConditionsObjectCollection<?>> CollectionType newCollection(
+ final Class<CollectionType> collectionType) {
+ final List<TableMetaData> tableMetaDataList = TableRegistry.getTableRegistry().findByCollectionType(
+ collectionType);
+ if (tableMetaDataList.size() > 1) {
+ throw new RuntimeException("More than one table meta data object returned for type: "
+ + collectionType.getName());
+ }
+ final TableMetaData tableMetaData = tableMetaDataList.get(0);
+ CollectionType collection;
+ try {
+ collection = collectionType.newInstance();
+ } catch (InstantiationException | IllegalAccessException e) {
+ throw new RuntimeException("Error creating new collection.", e);
+ }
+ collection.setTableMetaData(tableMetaData);
+ collection.setConnection(this.getConnection());
+ return collection;
+ }
+
+ public <CollectionType extends ConditionsObjectCollection<?>> CollectionType newCollection(
+ final Class<CollectionType> collectionType, final String tableName) {
+ final TableMetaData tableMetaData = TableRegistry.getTableRegistry().findByTableName(tableName);
+ CollectionType collection;
+ try {
+ collection = collectionType.newInstance();
+ } catch (InstantiationException | IllegalAccessException e) {
+ throw new RuntimeException("Error creating new collection.", e);
+ }
+ collection.setTableMetaData(tableMetaData);
+ collection.setConnection(this.getConnection());
+ return collection;
+ }
+
+ /**
+ * Open the database connection.
+ *
+ * @return <code>true</code> if a connection was opened; <code>false</code> if using an existing connection.
+ */
+ public synchronized boolean openConnection() {
+ boolean openedConnection = false;
+ if (!this.isConnected) {
+ // Do the connection parameters need to be figured out automatically?
+ if (this.connectionParameters == null) {
+ // Setup the default read-only connection, which will choose a SLAC or JLab database.
+ this.connectionParameters = ConnectionParameters.fromResource(DEFAULT_CONNECTION_PROPERTIES_RESOURCE);
+ }
+
+ if (!this.loggedConnectionParameters) {
+ // Print out detailed info to the log on first connection within the job.
+ logger.info("opening connection ... " + '\n' + "connection: "
+ + this.connectionParameters.getConnectionString() + '\n' + "host: "
+ + this.connectionParameters.getHostname() + '\n' + "port: "
+ + this.connectionParameters.getPort() + '\n' + "user: " + this.connectionParameters.getUser()
+ + '\n' + "database: " + this.connectionParameters.getDatabase());
+ this.loggedConnectionParameters = true;
+ }
+
+ // Create the connection using the parameters.
+ this.connection = this.connectionParameters.createConnection();
+ this.isConnected = true;
+ openedConnection = true;
+ }
+ logger.fine("connection opened successfully");
+
+ // Flag to indicate whether an existing connection was used or not.
+ return openedConnection;
+ }
+
+ /**
+ * Register the conditions converters with the manager.
+ */
+ private void registerConverters() {
+ if (this.svtConverter != null) {
+ // Remove old SVT converter.
+ this.removeConditionsConverter(this.svtConverter);
+ }
+
+ if (this.ecalConverter != null) {
+ // Remove old ECAL converter.
+ this.registerConditionsConverter(this.ecalConverter);
+ }
+
+ // Is configured for TestRun?
+ if (this.isTestRun()) {
+ // Load Test Run specific converters.
+ this.svtConverter = new TestRunSvtConditionsConverter();
+ this.ecalConverter = new TestRunEcalConditionsConverter();
+ logger.config("registering Test Run conditions converters");
+ } else {
+ // Load the default converters.
+ this.svtConverter = new SvtConditionsConverter();
+ this.ecalConverter = new EcalConditionsConverter();
+ logger.config("registering default conditions converters");
+ }
+ this.registerConditionsConverter(this.svtConverter);
+ this.registerConditionsConverter(this.ecalConverter);
+ }
+
+ /**
+ * This method can be used to perform a database SELECT query.
+ *
+ * @param query the SQL query string
+ * @return the <code>ResultSet</code> from the query
+ * @throws RuntimeException if there is a query error
+ */
+ ResultSet selectQuery(final String query) {
+ logger.fine("executing SQL select query ..." + '\n' + query);
+ ResultSet result = null;
+ Statement statement = null;
+ try {
+ statement = this.connection.createStatement();
+ result = statement.executeQuery(query);
+ } catch (final SQLException x) {
+ throw new RuntimeException("Error in query: " + query, x);
+ }
+ return result;
+ }
+
+ /**
+ * Set the connection parameters of the conditions database.
+ *
+ * @param connectionParameters the connection parameters
+ */
+ public void setConnectionParameters(final ConnectionParameters connectionParameters) {
+ this.connectionParameters = connectionParameters;
+ }
+
+ /**
+ * Set the path to a properties file containing connection settings.
+ *
+ * @param file the properties file
+ */
+ public void setConnectionProperties(final File file) {
+ logger.config("setting connection properties file " + file.getPath());
+ if (!file.exists()) {
+ throw new IllegalArgumentException("The connection properties file does not exist: "
+ + this.connectionPropertiesFile.getPath());
+ }
+ this.connectionParameters = ConnectionParameters.fromProperties(file);
+ }
+
+ /**
+ * Set the connection parameters from an embedded resource location.
+ *
+ * @param resource the classpath resource location
+ */
+ public void setConnectionResource(final String resource) {
+ logger.config("setting connection resource " + resource);
+ this.connectionParameters = ConnectionParameters.fromResource(resource);
+ }
+
+ /**
+ * This method handles changes to the detector name and run number. It is called every time an LCSim event is
+ * created, and so it has internal logic to figure out if the conditions system actually needs to be updated.
+ */
+ @Override
+ public synchronized void setDetector(final String detectorName, final int runNumber)
+ throws ConditionsNotFoundException {
+
+ logger.finest("setDetector " + detectorName + " with run number " + runNumber);
+
+ if (detectorName == null) {
+ throw new IllegalArgumentException("The detectorName argument is null.");
+ }
+
+ if (!this.isInitialized || !detectorName.equals(this.getDetector()) || runNumber != this.getRun()) {
+ if (!this.isFrozen) {
+ logger.info("new detector " + detectorName + " and run #" + runNumber);
+ this.initialize(detectorName, runNumber);
+ } else {
+ logger.finest("Conditions changed but will be ignored because manager is frozen.");
+ }
+ }
+ }
+
+ /**
+ * Set the name of the ECAL sub-detector.
+ *
+ * @param ecalName the name of the ECAL subdetector
+ */
+ private void setEcalName(final String ecalName) {
+ if (ecalName == null) {
+ throw new IllegalArgumentException("The ecalName is null");
+ }
+ this.ecalName = ecalName;
+ logger.info("ECAL name set to " + ecalName);
+ }
+
+ /**
+ * Set the log level.
+ *
+ * @param level the new log level
+ */
+ public void setLogLevel(final Level level) {
+ logger.config("setting log level to " + level);
+ logger.setLevel(level);
+ logger.getHandlers()[0].setLevel(level);
+ this.svtSetup.setLogLevel(level);
+ }
+
+ /**
+ * Set the name of the SVT subdetector.
+ *
+ * @param svtName the name of the SVT subdetector
+ */
+ private void setSvtName(final String svtName) {
+ if (svtName == null) {
+ throw new IllegalArgumentException("The svtName is null");
+ }
+ this.svtName = svtName;
+ logger.info("SVT name set to " + this.ecalName);
+ }
+
+ /**
+ * Set a tag used to filter the accessible conditions records
+ *
+ * @param tag the tag value used to filter returned conditions records
+ */
+ public void setTag(final String tag) {
+ this.tag = tag;
+ logger.info("using conditions tag: " + tag);
+ }
+
+ /**
+ * Setup the database connection from a file specified by a Java system property setting. This could be overridden
+ * by subsequent API calls to {@link #setConnectionProperties(File)} or {@link #setConnectionResource(String)}.
+ */
+ private void setupConnectionFromSystemProperty() {
+ final String systemPropertiesConnectionPath = (String) System.getProperties().get(CONNECTION_PROPERTY);
+ if (systemPropertiesConnectionPath != null) {
+ final File f = new File(systemPropertiesConnectionPath);
+ if (!f.exists()) {
+ throw new RuntimeException("Connection properties file from " + CONNECTION_PROPERTY
+ + " does not exist.");
+ }
+ this.setConnectionProperties(f);
+ logger.info("connection setup from system property " + CONNECTION_PROPERTY + " = "
+ + systemPropertiesConnectionPath);
+ }
+ }
+
+ /**
+ * Configure some properties of this object from an XML file
+ *
+ * @param file the XML file
+ */
+ public void setXmlConfig(final File file) {
+ logger.config("setting XML config from file " + file.getPath());
+ if (!file.exists()) {
+ throw new IllegalArgumentException("The config file does not exist: " + file.getPath());
+ }
+ try {
+ this.configure(new FileInputStream(file));
+ } catch (final FileNotFoundException e) {
+ throw new RuntimeException(e);
+ }
+ }
+
+ /**
+ * Configure this object from an embedded XML resource.
+ *
+ * @param resource the embedded XML resource
+ */
+ public void setXmlConfig(final String resource) {
+ logger.config("setting XML config from resource " + resource);
+ final InputStream is = this.getClass().getResourceAsStream(resource);
+ this.configure(is);
+ }
+
+ /**
+ * Un-freeze the conditions system so that updates will be received again.
+ */
+ public synchronized void unfreeze() {
+ this.isFrozen = false;
+ logger.info("conditions system unfrozen");
+ }
+
+ /**
+ * Perform a SQL query with an update command like INSERT, DELETE or UPDATE.
+ *
+ * @param query the SQL query string
+ * @return the keys of the rows affected
+ */
+ public List<Integer> updateQuery(final String query) {
+ final boolean openedConnection = this.openConnection();
+ logger.fine("executing SQL update query ..." + '\n' + query);
+ final List<Integer> keys = new ArrayList<Integer>();
+ Statement statement = null;
+ ResultSet resultSet = null;
+ try {
+ statement = this.connection.createStatement();
+ statement.executeUpdate(query, Statement.RETURN_GENERATED_KEYS);
+ resultSet = statement.getGeneratedKeys();
+ while (resultSet.next()) {
+ final int key = resultSet.getInt(1);
+ keys.add(key);
+ }
+ } catch (final SQLException x) {
+ throw new RuntimeException("Error in SQL query: " + query, x);
+ }
+ DatabaseUtilities.cleanup(resultSet);
+ this.closeConnection(openedConnection);
+ return keys;
+ }
+}
Added: java/sandbox/HPSTrack.java
=============================================================================
--- java/sandbox/HPSTrack.java (added)
+++ java/sandbox/HPSTrack.java Wed Jan 6 08:21:16 2016
@@ -0,0 +1,393 @@
+package org.hps.recon.tracking;
+
+import static org.lcsim.constants.Constants.fieldConversion;
+import hep.physics.matrix.SymmetricMatrix;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+
+import org.hps.util.Pair;
+import org.lcsim.event.MCParticle;
+import org.lcsim.fit.helicaltrack.HelicalTrackFit;
+import org.lcsim.fit.helicaltrack.HelicalTrackHit;
+import org.lcsim.fit.helicaltrack.HelixUtils;
+import org.lcsim.fit.helicaltrack.MultipleScatter;
+import org.lcsim.geometry.FieldMap;
+import org.lcsim.spacegeom.CartesianVector;
+import org.lcsim.spacegeom.SpacePoint;
+import org.lcsim.spacegeom.SpaceVector;
+import org.lcsim.util.swim.Helix;
+import org.lcsim.util.swim.Line;
+import org.lcsim.util.swim.Trajectory;
+
+/**
+ * Extension of HelicalTrackFit to include HPS-specific information and utilities.
+ *
+ * @author mgraham <[log in to unmask]>
+ * @author phansson <[log in to unmask]>
+ *
+ */
+
+public class HPSTrack extends HelicalTrackFit {
+
+ private BeamSpot _beam;
+ // all of the variables defined below are in the jlab (detector) frame
+ // this position & momentum are measured at the DOCA to the Y-axis,
+ // which is where the tracking returns it's parameters by default
+ private Hep3Vector _pDocaY;
+ private Hep3Vector _posDocaY;
+ // the position & momentum of the track at the intersection of the target (z=0)
+ private Hep3Vector _pTarget;
+ private Hep3Vector _posTarget;
+ // the position & momentum of the track at DOCA to the beam axis (z)
+ private Hep3Vector _pDocaZ;
+ private Hep3Vector _posDocaZ;
+ private double bField = 0.491; // make this set-able
+ private boolean _debug = false;
+ private boolean _debugForward = false;
+ private Trajectory _trajectory;
+ Map<Pair<Integer, Integer>, Double> _fieldMap;
+ Map<Pair<Integer, Integer>, Pair<Double, Double>> _fieldBins;
+ private MCParticle mcParticle;
+
+ public HPSTrack(double[] pars, SymmetricMatrix cov, double[] chisq, int[] ndf, Map<HelicalTrackHit, Double> smap, Map<HelicalTrackHit, MultipleScatter> msmap, BeamSpot beam) {
+ super(pars, cov, chisq, ndf, smap, msmap);
+ _beam = beam;
+ calculateParametersAtTarget();
+ calculateParametersAtDocaY();
+ calculateParametersAtDocaZ();
+ }
+
+ public HPSTrack(double[] pars, SymmetricMatrix cov, double[] chisq, int[] ndf, Map<HelicalTrackHit, Double> smap, Map<HelicalTrackHit, MultipleScatter> msmap) {
+ super(pars, cov, chisq, ndf, smap, msmap);
+ calculateParametersAtTarget();
+ calculateParametersAtDocaY();
+ calculateParametersAtDocaZ();
+ }
+
+ public HPSTrack(HelicalTrackFit htf, BeamSpot beam) {
+ super(htf.parameters(), htf.covariance(), htf.chisq(), htf.ndf(), htf.PathMap(), htf.ScatterMap());
+ _beam = beam;
+ calculateParametersAtTarget();
+ calculateParametersAtDocaY();
+ calculateParametersAtDocaZ();
+ }
+
+ public HPSTrack(HelicalTrackFit htf) {
+ super(htf.parameters(), htf.covariance(), htf.chisq(), htf.ndf(), htf.PathMap(), htf.ScatterMap());
+ calculateParametersAtTarget();
+ calculateParametersAtDocaY();
+ calculateParametersAtDocaZ();
+ }
+
+ public HPSTrack(double[] parameters, SymmetricMatrix covariance, double[] chiSquared, int[] ndf, Map<HelicalTrackHit, Double> sMap, Map<HelicalTrackHit, MultipleScatter> msMap, MCParticle mcParticle) {
+ super(parameters, covariance, chiSquared, ndf, sMap, msMap);
+
+ // Set the MC particle associated with this fit
+ this.mcParticle = mcParticle;
+ }
+
+ /**
+ * Get map of the the track trajectory within the uniform bfield
+ */
+ public Map<Integer, Double[]> trackTrajectory(double zStart, double zStop, int nSteps) {
+ Map<Integer, Double[]> traj = new HashMap<Integer, Double[]>();
+ double step = (zStop - zStart) / nSteps;
+ Double zVal = zStart;
+ Integer nstep = 0;
+ while (zVal < zStop) {
+ Double[] xyz = {0.0, 0.0, 0.0};
+ double s = HelixUtils.PathToXPlane(this, zVal, 1000.0, 1).get(0);
+ xyz[0] = HelixUtils.PointOnHelix(this, s).y();
+ xyz[1] = HelixUtils.PointOnHelix(this, s).z();
+ xyz[2] = zVal;
+ traj.put(nstep, xyz);
+ zVal += step;
+ nstep++;
+ }
+ return traj;
+ }
+
+ /**
+ * Get map of the the track direction within the uniform bfield
+ */
+ public Map<Integer, Double[]> trackDirection(double zStart, double zStop, int nSteps) {
+ Map<Integer, Double[]> traj = new HashMap<Integer, Double[]>();
+ double step = (zStop - zStart) / nSteps;
+ Double zVal = zStart;
+
+ Integer nstep = 0;
+ while (zVal < zStop) {
+ Double[] xyz = {0.0, 0.0, 0.0};
+ double s = HelixUtils.PathToXPlane(this, zVal, 1000.0, 1).get(0);
+ xyz[0] = HelixUtils.Direction(this, s).y();
+ xyz[1] = HelixUtils.Direction(this, s).z();
+ xyz[2] = zVal;
+ traj.put(nstep, xyz);
+ zVal += step;
+ nstep++;
+ }
+ return traj;
+ }
+
+ private void calculateParametersAtTarget() {
+ double pTot = this.p(bField);
+ // currently, PathToXPlane only returns a single path (no loopers!)
+ List<Double> paths = HelixUtils.PathToXPlane(this, 0.0, 100.0, 1);
+ Hep3Vector posTargetTrkSystem = HelixUtils.PointOnHelix(this, paths.get(0));
+ Hep3Vector dirTargetTrkSystem = HelixUtils.Direction(this, paths.get(0));
+ _posTarget = CoordinateTransformations.transformVectorToDetector(posTargetTrkSystem);
+ _pTarget = VecOp.mult(pTot, CoordinateTransformations.transformVectorToDetector(dirTargetTrkSystem));
+
+ }
+
+ private void calculateParametersAtDocaY() {
+ double pTot = this.p(bField);
+ Hep3Vector posDocaYTrkSystem = HelixUtils.PointOnHelix(this, 0);
+ Hep3Vector dirDocaYTrkSystem = HelixUtils.Direction(this, 0);
+ _posDocaY = CoordinateTransformations.transformVectorToDetector(posDocaYTrkSystem);
+ _pDocaY = VecOp.mult(pTot, CoordinateTransformations.transformVectorToDetector(dirDocaYTrkSystem));
+
+ }
+
+ private void calculateParametersAtDocaZ() {
+ double pTot = this.p(bField);
+ double sAtDocaZ = findPathToDocaZ();
+ Hep3Vector posDocaZTrkSystem = HelixUtils.PointOnHelix(this, sAtDocaZ);
+ Hep3Vector dirDocaZTrkSystem = HelixUtils.Direction(this, sAtDocaZ);
+ _posDocaZ = CoordinateTransformations.transformVectorToDetector(posDocaZTrkSystem);
+ _pDocaZ = VecOp.mult(pTot, CoordinateTransformations.transformVectorToDetector(dirDocaZTrkSystem));
+ }
+
+ public Hep3Vector[] getPositionAtZMap(double start, double xFinal, double step) {
+ System.out.println("Don't use this anymore! Use the contructor with FieldMap provided");
+ Hep3Vector out= new BasicHep3Vector(666,666,666);
+ Hep3Vector[] ret={out};
+ return ret;
+ }
+
+
+ public Hep3Vector[] getPositionAtZMap(double start, double xFinal, double step, boolean debugOK) {
+ System.out.println("Don't use this anymore! Use the contructor with FieldMap provided");
+ Hep3Vector out= new BasicHep3Vector(666,666,666);
+ Hep3Vector[] ret={out};
+ return ret;
+ }
+
+ /**
+ * Get the position and direction on the track using B-field map for
+ * extrapolation
+ *
+ * @param start = starting z-position of extrapolation
+ * @param zFinal = final z-position
+ * @param step = step size
+ * @param bmap = the B-Field map
+ * @return position[0] and direction[1] at Z=zfinal
+ */
+ public Hep3Vector[] getPositionAtZMap(double start, double xFinal, double step, FieldMap bmap) {
+ return this.getPositionAtZMap(start, xFinal, step, bmap, false);
+ }
+
+ public Hep3Vector[] getPositionAtZMap(double start, double xFinal, double step, FieldMap bmap, boolean debugOk) {
+
+ double startFringe = start;
+
+ _debugForward = false;
+ if (xFinal > 900)
+ _debugForward = debugOk ? true : false;
+ // if looking upstream, we'll be propagating backwards
+ if (xFinal < 0)
+ step = -step;
+ if (_debugForward)
+ System.out.println(this.toString());
+
+ double sStartFringe = HelixUtils.PathToXPlane(this, startFringe, 1000.0, 1).get(0);
+ if (_debugForward)
+ System.out.println("path to end of fringe = " + sStartFringe + "; xFinal = " + xFinal);
+ double xtmp = startFringe;
+ double ytmp = HelixUtils.PointOnHelix(this, sStartFringe).y();
+ double ztmp = HelixUtils.PointOnHelix(this, sStartFringe).z();
+ double Rorig = this.R();
+ double xCtmp = this.xc();
+ double yCtmp = this.yc();
+ double q = Math.signum(this.curvature());
+ double phitmp = getPhi(xtmp, ytmp, xCtmp, yCtmp, q);
+ if (_debugForward) {
+ System.out.println("\nOriginal xtmp = " + xtmp + "; ytmp = " + ytmp + "; ztmp = " + ztmp + "; phitmp = " + phitmp);
+
+ System.out.println("nOriginal Rorig = " + Rorig + "; xCtmp = " + xCtmp + "; yCtmp = " + yCtmp);
+ }
+ if (_debugForward)
+ System.out.println("Original Direction at Fringe: " + HelixUtils.Direction(this, startFringe).toString());
+ double Rtmp = Rorig;
+ // now start stepping through the fringe field
+ Hep3Vector r0Tmp = HelixUtils.PointOnHelix(this, sStartFringe);
+ SpacePoint r0 = new SpacePoint(r0Tmp);
+ double pTot = this.p(bField);
+ Hep3Vector dirOrig = HelixUtils.Direction(this, sStartFringe);
+ Hep3Vector p0 = VecOp.mult(pTot, dirOrig);
+ Hep3Vector dirTmp = dirOrig;
+ SpacePoint rTmp = r0;
+ Hep3Vector pTmp = p0;
+ double pXOrig = p0.x();
+ double pXTmp = pXOrig;
+ double totalS = sStartFringe;
+ if (_debugForward) {
+ double tmpdX = xFinal - startFringe;
+ Hep3Vector fooExt = extrapolateStraight(dirOrig, tmpdX);
+ Hep3Vector fooFinal = VecOp.add(fooExt, r0Tmp);
+ System.out.println("Extrpolating straight back from startFringe = (" + fooFinal.x() + "," + fooFinal.y() + "," + fooFinal.z() + ")");
+
+ }
+ // follow trajectory while: in fringe field, before end point, and we don't have a looper
+ while (Math.signum(step) * xtmp < Math.signum(step) * xFinal && Math.signum(pXOrig * pXTmp) > 0) {
+ if (_debugForward) {
+ System.out.println("New step in Fringe Field");
+ System.out.println("rTmp = " + rTmp.toString());
+ System.out.println("pTmp = " + pTmp.toString());
+ System.out.println("OriginalHelix pos = " + HelixUtils.PointOnHelix(this, totalS));
+ System.out.println("OriginalHelix Momentum = " + VecOp.mult(pTot, HelixUtils.Direction(this, totalS)));
+ }
+
+ double myBField = bmap.getField(new BasicHep3Vector(rTmp.z(), rTmp.z(), rTmp.x())).z(); // note weird co ordinates for track vs det
+ if (_debugForward)
+ System.out.println("rTmp.x() = " + rTmp.x() + " field = " + myBField);
+ setTrack(pTmp, rTmp, q, myBField);
+ rTmp = _trajectory.getPointAtDistance(step);
+ pTmp = VecOp.mult(pTot, _trajectory.getUnitTangentAtLength(step));
+ pXTmp = pTmp.x();
+ xtmp = rTmp.x();
+ if (_debugForward) {
+ System.out.println("############## done... #############");
+
+ System.out.println("\n");
+ }
+ totalS += step;
+ }
+
+ // Go with finer granularity in the end
+ rTmp = _trajectory.getPointAtDistance(0);
+ xtmp = rTmp.x();
+ pTmp = VecOp.mult(pTot, _trajectory.getUnitTangentAtLength(step));
+ pXTmp = pTmp.x();
+ step = step / 10.0;
+
+ while (Math.signum(step) * xtmp < Math.signum(step) * xFinal && Math.signum(pXOrig * pXTmp) > 0) {
+ if (_debugForward) {
+ System.out.println("New step in Fringe Field");
+ System.out.println("rTmp = " + rTmp.toString());
+ System.out.println("pTmp = " + pTmp.toString());
+ System.out.println("OriginalHelix pos = " + HelixUtils.PointOnHelix(this, totalS));
+ System.out.println("OriginalHelix Momentum = " + VecOp.mult(pTot, HelixUtils.Direction(this, totalS)));
+ }
+
+ double myBField = bmap.getField(new BasicHep3Vector(rTmp.z(), rTmp.z(), rTmp.x())).z(); // note weird co ordinates for track vs det
+
+ if (_debugForward)
+ System.out.println("rTmp.x() = " + rTmp.x() + " field = " + myBField);
+ setTrack(pTmp, rTmp, q, myBField);
+ rTmp = _trajectory.getPointAtDistance(step);
+ pTmp = VecOp.mult(pTot, _trajectory.getUnitTangentAtLength(step));
+ pXTmp = pTmp.x();
+ xtmp = rTmp.x();
+ if (_debugForward) {
+ System.out.println("############## done... #############");
+ System.out.println("\n");
+ }
+ totalS += step;
+ }
+
+ // ok, done with field.
+ Hep3Vector pointInTrking = new BasicHep3Vector(rTmp.x(), rTmp.y(), rTmp.z());
+ if (_debugForward)
+ System.out.println("Position xfinal (tracking) : x = " + xFinal + "; y = " + pointInTrking.y() + "; z = " + pointInTrking.z());
+ Hep3Vector[] out = {CoordinateTransformations.transformVectorToDetector(pointInTrking), CoordinateTransformations.transformVectorToDetector(pTmp)};
+
+ return out;
+ }
+
+ private double getPhi(double x, double y, double xc, double yc, double sign) {
+ // System.out.println("Math.atan2(y - yc, x - xc)="+Math.atan2(y - yc, x - xc));
+ return Math.atan2(y - yc, x - xc) - sign * Math.PI / 2;
+ }
+
+ private Hep3Vector extrapolateStraight(Hep3Vector dir, double deltaX) {
+ double deltaY = deltaX * dir.y();
+ double deltaZ = deltaX * dir.z();
+ return new BasicHep3Vector(deltaX, deltaY, deltaZ);
+ }
+
+ // field that changes linearly from Bmax->0
+ private double getFringe(double x, double halfWidth) {
+ // System.out.println("x = " + x + "; halfWidth = " + halfWidth);
+ if (x / halfWidth > 1)
+ return 1;
+ if (x / halfWidth < -1)
+ return 0;
+
+ return (1.0 / 2.0) * (1 + x / halfWidth);
+ }
+
+ private Hep3Vector getDirection(double phi, double sign) {
+ double ux = Math.cos(phi) * this.sth();
+ double uy = Math.sin(phi) * this.sth();
+ double uz = this.cth();
+ // Return the direction unit vector
+ return new BasicHep3Vector(ux, uy, uz);
+ }
+
+ private double findPathToDocaZ() {
+ double step = 0.1;// 100 um step size
+ double maxS = 100.0; // go to 10cm
+ double s = -30;
+ double minDist = 999999;
+ double minS = s;
+ double dist = 999998;
+ // once the distance starts increasing, quit and return
+ while (dist < minDist) {
+ Hep3Vector pos = HelixUtils.PointOnHelix(this, s);
+ dist = pos.y() * pos.y() + pos.z() * pos.z();
+ s += step;
+ }
+ return minS;
+ }
+
+ private void setTrack(Hep3Vector p0, SpacePoint r0, double q, double B) {
+ SpaceVector p = new CartesianVector(p0.v());
+ double phi = Math.atan2(p.y(), p.x());
+ double lambda = Math.atan2(p.z(), p.rxy());
+ double field = B * fieldConversion;
+
+ if (q != 0 && field != 0) {
+ double radius = p.rxy() / (q * field);
+ _trajectory = new Helix(r0, radius, phi, lambda);
+ } else
+ _trajectory = new Line(r0, phi, lambda);
+ }
+
+ public Trajectory getTrajectory() {
+ return this._trajectory;
+ }
+
+ /**
+ * Get the MC Particle associated with the HelicalTrackFit
+ *
+ * @return mcParticle :
+ */
+ public MCParticle getMCParticle() {
+ return this.mcParticle;
+ }
+
+ /**
+ * Set the MC Particle associated with the HelicalTrackFit
+ *
+ * @param mcParticle :
+ */
+ public void setMCParticle(MCParticle mcParticle) {
+ this.mcParticle = mcParticle;
+ }
+}
Added: java/sandbox/MultipleScattering.java
=============================================================================
--- java/sandbox/MultipleScattering.java (added)
+++ java/sandbox/MultipleScattering.java Wed Jan 6 08:21:16 2016
@@ -0,0 +1,471 @@
+package org.hps.recon.tracking;
+
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+
+import java.util.ArrayList;
+import java.util.Collections;
+import java.util.List;
+
+import org.hps.recon.tracking.MaterialSupervisor.ScatteringDetectorVolume;
+import org.hps.recon.tracking.MaterialSupervisor.SiStripPlane;
+import org.lcsim.detector.IDetectorElement;
+import org.lcsim.detector.solids.Inside;
+import org.lcsim.fit.helicaltrack.HelicalTrackFit;
+import org.lcsim.fit.helicaltrack.HelixUtils;
+import org.lcsim.recon.tracking.seedtracker.ScatterAngle;
+
+/**
+ * Extention of lcsim class to allow use of local classes. Finds scatter points
+ * and magnitude from detector geometry directly.
+ *
+ * @author Per Hansson <[log in to unmask]>
+ */
+public class MultipleScattering extends org.lcsim.recon.tracking.seedtracker.MultipleScattering {
+
+ private boolean _fixTrackMomentum = false;
+ private double _momentum = -99;//dummy
+ private static final double inside_tolerance = 1.0;//tolerance for first (approximate) test of track intersection with sensor
+
+ public MultipleScattering(MaterialManager materialmanager) {
+ super(materialmanager);
+ }
+
+ /**
+ * Override lcsim version and select material manager depending on object
+ * type. This allows to use a local extension of the material manager in the
+ * lcsim track fitting code.
+ *
+ * @param helix
+ * @return a list of ScatterAngle.
+ */
+ @Override
+ public List<ScatterAngle> FindScatters(HelicalTrackFit helix) {
+ if (_debug) {
+ System.out.printf("\n%s: FindScatters() for helix:\n%s\n", this.getClass().getSimpleName(), helix.toString());
+ }
+
+ if (MaterialSupervisor.class.isInstance(this._materialmanager)) {
+ if (_debug) {
+ System.out.printf("%s: use HPS scattering model", this.getClass().getSimpleName());
+ }
+ return this.FindHPSScatters(helix);
+ } else {
+ if (_debug) {
+ System.out.printf("%s: use default lcsim material manager to find scatters\n", this.getClass().getSimpleName());
+ }
+ return super.FindScatters(helix);
+ }
+ }
+
+ /**
+ * Extra interface to keep a function returning the same type as the lcsim
+ * version
+ *
+ * @param helix
+ * @return a list of ScatterAngle.
+ */
+ private List<ScatterAngle> FindHPSScatters(HelicalTrackFit helix) {
+ ScatterPoints scatterPoints = this.FindHPSScatterPoints(helix);
+ return scatterPoints.getScatterAngleList();
+ }
+
+ /**
+ * Find scatter points along helix using the local material manager
+ *
+ * @param helix
+ * @return the points of scatter along the helix
+ */
+ public ScatterPoints FindHPSScatterPoints(HelicalTrackFit helix) {
+ if (_debug) {
+ System.out.printf("\n%s: FindHPSScatters() for helix:\n%s\n", this.getClass().getSimpleName(), helix.toString());
+ System.out.printf("%s: momentum is p=%f,R=%f,B=%f \n", this.getClass().getSimpleName(), helix.p(Math.abs(_bfield)), helix.R(), _bfield);
+ }
+// MG TURN THIS OFF SO IT DOESN'T ABORT STRAIGHT TRACKS
+ // Check that B Field is set
+ if (_bfield == 0. && !_fixTrackMomentum) {
+ throw new RuntimeException("B Field or fixed momentum must be set before calling FindScatters method");
+ }
+
+ // Create a new list to contain the mutliple scatters
+ // List<ScatterAngle> scatters = new ArrayList<ScatterAngle>();
+ ScatterPoints scatters = new ScatterPoints();
+
+ MaterialSupervisor materialSupervisor = (MaterialSupervisor) this._materialmanager;
+
+ List<ScatteringDetectorVolume> materialVols = materialSupervisor.getMaterialVolumes();
+
+ if (_debug) {
+ System.out.printf("%s: there are %d detector volumes in the model\n", this.getClass().getSimpleName(), materialVols.size());
+ }
+
+ for (ScatteringDetectorVolume vol : materialVols) {
+
+ if (_debug) {
+ System.out.printf("\n%s: found detector volume \"%s\"\n", this.getClass().getSimpleName(), vol.getName());
+ }
+
+ // find intersection pathpoint with helix
+ Hep3Vector pos = getHelixIntersection(helix, vol);
+
+ if (pos != null) {
+
+ if (_debug) {
+ System.out.printf("%s: intersection position %s\n", this.getClass().getSimpleName(), pos.toString());
+ }
+
+ // find the track direction at the plane
+ double s = HelixUtils.PathToXPlane(helix, pos.x(), 0., 0).get(0);
+
+ if (_debug) {
+ System.out.printf("%s: path length %f\n", this.getClass().getSimpleName(), s);
+ }
+
+ Hep3Vector dir = HelixUtils.Direction(helix, s);
+
+ if (_debug) {
+ System.out.printf("%s: track dir %s\n", this.getClass().getSimpleName(), dir.toString());
+ }
+
+ // Calculate the material the track will traverse
+ double radlen = vol.getMaterialTraversedInRL(dir);
+
+ if (_debug) {
+ System.out.printf("%s: material traversed: %f R.L. (%fmm) \n", this.getClass().getSimpleName(), radlen, vol.getMaterialTraversed(dir));
+ }
+
+ double p;
+ if (_fixTrackMomentum) {
+ p = _momentum;
+ } else {
+ p = helix.p(this._bfield);
+ }
+ double msangle = this.msangle(p, radlen);
+
+ ScatterAngle scat = new ScatterAngle(s, msangle);
+
+ if (_debug) {
+ System.out.printf("%s: scatter angle %f rad for p %f GeV at path length %f\n", this.getClass().getSimpleName(), scat.Angle(), p, scat.PathLen());
+ }
+
+ ScatterPoint scatterPoint = new ScatterPoint(vol.getDetectorElement(), scat);
+ scatters.addPoint(scatterPoint);
+
+ } else if (_debug) {
+ System.out.printf("\n%s: helix did not intersect this volume \n", this.getClass().getSimpleName());
+ }
+
+ }
+
+ // Sort the multiple scatters by their path length
+ Collections.sort(scatters._points);
+
+ if (_debug) {
+ System.out.printf("\n%s: found %d scatters for this helix:\n", this.getClass().getSimpleName(), scatters.getPoints().size());
+ System.out.printf("%s: %10s %10s\n", this.getClass().getSimpleName(), "s (mm)", "theta(rad)");
+ for (ScatterPoint p : scatters.getPoints()) {
+ System.out.printf("%s: %10.2f %10f\n", this.getClass().getSimpleName(), p.getScatterAngle().PathLen(), p.getScatterAngle().Angle());
+ }
+ }
+ return scatters;
+ }
+
+ public Hep3Vector getHelixIntersection(HelicalTrackFit helix, ScatteringDetectorVolume plane) {
+
+ if (SiStripPlane.class.isInstance(plane)) {
+ return getHelixIntersection(helix, (SiStripPlane) plane);
+ } else {
+ throw new UnsupportedOperationException("This det volume type is not supported yet.");
+ }
+ }
+
+ /*
+ * Returns interception between helix and plane Uses the origin x posiution of the plane and
+ * extrapolates linearly to find teh intersection If inside use an iterative "exact" way to
+ * determine the final position
+ */
+ public Hep3Vector getHelixIntersection(HelicalTrackFit helix, SiStripPlane plane) {
+
+ if (_debug) {
+ System.out.printf("%s: calculate simple helix intercept\n", this.getClass().getSimpleName());
+ System.out.printf("%s: StripSensorPlane:\n", this.getClass().getSimpleName());
+ plane.print();
+ }
+
+ double s_origin = HelixUtils.PathToXPlane(helix, plane.origin().x(), 0., 0).get(0);
+
+ if (Double.isNaN(s_origin)) {
+ if (_debug) {
+ System.out.printf("%s: could not extrapolate to XPlane, too large curvature: origin is at %s \n", this.getClass().getSimpleName(), plane.origin().toString());
+ }
+ return null;
+ }
+
+ Hep3Vector pos = HelixUtils.PointOnHelix(helix, s_origin);
+ Hep3Vector direction = HelixUtils.Direction(helix, s_origin);
+
+ if (_debug) {
+ System.out.printf("%s: position at x=origin is %s with path length %f and direction %s\n", this.getClass().getSimpleName(), pos.toString(), s_origin, direction.toString());
+ }
+
+ // Use this approximate position to get a first estimate if the helix intercepted the plane
+ // This is only because the real intercept position is an iterative procedure and we'd
+ // like to avoid using it if possible
+ // Consider the plane as pure x-plane i.e. no rotations
+ // -> this is not very general, as it assumes that strips are (mostly) along y -> FIX
+ // THIS!?
+ // Transformation from tracking to detector frame
+ Hep3Vector pos_det = VecOp.mult(VecOp.inverse(CoordinateTransformations.getMatrix()), pos);
+ Hep3Vector direction_det = VecOp.mult(VecOp.inverse(CoordinateTransformations.getMatrix()), direction);
+
+ if (_debug) {
+ System.out.printf("%s: position in det frame %s and direction %s\n", this.getClass().getSimpleName(), pos_det.toString(), direction_det.toString());
+ }
+
+ // Transformation from detector frame to sensor frame
+ Hep3Vector pos_sensor = plane.getSensor().getGeometry().getGlobalToLocal().transformed(pos_det);
+ Hep3Vector direction_sensor = plane.getSensor().getGeometry().getGlobalToLocal().rotated(direction_det);
+
+ if (_debug) {
+ System.out.printf("%s: position in sensor frame %s and direction %s\n", this.getClass().getSimpleName(), pos_sensor.toString(), direction_sensor.toString());
+ }
+
+ // find step in w to cross sensor plane
+ double delta_w = -1.0 * pos_sensor.z() / direction_sensor.z();
+
+ // find the point where it crossed the plane
+ Hep3Vector pos_int = VecOp.add(pos_sensor, VecOp.mult(delta_w, direction_sensor));
+ Hep3Vector pos_int_det = plane.getSensor().getGeometry().getLocalToGlobal().transformed(pos_int);
+ // find the intercept in the tracking frame
+ Hep3Vector pos_int_trk = VecOp.mult(CoordinateTransformations.getMatrix(), pos_int_det);
+
+ if (_debug) {
+ System.out.printf("%s: take step %f to get intercept position in sensor frame %s (det: %s trk: %s)\n", this.getClass().getSimpleName(), delta_w, pos_int, pos_int_det.toString(), pos_int_trk.toString());
+ }
+
+ boolean isInside = true;
+ if (Math.abs(pos_int.x()) > plane.getMeasuredDimension() / 2.0 + inside_tolerance) {
+ if (_debug) {
+ System.out.printf("%s: intercept is outside in u\n", this.getClass().getSimpleName());
+ }
+ isInside = false;
+ }
+
+ if (Math.abs(pos_int.y()) > plane.getUnmeasuredDimension() / 2.0 + inside_tolerance) {
+ if (_debug) {
+ System.out.printf("%s: intercept is outside in v\n", this.getClass().getSimpleName());
+ }
+ isInside = false;
+ }
+
+ // Check if it's inside sensor and module and if it contradicts the manual calculation
+ // For now: trust manual calculation and output warning if it's outside BOTH sensor AND module
+ if (_debug) {
+ // check if it's inside the sensor
+ Inside result_inside = plane.getDetectorElement().getGeometry().getPhysicalVolume().getMotherLogicalVolume().getSolid().inside(pos_int);
+ Inside result_inside_module = plane.getSensor().getGeometry().getDetectorElement().getParent().getGeometry().inside(pos_int_det);
+ System.out.printf("%s: Inside result sensor: %s module: %s\n", this.getClass().getSimpleName(), result_inside.toString(), result_inside_module.toString());
+
+ boolean isInsideSolid = false;
+ if (result_inside.equals(Inside.INSIDE) || result_inside.equals(Inside.SURFACE)) {
+ isInsideSolid = true;
+ }
+
+ boolean isInsideSolidModule = false;
+ if (result_inside_module.equals(Inside.INSIDE) || result_inside_module.equals(Inside.SURFACE)) {
+ isInsideSolidModule = true;
+ }
+
+ if (isInside && !isInsideSolid) {
+ System.out.printf("%s: manual calculation says inside sensor, inside solid says outside -> contradiction \n", this.getClass().getSimpleName());
+ if (isInsideSolidModule) {
+ System.out.printf("%s: this intercept is outside sensor but inside module\n", this.getClass().getSimpleName());
+ } else {
+ System.out.printf("%s: warning: this intercept at %s, in sensor frame %s, (sensor origin at %s ) is outside sensor and module!\n", this.getClass().getSimpleName(), pos_int_trk.toString(), pos_int.toString(), plane.origin().toString());
+ }
+ }
+ }
+
+ if (!isInside) {
+ return null;
+ }
+
+ if (_debug) {
+ System.out.printf("%s: found simple intercept at %s \n", this.getClass().getSimpleName(), pos_int_trk.toString());
+ }
+
+ // TODO Catch special cases where the incidental iteration procedure seems to fail
+ if (Math.abs(helix.R()) < 2000 && Math.abs(helix.dca()) > 10.0) {
+ if (_debug) {
+ System.out.printf("%s: momentum is low (p=%f,R=%f,B=%f) and d0 is big (d0=%f), skip the iterative calculation\n", this.getClass().getSimpleName(), helix.p(Math.abs(_bfield)), helix.R(), _bfield, helix.dca());
+ }
+ return pos_int_trk;
+ }
+
+ if (_debug) {
+ System.out.printf("%s: calculate iterative helix intercept\n", this.getClass().getSimpleName());
+ }
+
+ //Hep3Vector pos_iter_trk = TrackUtils.getHelixPlaneIntercept(helix, plane.normal(), plane.origin(), _bfield, s_origin);
+ Hep3Vector pos_iter_trk = TrackUtils.getHelixPlaneIntercept(helix, plane.normal(), plane.origin(), _bfield);
+
+ if (pos_iter_trk == null) {
+ System.out.printf("%s: iterative intercept failed for helix \n%s\n at sensor with org=%s, unit w=%s\n", this.getClass().getSimpleName(), helix.toString(), plane.origin().toString(), plane.normal().toString());
+ System.out.printf("%s: => use simple intercept pos=%s\n", this.getClass().getSimpleName(), pos_int_trk);
+ System.out.printf("helix pos=%s dir=%s\n", pos, direction);
+ return pos_int_trk;
+ }
+
+ if (_debug) {
+// if (VecOp.sub(pos_iter_trk, pos_int_trk).magnitude()>1e-4)
+ System.out.printf("%s: iterative helix intercept point at %s (diff to approx: %s) \n", this.getClass().getSimpleName(), pos_iter_trk.toString(), VecOp.sub(pos_iter_trk, pos_int_trk).toString());
+ }
+
+ // find position in sensor frame
+ Hep3Vector pos_iter_sensor = plane.getSensor().getGeometry().getGlobalToLocal().transformed(VecOp.mult(VecOp.inverse(CoordinateTransformations.getMatrix()), pos_iter_trk));
+
+ if (_debug) {
+ System.out.printf("%s: found iterative helix intercept in sensor coordinates at %s\n", this.getClass().getSimpleName(), pos_iter_sensor.toString());
+ }
+
+ isInside = true;
+ if (Math.abs(pos_iter_sensor.x()) > plane.getMeasuredDimension() / 2.0) {
+ if (this._debug) {
+ System.out.printf("%s: intercept is outside in u\n", this.getClass().getSimpleName());
+ }
+ isInside = false;
+ }
+
+ if (Math.abs(pos_iter_sensor.y()) > plane.getUnmeasuredDimension() / 2.0) {
+ if (this._debug) {
+ System.out.printf("%s: intercept is outside in v\n", this.getClass().getSimpleName());
+ }
+ isInside = false;
+ }
+
+ if (_debug) {
+ Hep3Vector pos_iter_det = VecOp.mult(VecOp.inverse(CoordinateTransformations.getMatrix()), pos_iter_trk);
+ Inside result_inside = plane.getDetectorElement().getGeometry().getPhysicalVolume().getMotherLogicalVolume().getSolid().inside(pos_iter_sensor);
+ Inside result_inside_module = plane.getSensor().getGeometry().getDetectorElement().getParent().getGeometry().inside(pos_iter_det);
+ System.out.printf("%s: Inside result sensor: %s module: %s\n", this.getClass().getSimpleName(), result_inside.toString(), result_inside_module.toString());
+
+ boolean isInsideSolid = false;
+ if (result_inside.equals(Inside.INSIDE) || result_inside.equals(Inside.SURFACE)) {
+ isInsideSolid = true;
+ }
+
+ boolean isInsideSolidModule = false;
+ if (result_inside_module.equals(Inside.INSIDE) || result_inside_module.equals(Inside.SURFACE)) {
+ isInsideSolidModule = true;
+ }
+
+ // Check if it's inside sensor and module and if it contradicts the manual calculation
+ // For now: trust manual calculation and output warning if it's outside BOTH sensor AND
+ // module -> FIX THIS!?
+ if (isInside && !isInsideSolid) {
+ System.out.printf("%s: manual iterative calculation says inside sensor, inside solid says outside -> contradiction \n", this.getClass().getSimpleName());
+ if (isInsideSolidModule) {
+ System.out.printf("%s: this iterative intercept is outside sensor but inside module\n", this.getClass().getSimpleName());
+ } else {
+ System.out.printf("%s: warning: this iterative intercept %s, sensor frame %s, (sensor origin %s ) is outside sensor and module!\n", this.getClass().getSimpleName(), pos_iter_trk.toString(), pos_iter_sensor.toString(), plane.origin().toString());
+ }
+ }
+ }
+
+ if (!isInside) {
+ return null;
+ }
+
+ if (_debug) {
+ System.out.printf("%s: found intercept at %s \n", this.getClass().getSimpleName(), pos_iter_trk.toString());
+ }
+
+ return pos_iter_trk;
+ }
+
+ @Override
+ public void setDebug(boolean debug) {
+ _debug = debug;
+ }
+
+ public MaterialManager getMaterialManager() {
+ // Should be safe to cast here
+ return (MaterialManager) _materialmanager;
+ }
+
+ public void fixTrackMomentum(double mom) {
+ _fixTrackMomentum = true;
+ _momentum = mom;
+ }
+
+ /**
+ * Nested class to encapsulate the scatter angles and which detector element
+ * it is related to
+ */
+ public class ScatterPoint implements Comparable<ScatterPoint> {
+
+ IDetectorElement _det;
+ ScatterAngle _scatterAngle;
+
+ public ScatterPoint(IDetectorElement det, ScatterAngle scatterAngle) {
+ _det = det;
+ _scatterAngle = scatterAngle;
+ }
+
+ public IDetectorElement getDet() {
+ return _det;
+ }
+
+ public ScatterAngle getScatterAngle() {
+ return _scatterAngle;
+ }
+
+ @Override
+ public int compareTo(ScatterPoint p) {
+ return p.getScatterAngle().PathLen() > this._scatterAngle.PathLen() ? -1 : 1;
+ }
+ }
+
+ /**
+ * Nested class to encapsulate a list of scatters
+ *
+ */
+ public class ScatterPoints {
+
+ private List<ScatterPoint> _points;
+
+ public ScatterPoints(List<ScatterPoint> _points) {
+ this._points = _points;
+ }
+
+ private ScatterPoints() {
+ _points = new ArrayList<ScatterPoint>();
+ }
+
+ public List<ScatterPoint> getPoints() {
+ return _points;
+ }
+
+ public void addPoint(ScatterPoint point) {
+ _points.add(point);
+ }
+
+ private List<ScatterAngle> getScatterAngleList() {
+ List<ScatterAngle> scatters = new ArrayList<ScatterAngle>();
+ for (ScatterPoint p : _points) {
+ scatters.add(p._scatterAngle);
+ }
+ return scatters;
+ }
+
+ public ScatterPoint getScatterPoint(IDetectorElement detectorElement) {
+ for (ScatterPoint p : _points) {
+ if (p.getDet().equals(detectorElement)) {
+ return p;
+ }
+ }
+ return null;
+ }
+ }
+
+}
Added: java/sandbox/SvtBiasConstant.java
=============================================================================
--- java/sandbox/SvtBiasConstant.java (added)
+++ java/sandbox/SvtBiasConstant.java Wed Jan 6 08:21:16 2016
@@ -0,0 +1,124 @@
+package org.hps.conditions.svt;
+
+import java.sql.Timestamp;
+import java.util.Date;
+
+import org.hps.conditions.api.BaseConditionsObject;
+import org.hps.conditions.api.BaseConditionsObjectCollection;
+import org.hps.conditions.database.Converter;
+import org.hps.conditions.database.Field;
+import org.hps.conditions.database.MultipleCollectionsAction;
+import org.hps.conditions.database.Table;
+
+/**
+ *
+ * Encapsulates an SVT bias constant, which is range in time where bias was ON.
+ *
+ * @author Per Hansson Adrian, SLAC
+ */
+@Table(names = "svt_bias_constants")
+@Converter(multipleCollectionsAction = MultipleCollectionsAction.LAST_CREATED)
+public final class SvtBiasConstant extends BaseConditionsObject {
+
+ /**
+ * The collection implementation for {@link SvtBiasConstant}.
+ */
+ @SuppressWarnings("serial")
+ public static class SvtBiasConstantCollection extends BaseConditionsObjectCollection<SvtBiasConstant> {
+
+ /**
+ * Find bias constant by {@link Timestamp}.
+ *
+ * @param date the {@link Timestamp}
+ * @return the constant containing the date or <code>null</code> otherwise.
+ *
+ */
+ public SvtBiasConstant find(Timestamp date) {
+ for (SvtBiasConstant constant : this) {
+ if(date.getTime() >= constant.getStart() && date.getTime() <= constant.getEnd()) {
+ return constant;
+ }
+ }
+ return null;
+ }
+
+ /**
+ * Find bias constant by {@link Date} type. Converts date to {@link Timestamp} internally.
+ *
+ * @param date the {@link Date}.
+ * @return the constant containing the date or <code>null</code> otherwise.
+ *
+ */
+ public SvtBiasConstant find(Date date) {
+ return find(new Timestamp(date.getTime()));
+ }
+ }
+
+ private Timestamp getTimeStampFieldValue(String field) {
+ Object o = fieldValues.get(field);
+ Timestamp stamp = null;
+ if(o instanceof Timestamp)
+ stamp = (Timestamp)o;
+ else
+ throw new RuntimeException("the value for field " + field + " is not a Timestamp");
+ System.out.println("getTimeStampFieldValue for " + field + " is " + stamp.getTime());
+ return stamp;
+ }
+
+
+
+ /**
+ * The start date as a Unix timestamp in milliseconds (GMT).
+ *
+ * @return the start date as a Unix timestamp
+ */
+ @Field(names = {"start"})
+ public Long getStart() {
+ return getFieldValue("start");
+ }
+
+ /**
+ * The end date as a Unix timestamp in milliseconds (GMT).
+ *
+ * @return the end date as a Unix timestamp
+ */
+ @Field(names = {"end"})
+ public Long getEnd() {
+ return getFieldValue("end");
+ }
+
+ /**
+ * The start {@link Date} based on the {@link Timestamp}.
+ * This takes care of interfacing to the more used {@link Date} object.
+ * It is discouraged to treat {@link Timestamp} as type inheritance.
+ *
+ * @return the start date.
+ */
+ public Date getStart() {
+ return new Date(getStartTimestamp().getTime()*1000);
+ }
+
+
+ /**
+ * The end {@link Date} based on the {@link Timestamp}.
+ * This takes care of interfacing to the more used {@link Date} object.
+ * It is discouraged to treat {@link Timestamp} as type inheritance.
+ *
+ * @return the start date.
+ */
+ public Date getEnd() {
+ return new Date(getEndTimestamp().getTime()*1000);
+ }
+
+
+
+ /**
+ * The bias value.
+ *
+ * @return the bias value
+ */
+ @Field(names = {"value"})
+ public Double getValue() {
+ return getFieldValue("value");
+ }
+}
Added: java/sandbox/SvtOldHeaderAnalysisDriver.java
=============================================================================
--- java/sandbox/SvtOldHeaderAnalysisDriver.java (added)
+++ java/sandbox/SvtOldHeaderAnalysisDriver.java Wed Jan 6 08:21:16 2016
@@ -0,0 +1,217 @@
+/**
+ *
+ */
+package org.hps.users.phansson;
+
+import hep.aida.IHistogram2D;
+import hep.aida.IPlotter;
+
+import java.io.FileWriter;
+import java.io.IOException;
+import java.io.PrintWriter;
+import java.util.Date;
+import java.util.List;
+import java.util.logging.FileHandler;
+import java.util.logging.Level;
+import java.util.logging.Logger;
+
+import org.hps.analysis.trigger.util.TriggerDataUtils;
+import org.hps.evio.SvtEvioReader;
+import org.hps.evio.SvtEvioUtils;
+import org.hps.readout.svt.SvtHeaderDataInfo;
+import org.hps.record.triggerbank.AbstractIntData;
+import org.hps.record.triggerbank.HeadBankData;
+import org.hps.util.BasicLogFormatter;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.GenericObject;
+import org.lcsim.geometry.Detector;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+import org.lcsim.util.log.LogUtil;
+
+/**
+ * @author Per Hansson Adrian <[log in to unmask]>
+ *
+ */
+public class SvtOldHeaderAnalysisDriver extends Driver {
+
+ private final AIDA aida = AIDA.defaultInstance();
+
+ private final String HeaderCollectionName = "SvtHeaders";
+ private final Logger logger = LogUtil.create(SvtOldHeaderAnalysisDriver.class.getSimpleName(), new BasicLogFormatter(), Level.INFO);
+ private int nEventsProcessed = 0;
+ private Date eventDate = new Date(0);
+ private IHistogram2D rceSyncErrorCount;
+ private IHistogram2D rceOFErrorCount;
+ private IHistogram2D rceSkipCount;
+ private IHistogram2D rceMultisampleErrorCount;
+ //private Map< Integer, Map< String, IHistogram1D> > histError = new HashMap< Integer, Map<String, IHistogram1D >>();
+ private int nRceSyncErrorCountN = 0;
+ private int nRceOFErrorCount = 0;
+ private int nRceSkipCount = 0;
+ private int nRceMultisampleErrorCount = 0;
+ private int nRceSvtHeaders = 0;
+ IPlotter plotter;
+
+ private String logFileName = "dummy.log";
+ FileWriter fileWriter;
+ PrintWriter printWriter;
+
+ private boolean showPlots = false;
+ private final String triggerBankCollectionName = "TriggerBank";
+
+ /**
+ *
+ */
+ public SvtOldHeaderAnalysisDriver() {
+
+ }
+
+ public void setLogFileName(String name) {
+ this.logFileName = name;
+ }
+
+ public void setShowPlots(boolean show) {
+ this.showPlots = show;
+ }
+
+ @Override
+ protected void detectorChanged(Detector detector) {
+
+ FileHandler fh;
+ try {
+ fh = new FileHandler(logFileName);
+ logger.addHandler(fh);
+ fh.setFormatter(new BasicLogFormatter());
+ } catch (SecurityException | IOException e1) {
+ e1.printStackTrace();
+ }
+
+ plotter = aida.analysisFactory().createPlotterFactory().create("rceSyncError");
+ plotter.createRegions(2, 2);
+ final int n = SvtEvioReader.MAX_ROC_BANK_TAG+1 - SvtEvioReader.MIN_ROC_BANK_TAG;
+ rceSyncErrorCount = aida.histogram2D("rceSyncError", n, SvtEvioReader.MIN_ROC_BANK_TAG, SvtEvioReader.MAX_ROC_BANK_TAG+1, 2, 0, 1);
+ rceOFErrorCount = aida.histogram2D("rceOFError", n, SvtEvioReader.MIN_ROC_BANK_TAG, SvtEvioReader.MAX_ROC_BANK_TAG+1, 2, 0, 1);
+ rceSkipCount = aida.histogram2D("rceSkipCount", n, SvtEvioReader.MIN_ROC_BANK_TAG, SvtEvioReader.MAX_ROC_BANK_TAG+1, 2, 0, 1);
+ rceMultisampleErrorCount = aida.histogram2D("rceMultisampleError", n, SvtEvioReader.MIN_ROC_BANK_TAG, SvtEvioReader.MAX_ROC_BANK_TAG+1, 2, 0, 1);
+
+ plotter.region(0).plot(rceSyncErrorCount);
+ plotter.region(1).plot(rceOFErrorCount);
+ plotter.region(2).plot(rceSkipCount);
+ plotter.region(3).plot(rceMultisampleErrorCount);
+ if(showPlots ) plotter.show();
+
+
+ }
+
+
+ @Override
+ protected void process(EventHeader event) {
+
+ if(event.hasCollection(GenericObject.class, triggerBankCollectionName)) {
+ Date currentEventDate = TriggerDataUtils.getEventTimeStamp(event, triggerBankCollectionName);
+ if( currentEventDate == null) {
+ logger.info("Couldn't get event date from trigger bank for processed " + nEventsProcessed);
+ // throw new RuntimeException("Couldn't get event date from trigger bank!");
+ } else {
+ eventDate = currentEventDate;
+ }
+ }
+ // log start of run
+ if( nEventsProcessed == 0 )
+ logger.info("startOfRun: run " + event.getRunNumber() + " event " + event.getEventNumber() + " processed " + nEventsProcessed + " date " + eventDate.toString());
+
+ if( !event.hasCollection(GenericObject.class, HeaderCollectionName) )
+ return;
+
+ List<GenericObject> headers = event.get(GenericObject.class, HeaderCollectionName);
+
+ logger.fine("Found " + headers.size() + " SvtHeaders in event " + event.getEventNumber() + " run " + event.getRunNumber());
+
+ for(GenericObject header : headers ) {
+ logger.fine("nint " + header.getNInt());
+ int roc = SvtHeaderDataInfo.getNum(header);
+
+ // find the errors in the header
+ int syncError = SvtEvioUtils.getSvtTailSyncErrorBit(SvtHeaderDataInfo.getTail(header));
+ int oFError = SvtEvioUtils.getSvtTailOFErrorBit(SvtHeaderDataInfo.getTail(header));
+ int skipCount = SvtEvioUtils.getSvtTailMultisampleSkipCount(SvtHeaderDataInfo.getTail(header));
+
+ // check bits
+ checkBitValueRange(oFError);
+ checkBitValueRange(syncError);
+
+ // print header errors to log
+ if( syncError != 0) {
+ logger.info("syncError: run " + event.getRunNumber() + " event " + event.getEventNumber() + " date " + " processed " + nEventsProcessed + " date " + eventDate.toString()
+ + " roc " + roc);
+ }
+ if( oFError != 0) {
+ logger.info("oFError: run " + event.getRunNumber() + " event " + event.getEventNumber() + " date " + " processed " + nEventsProcessed + " date " + eventDate.toString()
+ + " roc " + roc);
+ }
+ for(int i=0; i < skipCount; ++i) {
+ if( oFError != 0) {
+ logger.info("skipCount: run " + event.getRunNumber() + " event " + event.getEventNumber() + " date " + " processed " + nEventsProcessed + " date " + eventDate.toString()
+ + " roc " + roc);
+ }
+ }
+
+
+
+ // Check for multisample tail error bit
+ int nMultisamples = SvtEvioUtils.getSvtTailMultisampleCount(SvtHeaderDataInfo.getTail(header));
+ logger.fine(nMultisamples + " multisamples");
+ int multisampleErrorBits = 0;
+ for(int iMultisample = 0; iMultisample != nMultisamples; ++iMultisample) {
+ int multisampleHeader = SvtHeaderDataInfo.getMultisample(iMultisample, header);
+ logger.fine("found multisample tail: " + Integer.toHexString(multisampleHeader));
+ int multisampleErrorBit = SvtEvioUtils.getErrorBitFromMultisampleHeader(multisampleHeader);
+ checkBitValueRange(multisampleErrorBit);
+ logger.fine("found multisample tail error bit: " + multisampleErrorBit);
+ if( multisampleErrorBit != 0) {
+ multisampleErrorBits++;
+ logger.info("multisample tail error: run " + event.getRunNumber() + " event " + event.getEventNumber() + " date " + eventDate.toString()
+ + " roc " + roc + " feb " + SvtEvioUtils.getFebIDFromMultisampleTail(multisampleHeader)
+ + " hybrid " + SvtEvioUtils.getFebHybridIDFromMultisampleTail(multisampleHeader)
+ + " apv " + SvtEvioUtils.getApvFromMultisampleTail(multisampleHeader));
+ }
+ }
+
+ // keep track how many headers have errors
+ rceSyncErrorCount.fill(roc, syncError);
+ rceOFErrorCount.fill(roc, oFError);
+ rceSkipCount.fill(roc, skipCount);
+ rceMultisampleErrorCount.fill(roc, multisampleErrorBits > 0 ? 1 : 0);
+ if( syncError > 0) nRceSyncErrorCountN++;
+ if( oFError > 0 ) nRceOFErrorCount++;
+ if( skipCount > 0 ) nRceSkipCount++;
+ if( multisampleErrorBits > 0 ) nRceMultisampleErrorCount++;
+ nRceSvtHeaders++;
+
+
+
+ }
+
+
+ nEventsProcessed++;
+ }
+
+ private void checkBitValueRange(int val) {
+ if( val != 0 && val != 1)
+ throw new RuntimeException("invalid value for error bit " + val);
+ }
+
+ @Override
+ protected void endOfData() {
+ logger.info("endOfData: processed " + nEventsProcessed + " date " + eventDate.toString());
+ logger.info("nRceSvtHeaders " + nRceSvtHeaders);
+ logger.info("nRceSyncErrorCountN " + nRceSyncErrorCountN);
+ logger.info("nRceOFErrorCount " + nRceOFErrorCount);
+ logger.info("nRceSkipCount " + nRceSkipCount);
+ logger.info("nRceMultisampleErrorCount " + nRceMultisampleErrorCount);
+
+ }
+
+
+}
Added: java/sandbox/SvtOldHeaderDataInfo.java
=============================================================================
--- java/sandbox/SvtOldHeaderDataInfo.java (added)
+++ java/sandbox/SvtOldHeaderDataInfo.java Wed Jan 6 08:21:16 2016
@@ -0,0 +1,119 @@
+/**
+ *
+ */
+package org.hps.readout.svt;
+
+import org.freehep.graphicsio.swf.SWFAction.GetMember;
+import org.lcsim.event.GenericObject;
+
+/**
+ * @author Per Hansson Adrian <[log in to unmask]>
+ *
+ */
+public class SvtOldHeaderDataInfo implements GenericObject {
+
+ private final int num;
+ private final int header;
+ private final int tail;
+ private int[] multisampleheader;
+
+
+
+ public SvtOldHeaderDataInfo(int num, int header, int tail) {
+ this.num = num;
+ this.header = header;
+ this.tail = tail;
+ }
+
+ public SvtOldHeaderDataInfo(int num, int header, int tail, int[] multisampleheaders) {
+ this.num = num;
+ this.header = header;
+ this.tail = tail;
+ this.multisampleheader = multisampleheaders;
+ }
+
+ public void setMultisampleHeaders(int[] multisampleheaders) {
+ this.multisampleheader = multisampleheaders;
+ }
+
+ public int[] getMultisampleHeaders() {
+ return this.multisampleheader;
+ }
+
+ public int getMultisampleHeader(int index) {
+ if( index >= getMultisampleHeaders().length || index < 0)
+ throw new ArrayIndexOutOfBoundsException(index);
+ return this.multisampleheader[index];
+ }
+
+ public int getNum() {
+ return this.num;
+ }
+
+ public int getTail() {
+ return this.tail;
+ }
+
+ public int getHeader() {
+ return this.header;
+ }
+
+ @Override
+ public int getNInt() {
+ return 3 + multisampleheader.length;
+ }
+
+ @Override
+ public int getNFloat() {
+ return 0;
+ }
+
+ @Override
+ public int getNDouble() {
+ return 0;
+ }
+
+ @Override
+ public int getIntVal(int index) {
+ int value;
+ switch (index) {
+ case 0:
+ value = getNum();
+ break;
+ case 1:
+ value = getHeader();
+ break;
+ case 2:
+ value = getTail();
+ break;
+ default:
+ if( (index-3) >= getMultisampleHeaders().length )
+ throw new RuntimeException("Invalid index " + Integer.toString(index));
+ else
+ value = getMultisampleHeader(index -3);
+ break;
+ }
+ return value;
+ }
+
+
+ @Override
+ public float getFloatVal(int index) {
+ throw new ArrayIndexOutOfBoundsException();
+ }
+
+
+ @Override
+ public double getDoubleVal(int index) {
+ throw new ArrayIndexOutOfBoundsException();
+ }
+
+
+ @Override
+ public boolean isFixedSize() {
+ return true;
+ }
+
+
+
+}
Added: java/sandbox/TrackUtils.java
=============================================================================
--- java/sandbox/TrackUtils.java (added)
+++ java/sandbox/TrackUtils.java Wed Jan 6 08:21:16 2016
@@ -0,0 +1,1197 @@
+package org.hps.recon.tracking;
+
+import hep.physics.matrix.SymmetricMatrix;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Matrix;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.SpacePoint;
+import hep.physics.vec.VecOp;
+import java.util.ArrayList;
+import java.util.HashMap;
+import java.util.HashSet;
+import java.util.List;
+import java.util.Map;
+import java.util.Set;
+import org.hps.recon.tracking.EventQuality.Quality;
+import org.hps.recon.tracking.gbl.HelicalTrackStripGbl;
+import static org.lcsim.constants.Constants.fieldConversion;
+import org.lcsim.detector.ITransform3D;
+import org.lcsim.detector.solids.Box;
+import org.lcsim.detector.solids.Point3D;
+import org.lcsim.detector.solids.Polygon3D;
+import org.lcsim.detector.tracker.silicon.HpsSiSensor;
+import org.lcsim.detector.tracker.silicon.SiSensor;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.LCRelation;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.RawTrackerHit;
+import org.lcsim.event.RelationalTable;
+import org.lcsim.event.Track;
+import org.lcsim.event.TrackState;
+import org.lcsim.event.TrackerHit;
+import org.lcsim.event.base.BaseRelationalTable;
+import org.lcsim.event.base.BaseTrackState;
+import org.lcsim.fit.helicaltrack.HelicalTrackFit;
+import org.lcsim.fit.helicaltrack.HelicalTrackHit;
+import org.lcsim.fit.helicaltrack.HelicalTrackStrip;
+import org.lcsim.fit.helicaltrack.HelixParamCalculator;
+import org.lcsim.fit.helicaltrack.HelixUtils;
+import org.lcsim.fit.helicaltrack.HitUtils;
+import org.lcsim.fit.helicaltrack.MultipleScatter;
+import org.lcsim.geometry.Detector;
+import org.lcsim.geometry.FieldMap;
+import org.lcsim.geometry.subdetector.BarrelEndcapFlag;
+import org.lcsim.recon.tracking.seedtracker.SeedCandidate;
+import org.lcsim.recon.tracking.seedtracker.SeedTrack;
+import org.lcsim.spacegeom.CartesianVector;
+import org.lcsim.spacegeom.SpaceVector;
+import org.lcsim.util.swim.Helix;
+import org.lcsim.util.swim.Line;
+import org.lcsim.util.swim.Trajectory;
+
+/**
+ * Assorted helper functions for the track and helix objects in lcsim. Re-use as
+ * much of HelixUtils as possible.
+ *
+ * @author Omar Moreno <[log in to unmask]>
+ */
+// TODO: Switch to tracking/LCsim coordinates for the extrapolation output!
+// FIXME: This class should probably be broken up into several different sets of utilities by type. --JM
+public class TrackUtils {
+
+ /**
+ * Private constructor to make class static only
+ */
+ private TrackUtils() {
+ }
+
+ /**
+ * Extrapolate track to a position along the x-axis. Turn the track into a
+ * helix object in order to use HelixUtils.
+ *
+ * @param track
+ * @param x
+ * @return
+ */
+ public static Hep3Vector extrapolateHelixToXPlane(Track track, double x) {
+ return extrapolateHelixToXPlane(getHTF(track), x);
+ }
+
+ /**
+ * Extrapolate helix to a position along the x-axis. Re-use HelixUtils.
+ *
+ * @param track
+ * @param x
+ * @return
+ */
+ public static Hep3Vector extrapolateHelixToXPlane(HelicalTrackFit htf, double x) {
+ double s = HelixUtils.PathToXPlane(htf, x, 0., 0).get(0);
+ return HelixUtils.PointOnHelix(htf, s);
+ }
+
+ // ==========================================================================
+ // Helper functions for track parameters and commonly used derived variables
+ public static double getPhi(Track track, Hep3Vector position) {
+ double x = Math.sin(getPhi0(track)) - (1 / getR(track)) * (position.x() - getX0(track));
+ double y = Math.cos(getPhi0(track)) + (1 / getR(track)) * (position.y() - getY0(track));
+ return Math.atan2(x, y);
+ }
+
+ public static double getX0(Track track) {
+ return -1 * getDoca(track) * Math.sin(getPhi0(track));
+ }
+
+ public static double getR(Track track) {
+ return 1.0 / track.getTrackStates().get(0).getOmega();
+ }
+
+ public static double getY0(Track track) {
+ return getDoca(track) * Math.cos(getPhi0(track));
+ }
+
+ public static double getDoca(Track track) {
+ return track.getTrackStates().get(0).getD0();
+ }
+
+ public static double getPhi0(Track track) {
+ return track.getTrackStates().get(0).getPhi();
+ }
+
+ public static double getZ0(Track track) {
+ return track.getTrackStates().get(0).getZ0();
+ }
+
+ public static double getTanLambda(Track track) {
+ return track.getTrackStates().get(0).getTanLambda();
+ }
+
+ public static double getSinTheta(Track track) {
+ return 1 / Math.sqrt(1 + Math.pow(getTanLambda(track), 2));
+ }
+
+ public static double getCosTheta(Track track) {
+ return getTanLambda(track) / Math.sqrt(1 + Math.pow(getTanLambda(track), 2));
+ }
+
+ // ==========================================================================
+ /**
+ * Calculate the point of interception between the helix and a plane in
+ * space. Uses an iterative procedure. This function makes assumptions on
+ * the sign and convecntion of the B-field. Be careful.
+ *
+ * @param helfit - helix
+ * @param unit_vec_normal_to_plane - unit vector normal to the plane
+ * @param point_on_plane - point on the plane
+ * @param bfield - magnetic field value
+ * @return point at intercept
+ */
+ public static Hep3Vector getHelixPlaneIntercept(HelicalTrackFit helfit, Hep3Vector unit_vec_normal_to_plane, Hep3Vector point_on_plane, double bfield) {
+ boolean debug = false;
+ //Hep3Vector B = new BasicHep3Vector(0, 0, -1);
+ //WTrack wtrack = new WTrack(helfit, -1.0*bfield); //
+ Hep3Vector B = new BasicHep3Vector(0, 0, 1);
+ WTrack wtrack = new WTrack(helfit, bfield); //
+ if (debug) {
+ System.out.printf("getHelixPlaneIntercept:find intercept between plane defined by point on plane %s, unit vec %s, bfield %.3f, h=%s and WTrack \n%s \n", point_on_plane.toString(), unit_vec_normal_to_plane.toString(), bfield, B.toString(), wtrack.toString());
+ }
+ Hep3Vector intercept_point = wtrack.getHelixAndPlaneIntercept(point_on_plane, unit_vec_normal_to_plane, B);
+ if (debug) {
+ System.out.printf("getHelixPlaneIntercept: found intercept point at %s\n", intercept_point.toString());
+ }
+ return intercept_point;
+ }
+
+
+ /**
+ * Calculate the point of interception between the helix and a plane in
+ * space. Uses an iterative procedure.
+ *
+ * @param helfit - helix
+ * @param strip - strip cluster that will define the plane
+ * @param bfield - magnetic field value
+ * @return point at intercept
+ */
+ public static Hep3Vector getHelixPlaneIntercept(HelicalTrackFit helfit, HelicalTrackStripGbl strip, double bfield) {
+ Hep3Vector point_on_plane = strip.origin();
+ Hep3Vector unit_vec_normal_to_plane = VecOp.cross(strip.u(), strip.v());// strip.w();
+ Hep3Vector intercept_point = getHelixPlaneIntercept(helfit, unit_vec_normal_to_plane, point_on_plane, bfield);
+ //double s_origin = HelixUtils.PathToXPlane(helfit, strip.origin().x(), 0., 0).get(0);
+ //Hep3Vector intercept_point = getHelixPlaneIntercept(helfit, unit_vec_normal_to_plane, point_on_plane, bfield, s_origin);
+ return intercept_point;
+ }
+
+ /*
+ * Calculates the point on the helix in the x-y plane at the intercept with plane The normal of
+ * the plane is in the same x-y plane as the circle.
+ * @param helix
+ * @param vector normal to plane
+ * @param origin of plane
+ * @return point in the x-y plane of the intercept
+ */
+ public Hep3Vector getHelixXPlaneIntercept(HelicalTrackFit helix, Hep3Vector w, Hep3Vector origin) {
+ throw new RuntimeException("this function is not working properly; don't use it");
+
+ // FInd the intercept point x_int,y_int, between the circle and sensor, which becomes a
+ // line in the x-y plane in this case.
+ // y_int = k*x_int + m
+ // R^2 = (y_int-y_c)^2 + (x_int-x_c)^2
+ // solve for x_int
+ }
+
+ /**
+ * Get position of a track extrapolated to the HARP in the HPS test run 2012
+ *
+ * @param track
+ * @return position at HARP
+ */
+ public static Hep3Vector getTrackPositionAtHarp(Track track) {
+ return extrapolateTrack(track, BeamlineConstants.HARP_POSITION_TESTRUN);
+ }
+
+ /**
+ * Get position of a track extrapolated to the ECAL face in the HPS test run
+ * 2012
+ *
+ * @param track
+ * @return position at ECAL
+ */
+ public static Hep3Vector getTrackPositionAtEcal(Track track) {
+ return extrapolateTrack(track, BeamlineConstants.ECAL_FACE);
+ }
+
+ /**
+ * Extrapolate track to given position.
+ *
+ * @param helix - to be extrapolated
+ * @param z
+ * @param useMap
+ * @param track - position along the x-axis of the helix in lcsim
+ * coordiantes
+ * @return
+ */
+ public static Hep3Vector extrapolateTrack(Track track, double z) {
+
+ Hep3Vector trackPosition;
+ double dz;
+ if (z >= BeamlineConstants.DIPOLE_EDGE_TESTRUN) {
+ trackPosition = extrapolateHelixToXPlane(track, BeamlineConstants.DIPOLE_EDGE_TESTRUN);
+ dz = z - BeamlineConstants.DIPOLE_EDGE_TESTRUN;
+ } else if (z <= BeamlineConstants.DIPOLE_EDGELOW_TESTRUN) {
+ trackPosition = extrapolateHelixToXPlane(track, BeamlineConstants.DIPOLE_EDGELOW_TESTRUN);
+ dz = z - trackPosition.x();
+ } else {
+ Hep3Vector detVecTracking = extrapolateHelixToXPlane(track, z);
+ // System.out.printf("detVec %s\n", detVecTracking.toString());
+ return new BasicHep3Vector(detVecTracking.y(), detVecTracking.z(), detVecTracking.x());
+ }
+
+ // Get the track azimuthal angle
+ double phi = getPhi(track, trackPosition);
+
+ // Find the distance to the point of interest
+ double r = dz / (getSinTheta(track) * Math.cos(phi));
+ double dx = r * getSinTheta(track) * Math.sin(phi);
+ double dy = r * getCosTheta(track);
+
+ // Find the track position at the point of interest
+ double x = trackPosition.y() + dx;
+ double y = trackPosition.z() + dy;
+
+ return new BasicHep3Vector(x, y, z);
+ }
+
+ /**
+ * Extrapolate track to given position, using dipole position from geometry.
+ *
+ * @param helix - to be extrapolated
+ * @param track - position along the x-axis of the helix in lcsim
+ * coordiantes
+ * @return
+ */
+ public static Hep3Vector extrapolateTrack(Track track, double z, Detector detector) {
+
+ Hep3Vector trackPosition;
+// <constant name="dipoleMagnetPositionZ" value="45.72*cm"/>
+// <constant name="dipoleMagnetLength" value="108*cm"/>
+
+ double magnetLength = detector.getConstants().get("dipoleMagnetLength").getValue();
+ double magnetZ = detector.getConstants().get("dipoleMagnetPositionZ").getValue();
+ double magnetDownstreamEdge = magnetZ + magnetLength / 2;
+ double magnetUpstreamEdge = magnetZ - magnetLength / 2;
+ if (z >= magnetDownstreamEdge)
+ trackPosition = extrapolateHelixToXPlane(track, magnetDownstreamEdge);
+ else if (z <= magnetUpstreamEdge)
+ trackPosition = extrapolateHelixToXPlane(track, magnetUpstreamEdge);
+ else {
+ Hep3Vector detVecTracking = extrapolateHelixToXPlane(track, z);
+ // System.out.printf("detVec %s\n", detVecTracking.toString());
+ return new BasicHep3Vector(detVecTracking.y(), detVecTracking.z(), detVecTracking.x());
+ }
+ double dz = z - trackPosition.x();
+
+ // Get the track azimuthal angle
+ double phi = getPhi(track, trackPosition);
+
+ // Find the distance to the point of interest
+ double r = dz / (getSinTheta(track) * Math.cos(phi));
+ double dx = r * getSinTheta(track) * Math.sin(phi);
+ double dy = r * getCosTheta(track);
+
+ // Find the track position at the point of interest
+ double x = trackPosition.y() + dx;
+ double y = trackPosition.z() + dy;
+
+ return new BasicHep3Vector(x, y, z);
+ }
+
+ /**
+ * Extrapolate helix to given position
+ *
+ * @param helix - to be extrapolated
+ * @param z - position along the x-axis of the helix in lcsim coordiantes
+ * @return
+ */
+ public static Hep3Vector extrapolateTrack(HelicalTrackFit helix, double z) {
+ SeedTrack trk = new SeedTrack();
+ // bfield = Math.abs((detector.getFieldMap().getField(new BasicHep3Vector(0, 0, 0)).y()));
+ double bfield = 0.;
+ // Here we aren't really using anything related to momentum so B-field is not important
+ trk.setTrackParameters(helix.parameters(), bfield); // Sets first TrackState.
+ trk.setCovarianceMatrix(helix.covariance()); // Modifies first TrackState.
+ trk.setChisq(helix.chisqtot());
+ trk.setNDF(helix.ndf()[0] + helix.ndf()[1]);
+ return TrackUtils.extrapolateTrack(trk, z);
+ }
+
+ /**
+ * @param helix input helix object
+ * @param origin of the plane to intercept
+ * @param normal of the plane to intercept
+ * @param eps criteria on the distance to the plane before stopping
+ * iteration
+ * @return position in space at the intercept of the plane
+ */
+ public static Hep3Vector getHelixPlanePositionIter(HelicalTrackFit helix, Hep3Vector origin, Hep3Vector normal, double eps) {
+ boolean debug = false;
+ if (debug) {
+ System.out.printf("--- getHelixPlanePositionIter ---\n");
+ System.out.printf("Target origin [%.10f %.10f %.10f] normal [%.10f %.10f %.10f]\n", origin.x(), origin.y(), origin.z(), normal.x(), normal.y(), normal.z());
+ System.out.printf("%.10f %.10f %.10f %.10f %.10f\n", helix.dca(), helix.z0(), helix.phi0(), helix.slope(), helix.R());
+ }
+ double x = origin.x();
+ double d = 9999.9;
+ double dx = 0.0;
+ int nIter = 0;
+ Hep3Vector pos = null;
+ while (Math.abs(d) > eps && nIter < 50) {
+ // Calculate position on helix at x
+ pos = getHelixPosAtX(helix, x + dx);
+ // Check if we are on the plane
+ d = VecOp.dot(VecOp.sub(pos, origin), normal);
+ dx += -1.0 * d / 2.0;
+ if (debug)
+ System.out.printf("%d d %.10f pos [%.10f %.10f %.10f] dx %.10f\n", nIter, d, pos.x(), pos.y(), pos.z(), dx);
+ nIter += 1;
+ }
+ return pos;
+ }
+
+ /*
+ * Calculates the point on the helix at a given point along the x-axis The normal of the plane
+ * is in the same x-y plane as the circle.
+ * @param helix
+ * @param x point along x-axis
+ * @return point on helix at x-coordinate
+ */
+ private static Hep3Vector getHelixPosAtX(HelicalTrackFit helix, double x) {
+ // double C = (double)Math.round(helix.curvature()*1000000)/1000000;
+ // double R = 1.0/C;
+ double R = helix.R();
+ double dca = helix.dca();
+ double z0 = helix.z0();
+ double phi0 = helix.phi0();
+ double slope = helix.slope();
+ // System.out.printf("%.10f %.10f %.10f %.10f %.10f\n",dca,z0,phi0,slope,R);
+
+ double xc = (R - dca) * Math.sin(phi0);
+ double sinPhi = (xc - x) / R;
+ double phi_at_x = Math.asin(sinPhi);
+ double dphi_at_x = phi_at_x - phi0;
+ if (dphi_at_x > Math.PI)
+ dphi_at_x -= 2.0 * Math.PI;
+ if (dphi_at_x < -Math.PI)
+ dphi_at_x += 2.0 * Math.PI;
+ double s_at_x = -1.0 * dphi_at_x * R;
+ double y = dca * Math.cos(phi0) - R * Math.cos(phi0) + R * Math.cos(phi_at_x);
+ double z = z0 + s_at_x * slope;
+ BasicHep3Vector pos = new BasicHep3Vector(x, y, z);
+ // System.out.printf("pos %s xc %f phi_at_x %f dphi_at_x %f s_at_x %f\n",
+ // pos.toString(),xc,phi_at_x,dphi_at_x,s_at_x);
+ Hep3Vector posXCheck = TrackUtils.extrapolateHelixToXPlane(helix, x);
+ if (VecOp.sub(pos, posXCheck).magnitude() > 0.0000001)
+ throw new RuntimeException(String.format("ERROR the helix propagation equations do not agree? (%f,%f,%f) vs (%f,%f,%f) in HelixUtils", pos.x(), pos.y(), pos.z(), posXCheck.x(), posXCheck.y(), posXCheck.z()));
+ return pos;
+ }
+
+ /**
+ *
+ */
+ public static double findTriangleArea(double x0, double y0, double x1, double y1, double x2, double y2) {
+ return .5 * (x1 * y2 - y1 * x2 - x0 * y2 + y0 * x2 + x0 * y1 - y0 * x1);
+ }
+
+ /**
+ *
+ */
+ public static boolean sensorContainsTrack(Hep3Vector trackPosition, SiSensor sensor) {
+ boolean debug = false;
+ ITransform3D localToGlobal = sensor.getGeometry().getLocalToGlobal();
+
+ Box sensorSolid = (Box) sensor.getGeometry().getLogicalVolume().getSolid();
+ Polygon3D sensorFace = sensorSolid.getFacesNormalTo(new BasicHep3Vector(0, 0, 1)).get(0);
+ if (debug)
+ System.out.println("sensorContainsTrack: Track Position: " + trackPosition.toString());
+
+ List<Point3D> vertices = new ArrayList<Point3D>();
+ for (int index = 0; index < 4; index++)
+ vertices.add(new Point3D());
+ for (Point3D vertex : sensorFace.getVertices())
+ if (vertex.y() < 0 && vertex.x() > 0) {
+ localToGlobal.transform(vertex);
+ // vertices.set(0, new Point3D(vertex.y() + sensorPos.x(), vertex.x() +
+ // sensorPos.y(), vertex.z() + sensorPos.z()));
+ vertices.set(0, new Point3D(vertex.x(), vertex.y(), vertex.z()));
+ if (debug)
+ System.out.println("sensorContainsTrack: Vertex 1 Position: " + vertices.get(0).toString()); // System.out.println("sensorContainsTrack: Transformed Vertex 1 Position: " + // localToGlobal.transformed(vertex).toString());
+ } else if (vertex.y() > 0 && vertex.x() > 0) {
+ localToGlobal.transform(vertex);
+ // vertices.set(1, new Point3D(vertex.y() + sensorPos.x(), vertex.x() +
+ // sensorPos.y(), vertex.z() + sensorPos.z()));
+ vertices.set(1, new Point3D(vertex.x(), vertex.y(), vertex.z()));
+ if (debug)
+ System.out.println("sensorContainsTrack: Vertex 2 Position: " + vertices.get(1).toString()); // System.out.println("sensorContainsTrack: Transformed Vertex 2 Position: " + // localToGlobal.transformed(vertex).toString());
+ } else if (vertex.y() > 0 && vertex.x() < 0) {
+ localToGlobal.transform(vertex);
+ // vertices.set(2, new Point3D(vertex.y() + sensorPos.x(), vertex.x() +
+ // sensorPos.y(), vertex.z() + sensorPos.z()));
+ vertices.set(2, new Point3D(vertex.x(), vertex.y(), vertex.z()));
+ if (debug)
+ System.out.println("sensorContainsTrack: Vertex 3 Position: " + vertices.get(2).toString()); // System.out.println("sensorContainsTrack: Transformed Vertex 3 Position: " + // localToGlobal.transformed(vertex).toString());
+ } else if (vertex.y() < 0 && vertex.x() < 0) {
+ localToGlobal.transform(vertex);
+ // vertices.set(3, new Point3D(vertex.y() + sensorPos.x(), vertex.x() +
+ // sensorPos.y(), vertex.z() + sensorPos.z()));
+ vertices.set(3, new Point3D(vertex.x(), vertex.y(), vertex.z()));
+ if (debug)
+ System.out.println("sensorContainsTrack: Vertex 4 Position: " + vertices.get(3).toString()); // System.out.println("sensorContainsTrack: Transformed Vertex 4 Position: " + // localToGlobal.transformed(vertex).toString());
+ }
+
+ double area1 = TrackUtils.findTriangleArea(vertices.get(0).x(), vertices.get(0).y(), vertices.get(1).x(), vertices.get(1).y(), trackPosition.y(), trackPosition.z());
+ double area2 = TrackUtils.findTriangleArea(vertices.get(1).x(), vertices.get(1).y(), vertices.get(2).x(), vertices.get(2).y(), trackPosition.y(), trackPosition.z());
+ double area3 = TrackUtils.findTriangleArea(vertices.get(2).x(), vertices.get(2).y(), vertices.get(3).x(), vertices.get(3).y(), trackPosition.y(), trackPosition.z());
+ double area4 = TrackUtils.findTriangleArea(vertices.get(3).x(), vertices.get(3).y(), vertices.get(0).x(), vertices.get(0).y(), trackPosition.y(), trackPosition.z());
+
+ return (area1 > 0 && area2 > 0 && area3 > 0 && area4 > 0) || (area1 < 0 && area2 < 0 && area3 < 0 && area4 < 0);
+ }
+
+ public static Map<String, Double> calculateTrackHitResidual(HelicalTrackHit hth, HelicalTrackFit track, boolean includeMS) {
+
+ boolean debug = false;
+ Map<String, Double> residuals = new HashMap<String, Double>();
+
+ Map<HelicalTrackHit, MultipleScatter> msmap = track.ScatterMap();
+ double msdrphi = 0;
+ double msdz = 0;
+
+ if (includeMS) {
+ msdrphi = msmap.get(hth).drphi();
+ msdz = msmap.get(hth).dz();
+ }
+
+ // Calculate the residuals that are being used in the track fit
+ // Start with the bendplane y
+ double drphi_res = hth.drphi();
+ double wrphi = Math.sqrt(drphi_res * drphi_res + msdrphi * msdrphi);
+ // This is the normal way to get s
+ double s_wrong = track.PathMap().get(hth);
+ // This is how I do it with HelicalTrackFits
+ double s = HelixUtils.PathToXPlane(track, hth.x(), 0, 0).get(0);
+ // System.out.printf("x %f s %f smap %f\n",hth.x(),s,s_wrong);
+ if (Double.isNaN(s)) {
+ double xc = track.xc();
+ double RC = track.R();
+ System.out.printf("calculateTrackHitResidual: s is NaN. p=%.3f RC=%.3f, x=%.3f, xc=%.3f\n", track.p(-0.491), RC, hth.x(), xc);
+ return residuals;
+ }
+
+ Hep3Vector posOnHelix = HelixUtils.PointOnHelix(track, s);
+ double resy = hth.y() - posOnHelix.y();
+ double erry = includeMS ? wrphi : drphi_res;
+
+ // Now the residual for the "measurement" direction z
+ double resz = hth.z() - posOnHelix.z();
+ double dz_res = HitUtils.zres(hth, msmap, track);
+ double dz_res2 = hth.getCorrectedCovMatrix().diagonal(2);
+
+ if (Double.isNaN(resy)) {
+ System.out.printf("calculateTrackHitResidual: resy is NaN. hit at %s posOnHelix=%s path=%.3f wrong_path=%.3f helix:\n%s\n", hth.getCorrectedPosition().toString(), posOnHelix.toString(), s, s_wrong, track.toString());
+ return residuals;
+ }
+
+ residuals.put("resy", resy);
+ residuals.put("erry", erry);
+ residuals.put("drphi", drphi_res);
+ residuals.put("msdrphi", msdrphi);
+
+ residuals.put("resz", resz);
+ residuals.put("errz", dz_res);
+ residuals.put("dz_res", Math.sqrt(dz_res2));
+ residuals.put("msdz", msdz);
+
+ if (debug) {
+ System.out.printf("calculateTrackHitResidual: HTH hit at (%f,%f,%f)\n", hth.x(), hth.y(), hth.z());
+ System.out.printf("calculateTrackHitResidual: helix params d0=%f phi0=%f R=%f z0=%f slope=%f chi2=%f/%f chi2tot=%f\n", track.dca(), track.phi0(), track.R(), track.z0(), track.slope(), track.chisq()[0], track.chisq()[1], track.chisqtot());
+ System.out.printf("calculateTrackHitResidual: => resz=%f resy=%f at s=%f\n", resz, resy, s);
+ // System.out.printf("calculateTrackHitResidual: resy=%f eresy=%f drphi=%f msdrphi=%f \n",resy,erry,drphi_res,msdrphi);
+ // System.out.printf("calculateTrackHitResidual: resz=%f eresz=%f dz_res=%f msdz=%f \n",resz,dz_res,Math.sqrt(dz_res2),msdz);
+ }
+
+ return residuals;
+ }
+
+ public static Map<String, Double> calculateLocalTrackHitResiduals(Track track, HelicalTrackHit hth, HelicalTrackStrip strip, double bFieldInZ) {
+ HelicalTrackStripGbl stripGbl = new HelicalTrackStripGbl(strip, true);
+ return calculateLocalTrackHitResiduals(track, hth, stripGbl, bFieldInZ);
+ }
+
+ public static Map<String, Double> calculateLocalTrackHitResiduals(Track track, HelicalTrackHit hth, HelicalTrackStripGbl strip, double bFieldInZ) {
+
+ SeedTrack st = (SeedTrack) track;
+ SeedCandidate seed = st.getSeedCandidate();
+ HelicalTrackFit _trk = seed.getHelix();
+ Map<HelicalTrackHit, MultipleScatter> msmap = seed.getMSMap();
+ double msdrdphi = msmap.get(hth).drphi();
+ double msdz = msmap.get(hth).dz();
+ return calculateLocalTrackHitResiduals(_trk, strip, msdrdphi, msdz, bFieldInZ);
+ }
+
+ public static Map<String, Double> calculateLocalTrackHitResiduals(HelicalTrackFit _trk, HelicalTrackStrip strip, double bFieldInZ) {
+ HelicalTrackStripGbl stripGbl = new HelicalTrackStripGbl(strip, true);
+ return calculateLocalTrackHitResiduals(_trk, stripGbl, 0.0, 0.0, bFieldInZ);
+ }
+
+ public static Map<String, Double> calculateLocalTrackHitResiduals(HelicalTrackFit _trk, HelicalTrackStripGbl strip, double msdrdphi, double msdz, double bFieldInZ) {
+
+ boolean debug = false;
+ boolean includeMS = true;
+
+ if (debug)
+ System.out.printf("calculateLocalTrackHitResiduals: for strip on sensor %s \n",
+ ((RawTrackerHit) strip.getStrip().rawhits().get(0)).getDetectorElement().getName());
+
+ Hep3Vector u = strip.u();
+ Hep3Vector corigin = strip.origin();
+
+ // Find interception with plane that the strips belongs to
+ Hep3Vector trkpos = TrackUtils.getHelixPlaneIntercept(_trk, strip, Math.abs(bFieldInZ));
+
+ if (debug) {
+ System.out.printf("calculateLocalTrackHitResiduals: strip u %s origin %s \n", u.toString(), corigin.toString());
+ System.out.printf("calculateLocalTrackHitResiduals: found interception point with sensor at %s \n", trkpos.toString());
+ }
+
+ if (Double.isNaN(trkpos.x()) || Double.isNaN(trkpos.y()) || Double.isNaN(trkpos.z())) {
+ System.out.printf("calculateLocalTrackHitResiduals: failed to get interception point (%s) \n", trkpos.toString());
+ System.out.printf("calculateLocalTrackHitResiduals: track params\n%s\n", _trk.toString());
+ System.out.printf("calculateLocalTrackHitResiduals: track pT=%.3f chi2=[%.3f][%.3f] \n", _trk.pT(bFieldInZ), _trk.chisq()[0], _trk.chisq()[1]);
+// trkpos = TrackUtils.getHelixPlaneIntercept(_trk, strip, bFieldInZ);
+ throw new RuntimeException();
+ }
+
+ double xint = trkpos.x();
+ double phi0 = _trk.phi0();
+ double R = _trk.R();
+ double s = HelixUtils.PathToXPlane(_trk, xint, 0, 0).get(0);
+ double phi = -s / R + phi0;
+
+ Hep3Vector mserr = new BasicHep3Vector(msdrdphi * Math.sin(phi), msdrdphi * Math.sin(phi), msdz);
+ double msuError = VecOp.dot(mserr, u);
+
+ Hep3Vector vdiffTrk = VecOp.sub(trkpos, corigin);
+ TrackerHitUtils thu = new TrackerHitUtils(debug);
+ Hep3Matrix trkToStrip = thu.getTrackToStripRotation(strip.getStrip());
+ Hep3Vector vdiff = VecOp.mult(trkToStrip, vdiffTrk);
+
+ double umc = vdiff.x();
+ double vmc = vdiff.y();
+ double wmc = vdiff.z();
+ double umeas = strip.umeas();
+ double uError = strip.du();
+ double vmeas = 0;
+ double vError = (strip.vmax() - strip.vmin()) / Math.sqrt(12);
+ double wmeas = 0;
+ double wError = 10.0 / Math.sqrt(12); // 0.001;
+
+ if (debug)
+ System.out.printf("calculateLocalTrackHitResiduals: vdiffTrk %s vdiff %s umc %f umeas %f du %f\n",
+ vdiffTrk.toString(), vdiff.toString(), umc, umeas, umeas - umc);
+
+ Map<String, Double> res = new HashMap<String, Double>();
+ res.put("ures", umeas - umc);
+ res.put("ureserr", includeMS ? Math.sqrt(uError * uError + msuError * msuError) : uError);
+ res.put("vres", vmeas - vmc);
+ res.put("vreserr", vError);
+ res.put("wres", wmeas - wmc);
+ res.put("wreserr", wError);
+
+ res.put("vdiffTrky", vdiffTrk.y());
+
+ return res;
+ }
+
+ public static int[] getHitsInTopBottom(Track track) {
+ int n[] = {0, 0};
+ List<TrackerHit> hitsOnTrack = track.getTrackerHits();
+ for (TrackerHit hit : hitsOnTrack) {
+ HelicalTrackHit hth = (HelicalTrackHit) hit;
+ //===> if (SvtUtils.getInstance().isTopLayer((SiSensor) ((RawTrackerHit) hth.getRawHits().get(0)).getDetectorElement())) {
+ HpsSiSensor sensor = ((HpsSiSensor) ((RawTrackerHit) hth.getRawHits().get(0)).getDetectorElement());
+ if (sensor.isTopLayer())
+ n[0] = n[0] + 1;
+ else
+ n[1] = n[1] + 1;
+ }
+ return n;
+ }
+
+ public static boolean isTopTrack(Track track, int minhits) {
+ return isTopOrBottomTrack(track, minhits) == 1;
+ }
+
+ public static boolean isBottomTrack(Track track, int minhits) {
+ return isTopOrBottomTrack(track, minhits) == 0;
+ }
+
+ public static int isTopOrBottomTrack(Track track, int minhits) {
+ int nhits[] = getHitsInTopBottom(track);
+ if (nhits[0] >= minhits && nhits[1] == 0)
+ return 1;
+ else if (nhits[1] >= minhits && nhits[0] == 0)
+ return 0;
+ else
+ return -1;
+ }
+
+ public static boolean hasTopBotHit(Track track) {
+ int nhits[] = getHitsInTopBottom(track);
+ return nhits[0] > 0 && nhits[1] > 0;
+ }
+
+ public static boolean isSharedHit(TrackerHit hit, List<Track> othertracks) {
+ HelicalTrackHit hth = (HelicalTrackHit) hit;
+ for (Track track : othertracks) {
+ List<TrackerHit> hitsOnTrack = track.getTrackerHits();
+ for (TrackerHit loop_hit : hitsOnTrack) {
+ HelicalTrackHit loop_hth = (HelicalTrackHit) loop_hit;
+ if (hth.equals(loop_hth)) // System.out.printf("share hit at layer %d at %s (%s) with track w/ chi2=%f\n",hth.Layer(),hth.getCorrectedPosition().toString(),loop_hth.getCorrectedPosition().toString(),track.getChi2());
+
+ return true;
+ }
+ }
+ return false;
+ }
+
+ public static int numberOfSharedHits(Track track, List<Track> tracklist) {
+ List<Track> tracks = new ArrayList<Track>();
+ // System.out.printf("%d tracks in event\n",tracklist.size());
+ // System.out.printf("look for another track with chi2=%f and px=%f \n",track.getChi2(),track.getTrackStates().get(0).getMomentum()[0]);
+ for (Track t : tracklist) {
+ // System.out.printf("add track with chi2=%f and px=%f ?\n",t.getChi2(),t.getTrackStates().get(0).getMomentum()[0]);
+ if (t.equals(track)) // System.out.printf("NOPE\n");
+
+ continue;
+ // System.out.printf("YEPP\n");
+ tracks.add(t);
+ }
+ List<TrackerHit> hitsOnTrack = track.getTrackerHits();
+ int n_shared = 0;
+ for (TrackerHit hit : hitsOnTrack)
+ if (isSharedHit(hit, tracks))
+ ++n_shared;
+ return n_shared;
+ }
+
+ public static boolean hasSharedHits(Track track, List<Track> tracklist) {
+ return numberOfSharedHits(track, tracklist) != 0;
+ }
+
+ public static void cut(int cuts[], EventQuality.Cut bit) {
+ cuts[0] = cuts[0] | (1 << bit.getValue());
+ }
+
+ public static boolean isGoodTrack(Track track, List<Track> tracklist, EventQuality.Quality trk_quality) {
+ int cuts = passTrackSelections(track, tracklist, trk_quality);
+ return cuts == 0;
+ }
+
+ public static int passTrackSelections(Track track, List<Track> tracklist, EventQuality.Quality trk_quality) {
+ int cuts[] = {0};
+ if (trk_quality.compareTo(Quality.NONE) != 0) {
+ if (track.getTrackStates().get(0).getMomentum()[0] < EventQuality.instance().getCutValue(EventQuality.Cut.PZ, trk_quality))
+ cut(cuts, EventQuality.Cut.PZ);
+ if (track.getChi2() >= EventQuality.instance().getCutValue(EventQuality.Cut.CHI2, trk_quality))
+ cut(cuts, EventQuality.Cut.CHI2);
+ if (numberOfSharedHits(track, tracklist) > ((int) Math.round(EventQuality.instance().getCutValue(EventQuality.Cut.SHAREDHIT, trk_quality))))
+ cut(cuts, EventQuality.Cut.SHAREDHIT);
+ if (hasTopBotHit(track))
+ cut(cuts, EventQuality.Cut.TOPBOTHIT);
+ if (track.getTrackerHits().size() < ((int) Math.round(EventQuality.instance().getCutValue(EventQuality.Cut.NHITS, trk_quality))))
+ cut(cuts, EventQuality.Cut.NHITS);
+ }
+ return cuts[0];
+ }
+
+ public static boolean isTopTrack(HelicalTrackFit htf) {
+ return htf.slope() > 0;
+ }
+
+ public static boolean isBottomTrack(HelicalTrackFit htf) {
+ return !isTopTrack(htf);
+ }
+
+ /**
+ * Transform MCParticle into a Helix object. Note that it produces the helix
+ * parameters at nominal x=0 and assumes that there is no field at x<0
+ *
+ * @param mcp MC particle to be transformed
+ * @return helix object based on the MC particle
+ */
+ public static HelicalTrackFit getHTF(MCParticle mcp, double Bz) {
+ boolean debug = true;
+ if (debug) {
+ System.out.printf("getHTF\n");
+ System.out.printf("mcp org %s mc p %s\n", mcp.getOrigin().toString(), mcp.getMomentum().toString());
+ }
+ Hep3Vector org = CoordinateTransformations.transformVectorToTracking(mcp.getOrigin());
+ Hep3Vector p = CoordinateTransformations.transformVectorToTracking(mcp.getMomentum());
+
+ if (debug)
+ System.out.printf("mcp org %s mc p %s (trans)\n", org.toString(), p.toString());
+
+ // Move to x=0 if needed
+ double targetX = BeamlineConstants.DIPOLE_EDGELOW_TESTRUN;
+ if (org.x() < targetX) {
+ double dydx = p.y() / p.x();
+ double dzdx = p.z() / p.x();
+ double delta_x = targetX - org.x();
+ double y = delta_x * dydx + org.y();
+ double z = delta_x * dzdx + org.z();
+ double x = org.x() + delta_x;
+ if (Math.abs(x - targetX) > 1e-8)
+ throw new RuntimeException("Error: origin is not zero!");
+ org = new BasicHep3Vector(x, y, z);
+ // System.out.printf("org %s p %s -> org %s\n",
+ // old.toString(),p.toString(),org.toString());
+ }
+
+ if (debug)
+ System.out.printf("mcp org %s mc p %s (trans2)\n", org.toString(), p.toString());
+
+ HelixParamCalculator helixParamCalculator = new HelixParamCalculator(p, org, -1 * ((int) mcp.getCharge()), Bz);
+ double par[] = new double[5];
+ par[HelicalTrackFit.dcaIndex] = helixParamCalculator.getDCA();
+ par[HelicalTrackFit.slopeIndex] = helixParamCalculator.getSlopeSZPlane();
+ par[HelicalTrackFit.phi0Index] = helixParamCalculator.getPhi0();
+ par[HelicalTrackFit.curvatureIndex] = 1.0 / helixParamCalculator.getRadius();
+ par[HelicalTrackFit.z0Index] = helixParamCalculator.getZ0();
+ HelicalTrackFit htf = getHTF(par);
+ System.out.printf("d0 %f z0 %f R %f phi %f lambda %s\n",
+ htf.dca(), htf.z0(), htf.R(), htf.phi0(), htf.slope());
+ return htf;
+ }
+
+ public static HelicalTrackFit getHTF(Track track) {
+ if (track.getClass().isInstance(SeedTrack.class))
+ return ((SeedTrack) track).getSeedCandidate().getHelix();
+ else
+ return getHTF(track.getTrackStates().get(0));
+ }
+
+ public static HelicalTrackFit getHTF(double par[]) {
+ // need to have matrix that makes sense? Really?
+ SymmetricMatrix cov = new SymmetricMatrix(5);
+ for (int i = 0; i < cov.getNRows(); ++i)
+ cov.setElement(i, i, 1.);
+ HelicalTrackFit htf = new HelicalTrackFit(par, cov, new double[2], new int[2], null, null);
+ return htf;
+ }
+
+ public static HelicalTrackFit getHTF(TrackState state) {
+ double par[] = state.getParameters();
+ SymmetricMatrix cov = new SymmetricMatrix(5, state.getCovMatrix(), true);
+ HelicalTrackFit htf = new HelicalTrackFit(par, cov, new double[2], new int[2], null, null);
+ return htf;
+ }
+
+ public static StraightLineTrack findSLTAtZ(Track trk1, double zVal, boolean useFringe) {
+ SeedTrack s1 = (SeedTrack) trk1;
+ HelicalTrackFit htf1 = s1.getSeedCandidate().getHelix();
+ HPSTrack hpstrk1 = new HPSTrack(htf1);
+ Hep3Vector pos1;
+ if (useFringe) {
+ // broken because you need ot provide the Field Map to get this...
+// pos1 = hpstrk1.getPositionAtZMap(100.0, zVal, 5.0)[0];
+ } else
+ pos1 = TrackUtils.extrapolateTrack(trk1, zVal);
+ // System.out.printf("%s: Position1 at edge of fringe %s\n",this.getClass().getSimpleName(),pos1.toString());
+ Helix traj = (Helix) hpstrk1.getTrajectory();
+ if (traj == null) {
+ SpacePoint r0 = new SpacePoint(HelixUtils.PointOnHelix(htf1, 0));
+ traj = new Helix(r0, htf1.R(), htf1.phi0(), Math.atan(htf1.slope()));
+ }
+ HelixConverter converter = new HelixConverter(0.);
+ StraightLineTrack slt1 = converter.Convert(traj);
+ // System.out.printf("%s: straight line track: x0=%f,y0=%f,z0=%f dz/dx=%f dydx=%f targetY=%f targetZ=%f \n",this.getClass().getSimpleName(),slt1.x0(),slt1.y0(),slt1.z0(),slt1.dzdx(),slt1.dydx(),slt1.TargetYZ()[0],slt1.TargetYZ()[1]);
+ return slt1;
+ }
+
+ public static MCParticle getMatchedTruthParticle(Track track) {
+ boolean debug = false;
+
+ Map<MCParticle, Integer> particlesOnTrack = new HashMap<MCParticle, Integer>();
+
+ if (debug)
+ System.out.printf("getMatchedTruthParticle: getmatched mc particle from %d tracker hits on the track \n", track.getTrackerHits().size());
+
+ for (TrackerHit hit : track.getTrackerHits()) {
+ List<MCParticle> mcps = ((HelicalTrackHit) hit).getMCParticles();
+ if (mcps == null)
+ System.out.printf("getMatchedTruthParticle: warning, this hit (layer %d pos=%s) has no mc particles.\n", ((HelicalTrackHit) hit).Layer(), ((HelicalTrackHit) hit).getCorrectedPosition().toString());
+ else {
+ if (debug)
+ System.out.printf("getMatchedTruthParticle: this hit (layer %d pos=%s) has %d mc particles.\n", ((HelicalTrackHit) hit).Layer(), ((HelicalTrackHit) hit).getCorrectedPosition().toString(), mcps.size());
+ for (MCParticle mcp : mcps) {
+ if (!particlesOnTrack.containsKey(mcp))
+ particlesOnTrack.put(mcp, 0);
+ int c = particlesOnTrack.get(mcp);
+ particlesOnTrack.put(mcp, c + 1);
+ }
+ }
+ }
+ if (debug) {
+ System.out.printf("Track p=[ %f, %f, %f] \n", track.getTrackStates().get(0).getMomentum()[0], track.getTrackStates().get(0).getMomentum()[1], track.getTrackStates().get(0).getMomentum()[1]);
+ System.out.printf("Found %d particles\n", particlesOnTrack.size());
+ for (Map.Entry<MCParticle, Integer> entry : particlesOnTrack.entrySet())
+ System.out.printf("%d hits assigned to %d p=%s \n", entry.getValue(), entry.getKey().getPDGID(), entry.getKey().getMomentum().toString());
+ }
+ Map.Entry<MCParticle, Integer> maxEntry = null;
+ for (Map.Entry<MCParticle, Integer> entry : particlesOnTrack.entrySet())
+ if (maxEntry == null || entry.getValue().compareTo(maxEntry.getValue()) > 0)
+ maxEntry = entry; //if ( maxEntry != null ) { // if(entry.getValue().compareTo(maxEntry.getValue()) < 0) continue; //}
+ //maxEntry = entry;
+ if (debug)
+ if (maxEntry != null)
+ System.out.printf("Matched particle with pdgId=%d and mom %s to track with charge %d and momentum [%f %f %f]\n",
+ maxEntry.getKey().getPDGID(), maxEntry.getKey().getMomentum().toString(),
+ track.getCharge(), track.getTrackStates().get(0).getMomentum()[0], track.getTrackStates().get(0).getMomentum()[1], track.getTrackStates().get(0).getMomentum()[2]);
+ else
+ System.out.printf("No truth particle found on this track\n");
+ return maxEntry == null ? null : maxEntry.getKey();
+ }
+
+ /*
+ try to make a seed track from a base track
+ ...some things are irretrivable (tracking strategy)..
+ and some things I just don't care much about to dig out (det name)
+ */
+ public static SeedTrack makeSeedTrackFromBaseTrack(Track track) {
+
+ TrackState trkState = track.getTrackStates().get(0);
+ // first make the HelicalTrackFit Object
+ double[] covArray = trkState.getCovMatrix();
+ SymmetricMatrix cov = new SymmetricMatrix(5, covArray, true);
+ double[] chisq = {track.getChi2(), 0};
+ int[] ndf = {track.getNDF(), 0};
+ Map<HelicalTrackHit, Double> smap = new HashMap<>(); // just leave these empty
+ Map<HelicalTrackHit, MultipleScatter> msmap = new HashMap<>();// just leave these empty
+ double[] pars = {trkState.getD0(), trkState.getPhi(), trkState.getOmega(), trkState.getZ0(), trkState.getTanLambda()};
+ HelicalTrackFit htf = new HelicalTrackFit(pars, cov, chisq, ndf, smap, msmap);
+ // now get the hits and make them helicaltrackhits
+ List<TrackerHit> rth = track.getTrackerHits();
+ List<HelicalTrackHit> hth = new ArrayList<>();
+ for (TrackerHit hit : rth)
+ hth.add(makeHelicalTrackHitFromTrackerHit(hit));
+// SeedCandidate(List<HelicalTrackHit> , SeedStrategy strategy, HelicalTrackFit helix, double bfield) ;
+ SeedCandidate scand = new SeedCandidate(hth, null, htf, 0.24);
+ SeedTrack st = new SeedTrack();
+ st.setSeedCandidate(scand);
+ return st;
+ }
+
+ /*
+ cast TrackerHit as HTH...this is pretty dumb; just a rearrangment of information
+ in TrackerHit. The important information that's in HTH but not
+ in Tracker hit is the HelicalTrackCrosses (and from there the individual strip clusters)
+ is lost; some work to get them back.
+ */
+ public static HelicalTrackHit makeHelicalTrackHitFromTrackerHit(TrackerHit hit) {
+ Hep3Vector pos = new BasicHep3Vector(hit.getPosition());
+ SymmetricMatrix hitcov = new SymmetricMatrix(3, hit.getCovMatrix(), true);
+ double dedx = hit.getdEdx();
+ double time = hit.getTime();
+ int type = hit.getType();
+ List rhits = hit.getRawHits();
+ String detname = "Foobar";
+ int layer = 666;
+ BarrelEndcapFlag beflag = BarrelEndcapFlag.BARREL;
+ return new HelicalTrackHit(pos, hitcov, dedx, time, type, rhits, detname, layer, beflag);
+ }
+
+ public static RelationalTable getHitToStripsTable(EventHeader event) {
+ RelationalTable hitToStrips = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
+ List<LCRelation> hitrelations = event.get(LCRelation.class, "HelicalTrackHitRelations");
+ for (LCRelation relation : hitrelations)
+ if (relation != null && relation.getFrom() != null && relation.getTo() != null)
+ hitToStrips.add(relation.getFrom(), relation.getTo());
+
+ return hitToStrips;
+ }
+
+ public static RelationalTable getHitToRotatedTable(EventHeader event) {
+
+ RelationalTable hitToRotated = new BaseRelationalTable(RelationalTable.Mode.ONE_TO_ONE, RelationalTable.Weighting.UNWEIGHTED);
+ List<LCRelation> rotaterelations = event.get(LCRelation.class, "RotatedHelicalTrackHitRelations");
+ for (LCRelation relation : rotaterelations)
+ if (relation != null && relation.getFrom() != null && relation.getTo() != null)
+ hitToRotated.add(relation.getFrom(), relation.getTo());
+ return hitToRotated;
+ }
+
+ public static double getTrackTime(Track track, RelationalTable hitToStrips, RelationalTable hitToRotated) {
+ int nStrips = 0;
+ double meanTime = 0;
+ for (TrackerHit hit : track.getTrackerHits()) {
+ Set<TrackerHit> htsList = hitToStrips.allFrom(hitToRotated.from(hit));
+ for (TrackerHit hts : htsList) {
+ nStrips++;
+ meanTime += hts.getTime();
+ }
+ }
+ meanTime /= nStrips;
+ return meanTime;
+ }
+
+ public static List<TrackerHit> getStripHits(Track track, RelationalTable hitToStrips, RelationalTable hitToRotated) {
+ List<TrackerHit> hits = new ArrayList<TrackerHit>();
+ for (TrackerHit hit : track.getTrackerHits())
+ hits.addAll(hitToStrips.allFrom(hitToRotated.from(hit)));
+ return hits;
+ }
+
+ public static boolean hasSharedStrips(Track track1, Track track2, RelationalTable hitToStrips, RelationalTable hitToRotated) {
+ Set<TrackerHit> track1hits = new HashSet<TrackerHit>(getStripHits(track1, hitToStrips, hitToRotated));
+ for (TrackerHit hit : track2.getTrackerHits())
+ for (TrackerHit hts : (Set<TrackerHit>) hitToStrips.allFrom(hitToRotated.from(hit)))
+ if (track1hits.contains(hts))
+ return true;
+ return false;
+ }
+
+ public static int getLayer(TrackerHit strip) {
+ return ((RawTrackerHit) strip.getRawHits().get(0)).getLayerNumber();
+ }
+
+ /**
+ * Compute strip isolation, defined as the minimum distance to another strip
+ * in the same sensor. Strips are only checked if they formed valid crosses
+ * with the other strip in the cross (passing time and tolerance cuts).
+ *
+ * @param strip The strip whose isolation is being calculated.
+ * @param otherStrip The other strip in the stereo hit.
+ * @param hitToStrips
+ * @param hitToRotated
+ * @return Double_MAX_VALUE if no other strips found.
+ */
+ public static double getIsolation(TrackerHit strip, TrackerHit otherStrip, RelationalTable hitToStrips, RelationalTable hitToRotated) {
+ double nearestDistance = 99999999.0;
+ for (TrackerHit cross : (Set<TrackerHit>) hitToStrips.allTo(otherStrip))
+ for (TrackerHit crossStrip : (Set<TrackerHit>) hitToStrips.allFrom(cross))
+ if (crossStrip != strip && crossStrip != otherStrip) {
+ int stripMin = Integer.MAX_VALUE;
+ int stripMax = Integer.MIN_VALUE;
+ int crossMin = Integer.MAX_VALUE;
+ int crossMax = Integer.MIN_VALUE;
+ for (Object rth : strip.getRawHits()) {
+ int hitStrip = ((RawTrackerHit) rth).getIdentifierFieldValue("strip");
+ stripMin = Math.min(stripMin, hitStrip);
+ stripMax = Math.max(stripMax, hitStrip);
+ }
+ for (Object rth : crossStrip.getRawHits()) {
+ int hitStrip = ((RawTrackerHit) rth).getIdentifierFieldValue("strip");
+ crossMin = Math.min(crossMin, hitStrip);
+ crossMax = Math.max(crossMax, hitStrip);
+ }
+ if (stripMin - crossMax <= 1 && crossMin - stripMax <= 1)
+ continue; //adjacent strips don't count
+ Hep3Vector stripPosition = new BasicHep3Vector(strip.getPosition());
+ Hep3Vector crossStripPosition = new BasicHep3Vector(crossStrip.getPosition());
+ double distance = VecOp.sub(stripPosition, crossStripPosition).magnitude();
+ if (Math.abs(stripPosition.y()) > Math.abs(crossStripPosition.y()))
+ distance = -distance;
+// System.out.format("%s, %s, %s, %f\n", stripPosition, crossStripPosition, VecOp.sub(stripPosition, crossStripPosition), distance);
+// if (distance<=0.0601) continue; //hack to avoid counting adjacent strips that didn't get clustered together
+// if (distance<0.1) System.out.format("%d, %d, %f\n",strip.getRawHits().size(),crossStrip.getRawHits().size(), distance);
+// if (distance < 0.1) {
+// System.out.format("%s, %s, %s, %f\n", stripPosition, crossStripPosition, VecOp.sub(stripPosition, crossStripPosition), distance);
+// }
+ if (Math.abs(distance) < Math.abs(nearestDistance))
+ nearestDistance = distance;
+ }
+ return nearestDistance;
+ }
+
+ /**
+ * Make an array with isolations for all 12 strip layers. If the track
+ * doesn't use a layer, the isolation is null; if there is no other strip
+ * hit in a layer, the isolation is Double.MAX_VALUE.
+ *
+ * @param trk
+ * @param hitToStrips
+ * @param hitToRotated
+ * @return
+ */
+ public static Double[] getIsolations(Track trk, RelationalTable hitToStrips, RelationalTable hitToRotated) {
+ Double[] isolations = new Double[12];
+ for (TrackerHit hit : trk.getTrackerHits()) {
+ Set<TrackerHit> htsList = hitToStrips.allFrom(hitToRotated.from(hit));
+ TrackerHit[] strips = new TrackerHit[2];
+ htsList.toArray(strips);
+ isolations[TrackUtils.getLayer(strips[0]) - 1] = TrackUtils.getIsolation(strips[0], strips[1], hitToStrips, hitToRotated);
+ isolations[TrackUtils.getLayer(strips[1]) - 1] = TrackUtils.getIsolation(strips[1], strips[0], hitToStrips, hitToRotated);
+ }
+ return isolations;
+ }
+
+ /**
+ * iteratively extrapolates a track to a specified value of z using the full
+ * field map.
+ *
+ * @param track
+ * @param start = value to start extrapolation (use analytic=constant field
+ * calculation up to here)
+ * @param zFinal = value to extrapolate to
+ * @param step = step size
+ * @param bmap = the field map
+ * @param debugOk
+ * @return
+ */
+ public static TrackState extrapolateTrackUsingFieldMap(Track track, double start, double zFinal, double step, FieldMap bmap, boolean debugOk) {
+ Trajectory _trajectory;
+ double startFringe = start;
+ HelicalTrackFit helix = getHTF(track);
+
+ // if looking upstream, we'll be propagating backwards
+ if (zFinal < 0)
+ step = -step;
+ if (debugOk)
+ System.out.println(track.toString());
+
+ double sStartFringe = HelixUtils.PathToXPlane(helix, startFringe, 1000.0, 1).get(0);
+ if (debugOk)
+ System.out.println("path to end of fringe = " + sStartFringe + "; zFinal = " + zFinal);
+ double xtmp = startFringe;
+ double ytmp = HelixUtils.PointOnHelix(helix, sStartFringe).y();
+ double ztmp = HelixUtils.PointOnHelix(helix, sStartFringe).z();
+ double Rorig = helix.R();
+ double xCtmp = helix.xc();
+ double yCtmp = helix.yc();
+ double q = Math.signum(helix.curvature());
+ double phitmp = calculatePhi(xtmp, ytmp, xCtmp, yCtmp, q);
+ if (debugOk) {
+ System.out.println("\nOriginal xtmp = " + xtmp + "; ytmp = " + ytmp + "; ztmp = " + ztmp + "; phitmp = " + phitmp);
+
+ System.out.println("nOriginal Rorig = " + Rorig + "; xCtmp = " + xCtmp + "; yCtmp = " + yCtmp);
+ }
+ if (debugOk)
+ System.out.println("Original Direction at Fringe: " + HelixUtils.Direction(helix, startFringe).toString());
+ double Rtmp = Rorig;
+ // now start stepping through the fringe field
+ Hep3Vector r0Tmp = HelixUtils.PointOnHelix(helix, sStartFringe);
+ org.lcsim.spacegeom.SpacePoint r0 = new org.lcsim.spacegeom.SpacePoint(r0Tmp);
+ Hep3Vector bMid = bmap.getField(new BasicHep3Vector(0, 0, 500.0));
+ double pTot = helix.p(bMid.y());//50 cm ~ middle of dipole
+ Hep3Vector dirOrig = HelixUtils.Direction(helix, sStartFringe);
+ Hep3Vector p0 = VecOp.mult(pTot, dirOrig);
+ Hep3Vector dirTmp = dirOrig;
+ org.lcsim.spacegeom.SpacePoint rTmp = r0;
+ Hep3Vector pTmp = p0;
+ double pXOrig = p0.x();
+ double pXTmp = pXOrig;
+ double totalS = sStartFringe;
+ double myBField = bMid.y();
+ // follow trajectory while: in fringe field, before end point, and we don't have a looper
+ while (Math.signum(step) * xtmp < Math.signum(step) * zFinal && Math.signum(pXOrig * pXTmp) > 0) {
+ if (debugOk) {
+ System.out.println("New step in Fringe Field");
+ System.out.println("rTmp = " + rTmp.toString());
+ System.out.println("pTmp = " + pTmp.toString());
+ System.out.println("OriginalHelix pos = " + HelixUtils.PointOnHelix(helix, totalS));
+ System.out.println("OriginalHelix Momentum = " + VecOp.mult(pTot, HelixUtils.Direction(helix, totalS)));
+ }
+
+ myBField = bmap.getField(new BasicHep3Vector(rTmp.y(), rTmp.z(), rTmp.x())).y(); // note weird co ordinates for track vs det
+ if (debugOk)
+ System.out.println("rTmp.x() = " + rTmp.x() + " field = " + myBField);
+ _trajectory = getTrajectory(pTmp, rTmp, q, myBField);
+ rTmp = _trajectory.getPointAtDistance(step);
+ pTmp = VecOp.mult(pTot, _trajectory.getUnitTangentAtLength(step));
+ pXTmp = pTmp.x();
+ xtmp = rTmp.x();
+ if (debugOk) {
+ System.out.println("############## done... #############");
+
+ System.out.println("\n");
+ }
+ totalS += step;
+ }
+ myBField = bmap.getField(new BasicHep3Vector(rTmp.y(), rTmp.z(), rTmp.x())).y(); // note weird co ordinates for track vs det
+ _trajectory = getTrajectory(pTmp, rTmp, q, myBField);
+ // Go with finer granularity in the end
+ rTmp = _trajectory.getPointAtDistance(0);
+ xtmp = rTmp.x();
+ pTmp = VecOp.mult(pTot, _trajectory.getUnitTangentAtLength(step));
+ pXTmp = pTmp.x();
+ step = step / 10.0;
+
+ while (Math.signum(step) * xtmp < Math.signum(step) * zFinal && Math.signum(pXOrig * pXTmp) > 0) {
+ if (debugOk) {
+ System.out.println("New step in Fringe Field");
+ System.out.println("rTmp = " + rTmp.toString());
+ System.out.println("pTmp = " + pTmp.toString());
+ System.out.println("OriginalHelix pos = " + HelixUtils.PointOnHelix(helix, totalS));
+ System.out.println("OriginalHelix Momentum = " + VecOp.mult(pTot, HelixUtils.Direction(helix, totalS)));
+ }
+
+ myBField = bmap.getField(new BasicHep3Vector(rTmp.y(), rTmp.z(), rTmp.x())).y(); // note weird co ordinates for track vs det
+
+ if (debugOk)
+ System.out.println("rTmp.x() = " + rTmp.x() + " field = " + myBField);
+ _trajectory = getTrajectory(pTmp, rTmp, q, myBField);
+ rTmp = _trajectory.getPointAtDistance(step);
+ pTmp = VecOp.mult(pTot, _trajectory.getUnitTangentAtLength(step));
+ pXTmp = pTmp.x();
+ xtmp = rTmp.x();
+ if (debugOk) {
+ System.out.println("############## done... #############");
+ System.out.println("\n");
+ }
+ totalS += step;
+ }
+
+ // ok, done with field.
+ //store the helical track parameters at rTmp
+ double pars[] = {Math.sqrt(rTmp.x() * rTmp.x() + rTmp.y() * rTmp.y()),
+ calculatePhi(pTmp.x(), pTmp.y()),
+ calculateCurvature(pTmp.magnitude(), q, myBField),
+ calculateLambda(pTmp.z(), pTmp.magnitude()),
+ rTmp.z()};
+ Hep3Vector pointInTrking = new BasicHep3Vector(rTmp.x(), rTmp.y(), rTmp.z());
+
+ if (debugOk)
+ System.out.println("Position xfinal (tracking) : x = " + zFinal + "; y = " + pointInTrking.y() + "; z = " + pointInTrking.z());
+ BaseTrackState ts = new BaseTrackState(pars, myBField);
+ ts.setReferencePoint(pointInTrking.v());
+ ts.setLocation(TrackState.AtCalorimeter);
+ return ts;
+ }
+
+ public static double calculatePhi(double x, double y, double xc, double yc, double sign) {
+ return Math.atan2(y - yc, x - xc) - sign * Math.PI / 2;
+ }
+
+ public static double calculatePhi(double px, double py) {
+ return Math.atan2(py, px);
+ }
+
+ public static double calculateLambda(double pz, double p) {
+ return Math.atan2(pz, p);
+ }
+
+ public static double calculateCurvature(double p, double q, double B) {
+ return q * B / p;
+ }
+
+ public static Trajectory getTrajectory(Hep3Vector p0, org.lcsim.spacegeom.SpacePoint r0, double q, double B) {
+ SpaceVector p = new CartesianVector(p0.v());
+ double phi = Math.atan2(p.y(), p.x());
+ double lambda = Math.atan2(p.z(), p.rxy());
+ double field = B * fieldConversion;
+
+ if (q != 0 && field != 0) {
+ double radius = p.rxy() / (q * field);
+ return new Helix(r0, radius, phi, lambda);
+ } else
+ return new Line(r0, phi, lambda);
+ }
+
+}
