diff --git a/common-tools/clas-analysis/src/main/java/org/jlab/analysis/eventmerger/ADCTDCMerger.java b/common-tools/clas-analysis/src/main/java/org/jlab/analysis/eventmerger/ADCTDCMerger.java index 290bf12cfb..11788285c3 100644 --- a/common-tools/clas-analysis/src/main/java/org/jlab/analysis/eventmerger/ADCTDCMerger.java +++ b/common-tools/clas-analysis/src/main/java/org/jlab/analysis/eventmerger/ADCTDCMerger.java @@ -86,7 +86,8 @@ public List readADCs(DetectorType detector,DataBank bank) { ADC adcData = new ADC(detector); adcData.readFromBank(bank, i); - if(!adcData.isGood()) adcData.skip(); + if(!adcData.isGood() || + (detector==DetectorType.FTCAL && adcData.getAdc()<200)) adcData.skip(); adcStore.add(adcData); } diff --git a/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/DCGeant4Factory.java b/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/DCGeant4Factory.java index 28e5782faf..236af2503a 100644 --- a/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/DCGeant4Factory.java +++ b/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/DCGeant4Factory.java @@ -580,6 +580,7 @@ public final class DCGeant4Factory extends Geant4Factory { private final HashMap properties = new HashMap<>(); private int nsgwires; + private final double x_enlargement = 1.00; private final double y_enlargement = 3.65; private final double z_enlargement = -2.46; private final double microgap = 0.01; @@ -825,19 +826,20 @@ private Geant4Basic getRegion(int isec, int ireg) { } private Geant4Basic getLayer(int isec, int isuper, int ilayer) { - return getRegion(isec, isuper/2).getChildren().get((isuper%2)*6 + ilayer); + return getRegion(isec, isuper/2).getChildren().get(0).getChildren().get((isuper%2)*6 + ilayer); } /////////////////////////////////////////////////// public Geant4Basic createRegion(int isector, int iregion) { Wire regw0 = new Wire(isector+1, iregion * 2, 0, 0); Wire regw1 = new Wire(isector+1, iregion * 2 + 1, 7, nsgwires - 1); - double dx_shift = y_enlargement * Math.tan(Math.toRadians(29.5)); - + double dx_shift = y_enlargement * Math.tan(dbref.thopen(iregion)/2); + double extra = 0.2; + double reg_dz = (dbref.frontgap(iregion) + dbref.backgap(iregion) + dbref.midgap(iregion) + dbref.superwidth(iregion * 2) + dbref.superwidth(iregion * 2 + 1)) / 2.0 + z_enlargement; - double reg_dx0 = Math.abs(regw0.bottom().x) - dx_shift + 1.0; - double reg_dx1 = Math.abs(regw1.top().x) + dx_shift + 1.0; - double reg_dy = regw1.top().minus(regw0.bottom()).y / Math.cos(dbref.thtilt(iregion)) / 2.0 + y_enlargement + 1.0; + double reg_dx0 = Math.abs(regw0.bottom().x) - dx_shift + x_enlargement; + double reg_dx1 = Math.abs(regw1.top().x) + dx_shift + x_enlargement; + double reg_dy = regw1.top().minus(regw0.bottom()).y / Math.cos(dbref.thtilt(iregion)) / 2.0 + y_enlargement; double reg_skew = 0.0; double reg_thtilt = dbref.thtilt(iregion); @@ -854,6 +856,17 @@ public Geant4Basic createRegion(int isector, int iregion) { regionVolume.translate(vcenter.x, vcenter.y, vcenter.z); regionVolume.setId(isector + 1, iregion + 1, 0, 0); + double gas_dx0 = reg_dx0 - x_enlargement + dx_shift + extra; + double gas_dx1 = reg_dx1 - x_enlargement - dx_shift + extra; + double gas_dy = reg_dy - y_enlargement + extra; + Geant4Basic regionGas = new G4Trap("regionGas" + (iregion + 1) + "_s" + (isector + 1), + reg_dz, -reg_thtilt, Math.toRadians(90.0), + gas_dy, gas_dx0, gas_dx1, 0.0, + gas_dy, gas_dx0, gas_dx1, 0.0); + regionGas.setPosition(0, 0, 0); + regionGas.setMother(regionVolume); + regionGas.setId(isector + 1, iregion + 1, 0, 0); + for (int isup = 0; isup < 2; isup++) { int isuper = iregion * 2 + isup; Geant4Basic superlayerVolume = this.createSuperlayer(isuper); @@ -866,7 +879,7 @@ public Geant4Basic createRegion(int isector, int iregion) { superlayerVolume.rotate("zxy", -dbref.thster(isuper), 0.0, 0.0); superlayerVolume.setPosition(slshift.x, slshift.y, slshift.z); - superlayerVolume.setMother(regionVolume); + superlayerVolume.setMother(regionGas); superlayerVolume.setId(isector + 1, iregion + 1, isuper + 1); int nsglayers = dbref.nsenselayers(isuper) + dbref.nguardlayers(isuper); @@ -881,7 +894,7 @@ public Geant4Basic createRegion(int isector, int iregion) { layerVolume.rotate("zxy", -dbref.thster(isuper), 0.0, 0.0); layerVolume.setPosition(lshift.x, lshift.y, lshift.z); - layerVolume.setMother(regionVolume); + layerVolume.setMother(regionGas); layerVolume.setId(isector + 1, iregion + 1, isuper + 1, ilayer); } } diff --git a/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellConstants.java b/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellConstants.java index 81b08bf569..c9ef54163d 100644 --- a/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellConstants.java +++ b/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellConstants.java @@ -9,23 +9,30 @@ public class URWellConstants { private final static String CCDBPATH = "/geometry/urwell/"; - public final static int NMAXREGIONS = 2; //max number of regions - public final static int NREGIONS = 1; //number of regions + public final static int NMAXREGIONS = 6; //max number of regions + public final static int NREGIONS = 6; //number of regions public final static int NSECTORS = 6; //number of sectors public final static int NLAYERS = 2; //number of layers public final static int NCHAMBERS = 3; //number of chambers in a sector + public final static int NCHAMBERS_ddvcs = 1; //number of chambers in a sector - public final static double XENLARGEMENT = 0.5; // cm - public final static double YENLARGEMENT = 1.; // cm + public final static double XENLARGEMENT = 0.1; // cm + public final static double YENLARGEMENT = 0.1; // cm public final static double ZENLARGEMENT = 0.1; // cm // Sector geometrical parameters public final static double THOPEN = 54.; // opening angle between endplate planes (deg) public final static double THTILT = 25; // theta tilt (deg) public final static double THMIN = 4.694; // polar angle to the base of first chamber (deg) + public final static double THMAX = 40; public final static double SECTORHEIGHT = 146.21; //height of each sector (cm) public final static double DX0CHAMBER0 = 5.197; // halfbase of chamber 1 (cm) + public final static double SECTORHEIGHT_ddvcs = 22; + public final static double THMIN_ddvcs = 6.8; + public final static double DX0CHAMBER0_ddvcs = 3.207; + public final static double THMAX_ddvcs = 35; + // Chamber volumes and materials (units are cm) public final static double[] CHAMBERVOLUMESTHICKNESS = {0.0025, 0.0005,0.3, // window 0.0025, 0.0005,0.4, // cathode @@ -43,19 +50,32 @@ public class URWellConstants { "support_skin1_g10", "support_honeycomb_nomex", "support_skin2_g10"}; // URWELL position in the CLAS12 frame - public final static double TGT2DC0 = 228.078; // cm + public final static double TGT2DC0 = 228.078; // cm // public final static double URWELL2DC0 = 2; // cm public final static double URWELL2DC0[] = new double[NMAXREGIONS]; public final static double DIST2TGT[] = new double[NMAXREGIONS]; public final static double W2TGT[] = new double[NMAXREGIONS];; public final static double YMIN[] = new double[NMAXREGIONS]; public final static double ZMIN[] = new double[NMAXREGIONS]; + public final static double HSECTOR[] = new double[NMAXREGIONS]; + + + public final static double TGT2DC0_ddvcs = 48.9; // cm + // public final static double URWELL2DC0 = 2; // cm + public final static double URWELL2DC0_ddvcs[] = new double[NMAXREGIONS]; + public final static double DIST2TGT_ddvcs[] = new double[NMAXREGIONS]; + public final static double W2TGT_ddvcs[] = new double[NMAXREGIONS]; + ; + public final static double YMIN_ddvcs[] = new double[NMAXREGIONS]; + public final static double ZMIN_ddvcs[] = new double[NMAXREGIONS]; + public final static double HSECTOR_ddvcs[] = new double[NMAXREGIONS]; // public final static double DIST2TGT = (TGT2DC0-URWELL2DC0); // public final static double W2TGT = DIST2TGT/Math.cos(Math.toRadians(THTILT-THMIN)); // public final static double YMIN = W2TGT*Math.sin(Math.toRadians(THMIN)); // distance from the base chamber1 and beamline // public final static double ZMIN = W2TGT*Math.cos(Math.toRadians(THMIN)); public final static double PITCH = 0.1 ; // cm + public final static double PITCH_ddvcs = 0.05 ; // cm public final static double STEREOANGLE = 10; // deg @@ -121,8 +141,21 @@ public static synchronized void load( DatabaseConstantProvider cp ) W2TGT[i] = DIST2TGT[i]/Math.cos(Math.toRadians(THTILT-THMIN)); YMIN[i]= W2TGT[i]*Math.sin(Math.toRadians(THMIN)); // distance from the base chamber1 and beamline ZMIN[i] = W2TGT[i]*Math.cos(Math.toRadians(THMIN)); - + + double h1 = DIST2TGT[i]*Math.tan((Math.toRadians(THMAX -THTILT ))); + double h2 = DIST2TGT[i]*Math.tan((Math.toRadians(THTILT - THMIN))); + HSECTOR[i] = h1 +h2; + + URWELL2DC0_ddvcs[i] = -2. + i * 2; + DIST2TGT_ddvcs[i] = (TGT2DC0_ddvcs + URWELL2DC0_ddvcs[i]); + W2TGT_ddvcs[i] = DIST2TGT_ddvcs[i] / Math.cos(Math.toRadians(THTILT - THMIN_ddvcs)); + YMIN_ddvcs[i] = W2TGT_ddvcs[i] * Math.sin(Math.toRadians(THMIN_ddvcs)); // distance from the base chamber1 and beamline + ZMIN_ddvcs[i] = W2TGT_ddvcs[i] * Math.cos(Math.toRadians(THMIN_ddvcs)); + double h1_ddvcs = DIST2TGT_ddvcs[i]*Math.tan((Math.toRadians(THMAX_ddvcs -THTILT ))); + double h2_ddvcs = DIST2TGT_ddvcs[i]*Math.tan((Math.toRadians(THTILT - THMIN_ddvcs))); + HSECTOR_ddvcs[i] = h1_ddvcs +h2_ddvcs; + } diff --git a/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellGeant4Factory.java b/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellGeant4Factory.java index deaac1cf2d..66bc6f43ab 100644 --- a/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellGeant4Factory.java +++ b/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellGeant4Factory.java @@ -8,6 +8,9 @@ import org.jlab.detector.calib.utils.DatabaseConstantProvider; import org.jlab.geom.prim.Line3D; import org.jlab.geom.prim.Point3D; +import java.util.Arrays; + + /** @@ -20,8 +23,15 @@ public final class URWellGeant4Factory extends Geant4Factory { private int nRegions = URWellConstants.NREGIONS; private int nSectors = URWellConstants.NSECTORS; private int nChambers = URWellConstants.NCHAMBERS; + private double sectorheight = URWellConstants.SECTORHEIGHT; + private double thilt = URWellConstants.THTILT; + private double dx0chamber0 = URWellConstants.DX0CHAMBER0; + private double thopen = URWellConstants.THOPEN; private boolean isProto = false; - + private double[] ymin = new double[URWellConstants.NMAXREGIONS]; + private double[] zmin = new double[URWellConstants.NMAXREGIONS]; + private double[] Hsector = new double[URWellConstants.NMAXREGIONS]; + private String name = ""; /** * Create the URWELL full geometry @@ -34,19 +44,106 @@ public URWellGeant4Factory( DatabaseConstantProvider cp, boolean prototype, int this.init(cp, prototype, nRegions); } + /** + * + * @param config is 0 standard, 1 ddvcs proposal, 2 prototype + */ + + public URWellGeant4Factory( DatabaseConstantProvider cp, int config, int nRegions) { + URWellConstants.connect(cp ); + this.init(cp, config, nRegions); + } + + + public void init(DatabaseConstantProvider cp, int config, int regions ) { + + motherVolume = new G4World("root"); + + switch (config) { + case 0 -> { + + nSectors = URWellConstants.NSECTORS; + nChambers = URWellConstants.NCHAMBERS; + nRegions = Math.min(URWellConstants.NMAXREGIONS, regions); + isProto = false; + sectorheight = URWellConstants.SECTORHEIGHT; + Arrays.fill(Hsector, URWellConstants.SECTORHEIGHT); + thilt = URWellConstants.THTILT; + dx0chamber0 = URWellConstants.DX0CHAMBER0; + thopen = URWellConstants.THOPEN; + ymin =Arrays.copyOf(URWellConstants.YMIN, URWellConstants.YMIN.length); + zmin =Arrays.copyOf(URWellConstants.ZMIN, URWellConstants.ZMIN.length); + name =""; + + } + case 1 -> { + nRegions = URWellConstants.NREGIONS_PROTO; + nSectors = URWellConstants.NSECTORS_PROTO; + nChambers = URWellConstants.NCHAMBERS_PROTO; + isProto = true; + name = "_proto"; + + } + case 2 -> { + + nRegions = Math.min(URWellConstants.NMAXREGIONS, regions); + nSectors = URWellConstants.NSECTORS; + nChambers = URWellConstants.NCHAMBERS_ddvcs; + Hsector = Arrays.copyOf(URWellConstants.HSECTOR_ddvcs, URWellConstants.HSECTOR_ddvcs.length); + sectorheight = URWellConstants.SECTORHEIGHT_ddvcs; + thilt = URWellConstants.THTILT; + dx0chamber0 = URWellConstants.DX0CHAMBER0_ddvcs; + thopen = URWellConstants.THOPEN; + ymin = Arrays.copyOf(URWellConstants.YMIN_ddvcs, URWellConstants.YMIN_ddvcs.length); + zmin = Arrays.copyOf(URWellConstants.ZMIN_ddvcs, URWellConstants.ZMIN_ddvcs.length); + name = "_ddvcs"; + isProto = false; + } + default -> { + } + } + + for (int iregion = 0; iregion (volume.getName().contains(volumeName))) - .findAny() - .orElse(null); + + volumeName = "rg" + r + "_s" + s + "_c" + c + "_cathode_gas" + name; + + return this.getAllVolumes().stream() + .filter(volume -> (volume.getName().contains(volumeName))) + .findAny() + .orElse(null); } /** @@ -399,7 +476,8 @@ public Geant4Basic getSectorVolume(int region, int sector) { int r = region; int s = sector; - String volName = "region_uRwell_" + r + "_s" + s; + String volName = "region_uRwell_" + r + "_s" + s +name; + return this.getAllVolumes().stream() .filter(volume -> (volume.getName().contains(volName))) .findAny() @@ -415,8 +493,8 @@ public static void main(String[] args) { URWellConstants.connect(cp); - URWellGeant4Factory factory = new URWellGeant4Factory(cp, false, 2); - + // URWellGeant4Factory factory = new URWellGeant4Factory(cp, false, 2); + URWellGeant4Factory factory = new URWellGeant4Factory(cp, 2, 6); factory.getAllVolumes().forEach(volume -> { System.out.println(volume.gemcString()); }); diff --git a/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellStripFactory.java b/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellStripFactory.java index a7ec2edb5f..8cf110c400 100644 --- a/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellStripFactory.java +++ b/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/URWELL/URWellStripFactory.java @@ -1,6 +1,6 @@ package org.jlab.detector.geant4.v2.URWELL; - +import eu.mihosoft.vrl.v3d.Transform; import eu.mihosoft.vrl.v3d.Vector3d; import java.util.List; import org.jlab.detector.calib.utils.DatabaseConstantProvider; @@ -11,51 +11,60 @@ import org.jlab.geom.prim.Point3D; import org.jlab.geom.prim.Vector3D; import org.jlab.geometry.prim.Line3d; +import org.jlab.geom.prim.Trap3D; import org.jlab.geometry.prim.Straight; import org.jlab.utils.groups.IndexedList; /** * Creates and handles the URWELL detector strips as 3D lines - * + * * @author bondi */ public final class URWellStripFactory { private URWellGeant4Factory factory; - private IndexedList globalStrips = new IndexedList(3); - private IndexedList localStrips = new IndexedList(3); - private IndexedList planeStrips = new IndexedList(3); + private IndexedList globalStrips = new IndexedList(3); + private IndexedList localStrips = new IndexedList(3); + private IndexedList planeStrips = new IndexedList(3); + private IndexedList surfaceLayers = new IndexedList(2); + + private int nRegions; private int nSectors; private int nChambers; private int nLayers; + private int Nconfig; + private double pitch; private boolean isProto; - + public URWellStripFactory() { } - + /** - * Create the strip factory based on constants from CCDB. - * Currently constants are defined in the URWellConstants class. - * They will be moved to CCDB when finalized). + * Create the strip factory based on constants from CCDB. Currently + * constants are defined in the URWellConstants class. They will be moved to + * CCDB when finalized). + * * @param cp database provide */ public URWellStripFactory(DatabaseConstantProvider cp) { this.init(cp); } - + /** * Initialize the factory by the strip maps + * * @param cp */ public void init(DatabaseConstantProvider cp) { this.init(cp, false, 1); } - + /** - * Create the strip factory based on constants from CCDB. - * Currently constants are defined in the URWellConstants class. - * They will be moved to CCDB when finalized). + * Create the strip factory based on constants from CCDB. Currently + * constants are defined in the URWellConstants class. They will be moved to + * CCDB when finalized). + * * @param cp database provide * @param prototype * @param regions @@ -63,9 +72,23 @@ public void init(DatabaseConstantProvider cp) { public URWellStripFactory(DatabaseConstantProvider cp, boolean prototype, int regions) { this.init(cp, prototype, regions); } - + + /** + * Create the strip factory based on constants from CCDB. Currently + * constants are defined in the URWellConstants class. They will be moved to + * CCDB when finalized). + * + * @param cp database provide + * @param config 0 standard, 1 ddvcs proposal, 2 prototype + * @param regions + */ + public URWellStripFactory(DatabaseConstantProvider cp, int config, int regions) { + this.init(cp, config, regions); + } + /** * Initialize the factory by the strip maps + * * @param cp * @param prototype * @param regions @@ -73,49 +96,134 @@ public URWellStripFactory(DatabaseConstantProvider cp, boolean prototype, int re public void init(DatabaseConstantProvider cp, boolean prototype, int regions) { factory = new URWellGeant4Factory(cp, prototype, regions); isProto = prototype; - if(!isProto){ - nRegions = Math.min(URWellConstants.NMAXREGIONS, regions); - nSectors = URWellConstants.NSECTORS; + if (!isProto) { + nRegions = Math.min(URWellConstants.NMAXREGIONS, regions); + nSectors = URWellConstants.NSECTORS; nChambers = URWellConstants.NCHAMBERS; - nLayers = URWellConstants.NLAYERS; - } - else { - nRegions = URWellConstants.NREGIONS_PROTO; - nSectors = URWellConstants.NSECTORS_PROTO; + nLayers = URWellConstants.NLAYERS; + pitch = URWellConstants.PITCH; + } else { + nRegions = URWellConstants.NREGIONS_PROTO; + nSectors = URWellConstants.NSECTORS_PROTO; nChambers = URWellConstants.NCHAMBERS_PROTO; - nLayers = URWellConstants.NLAYERS; + nLayers = URWellConstants.NLAYERS; + pitch = URWellConstants.PITCH; } this.fillStripLists(); this.fillPlaneLists(); + this.fillSurfaceLists(); + } + + /** + * Initialize the factory by the strip maps + * + * @param cp + * @param config 0 standard, 1 ddvcs proposal, 2 prototype + * @param regions + */ + public void init(DatabaseConstantProvider cp, int config, int regions) { + factory = new URWellGeant4Factory(cp, config, regions); + + switch (config) { + case 0 -> { + + nSectors = URWellConstants.NSECTORS; + nChambers = URWellConstants.NCHAMBERS; + nRegions = Math.min(URWellConstants.NMAXREGIONS, regions); + nLayers = URWellConstants.NLAYERS; + isProto = false; + Nconfig = config; + pitch = URWellConstants.PITCH; + } + case 1 -> { + + nRegions = URWellConstants.NREGIONS_PROTO; + nSectors = URWellConstants.NSECTORS_PROTO; + nChambers = URWellConstants.NCHAMBERS_PROTO; + nLayers = URWellConstants.NLAYERS; + isProto = true; + Nconfig = config; + pitch = URWellConstants.PITCH; + } + case 2 -> { + nRegions = Math.min(URWellConstants.NMAXREGIONS, regions); + nSectors = URWellConstants.NSECTORS; + nChambers = URWellConstants.NCHAMBERS_ddvcs; + nLayers = URWellConstants.NLAYERS; + isProto = false; + Nconfig = config; + pitch = URWellConstants.PITCH_ddvcs; + + } + default -> { + } + } + + this.fillStripLists(); + this.fillPlaneLists(); + this.fillSurfaceLists(); + } + + public double getStereoAngle(int layer, int ichamber) { + + double stereoAngle = Math.toRadians(URWellConstants.STEREOANGLE); + + double[] dim = factory.getChamber_daughter_Dimensions(layer, ichamber); + + double yHalf = dim[0]; + double xHalfSmallBase = dim[1]; + double xHalfLargeBase = dim[2]; + + if (Nconfig == 2) { + + stereoAngle = Math.abs(Math.atan2(2 * yHalf, (-xHalfSmallBase + xHalfLargeBase))); + + } + + + if (layer % 2 != 0) { + stereoAngle = -stereoAngle; + } + + return stereoAngle; + } /** * Calculates the total number of strips in a sector - * + * + * @param layer * @return the strip number */ - public int getNStripSector() { + public int getNStripSector(int layer) { int nStrips = 0; for (int i = 0; i < nChambers; i++) { - nStrips += getNStripChamber(i); + // System.out.println("indice camera "+i+" numero strip "+getNStripChamber(i)); + + nStrips += getNStripChamber(layer, i); + } return nStrips; } /** * Calculates the number of strips in the given chamber - * + * + * @param layer (1 to N) * @param ichamber (0, 1, 2) * @return the strip number (1-N) */ - public int getNStripChamber(int ichamber) { - - double[] dim = factory.getChamber_daughter_Dimensions(ichamber); + public int getNStripChamber(int layer, int ichamber) { + + + double[] dim = factory.getChamber_daughter_Dimensions(layer, ichamber); - double yHalf = dim[0]; + double yHalf = dim[0]; double xHalfSmallBase = dim[1]; double xHalfLargeBase = dim[2]; + double stereoAngle = Math.abs(getStereoAngle(layer, ichamber)); + // C-------------D // // ------------- // // ----------- // @@ -123,59 +231,66 @@ public int getNStripChamber(int ichamber) { /** * * number of strip in AB** */ - int nAB = (int) (2 * xHalfSmallBase / (URWellConstants.PITCH - / Math.sin(Math.toRadians(URWellConstants.STEREOANGLE)))); + + int nAB = (int) Math.round(2 * xHalfSmallBase / (pitch / Math.sin(stereoAngle))); + + + double AC = Math.sqrt((Math.pow((xHalfSmallBase - xHalfLargeBase), 2) + Math.pow((2 * yHalf), 2))); double theta = Math.acos(2 * yHalf / AC); - int nAC = (int) (AC / (URWellConstants.PITCH - / Math.cos(theta - Math.toRadians(URWellConstants.STEREOANGLE)))); - int nStrips = nAB + nAC +1 ; + int nAC = (int) Math.round(AC / (pitch / Math.cos(theta - stereoAngle))); + + int nStrips = nAB + nAC; return nStrips; } /** * Provides the index of the chamber containing the strip with the given ID - * + * + * @param layer (1 to N) * @param strip (1 to N) * @return the chamber index (0, 1, 2) */ - public int getChamberIndex(int strip) { + public int getChamberIndex(int layer, int strip) { int nStripTotal = 0; - for(int i=0; i strip ID chamber (from 1 to getNStripChamber) int nStripTotal = 0; - if (chamberIndex > 0) { + for (int i = 0; i < chamberIndex; i++) { - nStripTotal += this.getNStripChamber(i); + + nStripTotal += this.getNStripChamber(layer, i); } - } - + + //Strip ID: from 1 to getNStripChamber int cStrip = strip - nStripTotal; - + return cStrip; } - + /** * Builds the given strip line in the CLAS12 frame + * * @param sector (1-6) * @param layer (1-2) * @param strip (1-N) @@ -183,67 +298,60 @@ private int getLocalStripId(int strip) { */ private Line3d createStrip(int sector, int layer, int strip) { - int chamberIndex = getChamberIndex(strip); - - int cStrip = this.getLocalStripId(strip); + int chamberIndex = getChamberIndex(layer, strip); + double stereoAngle = getStereoAngle(layer,chamberIndex); - + int cStrip = this.getLocalStripId(layer, strip); // CHAMBER reference frame // new numeration with stri ID_strip=0 crossing (0,0,0) of chamber - double[] dim = factory.getChamber_daughter_Dimensions(chamberIndex); - - double yHalf = dim[0]; + double[] dim = factory.getChamber_daughter_Dimensions(layer, chamberIndex); + + double yHalf = dim[0]; double xHalfSmallBase = dim[1]; double xHalfLargeBase = dim[2]; - - + // Y coordinate of the intersection point between the x=0 and the strip line crossing for B + double DY = -yHalf - Math.tan(Math.abs(stereoAngle)) * xHalfSmallBase; - double DY = -yHalf - Math.tan(Math.toRadians(URWellConstants.STEREOANGLE)) *xHalfSmallBase; - // ID of the strip - int nS = (int) (DY * Math.cos(Math.toRadians(URWellConstants.STEREOANGLE)) / URWellConstants.PITCH); + int nS = (int) (DY * Math.cos(Math.abs(stereoAngle)) / pitch); int nCStrip = nS + (cStrip - 1); - - //strip straight line chamber reference frame -> y = mx +c; - double stereoAngle = URWellConstants.STEREOANGLE; - if (layer % 2 != 0) { - stereoAngle = -URWellConstants.STEREOANGLE; - } - double m = Math.tan(Math.toRadians(stereoAngle)); - double c = nCStrip * URWellConstants.PITCH / Math.cos(Math.toRadians(stereoAngle)); - + + + + double m = Math.tan(stereoAngle); + double c = nCStrip * pitch / Math.cos(stereoAngle); + // Take 2 points in the strip straight line. They needs to define Line object double oX = -xHalfLargeBase; double oY = -xHalfLargeBase * m + c; double oZ = 0; Vector3d origin = new Vector3d(oX, oY, oZ); - + double eX = xHalfLargeBase; double eY = xHalfLargeBase * m + c; double eZ = 0; Vector3d end = new Vector3d(eX, eY, eZ); - + // Get Chamber Volume - Geant4Basic chamberVolume = factory.getChamberVolume(sector, chamberIndex+1, layer, isProto); - + Geant4Basic chamberVolume = factory.getChamberVolume(sector, chamberIndex + 1, layer); + // 2 point defined before wrt the GLOBAL frame Vector3d globalOrigin = chamberVolume.getGlobalTransform().transform(origin); - - Vector3d globalEnd = chamberVolume.getGlobalTransform().transform(end); + Vector3d globalEnd = chamberVolume.getGlobalTransform().transform(end); Straight line = new Line3d(globalOrigin, globalEnd); - + // CHECK intersections between line and volume chamberVolume.makeSensitive(); List Hits = chamberVolume.getIntersections(line); - + if (Hits.size() >= 1) { - - Vector3d TestOrigin = Hits.get(0).origin(); - Vector3d TestEnd = Hits.get(0).end(); + + Vector3d TestOrigin = Hits.get(0).origin(); + Vector3d TestEnd = Hits.get(0).end(); return new Line3d(Hits.get(0).origin(), Hits.get(0).end()); @@ -252,75 +360,80 @@ private Line3d createStrip(int sector, int layer, int strip) { } } - /** + /** * Provides the given strip line in the Chamber local frame - * @param region (1-2) + * + * * @param sector (1-6) * @param layer (1-4) * @param strip (1-N) * @return the 3D strip line as a Line3d */ - - private Line3d getChamberStrip(int region, int sector, int chamber, int layer, int strip) { + private Line3d getChamberStrip(int sector, int layer, int strip) { - Line3d globalStrip = createStrip(sector, layer, strip); - Geant4Basic chamberVolume = factory.getChamberVolume(sector, chamber, layer, isProto); + int chamber = getChamberIndex(layer, strip); + + Geant4Basic chamberVolume = factory.getChamberVolume(sector, chamber+1, layer); Vector3d origin = chamberVolume.getGlobalTransform().invert().transform(globalStrip.origin()); - Vector3d end = chamberVolume.getGlobalTransform().invert().transform(globalStrip.end()); + Vector3d end = chamberVolume.getGlobalTransform().invert().transform(globalStrip.end()); Line3d localStrip = new Line3d(origin, end); - return localStrip; } - - /** * Provides the given strip line in the sector local frame + * * @param sector (1-6) * @param layer (1-2) * @param strip (1-N) * @return the 3D strip line as a Line3d */ - private Line3d getLocalStrip(int region, int sector, int layer, int strip) { + private Line3d getLocalStrip(int sector, int layer, int strip) { + + int region = (layer - 1) / 2 + 1; Line3d globalStrip = createStrip(sector, layer, strip); Geant4Basic sVolume = factory.getSectorVolume(region, sector); Vector3d origin = sVolume.getGlobalTransform().invert().transform(globalStrip.origin()); - Vector3d end = sVolume.getGlobalTransform().invert().transform(globalStrip.end()); + Vector3d end = sVolume.getGlobalTransform().invert().transform(globalStrip.end()); Line3d localStrip = new Line3d(origin, end); return localStrip; } - - + private void fillStripLists() { - - for(int ir=0; irnull will be used to exit the integration + * early because some condition has been reached. + * @return the number of steps used--may be less than the space provided if + * the integration ended early as a result of an exit condition. + */ + public int uniformStepWithEnergyLoss(double yo[], + double to, + double tf, + final double y[][], + final double t[], + IDerivative deriv, + IStopper stopper) { + int nstep = t.length; // the number of steps to store + + // the dimensionality of the problem, e.g. six if (x, y, z, vx, vy, vz) + final int nDim = yo.length; + + // uniform step size + double h = (tf - to) / (nstep - 1); + + // put starting step in + t[0] = to; + for (int i = 0; i < nDim; i++) { + y[i][0] = yo[i]; + } + + IRkListener listener = new IRkListener() { + + int step = 1; + + @Override + public void nextStep(double tNext, double yNext[], double h) { + + if (step < t.length) { + t[step] = tNext; + for (int i = 0; i < nDim; i++) { + // store results for this step + y[i][step] = yNext[i]; + } + } + step++; + } + + }; + + return uniformStepWithEnergyLoss(yo, to, tf, h, deriv, stopper, listener); + + } + /** * Integrator that uses the standard RK4 advance with a uniform step size. @@ -155,6 +238,46 @@ public int uniformStep(double yo[], UniformAdvance advancer = new UniformAdvance(); return driver(yo, to, tf, h, deriv, stopper, listener, advancer); } + + /** + * Integrator that uses the standard RK4 advance with a uniform step size. + * (i.e., this does NOT use an adaptive step size.) + * + * This version uses an IRk4Listener to notify the listener that the next + * step has been advanced. + * + * A very typical case is a 2nd order ODE converted to a 1st order where the + * dependent variables are x, y, z, vx, vy, vz and the independent variable + * is time. + * + * @param yo + * initial values. Probably something like (xo, yo, zo, vxo, vyo, + * vzo). + * @param to + * the initial value of the independent variable, e.g., time. + * @param tf + * the maximum value of the independent variable. + * @param deriv + * the derivative computer (interface). This is where the problem + * specificity resides. + * @param stopper + * if not null will be used to exit the integration + * early because some condition has been reached. + * @param listener + * listens for each step + * @return the number of steps used. + */ + public int uniformStepWithEnergyLoss(double yo[], + double to, + double tf, + double h, + IDerivative deriv, + IStopper stopper, + IRkListener listener) { + + UniformAdvance advancer = new UniformAdvance(); + return driverWithEnergyLoss(yo, to, tf, h, deriv, stopper, listener, advancer); + } /** * Integrator that uses the RungeKutta advance with a Butcher Tableau and @@ -563,6 +686,7 @@ private int driver(double yo[], // typically [x, y, z, vx, vy, vz] and derivative double yt[] = new double[nDim]; double dydt[] = new double[nDim]; + double yttemp[] = new double[nDim]; double t = to; for (int i = 0; i < nDim; i++) { @@ -570,10 +694,15 @@ private int driver(double yo[], } for (int k = 1; k < nstep; k++) { + for (int i = 0; i < nDim; i++) { + yttemp[i] = yt[i]; + } + // use derivs at previous t deriv.derivative(t, yt, dydt); - - advancer.advance(t, yt, dydt, h, deriv, yt, null); + + advancer.advance(t, yt, dydt, h, deriv, yt, null); // yt is updated + t += h; // someone listening? @@ -592,6 +721,89 @@ private int driver(double yo[], return nstep; } + + + /** + * Driver that uses the RungeKutta advance with a uniform step size. (I.e., + * this does NOT use an adaptive step size.) + * + * This version uses an IRk4Listener to notify the listener that the next + * step has been advanced. + * + * A very typical case is a 2nd order ODE converted to a 1st order where the + * dependent variables are x, y, z, vx, vy, vz and the independent variable + * is time. + * + * @param yo + * initial values. Probably something like (xo, yo, zo, vxo, vyo, + * vzo). + * @param to + * the initial value of the independent variable, e.g., time. + * @param tf + * the maximum value of the independent variable. + * @param deriv + * the derivative computer (interface). This is where the problem + * specificity resides. + * @param stopper + * if not null will be used to exit the integration + * early because some condition has been reached. + * @return the number of steps used. + */ + private int driverWithEnergyLoss(double yo[], + double to, + double tf, + double h, + IDerivative deriv, + IStopper stopper, + IRkListener listener, + IAdvance advancer) { + int nstep = (int) (1 + (tf - to) / h); // the number of steps to store + + // the dimensionality of the problem. E.., 6 if (x, y, z, vx, vy, vz) + int nDim = yo.length; + + // yt is the current value of the state vector, + // typically [x, y, z, vx, vy, vz] and derivative + double yt[] = new double[nDim]; + double dydt[] = new double[nDim]; + double yttemp[] = new double[nDim]; + + double t = to; + for (int i = 0; i < nDim; i++) { + yt[i] = yo[i]; + } + + for (int k = 1; k < nstep; k++) { + for (int i = 0; i < nDim; i++) { + yttemp[i] = yt[i]; + } + + // use derivs at previous t + deriv.derivative(t, yt, dydt); + + advancer.advance(t, yt, dydt, h, deriv, yt, null); // yt is updated + + // Calculate energy loss, update momentum and alpha, and accumulate energy loss into totalEnergyLoss + deriv.energyLossUpdate(yttemp, h); + + t += h; + + // someone listening? + if (listener != null) { + listener.nextStep(t, yt, h); + } + + // premature termination? Skip if stopper is null. + if (stopper != null) { + stopper.setFinalT(t); + if (stopper.stopIntegration(t, yt)) { + return k + 1; // actual number of steps taken + } + } + } + + return nstep; + } @@ -1018,7 +1230,7 @@ private int driver(double yo[], hdata[1] += h; hdata[2] = Math.max(hdata[2], h); } - + for (int i = 0; i < nDim; i++) { yt[i] = yt2[i]; } @@ -1303,6 +1515,5 @@ public double getMaxStepSize() { public double getMinStepSize() { return _minStepSize; } - - + } diff --git a/common-tools/cnuphys/swimmer/src/main/java/cnuphys/swim/DefaultDerivative.java b/common-tools/cnuphys/swimmer/src/main/java/cnuphys/swim/DefaultDerivative.java index 165a32eddb..c09040d9b4 100644 --- a/common-tools/cnuphys/swimmer/src/main/java/cnuphys/swim/DefaultDerivative.java +++ b/common-tools/cnuphys/swimmer/src/main/java/cnuphys/swim/DefaultDerivative.java @@ -5,7 +5,26 @@ public class DefaultDerivative implements IDerivative { + // ddvcs parameters + //geometry + private static double[] THETA_SHIELD = {6.64, 41.3}; // min,max polar angle (deg) + private static double[] THETA_ECAL = {6.64, 30}; // min,max polar angle (deg) + private static double[] RHO_SHIELD = {60.0, 140}; // min,max distance from (0,0,0) at theta=25 deg + private static double[] RHO_ECAL = {60.0, 80.0}; // min,max distance from (0,0,0) at theta=25 deg + private static double TAN25 = Math.tan(Math.toRadians(25)); + private static double COS25 = Math.cos(Math.toRadians(25)); + // eloss + private static double EMASS = 0.511E-3; + private static double MUMASS = 105.66E-3; + private static double K = 0.000307075; // GeV mol-1 cm2 + // materials + private static double[] IeV = {823E-9, 600.7E-9}; // GeV + private static double[] DENSITY = {11.35, 8.3}; // g/cm3 + private static double[] ZOVERA = {0.39573, 0.41315}; + private FieldProbe _probe; + + private int _charge; private double _momentum; //GeV/c @@ -13,6 +32,9 @@ public class DefaultDerivative implements IDerivative { // e is the electron charge = (10^-9) in GeV/(T*m) // p is in GeV/c private double _alpha; + + private double _totalEnergyLoss = 0; + private double _energyLoss = 0; //for mag field result float b[] = new float[3]; @@ -29,6 +51,7 @@ public class DefaultDerivative implements IDerivative { */ public DefaultDerivative(int charge, double momentum, FieldProbe field) { _probe = field; + _charge = charge; _momentum = momentum; //units of this alpha are 1/(T*m) _alpha = 1.0e-9 * charge * Swimmer.C / _momentum; @@ -36,6 +59,7 @@ public DefaultDerivative(int charge, double momentum, FieldProbe field) { public void set(int charge, double momentum, FieldProbe field) { _probe = field; + _charge = charge; _momentum = momentum; //units of this alpha are 1/(T*m) _alpha = 1.0e-9 * charge * Swimmer.C / _momentum; @@ -79,5 +103,65 @@ public void derivative(double s, double[] u, double[] du) { du[1] = u[4]; du[2] = u[5]; } - + + /** + * Calculate energy loss, update momentum and alpha for each step, and accumulate energy loss into totalEnergyLoss + * + * @param u state vector ([x,y,z, tx = px/p, ty = py/p, tz = pz/p]). + * @param dx path length for the step; (= step size if step size is small enough) + */ + @Override + public void energyLossUpdate(double[] u, double dx){ + _energyLoss = getEloss(_momentum, u[0], u[1], u[2], dx); //todo: energy loss function + double energyOrigin = Math.sqrt(_momentum*_momentum + MUMASS*MUMASS); + double energyFinal = energyOrigin + _energyLoss; + _momentum = Math.sqrt(energyFinal*energyFinal - MUMASS*MUMASS); + _alpha = 1.0e-9 * _charge * Swimmer.C / _momentum; + _totalEnergyLoss += _energyLoss; + } + + + /** + * Get total energy loss + * @return total energy loss + */ + public double getTotalEnergyLoss(){ + return _totalEnergyLoss; + } + + /** + * Calculate energy loss + * @param p momentum in GeV + * @param x position in meters + * @param y position in meters + * @param z position in meters + * @param dx step pathlength in meters + * @return energy loss + */ + public double getEloss(double p, double x, double y, double z, double dx) { + + double dz = Math.sqrt(x*x+y*y)*TAN25; + double r = (z+dz)*COS25; + double r_cm = r * 100; // convert to cm + double theta = Math.toDegrees(Math.acos(z/Math.sqrt(x*x+y*y+z*z))); // degree + + if(r_cmRHO_SHIELD[1] || thetaTHETA_SHIELD[1]) + return 0; + + int imat = 0; + if(r_cm implements Serializable /** The source of the trajectory e.g. hbtracking */ private String _source = "???"; + + private double totalEnergyLoss = 0; + + /** + * Set total energy loss + * @param totalEnergyLoss total energy loss in GeV + */ + public void setTotalEnergyLoss(double totalEnergyLoss){ + this.totalEnergyLoss = totalEnergyLoss; + } + + /** + * Get total energy loss + * @return total energy loss + */ + public double getTotalEnergyLoss(){ + return totalEnergyLoss; + } /** * Create a swim trajectory with no initial content diff --git a/common-tools/cnuphys/swimmer/src/main/java/cnuphys/swim/Swimmer.java b/common-tools/cnuphys/swimmer/src/main/java/cnuphys/swim/Swimmer.java index c170bebcf8..a45afd3ea7 100644 --- a/common-tools/cnuphys/swimmer/src/main/java/cnuphys/swim/Swimmer.java +++ b/common-tools/cnuphys/swimmer/src/main/java/cnuphys/swim/Swimmer.java @@ -845,6 +845,94 @@ public SwimTrajectory swim(int charge, double xo, double yo, double zo, double m // Integrate DefaultDerivative deriv = new DefaultDerivative(charge, momentum, _probe); ntotal = (new RungeKutta()).uniformStep(uo, 0, maxPathLength, u, s, deriv, stopper); + + trajectory.setTotalEnergyLoss(deriv.getTotalEnergyLoss()); + + // now cycle through and get the save points + for (int i = 0; i < ntotal; i++) { + if (((i % cycle) == 0) || (i == (ntotal - 1))) { + double v[] = makeVector(u[0][i], u[1][i], u[2][i], u[3][i], u[4][i], u[5][i]); + trajectory.add(v); + } + } + + return trajectory; + } + + + /** + * Swims a Lund particle with a built it stopper for the maximum value of + * the radial coordinate. This is for the trajectory mode, where you want to + * cache steps along the path. Uses a fixed stepsize algorithm. + * + * @param charge + * the charge: -1 for electron, 1 for proton, etc + * @param xo + * the x vertex position in meters + * @param yo + * the y vertex position in meters + * @param zo + * the z vertex position in meters + * @param momentum + * initial momentum in GeV/c + * @param theta + * initial polar angle in degrees + * @param phi + * initial azimuthal angle in degrees + * @param stopper + * an optional object that can terminate the swimming based on + * some condition + * @param maxPathLength + * in meters. This determines the max number of steps based on + * the step size. If a stopper is used, the integration might + * terminate before all the steps are taken. A reasonable value + * for CLAS is 8. meters + * @param stepSize + * the uniform step size in meters. + * @param distanceBetweenSaves + * this distance is in meters. It should be bigger than stepSize. + * It is approximately the distance between "saves" where the + * point is saved in a trajectory for later drawing. + * @return the trajectory of the particle + */ + public SwimTrajectory swimWithEnergyLoss(int charge, double xo, double yo, double zo, double momentum, double theta, double phi, + IStopper stopper, double maxPathLength, double stepSize, double distanceBetweenSaves) { + + // if no magnetic field or no charge, then simple straight line tracks. + // the path will consist of just two points + if ((_probe == null) || (charge == 0)) { + GeneratedParticleRecord genPartRec = new GeneratedParticleRecord(charge, xo, yo, zo, momentum, theta, phi); + return straightLineTrajectory(genPartRec, maxPathLength); + } + + if (momentum < MINMOMENTUM) { + //System.err.println("Skipping low momentum swim (A)"); + return new SwimTrajectory(charge, xo, yo, zo, momentum, theta, phi); + } + + // cycle is the number of advances per save + int cycle = (int) (distanceBetweenSaves / stepSize); + cycle = Math.max(2, cycle); + + // max number of possible steps--may not use all of them + int ntotal = (int) (maxPathLength / stepSize); // number steps + int nsave = ntotal / cycle; // aprox number saves + + // the the initial six vector + double uo[] = initialState(xo, yo, zo, theta, phi); + + // storage for time and state + double s[] = new double[ntotal]; + double u[][] = new double[6][ntotal]; + + // create the trajectory container + SwimTrajectory trajectory = new SwimTrajectory(charge, xo, yo, zo, momentum, theta, phi, nsave); + + // Integrate + DefaultDerivative deriv = new DefaultDerivative(charge, momentum, _probe); + ntotal = (new RungeKutta()).uniformStepWithEnergyLoss(uo, 0, maxPathLength, u, s, deriv, stopper); + + trajectory.setTotalEnergyLoss(deriv.getTotalEnergyLoss()); // now cycle through and get the save points for (int i = 0; i < ntotal; i++) { diff --git a/common-tools/cnuphys/swimmer/src/main/java/cnuphys/swimZ/SwimZDerivative.java b/common-tools/cnuphys/swimmer/src/main/java/cnuphys/swimZ/SwimZDerivative.java index b0cc781def..6978589066 100644 --- a/common-tools/cnuphys/swimmer/src/main/java/cnuphys/swimZ/SwimZDerivative.java +++ b/common-tools/cnuphys/swimmer/src/main/java/cnuphys/swimZ/SwimZDerivative.java @@ -89,5 +89,8 @@ public void derivative(double z, double[] x, double[] dxdz) { dxdz[2] = _qv * Ax; dxdz[3] = _qv * Ay; } + + @Override + public void energyLossUpdate(double[] u, double dx){} } diff --git a/common-tools/swim-tools/src/main/java/org/jlab/clas/swimtools/Swim.java b/common-tools/swim-tools/src/main/java/org/jlab/clas/swimtools/Swim.java index b0e45126e9..68995a0c70 100644 --- a/common-tools/swim-tools/src/main/java/org/jlab/clas/swimtools/Swim.java +++ b/common-tools/swim-tools/src/main/java/org/jlab/clas/swimtools/Swim.java @@ -52,6 +52,8 @@ public class Swim { private ProbeCollection PC; + private static double MUMASS = 105.66E-3; + /** * Class for swimming to various surfaces. The input and output units are cm and GeV/c */ @@ -963,12 +965,18 @@ public double[] SwimToBeamLine(double xB, double yB) { return null; BeamLineSwimStopper stopper = new BeamLineSwimStopper(xB, yB); - SwimTrajectory st = PC.CF.swim(_charge, _x0, _y0, _z0, _pTot, _theta, _phi, stopper, _maxPathLength, stepSize, + SwimTrajectory st = PC.CF.swimWithEnergyLoss(_charge, _x0, _y0, _z0, _pTot, _theta, _phi, stopper, _maxPathLength, stepSize, 0.0005); if(st==null) return null; st.computeBDL(PC.CP); // st.computeBDL(compositeField); + + // Get total energy loss and update final momentum + double totalEnergyLoss = st.getTotalEnergyLoss(); + double energyOrigin = Math.sqrt(_pTot*_pTot + MUMASS*MUMASS); + double energyFinal = energyOrigin + totalEnergyLoss; + _pTot = Math.sqrt(energyFinal*energyFinal - MUMASS*MUMASS); double[] lastY = st.lastElement(); diff --git a/reconstruction/dc/src/main/java/org/jlab/rec/dc/Constants.java b/reconstruction/dc/src/main/java/org/jlab/rec/dc/Constants.java index 25473ca915..d66f313a5e 100644 --- a/reconstruction/dc/src/main/java/org/jlab/rec/dc/Constants.java +++ b/reconstruction/dc/src/main/java/org/jlab/rec/dc/Constants.java @@ -14,6 +14,8 @@ import org.jlab.detector.geom.RICH.RICHGeoFactory; import org.jlab.geom.base.ConstantProvider; import org.jlab.detector.calib.utils.ConstantsManager; +import org.jlab.detector.calib.utils.DatabaseConstantProvider; +import org.jlab.detector.geant4.v2.URWELL.URWellStripFactory; import org.jlab.geom.base.Detector; import org.jlab.rec.dc.trajectory.TrajectorySurfaces; import org.jlab.utils.groups.IndexedTable; @@ -122,7 +124,7 @@ public static Constants getInstance() { public DCGeant4Factory dcDetector = null; public FTOFGeant4Factory ftofDetector = null; public Detector ecalDetector = null; - public Detector fmtDetector = null; + public URWellStripFactory fmtDetector = null; public RICHGeoFactory richDetector = null; public TrajectorySurfaces trajSurfaces = null; @@ -546,7 +548,7 @@ private synchronized void LoadGeometry(String geoVariation, double[][] shifts) { ConstantProvider providerFTOF = GeometryFactory.getConstants(DetectorType.FTOF, 11, geoVariation); ftofDetector = new FTOFGeant4Factory(providerFTOF); ecalDetector = GeometryFactory.getDetector(DetectorType.ECAL, 11, geoVariation); - fmtDetector = GeometryFactory.getDetector(DetectorType.FMT, 11, geoVariation); + fmtDetector = new URWellStripFactory(new DatabaseConstantProvider(),2, 6); ConstantsManager managerRICH = new ConstantsManager(geoVariation);; richDetector = new RICHGeoFactory(0, managerRICH, 11, false); // create the surfaces diff --git a/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/TrackCandListFinder.java b/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/TrackCandListFinder.java index c3d248a1bf..f2521e58d8 100644 --- a/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/TrackCandListFinder.java +++ b/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/TrackCandListFinder.java @@ -561,7 +561,7 @@ public void setTrackPars(Track cand, cand.fit_Successful = false; return; } - + // recalc new vertex using plane stopper //int sector = cand.get(2).getSector(); //double[] Vt = null; @@ -592,6 +592,7 @@ public void setTrackPars(Track cand, cand.fit_Successful = true; cand.set_TrackingInfoString(trking); + /* dcSwim.SetSwimParameters(xOrFix, yOrFix, zOrFix, pxOrFix, pyOrFix, pzOrFix, cand.get_Q()); @@ -610,6 +611,7 @@ public void setTrackPars(Track cand, //set the pseudocross at extrapolated position cand.set_PreRegion1CrossPoint(new Point3D(xInner, yInner, zInner)); cand.set_PreRegion1CrossDir(new Point3D(uxInner, uyInner, uzInner)); + */ } private Integer getKey(Track trk) { @@ -755,7 +757,7 @@ public void matchHits(List stateVecAtPlanesList, Track trk, } List fhits = new ArrayList<>(); - + /* dcSwim.SetSwimParameters(trk.get_Vtx0().x(), trk.get_Vtx0().y(), trk.get_Vtx0().z(), trk.get_pAtOrig().x(), trk.get_pAtOrig().y(), trk.get_pAtOrig().z(), trk.get_Q()); @@ -764,6 +766,7 @@ public void matchHits(List stateVecAtPlanesList, Track trk, if (ToFirstMeas == null) { return; } + */ for (StateVec st : stateVecAtPlanesList) { if (st == null) { diff --git a/reconstruction/dc/src/main/java/org/jlab/rec/dc/trajectory/TrajectorySurfaces.java b/reconstruction/dc/src/main/java/org/jlab/rec/dc/trajectory/TrajectorySurfaces.java index 06241218dc..5d6b9ceb8e 100644 --- a/reconstruction/dc/src/main/java/org/jlab/rec/dc/trajectory/TrajectorySurfaces.java +++ b/reconstruction/dc/src/main/java/org/jlab/rec/dc/trajectory/TrajectorySurfaces.java @@ -17,6 +17,7 @@ import java.io.PrintWriter; import java.util.logging.Logger; +import org.jlab.detector.geant4.v2.URWELL.URWellStripFactory; import org.jlab.geom.detector.ec.ECLayer; import org.jlab.geom.detector.ec.ECSuperlayer; import org.jlab.geom.detector.fmt.FMTLayer; @@ -46,7 +47,7 @@ public void setDetectorPlanes(List> planes) { } public void loadSurface(double targetPosition, double targetLength, DCGeant4Factory dcDetector, - FTOFGeant4Factory ftofDetector, Detector ecalDetector, Detector fmtDetector, RICHGeoFactory richDetector) { + FTOFGeant4Factory ftofDetector, Detector ecalDetector, URWellStripFactory fmtDetector, RICHGeoFactory richDetector) { // creating Boundaries for MS Constants.getInstance().Z[0]= targetPosition; Constants.getInstance().Z[1]= dcDetector.getWireMidpoint(0, 0, 0, 0).z; @@ -76,9 +77,9 @@ public void loadSurface(double targetPosition, double targetLength, DCGeant4Fact this.detectorPlanes.get(isector).add(new Surface(DetectorType.TARGET, sector, DetectorLayer.TARGET_CENTER, new Plane3D(0, 0, targetPosition, 0, 0, 1))); // Add FMT layers - for (int ilayer=0; ilayer<6; ++ilayer) { - FMTLayer fmtLayer = (FMTLayer) fmtDetector.getSector(0).getSuperlayer(0).getLayer(ilayer); - this.detectorPlanes.get(isector).add(new Surface(DetectorType.FMT, sector, ilayer+1, fmtLayer.getTrajectorySurface(), 0)); + for (int ilayer=0; ilayer<12; ++ilayer) { +// FMTLayer fmtLayer = (FMTLayer) fmtDetector.getSector(0).getSuperlayer(0).getLayer(ilayer); + this.detectorPlanes.get(isector).add(new Surface(DetectorType.FMT, sector, ilayer+1, fmtDetector.getSurface(sector, ilayer+1), 0)); } // Add DC diff --git a/reconstruction/fmt/pom.xml b/reconstruction/fmt/pom.xml index d550549ccf..6dfa6d77a1 100644 --- a/reconstruction/fmt/pom.xml +++ b/reconstruction/fmt/pom.xml @@ -14,6 +14,11 @@ + + org.jlab.clas12.detector + clas12detector-urwell + 13.0.2-SNAPSHOT + org.jlab.jnp jnp-hipo diff --git a/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/Constants.java b/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/Constants.java index 9e38a007d5..b05a2a809b 100644 --- a/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/Constants.java +++ b/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/Constants.java @@ -13,10 +13,13 @@ */ public class Constants { - public static final int NLAYERS = 6; + public static final int NLAYERS = 12; + // rMax difference between cross's clusters energy + public static double CROSSDELTAE = 200; + // DC-tracks to FMT-clusters matching parameter - public static double CIRCLECONFUSION = 1.2; // cm + public static double CIRCLECONFUSION = 12; // cm // min path for final swimming to beamline to reject failed swimming public static double MIN_SWIM_PATH = 0.2; diff --git a/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/banks/RecoBankWriter.java b/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/banks/RecoBankWriter.java index 9a0ac60374..87aed80467 100644 --- a/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/banks/RecoBankWriter.java +++ b/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/banks/RecoBankWriter.java @@ -113,6 +113,7 @@ private static DataBank fillFMTTracksBank(DataEvent event,List candlist) return bank; } + /* private static DataBank fillFMTTrajectoryBank(DataEvent event,List candlist) { DataBank bank = event.createBank("FMT::Trajectory", candlist.size()*Constants.NLAYERS); @@ -147,27 +148,14 @@ private static DataBank fillFMTTrajectoryBank(DataEvent event,List candli } return bank; } + */ - public static void appendFMTBanks(DataEvent event, List fhits, List clusters, - List tracks) { + public static void appendFMTBanks(DataEvent event, List tracks) { if (event == null) return; - if (fhits != null) { - event.appendBanks(fillFMTHitsBank(event, fhits)); - } - - if (clusters != null && clusters.size()>0) { - event.appendBanks(fillFMTClustersBank(event, clusters)); - } - -// if (crosses != null && crosses.size() > 0) { -// event.appendBanks(this.fillFMTCrossesBank(event, crosses)); -// } - if (tracks != null && tracks.size() > 0) { event.appendBanks(fillFMTTracksBank(event, tracks)); - event.appendBanks(fillFMTTrajectoryBank(event, tracks)); } } diff --git a/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/Track.java b/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/Track.java index 4b6bb58d2a..56e662cd67 100644 --- a/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/Track.java +++ b/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/Track.java @@ -6,12 +6,15 @@ import java.util.Map; import java.util.Map.Entry; import org.jlab.clas.swimtools.Swim; +import org.jlab.detector.base.DetectorType; import org.jlab.geom.prim.Line3D; +import org.jlab.geom.prim.Point3D; import org.jlab.geom.prim.Vector3D; import org.jlab.io.base.DataBank; import org.jlab.io.base.DataEvent; import org.jlab.rec.fmt.Constants; -import org.jlab.rec.fmt.cluster.Cluster; + +import org.jlab.rec.urwell.reader.URWellCluster; /** * @@ -45,7 +48,7 @@ public class Track { private int _NDF; private final Trajectory[] _DCtrajs = new Trajectory[Constants.NLAYERS]; - private final List[] _clusters = new ArrayList[Constants.NLAYERS]; + private final List[] _clusters = new ArrayList[Constants.NLAYERS]; private final Trajectory[] _FMTtrajs = new Trajectory[Constants.NLAYERS]; public Track() { @@ -53,7 +56,7 @@ public Track() { public Track(int _index, int _sector, int _q, double _x, double _y, double _z, - double _px, double _py, double _pz, List clusters) { + double _px, double _py, double _pz, List clusters) { this._index = _index; this._sector = _sector; this._q = _q; @@ -63,10 +66,10 @@ public Track(int _index, int _sector, int _q, double _x, double _y, double _z, this._px = _px; this._py = _py; this._pz = _pz; - for(Cluster cluster : clusters) this.addCluster(cluster); + for(URWellCluster cluster : clusters) this.addCluster(cluster); } - + /** * @param layer @@ -84,20 +87,20 @@ public void setDCTraj(Trajectory trj) { this._DCtrajs[trj.getLayer()-1] = trj; } - public List getClusters() { - List clusters = new ArrayList<>(); + public List getClusters() { + List clusters = new ArrayList<>(); for(int i=0; i getClusters(int layer) { + public List getClusters(int layer) { if(layer<=0 || layer>Constants.NLAYERS) return null; else return _clusters[layer-1]; } - public Cluster getCluster(int layer) { + public URWellCluster getCluster(int layer) { if(layer<=0 || layer>Constants.NLAYERS) return null; else if(_clusters[layer-1]== null || _clusters[layer-1].size()==0) return null; else return _clusters[layer-1].get(0); @@ -115,12 +118,12 @@ public int getClusterLayer(int layer) { if(_clusters[layer-1]!=null) return _clusters[layer-1].size(); else return 0; } - - public final void addCluster(Cluster cluster) { - if(this._clusters[cluster.getLayer()-1]==null) - this._clusters[cluster.getLayer()-1] = new ArrayList<>(); - this._clusters[cluster.getLayer()-1].add(cluster); - } + + public final void addCluster(URWellCluster cluster) { + if(this._clusters[cluster.layer()-1]==null) + this._clusters[cluster.layer()-1] = new ArrayList<>(); + this._clusters[cluster.layer()-1].add(cluster); + } public void clearClusters(int layer) { this._clusters[layer-1].clear(); @@ -305,71 +308,71 @@ public double getPz() { public void setPz(double _pz) { this._pz = _pz; } - - // FIXME: THIS METHOD SHOULD BE GENERALIZED - public double getSeedQuality() { - double quality=99; - if(this.getClusterLayer(1)>0 && this.getClusterLayer(2)>0 && this.getClusterLayer(3)>0) { - Line3D seg1 = this.getClusters(1).get(0).getGlobalSegment(); - Line3D seg2 = this.getClusters(2).get(0).getGlobalSegment(); - Line3D seg3 = this.getClusters(3).get(0).getGlobalSegment(); - quality = seg1.distanceSegments(seg3).distanceSegments(seg2).length(); - } - return quality; - } - - // FIXME: THIS METHOD SHOULD BE GENERALIZED - public void filterClusters(int mode) { - - double dref = Double.POSITIVE_INFINITY; - - int[] ibest = new int[Constants.NLAYERS]; - for(int i=0; i0 && this.getClusterLayer(2)>0 && this.getClusterLayer(3)>0) { - for(int i1=0; i1=0) { - List cls = new ArrayList<>(); - cls.add(_clusters[i].get(ibest[i])); - _clusters[i] = cls; - } - } - } +// // FIXME: THIS METHOD SHOULD BE GENERALIZED +// public double getSeedQuality() { +// double quality=99; +// if(this.getClusterLayer(1)>0 && this.getClusterLayer(2)>0 && this.getClusterLayer(3)>0) { +// Line3D seg1 = this.getClusters(1).get(0).getGlobalSegment(); +// Line3D seg2 = this.getClusters(2).get(0).getGlobalSegment(); +// Line3D seg3 = this.getClusters(3).get(0).getGlobalSegment(); +// quality = seg1.distanceSegments(seg3).distanceSegments(seg2).length(); +// } +// return quality; +// } + + +// // FIXME: THIS METHOD SHOULD BE GENERALIZED +// public void filterClusters(int mode) { +// +// double dref = Double.POSITIVE_INFINITY; +// +// int[] ibest = new int[Constants.NLAYERS]; +// for(int i=0; i0 && this.getClusterLayer(2)>0 && this.getClusterLayer(3)>0) { +// for(int i1=0; i1=0) { +// List cls = new ArrayList<>(); +// cls.add(_clusters[i].get(ibest[i])); +// _clusters[i] = cls; +// } +// } +// } - public static List getDCTracks(DataEvent event, Swim swimmer) { + public List getDCTracks(DataEvent event, Swim swimmer) { Map trackmap = new LinkedHashMap(); @@ -386,48 +389,44 @@ public static List getDCTracks(DataEvent event, Swim swimmer) { trk.setIndex(i); trk.setSector(trackBank.getByte("sector", i)); trk.setQ(trackBank.getByte("q", i)); - trk.setX(trackBank.getFloat("Vtx0_x", i)); - trk.setY(trackBank.getFloat("Vtx0_y", i)); - trk.setZ(trackBank.getFloat("Vtx0_z", i)); - trk.setPx(trackBank.getFloat("p0_x", i)); - trk.setPy(trackBank.getFloat("p0_y", i)); - trk.setPz(trackBank.getFloat("p0_z", i)); + + double vx = trackBank.getFloat("Vtx0_x", i); + double vy = trackBank.getFloat("Vtx0_y", i); + double vz = trackBank.getFloat("Vtx0_z", i); + + double px = trackBank.getFloat("p0_x", i); + double py = trackBank.getFloat("p0_y", i); + double pz = trackBank.getFloat("p0_z", i); + + + Point3D vertexLocal = transGlobaltoLocal(trk.getSector(), vx, vy, vz); + + Point3D momLocal = transGlobaltoLocal(trk.getSector(), px, py, pz); + + trk.setX(vertexLocal.x()); + trk.setY(vertexLocal.y()); + trk.setZ(vertexLocal.z()); + trk.setPx(momLocal.x()); + trk.setPy(momLocal.y()); + trk.setPz(momLocal.z()); trk.setStatus(1); - trackmap.put(id,trk); - - for(int j=0; j tracks = new ArrayList<>(); for(Entry entry: trackmap.entrySet()) { @@ -436,39 +435,33 @@ public static List getDCTracks(DataEvent event, Swim swimmer) { return tracks; } - private static double[] getTrajectory(double x, double y, double z, double px, double py, double pz, - int q, int layer, Swim swim) { - Vector3D p = Constants.getLayer(layer).getPlane().point().toVector3D(); - Vector3D n = Constants.getLayer(layer).getPlane().normal(); - Vector3D v = new Vector3D(x, y, z); - double d = p.dot(n); - if(v.dot(n) clusters; + private List clusters; public KFitter(Track track, Swim swimmer, int c) { sv = new StateVecs(swimmer); @@ -46,7 +47,7 @@ public KFitter(Track track, Swim swimmer, int c) { track.getPx(), track.getPy(), track.getPz(), track.getQ(), c); } - public KFitter(List clusters, int sector, double xVtx, double yVtx, double zVtx, + public KFitter(List clusters, int sector, double xVtx, double yVtx, double zVtx, double pxVtx, double pyVtx, double pzVtx, int q, Swim swimmer, int c) { sv = new StateVecs(swimmer); this.track = new Track(0,sector, q, xVtx, yVtx, zVtx, pxVtx, pyVtx, pzVtx, clusters); @@ -54,7 +55,7 @@ public KFitter(List clusters, int sector, double xVtx, double yVtx, dou this.init(clusters, sector, xVtx, yVtx, zVtx, pxVtx, pyVtx, pzVtx, q, c); } - private void init(List clusters, int sector, double xVtx, double yVtx, double zVtx, + private void init(List clusters, int sector, double xVtx, double yVtx, double zVtx, double pxVtx, double pyVtx, double pzVtx, int q, int c) { // initialize measVecs mv.setMeasVecs(clusters); @@ -65,7 +66,7 @@ private void init(List clusters, int sector, double xVtx, double yVtx, for (int i = 0; i < mSize; i++) sv.Z[i] = mv.measurements.get(i).z; // initialize stateVecs - sv.init(sector, xVtx, yVtx, zVtx, pxVtx, pyVtx, pzVtx, q, sv.Z[0], this, c); + sv.init(sector, xVtx, yVtx, zVtx, pxVtx, pyVtx, pzVtx, q, sv.Z[0], this, c); } public void runFitter(int sector) { @@ -115,6 +116,7 @@ public void runFitter(int sector) { // save final trajectory points + /* if(this.finalStateVec!=null) { for (int k = 0; k < svzLength; ++k) { Trajectory trj = new Trajectory(mv.measurements.get(k).layer, @@ -128,51 +130,29 @@ public void runFitter(int sector) { track.setFMTtraj(trj); } } + */ -// for (int li = 1; li <= 6; ++li) { -// // Get the state vector closest in z to the FMT layer. -// // NOTE: A simple optimization would be to do this on another for loop with only the -// // state vectors, saving the ones closest to the FMT layers to use them later -// // instead of looping through them 6 times. -// int closestSVID = -1; -// double closestSVDistance = Double.POSITIVE_INFINITY; -// for (int si = 0; si < sv.trackTraj.size(); ++si) { -// double svDistance = Math.abs(sv.trackTraj.get(si).z - Geometry.getLayerZ(li-1)); -// if (svDistance < closestSVDistance) { -// closestSVID = si; -// closestSVDistance = svDistance; -// } -// } -// -// // Get the state vector's y position in the layer's local coordinates. -// Point3D Pos = new Point3D(sv.trackTraj.get(closestSVID).x, sv.trackTraj.get(closestSVID).y, 0); -// Point3D locPos = Geometry.globalToLocal(Pos,li); -// -// // Store the cluster's residual. -// for (Cluster cl : this.clusters) { -// if (cl.get_Layer() == li) { -// cl.set_CentroidResidual(cl.get_Centroid() - locPos.y()); -// } -// } -// } } - + public Matrix filterCovMat(double[] H, Matrix Ci, double V) { + double det = Matrix5x5.inverse(Ci, first_inverse, adj); - if (Math.abs(det) < 1.e-30) return null; + if (Math.abs(det) < 1.e-60) { + return null; + } addition.set( - H[0] * H[0] / V, H[0] * H[1] / V, H[0] * H[2] / V, H[0] * H[3] / V, H[0] * H[4] / V, - H[1] * H[0] / V, H[1] * H[1] / V, H[1] * H[2] / V, H[1] * H[3] / V, H[1] * H[4] / V, - H[2] * H[0] / V, H[2] * H[1] / V, H[2] * H[2] / V, H[2] * H[3] / V, H[2] * H[4] / V, - H[3] * H[0] / V, H[3] * H[1] / V, H[3] * H[2] / V, H[3] * H[3] / V, H[3] * H[4] / V, - H[4] * H[0] / V, H[4] * H[1] / V, H[4] * H[2] / V, H[4] * H[3] / V, H[4] * H[4] / V - ); + H[0] * H[0] / V, H[0] * H[1] / V, 0, 0, 0, + H[0] * H[1] / V, H[1] * H[1] / V, 0, 0, 0, + 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0); Matrix5x5.add(first_inverse, addition, result); double det2 = Matrix5x5.inverse(result, result_inv, adj); - - if (Math.abs(det2) < 1.e-30) return null; + if (Math.abs(det2) < 1.e-60) { + return null; + } return result_inv; } @@ -180,29 +160,41 @@ public Matrix filterCovMat(double[] H, Matrix Ci, double V) { private void filter(int k) { if (sv.trackTraj.get(k) != null && sv.trackCov.get(k).covMat != null && k < sv.Z.length ) { double[] K = new double[5]; - double V = mv.measurements.get(k).error; + double V = mv.measurements.get(k).error * mv.measurements.get(k).error; - double[] H = mv.H(sv.trackTraj.get(k),sv); + double[] H = mv.HURWell(sv.trackTraj.get(k),sv); Matrix CaInv = this.filterCovMat(H, sv.trackCov.get(k).covMat, V); if (CaInv != null) sv.trackCov.get(k).covMat = CaInv; else return; + Matrix cMat = new Matrix(); + if (CaInv != null) { + Matrix5x5.copy(CaInv, cMat); + } else { + return; + } + // Calculate the gain matrix. for (int j = 0; j < 5; j++) { - K[j] = 0; - for (int i = 0; i < 5; i++) K[j] += H[i] * sv.trackCov.get(k).covMat.get(j, i) / V ; + // the gain matrix + K[j] = (H[0] * cMat.get(j, 0) + + H[1] * cMat.get(j, 1)) / V; } - // Update Chi^2 and filtered state vector. - double h = mv.h(sv.trackTraj.get(k)); - double m = mv.measurements.get(k).centroid; - this.chi2 += ((m-h)*(m-h)/mv.measurements.get(k).error/mv.measurements.get(k).error); - - double x_filt = sv.trackTraj.get(k).x + K[0] * (m-h); - double y_filt = sv.trackTraj.get(k).y + K[1] * (m-h); - double tx_filt = sv.trackTraj.get(k).tx + 0*K[2] * (m-h); - double ty_filt = sv.trackTraj.get(k).ty + 0*K[3] * (m-h); - double Q_filt = sv.trackTraj.get(k).Q + 0*K[4] * (m-h); + // Update Chi^2 and filtered state vector. + double res = mv.dhURWell(sv.trackTraj.get(k)); + double filt[] = new double[5]; + for(int j = 0; j < 5; j ++){ + filt[j] += K[j]*res; + } + + this.chi2 += (res*res/mv.measurements.get(k).error/mv.measurements.get(k).error); + + double x_filt = sv.trackTraj.get(k).x + filt[0]; + double y_filt = sv.trackTraj.get(k).y + filt[1]; + double tx_filt = sv.trackTraj.get(k).tx + filt[2]; + double ty_filt = sv.trackTraj.get(k).ty + filt[3]; + double Q_filt = sv.trackTraj.get(k).Q + filt[4]; if (filterOn) { sv.trackTraj.get(k).x = x_filt; diff --git a/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/fit/MeasVecs.java b/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/fit/MeasVecs.java index dba901877d..381b659b9b 100644 --- a/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/fit/MeasVecs.java +++ b/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/fit/MeasVecs.java @@ -4,9 +4,13 @@ import java.util.List; import org.jlab.geom.prim.Vector3D; import org.jlab.rec.fmt.Constants; -import org.jlab.rec.fmt.cluster.Cluster; +import org.jlab.rec.urwell.reader.URWellCluster; import org.jlab.rec.fmt.track.fit.StateVecs.StateVec; +import org.jlab.geom.prim.Point3D; +import org.jlab.geom.prim.Line3D; +import org.jlab.geom.prim.Plane3D; + /** * @author ziegler */ @@ -17,7 +21,8 @@ public class MeasVecs { public class MeasVec implements Comparable { public double z = Double.NaN; - public double centroid; + public Point3D lineEndPoint1 = null; + public Point3D lineEndPoint2 = null; public double seed; public double error; public int layer; @@ -33,47 +38,84 @@ public int compareTo(MeasVec arg) { } } - public void setMeasVecs(List clusters) { + public void setMeasVecs(List clusters) { measurements = new ArrayList<>(); for (int i = 0; i < clusters.size(); i++) { - int l = clusters.get(i).getLayer()-1; - double cent = clusters.get(i).getCentroid(); - double error = clusters.get(i).getCentroidError(); - double z = clusters.get(i).getGlobalSegment().origin().z(); - int seed = clusters.get(i).getSeedStrip(); - MeasVec meas = this.setMeasVec(l, cent, error, z, seed); + int l = clusters.get(i).layer()-1; + Point3D lineEndPoint1 = clusters.get(i).getLineLocal().origin(); + Point3D lineEndPoint2 = clusters.get(i).getLineLocal().end(); + double error = 0.0144; // = pitch/sqrt(12), where pitch is 500 um + double z = lineEndPoint1.z(); + int seed = clusters.get(i).strip(); + MeasVec meas = this.setMeasVec(l, lineEndPoint1, lineEndPoint2, error, z, seed); measurements.add(meas); } } - public MeasVec setMeasVec(int l, double cent, double error, double z, int seed) { + public MeasVec setMeasVec(int l, Point3D lineEndPoint1, Point3D lineEndPoint2, double error, double z, int seed) { MeasVec meas = new MeasVec(); meas.layer = l+1; - meas.centroid = cent; meas.error = error; meas.z = z; meas.seed = seed; + meas.lineEndPoint1 = lineEndPoint1; + meas.lineEndPoint2 = lineEndPoint2; + return meas; } - public double h(StateVec stateVec) { - if (stateVec == null) return 0; - if (this.measurements.get(stateVec.k) == null) return 0; - - int layer = this.measurements.get(stateVec.k).layer; - - return Constants.toLocal(layer, stateVec.x, stateVec.y, stateVec.z).y(); + public double dhURWell(StateVec stateVec) { + double value = Double.NaN; + if (stateVec == null|| this.measurements.get(stateVec.k) == null) { + return value; + } + + Line3D l = new Line3D(this.measurements.get(stateVec.k).lineEndPoint1, + this.measurements.get(stateVec.k).lineEndPoint2); + Line3D WL = new Line3D(); + WL.copy(l); + Point3D svP = new Point3D(stateVec.x, stateVec.y, stateVec.z); + WL.copy(WL.distance(svP)); + Plane3D plane = new Plane3D(0, 0, this.measurements.get(stateVec.k).z, 0, 0, 1); // plan perpenticular to z axis in TSC + double sideStrip = -Math.signum(l.direction().cross(WL.direction()). + dot(plane.normal())); + value = WL.length()*sideStrip; + + return value; } + + public double[] HURWell(StateVec stateVec, StateVecs sv) { + double[] hMatrix = new double[5]; + double Err = 0.01; + double[][] Result = new double[2][2]; + for (int i = 0; i < 2; i++) { + StateVec svc = sv.new StateVec(stateVec.k); + svc.x = stateVec.x; + svc.y = stateVec.y; + svc.z = stateVec.z; + svc.x = stateVec.x + (double) Math.pow(-1, i) * Err; + Result[i][0] = dhURWell(svc); + } + for (int i = 0; i < 2; i++) { + StateVec svc = sv.new StateVec(stateVec.k); + svc.x = stateVec.x; + svc.y = stateVec.y; + svc.z = stateVec.z; + svc.y = stateVec.y + (double) Math.pow(-1, i) * Err; + Result[i][1] = dhURWell(svc); + } - public double[] H(StateVec stateVec, StateVecs sv) { - int layer = this.measurements.get(stateVec.k).layer; - Vector3D derivatives = Constants.getDerivatives(layer, stateVec.x, stateVec.y, stateVec.z); - double[] H = new double[]{derivatives.x(), derivatives.y(), 0, 0, 0}; - return H; - } + hMatrix[0] = -(Result[0][0] - Result[1][0]) / (2. * Err); // Add negative sign since dh = meas - h; here use dh to replace h since meas is cancelled when derivative + hMatrix[1] = -(Result[0][1] - Result[1][1]) / (2. * Err); // Add negative sign since dh = meas - h; here use dh to replace h since meas is cancelled when derivative + hMatrix[2] = 0; + hMatrix[3] = 0; + hMatrix[4] = 0; + + return hMatrix; + } private StateVec reset(StateVec SVplus, StateVec stateVec, StateVecs sv) { SVplus = sv.new StateVec(stateVec.k); diff --git a/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/fit/RungeKutta.java b/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/fit/RungeKutta.java index 809acc4fde..21162bd99d 100644 --- a/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/fit/RungeKutta.java +++ b/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/fit/RungeKutta.java @@ -34,7 +34,7 @@ public RungeKutta() { public void SwimToZ(int sector, StateVecs.StateVec fVec, Swim dcSwim, double z0, float[] bf){ double stepSize = 0.5; - dcSwim.BfieldLab(fVec.x, fVec.y, fVec.z, bf); + dcSwim.Bfield(sector, fVec.x, fVec.y, fVec.z, bf); fVec.B = Math.sqrt(bf[0]*bf[0]+bf[1]*bf[1]+bf[2]*bf[2]); double s = fVec.B; @@ -67,25 +67,25 @@ public void SwimToZ(int sector, StateVecs.StateVec fVec, Swim dcSwim, double z0, void RK4transport(int sector, double q, double x0, double y0, double z0, double tx0, double ty0, double h, Swim swimmer, double dPath, StateVecs.StateVec fVec) { // lab system = 1, TSC =0 - swimmer.BfieldLab(x0, y0, z0, _b); + swimmer.Bfield(sector, x0, y0, z0, _b); double x1 = tx0; double y1 = ty0; double tx1=q*v*Ax(tx0, ty0, _b[0], _b[1], _b[2]); double ty1=q*v*Ay(tx0, ty0, _b[0], _b[1], _b[2]); - swimmer.BfieldLab(x0+0.5*h*x1, y0+0.5*h*y1, z0+0.5*h, _b); + swimmer.Bfield(sector, x0+0.5*h*x1, y0+0.5*h*y1, z0+0.5*h, _b); double x2 = tx0+0.5*h*tx1; double y2 = ty0+0.5*h*ty1; double tx2=q*v*Ax((tx0+0.5*h*tx1), (ty0+0.5*h*ty1), _b[0], _b[1], _b[2]); double ty2=q*v*Ay((tx0+0.5*h*tx1), (ty0+0.5*h*ty1), _b[0], _b[1], _b[2]); - swimmer.BfieldLab(x0+0.5*h*x2, y0+0.5*h*y2, z0+0.5*h, _b); + swimmer.Bfield(sector, x0+0.5*h*x2, y0+0.5*h*y2, z0+0.5*h, _b); double x3 = tx0+0.5*h*tx2; double y3 = ty0+0.5*h*ty2; double tx3=q*v*Ax((tx0+0.5*h*tx2), (ty0+0.5*h*ty2), _b[0], _b[1], _b[2]); double ty3=q*v*Ay((tx0+0.5*h*tx2), (ty0+0.5*h*ty2), _b[0], _b[1], _b[2]); - swimmer.BfieldLab(x0+h*x3, y0+h*y3, z0+h, _b); + swimmer.Bfield(sector, x0+h*x3, y0+h*y3, z0+h, _b); double x4 = tx0+h*tx3; double y4 = ty0+h*ty3; double tx4=q*v*Ax((tx0+h*tx3), (ty0+h*ty3), _b[0], _b[1], _b[2]); @@ -125,7 +125,7 @@ void RK4transport(int sector, double q, double x0, double y0, double z0, double double delty_delq0_0 =0; //System.out.println("RK0 "+x0+","+y0+","+z0+";"+tx0+","+ty0+","+" z0 "+z0+" h "+h); //State - swimmer.BfieldLab(x0, y0, z0, _b); + swimmer.Bfield(sector, x0, y0, z0, _b); double x1 = tx0; double y1 = ty0; double tx1=q*v*Ax(tx0, ty0, _b[0], _b[1], _b[2]); @@ -157,7 +157,7 @@ void RK4transport(int sector, double q, double x0, double y0, double z0, double + delAy_delty(tx0,ty0,_b[0],_b[1],_b[2])*delty_delq0_0); - swimmer.BfieldLab(x0+0.5*h*x1, y0+0.5*h*y1, z0+0.5*h, _b); + swimmer.Bfield(sector, x0+0.5*h*x1, y0+0.5*h*y1, z0+0.5*h, _b); double x2 = tx0+0.5*h*tx1; double y2 = ty0+0.5*h*ty1; double tx2=q*v*Ax((tx0+0.5*h*tx1), (ty0+0.5*h*ty1), _b[0], _b[1], _b[2]); @@ -186,7 +186,7 @@ void RK4transport(int sector, double q, double x0, double y0, double z0, double double delty_delq0_2 = this.delty_delq0_next(q,v,tx0+0.5*h*tx1,ty0+0.5*h*ty1,_b[0],_b[1],_b[2], deltx_delq0_0+0.5*h*deltx_delq0_1,delty_delq0_0+0.5*h*delty_delq0_1); - swimmer.BfieldLab(x0+0.5*h*x2, y0+0.5*h*y2, z0+0.5*h, _b); + swimmer.Bfield(sector, x0+0.5*h*x2, y0+0.5*h*y2, z0+0.5*h, _b); double x3 = tx0+0.5*h*tx2; double y3 = ty0+0.5*h*ty2; double tx3=q*v*Ax((tx0+0.5*h*tx2), (ty0+0.5*h*ty2), _b[0], _b[1], _b[2]); @@ -215,7 +215,7 @@ void RK4transport(int sector, double q, double x0, double y0, double z0, double double delty_delq0_3 = this.delty_delq0_next(q,v,tx0+0.5*h*tx2,ty0+0.5*h*ty2,_b[0],_b[1],_b[2], deltx_delq0_0+0.5*h*deltx_delq0_2,delty_delq0_0+0.5*h*delty_delq0_2); - swimmer.BfieldLab(x0+h*x3, y0+h*y3, z0+h, _b); + swimmer.Bfield(sector, x0+h*x3, y0+h*y3, z0+h, _b); double x4 = tx0+h*tx3; double y4 = ty0+h*ty3; double tx4=q*v*Ax((tx0+h*tx3), (ty0+h*ty3), _b[0], _b[1], _b[2]); diff --git a/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/fit/StateVecs.java b/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/fit/StateVecs.java index 620faf7214..dad52c309e 100644 --- a/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/fit/StateVecs.java +++ b/reconstruction/fmt/src/main/java/org/jlab/rec/fmt/track/fit/StateVecs.java @@ -12,6 +12,8 @@ */ public class StateVecs { private final double Bmax = 2.366498; // averaged + + private final double MuonMass = 0.105658; // GeV/c^2 final double speedLight = 0.002997924580; public double[] Z; @@ -29,7 +31,7 @@ public class StateVecs { private final float[] lbf = new float[3]; private final Swim dcSwim; private final RungeKutta rk; - + /** * State vector representing the track in the sector coordinate system at the measurement layer * @param swimmer @@ -90,7 +92,8 @@ public Matrix transport(int sector, int i, double Zf, StateVec iVec, CovMat covM double X0 = this.getX0(z); double t_ov_X0 = Math.abs(s) / X0;//path length in radiation length units = t/X0 [true path length/ X0] ; Ar radiation length = 14 cm - double beta = this.beta; + double energy = Math.sqrt(p*p + MuonMass * MuonMass); + double beta = p/energy; if(beta>1.0 || beta<=0) beta =1.0; @@ -361,6 +364,7 @@ public void init(int sector, double xVtx, double yVtx, double zVtx, rinitSV.tx = pxVtx/pzVtx; rinitSV.ty = pyVtx/pzVtx; rinitSV.Q = (double)q / p; + /* double[] FTF = new double[25]; double[] F = this.F(sector, z0, rinitSV); for(int i = 0; i<5; i++) { @@ -369,6 +373,16 @@ public void init(int sector, double xVtx, double yVtx, double zVtx, Matrix initCMatrix = new Matrix(); initCMatrix.set(FTF); initCM.covMat = initCMatrix; + */ + Matrix initCMatrix = new Matrix(); + initCMatrix.set(8 * 8, 0, 0, 0, 0, + 0, 8 * 8, 0, 0, 0, + 0, 0, 0.1 * 0.1, 0, 0, + 0, 0, 0, 0.1 * 0.1, 0, + 0, 0, 0, 0, 0.03 * 0.03 + ); + initCM.covMat = initCMatrix; + this.trackCov.put(0, initCM); } private StateVec reset(StateVec SVplus, StateVec stateVec) { diff --git a/reconstruction/fmt/src/main/java/org/jlab/service/fmt/FMTEngine.java b/reconstruction/fmt/src/main/java/org/jlab/service/fmt/FMTEngine.java index de6e67e4d5..0c41257436 100644 --- a/reconstruction/fmt/src/main/java/org/jlab/service/fmt/FMTEngine.java +++ b/reconstruction/fmt/src/main/java/org/jlab/service/fmt/FMTEngine.java @@ -2,23 +2,30 @@ import java.util.Arrays; import java.util.List; +import java.util.ArrayList; +import java.util.Comparator; import org.jlab.clas.reco.ReconstructionEngine; import org.jlab.clas.swimtools.Swim; import org.jlab.detector.base.DetectorType; import org.jlab.detector.base.GeometryFactory; +import org.jlab.geom.prim.Line3D; +import org.jlab.geom.prim.Point3D; +import org.jlab.geom.prim.Vector3D; import org.jlab.utils.groups.IndexedTable; import org.jlab.io.base.*; import org.jlab.rec.fmt.Constants; import org.jlab.rec.fmt.banks.RecoBankWriter; -import org.jlab.rec.fmt.cluster.Cluster; -import org.jlab.rec.fmt.hit.Hit; import org.jlab.rec.fmt.track.Track; import org.jlab.rec.fmt.track.Trajectory; import org.jlab.rec.fmt.track.fit.KFitter; import org.jlab.rec.fmt.track.fit.StateVecs.StateVec; +import org.jlab.rec.urwell.reader.URWellCluster; +import org.jlab.rec.urwell.reader.URWellCross; +import org.jlab.rec.urwell.reader.URWellReader; + /** * Service to return reconstructed track candidates - the output is in hipo format * @@ -46,7 +53,6 @@ public boolean init() { System.out.println("["+this.getName()+"] run with FMT default geometry"); } - String[] tables = new String[]{ "/geometry/beam/position", "/calibration/mvt/fmt_time", @@ -57,12 +63,9 @@ public boolean init() { // Load the geometry int run = 10; - Constants.setDetector(GeometryFactory.getDetector(DetectorType.FMT,run, variation)); - + Constants.setDetector(GeometryFactory.getDetector(DetectorType.FMT,run, variation)); + // Register output banks - super.registerOutputBank("FMT::Hits"); - super.registerOutputBank("FMT::Clusters"); - super.registerOutputBank("FMT::Crosses"); super.registerOutputBank("FMT::Tracks"); super.registerOutputBank("FMT::Trajectory"); @@ -78,7 +81,7 @@ public boolean processDataEvent(DataEvent event) { DataBank runConfig = event.getBank("RUN::config"); if (runConfig == null || runConfig.rows() == 0) return true; int run = runConfig.getInt("run", 0); - + // Set swimmer. Swim swimmer = new Swim(); @@ -86,107 +89,105 @@ public boolean processDataEvent(DataEvent event) { // uncomment this code. IndexedTable beamOffset = this.getConstantsManager().getConstants(run, "/geometry/beam/position"); double xB = beamOffset.getDoubleValue("x_offset", 0,0,0); - double yB = beamOffset.getDoubleValue("y_offset", 0,0,0); - - // get status table - IndexedTable status = this.getConstantsManager().getConstants(run, "/calibration/mvt/fmt_status"); - - // get time table - IndexedTable timecuts = this.getConstantsManager().getConstants(run, "/calibration/mvt/fmt_time"); - - // === HITS ================================================================================ - List hits = Hit.fetchHits(event, timecuts, status); - if (hits.isEmpty()) return true; + double yB = beamOffset.getDoubleValue("y_offset", 0,0,0); + // === CLUSTERS ============================================================================ - List clusters = Cluster.findClusters(hits); - if(debug) for (int i = 0; i < clusters.size(); i++) System.out.println(clusters.get(i).toString()); - + URWellReader uRWellReader = new URWellReader(event); + List clusters = uRWellReader.getUrwellClusters(); + List crosses = uRWellReader.getUrwellCrosses(); + //System.out.println(clusters.size()); + //for (int i = 0; i < clusters.size(); i++) System.out.println(clusters.get(i).toString()); + // === DC TRACKS =========================================================================== - List tracks = Track.getDCTracks(event, swimmer); + Track trk = new Track(); + List tracks = trk.getDCTracks(event, swimmer); if(tracks.isEmpty()) return true; - + // === SEEDS ============================================================================= + Point3D target = new Point3D(0,0,0); for(int i=0; i[] trackCrosses = new ArrayList[Constants.NLAYERS/2]; + for(int j=0; j(); + for(int j=0; jother.getSeedQuality()) { - other.clearClusters(layer); - cluster.setTrackIndex(track.getIndex()); - cluster.setDoca(doca); + System.out.println(); + for(int j=0; j> segments = new ArrayList<>(); + List end = new ArrayList<>(); + for(int r=6; r>3; r--) + end.addAll(trackCrosses[r-1]); + for(URWellCross ce : end) { + Line3D ve = new Line3D(ce.position(),target); + for(int ro=1; ro segment = new ArrayList(); + segment.add(ce); + for(int r=co.region()+1; r=4) + segments.add(segment); } } } } + System.out.print(segments.size()); + if(!segments.isEmpty()) { + segments.sort(Comparator.comparingInt(List::size).reversed()); + for(URWellCross cross : segments.get(0)) { + track.addCluster(cross.getCluster1()); + track.addCluster(cross.getCluster2()); + } + } } - + + List filtedTracks = new ArrayList(); + for(Track track : tracks){ + if (track.getClusters().size() > 4) filtedTracks.add(track); + } + + // === TRACKS ============================================================================== KFitter kf = null; // Iterate on list to run the fit. - for(int i=0; i crs = crossesMap.get(id); -// for(Cross c : crs) { -// // Filter clusters to use only the best cluster (minimum Tmin) per FMT layer. -// int lyr = c.getCluster1().getLayer(); -//// System.out.println(c.getCluster1().getDoca()); -// if (cls.get(lyr) == null || c.getCluster1().getDoca()< cls.get(lyr).getDoca()) -// cls.put(lyr, c.getCluster1()); -// } - // Set status and stop if there are not at least two measurements to fit against. - List trackClusters = track.getClusters(); - - if (trackClusters.size()<2) continue; + List trackClusters = track.getClusters(); - kf = new KFitter(track, swimmer, 0); + kf = new KFitter(track, swimmer, 0); + kf.runFitter(track.getSector()); + // Do one last KF pass with filtering off to get the final Chi^2. kf.totNumIter = 1; kf.filterOn = false; @@ -197,11 +198,15 @@ public boolean processDataEvent(DataEvent event) { // swim to beamline to get vertex parameters int charge = (int)Math.signum(sv.Q); - swimmer.SetSwimParameters(sv.x,sv.y,sv.z,-sv.getPx(),-sv.getPy(),-sv.getPz(),-charge); + + Point3D posGlobal = track.transLocaltoGlobal(track.getSector(), sv.x, sv.y, sv.z); + Point3D momGlobal = track.transLocaltoGlobal(track.getSector(), sv.getPx()/sv.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 || Vt[6]thresholds.getDoubleValue("thresholdCluster",1,1,hit.get_COMPONENT())) { + if(hit.get_Edep()>thresholds.getDoubleValue("thresholdCluster",1,1,FTCALHit.REFCOMPONENT)) { for(int j = 0; j< this.size(); j++) { double tDiff = Math.abs(hit.get_Time() - this.get(j).get_Time()); double xDiff = Math.abs(hit.get_IDX() - this.get(j).get_IDX()); diff --git a/reconstruction/ft/src/main/java/org/jlab/rec/ft/cal/FTCALHit.java b/reconstruction/ft/src/main/java/org/jlab/rec/ft/cal/FTCALHit.java index 507bdf089c..de425fe10c 100644 --- a/reconstruction/ft/src/main/java/org/jlab/rec/ft/cal/FTCALHit.java +++ b/reconstruction/ft/src/main/java/org/jlab/rec/ft/cal/FTCALHit.java @@ -6,25 +6,26 @@ public class FTCALHit implements Comparable{ // class implements Comparable interface to allow for sorting a collection of hits by Edep values - + public static final int REFCOMPONENT =245; + // constructor public FTCALHit(int i, int ICOMPONENT, int ADC, int TDC, IndexedTable charge2Energy, IndexedTable timeOffsets, IndexedTable timeWalk, IndexedTable cluster) { this._COMPONENT = ICOMPONENT; - this._IDY = ((int) ICOMPONENT/22) + 1; - this._IDX = ICOMPONENT + 1 - (this._IDY-1)*22; + this._IDY = ((int) ICOMPONENT/44) + 1; + this._IDX = ICOMPONENT + 1 - (this._IDY-1)*44; this._ADC = ADC; this._TDC = TDC; - this._Charge = ((double) this._ADC)*charge2Energy.getDoubleValue("fadc_to_charge", 1,1,ICOMPONENT); - this.set_Edep(this._Charge*charge2Energy.getDoubleValue("mips_energy", 1,1,ICOMPONENT) - /charge2Energy.getDoubleValue("mips_charge", 1,1,ICOMPONENT)/1000.); + this._Charge = ((double) this._ADC)*charge2Energy.getDoubleValue("fadc_to_charge", 1,1,REFCOMPONENT); + this.set_Edep(this._Charge*charge2Energy.getDoubleValue("mips_energy", 1,1,REFCOMPONENT) + /charge2Energy.getDoubleValue("mips_charge", 1,1,REFCOMPONENT)/1000.); double twCorr=0; if(this._Charge>0) { - twCorr = timeWalk.getDoubleValue("amplitude", 1,1,ICOMPONENT)*Math.exp(-this._Charge*timeWalk.getDoubleValue("lambda", 1,1,ICOMPONENT)); + twCorr = timeWalk.getDoubleValue("amplitude", 1,1,REFCOMPONENT)*Math.exp(-this._Charge*timeWalk.getDoubleValue("lambda", 1,1,REFCOMPONENT)); } this.set_Time(((double) this._TDC)/FTCALConstantsLoader.TIMECONVFAC -(FTCALConstantsLoader.CRYS_LENGTH-cluster.getDoubleValue("depth_z", 1,1,0))/FTCALConstantsLoader.VEFF - -timeOffsets.getDoubleValue("time_offset", 1,1,ICOMPONENT)-twCorr); + -timeOffsets.getDoubleValue("time_offset", 1,1,REFCOMPONENT)-twCorr); // if(this.get_Edep()>0.1) System.out.println(ICOMPONENT + " " + this._TDC + " " + // FTCALConstantsLoader.TIMECONVFAC + " " + FTCALConstantsLoader.time_offset[0][0][ICOMPONENT-1] + " " + // this.get_Time()); @@ -37,21 +38,21 @@ public FTCALHit(int i, int ICOMPONENT, int ADC, int TDC, IndexedTable charge2Ene public FTCALHit(int i, int ICOMPONENT, int ADC, float time, IndexedTable charge2Energy, IndexedTable timeOffsets, IndexedTable timeWalk, IndexedTable cluster) { this._COMPONENT = ICOMPONENT; - this._IDY = ((int) ICOMPONENT/22) + 1; - this._IDX = ICOMPONENT + 1 - (this._IDY-1)*22; + this._IDY = ((int) ICOMPONENT/44) + 1; + this._IDX = ICOMPONENT + 1 - (this._IDY-1)*44; this._ADC = ADC; - this._Charge = ((double) this._ADC)*charge2Energy.getDoubleValue("fadc_to_charge", 1,1,ICOMPONENT); - this.set_Edep(this._Charge*charge2Energy.getDoubleValue("mips_energy", 1,1,ICOMPONENT) - /charge2Energy.getDoubleValue("mips_charge", 1,1,ICOMPONENT)/1000.); + this._Charge = ((double) this._ADC)*charge2Energy.getDoubleValue("fadc_to_charge", 1,1,REFCOMPONENT); + this.set_Edep(this._Charge*charge2Energy.getDoubleValue("mips_energy", 1,1,REFCOMPONENT) + /charge2Energy.getDoubleValue("mips_charge", 1,1,REFCOMPONENT)/1000.); double twCorr=0; if(this._Charge>0) { - twCorr = timeWalk.getDoubleValue("amplitude", 1,1,ICOMPONENT)*Math.exp(-this._Charge*timeWalk.getDoubleValue("lambda", 1,1,ICOMPONENT)); + twCorr = timeWalk.getDoubleValue("amplitude", 1,1,REFCOMPONENT)*Math.exp(-this._Charge*timeWalk.getDoubleValue("lambda", 1,1,REFCOMPONENT)); } this.set_Time(time -(FTCALConstantsLoader.CRYS_LENGTH-cluster.getDoubleValue("depth_z", 1,1,0))/FTCALConstantsLoader.VEFF - -timeOffsets.getDoubleValue("time_offset", 1,1,ICOMPONENT)-twCorr); + -timeOffsets.getDoubleValue("time_offset", 1,1,REFCOMPONENT)-twCorr); // if(this.get_Edep()>0.1) System.out.println(ICOMPONENT + " " + this._TDC + " " + // FTCALConstantsLoader.TIMECONVFAC + " " + FTCALConstantsLoader.time_offset[0][0][ICOMPONENT-1] + " " + // this.get_Time()); @@ -208,7 +209,7 @@ public final void set_ClusID(int _ClusIndex) { public static boolean passHitSelection(FTCALHit hit, IndexedTable thresholds) { // a selection cut to pass the hit. - if(hit.get_Edep() > thresholds.getDoubleValue("thresholdHit", 1,1,hit.get_COMPONENT())) { + if(hit.get_Edep() > thresholds.getDoubleValue("thresholdHit", 1,1,REFCOMPONENT)) { return true; } else { return false; diff --git a/reconstruction/ft/src/main/java/org/jlab/rec/ft/cal/FTCALReconstruction.java b/reconstruction/ft/src/main/java/org/jlab/rec/ft/cal/FTCALReconstruction.java index fd9f113a1d..c4b4ad4f5e 100644 --- a/reconstruction/ft/src/main/java/org/jlab/rec/ft/cal/FTCALReconstruction.java +++ b/reconstruction/ft/src/main/java/org/jlab/rec/ft/cal/FTCALReconstruction.java @@ -216,7 +216,7 @@ public List readRawHitsHipo(DataEvent event, IndexedTable charge2Energ int adc = bankDGTZ.getInt("ADC",row); float time = bankDGTZ.getFloat("time",row); if(ilayer==0) ilayer=1; // fix for wrong layer in TT - if(adc!=-1 && time!=-1 && status.getIntValue("status", isector, ilayer, icomponent)==0){ + if(adc!=-1 && time!=-1 && status.getIntValue("status", isector, ilayer, FTCALHit.REFCOMPONENT)==0){ FTCALHit hit = new FTCALHit(bankDGTZ.trueIndex(row),icomponent, adc, time, charge2Energy, timeOffsets, timeWalk, cluster); //////////////////////////////////////////////////////////////////////////////////////////////////// diff --git a/reconstruction/urwell/src/main/java/org/jlab/rec/urwell/reader/URWellCluster.java b/reconstruction/urwell/src/main/java/org/jlab/rec/urwell/reader/URWellCluster.java new file mode 100644 index 0000000000..bbba97d199 --- /dev/null +++ b/reconstruction/urwell/src/main/java/org/jlab/rec/urwell/reader/URWellCluster.java @@ -0,0 +1,111 @@ +package org.jlab.rec.urwell.reader; + +import org.jlab.detector.base.DetectorDescriptor; +import org.jlab.detector.base.DetectorType; +import org.jlab.geom.prim.Point3D; +import org.jlab.geom.prim.Line3D; +import org.jlab.service.urwell.URWellConstants; + +/** + * + * @author Tongtong Cao + */ + +public class URWellCluster { + + private int id = -1; + private DetectorDescriptor desc = new DetectorDescriptor(DetectorType.URWELL); + private int size = 0; + private double energy = 0; + private double time = 0; + private int crossIndex = -1; + private double stereo = 10.; + private Line3D line = new Line3D(); + private Line3D lineLocal = new Line3D(); + + public URWellCluster(int id, int sector, int layer, int component, int size, double energy, double time, Point3D pointOrigin, Point3D pointEnd) { + this.id = id; + this.desc.setSectorLayerComponent(sector, layer, component); + this.size = size; + this.energy = energy; + this.time = time; + this.line.set(pointOrigin, pointEnd); + + Point3D pointOriginLocal = new Point3D(); + pointOriginLocal.copy(pointOrigin); + pointOriginLocal.rotateZ(Math.toRadians(-60 * (sector - 1))); + pointOriginLocal.rotateY(Math.toRadians(-25)); + Point3D pointEndLocal = new Point3D(); + pointEndLocal.copy(pointEnd); + pointEndLocal.rotateZ(Math.toRadians(-60 * (sector - 1))); + pointEndLocal.rotateY(Math.toRadians(-25)); + lineLocal.set(pointOriginLocal, pointEndLocal); + } + + public int id() { + return this.id; + } + + public int sector() { + return this.desc.getSector(); + } + + public int layer() { + return this.desc.getLayer(); + } + + public int strip() { + return this.desc.getComponent(); + } + + public int size() { + return size; + } + + public double energy() { + return energy; + } + + public double time() { + return time; + } + + public int getCrossIndex() { + return crossIndex; + } + + public void setCrossIndex(int crossIndex) { + this.crossIndex = crossIndex; + } + + public double getStereo(){ + return stereo; + } + + public Line3D getLineLocal(){ + return lineLocal; + } + + public Line3D getLine(){ + return line; + } + + /** + * + * @return cluster info. about location and number of hits contained in it + */ + @Override + public String toString() { + String str = "Cluster: Index " + this.id() + + ", Layer " + this.layer() + + ", Seed " + this.strip() + + ", Size " + this.size() + + ", originLocX " + String.format("%.4f", this.getLineLocal().origin().x()) + + ", originLocY " + String.format("%.4f", this.getLineLocal().origin().y()) + + ", originLocZ " + String.format("%.4f", this.getLineLocal().origin().z()) + + ", endLocX " + String.format("%.4f", this.getLineLocal().end().x()) + + ", endLocY " + String.format("%.4f", this.getLineLocal().end().y()) + + ", endLocZ " + String.format("%.4f", this.getLineLocal().end().z()); + return str; + } +} diff --git a/reconstruction/urwell/src/main/java/org/jlab/rec/urwell/reader/URWellCross.java b/reconstruction/urwell/src/main/java/org/jlab/rec/urwell/reader/URWellCross.java new file mode 100644 index 0000000000..1783aba777 --- /dev/null +++ b/reconstruction/urwell/src/main/java/org/jlab/rec/urwell/reader/URWellCross.java @@ -0,0 +1,189 @@ +package org.jlab.rec.urwell.reader; + +import java.util.List; +import org.jlab.detector.base.DetectorDescriptor; +import org.jlab.detector.base.DetectorType; +import org.jlab.geom.prim.Point3D; +import org.jlab.service.urwell.URWellConstants; +/** + * + * @author Tongtong Cao + */ + +public class URWellCross { + + private DetectorDescriptor desc = new DetectorDescriptor(DetectorType.URWELL); + private Point3D global; + private Point3D local; // Points in local coordinates + private int region; + private double energy = 0; + private double time = 0; + private int id = -1; + private int cluster1 = -1; + private int cluster2 = -1; + private int status = -1; + private int tid = -1; // Track id; + private URWellCluster cls1 = null; + private URWellCluster cls2 = null; + + public URWellCross(int id, int sector, int region, double x, double y, double z, double energy, double time, int cluster1, int cluster2, int status) { + this.id = id; + this.desc.setSectorLayerComponent(sector, 0, 0); + this.region = region; + this.global = new Point3D(x, y, z); + this.local = new Point3D(x, y, z); + local.rotateZ(Math.toRadians(-60 * (sector - 1))); + local.rotateY(Math.toRadians(-25)); + this.energy = energy; + this.time = time; + this.cluster1 = cluster1; + this.cluster2 = cluster2; + this.status = status; + } + + public URWellCross(int id, int tid, int sector, int region, double x, double y, double z, double x_local, double y_local, double z_local, double energy, double time, int cluster1, int cluster2, int status) { + this.id = id; + this.tid = tid; + this.desc.setSectorLayerComponent(sector, 0, 0); + this.region = region; + this.global = new Point3D(x, y, z); + this.local = new Point3D(x_local, y_local, z_local); + this.energy = energy; + this.time = time; + this.cluster1 = cluster1; + this.cluster2 = cluster2; + this.status = status; + } + + /** + * return track id + */ + + public int get_tid(){ + return tid; + } + + /** + * @param tid track id + */ + + public void set_tid(int tid){ + this.tid = tid; + } + + public int id() { + return this.id; + } + + public int sector() { + return this.desc.getSector(); + } + + public int region() { + return this.region; + } + + public Point3D position() { + return this.global; + } + + public Point3D local() { + return this.local; + } + + public double energy() { + return energy; + } + + public double time() { + return time; + } + + public void setClusterIndex1(int cluster) { + this.cluster1 = cluster; + } + + public void setClusterIndex2(int cluster) { + this.cluster2 = cluster; + } + + public int cluster1() { + return this.cluster1; + } + + public int cluster2() { + return this.cluster2; + } + + public int status() { + return this.status; + } + + public void setCluster1(URWellCluster cluster) { + cls1 = cluster; + } + + public void setCluster2(URWellCluster cluster) { + cls2 = cluster; + } + + public URWellCluster getCluster1() { + return cls1; + } + + public URWellCluster getCluster2() { + return cls2; + } + + public void setCluster1(List urClusters) { + URWellCluster cluster = null; + + for (URWellCluster cl : urClusters) { + if (cl.id() == cluster1) { + cluster = cl; + break; + } + } + + cls1 = cluster; + } + + public void setCluster2(List urClusters) { + URWellCluster cluster = null; + + for (URWellCluster cl : urClusters) { + if (cl.id() == cluster2) { + cluster = cl; + break; + } + } + + cls2 = cluster; + } + + public URWellCluster getCluster1(List urClusters) { + URWellCluster cluster = null; + + for (URWellCluster cl : urClusters) { + if (cl.id() == cluster1) { + cluster = cl; + break; + } + } + + return cluster; + } + + public URWellCluster getCluster2(List urClusters) { + URWellCluster cluster = null; + + for (URWellCluster cl : urClusters) { + if (cl.id() == cluster2) { + cluster = cl; + break; + } + } + + return cluster; + } +} diff --git a/reconstruction/urwell/src/main/java/org/jlab/rec/urwell/reader/URWellHit.java b/reconstruction/urwell/src/main/java/org/jlab/rec/urwell/reader/URWellHit.java new file mode 100644 index 0000000000..33fa88e9f4 --- /dev/null +++ b/reconstruction/urwell/src/main/java/org/jlab/rec/urwell/reader/URWellHit.java @@ -0,0 +1,43 @@ +package org.jlab.rec.urwell.reader; + +import org.jlab.detector.base.DetectorDescriptor; +import org.jlab.detector.base.DetectorType; + +/** + * + * @author Tongtong Cao + */ + +public class URWellHit { + + private DetectorDescriptor desc = new DetectorDescriptor(DetectorType.URWELL); + private double energy = 0; + private double time = 0; + + public URWellHit(int sector, int layer, int component, double energy, double time) { + this.desc.setSectorLayerComponent(sector, layer, component); + this.energy = energy; + this.time = time; + } + + public int sector() { + return this.desc.getSector(); + } + + public int layer() { + return this.desc.getLayer(); + } + + public int strip() { + return this.desc.getComponent(); + } + + public double energy() { + return energy; + } + + public double time() { + return time; + } + +} diff --git a/reconstruction/urwell/src/main/java/org/jlab/rec/urwell/reader/URWellReader.java b/reconstruction/urwell/src/main/java/org/jlab/rec/urwell/reader/URWellReader.java new file mode 100644 index 0000000000..2723cd7032 --- /dev/null +++ b/reconstruction/urwell/src/main/java/org/jlab/rec/urwell/reader/URWellReader.java @@ -0,0 +1,144 @@ +package org.jlab.rec.urwell.reader; + +import java.util.ArrayList; +import java.util.List; +import org.jlab.geom.prim.Point3D; +import org.jlab.io.base.DataBank; +import org.jlab.io.base.DataEvent; + +/** + * + * @author Tongtong Cao + */ +public class URWellReader{ + + private final List urHits = new ArrayList<>(); + private final List urClusters = new ArrayList<>(); + private final List urCrosses = new ArrayList<>(); + + + public URWellReader(DataEvent event) { + + if(event.hasBank("URWELL::hits")) + this.readHits(event.getBank("URWELL::hits")); + if(event.hasBank("URWELL::clusters")) + this.readClusters(event.getBank("URWELL::clusters")); + if(event.hasBank("URWELL::crosses")) + this.readCrosses(event.getBank("URWELL::crosses")); + } + + public List getUrwellHits() { + return urHits; + } + + public List getUrwellClusters() { + return urClusters; + } + + public List getUrwellCrosses() { + return urCrosses; + } + + public List getUrwellR1Crosses() { + List urCrossesR1 = new ArrayList<>(); + for (URWellCross crs : urCrosses) { + if (crs.region() == 1) { + urCrossesR1.add(crs); + } + } + return urCrossesR1; + } + + public final void readHits(DataBank bank) { + + for(int i=0; i getStrips(DataEvent event, URWellStripFactory fa strip.setEnergy(strip.ADC*URWellConstants.ADCTOENERGY); strip.setTime(strip.TDC*URWellConstants.TDCTOTIME); strip.setLine(factory.getStrip(sector, layer, comp)); - strip.setChamber(factory.getChamberIndex(comp)+1); + strip.setChamber(factory.getChamberIndex(layer, comp)+1); strip.setStatus(0); if(strip.getEnergy()>URWellConstants.THRESHOLD) strips.add(strip);