diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Hit/Hit.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Hit/Hit.java index 4a04862c2d..008e892104 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Hit/Hit.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Hit/Hit.java @@ -179,7 +179,13 @@ public RealVector get_Vector() { } public RealMatrix get_MeasurementNoise() { - return new Array2DRowRealMatrix(new double[][]{{0.09}}); + double mean_error = 0.471; // mm (no difference between adc and time) + double error_on_adc = (1.15146*raw_adc + 437.63)/(3.21187*raw_adc + 878.855); // mm + double error_on_time = (0.4423*time + 13.7215)/(0.846038*time + 31.9867); // mm + double error = error_on_adc*error_on_time/mean_error; // mm + + return new Array2DRowRealMatrix(new double[][]{{Math.pow(error, 2)}}); // mm^2 + //return new Array2DRowRealMatrix(new double[][]{{0.09}}); } // a signature for KalmanFilter.Hit_beam @@ -206,6 +212,32 @@ public static void main(String[] args) { System.out.println("h1 compare to h2 : " + h1.compareTo(h2)); System.out.println("h2 compare to h1 : " + h2.compareTo(h1)); System.out.println("h1 compare to h3 : " + h1.compareTo(h3)); + + System.out.println("/////////////////////////"); + System.out.println("Test AHDC geometry"); + System.out.println(""); + System.out.println("s : sector"); + System.out.println("sl : super layer"); + System.out.println("l : layer"); + System.out.println("c : component"); + System.out.println("/////////////////////////"); + System.out.println("------------------------------------------------------------------------------"); + System.out.println(" | origin | end"); + System.out.println("------------------------------------------------------------------------------"); + System.out.println("s sl l c | x y z | x y z"); + System.out.println("------------------------------------------------------------------------------"); + for (int s = 1; s <= factory.getNumSectors(); s++) { + for (int sl = 1; sl <= factory.getSector(s).getNumSuperlayers(); sl++) { + for (int l = 1; l <= factory.getSector(s).getSuperlayer(sl).getNumLayers(); l++) { + for (int c = 1; c <= factory.getSector(s).getSuperlayer(sl).getLayer(l).getNumComponents(); c++) { + Line3D line = factory.getSector(s).getSuperlayer(sl).getLayer(l).getComponent(c).getLine(); + Point3D end = line.end(); + Point3D origin = line.origin(); + System.out.printf("%2d %2d %2d %2d | %7.3f %7.3f %7.3f | %7.3f %7.3f %7.3f\n", s, sl, l, c, origin.x(), origin.y(), origin.z(), end.x(), end.y(), end.z()); + } + } + } + } } } diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFitter.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFitter.java index 3f19dc0305..7263d6638f 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFitter.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFitter.java @@ -29,7 +29,7 @@ public class KFitter { // masses/energies in MeV private final double electron_mass_c2 = PhysicsConstants.massElectron() * 1000; private final double proton_mass_c2 = PhysicsConstants.massProton() * 1000; - private boolean isvertexdefined = false; + private double[] vertex_resolutions = {0.09,1e10}; // default values // dr^2 and dz^2 in mm^2 public KFitter(final RealVector initialStateEstimate, final RealMatrix initialErrorCovariance, final Stepper stepper, final Propagator propagator, final HashMap materialHashMap) { this.stateEstimation = initialStateEstimate; @@ -93,15 +93,14 @@ public void correct(Hit hit) { RealVector h; // check if the hit is the beamline if (hit.getRadius() < 1) { - double z_beam_res_sq = 1.e10;//in mm - if(isvertexdefined)z_beam_res_sq = 4.0;//assuming 2. mm resolution + // the diagonal elements are the squared errors in r, phi, z measurementNoise = - new Array2DRowRealMatrix( - new double[][]{ - {0.09, 0.0000, 0.0000}, - {0.00, 1e10, 0.0000}, - {0.00, 0.0000, z_beam_res_sq} - });//3x3 + new Array2DRowRealMatrix( + new double[][]{ + {vertex_resolutions[0], 0.0000, 0.0000}, + {0.00, 1e10, 0.0000}, + {0.00, 0.0000, vertex_resolutions[1]} + });//3x3 measurementMatrix = H_beam(stateEstimation);//6x3 h = h_beam(stateEstimation);//3x1 z = hit.get_Vector_beam();//0! @@ -274,7 +273,10 @@ public RealVector getStateEstimationVector() { public RealMatrix getErrorCovarianceMatrix() { return errorCovariance.copy(); } - - public void setVertexDefined(boolean isvtxdef) {isvertexdefined = isvtxdef;} + + public void setVertexResolution(double[] res) { + vertex_resolutions[0] = res[0]; + vertex_resolutions[1] = res[1]; + } } diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KalmanFilter.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KalmanFilter.java index 805a736558..4d2840ee73 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KalmanFilter.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KalmanFilter.java @@ -11,6 +11,7 @@ import org.jlab.clas.pdg.PDGDatabase; import org.jlab.clas.pdg.PDGParticle; import org.jlab.clas.tracking.kalmanfilter.Material; +import org.jlab.io.base.DataBank; import org.jlab.io.base.DataEvent; import org.jlab.rec.ahdc.Hit.Hit; import org.jlab.rec.ahdc.Track.Track; @@ -32,6 +33,9 @@ public class KalmanFilter { private final int Niter = 40; // number of iterations for the Kalman Filter private boolean IsVtxDefined = false; // implemented but not used yet + private double[] vertex_resolutions = {0.09, 1e10}; // {error in r squared in mm^2, error in z squared in mm^2} + // mm, CLAS and AHDC don't necessary have the same alignement (ZERO), this parameter may be subject to calibration + private double clas_alignement = -54; private void propagation(ArrayList tracks, DataEvent event, final double magfield, boolean IsMC) { @@ -44,6 +48,30 @@ private void propagation(ArrayList tracks, DataEvent event, final double final double tesla = 0.001; final double[] B = {0.0, 0.0, magfield / 10 * tesla}; HashMap materialHashMap = MaterialMap.generateMaterials(); + // Recover the vertex of the electron + if (event.hasBank("REC::Particle")) { + DataBank recBank = event.getBank("REC::Particle"); + int row = 0; + while ((!IsVtxDefined) && row < recBank.rows()) { + if (recBank.getInt("pid", row) == 11) { + IsVtxDefined = true; + vz_constraint = 10*recBank.getFloat("vz",row) - (IsMC ? 0 : clas_alignement); // mm + //////////////////////////////////////// + /// compute electron resolution here + /// it depends en p and theta + /// the fine tuning will be done later + /// //////////////////////////////////// + //double px = recBank.getFloat("px",row); + //double py = recBank.getFloat("py",row); + //double pz = recBank.getFloat("pz",row); + //double p = Math.sqrt(px*px+py*py+pz*pz); + //double theta = Math.acos(pz/p); + vertex_resolutions[0] = 0.09; + vertex_resolutions[1] = 64;//4 + 1e10*theta + 1e10*p; + } + row++; + } + } // Loop over tracks int trackId = 0; @@ -53,7 +81,7 @@ private void propagation(ArrayList tracks, DataEvent event, final double // Initialize state vector double x0 = 0.0; double y0 = 0.0; - double z0 = track.get_Z0(); + double z0 = IsVtxDefined ? vz_constraint : track.get_Z0(); double px0 = track.get_px(); double py0 = track.get_py(); double pz0 = track.get_pz(); @@ -62,9 +90,6 @@ private void propagation(ArrayList tracks, DataEvent event, final double // Read list of hits ArrayList AHDC_hits = track.getHits(); Collections.sort(AHDC_hits); // sorted following the compareTo() method in Hit.java - - double zbeam = 0; - if(IsVtxDefined)zbeam = vz_constraint; // Start propagation Stepper stepper = new Stepper(y); @@ -76,7 +101,7 @@ private void propagation(ArrayList tracks, DataEvent event, final double RealVector initialStateEstimate = new ArrayRealVector(stepper.y); RealMatrix initialErrorCovariance = MatrixUtils.createRealMatrix(new double[][]{{50.0, 0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 50.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 900.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 100.00, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 100.00, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0, 900.0}}); KFitter TrackFitter = new KFitter(initialStateEstimate, initialErrorCovariance, stepper, propagator, materialHashMap); - TrackFitter.setVertexDefined(IsVtxDefined); + if (IsVtxDefined) TrackFitter.setVertexResolution(vertex_resolutions); // Loop over number of iterations for (int k = 0; k < Niter; k++) { @@ -94,7 +119,7 @@ private void propagation(ArrayList tracks, DataEvent event, final double } // Backward propagation (first layer to beamline) { - Hit hit = new Hit_beam(0, 0, zbeam); + Hit hit = new Hit_beam(0, 0, vz_constraint); TrackFitter.predict(hit, false); TrackFitter.correct(hit); } @@ -136,11 +161,11 @@ private void propagation(ArrayList tracks, DataEvent event, final double hit.setTrackId(trackId); sum_adc += hit.getADC(); sum_residuals += hit.getResidual(); - chi2 += Math.pow(hit.getResidual(),2.0); + chi2 += Math.pow(hit.getResidual(),2)/hit.get_MeasurementNoise().getEntry(0,0); } track.set_sum_adc(sum_adc); track.set_sum_residuals(sum_residuals); - track.set_chi2(chi2); + track.set_chi2(chi2/(AHDC_hits.size()-3)); track.set_p_drift_kf(p_drift); track.set_dEdx_kf(sum_adc/s); track.set_path_kf(s); diff --git a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java index b0c2a4992a..a5cdc2518e 100644 --- a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java +++ b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java @@ -187,11 +187,11 @@ public boolean processDataEvent(DataEvent event) { int matchHitId = -1; for (int k = 0; k < bank_ATOFHits.rows(); k++) { - int component = bank.getInt("component", k); + int component = bank_ATOFHits.getInt("component", k); if (component == 10) continue; - int sector = bank.getInt("sector", k); - int layer = bank.getInt("layer", k); + int sector = bank_ATOFHits.getInt("sector", k); + int layer = bank_ATOFHits.getInt("layer", k); ATOFHit hit = new ATOFHit(sector, layer, component, 0, 0, 0, 0f, ATOF); @@ -213,6 +213,7 @@ public boolean processDataEvent(DataEvent event) { System.out.println("Exception in ALERTEngine processDataEvent: " + ex); // TODO: proper logging } } + rbc.appendTrackMatchingAIBank(event, matched_ATOF_hit_id); return true; }