Skip to content

Transformation for covaraince matrix of DC TB tracks at vertex from local to global #1059

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
tongtongcao wants to merge 1 commit into
base: development
Choose a base branch
from
Open

Transformation for covaraince matrix of DC TB tracks at vertex from local to global #1059

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/utilities/MatrixOps.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -128,7 +128,7 @@ public double[][] MatrixTranspose(Object obj1) {

for(int i = 0; i < r; i++) {
for (int j = 0; j < c; j++) {
result[i][j] = arr1[j][i];
result[j][i] = arr1[i][j];
}
}
arr1 = null;
Expand Down
45 changes: 45 additions & 0 deletions etc/bankdefs/hipo4/dc.json
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -559,6 +559,51 @@
{"name":"C55", "type":"F", "info":"C55 covariance matrix element at last superlayer used in the fit"}
]
},
{
"name": "TimeBasedTrkg::TBCovMatLab",
"group": 20600,
"item" : 38,
"info": "reconstructed track covariance matrix in lab frame",
"entries": [
{"name":"id", "type":"S", "info":"id of the track"},
{"name":"C11", "type":"F", "info":"var(x) at vertex in lab frame"},
{"name":"C12", "type":"F", "info":"cov(x, y) at vertex in lab frame"},
{"name":"C13", "type":"F", "info":"cov(x, z) at vertex in lab frame"},
{"name":"C14", "type":"F", "info":"cov(x, p_theta) at vertex in lab frame"},
{"name":"C15", "type":"F", "info":"cov(x, p_phi) at vertex in lab frame"},
{"name":"C16", "type":"F", "info":"cov(x, p) at vertex in lab frame"},
{"name":"C21", "type":"F", "info":"cov(y, x) at vertex in lab frame"},
{"name":"C22", "type":"F", "info":"var(y) at vertex in lab frame"},
{"name":"C23", "type":"F", "info":"cov(y, z) at vertex in lab frame"},
{"name":"C24", "type":"F", "info":"cov(y, p_theta) at vertex in lab frame"},
{"name":"C25", "type":"F", "info":"cov(y, p_phi) at vertex in lab frame"},
{"name":"C26", "type":"F", "info":"cov(y, p) at vertex in lab frame"},
{"name":"C31", "type":"F", "info":"cov(z, x) at vertex in lab frame"},
{"name":"C32", "type":"F", "info":"cov(z, y) at vertex in lab frame"},
{"name":"C33", "type":"F", "info":"var(z) at vertex in lab frame"},
{"name":"C34", "type":"F", "info":"cov(z, p_theta) at vertex in lab frame"},
{"name":"C35", "type":"F", "info":"cov(z, p_phi) at vertex in lab frame"},
{"name":"C36", "type":"F", "info":"cov(z, p) at vertex in lab frame"},
{"name":"C41", "type":"F", "info":"cov(p_theta, x) at vertex in lab frame"},
{"name":"C42", "type":"F", "info":"cov(p_theta, y) at vertex in lab frame"},
{"name":"C43", "type":"F", "info":"cov(p_theta, z) at vertex in lab frame"},
{"name":"C44", "type":"F", "info":"var(p_theta) at vertex in lab frame"},
{"name":"C45", "type":"F", "info":"cov(p_theta, p_phi) at vertex in lab frame"},
{"name":"C46", "type":"F", "info":"cov(p_theta, p) at vertex in lab frame"},
{"name":"C51", "type":"F", "info":"cov(p_phi, x) at vertex in lab frame"},
{"name":"C52", "type":"F", "info":"cov(p_phi, y) at vertex in lab frame"},
{"name":"C53", "type":"F", "info":"cov(p_phi, z) at vertex in lab frame"},
{"name":"C54", "type":"F", "info":"cov(p_phi, p_theta) at vertex in lab frame"},
{"name":"C55", "type":"F", "info":"var(p_phi) at vertex in lab frame"},
{"name":"C56", "type":"F", "info":"cov(p_phi, p) at vertex in lab frame"},
{"name":"C61", "type":"F", "info":"cov(p, x) at vertex in lab frame"},
{"name":"C62", "type":"F", "info":"cov(p, y) at vertex in lab frame"},
{"name":"C63", "type":"F", "info":"cov(p, z) at vertex in lab frame"},
{"name":"C64", "type":"F", "info":"cov(p, p_theta) at vertex in lab frame"},
{"name":"C65", "type":"F", "info":"cov(p, p_phi) at vertex in lab frame"},
{"name":"C66", "type":"F", "info":"var(p) at vertex in lab frame"}
]
},
{
"name": "TimeBasedTrkg::Trajectory",
"group": 20600,
Expand Down
45 changes: 45 additions & 0 deletions etc/bankdefs/hipo4/dcnn.json
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -420,6 +420,51 @@
{"name":"C55", "type":"F", "info":"C55 covariance matrix element at last superlayer used in the fit"}
]
},
{
"name": "TimeBasedTrkg::AICovMatLab",
"group": 20600,
"item" : 78,
"info": "reconstructed track covariance matrix in lab frame",
"entries": [
{"name":"id", "type":"S", "info":"id of the track"},
{"name":"C11", "type":"F", "info":"var(x) at vertex in lab frame"},
{"name":"C12", "type":"F", "info":"cov(x, y) at vertex in lab frame"},
{"name":"C13", "type":"F", "info":"cov(x, z) at vertex in lab frame"},
{"name":"C14", "type":"F", "info":"cov(x, p_theta) at vertex in lab frame"},
{"name":"C15", "type":"F", "info":"cov(x, p_phi) at vertex in lab frame"},
{"name":"C16", "type":"F", "info":"cov(x, p) at vertex in lab frame"},
{"name":"C21", "type":"F", "info":"cov(y, x) at vertex in lab frame"},
{"name":"C22", "type":"F", "info":"var(y) at vertex in lab frame"},
{"name":"C23", "type":"F", "info":"cov(y, z) at vertex in lab frame"},
{"name":"C24", "type":"F", "info":"cov(y, p_theta) at vertex in lab frame"},
{"name":"C25", "type":"F", "info":"cov(y, p_phi) at vertex in lab frame"},
{"name":"C26", "type":"F", "info":"cov(y, p) at vertex in lab frame"},
{"name":"C31", "type":"F", "info":"cov(z, x) at vertex in lab frame"},
{"name":"C32", "type":"F", "info":"cov(z, y) at vertex in lab frame"},
{"name":"C33", "type":"F", "info":"var(z) at vertex in lab frame"},
{"name":"C34", "type":"F", "info":"cov(z, p_theta) at vertex in lab frame"},
{"name":"C35", "type":"F", "info":"cov(z, p_phi) at vertex in lab frame"},
{"name":"C36", "type":"F", "info":"cov(z, p) at vertex in lab frame"},
{"name":"C41", "type":"F", "info":"cov(p_theta, x) at vertex in lab frame"},
{"name":"C42", "type":"F", "info":"cov(p_theta, y) at vertex in lab frame"},
{"name":"C43", "type":"F", "info":"cov(p_theta, z) at vertex in lab frame"},
{"name":"C44", "type":"F", "info":"var(p_theta) at vertex in lab frame"},
{"name":"C45", "type":"F", "info":"cov(p_theta, p_phi) at vertex in lab frame"},
{"name":"C46", "type":"F", "info":"cov(p_theta, p) at vertex in lab frame"},
{"name":"C51", "type":"F", "info":"cov(p_phi, x) at vertex in lab frame"},
{"name":"C52", "type":"F", "info":"cov(p_phi, y) at vertex in lab frame"},
{"name":"C53", "type":"F", "info":"cov(p_phi, z) at vertex in lab frame"},
{"name":"C54", "type":"F", "info":"cov(p_phi, p_theta) at vertex in lab frame"},
{"name":"C55", "type":"F", "info":"var(p_phi) at vertex in lab frame"},
{"name":"C56", "type":"F", "info":"cov(p_phi, p) at vertex in lab frame"},
{"name":"C61", "type":"F", "info":"cov(p, x) at vertex in lab frame"},
{"name":"C62", "type":"F", "info":"cov(p, y) at vertex in lab frame"},
{"name":"C63", "type":"F", "info":"cov(p, z) at vertex in lab frame"},
{"name":"C64", "type":"F", "info":"cov(p, p_theta) at vertex in lab frame"},
{"name":"C65", "type":"F", "info":"cov(p, p_phi) at vertex in lab frame"},
{"name":"C66", "type":"F", "info":"var(p) at vertex in lab frame"}
]
},
{
"name": "TimeBasedTrkg::AITrajectory",
"group": 20600,
Expand Down
4 changes: 4 additions & 0 deletions reconstruction/dc/src/main/java/org/jlab/rec/dc/banks/Banks.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -138,6 +138,10 @@ public String getCovmatBank() {
return this.getOutputBank("CovMat");
}

public String getCovmatLabBank() {
return this.getOutputBank("CovMatLab");
}

public String getRecEventBank() {
return this.getRecBank("Event");
}
Expand Down
66 changes: 65 additions & 1 deletion reconstruction/dc/src/main/java/org/jlab/rec/dc/banks/RecoBankWriter.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -490,6 +490,69 @@ private DataBank fillTrackCovMatBank(DataEvent event, List<Track> candlist) {
//bank.show();
return bank;
}

/**
*
* @param event hipo event
* @param candlist tracks
* @return covariance matrix for momentum and vertex in lab frame
*/
private DataBank fillTrackCovMatLabBank(DataEvent event, List<Track> candlist) {

DataBank bank = event.createBank(bankNames.getCovmatLabBank(), candlist.size());

for (int i = 0; i < candlist.size(); i++) {
bank.setShort("id", i, (short) candlist.get(i).get_Id());
if(candlist.get(i).get_CMInLab()!=null) {
double[][] CM = candlist.get(i).get_CMInLab();
bank.setFloat("C11", i, (float) CM[0][0]);
bank.setFloat("C12", i, (float) CM[0][1]);
bank.setFloat("C13", i, (float) CM[0][2]);
bank.setFloat("C14", i, (float) CM[0][3]);
bank.setFloat("C15", i, (float) CM[0][4]);
bank.setFloat("C16", i, (float) CM[0][5]);

bank.setFloat("C21", i, (float) CM[1][0]);
bank.setFloat("C22", i, (float) CM[1][1]);
bank.setFloat("C23", i, (float) CM[1][2]);
bank.setFloat("C24", i, (float) CM[1][3]);
bank.setFloat("C25", i, (float) CM[1][4]);
bank.setFloat("C26", i, (float) CM[1][5]);

bank.setFloat("C31", i, (float) CM[2][0]);
bank.setFloat("C32", i, (float) CM[2][1]);
bank.setFloat("C33", i, (float) CM[2][2]);
bank.setFloat("C34", i, (float) CM[2][3]);
bank.setFloat("C35", i, (float) CM[2][4]);
bank.setFloat("C36", i, (float) CM[2][5]);

bank.setFloat("C41", i, (float) CM[3][0]);
bank.setFloat("C42", i, (float) CM[3][1]);
bank.setFloat("C43", i, (float) CM[3][2]);
bank.setFloat("C44", i, (float) CM[3][3]);
bank.setFloat("C45", i, (float) CM[3][4]);
bank.setFloat("C46", i, (float) CM[3][5]);

bank.setFloat("C51", i, (float) CM[4][0]);
bank.setFloat("C52", i, (float) CM[4][1]);
bank.setFloat("C53", i, (float) CM[4][2]);
bank.setFloat("C54", i, (float) CM[4][3]);
bank.setFloat("C55", i, (float) CM[4][4]);
bank.setFloat("C56", i, (float) CM[4][5]);

bank.setFloat("C61", i, (float) CM[5][0]);
bank.setFloat("C62", i, (float) CM[5][1]);
bank.setFloat("C63", i, (float) CM[5][2]);
bank.setFloat("C64", i, (float) CM[5][3]);
bank.setFloat("C65", i, (float) CM[5][4]);
bank.setFloat("C66", i, (float) CM[5][5]);

}
}
//bank.show();
return bank;
}

/**
*
* @param event the EvioEvent
Expand Down Expand Up @@ -939,7 +1002,8 @@ public void fillAllTBBanks(DataEvent event, List<FittedHit> fhits, List<FittedCl
this.fillTBCrossesBank(event, crosses),
this.fillTBTracksBank(event, trkcands),
this.fillTrajectoryBank(event, trkcands),
this.fillTrackCovMatBank(event, trkcands)
this.fillTrackCovMatBank(event, trkcands),
this.fillTrackCovMatLabBank(event, trkcands)
);
}
if (crosses != null && trkcands == null) {
Expand Down
87 changes: 87 additions & 0 deletions reconstruction/dc/src/main/java/org/jlab/rec/dc/track/Track.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,8 @@
import org.jlab.rec.dc.segment.Segment;
import org.jlab.rec.dc.trajectory.StateVec;
import org.jlab.rec.dc.trajectory.Trajectory;
import org.jlab.clas.tracking.utilities.MatrixOps;
import org.jlab.clas.tracking.utilities.MatrixOps.Libr;

/**
* A class representing track candidates in the DC. A track has a trajectory represented by an ensemble of geometrical state vectors along its path,
Expand Down Expand Up @@ -42,6 +44,7 @@ public void setFinalStateVec(StateVec finalStateVec) {
private int _Q;
private double _P;
private Matrix _CovMat;
private double[][] _CMInLab = new double[6][6];

private Point3D _Region3CrossPoint;
private Point3D _Region3CrossDir;
Expand Down Expand Up @@ -396,6 +399,90 @@ public void set_CovMat(Matrix _CovMat) {
this._CovMat = _CovMat;
}

/**
*
* @return covariance matrix in lab frame
*/
public double[][] get_CMInLab() {
return _CMInLab;
}

/**
*
* @param _CMInLab covariance matrix in lab frame
*/
public void set_CMInLab(double[][] _CMInLab) {
this._CMInLab = _CMInLab;
}

public void transCMToGlobal(){
MatrixOps mo = new MatrixOps(Libr.EJML);

int sector = this.getSector();

int charge = this.get_Q();
double p = this.get_P();

Vector3D momGlobal = this.get_pAtOrig();
Cross C = this.get(this.size() - 1);
Point3D momLocal = C.getCoordsInTiltedSector(momGlobal.x(), momGlobal.y(), momGlobal.z());
if(momLocal.z() == 0) set_CMInLab(new double[6][6]);
double tx = momLocal.x()/momLocal.z();
double ty = momLocal.y()/momLocal.z();
double theta = Math.atan(Math.sqrt(tx*tx + ty*ty));
double phi = Math.atan2(ty, tx);

//jm: Jacobi matrix for transformation from (x, y, tx, ty, Q) to (x, y, z, θ, φ, p) in the tilted sector frame
double dpdQ = -charge * p * p;
double dphidtx = -ty / (tx*tx + ty*ty);
double dphidty = tx / (tx*tx + ty*ty);
double dthetadtx = tx/(1 + tx*tx + ty*ty)/Math.sqrt(tx*tx + ty*ty);
double dthetadty = ty/(1 + tx*tx + ty*ty)/Math.sqrt(tx*tx + ty*ty);
double[][] jm = {{1, 0, 0, 0, 0}, {0, 1, 0, 0, 0}, {0, 0, 0, 0, 0},
{0, 0, dthetadtx, dthetadty, 0}, {0, 0, dphidtx, dphidty, 0}, {0, 0, 0, 0, dpdQ}
};
double[][] jmT = mo.MatrixTranspose(jm);

//jt: Jacobi matrix for transformation from (x, y, z, θ, φ, p) in the tilted sector frame to (x’, y’, z’, θ’, φ’, p’) in the sector frame (rotate 250 around y)
double numerator11 = Constants.SIN25 * Math.cos(theta) * Math.cos(phi) + Constants.COS25 * Math.sin(theta);
double numerator12 = -Constants.SIN25 * Math.sin(theta) * Math.sin(phi);
double denominator1 = Math.sqrt(1 - Math.pow(Constants.SIN25*Math.sin(theta)*Math.cos(phi) - Constants.COS25*Math.cos(theta), 2));
double numerator21 = Constants.SIN25 * Math.sin(phi);
double numerator22 = Constants.COS25 * Math.pow(Math.sin(theta),2) + Constants.SIN25 * Math.sin(theta) * Math.cos(theta) * Math.cos(phi);
double denoninator2 = Math.pow(Math.sin(theta) * Math.sin(phi), 2) + Math.pow(Constants.COS25 * Math.sin(theta) * Math.cos(phi) + Constants.SIN25 * Math.cos(theta),2);
double[][] jt = {{Constants.COS25, 0, Constants.SIN25, 0, 0, 0}, {0, 1, 0, 0, 0, 0}, {-Constants.SIN25, 0, Constants.COS25, 0, 0, 0},
{0, 0, 0, numerator11/denominator1, numerator12/denominator1, 0}, {0, 0, 0, numerator21/denoninator2, numerator22/denoninator2, 0}, {0, 0, 0, 0, 0, 1}
};
double[][] jtT = mo.MatrixTranspose(jt);

//js: Jacobi matrix for transformation from (x’, y’, z’, θ’, φ’, p’) to (x’’, y’’, z’’, θ’’, φ’’, p’’) in the lab frame (rotate around z’)
double[][] js = {{Constants.COSSECTOR60[sector - 1], -Constants.SINSECTOR60[sector - 1], 0, 0, 0, 0}, {Constants.SINSECTOR60[sector - 1], Constants.COSSECTOR60[sector - 1], 0, 0, 0, 0}, {0, 0, 1, 0, 0, 0},
{0, 0, 0, 1, 0, 0}, {0, 0, 0, 0, 1, 0}, {0, 0, 0, 0, 0, 1}
};
double[][] jsT = mo.MatrixTranspose(js);

Matrix CM = this.get_CovMat();
double[][] CMOrig = {
{CM.get(0, 0), CM.get(0, 1), CM.get(0, 2), CM.get(0, 3), CM.get(0, 4)},
{CM.get(1, 0), CM.get(1, 1), CM.get(1, 2), CM.get(1, 3), CM.get(1, 4)},
{CM.get(2, 0), CM.get(2, 1), CM.get(2, 2), CM.get(2, 3), CM.get(2, 4)},
{CM.get(3, 0), CM.get(3, 1), CM.get(3, 2), CM.get(3, 3), CM.get(3, 4)},
{CM.get(4, 0), CM.get(4, 1), CM.get(4, 2), CM.get(4, 3), CM.get(4, 4)}

};

double[][] transCM1 = mo.MatrixMultiplication(jm, CMOrig);
double[][] transCM2 = mo.MatrixMultiplication(transCM1, jmT);

double[][] transCM3 = mo.MatrixMultiplication(jt, transCM2);
double[][] transCM4 = mo.MatrixMultiplication(transCM3, jtT);

double[][] transCM5 = mo.MatrixMultiplication(js, transCM4);
double[][] transCM6 = mo.MatrixMultiplication(transCM5, jsT);

set_CMInLab(transCM6);
}

/**
*
* @param fitConvergenceStatus fit convergence status 0 if OK, 1 if the fit exits before converging
Expand Down
12 changes: 7 additions & 5 deletions reconstruction/dc/src/main/java/org/jlab/service/dc/DCTBEngine.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -263,8 +263,9 @@ public boolean processDataEvent(DataEvent event) {
}

// get CovMat at vertex
Point3D VTCS = crosses.get(0).getCoordsInTiltedSector(TrackArray1.get_Vtx0().x(), TrackArray1.get_Vtx0().y(), TrackArray1.get_Vtx0().z());
TrackArray1.set_CovMat(kFZRef.propagateToVtx(crosses.get(0).get_Sector(), VTCS.z()));
Point3D VTCS = TrackArray1.get(TrackArray1.size()-1).getCoordsInTiltedSector(TrackArray1.get_Vtx0().x(), TrackArray1.get_Vtx0().y(), TrackArray1.get_Vtx0().z());
TrackArray1.set_CovMat(kFZRef.propagateToVtx(TrackArray1.get(TrackArray1.size()-1).get_Sector(), VTCS.z()));
TrackArray1.transCMToGlobal();

double deltaPathToVtx = kFZRef.getDeltaPathToVtx(TrackArray1.get(TrackArray1.size()-1).get_Sector(), VTCS.z());
List<org.jlab.rec.dc.trajectory.StateVec> kfStateVecsAlongTrajectory = setKFStateVecsAlongTrajectory(kFZRef, deltaPathToVtx);
Expand Down Expand Up @@ -306,9 +307,10 @@ public boolean processDataEvent(DataEvent event) {
continue;
}

// get CovMat at vertex
Point3D VTCS = crosses.get(0).getCoordsInTiltedSector(TrackArray1.get_Vtx0().x(), TrackArray1.get_Vtx0().y(), TrackArray1.get_Vtx0().z());
TrackArray1.set_CovMat(kFZRef.propagateToVtx(crosses.get(0).get_Sector(), VTCS.z()));
// get CovMat at vertex
Point3D VTCS = TrackArray1.get(TrackArray1.size()-1).getCoordsInTiltedSector(TrackArray1.get_Vtx0().x(), TrackArray1.get_Vtx0().y(), TrackArray1.get_Vtx0().z());
TrackArray1.set_CovMat(kFZRef.propagateToVtx(TrackArray1.get(TrackArray1.size()-1).get_Sector(), VTCS.z()));
TrackArray1.transCMToGlobal();

double deltaPathToVtx = kFZRef.getDeltaPathToVtx(TrackArray1.get(TrackArray1.size()-1).get_Sector(), VTCS.z());
List<org.jlab.rec.dc.trajectory.StateVec> kfStateVecsAlongTrajectory = setKFStateVecsAlongTrajectory(kFZRef, deltaPathToVtx);
Expand Down