Updated target geometry and material handling in CVT Eloss

common-tools/clas-detector/src/main/java/org/jlab/detector/base/GeometryFactory.java

@@ -108,6 +108,7 @@ public static ConstantProvider getConstants(DetectorType type, int run, String v
 
         if(type==DetectorType.TARGET){
             provider.loadTable("/geometry/target");
+            provider.loadTable("/geometry/materials/target");
         }
 
         if(type==DetectorType.FMT){

common-tools/clas-geometry/src/main/java/org/jlab/geom/prim/Cylindrical3D.java

@@ -403,6 +403,46 @@ public int intersectionSegment(Line3D line, List<Point3D> intersections) {
         return 0; // no intersections, disjoint
     }
 
+    /**
+     * Compute the intersections of a ray within the 3D volume
+     * @param line
+     * @param intersections
+     * @return ray path length within the volume
+     */
+    public double intersectionLength(Line3D line, List<Point3D> intersections) {
+        if(this.baseArc.theta()<2*Math.PI) {
+            throw new UnsupportedOperationException("Not supported yet."); 
+        }
+        else {
+            int count = this.intersectionRay(line, intersections);
+            // add intersections with bases
+            Plane3D base0 = new Plane3D(this.baseArc.center(), this.baseArc.normal());
+            Plane3D base1 = new Plane3D(this.highArc().center(), this.baseArc.normal());
+
+            Point3D intersect0 = new Point3D();
+            if(base0.intersectionRay(line, intersect0)>0 &&
+               intersect0.distance(base0.point())<this.baseArc.radius()) {
+                count++;
+                intersections.add(intersect0);
+            }
+
+            Point3D intersect1 = new Point3D();
+            if(base1.intersectionRay(line, intersect1)>0 &&
+               intersect1.distance(base1.point())<this.baseArc.radius()) {
+                count++;
+                intersections.add(intersect1);
+            }
+
+            if(count==2) {
+                return intersections.get(0).distance(intersections.get(1));
+            }
+            else if(count==1) {
+                return line.origin().distance(intersections.get(0));
+            }
+        }
+        return 0;
+    }
+
     /**
      * Returns true if the given point is on the surface of this cylindrical 
      * segment.

common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/SVT/SVTConstants.java

@@ -108,17 +108,7 @@ public class SVTConstants {
 	public static double SECTORLEN;
          public static final double SIDETOL   = 1.0;  // extra width used in defining the module corners
 	public static final double LENGTHTOL = 10.0; // extra length for track intersection
-
-        // tungsten shield
-        public static double TSHIELDRMIN    = 51;
-        public static double TSHIELDRMAX    = 51.051;
-        public static double TSHIELDLENGTH  = 360;
-        public static double TSHIELDZPOS    = -50;
-        public static double TSHIELDRADLEN  = 6.76/19.3 *10; // X0(g/cm2) / density(g/cm3) * 10; 
-        public static double TSHIELDZOVERA  = 0.40252;
-        public static double TSHIELDRHO     = 19.3E-3; // g/mm3
-        public static double TSHIELDI       = 727;     // eV
-
+
          // faraday cup cage
          public static double[] FARADAYCAGERMIN    = new double[4];
          public static double[] FARADAYCAGERMAX    = new double[4];

common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/Material.java

@@ -74,5 +74,12 @@ public double getEloss(double p, double mass) {
         return dE;
     }
 
-
+    @Override
+    public String toString() {
+        String s = "Material: ";
+        s = s + String.format("Name=%s thickness=%.9f density=%.9f Z/A=%.6f x0=%9f IeV=%.9f",
+                               this.getName(), this.getThickness(),this.getDensity(),
+                               this.getZoverA(),this.getX0(), this.getIeV());
+        return s;
+    }
 }

common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/Surface.java

@@ -7,6 +7,7 @@
 import org.jlab.geom.prim.Point3D;
 import org.jlab.geom.prim.Arc3D;
 import org.jlab.geom.prim.Cylindrical3D;
+import org.jlab.geom.prim.Line3D;
 import org.jlab.geom.prim.Transformation3D;
 import org.jlab.geom.prim.Vector3D;
 
@@ -24,6 +25,7 @@ public class Surface implements Comparable<Surface> {
     public Point3D finitePlaneCorner1;
     public Point3D finitePlaneCorner2;
     public Cylindrical3D cylinder;
+    public Cylindrical3D lineVolume;
     private Transformation3D toGlobal = new Transformation3D();
     private Transformation3D toLocal  = new Transformation3D();
     public Arc3D arc;
@@ -111,11 +113,19 @@ public Surface(Point3D endPoint1, Point3D endPoint2, double accuracy) {
         swimAccuracy = accuracy;
     }
 
+    public Surface(Point3D endPoint1, Point3D endPoint2, Cylindrical3D volume, double accuracy) {
+        type = Type.LINE;
+        lineEndPoint1 = endPoint1;
+        lineEndPoint2 = endPoint2;
+        lineVolume = volume;
+        swimAccuracy = accuracy;
+    }
+
     @Override
     public String toString() {
         String s = "Surface: ";
-        s = s + String.format("Type=%s Index=%d  Layer=%d  Sector=%d  Emisphere=%.1f X0=%.4f  Z/A=%.4f  Error=%.4f Passive=%b",
-                               this.type.name(), this.getIndex(),this.getLayer(),this.getSector(),this.hemisphere,this.getToverX0(),
+        s = s + String.format("Type=%s Index=%d  Layer=%d  Sector=%d  Emisphere=%.1f thickness=%.4f X0=%.4f Z/A=%.4f Error=%.4f Passive=%b",
+                               this.type.name(), this.getIndex(),this.getLayer(),this.getSector(),this.hemisphere,this.getThickness(),this.getToverX0(),
                                this.getZoverA(),this.getError(), this.passive);
         if(type==Type.PLANEWITHSTRIP) {
             s = s + "\n\t" + this.plane.toString();
@@ -127,9 +137,15 @@ else if(type==Type.CYLINDERWITHSTRIP) {
             s = s + "\n\t" + this.cylinder.toString();
             s = s + "\n\t" + this.strip.toString();
         }
+        else if(type==Type.CYLINDERWITHLINE) {
+            s = s + "\n\t" + this.cylinder.toString();
+            s = s + "\n\t" + this.lineEndPoint1.toString();
+            s = s + "\n\t" + this.lineEndPoint2.toString();
+        }
         else if(type==Type.LINE) {
             s = s + "\n\t" + this.lineEndPoint1.toString();
             s = s + "\n\t" + this.lineEndPoint2.toString();
+            if(this.lineVolume!=null) s = s + "\n\t" + this.lineVolume.toString();
         }
         return s;
     }
@@ -195,6 +211,19 @@ public void addMaterial(String name, double thickness, double density, double Zo
         this.materials.add(new Material(name, thickness, density, ZoverA, X0, IeV, unit));
     }
 
+    public double getRadius() {
+        if(this.cylinder!=null)
+            return this.cylinder.baseArc().radius();
+        else if(this.lineVolume!=null) 
+            return this.lineVolume.baseArc().radius();
+        else
+            return 0;
+    }
+
+    public double getLength() {
+        return this.cylinder.height();
+    }
+
     public double getThickness() {
         double t = 0;
         for(Material m : this.materials) {
@@ -221,56 +250,89 @@ public double getZoverA() {
         return ZA/RhoX;
     }
 
-    public double getLocalDir(Vector3D dir) {
+    public double getTrackLength(Point3D pos, Vector3D dir) {
+        return this.getTrackLength(pos, dir, 0);
+    }
+
+    public double getTrackLength(Point3D pos, Vector3D dir, int materialIndex) {
         if(this.type!=Type.PLANEWITHSTRIP && 
            this.type!=Type.CYLINDERWITHSTRIP && 
            this.type!=Type.LINE) 
-           return 1;
+           return this.getMaterials().get(materialIndex).getThickness();
         else {
             if(this.type==Type.PLANEWITHSTRIP) {
                 Vector3D norm = this.plane.normal();
-                return Math.abs(norm.dot(dir));
+                return this.getMaterials().get(materialIndex).getThickness()/Math.abs(norm.dot(dir));
             }
             else if(this.type==Type.CYLINDERWITHSTRIP) {
                 Vector3D axis = this.cylinder.getAxis().direction().asUnit();
                 dir.sub(dir.projection(axis));
-                return Math.abs(dir.mag());
+                return this.getMaterials().get(materialIndex).getThickness()/Math.abs(dir.mag());
             }
             else if(this.type==Type.LINE) {
-                Vector3D axis = this.lineEndPoint1.vectorTo(this.lineEndPoint2).asUnit();
-                dir.sub(dir.projection(axis));
-                return Math.abs(dir.mag());
+                if(this.lineVolume==null) {
+                    Vector3D line = this.lineEndPoint1.vectorTo(this.lineEndPoint2).asUnit();
+                    dir.sub(dir.projection(line));
+                    return this.getMaterials().get(materialIndex).getThickness()/Math.abs(dir.mag());
+                }
+                else {  
+                    Line3D track = new Line3D(pos, dir);
+                    List<Point3D> intersections = new ArrayList<>();
+                    double trackLength = this.lineVolume.intersectionLength(track, intersections);
+                    if(materialIndex==0) { 
+                        return trackLength;
+                    }
+                    else {
+                        for(Point3D p : intersections) {
+                            if(this.lineVolume.isOnSurface(p)){
+                                Vector3D axis = this.lineVolume.getAxis().direction().asUnit();
+                                dir.sub(dir.projection(axis));
+                                return this.getMaterials().get(materialIndex).getThickness()/Math.abs(dir.mag());
+                            }
+                        }
+                        return 0;
+                    }
+                }
             }
             return 0;
         }
     }    
 
     public double getEloss(double p, double mass) {
-        double dE=0;
-        for(Material m : this.materials) {
+        double dE = 0;
+        for(int im=0; im<this.getMaterials().size(); im++) {
+            Material m = this.getMaterials().get(im);
             dE += m.getEloss(p, mass);
-        }
+        }    
         return dE;
     }
 
-    public double getEloss(Vector3D mom, double mass, int dir) {
-        double cosDir = this.getLocalDir(mom.asUnit());
+    public double getEloss(Point3D pos, Vector3D mom, double mass) {
+        double dE = 0;
+        double p = mom.mag();
+        for(int im=0; im<this.getMaterials().size(); im++) {
+            Material m = this.getMaterials().get(im);
+            dE += m.getEloss(p, mass)*this.getTrackLength(pos, mom.asUnit(), im)/m.getThickness();
+        }    
+        return dE;
+    }
+
+    public double getElossScale(Point3D pos, Vector3D mom, double mass, int dir) {
         double scale = 0;
-        if(cosDir!=0) {
-            double dE = -dir*this.getEloss(mom.mag(), mass)/cosDir;
-            double Ecorr = Math.sqrt(mom.mag2() + mass*mass) + dE;
-            if(Ecorr>mass) scale = Math.sqrt(Ecorr*Ecorr - mass*mass)/mom.mag();
-            mom.scale(scale);
-        }
+        double dE = -dir*this.getEloss(pos, mom, mass);
+        double Ecorr = Math.sqrt(mom.mag2() + mass*mass) + dE;
+        if(Ecorr>mass) scale = Math.sqrt(Ecorr*Ecorr - mass*mass)/mom.mag();
+        mom.scale(scale);
         return scale;
     }
 
-    public double getDx(Vector3D mom) {
-        double cosDir = this.getLocalDir(mom.asUnit());
-        if(cosDir!=0)
-          return this.getThickness()/cosDir;
-        else
-            return 0;
+    public double getDx(Point3D pos, Vector3D mom) {
+        double dx = 0;
+        for(int im=0; im<this.getMaterials().size(); im ++) {
+            Material m = this.getMaterials().get(im);
+            dx += this.getTrackLength(pos, mom.asUnit(), im);
+        }
+        return dx;
     }
 
     public double getThetaMS(double p, double mass, double cosEntranceAngle) {

common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/helical/StateVecs.java

@@ -166,11 +166,11 @@ public void corrForEloss(int dir, StateVec vec, AMeasVecs mv) {
         if(this.straight || mass<0) return;
 
         Surface  surf = mv.measurements.get(vec.k).surface;
-        double pScale = surf.getEloss(vec.getMomentum(), mass, dir);
+        double pScale = surf.getElossScale(vec.getPosition(), vec.getMomentum(), mass, dir);
         if(pScale>0) {
             vec.kappa = vec.kappa/pScale;
-            vec.energyLoss = surf.getEloss(vec.getMomentum().mag(), mass);
-            vec.dx = surf.getDx(vec.getMomentum());
+            vec.energyLoss = surf.getEloss(vec.getPosition(), vec.getMomentum(), mass);
+            vec.dx = surf.getDx(vec.getPosition(), vec.getMomentum());
             vec.updateFromHelix();
         }
         else {

reconstruction/cvt/src/main/java/org/jlab/rec/cvt/Constants.java

@@ -6,8 +6,6 @@
 import java.util.logging.Logger;
 import org.jlab.clas.pdg.PhysicsConstants;
 import org.jlab.clas.swimtools.Swim;
-import org.jlab.clas.tracking.kalmanfilter.Material;
-import org.jlab.clas.tracking.kalmanfilter.Units;
 
 import org.jlab.clas.tracking.utilities.MatrixOps.Libr;
 import org.jlab.rec.cvt.svt.SVTParameters;
@@ -63,7 +61,7 @@ public static Constants getInstance() {
     public boolean   useOnlyBMTC50PercTruthHits = false;
     public boolean   useOnlyBMTZ50PercTruthHits = false;
     public boolean   preElossCorrection = true;
-    private Material targetMaterial = LH2;
+    private String   targetType = "";
     public Libr      KFMatrixLibrary;
     private int svtmaxclussize = 30;
     private int bmtcmaxclussize =30;
@@ -100,14 +98,7 @@ public static Constants getInstance() {
     public static final double COSMICSMINRESIDUALZ =  12; // in mm
 
     public static final int SEEDFITITERATIONS = 5;
-
-    private static final Material LH2 = new Material("LH2", 8.85, 0.0708E-3, 0.99212, 8904.0, 21.8, Units.MM);
-    private static final Material LD2 = new Material("LD2", 8.85, 0.1638E-3, 0.49650, 7691.0, 21.8, Units.MM);
-    public  static final Material TARGETKAPTON = new Material("Kapton", 50E-3, 1.42E-3, 0.501, 285.7, 79.6, Units.MM);
-    public  static final Material TARGETRHOACELL = new Material("Rhoacell", 10.4, 0.1E-3, 0.5392, 1000, 93.0, Units.MM);
-    public  static final Material SCINTILLATOR = new Material("Scintillator", 1, 1.03E-3, 0.54141, 439.0, 64.7, Units.MM);
-    public  static final Material VACUUM = new Material("Vacuum", 1, 0, 1, Double.POSITIVE_INFINITY, 100, Units.MM);
-
+
     public static boolean KFFILTERON = true;
     public static boolean INITFROMMC = false;
     public static int     KFITERATIONS = 5;
@@ -133,16 +124,15 @@ public boolean kfBeamSpotConstraint() {
     }
 
     public void setTargetMaterial(String material) {
-        if(material.equalsIgnoreCase("LH2"))
-            targetMaterial = LH2;
-        else if(material.equalsIgnoreCase("LD2") )
-            targetMaterial = LD2;
+        if(!material.equalsIgnoreCase("LH2") &&
+           !material.equalsIgnoreCase("LD2") )
+            System.out.println("Unknown target material " + material + ", keeping current setting " + targetType);
         else
-            System.out.println("Unknown target material " + material + ", keeping current setting " + targetMaterial.getName());
+            targetType = material;
     }
 
-    public Material getTargetMaterial() {
-        return targetMaterial;
+    public String getTargetType() {
+        return targetType;
     }
 
     public int getRmReg() {