Implement the electron vertex for AHDC (#1077)

* implement electron vertex for ahdc

* convert vz_constraint to mm

* test alignement for clas and ahdc

* restore alignement at 0

* fix alert engine for ai track and atof matching

* fix alignement sign

* do not apply vertex shift for simulation

* add ahdc geometry verification in the test

* make the error on the ahdc distance dependent on adc and time

* correct the error formula

* correct the chi2 calculation

* fix chi2 computation

* remove redundancy

Author: ftouchte

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());
+					}
+				}
+			}
+		}
     }
 
 }

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<String, Material> 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];
+	}
 
 }

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<Track> tracks, DataEvent event, final double magfield, boolean IsMC) {
 
@@ -44,6 +48,30 @@ private void propagation(ArrayList<Track> tracks, DataEvent event, final double
 			final double      tesla             = 0.001;
 			final double[]    B                 = {0.0, 0.0, magfield / 10 * tesla};
 			HashMap<String, Material> 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<Track> 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<Track> tracks, DataEvent event, final double
 			    // Read list of hits
 			    ArrayList<Hit> 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<Track> 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<Track> 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<Track> 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);

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;
     }