Skip to content
This repository has been archived by the owner on May 18, 2023. It is now read-only.

Commit

Permalink
Usage of sector for swimmer and Bfield estimate in TSC
Browse files Browse the repository at this point in the history
  • Loading branch information
@zieglerv
zieglerv committed May 31, 2018
1 parent 7c1c879 commit c3a428b
Show file tree
Hide file tree
Showing 12 changed files with 61 additions and 59 deletions.
5 changes: 3 additions & 2 deletions reconstruction/dc/src/main/java/org/jlab/rec/dc/banks/HitReader.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@
import cnuphys.snr.clas12.Clas12NoiseResult;
import org.jlab.detector.geant4.v2.DCGeant4Factory;
import org.jlab.rec.dc.Constants;
import org.jlab.rec.dc.trajectory.DCSwimmer;
import org.jlab.utils.groups.IndexedTable;

/**
Expand Down Expand Up @@ -215,10 +216,10 @@ public void fetch_DCHits(DataEvent event, Clas12NoiseAnalysis noiseAnalysis, Noi
passTimingCut=true;
if(region ==2) {
if(wire[i]>=56) {
if(T0Sub>timeCutMin && T0Sub<timeCutMax+timeCutLC*(float)(112-wire[i]/56))
if(T0Sub>timeCutMin && T0Sub<timeCutMax + timeCutLC*(float)(112-wire[i]/56) )
passTimingCut=true;
} else {
if(T0Sub>timeCutMin && T0Sub<timeCutMax+timeCutLC*(float)(56-wire[i]/56))
if(T0Sub>timeCutMin && T0Sub<timeCutMax + (200 + timeCutLC*(float)(56-wire[i]/56)) )
passTimingCut=true;
}
}
Expand Down
4 changes: 2 additions & 2 deletions reconstruction/dc/src/main/java/org/jlab/rec/dc/cross/CrossListFinder.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -238,14 +238,14 @@ private void updateBFittedHits(Cross c, IndexedTable tab, DCGeant4Factory DcDete
for(int i =0; i<c.get_Segment1().size(); i++) {
Point3D ref =c.get_Segment1().get(i).getCrossDirIntersWire();
float[] result = new float[3];
swimmer.Bfield(ref.x(), ref.y(), ref.z(), result);
swimmer.Bfield(c.get_Sector(), ref.x(), ref.y(), ref.z(), result);
c.get_Segment1().get(i).setB(Math.sqrt(result[0]*result[0]+result[1]*result[1]+result[2]*result[2]) );

}
for(int i =0; i<c.get_Segment2().size(); i++) {
Point3D ref =c.get_Segment2().get(i).getCrossDirIntersWire();
float[] result = new float[3];
swimmer.Bfield(ref.x(), ref.y(), ref.z(), result);
swimmer.Bfield(c.get_Sector(), ref.x(), ref.y(), ref.z(), result);
c.get_Segment2().get(i).setB(Math.sqrt(result[0]*result[0]+result[1]*result[1]+result[2]*result[2]) );
}
if(tde!=null) {
Expand Down
30 changes: 15 additions & 15 deletions reconstruction/dc/src/main/java/org/jlab/rec/dc/track/TrackCandListFinder.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,34 +71,34 @@ public boolean PassNSuperlayerTracking(List<Cross> crossesInTrk, Track cand) {
return pass;
}

public double getHitBasedFitChi2ToCrosses(double x1, double y1, double z1, double x2, double y2, double z2, double x3, double y3, double z3,
public double getHitBasedFitChi2ToCrosses(int sector, double x1, double y1, double z1, double x2, double y2, double z2, double x3, double y3, double z3,
double p, int q, double x, double y, double z, double tanThX, double tanThY, DCSwimmer dcSwim) {
double pz = p / Math.sqrt(tanThX*tanThX + tanThY*tanThY + 1);

dcSwim.SetSwimParameters(x,y,z,
-pz*tanThX,-pz*tanThY,-pz,
-q);
double chi2 = 0; // assume err =1 on points
double[] R = dcSwim.SwimToPlane(z3);
double[] R = dcSwim.SwimToPlane(sector, z3);
chi2+= (R[0] - x3)*(R[0] - x3) + (R[1] - y3)*(R[1] - y3);
dcSwim.SetSwimParameters(R[0],R[1],R[2],
R[3],R[4],R[5],
-q);
R = dcSwim.SwimToPlane(z2);
R = dcSwim.SwimToPlane(sector, z2);
dcSwim.SetSwimParameters(R[0],R[1],R[2],
R[3],R[4],R[5],
-q);
chi2+= (R[0] - x2)*(R[0] - x2) + (R[1] - y2)*(R[1] - y2);
dcSwim.SetSwimParameters(R[0],R[1],R[2],
R[3],R[4],R[5],
-q);
R = dcSwim.SwimToPlane(z1);
R = dcSwim.SwimToPlane(sector, z1);
chi2+= (R[0] - x1)*(R[0] - x1) + (R[1] - y1)*(R[1] - y1);

return chi2;
}

private double[] getTrackInitFit(double x1, double y1, double z1, double x2, double y2, double z2, double x3, double y3, double z3,
private double[] getTrackInitFit(int sector, double x1, double y1, double z1, double x2, double y2, double z2, double x3, double y3, double z3,
double ux, double uy, double uz, double thX, double thY,
double theta1, double theta3,
double iBdl, double TORSCALE, DCSwimmer dcSwim) {
Expand All @@ -120,13 +120,13 @@ private double[] getTrackInitFit(double x1, double y1, double z1, double x2, dou
int q= calcInitTrkQ(theta1, theta3,TORSCALE);

dcSwim.SetSwimParameters(x1,y1,z1,p_x, p_y, p_z, q);
double[] R = dcSwim.SwimToPlane(z2);
double[] R = dcSwim.SwimToPlane(sector, z2);
chi2+= (R[0] - x2)*(R[0] - x2) + (R[1] - y2)*(R[1] - y2);
intBdl+=R[7];
dcSwim.SetSwimParameters(R[0],R[1],R[2],
R[3],R[4],R[5],
q);
R = dcSwim.SwimToPlane(z3);
R = dcSwim.SwimToPlane(sector, z3);
chi2+= (R[0] - x3)*(R[0] - x3) + (R[1] - y3)*(R[1] - y3);
intBdl+=R[7];

Expand Down Expand Up @@ -255,7 +255,7 @@ public List<Track> getTrackCands(CrossList crossList, DCGeant4Factory DcDetector
double iBdl = traj.get_IntegralBdl();
double[] pars;

pars = getTrackInitFit(x1, y1, z1, x2, y2, z2, x3, y3, z3,
pars = getTrackInitFit(cand.get_Sector(), x1, y1, z1, x2, y2, z2, x3, y3, z3,
ux, uy, uz, thX, thY,
theta1s1, theta3s1,
traj.get_IntegralBdl(), TORSCALE, dcSwim);
Expand All @@ -267,7 +267,7 @@ public List<Track> getTrackCands(CrossList crossList, DCGeant4Factory DcDetector
iBdl = pars[1];
}

pars = getTrackInitFit(x1, y1, z1, x2, y2, z2, x3, y3, z3,
pars = getTrackInitFit(cand.get_Sector(), x1, y1, z1, x2, y2, z2, x3, y3, z3,
ux, uy, uz, thX, thY,
theta1s1, theta3s2,
traj.get_IntegralBdl(), TORSCALE, dcSwim);
Expand All @@ -279,7 +279,7 @@ public List<Track> getTrackCands(CrossList crossList, DCGeant4Factory DcDetector
iBdl = pars[1];
}

pars = getTrackInitFit(x1, y1, z1, x2, y2, z2, x3, y3, z3,
pars = getTrackInitFit(cand.get_Sector(), x1, y1, z1, x2, y2, z2, x3, y3, z3,
ux, uy, uz, thX, thY,
theta1s2, theta3s1,
traj.get_IntegralBdl(), TORSCALE, dcSwim);
Expand All @@ -291,7 +291,7 @@ public List<Track> getTrackCands(CrossList crossList, DCGeant4Factory DcDetector
iBdl = pars[1];
}

pars = getTrackInitFit(x1, y1, z1, x2, y2, z2, x3, y3, z3,
pars = getTrackInitFit(cand.get_Sector(), x1, y1, z1, x2, y2, z2, x3, y3, z3,
ux, uy, uz, thX, thY,
theta1s2, theta3s2,
traj.get_IntegralBdl(), TORSCALE, dcSwim);
Expand Down Expand Up @@ -353,7 +353,7 @@ public List<Track> getTrackCands(CrossList crossList, DCGeant4Factory DcDetector
System.out.println("x "+ kFit.finalStateVec.x);
System.out.println("y "+ kFit.finalStateVec.y);
System.out.println("z "+ kFit.finalStateVec.z); */
double HBc2 = getHitBasedFitChi2ToCrosses( x1, y1, z1, x2, y2, z2, x3, y3, z3,
double HBc2 = getHitBasedFitChi2ToCrosses(cand.get_Sector(), x1, y1, z1, x2, y2, z2, x3, y3, z3,
1./Math.abs(kFit.finalStateVec.Q), (int)Math.signum(kFit.finalStateVec.Q),
kFit.finalStateVec.x, kFit.finalStateVec.y, kFit.finalStateVec.z, kFit.finalStateVec.tx, kFit.finalStateVec.ty, dcSwim);
//System.out.println(cand.get(0).get_Sector()+" HB fit to crosses c2 "+HBc2);
Expand Down Expand Up @@ -506,7 +506,7 @@ public void setTrackPars(Track cand, Trajectory traj, TrajectoryFinder trjFind,
cand.get_Q());

// swimming to a ref points outside of the last DC region
double[] VecAtTarOut = dcSwim.SwimToPlane(592);
double[] VecAtTarOut = dcSwim.SwimToPlane(cand.get_Sector(), 592);
double xOuter = VecAtTarOut[0];
double yOuter = VecAtTarOut[1];
double zOuter = VecAtTarOut[2];
Expand All @@ -527,7 +527,7 @@ public void setTrackPars(Track cand, Trajectory traj, TrajectoryFinder trjFind,
-cand.get_Q());

//swimming to a ref point upstream of the first DC region
double[] VecAtTarIn = dcSwim.SwimToPlane(180);
double[] VecAtTarIn = dcSwim.SwimToPlane(cand.get_Sector(), 180);

if(VecAtTarIn==null) {
cand.fit_Successful=false;
Expand All @@ -547,7 +547,7 @@ public void setTrackPars(Track cand, Trajectory traj, TrajectoryFinder trjFind,
double pzOr = -VecAtTarIn[5];

if(traj!=null && trjFind!=null) {
traj.set_Trajectory(trjFind.getStateVecsAlongTrajectory(xOr, yOr, pxOr/pzOr, pyOr/pzOr, cand.get_P(),cand.get_Q(), getDcDetector));
traj.set_Trajectory(trjFind.getStateVecsAlongTrajectory(cand.get_Sector(), xOr, yOr, pxOr/pzOr, pyOr/pzOr, cand.get_P(),cand.get_Q(), getDcDetector));
cand.set_Trajectory(traj.get_Trajectory());
}
//cand.set_Vtx0_TiltedCS(trakOrigTiltSec);
Expand Down
23 changes: 12 additions & 11 deletions reconstruction/dc/src/main/java/org/jlab/rec/dc/track/fit/KFitter.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,10 +24,10 @@ public class KFitter {
public double chi2kf = 0;
public int NDF = 0;
public int ConvStatus = 1;

public int sector;
public KFitter(Track trk, DCGeant4Factory DcDetector, boolean TimeBasedUsingHBtrack, DCSwimmer dcSwim) {
sv = new StateVecs(dcSwim);

sector = trk.get_Sector();
if(TimeBasedUsingHBtrack==true) {
this.initFromHB(trk, DcDetector);
} else {
Expand Down Expand Up @@ -65,10 +65,10 @@ public void runFitter() {
for (int i = 1; i <= totNumIter; i++) {
this.chi2kf = 0;
if (i > 1) {
sv.transport(sv.Z.length - 1, 0, sv.trackTraj.get(sv.Z.length - 1), sv.trackCov.get(sv.Z.length - 1)); //get new state vec at 1st measurement after propagating back from the last filtered state
sv.transport(sector, sv.Z.length - 1, 0, sv.trackTraj.get(sv.Z.length - 1), sv.trackCov.get(sv.Z.length - 1)); //get new state vec at 1st measurement after propagating back from the last filtered state
}
for (int k = 0; k < sv.Z.length - 1; k++) {
sv.transport(k, k + 1, sv.trackTraj.get(k), sv.trackCov.get(k));
sv.transport(sector, k, k + 1, sv.trackTraj.get(k), sv.trackCov.get(k));
//sv.trackTraj.add(k+1, sv.StateVec);
//sv.trackCov.add(k+1, sv.CovMat);
// System.out.println((k)+"] trans "+sv.trackTraj.get(k).x+","+sv.trackTraj.get(k).y+","+
Expand Down Expand Up @@ -105,7 +105,7 @@ public void runFitter() {
}
}
}
this.calcFinalChisq();
this.calcFinalChisq(sector);

}

Expand Down Expand Up @@ -187,24 +187,25 @@ private void filter(int k) {
}

@SuppressWarnings("unused")
private void smooth(int k) {
private void smooth(int sector, int k) {
this.chi2 = 0;
if (sv.trackTraj.get(k) != null && sv.trackCov.get(k).covMat != null) {
sv.transport(k, 0, sv.trackTraj.get(k), sv.trackCov.get(k));
sv.transport(sector, k, 0, sv.trackTraj.get(k), sv.trackCov.get(k));
for (int k1 = 0; k1 < k; k1++) {
sv.transport(k1, k1 + 1, sv.trackTraj.get(k1), sv.trackCov.get(k1));
sv.transport(sector, k1, k1 + 1, sv.trackTraj.get(k1), sv.trackCov.get(k1));
this.filter(k1 + 1);
}
}
}

private void calcFinalChisq() {
private void calcFinalChisq(int sector) {
int k = sv.Z.length - 1;
this.chi2 = 0;
if (sv.trackTraj.get(k) != null && sv.trackCov.get(k).covMat != null) {
sv.rinit(sv.Z[0], k);
sv.transport(sector, sv.Z.length - 1, 0, sv.trackTraj.get(sv.Z.length - 1), sv.trackCov.get(sv.Z.length - 1)); //get new state vec at 1st measurement after propagating back from the last filtered state

for (int k1 = 0; k1 < k; k1++) {
sv.transport(k1, k1 + 1, sv.trackTraj.get(k1), sv.trackCov.get(k1));
sv.transport(sector, k1, k1 + 1, sv.trackTraj.get(k1), sv.trackCov.get(k1));
double V = mv.measurements.get(k1 + 1).error;
double h = mv.h(new double[]{sv.trackTraj.get(k1 + 1).x, sv.trackTraj.get(k1 + 1).y},
(int) mv.measurements.get(k1 + 1).tilt, mv.measurements.get(k1 + 1).wireMaxSag, mv.measurements.get(k1 + 1).wireLen);
Expand Down
18 changes: 9 additions & 9 deletions reconstruction/dc/src/main/java/org/jlab/rec/dc/track/fit/StateVecs.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,7 @@ public StateVecs(DCSwimmer dcSwimmer) {
* @param iVec state vector at the initial state index
* @return state vector at the final state index
*/
public StateVec f(int i, int f, StateVec iVec) {
public StateVec f(int sector, int i, int f, StateVec iVec) {

double x = iVec.x;
double y = iVec.y;
Expand All @@ -68,7 +68,7 @@ public StateVec f(int i, int f, StateVec iVec) {
}

// propagate the state vector a next step
dcSwim.Bfield(x, y, z, bf);
dcSwim.Bfield(sector, x, y, z, bf);
A(tx, ty, bf[0], bf[1], bf[2], A);
// transport stateVec
x += tx * s + 0.5 * Q * speedLight * A[0] * s * s;
Expand Down Expand Up @@ -98,7 +98,7 @@ public StateVec f(int i, int f, StateVec iVec) {
* @param iVec state vector at the initial index
* @param covMat state covariance matrix at the initial index
*/
public void transport(int i, int f, StateVec iVec, CovMat covMat) { // s = signed step-size
public void transport(int sector, int i, int f, StateVec iVec, CovMat covMat) { // s = signed step-size
if(iVec==null)
return;
//StateVec iVec = trackTraj.get(i);
Expand All @@ -118,7 +118,7 @@ public void transport(int i, int f, StateVec iVec, CovMat covMat) { // s = signe
double Q = iVec.Q;

// B-field components at state vector coordinates
dcSwim.Bfield(x, y, Z[i], bf);
dcSwim.Bfield(sector, x, y, Z[i], bf);

// if (bfieldPoints.size() > 0) {
// double B = new Vector3D(bfieldPoints.get(bfieldPoints.size() - 1).Bx, bfieldPoints.get(bfieldPoints.size() - 1).By, bfieldPoints.get(bfieldPoints.size() - 1).Bz).mag();
Expand Down Expand Up @@ -158,7 +158,7 @@ public void transport(int i, int f, StateVec iVec, CovMat covMat) { // s = signe

//B bf = new B(i, z, x, y, tx, ty, s);
//bfieldPoints.add(bf);
dcSwim.Bfield(x, y, z, bf);
dcSwim.Bfield(sector, x, y, z, bf);

A(tx, ty, bf[0], bf[1], bf[2], A);
delA_delt(tx, ty, bf[0], bf[1], bf[2], dA);
Expand Down Expand Up @@ -372,7 +372,7 @@ private void setMass(int hypo, double mass) {
* @param z0 the value at which the state vector needs to be reinitialized
* @param kf the final state measurement index
*/
public void rinit(double z0, int kf) {
/*public void rinit(double z0, int kf) {
if (this.trackTraj.get(kf) != null) {
double x = this.trackTraj.get(kf).x;
double y = this.trackTraj.get(kf).y;
Expand All @@ -398,7 +398,7 @@ public void rinit(double z0, int kf) {
} else {
}
}

*/
/**
*
* @param trkcand the track candidate
Expand All @@ -412,7 +412,7 @@ public void init(Track trkcand, double z0, KFitter kf) {
trkcand.get_StateVecAtReg1MiddlePlane().tanThetaX(), trkcand.get_StateVecAtReg1MiddlePlane().tanThetaY(), trkcand.get_P(),
trkcand.get_Q());

double[] VecAtFirstMeasSite = dcSwim.SwimToPlane(z0);
double[] VecAtFirstMeasSite = dcSwim.SwimToPlane(trkcand.get_Sector(), z0);
StateVec initSV = new StateVec(0);
initSV.x = VecAtFirstMeasSite[0];
initSV.y = VecAtFirstMeasSite[1];
Expand Down Expand Up @@ -475,7 +475,7 @@ void initFromHB(Track trkcand, double z0, KFitter kf) {
dcSwim.SetSwimParameters(trkR1X.x(), trkR1X.y(), trkR1X.z(),
trkR1P.x(), trkR1P.y(), trkR1P.z(), trkcand.get_Q());

double[] VecAtFirstMeasSite = dcSwim.SwimToPlane(z0);
double[] VecAtFirstMeasSite = dcSwim.SwimToPlane(trkcand.get_Sector(), z0);
StateVec initSV = new StateVec(0);
initSV.x = VecAtFirstMeasSite[0];
initSV.y = VecAtFirstMeasSite[1];
Expand Down
2 changes: 1 addition & 1 deletion reconstruction/dc/src/main/java/org/jlab/rec/dc/track/fit/basefit/CallB.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,7 @@ public static void main(String[] args) {
DCSwimmer.setMagneticFieldsScales(0.0, -1.0, 0);
DCSwimmer swim2 = new DCSwimmer();
swim2.SetSwimParameters(0, 0, 0, 2*Math.sin(Math.toRadians(20.-25.)), 0, 2*Math.cos(Math.toRadians(20.-25)), -1);
double[]swimVal =swim2.SwimToPlane(500);
double[]swimVal =swim2.SwimToPlane(2,500);
for(int i = 0; i<swimVal.length; i++)
System.out.println("swimVal["+i+"]= "+swimVal[i]);
//CallB en = new CallB();
Expand Down
10 changes: 5 additions & 5 deletions reconstruction/dc/src/main/java/org/jlab/rec/dc/trajectory/DCSwimmer.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -169,7 +169,7 @@ public void SetSwimParameters(double x0, double y0, double z0, double px, double

}

public double[] SwimToPlane(double z_cm) {
public double[] SwimToPlane(int sector, double z_cm) {
double z = z_cm / 100; // the magfield method uses meters
double[] value = new double[8];
double accuracy = 20e-6; //20 microns
Expand All @@ -183,7 +183,7 @@ public double[] SwimToPlane(double z_cm) {
double hdata[] = new double[3];

try {
traj = swimmer.swim(_charge, _x0, _y0, _z0, _pTot,
traj = swimmer.sectorSwim(sector, _charge, _x0, _y0, _z0, _pTot,
_theta, _phi, z, accuracy, _rMax,
_maxPathLength, stepSize, Swimmer.CLAS_Tolerance, hdata);

Expand Down Expand Up @@ -495,9 +495,9 @@ public double[] SwimToPlaneBoundary(double d_cm, Vector3D n, int dir) {
}


public void Bfield(double x_cm, double y_cm, double z_cm, float[] result) {

rprob.field((float) x_cm, (float) y_cm, (float) z_cm, result);
public void Bfield(int sector, double x_cm, double y_cm, double z_cm, float[] result) {
rprob.field(sector, (float) x_cm, (float) y_cm, (float) z_cm, result);
//rcompositeField.field((float) x_cm, (float) y_cm, (float) z_cm, result);
result[0] =result[0]/10;
result[1] =result[1]/10;
Expand Down
Loading

0 comments on commit c3a428b

Please sign in to comment.