diff --git a/etc/bankdefs/hipo4/alert.json b/etc/bankdefs/hipo4/alert.json index 7c125eb9c..ea6f4e251 100644 --- a/etc/bankdefs/hipo4/alert.json +++ b/etc/bankdefs/hipo4/alert.json @@ -1,5 +1,105 @@ [ { + "name": "AHDC::Projections", + "group": 23000, + "item": 31, + "info": "Track Projections to ATOF", + "entries": [ + { + "name": "x_at_bar", + "type": "F", + "info": "x position at atof bar (middle surface) in mm" + }, { + "name": "y_at_bar", + "type": "F", + "info": "y position at atof bar (middle surface) in mm" + }, { + "name": "z_at_bar", + "type": "F", + "info": "z position at atof bar (middle surface) in mm" + },{ + "name": "L_at_bar", + "type": "F", + "info": "path length at atof bar (inner surface) in mm" + },{ + "name": "L_in_bar", + "type": "F", + "info": "path length inside atof bar in mm" + },{ + "name": "x_at_wedge", + "type": "F", + "info": "x position at atof wedge (middle surface) in mm" + }, { + "name": "y_at_wedge", + "type": "F", + "info": "y position at atof wedge (middle surface) in mm" + }, { + "name": "z_at_wedge", + "type": "F", + "info": "z position at atof wedge (middle surface) in mm" + },{ + "name": "L_at_wedge", + "type": "F", + "info": "path length at atof wedge (inner surface) in mm" + },{ + "name": "L_in_wedge", + "type": "F", + "info": "path length inside atof wedge in mm" + } + ] + },{ + "name": "ATOF::hits", + "group": 22500, + "item": 21, + "info": "Hits in ATOF", + "entries": [ + { + "name": "id", + "type": "S", + "info": "hit id" + }, { + "name": "sector", + "type": "I", + "info": "atof sector" + }, { + "name": "layer", + "type": "I", + "info": "atof layer" + },{ + "name": "component", + "type": "I", + "info": "atof component" + },{ + "name": "time", + "type": "F", + "info": "time in ns" + },{ + "name": "x", + "type": "F", + "info": "x position in mm" + }, { + "name": "y", + "type": "F", + "info": "y position in mm" + }, { + "name": "z", + "type": "F", + "info": "z position in mm" + },{ + "name": "energy", + "type": "F", + "info": "deposited energy in MeV" + },{ + "name": "inlength", + "type": "F", + "info": "path length inside the detector (from entrance to hit) in mm" + },{ + "name": "pathlength", + "type": "F", + "info": "path length to the hit in mm" + } + ] + },{ "name": "AHDC::Hits", "group": 23000, "item": 23, diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/atof/Cluster/AtofCluster.java b/reconstruction/alert/src/main/java/org/jlab/rec/atof/Cluster/AtofCluster.java new file mode 100644 index 000000000..87c895082 --- /dev/null +++ b/reconstruction/alert/src/main/java/org/jlab/rec/atof/Cluster/AtofCluster.java @@ -0,0 +1,145 @@ +package org.jlab.rec.atof.cluster; +import java.util.ArrayList; +import org.jlab.rec.atof.hit.AtofHit; +import org.jlab.rec.atof.hit.BarHit; + +/** + * + * @authors npilleux, churaman + */ + + +public class AtofCluster { + + ArrayList bar_hits; + ArrayList wedge_hits; + double x,y,z,time,energy; + double path_length; + + public ArrayList getBarHits() { + return bar_hits; + } + + public void setBarHits(ArrayList bar_hits) { + this.bar_hits = bar_hits; + } + + public ArrayList getWedgeHits() { + return wedge_hits; + } + + public void setWedgeHits(ArrayList wedge_hits) { + this.wedge_hits = wedge_hits; + } + + public double getX() { + return x; + } + + public void setX(double x) { + this.x = x; + } + + public double getY() { + return y; + } + + public void setY(double y) { + this.y = y; + } + + public double getZ() { + return z; + } + + public void setZ(double z) { + this.z = z; + } + + public double getTime() { + return time; + } + + public void setTime(double time) { + this.time = time; + } + + public double getEnergy() { + return energy; + } + + public void setEnergy(double energy) { + this.energy = energy; + } + + public double getPath_length() { + return path_length; + } + + public void setPath_length(double path_length) { + this.path_length = path_length; + } + + //Cluster coordinates and time are defined as the coordinates and time of the max energy hit + //Can be changed later + + public final void computeClusterProperties() { + this.energy=0; + double max_energy = -1; + AtofHit max_energy_hit = new AtofHit(); + BarHit max_energy_barhit = new BarHit(); + + for(int i_wedge = 0; i_wedgemax_energy){max_energy_hit = this_wedge_hit; max_energy = this_energy;} + } + + for(int i_bar = 0; i_barmax_energy){max_energy_barhit = this_bar_hit; max_energy = this_energy;} + } + + if(max_energy_hit.getEnergy() > max_energy_barhit.getEnergy()) + { + this.time = max_energy_hit.getTime(); + this.x = max_energy_hit.getX(); + this.y = max_energy_hit.getY(); + this.z = max_energy_hit.getZ(); + this.path_length = max_energy_hit.getPath_length(); + } + else + { + this.time = max_energy_barhit.getTime(); + this.x = max_energy_barhit.getX(); + this.y = max_energy_barhit.getY(); + this.z = max_energy_barhit.getZ(); + this.path_length = max_energy_barhit.getPath_length(); + } + + } + + public double getPhi() + { + return Math.atan2(this.y, this.x); + } + + public AtofCluster(ArrayList bar_hits, ArrayList wedge_hits) + { + this.bar_hits = bar_hits; + this.wedge_hits = wedge_hits; + this.computeClusterProperties(); + } + + /** + * @param args the command line arguments + */ + public static void main(String[] args) { + } + +} diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/atof/Cluster/ClusterFinder.java b/reconstruction/alert/src/main/java/org/jlab/rec/atof/Cluster/ClusterFinder.java new file mode 100644 index 000000000..28e6af5ce --- /dev/null +++ b/reconstruction/alert/src/main/java/org/jlab/rec/atof/Cluster/ClusterFinder.java @@ -0,0 +1,231 @@ +package org.jlab.rec.atof.cluster; + +import cnuphys.magfield.MagneticFields; +import java.util.ArrayList; +import org.jlab.clas.swimtools.Swim; +import org.jlab.detector.calib.utils.DatabaseConstantProvider; +import org.jlab.geom.base.Detector; +import org.jlab.geom.detector.alert.ATOF.AlertTOFFactory; +import org.jlab.io.base.DataEvent; +import org.jlab.io.hipo.HipoDataSource; +import org.jlab.rec.atof.hit.AtofHit; +import org.jlab.rec.atof.hit.BarHit; +import org.jlab.rec.atof.hit.HitFinder; +import org.jlab.rec.atof.trackMatch.TrackProjector; +import org.jlab.utils.CLASResources; + +/** + * + * @authors npilleux,churaman + */ + +public class ClusterFinder { + + private ArrayList clusters; + + public void setClusters(ArrayList clusters) { + this.clusters = clusters; + } + + public ArrayList getClusters() { + return clusters; + } + + public void MakeClusters(HitFinder hitfinder) { + + //A list of clusters is built for each event + clusters.clear(); + + //Getting the list of hits, they must have been ordered by energy already + ArrayList wedge_hits = hitfinder.getWedgeHits(); + ArrayList bar_hits = hitfinder.getBarHits(); + + double sigma_Phi = 6.0; //angle opening of a layer. to be read from DB in the future + double sigma_Z = 6000;//to be read from DB in the future + double sigma_T = 1000;//timing resolution to be read from DB in the future + int sigma_module = 1; //hits are always within +-1 phi module of the most energetic + int sigma_component = 1;//hits are always within +-1 z component of the most energetic + + + //Looping through wedge hits first + for (int i_wedge = 0; i_wedge < wedge_hits.size(); i_wedge++) { + AtofHit this_wedge_hit = wedge_hits.get(i_wedge); + //Make a cluster for each wedge hit that has not been previously clustered + if (this_wedge_hit.getIs_in_a_cluster()) { + continue; + } + //Holding onto the hits composing the cluster + ArrayList this_cluster_wedge_hits = new ArrayList<>(); + ArrayList this_cluster_bar_hits = new ArrayList<>(); + + //Indicate that this hit now is in a cluster + this_wedge_hit.setIs_in_a_cluster(true); + //And store it + this_cluster_wedge_hits.add(this_wedge_hit); + + //Check if other wedge hits should be clustered with the current one + //Start from the index of the current one and look at less energetic hits + for (int j_wedge = i_wedge + 1; j_wedge < wedge_hits.size(); j_wedge++) { + AtofHit other_wedge_hit = wedge_hits.get(j_wedge); + //If that other hit is already involved in a cluster, skip it + if (other_wedge_hit.getIs_in_a_cluster()) { + continue; + } + //Check the distance between the hits + //For now we use phi module and z component differences from what is observed in simu + int delta_module = Math.abs(this_wedge_hit.computeModule_index() - other_wedge_hit.computeModule_index()); + if (delta_module > 30) { + delta_module = 60 - delta_module; + } + int delta_component = Math.abs(this_wedge_hit.getComponent() - other_wedge_hit.getComponent()); + //Later we could use z and phi threshold + //double delta_Phi = Math.abs(this_wedge_hit.getPhi() - other_wedge_hit.getPhi()); + //double delta_Z = Math.abs(this_wedge_hit.getZ() - other_wedge_hit.getZ()); + //Time matching + double delta_T = Math.abs(this_wedge_hit.getTime() - other_wedge_hit.getTime()); + + if (delta_module <= sigma_module)//delta_Phi <= sigma_Phi) + { + if (delta_component <= sigma_component)//delta_Z <= sigma_Z) + { + if (delta_T < sigma_T) { + other_wedge_hit.setIs_in_a_cluster(true); + this_cluster_wedge_hits.add(other_wedge_hit); + } + } + } + } + + //After clustering wedge hits, check if bar hits should be clustered with them + for (int j_bar = 0; j_bar < bar_hits.size(); j_bar++) { + BarHit other_bar_hit = bar_hits.get(j_bar); + //Skip already clustered hits + if (other_bar_hit.getIs_in_a_cluster()) { + continue; + } + //Check the distance between the hits + //For now we use phi module difference from what is observed in simu + int delta_module = Math.abs(this_wedge_hit.computeModule_index() - other_bar_hit.computeModule_index()); + if (delta_module > 30) { + delta_module = 60 - delta_module; + } + //Later we could use phi threshold + //double delta_Phi = Math.abs(this_wedge_hit.getPhi() - other_wedge_hit.getPhi()); + double delta_Z = Math.abs(this_wedge_hit.getZ() - other_bar_hit.getZ()); + //Time matching + double delta_T = Math.abs(this_wedge_hit.getTime() - other_bar_hit.getTime()); + if (delta_module <= sigma_module)//delta_Phi < sigma_Phi) + { + if (delta_Z < sigma_Z) { + if (delta_T < sigma_T) { + other_bar_hit.setIs_in_a_cluster(true); + this_cluster_bar_hits.add(other_bar_hit); + } + } + } + }//End loop bar hits + //After all wedge and bar hits have been grouped, build the cluster + AtofCluster cluster = new AtofCluster(this_cluster_bar_hits, this_cluster_wedge_hits); + //And add it to the list of clusters + clusters.add(cluster); + }//End loop on all wedge hits + + //Now make clusters from bar hits that are not associated with wedge hits + //Loop through all bar hits + for (int i_bar = 0; i_bar < bar_hits.size(); i_bar++) { + BarHit this_bar_hit = bar_hits.get(i_bar); + //Skip hits that have already been clustered + if (this_bar_hit.getIs_in_a_cluster()) { + continue; + } + + ArrayList this_cluster_wedge_hits = new ArrayList(); + ArrayList this_cluster_bar_hits = new ArrayList(); + this_bar_hit.setIs_in_a_cluster(true); + this_cluster_bar_hits.add(this_bar_hit); + + //Loop through less energetic clusters + for (int j_bar = i_bar + 1; j_bar < bar_hits.size(); j_bar++) { + BarHit other_bar_hit = bar_hits.get(j_bar); + //Skip already clustered hits + if (other_bar_hit.getIs_in_a_cluster()) { + continue; + } + + //Check the distance between the hits + //For now we use phi module difference from what is observed in simu + int delta_module = Math.abs(this_bar_hit.computeModule_index() - other_bar_hit.computeModule_index()); + if (delta_module > 30) { + delta_module = 60 - delta_module; + } + //Later we could use phi threshold + //double delta_Phi = Math.abs(this_wedge_hit.getPhi() - other_wedge_hit.getPhi()); + double delta_Z = Math.abs(this_bar_hit.getZ() - other_bar_hit.getZ()); + //Time matching + double delta_T = Math.abs(this_bar_hit.getTime() - other_bar_hit.getTime()); + + if (delta_module <= sigma_module)//delta_Phi < sigma_Phi) + { + if (delta_Z < sigma_Z) { + if (delta_T < sigma_T) { + other_bar_hit.setIs_in_a_cluster(true); + this_cluster_bar_hits.add(other_bar_hit); + } + } + } + } + AtofCluster cluster = new AtofCluster(this_cluster_bar_hits, this_cluster_wedge_hits); + clusters.add(cluster); + } + } + + public ClusterFinder() { + clusters = new ArrayList(); + } + + /** + * @param args the command line arguments + */ + public static void main(String[] args) { + // TODO code application logic here + AlertTOFFactory factory = new AlertTOFFactory(); + DatabaseConstantProvider cp = new DatabaseConstantProvider(11, "default"); + Detector atof = factory.createDetectorCLAS(cp); + + //READING MAG FIELD MAP + System.setProperty("CLAS12DIR", "../../"); + String mapDir = CLASResources.getResourcePath("etc") + "/data/magfield"; + try { + MagneticFields.getInstance().initializeMagneticFields(mapDir, + "Symm_torus_r2501_phi16_z251_24Apr2018.dat", "Symm_solenoid_r601_phi1_z1201_13June2018.dat"); + } catch (Exception e) { + e.printStackTrace(); + } + float[] b = new float[3]; + Swim swimmer = new Swim(); + swimmer.BfieldLab(0, 0, 0, b); + double B = Math.abs(b[2]); + + //Track Projector Initialisation with B field + TrackProjector projector = new TrackProjector(); + projector.setB(B); + + //Input to be read + String input = "/Users/npilleux/Desktop/alert/atof-reconstruction/coatjava/reconstruction/alert/src/main/java/org/jlab/rec/atof/Hit/update_protons_to_test_with_tracks.hipo"; + HipoDataSource reader = new HipoDataSource(); + reader.open(input); + + HitFinder hitfinder = new HitFinder(); + + int event_number = 0; + while (reader.hasEvent()) { + DataEvent event = (DataEvent) reader.getNextEvent(); + event_number++; + projector.ProjectTracks(event); + hitfinder.FindHits(event, atof, projector); + ClusterFinder clusterfinder = new ClusterFinder(); + clusterfinder.MakeClusters(hitfinder); + } + } + +} diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/atof/Hit/AtofHit.java b/reconstruction/alert/src/main/java/org/jlab/rec/atof/Hit/AtofHit.java new file mode 100644 index 000000000..88d0a4e16 --- /dev/null +++ b/reconstruction/alert/src/main/java/org/jlab/rec/atof/Hit/AtofHit.java @@ -0,0 +1,526 @@ +package org.jlab.rec.atof.hit; + +import java.util.List; +import org.jlab.geom.base.*; +import org.jlab.geom.detector.alert.ATOF.*; +import org.jlab.detector.calib.utils.DatabaseConstantProvider; +import org.jlab.geom.prim.Point3D; +import org.jlab.io.base.DataBank; +import org.jlab.io.base.DataEvent; +import org.jlab.io.hipo.HipoDataSource; +import org.jlab.rec.atof.constants.Parameters; +import org.jlab.rec.atof.trackMatch.TrackProjection; +import org.jlab.rec.atof.trackMatch.TrackProjector; + +/** + * + * Represents a hit in the atof. Stores info about the sector, layer, component, + * order, TDC, ToT. Type is wedge/bar up/bar down/ bar Spatial coordinates are + * computed from atof detector object using the geometry service Stores whether + * the hit is part of a cluster. Calculates time, energy based on TDC/ToT. + * + * @authors npilleux, churaman + */ +public class AtofHit { + + private int sector, layer, component, order; + private int TDC, ToT; + private double time, energy, x, y, z; + private String type; + private boolean is_in_a_cluster; + private double path_length, inpath_length; + + public int getSector() { + return sector; + } + + public void setSector(int sector) { + this.sector = sector; + } + + public int getLayer() { + return layer; + } + + public void setLayer(int layer) { + this.layer = layer; + } + + public int getOrder() { + return order; + } + + public void setOrder(int order) { + this.order = order; + } + + public int getComponent() { + return component; + } + + public void setComponent(int component) { + this.component = component; + } + + public int getTDC() { + return TDC; + } + + public void setTDC(int tdc) { + this.TDC = tdc; + } + + public int getToT() { + return ToT; + } + + public void setToT(int tot) { + this.ToT = tot; + } + + public double getTime() { + return time; + } + + public void setTime(double time) { + this.time = time; + } + + public double getEnergy() { + return energy; + } + + public void setEnergy(double energy) { + this.energy = energy; + } + + public double getX() { + return x; + } + + public void setX(double x) { + this.x = x; + } + + public double getY() { + return y; + } + + public void setY(double y) { + this.y = y; + } + + public double getZ() { + return z; + } + + public void setZ(double z) { + this.z = z; + } + + public String getType() { + return type; + } + + public void setType(String type) { + this.type = type; + } + + public boolean getIs_in_a_cluster() { + return is_in_a_cluster; + } + + public void setIs_in_a_cluster(boolean is_in_a_cluster) { + this.is_in_a_cluster = is_in_a_cluster; + } + + public double getPath_length() { + return path_length; + } + + public void setPath_length(double path_length) { + this.path_length = path_length; + } + + public double getInpath_length() { + return inpath_length; + } + + public void setInpath_length(double inpath_length) { + this.inpath_length = inpath_length; + } + + public int computeModule_index() { + //Index ranging 0 to 60 for each wedge+bar module + return 4 * this.sector + this.layer; + } + + public final String makeType() { + String itype = "undefined"; + if (this.component == 10 && this.order == 0) { + itype = "bar down"; + } else if (this.component == 10 && this.order == 1) { + itype = "bar up"; + } else if (this.component < 10) { + itype = "wedge"; + } + this.type = itype; + return itype; + } + + public final int TDC_to_time() { + double tdc2time = 1; + double veff = 1; + if (null == this.type) { + return 1; + } else { + switch (this.type) { + case "wedge" -> { + tdc2time = Parameters.TDC2TIME; + veff = Parameters.VEFF; + } + case "bar up" -> { + tdc2time = Parameters.TDC2TIME; + veff = Parameters.VEFF; + } + case "bar down" -> { + tdc2time = Parameters.TDC2TIME; + veff = Parameters.VEFF; + } + case "bar" -> { + tdc2time = Parameters.TDC2TIME; + veff = Parameters.VEFF; + } + default -> { + return 1; + } + } + } + //Time at the inner surface + this.time = tdc2time * this.TDC - this.inpath_length / veff; + return 0; + } + + public final int ToT_to_energy() { + double tot2energy = 1; + if (null == this.type) { + return 1; + } else { + switch (this.type) { + case "wedge" -> { + tot2energy = Parameters.TOT2ENERGY_WEDGE; + //For now hits are considered in the middle of the wedge + //And the SiPM on top + double distance_hit_to_sipm = Parameters.WEDGE_THICKNESS/2.; + this.energy = tot2energy * this.ToT * Math.exp(distance_hit_to_sipm / Parameters.ATT_L); + } + case "bar up" -> { + tot2energy = Parameters.TOT2ENERGY_BAR; + //only half the information in the bar, + //the attenuation will be computed when the full hit is formed + this.energy = tot2energy * this.ToT; + } + case "bar down" -> { + tot2energy = Parameters.TOT2ENERGY_BAR; + //only half the information in the bar, + //the attenuation will be computed when the full hit is formed + this.energy = tot2energy * this.ToT; + } + default -> { + return 1; + } + } + } + return 0; + } + + /** + * Calculates spatial coordinates for the hit based on associated detector + * component. Retrieves the midpoint of the atof component to assign the + * corresponding x, y, z coordinates to the hit. + * + * + * @param atof The Detector object representing the atof. + * @return 0 if the coordinates were successfully set, or 1 if the hit type + * is undefined or unsupported. + */ + public final int slc_to_xyz(Detector atof) { + int sl = 999; + if (null == this.type) { + return 1; + } else { + switch (this.type) { + case "wedge" -> + sl = 1; + case "bar up", "bar down", "bar" -> + sl = 0; + default -> { + return 1; + } + } + } + Component comp = atof.getSector(this.sector).getSuperlayer(sl).getLayer(this.layer).getComponent(this.component); + Point3D midpoint = comp.getMidpoint(); + + //coordinates centered at the center of the atof in mm + this.x = midpoint.x() - Parameters.LENGTH_ATOF / 2.; + this.y = midpoint.y() - Parameters.LENGTH_ATOF / 2.; + this.z = midpoint.z() - Parameters.LENGTH_ATOF / 2.; + return 0; + } + + /** + * Compares two AtofHit objects to check if they match in the bar. + *
    + *
  • If the sector or layer of the two hits do not match, the method + * returns {@code false}.
    • + *
  • If either hit is not in the bar (component must be 10), the method + * returns {@code false}.
    • + *
  • If both hits are in the same SiPM (i.e., their order is the same), + * the method returns {@code false}.
    • + *
+ * If none of these conditions are violated, the method returns + * {@code true}, indicating the two hits match. + * + * @param hit2match The AtofHit object to compare with the current instance. + * @return {@code true} if the hits match; {@code false} otherwise. + */ + + public boolean matchBar(AtofHit hit2match) { + if (this.getSector() != hit2match.getSector()) { + return false; //System.out.print("Two hits in different sectors \n"); + } else if (this.getLayer() != hit2match.getLayer()) { + return false; //System.out.print("Two hits in different layers \n"); + } else if (this.getComponent() != 10 || hit2match.getComponent() != 10) { + return false; //System.out.print("At least one hit is not in the bar \n"); + } else if (this.getOrder() > 1 || hit2match.getComponent() > 1) { + return false; //System.out.print("At least one hit is not a downstram or upstream hit. It may be a full bar hit. \n"); + } else { + return this.getOrder() != hit2match.getOrder(); //System.out.print("Two hits in same SiPM \n"); + } + } + + /** + * Returns the azimuthal angle (phi) of the hit. + * + * @return The azimuthal angle (phi) in radians, in the range [-π, π]. + */ + public double getPhi() { + return Math.atan2(this.y, this.x); + } + + /** + * Constructor for a fit in the atof. Initializes the hit's sector, layer, + * component, order, TDC, ToT. Sets the hit's initial state regarding + * clustering. Set up the hit's type, time, energy, and spatial coordinates. + * + * @param sector The sector of the detector where the hit occurred. + * @param layer The layer of the detector where the hit was detected. + * @param component The component within the layer that registered the hit. + * @param order Order of the hit. + * @param tdc TDC value. + * @param tot ToT value. + * @param atof Detector object representing the atof, used to calculate + * spatial coordinates. + */ + public AtofHit(int sector, int layer, int component, int order, int tdc, int tot, Detector atof) { + this.sector = sector; + this.layer = layer; + this.component = component; + this.order = order; + this.TDC = tdc; + this.ToT = tot; + this.is_in_a_cluster = false; + + this.makeType(); + int is_ok = this.TDC_to_time(); + if (is_ok != 1) { + is_ok = this.ToT_to_energy(); + } + if (is_ok != 1) { + is_ok = this.slc_to_xyz(atof); + } + } + + /** + * Constructor for a fit in the atof. Initializes the hit's sector, layer, + * component, order, TDC, ToT. Sets the hit's initial state regarding + * clustering. Set up the hit's type, time, energy, and spatial coordinates. + * + * @param sector The sector of the detector where the hit occurred. + * @param layer The layer of the detector where the hit was detected. + * @param component The component within the layer that registered the hit. + * @param order Order of the hit. + * @param tdc TDC value. + * @param tot ToT value. + * @param atof Detector object representing the atof, used to calculate + * @param track_projector TrackProjector object with ahdc tracks to be + * matched to the hit. + */ + public AtofHit(int sector, int layer, int component, int order, int tdc, int tot, Detector atof, TrackProjector track_projector) { + this.sector = sector; + this.layer = layer; + this.component = component; + this.order = order; + this.TDC = tdc; + this.ToT = tot; + this.is_in_a_cluster = false; + + //First the type needs to be set + this.makeType(); + //From it the coordinates can be computed + this.slc_to_xyz(atof); + //From them tracks can be matched + this.matchTrack(track_projector); + //And energy and time can then be computed + this.ToT_to_energy(); + this.TDC_to_time(); + } + + /** + * Matches the current track with ahdc tracks projections. Calculates the + * match by comparing the hit's azimuthal angle and longitudinal position + * (z) with the track projection. If a match is found within defined + * tolerances for phi and z, the path length of the matched hit is updated. + * + * @param track_projector The TrackProjector object that provides a list of + * TrackProjections. + * + */ + + public final void matchTrack(TrackProjector track_projector) { + double sigma_phi = 0; + double sigma_z = 0; + + List Projections = track_projector.getProjections(); + for (int i_track = 0; i_track < Projections.size(); i_track++) { + Point3D projection_point = new Point3D(); + if (null == this.getType()) { + System.out.print("Impossible to match track and hit; hit type is null \n"); + } else { + switch (this.getType()) { + case "wedge" -> { + sigma_phi = Parameters.SIGMA_PHI_TRACK_MATCHING_WEDGE; + sigma_z = Parameters.SIGMA_Z_TRACK_MATCHING_WEDGE; + projection_point = Projections.get(i_track).get_WedgeIntersect(); + } + case "bar up", "bar down" -> { + System.out.print("WARNING : YOU ARE MATCHING A TRACK TO A SINGLE HIT IN THE BAR. \n"); + sigma_phi = Parameters.SIGMA_PHI_TRACK_MATCHING_BAR; + sigma_z = Parameters.SIGMA_Z_TRACK_MATCHING_BAR; + projection_point = Projections.get(i_track).get_BarIntersect(); + } + default -> + System.out.print("Impossible to match track and hit; hit type is undefined \n"); + } + } + if (Math.abs(this.getPhi() - projection_point.toVector3D().phi()) < sigma_phi) { + if (Math.abs(this.getZ() - projection_point.z()) < sigma_z) { + if ("wedge".equals(this.getType())) { + this.setPath_length(Projections.get(i_track).get_WedgePathLength()); + } else { + this.setPath_length(Projections.get(i_track).get_BarPathLength()); + } + } + } + } + } + + public int matchTrack(DataEvent event) { + + String track_bank_name = "AHDC::Projections"; + if (event == null) { // check if there is an event + //System.out.print(" no event \n"); + return 1; + } else if (event.hasBank(track_bank_name) == false) { + return 1; + // check if there are ahdc tracks in the event + //System.out.print("no tracks \n"); + } else { + DataBank track_bank = event.getBank(track_bank_name); + int nt = track_bank.rows(); // number of tracks + double sigma_phi = 0; + double sigma_z = 0; + + //Looping through all tracks + for (int i = 0; i < nt; i++) { + + Float xt = null, yt = null, zt = null, path = null, inpath = null; + + if (null == this.getType()) { + System.out.print("Impossible to match track and hit; hit type is null \n"); + } else { + switch (this.getType()) { + case "wedge" -> { + sigma_phi = Parameters.SIGMA_PHI_TRACK_MATCHING_WEDGE; + sigma_z = Parameters.SIGMA_Z_TRACK_MATCHING_WEDGE; + xt = track_bank.getFloat("x_at_wedge", i); + yt = track_bank.getFloat("y_at_wedge", i); + zt = track_bank.getFloat("z_at_wedge", i); + path = track_bank.getFloat("L_at_wedge", i); + inpath = track_bank.getFloat("L_in_wedge", i); + } + case "bar up", "bar down" -> { + System.out.print("WARNING : YOU ARE MATCHING A TRACK TO A SINGLE HIT IN THE BAR. \n"); + sigma_phi = Parameters.SIGMA_PHI_TRACK_MATCHING_BAR; + sigma_z = Parameters.SIGMA_Z_TRACK_MATCHING_BAR; + xt = track_bank.getFloat("x_at_bar", i); + yt = track_bank.getFloat("y_at_bar", i); + zt = track_bank.getFloat("z_at_bar", i); + path = track_bank.getFloat("L_at_bar", i); + inpath = track_bank.getFloat("L_in_bar", i); + } + default -> + System.out.print("Impossible to match track and hit; hit type is undefined \n"); + } + } + Point3D projection_point = new Point3D(xt, yt, zt); + if (Math.abs(this.getPhi() - projection_point.toVector3D().phi()) < sigma_phi) { + if (Math.abs(this.getZ() - projection_point.z()) < sigma_z) { + System.out.print("PASSED CUTS \n"); + this.setPath_length(path); + this.setInpath_length(inpath); + } + } + } + } + return 0; + } + + public AtofHit() { + } + + /** + * @param args the command line arguments + */ + public static void main(String[] args) { + // TODO code application logic here + AlertTOFFactory factory = new AlertTOFFactory(); + DatabaseConstantProvider cp = new DatabaseConstantProvider(11, "default"); + Detector atof = factory.createDetectorCLAS(cp); + + //Input to be read + String input = "/Users/npilleux/Desktop/alert/atof-reconstruction/coatjava/reconstruction/alert/src/main/java/org/jlab/rec/atof/Hit/test_tdc_atof.hipo"; + HipoDataSource reader = new HipoDataSource(); + reader.open(input); + + int event_number = 0; + while (reader.hasEvent()) { + DataEvent event = (DataEvent) reader.getNextEvent(); + event_number++; + DataBank bank = event.getBank("ATOF::tdc"); + int nt = bank.rows(); // number of tracks + + for (int i = 0; i < nt; i++) { + int sector = bank.getInt("sector", i); + int layer = bank.getInt("layer", i); + int component = bank.getInt("component", i); + int order = bank.getInt("order", i); + int tdc = bank.getInt("TDC", i); + int tot = bank.getInt("ToT", i); + AtofHit hit = new AtofHit(sector, layer, component, order, tdc, tot, atof); + System.out.print(hit.getX() + "\n"); + } + } + } +} diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/atof/Hit/BarHit.java b/reconstruction/alert/src/main/java/org/jlab/rec/atof/Hit/BarHit.java new file mode 100644 index 000000000..3edae6bd5 --- /dev/null +++ b/reconstruction/alert/src/main/java/org/jlab/rec/atof/Hit/BarHit.java @@ -0,0 +1,140 @@ +package org.jlab.rec.atof.hit; +import org.jlab.detector.calib.utils.DatabaseConstantProvider; +import org.jlab.geom.base.Detector; +import org.jlab.geom.detector.alert.ATOF.AlertTOFFactory; +import org.jlab.io.base.DataBank; +import org.jlab.io.base.DataEvent; +import org.jlab.io.hipo.HipoDataSource; +import java.util.ArrayList; +import org.jlab.rec.atof.constants.Parameters; + +/** + * + * Represents a hit in the atof bar. Extends class AtofHit. Is further defined + * by the two hits upstream and downstream composing a full bar hit. z position, + * time and energy are defined from the up/down hits + * + * @authors npilleux,churaman + */ + +public class BarHit extends AtofHit { + + //A bar hit is the combination of a downstream and upstream hits + private AtofHit hit_up, hit_down; + + public AtofHit getHitUp() { + return hit_up; + } + + public void setHitUp(AtofHit hit_up) { + this.hit_up = hit_up; + } + + public AtofHit getHitDown() { + return hit_down; + } + + public void setHitDown(AtofHit hit_down) { + this.hit_down = hit_down; + } + + public final void computeZ() { + double some_calibration = 10; //Here read the calibration DB + this.setZ(some_calibration * (hit_up.getTime() - hit_down.getTime())); + } + + public final void computeTime() { + double some_calibration = 10; //Here read the calibration DB + this.setTime(some_calibration + ((hit_up.getTime() + hit_down.getTime()) / 2.)); + } + + public final void computeEnergy() { + this.computeZ(); + double distance_hit_to_sipm_up = Parameters.LENGTH_ATOF / 2. + this.getZ(); + double distance_hit_to_sipm_down = Parameters.LENGTH_ATOF / 2. - this.getZ(); + double Edep_up = hit_up.getEnergy() * Math.exp(distance_hit_to_sipm_up / Parameters.ATT_L); + double Edep_down = hit_down.getEnergy() * Math.exp(distance_hit_to_sipm_down / Parameters.ATT_L); + this.setEnergy(Edep_up + Edep_down); + } + + public BarHit(AtofHit hit_down, AtofHit hit_up) { + boolean hits_match = hit_down.matchBar(hit_up); + if (!hits_match) { + throw new UnsupportedOperationException("Hits do not match \n"); + } + this.setType("bar"); + this.setOrder(2);//Fake order for bar hits + this.hit_up = hit_up; + this.hit_down = hit_down; + this.setLayer(hit_up.getLayer()); + this.setSector(hit_up.getSector()); + this.setX(hit_up.getX()); + this.setY(hit_up.getY()); + this.computeZ(); + this.computeTime(); + this.computeEnergy(); + } + + public BarHit() { + } + + /** + * @param args the command line arguments + */ + public static void main(String[] args) { + // TODO code application logic here + AlertTOFFactory factory = new AlertTOFFactory(); + DatabaseConstantProvider cp = new DatabaseConstantProvider(11, "default"); + Detector atof = factory.createDetectorCLAS(cp); + + //Input to be read + String input = "/Users/npilleux/Desktop/alert/atof-reconstruction/coatjava/reconstruction/alert/src/main/java/org/jlab/rec/atof/Hit/test_tdc_atof.hipo"; + HipoDataSource reader = new HipoDataSource(); + reader.open(input); + + int event_number = 0; + while (reader.hasEvent()) { + DataEvent event = (DataEvent) reader.getNextEvent(); + event_number++; + DataBank bank = event.getBank("ATOF::tdc"); + int nt = bank.rows(); // number of tracks + ArrayList hit_up = new ArrayList<>(); + ArrayList hit_down = new ArrayList<>(); + ArrayList hit_wedge = new ArrayList<>(); + for (int i = 0; i < nt; i++) { + int sector = bank.getInt("sector", i); + int layer = bank.getInt("layer", i); + int component = bank.getInt("component", i); + int order = bank.getInt("order", i); + int tdc = bank.getInt("TDC", i); + int tot = bank.getInt("ToT", i); + AtofHit hit = new AtofHit(sector, layer, component, order, tdc, tot, atof); + if (null == hit.getType()) { + System.out.print("Undefined hit type \n"); + } else { + switch (hit.getType()) { + case "bar up" -> + hit_up.add(hit); + case "bar down" -> + hit_down.add(hit); + case "wedge" -> + hit_wedge.add(hit); + default -> + System.out.print("Undefined hit type \n"); + } + } + } + ArrayList bar_hits = new ArrayList<>(); + for (int i_up = 0; i_up < hit_up.size(); i_up++) { + AtofHit this_hit_up = hit_up.get(i_up); + for (int i_down = 0; i_down < hit_down.size(); i_down++) { + AtofHit this_hit_down = hit_down.get(i_down); + if (this_hit_up.matchBar(this_hit_down)) { + BarHit this_bar_hit = new BarHit(this_hit_up, this_hit_down); + bar_hits.add(this_bar_hit); + } + } + } + } + } +} diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/atof/Hit/HitFinder.java b/reconstruction/alert/src/main/java/org/jlab/rec/atof/Hit/HitFinder.java new file mode 100644 index 000000000..8f33db7f6 --- /dev/null +++ b/reconstruction/alert/src/main/java/org/jlab/rec/atof/Hit/HitFinder.java @@ -0,0 +1,294 @@ +package org.jlab.rec.atof.hit; + +import cnuphys.magfield.MagneticFields; +import java.util.ArrayList; +import java.util.Collections; +import javax.swing.JFrame; +import org.jlab.clas.swimtools.Swim; +import org.jlab.detector.calib.utils.DatabaseConstantProvider; +import org.jlab.geom.base.Detector; +import org.jlab.geom.detector.alert.ATOF.AlertTOFFactory; +import org.jlab.groot.data.H1F; +import org.jlab.groot.graphics.EmbeddedCanvas; +import org.jlab.io.base.DataBank; +import org.jlab.io.base.DataEvent; +import org.jlab.io.hipo.HipoDataSource; +import static org.jlab.rec.atof.banks.RecoBankWriter.fillAtofHitBank; +import org.jlab.rec.atof.trackMatch.TrackProjector; +import org.jlab.utils.CLASResources; + +/** + * + * @authors npilleux,churaman + */ +public class HitFinder { + + private ArrayList bar_hits; + private ArrayList wedge_hits; + + public HitFinder() { + this.bar_hits = new ArrayList<>(); + this.wedge_hits = new ArrayList<>(); + } + + // Getter and Setter for bar_hits + public ArrayList getBarHits() { + return bar_hits; + } + + public void setBarHits(ArrayList bar_hits) { + this.bar_hits = bar_hits; + } + + // Getter and Setter for wedge_hits + public ArrayList getWedgeHits() { + return wedge_hits; + } + + public void setWedgeHits(ArrayList wedge_hits) { + this.wedge_hits = wedge_hits; + } + + public void FindHits(DataEvent event, Detector atof, TrackProjector track_projector) { + + //For each event a list of bar hits and a list of wedge hits are filled + this.bar_hits.clear(); + this.wedge_hits.clear(); + //They are read from the ATOF TDC bank + DataBank bank = event.getBank("ATOF::tdc"); + int nt = bank.rows(); // number of hits + //Hits in the bar downstream and upstream will be matched + ArrayList hit_up = new ArrayList<>(); + ArrayList hit_down = new ArrayList<>(); + //Looping through all hits + for (int i = 0; i < nt; i++) { + //Getting their properties + int sector = bank.getInt("sector", i); + int layer = bank.getInt("layer", i); + int component = bank.getInt("component", i); + int order = bank.getInt("order", i); + int tdc = bank.getInt("TDC", i); + int tot = bank.getInt("ToT", i); + //Building a Hit + AtofHit hit = new AtofHit(sector, layer, component, order, tdc, tot, atof); + if (hit.getEnergy() < 0.1) { + continue; //energy threshold + } //Sorting the hits into wedge, upstream and downstream bar hits + //Lists are built for up/down bar to match them after + //Wedge hits are mayched to ahdc tracks and listed + if (null == hit.getType()) { + System.out.print("Undefined hit type \n"); + } else { + switch (hit.getType()) { + case "bar up" -> + hit_up.add(hit); + case "bar down" -> + hit_down.add(hit); + case "wedge" -> { + hit.matchTrack(track_projector); + this.wedge_hits.add(hit); + } + default -> + System.out.print("Undefined hit type \n"); + } + } + }//End loop through all hits + + //Starting loop through up hits in the bar + for (int i_up = 0; i_up < hit_up.size(); i_up++) { + AtofHit this_hit_up = hit_up.get(i_up); + //Starting loop through down hits in the bar + for (int i_down = 0; i_down < hit_down.size(); i_down++) { + AtofHit this_hit_down = hit_down.get(i_down); + //Matching the hits: if same module and different order, they make up a bar hit + if (this_hit_up.matchBar(this_hit_down)) { + //Bar hits are matched to ahdc tracks and listed + BarHit this_bar_hit = new BarHit(this_hit_up, this_hit_down); + this_bar_hit.matchTrack(track_projector); + this.bar_hits.add(this_bar_hit); + } + } + } + //Once all has been listed, hits are sorted by energy + Collections.sort(this.bar_hits, (hit1, hit2) -> Double.compare(hit1.getEnergy(), hit2.getEnergy())); + Collections.sort(this.wedge_hits, (hit1, hit2) -> Double.compare(hit1.getEnergy(), hit2.getEnergy())); + } + + public void FindHits(DataEvent event, Detector atof) { + + //For each event a list of bar hits and a list of wedge hits are filled + this.bar_hits.clear(); + this.wedge_hits.clear(); + //They are read from the ATOF TDC bank + DataBank bank_atof_hits = event.getBank("ATOF::tdc"); + int nt = bank_atof_hits.rows(); // number of hits + //Hits in the bar downstream and upstream will be matched + ArrayList hit_up = new ArrayList<>(); + ArrayList hit_down = new ArrayList<>(); + //Looping through all hits + for (int i = 0; i < nt; i++) { + //Getting their properties + int sector = bank_atof_hits.getInt("sector", i); + int layer = bank_atof_hits.getInt("layer", i); + int component = bank_atof_hits.getInt("component", i); + int order = bank_atof_hits.getInt("order", i); + int tdc = bank_atof_hits.getInt("TDC", i); + int tot = bank_atof_hits.getInt("ToT", i); + //Building a Hit + AtofHit hit = new AtofHit(sector, layer, component, order, tdc, tot, atof); + if (hit.getEnergy() < 0.1) { + continue; //energy threshold + } //Sorting the hits into wedge, upstream and downstream bar hits + //Lists are built for up/down bar to match them after + //Wedge hits are mayched to ahdc tracks and listed + if (null == hit.getType()) { + System.out.print("Undefined hit type \n"); + } else { + switch (hit.getType()) { + case "bar up" -> + hit_up.add(hit); + case "bar down" -> + hit_down.add(hit); + case "wedge" -> { + hit.matchTrack(event); + this.wedge_hits.add(hit); + } + default -> + System.out.print("Undefined hit type \n"); + } + } + }//End loop through all hits + + //Starting loop through up hits in the bar + for (int i_up = 0; i_up < hit_up.size(); i_up++) { + AtofHit this_hit_up = hit_up.get(i_up); + //Starting loop through down hits in the bar + for (int i_down = 0; i_down < hit_down.size(); i_down++) { + AtofHit this_hit_down = hit_down.get(i_down); + //Matching the hits: if same module and different order, they make up a bar hit + if (this_hit_up.matchBar(this_hit_down)) { + //Bar hits are matched to ahdc tracks and listed + BarHit this_bar_hit = new BarHit(this_hit_up, this_hit_down); + this_bar_hit.matchTrack(event); + this.bar_hits.add(this_bar_hit); + } + } + } + //Once all has been listed, hits are sorted by energy + Collections.sort(this.bar_hits, (hit1, hit2) -> Double.compare(hit1.getEnergy(), hit2.getEnergy())); + Collections.sort(this.wedge_hits, (hit1, hit2) -> Double.compare(hit1.getEnergy(), hit2.getEnergy())); + ArrayList allhits = new ArrayList<>();; + allhits.addAll(this.wedge_hits); + allhits.addAll(this.bar_hits); + DataBank hitbank = fillAtofHitBank(event, allhits); + event.appendBank(hitbank); + } + + /** + * @param args the command line arguments + */ + public static void main(String[] args) { + + //Building ALERT geometry + AlertTOFFactory factory = new AlertTOFFactory(); + DatabaseConstantProvider cp = new DatabaseConstantProvider(11, "default"); + Detector atof = factory.createDetectorCLAS(cp); + + //READING MAG FIELD MAP + System.setProperty("CLAS12DIR", "../../"); + String mapDir = CLASResources.getResourcePath("etc") + "/data/magfield"; + try { + MagneticFields.getInstance().initializeMagneticFields(mapDir, + "Symm_torus_r2501_phi16_z251_24Apr2018.dat", "Symm_solenoid_r601_phi1_z1201_13June2018.dat"); + } catch (Exception e) { + e.printStackTrace(); + } + float[] b = new float[3]; + Swim swimmer = new Swim(); + swimmer.BfieldLab(0, 0, 0, b); + double B = Math.abs(b[2]); + + //Track Projector Initialisation with B field + TrackProjector projector = new TrackProjector(); + projector.setB(B); + + //Hit finder init + HitFinder hitfinder = new HitFinder(); + + //Input to be read + String input = "/Users/npilleux/Desktop/alert/atof-reconstruction/coatjava/reconstruction/alert/src/main/java/org/jlab/rec/atof/hit/updated_updated_rec_more_protons_50_to_650.hipo"; + HipoDataSource reader = new HipoDataSource(); + reader.open(input); + + H1F h_delta_energy = new H1F("Energy", "Energy", 100, -10, 10); + h_delta_energy.setTitleX("delta energy [GeV]"); + + int event_number = 0; + while (reader.hasEvent()) { + DataEvent event = (DataEvent) reader.getNextEvent(); + event_number++; + projector.ProjectTracks(event); + hitfinder.FindHits(event, atof); + DataBank MC_True = event.getBank("MC::True"); + DataBank tdc = event.getBank("ATOF::tdc"); + DataBank hits = event.getBank("ATOF::hits"); + double totEdep = 0; + + for (int i = 0; i < MC_True.rows(); i++) { + + if (MC_True.getByte("detector", i) != 24) { + continue; + } + Float true_energy = MC_True.getFloat("avgT", i); + if (true_energy < 5) { + continue; + } + System.out.print(true_energy + " TRUE \n"); + double min_diff = 9999.; + double energy_at_min = 9999.; + + for (int j = 0; j < hits.rows(); j++) { + Float hit_energy = hits.getFloat("time", j); + Float diff = true_energy - hit_energy; + if (diff < min_diff) { + min_diff = diff; + energy_at_min = true_energy; + } + } + System.out.print("ICI " + energy_at_min + " " + true_energy + " \n"); + h_delta_energy.fill(min_diff / energy_at_min); + } + System.out.print("------------------- \n"); + } + + System.out.print ( + + "Read " + event_number + " events"); + JFrame frame = new JFrame("Raster"); + + frame.setSize ( + + 2500,800); + EmbeddedCanvas canvas = new EmbeddedCanvas(); + + canvas.divide ( + + + 3,2); + canvas.cd ( + + + 3); canvas.draw (h_delta_energy); + + frame.add (canvas); + + frame.setLocationRelativeTo ( + + + null); + frame.setVisible ( + + +true); + } +} diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/atof/TrackMatch/TrackProjection.java b/reconstruction/alert/src/main/java/org/jlab/rec/atof/TrackMatch/TrackProjection.java new file mode 100644 index 000000000..5138d1508 --- /dev/null +++ b/reconstruction/alert/src/main/java/org/jlab/rec/atof/TrackMatch/TrackProjection.java @@ -0,0 +1,176 @@ +package org.jlab.rec.atof.trackMatch; + +import org.jlab.geom.prim.Point3D; + +/** + * The {@code TrackProjection} class holds ahdc track information relevant for atof analysis + * i.e projected to the surfaces of the bar and wedges + * @author pilleux + */ + +public class TrackProjection { + + /** + * Intersection point of the track with the middle surface of the bar. + */ + private Point3D _BarIntersect = new Point3D(); + + /** + * Intersection point of the track with the middle surface of the wedges. + */ + private Point3D _WedgeIntersect = new Point3D(); + + /** + * Path length of the track from the DOCA to the beam line + * to the entrance surface of the bar. + */ + Float _BarPathLength; + + /** + * Path length of the track from the DOCA to the beam line + * to the entrance surface of the wedges. + */ + Float _WedgePathLength; + + /** + * Path length inside the bar. + */ + Float _BarInPathLength; + + /** + * Path length inside the wedge. + */ + Float _WedgeInPathLength; + + + /** + * Default constructor that initializes the intersection points and path lengths to {@code NaN}. + */ + public TrackProjection() { + _BarIntersect = new Point3D(Double.NaN, Double.NaN, Double.NaN); + _WedgeIntersect = new Point3D(Double.NaN, Double.NaN, Double.NaN); + _BarPathLength = Float.NaN; + _WedgePathLength = Float.NaN; + _BarInPathLength = Float.NaN; + _WedgeInPathLength = Float.NaN; + } + + /** + * Gets the intersection point of the track with the middle surface of the bar. + * + * @return {@link Point3D} bar's intersection point. + */ + public Point3D get_BarIntersect() { + return _BarIntersect; + } + + /** + * Gets the intersection point of the track with the middle surface of the wedges. + * + * @return {@link Point3D} wedge's intersection point. + */ + public Point3D get_WedgeIntersect() { + return _WedgeIntersect; + } + + /** + * Gets the path length of the track from the DOCA to the beam line to the inner surface of the bar. + * + * @return {@code Float} path length to the bar's middle surface. + */ + public Float get_BarPathLength() { + return _BarPathLength; + } + + /** + * Gets the path length of the track from the inner surface of the bar + * to its middle surface. + * + * @return {@code Float} path length inside the bar. + */ + public Float get_BarInPathLength() { + return _BarInPathLength; + } + + /** + * Gets the path length of the track from the DOCA to the beam line to the inner surface of the wedges. + * + * @return {@code Float} path length to the wedge's middle surface. + */ + public Float get_WedgePathLength() { + return _WedgePathLength; + } + + /** + * Gets the path length of the track from the the inner surface of the wedge + * to its middle surface. + * + * @return {@code Float} path length inside the wedge. + */ + public Float get_WedgeInPathLength() { + return _WedgeInPathLength; + } + + /** + * Sets the intersection point of the track with the middle surface of the bar. + * + * @param BarIntersect {@link Point3D} intersection with the bar. + */ + public void set_BarIntersect(Point3D BarIntersect) { + this._BarIntersect = BarIntersect; + } + + /** + * Sets the intersection point of the track with the middle surface of the wedges. + * + * @param WedgeIntersect {@link Point3D} intersection with the wedge. + */ + public void set_WedgeIntersect(Point3D WedgeIntersect) { + this._WedgeIntersect = WedgeIntersect; + } + + /** + * Sets the path length of the track from the DOCA to the beam line to the inner surface of the bar. + * + * @param BarPathLength {@code Float} path length to the bar inner surface. + */ + public void set_BarPathLength(Float BarPathLength) { + this._BarPathLength = BarPathLength; + } + + /** + * Sets the path length of the track from the DOCA to the beam line to the inner surface of the wedges. + * + * @param WedgePathLength {@code Float} path length to the wedge inner surface. + */ + public void set_WedgePathLength(Float WedgePathLength) { + this._WedgePathLength = WedgePathLength; + } + + /** + * Sets the path length of the track inside the bar. + * + * @param BarInPathLength {@code Float} path length inside the bar. + */ + public void set_BarInPathLength(Float BarInPathLength) { + this._BarInPathLength = BarInPathLength; + } + + /** + * Sets the path length of the track inside the wedges. + * + * @param WedgeInPathLength {@code Float} path length inside the wedge. + */ + public void set_WedgeInPathLength(Float WedgeInPathLength) { + this._WedgeInPathLength = WedgeInPathLength; + } + + + /** + * testing purposes. + * + * @param arg command-line arguments. + */ + public static void main(String arg[]) { + } +} \ No newline at end of file diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/atof/TrackMatch/TrackProjector.java b/reconstruction/alert/src/main/java/org/jlab/rec/atof/TrackMatch/TrackProjector.java new file mode 100644 index 000000000..90d5936b8 --- /dev/null +++ b/reconstruction/alert/src/main/java/org/jlab/rec/atof/TrackMatch/TrackProjector.java @@ -0,0 +1,269 @@ +package org.jlab.rec.atof.trackMatch; + +import java.util.ArrayList; +import java.util.List; + +import org.jlab.io.base.DataBank; +import org.jlab.io.base.DataEvent; +import org.jlab.clas.tracking.trackrep.Helix; +import org.jlab.clas.tracking.kalmanfilter.Units; +import org.jlab.io.hipo.HipoDataSource; +import org.jlab.clas.swimtools.Swim; +import org.jlab.utils.CLASResources; +import cnuphys.magfield.MagneticFields; +import org.jlab.rec.atof.constants.Parameters; + +/** + * The {@code TrackProjector} class projects ahdc tracks to the inner surfaces + * of the bar and wedges of the atof + * + *

+ * Uses ahdc track bank information (for now position, momentum) Creates a + * {@link TrackProjection} for each track. + *

+ * + *

+ * TO DO: - replace hardcoded values with database values. - magnetic field ? + * use swimmer tools? - charge ? + *

+ * + * @author pilleux + */ +public class TrackProjector { + + /** + * projections of tracks. + */ + private List Projections; + + /** + * solenoid magnitude + */ + private Double B; + + /** + * Default constructor that initializes the list of projections as new empty + * list and the magnetic field to 5T. + */ + public TrackProjector() { + Projections = new ArrayList(); + B = 5.0; + } + + /** + * Gets the list of track projections. + * + * @return a {@link List} of {@link TrackProjection} objects representing + * the projections. + */ + public List getProjections() { + return Projections; + } + + /** + * Gets the solenoid magnitude + * + * @return solenoid magnitude + */ + public Double getB() { + return B; + } + + /** + * Sets the list of track projections. + * + * @param Projections a {@link List} of {@link TrackProjection}. + */ + public void setProjections(List Projections) { + this.Projections = Projections; + } + + /** + * Sets the solenoid magnitude. + * + * @param B a double. + */ + public void setB(Double B) { + this.B = B; + } + + /** + * Projects the ahdc tracks in the event onto the atof using a {@link Helix} + * model. + * + * @param event the {@link DataEvent} containing track data to be projected. + */ + public void ProjectTracks(DataEvent event) {//, CalibrationConstantsLoader ccdb) { + + Projections.clear(); + + String track_bank_name = "AHDC::Track"; + + if (event == null) { // check if there is an event + //System.out.print(" no event \n"); + } else if (event.hasBank(track_bank_name) == false) { + // check if there are ahdc tracks in the event + //System.out.print("no tracks \n"); + } else { + DataBank bank = event.getBank(track_bank_name); + int nt = bank.rows(); // number of tracks + TrackProjection projection = new TrackProjection(); + DataBank outputBank = event.createBank("AHDC::Projections", nt); + + for (int i = 0; i < nt; i++) { + + double x = bank.getFloat("x", i); + double y = bank.getFloat("y", i); + double z = bank.getFloat("z", i); + double px = bank.getFloat("px", i); + double py = bank.getFloat("py", i); + double pz = bank.getFloat("pz", i); + + int q = 1; //need the charge sign from tracking + + Units units = Units.MM; //can be MM or CM. + + double xb = 0; + double yb = 0; + + //momenta must be in GeV for the helix class + Helix helix = new Helix(x, y, z, px/1000., py/1000., pz/1000., q, B, xb, yb, units); + + //Intersection points with the middle of the bar or wedge + projection.set_BarIntersect(helix.getHelixPointAtR(Parameters.BAR_MIDDLE_RADIUS)); + projection.set_WedgeIntersect(helix.getHelixPointAtR(Parameters.WEDGE_MIDDLE_RADIUS)); + + //Path length to the middle of the bar or wedge + projection.set_BarPathLength((float) Math.abs(helix.getLAtR(Parameters.BAR_INNER_RADIUS))); + projection.set_WedgePathLength((float) Math.abs(helix.getLAtR(Parameters.WEDGE_INNER_RADIUS))); + + //Path length from the inner radius to the middle radius + projection.set_BarInPathLength((float) Math.abs(helix.getLAtR(Parameters.BAR_MIDDLE_RADIUS)) - projection.get_BarPathLength()); + projection.set_WedgeInPathLength((float) Math.abs(helix.getLAtR(Parameters.WEDGE_MIDDLE_RADIUS)) - projection.get_WedgePathLength()); + Projections.add(projection); + fill_out_bank(outputBank, projection, i); + } + event.appendBank(outputBank); + } + } + + /** + * Projects the MC particles onto the atof using a {@link Helix} + * model. + * + * @param event the {@link DataEvent} containing track data to be projected. + */ + public void projectMCTracks(DataEvent event) {//, CalibrationConstantsLoader ccdb) { + + Projections.clear(); + + String track_bank_name = "MC::Particle"; + + if (event == null) { // check if there is an event + //System.out.print(" no event \n"); + } else if (event.hasBank(track_bank_name) == false) { + // check if there are ahdc tracks in the event + //System.out.print("no tracks \n"); + } else { + DataBank bank = event.getBank(track_bank_name); + int nt = bank.rows(); // number of tracks + TrackProjection projection = new TrackProjection(); + DataBank outputBank = event.createBank("AHDC::Projections", nt); + + for (int i = 0; i < nt; i++) { + + double x = bank.getFloat("vx", i); + double y = bank.getFloat("vy", i); + double z = bank.getFloat("vz", i); + double px = bank.getFloat("px", i); + double py = bank.getFloat("py", i); + double pz = bank.getFloat("pz", i); + + //Put everything in MM + + x = x*10; + y = y*10; + z = z*10; + + Units units = Units.MM; + + int q = 1; //need the charge sign from tracking + + double xb = 0; + double yb = 0; + + //momenta must be in GeV for the helix class + Helix helix = new Helix(x, y, z, px, py, pz, q, B, xb, yb, units); + + //Intersection points with the middle of the bar or wedge + + projection.set_BarIntersect(helix.getHelixPointAtR(Parameters.BAR_MIDDLE_RADIUS)); + projection.set_WedgeIntersect(helix.getHelixPointAtR(Parameters.WEDGE_MIDDLE_RADIUS)); + + //Path length to the middle of the bar or wedge + + projection.set_BarPathLength((float) Math.abs(helix.getLAtR(Parameters.BAR_INNER_RADIUS))); + projection.set_WedgePathLength((float) Math.abs(helix.getLAtR(Parameters.WEDGE_INNER_RADIUS))); + + //Path length from the inner radius to the middle radius + + projection.set_BarInPathLength((float) Math.abs(helix.getLAtR(Parameters.BAR_MIDDLE_RADIUS)) - projection.get_BarPathLength()); + projection.set_WedgeInPathLength((float) Math.abs(helix.getLAtR(Parameters.WEDGE_MIDDLE_RADIUS)) - projection.get_WedgePathLength()); + Projections.add(projection); + fill_out_bank(outputBank, projection, i); + } + event.appendBank(outputBank); + } + } + + public static void fill_out_bank(DataBank outputBank, TrackProjection projection, int i) { + outputBank.setFloat("x_at_bar", i, (float) projection.get_BarIntersect().x()); + outputBank.setFloat("y_at_bar", i, (float) projection.get_BarIntersect().y()); + outputBank.setFloat("z_at_bar", i, (float) projection.get_BarIntersect().z()); + outputBank.setFloat("L_at_bar", i, (float) projection.get_BarPathLength()); + outputBank.setFloat("L_in_bar", i, (float) projection.get_BarInPathLength()); + outputBank.setFloat("x_at_wedge", i, (float) projection.get_WedgeIntersect().x()); + outputBank.setFloat("y_at_wedge", i, (float) projection.get_WedgeIntersect().y()); + outputBank.setFloat("z_at_wedge", i, (float) projection.get_WedgeIntersect().z()); + outputBank.setFloat("L_at_wedge", i, (float) projection.get_WedgePathLength()); + outputBank.setFloat("L_in_wedge", i, (float) projection.get_WedgeInPathLength()); + } + + public static void main(String arg[]) { + + //READING MAG FIELD MAP + System.setProperty("CLAS12DIR", "../../"); + + String mapDir = CLASResources.getResourcePath("etc") + "/data/magfield"; + try { + MagneticFields.getInstance().initializeMagneticFields(mapDir, + "Symm_torus_r2501_phi16_z251_24Apr2018.dat", "Symm_solenoid_r601_phi1_z1201_13June2018.dat"); + } catch (Exception e) { + e.printStackTrace(); + } + + float[] b = new float[3]; + Swim swimmer = new Swim(); + swimmer.BfieldLab(0, 0, 0, b); + double B = Math.abs(b[2]); + + //Track Projector Initialisation with B field + TrackProjector projector = new TrackProjector(); + projector.setB(B); + + //Input to be read + String input = "/Users/npilleux/Desktop/alert/atof-reconstruction/coatjava/reconstruction/alert/src/main/java/org/jlab/rec/atof/Hit/update_protons_to_test_with_tracks.hipo"; + HipoDataSource reader = new HipoDataSource(); + reader.open(input); + int event_number = 0; + while (reader.hasEvent()) { + DataEvent event = (DataEvent) reader.getNextEvent(); + event_number++; + + projector.ProjectTracks(event); + event.getBank("AHDC::Projections").show(); + } + + System.out.print("Read " + event_number + " events"); + } +} diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/atof/banks/RecoBankWriter.java b/reconstruction/alert/src/main/java/org/jlab/rec/atof/banks/RecoBankWriter.java new file mode 100644 index 000000000..c87040d17 --- /dev/null +++ b/reconstruction/alert/src/main/java/org/jlab/rec/atof/banks/RecoBankWriter.java @@ -0,0 +1,53 @@ +/* + * Click nbfs://nbhost/SystemFileSystem/Templates/Licenses/license-default.txt to change this license + * Click nbfs://nbhost/SystemFileSystem/Templates/Classes/Main.java to edit this template + */ +package org.jlab.rec.atof.banks; + +import java.util.ArrayList; +import org.jlab.io.base.DataBank; +import org.jlab.io.base.DataEvent; +import org.jlab.rec.atof.hit.AtofHit; + +/** + * + * @authors npilleux,churaman + */ +public class RecoBankWriter { + + // write useful information in the bank + public static DataBank fillAtofHitBank(DataEvent event, ArrayList hitlist) { + + DataBank bank = event.createBank("ATOF::hits", hitlist.size()); + + if (bank == null) { + System.err.println("COULD NOT CREATE A ATOF::Hits BANK!!!!!!"); + return null; + } + + for(int i =0; i< hitlist.size(); i++) { + bank.setShort("id",i, (short)(i+1)); + bank.setInt("sector",i, (int) hitlist.get(i).getSector()); + bank.setInt("layer",i, (int) hitlist.get(i).getLayer()); + bank.setInt("component",i, (int) hitlist.get(i).getComponent()); + //bank.setShort("trkID",i, (short) hitlist.get(i).get_AssociatedTrkId()); + //bank.setShort("clusterid", i, (short) hitlist.get(i).get_AssociatedClusterID()); + bank.setFloat("time",i, (float) hitlist.get(i).getTime()); + bank.setFloat("x",i, (float) (hitlist.get(i).getX())); + bank.setFloat("y",i, (float) (hitlist.get(i).getY())); + bank.setFloat("z",i, (float) (hitlist.get(i).getZ())); + bank.setFloat("energy",i, (float) hitlist.get(i).getEnergy()); + bank.setFloat("inlength",i, (float) (hitlist.get(i).getInpath_length())); + bank.setFloat("pathlength",i, (float) (hitlist.get(i).getPath_length())); + } + return bank; + } + + /** + * @param args the command line arguments + */ + public static void main(String[] args) { + // TODO code application logic here + } + +} diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/atof/constants/Parameters.java b/reconstruction/alert/src/main/java/org/jlab/rec/atof/constants/Parameters.java new file mode 100644 index 000000000..2c15a5200 --- /dev/null +++ b/reconstruction/alert/src/main/java/org/jlab/rec/atof/constants/Parameters.java @@ -0,0 +1,44 @@ +package org.jlab.rec.atof.constants; + +/** + * + * @authors npilleux,churaman + */ +public class Parameters { + + Parameters() { + } + + //In millimiters + public static final double BAR_INNER_RADIUS = 77; + public static final double WEDGE_INNER_RADIUS = 80; + public static final double BAR_THICKNESS = 3; + public static final double WEDGE_THICKNESS = 20; + public static final double BAR_MIDDLE_RADIUS = BAR_INNER_RADIUS + BAR_THICKNESS / 2; + public static final double WEDGE_MIDDLE_RADIUS = WEDGE_INNER_RADIUS + WEDGE_THICKNESS / 2; + + + public static final double VEFF = 200.0;//mm/ns + public static final double TDC2TIME = 0.015625;//ns per channel bin + public static final double ATT_L = 1600.0;//mm + public static final double TOT2ENERGY_BAR = 1.956 * 0.3 /1000;//to MeV + public static final double TOT2ENERGY_WEDGE = 1.956 * 2.0 /1000;//to MeV + + public static double LENGTH_ATOF = 279.7; //detector length in mm + + public static double SIGMA_PHI_TRACK_MATCHING_BAR = 180 * Math.PI/180.;//in rad + public static double SIGMA_PHI_TRACK_MATCHING_WEDGE = 180 * Math.PI/180.;//in rad + + public static double SIGMA_Z_TRACK_MATCHING_BAR = 200;//in mm + public static double SIGMA_Z_TRACK_MATCHING_WEDGE = 200;//in mm + + public static double SIGMA_PHI_CLUSTERING = 6;//in deg + + /** + * @param args the command line arguments + */ + public static void main(String[] args) { + // TODO code application logic here + } + +}