Skip to content

Commit bf360d3

Browse files
tongtongcaotongtongcao
committed
transformation for covaraince matrix of DC TB tracks from local to global
1 parent d59c682 commit bf360d3

7 files changed

Lines changed: 254 additions & 7 deletions

File tree

‎common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/utilities/MatrixOps.java‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -128,7 +128,7 @@ public double[][] MatrixTranspose(Object obj1) {
128128

129129
for(int i = 0; i < r; i++) {
130130
for (int j = 0; j < c; j++) {
131-
result[i][j] = arr1[j][i];
131+
result[j][i] = arr1[i][j];
132132
}
133133
}
134134
arr1 = null;

‎etc/bankdefs/hipo4/dc.json‎

Lines changed: 45 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -559,6 +559,51 @@
559559
{"name":"C55", "type":"F", "info":"C55 covariance matrix element at last superlayer used in the fit"}
560560
]
561561
},
562+
{
563+
"name": "TimeBasedTrkg::TBCovMatLab",
564+
"group": 20600,
565+
"item" : 38,
566+
"info": "reconstructed track covariance matrix in lab frame",
567+
"entries": [
568+
{"name":"id", "type":"S", "info":"id of the track"},
569+
{"name":"C11", "type":"F", "info":"var(x) at vertex in lab frame"},
570+
{"name":"C12", "type":"F", "info":"cov(x, y) at vertex in lab frame"},
571+
{"name":"C13", "type":"F", "info":"cov(x, z) at vertex in lab frame"},
572+
{"name":"C14", "type":"F", "info":"cov(x, p_theta) at vertex in lab frame"},
573+
{"name":"C15", "type":"F", "info":"cov(x, p_phi) at vertex in lab frame"},
574+
{"name":"C16", "type":"F", "info":"cov(x, p) at vertex in lab frame"},
575+
{"name":"C21", "type":"F", "info":"cov(y, x) at vertex in lab frame"},
576+
{"name":"C22", "type":"F", "info":"var(y) at vertex in lab frame"},
577+
{"name":"C23", "type":"F", "info":"cov(y, z) at vertex in lab frame"},
578+
{"name":"C24", "type":"F", "info":"cov(y, p_theta) at vertex in lab frame"},
579+
{"name":"C25", "type":"F", "info":"cov(y, p_phi) at vertex in lab frame"},
580+
{"name":"C26", "type":"F", "info":"cov(y, p) at vertex in lab frame"},
581+
{"name":"C31", "type":"F", "info":"cov(z, x) at vertex in lab frame"},
582+
{"name":"C32", "type":"F", "info":"cov(z, y) at vertex in lab frame"},
583+
{"name":"C33", "type":"F", "info":"var(z) at vertex in lab frame"},
584+
{"name":"C34", "type":"F", "info":"cov(z, p_theta) at vertex in lab frame"},
585+
{"name":"C35", "type":"F", "info":"cov(z, p_phi) at vertex in lab frame"},
586+
{"name":"C36", "type":"F", "info":"cov(z, p) at vertex in lab frame"},
587+
{"name":"C41", "type":"F", "info":"cov(p_theta, x) at vertex in lab frame"},
588+
{"name":"C42", "type":"F", "info":"cov(p_theta, y) at vertex in lab frame"},
589+
{"name":"C43", "type":"F", "info":"cov(p_theta, z) at vertex in lab frame"},
590+
{"name":"C44", "type":"F", "info":"var(p_theta) at vertex in lab frame"},
591+
{"name":"C45", "type":"F", "info":"cov(p_theta, p_phi) at vertex in lab frame"},
592+
{"name":"C46", "type":"F", "info":"cov(p_theta, p) at vertex in lab frame"},
593+
{"name":"C51", "type":"F", "info":"cov(p_phi, x) at vertex in lab frame"},
594+
{"name":"C52", "type":"F", "info":"cov(p_phi, y) at vertex in lab frame"},
595+
{"name":"C53", "type":"F", "info":"cov(p_phi, z) at vertex in lab frame"},
596+
{"name":"C54", "type":"F", "info":"cov(p_phi, p_theta) at vertex in lab frame"},
597+
{"name":"C55", "type":"F", "info":"var(p_phi) at vertex in lab frame"},
598+
{"name":"C56", "type":"F", "info":"cov(p_phi, p) at vertex in lab frame"},
599+
{"name":"C61", "type":"F", "info":"cov(p, x) at vertex in lab frame"},
600+
{"name":"C62", "type":"F", "info":"cov(p, y) at vertex in lab frame"},
601+
{"name":"C63", "type":"F", "info":"cov(p, z) at vertex in lab frame"},
602+
{"name":"C64", "type":"F", "info":"cov(p, p_theta) at vertex in lab frame"},
603+
{"name":"C65", "type":"F", "info":"cov(p, p_phi) at vertex in lab frame"},
604+
{"name":"C66", "type":"F", "info":"var(p) at vertex in lab frame"}
605+
]
606+
},
562607
{
563608
"name": "TimeBasedTrkg::Trajectory",
564609
"group": 20600,

‎etc/bankdefs/hipo4/dcnn.json‎

Lines changed: 45 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -420,6 +420,51 @@
420420
{"name":"C55", "type":"F", "info":"C55 covariance matrix element at last superlayer used in the fit"}
421421
]
422422
},
423+
{
424+
"name": "TimeBasedTrkg::AICovMatLab",
425+
"group": 20600,
426+
"item" : 78,
427+
"info": "reconstructed track covariance matrix in lab frame",
428+
"entries": [
429+
{"name":"id", "type":"S", "info":"id of the track"},
430+
{"name":"C11", "type":"F", "info":"var(x) at vertex in lab frame"},
431+
{"name":"C12", "type":"F", "info":"cov(x, y) at vertex in lab frame"},
432+
{"name":"C13", "type":"F", "info":"cov(x, z) at vertex in lab frame"},
433+
{"name":"C14", "type":"F", "info":"cov(x, p_theta) at vertex in lab frame"},
434+
{"name":"C15", "type":"F", "info":"cov(x, p_phi) at vertex in lab frame"},
435+
{"name":"C16", "type":"F", "info":"cov(x, p) at vertex in lab frame"},
436+
{"name":"C21", "type":"F", "info":"cov(y, x) at vertex in lab frame"},
437+
{"name":"C22", "type":"F", "info":"var(y) at vertex in lab frame"},
438+
{"name":"C23", "type":"F", "info":"cov(y, z) at vertex in lab frame"},
439+
{"name":"C24", "type":"F", "info":"cov(y, p_theta) at vertex in lab frame"},
440+
{"name":"C25", "type":"F", "info":"cov(y, p_phi) at vertex in lab frame"},
441+
{"name":"C26", "type":"F", "info":"cov(y, p) at vertex in lab frame"},
442+
{"name":"C31", "type":"F", "info":"cov(z, x) at vertex in lab frame"},
443+
{"name":"C32", "type":"F", "info":"cov(z, y) at vertex in lab frame"},
444+
{"name":"C33", "type":"F", "info":"var(z) at vertex in lab frame"},
445+
{"name":"C34", "type":"F", "info":"cov(z, p_theta) at vertex in lab frame"},
446+
{"name":"C35", "type":"F", "info":"cov(z, p_phi) at vertex in lab frame"},
447+
{"name":"C36", "type":"F", "info":"cov(z, p) at vertex in lab frame"},
448+
{"name":"C41", "type":"F", "info":"cov(p_theta, x) at vertex in lab frame"},
449+
{"name":"C42", "type":"F", "info":"cov(p_theta, y) at vertex in lab frame"},
450+
{"name":"C43", "type":"F", "info":"cov(p_theta, z) at vertex in lab frame"},
451+
{"name":"C44", "type":"F", "info":"var(p_theta) at vertex in lab frame"},
452+
{"name":"C45", "type":"F", "info":"cov(p_theta, p_phi) at vertex in lab frame"},
453+
{"name":"C46", "type":"F", "info":"cov(p_theta, p) at vertex in lab frame"},
454+
{"name":"C51", "type":"F", "info":"cov(p_phi, x) at vertex in lab frame"},
455+
{"name":"C52", "type":"F", "info":"cov(p_phi, y) at vertex in lab frame"},
456+
{"name":"C53", "type":"F", "info":"cov(p_phi, z) at vertex in lab frame"},
457+
{"name":"C54", "type":"F", "info":"cov(p_phi, p_theta) at vertex in lab frame"},
458+
{"name":"C55", "type":"F", "info":"var(p_phi) at vertex in lab frame"},
459+
{"name":"C56", "type":"F", "info":"cov(p_phi, p) at vertex in lab frame"},
460+
{"name":"C61", "type":"F", "info":"cov(p, x) at vertex in lab frame"},
461+
{"name":"C62", "type":"F", "info":"cov(p, y) at vertex in lab frame"},
462+
{"name":"C63", "type":"F", "info":"cov(p, z) at vertex in lab frame"},
463+
{"name":"C64", "type":"F", "info":"cov(p, p_theta) at vertex in lab frame"},
464+
{"name":"C65", "type":"F", "info":"cov(p, p_phi) at vertex in lab frame"},
465+
{"name":"C66", "type":"F", "info":"var(p) at vertex in lab frame"}
466+
]
467+
},
423468
{
424469
"name": "TimeBasedTrkg::AITrajectory",
425470
"group": 20600,

‎reconstruction/dc/src/main/java/org/jlab/rec/dc/banks/Banks.java‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -138,6 +138,10 @@ public String getCovmatBank() {
138138
return this.getOutputBank("CovMat");
139139
}
140140

141+
public String getCovmatLabBank() {
142+
return this.getOutputBank("CovMatLab");
143+
}
144+
141145
public String getRecEventBank() {
142146
return this.getRecBank("Event");
143147
}

‎reconstruction/dc/src/main/java/org/jlab/rec/dc/banks/RecoBankWriter.java‎

Lines changed: 65 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -490,6 +490,69 @@ private DataBank fillTrackCovMatBank(DataEvent event, List<Track> candlist) {
490490
//bank.show();
491491
return bank;
492492
}
493+
494+
/**
495+
*
496+
* @param event hipo event
497+
* @param candlist tracks
498+
* @return covariance matrix for momentum and vertex in lab frame
499+
*/
500+
private DataBank fillTrackCovMatLabBank(DataEvent event, List<Track> candlist) {
501+
502+
DataBank bank = event.createBank(bankNames.getCovmatLabBank(), candlist.size());
503+
504+
for (int i = 0; i < candlist.size(); i++) {
505+
bank.setShort("id", i, (short) candlist.get(i).get_Id());
506+
if(candlist.get(i).get_CMInLab()!=null) {
507+
double[][] CM = candlist.get(i).get_CMInLab();
508+
bank.setFloat("C11", i, (float) CM[0][0]);
509+
bank.setFloat("C12", i, (float) CM[0][1]);
510+
bank.setFloat("C13", i, (float) CM[0][2]);
511+
bank.setFloat("C14", i, (float) CM[0][3]);
512+
bank.setFloat("C15", i, (float) CM[0][4]);
513+
bank.setFloat("C16", i, (float) CM[0][5]);
514+
515+
bank.setFloat("C21", i, (float) CM[1][0]);
516+
bank.setFloat("C22", i, (float) CM[1][1]);
517+
bank.setFloat("C23", i, (float) CM[1][2]);
518+
bank.setFloat("C24", i, (float) CM[1][3]);
519+
bank.setFloat("C25", i, (float) CM[1][4]);
520+
bank.setFloat("C26", i, (float) CM[1][5]);
521+
522+
bank.setFloat("C31", i, (float) CM[2][0]);
523+
bank.setFloat("C32", i, (float) CM[2][1]);
524+
bank.setFloat("C33", i, (float) CM[2][2]);
525+
bank.setFloat("C34", i, (float) CM[2][3]);
526+
bank.setFloat("C35", i, (float) CM[2][4]);
527+
bank.setFloat("C36", i, (float) CM[2][5]);
528+
529+
bank.setFloat("C41", i, (float) CM[3][0]);
530+
bank.setFloat("C42", i, (float) CM[3][1]);
531+
bank.setFloat("C43", i, (float) CM[3][2]);
532+
bank.setFloat("C44", i, (float) CM[3][3]);
533+
bank.setFloat("C45", i, (float) CM[3][4]);
534+
bank.setFloat("C46", i, (float) CM[3][5]);
535+
536+
bank.setFloat("C51", i, (float) CM[4][0]);
537+
bank.setFloat("C52", i, (float) CM[4][1]);
538+
bank.setFloat("C53", i, (float) CM[4][2]);
539+
bank.setFloat("C54", i, (float) CM[4][3]);
540+
bank.setFloat("C55", i, (float) CM[4][4]);
541+
bank.setFloat("C56", i, (float) CM[4][5]);
542+
543+
bank.setFloat("C61", i, (float) CM[5][0]);
544+
bank.setFloat("C62", i, (float) CM[5][1]);
545+
bank.setFloat("C63", i, (float) CM[5][2]);
546+
bank.setFloat("C64", i, (float) CM[5][3]);
547+
bank.setFloat("C65", i, (float) CM[5][4]);
548+
bank.setFloat("C66", i, (float) CM[5][5]);
549+
550+
}
551+
}
552+
//bank.show();
553+
return bank;
554+
}
555+
493556
/**
494557
*
495558
* @param event the EvioEvent
@@ -939,7 +1002,8 @@ public void fillAllTBBanks(DataEvent event, List<FittedHit> fhits, List<FittedCl
9391002
this.fillTBCrossesBank(event, crosses),
9401003
this.fillTBTracksBank(event, trkcands),
9411004
this.fillTrajectoryBank(event, trkcands),
942-
this.fillTrackCovMatBank(event, trkcands)
1005+
this.fillTrackCovMatBank(event, trkcands),
1006+
this.fillTrackCovMatLabBank(event, trkcands)
9431007
);
9441008
}
9451009
if (crosses != null && trkcands == null) {

‎reconstruction/dc/src/main/java/org/jlab/rec/dc/track/Track.java‎

Lines changed: 87 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,8 @@
1111
import org.jlab.rec.dc.segment.Segment;
1212
import org.jlab.rec.dc.trajectory.StateVec;
1313
import org.jlab.rec.dc.trajectory.Trajectory;
14+
import org.jlab.clas.tracking.utilities.MatrixOps;
15+
import org.jlab.clas.tracking.utilities.MatrixOps.Libr;
1416

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

4649
private Point3D _Region3CrossPoint;
4750
private Point3D _Region3CrossDir;
@@ -396,6 +399,90 @@ public void set_CovMat(Matrix _CovMat) {
396399
this._CovMat = _CovMat;
397400
}
398401

402+
/**
403+
*
404+
* @return covariance matrix in lab frame
405+
*/
406+
public double[][] get_CMInLab() {
407+
return _CMInLab;
408+
}
409+
410+
/**
411+
*
412+
* @param _CMInLab covariance matrix in lab frame
413+
*/
414+
public void set_CMInLab(double[][] _CMInLab) {
415+
this._CMInLab = _CMInLab;
416+
}
417+
418+
public void transCMToGlobal(){
419+
MatrixOps mo = new MatrixOps(Libr.EJML);
420+
421+
int sector = this.getSector();
422+
423+
int charge = this.get_Q();
424+
double p = this.get_P();
425+
426+
Vector3D momGlobal = this.get_pAtOrig();
427+
Cross C = this.get(this.size() - 1);
428+
Point3D momLocal = C.getCoordsInTiltedSector(momGlobal.x(), momGlobal.y(), momGlobal.z());
429+
if(momLocal.z() == 0) set_CMInLab(new double[6][6]);
430+
double tx = momLocal.x()/momLocal.z();
431+
double ty = momLocal.y()/momLocal.z();
432+
double theta = Math.atan(Math.sqrt(tx*tx + ty*ty));
433+
double phi = Math.atan2(ty, tx);
434+
435+
//jm: Jacobi matrix for transformation from (x, y, tx, ty, Q) to (x, y, z, θ, φ, p) in the tilted sector frame
436+
double dpdQ = -charge * p * p;
437+
double dphidtx = -ty / (tx*tx + ty*ty);
438+
double dphidty = tx / (tx*tx + ty*ty);
439+
double dthetadtx = tx/(1 + tx*tx + ty*ty)/Math.sqrt(tx*tx + ty*ty);
440+
double dthetadty = ty/(1 + tx*tx + ty*ty)/Math.sqrt(tx*tx + ty*ty);
441+
double[][] jm = {{1, 0, 0, 0, 0}, {0, 1, 0, 0, 0}, {0, 0, 0, 0, 0},
442+
{0, 0, dthetadtx, dthetadty, 0}, {0, 0, dphidtx, dphidty, 0}, {0, 0, 0, 0, dpdQ}
443+
};
444+
double[][] jmT = mo.MatrixTranspose(jm);
445+
446+
//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)
447+
double numerator11 = Constants.SIN25 * Math.cos(theta) * Math.cos(phi) + Constants.COS25 * Math.sin(theta);
448+
double numerator12 = -Constants.SIN25 * Math.sin(theta) * Math.sin(phi);
449+
double denominator1 = Math.sqrt(1 - Math.pow(Constants.SIN25*Math.sin(theta)*Math.cos(phi) - Constants.COS25*Math.cos(theta), 2));
450+
double numerator21 = Constants.SIN25 * Math.sin(phi);
451+
double numerator22 = Constants.COS25 * Math.pow(Math.sin(theta),2) + Constants.SIN25 * Math.sin(theta) * Math.cos(theta) * Math.cos(phi);
452+
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);
453+
double[][] jt = {{Constants.COS25, 0, Constants.SIN25, 0, 0, 0}, {0, 1, 0, 0, 0, 0}, {-Constants.SIN25, 0, Constants.COS25, 0, 0, 0},
454+
{0, 0, 0, numerator11/denominator1, numerator12/denominator1, 0}, {0, 0, 0, numerator21/denoninator2, numerator22/denoninator2, 0}, {0, 0, 0, 0, 0, 1}
455+
};
456+
double[][] jtT = mo.MatrixTranspose(jt);
457+
458+
//js: Jacobi matrix for transformation from (x’, y’, z’, θ’, φ’, p’) to (x’’, y’’, z’’, θ’’, φ’’, p’’) in the lab frame (rotate around z’)
459+
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},
460+
{0, 0, 0, 1, 0, 0}, {0, 0, 0, 0, 1, 0}, {0, 0, 0, 0, 0, 1}
461+
};
462+
double[][] jsT = mo.MatrixTranspose(js);
463+
464+
Matrix CM = this.get_CovMat();
465+
double[][] CMOrig = {
466+
{CM.get(0, 0), CM.get(0, 1), CM.get(0, 2), CM.get(0, 3), CM.get(0, 4)},
467+
{CM.get(1, 0), CM.get(1, 1), CM.get(1, 2), CM.get(1, 3), CM.get(1, 4)},
468+
{CM.get(2, 0), CM.get(2, 1), CM.get(2, 2), CM.get(2, 3), CM.get(2, 4)},
469+
{CM.get(3, 0), CM.get(3, 1), CM.get(3, 2), CM.get(3, 3), CM.get(3, 4)},
470+
{CM.get(4, 0), CM.get(4, 1), CM.get(4, 2), CM.get(4, 3), CM.get(4, 4)}
471+
472+
};
473+
474+
double[][] transCM1 = mo.MatrixMultiplication(jm, CMOrig);
475+
double[][] transCM2 = mo.MatrixMultiplication(transCM1, jmT);
476+
477+
double[][] transCM3 = mo.MatrixMultiplication(jt, transCM2);
478+
double[][] transCM4 = mo.MatrixMultiplication(transCM3, jtT);
479+
480+
double[][] transCM5 = mo.MatrixMultiplication(js, transCM4);
481+
double[][] transCM6 = mo.MatrixMultiplication(transCM5, jsT);
482+
483+
set_CMInLab(transCM6);
484+
}
485+
399486
/**
400487
*
401488
* @param fitConvergenceStatus fit convergence status 0 if OK, 1 if the fit exits before converging

‎reconstruction/dc/src/main/java/org/jlab/service/dc/DCTBEngine.java‎

Lines changed: 7 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -263,8 +263,9 @@ public boolean processDataEvent(DataEvent event) {
263263
}
264264

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

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

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

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

0 commit comments

Comments
 (0)