From 2cb0c17d0b91bc0a129a4bc98699a422f28c02eb Mon Sep 17 00:00:00 2001 From: Raffaella De Vita Date: Mon, 16 Jun 2025 10:50:39 -0400 Subject: [PATCH 01/16] add Eloss for DDVCS --- .../src/main/java/cnuphys/rk4/RungeKutta.java | 54 ++++++++++++++++++- 1 file changed, 53 insertions(+), 1 deletion(-) diff --git a/common-tools/cnuphys/swimmer/src/main/java/cnuphys/rk4/RungeKutta.java b/common-tools/cnuphys/swimmer/src/main/java/cnuphys/rk4/RungeKutta.java index d1fb8ca365..5e2e2e12b8 100644 --- a/common-tools/cnuphys/swimmer/src/main/java/cnuphys/rk4/RungeKutta.java +++ b/common-tools/cnuphys/swimmer/src/main/java/cnuphys/rk4/RungeKutta.java @@ -24,7 +24,24 @@ public class RungeKutta { // the max dimension we'll use is probably 6, for state vectors // e.g., [x, y, z, px/p, py/p, pz/p]. private static int MAXDIM = 6; // we'll know if this fails! - + + // ddvcs parameters + //geometry + private static double[] THETA_SHIELD = {6.64, 41.3}; // min,max polar angle (deg) + private static double[] THETA_ECAL = {6.64, 30}; // min,max polar angle (deg) + private static double[] RHO_SHIELD = {60.0, 140}; // min,max distance from (0,0,0) at theta=25 deg + private static double[] RHO_ECAL = {60.0, 80.0}; // min,max distance from (0,0,0) at theta=25 deg + private static double TAN25 = Math.tan(Math.toRadians(25)); + private static double COS25 = Math.cos(Math.toRadians(25)); + // eloss + private static double EMASS = 0.511E-3; + private static double MUMASS = 105.66E-3; + private static double K = 0.000307075; // GeV mol-1 cm2 + // materials + private static double[] IeV = {823E-9, 600.7E-9}; // GeV + private static double[] DENSITY = {11.35, 8.3}; // g/cm3 + private static double[] ZOVERA = {0.39573, 0.41315}; + /** * Create a RungeKutta object that can be used for integration */ @@ -1304,5 +1321,40 @@ public double getMinStepSize() { return _minStepSize; } + + /** + * Calculate energy loss + * @param p momentum in GeV + * @param x position in meters + * @param y position in meters + * @param z position in meters + * @param dx step pathlength in meters + * @return + */ + public double getEloss(double p, double x, double y, double z, double dx) { + + double dz = Math.sqrt(x*x+y*y)*TAN25; + double r = (z+dz)*COS25; + if(rRHO_SHIELD[1]) + return 0; + + int imat = 0; + if(r Date: Tue, 17 Jun 2025 17:22:18 -0400 Subject: [PATCH 02/16] update swimmer with energy loss --- .../main/java/cnuphys/rk4/IDerivative.java | 5 +- .../src/main/java/cnuphys/rk4/RungeKutta.java | 68 +++------------ .../java/cnuphys/swim/DefaultDerivative.java | 86 ++++++++++++++++++- .../java/cnuphys/swim/SwimTrajectory.java | 18 ++++ .../src/main/java/cnuphys/swim/Swimmer.java | 2 + .../java/org/jlab/clas/swimtools/Swim.java | 8 ++ 6 files changed, 126 insertions(+), 61 deletions(-) diff --git a/common-tools/cnuphys/swimmer/src/main/java/cnuphys/rk4/IDerivative.java b/common-tools/cnuphys/swimmer/src/main/java/cnuphys/rk4/IDerivative.java index c8c0f40036..5b9c72e738 100644 --- a/common-tools/cnuphys/swimmer/src/main/java/cnuphys/rk4/IDerivative.java +++ b/common-tools/cnuphys/swimmer/src/main/java/cnuphys/rk4/IDerivative.java @@ -27,8 +27,5 @@ public interface IDerivative { */ public void derivative(double t, double y[], double dydt[]); - - - - + public void energyLossUpdate(double[] u, double dx); } diff --git a/common-tools/cnuphys/swimmer/src/main/java/cnuphys/rk4/RungeKutta.java b/common-tools/cnuphys/swimmer/src/main/java/cnuphys/rk4/RungeKutta.java index 5e2e2e12b8..9484f9d9cc 100644 --- a/common-tools/cnuphys/swimmer/src/main/java/cnuphys/rk4/RungeKutta.java +++ b/common-tools/cnuphys/swimmer/src/main/java/cnuphys/rk4/RungeKutta.java @@ -25,23 +25,6 @@ public class RungeKutta { // e.g., [x, y, z, px/p, py/p, pz/p]. private static int MAXDIM = 6; // we'll know if this fails! - // ddvcs parameters - //geometry - private static double[] THETA_SHIELD = {6.64, 41.3}; // min,max polar angle (deg) - private static double[] THETA_ECAL = {6.64, 30}; // min,max polar angle (deg) - private static double[] RHO_SHIELD = {60.0, 140}; // min,max distance from (0,0,0) at theta=25 deg - private static double[] RHO_ECAL = {60.0, 80.0}; // min,max distance from (0,0,0) at theta=25 deg - private static double TAN25 = Math.tan(Math.toRadians(25)); - private static double COS25 = Math.cos(Math.toRadians(25)); - // eloss - private static double EMASS = 0.511E-3; - private static double MUMASS = 105.66E-3; - private static double K = 0.000307075; // GeV mol-1 cm2 - // materials - private static double[] IeV = {823E-9, 600.7E-9}; // GeV - private static double[] DENSITY = {11.35, 8.3}; // g/cm3 - private static double[] ZOVERA = {0.39573, 0.41315}; - /** * Create a RungeKutta object that can be used for integration */ @@ -580,6 +563,7 @@ private int driver(double yo[], // typically [x, y, z, vx, vy, vz] and derivative double yt[] = new double[nDim]; double dydt[] = new double[nDim]; + double yttemp[] = new double[nDim]; double t = to; for (int i = 0; i < nDim; i++) { @@ -587,10 +571,18 @@ private int driver(double yo[], } for (int k = 1; k < nstep; k++) { + for (int i = 0; i < nDim; i++) { + yttemp[i] = yt[i]; + } + // use derivs at previous t deriv.derivative(t, yt, dydt); - - advancer.advance(t, yt, dydt, h, deriv, yt, null); + + advancer.advance(t, yt, dydt, h, deriv, yt, null); // yt is updated + + // Calculate energy loss, update momentum and alpha, and accumulate energy loss into totalEnergyLoss + deriv.energyLossUpdate(yttemp, h); + t += h; // someone listening? @@ -1035,7 +1027,7 @@ private int driver(double yo[], hdata[1] += h; hdata[2] = Math.max(hdata[2], h); } - + for (int i = 0; i < nDim; i++) { yt[i] = yt2[i]; } @@ -1320,41 +1312,5 @@ public double getMaxStepSize() { public double getMinStepSize() { return _minStepSize; } - - /** - * Calculate energy loss - * @param p momentum in GeV - * @param x position in meters - * @param y position in meters - * @param z position in meters - * @param dx step pathlength in meters - * @return - */ - public double getEloss(double p, double x, double y, double z, double dx) { - - double dz = Math.sqrt(x*x+y*y)*TAN25; - double r = (z+dz)*COS25; - if(rRHO_SHIELD[1]) - return 0; - - int imat = 0; - if(rRHO_SHIELD[1] || thetaTHETA_SHIELD[1]) + return 0; + + int imat = 0; + if(r_cm implements Serializable /** The source of the trajectory e.g. hbtracking */ private String _source = "???"; + + private double totalEnergyLoss = 0; + + /** + * Set total energy loss + * @param totalEnergyLoss total energy loss in GeV + */ + public void setTotalEnergyLoss(double totalEnergyLoss){ + this.totalEnergyLoss = totalEnergyLoss; + } + + /** + * Get total energy loss + * @return total energy loss + */ + public double getTotalEnergyLoss(){ + return totalEnergyLoss; + } /** * Create a swim trajectory with no initial content diff --git a/common-tools/cnuphys/swimmer/src/main/java/cnuphys/swim/Swimmer.java b/common-tools/cnuphys/swimmer/src/main/java/cnuphys/swim/Swimmer.java index c170bebcf8..733c64dc45 100644 --- a/common-tools/cnuphys/swimmer/src/main/java/cnuphys/swim/Swimmer.java +++ b/common-tools/cnuphys/swimmer/src/main/java/cnuphys/swim/Swimmer.java @@ -845,6 +845,8 @@ public SwimTrajectory swim(int charge, double xo, double yo, double zo, double m // Integrate DefaultDerivative deriv = new DefaultDerivative(charge, momentum, _probe); ntotal = (new RungeKutta()).uniformStep(uo, 0, maxPathLength, u, s, deriv, stopper); + + trajectory.setTotalEnergyLoss(deriv.getTotalEnergyLoss()); // now cycle through and get the save points for (int i = 0; i < ntotal; i++) { diff --git a/common-tools/swim-tools/src/main/java/org/jlab/clas/swimtools/Swim.java b/common-tools/swim-tools/src/main/java/org/jlab/clas/swimtools/Swim.java index b0e45126e9..e0252320c2 100644 --- a/common-tools/swim-tools/src/main/java/org/jlab/clas/swimtools/Swim.java +++ b/common-tools/swim-tools/src/main/java/org/jlab/clas/swimtools/Swim.java @@ -52,6 +52,8 @@ public class Swim { private ProbeCollection PC; + private static double MUMASS = 105.66E-3; + /** * Class for swimming to various surfaces. The input and output units are cm and GeV/c */ @@ -969,6 +971,12 @@ public double[] SwimToBeamLine(double xB, double yB) { return null; st.computeBDL(PC.CP); // st.computeBDL(compositeField); + + // Get total energy loss and update final momentum + double totalEnergyLoss = st.getTotalEnergyLoss(); + double energyOrigin = Math.sqrt(_pTot*_pTot + MUMASS*MUMASS); + double energyFinal = energyOrigin + totalEnergyLoss; + _pTot = Math.sqrt(energyFinal*energyFinal - MUMASS*MUMASS); double[] lastY = st.lastElement(); From 3d7524658eb3323203a67fd69eeb4a2aeb1655f9 Mon Sep 17 00:00:00 2001 From: tongtongcao Date: Wed, 18 Jun 2025 13:00:36 -0400 Subject: [PATCH 03/16] override energyLossUpdate() in daughter classes --- .../cnuphys/swimmer/src/main/java/cnuphys/rk4/RkTest.java | 3 +++ .../swimmer/src/main/java/cnuphys/swim/SectorDerivative.java | 3 +++ .../swimmer/src/main/java/cnuphys/swimZ/SwimZDerivative.java | 3 +++ 3 files changed, 9 insertions(+) diff --git a/common-tools/cnuphys/swimmer/src/main/java/cnuphys/rk4/RkTest.java b/common-tools/cnuphys/swimmer/src/main/java/cnuphys/rk4/RkTest.java index 0fd048169d..390097fd0c 100644 --- a/common-tools/cnuphys/swimmer/src/main/java/cnuphys/rk4/RkTest.java +++ b/common-tools/cnuphys/swimmer/src/main/java/cnuphys/rk4/RkTest.java @@ -14,6 +14,9 @@ public void derivative(double t, double[] y, double[] dydt) { dydt[0] = y[1]; dydt[1] = y[0]; } + + @Override + public void energyLossUpdate(double[] u, double dx){} }; diff --git a/common-tools/cnuphys/swimmer/src/main/java/cnuphys/swim/SectorDerivative.java b/common-tools/cnuphys/swimmer/src/main/java/cnuphys/swim/SectorDerivative.java index 01d9986cb7..e99e09c6f6 100644 --- a/common-tools/cnuphys/swimmer/src/main/java/cnuphys/swim/SectorDerivative.java +++ b/common-tools/cnuphys/swimmer/src/main/java/cnuphys/swim/SectorDerivative.java @@ -97,5 +97,8 @@ public void derivative(double s, double[] Q, double[] dQds) { dQds[1] = Q[4]; dQds[2] = Q[5]; } + + @Override + public void energyLossUpdate(double[] u, double dx){} } diff --git a/common-tools/cnuphys/swimmer/src/main/java/cnuphys/swimZ/SwimZDerivative.java b/common-tools/cnuphys/swimmer/src/main/java/cnuphys/swimZ/SwimZDerivative.java index b0cc781def..6978589066 100644 --- a/common-tools/cnuphys/swimmer/src/main/java/cnuphys/swimZ/SwimZDerivative.java +++ b/common-tools/cnuphys/swimmer/src/main/java/cnuphys/swimZ/SwimZDerivative.java @@ -89,5 +89,8 @@ public void derivative(double z, double[] x, double[] dxdz) { dxdz[2] = _qv * Ax; dxdz[3] = _qv * Ay; } + + @Override + public void energyLossUpdate(double[] u, double dx){} } From bf47759949865492cac8cb82e2cb15c6b988a08d Mon Sep 17 00:00:00 2001 From: tongtongcao Date: Wed, 18 Jun 2025 15:23:02 -0400 Subject: [PATCH 04/16] let energy loss only apply into swimming to beam line after TB tracking --- .../src/main/java/cnuphys/rk4/RungeKutta.java | 205 +++++++++++++++++- .../src/main/java/cnuphys/swim/Swimmer.java | 86 ++++++++ .../java/org/jlab/clas/swimtools/Swim.java | 2 +- .../rec/dc/track/TrackCandListFinder.java | 4 +- 4 files changed, 294 insertions(+), 3 deletions(-) diff --git a/common-tools/cnuphys/swimmer/src/main/java/cnuphys/rk4/RungeKutta.java b/common-tools/cnuphys/swimmer/src/main/java/cnuphys/rk4/RungeKutta.java index 9484f9d9cc..99764fb350 100644 --- a/common-tools/cnuphys/swimmer/src/main/java/cnuphys/rk4/RungeKutta.java +++ b/common-tools/cnuphys/swimmer/src/main/java/cnuphys/rk4/RungeKutta.java @@ -115,6 +115,89 @@ public void nextStep(double tNext, double yNext[], double h) { } + /** + * Driver that uses the RungeKutta advance with a uniform step size. (i.e., + * this does NOT use an adaptive step size.) + * + * This version stores each step into the arrays t[] and y[][]. An + * alternative does not store the results but instead uses an IRk4Listener + * to notify the listener that the next step has been advanced. + * + * A very typical case is a 2nd order ODE converted to a 1st order where the + * dependent variables are x, y, z, vx, vy, vz and the independent variable + * is time. + * + * @param yo + * initial values. Probably something like (xo, yo, zo, vxo, vyo, + * vzo). + * @param to + * the initial value of the independent variable, e.g., time. + * @param tf + * the maximum value of the independent variable. + * @param y + * will be filled with results. The first index is small-- the + * dimensionality of the problem-- i.e., often it is 6 for (xo, + * yo, zo, vxo, vyo, vzo). The second dimension is for storing + * results and determining stepsize. If it is 1000, we will have + * a thousand steps and the stepsize will be (tf-to)/1000 + * (actually 999). + * @param t + * will filled with the locations of t--should have the exact + * same large dimension as the second index of y--i.e., something + * like 1000. + * @param deriv + * the derivative computer (interface). This is where the problem + * specificity resides. + * @param stopper + * if not null will be used to exit the integration + * early because some condition has been reached. + * @return the number of steps used--may be less than the space provided if + * the integration ended early as a result of an exit condition. + */ + public int uniformStepWithEnergyLoss(double yo[], + double to, + double tf, + final double y[][], + final double t[], + IDerivative deriv, + IStopper stopper) { + int nstep = t.length; // the number of steps to store + + // the dimensionality of the problem, e.g. six if (x, y, z, vx, vy, vz) + final int nDim = yo.length; + + // uniform step size + double h = (tf - to) / (nstep - 1); + + // put starting step in + t[0] = to; + for (int i = 0; i < nDim; i++) { + y[i][0] = yo[i]; + } + + IRkListener listener = new IRkListener() { + + int step = 1; + + @Override + public void nextStep(double tNext, double yNext[], double h) { + + if (step < t.length) { + t[step] = tNext; + for (int i = 0; i < nDim; i++) { + // store results for this step + y[i][step] = yNext[i]; + } + } + step++; + } + + }; + + return uniformStepWithEnergyLoss(yo, to, tf, h, deriv, stopper, listener); + + } + /** * Integrator that uses the standard RK4 advance with a uniform step size. @@ -155,6 +238,46 @@ public int uniformStep(double yo[], UniformAdvance advancer = new UniformAdvance(); return driver(yo, to, tf, h, deriv, stopper, listener, advancer); } + + /** + * Integrator that uses the standard RK4 advance with a uniform step size. + * (i.e., this does NOT use an adaptive step size.) + * + * This version uses an IRk4Listener to notify the listener that the next + * step has been advanced. + * + * A very typical case is a 2nd order ODE converted to a 1st order where the + * dependent variables are x, y, z, vx, vy, vz and the independent variable + * is time. + * + * @param yo + * initial values. Probably something like (xo, yo, zo, vxo, vyo, + * vzo). + * @param to + * the initial value of the independent variable, e.g., time. + * @param tf + * the maximum value of the independent variable. + * @param deriv + * the derivative computer (interface). This is where the problem + * specificity resides. + * @param stopper + * if not null will be used to exit the integration + * early because some condition has been reached. + * @param listener + * listens for each step + * @return the number of steps used. + */ + public int uniformStepWithEnergyLoss(double yo[], + double to, + double tf, + double h, + IDerivative deriv, + IStopper stopper, + IRkListener listener) { + + UniformAdvance advancer = new UniformAdvance(); + return driverWithEnergyLoss(yo, to, tf, h, deriv, stopper, listener, advancer); + } /** * Integrator that uses the RungeKutta advance with a Butcher Tableau and @@ -570,6 +693,86 @@ private int driver(double yo[], yt[i] = yo[i]; } + for (int k = 1; k < nstep; k++) { + for (int i = 0; i < nDim; i++) { + yttemp[i] = yt[i]; + } + + // use derivs at previous t + deriv.derivative(t, yt, dydt); + + advancer.advance(t, yt, dydt, h, deriv, yt, null); // yt is updated + + t += h; + + // someone listening? + if (listener != null) { + listener.nextStep(t, yt, h); + } + + // premature termination? Skip if stopper is null. + if (stopper != null) { + stopper.setFinalT(t); + if (stopper.stopIntegration(t, yt)) { + return k + 1; // actual number of steps taken + } + } + } + + return nstep; + } + + + /** + * Driver that uses the RungeKutta advance with a uniform step size. (I.e., + * this does NOT use an adaptive step size.) + * + * This version uses an IRk4Listener to notify the listener that the next + * step has been advanced. + * + * A very typical case is a 2nd order ODE converted to a 1st order where the + * dependent variables are x, y, z, vx, vy, vz and the independent variable + * is time. + * + * @param yo + * initial values. Probably something like (xo, yo, zo, vxo, vyo, + * vzo). + * @param to + * the initial value of the independent variable, e.g., time. + * @param tf + * the maximum value of the independent variable. + * @param deriv + * the derivative computer (interface). This is where the problem + * specificity resides. + * @param stopper + * if not null will be used to exit the integration + * early because some condition has been reached. + * @return the number of steps used. + */ + private int driverWithEnergyLoss(double yo[], + double to, + double tf, + double h, + IDerivative deriv, + IStopper stopper, + IRkListener listener, + IAdvance advancer) { + int nstep = (int) (1 + (tf - to) / h); // the number of steps to store + + // the dimensionality of the problem. E.., 6 if (x, y, z, vx, vy, vz) + int nDim = yo.length; + + // yt is the current value of the state vector, + // typically [x, y, z, vx, vy, vz] and derivative + double yt[] = new double[nDim]; + double dydt[] = new double[nDim]; + double yttemp[] = new double[nDim]; + + double t = to; + for (int i = 0; i < nDim; i++) { + yt[i] = yo[i]; + } + for (int k = 1; k < nstep; k++) { for (int i = 0; i < nDim; i++) { yttemp[i] = yt[i]; @@ -600,7 +803,7 @@ private int driver(double yo[], } return nstep; - } + } diff --git a/common-tools/cnuphys/swimmer/src/main/java/cnuphys/swim/Swimmer.java b/common-tools/cnuphys/swimmer/src/main/java/cnuphys/swim/Swimmer.java index 733c64dc45..a45afd3ea7 100644 --- a/common-tools/cnuphys/swimmer/src/main/java/cnuphys/swim/Swimmer.java +++ b/common-tools/cnuphys/swimmer/src/main/java/cnuphys/swim/Swimmer.java @@ -858,6 +858,92 @@ public SwimTrajectory swim(int charge, double xo, double yo, double zo, double m return trajectory; } + + + /** + * Swims a Lund particle with a built it stopper for the maximum value of + * the radial coordinate. This is for the trajectory mode, where you want to + * cache steps along the path. Uses a fixed stepsize algorithm. + * + * @param charge + * the charge: -1 for electron, 1 for proton, etc + * @param xo + * the x vertex position in meters + * @param yo + * the y vertex position in meters + * @param zo + * the z vertex position in meters + * @param momentum + * initial momentum in GeV/c + * @param theta + * initial polar angle in degrees + * @param phi + * initial azimuthal angle in degrees + * @param stopper + * an optional object that can terminate the swimming based on + * some condition + * @param maxPathLength + * in meters. This determines the max number of steps based on + * the step size. If a stopper is used, the integration might + * terminate before all the steps are taken. A reasonable value + * for CLAS is 8. meters + * @param stepSize + * the uniform step size in meters. + * @param distanceBetweenSaves + * this distance is in meters. It should be bigger than stepSize. + * It is approximately the distance between "saves" where the + * point is saved in a trajectory for later drawing. + * @return the trajectory of the particle + */ + public SwimTrajectory swimWithEnergyLoss(int charge, double xo, double yo, double zo, double momentum, double theta, double phi, + IStopper stopper, double maxPathLength, double stepSize, double distanceBetweenSaves) { + + // if no magnetic field or no charge, then simple straight line tracks. + // the path will consist of just two points + if ((_probe == null) || (charge == 0)) { + GeneratedParticleRecord genPartRec = new GeneratedParticleRecord(charge, xo, yo, zo, momentum, theta, phi); + return straightLineTrajectory(genPartRec, maxPathLength); + } + + if (momentum < MINMOMENTUM) { + //System.err.println("Skipping low momentum swim (A)"); + return new SwimTrajectory(charge, xo, yo, zo, momentum, theta, phi); + } + + // cycle is the number of advances per save + int cycle = (int) (distanceBetweenSaves / stepSize); + cycle = Math.max(2, cycle); + + // max number of possible steps--may not use all of them + int ntotal = (int) (maxPathLength / stepSize); // number steps + int nsave = ntotal / cycle; // aprox number saves + + // the the initial six vector + double uo[] = initialState(xo, yo, zo, theta, phi); + + // storage for time and state + double s[] = new double[ntotal]; + double u[][] = new double[6][ntotal]; + + // create the trajectory container + SwimTrajectory trajectory = new SwimTrajectory(charge, xo, yo, zo, momentum, theta, phi, nsave); + + // Integrate + DefaultDerivative deriv = new DefaultDerivative(charge, momentum, _probe); + ntotal = (new RungeKutta()).uniformStepWithEnergyLoss(uo, 0, maxPathLength, u, s, deriv, stopper); + + trajectory.setTotalEnergyLoss(deriv.getTotalEnergyLoss()); + + // now cycle through and get the save points + for (int i = 0; i < ntotal; i++) { + if (((i % cycle) == 0) || (i == (ntotal - 1))) { + double v[] = makeVector(u[0][i], u[1][i], u[2][i], u[3][i], u[4][i], u[5][i]); + trajectory.add(v); + } + } + + return trajectory; + } /** * Swims a Lund particle with a built in stopper for the maximum value of diff --git a/common-tools/swim-tools/src/main/java/org/jlab/clas/swimtools/Swim.java b/common-tools/swim-tools/src/main/java/org/jlab/clas/swimtools/Swim.java index e0252320c2..68995a0c70 100644 --- a/common-tools/swim-tools/src/main/java/org/jlab/clas/swimtools/Swim.java +++ b/common-tools/swim-tools/src/main/java/org/jlab/clas/swimtools/Swim.java @@ -965,7 +965,7 @@ public double[] SwimToBeamLine(double xB, double yB) { return null; BeamLineSwimStopper stopper = new BeamLineSwimStopper(xB, yB); - SwimTrajectory st = PC.CF.swim(_charge, _x0, _y0, _z0, _pTot, _theta, _phi, stopper, _maxPathLength, stepSize, + SwimTrajectory st = PC.CF.swimWithEnergyLoss(_charge, _x0, _y0, _z0, _pTot, _theta, _phi, stopper, _maxPathLength, stepSize, 0.0005); if(st==null) return null; diff --git a/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/TrackCandListFinder.java b/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/TrackCandListFinder.java index 6e8ca4cb76..8ec3b8a7d1 100644 --- a/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/TrackCandListFinder.java +++ b/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/TrackCandListFinder.java @@ -561,7 +561,7 @@ public void setTrackPars(Track cand, cand.fit_Successful = false; return; } - + // recalc new vertex using plane stopper //int sector = cand.get(2).getSector(); //double[] Vt = null; @@ -592,6 +592,7 @@ public void setTrackPars(Track cand, cand.fit_Successful = true; cand.set_TrackingInfoString(trking); + /* dcSwim.SetSwimParameters(xOrFix, yOrFix, zOrFix, pxOrFix, pyOrFix, pzOrFix, cand.get_Q()); @@ -610,6 +611,7 @@ public void setTrackPars(Track cand, //set the pseudocross at extrapolated position cand.set_PreRegion1CrossPoint(new Point3D(xInner, yInner, zInner)); cand.set_PreRegion1CrossDir(new Point3D(uxInner, uyInner, uzInner)); + */ } private Integer getKey(Track trk) { From 20e2ba4ddd83b82f710626ff3237d5a61b73e50e Mon Sep 17 00:00:00 2001 From: tongtongcao Date: Wed, 18 Jun 2025 15:32:57 -0400 Subject: [PATCH 05/16] cancel swimming from vertex to first layer of DC in TrackCandListFinder::matchHits() --- .../main/java/org/jlab/rec/dc/track/TrackCandListFinder.java | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/TrackCandListFinder.java b/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/TrackCandListFinder.java index 8ec3b8a7d1..0fef7c1f76 100644 --- a/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/TrackCandListFinder.java +++ b/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/TrackCandListFinder.java @@ -757,7 +757,7 @@ public void matchHits(List stateVecAtPlanesList, Track trk, } List fhits = new ArrayList<>(); - + /* dcSwim.SetSwimParameters(trk.get_Vtx0().x(), trk.get_Vtx0().y(), trk.get_Vtx0().z(), trk.get_pAtOrig().x(), trk.get_pAtOrig().y(), trk.get_pAtOrig().z(), trk.get_Q()); @@ -766,6 +766,7 @@ public void matchHits(List stateVecAtPlanesList, Track trk, if (ToFirstMeas == null) { return; } + */ for (StateVec st : stateVecAtPlanesList) { if (st == null) { From 26bf391fb32adc53330c332ef2c95cc2549b068d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mariangela=20Bond=C3=AD?= Date: Thu, 12 Jun 2025 15:39:31 +0200 Subject: [PATCH 06/16] new branch for urwell ddvcs proposal --- .../geant4/v2/URWELL/URWellConstants.java | 29 +++- .../geant4/v2/URWELL/URWellGeant4Factory.java | 156 +++++++++++++----- 2 files changed, 141 insertions(+), 44 deletions(-) diff --git a/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellConstants.java b/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellConstants.java index 81b08bf569..c41492485f 100644 --- a/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellConstants.java +++ b/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellConstants.java @@ -9,11 +9,12 @@ public class URWellConstants { private final static String CCDBPATH = "/geometry/urwell/"; - public final static int NMAXREGIONS = 2; //max number of regions + public final static int NMAXREGIONS = 6; //max number of regions public final static int NREGIONS = 1; //number of regions public final static int NSECTORS = 6; //number of sectors public final static int NLAYERS = 2; //number of layers public final static int NCHAMBERS = 3; //number of chambers in a sector + public final static int NCHAMBERS_ddvcs = 1; //number of chambers in a sector public final static double XENLARGEMENT = 0.5; // cm public final static double YENLARGEMENT = 1.; // cm @@ -26,6 +27,10 @@ public class URWellConstants { public final static double SECTORHEIGHT = 146.21; //height of each sector (cm) public final static double DX0CHAMBER0 = 5.197; // halfbase of chamber 1 (cm) + public final static double SECTORHEIGHT_ddvcs = 22; + public final static double THMIN_ddvcs = 7; + public final static double DX0CHAMBER0_ddvcs = 3.597; + // Chamber volumes and materials (units are cm) public final static double[] CHAMBERVOLUMESTHICKNESS = {0.0025, 0.0005,0.3, // window 0.0025, 0.0005,0.4, // cathode @@ -43,7 +48,7 @@ public class URWellConstants { "support_skin1_g10", "support_honeycomb_nomex", "support_skin2_g10"}; // URWELL position in the CLAS12 frame - public final static double TGT2DC0 = 228.078; // cm + public final static double TGT2DC0 = 228.078; // cm // public final static double URWELL2DC0 = 2; // cm public final static double URWELL2DC0[] = new double[NMAXREGIONS]; public final static double DIST2TGT[] = new double[NMAXREGIONS]; @@ -51,12 +56,22 @@ public class URWellConstants { public final static double YMIN[] = new double[NMAXREGIONS]; public final static double ZMIN[] = new double[NMAXREGIONS]; + + public final static double TGT2DC0_ddvcs = 52.9; // cm + // public final static double URWELL2DC0 = 2; // cm + public final static double URWELL2DC0_ddvcs[] = new double[NMAXREGIONS]; + public final static double DIST2TGT_ddvcs[] = new double[NMAXREGIONS]; + public final static double W2TGT_ddvcs[] = new double[NMAXREGIONS]; + ; + public final static double YMIN_ddvcs[] = new double[NMAXREGIONS]; + public final static double ZMIN_ddvcs[] = new double[NMAXREGIONS]; + // public final static double DIST2TGT = (TGT2DC0-URWELL2DC0); // public final static double W2TGT = DIST2TGT/Math.cos(Math.toRadians(THTILT-THMIN)); // public final static double YMIN = W2TGT*Math.sin(Math.toRadians(THMIN)); // distance from the base chamber1 and beamline // public final static double ZMIN = W2TGT*Math.cos(Math.toRadians(THMIN)); public final static double PITCH = 0.1 ; // cm - public final static double STEREOANGLE = 10; // deg + public final static double STEREOANGLE = 6; // deg @@ -121,7 +136,13 @@ public static synchronized void load( DatabaseConstantProvider cp ) W2TGT[i] = DIST2TGT[i]/Math.cos(Math.toRadians(THTILT-THMIN)); YMIN[i]= W2TGT[i]*Math.sin(Math.toRadians(THMIN)); // distance from the base chamber1 and beamline ZMIN[i] = W2TGT[i]*Math.cos(Math.toRadians(THMIN)); - + + URWELL2DC0_ddvcs[i] = -2. + i * 1.3; + DIST2TGT_ddvcs[i] = (TGT2DC0_ddvcs + URWELL2DC0_ddvcs[i]); + W2TGT_ddvcs[i] = DIST2TGT_ddvcs[i] / Math.cos(Math.toRadians(THTILT - THMIN_ddvcs)); + YMIN_ddvcs[i] = W2TGT_ddvcs[i] * Math.sin(Math.toRadians(THMIN_ddvcs)); // distance from the base chamber1 and beamline + ZMIN_ddvcs[i] = W2TGT_ddvcs[i] * Math.cos(Math.toRadians(THMIN_ddvcs)); + } diff --git a/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellGeant4Factory.java b/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellGeant4Factory.java index deaac1cf2d..6d128ded5d 100644 --- a/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellGeant4Factory.java +++ b/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellGeant4Factory.java @@ -8,6 +8,8 @@ import org.jlab.detector.calib.utils.DatabaseConstantProvider; import org.jlab.geom.prim.Line3D; import org.jlab.geom.prim.Point3D; +import java.util.Arrays; + /** @@ -20,7 +22,13 @@ public final class URWellGeant4Factory extends Geant4Factory { private int nRegions = URWellConstants.NREGIONS; private int nSectors = URWellConstants.NSECTORS; private int nChambers = URWellConstants.NCHAMBERS; + private double sectorheight = URWellConstants.SECTORHEIGHT; + private double thilt = URWellConstants.THTILT; + private double dx0chamber0 = URWellConstants.DX0CHAMBER0; + private double thopen = URWellConstants.THOPEN; private boolean isProto = false; + private double[] ymin = new double[URWellConstants.NMAXREGIONS]; + private double[] zmin = new double[URWellConstants.NMAXREGIONS]; /** @@ -34,12 +42,91 @@ public URWellGeant4Factory( DatabaseConstantProvider cp, boolean prototype, int this.init(cp, prototype, nRegions); } + /** + * + * @param config is 0 standard, 1 ddvcs proposal, 2 prototype + */ + + public URWellGeant4Factory( DatabaseConstantProvider cp, int config, int nRegions) { + URWellConstants.connect(cp ); + this.init(cp, config, nRegions); + } + + + public void init(DatabaseConstantProvider cp, int config, int regions ) { + + motherVolume = new G4World("root"); + + switch (config) { + case 0 -> { + + nSectors = URWellConstants.NSECTORS; + nChambers = URWellConstants.NCHAMBERS; + nRegions = Math.min(URWellConstants.NMAXREGIONS, regions); + isProto = false; + sectorheight = URWellConstants.SECTORHEIGHT; + thilt = URWellConstants.THTILT; + dx0chamber0 = URWellConstants.DX0CHAMBER0; + thopen = URWellConstants.THOPEN; + ymin =Arrays.copyOf(URWellConstants.YMIN, URWellConstants.YMIN.length); + zmin =Arrays.copyOf(URWellConstants.ZMIN, URWellConstants.ZMIN.length); + + } + case 1 -> { + nRegions = Math.min(URWellConstants.NMAXREGIONS, regions); + nSectors = URWellConstants.NSECTORS; + nChambers = URWellConstants.NCHAMBERS_ddvcs; + sectorheight = URWellConstants.SECTORHEIGHT_ddvcs; + thilt = URWellConstants.THTILT; + dx0chamber0 = URWellConstants.DX0CHAMBER0_ddvcs; + thopen = URWellConstants.THOPEN; + ymin = Arrays.copyOf(URWellConstants.YMIN_ddvcs, URWellConstants.YMIN_ddvcs.length); + zmin = Arrays.copyOf(URWellConstants.ZMIN_ddvcs, URWellConstants.ZMIN_ddvcs.length); + + isProto = false; + } + case 2 -> { + nRegions = URWellConstants.NREGIONS_PROTO; + nSectors = URWellConstants.NSECTORS_PROTO; + nChambers = URWellConstants.NCHAMBERS_PROTO; + isProto = true; + + } + default -> { + } + } + + for (int iregion = 0; iregion { System.out.println(volume.gemcString()); }); From ea0419a3128241e99cb32e9c344b8a6dd760239f Mon Sep 17 00:00:00 2001 From: tongtongcao Date: Fri, 20 Jun 2025 14:52:05 -0400 Subject: [PATCH 07/16] Change logger level from FINE to FINEST for DC codes at the event level --- .../kalmanfilter/zReference/StateVecs.java | 6 +- .../java/org/jlab/rec/dc/banks/HitReader.java | 14 ++-- .../java/org/jlab/rec/dc/cross/Cross.java | 10 +-- .../java/org/jlab/rec/dc/nn/PatternRec.java | 2 +- .../java/org/jlab/rec/dc/segment/Segment.java | 4 +- .../rec/dc/track/TrackCandListFinder.java | 64 +++++++++---------- .../jlab/rec/dc/track/fit/KFitterDoca.java | 16 ++--- .../jlab/rec/dc/track/fit/MeasVecsDoca.java | 8 +-- .../jlab/rec/dc/track/fit/StateVecsDoca.java | 12 ++-- .../java/org/jlab/service/dc/DCEngine.java | 2 +- .../jlab/service/dc/DCHBPostClusterAI.java | 6 +- .../jlab/service/dc/DCHBPostClusterConv.java | 8 +-- .../java/org/jlab/service/dc/DCTBEngine.java | 10 +-- 13 files changed, 81 insertions(+), 81 deletions(-) diff --git a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/zReference/StateVecs.java b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/zReference/StateVecs.java index 269e966c21..1fe0751241 100644 --- a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/zReference/StateVecs.java +++ b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/zReference/StateVecs.java @@ -97,7 +97,7 @@ public Matrix transport(int sector, int i, double Zf, StateVec iVec, AMeasVecs m Matrix5x5.copy(fVec.CM, cMat); s = Math.signum(Zf - zInit) * stepSize; - // LOGGER.log(Level.FINE, " from "+(float)Z[i]+" to "+(float)Z[f]+" at "+(float)z+" By is "+bf[1]+" B is "+Math.sqrt(bf[0]*bf[0]+bf[1]*bf[1]+bf[2]*bf[2])/Bmax+" stepSize is "+s); + // LOGGER.log(Level.FINEST, " from "+(float)Z[i]+" to "+(float)Z[f]+" at "+(float)z+" By is "+bf[1]+" B is "+Math.sqrt(bf[0]*bf[0]+bf[1]*bf[1]+bf[2]*bf[2])/Bmax+" stepSize is "+s); if (Math.signum(Zf - zInit) * (z + s) > Math.signum(Zf - zInit) * Zf) { s = Math.signum(Zf - zInit) * Math.abs(Zf - z); } @@ -195,7 +195,7 @@ public double getDeltaPath(int sector, int i, double Zf, StateVec iVec, AMeasVec Matrix5x5.copy(fVec.CM, cMat); s = Math.signum(Zf - zInit) * stepSize; - // LOGGER.log(Level.FINE, " from "+(float)Z[i]+" to "+(float)Z[f]+" at "+(float)z+" By is "+bf[1]+" B is "+Math.sqrt(bf[0]*bf[0]+bf[1]*bf[1]+bf[2]*bf[2])/Bmax+" stepSize is "+s); + // LOGGER.log(Level.FINEST, " from "+(float)Z[i]+" to "+(float)Z[f]+" at "+(float)z+" By is "+bf[1]+" B is "+Math.sqrt(bf[0]*bf[0]+bf[1]*bf[1]+bf[2]*bf[2])/Bmax+" stepSize is "+s); if (Math.signum(Zf - zInit) * (z + s) > Math.signum(Zf - zInit) * Zf) { s = Math.signum(Zf - zInit) * Math.abs(Zf - z); } @@ -296,7 +296,7 @@ public boolean transport(int sector, int i, int f, StateVec iVec, AMeasVecs mv, Matrix5x5.copy(fVec.CM, cMat); s = Math.signum(zFinal - zInit) * stepSize; - // LOGGER.log(Level.FINE, " from "+(float)Z[i]+" to "+(float)Z[f]+" at "+(float)z+" By is "+bf[1]+" B is "+Math.sqrt(bf[0]*bf[0]+bf[1]*bf[1]+bf[2]*bf[2])/Bmax+" stepSize is "+s); + // LOGGER.log(Level.FINEST, " from "+(float)Z[i]+" to "+(float)Z[f]+" at "+(float)z+" By is "+bf[1]+" B is "+Math.sqrt(bf[0]*bf[0]+bf[1]*bf[1]+bf[2]*bf[2])/Bmax+" stepSize is "+s); if (Math.signum(zFinal - zInit) * (z + s) > Math.signum(zFinal - zInit) * zFinal) { s = Math.signum(zFinal - zInit) * Math.abs(zFinal - z); } diff --git a/reconstruction/dc/src/main/java/org/jlab/rec/dc/banks/HitReader.java b/reconstruction/dc/src/main/java/org/jlab/rec/dc/banks/HitReader.java index b2b4c30574..2dbac67675 100644 --- a/reconstruction/dc/src/main/java/org/jlab/rec/dc/banks/HitReader.java +++ b/reconstruction/dc/src/main/java/org/jlab/rec/dc/banks/HitReader.java @@ -352,7 +352,7 @@ private void read_HBHits(TimeToDistanceEstimator tde) { String pointName = bankNames.getInputIdsBank(); String recBankName = bankNames.getRecEventBank(); - LOGGER.log(Level.FINE,"Reading hb banks for "+ bankName + ", " + pointName + " " + recBankName); + LOGGER.log(Level.FINEST,"Reading hb banks for "+ bankName + ", " + pointName + " " + recBankName); _HBHits = new ArrayList<>(); @@ -489,7 +489,7 @@ private void read_HBHits(TimeToDistanceEstimator tde) { //if(hit.betaFlag == 0) if(passHit(hit.betaFlag)) { this._HBHits.add(hit); - LOGGER.log(Level.FINE, "Passing "+hit.printInfo()+" for "+ bankNames.getHitsBank()); + LOGGER.log(Level.FINEST, "Passing "+hit.printInfo()+" for "+ bankNames.getHitsBank()); } } } @@ -562,7 +562,7 @@ private void read_NNHits() { hit.NNTrkP = this.aimatch.get(clusterID)[0]; hit.NNTrkTheta = this.aimatch.get(clusterID)[1]; hit.NNTrkPhi = this.aimatch.get(clusterID)[2]; - LOGGER.log(Level.FINE, "NN"+hit.printInfo()); + LOGGER.log(Level.FINEST, "NN"+hit.printInfo()); this._DCHits.add(hit); } } @@ -788,11 +788,11 @@ private void swapWires(DataEvent event, int sector, int layer, int wire) { } else { for(int i = 0; i clusters, if(entry.getValue().size()==3) crossList.add(entry.getValue()); for(Cross c : entry.getValue()) - LOGGER.log(Level.FINE, "AI"+c.printInfo()+c.get_Segment1().printInfo()+c.get_Segment2().printInfo()); + LOGGER.log(Level.FINEST, "AI"+c.printInfo()+c.get_Segment1().printInfo()+c.get_Segment2().printInfo()); } return crossList; } diff --git a/reconstruction/dc/src/main/java/org/jlab/rec/dc/segment/Segment.java b/reconstruction/dc/src/main/java/org/jlab/rec/dc/segment/Segment.java index 793dfe0e5f..83f85af139 100644 --- a/reconstruction/dc/src/main/java/org/jlab/rec/dc/segment/Segment.java +++ b/reconstruction/dc/src/main/java/org/jlab/rec/dc/segment/Segment.java @@ -212,7 +212,7 @@ public boolean isCloseTo(Segment otherseg) { ///
|Xwires2-Xwires1| = a*Xwires1 + b
\n /// where a and b are DC parameters set by DC_RSEG_a and DC_RSEG_B .\n\n boolean value = false; - //LOGGER.log(Level.FINE, "in Segment DeltaW "+Math.abs(this.getAvgwire()-otherseg.getAvgwire() )+ + //LOGGER.log(Level.FINEST, "in Segment DeltaW "+Math.abs(this.getAvgwire()-otherseg.getAvgwire() )+ // " < ? "+(Constants.DC_RSEG_A * this.getAvgwire() + Constants.DC_RSEG_B)); if (Math.abs(this.getAvgwire() - otherseg.getAvgwire()) < Constants.DC_RSEG_A * this.getAvgwire() + Constants.DC_RSEG_B) { value = true; @@ -229,7 +229,7 @@ public boolean hasNoMatchingSegment(List othersegs) { ///
|Xwires2-Xwires1| = a*Xwires1 + b
\n /// where a and b are DC parameters set by DC_RSEG_a and DC_RSEG_B .\n\n boolean value = true; - //LOGGER.log(Level.FINE, "in Segment DeltaW "+Math.abs(this.getAvgwire()-otherseg.getAvgwire() )+ + //LOGGER.log(Level.FINEST, "in Segment DeltaW "+Math.abs(this.getAvgwire()-otherseg.getAvgwire() )+ // " < ? "+(Constants.DC_RSEG_A * this.getAvgwire() + Constants.DC_RSEG_B)); for(Segment otherseg : othersegs) { if (Math.abs(this.getAvgwire() - otherseg.getAvgwire()) < Constants.DC_RSEG_A * this.getAvgwire() + Constants.DC_RSEG_B) { diff --git a/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/TrackCandListFinder.java b/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/TrackCandListFinder.java index 6e8ca4cb76..c3d248a1bf 100644 --- a/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/TrackCandListFinder.java +++ b/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/TrackCandListFinder.java @@ -351,7 +351,7 @@ public void setTrackPars(Track cand, double pz = cand.get_P() / Math.sqrt(stateVec.tanThetaX() * stateVec.tanThetaX() + stateVec.tanThetaY() * stateVec.tanThetaY() + 1); - //LOGGER.log(Level.FINE, "Setting track params for ");stateVec.printInfo(); + //LOGGER.log(Level.FINEST, "Setting track params for ");stateVec.printInfo(); dcSwim.SetSwimParameters(stateVec.x(), stateVec.y(), z, pz * stateVec.tanThetaX(), pz * stateVec.tanThetaY(), pz, cand.get_Q()); @@ -489,7 +489,7 @@ public void setTrackPars(Track cand, double pz = cand.get_P() / Math.sqrt(stateVec.tanThetaX() * stateVec.tanThetaX() + stateVec.tanThetaY() * stateVec.tanThetaY() + 1); - //LOGGER.log(Level.FINE, "Setting track params for ");stateVec.printInfo(); + //LOGGER.log(Level.FINEST, "Setting track params for ");stateVec.printInfo(); dcSwim.SetSwimParameters(stateVec.x(), stateVec.y(), z, pz * stateVec.tanThetaX(), pz * stateVec.tanThetaY(), pz, cand.get_Q()); @@ -620,13 +620,13 @@ private Integer getKey(Track trk) { public void removeOverlappingTracksOld(List trkcands) { if(Constants.DEBUG) { - LOGGER.log(Level.FINE, "Found "+trkcands.size()+" HB seeds "); + LOGGER.log(Level.FINEST, "Found "+trkcands.size()+" HB seeds "); for(int i = 0; i< trkcands.size(); i++) { - LOGGER.log(Level.FINE, "cand "+i); + LOGGER.log(Level.FINEST, "cand "+i); for(Cross c : trkcands.get(i)) { - LOGGER.log(Level.FINE, c.printInfo()); + LOGGER.log(Level.FINEST, c.printInfo()); } - LOGGER.log(Level.FINE, "------------------------------------------------------------------ "); + LOGGER.log(Level.FINEST, "------------------------------------------------------------------ "); } } Map selectedTracksMap = new HashMap<>(); @@ -648,26 +648,26 @@ public void removeOverlappingTracksOld(List trkcands) { trkcands.add(entry.getValue()); }); if(Constants.DEBUG) { - LOGGER.log(Level.FINE, "After Overlap Remvr "+trkcands.size()+" HB seeds "); + LOGGER.log(Level.FINEST, "After Overlap Remvr "+trkcands.size()+" HB seeds "); for(int i = 0; i< trkcands.size(); i++) { - LOGGER.log(Level.FINE, "cand "+i); + LOGGER.log(Level.FINEST, "cand "+i); for(Cross c : trkcands.get(i)) { - LOGGER.log(Level.FINE, c.printInfo()); + LOGGER.log(Level.FINEST, c.printInfo()); } - LOGGER.log(Level.FINE, "------------------------------------------------------------------ "); + LOGGER.log(Level.FINEST, "------------------------------------------------------------------ "); } } } public void removeOverlappingTracks(List trkcands) { if(Constants.DEBUG) { - LOGGER.log(Level.FINE, "Found "+trkcands.size()+" HB seeds "); + LOGGER.log(Level.FINEST, "Found "+trkcands.size()+" HB seeds "); for(int i = 0; i< trkcands.size(); i++) { - LOGGER.log(Level.FINE, "cand "+i); + LOGGER.log(Level.FINEST, "cand "+i); for(Cross c : trkcands.get(i)) { - LOGGER.log(Level.FINE, c.printInfo()); + LOGGER.log(Level.FINEST, c.printInfo()); } - LOGGER.log(Level.FINE, "------------------------------------------------------------------ "); + LOGGER.log(Level.FINEST, "------------------------------------------------------------------ "); } } List badTracks = new ArrayList<>(); @@ -680,7 +680,7 @@ public void removeOverlappingTracks(List trkcands) { Track t1 = trkcands.get(i); for(int j=0; jt2.get_FitChi2()/t2.get_FitNDF()) @@ -688,7 +688,7 @@ public void removeOverlappingTracks(List trkcands) { else if(t1.get_FitChi2()/t1.get_FitNDF()==t2.get_FitChi2()/t2.get_FitNDF() && i>j) overlap=true; } -// LOGGER.log(Level.FINE, overlap); +// LOGGER.log(Level.FINEST, overlap); } if(!overlap) selectedTracks.add(t1); } @@ -852,7 +852,7 @@ private double calcCurvSign(Track cand) { private List findStraightTracks(CrossList crossList, DCGeant4Factory DcDetector, double TORSCALE, Swim dcSwim) { - if(LOGGER.getLevel()==Level.FINE) { + if(LOGGER.getLevel()==Level.FINEST) { startTime2 = System.currentTimeMillis(); } @@ -866,11 +866,11 @@ private List findStraightTracks(CrossList crossList, DCGeant4Factory DcDe Track cand = new Track(); TrajectoryFinder trjFind = new TrajectoryFinder(); - if(LOGGER.getLevel()==Level.FINE) { + if(LOGGER.getLevel()==Level.FINEST) { startTime = System.currentTimeMillis(); } Trajectory traj = trjFind.findTrajectory(aCrossList, DcDetector, dcSwim); - LOGGER.log(Level.FINE, "Trajectory finding = " + (System.currentTimeMillis() - startTime)); + LOGGER.log(Level.FINEST, "Trajectory finding = " + (System.currentTimeMillis() - startTime)); if (traj == null) { @@ -892,7 +892,7 @@ private List findStraightTracks(CrossList crossList, DCGeant4Factory DcDe cand.get(0).get_Dir().y() / cand.get(0).get_Dir().z()); cand.set_StateVecAtReg1MiddlePlane(VecAtReg1MiddlePlane); - LOGGER.log(Level.FINE, "Kalman fitter - 2 = " + (System.currentTimeMillis() - startTime)); + LOGGER.log(Level.FINEST, "Kalman fitter - 2 = " + (System.currentTimeMillis() - startTime)); KFitterStraight kFZRef = new KFitterStraight(true, 1, 1, dcSwim, Constants.getInstance().Z, Libr.JNP); List measSurfaces = getMeasSurfaces(cand, DcDetector); @@ -946,7 +946,7 @@ private List findCurvedTracks(CrossList crossList, DCGeant4Factory DcDete boolean donotapplyCuts) { - if(LOGGER.getLevel()==Level.FINE) { + if(LOGGER.getLevel()==Level.FINEST) { startTime2 = System.currentTimeMillis(); } @@ -964,12 +964,12 @@ private List findCurvedTracks(CrossList crossList, DCGeant4Factory DcDete Track cand = new Track(); TrajectoryFinder trjFind = new TrajectoryFinder(); - if(LOGGER.getLevel()==Level.FINE) { + if(LOGGER.getLevel()==Level.FINEST) { startTime = System.currentTimeMillis(); } Trajectory traj = trjFind.findTrajectory(aCrossList, DcDetector, dcSwim); - LOGGER.log(Level.FINE, "Trajectory finding = " + (System.currentTimeMillis() - startTime)); + LOGGER.log(Level.FINEST, "Trajectory finding = " + (System.currentTimeMillis() - startTime)); if (traj == null) { @@ -987,11 +987,11 @@ private List findCurvedTracks(CrossList crossList, DCGeant4Factory DcDete //require 3 crosses to make a track (allows for 1 pseudo-cross) if (cand.size() == 3) { - // LOGGER.log(Level.FINE, "---- cand in sector " + aCrossList.get(0).getSector()); - // LOGGER.log(Level.FINE, aCrossList.get(0).printInfo()); - // LOGGER.log(Level.FINE, aCrossList.get(1).printInfo()); - // LOGGER.log(Level.FINE, aCrossList.get(2).printInfo()); - // LOGGER.log(Level.FINE, "---------------"); + // LOGGER.log(Level.FINEST, "---- cand in sector " + aCrossList.get(0).getSector()); + // LOGGER.log(Level.FINEST, aCrossList.get(0).printInfo()); + // LOGGER.log(Level.FINEST, aCrossList.get(1).printInfo()); + // LOGGER.log(Level.FINEST, aCrossList.get(2).printInfo()); + // LOGGER.log(Level.FINEST, "---------------"); double x1 = aCrossList.get(0).get_Point().x(); double y1 = aCrossList.get(0).get_Point().y(); double z1 = aCrossList.get(0).get_Point().z(); @@ -1032,12 +1032,12 @@ private List findCurvedTracks(CrossList crossList, DCGeant4Factory DcDete if (iBdl != 0) { // momentum estimate if Bdl is non zero and the track has curvature double p = calcInitTrkP(thX, thY, theta1, theta3, iBdl); - if(LOGGER.getLevel()==Level.FINE) { + if(LOGGER.getLevel()==Level.FINEST) { startTime = System.currentTimeMillis(); } int q = this.calcInitTrkQ(traj.getA(), TORSCALE); - LOGGER.log(Level.FINE, "calcInitTrkQ = " + (System.currentTimeMillis() - startTime)); + LOGGER.log(Level.FINEST, "calcInitTrkQ = " + (System.currentTimeMillis() - startTime)); if (p > 11) { p = 11; @@ -1064,7 +1064,7 @@ private List findCurvedTracks(CrossList crossList, DCGeant4Factory DcDete crossIdxinList = 0; } - LOGGER.log(Level.FINE, "Kalman fitter - 2 = " + (System.currentTimeMillis() - startTime)); + LOGGER.log(Level.FINEST, "Kalman fitter - 2 = " + (System.currentTimeMillis() - startTime)); KFitter kFZRef = new KFitter(true, 10, 1, dcSwim, Constants.getInstance().Z, Libr.JNP); List measSurfaces = getMeasSurfaces(cand, DcDetector); @@ -1236,7 +1236,7 @@ public List getMeasSurfaces(Track trkcand, DCGeant4Factory DcDetector) hot._doca[0]*=-LR; hot._hitError = trkcand.get(c).get(s).get(h).get_DocaErr()*trkcand.get(c).get(s).get(h).get_DocaErr(); - //LOGGER.log(Level.FINE, " Z "+Z+" ferr "+(float)(hot._Unc /(hot._hitError/4.))); + //LOGGER.log(Level.FINEST, " Z "+Z+" ferr "+(float)(hot._Unc /(hot._hitError/4.))); hot._Unc[0] = hot._hitError; hOTS.add(hot); diff --git a/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/fit/KFitterDoca.java b/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/fit/KFitterDoca.java index 11d79bde7c..dc93e2e965 100644 --- a/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/fit/KFitterDoca.java +++ b/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/fit/KFitterDoca.java @@ -117,7 +117,7 @@ public void runFitter(int sector) { sv.trackCov.get(svzLength- 1)); */ for (int k = svzLength - 1; k >0; k--) { //if(i==2 && this.totNumIter==30) - //LOGGER.log(Level.FINE, "sector " +sector+"stateVec "+sv.trackTraj.get(k).printInfo()); + //LOGGER.log(Level.FINEST, "sector " +sector+"stateVec "+sv.trackTraj.get(k).printInfo()); if(k>=2) { sv.transport(sector, k, k - 2, sv.trackTraj.get(k), @@ -137,7 +137,7 @@ public void runFitter(int sector) { } for (int k = 0; k < svzLength - 1; k++) { //if(i==2 && this.totNumIter==30) - //LOGGER.log(Level.FINE, "stateVec "+sv.trackTraj.get(k).printInfo()); + //LOGGER.log(Level.FINEST, "stateVec "+sv.trackTraj.get(k).printInfo()); sv.transport(sector, k, k + 1, sv.trackTraj.get(k), sv.trackCov.get(k)); @@ -208,18 +208,18 @@ public Matrix filterCovMat(double[] H, Matrix Ci, double V) { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0); - //LOGGER.log(Level.FINE, "Ci "); + //LOGGER.log(Level.FINEST, "Ci "); //Matrix5x5.show(Ci); - //LOGGER.log(Level.FINE, "Cinv "); + //LOGGER.log(Level.FINEST, "Cinv "); //Matrix5x5.show(first_inverse); - //LOGGER.log(Level.FINE, "addition "); + //LOGGER.log(Level.FINEST, "addition "); //Matrix5x5.show(addition); Matrix5x5.add(first_inverse, addition, result); double det2 = Matrix5x5.inverse(result, result_inv, adj); - //LOGGER.log(Level.FINE, "addition result"); + //LOGGER.log(Level.FINEST, "addition result"); //Matrix5x5.show(result); - //LOGGER.log(Level.FINE, "inv result"); + //LOGGER.log(Level.FINEST, "inv result"); //Matrix5x5.show(result_inv); if(Math.abs(det2)<1.e-30) return null; @@ -276,7 +276,7 @@ private void filter(int k) { // signMeas = Math.signum(h); double c2 = ((signMeas*Math.abs(mv.measurements.get(k).doca[0]) - sign*Math.abs(h)) * (signMeas*Math.abs(mv.measurements.get(k).doca[0]) - sign*Math.abs(h)) / V); - //if(signMeas!=Math.signum(h) && this.interNum>1) LOGGER.log(Level.FINE, sv.trackTraj.get(k).printInfo()+" h "+(float)h); + //if(signMeas!=Math.signum(h) && this.interNum>1) LOGGER.log(Level.FINEST, sv.trackTraj.get(k).printInfo()+" h "+(float)h); double x_filt = sv.trackTraj.get(k).x + K[0] * (signMeas*Math.abs(mv.measurements.get(k).doca[0]) - sign*Math.abs(h)); double y_filt = sv.trackTraj.get(k).y + K[1] * (signMeas*Math.abs(mv.measurements.get(k).doca[0]) - sign*Math.abs(h)); double tx_filt = sv.trackTraj.get(k).tx + K[2] * (signMeas*Math.abs(mv.measurements.get(k).doca[0]) - sign*Math.abs(h)); diff --git a/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/fit/MeasVecsDoca.java b/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/fit/MeasVecsDoca.java index 6dba453a00..59e3388e1c 100644 --- a/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/fit/MeasVecsDoca.java +++ b/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/fit/MeasVecsDoca.java @@ -46,7 +46,7 @@ public double h(double[] stateV, double Z, Line3D wireLine) { WL.copy(wireLine); WL.copy(WL.distance(new Point3D(stateV[0], stateV[1], Z))); - //LOGGER.log(Level.FINE, Math.signum(-WL.direction().x())+ + //LOGGER.log(Level.FINEST, Math.signum(-WL.direction().x())+ // wireLine.origin().toString()+WL.toString()+" "+stateV[0]+" ,"+stateV[1]); return WL.length()*Math.signum(WL.direction().x()); } @@ -100,7 +100,7 @@ public void setMeasVecs(Track trkcand, DCGeant4Factory DcDetector) { hot._doca[0]*=LR; //hot._hitError = err_sl1 * err_sl1 * Z * Z + err_it1 * err_it1 + 2 * Z * err_cov1 + trk.get_ListOfHBSegments().get(s).get(h).get_DocaErr()*trk.get_ListOfHBSegments().get(s).get(h).get_DocaErr(); hot._hitError = trkcand.get(c).get(s).get(h).get_DocaErr()*trkcand.get(c).get(s).get(h).get_DocaErr(); - //LOGGER.log(Level.FINE, " Z "+Z+" ferr "+(float)(hot._Unc /(hot._hitError/4.))); + //LOGGER.log(Level.FINEST, " Z "+Z+" ferr "+(float)(hot._Unc /(hot._hitError/4.))); hot._Unc[0] = hot._hitError; hOTS.add(hot); @@ -162,7 +162,7 @@ void setMeasVecsFromHB(Track trk, DCGeant4Factory DcDetector) { double LR = Math.signum(trk.get_ListOfHBSegments().get(s).get(h).get_XWire()-trk.get_ListOfHBSegments().get(s).get(h).get_X()); hot._doca[0]*=LR; hot._hitError = trk.get_ListOfHBSegments().get(s).get(h).get_DocaErr()*trk.get_ListOfHBSegments().get(s).get(h).get_DocaErr(); - //LOGGER.log(Level.FINE, " Z "+Z+" ferr "+(float)(hot._Unc /(hot._hitError/4.))); + //LOGGER.log(Level.FINEST, " Z "+Z+" ferr "+(float)(hot._Unc /(hot._hitError/4.))); hot._Unc[0] = hot._hitError; hot.region = trk.get_ListOfHBSegments().get(s).get(h).get_Region(); hOTS.add(hot); @@ -199,7 +199,7 @@ void setMeasVecsFromHB(Track trk, DCGeant4Factory DcDetector) { meas.wireLine[0] = hOTS.get(i)._wireLine[0]; meas.wireLine[1] = hOTS.get(i)._wireLine[1]; this.measurements.add(i, meas); - //LOGGER.log(Level.FINE, " measurement "+i+" = "+meas.x+" at "+meas.z); + //LOGGER.log(Level.FINEST, " measurement "+i+" = "+meas.x+" at "+meas.z); } } diff --git a/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/fit/StateVecsDoca.java b/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/fit/StateVecsDoca.java index 25fe85a220..70d56dd75c 100644 --- a/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/fit/StateVecsDoca.java +++ b/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/fit/StateVecsDoca.java @@ -78,7 +78,7 @@ public Matrix transport(int sector, int i, double Zf, StateVec iVec, CovMat covM double BatMeas = iVec.B; while(Math.signum(Zf - Z[i]) *zMath.signum(Zf - Z[i]) *Zf) s=Math.signum(Zf - Z[i]) *Math.abs(Zf-z); rk.RK4transport(sector, s, dcSwim, covMat, fVec, fCov); @@ -169,7 +169,7 @@ public void transport(int sector, int i, int f, StateVec iVec, CovMat covMat) { double BatMeas = iVec.B; while(Math.signum(Z[f] - Z[i]) *zMath.signum(Z[f] - Z[i]) *Z[f]) s=Math.signum(Z[f] - Z[i]) *Math.abs(Z[f]-z); @@ -277,7 +277,7 @@ public void transportFixed(int sector, int i, int f, StateVec iVec, CovMat covMa if (j == nSteps - 1) { s = Math.signum(Z[f] - Z[i]) * Math.abs(z - Z[f]); } - //LOGGER.log(Level.FINE, " RK step num "+(j+1)+" = "+(float)s+" nSteps = "+nSteps); + //LOGGER.log(Level.FINEST, " RK step num "+(j+1)+" = "+(float)s+" nSteps = "+nSteps); double x = fVec.x; double y = fVec.y; z = fVec.z; @@ -368,7 +368,7 @@ public void init(Track trkcand, double z0, KFitterDoca kf, int c) { kf.setFitFailed = true; return; } - //LOGGER.log(Level.FINE, (0)+"] init "+this.trackTraj.get(0).x+","+this.trackTraj.get(0).y+","+ + //LOGGER.log(Level.FINEST, (0)+"] init "+this.trackTraj.get(0).x+","+this.trackTraj.get(0).y+","+ // this.trackTraj.get(0).z+","+this.trackTraj.get(0).tx+","+this.trackTraj.get(0).ty+" "+1/this.trackTraj.get(0).Q); double err_sl1 = trkcand.get(0).get_Segment1().get_fittedCluster().get_clusterLineFitSlopeErr(); double err_sl2 = trkcand.get(0).get_Segment2().get_fittedCluster().get_clusterLineFitSlopeErr(); diff --git a/reconstruction/dc/src/main/java/org/jlab/service/dc/DCEngine.java b/reconstruction/dc/src/main/java/org/jlab/service/dc/DCEngine.java index cdc44a3724..5bfba7b14e 100644 --- a/reconstruction/dc/src/main/java/org/jlab/service/dc/DCEngine.java +++ b/reconstruction/dc/src/main/java/org/jlab/service/dc/DCEngine.java @@ -222,7 +222,7 @@ public int getRun(DataEvent event) { return 0; } DataBank bank = event.getBank("RUN::config"); - LOGGER.log(Level.FINE,"["+this.getName()+"] EVENT "+bank.getInt("event", 0)); + LOGGER.log(Level.FINEST,"["+this.getName()+"] EVENT "+bank.getInt("event", 0)); int run = bank.getInt("run", 0); return run; diff --git a/reconstruction/dc/src/main/java/org/jlab/service/dc/DCHBPostClusterAI.java b/reconstruction/dc/src/main/java/org/jlab/service/dc/DCHBPostClusterAI.java index 7ac9ed6969..db268e567e 100644 --- a/reconstruction/dc/src/main/java/org/jlab/service/dc/DCHBPostClusterAI.java +++ b/reconstruction/dc/src/main/java/org/jlab/service/dc/DCHBPostClusterAI.java @@ -66,7 +66,7 @@ public boolean processDataEvent(DataEvent event) { RecoBankWriter writer = new RecoBankWriter(this.getBanks()); // get Field Swim dcSwim = new Swim(); - LOGGER.log(Level.FINE, "HB AI process event"); + LOGGER.log(Level.FINEST, "HB AI process event"); //AI List trkcands = null; @@ -96,11 +96,11 @@ public boolean processDataEvent(DataEvent event) { CrossList crosslist = pr.RecomposeCrossList(segments, Constants.getInstance().dcDetector); crosses = new ArrayList<>(); - LOGGER.log(Level.FINE, "num cands = "+crosslist.size()); + LOGGER.log(Level.FINEST, "num cands = "+crosslist.size()); for (List clist : crosslist) { crosses.addAll(clist); for(Cross c : clist) - LOGGER.log(Level.FINE, "Pass Cross"+c.printInfo()); + LOGGER.log(Level.FINEST, "Pass Cross"+c.printInfo()); } if (crosses.isEmpty()) { for(Segment seg : segments) { diff --git a/reconstruction/dc/src/main/java/org/jlab/service/dc/DCHBPostClusterConv.java b/reconstruction/dc/src/main/java/org/jlab/service/dc/DCHBPostClusterConv.java index 2cf8c6d5b2..1c8a2bcbdf 100644 --- a/reconstruction/dc/src/main/java/org/jlab/service/dc/DCHBPostClusterConv.java +++ b/reconstruction/dc/src/main/java/org/jlab/service/dc/DCHBPostClusterConv.java @@ -228,13 +228,13 @@ public boolean processDataEvent(DataEvent event) { trkcands.addAll(mistrkcands); - LOGGER.log(Level.FINE, "Found after 5STg "+mistrkcands.size()+" HB seeds "); + LOGGER.log(Level.FINEST, "Found after 5STg "+mistrkcands.size()+" HB seeds "); for(int i = 0; i< trkcands.size(); i++) { - LOGGER.log(Level.FINE, "cand "+i); + LOGGER.log(Level.FINEST, "cand "+i); for(Cross c : trkcands.get(i)) { - LOGGER.log(Level.FINE, c.printInfo()); + LOGGER.log(Level.FINEST, c.printInfo()); } - LOGGER.log(Level.FINE, "------------------------------------------------------------------ "); + LOGGER.log(Level.FINEST, "------------------------------------------------------------------ "); } //gather all the hits for pointer bank creation diff --git a/reconstruction/dc/src/main/java/org/jlab/service/dc/DCTBEngine.java b/reconstruction/dc/src/main/java/org/jlab/service/dc/DCTBEngine.java index 096a1eee80..93abb7cf23 100644 --- a/reconstruction/dc/src/main/java/org/jlab/service/dc/DCTBEngine.java +++ b/reconstruction/dc/src/main/java/org/jlab/service/dc/DCTBEngine.java @@ -106,7 +106,7 @@ public boolean processDataEvent(DataEvent event) { List crosses = new ArrayList<>(); List trkcands = new ArrayList<>(); - LOGGER.log(Level.FINE, "TB AI "+ this.getName()); + LOGGER.log(Level.FINEST, "TB AI "+ this.getName()); //instantiate bank writer RecoBankWriter rbc = new RecoBankWriter(this.getBanks()); @@ -331,7 +331,7 @@ public boolean processDataEvent(DataEvent event) { //trk.set_Id(trkId); trkcandFinder.matchHits(trk.getStateVecs(), trk, Constants.getInstance().dcDetector, dcSwim); trk.calcTrajectory(trkId, dcSwim, trk.get_Vtx0(), trk.get_pAtOrig(), trk.get_Q()); - LOGGER.log(Level.FINE, trk.toString()); + LOGGER.log(Level.FINEST, trk.toString()); for(Cross c : trk) { c.set_CrossDirIntersSegWires(); @@ -574,7 +574,7 @@ public List getMeasSurfaces(Track trk, DCGeant4Factory DcDetector) { double LR = Math.signum(trk.get_ListOfHBSegments().get(s).get(h).get_XWire()-trk.get_ListOfHBSegments().get(s).get(h).get_X()); hot._doca[0]*=-LR; hot._hitError = trk.get_ListOfHBSegments().get(s).get(h).get_DocaErr()*trk.get_ListOfHBSegments().get(s).get(h).get_DocaErr(); - //LOGGER.log(Level.FINE, " Z "+Z+" ferr "+(float)(hot._Unc /(hot._hitError/4.))); + //LOGGER.log(Level.FINEST, " Z "+Z+" ferr "+(float)(hot._Unc /(hot._hitError/4.))); hot._Unc[0] = hot._hitError; hot.region = trk.get_ListOfHBSegments().get(s).get(h).get_Region(); hot.sector = trk.get_ListOfHBSegments().get(s).get(h).get_Sector(); @@ -671,11 +671,11 @@ private int get_Status(Track track) { miss=l+1; if(miss%2==0 && SegMap.containsKey(l)) { //missing sl in 2,4,6 track.setSingleSuperlayer(SegMap.get(l)); //isolated sl in 1,3,5 - LOGGER.log(Level.FINE, "Missing superlayer "+miss+" seg "+SegMap.get(l).printInfo()); + LOGGER.log(Level.FINEST, "Missing superlayer "+miss+" seg "+SegMap.get(l).printInfo()); } else if(miss%2==1 && SegMap.containsKey(l+2)) { //missing sl in 1,3,5 track.setSingleSuperlayer(SegMap.get(l+2)); //isolated sl in 2,4,6 - LOGGER.log(Level.FINE, "Missing superlayer "+miss+" seg "+track.getSingleSuperlayer().printInfo()); + LOGGER.log(Level.FINEST, "Missing superlayer "+miss+" seg "+track.getSingleSuperlayer().printInfo()); } } } From 601caf728d42e9ea1059f720bd348bcd8abd5d6f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mariangela=20Bond=C3=AD?= Date: Thu, 3 Jul 2025 12:45:03 +0200 Subject: [PATCH 08/16] new FVT detector for ddvcs proposal; Improved strip geometry for URWell detector --- .../geant4/v2/URWELL/URWellConstants.java | 34 +- .../geant4/v2/URWELL/URWellGeant4Factory.java | 218 ++++----- .../geant4/v2/URWELL/URWellStripFactory.java | 450 +++++++++++------- 3 files changed, 418 insertions(+), 284 deletions(-) diff --git a/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellConstants.java b/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellConstants.java index c41492485f..c9ef54163d 100644 --- a/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellConstants.java +++ b/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellConstants.java @@ -10,26 +10,28 @@ public class URWellConstants { private final static String CCDBPATH = "/geometry/urwell/"; public final static int NMAXREGIONS = 6; //max number of regions - public final static int NREGIONS = 1; //number of regions + public final static int NREGIONS = 6; //number of regions public final static int NSECTORS = 6; //number of sectors public final static int NLAYERS = 2; //number of layers public final static int NCHAMBERS = 3; //number of chambers in a sector public final static int NCHAMBERS_ddvcs = 1; //number of chambers in a sector - public final static double XENLARGEMENT = 0.5; // cm - public final static double YENLARGEMENT = 1.; // cm + public final static double XENLARGEMENT = 0.1; // cm + public final static double YENLARGEMENT = 0.1; // cm public final static double ZENLARGEMENT = 0.1; // cm // Sector geometrical parameters public final static double THOPEN = 54.; // opening angle between endplate planes (deg) public final static double THTILT = 25; // theta tilt (deg) public final static double THMIN = 4.694; // polar angle to the base of first chamber (deg) + public final static double THMAX = 40; public final static double SECTORHEIGHT = 146.21; //height of each sector (cm) public final static double DX0CHAMBER0 = 5.197; // halfbase of chamber 1 (cm) public final static double SECTORHEIGHT_ddvcs = 22; - public final static double THMIN_ddvcs = 7; - public final static double DX0CHAMBER0_ddvcs = 3.597; + public final static double THMIN_ddvcs = 6.8; + public final static double DX0CHAMBER0_ddvcs = 3.207; + public final static double THMAX_ddvcs = 35; // Chamber volumes and materials (units are cm) public final static double[] CHAMBERVOLUMESTHICKNESS = {0.0025, 0.0005,0.3, // window @@ -55,9 +57,10 @@ public class URWellConstants { public final static double W2TGT[] = new double[NMAXREGIONS];; public final static double YMIN[] = new double[NMAXREGIONS]; public final static double ZMIN[] = new double[NMAXREGIONS]; + public final static double HSECTOR[] = new double[NMAXREGIONS]; - public final static double TGT2DC0_ddvcs = 52.9; // cm + public final static double TGT2DC0_ddvcs = 48.9; // cm // public final static double URWELL2DC0 = 2; // cm public final static double URWELL2DC0_ddvcs[] = new double[NMAXREGIONS]; public final static double DIST2TGT_ddvcs[] = new double[NMAXREGIONS]; @@ -65,13 +68,15 @@ public class URWellConstants { ; public final static double YMIN_ddvcs[] = new double[NMAXREGIONS]; public final static double ZMIN_ddvcs[] = new double[NMAXREGIONS]; + public final static double HSECTOR_ddvcs[] = new double[NMAXREGIONS]; // public final static double DIST2TGT = (TGT2DC0-URWELL2DC0); // public final static double W2TGT = DIST2TGT/Math.cos(Math.toRadians(THTILT-THMIN)); // public final static double YMIN = W2TGT*Math.sin(Math.toRadians(THMIN)); // distance from the base chamber1 and beamline // public final static double ZMIN = W2TGT*Math.cos(Math.toRadians(THMIN)); public final static double PITCH = 0.1 ; // cm - public final static double STEREOANGLE = 6; // deg + public final static double PITCH_ddvcs = 0.05 ; // cm + public final static double STEREOANGLE = 10; // deg @@ -136,14 +141,21 @@ public static synchronized void load( DatabaseConstantProvider cp ) W2TGT[i] = DIST2TGT[i]/Math.cos(Math.toRadians(THTILT-THMIN)); YMIN[i]= W2TGT[i]*Math.sin(Math.toRadians(THMIN)); // distance from the base chamber1 and beamline ZMIN[i] = W2TGT[i]*Math.cos(Math.toRadians(THMIN)); - - URWELL2DC0_ddvcs[i] = -2. + i * 1.3; + + double h1 = DIST2TGT[i]*Math.tan((Math.toRadians(THMAX -THTILT ))); + double h2 = DIST2TGT[i]*Math.tan((Math.toRadians(THTILT - THMIN))); + HSECTOR[i] = h1 +h2; + + + URWELL2DC0_ddvcs[i] = -2. + i * 2; DIST2TGT_ddvcs[i] = (TGT2DC0_ddvcs + URWELL2DC0_ddvcs[i]); W2TGT_ddvcs[i] = DIST2TGT_ddvcs[i] / Math.cos(Math.toRadians(THTILT - THMIN_ddvcs)); YMIN_ddvcs[i] = W2TGT_ddvcs[i] * Math.sin(Math.toRadians(THMIN_ddvcs)); // distance from the base chamber1 and beamline ZMIN_ddvcs[i] = W2TGT_ddvcs[i] * Math.cos(Math.toRadians(THMIN_ddvcs)); - - + double h1_ddvcs = DIST2TGT_ddvcs[i]*Math.tan((Math.toRadians(THMAX_ddvcs -THTILT ))); + double h2_ddvcs = DIST2TGT_ddvcs[i]*Math.tan((Math.toRadians(THTILT - THMIN_ddvcs))); + HSECTOR_ddvcs[i] = h1_ddvcs +h2_ddvcs; + } diff --git a/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellGeant4Factory.java b/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellGeant4Factory.java index 6d128ded5d..a20cb42604 100644 --- a/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellGeant4Factory.java +++ b/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellGeant4Factory.java @@ -12,6 +12,7 @@ + /** * Generate GEANT4 volume for the URWELL detector * @@ -29,7 +30,8 @@ public final class URWellGeant4Factory extends Geant4Factory { private boolean isProto = false; private double[] ymin = new double[URWellConstants.NMAXREGIONS]; private double[] zmin = new double[URWellConstants.NMAXREGIONS]; - + private double[] Hsector = new double[URWellConstants.NMAXREGIONS]; + private String name = ""; /** * Create the URWELL full geometry @@ -59,44 +61,51 @@ public void init(DatabaseConstantProvider cp, int config, int regions ) { switch (config) { case 0 -> { - + nSectors = URWellConstants.NSECTORS; nChambers = URWellConstants.NCHAMBERS; nRegions = Math.min(URWellConstants.NMAXREGIONS, regions); isProto = false; sectorheight = URWellConstants.SECTORHEIGHT; + Arrays.fill(Hsector, URWellConstants.SECTORHEIGHT); thilt = URWellConstants.THTILT; dx0chamber0 = URWellConstants.DX0CHAMBER0; thopen = URWellConstants.THOPEN; ymin =Arrays.copyOf(URWellConstants.YMIN, URWellConstants.YMIN.length); zmin =Arrays.copyOf(URWellConstants.ZMIN, URWellConstants.ZMIN.length); + name =""; } case 1 -> { + nRegions = URWellConstants.NREGIONS_PROTO; + nSectors = URWellConstants.NSECTORS_PROTO; + nChambers = URWellConstants.NCHAMBERS_PROTO; + isProto = true; + name = "_proto"; + + } + case 2 -> { + nRegions = Math.min(URWellConstants.NMAXREGIONS, regions); nSectors = URWellConstants.NSECTORS; nChambers = URWellConstants.NCHAMBERS_ddvcs; + Hsector = Arrays.copyOf(URWellConstants.HSECTOR_ddvcs, URWellConstants.HSECTOR_ddvcs.length); sectorheight = URWellConstants.SECTORHEIGHT_ddvcs; thilt = URWellConstants.THTILT; dx0chamber0 = URWellConstants.DX0CHAMBER0_ddvcs; thopen = URWellConstants.THOPEN; ymin = Arrays.copyOf(URWellConstants.YMIN_ddvcs, URWellConstants.YMIN_ddvcs.length); zmin = Arrays.copyOf(URWellConstants.ZMIN_ddvcs, URWellConstants.ZMIN_ddvcs.length); - + name = "_ddvcs"; isProto = false; } - case 2 -> { - nRegions = URWellConstants.NREGIONS_PROTO; - nSectors = URWellConstants.NSECTORS_PROTO; - nChambers = URWellConstants.NCHAMBERS_PROTO; - isProto = true; - - } default -> { } } for (int iregion = 0; iregion (volume.getName().contains(volumeName))) - .findAny() - .orElse(null); + + volumeName = "rg" + r + "_s" + s + "_c" + c + "_cathode_gas" + name; + + return this.getAllVolumes().stream() + .filter(volume -> (volume.getName().contains(volumeName))) + .findAny() + .orElse(null); } /** @@ -475,7 +476,8 @@ public Geant4Basic getSectorVolume(int region, int sector) { int r = region; int s = sector; - String volName = "region_uRwell_" + r + "_s" + s; + String volName = "region_uRwell_" + r + "_s" + s +name; + return this.getAllVolumes().stream() .filter(volume -> (volume.getName().contains(volName))) .findAny() @@ -492,7 +494,7 @@ public static void main(String[] args) { URWellConstants.connect(cp); // URWellGeant4Factory factory = new URWellGeant4Factory(cp, false, 2); - URWellGeant4Factory factory = new URWellGeant4Factory(cp, 1, 2); + URWellGeant4Factory factory = new URWellGeant4Factory(cp, 2, 6); factory.getAllVolumes().forEach(volume -> { System.out.println(volume.gemcString()); }); diff --git a/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellStripFactory.java b/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellStripFactory.java index a7ec2edb5f..469c1f302f 100644 --- a/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellStripFactory.java +++ b/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellStripFactory.java @@ -1,6 +1,6 @@ package org.jlab.detector.geant4.v2.URWELL; - +import eu.mihosoft.vrl.v3d.Transform; import eu.mihosoft.vrl.v3d.Vector3d; import java.util.List; import org.jlab.detector.calib.utils.DatabaseConstantProvider; @@ -16,46 +16,51 @@ /** * Creates and handles the URWELL detector strips as 3D lines - * + * * @author bondi */ public final class URWellStripFactory { private URWellGeant4Factory factory; - private IndexedList globalStrips = new IndexedList(3); - private IndexedList localStrips = new IndexedList(3); - private IndexedList planeStrips = new IndexedList(3); + private IndexedList globalStrips = new IndexedList(3); + private IndexedList localStrips = new IndexedList(3); + private IndexedList planeStrips = new IndexedList(3); private int nRegions; private int nSectors; private int nChambers; private int nLayers; + private int Nconfig; + private double pitch; private boolean isProto; - + public URWellStripFactory() { } - + /** - * Create the strip factory based on constants from CCDB. - * Currently constants are defined in the URWellConstants class. - * They will be moved to CCDB when finalized). + * Create the strip factory based on constants from CCDB. Currently + * constants are defined in the URWellConstants class. They will be moved to + * CCDB when finalized). + * * @param cp database provide */ public URWellStripFactory(DatabaseConstantProvider cp) { this.init(cp); } - + /** * Initialize the factory by the strip maps + * * @param cp */ public void init(DatabaseConstantProvider cp) { this.init(cp, false, 1); } - + /** - * Create the strip factory based on constants from CCDB. - * Currently constants are defined in the URWellConstants class. - * They will be moved to CCDB when finalized). + * Create the strip factory based on constants from CCDB. Currently + * constants are defined in the URWellConstants class. They will be moved to + * CCDB when finalized). + * * @param cp database provide * @param prototype * @param regions @@ -63,9 +68,23 @@ public void init(DatabaseConstantProvider cp) { public URWellStripFactory(DatabaseConstantProvider cp, boolean prototype, int regions) { this.init(cp, prototype, regions); } - + + /** + * Create the strip factory based on constants from CCDB. Currently + * constants are defined in the URWellConstants class. They will be moved to + * CCDB when finalized). + * + * @param cp database provide + * @param config 0 standard, 1 ddvcs proposal, 2 prototype + * @param regions + */ + public URWellStripFactory(DatabaseConstantProvider cp, int config, int regions) { + this.init(cp, config, regions); + } + /** * Initialize the factory by the strip maps + * * @param cp * @param prototype * @param regions @@ -73,49 +92,132 @@ public URWellStripFactory(DatabaseConstantProvider cp, boolean prototype, int re public void init(DatabaseConstantProvider cp, boolean prototype, int regions) { factory = new URWellGeant4Factory(cp, prototype, regions); isProto = prototype; - if(!isProto){ - nRegions = Math.min(URWellConstants.NMAXREGIONS, regions); - nSectors = URWellConstants.NSECTORS; + if (!isProto) { + nRegions = Math.min(URWellConstants.NMAXREGIONS, regions); + nSectors = URWellConstants.NSECTORS; nChambers = URWellConstants.NCHAMBERS; - nLayers = URWellConstants.NLAYERS; - } - else { - nRegions = URWellConstants.NREGIONS_PROTO; - nSectors = URWellConstants.NSECTORS_PROTO; + nLayers = URWellConstants.NLAYERS; + pitch = URWellConstants.PITCH; + } else { + nRegions = URWellConstants.NREGIONS_PROTO; + nSectors = URWellConstants.NSECTORS_PROTO; nChambers = URWellConstants.NCHAMBERS_PROTO; - nLayers = URWellConstants.NLAYERS; + nLayers = URWellConstants.NLAYERS; + pitch = URWellConstants.PITCH; } this.fillStripLists(); this.fillPlaneLists(); } + /** + * Initialize the factory by the strip maps + * + * @param cp + * @param config 0 standard, 1 ddvcs proposal, 2 prototype + * @param regions + */ + public void init(DatabaseConstantProvider cp, int config, int regions) { + factory = new URWellGeant4Factory(cp, config, regions); + + switch (config) { + case 0 -> { + + nSectors = URWellConstants.NSECTORS; + nChambers = URWellConstants.NCHAMBERS; + nRegions = Math.min(URWellConstants.NMAXREGIONS, regions); + nLayers = URWellConstants.NLAYERS; + isProto = false; + Nconfig = config; + pitch = URWellConstants.PITCH; + } + case 1 -> { + + nRegions = URWellConstants.NREGIONS_PROTO; + nSectors = URWellConstants.NSECTORS_PROTO; + nChambers = URWellConstants.NCHAMBERS_PROTO; + nLayers = URWellConstants.NLAYERS; + isProto = true; + Nconfig = config; + pitch = URWellConstants.PITCH; + } + case 2 -> { + nRegions = Math.min(URWellConstants.NMAXREGIONS, regions); + nSectors = URWellConstants.NSECTORS; + nChambers = URWellConstants.NCHAMBERS_ddvcs; + nLayers = URWellConstants.NLAYERS; + isProto = false; + Nconfig = config; + pitch = URWellConstants.PITCH_ddvcs; + + } + default -> { + } + } + + this.fillStripLists(); + this.fillPlaneLists(); + } + + public double getStereoAngle(int layer, int ichamber) { + + double stereoAngle = Math.toRadians(URWellConstants.STEREOANGLE); + + double[] dim = factory.getChamber_daughter_Dimensions(layer, ichamber); + + double yHalf = dim[0]; + double xHalfSmallBase = dim[1]; + double xHalfLargeBase = dim[2]; + + if (Nconfig == 2) { + + stereoAngle = Math.abs(Math.atan2(2 * yHalf, (-xHalfSmallBase + xHalfLargeBase))); + + } + + + if (layer % 2 != 0) { + stereoAngle = -stereoAngle; + } + + return stereoAngle; + + } + /** * Calculates the total number of strips in a sector - * + * + * @param layer * @return the strip number */ - public int getNStripSector() { + public int getNStripSector(int layer) { int nStrips = 0; for (int i = 0; i < nChambers; i++) { - nStrips += getNStripChamber(i); + // System.out.println("indice camera "+i+" numero strip "+getNStripChamber(i)); + + nStrips += getNStripChamber(layer, i); + } return nStrips; } /** * Calculates the number of strips in the given chamber - * + * + * @param layer (1 to N) * @param ichamber (0, 1, 2) * @return the strip number (1-N) */ - public int getNStripChamber(int ichamber) { - - double[] dim = factory.getChamber_daughter_Dimensions(ichamber); + public int getNStripChamber(int layer, int ichamber) { + + + double[] dim = factory.getChamber_daughter_Dimensions(layer, ichamber); - double yHalf = dim[0]; + double yHalf = dim[0]; double xHalfSmallBase = dim[1]; double xHalfLargeBase = dim[2]; + double stereoAngle = Math.abs(getStereoAngle(layer, ichamber)); + // C-------------D // // ------------- // // ----------- // @@ -123,59 +225,66 @@ public int getNStripChamber(int ichamber) { /** * * number of strip in AB** */ - int nAB = (int) (2 * xHalfSmallBase / (URWellConstants.PITCH - / Math.sin(Math.toRadians(URWellConstants.STEREOANGLE)))); + + int nAB = (int) Math.round(2 * xHalfSmallBase / (pitch / Math.sin(stereoAngle))); + + + double AC = Math.sqrt((Math.pow((xHalfSmallBase - xHalfLargeBase), 2) + Math.pow((2 * yHalf), 2))); double theta = Math.acos(2 * yHalf / AC); - int nAC = (int) (AC / (URWellConstants.PITCH - / Math.cos(theta - Math.toRadians(URWellConstants.STEREOANGLE)))); - int nStrips = nAB + nAC +1 ; + int nAC = (int) Math.round(AC / (pitch / Math.cos(theta - stereoAngle))); + + int nStrips = nAB + nAC; return nStrips; } /** * Provides the index of the chamber containing the strip with the given ID - * + * + * @param layer (1 to N) * @param strip (1 to N) * @return the chamber index (0, 1, 2) */ - public int getChamberIndex(int strip) { + public int getChamberIndex(int layer, int strip) { int nStripTotal = 0; - for(int i=0; i strip ID chamber (from 1 to getNStripChamber) int nStripTotal = 0; - if (chamberIndex > 0) { + for (int i = 0; i < chamberIndex; i++) { - nStripTotal += this.getNStripChamber(i); + + nStripTotal += this.getNStripChamber(layer, i); } - } - + + //Strip ID: from 1 to getNStripChamber int cStrip = strip - nStripTotal; - + return cStrip; } - + /** * Builds the given strip line in the CLAS12 frame + * * @param sector (1-6) * @param layer (1-2) * @param strip (1-N) @@ -183,67 +292,60 @@ private int getLocalStripId(int strip) { */ private Line3d createStrip(int sector, int layer, int strip) { - int chamberIndex = getChamberIndex(strip); - - int cStrip = this.getLocalStripId(strip); + int chamberIndex = getChamberIndex(layer, strip); + double stereoAngle = getStereoAngle(layer,chamberIndex); - + int cStrip = this.getLocalStripId(layer, strip); // CHAMBER reference frame // new numeration with stri ID_strip=0 crossing (0,0,0) of chamber - double[] dim = factory.getChamber_daughter_Dimensions(chamberIndex); - - double yHalf = dim[0]; + double[] dim = factory.getChamber_daughter_Dimensions(layer, chamberIndex); + + double yHalf = dim[0]; double xHalfSmallBase = dim[1]; double xHalfLargeBase = dim[2]; - - + // Y coordinate of the intersection point between the x=0 and the strip line crossing for B + double DY = -yHalf - Math.tan(Math.abs(stereoAngle)) * xHalfSmallBase; - double DY = -yHalf - Math.tan(Math.toRadians(URWellConstants.STEREOANGLE)) *xHalfSmallBase; - // ID of the strip - int nS = (int) (DY * Math.cos(Math.toRadians(URWellConstants.STEREOANGLE)) / URWellConstants.PITCH); + int nS = (int) (DY * Math.cos(Math.abs(stereoAngle)) / pitch); int nCStrip = nS + (cStrip - 1); - - //strip straight line chamber reference frame -> y = mx +c; - double stereoAngle = URWellConstants.STEREOANGLE; - if (layer % 2 != 0) { - stereoAngle = -URWellConstants.STEREOANGLE; - } - double m = Math.tan(Math.toRadians(stereoAngle)); - double c = nCStrip * URWellConstants.PITCH / Math.cos(Math.toRadians(stereoAngle)); - + + + + double m = Math.tan(stereoAngle); + double c = nCStrip * pitch / Math.cos(stereoAngle); + // Take 2 points in the strip straight line. They needs to define Line object double oX = -xHalfLargeBase; double oY = -xHalfLargeBase * m + c; double oZ = 0; Vector3d origin = new Vector3d(oX, oY, oZ); - + double eX = xHalfLargeBase; double eY = xHalfLargeBase * m + c; double eZ = 0; Vector3d end = new Vector3d(eX, eY, eZ); - + // Get Chamber Volume - Geant4Basic chamberVolume = factory.getChamberVolume(sector, chamberIndex+1, layer, isProto); - + Geant4Basic chamberVolume = factory.getChamberVolume(sector, chamberIndex + 1, layer); + // 2 point defined before wrt the GLOBAL frame Vector3d globalOrigin = chamberVolume.getGlobalTransform().transform(origin); - - Vector3d globalEnd = chamberVolume.getGlobalTransform().transform(end); + Vector3d globalEnd = chamberVolume.getGlobalTransform().transform(end); Straight line = new Line3d(globalOrigin, globalEnd); - + // CHECK intersections between line and volume chamberVolume.makeSensitive(); List Hits = chamberVolume.getIntersections(line); - + if (Hits.size() >= 1) { - - Vector3d TestOrigin = Hits.get(0).origin(); - Vector3d TestEnd = Hits.get(0).end(); + + Vector3d TestOrigin = Hits.get(0).origin(); + Vector3d TestEnd = Hits.get(0).end(); return new Line3d(Hits.get(0).origin(), Hits.get(0).end()); @@ -252,75 +354,80 @@ private Line3d createStrip(int sector, int layer, int strip) { } } - /** + /** * Provides the given strip line in the Chamber local frame - * @param region (1-2) + * + * * @param sector (1-6) * @param layer (1-4) * @param strip (1-N) * @return the 3D strip line as a Line3d */ - - private Line3d getChamberStrip(int region, int sector, int chamber, int layer, int strip) { + private Line3d getChamberStrip(int sector, int layer, int strip) { - Line3d globalStrip = createStrip(sector, layer, strip); - Geant4Basic chamberVolume = factory.getChamberVolume(sector, chamber, layer, isProto); + int chamber = getChamberIndex(layer, strip); + + Geant4Basic chamberVolume = factory.getChamberVolume(sector, chamber+1, layer); Vector3d origin = chamberVolume.getGlobalTransform().invert().transform(globalStrip.origin()); - Vector3d end = chamberVolume.getGlobalTransform().invert().transform(globalStrip.end()); + Vector3d end = chamberVolume.getGlobalTransform().invert().transform(globalStrip.end()); Line3d localStrip = new Line3d(origin, end); - return localStrip; } - - /** * Provides the given strip line in the sector local frame + * * @param sector (1-6) * @param layer (1-2) * @param strip (1-N) * @return the 3D strip line as a Line3d */ - private Line3d getLocalStrip(int region, int sector, int layer, int strip) { + private Line3d getLocalStrip(int sector, int layer, int strip) { + + int region = (layer - 1) / 2 + 1; Line3d globalStrip = createStrip(sector, layer, strip); Geant4Basic sVolume = factory.getSectorVolume(region, sector); Vector3d origin = sVolume.getGlobalTransform().invert().transform(globalStrip.origin()); - Vector3d end = sVolume.getGlobalTransform().invert().transform(globalStrip.end()); + Vector3d end = sVolume.getGlobalTransform().invert().transform(globalStrip.end()); Line3d localStrip = new Line3d(origin, end); return localStrip; } - - + private void fillStripLists() { - - for(int ir=0; ir Date: Thu, 3 Jul 2025 13:13:34 +0200 Subject: [PATCH 09/16] fix a bug on URWellStrip in the recostruction part --- .../org/jlab/detector/geant4/v2/URWELL/URWellStripFactory.java | 2 +- .../src/main/java/org/jlab/service/urwell/URWellStrip.java | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellStripFactory.java b/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellStripFactory.java index 469c1f302f..f78c8c40b2 100644 --- a/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellStripFactory.java +++ b/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellStripFactory.java @@ -551,7 +551,7 @@ public static void main(String[] args) { double angle = factory2.getStereoAngle(0,0); int strip = 1418; - int layer = 2; + int layer = 1; int sector =2; System.out.println("stereo angle"); diff --git a/reconstruction/urwell/src/main/java/org/jlab/service/urwell/URWellStrip.java b/reconstruction/urwell/src/main/java/org/jlab/service/urwell/URWellStrip.java index cbfbfdf4fb..d9d2575999 100644 --- a/reconstruction/urwell/src/main/java/org/jlab/service/urwell/URWellStrip.java +++ b/reconstruction/urwell/src/main/java/org/jlab/service/urwell/URWellStrip.java @@ -185,7 +185,7 @@ public static List getStrips(DataEvent event, URWellStripFactory fa strip.setEnergy(strip.ADC*URWellConstants.ADCTOENERGY); strip.setTime(strip.TDC*URWellConstants.TDCTOTIME); strip.setLine(factory.getStrip(sector, layer, comp)); - strip.setChamber(factory.getChamberIndex(comp)+1); + strip.setChamber(factory.getChamberIndex(layer, comp)+1); strip.setStatus(0); if(strip.getEnergy()>URWellConstants.THRESHOLD) strips.add(strip); From 4e6fc6a4fb77c669112efeb33e6c606e38304c9a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mariangela=20Bond=C3=AD?= Date: Mon, 14 Jul 2025 17:31:42 -0400 Subject: [PATCH 10/16] fix issue on the urwell geometry --- .../jlab/detector/geant4/v2/URWELL/URWellGeant4Factory.java | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellGeant4Factory.java b/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellGeant4Factory.java index a20cb42604..66bc6f43ab 100644 --- a/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellGeant4Factory.java +++ b/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellGeant4Factory.java @@ -403,9 +403,9 @@ private double[] getNonProtoDimensions(int ichamber, double h, boolean full) { if (full) { return new double[]{ this.getChamberThickness() / 2. + URWellConstants.ZENLARGEMENT / 2., - halfHeight + 0.05, - dx0 + 0.1, - dx1, + halfHeight + URWellConstants.YENLARGEMENT / 2, + dx0 + URWellConstants.XENLARGEMENT / 2., + dx1 + URWellConstants.XENLARGEMENT / 2., Math.toRadians(thilt) }; } else { From 995883cccf7509ef6d0efa751167ad551f9d0388 Mon Sep 17 00:00:00 2001 From: Raffaella De Vita Date: Mon, 14 Jul 2025 18:48:44 -0400 Subject: [PATCH 11/16] initializing the urwell service for ddvcs --- .../org/jlab/service/urwell/URWellConstants.java | 16 ++++++++-------- .../org/jlab/service/urwell/URWellEngine.java | 2 +- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/reconstruction/urwell/src/main/java/org/jlab/service/urwell/URWellConstants.java b/reconstruction/urwell/src/main/java/org/jlab/service/urwell/URWellConstants.java index a825418809..12237f2d21 100644 --- a/reconstruction/urwell/src/main/java/org/jlab/service/urwell/URWellConstants.java +++ b/reconstruction/urwell/src/main/java/org/jlab/service/urwell/URWellConstants.java @@ -9,14 +9,14 @@ public class URWellConstants { // geometry public final static int NSECTOR = 6; - public final static int NLAYER = 4; - public final static int NREGION = 2; - public final static int NCHAMBER = 3; - public final static int[] NSTRIPS = { 542, 628, 714}; // number of strips for the three chambers - public final static int[] STRIPMIN = { 1, 543, 1171}; // lower strip number - public final static int[] STRIPMAX = { 542, 1170, 1884}; // higher strip number - public final static double PITCH = 0.1; // mm - public final static double[] STEREO = { 10.0, 10.0 }; + public final static int NLAYER = 12; + public final static int NREGION = 6; + public final static int NCHAMBER = 1; +// public final static int[] NSTRIPS = { 542, 628, 714}; // number of strips for the three chambers +// public final static int[] STRIPMIN = { 1, 543, 1171}; // lower strip number +// public final static int[] STRIPMAX = { 542, 1170, 1884}; // higher strip number + public final static double PITCH = 0.5; // mm +// public final static double[] STEREO = { 10.0, 10.0 }; // strips public final static double THRESHOLD = 0; diff --git a/reconstruction/urwell/src/main/java/org/jlab/service/urwell/URWellEngine.java b/reconstruction/urwell/src/main/java/org/jlab/service/urwell/URWellEngine.java index fcb53c5ad1..3ec66baf51 100644 --- a/reconstruction/urwell/src/main/java/org/jlab/service/urwell/URWellEngine.java +++ b/reconstruction/urwell/src/main/java/org/jlab/service/urwell/URWellEngine.java @@ -41,7 +41,7 @@ public boolean init() { // init ConstantsManager to read constants from CCDB String variationName = Optional.ofNullable(this.getEngineConfigString("variation")).orElse("default"); DatabaseConstantProvider cp = new DatabaseConstantProvider(11, variationName); - factory.init(cp, false, URWellConstants.NREGION); + factory.init(cp, 2, URWellConstants.NREGION); // register output banks for drop option this.registerOutputBank("URWELL::hits"); this.registerOutputBank("URWELL::clusters"); From 51b526938c3cbbb6d2b7803bff847dbb361e33b5 Mon Sep 17 00:00:00 2001 From: Raffaella De Vita Date: Mon, 14 Jul 2025 18:50:12 -0400 Subject: [PATCH 12/16] added dc frame --- .../detector/geant4/v2/DCGeant4Factory.java | 29 ++++++++++++++----- 1 file changed, 21 insertions(+), 8 deletions(-) diff --git a/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/DCGeant4Factory.java b/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/DCGeant4Factory.java index 28e5782faf..236af2503a 100644 --- a/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/DCGeant4Factory.java +++ b/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/DCGeant4Factory.java @@ -580,6 +580,7 @@ public final class DCGeant4Factory extends Geant4Factory { private final HashMap properties = new HashMap<>(); private int nsgwires; + private final double x_enlargement = 1.00; private final double y_enlargement = 3.65; private final double z_enlargement = -2.46; private final double microgap = 0.01; @@ -825,19 +826,20 @@ private Geant4Basic getRegion(int isec, int ireg) { } private Geant4Basic getLayer(int isec, int isuper, int ilayer) { - return getRegion(isec, isuper/2).getChildren().get((isuper%2)*6 + ilayer); + return getRegion(isec, isuper/2).getChildren().get(0).getChildren().get((isuper%2)*6 + ilayer); } /////////////////////////////////////////////////// public Geant4Basic createRegion(int isector, int iregion) { Wire regw0 = new Wire(isector+1, iregion * 2, 0, 0); Wire regw1 = new Wire(isector+1, iregion * 2 + 1, 7, nsgwires - 1); - double dx_shift = y_enlargement * Math.tan(Math.toRadians(29.5)); - + double dx_shift = y_enlargement * Math.tan(dbref.thopen(iregion)/2); + double extra = 0.2; + double reg_dz = (dbref.frontgap(iregion) + dbref.backgap(iregion) + dbref.midgap(iregion) + dbref.superwidth(iregion * 2) + dbref.superwidth(iregion * 2 + 1)) / 2.0 + z_enlargement; - double reg_dx0 = Math.abs(regw0.bottom().x) - dx_shift + 1.0; - double reg_dx1 = Math.abs(regw1.top().x) + dx_shift + 1.0; - double reg_dy = regw1.top().minus(regw0.bottom()).y / Math.cos(dbref.thtilt(iregion)) / 2.0 + y_enlargement + 1.0; + double reg_dx0 = Math.abs(regw0.bottom().x) - dx_shift + x_enlargement; + double reg_dx1 = Math.abs(regw1.top().x) + dx_shift + x_enlargement; + double reg_dy = regw1.top().minus(regw0.bottom()).y / Math.cos(dbref.thtilt(iregion)) / 2.0 + y_enlargement; double reg_skew = 0.0; double reg_thtilt = dbref.thtilt(iregion); @@ -854,6 +856,17 @@ public Geant4Basic createRegion(int isector, int iregion) { regionVolume.translate(vcenter.x, vcenter.y, vcenter.z); regionVolume.setId(isector + 1, iregion + 1, 0, 0); + double gas_dx0 = reg_dx0 - x_enlargement + dx_shift + extra; + double gas_dx1 = reg_dx1 - x_enlargement - dx_shift + extra; + double gas_dy = reg_dy - y_enlargement + extra; + Geant4Basic regionGas = new G4Trap("regionGas" + (iregion + 1) + "_s" + (isector + 1), + reg_dz, -reg_thtilt, Math.toRadians(90.0), + gas_dy, gas_dx0, gas_dx1, 0.0, + gas_dy, gas_dx0, gas_dx1, 0.0); + regionGas.setPosition(0, 0, 0); + regionGas.setMother(regionVolume); + regionGas.setId(isector + 1, iregion + 1, 0, 0); + for (int isup = 0; isup < 2; isup++) { int isuper = iregion * 2 + isup; Geant4Basic superlayerVolume = this.createSuperlayer(isuper); @@ -866,7 +879,7 @@ public Geant4Basic createRegion(int isector, int iregion) { superlayerVolume.rotate("zxy", -dbref.thster(isuper), 0.0, 0.0); superlayerVolume.setPosition(slshift.x, slshift.y, slshift.z); - superlayerVolume.setMother(regionVolume); + superlayerVolume.setMother(regionGas); superlayerVolume.setId(isector + 1, iregion + 1, isuper + 1); int nsglayers = dbref.nsenselayers(isuper) + dbref.nguardlayers(isuper); @@ -881,7 +894,7 @@ public Geant4Basic createRegion(int isector, int iregion) { layerVolume.rotate("zxy", -dbref.thster(isuper), 0.0, 0.0); layerVolume.setPosition(lshift.x, lshift.y, lshift.z); - layerVolume.setMother(regionVolume); + layerVolume.setMother(regionGas); layerVolume.setId(isector + 1, iregion + 1, isuper + 1, ilayer); } } From a89eb257501d73b5c36f7d8f17c4ced131cfd84e Mon Sep 17 00:00:00 2001 From: tongtongcao Date: Wed, 16 Jul 2025 14:55:14 -0400 Subject: [PATCH 13/16] update FMT tracking to work on FVT --- reconstruction/fmt/pom.xml | 5 + .../main/java/org/jlab/rec/fmt/Constants.java | 2 +- .../jlab/rec/fmt/banks/RecoBankWriter.java | 18 +- .../java/org/jlab/rec/fmt/track/Track.java | 212 +++++------------- .../org/jlab/rec/fmt/track/fit/KFitter.java | 104 ++++----- .../org/jlab/rec/fmt/track/fit/MeasVecs.java | 90 ++++++-- .../jlab/rec/fmt/track/fit/RungeKutta.java | 18 +- .../org/jlab/rec/fmt/track/fit/StateVecs.java | 18 +- .../java/org/jlab/service/fmt/FMTEngine.java | 147 ++++-------- .../jlab/rec/urwell/reader/URWellCluster.java | 105 +++++++++ .../jlab/rec/urwell/reader/URWellCross.java | 189 ++++++++++++++++ .../org/jlab/rec/urwell/reader/URWellHit.java | 43 ++++ .../jlab/rec/urwell/reader/URWellReader.java | 144 ++++++++++++ 13 files changed, 724 insertions(+), 371 deletions(-) create mode 100644 reconstruction/urwell/src/main/java/org/jlab/rec/urwell/reader/URWellCluster.java create mode 100644 reconstruction/urwell/src/main/java/org/jlab/rec/urwell/reader/URWellCross.java create mode 100644 reconstruction/urwell/src/main/java/org/jlab/rec/urwell/reader/URWellHit.java create mode 100644 reconstruction/urwell/src/main/java/org/jlab/rec/urwell/reader/URWellReader.java diff --git a/reconstruction/fmt/pom.xml b/reconstruction/fmt/pom.xml index d550549ccf..6dfa6d77a1 100644 --- a/reconstruction/fmt/pom.xml +++ b/reconstruction/fmt/pom.xml @@ -14,6 +14,11 @@ + + org.jlab.clas12.detector + clas12detector-urwell + 13.0.2-SNAPSHOT + org.jlab.jnp jnp-hipo diff --git a/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/Constants.java b/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/Constants.java index 9e38a007d5..08969f9af1 100644 --- a/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/Constants.java +++ b/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/Constants.java @@ -13,7 +13,7 @@ */ public class Constants { - public static final int NLAYERS = 6; + public static final int NLAYERS = 12; // DC-tracks to FMT-clusters matching parameter public static double CIRCLECONFUSION = 1.2; // cm diff --git a/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/banks/RecoBankWriter.java b/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/banks/RecoBankWriter.java index 9a0ac60374..87aed80467 100644 --- a/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/banks/RecoBankWriter.java +++ b/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/banks/RecoBankWriter.java @@ -113,6 +113,7 @@ private static DataBank fillFMTTracksBank(DataEvent event,List candlist) return bank; } + /* private static DataBank fillFMTTrajectoryBank(DataEvent event,List candlist) { DataBank bank = event.createBank("FMT::Trajectory", candlist.size()*Constants.NLAYERS); @@ -147,27 +148,14 @@ private static DataBank fillFMTTrajectoryBank(DataEvent event,List candli } return bank; } + */ - public static void appendFMTBanks(DataEvent event, List fhits, List clusters, - List tracks) { + public static void appendFMTBanks(DataEvent event, List tracks) { if (event == null) return; - if (fhits != null) { - event.appendBanks(fillFMTHitsBank(event, fhits)); - } - - if (clusters != null && clusters.size()>0) { - event.appendBanks(fillFMTClustersBank(event, clusters)); - } - -// if (crosses != null && crosses.size() > 0) { -// event.appendBanks(this.fillFMTCrossesBank(event, crosses)); -// } - if (tracks != null && tracks.size() > 0) { event.appendBanks(fillFMTTracksBank(event, tracks)); - event.appendBanks(fillFMTTrajectoryBank(event, tracks)); } } diff --git a/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/Track.java b/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/Track.java index 4b6bb58d2a..6167ba64f7 100644 --- a/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/Track.java +++ b/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/Track.java @@ -7,11 +7,13 @@ import java.util.Map.Entry; import org.jlab.clas.swimtools.Swim; import org.jlab.geom.prim.Line3D; +import org.jlab.geom.prim.Point3D; import org.jlab.geom.prim.Vector3D; import org.jlab.io.base.DataBank; import org.jlab.io.base.DataEvent; import org.jlab.rec.fmt.Constants; -import org.jlab.rec.fmt.cluster.Cluster; + +import org.jlab.rec.urwell.reader.URWellCluster; /** * @@ -44,8 +46,7 @@ public class Track { private double _pz; private int _NDF; - private final Trajectory[] _DCtrajs = new Trajectory[Constants.NLAYERS]; - private final List[] _clusters = new ArrayList[Constants.NLAYERS]; + private final List[] _clusters = new ArrayList[Constants.NLAYERS]; private final Trajectory[] _FMTtrajs = new Trajectory[Constants.NLAYERS]; public Track() { @@ -53,7 +54,7 @@ public Track() { public Track(int _index, int _sector, int _q, double _x, double _y, double _z, - double _px, double _py, double _pz, List clusters) { + double _px, double _py, double _pz, List clusters) { this._index = _index; this._sector = _sector; this._q = _q; @@ -63,41 +64,24 @@ public Track(int _index, int _sector, int _q, double _x, double _y, double _z, this._px = _px; this._py = _py; this._pz = _pz; - for(Cluster cluster : clusters) this.addCluster(cluster); + for(URWellCluster cluster : clusters) this.addCluster(cluster); } - - - /** - * @param layer - * @return the _traj - */ - public Trajectory getDCTraj(int layer) { - if(layer<=0 || layer>Constants.NLAYERS) return null; - else return _DCtrajs[layer-1]; - } - /** - * @param trj - */ - public void setDCTraj(Trajectory trj) { - this._DCtrajs[trj.getLayer()-1] = trj; - } - - public List getClusters() { - List clusters = new ArrayList<>(); + public List getClusters() { + List clusters = new ArrayList<>(); for(int i=0; i getClusters(int layer) { + public List getClusters(int layer) { if(layer<=0 || layer>Constants.NLAYERS) return null; else return _clusters[layer-1]; } - public Cluster getCluster(int layer) { + public URWellCluster getCluster(int layer) { if(layer<=0 || layer>Constants.NLAYERS) return null; else if(_clusters[layer-1]== null || _clusters[layer-1].size()==0) return null; else return _clusters[layer-1].get(0); @@ -115,12 +99,12 @@ public int getClusterLayer(int layer) { if(_clusters[layer-1]!=null) return _clusters[layer-1].size(); else return 0; } - - public final void addCluster(Cluster cluster) { - if(this._clusters[cluster.getLayer()-1]==null) - this._clusters[cluster.getLayer()-1] = new ArrayList<>(); - this._clusters[cluster.getLayer()-1].add(cluster); - } + + public final void addCluster(URWellCluster cluster) { + if(this._clusters[cluster.layer()-1]==null) + this._clusters[cluster.layer()-1] = new ArrayList<>(); + this._clusters[cluster.layer()-1].add(cluster); + } public void clearClusters(int layer) { this._clusters[layer-1].clear(); @@ -305,71 +289,8 @@ public double getPz() { public void setPz(double _pz) { this._pz = _pz; } - - // FIXME: THIS METHOD SHOULD BE GENERALIZED - public double getSeedQuality() { - double quality=99; - if(this.getClusterLayer(1)>0 && this.getClusterLayer(2)>0 && this.getClusterLayer(3)>0) { - Line3D seg1 = this.getClusters(1).get(0).getGlobalSegment(); - Line3D seg2 = this.getClusters(2).get(0).getGlobalSegment(); - Line3D seg3 = this.getClusters(3).get(0).getGlobalSegment(); - quality = seg1.distanceSegments(seg3).distanceSegments(seg2).length(); - } - return quality; - } - - // FIXME: THIS METHOD SHOULD BE GENERALIZED - public void filterClusters(int mode) { - - double dref = Double.POSITIVE_INFINITY; - - int[] ibest = new int[Constants.NLAYERS]; - for(int i=0; i0 && this.getClusterLayer(2)>0 && this.getClusterLayer(3)>0) { - for(int i1=0; i1=0) { - List cls = new ArrayList<>(); - cls.add(_clusters[i].get(ibest[i])); - _clusters[i] = cls; - } - } - } - - public static List getDCTracks(DataEvent event, Swim swimmer) { + public List getDCTracks(DataEvent event, Swim swimmer) { Map trackmap = new LinkedHashMap(); @@ -386,48 +307,29 @@ public static List getDCTracks(DataEvent event, Swim swimmer) { trk.setIndex(i); trk.setSector(trackBank.getByte("sector", i)); trk.setQ(trackBank.getByte("q", i)); - trk.setX(trackBank.getFloat("Vtx0_x", i)); - trk.setY(trackBank.getFloat("Vtx0_y", i)); - trk.setZ(trackBank.getFloat("Vtx0_z", i)); - trk.setPx(trackBank.getFloat("p0_x", i)); - trk.setPy(trackBank.getFloat("p0_y", i)); - trk.setPz(trackBank.getFloat("p0_z", i)); + + double vx = trackBank.getFloat("Vtx0_x", i); + double vy = trackBank.getFloat("Vtx0_y", i); + double vz = trackBank.getFloat("Vtx0_z", i); + + double px = trackBank.getFloat("p0_x", i); + double py = trackBank.getFloat("p0_y", i); + double pz = trackBank.getFloat("p0_z", i); + + + Point3D vertexLocal = transGlobaltoLocal(trk.getSector(), vx, vy, vz); + + Point3D momLocal = transGlobaltoLocal(trk.getSector(), px, py, pz); + + trk.setX(vertexLocal.x()); + trk.setY(vertexLocal.y()); + trk.setZ(vertexLocal.z()); + trk.setPx(momLocal.x()); + trk.setPy(momLocal.y()); + trk.setPz(momLocal.z()); trk.setStatus(1); - trackmap.put(id,trk); - - for(int j=0; j tracks = new ArrayList<>(); for(Entry entry: trackmap.entrySet()) { @@ -436,39 +338,33 @@ public static List getDCTracks(DataEvent event, Swim swimmer) { return tracks; } - private static double[] getTrajectory(double x, double y, double z, double px, double py, double pz, - int q, int layer, Swim swim) { - Vector3D p = Constants.getLayer(layer).getPlane().point().toVector3D(); - Vector3D n = Constants.getLayer(layer).getPlane().normal(); - Vector3D v = new Vector3D(x, y, z); - double d = p.dot(n); - if(v.dot(n) clusters; + private List clusters; public KFitter(Track track, Swim swimmer, int c) { sv = new StateVecs(swimmer); @@ -46,7 +47,7 @@ public KFitter(Track track, Swim swimmer, int c) { track.getPx(), track.getPy(), track.getPz(), track.getQ(), c); } - public KFitter(List clusters, int sector, double xVtx, double yVtx, double zVtx, + public KFitter(List clusters, int sector, double xVtx, double yVtx, double zVtx, double pxVtx, double pyVtx, double pzVtx, int q, Swim swimmer, int c) { sv = new StateVecs(swimmer); this.track = new Track(0,sector, q, xVtx, yVtx, zVtx, pxVtx, pyVtx, pzVtx, clusters); @@ -54,7 +55,7 @@ public KFitter(List clusters, int sector, double xVtx, double yVtx, dou this.init(clusters, sector, xVtx, yVtx, zVtx, pxVtx, pyVtx, pzVtx, q, c); } - private void init(List clusters, int sector, double xVtx, double yVtx, double zVtx, + private void init(List clusters, int sector, double xVtx, double yVtx, double zVtx, double pxVtx, double pyVtx, double pzVtx, int q, int c) { // initialize measVecs mv.setMeasVecs(clusters); @@ -65,7 +66,7 @@ private void init(List clusters, int sector, double xVtx, double yVtx, for (int i = 0; i < mSize; i++) sv.Z[i] = mv.measurements.get(i).z; // initialize stateVecs - sv.init(sector, xVtx, yVtx, zVtx, pxVtx, pyVtx, pzVtx, q, sv.Z[0], this, c); + sv.init(sector, xVtx, yVtx, zVtx, pxVtx, pyVtx, pzVtx, q, sv.Z[0], this, c); } public void runFitter(int sector) { @@ -115,6 +116,7 @@ public void runFitter(int sector) { // save final trajectory points + /* if(this.finalStateVec!=null) { for (int k = 0; k < svzLength; ++k) { Trajectory trj = new Trajectory(mv.measurements.get(k).layer, @@ -128,51 +130,29 @@ public void runFitter(int sector) { track.setFMTtraj(trj); } } + */ -// for (int li = 1; li <= 6; ++li) { -// // Get the state vector closest in z to the FMT layer. -// // NOTE: A simple optimization would be to do this on another for loop with only the -// // state vectors, saving the ones closest to the FMT layers to use them later -// // instead of looping through them 6 times. -// int closestSVID = -1; -// double closestSVDistance = Double.POSITIVE_INFINITY; -// for (int si = 0; si < sv.trackTraj.size(); ++si) { -// double svDistance = Math.abs(sv.trackTraj.get(si).z - Geometry.getLayerZ(li-1)); -// if (svDistance < closestSVDistance) { -// closestSVID = si; -// closestSVDistance = svDistance; -// } -// } -// -// // Get the state vector's y position in the layer's local coordinates. -// Point3D Pos = new Point3D(sv.trackTraj.get(closestSVID).x, sv.trackTraj.get(closestSVID).y, 0); -// Point3D locPos = Geometry.globalToLocal(Pos,li); -// -// // Store the cluster's residual. -// for (Cluster cl : this.clusters) { -// if (cl.get_Layer() == li) { -// cl.set_CentroidResidual(cl.get_Centroid() - locPos.y()); -// } -// } -// } } - + public Matrix filterCovMat(double[] H, Matrix Ci, double V) { + double det = Matrix5x5.inverse(Ci, first_inverse, adj); - if (Math.abs(det) < 1.e-30) return null; + if (Math.abs(det) < 1.e-60) { + return null; + } addition.set( - H[0] * H[0] / V, H[0] * H[1] / V, H[0] * H[2] / V, H[0] * H[3] / V, H[0] * H[4] / V, - H[1] * H[0] / V, H[1] * H[1] / V, H[1] * H[2] / V, H[1] * H[3] / V, H[1] * H[4] / V, - H[2] * H[0] / V, H[2] * H[1] / V, H[2] * H[2] / V, H[2] * H[3] / V, H[2] * H[4] / V, - H[3] * H[0] / V, H[3] * H[1] / V, H[3] * H[2] / V, H[3] * H[3] / V, H[3] * H[4] / V, - H[4] * H[0] / V, H[4] * H[1] / V, H[4] * H[2] / V, H[4] * H[3] / V, H[4] * H[4] / V - ); + H[0] * H[0] / V, H[0] * H[1] / V, 0, 0, 0, + H[0] * H[1] / V, H[1] * H[1] / V, 0, 0, 0, + 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0); Matrix5x5.add(first_inverse, addition, result); double det2 = Matrix5x5.inverse(result, result_inv, adj); - - if (Math.abs(det2) < 1.e-30) return null; + if (Math.abs(det2) < 1.e-60) { + return null; + } return result_inv; } @@ -180,29 +160,41 @@ public Matrix filterCovMat(double[] H, Matrix Ci, double V) { private void filter(int k) { if (sv.trackTraj.get(k) != null && sv.trackCov.get(k).covMat != null && k < sv.Z.length ) { double[] K = new double[5]; - double V = mv.measurements.get(k).error; + double V = mv.measurements.get(k).error * mv.measurements.get(k).error; - double[] H = mv.H(sv.trackTraj.get(k),sv); + double[] H = mv.HURWell(sv.trackTraj.get(k),sv); Matrix CaInv = this.filterCovMat(H, sv.trackCov.get(k).covMat, V); if (CaInv != null) sv.trackCov.get(k).covMat = CaInv; else return; + Matrix cMat = new Matrix(); + if (CaInv != null) { + Matrix5x5.copy(CaInv, cMat); + } else { + return; + } + // Calculate the gain matrix. for (int j = 0; j < 5; j++) { - K[j] = 0; - for (int i = 0; i < 5; i++) K[j] += H[i] * sv.trackCov.get(k).covMat.get(j, i) / V ; + // the gain matrix + K[j] = (H[0] * cMat.get(j, 0) + + H[1] * cMat.get(j, 1)) / V; } - // Update Chi^2 and filtered state vector. - double h = mv.h(sv.trackTraj.get(k)); - double m = mv.measurements.get(k).centroid; - this.chi2 += ((m-h)*(m-h)/mv.measurements.get(k).error/mv.measurements.get(k).error); - - double x_filt = sv.trackTraj.get(k).x + K[0] * (m-h); - double y_filt = sv.trackTraj.get(k).y + K[1] * (m-h); - double tx_filt = sv.trackTraj.get(k).tx + 0*K[2] * (m-h); - double ty_filt = sv.trackTraj.get(k).ty + 0*K[3] * (m-h); - double Q_filt = sv.trackTraj.get(k).Q + 0*K[4] * (m-h); + // Update Chi^2 and filtered state vector. + double res = mv.dhURWell(sv.trackTraj.get(k)); + double filt[] = new double[5]; + for(int j = 0; j < 5; j ++){ + filt[j] += K[j]*res; + } + + this.chi2 += (res*res/mv.measurements.get(k).error/mv.measurements.get(k).error); + + double x_filt = sv.trackTraj.get(k).x + filt[0]; + double y_filt = sv.trackTraj.get(k).y + filt[1]; + double tx_filt = sv.trackTraj.get(k).tx + filt[2]; + double ty_filt = sv.trackTraj.get(k).ty + filt[3]; + double Q_filt = sv.trackTraj.get(k).Q + filt[4]; if (filterOn) { sv.trackTraj.get(k).x = x_filt; diff --git a/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/fit/MeasVecs.java b/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/fit/MeasVecs.java index dba901877d..381b659b9b 100644 --- a/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/fit/MeasVecs.java +++ b/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/fit/MeasVecs.java @@ -4,9 +4,13 @@ import java.util.List; import org.jlab.geom.prim.Vector3D; import org.jlab.rec.fmt.Constants; -import org.jlab.rec.fmt.cluster.Cluster; +import org.jlab.rec.urwell.reader.URWellCluster; import org.jlab.rec.fmt.track.fit.StateVecs.StateVec; +import org.jlab.geom.prim.Point3D; +import org.jlab.geom.prim.Line3D; +import org.jlab.geom.prim.Plane3D; + /** * @author ziegler */ @@ -17,7 +21,8 @@ public class MeasVecs { public class MeasVec implements Comparable { public double z = Double.NaN; - public double centroid; + public Point3D lineEndPoint1 = null; + public Point3D lineEndPoint2 = null; public double seed; public double error; public int layer; @@ -33,47 +38,84 @@ public int compareTo(MeasVec arg) { } } - public void setMeasVecs(List clusters) { + public void setMeasVecs(List clusters) { measurements = new ArrayList<>(); for (int i = 0; i < clusters.size(); i++) { - int l = clusters.get(i).getLayer()-1; - double cent = clusters.get(i).getCentroid(); - double error = clusters.get(i).getCentroidError(); - double z = clusters.get(i).getGlobalSegment().origin().z(); - int seed = clusters.get(i).getSeedStrip(); - MeasVec meas = this.setMeasVec(l, cent, error, z, seed); + int l = clusters.get(i).layer()-1; + Point3D lineEndPoint1 = clusters.get(i).getLineLocal().origin(); + Point3D lineEndPoint2 = clusters.get(i).getLineLocal().end(); + double error = 0.0144; // = pitch/sqrt(12), where pitch is 500 um + double z = lineEndPoint1.z(); + int seed = clusters.get(i).strip(); + MeasVec meas = this.setMeasVec(l, lineEndPoint1, lineEndPoint2, error, z, seed); measurements.add(meas); } } - public MeasVec setMeasVec(int l, double cent, double error, double z, int seed) { + public MeasVec setMeasVec(int l, Point3D lineEndPoint1, Point3D lineEndPoint2, double error, double z, int seed) { MeasVec meas = new MeasVec(); meas.layer = l+1; - meas.centroid = cent; meas.error = error; meas.z = z; meas.seed = seed; + meas.lineEndPoint1 = lineEndPoint1; + meas.lineEndPoint2 = lineEndPoint2; + return meas; } - public double h(StateVec stateVec) { - if (stateVec == null) return 0; - if (this.measurements.get(stateVec.k) == null) return 0; - - int layer = this.measurements.get(stateVec.k).layer; - - return Constants.toLocal(layer, stateVec.x, stateVec.y, stateVec.z).y(); + public double dhURWell(StateVec stateVec) { + double value = Double.NaN; + if (stateVec == null|| this.measurements.get(stateVec.k) == null) { + return value; + } + + Line3D l = new Line3D(this.measurements.get(stateVec.k).lineEndPoint1, + this.measurements.get(stateVec.k).lineEndPoint2); + Line3D WL = new Line3D(); + WL.copy(l); + Point3D svP = new Point3D(stateVec.x, stateVec.y, stateVec.z); + WL.copy(WL.distance(svP)); + Plane3D plane = new Plane3D(0, 0, this.measurements.get(stateVec.k).z, 0, 0, 1); // plan perpenticular to z axis in TSC + double sideStrip = -Math.signum(l.direction().cross(WL.direction()). + dot(plane.normal())); + value = WL.length()*sideStrip; + + return value; } + + public double[] HURWell(StateVec stateVec, StateVecs sv) { + double[] hMatrix = new double[5]; + double Err = 0.01; + double[][] Result = new double[2][2]; + for (int i = 0; i < 2; i++) { + StateVec svc = sv.new StateVec(stateVec.k); + svc.x = stateVec.x; + svc.y = stateVec.y; + svc.z = stateVec.z; + svc.x = stateVec.x + (double) Math.pow(-1, i) * Err; + Result[i][0] = dhURWell(svc); + } + for (int i = 0; i < 2; i++) { + StateVec svc = sv.new StateVec(stateVec.k); + svc.x = stateVec.x; + svc.y = stateVec.y; + svc.z = stateVec.z; + svc.y = stateVec.y + (double) Math.pow(-1, i) * Err; + Result[i][1] = dhURWell(svc); + } - public double[] H(StateVec stateVec, StateVecs sv) { - int layer = this.measurements.get(stateVec.k).layer; - Vector3D derivatives = Constants.getDerivatives(layer, stateVec.x, stateVec.y, stateVec.z); - double[] H = new double[]{derivatives.x(), derivatives.y(), 0, 0, 0}; - return H; - } + hMatrix[0] = -(Result[0][0] - Result[1][0]) / (2. * Err); // Add negative sign since dh = meas - h; here use dh to replace h since meas is cancelled when derivative + hMatrix[1] = -(Result[0][1] - Result[1][1]) / (2. * Err); // Add negative sign since dh = meas - h; here use dh to replace h since meas is cancelled when derivative + hMatrix[2] = 0; + hMatrix[3] = 0; + hMatrix[4] = 0; + + return hMatrix; + } private StateVec reset(StateVec SVplus, StateVec stateVec, StateVecs sv) { SVplus = sv.new StateVec(stateVec.k); diff --git a/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/fit/RungeKutta.java b/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/fit/RungeKutta.java index 809acc4fde..21162bd99d 100644 --- a/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/fit/RungeKutta.java +++ b/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/fit/RungeKutta.java @@ -34,7 +34,7 @@ public RungeKutta() { public void SwimToZ(int sector, StateVecs.StateVec fVec, Swim dcSwim, double z0, float[] bf){ double stepSize = 0.5; - dcSwim.BfieldLab(fVec.x, fVec.y, fVec.z, bf); + dcSwim.Bfield(sector, fVec.x, fVec.y, fVec.z, bf); fVec.B = Math.sqrt(bf[0]*bf[0]+bf[1]*bf[1]+bf[2]*bf[2]); double s = fVec.B; @@ -67,25 +67,25 @@ public void SwimToZ(int sector, StateVecs.StateVec fVec, Swim dcSwim, double z0, void RK4transport(int sector, double q, double x0, double y0, double z0, double tx0, double ty0, double h, Swim swimmer, double dPath, StateVecs.StateVec fVec) { // lab system = 1, TSC =0 - swimmer.BfieldLab(x0, y0, z0, _b); + swimmer.Bfield(sector, x0, y0, z0, _b); double x1 = tx0; double y1 = ty0; double tx1=q*v*Ax(tx0, ty0, _b[0], _b[1], _b[2]); double ty1=q*v*Ay(tx0, ty0, _b[0], _b[1], _b[2]); - swimmer.BfieldLab(x0+0.5*h*x1, y0+0.5*h*y1, z0+0.5*h, _b); + swimmer.Bfield(sector, x0+0.5*h*x1, y0+0.5*h*y1, z0+0.5*h, _b); double x2 = tx0+0.5*h*tx1; double y2 = ty0+0.5*h*ty1; double tx2=q*v*Ax((tx0+0.5*h*tx1), (ty0+0.5*h*ty1), _b[0], _b[1], _b[2]); double ty2=q*v*Ay((tx0+0.5*h*tx1), (ty0+0.5*h*ty1), _b[0], _b[1], _b[2]); - swimmer.BfieldLab(x0+0.5*h*x2, y0+0.5*h*y2, z0+0.5*h, _b); + swimmer.Bfield(sector, x0+0.5*h*x2, y0+0.5*h*y2, z0+0.5*h, _b); double x3 = tx0+0.5*h*tx2; double y3 = ty0+0.5*h*ty2; double tx3=q*v*Ax((tx0+0.5*h*tx2), (ty0+0.5*h*ty2), _b[0], _b[1], _b[2]); double ty3=q*v*Ay((tx0+0.5*h*tx2), (ty0+0.5*h*ty2), _b[0], _b[1], _b[2]); - swimmer.BfieldLab(x0+h*x3, y0+h*y3, z0+h, _b); + swimmer.Bfield(sector, x0+h*x3, y0+h*y3, z0+h, _b); double x4 = tx0+h*tx3; double y4 = ty0+h*ty3; double tx4=q*v*Ax((tx0+h*tx3), (ty0+h*ty3), _b[0], _b[1], _b[2]); @@ -125,7 +125,7 @@ void RK4transport(int sector, double q, double x0, double y0, double z0, double double delty_delq0_0 =0; //System.out.println("RK0 "+x0+","+y0+","+z0+";"+tx0+","+ty0+","+" z0 "+z0+" h "+h); //State - swimmer.BfieldLab(x0, y0, z0, _b); + swimmer.Bfield(sector, x0, y0, z0, _b); double x1 = tx0; double y1 = ty0; double tx1=q*v*Ax(tx0, ty0, _b[0], _b[1], _b[2]); @@ -157,7 +157,7 @@ void RK4transport(int sector, double q, double x0, double y0, double z0, double + delAy_delty(tx0,ty0,_b[0],_b[1],_b[2])*delty_delq0_0); - swimmer.BfieldLab(x0+0.5*h*x1, y0+0.5*h*y1, z0+0.5*h, _b); + swimmer.Bfield(sector, x0+0.5*h*x1, y0+0.5*h*y1, z0+0.5*h, _b); double x2 = tx0+0.5*h*tx1; double y2 = ty0+0.5*h*ty1; double tx2=q*v*Ax((tx0+0.5*h*tx1), (ty0+0.5*h*ty1), _b[0], _b[1], _b[2]); @@ -186,7 +186,7 @@ void RK4transport(int sector, double q, double x0, double y0, double z0, double double delty_delq0_2 = this.delty_delq0_next(q,v,tx0+0.5*h*tx1,ty0+0.5*h*ty1,_b[0],_b[1],_b[2], deltx_delq0_0+0.5*h*deltx_delq0_1,delty_delq0_0+0.5*h*delty_delq0_1); - swimmer.BfieldLab(x0+0.5*h*x2, y0+0.5*h*y2, z0+0.5*h, _b); + swimmer.Bfield(sector, x0+0.5*h*x2, y0+0.5*h*y2, z0+0.5*h, _b); double x3 = tx0+0.5*h*tx2; double y3 = ty0+0.5*h*ty2; double tx3=q*v*Ax((tx0+0.5*h*tx2), (ty0+0.5*h*ty2), _b[0], _b[1], _b[2]); @@ -215,7 +215,7 @@ void RK4transport(int sector, double q, double x0, double y0, double z0, double double delty_delq0_3 = this.delty_delq0_next(q,v,tx0+0.5*h*tx2,ty0+0.5*h*ty2,_b[0],_b[1],_b[2], deltx_delq0_0+0.5*h*deltx_delq0_2,delty_delq0_0+0.5*h*delty_delq0_2); - swimmer.BfieldLab(x0+h*x3, y0+h*y3, z0+h, _b); + swimmer.Bfield(sector, x0+h*x3, y0+h*y3, z0+h, _b); double x4 = tx0+h*tx3; double y4 = ty0+h*ty3; double tx4=q*v*Ax((tx0+h*tx3), (ty0+h*ty3), _b[0], _b[1], _b[2]); diff --git a/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/fit/StateVecs.java b/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/fit/StateVecs.java index 620faf7214..dad52c309e 100644 --- a/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/fit/StateVecs.java +++ b/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/fit/StateVecs.java @@ -12,6 +12,8 @@ */ public class StateVecs { private final double Bmax = 2.366498; // averaged + + private final double MuonMass = 0.105658; // GeV/c^2 final double speedLight = 0.002997924580; public double[] Z; @@ -29,7 +31,7 @@ public class StateVecs { private final float[] lbf = new float[3]; private final Swim dcSwim; private final RungeKutta rk; - + /** * State vector representing the track in the sector coordinate system at the measurement layer * @param swimmer @@ -90,7 +92,8 @@ public Matrix transport(int sector, int i, double Zf, StateVec iVec, CovMat covM double X0 = this.getX0(z); double t_ov_X0 = Math.abs(s) / X0;//path length in radiation length units = t/X0 [true path length/ X0] ; Ar radiation length = 14 cm - double beta = this.beta; + double energy = Math.sqrt(p*p + MuonMass * MuonMass); + double beta = p/energy; if(beta>1.0 || beta<=0) beta =1.0; @@ -361,6 +364,7 @@ public void init(int sector, double xVtx, double yVtx, double zVtx, rinitSV.tx = pxVtx/pzVtx; rinitSV.ty = pyVtx/pzVtx; rinitSV.Q = (double)q / p; + /* double[] FTF = new double[25]; double[] F = this.F(sector, z0, rinitSV); for(int i = 0; i<5; i++) { @@ -369,6 +373,16 @@ public void init(int sector, double xVtx, double yVtx, double zVtx, Matrix initCMatrix = new Matrix(); initCMatrix.set(FTF); initCM.covMat = initCMatrix; + */ + Matrix initCMatrix = new Matrix(); + initCMatrix.set(8 * 8, 0, 0, 0, 0, + 0, 8 * 8, 0, 0, 0, + 0, 0, 0.1 * 0.1, 0, 0, + 0, 0, 0, 0.1 * 0.1, 0, + 0, 0, 0, 0, 0.03 * 0.03 + ); + initCM.covMat = initCMatrix; + this.trackCov.put(0, initCM); } private StateVec reset(StateVec SVplus, StateVec stateVec) { diff --git a/reconstruction/fmt/src/main/java/org/jlab/service/fmt/FMTEngine.java b/reconstruction/fmt/src/main/java/org/jlab/service/fmt/FMTEngine.java index de6e67e4d5..9b19bf2344 100644 --- a/reconstruction/fmt/src/main/java/org/jlab/service/fmt/FMTEngine.java +++ b/reconstruction/fmt/src/main/java/org/jlab/service/fmt/FMTEngine.java @@ -2,23 +2,26 @@ import java.util.Arrays; import java.util.List; +import java.util.ArrayList; import org.jlab.clas.reco.ReconstructionEngine; import org.jlab.clas.swimtools.Swim; import org.jlab.detector.base.DetectorType; import org.jlab.detector.base.GeometryFactory; +import org.jlab.geom.prim.Point3D; import org.jlab.utils.groups.IndexedTable; import org.jlab.io.base.*; import org.jlab.rec.fmt.Constants; import org.jlab.rec.fmt.banks.RecoBankWriter; -import org.jlab.rec.fmt.cluster.Cluster; -import org.jlab.rec.fmt.hit.Hit; import org.jlab.rec.fmt.track.Track; import org.jlab.rec.fmt.track.Trajectory; import org.jlab.rec.fmt.track.fit.KFitter; import org.jlab.rec.fmt.track.fit.StateVecs.StateVec; +import org.jlab.rec.urwell.reader.URWellCluster; +import org.jlab.rec.urwell.reader.URWellReader; + /** * Service to return reconstructed track candidates - the output is in hipo format * @@ -46,7 +49,6 @@ public boolean init() { System.out.println("["+this.getName()+"] run with FMT default geometry"); } - String[] tables = new String[]{ "/geometry/beam/position", "/calibration/mvt/fmt_time", @@ -57,12 +59,9 @@ public boolean init() { // Load the geometry int run = 10; - Constants.setDetector(GeometryFactory.getDetector(DetectorType.FMT,run, variation)); - + Constants.setDetector(GeometryFactory.getDetector(DetectorType.FMT,run, variation)); + // Register output banks - super.registerOutputBank("FMT::Hits"); - super.registerOutputBank("FMT::Clusters"); - super.registerOutputBank("FMT::Crosses"); super.registerOutputBank("FMT::Tracks"); super.registerOutputBank("FMT::Trajectory"); @@ -78,7 +77,7 @@ public boolean processDataEvent(DataEvent event) { DataBank runConfig = event.getBank("RUN::config"); if (runConfig == null || runConfig.rows() == 0) return true; int run = runConfig.getInt("run", 0); - + // Set swimmer. Swim swimmer = new Swim(); @@ -86,107 +85,53 @@ public boolean processDataEvent(DataEvent event) { // uncomment this code. IndexedTable beamOffset = this.getConstantsManager().getConstants(run, "/geometry/beam/position"); double xB = beamOffset.getDoubleValue("x_offset", 0,0,0); - double yB = beamOffset.getDoubleValue("y_offset", 0,0,0); - - // get status table - IndexedTable status = this.getConstantsManager().getConstants(run, "/calibration/mvt/fmt_status"); - - // get time table - IndexedTable timecuts = this.getConstantsManager().getConstants(run, "/calibration/mvt/fmt_time"); - - // === HITS ================================================================================ - List hits = Hit.fetchHits(event, timecuts, status); - if (hits.isEmpty()) return true; + double yB = beamOffset.getDoubleValue("y_offset", 0,0,0); + // === CLUSTERS ============================================================================ - List clusters = Cluster.findClusters(hits); - if(debug) for (int i = 0; i < clusters.size(); i++) System.out.println(clusters.get(i).toString()); + URWellReader uRWellReader = new URWellReader(event); + List clusters = uRWellReader.getUrwellClusters(); + //System.out.println(clusters.size()); + //for (int i = 0; i < clusters.size(); i++) System.out.println(clusters.get(i).toString()); // === DC TRACKS =========================================================================== - List tracks = Track.getDCTracks(event, swimmer); + Track trk = new Track(); + List tracks = trk.getDCTracks(event, swimmer); if(tracks.isEmpty()) return true; - + // === SEEDS ============================================================================= for(int i=0; iother.getSeedQuality()) { - other.clearClusters(layer); - cluster.setTrackIndex(track.getIndex()); - cluster.setDoca(doca); - } - } - } - } + + List filtedTracks = new ArrayList(); + for(Track track : tracks){ + if (track.getClusters().size() > 4) filtedTracks.add(track); } - + + // === TRACKS ============================================================================== KFitter kf = null; // Iterate on list to run the fit. - for(int i=0; i crs = crossesMap.get(id); -// for(Cross c : crs) { -// // Filter clusters to use only the best cluster (minimum Tmin) per FMT layer. -// int lyr = c.getCluster1().getLayer(); -//// System.out.println(c.getCluster1().getDoca()); -// if (cls.get(lyr) == null || c.getCluster1().getDoca()< cls.get(lyr).getDoca()) -// cls.put(lyr, c.getCluster1()); -// } - // Set status and stop if there are not at least two measurements to fit against. - List trackClusters = track.getClusters(); - - if (trackClusters.size()<2) continue; + List trackClusters = track.getClusters(); - kf = new KFitter(track, swimmer, 0); + kf = new KFitter(track, swimmer, 0); + kf.runFitter(track.getSector()); + // Do one last KF pass with filtering off to get the final Chi^2. kf.totNumIter = 1; kf.filterOn = false; @@ -197,11 +142,15 @@ public boolean processDataEvent(DataEvent event) { // swim to beamline to get vertex parameters int charge = (int)Math.signum(sv.Q); - swimmer.SetSwimParameters(sv.x,sv.y,sv.z,-sv.getPx(),-sv.getPy(),-sv.getPz(),-charge); + + Point3D posGlobal = track.transLocaltoGlobal(track.getSector(), sv.x, sv.y, sv.z); + Point3D momGlobal = track.transLocaltoGlobal(track.getSector(), sv.getPx(), sv.getPy(), sv.getPz()); + + swimmer.SetSwimParameters(posGlobal.x(),posGlobal.y(),posGlobal.z(), -momGlobal.x(),-momGlobal.y(),-momGlobal.z(),-charge); double[] Vt = swimmer.SwimToBeamLine(xB, yB); // if successful, save track parameters - if(Vt == null || Vt[6] urClusters) { + URWellCluster cluster = null; + + for (URWellCluster cl : urClusters) { + if (cl.id() == cluster1) { + cluster = cl; + break; + } + } + + cls1 = cluster; + } + + public void setCluster2(List urClusters) { + URWellCluster cluster = null; + + for (URWellCluster cl : urClusters) { + if (cl.id() == cluster2) { + cluster = cl; + break; + } + } + + cls2 = cluster; + } + + public URWellCluster getCluster1(List urClusters) { + URWellCluster cluster = null; + + for (URWellCluster cl : urClusters) { + if (cl.id() == cluster1) { + cluster = cl; + break; + } + } + + return cluster; + } + + public URWellCluster getCluster2(List urClusters) { + URWellCluster cluster = null; + + for (URWellCluster cl : urClusters) { + if (cl.id() == cluster2) { + cluster = cl; + break; + } + } + + return cluster; + } +} diff --git a/reconstruction/urwell/src/main/java/org/jlab/rec/urwell/reader/URWellHit.java b/reconstruction/urwell/src/main/java/org/jlab/rec/urwell/reader/URWellHit.java new file mode 100644 index 0000000000..33fa88e9f4 --- /dev/null +++ b/reconstruction/urwell/src/main/java/org/jlab/rec/urwell/reader/URWellHit.java @@ -0,0 +1,43 @@ +package org.jlab.rec.urwell.reader; + +import org.jlab.detector.base.DetectorDescriptor; +import org.jlab.detector.base.DetectorType; + +/** + * + * @author Tongtong Cao + */ + +public class URWellHit { + + private DetectorDescriptor desc = new DetectorDescriptor(DetectorType.URWELL); + private double energy = 0; + private double time = 0; + + public URWellHit(int sector, int layer, int component, double energy, double time) { + this.desc.setSectorLayerComponent(sector, layer, component); + this.energy = energy; + this.time = time; + } + + public int sector() { + return this.desc.getSector(); + } + + public int layer() { + return this.desc.getLayer(); + } + + public int strip() { + return this.desc.getComponent(); + } + + public double energy() { + return energy; + } + + public double time() { + return time; + } + +} diff --git a/reconstruction/urwell/src/main/java/org/jlab/rec/urwell/reader/URWellReader.java b/reconstruction/urwell/src/main/java/org/jlab/rec/urwell/reader/URWellReader.java new file mode 100644 index 0000000000..2723cd7032 --- /dev/null +++ b/reconstruction/urwell/src/main/java/org/jlab/rec/urwell/reader/URWellReader.java @@ -0,0 +1,144 @@ +package org.jlab.rec.urwell.reader; + +import java.util.ArrayList; +import java.util.List; +import org.jlab.geom.prim.Point3D; +import org.jlab.io.base.DataBank; +import org.jlab.io.base.DataEvent; + +/** + * + * @author Tongtong Cao + */ +public class URWellReader{ + + private final List urHits = new ArrayList<>(); + private final List urClusters = new ArrayList<>(); + private final List urCrosses = new ArrayList<>(); + + + public URWellReader(DataEvent event) { + + if(event.hasBank("URWELL::hits")) + this.readHits(event.getBank("URWELL::hits")); + if(event.hasBank("URWELL::clusters")) + this.readClusters(event.getBank("URWELL::clusters")); + if(event.hasBank("URWELL::crosses")) + this.readCrosses(event.getBank("URWELL::crosses")); + } + + public List getUrwellHits() { + return urHits; + } + + public List getUrwellClusters() { + return urClusters; + } + + public List getUrwellCrosses() { + return urCrosses; + } + + public List getUrwellR1Crosses() { + List urCrossesR1 = new ArrayList<>(); + for (URWellCross crs : urCrosses) { + if (crs.region() == 1) { + urCrossesR1.add(crs); + } + } + return urCrossesR1; + } + + public final void readHits(DataBank bank) { + + for(int i=0; i Date: Wed, 16 Jul 2025 18:07:45 -0400 Subject: [PATCH 14/16] updated FT reconstruction for ddvcs calorimeter, added 20 MeV threshold for FTCAL bg merging --- .../analysis/eventmerger/ADCTDCMerger.java | 3 +- .../org/jlab/rec/ft/cal/FTCALCluster.java | 4 +-- .../java/org/jlab/rec/ft/cal/FTCALHit.java | 33 ++++++++++--------- .../jlab/rec/ft/cal/FTCALReconstruction.java | 2 +- 4 files changed, 22 insertions(+), 20 deletions(-) diff --git a/common-tools/clas-analysis/src/main/java/org/jlab/analysis/eventmerger/ADCTDCMerger.java b/common-tools/clas-analysis/src/main/java/org/jlab/analysis/eventmerger/ADCTDCMerger.java index 290bf12cfb..11788285c3 100644 --- a/common-tools/clas-analysis/src/main/java/org/jlab/analysis/eventmerger/ADCTDCMerger.java +++ b/common-tools/clas-analysis/src/main/java/org/jlab/analysis/eventmerger/ADCTDCMerger.java @@ -86,7 +86,8 @@ public List readADCs(DetectorType detector,DataBank bank) { ADC adcData = new ADC(detector); adcData.readFromBank(bank, i); - if(!adcData.isGood()) adcData.skip(); + if(!adcData.isGood() || + (detector==DetectorType.FTCAL && adcData.getAdc()<200)) adcData.skip(); adcStore.add(adcData); } diff --git a/reconstruction/ft/src/main/java/org/jlab/rec/ft/cal/FTCALCluster.java b/reconstruction/ft/src/main/java/org/jlab/rec/ft/cal/FTCALCluster.java index b7e09cbd53..dacbad6cd2 100644 --- a/reconstruction/ft/src/main/java/org/jlab/rec/ft/cal/FTCALCluster.java +++ b/reconstruction/ft/src/main/java/org/jlab/rec/ft/cal/FTCALCluster.java @@ -47,7 +47,7 @@ public double getEnergy() { public double getFullEnergy(IndexedTable energyTable) { // return energy corrected for leakage and threshold effects double clusterEnergy = this.getEnergy(); - int seedID = this.get(0).get_COMPONENT(); + int seedID = FTCALHit.REFCOMPONENT;//this.get(0).get_COMPONENT(); double energyCorr = (energyTable.getDoubleValue("c0",1,1,seedID) + energyTable.getDoubleValue("c1",1,1,seedID)*clusterEnergy + energyTable.getDoubleValue("c2",1,1,seedID)*clusterEnergy*clusterEnergy @@ -197,7 +197,7 @@ private double weight(FTCALHit hit, double clusEnergy) { public boolean containsHit(FTCALHit hit, IndexedTable thresholds, IndexedTable clusterTable) { boolean addFlag = false; - if(hit.get_Edep()>thresholds.getDoubleValue("thresholdCluster",1,1,hit.get_COMPONENT())) { + if(hit.get_Edep()>thresholds.getDoubleValue("thresholdCluster",1,1,FTCALHit.REFCOMPONENT)) { for(int j = 0; j< this.size(); j++) { double tDiff = Math.abs(hit.get_Time() - this.get(j).get_Time()); double xDiff = Math.abs(hit.get_IDX() - this.get(j).get_IDX()); diff --git a/reconstruction/ft/src/main/java/org/jlab/rec/ft/cal/FTCALHit.java b/reconstruction/ft/src/main/java/org/jlab/rec/ft/cal/FTCALHit.java index 507bdf089c..de425fe10c 100644 --- a/reconstruction/ft/src/main/java/org/jlab/rec/ft/cal/FTCALHit.java +++ b/reconstruction/ft/src/main/java/org/jlab/rec/ft/cal/FTCALHit.java @@ -6,25 +6,26 @@ public class FTCALHit implements Comparable{ // class implements Comparable interface to allow for sorting a collection of hits by Edep values - + public static final int REFCOMPONENT =245; + // constructor public FTCALHit(int i, int ICOMPONENT, int ADC, int TDC, IndexedTable charge2Energy, IndexedTable timeOffsets, IndexedTable timeWalk, IndexedTable cluster) { this._COMPONENT = ICOMPONENT; - this._IDY = ((int) ICOMPONENT/22) + 1; - this._IDX = ICOMPONENT + 1 - (this._IDY-1)*22; + this._IDY = ((int) ICOMPONENT/44) + 1; + this._IDX = ICOMPONENT + 1 - (this._IDY-1)*44; this._ADC = ADC; this._TDC = TDC; - this._Charge = ((double) this._ADC)*charge2Energy.getDoubleValue("fadc_to_charge", 1,1,ICOMPONENT); - this.set_Edep(this._Charge*charge2Energy.getDoubleValue("mips_energy", 1,1,ICOMPONENT) - /charge2Energy.getDoubleValue("mips_charge", 1,1,ICOMPONENT)/1000.); + this._Charge = ((double) this._ADC)*charge2Energy.getDoubleValue("fadc_to_charge", 1,1,REFCOMPONENT); + this.set_Edep(this._Charge*charge2Energy.getDoubleValue("mips_energy", 1,1,REFCOMPONENT) + /charge2Energy.getDoubleValue("mips_charge", 1,1,REFCOMPONENT)/1000.); double twCorr=0; if(this._Charge>0) { - twCorr = timeWalk.getDoubleValue("amplitude", 1,1,ICOMPONENT)*Math.exp(-this._Charge*timeWalk.getDoubleValue("lambda", 1,1,ICOMPONENT)); + twCorr = timeWalk.getDoubleValue("amplitude", 1,1,REFCOMPONENT)*Math.exp(-this._Charge*timeWalk.getDoubleValue("lambda", 1,1,REFCOMPONENT)); } this.set_Time(((double) this._TDC)/FTCALConstantsLoader.TIMECONVFAC -(FTCALConstantsLoader.CRYS_LENGTH-cluster.getDoubleValue("depth_z", 1,1,0))/FTCALConstantsLoader.VEFF - -timeOffsets.getDoubleValue("time_offset", 1,1,ICOMPONENT)-twCorr); + -timeOffsets.getDoubleValue("time_offset", 1,1,REFCOMPONENT)-twCorr); // if(this.get_Edep()>0.1) System.out.println(ICOMPONENT + " " + this._TDC + " " + // FTCALConstantsLoader.TIMECONVFAC + " " + FTCALConstantsLoader.time_offset[0][0][ICOMPONENT-1] + " " + // this.get_Time()); @@ -37,21 +38,21 @@ public FTCALHit(int i, int ICOMPONENT, int ADC, int TDC, IndexedTable charge2Ene public FTCALHit(int i, int ICOMPONENT, int ADC, float time, IndexedTable charge2Energy, IndexedTable timeOffsets, IndexedTable timeWalk, IndexedTable cluster) { this._COMPONENT = ICOMPONENT; - this._IDY = ((int) ICOMPONENT/22) + 1; - this._IDX = ICOMPONENT + 1 - (this._IDY-1)*22; + this._IDY = ((int) ICOMPONENT/44) + 1; + this._IDX = ICOMPONENT + 1 - (this._IDY-1)*44; this._ADC = ADC; - this._Charge = ((double) this._ADC)*charge2Energy.getDoubleValue("fadc_to_charge", 1,1,ICOMPONENT); - this.set_Edep(this._Charge*charge2Energy.getDoubleValue("mips_energy", 1,1,ICOMPONENT) - /charge2Energy.getDoubleValue("mips_charge", 1,1,ICOMPONENT)/1000.); + this._Charge = ((double) this._ADC)*charge2Energy.getDoubleValue("fadc_to_charge", 1,1,REFCOMPONENT); + this.set_Edep(this._Charge*charge2Energy.getDoubleValue("mips_energy", 1,1,REFCOMPONENT) + /charge2Energy.getDoubleValue("mips_charge", 1,1,REFCOMPONENT)/1000.); double twCorr=0; if(this._Charge>0) { - twCorr = timeWalk.getDoubleValue("amplitude", 1,1,ICOMPONENT)*Math.exp(-this._Charge*timeWalk.getDoubleValue("lambda", 1,1,ICOMPONENT)); + twCorr = timeWalk.getDoubleValue("amplitude", 1,1,REFCOMPONENT)*Math.exp(-this._Charge*timeWalk.getDoubleValue("lambda", 1,1,REFCOMPONENT)); } this.set_Time(time -(FTCALConstantsLoader.CRYS_LENGTH-cluster.getDoubleValue("depth_z", 1,1,0))/FTCALConstantsLoader.VEFF - -timeOffsets.getDoubleValue("time_offset", 1,1,ICOMPONENT)-twCorr); + -timeOffsets.getDoubleValue("time_offset", 1,1,REFCOMPONENT)-twCorr); // if(this.get_Edep()>0.1) System.out.println(ICOMPONENT + " " + this._TDC + " " + // FTCALConstantsLoader.TIMECONVFAC + " " + FTCALConstantsLoader.time_offset[0][0][ICOMPONENT-1] + " " + // this.get_Time()); @@ -208,7 +209,7 @@ public final void set_ClusID(int _ClusIndex) { public static boolean passHitSelection(FTCALHit hit, IndexedTable thresholds) { // a selection cut to pass the hit. - if(hit.get_Edep() > thresholds.getDoubleValue("thresholdHit", 1,1,hit.get_COMPONENT())) { + if(hit.get_Edep() > thresholds.getDoubleValue("thresholdHit", 1,1,REFCOMPONENT)) { return true; } else { return false; diff --git a/reconstruction/ft/src/main/java/org/jlab/rec/ft/cal/FTCALReconstruction.java b/reconstruction/ft/src/main/java/org/jlab/rec/ft/cal/FTCALReconstruction.java index fd9f113a1d..c4b4ad4f5e 100644 --- a/reconstruction/ft/src/main/java/org/jlab/rec/ft/cal/FTCALReconstruction.java +++ b/reconstruction/ft/src/main/java/org/jlab/rec/ft/cal/FTCALReconstruction.java @@ -216,7 +216,7 @@ public List readRawHitsHipo(DataEvent event, IndexedTable charge2Energ int adc = bankDGTZ.getInt("ADC",row); float time = bankDGTZ.getFloat("time",row); if(ilayer==0) ilayer=1; // fix for wrong layer in TT - if(adc!=-1 && time!=-1 && status.getIntValue("status", isector, ilayer, icomponent)==0){ + if(adc!=-1 && time!=-1 && status.getIntValue("status", isector, ilayer, FTCALHit.REFCOMPONENT)==0){ FTCALHit hit = new FTCALHit(bankDGTZ.trueIndex(row),icomponent, adc, time, charge2Energy, timeOffsets, timeWalk, cluster); //////////////////////////////////////////////////////////////////////////////////////////////////// From c8058bd08bb7c6d74b5051eb8daa8a6dfbd57b9e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mariangela=20Bond=C3=AD?= Date: Thu, 17 Jul 2025 17:46:14 -0400 Subject: [PATCH 15/16] implemented layer surface info for the FVT --- .../geant4/v2/URWELL/URWellStripFactory.java | 92 ++++++++++++++++++- 1 file changed, 89 insertions(+), 3 deletions(-) diff --git a/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellStripFactory.java b/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellStripFactory.java index f78c8c40b2..8cf110c400 100644 --- a/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellStripFactory.java +++ b/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellStripFactory.java @@ -11,6 +11,7 @@ import org.jlab.geom.prim.Point3D; import org.jlab.geom.prim.Vector3D; import org.jlab.geometry.prim.Line3d; +import org.jlab.geom.prim.Trap3D; import org.jlab.geometry.prim.Straight; import org.jlab.utils.groups.IndexedList; @@ -25,6 +26,9 @@ public final class URWellStripFactory { private IndexedList globalStrips = new IndexedList(3); private IndexedList localStrips = new IndexedList(3); private IndexedList planeStrips = new IndexedList(3); + private IndexedList surfaceLayers = new IndexedList(2); + + private int nRegions; private int nSectors; private int nChambers; @@ -107,6 +111,7 @@ public void init(DatabaseConstantProvider cp, boolean prototype, int regions) { } this.fillStripLists(); this.fillPlaneLists(); + this.fillSurfaceLists(); } /** @@ -156,6 +161,7 @@ public void init(DatabaseConstantProvider cp, int config, int regions) { this.fillStripLists(); this.fillPlaneLists(); + this.fillSurfaceLists(); } public double getStereoAngle(int layer, int ichamber) { @@ -535,9 +541,86 @@ private Plane3D createPLane(int sector, int layer, int strip) { } Plane3D plane = new Plane3D(First_strip.origin(), normal_plane); - + return plane; } + + + private void fillSurfaceLists() { + + for (int ir = 0; ir < nRegions; ir++) { + int region = ir + 1; + + for (int is = 0; is < nSectors; is++) { + int sector = is + 1; + if (isProto == true) { + sector = 6; + } + + for (int il = 0; il < nLayers; il++) { + + int layer = (2 * region - 1) + il; + + Trap3D surf = this.createSurface(sector, layer); + this.surfaceLayers.add(surf, sector, layer); + + + } + } + } + } + + /** + * Define the boundary points of the uRWELL surface using the first strip + * from the odd and even layers. This method assumes the strips are aligned + * with the trapezoid sides, and is only valid under that geometric + * condition. + * + * @param sector + * @param layer + * @return + */ + private Trap3D createSurface(int sector, int layer) { + + int layer2 = 0; + if (layer % 2 == 0) { + layer2 = layer - 1; + } else { + layer2 = layer + 1; + } + + Line3D first_strip = this.getStrip(sector, layer, 1); + Line3D second_strip = this.getStrip(sector, layer2, 1); + + Point3D p0 = first_strip.origin(); + Point3D p1 = first_strip.end(); + Point3D p2 = second_strip.origin(); + Point3D p3 = second_strip.end(); + + Trap3D trapezoid = new Trap3D(p0.x(), p0.y(), p0.z(), p1.x(), p1.y(), p1.z(), p2.x(), p2.y(), p2.z(), p3.x(), p3.y(), p3.z()); + + return trapezoid; + } + + /** + * + * @param sector + * @param layer + * @return + */ + public Trap3D getSurface(int sector, int layer) { + + return surfaceLayers.getItem(sector, layer); + } + + + + + + + + + public static void main(String[] args) { DatabaseConstantProvider cp = new DatabaseConstantProvider(11, "default"); @@ -546,7 +629,7 @@ public static void main(String[] args) { // URWellGeant4Factory factory2 = new URWellGeant4Factory(cp, true, 2); // URWellStripFactory factory2 = new URWellStripFactory(cp, false, 1); - URWellStripFactory factory2 = new URWellStripFactory(cp, 0, 2); + URWellStripFactory factory2 = new URWellStripFactory(cp, 2, 6); double angle = factory2.getStereoAngle(0,0); @@ -559,12 +642,15 @@ public static void main(String[] args) { System.out.println("Numero di strip "+factory2.getNStripSector(layer)); System.out.println("Numero di strip chamber "+factory2.getNStripChamber(layer,0)); + + Trap3D surf = factory2.getSurface(5,1); + System.out.println("norm " +surf.normal()); // Plane3D plane = factory2.getPlane(6, 1, 200); // System.out.println(plane.toString()); - System.out.println((strip) + " " + factory2.getLocalStripId(layer,strip) + "\n" + factory2.getChamberStrip(sector, layer, strip) + "\n" + factory2.getLocalStrip(sector, layer, strip) + "\n" + factory2.getStrip(sector, layer, strip)) ; +// System.out.println((strip) + " " + factory2.getLocalStripId(layer,strip) + "\n" + factory2.getChamberStrip(sector, layer, strip) + "\n" + factory2.getLocalStrip(sector, layer, strip) + "\n" + factory2.getStrip(sector, layer, strip)) ; /* for(int istrip=0; istrip Date: Fri, 15 Aug 2025 14:26:56 -0400 Subject: [PATCH 16/16] ddvcs fvt: implemented seeding based on crosses --- .../main/java/org/jlab/rec/dc/Constants.java | 6 +- .../rec/dc/trajectory/TrajectorySurfaces.java | 9 +- .../main/java/org/jlab/rec/fmt/Constants.java | 5 +- .../java/org/jlab/rec/fmt/track/Track.java | 97 +++++++++++++++++++ .../org/jlab/rec/fmt/track/Trajectory.java | 2 +- .../java/org/jlab/service/fmt/FMTEngine.java | 70 +++++++++++-- .../jlab/rec/urwell/reader/URWellCluster.java | 6 ++ 7 files changed, 180 insertions(+), 15 deletions(-) diff --git a/reconstruction/dc/src/main/java/org/jlab/rec/dc/Constants.java b/reconstruction/dc/src/main/java/org/jlab/rec/dc/Constants.java index 25473ca915..d66f313a5e 100644 --- a/reconstruction/dc/src/main/java/org/jlab/rec/dc/Constants.java +++ b/reconstruction/dc/src/main/java/org/jlab/rec/dc/Constants.java @@ -14,6 +14,8 @@ import org.jlab.detector.geom.RICH.RICHGeoFactory; import org.jlab.geom.base.ConstantProvider; import org.jlab.detector.calib.utils.ConstantsManager; +import org.jlab.detector.calib.utils.DatabaseConstantProvider; +import org.jlab.detector.geant4.v2.URWELL.URWellStripFactory; import org.jlab.geom.base.Detector; import org.jlab.rec.dc.trajectory.TrajectorySurfaces; import org.jlab.utils.groups.IndexedTable; @@ -122,7 +124,7 @@ public static Constants getInstance() { public DCGeant4Factory dcDetector = null; public FTOFGeant4Factory ftofDetector = null; public Detector ecalDetector = null; - public Detector fmtDetector = null; + public URWellStripFactory fmtDetector = null; public RICHGeoFactory richDetector = null; public TrajectorySurfaces trajSurfaces = null; @@ -546,7 +548,7 @@ private synchronized void LoadGeometry(String geoVariation, double[][] shifts) { ConstantProvider providerFTOF = GeometryFactory.getConstants(DetectorType.FTOF, 11, geoVariation); ftofDetector = new FTOFGeant4Factory(providerFTOF); ecalDetector = GeometryFactory.getDetector(DetectorType.ECAL, 11, geoVariation); - fmtDetector = GeometryFactory.getDetector(DetectorType.FMT, 11, geoVariation); + fmtDetector = new URWellStripFactory(new DatabaseConstantProvider(),2, 6); ConstantsManager managerRICH = new ConstantsManager(geoVariation);; richDetector = new RICHGeoFactory(0, managerRICH, 11, false); // create the surfaces diff --git a/reconstruction/dc/src/main/java/org/jlab/rec/dc/trajectory/TrajectorySurfaces.java b/reconstruction/dc/src/main/java/org/jlab/rec/dc/trajectory/TrajectorySurfaces.java index 06241218dc..5d6b9ceb8e 100644 --- a/reconstruction/dc/src/main/java/org/jlab/rec/dc/trajectory/TrajectorySurfaces.java +++ b/reconstruction/dc/src/main/java/org/jlab/rec/dc/trajectory/TrajectorySurfaces.java @@ -17,6 +17,7 @@ import java.io.PrintWriter; import java.util.logging.Logger; +import org.jlab.detector.geant4.v2.URWELL.URWellStripFactory; import org.jlab.geom.detector.ec.ECLayer; import org.jlab.geom.detector.ec.ECSuperlayer; import org.jlab.geom.detector.fmt.FMTLayer; @@ -46,7 +47,7 @@ public void setDetectorPlanes(List> planes) { } public void loadSurface(double targetPosition, double targetLength, DCGeant4Factory dcDetector, - FTOFGeant4Factory ftofDetector, Detector ecalDetector, Detector fmtDetector, RICHGeoFactory richDetector) { + FTOFGeant4Factory ftofDetector, Detector ecalDetector, URWellStripFactory fmtDetector, RICHGeoFactory richDetector) { // creating Boundaries for MS Constants.getInstance().Z[0]= targetPosition; Constants.getInstance().Z[1]= dcDetector.getWireMidpoint(0, 0, 0, 0).z; @@ -76,9 +77,9 @@ public void loadSurface(double targetPosition, double targetLength, DCGeant4Fact this.detectorPlanes.get(isector).add(new Surface(DetectorType.TARGET, sector, DetectorLayer.TARGET_CENTER, new Plane3D(0, 0, targetPosition, 0, 0, 1))); // Add FMT layers - for (int ilayer=0; ilayer<6; ++ilayer) { - FMTLayer fmtLayer = (FMTLayer) fmtDetector.getSector(0).getSuperlayer(0).getLayer(ilayer); - this.detectorPlanes.get(isector).add(new Surface(DetectorType.FMT, sector, ilayer+1, fmtLayer.getTrajectorySurface(), 0)); + for (int ilayer=0; ilayer<12; ++ilayer) { +// FMTLayer fmtLayer = (FMTLayer) fmtDetector.getSector(0).getSuperlayer(0).getLayer(ilayer); + this.detectorPlanes.get(isector).add(new Surface(DetectorType.FMT, sector, ilayer+1, fmtDetector.getSurface(sector, ilayer+1), 0)); } // Add DC diff --git a/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/Constants.java b/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/Constants.java index 08969f9af1..b05a2a809b 100644 --- a/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/Constants.java +++ b/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/Constants.java @@ -15,8 +15,11 @@ public class Constants { public static final int NLAYERS = 12; + // rMax difference between cross's clusters energy + public static double CROSSDELTAE = 200; + // DC-tracks to FMT-clusters matching parameter - public static double CIRCLECONFUSION = 1.2; // cm + public static double CIRCLECONFUSION = 12; // cm // min path for final swimming to beamline to reject failed swimming public static double MIN_SWIM_PATH = 0.2; diff --git a/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/Track.java b/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/Track.java index 6167ba64f7..56e662cd67 100644 --- a/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/Track.java +++ b/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/Track.java @@ -6,6 +6,7 @@ import java.util.Map; import java.util.Map.Entry; import org.jlab.clas.swimtools.Swim; +import org.jlab.detector.base.DetectorType; import org.jlab.geom.prim.Line3D; import org.jlab.geom.prim.Point3D; import org.jlab.geom.prim.Vector3D; @@ -46,6 +47,7 @@ public class Track { private double _pz; private int _NDF; + private final Trajectory[] _DCtrajs = new Trajectory[Constants.NLAYERS]; private final List[] _clusters = new ArrayList[Constants.NLAYERS]; private final Trajectory[] _FMTtrajs = new Trajectory[Constants.NLAYERS]; @@ -68,6 +70,23 @@ public Track(int _index, int _sector, int _q, double _x, double _y, double _z, } + + /** + * @param layer + * @return the _traj + */ + public Trajectory getDCTraj(int layer) { + if(layer<=0 || layer>Constants.NLAYERS) return null; + else return _DCtrajs[layer-1]; + } + + /** + * @param trj + */ + public void setDCTraj(Trajectory trj) { + this._DCtrajs[trj.getLayer()-1] = trj; + } + public List getClusters() { List clusters = new ArrayList<>(); for(int i=0; i0 && this.getClusterLayer(2)>0 && this.getClusterLayer(3)>0) { +// Line3D seg1 = this.getClusters(1).get(0).getGlobalSegment(); +// Line3D seg2 = this.getClusters(2).get(0).getGlobalSegment(); +// Line3D seg3 = this.getClusters(3).get(0).getGlobalSegment(); +// quality = seg1.distanceSegments(seg3).distanceSegments(seg2).length(); +// } +// return quality; +// } + + +// // FIXME: THIS METHOD SHOULD BE GENERALIZED +// public void filterClusters(int mode) { +// +// double dref = Double.POSITIVE_INFINITY; +// +// int[] ibest = new int[Constants.NLAYERS]; +// for(int i=0; i0 && this.getClusterLayer(2)>0 && this.getClusterLayer(3)>0) { +// for(int i1=0; i1=0) { +// List cls = new ArrayList<>(); +// cls.add(_clusters[i].get(ibest[i])); +// _clusters[i] = cls; +// } +// } +// } + public List getDCTracks(DataEvent event, Swim swimmer) { Map trackmap = new LinkedHashMap(); @@ -330,6 +412,21 @@ public List getDCTracks(DataEvent event, Swim swimmer) { trk.setStatus(1); trackmap.put(id,trk); } + for (int i = 0; i < trajBank.rows(); i++) { + if (trajBank.getByte("detector", i) == DetectorType.FMT.getDetectorId()) { + int id = trajBank.getShort("id", i); + int layer = trajBank.getByte("layer", i); + Trajectory trj = new Trajectory(layer, + trajBank.getFloat("x", i), + trajBank.getFloat("y", i), + trajBank.getFloat("z", i), + trajBank.getFloat("tx", i), + trajBank.getFloat("ty", i), + trajBank.getFloat("tz", i), + trajBank.getFloat("path", i)); + trackmap.get(id).setDCTraj(trj); + } + } } List tracks = new ArrayList<>(); for(Entry entry: trackmap.entrySet()) { diff --git a/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/Trajectory.java b/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/Trajectory.java index 233bc6692e..d5ea451410 100644 --- a/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/Trajectory.java +++ b/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/Trajectory.java @@ -26,7 +26,7 @@ public Trajectory(int layer, Point3D p, Vector3D d, double path) { public Trajectory(int layer, double x, double y, double z, double tx, double ty, double tz, double path) { this.layer = layer; this.position = new Point3D(x,y,z); - this.localPosition = Constants.toLocal(layer, position); +// this.localPosition = Constants.toLocal(layer, position); this.direction = new Vector3D(tx,ty,tz); this.path = path; } diff --git a/reconstruction/fmt/src/main/java/org/jlab/service/fmt/FMTEngine.java b/reconstruction/fmt/src/main/java/org/jlab/service/fmt/FMTEngine.java index 9b19bf2344..0c41257436 100644 --- a/reconstruction/fmt/src/main/java/org/jlab/service/fmt/FMTEngine.java +++ b/reconstruction/fmt/src/main/java/org/jlab/service/fmt/FMTEngine.java @@ -3,12 +3,15 @@ import java.util.Arrays; import java.util.List; import java.util.ArrayList; +import java.util.Comparator; import org.jlab.clas.reco.ReconstructionEngine; import org.jlab.clas.swimtools.Swim; import org.jlab.detector.base.DetectorType; import org.jlab.detector.base.GeometryFactory; +import org.jlab.geom.prim.Line3D; import org.jlab.geom.prim.Point3D; +import org.jlab.geom.prim.Vector3D; import org.jlab.utils.groups.IndexedTable; import org.jlab.io.base.*; import org.jlab.rec.fmt.Constants; @@ -20,6 +23,7 @@ import org.jlab.rec.fmt.track.fit.StateVecs.StateVec; import org.jlab.rec.urwell.reader.URWellCluster; +import org.jlab.rec.urwell.reader.URWellCross; import org.jlab.rec.urwell.reader.URWellReader; /** @@ -91,22 +95,74 @@ public boolean processDataEvent(DataEvent event) { // === CLUSTERS ============================================================================ URWellReader uRWellReader = new URWellReader(event); List clusters = uRWellReader.getUrwellClusters(); + List crosses = uRWellReader.getUrwellCrosses(); //System.out.println(clusters.size()); //for (int i = 0; i < clusters.size(); i++) System.out.println(clusters.get(i).toString()); - + // === DC TRACKS =========================================================================== Track trk = new Track(); List tracks = trk.getDCTracks(event, swimmer); if(tracks.isEmpty()) return true; // === SEEDS ============================================================================= + Point3D target = new Point3D(0,0,0); for(int i=0; i[] trackCrosses = new ArrayList[Constants.NLAYERS/2]; + for(int j=0; j(); + for(int j=0; j> segments = new ArrayList<>(); + List end = new ArrayList<>(); + for(int r=6; r>3; r--) + end.addAll(trackCrosses[r-1]); + for(URWellCross ce : end) { + Line3D ve = new Line3D(ce.position(),target); + for(int ro=1; ro segment = new ArrayList(); + segment.add(ce); + for(int r=co.region()+1; r=4) + segments.add(segment); + } + } + } + } + System.out.print(segments.size()); + if(!segments.isEmpty()) { + segments.sort(Comparator.comparingInt(List::size).reversed()); + for(URWellCross cross : segments.get(0)) { + track.addCluster(cross.getCluster1()); + track.addCluster(cross.getCluster2()); } } } @@ -144,7 +200,7 @@ public boolean processDataEvent(DataEvent event) { int charge = (int)Math.signum(sv.Q); Point3D posGlobal = track.transLocaltoGlobal(track.getSector(), sv.x, sv.y, sv.z); - Point3D momGlobal = track.transLocaltoGlobal(track.getSector(), sv.getPx(), sv.getPy(), sv.getPz()); + Point3D momGlobal = track.transLocaltoGlobal(track.getSector(), sv.getPx()/sv.getP()*track.getP(), sv.getPy()/sv.getP()*track.getP(), sv.getPz()/sv.getP()*track.getP()); swimmer.SetSwimParameters(posGlobal.x(),posGlobal.y(),posGlobal.z(), -momGlobal.x(),-momGlobal.y(),-momGlobal.z(),-charge); double[] Vt = swimmer.SwimToBeamLine(xB, yB); diff --git a/reconstruction/urwell/src/main/java/org/jlab/rec/urwell/reader/URWellCluster.java b/reconstruction/urwell/src/main/java/org/jlab/rec/urwell/reader/URWellCluster.java index 7a74e514db..bbba97d199 100644 --- a/reconstruction/urwell/src/main/java/org/jlab/rec/urwell/reader/URWellCluster.java +++ b/reconstruction/urwell/src/main/java/org/jlab/rec/urwell/reader/URWellCluster.java @@ -20,6 +20,7 @@ public class URWellCluster { private double time = 0; private int crossIndex = -1; private double stereo = 10.; + private Line3D line = new Line3D(); private Line3D lineLocal = new Line3D(); public URWellCluster(int id, int sector, int layer, int component, int size, double energy, double time, Point3D pointOrigin, Point3D pointEnd) { @@ -28,6 +29,7 @@ public URWellCluster(int id, int sector, int layer, int component, int size, dou this.size = size; this.energy = energy; this.time = time; + this.line.set(pointOrigin, pointEnd); Point3D pointOriginLocal = new Point3D(); pointOriginLocal.copy(pointOrigin); @@ -84,6 +86,10 @@ public Line3D getLineLocal(){ return lineLocal; } + public Line3D getLine(){ + return line; + } + /** * * @return cluster info. about location and number of hits contained in it