diff --git a/common-tools/clas-detector/src/main/java/org/jlab/detector/base/DetectorType.java b/common-tools/clas-detector/src/main/java/org/jlab/detector/base/DetectorType.java index ff90cc4a11..dfebf090e0 100644 --- a/common-tools/clas-detector/src/main/java/org/jlab/detector/base/DetectorType.java +++ b/common-tools/clas-detector/src/main/java/org/jlab/detector/base/DetectorType.java @@ -34,6 +34,10 @@ public enum DetectorType { AHDC (24, "AHDC"), ATOF (25, "ATOF"), RECOIL (26, "RECOIL"), + MUCAL (28, "MUCAL"), + MUVT (29, "MUVT"), + MURT (30, "MURT"), + MURH (31, "MURH"), TARGET (100, "TARGET"), MAGNETS (101, "MAGNETS"); 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 85ba18104d..ab8564fb12 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 @@ -579,6 +579,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; @@ -823,19 +824,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); @@ -852,6 +854,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); @@ -864,7 +877,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); @@ -879,7 +892,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); } } diff --git a/etc/bankdefs/hipo4/data.json b/etc/bankdefs/hipo4/data.json index c76de0015b..83383bf30c 100644 --- a/etc/bankdefs/hipo4/data.json +++ b/etc/bankdefs/hipo4/data.json @@ -396,6 +396,66 @@ { "name":"TDC" , "type":"I", "info":"TDC value"} ] }, + { + "name" : "MUCAL::adc", + "group": 22800, + "item" : 11, + "info": "ADC bank for the muCLAS12 Electromagnetic Calorimeter", + "entries":[ + { "name":"sector" , "type":"B", "info":"sector (=1)"}, + { "name":"layer" , "type":"B", "info":"layer (1-6)"}, + { "name":"component" , "type":"S", "info":"strip"}, + { "name":"order" , "type":"B", "info":"order: 0 - ADCL , 1 - ADCR"}, + { "name":"ADC" , "type":"I", "info":"ADC maximum"}, + { "name":"time" , "type":"F", "info":"time from pulse fit"}, + { "name":"ped" , "type":"S", "info":"pedestal from pulse analysis"} + ] + }, + { + "name" : "MUVT::adc", + "group": 22900, + "item" : 11, + "info": "ADC bank for the muCLAS12 Vertex Tracker", + "entries":[ + {"name":"sector" , "type":"B", "info":"sector (=1)"}, + {"name":"layer" , "type":"B", "info":"layer (=1)"}, + {"name":"component" , "type":"S", "info":"crystal id"}, + { "name":"order" , "type":"B", "info":"order: 0 - ADCL , 1 - ADCR"}, + { "name":"ADC" , "type":"I", "info":"ADC integral from pulse fit"}, + { "name":"time" , "type":"F", "info":"time from pulse fit"}, + { "name":"ped" , "type":"S", "info":"pedestal from pulse analysis"} + ] + }, + { + "name" : "MURT::adc", + "group": 23000, + "item" : 11, + "info": "ADC bank for the muCLAS12 Recoil Tracker", + "entries":[ + {"name":"sector" , "type":"B", "info":"sector (=1)"}, + {"name":"layer" , "type":"B", "info":"layer (=1)"}, + {"name":"component" , "type":"S", "info":"crystal id"}, + { "name":"order" , "type":"B", "info":"order: 0 - ADCL , 1 - ADCR"}, + { "name":"ADC" , "type":"I", "info":"ADC integral from pulse fit"}, + { "name":"time" , "type":"F", "info":"time from pulse fit"}, + { "name":"ped" , "type":"S", "info":"pedestal from pulse analysis"} + ] + }, + { + "name" : "MURH::adc", + "group": 23100, + "item" : 11, + "info": "ADC bank for the muCLAS12 Recoil Hodoscope", + "entries":[ + {"name":"sector" , "type":"B", "info":"sector (1-8)"}, + {"name":"layer" , "type":"B", "info":"layer (1-2)"}, + {"name":"component" , "type":"S", "info":"tile id"}, + { "name":"order" , "type":"B", "info":"order: 0 - ADCL , 1 - ADCR"}, + { "name":"ADC" , "type":"I", "info":"ADC integral from pulse fit"}, + { "name":"time" , "type":"F", "info":"time from pulse fit"}, + { "name":"ped" , "type":"S", "info":"pedestal from pulse analysis"} + ] + }, { "name" : "RF::adc", "group": 21700, diff --git a/etc/bankdefs/hipo4/mu.json b/etc/bankdefs/hipo4/mu.json new file mode 100644 index 0000000000..6d171f0207 --- /dev/null +++ b/etc/bankdefs/hipo4/mu.json @@ -0,0 +1,176 @@ +[ + { + "name": "MUCAL::hits", + "group": 22800, + "item" : 21, + "info": "Reconstructed Hits in muCLAS12 calorimeter", + "entries": [ + {"name":"idx", "type":"B", "info":"x id of grid"}, + {"name":"idy", "type":"B", "info":"y id of grid"}, + {"name":"x", "type":"F", "info":"Hit X position (cm)" }, + {"name":"y", "type":"F", "info":"Hit Y position (cm)" }, + {"name":"z", "type":"F", "info":"Hit Z position (cm)" }, + {"name":"energy", "type":"F", "info":"Hit Energy" }, + {"name":"time", "type":"F", "info":"Hit Time" }, + {"name":"hitID", "type":"S", "info":"Hit Pointer to ADC bank"}, + {"name":"clusterID", "type":"S", "info":"Hit Pointer to Cluster Bank"} + ] + }, + { + "name": "MUCAL::clusters", + "group": 22800, + "item" : 22, + "info": "Reconstructed Clusters in muCLAS12 calorimeter", + "entries": [ + {"name":"size", "type":"S", "info":"Cluster size" }, + {"name":"id", "type":"S", "info":"Cluster ID" }, + {"name":"x", "type":"F", "info":"Cluster centroid X moment (cm)" }, + {"name":"y", "type":"F", "info":"Cluster centroid Y moment (cm)" }, + {"name":"z", "type":"F", "info":"Cluster centroid Z moment (cm)" }, + {"name":"widthX", "type":"F", "info":"Cluster width in x (cm)" }, + {"name":"widthY", "type":"F", "info":"Cluster width in y (cm)" }, + {"name":"radius", "type":"F", "info":"Cluster radius (cm)" }, + {"name":"time", "type":"F", "info":"Cluster timing information" }, + {"name":"energy", "type":"F", "info":"Cluster total energy" }, + {"name":"maxEnergy", "type":"F", "info":"Cluster maximum energy" }, + {"name":"recEnergy", "type":"F", "info":"Cluster reconstructed energy" } + ] + }, + { + "name": "MUVT::hits", + "group": 22900, + "item" : 21, + "info": "MUVT hits", + "entries": [ + {"name":"id", "type":"S", "info":"id of the hit"}, + {"name":"sector", "type":"B", "info":"sector number"}, + {"name":"layer", "type":"B", "info":"layer number"}, + {"name":"strip", "type":"S", "info":"strip number"}, + {"name":"energy", "type":"F", "info":"energy of the hit (eV)"}, + {"name":"time", "type":"F", "info":"time of the hit (ns)"}, + {"name":"clusterId", "type":"S", "info":"id of the cluster the hit belongs to"}, + {"name":"status", "type":"S", "info":"status of the hit"} + ] + }, + { + "name": "MUVT::clusters", + "group": 22900, + "item" : 22, + "info": "reconstructed clusters from MUVT", + "entries": [ + {"name":"id", "type":"S", "info":"id of the cluster"}, + {"name":"sector", "type":"B", "info":"sector"}, + {"name":"layer", "type":"B", "info":"layer"}, + {"name":"strip", "type":"S", "info":"seed strip"}, + {"name":"energy", "type":"F", "info":"energy of the cluster (eV)"}, + {"name":"time", "type":"F", "info":"time of the cluster (ns)"}, + {"name":"xo", "type":"F", "info":"strip origin X coordinate (cm)"}, + {"name":"yo", "type":"F", "info":"strip origin Y coordinate (cm)"}, + {"name":"zo", "type":"F", "info":"strip origin Z coordinate (cm)"}, + {"name":"xe", "type":"F", "info":"strip end X coordinate (cm)"}, + {"name":"ye", "type":"F", "info":"strip end Y coordinate (cm)"}, + {"name":"ze", "type":"F", "info":"strip end Z coordinate (cm)"}, + {"name":"size", "type":"S", "info":"size of the cluster"}, + {"name":"status", "type":"S", "info":"status of the cluster"} + ] + }, + { + "name": "MUVT::crosses", + "group": 22900, + "item" : 23, + "info": "reconstructed crosses from MUVT", + "entries": [ + {"name":"id", "type":"S", "info":"id of the cross"}, + {"name":"sector", "type":"B", "info":"sector"}, + {"name":"region", "type":"B", "info":"region"}, + {"name":"energy", "type":"F", "info":"energy of the cross (eV)"}, + {"name":"time", "type":"F", "info":"time of the cross (ns)"}, + {"name":"x", "type":"F", "info":"x coordinate (cm)"}, + {"name":"y", "type":"F", "info":"y coordinate (cm)"}, + {"name":"z", "type":"F", "info":"z coordinate (cm)"}, + {"name":"cluster1", "type":"S", "info":"id of the cluster in the V layer"}, + {"name":"cluster2", "type":"S", "info":"id of the cluster in the W layer"}, + {"name":"status", "type":"S", "info":"status of the cluster"} + ] + }, + { + "name": "MUVT::tracks", + "group": 22900, + "item" : 36, + "info": "reconstructed tracks using MUVT information", + "entries": [ + {"name":"index", "type":"S", "info":"index of the track in the DC bank"}, + {"name":"status", "type":"B", "info":"status of the track (0: refitted using FMT, 1: original DC track)"}, + {"name":"sector", "type":"B", "info":"sector of the track in DC"}, + {"name":"vx", "type":"F", "info":"Vertex x-position of the swam track to the DOCA to the beamline (in cm)"}, + {"name":"vy", "type":"F", "info":"Vertex y-position of the swam track to the DOCA to the beamline (in cm)"}, + {"name":"vz", "type":"F", "info":"Vertex z-position of the swam track to the DOCA to the beamline (in cm)"}, + {"name":"px", "type":"F", "info":"3-momentum x-coordinate of the swam track to the DOCA to the beamline (in GeV)"}, + {"name":"py", "type":"F", "info":"3-momentum y-coordinate of the swam track to the DOCA to the beamline (in GeV)"}, + {"name":"pz", "type":"F", "info":"3-momentum z-coordinate of the swam track to the DOCA to the beamline (in GeV)"}, + {"name":"charge", "type":"B", "info":"charge of the track"}, + {"name":"chi2", "type":"F", "info":"chi^2 of the fit"}, + {"name":"NDF", "type":"B", "info":"number of degrees of freedom of the fit"} + ] + }, + { + "name": "MUVT::trajectory", + "group": 22900, + "item" : 37, + "info": "MUVT tracks trajectory bank", + "entries": [ + {"name":"index", "type":"S", "info":"index of the track in the DC bank"}, + {"name":"detector", "type":"B", "info":"id of the detector"}, + {"name":"layer", "type":"B", "info":"id of the layer"}, + {"name":"x", "type":"F", "info":"track x position at detector surface (cm)"}, + {"name":"y", "type":"F", "info":"track y position at detector surface (cm)"}, + {"name":"z", "type":"F", "info":"track z position at detector surface (cm)"}, + {"name":"tx", "type":"F", "info":"track unit direction vector x component at detector surface"}, + {"name":"ty", "type":"F", "info":"track unit direction vector y component at detector surface"}, + {"name":"tz", "type":"F", "info":"track unit direction vector z component at detector surface"}, + {"name":"lx", "type":"F", "info":"track x position in local coordinates (cm)"}, + {"name":"ly", "type":"F", "info":"track y position in local coordinates (cm)"}, + {"name":"lz", "type":"F", "info":"track z position in local coordinates (cm)"}, + {"name":"dx", "type":"F", "info":"DC track x position in local coordinates (cm)"}, + {"name":"dy", "type":"F", "info":"DC track y position in local coordinates (cm)"}, + {"name":"dz", "type":"F", "info":"DC track z position in local coordinates (cm)"}, + {"name":"path", "type":"F", "info":"pathlength of the track from the track vertex to the detector surface (cm)"} + ] + }, + { + "name": "MURH::hits", + "group": 23100, + "item" : 21, + "info": "Reconstructed Hits in muCLAS12 hodoscope", + "entries": [ + {"name":"sector", "type":"B", "info":"sector number"}, + {"name":"layer", "type":"B", "info":"layer number"}, + {"name":"component", "type":"S", "info":"component number"}, + {"name":"x", "type":"F", "info":"Hit X position (cm)" }, + {"name":"y", "type":"F", "info":"Hit Y position (cm)" }, + {"name":"z", "type":"F", "info":"Hit Z position (cm)" }, + {"name":"energy", "type":"F", "info":"Hit Energy" }, + {"name":"time", "type":"F", "info":"Hit Time" }, + {"name":"hitID", "type":"S", "info":"Hit Pointer to ADC bank"}, + {"name":"clusterID", "type":"S", "info":"Hit Pointer to Cluster Bank"} + ] + }, + { + "name": "MURH::clusters", + "group": 23100, + "item" : 22, + "info": "Reconstructed clusters in muCLAS12 hodoscope", + "entries": [ + {"name":"size", "type":"S", "info":"Cluster size" }, + {"name":"id", "type":"S", "info":"Cluster ID" }, + {"name":"x", "type":"F", "info":"Cluster centroid X moment (cm)" }, + {"name":"y", "type":"F", "info":"Cluster centroid Y moment (cm)" }, + {"name":"z", "type":"F", "info":"Cluster centroid Z moment (cm)" }, + {"name":"widthX", "type":"F", "info":"Cluster width in x (cm)" }, + {"name":"widthY", "type":"F", "info":"Cluster width in y (cm)" }, + {"name":"radius", "type":"F", "info":"Cluster radius (cm)" }, + {"name":"time", "type":"F", "info":"Cluster timing information" }, + {"name":"energy", "type":"F", "info":"Cluster total energy" } + ] + } +] diff --git a/reconstruction/ft/src/main/java/org/jlab/rec/ft/cal/FTCALEngine.java b/reconstruction/ft/src/main/java/org/jlab/rec/ft/cal/FTCALEngine.java index eb9bc1cdfa..454fc11b83 100644 --- a/reconstruction/ft/src/main/java/org/jlab/rec/ft/cal/FTCALEngine.java +++ b/reconstruction/ft/src/main/java/org/jlab/rec/ft/cal/FTCALEngine.java @@ -6,6 +6,7 @@ import javax.swing.JFrame; import org.jlab.clas.detector.DetectorData; import org.jlab.clas.detector.DetectorEvent; +import org.jlab.clas.physics.PhysicsEvent; import org.jlab.clas.reco.ReconstructionEngine; import org.jlab.geom.prim.Vector3D; import org.jlab.groot.data.H1F; @@ -89,4 +90,75 @@ public int setRunConditionsParameters(DataEvent event) { return run; } + + + public static void main (String arg[]) { + FTCALEngine cal = new FTCALEngine(); + cal.init(); + // String input = "/Users/devita/Work/clas12/simulations/clas12Tags/4.4.0/out.hipo"; +// String input = "/home/filippi/clas/ForwardTracker/DATA/out_realGeo_noMagField.data"; +// String input = "/home/filippi/clas/gemc/electronGun/gemc.hipo"; + String input = "/home/filippi/clas12/fttrkDev/clas12-offline-software-6.5.13-fttrkDev/ft_005038.evio.01231.hipo"; + HipoDataSource reader = new HipoDataSource(); +// String input = "/Users/devita/Work/clas12/simulations/tests/detectors/clas12/ft/out_header.ev"; +// EvioSource reader = new EvioSource(); + reader.open(input); + + // initialize histos + H1F h1 = new H1F("Cluster Energy",100, 0.,5.); + h1.setOptStat(Integer.parseInt("1111")); h1.setTitleX("Cluster Energy (GeV)"); + H1F h2 = new H1F("Energy Resolution",100, -1, 1); + h2.setOptStat(Integer.parseInt("1111")); h2.setTitleX("Energy Resolution(GeV)"); + H1F h3 = new H1F("Theta Resolution",100, -2, 2); + h3.setOptStat(Integer.parseInt("1111")); h3.setTitleX("Theta Resolution(deg)"); + H1F h4 = new H1F("Phi Resolution",100, -10, 10); + h4.setOptStat(Integer.parseInt("1111")); h4.setTitleX("Phi Resolution(deg)"); +// H1F h5 = new H1F("Time Resolution",100, -10, 10); + H1F h5 = new H1F("Time Resolution",100, -100, 300); + h5.setOptStat(Integer.parseInt("1111")); h5.setTitleX("Time Resolution(ns)"); + H2F h6 = new H2F("cluster xy", 100, -15., 15., 100, -15., 15.); + h6.setTitleX("cluster x"); h6.setTitleY("cluster y"); + + while(reader.hasEvent()){ +// for(int nev=0; nev<2; nev++){ + DataEvent event = (DataEvent) reader.getNextEvent(); + cal.processDataEvent(event); + + DetectorEvent detectorEvent = DetectorData.readDetectorEvent(event); + PhysicsEvent gen = detectorEvent.getGeneratedEvent(); + if(event.hasBank("FTCAL::clusters")) { + DataBank bank = event.getBank("FTCAL::clusters"); + int nrows = bank.rows(); + for(int i=0; i + + 4.0.0 + + org.jlab.clas12.detector + clas12detector-mu + 13.5.0-SNAPSHOT + jar + + + org.jlab.clas12 + reconstruction + 13.5.0-SNAPSHOT + + + + + org.jlab.jnp + jnp-hipo + + + org.jlab + groot + + + org.jlab.clas + swim-tools + 13.5.0-SNAPSHOT + + + org.jlab.clas + clas-reco + 13.5.0-SNAPSHOT + + + org.jlab.clas + clas-utils + 13.5.0-SNAPSHOT + + + org.jlab.clas + clas-detector + 13.5.0-SNAPSHOT + + + org.jlab.clas + clas-io + 13.5.0-SNAPSHOT + + + org.jlab.clas + clas-physics + 13.5.0-SNAPSHOT + + + org.jlab.clas + clas-geometry + 13.5.0-SNAPSHOT + + + org.jlab.clas + clas-jcsg + 13.5.0-SNAPSHOT + + + + diff --git a/reconstruction/mu/src/main/java/org/jlab/rec/mucal/MUCALCluster.java b/reconstruction/mu/src/main/java/org/jlab/rec/mucal/MUCALCluster.java new file mode 100644 index 0000000000..2741127d21 --- /dev/null +++ b/reconstruction/mu/src/main/java/org/jlab/rec/mucal/MUCALCluster.java @@ -0,0 +1,234 @@ +package org.jlab.rec.mucal; + +import java.util.ArrayList; +import org.jlab.geom.prim.Point3D; +import org.jlab.utils.groups.IndexedTable; + + + +public class MUCALCluster extends ArrayList { + + /** + * A cluster in the calorimeter consists of an array of crystals that are grouped together according to + * the algorithm of the CalClusterFinder class + */ + private static final long serialVersionUID = 1L; + + private int _clusID; // cluster ID + private boolean _clusStat=true; // cluster status flag (true==good, false==bad) + + // constructor + public MUCALCluster(int cid) { + this.setID(cid); + } + + public int getID() { + return _clusID; + } + + public void setID(int _clusID) { + this._clusID = _clusID; + } + + public int getSize() { + // return number of crystals in cluster + return this.size(); + } + + public double getEnergy() { + // return measured energy + double clusterEnergy = 0; + for(int i=0; i clusterTable.getDoubleValue("cluster_min_size", 1,1,0) && + this.getEnergy()> clusterTable.getDoubleValue("cluster_min_energy", 1,1,0)/1000) { + this._clusStat=true; + } + else { + this._clusStat=false; + } + } + + public boolean getStatus() { + return this._clusStat; + } + + private double weight(MUCALHit hit, double clusEnergy) { + return Math.max(0., (3.45+Math.log(hit.get_Edep()/clusEnergy))); + } + + public boolean containsHit(MUCALHit hit, IndexedTable thresholds, IndexedTable clusterTable) { + boolean addFlag = false; + if(hit.get_Edep()>thresholds.getDoubleValue("thresholdCluster",1,1,MUCALHit.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()); + double yDiff = Math.abs(hit.get_IDY() - this.get(j).get_IDY()); + if(tDiff <= clusterTable.getDoubleValue("time_window", 1,1,0) && xDiff <= 1 && yDiff <= 1 && (xDiff + yDiff) >0) addFlag = true; + } + } + return addFlag; + } + + @Override + public String toString(){ + StringBuilder str = new StringBuilder(); + + str.append(String.format("Cluster: %4d\n", this.getID())); + str.append(String.format("\tStatus: %b", this.getStatus())); + str.append(String.format("\tSize: %4d", this.getSize())); + str.append(String.format("\tE: %7.3f", this.getEnergy())); + str.append(String.format("\tTime: %7.3f", this.getTime())); + str.append(String.format("\tTheta: %7.3f", this.getTheta())); + str.append(String.format("\tPhi: %7.3f\n", this.getPhi())); + for(int j = 0; j< this.size(); j++) { + str.append(String.format("\thit #%d\t%d\t%d\t%7.3f\n", j, this.get(j).get_IDX(), this.get(j).get_IDY(), this.get(j).get_Edep())); + } + return str.toString(); + } + + public void show() { + System.out.println(this.toString()); + } +} + + + diff --git a/reconstruction/mu/src/main/java/org/jlab/rec/mucal/MUCALConstants.java b/reconstruction/mu/src/main/java/org/jlab/rec/mucal/MUCALConstants.java new file mode 100644 index 0000000000..5917c9770b --- /dev/null +++ b/reconstruction/mu/src/main/java/org/jlab/rec/mucal/MUCALConstants.java @@ -0,0 +1,21 @@ +package org.jlab.rec.mucal; + +public class MUCALConstants { + + public MUCALConstants() { + } + + public static int DEBUGMODE = 0; + + // RECONSTRUCTION CONSTANTS + public static final double TIMECONVFAC = 100./4.; // conversion factor from TDC channel to time (ns^-1) + public static final double VEFF = 150.; // speed of light in the scintillator mm/ns + + // energy threshold in GeV + // GEOMETRY PARAMETERS + public static double CRYS_DELTA = 11.5; + public static double CRYS_WIDTH = 15.3; // crystal width in mm + public static double CRYS_LENGTH = 200.; // crystal length in mm + public static double CRYS_ZPOS = 1898.; // position of the crystal front face + +} diff --git a/reconstruction/mu/src/main/java/org/jlab/rec/mucal/MUCALEngine.java b/reconstruction/mu/src/main/java/org/jlab/rec/mucal/MUCALEngine.java new file mode 100644 index 0000000000..0d9dc58efb --- /dev/null +++ b/reconstruction/mu/src/main/java/org/jlab/rec/mucal/MUCALEngine.java @@ -0,0 +1,164 @@ +package org.jlab.rec.mucal; + +import java.util.ArrayList; +import java.util.Arrays; +import java.util.List; +import javax.swing.JFrame; +import org.jlab.clas.detector.DetectorData; +import org.jlab.clas.detector.DetectorEvent; +import org.jlab.clas.physics.PhysicsEvent; +import org.jlab.clas.reco.ReconstructionEngine; +import org.jlab.geom.prim.Vector3D; +import org.jlab.groot.data.H1F; +import org.jlab.groot.data.H2F; +import org.jlab.groot.graphics.EmbeddedCanvas; +import org.jlab.io.base.DataBank; +import org.jlab.io.base.DataEvent; +import org.jlab.io.evio.EvioDataBank; +import org.jlab.io.evio.EvioDataEvent; +import org.jlab.io.hipo.HipoDataSource; + + +public class MUCALEngine extends ReconstructionEngine { + + public MUCALEngine() { + super("MUCAL", "devita", "3.0"); + } + + MUCALReconstruction reco; + + @Override + public boolean init() { + reco = new MUCALReconstruction(); + reco.debugMode=0; + + String[] tables = new String[]{ + "/calibration/ft/ftcal/charge_to_energy", + "/calibration/ft/ftcal/time_offsets", + "/calibration/ft/ftcal/time_walk", + "/calibration/ft/ftcal/status", + "/calibration/ft/ftcal/thresholds", + "/calibration/ft/ftcal/cluster", + "/calibration/ft/ftcal/energycorr" + }; + requireConstants(Arrays.asList(tables)); + this.getConstantsManager().setVariation("default"); + + this.registerOutputBank("MUCAL::hits","MUCAL::clusters"); + + return true; + } + + @Override + public boolean processDataEvent(DataEvent event) { + List allHits = new ArrayList(); + List selectedHits = new ArrayList(); + List clusters = new ArrayList(); + + // update calibration constants based on run number if changed + int run = setRunConditionsParameters(event); + + if(run>=0) { + // get hits fron banks + allHits = reco.initMUCAL(event,this.getConstantsManager(), run); + // select good hits and order them by energy + selectedHits = reco.selectHits(allHits,this.getConstantsManager(), run); + // create clusters + clusters = reco.findClusters(selectedHits, this.getConstantsManager(), run); + // set cluster status + reco.selectClusters(clusters, this.getConstantsManager(), run); + // write output banks + reco.writeBanks(event, selectedHits, clusters, this.getConstantsManager(), run); + } + return true; + } + + public int setRunConditionsParameters(DataEvent event) { + int run = -1; + if(event.hasBank("RUN::config")==false) { + System.out.println("RUN CONDITIONS NOT READ!"); + } + + if(event instanceof EvioDataEvent) { + EvioDataBank bank = (EvioDataBank) event.getBank("RUN::config"); + run = bank.getInt("Run",0); + } + else { + DataBank bank = event.getBank("RUN::config"); + run = bank.getInt("run",0); + } + + return run; + } + + + public static void main (String arg[]) { + MUCALEngine cal = new MUCALEngine(); + cal.init(); + // String input = "/Users/devita/Work/clas12/simulations/clas12Tags/4.4.0/out.hipo"; +// String input = "/home/filippi/clas/ForwardTracker/DATA/out_realGeo_noMagField.data"; +// String input = "/home/filippi/clas/gemc/electronGun/gemc.hipo"; + String input = "/home/filippi/clas12/fttrkDev/clas12-offline-software-6.5.13-fttrkDev/ft_005038.evio.01231.hipo"; + HipoDataSource reader = new HipoDataSource(); +// String input = "/Users/devita/Work/clas12/simulations/tests/detectors/clas12/ft/out_header.ev"; +// EvioSource reader = new EvioSource(); + reader.open(input); + + // initialize histos + H1F h1 = new H1F("Cluster Energy",100, 0.,5.); + h1.setOptStat(Integer.parseInt("1111")); h1.setTitleX("Cluster Energy (GeV)"); + H1F h2 = new H1F("Energy Resolution",100, -1, 1); + h2.setOptStat(Integer.parseInt("1111")); h2.setTitleX("Energy Resolution(GeV)"); + H1F h3 = new H1F("Theta Resolution",100, -2, 2); + h3.setOptStat(Integer.parseInt("1111")); h3.setTitleX("Theta Resolution(deg)"); + H1F h4 = new H1F("Phi Resolution",100, -10, 10); + h4.setOptStat(Integer.parseInt("1111")); h4.setTitleX("Phi Resolution(deg)"); +// H1F h5 = new H1F("Time Resolution",100, -10, 10); + H1F h5 = new H1F("Time Resolution",100, -100, 300); + h5.setOptStat(Integer.parseInt("1111")); h5.setTitleX("Time Resolution(ns)"); + H2F h6 = new H2F("cluster xy", 100, -15., 15., 100, -15., 15.); + h6.setTitleX("cluster x"); h6.setTitleY("cluster y"); + + while(reader.hasEvent()){ +// for(int nev=0; nev<2; nev++){ + DataEvent event = (DataEvent) reader.getNextEvent(); + cal.processDataEvent(event); + + DetectorEvent detectorEvent = DetectorData.readDetectorEvent(event); + PhysicsEvent gen = detectorEvent.getGeneratedEvent(); + if(event.hasBank("MUCAL::clusters")) { + DataBank bank = event.getBank("MUCAL::clusters"); + int nrows = bank.rows(); + for(int i=0; i{ + // class implements Comparable interface to allow for sorting a collection of hits by Edep values + public static final int REFCOMPONENT =245; + + // constructor + public MUCALHit(int i, int ICOMPONENT, int ADC, int TDC, IndexedTable charge2Energy, IndexedTable timeOffsets, IndexedTable timeWalk, IndexedTable cluster) { + this._COMPONENT = ICOMPONENT; + 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,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,REFCOMPONENT)*Math.exp(-this._Charge*timeWalk.getDoubleValue("lambda", 1,1,REFCOMPONENT)); + } + this.set_Time(((double) this._TDC)/MUCALConstants.TIMECONVFAC + -(MUCALConstants.CRYS_LENGTH-cluster.getDoubleValue("depth_z", 1,1,0))/MUCALConstants.VEFF + -timeOffsets.getDoubleValue("time_offset", 1,1,REFCOMPONENT)-twCorr); +// if(this.get_Edep()>0.1) System.out.println(ICOMPONENT + " " + this._TDC + " " + +// MUCALConstantsLoader.TIMECONVFAC + " " + MUCALConstantsLoader.time_offset[0][0][ICOMPONENT-1] + " " + +// this.get_Time()); + this.set_Dx( (this._IDX-MUCALConstants.CRYS_DELTA )* MUCALConstants.CRYS_WIDTH); + this.set_Dy( (this._IDY-MUCALConstants.CRYS_DELTA )* MUCALConstants.CRYS_WIDTH); + this.set_Dz(MUCALConstants.CRYS_ZPOS+cluster.getDoubleValue("depth_z", 1,1,0)); + this.set_DGTZIndex(i); + this.set_ClusID(0); + } + + public MUCALHit(int i, int ICOMPONENT, int ADC, float time, IndexedTable charge2Energy, IndexedTable timeOffsets, IndexedTable timeWalk, IndexedTable cluster) { + this._COMPONENT = ICOMPONENT; + 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,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,REFCOMPONENT)*Math.exp(-this._Charge*timeWalk.getDoubleValue("lambda", 1,1,REFCOMPONENT)); + } + + this.set_Time(time -(MUCALConstants.CRYS_LENGTH-cluster.getDoubleValue("depth_z", 1,1,0))/MUCALConstants.VEFF + -timeOffsets.getDoubleValue("time_offset", 1,1,REFCOMPONENT)-twCorr); +// if(this.get_Edep()>0.1) System.out.println(ICOMPONENT + " " + this._TDC + " " + +// MUCALConstantsLoader.TIMECONVFAC + " " + MUCALConstantsLoader.time_offset[0][0][ICOMPONENT-1] + " " + +// this.get_Time()); + this.set_Dx( (this._IDX-MUCALConstants.CRYS_DELTA )* MUCALConstants.CRYS_WIDTH); + this.set_Dy( (this._IDY-MUCALConstants.CRYS_DELTA )* MUCALConstants.CRYS_WIDTH); + this.set_Dz(MUCALConstants.CRYS_ZPOS+cluster.getDoubleValue("depth_z", 1,1,0)); + this.set_DGTZIndex(i); + this.set_ClusID(0); + } + + private int _COMPONENT; // Component number + private int _IDX; // Crystal ID: X + private int _IDY; // Crystal ID: Y + private int _ADC; // ADC + private int _TDC; // TDC + + private double _Charge; // Reconstructed energy deposited by the hit in the crystal + private double _Edep; // Reconstructed energy deposited by the hit in the crystal + private double _Time; // Reconstructed time, for now it is the gemc time + private double _Dx; + private double _Dy; + private double _Dz; + private int _DGTZIndex; // Pointer to cluster + private int _ClusIndex; // Pointer to cluster + + + public int get_COMPONENT() { + return _COMPONENT; + } + + + + public void set_COMPONENT(int COMPONENT) { + this._COMPONENT = COMPONENT; + } + + + public int get_IDX() { + return _IDX; + } + + + + public void set_IDX(int IDX) { + this._IDX = IDX; + } + + + + public int get_IDY() { + return _IDY; + } + + + + public void set_IDY(int IDY) { + this._IDY = IDY; + } + + + + public int get_ADC() { + return _ADC; + } + + + + public void set_ADC(int ADC) { + this._ADC = ADC; + } + + + + public int get_TDC() { + return _TDC; + } + + + + public void set_TDC(int TDC) { + this._TDC = TDC; + } + + + public double get_Edep() { + return _Edep; + } + + + public final void set_Edep(double edep) { + this._Edep = edep; + } + + + + public double get_Time() { + return _Time; + } + + + public final void set_Time(double Time) { + this._Time = Time; + } + + + public double get_Dx() { + return _Dx; + } + + + public final void set_Dx(double Dx) { + this._Dx = Dx; + } + + + public double get_Dy() { + return _Dy; + } + + + public final void set_Dy(double Dy) { + this._Dy = Dy; + } + + + public double get_Dz() { + return _Dz; + } + + + public final void set_Dz(double Dz) { + this._Dz = Dz; + } + + + public int get_DGTZIndex() { + return _DGTZIndex; + } + + + public final void set_DGTZIndex(int _DGTZIndex) { + this._DGTZIndex = _DGTZIndex; + } + + + public int get_ClusID() { + return _ClusIndex; + } + + + public final void set_ClusID(int _ClusIndex) { + this._ClusIndex = _ClusIndex; + } + + public static boolean passHitSelection(MUCALHit hit, IndexedTable thresholds) { + // a selection cut to pass the hit. + if(hit.get_Edep() > thresholds.getDoubleValue("thresholdHit", 1,1,REFCOMPONENT)) { + return true; + } else { + return false; + } + } + + @Override + public int compareTo(MUCALHit arg0) { + if(this.get_Edep() initMUCAL(DataEvent event, ConstantsManager manager, int run) { + + IndexedTable charge2Energy = manager.getConstants(run, "/calibration/ft/ftcal/charge_to_energy"); + IndexedTable timeOffsets = manager.getConstants(run, "/calibration/ft/ftcal/time_offsets"); + IndexedTable timeWalk = manager.getConstants(run, "/calibration/ft/ftcal/time_walk"); + IndexedTable cluster = manager.getConstants(run, "/calibration/ft/ftcal/cluster"); + IndexedTable status = manager.getConstants(run, "/calibration/ft/ftcal/status"); + + if(this.debugMode>=1) System.out.println("\nAnalyzing new event"); + + List allhits = this.readRawHits(event,charge2Energy,timeOffsets,timeWalk,cluster,status); + + if(debugMode>=1) { + System.out.println("Found " + allhits.size() + " hits"); + for(int i = 0; i < allhits.size(); i++) { + System.out.print(i + "\t"); + allhits.get(i).show(); + } + } + return allhits; + } + + public List selectHits(List allhits, ConstantsManager manager, int run) { + + if(debugMode>=1) System.out.println("\nSelecting hits"); + ArrayList hits = new ArrayList<>(); + + IndexedTable thresholds = manager.getConstants(run, "/calibration/ft/ftcal/thresholds"); + + for(int i = 0; i < allhits.size(); i++) + { + if(MUCALHit.passHitSelection(allhits.get(i), thresholds)) { + hits.add(allhits.get(i)); + } + } + Collections.sort(hits); + if(debugMode>=1) { + System.out.println("List of selected hits"); + for(int i = 0; i < hits.size(); i++) + { + System.out.print(i + "\t"); + hits.get(i).show(); + } + } + return hits; + } + + public List findClusters(List hits, ConstantsManager manager, int run) { + + List clusters = new ArrayList(); + + IndexedTable thresholds = manager.getConstants(run, "/calibration/ft/ftcal/thresholds"); + IndexedTable clusterTable = manager.getConstants(run, "/calibration/ft/ftcal/cluster"); + + if(debugMode>=1) System.out.println("\nBuilding clusters"); + for(int ihit=0; ihit=1) System.out.println("Attaching hit " + ihit + " to cluster " + cluster.getID()); + break; + } + } + } + if(hit.get_ClusID()==0) { // new cluster found + MUCALCluster cluster = new MUCALCluster(clusters.size()+1); + hit.set_ClusID(cluster.getID()); + cluster.add(hit); + clusters.add(cluster); + if(debugMode>=1) System.out.println("Creating new cluster with ID " + cluster.getID()); + } + } + return clusters; + } + + public void selectClusters(List clusters, ConstantsManager manager, int run) { + + IndexedTable clusterTable = manager.getConstants(run, "/calibration/ft/ftcal/cluster"); + + for(int i=0; i=1) System.out.println("Setting status for cluster " + i + " " + clusters.get(i).toString()); + } + } + + + public void writeBanks(DataEvent event, List hits, List clusters, ConstantsManager manager, int run){ + + IndexedTable energyTable = manager.getConstants(run, "/calibration/ft/ftcal/energycorr"); + + // hits banks + if(!hits.isEmpty()) { + DataBank bankHits = event.createBank("MUCAL::hits", hits.size()); + if(bankHits==null){ + System.out.println("ERROR CREATING BANK : MUCAL::hits"); + return; + } + for(int i = 0; i < hits.size(); i++){ + bankHits.setByte("idx",i,(byte) hits.get(i).get_IDX()); + bankHits.setByte("idy",i,(byte) hits.get(i).get_IDY()); + bankHits.setFloat("x",i,(float) (hits.get(i).get_Dx()/10.0)); + bankHits.setFloat("y",i,(float) (hits.get(i).get_Dy()/10.0)); + bankHits.setFloat("z",i,(float) (hits.get(i).get_Dz()/10.0)); + bankHits.setFloat("energy",i,(float) hits.get(i).get_Edep()); + bankHits.setFloat("time",i,(float) hits.get(i).get_Time()); + bankHits.setShort("hitID",i,(short) hits.get(i).get_DGTZIndex()); + if(!clusters.isEmpty() && clusters.get(hits.get(i).get_ClusID()-1).getStatus()) { + bankHits.setShort("clusterID",i,(short) hits.get(i).get_ClusID()); + } + else { + bankHits.setShort("clusterID",i,(short) 0); + } + } + if(debugMode>=1) bankHits.show(); + event.appendBanks(bankHits); + } + // cluster bank + if(!clusters.isEmpty()) { + List selectedClusters = new ArrayList(); + for(int i =0; i< clusters.size(); i++) { + if(clusters.get(i).getStatus()) selectedClusters.add(clusters.get(i)); + } + if(!selectedClusters.isEmpty()) { + DataBank bankCluster = event.createBank("MUCAL::clusters", selectedClusters.size()); + if(bankCluster==null){ + System.out.println("ERROR CREATING BANK : MUCAL::clusters"); + return; + } + for(int i = 0; i < selectedClusters.size(); i++){ + bankCluster.setShort("id", i,(short) selectedClusters.get(i).getID()); + bankCluster.setShort("size", i,(short) selectedClusters.get(i).getSize()); + bankCluster.setFloat("x",i,(float) (selectedClusters.get(i).getX()/10.0)); + bankCluster.setFloat("y",i, (float) (selectedClusters.get(i).getY()/10.0)); + bankCluster.setFloat("z",i, (float) (selectedClusters.get(i).getZ()/10.0)); + bankCluster.setFloat("widthX",i, (float) (selectedClusters.get(i).getWidthX()/10.0)); + bankCluster.setFloat("widthY",i, (float) (selectedClusters.get(i).getWidthY()/10.0)); + bankCluster.setFloat("radius",i, (float) (selectedClusters.get(i).getRadius()/10.0)); + bankCluster.setFloat("time",i, (float) selectedClusters.get(i).getTime()); + bankCluster.setFloat("energy",i, (float) selectedClusters.get(i).getFullEnergy(energyTable)); + bankCluster.setFloat("recEnergy",i, (float) selectedClusters.get(i).getEnergy()); + bankCluster.setFloat("maxEnergy",i, (float) selectedClusters.get(i).getSeedEnergy()); + } + if(debugMode>=1) bankCluster.show(); + event.appendBanks(bankCluster); + } + } + } + + public List readRawHits(DataEvent event, IndexedTable charge2Energy, IndexedTable timeOffsets, IndexedTable timeWalk, IndexedTable cluster, IndexedTable status) { + // getting raw data bank + if(debugMode>=1) System.out.println("Getting raw hits from MUCAL:adc bank"); + + List hits = new ArrayList<>(); + if(event.hasBank("MUCAL::adc")==true) { + RawDataBank bankDGTZ = new RawDataBank("MUCAL::adc"); + bankDGTZ.read(event); + int nrows = bankDGTZ.rows(); + for(int row = 0; row < nrows; row++){ + int isector = bankDGTZ.getByte("sector",row); + int ilayer = bankDGTZ.getByte("layer",row); + int icomponent = bankDGTZ.getShort("component",row); + 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, MUCALHit.REFCOMPONENT)==0){ + MUCALHit hit = new MUCALHit(bankDGTZ.trueIndex(row),icomponent, adc, time, charge2Energy, timeOffsets, timeWalk, cluster); + hits.add(hit); + } + } + } + return hits; + } + +} diff --git a/reconstruction/mu/src/main/java/org/jlab/rec/muhd/MUHDCluster.java b/reconstruction/mu/src/main/java/org/jlab/rec/muhd/MUHDCluster.java new file mode 100644 index 0000000000..4abc165684 --- /dev/null +++ b/reconstruction/mu/src/main/java/org/jlab/rec/muhd/MUHDCluster.java @@ -0,0 +1,235 @@ +package org.jlab.rec.muhd; + + +import java.util.ArrayList; +import org.jlab.geom.prim.Point3D; + +import org.jlab.io.evio.EvioDataBank; + + + +public class MUHDCluster extends ArrayList { + + /** + * A signal in the hodoscope is a pair of hits above threshold in the two layers. + * The hit pair is define according o the algorithm in HodoSignalFinder + */ + private static final long serialVersionUID = 1L; + + private int _clusID; // signal ID + + + // constructor + public MUHDCluster(int cid) { + this.setID(cid); + } + + public int getID() { + return _clusID; + } + + public void setID(int _clusterID) { + this._clusID = _clusterID; + } + + public int getSize() { + return this.size(); + } + + public double getEnergy() { + // return measured energy + double clusterEnergy = 0; + for(int i=0; i MUHDConstants.CLUSTER_MIN_SIZE && + this.getEnergy()> MUHDConstants.CLUSTER_MIN_ENERGY && + this.getLayerMultiplicity(1)>0 && this.getLayerMultiplicity(2)>0) { + return true; + } + else { + return false; + } + } + + public boolean containsHit(MUHDHit hit) { + boolean addFlag = false; + 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_Dx() - this.get(j).get_Dx()); + double yDiff = Math.abs(hit.get_Dy() - this.get(j).get_Dy()); +// System.out.println("DT: " + tDiff + "(" + MUHDConstantsLoader.time_window +// + ")\t DX: " + xDiff + "(" + MUHDConstantsLoader.hit_distance +// + ")\t DY: " + yDiff + "(" + MUHDConstantsLoader.hit_distance + ")"); + if(tDiff <= MUHDConstants.TIME_WINDOW && + xDiff <= MUHDConstants.HIT_DISTANCE && + yDiff <= MUHDConstants.HIT_DISTANCE) addFlag = true; + } + return addFlag; + } + + public void showCluster() { + System.out.println("\nCluster = " + this._clusID); + System.out.println("Size = " + this.getSize() + + "\tE = " + this.getEnergy() + + "\tTime = " + this.getTime() + + "\tTheta = " + this.getTheta() + + "\tPhi = " + this.getPhi()); + for(int j = 0; j< this.size(); j++) { + System.out.println("hit # " + j + "\t" + this.get(j).get_Layer() + "\t" + this.get(j).get_Sector() + "\t" + this.get(j).get_ID() + "\t" + this.get(j).get_Edep()); + } + } + + + public static ArrayList getSignals(EvioDataBank bank) { + + int[] signalID = bank.getInt("signalID"); + int[] signalSize = bank.getInt("signalSize"); + double[] signalX = bank.getDouble("signalX"); + double[] signalY = bank.getDouble("signalY"); + double[] signalDX = bank.getDouble("signalDX"); + double[] signalDY = bank.getDouble("signalDY"); + double[] signalTime = bank.getDouble("signalTime"); + double[] signalEnergy = bank.getDouble("signalEnergy"); + double[] signalTheta = bank.getDouble("signalTheta"); + double[] signalPhi = bank.getDouble("signalPhi"); + + + int size = signalID.length; + ArrayList signals = new ArrayList(); + +// for(int i = 0; i=0) { + // get hits fron banks + List allHits = reco.initMUHD(event,this.getConstantsManager(), run); + // select good hits and order them by energy + List selectedHits = reco.selectHits(allHits); + // create clusters + List clusters = reco.findClusters(selectedHits); + // write output banks + reco.writeBanks(event, selectedHits, clusters); + } + return true; + } + + public int setRunConditionsParameters(DataEvent event) { + int run = -1; + if(event.hasBank("RUN::config")==false) { + System.out.println("RUN CONDITIONS NOT READ!"); + } + + if(event instanceof EvioDataEvent) { + EvioDataBank bank = (EvioDataBank) event.getBank("RUN::config"); + run = bank.getInt("Run",0); + } + else { + DataBank bank = event.getBank("RUN::config"); + run = bank.getInt("run",0); + } + return run; + } + + + public static void main (String arg[]) { + MUHDEngine cal = new MUHDEngine(); + cal.init(); +// String input = "/Users/devita/data/out_clasdispr.00.e11.000.emn0.75tmn.09.xs65.61nb.dis.1.V5.hipo"; +// HipoDataSource reader = new HipoDataSource(); + String input = "/Users/devita/test4.hipo"; + HipoDataSource reader = new HipoDataSource(); + reader.open(input); + + // initialize histos + H1F h1 = new H1F("Cluster Energy",100, 0.,5.); + h1.setOptStat(Integer.parseInt("1111")); h1.setTitleX("Cluster Energy (GeV)"); + H1F h2 = new H1F("Energy Resolution",100, -1, 1); + h2.setOptStat(Integer.parseInt("1111")); h2.setTitleX("Energy Resolution(GeV)"); + H1F h3 = new H1F("Theta Resolution",100, -2, 2); + h3.setOptStat(Integer.parseInt("1111")); h3.setTitleX("Theta Resolution(deg)"); + H1F h4 = new H1F("Phi Resolution",100, -10, 10); + h4.setOptStat(Integer.parseInt("1111")); h4.setTitleX("Phi Resolution(deg)"); + H1F h5 = new H1F("Time Resolution",100, -10, 10); + h5.setOptStat(Integer.parseInt("1111")); h5.setTitleX("Time Resolution(ns)"); + + while(reader.hasEvent()){ + DataEvent event = (DataEvent) reader.getNextEvent(); + cal.processDataEvent(event); + + DetectorEvent detectorEvent = DetectorData.readDetectorEvent(event); + PhysicsEvent gen = detectorEvent.getGeneratedEvent(); + if(event.hasBank("MUHD::clusters")) { + DataBank bank = event.getBank("MUHD::clusters"); + int nrows = bank.rows(); + for(int i=0; i{ + // class implements Comparable interface to allow for sorting a collection of hits by Edep values + + // constructor + public MUHDHit(int i, int Sector, int Layer, int ID, int ADC, int TDC, IndexedTable charge2Energy, IndexedTable timeOffsets, IndexedTable geometry) { + this._Sector = Sector; + this._Layer = Layer; + this._ID = ID; + this._ADC = ADC; + this._TDC = TDC; + + this.set_Edep(((double) this._ADC)*MUHDConstants.FADC_TO_CHARGE + *charge2Energy.getDoubleValue("mips_energy", Sector, Layer, ID) + /charge2Energy.getDoubleValue("mips_charge", Sector, Layer, ID)); + this.set_Time(((double) this._TDC)/MUHDConstants.TIMECONVFAC + -timeOffsets.getDoubleValue("time_offset", Sector, Layer, ID)); // Time set to gemc value + this.set_Dx(geometry.getDoubleValue("x", Sector, Layer, ID)); + this.set_Dy(geometry.getDoubleValue("y", Sector, Layer, ID)); + this.set_Dz(geometry.getDoubleValue("z", Sector, Layer, ID)); + this.set_DGTZIndex(i); + this.set_ClusterIndex(0); +// System.out.println(this._Dx + " " + this._Dy); + } + + public MUHDHit(int i, int Sector, int Layer, int ID, int ADC, float time, IndexedTable charge2Energy, IndexedTable timeOffsets, IndexedTable geometry) { + this._Sector = Sector; + this._Layer = Layer; + this._ID = ID; + this._ADC = ADC; + + this.set_Edep(((double) this._ADC)*MUHDConstants.FADC_TO_CHARGE + *charge2Energy.getDoubleValue("mips_energy", Sector, Layer, ID) + /charge2Energy.getDoubleValue("mips_charge", Sector, Layer, ID)); + this.set_Time(time-timeOffsets.getDoubleValue("time_offset", Sector, Layer, ID)); // Time set to gemc value + this.set_Dx(geometry.getDoubleValue("x", Sector, Layer, ID)); + this.set_Dy(geometry.getDoubleValue("y", Sector, Layer, ID)); + this.set_Dz(geometry.getDoubleValue("z", Sector, Layer, ID)); + this.set_DGTZIndex(i); + this.set_ClusterIndex(0); +// System.out.println(this._Dx + " " + this._Dy); + } + + + private int _ID; // ID + private int _Sector; // Sector + private int _Layer; // Layer + private int _ADC; // ADC + private int _TDC; // TDC + + private double _Edep; // Reconstructed energy deposited by the hit in the tile + private double _Time; // Reconstructed time, for now it is the gemc time + private double _Dx; // X position (corresponds to tile centroid) + private double _Dy; // Y position (corresponds to tile centroid) + private double _Dz; // Z position (corresponds to tile centroid) + private int _DGTZIndex; // Pointer to cluster + private int _ClusterIndex; // Pointer to cluster + + + + public int get_ID() { + return _ID; + } + + + + public void set_ID(int ID) { + this._ID = ID; + } + + + + public int get_Sector() { + return _Sector; + } + + + + public void set_Sector(int Sector) { + this._Sector = Sector; + } + + + + public int get_Layer() { + return _Layer; + } + + + + public void set_Layer(int Layer) { + this._Layer = Layer; + } + + + + public int get_ADC() { + return _ADC; + } + + + + public void set_ADC(int ADC) { + this._ADC = ADC; + } + + + + public int get_TDC() { + return _TDC; + } + + + + public void set_TDC(int TDC) { + this._TDC = TDC; + } + + + public double get_Edep() { + return _Edep; + } + + + public void set_Edep(double edep) { + this._Edep = edep; + } + + + + public double get_Time() { + return _Time; + } + + + public void set_Time(double Time) { + this._Time = Time; + } + + + public double get_Dx() { + return _Dx; + } + + + + public void set_Dx(double Dx) { + this._Dx = Dx; + } + + + + public double get_Dy() { + return _Dy; + } + + + + public void set_Dy(double Dy) { + this._Dy = Dy; + } + + + public double get_Dz() { + return _Dz; + } + + + + public void set_Dz(double Dz) { + this._Dz = Dz; + } + + + public int get_DGTZIndex() { + return _DGTZIndex; + } + + + + public void set_DGTZIndex(int _DGTZIndex) { + this._DGTZIndex = _DGTZIndex; + } + + + + public int get_ClusterIndex() { + return _ClusterIndex; + } + + + + public void set_ClusterIndex(int _SignalIndex) { + this._ClusterIndex = _SignalIndex; + } + + + +// public static ArrayList getRawHits(EvioDataBank bank) { +// +// int[] id = bank.getInt("component"); +// int[] sector = bank.getInt("sector"); +// int[] layer = bank.getInt("layer"); +// int[] adc = bank.getInt("ADC"); +// int[] tdc = bank.getInt("TDC"); +// +// int size = id.length; +// +// ArrayList hits = new ArrayList(); +// +// for(int i = 0; i MUHDConstants.EN_THRES) { + return true; + } else { + return false; + } + } + + + + @Override + public int compareTo(MUHDHit arg0) { + if(this.get_Edep() initMUHD(DataEvent event, ConstantsManager manager, int run) { + + IndexedTable charge2Energy = manager.getConstants(run, "/calibration/ft/fthodo/charge_to_energy"); + IndexedTable timeOffsets = manager.getConstants(run, "/calibration/ft/fthodo/time_offsets"); + IndexedTable status = manager.getConstants(run, "/calibration/ft/fthodo/status"); + IndexedTable geometry = manager.getConstants(run, "/geometry/ft/fthodo"); + + if(debugMode>=1) System.out.println("\nAnalyzing new event"); + List allhits = null; + + if(event instanceof HipoDataEvent) { + allhits = this.readRawHits(event,charge2Energy,timeOffsets,status,geometry); + } + if(debugMode>=1) { + System.out.println("Found " + allhits.size() + " hits"); + for(int i = 0; i < allhits.size(); i++) { + System.out.print(i + "\t"); + allhits.get(i).showHit(); + } + } + return allhits; + } + + public List selectHits(List allhits) { + + if(debugMode>=1) System.out.println("\nSelecting hits"); + ArrayList hits = new ArrayList<>(); + + for(int i = 0; i < allhits.size(); i++) + { + if(MUHDHit.passHitSelection(allhits.get(i))) { + hits.add(allhits.get(i)); + } + } + Collections.sort(hits); + if(debugMode>=1) { + System.out.println("List of selected hits"); + for(int i = 0; i < hits.size(); i++) + { + System.out.print(i + "\t"); + hits.get(i).showHit(); + } + } + return hits; + } + + public List findClusters(List hits) { + + List clusters = new ArrayList(); + + if(debugMode>=1) System.out.println("\nBuilding clusters"); + for(int ihit=0; ihit=1) System.out.println("Attaching hit " + ihit + " to cluster " + cluster.getID()); + } + } + } + if(hit.get_ClusterIndex()==0) { // new cluster found + MUHDCluster cluster = new MUHDCluster(clusters.size()+1); + hit.set_ClusterIndex(cluster.getID()); + cluster.add(hit); + clusters.add(cluster); + if(debugMode>=1) System.out.println("Creating new cluster with ID " + cluster.getID()); + } + } + return clusters; + } + + public void writeBanks(DataEvent event, List hits, List clusters){ + // hits banks + if(!hits.isEmpty()) { + DataBank bankHits = event.createBank("MUHD::hits", hits.size()); + if(bankHits==null){ + System.out.println("ERROR CREATING BANK : MUHD::hits"); + return; + } + for(int i = 0; i < hits.size(); i++){ + bankHits.setByte("sector",i,(byte) hits.get(i).get_Sector()); + bankHits.setByte("layer",i,(byte) hits.get(i).get_Layer()); + bankHits.setShort("component",i,(short) hits.get(i).get_ID()); + bankHits.setFloat("x",i,(float) (hits.get(i).get_Dx()/10.0)); + bankHits.setFloat("y",i,(float) (hits.get(i).get_Dy()/10.0)); + bankHits.setFloat("z",i,(float) (hits.get(i).get_Dz()/10.0)); + bankHits.setFloat("energy",i,(float) hits.get(i).get_Edep()); + bankHits.setFloat("time",i,(float) hits.get(i).get_Time()); + bankHits.setShort("hitID",i,(short) hits.get(i).get_DGTZIndex()); + bankHits.setShort("clusterID",i,(short) hits.get(i).get_ClusterIndex()); + } + event.appendBanks(bankHits); + } + // cluster bank + if(!clusters.isEmpty()){ + DataBank bankCluster = event.createBank("MUHD::clusters", clusters.size()); + if(bankCluster==null){ + System.out.println("ERROR CREATING BANK : MUHD::clusters"); + return; + } + for(int i = 0; i < clusters.size(); i++){ + bankCluster.setShort("id", i,(short) clusters.get(i).getID()); + bankCluster.setShort("size", i,(short) clusters.get(i).getSize()); + bankCluster.setFloat("x",i,(float) (clusters.get(i).getX()/10.0)); + bankCluster.setFloat("y",i,(float) (clusters.get(i).getY()/10.0)); + bankCluster.setFloat("z",i,(float) (clusters.get(i).getZ()/10.0)); + bankCluster.setFloat("widthX",i,(float) (clusters.get(i).getWidthX()/10.0)); + bankCluster.setFloat("widthY",i,(float) (clusters.get(i).getWidthY()/10.0)); + bankCluster.setFloat("radius",i,(float) (clusters.get(i).getRadius()/10.0)); + bankCluster.setFloat("time",i,(float) clusters.get(i).getTime()); + bankCluster.setFloat("energy",i,(float) clusters.get(i).getEnergy()); + } + event.appendBanks(bankCluster); + } + } + + public List readRawHits(DataEvent event, IndexedTable charge2Energy, IndexedTable timeOffsets, IndexedTable status, IndexedTable geometry) { + // getting raw data bank + if(debugMode>=1) System.out.println("Getting raw hits from MUHD:adc bank"); + + List hits = new ArrayList<>(); + if(event.hasBank("MUHD::adc")==true) { + RawDataBank bankDGTZ = new RawDataBank("MUHD::adc"); + bankDGTZ.read(event); + int nrows = bankDGTZ.rows(); + for(int row = 0; row < nrows; row++){ + int isector = bankDGTZ.getByte("sector",row); + int ilayer = bankDGTZ.getByte("layer",row); + int icomponent = bankDGTZ.getShort("component",row); + int adc = bankDGTZ.getInt("ADC",row); + float time = bankDGTZ.getFloat("time",row); + if(adc!=-1 && time!=-1 && status.getIntValue("status", isector, ilayer, icomponent)==0){ + MUHDHit hit = new MUHDHit(bankDGTZ.trueIndex(row),isector,ilayer,icomponent, adc, time, charge2Energy,timeOffsets,geometry); + hits.add(hit); + } + } + } + return hits; + } +} diff --git a/reconstruction/mu/src/main/java/org/jlab/rec/muvt/MUVTCluster.java b/reconstruction/mu/src/main/java/org/jlab/rec/muvt/MUVTCluster.java new file mode 100644 index 0000000000..36df92ca51 --- /dev/null +++ b/reconstruction/mu/src/main/java/org/jlab/rec/muvt/MUVTCluster.java @@ -0,0 +1,221 @@ +package org.jlab.rec.muvt; + +import java.util.ArrayList; +import java.util.List; +import org.jlab.detector.base.DetectorDescriptor; +import org.jlab.detector.base.DetectorType; +import org.jlab.geom.prim.Line3D; +import org.jlab.geom.prim.Point3D; + +/** + * MUVT in-layer cluster + * + * @author bondi, devita, caot + */ +public class MUVTCluster extends ArrayList { + + + private DetectorDescriptor desc = new DetectorDescriptor(DetectorType.MUVT); + private int id; + private Line3D clusterLine = null; + private Line3D clusterLineTS = null; + public int indexMaxStrip = -1; + private byte clusterStatus = 1; + + public MUVTCluster(MUVTStrip strip){ + this.desc.setSectorLayerComponent(strip.getDescriptor().getSector(), + strip.getDescriptor().getLayer(), 0); + this.add(strip); + this.clusterLine = new Line3D(strip.getLine()); + this.clusterLineTS = MUVTConstants.toTiltedSectorFrame(this.desc.getSector(), clusterLine); + this.indexMaxStrip = 0; + } + + public int getId() { + return id; + } + + public DetectorDescriptor getDescriptor(){ + return this.desc; + } + + public int getSector() { + return this.desc.getSector(); + } + + public int getLayer() { + return this.desc.getLayer(); + } + + public int getChamber() { + return this.get(0).getChamber(); + } + + public Line3D getLine() { + return this.clusterLine; + } + + public Line3D getLineTS() { + return this.clusterLineTS; + } + + public double getEnergy(){ + double energy = 0.0; + for(MUVTStrip strip : this){ + energy += strip.getEnergy(); + } + return energy; + } + + public double getTime(){ + double time = 0.0; + for(MUVTStrip strip : this){ + time += strip.getTime()*strip.getEnergy(); + } + time /= this.getEnergy(); + return time; + } + + public double getSeedTime(){ + if(this.indexMaxStrip >= 0 && this.indexMaxStrip < this.size()){ + return this.get(indexMaxStrip).getTime(); + } + return 0.0; + } + + public MUVTStrip getSeedStrip() { + return this.get(this.indexMaxStrip); + } + + public int getMaxStrip(){ + return this.get(this.indexMaxStrip).getDescriptor().getComponent(); + } + + public boolean addStrip(MUVTStrip strip){ + for(MUVTStrip s : this){ + if(s.isNeighbour(strip)){ + this.add(strip); + if(strip.getEnergy()>this.get(indexMaxStrip).getEnergy()){ + this.indexMaxStrip = this.size()-1; + this.clusterLine.copy(strip.getLine()); + } + return true; + } + } + return false; + } + + public int getADC(){ + int adc = 0; + for(MUVTStrip s : this){ + adc+= s.getADC(); + } + return adc; + } + + public void setStatus(int val) {this.clusterStatus+=val;} + + public byte getStatus() {return clusterStatus;} + + public void setClusterId(int id){ + this.id = id; + for(MUVTStrip strip : this){ + strip.setClusterId(id); + } + } + + public void redoClusterLine(){ + + Point3D pointOrigin = new Point3D(0.0,0.0,0.0); + Point3D pointEnd = new Point3D(0.0,0.0,0.0); + + double logSumm = 0.0; + double summE = 0.0; + + for(int i = 0; i < this.size(); i++){ + Line3D line = this.get(i).getLine(); + + double energy = this.get(i).getEnergy(); + double energymev = energy*1000.0; + double le = Math.log(energymev); + + pointOrigin.setX(pointOrigin.x()+line.origin().x()*le); + pointOrigin.setY(pointOrigin.y()+line.origin().y()*le); + pointOrigin.setZ(pointOrigin.z()+line.origin().z()*le); + + pointEnd.setX(pointEnd.x()+line.end().x()*le); + pointEnd.setY(pointEnd.y()+line.end().y()*le); + pointEnd.setZ(pointEnd.z()+line.end().z()*le); + + logSumm += le; + + summE += energy; + } + + this.clusterLine.set( + pointOrigin.x()/logSumm, + pointOrigin.y()/logSumm, + pointOrigin.z()/logSumm, + pointEnd.x()/logSumm, + pointEnd.y()/logSumm, + pointEnd.z()/logSumm + ); + this.clusterLineTS = MUVTConstants.toTiltedSectorFrame(this.desc.getSector(), clusterLine); + } + + + public static List createClusters(List stripList){ + + List clusterList = new ArrayList<>(); + + if(!stripList.isEmpty()){ + for(int loop = 0; loop < stripList.size(); loop++){ //Loop over all strips + boolean stripAdded = false; + for(MUVTCluster cluster : clusterList) { + if(cluster.addStrip(stripList.get(loop))){ //Add adjacent strip to newly seeded peak + stripAdded = true; + } + } + if(!stripAdded){ + MUVTCluster newPeak = new MUVTCluster(stripList.get(loop)); //Non-adjacent strip seeds new peak + clusterList.add(newPeak); + } + } + } + for(int loop = 0; loop < clusterList.size(); loop++){ + clusterList.get(loop).setClusterId(loop+1); + clusterList.get(loop).redoClusterLine(); + } + return clusterList; + } + + public static List getClusters(List clusters, int sector, int layer) { + List selectedClusters = new ArrayList<>(); + for(MUVTCluster cluster : clusters) { + if(cluster.getSector()==sector && cluster.getLayer()==layer) + selectedClusters.add(cluster); + } + return selectedClusters; + } + + @Override + public String toString(){ + StringBuilder str = new StringBuilder(); + str.append(String.format("----> cluster ( %3d %3d ) ENERGY = %12.5f\n", + this.desc.getSector(),this.desc.getLayer(), this.getEnergy())); + str.append(this.clusterLine.toString()); + str.append("\n"); + for(MUVTStrip strip : this){ + str.append("\t\t"); + str.append(strip.toString()); + str.append("\n"); + } + + return str.toString(); + } + + + + + +} diff --git a/reconstruction/mu/src/main/java/org/jlab/rec/muvt/MUVTConstants.java b/reconstruction/mu/src/main/java/org/jlab/rec/muvt/MUVTConstants.java new file mode 100644 index 0000000000..863236df8c --- /dev/null +++ b/reconstruction/mu/src/main/java/org/jlab/rec/muvt/MUVTConstants.java @@ -0,0 +1,118 @@ +package org.jlab.rec.muvt; + +import org.jlab.geom.base.Detector; +import org.jlab.geom.detector.fmt.FMTLayer; +import org.jlab.geom.prim.Line3D; +import org.jlab.geom.prim.Point3D; +import org.jlab.geom.prim.Transformation3D; +import org.jlab.geom.prim.Vector3D; + +/** + * + * @author devita + */ +public class MUVTConstants { + + public final static int NSECTOR = 6; + public final static int NLAYER = 12; + public final static int NREGION = 6; + public final static int NCHAMBER = 1; + public final static double PITCH = 0.5; // mm + public final static double TILT = 25; // deg + public final static double[] STEREO = { 10.0, 10.0 }; + + // strips + public final static double THRESHOLD = 0; + public final static double ADCTOENERGY = 25/1E4; // in eV, values from gemc ADC = (uRwellC.gain*1e6*tInfos.eTot/uRwellC.w_i); with gain = 10^4 and w_i = 25 eV + public final static double TDCTOTIME = 1; + + // cluster + public final static double COINCTIME = 100; + + // rMax difference between cross's clusters energy + public static double CROSSDELTAE = 200; + + // DC-tracks to FMT-clusters matching parameter + 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; + + // small distance (cm) for derivatives calculations + public static double EPSILON = 1e-4; + + public static int MAX_NB_CROSSES = 30; + +// private static Detector fmtDetector = null; + + + +// public static void setDetector(Detector detector) { +// MUVTConstants.fmtDetector = detector; +// } + +// public static Vector3D getDerivatives(int layer, double x, double y, double z) { +// Vector3D p0 = new Vector3D(x,y,z); +// Vector3D p1 = new Vector3D(x,y+MUVTConstants.EPSILON,z); +// MUVTConstants.getInverseTransform(layer).apply(p0); +// MUVTConstants.getInverseTransform(layer).apply(p1); +// return p1.sub(p0).divide(MUVTConstants.EPSILON); +// } + +// public static FMTLayer getLayer(int layer) { +// if(layer<1 || layer>NLAYER) +// throw new IllegalArgumentException("Error: invalid layer="+layer); +// return (FMTLayer) MUVTConstants.fmtDetector.getSector(0).getSuperlayer(0).getLayer(layer-1); +// } +// +// public static double getPitch() { +// return MUVTConstants.getLayer(1).getComponent(0).getWidth(); +// } +// +// public static Line3D getStrip(int layer, int strip) { +// if(strip<1 || strip>MUVTConstants.getLayer(layer).getNumComponents()) +// throw new IllegalArgumentException("Error: invalid strip="+strip); +// return MUVTConstants.getLayer(layer).getComponent(strip-1).getLine(); +// } + +// public static Line3D getLocalStrip(int layer, int strip) { +// Line3D local = new Line3D(MUVTConstants.getStrip(layer, strip)); +// MUVTConstants.getInverseTransform(layer).apply(local); +// return local; +// } +// +// public static double getThickness() { +// return MUVTConstants.getLayer(1).getComponent(0).getThickness(); +// } + + public static Transformation3D toTiltedSectorFrame(int sector) { + Transformation3D t = new Transformation3D(); + t.rotateZ(-60 * (sector-1)); + t.rotateZ(-TILT); + return t; + } + + public static Transformation3D toLab(int sector) { + Transformation3D inverse = new Transformation3D(MUVTConstants.toTiltedSectorFrame(sector)); + return inverse.inverse(); + } + + public static Point3D toTiltedSectorFrame(int sector, Point3D p) { + Point3D local = new Point3D(p); + MUVTConstants.toTiltedSectorFrame(sector).apply(local); + return local; + } + + public static Point3D toTiltedSectorFrame(int sector, double x, double y, double z) { + Point3D local = new Point3D(x,y,z); + MUVTConstants.toTiltedSectorFrame(sector).apply(local); + return local; + } + + public static Line3D toTiltedSectorFrame(int sector, Line3D l) { + Line3D local = new Line3D(l); + MUVTConstants.toTiltedSectorFrame(sector).apply(local); + return local; + } + +} diff --git a/reconstruction/mu/src/main/java/org/jlab/rec/muvt/MUVTCross.java b/reconstruction/mu/src/main/java/org/jlab/rec/muvt/MUVTCross.java new file mode 100644 index 0000000000..4fde7af6b2 --- /dev/null +++ b/reconstruction/mu/src/main/java/org/jlab/rec/muvt/MUVTCross.java @@ -0,0 +1,129 @@ +package org.jlab.rec.muvt; + +import java.util.ArrayList; +import java.util.List; +import org.jlab.geom.prim.Plane3D; +import org.jlab.geom.prim.Point3D; +import org.jlab.geom.prim.Vector3D; + +/** + * MUVT V-W clusters + * @author devita + */ +public class MUVTCross { + + private int id; + + private int sector; + private int region; + private int chamber; + + private int cluster1; + private int cluster2; + + private Point3D cross; + private double energy; + private double deltaE; + private double time; + private int status; + + + + public MUVTCross(MUVTCluster c1, MUVTCluster c2) { + + Vector3D dir = c1.getLine().direction().cross(c2.getLine().direction()); + Plane3D plane = new Plane3D(c1.getLine().origin(), c1.getLine().direction().cross(dir)); + Point3D point = new Point3D(); + int nint = plane.intersectionSegment(c2.getLine(), point); + if(nint==1) { + this.sector = c1.getSector(); + this.region = (c1.getLayer()-1)/(MUVTConstants.NLAYER/MUVTConstants.NREGION)+1; + this.cross = point; + this.energy = c1.getEnergy() + c2.getEnergy(); + this.deltaE = c1.getEnergy() - c2.getEnergy(); + this.time = (c1.getTime() + c2.getTime())/2; + this.cluster1 = c1.getId(); + this.cluster2 = c2.getId(); + } + } + + public void setId(int id) { + this.id = id; + } + + public int getId() { + return id; + } + + public int getSector() { + return this.sector; + } + + public int getRegion() { + return this.region; + } + + public int getChamber() { + return this.chamber; + } + + public int getCluster1() { + return cluster1; + } + + public int getCluster2() { + return cluster2; + } + + public Point3D point() { + return cross; + } + + public double getEnergy() { + return energy; + } + + public double getDeltaEnergy() { + return deltaE; + } + + public double getTime() { + return time; + } + + public int getStatus() { + return status; + } + + public static List createCrosses(List clusters) { + + List crosses = new ArrayList<>(); + + for(int is=0; is clustersV = MUVTCluster.getClusters(clusters, is+1, (MUVTConstants.NLAYER/MUVTConstants.NREGION)*ir+1); + List clustersW = MUVTCluster.getClusters(clusters, is+1, (MUVTConstants.NLAYER/MUVTConstants.NREGION)*ir+2); + + for(MUVTCluster v : clustersV) { + for(MUVTCluster w : clustersW) { + + if(v.getChamber()==w.getChamber()) { + MUVTCross cross = new MUVTCross(v, w); + if(cross.point()!=null) crosses.add(cross); + cross.setId(crosses.size()); + } + } + } + } + } + return crosses; + } + + @Override + public String toString(){ + StringBuilder str = new StringBuilder(); + str.append(String.format("----> cross ( %3d %3d )\n", this.getSector(),this)); + str.append(this.point().toString()); + return str.toString(); + } +} diff --git a/reconstruction/mu/src/main/java/org/jlab/rec/muvt/MUVTEngine.java b/reconstruction/mu/src/main/java/org/jlab/rec/muvt/MUVTEngine.java new file mode 100644 index 0000000000..4823e2c4b7 --- /dev/null +++ b/reconstruction/mu/src/main/java/org/jlab/rec/muvt/MUVTEngine.java @@ -0,0 +1,288 @@ +package org.jlab.rec.muvt; + +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.geant4.v2.URWELL.URWellStripFactory; +import org.jlab.geom.prim.Line3D; +import org.jlab.geom.prim.Point3D; +import org.jlab.utils.groups.IndexedTable; +import org.jlab.io.base.*; + +import org.jlab.rec.muvt.track.fit.KFitter; +import org.jlab.rec.muvt.track.fit.StateVecs.StateVec; + +/** + * Service to return reconstructed track candidates - the output is in hipo format + * + * @author ziegler, benkel, devita + */ +public class MUVTEngine extends ReconstructionEngine { + + boolean debug = false; + + public static URWellStripFactory factory = new URWellStripFactory(); + + public MUVTEngine() { + super("MUVT", "ziegler", "5.0"); + } + + @Override + public boolean init() { + + // Get the constants for the correct variation + String variation = this.getEngineConfigString("variation"); + if (variation!=null) { + System.out.println("["+this.getName()+"] " + + "run with MUVT geometry variation based on yaml = " + variation); + } + else { + variation = "default"; + System.out.println("["+this.getName()+"] run with MUVT default geometry"); + } + + String[] tables = new String[]{ + "/geometry/beam/position" + }; + requireConstants(Arrays.asList(tables)); + this.getConstantsManager().setVariation(variation); + // Register output banks + super.registerOutputBank("MUVT::hits"); + super.registerOutputBank("MUVT::clusters"); + super.registerOutputBank("MUVT::crosses"); + super.registerOutputBank("MUVT::tracks"); + super.registerOutputBank("MUVT::trajectory"); + + return true; + } + + @Override + public boolean processDataEvent(DataEvent event) { + // Initial setup. + if(debug) System.out.println("\nNew event"); + + // Set run number. + 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(); + + // Set beam shift. NOTE: Set to zero for the time being, when beam alignment is done + // 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); + + + // === from STRIPS to CROSSES ============================================================= + List strips = MUVTStrip.getStrips(event, factory, this.getConstantsManager()); + List clusters = MUVTCluster.createClusters(strips); + List crosses = MUVTCross.createCrosses(clusters); + + //System.out.println(clusters.size()); + //for (int i = 0; i < clusters.size(); i++) System.out.println(clusters.get(i).toString()); + + // === DC TRACKS =========================================================================== + MUVTTrack trk = new MUVTTrack(); + 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[MUVTConstants.NLAYER/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(MUVTCross ce : end) { + Line3D ve = new Line3D(ce.point(),target); + for(int ro=1; ro segment = new ArrayList(); + segment.add(ce); + for(int r=co.getRegion()+1; r=4) + segments.add(segment); + } + } + } + } + System.out.print(segments.size()); + if(!segments.isEmpty()) { + segments.sort(Comparator.comparingInt(List::size).reversed()); + for(MUVTCross cross : segments.get(0)) { + track.addCluster(clusters.get(cross.getCluster1()-1)); + track.addCluster(clusters.get(cross.getCluster2()-1)); + } + } + } + + List filtedTracks = new ArrayList(); + for(MUVTTrack 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 trackClusters = track.getClusters(); + + 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; + kf.runFitter(track.getSector()); + + if (kf.finalStateVec != null) { + StateVec sv = kf.finalStateVec; + + // swim to beamline to get vertex parameters + 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.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); + + // if successful, save track parameters + if(Vt == null) continue; + track.setStatus(0); + track.setNDF(trackClusters.size()); + track.setQ(charge); + track.setChi2(kf.chi2); + track.setX(Vt[0]); + track.setY(Vt[1]); + track.setZ(Vt[2]); + track.setPx(-Vt[3]); + track.setPy(-Vt[4]); + track.setPz(-Vt[5]); + } + } + + this.writeHipoBanks(event, strips, clusters, crosses, tracks); + return true; + } + + + + private void writeHipoBanks(DataEvent de, + List strips, + List clusters, + List crosses, + List tracks){ + + DataBank bankS = de.createBank("MUVT::hits", strips.size()); + for(int h = 0; h < strips.size(); h++){ + bankS.setShort("id", h, (short) strips.get(h).getId()); + bankS.setByte("sector", h, (byte) strips.get(h).getDescriptor().getSector()); + bankS.setByte("layer", h, (byte) strips.get(h).getDescriptor().getLayer()); + bankS.setShort("strip", h, (short) strips.get(h).getDescriptor().getComponent()); + bankS.setFloat("energy", h, (float) strips.get(h).getEnergy()); + bankS.setFloat("time", h, (float) strips.get(h).getTime()); + bankS.setShort("status", h, (short) strips.get(h).getStatus()); + bankS.setShort("clusterId", h, (short) strips.get(h).getClusterId()); + } + + DataBank bankC = de.createBank("MUVT::clusters", clusters.size()); + for(int c = 0; c < clusters.size(); c++){ + bankC.setShort("id", c, (short) clusters.get(c).getId()); + bankC.setByte("sector", c, (byte) clusters.get(c).get(0).getDescriptor().getSector()); + bankC.setByte("layer", c, (byte) clusters.get(c).get(0).getDescriptor().getLayer()); + bankC.setShort("strip", c, (short) clusters.get(c).getMaxStrip()); + bankC.setFloat("energy", c, (float) clusters.get(c).getEnergy()); + bankC.setFloat("time", c, (float) clusters.get(c).getTime()); + bankC.setFloat("xo", c, (float) clusters.get(c).getLine().origin().x()); + bankC.setFloat("yo", c, (float) clusters.get(c).getLine().origin().y()); + bankC.setFloat("zo", c, (float) clusters.get(c).getLine().origin().z()); + bankC.setFloat("xe", c, (float) clusters.get(c).getLine().end().x()); + bankC.setFloat("ye", c, (float) clusters.get(c).getLine().end().y()); + bankC.setFloat("ze", c, (float) clusters.get(c).getLine().end().z()); + bankC.setShort("size", c, (short) clusters.get(c).size()); + bankC.setShort("status", c, (short) clusters.get(c).getStatus()); + } + + DataBank bankX = de.createBank("MUVT::crosses", crosses.size()); + for(int c = 0; c < crosses.size(); c++){ + bankX.setShort("id", c, (short) crosses.get(c).getId()); + bankX.setByte("sector", c, (byte) crosses.get(c).getSector()); + bankX.setByte("region", c, (byte) crosses.get(c).getRegion()); + bankX.setFloat("energy", c, (float) crosses.get(c).getEnergy()); + bankX.setFloat("time", c, (float) crosses.get(c).getTime()); + bankX.setFloat("x", c, (float) crosses.get(c).point().x()); + bankX.setFloat("y", c, (float) crosses.get(c).point().y()); + bankX.setFloat("z", c, (float) crosses.get(c).point().z()); + bankX.setShort("cluster1", c, (short) crosses.get(c).getCluster1()); + bankX.setShort("cluster2", c, (short) crosses.get(c).getCluster2()); + bankX.setShort("status", c, (short) crosses.get(c).getStatus()); + } + + + DataBank bankT = de.createBank("MUVT::tracks", tracks.size()); + for (int i=0; i this.desc.getSector()) return -1; + if(ob.getDescriptor().getLayer() < this.desc.getLayer()) return 1; + if(ob.getDescriptor().getLayer() > this.desc.getLayer()) return -1; + if(ob.getDescriptor().getComponent() < this.desc.getComponent()) return 1; + if(ob.getDescriptor().getComponent() == this.desc.getComponent()) return 0; + return -1; + } + + public static List getStrips(DataEvent event, URWellStripFactory factory, ConstantsManager ccdb) { + + List strips = new ArrayList<>(); + + if(event.hasBank("URWELL::adc")){ + RawDataBank bank = new RawDataBank("URWELL::adc"); + bank.read(event); + //DataBank bank = event.getBank("URWELL::adc"); + for(int i = 0; i < bank.rows(); i++){ + int sector = bank.getByte("sector", i); + int layer = bank.getByte("layer", i); + int comp = bank.getShort("component", i); + int adc = bank.getInt("ADC", i); + double time = bank.getFloat("time", i); + + MUVTStrip strip = new MUVTStrip(sector, layer, comp); + +// strip.setTriggerPhase(triggerPhase); + strip.setId(bank.trueIndex(i)+1); + strip.setADC(adc); + strip.setTDC((int) time); + strip.setEnergy(strip.ADC*MUVTConstants.ADCTOENERGY); + strip.setTime(strip.TDC*MUVTConstants.TDCTOTIME); + strip.setLine(factory.getStrip(sector, layer, comp)); + strip.setChamber(factory.getChamberIndex(comp)+1); + strip.setStatus(0); + + if(strip.getEnergy()>MUVTConstants.THRESHOLD) strips.add(strip); + + } + } + return strips; + } + + @Override + public String toString(){ + StringBuilder str = new StringBuilder(); + str.append(String.format("----> strip (%3d %3d %3d) ADC/TDC %5d %5d ENERGY = %8.5f TIME = %8.5f ", + this.desc.getSector(),this.desc.getLayer(),this.desc.getComponent(), + this.ADC,this.TDC,this.getEnergy(),this.getTime())); + return str.toString(); + } +} diff --git a/reconstruction/mu/src/main/java/org/jlab/rec/muvt/MUVTTrack.java b/reconstruction/mu/src/main/java/org/jlab/rec/muvt/MUVTTrack.java new file mode 100644 index 0000000000..b07edb79ec --- /dev/null +++ b/reconstruction/mu/src/main/java/org/jlab/rec/muvt/MUVTTrack.java @@ -0,0 +1,464 @@ +package org.jlab.rec.muvt; + +import java.util.ArrayList; +import java.util.LinkedHashMap; +import java.util.List; +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.Point3D; +import org.jlab.io.base.DataBank; +import org.jlab.io.base.DataEvent; + + +/** + * + * @author ziegler + * @author devita + */ +public class MUVTTrack { + + /** + * The status variable explains the number of tracks and the quality of the reconstruction. + * + * Its last digit is the number of FMT layers used in FVT tracking, so it can be any number + * from 0 to 3. If it's 0, it means that no FMT layers were used and the FVT track should be + * the same as the DC track. + * + * If there was an error in swimming due to an odd track shape or anything, a 100 is added to + * the variable to denote that. + */ + private int status; + private int _id; + private int _index; + private int _sector; + private int _q; + private double _chi2; + private double _x; + private double _y; + private double _z; + private double _px; + private double _py; + private double _pz; + private int _NDF; + + private final MUVTTrajectory[] _DCtrajs = new MUVTTrajectory[MUVTConstants.NLAYER]; + private final List[] _clusters = new ArrayList[MUVTConstants.NLAYER]; + private final MUVTTrajectory[] _FMTtrajs = new MUVTTrajectory[MUVTConstants.NLAYER]; + + public MUVTTrack() { + } + + + public MUVTTrack(int _index, int _sector, int _q, double _x, double _y, double _z, + double _px, double _py, double _pz, List clusters) { + this._index = _index; + this._sector = _sector; + this._q = _q; + this._x = _x; + this._y = _y; + this._z = _z; + this._px = _px; + this._py = _py; + this._pz = _pz; + for(MUVTCluster cluster : clusters) this.addCluster(cluster); + } + + + + /** + * @param layer + * @return the _traj + */ + public MUVTTrajectory getDCTraj(int layer) { + if(layer<=0 || layer>MUVTConstants.NLAYER) return null; + else return _DCtrajs[layer-1]; + } + + /** + * @param trj + */ + public void setDCTraj(MUVTTrajectory trj) { + this._DCtrajs[trj.getLayer()-1] = trj; + } + + public List getClusters() { + List clusters = new ArrayList<>(); + for(int i=0; i getClusters(int layer) { + if(layer<=0 || layer>MUVTConstants.NLAYER) return null; + else return _clusters[layer-1]; + } + + public MUVTCluster getCluster(int layer) { + if(layer<=0 || layer>MUVTConstants.NLAYER) return null; + else if(_clusters[layer-1]== null || _clusters[layer-1].size()==0) return null; + else return _clusters[layer-1].get(0); + } + + public int getClusterLayers() { + int n = 0; + for(int i=0; i(); + this._clusters[cluster.getLayer()-1].add(cluster); + } + + public void clearClusters(int layer) { + this._clusters[layer-1].clear(); + } + + public MUVTTrajectory getFMTTraj(int layer) { + if(layer<=0 || layer>MUVTConstants.NLAYER) return null; + return _FMTtrajs[layer-1]; + } + + public void setFMTtraj(MUVTTrajectory trj) { + this._FMTtrajs[trj.getLayer()-1] = trj; + } + + public int getId() { + return _id; + } + + public void setId(int _id) { + this._id = _id; + } + + /** + * @return the _id + */ + public int getIndex() { + return _index; + } + + /** + * @param _id the _id to set + */ + public void setIndex(int _id) { + this._index = _id; + } + + /** + * @return the sector + */ + public int getSector() { + return _sector; + } + + /** + * @param _sector the sector to set + */ + public void setSector(int _sector) { + this._sector = _sector; + } + + /** + * @return the _q + */ + public int getQ() { + return _q; + } + + /** + * @param _q the _q to set + */ + public void setQ(int _q) { + this._q = _q; + } + + /** + * @return the _chi^2. + */ + public double getChi2() { + return _chi2; + } + + /** + * @param _chi2 the _chi2 to set + */ + public void setChi2(double _chi2) { + this._chi2 = _chi2; + } + + public int getStatus() { + return status; + } + + public void setStatus(int status) { + this.status = status; + } + + public int getNDF() { + return _NDF; + } + + public void setNDF(int _NDF) { + this._NDF = _NDF; + } + + /** + * @return the _x + */ + public double getX() { + return _x; + } + + /** + * @param _x the _x to set + */ + public void setX(double _x) { + this._x = _x; + } + + /** + * @return the _y + */ + public double getY() { + return _y; + } + + /** + * @param _y the _y to set + */ + public void setY(double _y) { + this._y = _y; + } + + /** + * @return the _z + */ + public double getZ() { + return _z; + } + + /** + * @param _z the _z to set + */ + public void setZ(double _z) { + this._z = _z; + } + + /** + * @return the the tracke momentum + */ + public double getP() { + return Math.sqrt(_px*_px+_py*_py+_pz*_pz); + } + + /** + * @return the _px + */ + public double getPx() { + return _px; + } + + /** + * @param _px the _px to set + */ + public void setPx(double _px) { + this._px = _px; + } + + /** + * @return the _py + */ + public double getPy() { + return _py; + } + + /** + * @param _py the _py to set + */ + public void setPy(double _py) { + this._py = _py; + } + + /** + * @return the _pz + */ + public double getPz() { + return _pz; + } + + /** + * @param _pz the _pz to set + */ + 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 List getDCTracks(DataEvent event, Swim swimmer) { + + Map trackmap = new LinkedHashMap<>(); + + DataBank trackBank = null; + DataBank trajBank = null; + if(event.hasBank("TimeBasedTrkg::TBTracks")) trackBank = event.getBank("TimeBasedTrkg::TBTracks"); + if(event.hasBank("TimeBasedTrkg::Trajectory")) trajBank = event.getBank("TimeBasedTrkg::Trajectory"); + if (trackBank!=null) { + + for (int i = 0; i < trackBank.rows(); i++) { + MUVTTrack trk = new MUVTTrack(); + int id = trackBank.getShort("id", i); + trk.setId(id); + trk.setIndex(i); + trk.setSector(trackBank.getByte("sector", i)); + trk.setQ(trackBank.getByte("q", 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 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); + MUVTTrajectory trj = new MUVTTrajectory(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()) { + tracks.add(entry.getValue()); + } + return tracks; + } + + public Point3D transGlobaltoLocal(int sector, double x, double y, double z){ + Point3D point = new Point3D(x, y, z); + point.rotateZ(Math.toRadians(-60 * (sector - 1))); + point.rotateY(Math.toRadians(-25)); + + return point; + } + + public Point3D transLocaltoGlobal(int sector, double x, double y, double z){ + Point3D point = new Point3D(x, y, z); + point.rotateY(Math.toRadians(25)); + point.rotateZ(Math.toRadians(60 * (sector - 1))); + + return point; + } + + @Override + public String toString() { + String str = "FMT track :" + " Index " + this._index + + " Q " + this._q + + String.format(" P (%.4f,%.4f,%.4f)", this._px, this._py, this._pz) + + String.format(" D (%.4f,%.4f,%.4f)", this._x, this._y, this._z); + for(int i=0; i clusters; + + public KFitter(MUVTTrack track, Swim swimmer, int c) { + sv = new StateVecs(swimmer); + this.track = track; + this.clusters = track.getClusters(); + this.init(clusters, track.getSector(), track.getX(), track.getY(), track.getZ(), + track.getPx(), track.getPy(), track.getPz(), track.getQ(), c); + } + + 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 MUVTTrack(0,sector, q, xVtx, yVtx, zVtx, pxVtx, pyVtx, pzVtx, clusters); + this.clusters = clusters; + 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, + double pxVtx, double pyVtx, double pzVtx, int q, int c) { + // initialize measVecs + mv.setMeasVecs(clusters); + int mSize = mv.measurements.size(); + + sv.Z = new double[mSize]; + + 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); + } + + public void runFitter(int sector) { + int svzLength = sv.Z.length; + + for (int i = 1; i <= totNumIter; i++) { + interNum = i; + this.chi2 = 0; + if (i > 1) { + for (int k = svzLength - 1; k > 0; k--) { + if (k >= 1) { + sv.transport(sector, k, k - 1, sv.trackTraj.get(k), sv.trackCov.get(k)); + this.filter(k - 1); + } + } + } + for (int k = 0; k < svzLength - 1; k++) { + sv.transport(sector, k, k + 1, sv.trackTraj.get(k), sv.trackCov.get(k)); + this.filter(k + 1); + } + if (i > 1) { + if (this.setFitFailed) i = totNumIter; + if (!this.setFitFailed) { + this.finalStateVec = sv.trackTraj.get(svzLength - 1); + this.finalCovMat = sv.trackCov.get(svzLength - 1); + } else { + this.ConvStatus = 1; + } + } + } + if (totNumIter == 1) { + this.finalStateVec = sv.trackTraj.get(svzLength - 1); + this.finalCovMat = sv.trackCov.get(svzLength - 1); + } + + // Do one final pass to get the final chi^2 and the corresponding centroid residuals. + this.chi2 = 0; + for (int k = svzLength - 1; k > 0; --k) { + if (k >= 1) { + sv.transport(sector, k, k-1, sv.trackTraj.get(k), sv.trackCov.get(k)); + this.filter(k - 1); + } + } + for (int k = 0; k < svzLength - 1; ++k) { + sv.transport(sector, k, k+1, sv.trackTraj.get(k), sv.trackCov.get(k)); + } + + + // save final trajectory points + /* + if(this.finalStateVec!=null) { + for (int k = 0; k < svzLength; ++k) { + Trajectory trj = new Trajectory(mv.measurements.get(k).layer, + sv.trackTraj.get(k).x, + sv.trackTraj.get(k).y, + sv.trackTraj.get(k).z, + sv.trackTraj.get(k).tx, + sv.trackTraj.get(k).ty, + 0, + sv.trackTraj.get(k).deltaPath); + track.setFMTtraj(trj); + } + } + */ + + } + + public Matrix filterCovMat(double[] H, Matrix Ci, double V) { + + double det = Matrix5x5.inverse(Ci, first_inverse, adj); + if (Math.abs(det) < 1.e-60) { + return null; + } + + addition.set( + 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-60) { + return null; + } + + return result_inv; + } + + 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 * mv.measurements.get(k).error; + + double[] H = mv.HMUVT(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++) { + // 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 res = mv.dhMUVT(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; + sv.trackTraj.get(k).y = y_filt; + sv.trackTraj.get(k).tx = tx_filt; + sv.trackTraj.get(k).ty = ty_filt; + sv.trackTraj.get(k).Q = Q_filt; + } + } + } + + public Matrix propagateToVtx(int sector, double Zf) { + return sv.transport(sector, 0, Zf, sv.trackTraj.get(0), sv.trackCov.get(0)); + } +} diff --git a/reconstruction/mu/src/main/java/org/jlab/rec/muvt/track/fit/MeasVecs.java b/reconstruction/mu/src/main/java/org/jlab/rec/muvt/track/fit/MeasVecs.java new file mode 100644 index 0000000000..5cc75e942e --- /dev/null +++ b/reconstruction/mu/src/main/java/org/jlab/rec/muvt/track/fit/MeasVecs.java @@ -0,0 +1,129 @@ +package org.jlab.rec.muvt.track.fit; + +import java.util.ArrayList; +import java.util.List; +import org.jlab.rec.muvt.track.fit.StateVecs.StateVec; + +import org.jlab.geom.prim.Point3D; +import org.jlab.geom.prim.Line3D; +import org.jlab.geom.prim.Plane3D; +import org.jlab.rec.muvt.MUVTCluster; + +/** + * @author ziegler + */ + +public class MeasVecs { + + public List measurements ; + + public class MeasVec implements Comparable { + public double z = Double.NaN; + public Point3D lineEndPoint1 = null; + public Point3D lineEndPoint2 = null; + public double seed; + public double error; + public int layer; + public int k; + public int size; + + MeasVec() {} + + @Override + public int compareTo(MeasVec arg) { + int CompLay = this.layer < arg.layer ? -1 : this.layer == arg.layer ? 0 : 1; + return CompLay; + } + } + + public void setMeasVecs(List clusters) { + measurements = new ArrayList<>(); + + for (int i = 0; i < clusters.size(); i++) { + int l = clusters.get(i).getLayer()-1; + Point3D lineEndPoint1 = clusters.get(i).getLineTS().origin(); + Point3D lineEndPoint2 = clusters.get(i).getLineTS().end(); + double error = 0.0144; // = pitch/sqrt(12), where pitch is 500 um + double z = lineEndPoint1.z(); + int seed = clusters.get(i).getMaxStrip(); + MeasVec meas = this.setMeasVec(l, lineEndPoint1, lineEndPoint2, error, z, seed); + measurements.add(meas); + } + } + + public MeasVec setMeasVec(int l, Point3D lineEndPoint1, Point3D lineEndPoint2, double error, double z, int seed) { + + MeasVec meas = new MeasVec(); + meas.layer = l+1; + meas.error = error; + meas.z = z; + meas.seed = seed; + + meas.lineEndPoint1 = lineEndPoint1; + meas.lineEndPoint2 = lineEndPoint2; + + return meas; + } + + public double dhMUVT(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[] HMUVT(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] = dhMUVT(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] = dhMUVT(svc); + } + + 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); + SVplus.x = stateVec.x; + SVplus.y = stateVec.y; + SVplus.z = stateVec.z; + SVplus.tx = stateVec.tx; + SVplus.ty = stateVec.ty; + SVplus.Q = stateVec.Q; + + return SVplus; + } +} diff --git a/reconstruction/mu/src/main/java/org/jlab/rec/muvt/track/fit/RungeKutta.java b/reconstruction/mu/src/main/java/org/jlab/rec/muvt/track/fit/RungeKutta.java new file mode 100644 index 0000000000..3ccbf4e93d --- /dev/null +++ b/reconstruction/mu/src/main/java/org/jlab/rec/muvt/track/fit/RungeKutta.java @@ -0,0 +1,497 @@ +package org.jlab.rec.muvt.track.fit; + +import java.util.ArrayList; +import org.jlab.clas.swimtools.Swim; + +/** + * + * @author ziegler + */ +public class RungeKutta { + + private final float[] _b = new float[3]; + final double v = 0.0029979245; + private final ArrayList k1; + private final ArrayList k2; + private final ArrayList k3; + private final ArrayList k4; + private final ArrayList jk1; + private final ArrayList jk2; + private final ArrayList jk3; + private final ArrayList jk4; + + public RungeKutta() { + this.k1 = new ArrayList<>(4); + this.k2 = new ArrayList<>(4); + this.k3 = new ArrayList<>(4); + this.k4 = new ArrayList<>(4); + this.jk1 = new ArrayList<>(12); + this.jk2 = new ArrayList<>(12); + this.jk3 = new ArrayList<>(12); + this.jk4 = new ArrayList<>(12); + } + + public void SwimToZ(int sector, StateVecs.StateVec fVec, Swim dcSwim, double z0, float[] bf){ + + double stepSize = 0.5; + 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; + double z = fVec.z; + final double Zi = fVec.z; + double BatMeas = 0; + + while(Math.signum(z0 - Zi) *z 0) { +// C[2][2] += cov_txtx; +// C[2][3] += cov_txty; +// C[3][2] += cov_txty; +// C[3][3] += cov_tyty; +// } + + fVec.x = x; + fVec.y = y ; + fVec.z = z0+h; + fVec.tx = tx; + fVec.ty = ty; + fVec.Q = q; + fVec.B = Math.sqrt(_b[0]*_b[0]+_b[1]*_b[1]+_b[2]*_b[2]); + fVec.deltaPath = Math.sqrt((x0-x)*(x0-x)+(y0-y)*(y0-y)+h*h)+dPath; + fCov.covMat.set(C); + //System.out.println("Transported matrix"); + //Matrix5x5.show(fCov.covMat); + } + + + private double RK4(double k1, double k2, double k3, double k4, double h) { + return h/6*(k1 + 2*k2 +2*k3 + k4); + } + + private double Ax(double tx, double ty, double Bx, double By, double Bz) { + double C = Math.sqrt(1 + tx * tx + ty * ty); + return C * (ty * (tx * Bx + Bz) - (1 + tx * tx) * By); + } + private double Ay(double tx, double ty, double Bx, double By, double Bz) { + double C = Math.sqrt(1 + tx * tx + ty * ty); + return C * (-tx * (ty * By + Bz) + (1 + ty * ty) * Bx); + } + + private double delAx_deltx(double tx, double ty, double Bx, double By, double Bz) { + double C2 = 1 + tx * tx + ty * ty; + double C = Math.sqrt(1 + tx * tx + ty * ty); + double Ax = C * (ty * (tx * Bx + Bz) - (1 + tx * tx) * By); + double Ay = C * (-tx * (ty * By + Bz) + (1 + ty * ty) * Bx); + + return tx * Ax / C2 + C * (ty * Bx - 2 * tx * By); //delAx_deltx + } + private double delAx_delty(double tx, double ty, double Bx, double By, double Bz) { + double C2 = 1 + tx * tx + ty * ty; + double C = Math.sqrt(1 + tx * tx + ty * ty); + double Ax = C * (ty * (tx * Bx + Bz) - (1 + tx * tx) * By); + double Ay = C * (-tx * (ty * By + Bz) + (1 + ty * ty) * Bx); + + return ty * Ax / C2 + C * (tx * Bx + Bz); //delAx_delty + } + private double delAy_deltx(double tx, double ty, double Bx, double By, double Bz) { + double C2 = 1 + tx * tx + ty * ty; + double C = Math.sqrt(1 + tx * tx + ty * ty); + double Ax = C * (ty * (tx * Bx + Bz) - (1 + tx * tx) * By); + double Ay = C * (-tx * (ty * By + Bz) + (1 + ty * ty) * Bx); + + return tx * Ay / C2 + C * (-ty * By - Bz); //delAy_deltx + } + private double delAy_delty(double tx, double ty, double Bx, double By, double Bz) { + double C2 = 1 + tx * tx + ty * ty; + double C = Math.sqrt(1 + tx * tx + ty * ty); + double Ax = C * (ty * (tx * Bx + Bz) - (1 + tx * tx) * By); + double Ay = C * (-tx * (ty * By + Bz) + (1 + ty * ty) * Bx); + + return ty * Ay / C2 + C * (-tx * By + 2 * ty * Bx); //delAy_delty + } + + private void A(double tx, double ty, double Bx, double By, double Bz, double[] a) { + + double C = Math.sqrt(1 + tx * tx + ty * ty); + a[0] = C * (ty * (tx * Bx + Bz) - (1 + tx * tx) * By); + a[1] = C * (-tx * (ty * By + Bz) + (1 + ty * ty) * Bx); + } + + private void delA_delt(double tx, double ty, double Bx, double By, double Bz, double[] dela_delt) { + + double C2 = 1 + tx * tx + ty * ty; + double C = Math.sqrt(1 + tx * tx + ty * ty); + double Ax = C * (ty * (tx * Bx + Bz) - (1 + tx * tx) * By); + double Ay = C * (-tx * (ty * By + Bz) + (1 + ty * ty) * Bx); + + dela_delt[0] = tx * Ax / C2 + C * (ty * Bx - 2 * tx * By); //delAx_deltx + dela_delt[1] = ty * Ax / C2 + C * (tx * Bx + Bz); //delAx_delty + dela_delt[2] = tx * Ay / C2 + C * (-ty * By - Bz); //delAy_deltx + dela_delt[3] = ty * Ay / C2 + C * (-tx * By + 2 * ty * Bx); //delAy_delty + } + + private double deltx_deltx0_next(double q, double v, double tx1, double ty1, float b0, float b1, float b2, double deltx_deltx0_1, double delty_deltx0_1) { + return q*v*(delAx_deltx(tx1,ty1,b0,b1,b2)*(deltx_deltx0_1) + + delAx_delty(tx1,ty1,b0,b1,b2)*(delty_deltx0_1)); + } + + private double delty_deltx0_next(double q, double v, double tx1, double ty1, float b0, float b1, float b2, double deltx_deltx0_1, double delty_deltx0_1) { + return q*v*(delAy_deltx(tx1,ty1,b0,b1,b2)*(deltx_deltx0_1) + + delAy_delty(tx1,ty1,b0,b1,b2)*(delty_deltx0_1)); + } + + private double deltx_delty0_next(double q, double v, double tx1, double ty1, float b0, float b1, float b2, double deltx_delty0_1, double delty_delty0_1) { + return q*v*(delAx_delty(tx1,ty1,b0,b1,b2)*(deltx_delty0_1) + + delAx_delty(tx1,ty1,b0,b1,b2)*(delty_delty0_1)); + } + + private double delty_delty0_next(double q, double v, double tx1, double ty1, float b0, float b1, float b2, double deltx_delty0_1, double delty_delty0_1) { + return q*v*(delAy_delty(tx1,ty1,b0,b1,b2)*(deltx_delty0_1) + + delAy_delty(tx1,ty1,b0,b1,b2)*(delty_delty0_1)); + } + + private double deltx_delq0_next(double q, double v, double tx1, double ty1, float b0, float b1, float b2, double deltx_delq0_1, double delty_delq0_1) { + return v*Ax(tx1, ty1, b0, b1, b2) + + q*v*(delAx_deltx(tx1,ty1,b0,b1,b2)*(deltx_delq0_1) + + delAx_delty(tx1,ty1,b0,b1,b2)*(delty_delq0_1)); + } + + private double delty_delq0_next(double q, double v, double tx1, double ty1, float b0, float b1, float b2, double deltx_delq0_1, double delty_delq0_1) { + return v*Ay(tx1, ty1, b0, b1, b2) + + q*v*(delAy_deltx(tx1, ty1,b0,b1,b2)*(deltx_delq0_1) + + delAy_delty(tx1, ty1,b0,b1,b2)*(delty_delq0_1)); + } + + private void getRKn(int sector, ArrayList k1, ArrayList k2, double d, double x0, double y0, double z0, double tx0, double ty0, double q, float[] b) { + + double tx1 = k1.get(2); + double ty1 = k1.get(3); + + double x2 = tx0+d*tx1; + double y2 = ty0+d*ty1; + double tx2=q*v*Ax((tx0+d*tx1), (ty0+d*ty1), b[0], b[1], b[2]); + double ty2=q*v*Ay((tx0+d*tx1), (ty0+d*ty1), b[0], b[1], b[2]); + + k2.add(0, x2); + k2.add(1, y2); + k2.add(2, tx2); + k2.add(3, ty2); + } + + private void getjRKn(int sector, ArrayList k1, ArrayList jk1, ArrayList jk2, double d, double x0, double y0, double z0, double tx0, double ty0, double q, float[] _b, + double deltx_deltx0_0, double delty_deltx0_0, double deltx_delty0_0, double delty_delty0_0, double deltx_delq0_0, double delty_delq0_0) { + + double tx1 = k1.get(2); + double ty1 = k1.get(3); + + double delx_deltx0_1 = jk1.get(0); + double dely_deltx0_1 = jk1.get(1); + double delx_delty0_1 = jk1.get(2); + double dely_delty0_1 = jk1.get(3); + + double deltx_deltx0_1 = jk1.get(4); + double delty_deltx0_1 = jk1.get(5); + double deltx_delty0_1 = jk1.get(6); + double delty_delty0_1 = jk1.get(7); + + double delx_delq0_1 = jk1.get(8); + double dely_delq0_1 = jk1.get(9); + + double deltx_delq0_1 = jk1.get(10); + double delty_delq0_1 = jk1.get(11); + + double delx_deltx0_2 = deltx_deltx0_0+d*deltx_deltx0_1; + double dely_deltx0_2 = delty_deltx0_0+d*delty_deltx0_1; + double delx_delty0_2 = deltx_delty0_0+d*deltx_delty0_1; + double dely_delty0_2 = delty_delty0_0+d*delty_delty0_1; + + double deltx_deltx0_2 = this.deltx_deltx0_next(q,v,tx0+d*tx1,ty0+d*ty1,_b[0],_b[1],_b[2], + deltx_deltx0_0+d*deltx_deltx0_1,delty_deltx0_0+d*delty_deltx0_1); + double delty_deltx0_2 = this.delty_deltx0_next(q,v,tx0+d*tx1,ty0+d*ty1,_b[0],_b[1],_b[2], + deltx_deltx0_0+d*deltx_deltx0_1,delty_deltx0_0+d*delty_deltx0_1); + double deltx_delty0_2 = this.deltx_delty0_next(q,v,tx0+d*tx1,ty0+d*ty1,_b[0],_b[1],_b[2], + deltx_delty0_0+d*deltx_delty0_1,delty_delty0_0+d*delty_delty0_1); + double delty_delty0_2 = this.delty_delty0_next(q,v,tx0+d*tx1,ty0+d*ty1,_b[0],_b[1],_b[2], + deltx_delty0_0+d*deltx_delty0_1,delty_delty0_0+d*delty_delty0_1); + + double delx_delq0_2 = deltx_delq0_0+d*deltx_delq0_1; + double dely_delq0_2 = delty_delq0_0+d*delty_delq0_1; + + double deltx_delq0_2 = this.deltx_delq0_next(q,v,tx0+d*tx1,ty0+d*ty1,_b[0],_b[1],_b[2], + deltx_delq0_0+d*deltx_delq0_1,delty_delq0_0+d*delty_delq0_1); + double delty_delq0_2 = this.delty_delq0_next(q,v,tx0+d*tx1,ty0+d*ty1,_b[0],_b[1],_b[2], + deltx_delq0_0+d*deltx_delq0_1,delty_delq0_0+d*delty_delq0_1); + + jk2.add(0, delx_deltx0_2 ); + jk2.add(1, dely_deltx0_2 ); + jk2.add(2, delx_delty0_2 ); + jk2.add(3, dely_delty0_2 ); + + jk2.add(4, deltx_deltx0_2 ); + jk2.add(5, delty_deltx0_2 ); + jk2.add(6, deltx_delty0_2 ); + jk2.add(7, delty_delty0_2 ); + + jk2.add(8, delx_delq0_2 ); + jk2.add(9, dely_delq0_2 ); + + jk2.add(10, deltx_delq0_2 ); + jk2.add(11, delty_delq0_2 ); + } +} diff --git a/reconstruction/mu/src/main/java/org/jlab/rec/muvt/track/fit/StateVecs.java b/reconstruction/mu/src/main/java/org/jlab/rec/muvt/track/fit/StateVecs.java new file mode 100644 index 0000000000..d5155f2a6f --- /dev/null +++ b/reconstruction/mu/src/main/java/org/jlab/rec/muvt/track/fit/StateVecs.java @@ -0,0 +1,508 @@ +package org.jlab.rec.muvt.track.fit; + +import org.jlab.jnp.matrix.*; +import java.util.HashMap; +import java.util.Map; +import org.jlab.clas.pdg.PhysicsConstants; +import org.jlab.clas.swimtools.Swim; + +/** + * + * @author ziegler + */ +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; + public Map trackTraj = new HashMap(); + public Map trackCov = new HashMap(); + + public StateVec StateVec; + public CovMat CovMat; + public Matrix F = new Matrix(); + private final Matrix fMS = new Matrix(); + private final Matrix copyMatrix = new Matrix(); + private final double[] A = new double[2]; + private final double[] dA = new double[4]; + private final float[] bf = new float[3]; + 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 + */ + public StateVecs(Swim swimmer) { + dcSwim = swimmer; + rk = new RungeKutta(); + } + + /** + * + * @param sector + * @param i initial state vector index + * @param Zf + * @param iVec state vector at the initial index + * @param covMat state covariance matrix at the initial index + * @return + */ + public Matrix transport(int sector, int i, double Zf, StateVec iVec, CovMat covMat) { // s = signed step-size + double stepSize = 1.0; + StateVecs.StateVec fVec = new StateVec(0); + CovMat fCov = new CovMat(0); + fVec.x = iVec.x; + fVec.y = iVec.y; + fVec.z = iVec.z; + fVec.tx = iVec.tx; + fVec.ty = iVec.ty; + fVec.Q = iVec.Q; + fVec.B = iVec.B; + + Matrix5x5.copy(covMat.covMat, fCov.covMat); + double s = 0; + double z = Z[i]; + 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, Q, x, y, z, tx, ty, s, dcSwim, + covMat, fVec, fCov, dPath); + + // Q process noise matrix estimate + + double p = Math.abs(1. / iVec.Q); + + 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 energy = Math.sqrt(p*p + MuonMass * MuonMass); + double beta = p/energy; + if(beta>1.0 || beta<=0) + beta =1.0; + + double sctRMS = 0; + + if(Math.abs(s)>0) + sctRMS = ((0.0136)/(beta*PhysicsConstants.speedOfLight()*p))*Math.sqrt(t_ov_X0)* + (1 + 0.038 * Math.log(t_ov_X0)); + + double cov_txtx = (1 + tx * tx) * (1 + tx * tx + ty * ty) * sctRMS * sctRMS; + double cov_tyty = (1 + ty * ty) * (1 + tx * tx + ty * ty) * sctRMS * sctRMS; + double cov_txty = tx * ty * (1 + tx * tx + ty * ty) * sctRMS * sctRMS; + + fMS.set( + 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, + 0, 0, cov_txtx, cov_txty, 0, + 0, 0, cov_txty, cov_tyty, 0, + 0, 0, 0, 0, 0 + ); + + Matrix5x5.copy(fCov.covMat, copyMatrix); + Matrix5x5.add(copyMatrix, fMS, fCov.covMat); + + // end add process noise + if (Math.abs(fVec.B - BatMeas) < 0.0001) stepSize*=2; + BatMeas = fVec.B; + } + + return fCov.covMat; + } + + /** + * + * @param sector + * @param i initial state vector index + * @param f final state vector index + * @param iVec state vector at the initial index + * @param covMat state covariance matrix at the initial index + */ + public void transport(int sector, int i, int f, StateVec iVec, CovMat covMat) { // s = signed step-size + if(iVec==null) + return; + double stepSize = 1.0; + StateVecs.StateVec fVec = new StateVec(f); + CovMat fCov = new CovMat(f); + fVec.x = iVec.x; + fVec.y = iVec.y; + fVec.z = iVec.z; + fVec.tx = iVec.tx; + fVec.ty = iVec.ty; + fVec.Q = iVec.Q; + fVec.B = iVec.B; + //fCov.covMat = covMat.covMat; + Matrix5x5.copy(covMat.covMat, fCov.covMat); + double s = 0; + double z = Z[i]; + 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); + + rk.RK4transport(sector, Q, x, y, z, tx, ty, s, dcSwim, + covMat, fVec, fCov, dPath); + + // Q process noise matrix estimate + + double p = Math.abs(1. / iVec.Q); + + 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; + if(beta>1.0 || beta<=0) + beta =1.0; + + double sctRMS = 0; + + if(Math.abs(s)>0) + sctRMS = ((0.0136)/(beta*PhysicsConstants.speedOfLight()*p))*Math.sqrt(t_ov_X0)* + (1 + 0.038 * Math.log(t_ov_X0)); + + + double cov_txtx = (1 + tx * tx) * (1 + tx * tx + ty * ty) * sctRMS * sctRMS; + double cov_tyty = (1 + ty * ty) * (1 + tx * tx + ty * ty) * sctRMS * sctRMS; + double cov_txty = tx * ty * (1 + tx * tx + ty * ty) * sctRMS * sctRMS; + + + //if (Math.signum(Z[f] - Z[i]) > 0) { + fMS.set( + 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, + 0, 0, cov_txtx, cov_txty, 0, + 0, 0, cov_txty, cov_tyty, 0, + 0, 0, 0, 0, 0 + ); + + Matrix5x5.copy(fCov.covMat, copyMatrix); + Matrix5x5.add(copyMatrix, fMS, fCov.covMat); + + //} + // end add process noise + + if( Math.abs(fVec.B - BatMeas)<0.0001) + stepSize*=2; + + BatMeas = fVec.B; + } + this.trackTraj.put(f, fVec); + this.trackCov.put(f, fCov); + } + double AIRRADLEN = 30400; // radiation length in cm + public double getX0(double z) { + + return AIRRADLEN; + } + /** + * + * @param i initial state vector index + * @param f final state vector index + * @param iVec state vector at the initial index + * @param covMat state covariance matrix at the initial index + */ + public void transportFixed(int sector, int i, int f, StateVec iVec, CovMat covMat) { // s = signed step-size + if(iVec==null) + return; + double stepSize = 0.5; + + StateVecs.StateVec fVec = new StateVec(f); + StateVecs.CovMat fCov = new CovMat(f); + fVec.x = iVec.x; + fVec.y = iVec.y; + fVec.z = iVec.z; + fVec.tx = iVec.tx; + fVec.ty = iVec.ty; + fVec.Q = iVec.Q; + fCov.covMat = covMat.covMat; + int nSteps = (int) (Math.abs((Z[i] - Z[f]) / stepSize) + 1); + + double s = (Z[f] - Z[i]) / (double) nSteps; + double z = Z[i]; + + for (int j = 0; j < nSteps; j++) { + // get the sign of the step + if (j == nSteps - 1) { + s = Math.signum(Z[f] - Z[i]) * Math.abs(z - Z[f]); + } + //System.out.println(" RK step num "+(j+1)+" = "+(float)s+" nSteps = "+nSteps); + double x = fVec.x; + double y = fVec.y; + z = fVec.z; + double tx = fVec.tx; + double ty = fVec.ty; + double Q = fVec.Q; + double dPath = fVec.deltaPath; + covMat.covMat = fCov.covMat; + + rk.RK4transport(sector, Q, x, y, z, tx, ty, s, dcSwim, + covMat, fVec, fCov, dPath); + + } + + this.trackTraj.put(f, fVec); + this.trackCov.put(f, fCov); + } + + + private void A(double tx, double ty, double Bx, double By, double Bz, double[] a) { + + double C = Math.sqrt(1 + tx * tx + ty * ty); + a[0] = C * (ty * (tx * Bx + Bz) - (1 + tx * tx) * By); + a[1] = C * (-tx * (ty * By + Bz) + (1 + ty * ty) * Bx); + } + + private void delA_delt(double tx, double ty, double Bx, double By, double Bz, double[] dela_delt) { + + double C2 = 1 + tx * tx + ty * ty; + double C = Math.sqrt(1 + tx * tx + ty * ty); + double Ax = C * (ty * (tx * Bx + Bz) - (1 + tx * tx) * By); + double Ay = C * (-tx * (ty * By + Bz) + (1 + ty * ty) * Bx); + + dela_delt[0] = tx * Ax / C2 + C * (ty * Bx - 2 * tx * By); //delAx_deltx + dela_delt[1] = ty * Ax / C2 + C * (tx * Bx + Bz); //delAx_delty + dela_delt[2] = tx * Ay / C2 + C * (-ty * By - Bz); //delAy_deltx + dela_delt[3] = ty * Ay / C2 + C * (-tx * By + 2 * ty * Bx); //delAy_delty + } + + + private double beta = 1.0; + + /** + * + * @param sector + * @param xVtx + * @param yVtx + * @param zVtx + * @param z0 the value in z to which the track is swam back to + * @param pyVtx + * @param pzVtx + * @param q + * @param kf the final state measurement index + * @param pxVtx + * @param c + */ + public void init(int sector, double xVtx, double yVtx, double zVtx, + double pxVtx, double pyVtx, double pzVtx, + int q, + double z0, KFitter kf, int c) { + + StateVec initSV = new StateVec(0); + initSV.x = xVtx; + initSV.y = yVtx; + initSV.z = zVtx; + initSV.tx = pxVtx/pzVtx; + initSV.ty = pyVtx/pzVtx; + double p = Math.sqrt(pxVtx*pxVtx+pyVtx*pyVtx+pzVtx*pzVtx); + initSV.Q = (double)q / p; + + rk.SwimToZ(sector, initSV, dcSwim, z0, bf); + + if (initSV != null) { + + this.trackTraj.put(0, initSV); + } else { + kf.setFitFailed = true; + return; + } + + CovMat initCM = new CovMat(0); + StateVec rinitSV = new StateVec(0); + rinitSV.x = xVtx; + rinitSV.y = yVtx; + rinitSV.z = 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++) { + FTF[i*5+i]=F[i]*F[i]; + } + 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) { + SVplus = new StateVec(stateVec.k); + SVplus.x = stateVec.x; + SVplus.y = stateVec.y; + SVplus.tx = stateVec.tx; + SVplus.ty = stateVec.ty; + SVplus.z = stateVec.z; + SVplus.Q = stateVec.Q; + + return SVplus; + } + private void swimToSite(int sector, double z0, + StateVec SVplus, StateVec SVminus) { + + rk.SwimToZ(sector, SVplus, dcSwim, z0, bf); + rk.SwimToZ(sector, SVminus, dcSwim, z0, bf); + } + + double[] F(int sector, double z0, StateVec stateVec) { + double[] _F = new double[5]; + StateVec SVplus = null; + StateVec SVminus = null; + + SVplus = this.reset(SVplus, stateVec); + SVminus = this.reset(SVminus, stateVec); + + double delt_x = 0.05; + SVplus.x += delt_x/2.; + SVminus.x-= delt_x/2.; + + this.swimToSite(sector, z0, SVplus, SVminus); + + _F[0] = (SVplus.x - SVminus.x)/delt_x; + + SVplus = this.reset(SVplus, stateVec); + SVminus = this.reset(SVminus, stateVec); + + double delt_y = 0.05; + SVplus.y += delt_y/2.; + SVminus.y-= delt_y/2.; + + this.swimToSite(sector, z0, SVplus, SVminus); + + _F[1] = (SVplus.y - SVminus.y)/delt_y; + + SVplus = this.reset(SVplus, stateVec); + SVminus = this.reset(SVminus, stateVec); + + double delt_tx = 0.001; + SVplus.tx += delt_tx/2.; + SVminus.tx-= delt_tx/2.; + + this.swimToSite(sector, z0, SVplus, SVminus); + + _F[2] = (SVplus.tx - SVminus.tx)/delt_tx; + + SVplus = this.reset(SVplus, stateVec); + SVminus = this.reset(SVminus, stateVec); + + double delt_ty = 0.001; + SVplus.ty += delt_ty/2.; + SVminus.ty-= delt_ty/2.; + + this.swimToSite(sector, z0, SVplus, SVminus); + + _F[3] = (SVplus.ty - SVminus.ty)/delt_ty; + + SVplus = this.reset(SVplus, stateVec); + SVminus = this.reset(SVminus, stateVec); + + + _F[4] = 0.01/Math.abs(SVplus.Q); + + return _F; + + } + + public void printMatrix(Matrix C) { + for (int k = 0; k < 5; k++) { + for (int j = 0; j < 5; j++) { + System.out.println("C["+j+"]["+k+"] = "+C.get(j, k)); + } + } + } + + /** + * The state vector representing the track at a given measurement site + */ + public class StateVec { + + final int k; //index + public double z; //z (fixed measurement planes) + public double x; //track x in the tilted sector coordinate system at z + public double y; //track y in the tilted sector coordinate system at z + public double tx; //track px/pz in the tilted sector coordinate system at z + public double ty; //track py/pz in the tilted sector coordinate system at z + public double Q; //track q/p + double B; + double deltaPath; + + StateVec(int k) { + this.k = k; + } + + public double getPx() { + return this.getPz()*tx; + } + + public double getPy() { + return this.getPz()*ty; + } + + public double getPz() { + double pz = this.getP()/Math.sqrt(tx*tx+ty*ty+1); + return pz; + } + + public double getP() { + return 1./Math.abs(this.Q); + } + + + + String printInfo() { + return this.k+"] = "+(float)this.x+", "+(float)this.y+", "+(float)this.z+", " + +(float)this.tx+", "+(float)this.ty+", "+(float)this.Q+" B = "+(float)this.B; + } + } + /** + * The track covariance matrix + */ + public class CovMat { + + final int k; + public Matrix covMat = new Matrix(); + + CovMat(int k) { + this.k = k; + } + + } +} diff --git a/reconstruction/pom.xml b/reconstruction/pom.xml index 134042683f..645e0d8120 100644 --- a/reconstruction/pom.xml +++ b/reconstruction/pom.xml @@ -40,6 +40,7 @@ recoil calib uber + mu