lcsim-contrib/src/main/java/org/lcsim/contrib/Mbussonn
diff -N HelixUtils.java
--- HelixUtils.java 19 Jun 2009 19:11:19 -0000 1.1
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,430 +0,0 @@
-/*
- * HelixUtils.java
- *
- * Created on February 1, 2008, 5:48 PM
- *
- */
-
-package org.lcsim.contrib.Mbussonn;
-
-import hep.physics.vec.BasicHep3Vector;
-import hep.physics.vec.Hep3Vector;
-import hep.physics.matrix.BasicMatrix;
-import hep.physics.matrix.Matrix;
-
-import java.util.ArrayList;
-import java.util.List;
-
-import org.lcsim.fit.circle.CircleFit;
-import org.lcsim.fit.helicaltrack.HelicalTrackFit;
-import org.lcsim.fit.helicaltrack.HelicalTrackHit;
-import org.lcsim.fit.helicaltrack.TrackDirection;
-
-/**
- * Assorted helix utilities that operate on HelicalTrackFits.
- *
- * Since these utilities are typically called many times for a
- * given HelicalTrackFit (or in some cases for a given CircleFit),
- * a number of derived quantities used by the methods in this class
- * are cached to minimize re-computation. Since a HelicalTrackFit
- * (or CircleFit) is immutable, these quantities only need to be
- * calculated once.
- * @author Richard Partridge
- * @version 1.0
- */
-public class HelixUtils {
- private static double _minslope = 1.e-6;
- private static double _maxpath = 2400;
- private static double _epsilon = 0.01; // value below witch a z vector Or an angle is considered as null
-
- /**
- * Creates a new instance of HelixUtils.
- */
- public HelixUtils() {
- }
-
- /**
- * Return the x-y path length between two HelicalTrackHits.
- * @param cfit CircleFit to be used in calculating the path length
- * @param hit1 first hit
- * @param hit2 second hit
- * @return path length from hit1 to hit2
- */
- public static double PathLength(CircleFit cfit, HelicalTrackHit hit1, HelicalTrackHit hit2) {
- // Return the path length between hit1 and hit2
- return PathCalc(xc(cfit), yc(cfit), RC(cfit), hit1.x(), hit1.y(), hit2.x(), hit2.y());
- }
-
- /**
- * Return the x-y path length from the DCA.
- * @param cfit CircleFit to be used in calculating the path length
- * @param hit hit to be used for the path length calculation
- * @return path length from the DCA to the hit
- */
- public static double PathLength(CircleFit cfit, HelicalTrackHit hit) {
- // Return the path length from the DCA
- return PathCalc(xc(cfit), yc(cfit), RC(cfit), x0(cfit), y0(cfit), hit.x(), hit.y());
- }
-
- /**
- * Return the x-y path length from the DCA to a HelicalTrackHit.
- * @param helix HelicalTrackFit to be used in calculating the path length
- * @param hit hit to be used for the path length calculation
- * @return path length from the DCA to the hit
- */
- public static double PathLength(HelicalTrackFit helix, HelicalTrackHit hit) {
- // Return the path length from the DCA
- return PathCalc(helix.xc(), helix.yc(), helix.R(), helix.x0(), helix.y0(), hit.x(), hit.y());
- }
-
- /**
- * Return the x-y path length to a z-plane.
- * @param helix HelicalTrackFit to be used in calculating the path length
- * @param z location of z-plane
- * @return path length from the DCA to the z-plane
- */
- public static double PathToZPlane(HelicalTrackFit helix, double z) {
- // Find the z distace between the DCA and the z plane and calculate the x-y path length
- double zdist = z - helix.z0();
- double safeslope = helix.slope();
- if (Math.abs(safeslope) < _minslope) safeslope = _minslope * Math.signum(safeslope);
- return zdist / safeslope;
- }
-
- /**
- * Return a list of x-y path lengths to a cylinder centered on the z axis.
- * @param helix HelicalTrackFit to be used in calculating the path length
- * @param r desired radius
- * @param smax maximum path length to be considered
- * @param mxint Maximum number of intersections
- * @return list of path lengths
- */
- public static List<Double> PathToCylinder(HelicalTrackFit helix, double r, double smax, int mxint) {
- // Create a list to hold the path lengths
- List<Double> pathlist = new ArrayList<Double>();
- // Retrieve helix dca and RC
- double dca = helix.dca();
- double RC = helix.R();
- // Find cos(dphi) for dphi from the point of closest approach to the first crossing
- double cdphi = 1 - (r * r - dca * dca) / (2. * RC * (RC - dca));
- // Make sure we have a physical intersection
- if (Math.abs(cdphi) < 1.) {
- // Calculate dphi and the path length to the first crossing
- double dphi = Math.acos(cdphi);
- Double s = dphi * Math.abs(RC);
- // Loop over crossings until we exceed one of the limits
- while (s < smax && pathlist.size() < mxint) {
- // Add this odd-numbered crossing to the list
- pathlist.add(s);
- // Advance to the next even-numbered crossing
- s += 2. * (Math.PI - dphi) * Math.abs(RC);
- // Check to see if we should add it
- if (s < smax && pathlist.size() < mxint) pathlist.add(s);
- // Add this even-numbered crossing to the list
- s += 2. * dphi * Math.abs(RC);
- }
- }
- return pathlist;
- }
-
- /**
- * Return a unit vector giving the track direction at a given point on
- * the helix.
- * @param helix HelicalTrackFit to be used in calculating direction
- * @param s path length to the desired point on helix
- * @return direction unit vector
- */
- public static Hep3Vector Direction(HelicalTrackFit helix, double s) {
- // Calculate the azimuthal direction
- double phi = helix.phi0() - s / helix.R();
- double sth = helix.sth();
- // Calculate the components of the direction unit vector
- double ux = Math.cos(phi) * helix.sth();
- double uy = Math.sin(phi) * helix.sth();
- double uz = helix.cth();
- // Return the direction unit vector
- return new BasicHep3Vector(ux, uy, uz);
- }
-
- /**
- * Return the TrackDirection object for a given point on a helix.
- * This might seem like an odd place to put this code, but putting
- * it here allows the cached helix quantities to be used in constructin
- * the TrackDirection objects.
- * @param helix HelicalTrackFit to use in constructing the TrackDirection
- * @param s path length specifying location on helix
- * @return TrackDirection object for this point on the helix
- */
- public static TrackDirection CalculateTrackDirection(HelicalTrackFit helix, double s) {
- Hep3Vector dir = Direction(helix, s);
- Matrix deriv = DirectionDerivates(helix, dir, s);
- return new TrackDirection(dir, deriv);
- }
-
- /**
- * Return the location in space for a particular point on a helix.
- * @param helix HelicalTrackFit to be used
- * @param s path length
- * @return point in space corresponding to the given path length
- */
- public static Hep3Vector PointOnHelix(HelicalTrackFit helix, double s) {
- // Find the azimuthal direction at this path length
- double RC = helix.R();
- double phi = helix.phi0() - s / RC;
- // Calculate the position on the helix at this path length
- double x = helix.xc() - RC * Math.sin(phi);
- double y = helix.yc() + RC * Math.cos(phi);
- double z = helix.z0() + s * helix.slope();
- // Return the position as a Hep3Vector
- return new BasicHep3Vector(x, y, z);
- }
-
- private static double PathCalc(double xc, double yc, double RC, double x1, double y1, double x2, double y2) {
- // Find the angle between these points measured wrt the circle center
- double phi1 = Math.atan2(y1 - yc, x1 - xc);
- double phi2 = Math.atan2(y2 - yc, x2 - xc);
- double dphi = phi2 - phi1;
- // Make sure dphi is in the valid range (-pi, pi)
- if (dphi > Math.PI) dphi -= 2. * Math.PI;
- if (dphi < -Math.PI) dphi += 2. * Math.PI;
- // Return the arc length
- return -RC * dphi;
- }
-
- /**
- * Return the derivatives of the momentum unit vector with respect to the
- * helix parameters. The direction derivatives are returned in matrix
- * form, with the row giving the cartesian component of the momentum
- * vector and the column giving the helix parameter.
- * @param helix HelicalTrackFit to be used in calculating derivatives
- * @param s path length to the desired point on the helix
- * @return direction derivatives
- */
- private static Matrix DirectionDerivates(HelicalTrackFit helix, Hep3Vector u, double s) {
- // Create the matrix that will hold the derivatives
- BasicMatrix deriv = new BasicMatrix(3,5);
- // Retrieve some helix info
- double cphi0 = Math.cos(helix.phi0());
- double sphi0 = Math.sin(helix.phi0());
- double sth = helix.sth();
- double cth = helix.cth();
- double dca = helix.dca();
- double omega = helix.curvature();
- // Calculate the non-zero derivatives of the direction with respect to the helix parameters
- deriv.setElement(0, helix.curvatureIndex, (u.x() - cphi0 * sth) / omega); // du_x / domega
- deriv.setElement(1, helix.curvatureIndex, (u.y() - sphi0 * sth) / omega); // du_y / domega
- deriv.setElement(0, helix.dcaIndex, -omega * cphi0 * sth); // du_x / ddca
- deriv.setElement(1, helix.dcaIndex, -omega * sphi0 * sth); // du_y / ddca
- deriv.setElement(0, helix.phi0Index, -(1 - dca * omega) * sphi0 * sth); // du_x / dphi0
- deriv.setElement(1, helix.phi0Index, (1 - dca * omega) * cphi0 * sth); // du_y / dphi0
- deriv.setElement(0, helix.slopeIndex, -u.x() * sth * cth); // du_x / dslope
- deriv.setElement(1, helix.slopeIndex, -u.y() * sth * cth); // du_y / dslope
- deriv.setElement(2, helix.slopeIndex, sth*sth*sth); // du_z / dslope
- // Return the derivative matrix
- return deriv;
- }
-/**
- * return true if the helix is intercepting the bounded cylinder
- * @param helix
- * @param r rarius of the cylinder
- * @param zmin lower bound of the cylinder
- * @param zmax higher bound of the cylinder
- * @return
- */
- public static boolean isInterceptingBoundedCylinder(HelicalTrackFit helix,double r,double zmin,double zmax){
- double minpath = PathToZPlane(helix, zmin);//not sure it's very efficient to calculate the maximum path
- double maxpath = PathToZPlane(helix, zmax);
- double path = Math.max(minpath, maxpath);
- List<Double> pathlist = PathToCylinder(helix,r,path,10);
- for(double s:pathlist){
- double z = PointOnHelix(helix,s).z();
- if(z<zmax && z> zmin)
- return true;
- }
- return false;
- }
- /**
- * return true if the helix is intercepting the given disk
- * @param helix
- * @param rmin
- * @param rmax
- * @param z
- * @return
- */
- public static boolean isInterceptingZDisk(HelicalTrackFit helix,double rmin,double rmax, double z ){
- double s = PathToZPlane(helix,z);
- Hep3Vector point = PointOnHelix(helix,s);
- double x = point.x();
- double y = point.y();
- double r = Math.sqrt(x*x+y*y);
- if(r < rmax && r > rmin ){return true;}
- return false;
- }
-
- /**
- * return true if the poitn on the helix at coord z is intercepting a given ZPolygon, the
- * method don't check if the polygone is parallel to a z plane
- * @param helix
- * @param z
- * @param vertices
- * @return true OR false wether the helix intercept or not the polygone
- */
- public static boolean isInterceptingZpolygone(HelicalTrackFit helix,double z,List<Hep3Vector> vertices){
- //todo: check if all vertices are on same z
- double epsilon = Math.PI/360.;
- double s = PathToZPlane(helix,z);
- if(s<0){return false;}
- if(s>_maxpath){return false;}
- Hep3Vector point = PointOnHelix(helix,s);
- double angle = insidePolygone(point,vertices);
- if(Math.abs(angle-2*Math.PI)<epsilon)
- return true;
- return false;
- }
-
-
-
- /**
- * Check if the given helix is intercepting XY plan
- * note, the z coordinate of the XY plane should be 0, but due to numerical
- * approximation, it accept all z value < _epsilon
- * @param helix the Helix
- * @param normal a unitary normal vector to the plan
- * @param orig one point of the plane
- * @return true if the helix intersect at least once with the plane
- */
- public static boolean isInterceptingXYPlane(HelicalTrackFit helix, Hep3Vector normal,Hep3Vector orig){
- if(normal.z()>_epsilon)
- throw new UnsupportedOperationException("Not a XY plan !"+normal);
- double x = -helix.x0()+orig.x();
- double y = -helix.y0()+orig.y();
- double z = -helix.z0()+orig.z();
- double xn= normal.x();
- double yn= normal.y();
- double zn= normal.z();
- double Phip = Math.atan2(yn, xn);
- double dist = (xn*x+yn*y+zn*z);
- double verif = Math.sin(Phip-helix.phi0())+helix.curvature()*dist;
- if(normal.magnitude() < 0.99)
- throw new UnsupportedOperationException("normal vector not unitary :"+normal.magnitude());
- if(Math.abs(verif)>1)
- return false;
- else
- return true ;
- }
-
- /**
- * Check if the given helix is intercepting an XY plane bouded by a
- * list of nodes
- * @param helix
- * @param normal normal vector to the plane
- * @param nodes
- * @return true if th e helix intercept the bounded plane
- **/
- public static boolean isInterceptingBoundedXYPlane(HelicalTrackFit helix,Hep3Vector normal,List<Hep3Vector> nodes){
- Hep3Vector orig = nodes.get(0);
- if(!isInterceptingXYPlane(helix,normal,orig)){return false;}
- double s = PathToXYPlan(helix, normal, orig);
- if(s<0) {return false;}
- if(s>_maxpath){return false;}
- Hep3Vector point = PointOnHelix(helix, s);
- if(Math.abs(insidePolygone(point, nodes)-2*Math.PI)<_epsilon)
- return true;
- return false;
- }
- /**
- * return one path length (projection on the XY plane) to an XY plane
- * the methode is for now only used to check IF an helix is intercepting an
- * XY plane
- * note, the z coordinate of the XY plane should be 0, but due to numerical
- * approximation, it accept all z value < _epsilon
- * @param helix
- * @param normal a UNITARY vector NORMAL to the plan
- * @param origin one point of the plane
- * @return path length
- */
- public static double PathToXYPlan(HelicalTrackFit helix,Hep3Vector normal,Hep3Vector origin){
- if(normal.z()>_epsilon)
- throw new UnsupportedOperationException("Not a XY plan !"+normal);
- double x = -helix.x0()+origin.x();
- double y = -helix.y0()+origin.y();
- double z = -helix.z0()+origin.z();
- double xn= normal.x();
- double yn= normal.y();
- double zn= normal.z();
- double phip = Math.atan2(yn, xn);
- double dist = (xn*x+yn*y+zn*z);
- double sinPlan = Math.sin(phip-helix.phi0())+helix.curvature()*dist;
- double Cs = (Math.asin(sinPlan)-phip+helix.phi0());
- double s=dist/(Math.sqrt(xn*xn+yn*yn)*(2./Cs)*Math.sin(Cs/2.)*Math.cos(Cs/2+phip-helix.phi0())+zn*helix.slope());
- return s;
- }
-
-
- /**
- * adapted from http://local.wasp.uwa.edu.au/~pbourke/geometry/insidepoly/
- * find if a point is inside a given polygon (!!! not polyhedron !!!) in 3Dspace
- * if the point is in the polygone the returned value is 2*PI
- * if the point is outside of the polygone the returned value is 0
- * if the point is not of the plan of the polygone the value is between 0 and 2PI
- * The polygone HAVE to Be convex for the algorythme to work
- *
- * This fonction sum all the angle made by q and 2 consecutives point of the polygone
- *
- * @param q the point to determine if it is or not in the polygon
- * @param p List of edges of the polygone
- * @return 2PI if inside, 0 if outside,intermediate values if not on the plane
- */
- private static double insidePolygone(Hep3Vector q,List<Hep3Vector> p)
- {
- int i;
- double Epsilon = _epsilon;
- double m1,m2;
- double anglesum=0,costheta;
- Hep3Vector p1;
- Hep3Vector p2;
- int n = p.size();
- for (i=0;i<n;i++) {
- p1 = p.get(i);
- p2 = p.get((i+1)%n);
- double x1 = (p1.x()-q.x());
- double y1 = (p1.y()-q.y());
- double z1 = (p1.z()-q.z());
- double x2 = (p2.x()-q.x());
- double y2 = (p2.y()-q.y());
- double z2 = (p2.z()-q.z());
- m1 = Math.sqrt(x1*x1+y1*y1+z1*z1);
- m2 = Math.sqrt(x2*x2+y2*y2+z2*z2);
- if (m1*m2 <= Epsilon){
- return 2*Math.PI; /* We are on a node, consider this inside */
- }
- else{
- costheta = (x1*x2 + y1*y2 + z1*z2)/(m1*m2);
- anglesum += Math.acos(costheta);
- }
- }
- return anglesum;
- }
-
-
- private static double RC(CircleFit cfit) {
- return 1.0 / cfit.curvature();
- }
-
- private static double xc(CircleFit cfit) {
- return cfit.xref() + (RC(cfit) - cfit.dca()) * Math.sin(cfit.phi());
- }
-
- private static double yc(CircleFit cfit) {
- return cfit.yref() - (RC(cfit) - cfit.dca()) * Math.cos(cfit.phi());
- }
-
- private static double x0(CircleFit cfit) {
- return cfit.xref() - cfit.dca() * Math.sin(cfit.phi());
- }
-
- private static double y0(CircleFit cfit) {
- return cfit.yref() + cfit.dca() * Math.cos(cfit.phi());
- }
-
-}
lcsim-contrib/src/main/java/org/lcsim/contrib/Mbussonn/HelixTest
diff -u -r1.3 -r1.4
--- HelixTest.java 29 Jun 2009 18:13:08 -0000 1.3
+++ HelixTest.java 9 Jul 2009 00:43:49 -0000 1.4
@@ -19,15 +19,13 @@
import java.util.Set;
import org.freehep.application.studio.Studio;
import org.freehep.util.FreeHEPLookup;
-import org.lcsim.contrib.Mbussonn.HelixUtils;
-import org.lcsim.contrib.onoprien.util.swim.IntersectionNone;
+import org.lcsim.fit.helicaltrack.HelixUtils;
import org.lcsim.contrib.onoprien.util.swim.ZDisk;
import org.lcsim.contrib.onoprien.util.vector.ConstHep3Vector;
import org.lcsim.detector.*;
import org.lcsim.event.Track;
import org.lcsim.geometry.Detector;
import org.lcsim.util.Driver;
-import org.lcsim.contrib.onoprien.util.swim.Helix.CRep;
import org.lcsim.contrib.sATLAS.TrackReconstructionDriver;
import org.lcsim.detector.solids.*;
import org.lcsim.event.EventHeader;
@@ -37,53 +35,62 @@
import org.lcsim.util.heprep.LCSimHepRepConverter;
/**
- *
+ * Sand box for constructing a driver which find all the intersections between
+ * the material of a detector and an "HelicalTrackFit"
+ * The goal is to create all the "propagator" of this track to be feed in a kalmann filter
* @author matthiasbussonnier
*/
public class HelixTest extends Driver {
- private HepRepFactory factory;
- private HepRepInstanceTree instanceTree;
- private HepRepType typeM;
- private List<HelixPath> helixPath = new LinkedList<HelixPath>();
+ private static boolean alreadyloaded=false;
+ private boolean donothing;
+ private List<HelixPath> helixPath;
public HelixTest()
{
+ if(!alreadyloaded){
+ alreadyloaded=true;
+ System.out.println("loadded one time");
add(new TrackReconstructionDriver());
- Studio app = (Studio) Studio.getApplication();
+ ////////////////////////////////////////////////////////////////
+ //nexts lines are used to "draw" what the helix is intecepting//
+ ////////////////////////////////////////////////////////////////
+ Studio app = (Studio) Studio.getApplication(); //
FreeHEPLookup lookup = app.getLookup();
LCSimHepRepConverter converter = (LCSimHepRepConverter) lookup.lookup(LCSimHepRepConverter.class);
-
- converter.register(new HelixPathConverter());
+ converter.register(new HelixPathConverter()); //
+ ////////////////////////////////////////////////////////////////
+ }
+ else
+ {
+ System.out.println("will do nothing in here");
+ donothing=true;
+ }
}
+ @Override
protected void process(EventHeader event)
{
+ if(donothing){
+ System.out.println("doing nothing in here");
+ return;}
+ helixPath = new LinkedList<HelixPath>();
helixPath.clear();
super.process(event);
- List<Track> tracks = event.getTracks();
- Detector det = event.getDetector();
- Track track = null;
- double xyz[] = {0,0,0};
- Set<Helix> helixlist = new HashSet<Helix>();
- ConstHep3Vector org = new ConstHep3Vector(0.0,0.0,0.0);
- ConstHep3Vector ref =null;
- ConstHep3Vector vec =null;
- ref = new ConstHep3Vector(0,0,0);
+ List<Track> tracks = event.getTracks();
+ Detector det = event.getDetector();
+ Track track = null;
- Helix helix = null;
+
+ HelicalTrackFit hf = null;
if(tracks.size()>0){
- track = tracks.get(0);
- xyz = track.getReferencePoint();
- ref = new ConstHep3Vector(xyz[0],xyz[1],xyz[2]);
- vec = new ConstHep3Vector(track.getPX(),track.getPY(),track.getPZ());
- helix = new Helix(org , vec,track.getTrackParameter(2));
- helix.setReferencePoint(ref);
+ track=tracks.get(0);
+ hf = TrackToHTF(track);
}else{
throw new UnsupportedOperationException("no tracks founded");
}
- helixlist.add(helix);
+
//////////////////////////////////////////////////////////////
//**********************************************************//
//************** some helix for test ***********************//Ê
@@ -110,32 +117,33 @@
//////////////////////////////////////////////////////////////
-// System.out.println(helix);
-// System.out.println("==== Helix ========");
-// System.out.println("curvature = "+helix.getParameter(CRep.C));
-// System.out.println("dca = "+helix.getParameter(CRep.D0));
-// System.out.println("phi0 = "+helix.getParameter(CRep.PHI0));
-// System.out.println("Z0 = "+helix.getParameter(CRep.Z0));
-// System.out.println(" ");
-// System.out.println("path at 20 = "+helix.findLengthToZ(20));
-// System.out.println(" '' at -20 = "+helix.findLengthToZ(-20));
-// System.out.println("slope = "+helix.getParameter(CRep.TAN_LAMBDA));
-// ZDisk disq =new ZDisk(20.,0.,1000000.);
- //System.out.println("x at z:20 = "+helix.intersect(disq).getPosition().x());
- //System.out.println("y at z:20 = "+helix.intersect(disq).getPosition().y());
-// System.out.println("===== HelixFit ====");
- HelicalTrackFit hf = HelixToHTF(helix);
+
+
+ System.out.println("===== Track ===========");
+ System.out.println("curvature = "+track.getTrackParameter(HelicalTrackFit.curvatureIndex));
+ System.out.println("dca = "+track.getTrackParameter(HelicalTrackFit.dcaIndex));
+ System.out.println("phi0 = "+track.getTrackParameter(HelicalTrackFit.phi0Index));
+ System.out.println("Z0 = "+track.getTrackParameter(HelicalTrackFit.z0Index));
+
+ System.out.println("helical fit ");
System.out.println("curvature = "+hf.curvature());
System.out.println("dca = "+hf.dca());
System.out.println("phi0 = "+hf.phi0());
System.out.println("Z0 = "+hf.z0());
- System.out.println("path at 20 = "+HelixUtils.PathToZPlane(hf, 20));
+ System.out.println("======================");
+// System.out.println("path at 20 = "+HelixUtils.PathToZPlane(hf, 20));
//System.out.println("path at -20 = "+HelixUtils.PathToZPlane(hf, -20));
//System.out.println("x at z:20 = "+HelixUtils.isInterceptingZDisk(hf,0.0, 100.0, 20.0));
//System.out.println("slope = "+hf.slope());
//System.out.println("===================");
ISolid trkvol = det.getTrackingVolume().getLogicalVolume().getSolid();
+ ////////////////////////////////////////////////////////////////
+ //not in use yet... smin is the minimum path length to either://
+ // -go outside the dŽtector in the radial direction //
+ // -go outside the detector in the z direction //
+ ////////////////////////////////////////////////////////////////
+ //HelixUtils.setMaxPath(2512);
double smin;
if (trkvol instanceof Tube) {
Tube trktube = (Tube) trkvol;
@@ -143,9 +151,14 @@
double zmax = trktube.getZHalfLength();
System.out.println("Ecal radius = " + rmax);
System.out.println("ECal inner Z = " + zmax);
- double s0= HelixUtils.PathToZPlane(hf, zmax);
- List<Double> s1= HelixUtils.PathToCylinder(hf, s0, zmax, 1);
- //smin = Math.min(s0,s1.get(0));
+ double s0a= HelixUtils.PathToZPlane(hf, zmax);
+ double s0b= HelixUtils.PathToZPlane(hf, -zmax);
+ double s0 = Math.max(s0a, s0b);
+ List<Double> s1= HelixUtils.PathToCylinder(hf, s0, rmax, 1);
+ if(s1.size()!=0)
+ {smin = Math.min(s0,s1.get(0));}
+ else
+ {smin = s0;}
}
else
{
@@ -153,10 +166,15 @@
}
+
+
+
+
+ // Loop over the volumes defined at the compact.xml level
+ // and flatten everything
IPhysicalVolumeNavigator nav = det.getNavigator();
- ILogicalVolume ltrkr = det.getTrackingVolume().getLogicalVolume();
- // Loop over the volumes defined at the compact.xml level
+ ILogicalVolume ltrkr = det.getTrackingVolume().getLogicalVolume();
Set<UniquePV> allVolumes = new HashSet<UniquePV>();
for (IPhysicalVolume pvtree : ltrkr.getDaughters()) {
List<UniquePV> pvflat = Flatten(pvtree, nav);
@@ -165,25 +183,29 @@
}
Set<UniquePV> traversedVolumes = new HashSet<UniquePV>();
int i;
- ////////////////////////////
- ////////////////////////////
- ////////////////////////////
- List<Hep3Vector> thehelixpath = new LinkedList<Hep3Vector>();
- HelicalTrackFit tf = this.HelixToHTF(helix);
- HelixPath hp = new HelixPath();
- for(i=0;i<31;i++)
- thehelixpath.add(HelixUtils.PointOnHelix(tf, i*2000));
+ //////////////////////////////////////////////////////////////////////////
+ ///////////////////////////////////////////////////////////////////////////////
+ /// Graph the track path again, to be sure it's the good one ///
+ // ///
+ List<Hep3Vector> thehelixpath = new LinkedList<Hep3Vector>(); ///
+ HelicalTrackFit trackFit = this.TrackToHTF(track); ///
+ HelixPath hp = new HelixPath(); ///
+ for(i=0;i<31;i++) ///
+ thehelixpath.add(HelixUtils.PointOnHelix(trackFit, i*29)); ///
int[] x={0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30};
- hp.setHepRepVertexOrdering(x);
- hp.setVertices(thehelixpath);
- helixPath.add(hp);
- ////////////////////////////
- ////////////////////////////
- ////////////////////////////
- for(Helix he: helixlist){
- traversedVolumes = this.listOfTraverserVolumes(helix, allVolumes);
+ hp.setHepRepVertexOrdering(x); ///
+ hp.setVertices(thehelixpath); ///
+ helixPath.add(hp); ///
+ // ///
+ /// ///
+ //////////////////////////////////////////////////////////////////////////////
+ //////////////////////////////////////////////////////////////////////////
+
+
+ // try to loop through tracks
+ traversedVolumes = this.listOfTraverserVolumes(trackFit, allVolumes);
System.out.println("==========================================");
- }
+
event.put("HelixPath", helixPath, HelixPath.class, 0);
System.out.println("intersect with "+traversedVolumes.size()+" volumes among "+allVolumes.size());
return;
@@ -191,9 +213,20 @@
+ /////////////////////////////////////////////////////////////////////////
+ //////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+// Somes privates methodes //
+////////////////////////////////////////////////////////////////////////////////
+ //////////////////////////////////////////////////////////////////////////////
+ //////////////////////////////////////////////////////////////////////////
+
- private Set<UniquePV> listOfTraverserVolumes(Helix helix,Set<UniquePV> listeOfVolumes)
+////////////////////////////////////////////////////////////////////////////////
+// create a list of traversed volumes //
+////////////////////////////////////////////////////////////////////////////////
+ private Set<UniquePV> listOfTraverserVolumes(HelicalTrackFit helix,Set<UniquePV> listeOfVolumes)
{
Set<UniquePV> traversedVolumes = new HashSet<UniquePV>();
@@ -203,7 +236,8 @@
{traversedVolumes.add(pv);}
}
catch(UnsupportedOperationException e)
- {System.out.println(e);}
+ {//System.out.println(e);
+ }
}
return traversedVolumes;
}
@@ -212,7 +246,9 @@
-
+////////////////////////////////////////////////////////////////////////////////
+// flatten the tree of physical volumes //
+////////////////////////////////////////////////////////////////////////////////
private List<UniquePV> Flatten(IPhysicalVolume vol, IPhysicalVolumeNavigator nav) {
LinkedList<UniquePV> pvtree = new LinkedList<UniquePV>();
@@ -238,8 +274,10 @@
return pvflat;
}
+
+
/**
- * check i fan helix intersect with a given plan extending a faces ie don't
+ * check if an helix intersect with a given plan extending a faces ie don't
* check yet that the intersection is in the face
* @param hf
* @param faces
@@ -250,7 +288,7 @@
Hep3Vector norm = face.getPlane().getNormal();
List<Point3D> vertices = face.getVertices();
- Point3D vx = vertices.get(1);
+ Point3D vx = vertices.get(0);
double zmin = vx.z();
double zmax = vx.z();
List<Hep3Vector> verticesx = new LinkedList<Hep3Vector>();
@@ -258,18 +296,32 @@
zmin = Math.min(zmin, vert.z());
zmax = Math.max(zmax, vert.z());
verticesx.add(vert);
+ }
+ double zvalue = Math.abs(norm.z());
+
+ if(zvalue <0.02 )
+ {
+ if(HelixUtils.isInterceptingBoundedXYPlane(hf, norm, verticesx))
+ {return true;}
+ }
+ else if(zvalue >0.98){
+// //
+// if(Math.abs(norm.z())<0.98) // deal with vector who should be parralel to z but arn't because of numerical issues
+// {
+// System.out.println("not really z plane...continue ("+norm+")");
+// continue;
+// }//System.out.println("skip z!=0 : "+norm.z()+" "+norm);
+// if(HelixUtils.isInterceptingZpolygone(hf, zmax,verticesx)){//here, z should be close from 1 or -1
+// System.out.println("intercept z polygon...return true");
+// return true;
+// }//modifie that ...test for now
+// else
+// {continue;}
}
- if(Math.abs(norm.z())>0.001){
- if(Math.abs(norm.z())<0.98) // deal with vector who should be parralel to z but arn't because of numerical issues
- {continue;}//System.out.println("skip z!=0 : "+norm.z()+" "+norm);
- if(HelixUtils.isInterceptingZpolygone(hf, zmax,verticesx))//here, z should be close from 1 or -1
- {return true;}//modifie that ...test for now
- else
- {continue;}
- }
- else if(HelixUtils.isInterceptingBoundedXYPlane(hf, norm, verticesx))
- {return true;}
- }
+ else{
+ System.out.println("neither a z nor a xy plan");
+ }
+ }//endfor
return false;
}
@@ -283,38 +335,36 @@
* @return true if the helix intersect the Unique PV
*/
@SuppressWarnings({"unchecked"})
- private boolean intersect(Helix he, UniquePV pv) {
+ private boolean intersect(HelicalTrackFit helicalFit, UniquePV pv) {
ISolid solid = pv.getSolid();
- HelicalTrackFit hf = HelixToHTF(he);
- List<Point3D> vtx;
+ HelicalTrackFit hf = helicalFit;
List<Hep3Vector> vtxx = new LinkedList<Hep3Vector>();
- //System.out.println(solid.getClass());
if(solid instanceof Box||solid instanceof Trd||solid instanceof Trap)
{
AbstractPolyhedron box = (AbstractPolyhedron)solid;
- ITransform3D tr= pv.getLtoGTransform();
+ ITransform3D transformation= pv.getLtoGTransform();
List<Polygon3D> faces =box.getFaces();
- List<Polygon3D> facest = new LinkedList<Polygon3D>();
- for(Polygon3D fc : faces )
+ List<Polygon3D> transformedFaces = new LinkedList<Polygon3D>();
+ //get the "absolute position of the faces"
+ for(Polygon3D face : faces )
{
- Polygon3D pl = fc.transformed(tr);
- facest.add(pl);
- for(Point3D vt : pl.getVertices())// to be albe do draw each face
+ Polygon3D polygon = face.transformed(transformation);
+ transformedFaces.add(polygon);
+ for(Point3D vt : polygon.getVertices())// to be albe do draw each face
{vtxx.add(vt.getHep3Vector());}
-
}
- vtx = box.getVertices();
- //for(Point3D vt : vtx)
- // {vtxx.add(vt.transformed(tr));}
- faces.clear();
- faces=null;
- if(intersectWithFaces(hf,facest)){
+
+ if(intersectWithFaces(hf,transformedFaces)){
+ //////////////////////////////////////////////////////////////////////
+ // draw the face //
HelixPath hp = new HelixPath();
//hp.setHepRepVertexOrdering(box.getHepRepVertexOrdering());
int[] x = {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23};
hp.setHepRepVertexOrdering(x);
hp.setVertices(vtxx);
helixPath.add(hp);
+ // //
+ //////////////////////////////////////////////////////////////////////
return true;
}
return false;
@@ -322,57 +372,17 @@
else if(solid instanceof Tube)
{
Tube tube = (Tube) solid;
-// System.out.println("et un de rayon "+tube.getInnerRadius());
double z0 = pv.getLtoGTransform().getTranslation().z();
-
double zmax = z0 + tube.getZHalfLength();
double zmin = z0 - tube.getZHalfLength();
double rmin = tube.getInnerRadius();
double rmax = tube.getOuterRadius();
-
-
- ZDisk disque = new ZDisk(zmin,rmin,rmax);
- if(!(he.intersect(disque) instanceof IntersectionNone ))
- {
- if(!(HelixUtils.isInterceptingZDisk(hf, rmin, rmax, zmin))){
- System.out.println("disque zmin intersect in Helix, not in HelixUtil !"+he.intersect(disque).getClass());
-// HelixUtils.rayInterceptingZDisk(hf, rmin, rmax, zmin);
-// Hep3Vector pos = he.intersect(new ZDisk(zmin,0.,100000000.)).getPosition();
-// System.out.println("Helix : "+Math.sqrt(pos.x()*pos.x()+pos.y()*pos.y()));
- //System.out.println("----");
- }
- }
-
- if(HelixUtils.isInterceptingZDisk(hf, rmin, rmax, zmin))
- {
- if((he.intersect(disque) instanceof IntersectionNone ))
- System.out.println("disque zmin intersect in HelixUtil, not in Helix !");
- return true;
- }
-
-
-
- disque = new ZDisk(zmax,rmin,rmax);
if(HelixUtils.isInterceptingZDisk(hf, rmin, rmax, zmax))
- {
- if( (he.intersect(disque) instanceof IntersectionNone ))
- System.out.println("disque zmax intersect in HelixUtil, not in Helix !");
- return true;
- }
-
-
-
-
- if(HelixUtils.isInterceptingBoundedCylinder(hf, rmin, zmin, zmax)){
-
- return true;
- }
- else{
- //System.out.println("intercept nothing");
- return false;
- }
-
- //throw new UnsupportedOperationException("tube missed");
+ {return true;}
+ else if(HelixUtils.isInterceptingBoundedCylinder(hf, rmin, zmin, zmax))
+ {return true;}
+ else
+ {return false;}
}
else if(solid instanceof Polycone)
{
@@ -389,32 +399,71 @@
private boolean isCylinder(double rmin, double rmax, double zmin, double zmax) {
return (rmax - rmin) * Math.abs(zmax + zmin) < (zmax - zmin) * (rmax + rmin);
}
+////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+//// convert a track ton an Helical track fit ////
+////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+ private HelicalTrackFit TrackToHTF(Track trk){
- private HelicalTrackFit HelixToHTF(Helix he){
-
- int dcaIndex = 0;
- int phi0Index = 1;
- int curvatureIndex = 2;
- int z0Index = 3;
- int slopeIndex = 4;
- double[] pars={he.getParameter(CRep.D0),
- he.getParameter(CRep.PHI0),
- he.getParameter(CRep.C),
- he.getParameter(CRep.Z0),
- he.getParameter(CRep.TAN_LAMBDA)};
+ double[] pars={trk.getTrackParameter(HelicalTrackFit.dcaIndex),
+ trk.getTrackParameter(HelicalTrackFit.phi0Index),
+ trk.getTrackParameter(HelicalTrackFit.curvatureIndex),
+ trk.getTrackParameter(HelicalTrackFit.z0Index),
+ trk.getTrackParameter(HelicalTrackFit.slopeIndex)};
SymmetricMatrix cov = null;
double[] _chisq = new double[2];
- double _nhchisq = 0.;
int[] _ndf = new int[2];
Map<HelicalTrackHit, Double> _smap = null;
Map<HelicalTrackHit, MultipleScatter> _msmap = null;
HelicalTrackFit helixFakeFit=new HelicalTrackFit(pars, cov,_chisq,_ndf,_smap,_msmap);
return helixFakeFit;
- //HelicalTrackFit(pars, cov,_chisq, _ndf,smap, Map<HelicalTrackHit, MultipleScatter> msmap)
}
+ ////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+
+
/**
* A "struct" holding geometry information about a single physical volume
*/
@@ -430,9 +479,6 @@
double ox,oy,oz = 0;
double size = 0;
}
-
-
-
/**
* A UniquePV is a wrapper around IPhysicalVolumePath which provides
* some convenience methods and caches transformations.