Skip to content

Commit a854a60

Browse files
efucheybaltzell
efuchey
authored and
baltzell
committed
Looping on all track candidates in Kalman filter instead of the first track only, to address:
#553 Added simulation flag self-configuration in AHDC Engine processing; Initialized AHDC Kalman Filter vertex constraint to 0 to avoid a nonsensical value if the vertex constraint is enabled but no vertex constraint is available. Ensured this constraint will not be fetched from MC banks if MC banks are not available.
1 parent 6fbfa00 commit a854a60

File tree

2 files changed

+88
-84
lines changed

2 files changed

+88
-84
lines changed

‎reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KalmanFilter.java‎

Lines changed: 84 additions & 84 deletions
Original file line numberDiff line numberDiff line change
@@ -39,7 +39,7 @@ public class KalmanFilter {
3939
private void propagation(ArrayList<Track> tracks, DataEvent event, boolean IsMC) {
4040

4141
try {
42-
double vz_constraint;
42+
double vz_constraint = 0;
4343
if(IsMC) {//If simulation read MC::Particle Bank ------------------------------------------------
4444
DataBank bankParticle = event.getBank("MC::Particle");
4545
double vxmc = bankParticle.getFloat("vx", 0)*10;//mm
@@ -77,24 +77,23 @@ private void propagation(ArrayList<Track> tracks, DataEvent event, boolean IsMC)
7777
// Initialization material map
7878
HashMap<String, Material> materialHashMap = materialGeneration();
7979

80-
// Initialization State Vector
81-
final double x0 = 0.0;
82-
final double y0 = 0.0;
83-
final double z0 = tracks.get(0).get_Z0();
84-
//final
85-
double px0 = tracks.get(0).get_px();
86-
//final
87-
double py0 = tracks.get(0).get_py();
88-
final double pz0 = tracks.get(0).get_pz();
89-
//final double p_init = java.lang.Math.sqrt(px0*px0+py0*py0+pz0*pz0);
90-
double[] y = new double[]{x0, y0, z0, px0, py0, pz0};
91-
// EPAF: *the line below is for TEST ONLY!!!*
92-
//double[] y = new double[]{vxmc, vymc, vzmc, pxmc, pymc, pzmc};
93-
// Initialization hit
94-
ArrayList<org.jlab.rec.ahdc.Hit.Hit> AHDC_hits = tracks.get(0).getHits();
95-
ArrayList<Hit> KF_hits = new ArrayList<>();
96-
//System.out.println(" px " + y[3] + " py " + y[4] +" pz " + y[5] +" vz " + y[2] + " number of hits: " + AHDC_hits.size() + " MC hits? " + sim_hits.size());
97-
for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits) {
80+
for (Track track : tracks) {
81+
// Initialization State Vector
82+
final double x0 = 0.0;
83+
final double y0 = 0.0;
84+
final double z0 = track.get_Z0();
85+
//final
86+
double px0 = track.get_px();
87+
//final
88+
double py0 = track.get_py();
89+
final double pz0 = track.get_pz();
90+
//final double p_init = java.lang.Math.sqrt(px0*px0+py0*py0+pz0*pz0);
91+
double[] y = new double[]{x0, y0, z0, px0, py0, pz0};
92+
// Initialization hit
93+
ArrayList<org.jlab.rec.ahdc.Hit.Hit> AHDC_hits = track.getHits();
94+
ArrayList<Hit> KF_hits = new ArrayList<>();
95+
//System.out.println(" px " + y[3] + " py " + y[4] +" pz " + y[5] +" vz " + y[2] + " number of hits: " + AHDC_hits.size() + " MC hits? " + sim_hits.size());
96+
for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits) {
9897
Hit hit = new Hit(AHDC_hit.getSuperLayerId(), AHDC_hit.getLayerId(), AHDC_hit.getWireId(), AHDC_hit.getNbOfWires(), AHDC_hit.getRadius(), AHDC_hit.getDoca());
9998
hit.setADC(AHDC_hit.getADC());
10099
hit.setHitIdx(AHDC_hit.getId());
@@ -104,97 +103,98 @@ private void propagation(ArrayList<Track> tracks, DataEvent event, boolean IsMC)
104103
boolean phi_rollover = false;
105104
boolean aleardyHaveR = false;
106105
for (Hit o: KF_hits){
107-
if (o.r() == hit.r()){
108-
aleardyHaveR = true;
109-
// //sign+ means (phi track - phi wire) > 0
110-
// if(o.phi()>hit.phi()){
111-
// if(Math.abs(o.phi()-hit.phi())< 2*Math.toRadians(360./o.getNumWires()) ){
112-
// o.setSign(-1);
113-
// hit.setSign(+1);
114-
// }else{
115-
// phi_rollover = true;
116-
// hit.setSign(-1);
117-
// o.setSign(+1);
118-
// }
119-
// }else{
120-
// if(Math.abs(o.phi()-hit.phi())< 2*Math.toRadians(360./o.getNumWires()) ){
121-
// hit.setSign(-1);
122-
// o.setSign(+1);
123-
// }else{
124-
// phi_rollover = true;
125-
// o.setSign(-1);
126-
// hit.setSign(+1);
127-
// }
128-
// }
129-
// //System.out.println( " r = " + o.r() + " o.phi = " + o.phi() + " o.doca = " + o.getDoca()*o.getSign() + " hit.phi " + hit.phi() +" hit.doca = " + hit.getDoca()*hit.getSign() + " angle between wires: " + Math.toRadians(360./hit.getNumWires()) + " >= ? angle covered by docas: " + Math.atan( (o.getDoca()+hit.getDoca())/o.r() ) );
130-
}
106+
if (o.r() == hit.r()){
107+
aleardyHaveR = true;
108+
// //sign+ means (phi track - phi wire) > 0
109+
// if(o.phi()>hit.phi()){
110+
// if(Math.abs(o.phi()-hit.phi())< 2*Math.toRadians(360./o.getNumWires()) ){
111+
// o.setSign(-1);
112+
// hit.setSign(+1);
113+
// }else{
114+
// phi_rollover = true;
115+
// hit.setSign(-1);
116+
// o.setSign(+1);
117+
// }
118+
// }else{
119+
// if(Math.abs(o.phi()-hit.phi())< 2*Math.toRadians(360./o.getNumWires()) ){
120+
// hit.setSign(-1);
121+
// o.setSign(+1);
122+
// }else{
123+
// phi_rollover = true;
124+
// o.setSign(-1);
125+
// hit.setSign(+1);
126+
// }
127+
// }
128+
// //System.out.println( " r = " + o.r() + " o.phi = " + o.phi() + " o.doca = " + o.getDoca()*o.getSign() + " hit.phi " + hit.phi() +" hit.doca = " + hit.getDoca()*hit.getSign() + " angle between wires: " + Math.toRadians(360./hit.getNumWires()) + " >= ? angle covered by docas: " + Math.atan( (o.getDoca()+hit.getDoca())/o.r() ) );
129+
}
131130
}
132131
if(!aleardyHaveR)KF_hits.add(hit);
133132
// if (phi_rollover){
134133
// KF_hits.add(KF_hits.size()-1, hit);
135134
// }else{
136135
// KF_hits.add(hit);
137136
// }
138-
}
137+
}
139138

140-
double zbeam = 0;
141-
if(IsVtxDefined)zbeam = vz_constraint;//test
142-
final ArrayList<Indicator> forwardIndicators = forwardIndicators(KF_hits, materialHashMap);
143-
final ArrayList<Indicator> backwardIndicators = backwardIndicators(KF_hits, materialHashMap, zbeam);
139+
double zbeam = 0;
140+
if(IsVtxDefined)zbeam = vz_constraint;//test
141+
final ArrayList<Indicator> forwardIndicators = forwardIndicators(KF_hits, materialHashMap);
142+
final ArrayList<Indicator> backwardIndicators = backwardIndicators(KF_hits, materialHashMap, zbeam);
144143

145-
// Start propagation
146-
Stepper stepper = new Stepper(y);
147-
RungeKutta4 RK4 = new RungeKutta4(proton, numberOfVariables, B);
148-
Propagator propagator = new Propagator(RK4);
149-
150-
// ----------------------------------------------------------------------------------------
151-
152-
// Initialization of the Kalman Fitter
153-
RealVector initialStateEstimate = new ArrayRealVector(stepper.y);
154-
//first 3 lines in cm^2; last 3 lines in MeV^2
155-
RealMatrix initialErrorCovariance = MatrixUtils.createRealMatrix(new double[][]{{1.00, 0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 1.00, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 25.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 1.00, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 1.00, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0, 25.0}});
156-
KFitter kFitter = new KFitter(initialStateEstimate, initialErrorCovariance, stepper, propagator);
157-
kFitter.setVertexDefined(IsVtxDefined);
144+
// Start propagation
145+
Stepper stepper = new Stepper(y);
146+
RungeKutta4 RK4 = new RungeKutta4(proton, numberOfVariables, B);
147+
Propagator propagator = new Propagator(RK4);
148+
149+
// ----------------------------------------------------------------------------------------
150+
151+
// Initialization of the Kalman Fitter
152+
RealVector initialStateEstimate = new ArrayRealVector(stepper.y);
153+
//first 3 lines in cm^2; last 3 lines in MeV^2
154+
RealMatrix initialErrorCovariance = MatrixUtils.createRealMatrix(new double[][]{{1.00, 0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 1.00, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 25.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 1.00, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 1.00, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0, 25.0}});
155+
KFitter kFitter = new KFitter(initialStateEstimate, initialErrorCovariance, stepper, propagator);
156+
kFitter.setVertexDefined(IsVtxDefined);
158157

159-
for (int k = 0; k < Niter; k++) {
158+
for (int k = 0; k < Niter; k++) {
160159
//System.out.println("--------- ForWard propagation !! ---------");
161160
//Reset error covariance:
162161
//kFitter.ResetErrorCovariance(initialErrorCovariance);
163162
for (Indicator indicator : forwardIndicators) {
164-
kFitter.predict(indicator);
165-
if (indicator.haveAHit()) {
166-
if( k==0 && indicator.hit.getHitIdx()>0){
167-
for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits){
168-
if(AHDC_hit.getId()==indicator.hit.getHitIdx())AHDC_hit.setResidualPrefit(kFitter.residual(indicator));
169-
}
170-
}
171-
kFitter.correct(indicator);
163+
kFitter.predict(indicator);
164+
if (indicator.haveAHit()) {
165+
if( k==0 && indicator.hit.getHitIdx()>0){
166+
for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits){
167+
if(AHDC_hit.getId()==indicator.hit.getHitIdx())AHDC_hit.setResidualPrefit(kFitter.residual(indicator));
168+
}
172169
}
170+
kFitter.correct(indicator);
171+
}
173172
}
174173

175174
//System.out.println("--------- BackWard propagation !! ---------");
176175
for (Indicator indicator : backwardIndicators) {
177-
kFitter.predict(indicator);
178-
if (indicator.haveAHit()) {
179-
kFitter.correct(indicator);
180-
}
176+
kFitter.predict(indicator);
177+
if (indicator.haveAHit()) {
178+
kFitter.correct(indicator);
179+
}
181180
}
182-
}
181+
}
183182

184-
RealVector x_out = kFitter.getStateEstimationVector();
185-
tracks.get(0).setPositionAndMomentumForKF(x_out);
183+
RealVector x_out = kFitter.getStateEstimationVector();
184+
track.setPositionAndMomentumForKF(x_out);
186185

187-
//Residual calcuation post fit:
188-
for (Indicator indicator : forwardIndicators) {
186+
//Residual calcuation post fit:
187+
for (Indicator indicator : forwardIndicators) {
189188
kFitter.predict(indicator);
190189
if (indicator.haveAHit()) {
191-
if( indicator.hit.getHitIdx()>0){
192-
for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits){
193-
if(AHDC_hit.getId()==indicator.hit.getHitIdx())AHDC_hit.setResidual(kFitter.residual(indicator));
194-
}
190+
if( indicator.hit.getHitIdx()>0){
191+
for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits){
192+
if(AHDC_hit.getId()==indicator.hit.getHitIdx())AHDC_hit.setResidual(kFitter.residual(indicator));
195193
}
194+
}
196195
}
197-
}
196+
}
197+
}//end of loop on track candidates
198198
} catch (Exception e) {
199199
// e.printStackTrace();
200200
}

‎reconstruction/alert/src/main/java/org/jlab/service/ahdc/AHDCEngine.java‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -113,6 +113,10 @@ public boolean processDataEvent(DataEvent event) {
113113
double magfield = 50.0; // what is this?
114114
double magfieldfactor = 1; // why is this here?
115115

116+
if(event.hasBank("MC::Particle")){
117+
simulation = true;
118+
}
119+
116120
if (event.hasBank("RUN::config")) {
117121
DataBank bank = event.getBank("RUN::config");
118122
runNo = bank.getInt("run", 0);

0 commit comments

Comments
 (0)