Skip to content

Commit

Permalink
Merge pull request #996 from alspellm/master
Browse files Browse the repository at this point in the history
Add track state at target
  • Loading branch information
alspellm authored Jul 25, 2023
2 parents 2e3cf9e + a83f53d commit 7871a8e
Show file tree
Hide file tree
Showing 6 changed files with 177 additions and 20 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@
<driver name="RfFitter"/>

<!-- Ecal reconstruction drivers -->
<!--<driver name="EcalRunningPedestal"/>-->
<driver name="EcalRunningPedestal"/>
<driver name="EcalRawConverter" />
<driver name="HitTimeSmear"/>
<driver name="ReconClusterer" />
Expand All @@ -21,7 +21,7 @@
<driver name="RawTrackerHitSensorSetup"/>
<driver name="RawTrackerHitFitterDriver" />
<driver name="TrackerHitDriver"/>
<driver name="HelicalTrackHitDriver"/>
<!--<driver name="HelicalTrackHitDriver"/>-->

<!-- <driver name="MergeTrackCollections"/>
<driver name="GBLRefitterDriver" />
Expand All @@ -44,9 +44,9 @@
<driver name="RfFitter" type="org.hps.evio.RfFitterDriver"/>

<!-- Ecal reconstruction drivers -->
<!--<driver name="EcalRunningPedestal" type="org.hps.recon.ecal.EcalRunningPedestalDriver">
<driver name="EcalRunningPedestal" type="org.hps.recon.ecal.EcalRunningPedestalDriver">
<logLevel>CONFIG</logLevel>
</driver>-->
</driver>

<!--<driver name="EcalRawConverter" type="org.hps.recon.ecal.EcalRawConverter2Driver">
</driver>-->
Expand Down Expand Up @@ -94,7 +94,7 @@
<debug>false</debug>
</driver>
<driver name="TrackerHitDriver" type="org.hps.recon.tracking.DataTrackerHitDriver">
<neighborDeltaT>8.0</neighborDeltaT>
<neighborDeltaT>24.0</neighborDeltaT>
</driver>
<driver name="HelicalTrackHitDriver" type="org.hps.recon.tracking.HelicalTrackHitDriver">
<debug>false</debug>
Expand All @@ -111,6 +111,9 @@
<!-- <siHitsLimit>50</siHitsLimit> -->
<!-- <seedCompThr>0.05</seedCompThr> -->
<!-- <addResiduals>true</addResiduals> -->
<numEvtPlots>40</numEvtPlots>
<targetPosition>-4.3</targetPosition>
<addTrackStateAtTarget>true</addTrackStateAtTarge>
<verbose>false</verbose>
</driver>

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -83,7 +83,7 @@
<debug>false</debug>
</driver>
<driver name="TrackerHitDriver" type="org.hps.recon.tracking.DataTrackerHitDriver">
<neighborDeltaT>8.0</neighborDeltaT>
<neighborDeltaT>24.0</neighborDeltaT>
</driver>
<driver name="HelicalTrackHitDriver" type="org.hps.recon.tracking.HelicalTrackHitDriver">
<debug>false</debug>
Expand Down Expand Up @@ -126,6 +126,8 @@
<!-- <addResiduals>true</addResiduals> -->
<numEvtPlots>40</numEvtPlots>
<verbose>false</verbose>
<targetPosition>-4.3</targetPosition>
<addTrackStateAtTarget>true</addTrackStateAtTarget>
</driver>

<driver name="LCIOWriter" type="org.lcsim.util.loop.LCIODriver">
Expand Down
24 changes: 23 additions & 1 deletion tracking/src/main/java/org/hps/recon/tracking/TrackData.java
Original file line number Diff line number Diff line change
Expand Up @@ -18,11 +18,14 @@ public class TrackData implements GenericObject {
public static final int L1_ISOLATION_INDEX = 0;
public static final int L2_ISOLATION_INDEX = 1;
public static final int N_ISOLATIONS = 14; //Default
public static final int N_TRACK_PARAMS = 4;
public static final int N_TRACK_PARAMS = 7;
public static final int TRACK_TIME_INDEX = 0;
public static final int PX_INDEX = 1;
public static final int PY_INDEX = 2;
public static final int PZ_INDEX = 3;
public static final int ORIGIN_BFY_INDEX = 4; //BFieldY at Origin
public static final int TARGET_BFY_INDEX = 5; //BFieldY at Target TrackState
public static final int ECAL_BFY_INDEX = 6; //BFieldY at ECal TrackState
public static final int TRACK_VOLUME_INDEX = 0;
public static final String TRACK_DATA_COLLECTION = "TrackData";
public static final String TRACK_DATA_RELATION_COLLECTION = "TrackDataRelations";
Expand Down Expand Up @@ -81,6 +84,25 @@ public TrackData(int trackVolume, float trackTime, double[] isolations, float[]
this.floats = new float[]{trackTime,momentum[0],momentum[1],momentum[2]};
this.ints = new int[]{trackVolume};
}

/**
* Constructor
*
* @param trackVolume : SVT volume associated with the track
* @param trackTime : The track time
* @param isolations : an array of doubles containing isolations for every
* sensor layer
* @param momentum : an array of floats containing track momentum in the form (px,py,pz)
* @param origin_bfieldY : value of BfieldY at origin ref
* @param target_bfieldY : value of BfieldY at Target Track State
* @param ecal_bfieldY : value of BfieldY at ECal Track State
*/
public TrackData(int trackVolume, float trackTime, double[] isolations, float[] momentum, float origin_bfieldY, float target_bfieldY, float ecal_bfieldY) {

this.doubles = isolations;
this.floats = new float[]{trackTime,momentum[0],momentum[1],momentum[2],origin_bfieldY,target_bfieldY,ecal_bfieldY};
this.ints = new int[]{trackVolume};
}

/**
* Get isolation value for the hit in the given sensor layer.
Expand Down
58 changes: 58 additions & 0 deletions tracking/src/main/java/org/hps/recon/tracking/TrackUtils.java
Original file line number Diff line number Diff line change
Expand Up @@ -755,6 +755,60 @@ public static BaseTrackState getTrackExtrapAtHodoRK(Track trk, FieldMap fM, int
return getTrackExtrapAtHodoRK(trk, fM, 0, hodoLayer);
}

public static BaseTrackState getTrackExtrapAtTargetRK(Track track, double target_z, double[] beamPosition, FieldMap fM, double stepSize) {

TrackState ts = track.getTrackStates().get(0);

//if track passed to extrapolateHelixToXPlane, uses first track state
//by default, else if trackstate is passed, uses trackstate params.
//Forms HTF, projects to location, then returns point on helix
Hep3Vector startPos = extrapolateHelixToXPlane(ts, 0.0);
Hep3Vector startPosTrans = CoordinateTransformations.transformVectorToDetector(startPos);
double distanceZ = target_z;
double charge = -1.0 * Math.signum(getR(ts));

//extrapolateTrackUsingFieldMapRK gets HTF of input track/trackstate
//if track is passed, defaults to first track state
//if trackstate is passed, uses trackstate params
org.hps.util.Pair<Hep3Vector, Hep3Vector> RKresults = extrapolateTrackUsingFieldMapRK(ts, startPosTrans, distanceZ, stepSize, fM);
double bFieldY = fM.getField(RKresults.getFirstElement()).y();
Hep3Vector posTrans = CoordinateTransformations.transformVectorToTracking(RKresults.getFirstElement());
Hep3Vector momTrans = CoordinateTransformations.transformVectorToTracking(RKresults.getSecondElement());

Hep3Vector finalPos = posTrans;
if (RKresults.getFirstElement().z() != target_z) {
Hep3Vector mom = RKresults.getSecondElement();
double dz = target_z - RKresults.getFirstElement().z();
double dy = dz * mom.y() / mom.z();
double dx = dz * mom.x() / mom.z();
Hep3Vector dPos = new BasicHep3Vector(dx, dy, dz);
finalPos = CoordinateTransformations.transformVectorToTracking(VecOp.add(dPos, RKresults.getFirstElement()));
}
bFieldY = fM.getField(CoordinateTransformations.transformVectorToDetector(finalPos)).y();
//params are calculated with respect to ref = {trackX, trackY, z=0)
double[] params = getParametersFromPointAndMomentum(finalPos, momTrans, (int) charge, bFieldY);
BaseTrackState bts = new BaseTrackState(params, bFieldY);
//reference point is set to track position in X Y Z
bts.setReferencePoint(finalPos.v());
//Define new reference point, to which track parameters are calc wrt
double[] newRef = {target_z, beamPosition[0], beamPosition[1]};
params = getParametersAtNewRefPoint(newRef, bts);
bts.setParameters(params, bFieldY);
//Reference point records final position of track.
//This does not hold the reference point to which the track params are
//calculated from
bts.setReferencePoint(finalPos.v());
bts.setLocation(TrackState.LastLocation);

//Get covariance matrix at target. This does not use RK extrap, but
//simple change of reference point
SymmetricMatrix originCovMatrix = new SymmetricMatrix(5, ts.getCovMatrix(),true);
SymmetricMatrix covtrans = getCovarianceAtNewRefPoint(newRef, ts.getReferencePoint(), ts.getParameters(), originCovMatrix);
bts.setCovMatrix(covtrans.asPackedArray(true));

return bts;
}

/**
* Extrapolate track to given position. For backwards compatibility.
*
Expand Down Expand Up @@ -1869,6 +1923,10 @@ public static TrackState getTrackStateAtLocation(Track trk, int location) {
return null;
}

public static TrackState getTrackStateAtTarget(Track trk){
return getTrackStateAtLocation(trk, TrackState.LastLocation);
}

public static TrackState getTrackStateAtECal(Track trk) {
return getTrackStateAtLocation(trk, TrackState.AtCalorimeter);
}
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -79,6 +79,9 @@ public class KalmanInterface {
private int maxHits;
private int nBigEvents;
private int eventNumber;
private static double target_pos = -999.9;
private static boolean addTrackStateAtTarget = false;
private double[] beamPosition = null;

private static final boolean debug = false;
private static final double SVTcenter = 505.57;
Expand All @@ -97,6 +100,18 @@ public void setSiHitsLimit(int limit) {
public int getSiHitsLimit() {
return _siHitsLimit;
}

public void setTargetPosition(double target_pos){
this.target_pos = target_pos;
}

public void setAddTrackStateAtTarget(boolean input){
this.addTrackStateAtTarget = input;
}

public void setBeamPosition(double[] beamPosition){
this.beamPosition = beamPosition;
}

// Get the HPS tracker hit corresponding to a Kalman hit
public TrackerHit getHpsHit(Measurement km) {
Expand Down Expand Up @@ -520,7 +535,7 @@ public List<GBLStripClusterData> createGBLStripClusterData(KalTrack kT) {
public PropagatedTrackState propagateTrackState(TrackState stateHPS, double [] location, double [] direction) {
return new PropagatedTrackState(stateHPS, location, direction, detPlanes, fM);
}

// Create an HPS track from a Kalman track
public BaseTrack createTrack(KalTrack kT, boolean storeTrackStates) {
if (kT.SiteList == null) {
Expand Down Expand Up @@ -596,9 +611,18 @@ public BaseTrack createTrack(KalTrack kT, boolean storeTrackStates) {
}

// Extrapolate to the ECAL and make a new trackState there.

BaseTrackState ts_ecal = TrackUtils.getTrackExtrapAtEcalRK(newTrack, fM, runNumber);
BaseTrackState ts_ecal = new BaseTrackState();
ts_ecal = TrackUtils.getTrackExtrapAtEcalRK(newTrack, fM, runNumber);
newTrack.getTrackStates().add(ts_ecal);

// Extrapolate to Target and make a new trackState there.
BaseTrackState ts_target = new BaseTrackState();
if (target_pos != -999.9 && addTrackStateAtTarget){
ts_target = TrackUtils.getTrackExtrapAtTargetRK(newTrack, target_pos, beamPosition, fM, 0);
if (ts_target != null){
newTrack.getTrackStates().add(ts_target);
}
}

// other track properties
newTrack.setChisq(kT.chi2);
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@
import org.hps.recon.tracking.MaterialSupervisor.ScatteringDetectorVolume;
import org.hps.recon.tracking.MaterialSupervisor.SiStripPlane;
import org.hps.recon.tracking.gbl.GBLStripClusterData;
import org.hps.recon.tracking.TrackUtils;
import org.hps.util.Pair;
import org.lcsim.detector.tracker.silicon.HpsSiSensor;
import org.lcsim.event.EventHeader;
Expand Down Expand Up @@ -56,6 +57,8 @@ public class KalmanPatRecDriver extends Driver {
private KalmanParams kPar;
private KalmanPatRecPlots kPlot;
private static Logger logger;
private static double target_pos = -999.9;
private static boolean addTrackStateAtTarget = false;

// Parameters for the Kalman pattern recognition that can be set by the user in the steering file:
private ArrayList<String> strategies; // List of seed strategies for both top and bottom trackers, from steering
Expand Down Expand Up @@ -130,6 +133,14 @@ public void setSiHitsLimit(int input) {
public void setAddResiduals(boolean input) {
addResiduals = input;
}

public void setTargetPosition(double target_pos){
this.target_pos = target_pos;
}

public void setAddTrackStateAtTarget(boolean input){
this.addTrackStateAtTarget = input;
}

@Override
public void detectorChanged(Detector det) {
Expand Down Expand Up @@ -248,17 +259,28 @@ public void detectorChanged(Detector det) {

// Setup optional usage of beam positions from database.
final DatabaseConditionsManager mgr = DatabaseConditionsManager.getInstance();
if (useBeamPositionConditions && mgr.hasConditionsRecord("beam_positions")) {
logger.config("Using Kalman beam position from the conditions database");
double[] beamPositionArr = {beamPositionX, beamPositionY, beamPositionZ};
if (mgr.hasConditionsRecord("beam_positions")){
BeamPositionCollection beamPositions =
mgr.getCachedConditions(BeamPositionCollection.class, "beam_positions").getCachedData();
BeamPosition beamPositionCond = beamPositions.get(0);
if (!useFixedVertexZPosition) kPar.setBeamSpotY(beamPositionCond.getPositionZ());
else logger.config("Using fixed Kalman beam Z position: " + kPar.beamSpot[1]);
kPar.setBeamSpotX(beamPositionCond.getPositionX()); // Includes a transformation to Kalman coordinates
kPar.setBeamSpotZ(-beamPositionCond.getPositionY());
} else {
logger.config("Using Kalman beam position from the steering file or default");
beamPositionArr[0] = beamPositionCond.getPositionX();
beamPositionArr[1] = beamPositionCond.getPositionY();
beamPositionArr[2] = beamPositionCond.getPositionZ();
System.out.println("beamPosition[0]: "+beamPositionArr[0]);
System.out.println("beamPosition[1]: "+beamPositionArr[1]);
System.out.println("beamPosition[2]: "+beamPositionArr[2]);
if (useBeamPositionConditions) {
logger.config("Using Kalman beam position from the conditions database");
if (!useFixedVertexZPosition) kPar.setBeamSpotY(beamPositionCond.getPositionZ());
else logger.config("Using fixed Kalman beam Z position: " + kPar.beamSpot[1]);
kPar.setBeamSpotX(beamPositionCond.getPositionX()); // Includes a transformation to Kalman coordinates
kPar.setBeamSpotZ(-beamPositionCond.getPositionY());
}
else {
logger.config("Using Kalman beam position from the steering file or default");
}

}
logger.config("Using Kalman beam position [ Z, X, Y ]: " + String.format("[ %f, %f, %f ]",
kPar.beamSpot[0], -kPar.beamSpot[2], kPar.beamSpot[1]) + " in HPS coordinates.");
Expand All @@ -268,6 +290,12 @@ public void detectorChanged(Detector det) {
kPar.print();

KI = new KalmanInterface(uniformB, kPar, fm);
//Track State at Target uses beam position as track param reference
KI.setBeamPosition(beamPositionArr);
if (target_pos != -999.9 && addTrackStateAtTarget) {
KI.setAddTrackStateAtTarget(addTrackStateAtTarget);
KI.setTargetPosition(target_pos);
}
KI.setSiHitsLimit(siHitsLimit);
KI.createSiModules(detPlanes);
decoder = det.getSubdetector("Tracker").getIDDecoder();
Expand Down Expand Up @@ -398,6 +426,7 @@ private ArrayList<KalTrack>[] prepareTrackCollections(EventHeader event, List<Tr
//double pt = Math.abs(1./ptInv_check);

outputFullTracks.add(KalmanTrackHPS);

List<GBLStripClusterData> clstrs = KI.createGBLStripClusterData(kTk);
if (verbose) {
for (GBLStripClusterData clstr : clstrs) {
Expand Down Expand Up @@ -431,9 +460,28 @@ private ArrayList<KalTrack>[] prepareTrackCollections(EventHeader event, List<Tr
momentum_f[0] = (float) momentum.x();
momentum_f[1] = (float) momentum.y();
momentum_f[2] = (float) momentum.z();


//Get Bfield at origin
double origin_bFieldY = fm.getField((new BasicHep3Vector(0.0,0.0,0.0))).y();

//Get Bfield at target
double target_bFieldY = -999.9;
if (TrackUtils.getTrackStateAtTarget(KalmanTrackHPS) != null)
{
Hep3Vector target_pos = new BasicHep3Vector(TrackUtils.getTrackStateAtTarget(KalmanTrackHPS).getReferencePoint());
target_bFieldY = fm.getField(CoordinateTransformations.transformVectorToDetector(target_pos)).y();

}
//Get Bfield at ecal
double ecal_bFieldY = -999.9;
if (TrackUtils.getTrackStateAtECal(KalmanTrackHPS) != null)
{
Hep3Vector ecal_pos = new BasicHep3Vector(TrackUtils.getTrackStateAtECal(KalmanTrackHPS).getReferencePoint());
ecal_bFieldY = fm.getField(CoordinateTransformations.transformVectorToDetector(ecal_pos)).y();
}

//Add the Track Data
TrackData KFtrackData = new TrackData(trackerVolume, (float) kTk.getTime(), qualityArray, momentum_f);
TrackData KFtrackData = new TrackData(trackerVolume, (float) kTk.getTime(), qualityArray, momentum_f, (float) origin_bFieldY, (float) target_bFieldY, (float) ecal_bFieldY);
trackDataCollection.add(KFtrackData);
trackDataRelations.add(new BaseLCRelation(KFtrackData, KalmanTrackHPS));

Expand Down

0 comments on commit 7871a8e

Please sign in to comment.