From 5e60fa5c2d237f1b1e2571daa5c21b96e01517b8 Mon Sep 17 00:00:00 2001 From: Chengdong Fu <fucd@ihep.ac.cn> Date: Mon, 20 Jul 2020 14:20:35 +0800 Subject: [PATCH] first import from MarlinUtil --- Utilities/DataHelper/CMakeLists.txt | 14 + Utilities/DataHelper/CMakeLists.txt~ | 26 + .../DataHelper/DataHelper/CaloHitExtended.h | 72 ++ .../DataHelper/DataHelper/ClusterExtended.h | 105 +++ Utilities/DataHelper/DataHelper/GroupTracks.h | 46 ++ Utilities/DataHelper/DataHelper/HelixClass.h | 302 +++++++ Utilities/DataHelper/DataHelper/LineClass.h | 32 + .../DataHelper/DataHelper/TrackExtended.h | 130 +++ .../DataHelper/TrackerHitExtended.h | 84 ++ Utilities/DataHelper/src/CaloHitExtended.cc | 128 +++ Utilities/DataHelper/src/CaloHitExtended.cc~ | 128 +++ Utilities/DataHelper/src/CaloHitExtended.h~ | 72 ++ Utilities/DataHelper/src/ClusterExtended.cc | 232 ++++++ Utilities/DataHelper/src/ClusterExtended.cc~ | 232 ++++++ Utilities/DataHelper/src/ClusterExtended.h~ | 105 +++ Utilities/DataHelper/src/GroupTracks.cc | 40 + Utilities/DataHelper/src/GroupTracks.cc~ | 40 + Utilities/DataHelper/src/GroupTracks.h~ | 44 + Utilities/DataHelper/src/HelixClass.cc | 768 ++++++++++++++++++ Utilities/DataHelper/src/HelixClass.cc~ | 768 ++++++++++++++++++ Utilities/DataHelper/src/LineClass.cc | 73 ++ Utilities/DataHelper/src/LineClass.cc~ | 73 ++ Utilities/DataHelper/src/TrackExtended.cc | 229 ++++++ Utilities/DataHelper/src/TrackExtended.cc~ | 229 ++++++ Utilities/DataHelper/src/TrackExtended.h~ | 130 +++ .../DataHelper/src/TrackerHitExtended.cc | 133 +++ .../DataHelper/src/TrackerHitExtended.cc~ | 133 +++ .../DataHelper/src/TrackerHitExtended.h~ | 84 ++ 28 files changed, 4452 insertions(+) create mode 100644 Utilities/DataHelper/CMakeLists.txt create mode 100644 Utilities/DataHelper/CMakeLists.txt~ create mode 100644 Utilities/DataHelper/DataHelper/CaloHitExtended.h create mode 100644 Utilities/DataHelper/DataHelper/ClusterExtended.h create mode 100644 Utilities/DataHelper/DataHelper/GroupTracks.h create mode 100644 Utilities/DataHelper/DataHelper/HelixClass.h create mode 100644 Utilities/DataHelper/DataHelper/LineClass.h create mode 100644 Utilities/DataHelper/DataHelper/TrackExtended.h create mode 100644 Utilities/DataHelper/DataHelper/TrackerHitExtended.h create mode 100644 Utilities/DataHelper/src/CaloHitExtended.cc create mode 100644 Utilities/DataHelper/src/CaloHitExtended.cc~ create mode 100644 Utilities/DataHelper/src/CaloHitExtended.h~ create mode 100644 Utilities/DataHelper/src/ClusterExtended.cc create mode 100644 Utilities/DataHelper/src/ClusterExtended.cc~ create mode 100644 Utilities/DataHelper/src/ClusterExtended.h~ create mode 100644 Utilities/DataHelper/src/GroupTracks.cc create mode 100644 Utilities/DataHelper/src/GroupTracks.cc~ create mode 100644 Utilities/DataHelper/src/GroupTracks.h~ create mode 100644 Utilities/DataHelper/src/HelixClass.cc create mode 100644 Utilities/DataHelper/src/HelixClass.cc~ create mode 100644 Utilities/DataHelper/src/LineClass.cc create mode 100644 Utilities/DataHelper/src/LineClass.cc~ create mode 100644 Utilities/DataHelper/src/TrackExtended.cc create mode 100644 Utilities/DataHelper/src/TrackExtended.cc~ create mode 100644 Utilities/DataHelper/src/TrackExtended.h~ create mode 100644 Utilities/DataHelper/src/TrackerHitExtended.cc create mode 100644 Utilities/DataHelper/src/TrackerHitExtended.cc~ create mode 100644 Utilities/DataHelper/src/TrackerHitExtended.h~ diff --git a/Utilities/DataHelper/CMakeLists.txt b/Utilities/DataHelper/CMakeLists.txt new file mode 100644 index 00000000..b8cfb6ee --- /dev/null +++ b/Utilities/DataHelper/CMakeLists.txt @@ -0,0 +1,14 @@ +gaudi_subdir(DataHelper v0r0) + +find_package(EDM4HEP REQUIRED) + +gaudi_depends_on_subdirs() + +set(DataHelperLib_srcs src/*.cc) + +#gaudi_install_headers(DataHelper) + +gaudi_add_library(DataHelperLib ${DataHelperLib_srcs} + PUBLIC_HEADERS DataHelper + LINK_LIBRARIES EDM4HEP::edm4hep EDM4HEP::edm4hepDict +) diff --git a/Utilities/DataHelper/CMakeLists.txt~ b/Utilities/DataHelper/CMakeLists.txt~ new file mode 100644 index 00000000..ccbc04fc --- /dev/null +++ b/Utilities/DataHelper/CMakeLists.txt~ @@ -0,0 +1,26 @@ +gaudi_subdir(DataHelper v0r0) + +#find_package(CLHEP REQUIRED) +#find_package(GEAR REQUIRED) +#find_package(GSL REQUIRED ) +#find_package(LCIO REQUIRED ) +#find_package(podio REQUIRED ) +#find_package(K4FWCore REQUIRED) +find_package(EDM4HEP REQUIRED) + +gaudi_depends_on_subdirs() + +set(DataHelper_srcs src/*.cc) + +gaudi_install_headers(DataHelper) + +#gaudi_add_library(DataHelperLib ${DataHelperLib_srcs} +# PUBLIC_HEADERS DataHelper +# LINK_LIBRARIES EDM4HEP::edm4hep EDM4HEP::edm4hepDict +#) + +# Modules +gaudi_add_module(DataHelper ${DataHelper_srcs} + INCLUDE_DIRS K4FWCore GaudiKernel GaudiAlgLib CLHEP gear ${GSL_INCLUDE_DIRS} ${LCIO_INCLUDE_DIRS} + LINK_LIBRARIES K4FWCore GaudiKernel GaudiAlgLib CLHEP $ENV{GEAR}/lib/libgearsurf.so ${GSL_LIBRARIES} ${LCIO_LIBRARIES} EDM4HEP::edm4hep EDM4HEP::edm4hepDict +) diff --git a/Utilities/DataHelper/DataHelper/CaloHitExtended.h b/Utilities/DataHelper/DataHelper/CaloHitExtended.h new file mode 100644 index 00000000..07b717ad --- /dev/null +++ b/Utilities/DataHelper/DataHelper/CaloHitExtended.h @@ -0,0 +1,72 @@ +#ifndef CaloHitExtended_H +#define CaloHitExtended_H 1 + +//#include "lcio.h" +//#include "EVENT/LCIO.h" +#include "edm4hep/CalorimeterHit.h" +#include "ClusterExtended.h" + +using namespace edm4hep ; + +class ClusterExtended; + +/** + * Class extending native LCIO class CalorimeterHit. <br> + * Class CaloHitExtended is used in TrackwiseClustering <br> + * and Wolf processors. <br> + * @author A. Raspereza (DESY)<br> + * @version $Id: CaloHitExtended.h,v 1.4 2006-02-22 14:53:27 owendt Exp $<br> + */ +class CaloHitExtended { + + public: + + CaloHitExtended(const edm4hep::CalorimeterHit& calhit, int type); + + ~CaloHitExtended(); + + CalorimeterHit * getCalorimeterHit(); + CaloHitExtended * getCaloHitTo(); + CaloHitExtended * getCaloHitFrom(); + ClusterExtended * getClusterExtended(); + int getIndex(); + int getType(); + const float* getDirVec(); + float getYresTo(); + float getYresFrom(); + float getGenericDistance(); + + //void setCalorimeterHit(CalorimeterHit* calhit); + void setCaloHitTo(CaloHitExtended* calhitto); + void setCaloHitFrom(CaloHitExtended* calohitfrom); + void setClusterExtended(ClusterExtended* cluster); + void setIndex(int index); + void setType(int type); + void setDirVec(float* dirVec); + void setYresTo(float yresto); + void setYresFrom(float yresfrom); + void setGenericDistance(float genericDistance); + void setDistanceToCalo(float distanceToCalo); + float getDistanceToCalo(); + void setDistanceToNearestHit(float distanceToNearest); + float getDistanceToNearestHit(); + + + private: + + edm4hep::CalorimeterHit _calohit; + CaloHitExtended * _calohitTo; + CaloHitExtended * _calohitFrom; + ClusterExtended * _clusterAR; + int _index; + int _type; + float _dirVec[3]; + float _yresTo; + float _yresFrom; + float _genericDistance; + float _caloDistance; + float _distanceToNearestHit; + +}; + +#endif diff --git a/Utilities/DataHelper/DataHelper/ClusterExtended.h b/Utilities/DataHelper/DataHelper/ClusterExtended.h new file mode 100644 index 00000000..6b38d887 --- /dev/null +++ b/Utilities/DataHelper/DataHelper/ClusterExtended.h @@ -0,0 +1,105 @@ +#ifndef ClusterExtended_H +#define ClusterExtended_H 1 + +//#include "CaloHitExtended.h" +//#include "TrackExtended.h" +#include "edm4hep/Cluster.h" +#include "HelixClass.h" + +using namespace edm4hep; + +class TrackExtended; +typedef std::vector<TrackExtended*> TrackExtendedVec; + +class CaloHitExtended; +typedef std::vector<CaloHitExtended*> CaloHitExtendedVec; + +//class ClusterExtended; +//typedef std::vector<ClusterExtended*> ClusterExtendedVec; + +/** + * Class extending native LCIO class Cluster. <br> + * Class ClusterExtended is used in TrackwiseClustering <br> + * and Wolf processors. <br> + * @author A. Raspereza (DESY)<br> + * @version $Id: ClusterExtended.h,v 1.4 2006-02-22 14:41:41 owendt Exp $<br> + */ + +class ClusterExtended { + + public: + + ClusterExtended(); + ClusterExtended( Cluster * cluster ); + ClusterExtended(CaloHitExtended * calohit); + ClusterExtended(TrackExtended * track); + + ~ClusterExtended(); + + CaloHitExtendedVec & getCaloHitExtendedVec(); + TrackExtendedVec & getTrackExtendedVec(); + const float* getStartingPoint(); + const float* getDirection(); + void setStartingPoint(float* sPoint); + void setDirection(float* direct); + void addCaloHitExtended(CaloHitExtended * calohit); + void addTrackExtended(TrackExtended * track); + void setType( int type ); + int getType(); + + void Clear(); + + void setCluster(Cluster * cluster); + Cluster * getCluster(); + + void setAxis(float * axis); + float * getAxis(); + + void setEccentricity( float eccentricity); + float getEccentricity(); + + void setHelix(HelixClass helix); + HelixClass & getHelix(); + + void setHelixChi2R(float helixChi2); + float getHelixChi2R(); + + void setHelixChi2Z(float helixChi2); + float getHelixChi2Z(); + + void setPosition(float * position); + float * getPosition(); + + void setLowEdge(float * lowEdge); + float * getLowEdge(); + void setUpEdge(float * upEdge); + float * getUpEdge(); + + + private: + + TrackExtendedVec _trackVector; + CaloHitExtendedVec _hitVector; + float _startingPoint[3]; + float _direction[3]; + + int _type; + Cluster * _cluster; + + float _axis[3]; + float _position[3]; + float _eccentricity; + + HelixClass _helix; + float _helixChi2R; + float _helixChi2Z; + + float _lowEdge[3]; + float _upEdge[3]; + + +}; + +typedef std::vector<ClusterExtended*> ClusterExtendedVec; + +#endif diff --git a/Utilities/DataHelper/DataHelper/GroupTracks.h b/Utilities/DataHelper/DataHelper/GroupTracks.h new file mode 100644 index 00000000..6ae27f1c --- /dev/null +++ b/Utilities/DataHelper/DataHelper/GroupTracks.h @@ -0,0 +1,46 @@ +#ifndef GROUPTRACKS_H +#define GROUPTRACKS_H 1 + +//#include "TrackExtended.h" +//#include "ClusterExtended.h" +#include <vector> + +//fg : forwar declaration needed because of circular include .... +class TrackExtended; +typedef std::vector<TrackExtended*> TrackExtendedVec; +//fg : forward .... + +/** + * Class GroupTracks is needed to group track segments <br> + * with close helix parameters. Needed for Tracking. <br> + * * @author A. Raspereza (DESY)<br> + * @version $Id: GroupTracks.h,v 1.4 2007-09-05 09:39:49 rasp Exp $<br> + */ +//class GroupTracks; + +//typedef std::vector<GroupTracks*> GroupTracksVec; + +class GroupTracks { + + public: + GroupTracks(); + GroupTracks(TrackExtended * track ); + ~GroupTracks(); + + void addTrackExtended( TrackExtended * track ); + void ClearTrackExtendedVec(); + TrackExtendedVec & getTrackExtendedVec(); + void setEdges(float * edges); + float * getEdges(); + + private: + + TrackExtendedVec _trackARVec; + + float _edges[2]; + +}; + +typedef std::vector<GroupTracks*> GroupTracksVec; + +#endif diff --git a/Utilities/DataHelper/DataHelper/HelixClass.h b/Utilities/DataHelper/DataHelper/HelixClass.h new file mode 100644 index 00000000..794f8158 --- /dev/null +++ b/Utilities/DataHelper/DataHelper/HelixClass.h @@ -0,0 +1,302 @@ +#ifndef HELIXAR_H +#define HELIXAR_H 1 +#include <vector> +/** + * Utility class to manipulate with different parameterisations <br> + * of helix. Helix can be initialized in a three different <br> + * ways using the following public methods : <br> + * 1) Initialize_VP(float * pos, float * mom, float q, float B) : <br> + * initialization of helix is done using <br> + * - position of the reference point : pos[3]; <br> + * - momentum vector at the reference point : mom[3];<br> + * - particle charge : q;<br> + * - magnetic field (in Tesla) : B;<br> + * 2) Initialize_BZ(float xCentre, float yCentre, float radius, <br> + * float bZ, float phi0, float B, float signPz,<br> + * float zBegin):<br> + * initialization of helix is done according to the following<br> + * parameterization<br> + * x = xCentre + radius*cos(bZ*z + phi0)<br> + * y = yCentre + radius*sin(bZ*z + phi0)<br> + * where (x,y,z) - position of point on the helix<br> + * - (xCentre,yCentre) is the centre of circumference in R-Phi plane<br> + * - radius is the radius of circumference<br> + * - bZ is the helix slope parameter<br> + * - phi0 is the initial phase of circumference<br> + * - B is the magnetic field (in Tesla)<br> + * - signPz is the sign of the z component of momentum vector<br> + * - zBegin is the z coordinate of the reference point;<br> + * 3) void Initialize_Canonical(float phi0, float d0, float z0, float omega,<br> + * float tanlambda, float B) :<br> + * canonical (LEP-wise) parameterisation with the following parameters<br> + * - phi0 - Phi angle of momentum vector at the point of<br> + * closest approach to IP in R-Phi plane; + * - d0 - signed distance of closest approach to IP in R-Phi plane;<br> + * - z0 - z coordinate of the point of closest approach in R-Phi plane;<br> + * - omega - signed curvature;<br> + * - tanlambda - tangent of dip angle;<br> + * - B - magnetic field (in Tesla);<br> + * A set of public methods (getters) provide access to <br> + * various parameters associated with helix. Helix Class contains<br> + * also several utility methods, allowing for calculation of helix<br> + * intersection points with planes parallel and perpendicular to <br> + * z (beam) axis and determination of the distance of closest approach<br> + * from arbitrary 3D point to the helix. <br> + * @author A. Raspereza (DESY)<br> + * @version $Id: HelixClass.h,v 1.16 2008-06-05 13:47:18 rasp Exp $<br> + * + */ + +#include "LineClass.h" +class HelixClass; + +class HelixClass { + public: + +/** + * Constructor. Initializations of constants which are used + * to calculate various parameters associated with helix. + */ + HelixClass(); +/** + * Destructor. + */ + ~HelixClass(); +/** + * Initialization of helix using <br> + * - position of the reference point : pos[3]; <br> + * - momentum vector at the reference point : mom[3];<br> + * - particle charge : q;<br> + * - magnetic field (in Tesla) : B<br> + */ + void Initialize_VP(float * pos, float * mom, float q, float B); + +/** + * Initialization of helix according to the following<br> + * parameterization<br> + * x = xCentre + radius*cos(bZ*z + phi0)<br> + * y = yCentre + radius*sin(bZ*z + phi0)<br> + * where (x,y,z) - position of point on the helix<br> + * - (xCentre,yCentre) is the centre of circumference in R-Phi plane<br> + * - radius is the radius of circumference<br> + * - bZ is the helix slope parameter<br> + * - phi0 is the initial phase of circumference<br> + * - B is the magnetic field (in Tesla)<br> + * - signPz is the sign of the z component of momentum vector<br> + * - zBegin is the z coordinate of the reference point<br> + */ + void Initialize_BZ(float xCentre, float yCentre, float radius, + float bZ, float phi0, float B, float signPz, + float zBegin); +/** + * Canonical (LEP-wise) parameterisation with the following parameters<br> + * - phi0 - Phi angle of momentum vector at the point of<br> + * closest approach to IP in R-Phi plane; + * - d0 - signed distance of closest approach in R-Phi plane;<br> + * - z0 - z coordinate of the point of closest approach to IP + * in R-Phi plane;<br> + * - omega - signed curvature;<br> + * - tanlambda - tangent of dip angle;<br> + * - B - magnetic field (in Tesla)<br> + */ + void Initialize_Canonical(float phi0, float d0, float z0, float omega, + float tanlambda, float B); + /** + * Returns momentum of particle at the point of closest approach <br> + * to IP <br> + */ + const float * getMomentum(); + + /** + * Returns reference point of track <br> + */ + const float * getReferencePoint(); + + /** + * Returns Phi angle of the momentum vector <br> + * at the point of closest approach to IP <br> + */ + float getPhi0(); + + /** + * Returns signed distance of closest approach <br> + * to IP in the R-Phi plane <br> + */ + float getD0(); + + /** + * Returns z coordinate of the point of closest + * approach to IP in the R-Phi plane <br> + */ + float getZ0(); + + /** + * Returns signed curvature of the track <br> + */ + float getOmega(); + + /** + * Returns tangent of dip angle of the track <br> + */ + float getTanLambda(); + + /** + * Returns transverse momentum of the track <br> + */ + float getPXY(); + + + /** + * Returns x coordinate of circumference + */ + float getXC(); + + /** + * Returns y coordinate of circumference + */ + float getYC(); + + + /** + * Returns radius of circumference + */ + float getRadius(); + + + /** + * Returns helix intersection point with the plane <br> + * parallel to z axis. Plane is defined by two coordinates <br> + * of the point belonging to the plane (x0,y0) and normal <br> + * vector (ax,ay). ref[3] is the reference point of the helix. <br> + * point[3] - returned vector holding the coordinates of <br> + * intersection point <br> + */ + float getPointInXY(float x0, float y0, float ax, float ay, + float * ref , float * point); + + /** + * Returns helix intersection point with the plane <br> + * perpendicular to z axis. Plane is defined by z coordinate <br> + * of the plane. ref[3] is the reference point of the helix. <br> + * point[3] - returned vector holding the coordinates of <br> + * intersection point <br> + */ + float getPointInZ(float zLine, float * ref, float * point); + + /** + * Return distance of the closest approach of the helix to <br> + * arbitrary 3D point in space. xPoint[3] - coordinates of <br> + * space point. Distance[3] - vector of distances of helix to <br> + * a given point in various projections : <br> + * Distance[0] - distance in R-Phi plane <br> + * Distance[1] - distance along Z axis <br> + * Distance[2] - 3D distance <br> + */ + float getDistanceToPoint(float * xPoint, float * Distance); + + /** + * Return distance of the closest approach of the helix to <br> + * arbitrary 3D point in space. xPoint[3] - coordinates of <br> + * space point. distCut - limit on the distance between helix <br> + * and the point to reduce calculation time <br> + * If R-Phi is found to be greater than distCut, rPhi distance is returned <br> + * If the R-Phi distance is not too big, than the exact 3D distance is returned <br> + * This function can be used, if the exact distance is not always needed <br> + */ + float getDistanceToPoint(const float* xPoint, float distCut); + float getDistanceToPoint(const std::vector<float>& xPoint, float distCut); + + /** + * This method calculates coordinates of both intersection <br> + * of the helix with a cylinder. <br> + * Rotational symmetry with respect to z axis is assumed, <br> + * meaning that cylinder axis corresponds to the z axis <br> + * of reference frame. <br> + * Inputs : <br> + * Radius - radius of cylinder. <br> + * ref[3] - point of closest approach to the origin of the helix. <br> + * Output : <br> + * point[6] - coordinates of intersection point. <br> + * Method returns also generic time, defined as the <br> + * ratio of helix length from reference point to the intersection <br> + * point to the particle momentum <br> + */ + float getPointOnCircle(float Radius, float * ref, float * point); + + /** Returns distance between two helixes <br> + * Output : <br> + * pos[3] - position of the point of closest approach <br> + * mom[3] - momentum of V0 <br> + */ + float getDistanceToHelix(HelixClass * helix, float * pos, float * mom); + + /** + * Set Edges of helix + */ + void setHelixEdges(float * xStart, float * xEnd); + + /** + * Returns starting point of helix + */ + float * getStartingPoint() {return _xStart;} + + /** + * Returns endpoint of helix + */ + float * getEndPoint() {return _xEnd;} + + /** + * Returns BZ for the second parameterization + */ + float getBz(); + + /** + * Returns Phi for the second parameterization + */ + float getPhiZ(); + + /** + * Returns extrapolated momentum + */ + void getExtrapolatedMomentum(float * pos, float * momentum); + + /** + * Returns charge + */ + float getCharge(); + + private: + float _momentum[3]; // momentum @ ref point + float _referencePoint[3]; // coordinates @ ref point + float _phi0; // phi0 in canonical parameterization + float _d0; // d0 in canonical parameterisation + float _z0; // z0 in canonical parameterisation + float _omega; // signed curvuture in canonical parameterisation + float _tanLambda; // TanLambda + float _pxy; // Transverse momentum + float _charge; // Particle Charge + float _bField; // Magnetic field (assumed to point to Z>0) + float _radius; // radius of circle in XY plane + float _xCentre; // X of circle centre + float _yCentre; // Y of circle centre + float _phiRefPoint; // Phi w.r.t. (X0,Y0) of circle @ ref point + float _phiAtPCA; // Phi w.r.t. (X0,Y0) of circle @ PCA + float _xAtPCA; // X @ PCA + float _yAtPCA; // Y @ PCA + float _pxAtPCA; // PX @ PCA + float _pyAtPCA; // PY @ PCA + float _phiMomRefPoint; // Phi of Momentum vector @ ref point + float _const_pi; // PI + float _const_2pi; // 2*PI + float _const_pi2; // PI/2 + float _FCT; // 2.99792458E-4 + float _xStart[3]; // Starting point of track segment + float _xEnd[3]; // Ending point of track segment + + float _bZ; + float _phiZ; + +}; + + +#endif diff --git a/Utilities/DataHelper/DataHelper/LineClass.h b/Utilities/DataHelper/DataHelper/LineClass.h new file mode 100644 index 00000000..98970a26 --- /dev/null +++ b/Utilities/DataHelper/DataHelper/LineClass.h @@ -0,0 +1,32 @@ +#ifndef LINECLASS_H +#define LINECLASS_H +class LineClass { + + public: + LineClass(float x0, + float y0, + float z0, + float ax, + float ay, + float az); + + LineClass(float *x0, + float *ax); + + ~LineClass(); + + float * getReferencePoint(); + void setReferencePoint(float *x0); + float * getDirectionalVector(); + void setDirectionalVector(float *ax); + float getDistanceToPoint(float * xpoint, float * pos); + + private: + + float _x0[3]; + float _ax[3]; + + +}; + +#endif diff --git a/Utilities/DataHelper/DataHelper/TrackExtended.h b/Utilities/DataHelper/DataHelper/TrackExtended.h new file mode 100644 index 00000000..61ff5923 --- /dev/null +++ b/Utilities/DataHelper/DataHelper/TrackExtended.h @@ -0,0 +1,130 @@ +#ifndef TRACKAR_H +#define TRACKAR_H 1 + +//#include "lcio.h" +//#include "EVENT/LCIO.h" +#include "edm4hep/Track.h" +#include <vector> +//#include "ClusterExtended.h" +//#include "TrackerHitExtended.h" +#include "GroupTracks.h" + +//using namespace edm4hep; + +class ClusterExtended; +//class GroupTracks; +class TrackerHitExtended; +typedef std::vector<TrackerHitExtended*> TrackerHitExtendedVec; +typedef std::vector<ClusterExtended*> ClusterExtendedVec; + +/** + * Class extending native LCIO class Track. <br> + * Class TrackExtended is used in TrackwiseClustering <br> + * and Wolf processors. <br> + * @author A. Raspereza (DESY)<br> + * @version $Id: TrackExtended.h,v 1.8 2007-09-05 09:39:49 rasp Exp $<br> + */ + +class TrackExtended { + + public: + + + TrackExtended( ); + TrackExtended( TrackerHitExtended * trackerhit ); + TrackExtended( edm4hep::Track * track ); + ~TrackExtended(); + + edm4hep::Track * getTrack(); + const float * getSeedDirection(); + const float * getSeedPosition(); + ClusterExtendedVec & getClusterVec(); + ClusterExtended * getSuperCluster(); + TrackerHitExtendedVec & getTrackerHitExtendedVec(); + void addCluster(ClusterExtended * cluster); + void setSuperCluster(ClusterExtended * superCluster); + void setSeedDirection( float * seedDirection ); + void setSeedPosition( float * seedPosition); + void addTrackerHitExtended( TrackerHitExtended * trackerhit); + void ClearTrackerHitExtendedVec(); + + void setX0(float x0); + void setY0(float y0); + void setR0(float r0); + void setD0(float d0); + void setZ0(float z0); + void setBz(float bz); + void setPhi0(float phi0); + void setPhi(float phi); + void setOmega(float omega); + void setTanLambda(float tanLambda); + + + void setStart(float * xStart); + void setEnd(float * xEnd); + + + float getX0(); + float getY0(); + float getD0(); + float getZ0(); + float getOmega(); + float getTanLambda(); + float getPhi(); + float getR0(); + float getBz(); + float getPhi0(); + + float * getStart(); + float * getEnd(); + + void setGroupTracks( GroupTracks * group ); + GroupTracks * getGroupTracks(); + + float getChi2(); + void setChi2(float chi2); + + int getNDF(); + void setNDF(int ndf); + + float * getCovMatrix(); + void setCovMatrix( float * cov); + + private: + + ClusterExtended *_superCluster; + ClusterExtendedVec _clusterVec; + GroupTracks * _group; + edm4hep::Track * _track; + float _seedDirection[3]; + float _seedPosition[3]; + TrackerHitExtendedVec _trackerHitVector; + + float _x0; + float _y0; + float _r0; + float _bz; + float _phi0; + + float _xStart[3]; + float _xEnd[3]; + + float _d0; // d0 in canonical parameterisation + float _z0; // z0 in canonical parameterisation + float _omega; // omega in canonical parameterisation + float _tanLambda; // tanlambda in canonical parameterisation + float _phi; // phi in canonical parameterisation + + float _chi2; // chi2 of the fit + + float _cov[15]; // covariance matrix + + int _ndf; // NDF + + int _type; // Track type; + +}; + +typedef std::vector<TrackExtended*> TrackExtendedVec; + +#endif diff --git a/Utilities/DataHelper/DataHelper/TrackerHitExtended.h b/Utilities/DataHelper/DataHelper/TrackerHitExtended.h new file mode 100644 index 00000000..d625679e --- /dev/null +++ b/Utilities/DataHelper/DataHelper/TrackerHitExtended.h @@ -0,0 +1,84 @@ +#ifndef TRACKERHITAR_H +#define TRACKERHITAR_H 1 + +//#include "lcio.h" +//#include "EVENT/LCIO.h" +#include "edm4hep/TrackerHit.h" +//#include "TrackExtended.h" +#include <vector> + + +//using namespace edm4hep; + +class TrackExtended; +typedef std::vector<TrackExtended*> TrackExtendedVec; + +/** + * Class extending native LCIO class TrackerHit. <br> + * Class TrackerHitExtended is used in TrackwiseClustering <br> + * and Wolf processors. <br> + * @author A. Raspereza (DESY)<br> + * @version $Id: TrackerHitExtended.h,v 1.7 2007-09-05 09:39:49 rasp Exp $<br> + */ +class TrackerHitExtended { + + public: + + TrackerHitExtended(const edm4hep::TrackerHit& trackerhit); + ~TrackerHitExtended(); + void setTrackExtended(TrackExtended * trackAR); + void addTrackExtended(TrackExtended * trackAR); + void setTrackerHitTo(TrackerHitExtended * hitTo); + void setTrackerHitFrom(TrackerHitExtended * hitFrom); + void setGenericDistance(float genericDistance); + //void setTrackerHit(edm4hep::TrackerHit* hit); + void setYresTo(float yresTo); + void setYresFrom(float yresFrom); + void setDirVec(float * dirVec); + void clearTrackVec(); + void setResolutionRPhi(float rphiReso); + void setResolutionZ(float zReso); + void setType(int type); + void setDet(int idet); + void setUsedInFit(bool usedInFit); + + edm4hep::TrackerHit* getTrackerHit(); + TrackExtended * getTrackExtended(); + TrackExtendedVec & getTrackExtendedVec(); + TrackerHitExtended * getTrackerHitFrom(); + TrackerHitExtended * getTrackerHitTo(); + float getGenericDistance(); + float getYresTo(); + float getYresFrom(); + float * getDirVec(); + float getResolutionRPhi(); + float getResolutionZ(); + int getType(); + int getDet(); + bool getUsedInFit(); + + private: + + edm4hep::TrackerHit _trackerHit; + TrackExtended * _trackAR; + TrackerHitExtended * _hitTo; + TrackerHitExtended * _hitFrom; + TrackExtendedVec _trackVecAR; + + float _rphiReso; + float _zReso; + float _yresTo; + float _yresFrom; + float _genericDistance; + float _dirVec[3]; + + int _type; + int _idet; + + bool _usedInFit; + +}; + +typedef std::vector<TrackerHitExtended*> TrackerHitExtendedVec; + +#endif diff --git a/Utilities/DataHelper/src/CaloHitExtended.cc b/Utilities/DataHelper/src/CaloHitExtended.cc new file mode 100644 index 00000000..8c608d97 --- /dev/null +++ b/Utilities/DataHelper/src/CaloHitExtended.cc @@ -0,0 +1,128 @@ +#include "DataHelper/CaloHitExtended.h" +#include <math.h> + +CaloHitExtended::CaloHitExtended(const edm4hep::CalorimeterHit& calhit, int type) + :_calohit(calhit){ + _type = type; + _index = 0; + _calohitTo = NULL; + _calohitFrom = NULL; +} + +CaloHitExtended::~CaloHitExtended(){ +} + +CalorimeterHit* CaloHitExtended::getCalorimeterHit() { + return &_calohit; +} + +CaloHitExtended* CaloHitExtended::getCaloHitTo() { + return _calohitTo; +} + +CaloHitExtended* CaloHitExtended::getCaloHitFrom() { + return _calohitFrom; +} + +ClusterExtended* CaloHitExtended::getClusterExtended() { + return _clusterAR; +} + +int CaloHitExtended::getIndex() { + return _index; + +} + +int CaloHitExtended::getType() { + return _type; + +} +const float* CaloHitExtended::getDirVec() { + return _dirVec; +} + +float CaloHitExtended::getYresTo() { + return _yresTo; +} + +float CaloHitExtended::getYresFrom() { + return _yresFrom; +} + + + +float CaloHitExtended::getGenericDistance() { + return _genericDistance; +} + +//void CaloHitExtended::setCalorimeterHit(CalorimeterHit* calhit) { +// _calohit = calhit; +//} + +void CaloHitExtended::setCaloHitTo(CaloHitExtended* calhitto) { + _calohitTo = calhitto; +} + +void CaloHitExtended::setCaloHitFrom(CaloHitExtended* calhitfrom) { + _calohitFrom = calhitfrom; +} + +void CaloHitExtended::setClusterExtended(ClusterExtended* cluster) { + _clusterAR = cluster; +} + +void CaloHitExtended::setIndex(int index) { + _index = index; +} + +void CaloHitExtended::setType(int type) { + _type = type; +} + +void CaloHitExtended::setDirVec(float* dirVec) { + float modulus(0.); + for (int i(0); i < 3; ++i) { + modulus += dirVec[i]*dirVec[i]; + } + + modulus = sqrt(modulus); + + if (modulus <= 0.) { + _dirVec[0] = 0.; + _dirVec[1] = 0.; + _dirVec[2] = 1.; + } + else { + _dirVec[0] = dirVec[0] / modulus; + _dirVec[1] = dirVec[1] / modulus; + _dirVec[2] = dirVec[2] / modulus; + } +} + +void CaloHitExtended::setYresTo(float yresto) { + _yresTo = yresto; +} + +void CaloHitExtended::setYresFrom(float yresfrom) { + _yresFrom = yresfrom; +} + +void CaloHitExtended::setGenericDistance(float genericDistance) { + _genericDistance = genericDistance; +} + +void CaloHitExtended::setDistanceToCalo(float caloDistance) { + _caloDistance = caloDistance; +} + +float CaloHitExtended::getDistanceToCalo() { + return _caloDistance; +} + +void CaloHitExtended::setDistanceToNearestHit(float distanceToNearest) { + _distanceToNearestHit = distanceToNearest; +} + +float CaloHitExtended::getDistanceToNearestHit() { + return _distanceToNearestHit; +} diff --git a/Utilities/DataHelper/src/CaloHitExtended.cc~ b/Utilities/DataHelper/src/CaloHitExtended.cc~ new file mode 100644 index 00000000..3adab51d --- /dev/null +++ b/Utilities/DataHelper/src/CaloHitExtended.cc~ @@ -0,0 +1,128 @@ +#include "CaloHitExtended.h" +#include <math.h> + +CaloHitExtended::CaloHitExtended(const edm4hep::CalorimeterHit& calhit, int type) + :_calohit(calhit){ + _type = type; + _index = 0; + _calohitTo = NULL; + _calohitFrom = NULL; +} + +CaloHitExtended::~CaloHitExtended(){ +} + +CalorimeterHit* CaloHitExtended::getCalorimeterHit() { + return &_calohit; +} + +CaloHitExtended* CaloHitExtended::getCaloHitTo() { + return _calohitTo; +} + +CaloHitExtended* CaloHitExtended::getCaloHitFrom() { + return _calohitFrom; +} + +ClusterExtended* CaloHitExtended::getClusterExtended() { + return _clusterAR; +} + +int CaloHitExtended::getIndex() { + return _index; + +} + +int CaloHitExtended::getType() { + return _type; + +} +const float* CaloHitExtended::getDirVec() { + return _dirVec; +} + +float CaloHitExtended::getYresTo() { + return _yresTo; +} + +float CaloHitExtended::getYresFrom() { + return _yresFrom; +} + + + +float CaloHitExtended::getGenericDistance() { + return _genericDistance; +} + +//void CaloHitExtended::setCalorimeterHit(CalorimeterHit* calhit) { +// _calohit = calhit; +//} + +void CaloHitExtended::setCaloHitTo(CaloHitExtended* calhitto) { + _calohitTo = calhitto; +} + +void CaloHitExtended::setCaloHitFrom(CaloHitExtended* calhitfrom) { + _calohitFrom = calhitfrom; +} + +void CaloHitExtended::setClusterExtended(ClusterExtended* cluster) { + _clusterAR = cluster; +} + +void CaloHitExtended::setIndex(int index) { + _index = index; +} + +void CaloHitExtended::setType(int type) { + _type = type; +} + +void CaloHitExtended::setDirVec(float* dirVec) { + float modulus(0.); + for (int i(0); i < 3; ++i) { + modulus += dirVec[i]*dirVec[i]; + } + + modulus = sqrt(modulus); + + if (modulus <= 0.) { + _dirVec[0] = 0.; + _dirVec[1] = 0.; + _dirVec[2] = 1.; + } + else { + _dirVec[0] = dirVec[0] / modulus; + _dirVec[1] = dirVec[1] / modulus; + _dirVec[2] = dirVec[2] / modulus; + } +} + +void CaloHitExtended::setYresTo(float yresto) { + _yresTo = yresto; +} + +void CaloHitExtended::setYresFrom(float yresfrom) { + _yresFrom = yresfrom; +} + +void CaloHitExtended::setGenericDistance(float genericDistance) { + _genericDistance = genericDistance; +} + +void CaloHitExtended::setDistanceToCalo(float caloDistance) { + _caloDistance = caloDistance; +} + +float CaloHitExtended::getDistanceToCalo() { + return _caloDistance; +} + +void CaloHitExtended::setDistanceToNearestHit(float distanceToNearest) { + _distanceToNearestHit = distanceToNearest; +} + +float CaloHitExtended::getDistanceToNearestHit() { + return _distanceToNearestHit; +} diff --git a/Utilities/DataHelper/src/CaloHitExtended.h~ b/Utilities/DataHelper/src/CaloHitExtended.h~ new file mode 100644 index 00000000..65781a9d --- /dev/null +++ b/Utilities/DataHelper/src/CaloHitExtended.h~ @@ -0,0 +1,72 @@ +#ifndef CaloHitExtended_H +#define CaloHitExtended_H 1 + +//#include "lcio.h" +//#include "EVENT/LCIO.h" +#include "edm4hep/CalorimeterHit.h" +#include "ClusterExtended.h" + +using namespace edm4hep ; + +class ClusterExtended; + +/** + * Class extending native LCIO class CalorimeterHit. <br> + * Class CaloHitExtended is used in TrackwiseClustering <br> + * and Wolf processors. <br> + * @author A. Raspereza (DESY)<br> + * @version $Id: CaloHitExtended.h,v 1.4 2006-02-22 14:53:27 owendt Exp $<br> + */ +class CaloHitExtended { + + public: + + CaloHitExtended(CalorimeterHit* calhit, int type); + + ~CaloHitExtended(); + + CalorimeterHit * getCalorimeterHit(); + CaloHitExtended * getCaloHitTo(); + CaloHitExtended * getCaloHitFrom(); + ClusterExtended * getClusterExtended(); + int getIndex(); + int getType(); + const float* getDirVec(); + float getYresTo(); + float getYresFrom(); + float getGenericDistance(); + + void setCalorimeterHit(CalorimeterHit* calhit); + void setCaloHitTo(CaloHitExtended* calhitto); + void setCaloHitFrom(CaloHitExtended* calohitfrom); + void setClusterExtended(ClusterExtended* cluster); + void setIndex(int index); + void setType(int type); + void setDirVec(float* dirVec); + void setYresTo(float yresto); + void setYresFrom(float yresfrom); + void setGenericDistance(float genericDistance); + void setDistanceToCalo(float distanceToCalo); + float getDistanceToCalo(); + void setDistanceToNearestHit(float distanceToNearest); + float getDistanceToNearestHit(); + + + private: + + CalorimeterHit * _calohit; + CaloHitExtended * _calohitTo; + CaloHitExtended * _calohitFrom; + ClusterExtended * _clusterAR; + int _index; + int _type; + float _dirVec[3]; + float _yresTo; + float _yresFrom; + float _genericDistance; + float _caloDistance; + float _distanceToNearestHit; + +}; + +#endif diff --git a/Utilities/DataHelper/src/ClusterExtended.cc b/Utilities/DataHelper/src/ClusterExtended.cc new file mode 100644 index 00000000..beee779d --- /dev/null +++ b/Utilities/DataHelper/src/ClusterExtended.cc @@ -0,0 +1,232 @@ +#include "DataHelper/CaloHitExtended.h" +#include "DataHelper/TrackExtended.h" +#include "DataHelper/ClusterExtended.h" +#include <math.h> + +ClusterExtended::ClusterExtended() { + _hitVector.clear(); + _trackVector.clear(); + + _direction[0] = 0.; + _direction[1] = 0.; + _direction[2] = 0.; + + _startingPoint[0] = 0.; + _startingPoint[1] = 0.; + _startingPoint[2] = 0.; + +} + +ClusterExtended::ClusterExtended( Cluster * cluster ) { + _hitVector.clear(); + _trackVector.clear(); + + _direction[0] = 0.; + _direction[1] = 0.; + _direction[2] = 0.; + + _startingPoint[0] = 0.; + _startingPoint[1] = 0.; + _startingPoint[2] = 0.; + + _cluster = cluster; + +} + + + + +ClusterExtended::ClusterExtended(CaloHitExtended * calohit) { + + _hitVector.clear(); + _hitVector.push_back(calohit); + _trackVector.clear(); + + float rad(0); + + _startingPoint[0] = calohit->getCalorimeterHit()->getPosition().x; + _startingPoint[1] = calohit->getCalorimeterHit()->getPosition().y; + _startingPoint[2] = calohit->getCalorimeterHit()->getPosition().z; + + for (int i(0); i < 3; ++i) { + rad += _startingPoint[i]*_startingPoint[i]; + } + + rad = sqrt(rad); + + _direction[0] = _startingPoint[0]/rad; + _direction[1] = _startingPoint[1]/rad; + _direction[2] = _startingPoint[2]/rad; + +} + +ClusterExtended::ClusterExtended(TrackExtended * track) { + _hitVector.clear(); + _trackVector.clear(); + _trackVector.push_back(track); + for (int i(0); i < 3; ++i) { + _startingPoint[i] = track->getSeedPosition()[i]; + _direction[i] = track->getSeedDirection()[i]; + } +} + + +ClusterExtended::~ClusterExtended() { + _hitVector.clear(); + _trackVector.clear(); +} + +CaloHitExtendedVec& ClusterExtended::getCaloHitExtendedVec() { + return _hitVector; +} + +TrackExtendedVec& ClusterExtended::getTrackExtendedVec() { + return _trackVector; +} + +const float* ClusterExtended::getStartingPoint() { + return _startingPoint; +} + +const float* ClusterExtended::getDirection() { + return _direction; +} + +void ClusterExtended::setStartingPoint(float* sPoint) { + _startingPoint[0] = sPoint[0]; + _startingPoint[1] = sPoint[1]; + _startingPoint[2] = sPoint[2]; + +} + +void ClusterExtended::setDirection(float* direct) { + _direction[0] = direct[0]; + _direction[1] = direct[1]; + _direction[2] = direct[2]; +} + +void ClusterExtended::addCaloHitExtended(CaloHitExtended* calohit) { + _hitVector.push_back(calohit); +} + +void ClusterExtended::addTrackExtended(TrackExtended * track) { + _trackVector.push_back(track); +} + +void ClusterExtended::Clear() { + _hitVector.clear(); + _trackVector.clear(); + +} + +void ClusterExtended::setType( int type ) { + _type = type; +} + +int ClusterExtended::getType() { + return _type; +} + +void ClusterExtended::setCluster(Cluster * cluster) { + _cluster = cluster; +} + +Cluster * ClusterExtended::getCluster() { + return _cluster; +} + +void ClusterExtended::setAxis(float * axis) { + _axis[0] = axis[0]; + _axis[1] = axis[1]; + _axis[2] = axis[2]; +} +float * ClusterExtended::getAxis() { + return _axis; +} + +void ClusterExtended::setEccentricity( float eccentricity) { + _eccentricity = eccentricity; +} +float ClusterExtended::getEccentricity() { + return _eccentricity; +} + +void ClusterExtended::setHelix(HelixClass helix) { + _helix = helix; + int nHits = int(_hitVector.size()); + float timeMax = -1.0e+20; + float timeMin = 1.0e+20; + for (int ihit=0;ihit<nHits;++ihit) { + float pos[3]; + for (int i=0;i<3;++i) + pos[i]=_hitVector[ihit]->getCalorimeterHit()->getPosition()[i]; + + float distance[3]; + float time = _helix.getDistanceToPoint(pos,distance); + if (time > timeMax) { + timeMax = time; + _upEdge[0] = pos[0]; + _upEdge[1] = pos[1]; + _upEdge[2] = pos[2]; + } + if (time < timeMin) { + timeMin = time; + _lowEdge[0] = pos[0]; + _lowEdge[1] = pos[1]; + _lowEdge[2] = pos[2]; + } + + } + +} +HelixClass & ClusterExtended::getHelix() { + return _helix; +} + +void ClusterExtended::setHelixChi2R(float helixChi2) { + _helixChi2R = helixChi2; +} +float ClusterExtended::getHelixChi2R() { + return _helixChi2R; +} + +void ClusterExtended::setHelixChi2Z(float helixChi2) { + _helixChi2Z = helixChi2; +} +float ClusterExtended::getHelixChi2Z() { + return _helixChi2Z; +} + + + +void ClusterExtended::setPosition(float * position) { + _position[0] = position[0]; + _position[1] = position[1]; + _position[2] = position[2]; + +} + +float * ClusterExtended::getPosition() { + return _position; +} + +void ClusterExtended::setUpEdge(float * upEdge) { + _upEdge[0] = upEdge[0]; + _upEdge[1] = upEdge[1]; + _upEdge[2] = upEdge[2]; +} + +void ClusterExtended::setLowEdge(float * lowEdge) { + _lowEdge[0] = lowEdge[0]; + _lowEdge[1] = lowEdge[1]; + _lowEdge[2] = lowEdge[2]; +} + +float * ClusterExtended::getUpEdge() { + return _upEdge; +} + +float * ClusterExtended::getLowEdge() { + return _lowEdge; +} + diff --git a/Utilities/DataHelper/src/ClusterExtended.cc~ b/Utilities/DataHelper/src/ClusterExtended.cc~ new file mode 100644 index 00000000..a5f99aeb --- /dev/null +++ b/Utilities/DataHelper/src/ClusterExtended.cc~ @@ -0,0 +1,232 @@ +#include "CaloHitExtended.h" +#include "TrackExtended.h" +#include "ClusterExtended.h" +#include <math.h> + +ClusterExtended::ClusterExtended() { + _hitVector.clear(); + _trackVector.clear(); + + _direction[0] = 0.; + _direction[1] = 0.; + _direction[2] = 0.; + + _startingPoint[0] = 0.; + _startingPoint[1] = 0.; + _startingPoint[2] = 0.; + +} + +ClusterExtended::ClusterExtended( Cluster * cluster ) { + _hitVector.clear(); + _trackVector.clear(); + + _direction[0] = 0.; + _direction[1] = 0.; + _direction[2] = 0.; + + _startingPoint[0] = 0.; + _startingPoint[1] = 0.; + _startingPoint[2] = 0.; + + _cluster = cluster; + +} + + + + +ClusterExtended::ClusterExtended(CaloHitExtended * calohit) { + + _hitVector.clear(); + _hitVector.push_back(calohit); + _trackVector.clear(); + + float rad(0); + + _startingPoint[0] = calohit->getCalorimeterHit()->getPosition().x; + _startingPoint[1] = calohit->getCalorimeterHit()->getPosition().y; + _startingPoint[2] = calohit->getCalorimeterHit()->getPosition().z; + + for (int i(0); i < 3; ++i) { + rad += _startingPoint[i]*_startingPoint[i]; + } + + rad = sqrt(rad); + + _direction[0] = _startingPoint[0]/rad; + _direction[1] = _startingPoint[1]/rad; + _direction[2] = _startingPoint[2]/rad; + +} + +ClusterExtended::ClusterExtended(TrackExtended * track) { + _hitVector.clear(); + _trackVector.clear(); + _trackVector.push_back(track); + for (int i(0); i < 3; ++i) { + _startingPoint[i] = track->getSeedPosition()[i]; + _direction[i] = track->getSeedDirection()[i]; + } +} + + +ClusterExtended::~ClusterExtended() { + _hitVector.clear(); + _trackVector.clear(); +} + +CaloHitExtendedVec& ClusterExtended::getCaloHitExtendedVec() { + return _hitVector; +} + +TrackExtendedVec& ClusterExtended::getTrackExtendedVec() { + return _trackVector; +} + +const float* ClusterExtended::getStartingPoint() { + return _startingPoint; +} + +const float* ClusterExtended::getDirection() { + return _direction; +} + +void ClusterExtended::setStartingPoint(float* sPoint) { + _startingPoint[0] = sPoint[0]; + _startingPoint[1] = sPoint[1]; + _startingPoint[2] = sPoint[2]; + +} + +void ClusterExtended::setDirection(float* direct) { + _direction[0] = direct[0]; + _direction[1] = direct[1]; + _direction[2] = direct[2]; +} + +void ClusterExtended::addCaloHitExtended(CaloHitExtended* calohit) { + _hitVector.push_back(calohit); +} + +void ClusterExtended::addTrackExtended(TrackExtended * track) { + _trackVector.push_back(track); +} + +void ClusterExtended::Clear() { + _hitVector.clear(); + _trackVector.clear(); + +} + +void ClusterExtended::setType( int type ) { + _type = type; +} + +int ClusterExtended::getType() { + return _type; +} + +void ClusterExtended::setCluster(Cluster * cluster) { + _cluster = cluster; +} + +Cluster * ClusterExtended::getCluster() { + return _cluster; +} + +void ClusterExtended::setAxis(float * axis) { + _axis[0] = axis[0]; + _axis[1] = axis[1]; + _axis[2] = axis[2]; +} +float * ClusterExtended::getAxis() { + return _axis; +} + +void ClusterExtended::setEccentricity( float eccentricity) { + _eccentricity = eccentricity; +} +float ClusterExtended::getEccentricity() { + return _eccentricity; +} + +void ClusterExtended::setHelix(HelixClass helix) { + _helix = helix; + int nHits = int(_hitVector.size()); + float timeMax = -1.0e+20; + float timeMin = 1.0e+20; + for (int ihit=0;ihit<nHits;++ihit) { + float pos[3]; + for (int i=0;i<3;++i) + pos[i]=_hitVector[ihit]->getCalorimeterHit()->getPosition()[i]; + + float distance[3]; + float time = _helix.getDistanceToPoint(pos,distance); + if (time > timeMax) { + timeMax = time; + _upEdge[0] = pos[0]; + _upEdge[1] = pos[1]; + _upEdge[2] = pos[2]; + } + if (time < timeMin) { + timeMin = time; + _lowEdge[0] = pos[0]; + _lowEdge[1] = pos[1]; + _lowEdge[2] = pos[2]; + } + + } + +} +HelixClass & ClusterExtended::getHelix() { + return _helix; +} + +void ClusterExtended::setHelixChi2R(float helixChi2) { + _helixChi2R = helixChi2; +} +float ClusterExtended::getHelixChi2R() { + return _helixChi2R; +} + +void ClusterExtended::setHelixChi2Z(float helixChi2) { + _helixChi2Z = helixChi2; +} +float ClusterExtended::getHelixChi2Z() { + return _helixChi2Z; +} + + + +void ClusterExtended::setPosition(float * position) { + _position[0] = position[0]; + _position[1] = position[1]; + _position[2] = position[2]; + +} + +float * ClusterExtended::getPosition() { + return _position; +} + +void ClusterExtended::setUpEdge(float * upEdge) { + _upEdge[0] = upEdge[0]; + _upEdge[1] = upEdge[1]; + _upEdge[2] = upEdge[2]; +} + +void ClusterExtended::setLowEdge(float * lowEdge) { + _lowEdge[0] = lowEdge[0]; + _lowEdge[1] = lowEdge[1]; + _lowEdge[2] = lowEdge[2]; +} + +float * ClusterExtended::getUpEdge() { + return _upEdge; +} + +float * ClusterExtended::getLowEdge() { + return _lowEdge; +} + diff --git a/Utilities/DataHelper/src/ClusterExtended.h~ b/Utilities/DataHelper/src/ClusterExtended.h~ new file mode 100644 index 00000000..4a180240 --- /dev/null +++ b/Utilities/DataHelper/src/ClusterExtended.h~ @@ -0,0 +1,105 @@ +#ifndef ClusterExtended_H +#define ClusterExtended_H 1 + +#include "CaloHitExtended.h" +#include "TrackExtended.h" +#include "plcio/Cluster.h" +#include "HelixClass.h" + +using namespace plcio; + +//class TrackExtended; +//typedef std::vector<TrackExtended*> TrackExtendedVec; + +//class CaloHitExtended; +//typedef std::vector<CaloHitExtended*> CaloHitExtendedVec; + +//class ClusterExtended; +//typedef std::vector<ClusterExtended*> ClusterExtendedVec; + +/** + * Class extending native LCIO class Cluster. <br> + * Class ClusterExtended is used in TrackwiseClustering <br> + * and Wolf processors. <br> + * @author A. Raspereza (DESY)<br> + * @version $Id: ClusterExtended.h,v 1.4 2006-02-22 14:41:41 owendt Exp $<br> + */ + +class ClusterExtended { + + public: + + ClusterExtended(); + ClusterExtended( Cluster * cluster ); + ClusterExtended(CaloHitExtended * calohit); + ClusterExtended(TrackExtended * track); + + ~ClusterExtended(); + + CaloHitExtendedVec & getCaloHitExtendedVec(); + TrackExtendedVec & getTrackExtendedVec(); + const float* getStartingPoint(); + const float* getDirection(); + void setStartingPoint(float* sPoint); + void setDirection(float* direct); + void addCaloHitExtended(CaloHitExtended * calohit); + void addTrackExtended(TrackExtended * track); + void setType( int type ); + int getType(); + + void Clear(); + + void setCluster(Cluster * cluster); + Cluster * getCluster(); + + void setAxis(float * axis); + float * getAxis(); + + void setEccentricity( float eccentricity); + float getEccentricity(); + + void setHelix(HelixClass helix); + HelixClass & getHelix(); + + void setHelixChi2R(float helixChi2); + float getHelixChi2R(); + + void setHelixChi2Z(float helixChi2); + float getHelixChi2Z(); + + void setPosition(float * position); + float * getPosition(); + + void setLowEdge(float * lowEdge); + float * getLowEdge(); + void setUpEdge(float * upEdge); + float * getUpEdge(); + + + private: + + TrackExtendedVec _trackVector; + CaloHitExtendedVec _hitVector; + float _startingPoint[3]; + float _direction[3]; + + int _type; + Cluster * _cluster; + + float _axis[3]; + float _position[3]; + float _eccentricity; + + HelixClass _helix; + float _helixChi2R; + float _helixChi2Z; + + float _lowEdge[3]; + float _upEdge[3]; + + +}; + +typedef std::vector<ClusterExtended*> ClusterExtendedVec; + +#endif diff --git a/Utilities/DataHelper/src/GroupTracks.cc b/Utilities/DataHelper/src/GroupTracks.cc new file mode 100644 index 00000000..593fe02c --- /dev/null +++ b/Utilities/DataHelper/src/GroupTracks.cc @@ -0,0 +1,40 @@ +#include "DataHelper/TrackExtended.h" +#include "DataHelper/GroupTracks.h" + +GroupTracks::GroupTracks() { + _trackARVec.clear(); + _edges[0] = 0.0; + _edges[1] = 0.0; +} + +GroupTracks::GroupTracks( TrackExtended * track ) { + _trackARVec.clear(); + _trackARVec.push_back( track ); + _edges[0] = 0.0; + _edges[1] = 0.0; +} + +GroupTracks::~GroupTracks() {} + +void GroupTracks::addTrackExtended( TrackExtended * track ) { + _trackARVec.push_back( track ); +} + +void GroupTracks::ClearTrackExtendedVec() { + _trackARVec.clear(); +} + +TrackExtendedVec & GroupTracks::getTrackExtendedVec() { + return _trackARVec; +} + +void GroupTracks::setEdges(float * edges) { + + _edges[0] = edges[0]; + _edges[1] = edges[1]; + +} + +float * GroupTracks::getEdges() { + return _edges; +} diff --git a/Utilities/DataHelper/src/GroupTracks.cc~ b/Utilities/DataHelper/src/GroupTracks.cc~ new file mode 100644 index 00000000..fdbda40d --- /dev/null +++ b/Utilities/DataHelper/src/GroupTracks.cc~ @@ -0,0 +1,40 @@ +#include "TrackExtended.h" +#include "GroupTracks.h" + +GroupTracks::GroupTracks() { + _trackARVec.clear(); + _edges[0] = 0.0; + _edges[1] = 0.0; +} + +GroupTracks::GroupTracks( TrackExtended * track ) { + _trackARVec.clear(); + _trackARVec.push_back( track ); + _edges[0] = 0.0; + _edges[1] = 0.0; +} + +GroupTracks::~GroupTracks() {} + +void GroupTracks::addTrackExtended( TrackExtended * track ) { + _trackARVec.push_back( track ); +} + +void GroupTracks::ClearTrackExtendedVec() { + _trackARVec.clear(); +} + +TrackExtendedVec & GroupTracks::getTrackExtendedVec() { + return _trackARVec; +} + +void GroupTracks::setEdges(float * edges) { + + _edges[0] = edges[0]; + _edges[1] = edges[1]; + +} + +float * GroupTracks::getEdges() { + return _edges; +} diff --git a/Utilities/DataHelper/src/GroupTracks.h~ b/Utilities/DataHelper/src/GroupTracks.h~ new file mode 100644 index 00000000..681117bb --- /dev/null +++ b/Utilities/DataHelper/src/GroupTracks.h~ @@ -0,0 +1,44 @@ +#ifndef GROUPTRACKS_H +#define GROUPTRACKS_H 1 + +#include "TrackExtended.h" +#include "ClusterExtended.h" +#include <vector> + +//fg : forwar declaration needed because of circular include .... +class TrackExtended; +typedef std::vector<TrackExtended*> TrackExtendedVec; +//fg : forward .... + +/** + * Class GroupTracks is needed to group track segments <br> + * with close helix parameters. Needed for Tracking. <br> + * * @author A. Raspereza (DESY)<br> + * @version $Id: GroupTracks.h,v 1.4 2007-09-05 09:39:49 rasp Exp $<br> + */ +class GroupTracks; + +typedef std::vector<GroupTracks*> GroupTracksVec; + +class GroupTracks { + + public: + GroupTracks(); + GroupTracks(TrackExtended * track ); + ~GroupTracks(); + + void addTrackExtended( TrackExtended * track ); + void ClearTrackExtendedVec(); + TrackExtendedVec & getTrackExtendedVec(); + void setEdges(float * edges); + float * getEdges(); + + private: + + TrackExtendedVec _trackARVec; + + float _edges[2]; + +}; + +#endif diff --git a/Utilities/DataHelper/src/HelixClass.cc b/Utilities/DataHelper/src/HelixClass.cc new file mode 100644 index 00000000..591acc51 --- /dev/null +++ b/Utilities/DataHelper/src/HelixClass.cc @@ -0,0 +1,768 @@ +#include "DataHelper/HelixClass.h" +#include <math.h> +#include <stdlib.h> +#include <iostream> +//#include "ced_cli.h" + +HelixClass::HelixClass() { + _const_2pi = 2.0*M_PI; + _const_pi2 = 0.5*M_PI; + _FCT = 2.99792458E-4; +} + +HelixClass::~HelixClass() {} + +void HelixClass::Initialize_VP(float * pos, float * mom, float q, float B) { + _referencePoint[0] = pos[0]; + _referencePoint[1] = pos[1]; + _referencePoint[2] = pos[2]; + _momentum[0] = mom[0]; + _momentum[1] = mom[1]; + _momentum[2] = mom[2]; + _charge = q; + _bField = B; + _pxy = sqrt(mom[0]*mom[0]+mom[1]*mom[1]); + _radius = _pxy / (_FCT*B); + _omega = q/_radius; + _tanLambda = mom[2]/_pxy; + _phiMomRefPoint = atan2(mom[1],mom[0]); + _xCentre = pos[0] + _radius*cos(_phiMomRefPoint-_const_pi2*q); + _yCentre = pos[1] + _radius*sin(_phiMomRefPoint-_const_pi2*q); + _phiRefPoint = atan2(pos[1]-_yCentre,pos[0]-_xCentre); + _phiAtPCA = atan2(-_yCentre,-_xCentre); + _phi0 = -_const_pi2*q + _phiAtPCA; + while (_phi0<0) _phi0+=_const_2pi; + while (_phi0>=_const_2pi) _phi0-=_const_2pi; + _xAtPCA = _xCentre + _radius*cos(_phiAtPCA); + _yAtPCA = _yCentre + _radius*sin(_phiAtPCA); + // _d0 = -_xAtPCA*sin(_phi0) + _yAtPCA*cos(_phi0); + double pxy = double(_pxy); + double radius = pxy/double(_FCT*B); + double xCentre = double(pos[0]) + radius*double(cos(_phiMomRefPoint-_const_pi2*q)); + double yCentre = double(pos[1]) + radius*double(sin(_phiMomRefPoint-_const_pi2*q)); + + double d0; + + if (q>0) { + d0 = double(q)*radius - double(sqrt(xCentre*xCentre+yCentre*yCentre)); + } + else { + d0 = double(q)*radius + double(sqrt(xCentre*xCentre+yCentre*yCentre)); + } + + _d0 = float(d0); + +// if (fabs(_d0)>0.001 ) { +// std::cout << "New helix : " << std::endl; +// std::cout << " Position : " << pos[0] +// << " " << pos[1] +// << " " << pos[2] << std::endl; +// std::cout << " Radius = " << _radius << std::endl; +// std::cout << " RC = " << sqrt(_xCentre*_xCentre+_yCentre*_yCentre) << std::endl; +// std::cout << " D0 = " << _d0 << std::endl; +// } + + _pxAtPCA = _pxy*cos(_phi0); + _pyAtPCA = _pxy*sin(_phi0); + float deltaPhi = _phiRefPoint - _phiAtPCA; + float xCircles = -pos[2]*q/(_radius*_tanLambda) - deltaPhi; + xCircles = xCircles/_const_2pi; + int nCircles; + int n1,n2; + + if (xCircles >= 0.) { + n1 = int(xCircles); + n2 = n1 + 1; + } + else { + n1 = int(xCircles) - 1; + n2 = n1 + 1; + } + + if (fabs(n1-xCircles) < fabs(n2-xCircles)) { + nCircles = n1; + } + else { + nCircles = n2; + } + _z0 = pos[2] + _radius*_tanLambda*q*(deltaPhi + _const_2pi*nCircles); + +} + +void HelixClass::Initialize_Canonical(float phi0, float d0, float z0, + float omega, float tanLambda, float B) { + _omega = omega; + _d0 = d0; + _phi0 = phi0; + _z0 = z0; + _tanLambda = tanLambda; + _charge = omega/fabs(omega); + _radius = 1./fabs(omega); + _xAtPCA = -_d0*sin(_phi0); + _yAtPCA = _d0*cos(_phi0); + _referencePoint[0] = _xAtPCA; + _referencePoint[1] = _yAtPCA; + _referencePoint[2] = _z0; + _pxy = _FCT*B*_radius; + _momentum[0] = _pxy*cos(_phi0); + _momentum[1] = _pxy*sin(_phi0); + _momentum[2] = _tanLambda * _pxy; + _pxAtPCA = _momentum[0]; + _pyAtPCA = _momentum[1]; + _phiMomRefPoint = atan2(_momentum[1],_momentum[0]); + _xCentre = _referencePoint[0] + + _radius*cos(_phi0-_const_pi2*_charge); + _yCentre = _referencePoint[1] + + _radius*sin(_phi0-_const_pi2*_charge); + _phiAtPCA = atan2(-_yCentre,-_xCentre); + _phiRefPoint = _phiAtPCA ; + _bField = B; +} + + +void HelixClass::Initialize_BZ(float xCentre, float yCentre, float radius, + float bZ, float phi0, float B, float signPz, + float zBegin) { + + _radius = radius; + _pxy = _FCT*B*_radius; + _charge = -(bZ*signPz)/fabs(bZ*signPz); + _momentum[2] = -_charge*_pxy/(bZ*_radius); + _xCentre = xCentre; + _yCentre = yCentre; + _omega = _charge/radius; + _phiAtPCA = atan2(-_yCentre,-_xCentre); + _phi0 = -_const_pi2*_charge + _phiAtPCA; + while (_phi0<0) _phi0+=_const_2pi; + while (_phi0>=_const_2pi) _phi0-=_const_2pi; + _xAtPCA = _xCentre + _radius*cos(_phiAtPCA); + _yAtPCA = _yCentre + _radius*sin(_phiAtPCA); + _d0 = -_xAtPCA*sin(_phi0) + _yAtPCA*cos(_phi0); + _pxAtPCA = _pxy*cos(_phi0); + _pyAtPCA = _pxy*sin(_phi0); + _referencePoint[2] = zBegin; + _referencePoint[0] = xCentre + radius*cos(bZ*zBegin+phi0); + _referencePoint[1] = yCentre + radius*sin(bZ*zBegin+phi0); + _phiRefPoint = atan2(_referencePoint[1]-_yCentre,_referencePoint[0]-_xCentre); + _phiMomRefPoint = -_const_pi2*_charge + _phiRefPoint; + _tanLambda = _momentum[2]/_pxy; + _momentum[0] = _pxy*cos(_phiMomRefPoint); + _momentum[1] = _pxy*sin(_phiMomRefPoint); + + float deltaPhi = _phiRefPoint - _phiAtPCA; + float xCircles = bZ*_referencePoint[2] - deltaPhi; + xCircles = xCircles/_const_2pi; + int nCircles; + int n1,n2; + + if (xCircles >= 0.) { + n1 = int(xCircles); + n2 = n1 + 1; + } + else { + n1 = int(xCircles) - 1; + n2 = n1 + 1; + } + + if (fabs(n1-xCircles) < fabs(n2-xCircles)) { + nCircles = n1; + } + else { + nCircles = n2; + } + _z0 = _referencePoint[2] - (deltaPhi + _const_2pi*nCircles)/bZ; + _bField = B; + +} + +const float * HelixClass::getMomentum() { + return _momentum; +} +const float * HelixClass::getReferencePoint() { + return _referencePoint; +} +float HelixClass::getPhi0() { + if (_phi0<0.0) + _phi0 += 2*M_PI; + return _phi0; +} +float HelixClass::getD0() { + return _d0; +} +float HelixClass::getZ0() { + return _z0; +} +float HelixClass::getOmega() { + return _omega; +} +float HelixClass::getTanLambda() { + return _tanLambda; +} +float HelixClass::getPXY() { + return _pxy; +} +float HelixClass::getXC() { + return _xCentre; +} + +float HelixClass::getYC() { + return _yCentre; +} + +float HelixClass::getRadius() { + return _radius; +} + +float HelixClass::getBz() { + return _bZ; +} + +float HelixClass::getPhiZ() { + return _phiZ; +} + +float HelixClass::getCharge() { + return _charge; +} + +float HelixClass::getPointInXY(float x0, float y0, float ax, float ay, + float * ref , float * point) { + + float time; + + float AA = sqrt(ax*ax+ay*ay); + + + if (AA <= 0) { + time = -1.0e+20; + return time; + } + + + float BB = ax*(x0-_xCentre) + ay*(y0-_yCentre); + BB = BB / AA; + + float CC = (x0-_xCentre)*(x0-_xCentre) + + (y0-_yCentre)*(y0-_yCentre) - _radius*_radius; + + CC = CC / AA; + + float DET = BB*BB - CC; + float tt1 = 0.; + float tt2 = 0.; + float xx1,xx2,yy1,yy2; + + + if (DET < 0 ) { + time = 1.0e+10; + point[0]=0.0; + point[1]=0.0; + point[2]=0.0; + return time; + } + + + tt1 = - BB + sqrt(DET); + tt2 = - BB - sqrt(DET); + + xx1 = x0 + tt1*ax; + yy1 = y0 + tt1*ay; + xx2 = x0 + tt2*ax; + yy2 = y0 + tt2*ay; + + float phi1 = atan2(yy1-_yCentre,xx1-_xCentre); + float phi2 = atan2(yy2-_yCentre,xx2-_xCentre); + float phi0 = atan2(ref[1]-_yCentre,ref[0]-_xCentre); + + float dphi1 = phi1 - phi0; + float dphi2 = phi2 - phi0; + + if (dphi1 < 0 && _charge < 0) { + dphi1 = dphi1 + _const_2pi; + } + else if (dphi1 > 0 && _charge > 0) { + dphi1 = dphi1 - _const_2pi; + } + + if (dphi2 < 0 && _charge < 0) { + dphi2 = dphi2 + _const_2pi; + } + else if (dphi2 > 0 && _charge > 0) { + dphi2 = dphi2 - _const_2pi; + } + + // Times + tt1 = -_charge*dphi1*_radius/_pxy; + tt2 = -_charge*dphi2*_radius/_pxy; + + if (tt1 < 0. ) + std::cout << "WARNING " << tt1 << std::endl; + if (tt2 < 0. ) + std::cout << "WARNING " << tt2 << std::endl; + + + if (tt1 < tt2) { + point[0] = xx1; + point[1] = yy1; + time = tt1; + } + else { + point[0] = xx2; + point[1] = yy2; + time = tt2; + } + + point[2] = ref[2]+time*_momentum[2]; + + + + return time; + +} + + +float HelixClass::getPointOnCircle(float Radius, float * ref, float * point) { + + float distCenterToIP = sqrt(_xCentre*_xCentre + _yCentre*_yCentre); + + point[0] = 0.0; + point[1] = 0.0; + point[2] = 0.0; + + if ((distCenterToIP+_radius)<Radius) { + float xx = -1.0e+20; + return xx; + } + + if ((_radius+Radius)<distCenterToIP) { + float xx = -1.0e+20; + return xx; + } + + float phiCentre = atan2(_yCentre,_xCentre); + float phiStar = Radius*Radius + distCenterToIP*distCenterToIP + - _radius*_radius; + + phiStar = 0.5*phiStar/fmax(1.0e-20,Radius*distCenterToIP); + + if (phiStar > 1.0) + phiStar = 0.9999999; + + if (phiStar < -1.0) + phiStar = -0.9999999; + + phiStar = acos(phiStar); + + float tt1,tt2,time; + + float xx1 = Radius*cos(phiCentre+phiStar); + float yy1 = Radius*sin(phiCentre+phiStar); + + float xx2 = Radius*cos(phiCentre-phiStar); + float yy2 = Radius*sin(phiCentre-phiStar); + + + float phi1 = atan2(yy1-_yCentre,xx1-_xCentre); + float phi2 = atan2(yy2-_yCentre,xx2-_xCentre); + float phi0 = atan2(ref[1]-_yCentre,ref[0]-_xCentre); + + float dphi1 = phi1 - phi0; + float dphi2 = phi2 - phi0; + + if (dphi1 < 0 && _charge < 0) { + dphi1 = dphi1 + _const_2pi; + } + else if (dphi1 > 0 && _charge > 0) { + dphi1 = dphi1 - _const_2pi; + } + + if (dphi2 < 0 && _charge < 0) { + dphi2 = dphi2 + _const_2pi; + } + else if (dphi2 > 0 && _charge > 0) { + dphi2 = dphi2 - _const_2pi; + } + + // Times + tt1 = -_charge*dphi1*_radius/_pxy; + tt2 = -_charge*dphi2*_radius/_pxy; + + if (tt1 < 0. ) + std::cout << "WARNING " << tt1 << std::endl; + if (tt2 < 0. ) + std::cout << "WARNING " << tt2 << std::endl; + + + float time2; + if (tt1 < tt2) { + point[0] = xx1; + point[1] = yy1; + point[3] = xx2; + point[4] = yy2; + time = tt1; + time2 = tt2; + } + else { + point[0] = xx2; + point[1] = yy2; + point[3] = xx1; + point[4] = yy1; + time = tt2; + time2 = tt1; + } + + point[2] = ref[2]+time*_momentum[2]; + point[5] = ref[2]+time2*_momentum[2]; + + + return time; + +} + + +float HelixClass::getPointInZ(float zLine, float * ref, float * point) { + + float time = zLine - ref[2]; + + if (_momentum[2] == 0.) { + time = -1.0e+20; + point[0] = 0.; + point[1] = 0.; + point[2] = 0.; + return time; + } + + time = time/_momentum[2]; + + float phi0 = atan2(ref[1] - _yCentre , ref[0] - _xCentre); + float phi = phi0 - _charge*_pxy*time/_radius; + float xx = _xCentre + _radius*cos(phi); + float yy = _yCentre + _radius*sin(phi); + + point[0] = xx; + point[1] = yy; + point[2] = zLine; + + return time; + + +} + +float HelixClass::getDistanceToPoint(float * xPoint, float * Distance) { + + float zOnHelix; + float phi = atan2(xPoint[1]-_yCentre,xPoint[0]-_xCentre); + float phi0 = atan2(_referencePoint[1]-_yCentre,_referencePoint[0]-_xCentre); + //calculate distance to XYprojected centre of Helix, comparing this with distance to radius around centre gives DistXY + float DistXY = (_xCentre-xPoint[0])*(_xCentre-xPoint[0]) + (_yCentre-xPoint[1])*(_yCentre-xPoint[1]); + DistXY = sqrt(DistXY); + DistXY = fabs(DistXY - _radius); + + int nCircles = 0; + + if (fabs(_tanLambda*_radius)>1.0e-20) { + float xCircles = phi0 - phi -_charge*(xPoint[2]-_referencePoint[2])/(_tanLambda*_radius); + xCircles = xCircles/_const_2pi; + int n1,n2; + + if (xCircles >= 0.) { + n1 = int(xCircles); + n2 = n1 + 1; + } + else { + n1 = int(xCircles) - 1; + n2 = n1 + 1; + } + + if (fabs(n1-xCircles) < fabs(n2-xCircles)) { + nCircles = n1; + } + else { + nCircles = n2; + } + + } + + float DPhi = _const_2pi*((float)nCircles) + phi - phi0; + + zOnHelix = _referencePoint[2] - _charge*_radius*_tanLambda*DPhi; + + float DistZ = fabs(zOnHelix - xPoint[2]); + float Time; + + if (fabs(_momentum[2]) > 1.0e-20) { + Time = (zOnHelix - _referencePoint[2])/_momentum[2]; + } + else { + Time = _charge*_radius*DPhi/_pxy; + } + + Distance[0] = DistXY; + Distance[1] = DistZ; + Distance[2] = sqrt(DistXY*DistXY+DistZ*DistZ); + + return Time; + + +} + +//When we are not interested in the exact distance, we can check if we are +//already far enough away in XY, before we start calculating in Z as the +//distance will only increase +float HelixClass::getDistanceToPoint(const std::vector<float>& xPoint, float distCut) { + //calculate distance to XYprojected centre of Helix, comparing this with distance to radius around centre gives DistXY + float tempx = xPoint[0]-_xCentre; + float tempy = xPoint[1]-_yCentre; + float tempsq = sqrt(tempx*tempx + tempy*tempy); + float tempdf = tempsq - _radius; + float DistXY = fabs( tempdf ); + //If this is bigger than distCut, we dont have to know how much bigger this is + if( DistXY > distCut) { + return DistXY; + } + + int nCircles = 0; + float phi = atan2(tempy,tempx); + float phi0 = atan2(_referencePoint[1]-_yCentre,_referencePoint[0]-_xCentre); + float phidiff = phi0-phi; + float tempz = xPoint[2] - _referencePoint[2];//Yes referencePoint + float tanradius = _tanLambda*_radius; + if (fabs(tanradius)>1.0e-20) { + float xCircles = (phidiff -_charge*tempz/tanradius)/_const_2pi; + int n1,n2; + if (xCircles >= 0.) { + n1 = int(xCircles); + n2 = n1 + 1; + } + else { + n1 = int(xCircles) - 1; + n2 = n1 + 1; + } + if (fabs(n1-xCircles) < fabs(n2-xCircles)) { + nCircles = n1; + } + else { + nCircles = n2; + } + } + float DistZ = - tempz - _charge*tanradius*(_const_2pi*((float)nCircles) - phidiff); + return sqrt(DistXY*DistXY+DistZ*DistZ); +}//getDistanceToPoint(vector,float) + +float HelixClass::getDistanceToPoint(const float* xPoint, float distCut) { + std::vector<float> xPosition(xPoint, xPoint + 3 );//We are expecting three coordinates, must be +3, last element is excluded! + return getDistanceToPoint(xPosition, distCut); +}//getDistanceToPoint(float*,float) + + + +void HelixClass::setHelixEdges(float * xStart, float * xEnd) { + for (int i=0; i<3; ++i) { + _xStart[i] = xStart[i]; + _xEnd[i] = xEnd[i]; + } + +} + +float HelixClass::getDistanceToHelix(HelixClass * helix, float * pos, float * mom) { + + float x01 = getXC(); + float y01 = getYC(); + + float x02 = helix->getXC(); + float y02 = helix->getYC(); + + float rad1 = getRadius(); + float rad2 = helix->getRadius(); + + float distance = sqrt((x01-x02)*(x01-x02)+(y01-y02)*(y01-y02)); + + bool singlePoint = true; + + float phi1 = 0; + float phi2 = 0; + + if (rad1+rad2<distance) { + phi1 = atan2(y02-y01,x02-x01); + phi2 = atan2(y01-y02,x01-x02); + } + else if (distance+rad2<rad1) { + phi1 = atan2(y02-y01,x02-x01); + phi2 = phi1; + } + else if (distance+rad1<rad2) { + phi1 = atan2(y01-y02,x01-x02); + phi2 = phi1; + } + else { + singlePoint = false; + float cosAlpha = 0.5*(distance*distance+rad2*rad2-rad1*rad1)/(distance*rad2); + float alpha = acos(cosAlpha); + float phi0 = atan2(y01-y02,x01-x02); + phi1 = phi0 + alpha; + phi2 = phi0 - alpha; + } + + + float ref1[3]; + float ref2[3]; + for (int i=0;i<3;++i) { + ref1[i]=_referencePoint[i]; + ref2[i]=helix->getReferencePoint()[i]; + } + + float pos1[3]; + float pos2[3]; + float mom1[3]; + float mom2[3]; + + + if (singlePoint ) { + + float xSect1 = x01 + rad1*cos(phi1); + float ySect1 = y01 + rad1*sin(phi1); + float xSect2 = x02 + rad2*cos(phi2); + float ySect2 = y02 + rad2*sin(phi2); + float R1 = sqrt(xSect1*xSect1+ySect1*ySect1); + float R2 = sqrt(xSect2*xSect2+ySect2*ySect2); + + getPointOnCircle(R1,ref1,pos1); + helix->getPointOnCircle(R2,ref2,pos2); + + } + else { + + float xSect1 = x02 + rad2*cos(phi1); + float ySect1 = y02 + rad2*sin(phi1); + float xSect2 = x02 + rad2*cos(phi2); + float ySect2 = y02 + rad2*sin(phi2); + +// std::cout << "(xSect1,ySect1)=(" << xSect1 << "," << ySect1 << ")" << std::endl; +// std::cout << "(xSect2,ySect2)=(" << xSect2 << "," << ySect2 << ")" << std::endl; + + float temp21[3]; + float temp22[3]; + + float phiI2 = atan2(ref2[1]-y02,ref2[0]-x02); + float phiF21 = atan2(ySect1-y02,xSect1-x02); + float phiF22 = atan2(ySect2-y02,xSect2-x02); + float deltaPhi21 = phiF21 - phiI2; + float deltaPhi22 = phiF22 - phiI2; + float charge2 = helix->getCharge(); + float pxy2 = helix->getPXY(); + float pz2 = helix->getMomentum()[2]; + if (deltaPhi21 < 0 && charge2 < 0) { + deltaPhi21 += _const_2pi; + } + else if (deltaPhi21 > 0 && charge2 > 0) { + deltaPhi21 -= _const_2pi; + } + + if (deltaPhi22 < 0 && charge2 < 0) { + deltaPhi22 += _const_2pi; + } + else if (deltaPhi22 > 0 && charge2 > 0) { + deltaPhi22 -= _const_2pi; + } + + float time21 = -charge2*deltaPhi21*rad2/pxy2; + float time22 = -charge2*deltaPhi22*rad2/pxy2; + + float Z21 = ref2[2] + time21*pz2; + float Z22 = ref2[2] + time22*pz2; + + temp21[0] = xSect1; temp21[1] = ySect1; temp21[2] = Z21; + temp22[0] = xSect2; temp22[1] = ySect2; temp22[2] = Z22; + + +// std::cout << "temp21 = " << temp21[0] << " " << temp21[1] << " " << temp21[2] << std::endl; +// std::cout << "temp22 = " << temp22[0] << " " << temp22[1] << " " << temp22[2] << std::endl; + + + float temp11[3]; + float temp12[3]; + + float phiI1 = atan2(ref1[1]-y01,ref1[0]-x01); + float phiF11 = atan2(ySect1-y01,xSect1-x01); + float phiF12 = atan2(ySect2-y01,xSect2-x01); + float deltaPhi11 = phiF11 - phiI1; + float deltaPhi12 = phiF12 - phiI1; + float charge1 = _charge; + float pxy1 = _pxy; + float pz1 = _momentum[2]; + if (deltaPhi11 < 0 && charge1 < 0) { + deltaPhi11 += _const_2pi; + } + else if (deltaPhi11 > 0 && charge1 > 0) { + deltaPhi11 -= _const_2pi; + } + + if (deltaPhi12 < 0 && charge1 < 0) { + deltaPhi12 += _const_2pi; + } + else if (deltaPhi12 > 0 && charge1 > 0) { + deltaPhi12 -= _const_2pi; + } + + float time11 = -charge1*deltaPhi11*rad1/pxy1; + float time12 = -charge1*deltaPhi12*rad1/pxy1; + + float Z11 = ref1[2] + time11*pz1; + float Z12 = ref1[2] + time12*pz1; + + temp11[0] = xSect1; temp11[1] = ySect1; temp11[2] = Z11; + temp12[0] = xSect2; temp12[1] = ySect2; temp12[2] = Z12; + +// std::cout << "temp11 = " << temp11[0] << " " << temp11[1] << " " << temp11[2] << std::endl; +// std::cout << "temp12 = " << temp12[0] << " " << temp12[1] << " " << temp12[2] << std::endl; + + float Dist1 = 0; + float Dist2 = 0; + + for (int j=0;j<3;++j) { + Dist1 += (temp11[j]-temp21[j])*(temp11[j]-temp21[j]); + Dist2 += (temp12[j]-temp22[j])*(temp12[j]-temp22[j]); + } + + if (Dist1<Dist2) { + for (int l=0;l<3;++l) { + pos1[l] = temp11[l]; + pos2[l] = temp21[l]; + } + } + else { + for (int l=0;l<3;++l) { + pos1[l] = temp12[l]; + pos2[l] = temp22[l]; + } + } + + } + + getExtrapolatedMomentum(pos1,mom1); + helix->getExtrapolatedMomentum(pos2,mom2); + + float helixDistance = 0; + + for (int i=0;i<3;++i) { + helixDistance += (pos1[i] - pos2[i])*(pos1[i] - pos2[i]); + pos[i] = 0.5*(pos1[i]+pos2[i]); + mom[i] = mom1[i]+mom2[i]; + } + helixDistance = sqrt(helixDistance); + + return helixDistance; + +} + +void HelixClass::getExtrapolatedMomentum(float * pos, float * momentum) { + + float phi = atan2(pos[1]-_yCentre,pos[0]-_xCentre); + if (phi <0.) phi += _const_2pi; + phi = phi - _phiAtPCA + _phi0; + momentum[0] = _pxy*cos(phi); + momentum[1] = _pxy*sin(phi); + momentum[2] = _momentum[2]; + + +} diff --git a/Utilities/DataHelper/src/HelixClass.cc~ b/Utilities/DataHelper/src/HelixClass.cc~ new file mode 100644 index 00000000..d16ac955 --- /dev/null +++ b/Utilities/DataHelper/src/HelixClass.cc~ @@ -0,0 +1,768 @@ +#include "HelixClass.h" +#include <math.h> +#include <stdlib.h> +#include <iostream> +//#include "ced_cli.h" + +HelixClass::HelixClass() { + _const_2pi = 2.0*M_PI; + _const_pi2 = 0.5*M_PI; + _FCT = 2.99792458E-4; +} + +HelixClass::~HelixClass() {} + +void HelixClass::Initialize_VP(float * pos, float * mom, float q, float B) { + _referencePoint[0] = pos[0]; + _referencePoint[1] = pos[1]; + _referencePoint[2] = pos[2]; + _momentum[0] = mom[0]; + _momentum[1] = mom[1]; + _momentum[2] = mom[2]; + _charge = q; + _bField = B; + _pxy = sqrt(mom[0]*mom[0]+mom[1]*mom[1]); + _radius = _pxy / (_FCT*B); + _omega = q/_radius; + _tanLambda = mom[2]/_pxy; + _phiMomRefPoint = atan2(mom[1],mom[0]); + _xCentre = pos[0] + _radius*cos(_phiMomRefPoint-_const_pi2*q); + _yCentre = pos[1] + _radius*sin(_phiMomRefPoint-_const_pi2*q); + _phiRefPoint = atan2(pos[1]-_yCentre,pos[0]-_xCentre); + _phiAtPCA = atan2(-_yCentre,-_xCentre); + _phi0 = -_const_pi2*q + _phiAtPCA; + while (_phi0<0) _phi0+=_const_2pi; + while (_phi0>=_const_2pi) _phi0-=_const_2pi; + _xAtPCA = _xCentre + _radius*cos(_phiAtPCA); + _yAtPCA = _yCentre + _radius*sin(_phiAtPCA); + // _d0 = -_xAtPCA*sin(_phi0) + _yAtPCA*cos(_phi0); + double pxy = double(_pxy); + double radius = pxy/double(_FCT*B); + double xCentre = double(pos[0]) + radius*double(cos(_phiMomRefPoint-_const_pi2*q)); + double yCentre = double(pos[1]) + radius*double(sin(_phiMomRefPoint-_const_pi2*q)); + + double d0; + + if (q>0) { + d0 = double(q)*radius - double(sqrt(xCentre*xCentre+yCentre*yCentre)); + } + else { + d0 = double(q)*radius + double(sqrt(xCentre*xCentre+yCentre*yCentre)); + } + + _d0 = float(d0); + +// if (fabs(_d0)>0.001 ) { +// std::cout << "New helix : " << std::endl; +// std::cout << " Position : " << pos[0] +// << " " << pos[1] +// << " " << pos[2] << std::endl; +// std::cout << " Radius = " << _radius << std::endl; +// std::cout << " RC = " << sqrt(_xCentre*_xCentre+_yCentre*_yCentre) << std::endl; +// std::cout << " D0 = " << _d0 << std::endl; +// } + + _pxAtPCA = _pxy*cos(_phi0); + _pyAtPCA = _pxy*sin(_phi0); + float deltaPhi = _phiRefPoint - _phiAtPCA; + float xCircles = -pos[2]*q/(_radius*_tanLambda) - deltaPhi; + xCircles = xCircles/_const_2pi; + int nCircles; + int n1,n2; + + if (xCircles >= 0.) { + n1 = int(xCircles); + n2 = n1 + 1; + } + else { + n1 = int(xCircles) - 1; + n2 = n1 + 1; + } + + if (fabs(n1-xCircles) < fabs(n2-xCircles)) { + nCircles = n1; + } + else { + nCircles = n2; + } + _z0 = pos[2] + _radius*_tanLambda*q*(deltaPhi + _const_2pi*nCircles); + +} + +void HelixClass::Initialize_Canonical(float phi0, float d0, float z0, + float omega, float tanLambda, float B) { + _omega = omega; + _d0 = d0; + _phi0 = phi0; + _z0 = z0; + _tanLambda = tanLambda; + _charge = omega/fabs(omega); + _radius = 1./fabs(omega); + _xAtPCA = -_d0*sin(_phi0); + _yAtPCA = _d0*cos(_phi0); + _referencePoint[0] = _xAtPCA; + _referencePoint[1] = _yAtPCA; + _referencePoint[2] = _z0; + _pxy = _FCT*B*_radius; + _momentum[0] = _pxy*cos(_phi0); + _momentum[1] = _pxy*sin(_phi0); + _momentum[2] = _tanLambda * _pxy; + _pxAtPCA = _momentum[0]; + _pyAtPCA = _momentum[1]; + _phiMomRefPoint = atan2(_momentum[1],_momentum[0]); + _xCentre = _referencePoint[0] + + _radius*cos(_phi0-_const_pi2*_charge); + _yCentre = _referencePoint[1] + + _radius*sin(_phi0-_const_pi2*_charge); + _phiAtPCA = atan2(-_yCentre,-_xCentre); + _phiRefPoint = _phiAtPCA ; + _bField = B; +} + + +void HelixClass::Initialize_BZ(float xCentre, float yCentre, float radius, + float bZ, float phi0, float B, float signPz, + float zBegin) { + + _radius = radius; + _pxy = _FCT*B*_radius; + _charge = -(bZ*signPz)/fabs(bZ*signPz); + _momentum[2] = -_charge*_pxy/(bZ*_radius); + _xCentre = xCentre; + _yCentre = yCentre; + _omega = _charge/radius; + _phiAtPCA = atan2(-_yCentre,-_xCentre); + _phi0 = -_const_pi2*_charge + _phiAtPCA; + while (_phi0<0) _phi0+=_const_2pi; + while (_phi0>=_const_2pi) _phi0-=_const_2pi; + _xAtPCA = _xCentre + _radius*cos(_phiAtPCA); + _yAtPCA = _yCentre + _radius*sin(_phiAtPCA); + _d0 = -_xAtPCA*sin(_phi0) + _yAtPCA*cos(_phi0); + _pxAtPCA = _pxy*cos(_phi0); + _pyAtPCA = _pxy*sin(_phi0); + _referencePoint[2] = zBegin; + _referencePoint[0] = xCentre + radius*cos(bZ*zBegin+phi0); + _referencePoint[1] = yCentre + radius*sin(bZ*zBegin+phi0); + _phiRefPoint = atan2(_referencePoint[1]-_yCentre,_referencePoint[0]-_xCentre); + _phiMomRefPoint = -_const_pi2*_charge + _phiRefPoint; + _tanLambda = _momentum[2]/_pxy; + _momentum[0] = _pxy*cos(_phiMomRefPoint); + _momentum[1] = _pxy*sin(_phiMomRefPoint); + + float deltaPhi = _phiRefPoint - _phiAtPCA; + float xCircles = bZ*_referencePoint[2] - deltaPhi; + xCircles = xCircles/_const_2pi; + int nCircles; + int n1,n2; + + if (xCircles >= 0.) { + n1 = int(xCircles); + n2 = n1 + 1; + } + else { + n1 = int(xCircles) - 1; + n2 = n1 + 1; + } + + if (fabs(n1-xCircles) < fabs(n2-xCircles)) { + nCircles = n1; + } + else { + nCircles = n2; + } + _z0 = _referencePoint[2] - (deltaPhi + _const_2pi*nCircles)/bZ; + _bField = B; + +} + +const float * HelixClass::getMomentum() { + return _momentum; +} +const float * HelixClass::getReferencePoint() { + return _referencePoint; +} +float HelixClass::getPhi0() { + if (_phi0<0.0) + _phi0 += 2*M_PI; + return _phi0; +} +float HelixClass::getD0() { + return _d0; +} +float HelixClass::getZ0() { + return _z0; +} +float HelixClass::getOmega() { + return _omega; +} +float HelixClass::getTanLambda() { + return _tanLambda; +} +float HelixClass::getPXY() { + return _pxy; +} +float HelixClass::getXC() { + return _xCentre; +} + +float HelixClass::getYC() { + return _yCentre; +} + +float HelixClass::getRadius() { + return _radius; +} + +float HelixClass::getBz() { + return _bZ; +} + +float HelixClass::getPhiZ() { + return _phiZ; +} + +float HelixClass::getCharge() { + return _charge; +} + +float HelixClass::getPointInXY(float x0, float y0, float ax, float ay, + float * ref , float * point) { + + float time; + + float AA = sqrt(ax*ax+ay*ay); + + + if (AA <= 0) { + time = -1.0e+20; + return time; + } + + + float BB = ax*(x0-_xCentre) + ay*(y0-_yCentre); + BB = BB / AA; + + float CC = (x0-_xCentre)*(x0-_xCentre) + + (y0-_yCentre)*(y0-_yCentre) - _radius*_radius; + + CC = CC / AA; + + float DET = BB*BB - CC; + float tt1 = 0.; + float tt2 = 0.; + float xx1,xx2,yy1,yy2; + + + if (DET < 0 ) { + time = 1.0e+10; + point[0]=0.0; + point[1]=0.0; + point[2]=0.0; + return time; + } + + + tt1 = - BB + sqrt(DET); + tt2 = - BB - sqrt(DET); + + xx1 = x0 + tt1*ax; + yy1 = y0 + tt1*ay; + xx2 = x0 + tt2*ax; + yy2 = y0 + tt2*ay; + + float phi1 = atan2(yy1-_yCentre,xx1-_xCentre); + float phi2 = atan2(yy2-_yCentre,xx2-_xCentre); + float phi0 = atan2(ref[1]-_yCentre,ref[0]-_xCentre); + + float dphi1 = phi1 - phi0; + float dphi2 = phi2 - phi0; + + if (dphi1 < 0 && _charge < 0) { + dphi1 = dphi1 + _const_2pi; + } + else if (dphi1 > 0 && _charge > 0) { + dphi1 = dphi1 - _const_2pi; + } + + if (dphi2 < 0 && _charge < 0) { + dphi2 = dphi2 + _const_2pi; + } + else if (dphi2 > 0 && _charge > 0) { + dphi2 = dphi2 - _const_2pi; + } + + // Times + tt1 = -_charge*dphi1*_radius/_pxy; + tt2 = -_charge*dphi2*_radius/_pxy; + + if (tt1 < 0. ) + std::cout << "WARNING " << tt1 << std::endl; + if (tt2 < 0. ) + std::cout << "WARNING " << tt2 << std::endl; + + + if (tt1 < tt2) { + point[0] = xx1; + point[1] = yy1; + time = tt1; + } + else { + point[0] = xx2; + point[1] = yy2; + time = tt2; + } + + point[2] = ref[2]+time*_momentum[2]; + + + + return time; + +} + + +float HelixClass::getPointOnCircle(float Radius, float * ref, float * point) { + + float distCenterToIP = sqrt(_xCentre*_xCentre + _yCentre*_yCentre); + + point[0] = 0.0; + point[1] = 0.0; + point[2] = 0.0; + + if ((distCenterToIP+_radius)<Radius) { + float xx = -1.0e+20; + return xx; + } + + if ((_radius+Radius)<distCenterToIP) { + float xx = -1.0e+20; + return xx; + } + + float phiCentre = atan2(_yCentre,_xCentre); + float phiStar = Radius*Radius + distCenterToIP*distCenterToIP + - _radius*_radius; + + phiStar = 0.5*phiStar/fmax(1.0e-20,Radius*distCenterToIP); + + if (phiStar > 1.0) + phiStar = 0.9999999; + + if (phiStar < -1.0) + phiStar = -0.9999999; + + phiStar = acos(phiStar); + + float tt1,tt2,time; + + float xx1 = Radius*cos(phiCentre+phiStar); + float yy1 = Radius*sin(phiCentre+phiStar); + + float xx2 = Radius*cos(phiCentre-phiStar); + float yy2 = Radius*sin(phiCentre-phiStar); + + + float phi1 = atan2(yy1-_yCentre,xx1-_xCentre); + float phi2 = atan2(yy2-_yCentre,xx2-_xCentre); + float phi0 = atan2(ref[1]-_yCentre,ref[0]-_xCentre); + + float dphi1 = phi1 - phi0; + float dphi2 = phi2 - phi0; + + if (dphi1 < 0 && _charge < 0) { + dphi1 = dphi1 + _const_2pi; + } + else if (dphi1 > 0 && _charge > 0) { + dphi1 = dphi1 - _const_2pi; + } + + if (dphi2 < 0 && _charge < 0) { + dphi2 = dphi2 + _const_2pi; + } + else if (dphi2 > 0 && _charge > 0) { + dphi2 = dphi2 - _const_2pi; + } + + // Times + tt1 = -_charge*dphi1*_radius/_pxy; + tt2 = -_charge*dphi2*_radius/_pxy; + + if (tt1 < 0. ) + std::cout << "WARNING " << tt1 << std::endl; + if (tt2 < 0. ) + std::cout << "WARNING " << tt2 << std::endl; + + + float time2; + if (tt1 < tt2) { + point[0] = xx1; + point[1] = yy1; + point[3] = xx2; + point[4] = yy2; + time = tt1; + time2 = tt2; + } + else { + point[0] = xx2; + point[1] = yy2; + point[3] = xx1; + point[4] = yy1; + time = tt2; + time2 = tt1; + } + + point[2] = ref[2]+time*_momentum[2]; + point[5] = ref[2]+time2*_momentum[2]; + + + return time; + +} + + +float HelixClass::getPointInZ(float zLine, float * ref, float * point) { + + float time = zLine - ref[2]; + + if (_momentum[2] == 0.) { + time = -1.0e+20; + point[0] = 0.; + point[1] = 0.; + point[2] = 0.; + return time; + } + + time = time/_momentum[2]; + + float phi0 = atan2(ref[1] - _yCentre , ref[0] - _xCentre); + float phi = phi0 - _charge*_pxy*time/_radius; + float xx = _xCentre + _radius*cos(phi); + float yy = _yCentre + _radius*sin(phi); + + point[0] = xx; + point[1] = yy; + point[2] = zLine; + + return time; + + +} + +float HelixClass::getDistanceToPoint(float * xPoint, float * Distance) { + + float zOnHelix; + float phi = atan2(xPoint[1]-_yCentre,xPoint[0]-_xCentre); + float phi0 = atan2(_referencePoint[1]-_yCentre,_referencePoint[0]-_xCentre); + //calculate distance to XYprojected centre of Helix, comparing this with distance to radius around centre gives DistXY + float DistXY = (_xCentre-xPoint[0])*(_xCentre-xPoint[0]) + (_yCentre-xPoint[1])*(_yCentre-xPoint[1]); + DistXY = sqrt(DistXY); + DistXY = fabs(DistXY - _radius); + + int nCircles = 0; + + if (fabs(_tanLambda*_radius)>1.0e-20) { + float xCircles = phi0 - phi -_charge*(xPoint[2]-_referencePoint[2])/(_tanLambda*_radius); + xCircles = xCircles/_const_2pi; + int n1,n2; + + if (xCircles >= 0.) { + n1 = int(xCircles); + n2 = n1 + 1; + } + else { + n1 = int(xCircles) - 1; + n2 = n1 + 1; + } + + if (fabs(n1-xCircles) < fabs(n2-xCircles)) { + nCircles = n1; + } + else { + nCircles = n2; + } + + } + + float DPhi = _const_2pi*((float)nCircles) + phi - phi0; + + zOnHelix = _referencePoint[2] - _charge*_radius*_tanLambda*DPhi; + + float DistZ = fabs(zOnHelix - xPoint[2]); + float Time; + + if (fabs(_momentum[2]) > 1.0e-20) { + Time = (zOnHelix - _referencePoint[2])/_momentum[2]; + } + else { + Time = _charge*_radius*DPhi/_pxy; + } + + Distance[0] = DistXY; + Distance[1] = DistZ; + Distance[2] = sqrt(DistXY*DistXY+DistZ*DistZ); + + return Time; + + +} + +//When we are not interested in the exact distance, we can check if we are +//already far enough away in XY, before we start calculating in Z as the +//distance will only increase +float HelixClass::getDistanceToPoint(const std::vector<float>& xPoint, float distCut) { + //calculate distance to XYprojected centre of Helix, comparing this with distance to radius around centre gives DistXY + float tempx = xPoint[0]-_xCentre; + float tempy = xPoint[1]-_yCentre; + float tempsq = sqrt(tempx*tempx + tempy*tempy); + float tempdf = tempsq - _radius; + float DistXY = fabs( tempdf ); + //If this is bigger than distCut, we dont have to know how much bigger this is + if( DistXY > distCut) { + return DistXY; + } + + int nCircles = 0; + float phi = atan2(tempy,tempx); + float phi0 = atan2(_referencePoint[1]-_yCentre,_referencePoint[0]-_xCentre); + float phidiff = phi0-phi; + float tempz = xPoint[2] - _referencePoint[2];//Yes referencePoint + float tanradius = _tanLambda*_radius; + if (fabs(tanradius)>1.0e-20) { + float xCircles = (phidiff -_charge*tempz/tanradius)/_const_2pi; + int n1,n2; + if (xCircles >= 0.) { + n1 = int(xCircles); + n2 = n1 + 1; + } + else { + n1 = int(xCircles) - 1; + n2 = n1 + 1; + } + if (fabs(n1-xCircles) < fabs(n2-xCircles)) { + nCircles = n1; + } + else { + nCircles = n2; + } + } + float DistZ = - tempz - _charge*tanradius*(_const_2pi*((float)nCircles) - phidiff); + return sqrt(DistXY*DistXY+DistZ*DistZ); +}//getDistanceToPoint(vector,float) + +float HelixClass::getDistanceToPoint(const float* xPoint, float distCut) { + std::vector<float> xPosition(xPoint, xPoint + 3 );//We are expecting three coordinates, must be +3, last element is excluded! + return getDistanceToPoint(xPosition, distCut); +}//getDistanceToPoint(float*,float) + + + +void HelixClass::setHelixEdges(float * xStart, float * xEnd) { + for (int i=0; i<3; ++i) { + _xStart[i] = xStart[i]; + _xEnd[i] = xEnd[i]; + } + +} + +float HelixClass::getDistanceToHelix(HelixClass * helix, float * pos, float * mom) { + + float x01 = getXC(); + float y01 = getYC(); + + float x02 = helix->getXC(); + float y02 = helix->getYC(); + + float rad1 = getRadius(); + float rad2 = helix->getRadius(); + + float distance = sqrt((x01-x02)*(x01-x02)+(y01-y02)*(y01-y02)); + + bool singlePoint = true; + + float phi1 = 0; + float phi2 = 0; + + if (rad1+rad2<distance) { + phi1 = atan2(y02-y01,x02-x01); + phi2 = atan2(y01-y02,x01-x02); + } + else if (distance+rad2<rad1) { + phi1 = atan2(y02-y01,x02-x01); + phi2 = phi1; + } + else if (distance+rad1<rad2) { + phi1 = atan2(y01-y02,x01-x02); + phi2 = phi1; + } + else { + singlePoint = false; + float cosAlpha = 0.5*(distance*distance+rad2*rad2-rad1*rad1)/(distance*rad2); + float alpha = acos(cosAlpha); + float phi0 = atan2(y01-y02,x01-x02); + phi1 = phi0 + alpha; + phi2 = phi0 - alpha; + } + + + float ref1[3]; + float ref2[3]; + for (int i=0;i<3;++i) { + ref1[i]=_referencePoint[i]; + ref2[i]=helix->getReferencePoint()[i]; + } + + float pos1[3]; + float pos2[3]; + float mom1[3]; + float mom2[3]; + + + if (singlePoint ) { + + float xSect1 = x01 + rad1*cos(phi1); + float ySect1 = y01 + rad1*sin(phi1); + float xSect2 = x02 + rad2*cos(phi2); + float ySect2 = y02 + rad2*sin(phi2); + float R1 = sqrt(xSect1*xSect1+ySect1*ySect1); + float R2 = sqrt(xSect2*xSect2+ySect2*ySect2); + + getPointOnCircle(R1,ref1,pos1); + helix->getPointOnCircle(R2,ref2,pos2); + + } + else { + + float xSect1 = x02 + rad2*cos(phi1); + float ySect1 = y02 + rad2*sin(phi1); + float xSect2 = x02 + rad2*cos(phi2); + float ySect2 = y02 + rad2*sin(phi2); + +// std::cout << "(xSect1,ySect1)=(" << xSect1 << "," << ySect1 << ")" << std::endl; +// std::cout << "(xSect2,ySect2)=(" << xSect2 << "," << ySect2 << ")" << std::endl; + + float temp21[3]; + float temp22[3]; + + float phiI2 = atan2(ref2[1]-y02,ref2[0]-x02); + float phiF21 = atan2(ySect1-y02,xSect1-x02); + float phiF22 = atan2(ySect2-y02,xSect2-x02); + float deltaPhi21 = phiF21 - phiI2; + float deltaPhi22 = phiF22 - phiI2; + float charge2 = helix->getCharge(); + float pxy2 = helix->getPXY(); + float pz2 = helix->getMomentum()[2]; + if (deltaPhi21 < 0 && charge2 < 0) { + deltaPhi21 += _const_2pi; + } + else if (deltaPhi21 > 0 && charge2 > 0) { + deltaPhi21 -= _const_2pi; + } + + if (deltaPhi22 < 0 && charge2 < 0) { + deltaPhi22 += _const_2pi; + } + else if (deltaPhi22 > 0 && charge2 > 0) { + deltaPhi22 -= _const_2pi; + } + + float time21 = -charge2*deltaPhi21*rad2/pxy2; + float time22 = -charge2*deltaPhi22*rad2/pxy2; + + float Z21 = ref2[2] + time21*pz2; + float Z22 = ref2[2] + time22*pz2; + + temp21[0] = xSect1; temp21[1] = ySect1; temp21[2] = Z21; + temp22[0] = xSect2; temp22[1] = ySect2; temp22[2] = Z22; + + +// std::cout << "temp21 = " << temp21[0] << " " << temp21[1] << " " << temp21[2] << std::endl; +// std::cout << "temp22 = " << temp22[0] << " " << temp22[1] << " " << temp22[2] << std::endl; + + + float temp11[3]; + float temp12[3]; + + float phiI1 = atan2(ref1[1]-y01,ref1[0]-x01); + float phiF11 = atan2(ySect1-y01,xSect1-x01); + float phiF12 = atan2(ySect2-y01,xSect2-x01); + float deltaPhi11 = phiF11 - phiI1; + float deltaPhi12 = phiF12 - phiI1; + float charge1 = _charge; + float pxy1 = _pxy; + float pz1 = _momentum[2]; + if (deltaPhi11 < 0 && charge1 < 0) { + deltaPhi11 += _const_2pi; + } + else if (deltaPhi11 > 0 && charge1 > 0) { + deltaPhi11 -= _const_2pi; + } + + if (deltaPhi12 < 0 && charge1 < 0) { + deltaPhi12 += _const_2pi; + } + else if (deltaPhi12 > 0 && charge1 > 0) { + deltaPhi12 -= _const_2pi; + } + + float time11 = -charge1*deltaPhi11*rad1/pxy1; + float time12 = -charge1*deltaPhi12*rad1/pxy1; + + float Z11 = ref1[2] + time11*pz1; + float Z12 = ref1[2] + time12*pz1; + + temp11[0] = xSect1; temp11[1] = ySect1; temp11[2] = Z11; + temp12[0] = xSect2; temp12[1] = ySect2; temp12[2] = Z12; + +// std::cout << "temp11 = " << temp11[0] << " " << temp11[1] << " " << temp11[2] << std::endl; +// std::cout << "temp12 = " << temp12[0] << " " << temp12[1] << " " << temp12[2] << std::endl; + + float Dist1 = 0; + float Dist2 = 0; + + for (int j=0;j<3;++j) { + Dist1 += (temp11[j]-temp21[j])*(temp11[j]-temp21[j]); + Dist2 += (temp12[j]-temp22[j])*(temp12[j]-temp22[j]); + } + + if (Dist1<Dist2) { + for (int l=0;l<3;++l) { + pos1[l] = temp11[l]; + pos2[l] = temp21[l]; + } + } + else { + for (int l=0;l<3;++l) { + pos1[l] = temp12[l]; + pos2[l] = temp22[l]; + } + } + + } + + getExtrapolatedMomentum(pos1,mom1); + helix->getExtrapolatedMomentum(pos2,mom2); + + float helixDistance = 0; + + for (int i=0;i<3;++i) { + helixDistance += (pos1[i] - pos2[i])*(pos1[i] - pos2[i]); + pos[i] = 0.5*(pos1[i]+pos2[i]); + mom[i] = mom1[i]+mom2[i]; + } + helixDistance = sqrt(helixDistance); + + return helixDistance; + +} + +void HelixClass::getExtrapolatedMomentum(float * pos, float * momentum) { + + float phi = atan2(pos[1]-_yCentre,pos[0]-_xCentre); + if (phi <0.) phi += _const_2pi; + phi = phi - _phiAtPCA + _phi0; + momentum[0] = _pxy*cos(phi); + momentum[1] = _pxy*sin(phi); + momentum[2] = _momentum[2]; + + +} diff --git a/Utilities/DataHelper/src/LineClass.cc b/Utilities/DataHelper/src/LineClass.cc new file mode 100644 index 00000000..826a0d8f --- /dev/null +++ b/Utilities/DataHelper/src/LineClass.cc @@ -0,0 +1,73 @@ +#include "DataHelper/LineClass.h" +#include <math.h> + +/* + * Constructor + */ + +LineClass::LineClass(float x0, + float y0, + float z0, + float ax, + float ay, + float az) { + _x0[0] = x0; + _x0[1] = y0; + _x0[2] = z0; + _ax[0] = ax; + _ax[1] = ay; + _ax[2] = az; +} + +LineClass::LineClass(float *x0, + float *ax) { + + for (int i=0; i<3; ++i) { + _x0[i] = x0[i]; + _ax[i] = ax[i]; + } +} + +float * LineClass::getReferencePoint() { + return _x0; +} + +float * LineClass::getDirectionalVector() { + return _ax; +} + +void LineClass::setReferencePoint(float *x0) { + for (int i=0; i<3; ++i) + _x0[i] = x0[i]; +} + +void LineClass::setDirectionalVector(float *ax) { + for (int i=0; i<3; ++i) + _ax[i] = ax[i]; + +} + +float LineClass::getDistanceToPoint(float * xpoint, float * pos) { + + float dif[3]; + float prod = 0; + float den = 0; + for (int i=0; i<3; ++i) { + dif[i] = xpoint[i] - _x0[i]; + prod += _ax[i]*dif[i]; + den += _ax[i]*_ax[i]; + } + float time = prod/fmax(1e-10,den); + + float dist = 0.0; + for (int i=0; i<3; ++i) { + pos[i] = _x0[i] + _ax[i]*time; + dist += (xpoint[i]-pos[i])*(xpoint[i]-pos[i]); + } + dist = sqrt(dist); + + return dist; + + + +} diff --git a/Utilities/DataHelper/src/LineClass.cc~ b/Utilities/DataHelper/src/LineClass.cc~ new file mode 100644 index 00000000..e23db1fa --- /dev/null +++ b/Utilities/DataHelper/src/LineClass.cc~ @@ -0,0 +1,73 @@ +#include "LineClass.h" +#include <math.h> + +/* + * Constructor + */ + +LineClass::LineClass(float x0, + float y0, + float z0, + float ax, + float ay, + float az) { + _x0[0] = x0; + _x0[1] = y0; + _x0[2] = z0; + _ax[0] = ax; + _ax[1] = ay; + _ax[2] = az; +} + +LineClass::LineClass(float *x0, + float *ax) { + + for (int i=0; i<3; ++i) { + _x0[i] = x0[i]; + _ax[i] = ax[i]; + } +} + +float * LineClass::getReferencePoint() { + return _x0; +} + +float * LineClass::getDirectionalVector() { + return _ax; +} + +void LineClass::setReferencePoint(float *x0) { + for (int i=0; i<3; ++i) + _x0[i] = x0[i]; +} + +void LineClass::setDirectionalVector(float *ax) { + for (int i=0; i<3; ++i) + _ax[i] = ax[i]; + +} + +float LineClass::getDistanceToPoint(float * xpoint, float * pos) { + + float dif[3]; + float prod = 0; + float den = 0; + for (int i=0; i<3; ++i) { + dif[i] = xpoint[i] - _x0[i]; + prod += _ax[i]*dif[i]; + den += _ax[i]*_ax[i]; + } + float time = prod/fmax(1e-10,den); + + float dist = 0.0; + for (int i=0; i<3; ++i) { + pos[i] = _x0[i] + _ax[i]*time; + dist += (xpoint[i]-pos[i])*(xpoint[i]-pos[i]); + } + dist = sqrt(dist); + + return dist; + + + +} diff --git a/Utilities/DataHelper/src/TrackExtended.cc b/Utilities/DataHelper/src/TrackExtended.cc new file mode 100644 index 00000000..bd5444d7 --- /dev/null +++ b/Utilities/DataHelper/src/TrackExtended.cc @@ -0,0 +1,229 @@ +#include "DataHelper/ClusterExtended.h" +#include "DataHelper/TrackerHitExtended.h" +#include "DataHelper/TrackExtended.h" +#include <math.h> + +TrackExtended::TrackExtended( ) { + _track = NULL; + _superCluster = NULL; + _trackerHitVector.clear(); + _clusterVec.clear(); + _group = NULL; +} + +TrackExtended::TrackExtended( Track * track) { + _track = track; + _superCluster = NULL; + _trackerHitVector.clear(); + _clusterVec.clear(); + _group = NULL; + +} + +TrackExtended::TrackExtended( TrackerHitExtended * trackerhit) { + _trackerHitVector.clear(); + _trackerHitVector.push_back(trackerhit); + _track = NULL; + _superCluster = NULL; + _clusterVec.clear(); + _group = NULL; +} + + +TrackExtended::~TrackExtended() {} + +Track * TrackExtended::getTrack() { + return _track; +} + +const float * TrackExtended::getSeedPosition() { + return _seedPosition; +} + +const float * TrackExtended::getSeedDirection() { + return _seedDirection; +} + +ClusterExtendedVec & TrackExtended::getClusterVec() { + return _clusterVec; +} + +ClusterExtended * TrackExtended::getSuperCluster() { + return _superCluster; +} + +TrackerHitExtendedVec & TrackExtended::getTrackerHitExtendedVec() { + return _trackerHitVector; +} + +void TrackExtended::addCluster(ClusterExtended * cluster) { + _clusterVec.push_back(cluster); +} + +void TrackExtended::setSuperCluster(ClusterExtended * superCluster) { + _superCluster = superCluster; +} + +void TrackExtended::setSeedDirection( float * seedDirection ) { + _seedDirection[0] = seedDirection[0]; + _seedDirection[1] = seedDirection[1]; + _seedDirection[2] = seedDirection[2]; +} + +void TrackExtended::setSeedPosition( float * seedPosition) { + _seedPosition[0] = seedPosition[0]; + _seedPosition[1] = seedPosition[1]; + _seedPosition[2] = seedPosition[2]; +} + +void TrackExtended::addTrackerHitExtended( TrackerHitExtended * trackerhit) { + _trackerHitVector.push_back(trackerhit); + +} + +void TrackExtended::ClearTrackerHitExtendedVec() { + _trackerHitVector.clear(); +} + +void TrackExtended::setX0( float x0 ) { + _x0 = x0; +} + +void TrackExtended::setD0( float d0 ) { + _d0 = d0; +} + +void TrackExtended::setZ0( float z0 ) { + _z0 = z0; +} + +void TrackExtended::setY0( float y0 ) { + _y0 = y0; +} + +void TrackExtended::setR0( float r0 ) { + _r0 = r0; +} + +void TrackExtended::setBz( float bz ) { + _bz = bz; +} + +void TrackExtended::setPhi0( float phi0 ) { + _phi0 = phi0; +} + +void TrackExtended::setPhi( float phi ) { + _phi = phi; +} + +void TrackExtended::setOmega( float omega ) { + _omega = omega; +} + +void TrackExtended::setTanLambda( float tanLambda ) { + _tanLambda = tanLambda; +} +void TrackExtended::setNDF( int ndf) { + _ndf = ndf; +} + + +float TrackExtended::getX0() { + return _x0; +} + +float TrackExtended::getY0() { + return _y0; +} + +float TrackExtended::getR0() { + return _r0; +} + +float TrackExtended::getBz() { + return _bz; +} + +float TrackExtended::getD0() { + return _d0; +} + +float TrackExtended::getZ0() { + return _z0; +} + +float TrackExtended::getPhi0() { + return _phi0; +} + +float TrackExtended::getPhi() { + return _phi; +} + +float TrackExtended::getOmega() { + return _omega; +} + +float TrackExtended::getTanLambda() { + return _tanLambda; +} + +int TrackExtended::getNDF() { + return _ndf; +} + +void TrackExtended::setStart(float * xStart) { + _xStart[0] = xStart[0]; + _xStart[1] = xStart[1]; + _xStart[2] = xStart[2]; +} + +void TrackExtended::setEnd(float * xEnd) { + _xEnd[0] = xEnd[0]; + _xEnd[1] = xEnd[1]; + _xEnd[2] = xEnd[2]; +} + +float * TrackExtended::getStart() { + return _xStart; +} + +float * TrackExtended::getEnd() { + return _xEnd; +} + +void TrackExtended::setGroupTracks( GroupTracks * group ) { + _group = group; +} + +GroupTracks * TrackExtended::getGroupTracks() { + return _group; +} + +float TrackExtended::getChi2() { + return _chi2; +} + +void TrackExtended::setChi2(float chi2) { + _chi2 = chi2; +} + +void TrackExtended::setCovMatrix(float * cov) { + for (int i=0; i<15; ++i) + _cov[i] = cov[i]; +} + +float * TrackExtended::getCovMatrix() { + return _cov; +} + +// void TrackExtended::setType(int type) { +// _type = type + +// } + +// int TrackExtended::getType() { +// return _type; + +// } diff --git a/Utilities/DataHelper/src/TrackExtended.cc~ b/Utilities/DataHelper/src/TrackExtended.cc~ new file mode 100644 index 00000000..77393705 --- /dev/null +++ b/Utilities/DataHelper/src/TrackExtended.cc~ @@ -0,0 +1,229 @@ +#include "ClusterExtended.h" +#include "TrackerHitExtended.h" +#include "TrackExtended.h" +#include <math.h> + +TrackExtended::TrackExtended( ) { + _track = NULL; + _superCluster = NULL; + _trackerHitVector.clear(); + _clusterVec.clear(); + _group = NULL; +} + +TrackExtended::TrackExtended( Track * track) { + _track = track; + _superCluster = NULL; + _trackerHitVector.clear(); + _clusterVec.clear(); + _group = NULL; + +} + +TrackExtended::TrackExtended( TrackerHitExtended * trackerhit) { + _trackerHitVector.clear(); + _trackerHitVector.push_back(trackerhit); + _track = NULL; + _superCluster = NULL; + _clusterVec.clear(); + _group = NULL; +} + + +TrackExtended::~TrackExtended() {} + +Track * TrackExtended::getTrack() { + return _track; +} + +const float * TrackExtended::getSeedPosition() { + return _seedPosition; +} + +const float * TrackExtended::getSeedDirection() { + return _seedDirection; +} + +ClusterExtendedVec & TrackExtended::getClusterVec() { + return _clusterVec; +} + +ClusterExtended * TrackExtended::getSuperCluster() { + return _superCluster; +} + +TrackerHitExtendedVec & TrackExtended::getTrackerHitExtendedVec() { + return _trackerHitVector; +} + +void TrackExtended::addCluster(ClusterExtended * cluster) { + _clusterVec.push_back(cluster); +} + +void TrackExtended::setSuperCluster(ClusterExtended * superCluster) { + _superCluster = superCluster; +} + +void TrackExtended::setSeedDirection( float * seedDirection ) { + _seedDirection[0] = seedDirection[0]; + _seedDirection[1] = seedDirection[1]; + _seedDirection[2] = seedDirection[2]; +} + +void TrackExtended::setSeedPosition( float * seedPosition) { + _seedPosition[0] = seedPosition[0]; + _seedPosition[1] = seedPosition[1]; + _seedPosition[2] = seedPosition[2]; +} + +void TrackExtended::addTrackerHitExtended( TrackerHitExtended * trackerhit) { + _trackerHitVector.push_back(trackerhit); + +} + +void TrackExtended::ClearTrackerHitExtendedVec() { + _trackerHitVector.clear(); +} + +void TrackExtended::setX0( float x0 ) { + _x0 = x0; +} + +void TrackExtended::setD0( float d0 ) { + _d0 = d0; +} + +void TrackExtended::setZ0( float z0 ) { + _z0 = z0; +} + +void TrackExtended::setY0( float y0 ) { + _y0 = y0; +} + +void TrackExtended::setR0( float r0 ) { + _r0 = r0; +} + +void TrackExtended::setBz( float bz ) { + _bz = bz; +} + +void TrackExtended::setPhi0( float phi0 ) { + _phi0 = phi0; +} + +void TrackExtended::setPhi( float phi ) { + _phi = phi; +} + +void TrackExtended::setOmega( float omega ) { + _omega = omega; +} + +void TrackExtended::setTanLambda( float tanLambda ) { + _tanLambda = tanLambda; +} +void TrackExtended::setNDF( int ndf) { + _ndf = ndf; +} + + +float TrackExtended::getX0() { + return _x0; +} + +float TrackExtended::getY0() { + return _y0; +} + +float TrackExtended::getR0() { + return _r0; +} + +float TrackExtended::getBz() { + return _bz; +} + +float TrackExtended::getD0() { + return _d0; +} + +float TrackExtended::getZ0() { + return _z0; +} + +float TrackExtended::getPhi0() { + return _phi0; +} + +float TrackExtended::getPhi() { + return _phi; +} + +float TrackExtended::getOmega() { + return _omega; +} + +float TrackExtended::getTanLambda() { + return _tanLambda; +} + +int TrackExtended::getNDF() { + return _ndf; +} + +void TrackExtended::setStart(float * xStart) { + _xStart[0] = xStart[0]; + _xStart[1] = xStart[1]; + _xStart[2] = xStart[2]; +} + +void TrackExtended::setEnd(float * xEnd) { + _xEnd[0] = xEnd[0]; + _xEnd[1] = xEnd[1]; + _xEnd[2] = xEnd[2]; +} + +float * TrackExtended::getStart() { + return _xStart; +} + +float * TrackExtended::getEnd() { + return _xEnd; +} + +void TrackExtended::setGroupTracks( GroupTracks * group ) { + _group = group; +} + +GroupTracks * TrackExtended::getGroupTracks() { + return _group; +} + +float TrackExtended::getChi2() { + return _chi2; +} + +void TrackExtended::setChi2(float chi2) { + _chi2 = chi2; +} + +void TrackExtended::setCovMatrix(float * cov) { + for (int i=0; i<15; ++i) + _cov[i] = cov[i]; +} + +float * TrackExtended::getCovMatrix() { + return _cov; +} + +// void TrackExtended::setType(int type) { +// _type = type + +// } + +// int TrackExtended::getType() { +// return _type; + +// } diff --git a/Utilities/DataHelper/src/TrackExtended.h~ b/Utilities/DataHelper/src/TrackExtended.h~ new file mode 100644 index 00000000..e6cbd4a0 --- /dev/null +++ b/Utilities/DataHelper/src/TrackExtended.h~ @@ -0,0 +1,130 @@ +#ifndef TRACKAR_H +#define TRACKAR_H 1 + +//#include "lcio.h" +//#include "EVENT/LCIO.h" +#include "plcio/Track.h" +#include <vector> +//#include "ClusterExtended.h" +//#include "TrackerHitExtended.h" +#include "GroupTracks.h" + +using namespace plcio; + +class ClusterExtended; +//class GroupTracks; +class TrackerHitExtended; +typedef std::vector<TrackerHitExtended*> TrackerHitExtendedVec; +typedef std::vector<ClusterExtended*> ClusterExtendedVec; + +/** + * Class extending native LCIO class Track. <br> + * Class TrackExtended is used in TrackwiseClustering <br> + * and Wolf processors. <br> + * @author A. Raspereza (DESY)<br> + * @version $Id: TrackExtended.h,v 1.8 2007-09-05 09:39:49 rasp Exp $<br> + */ + +class TrackExtended { + + public: + + + TrackExtended( ); + TrackExtended( TrackerHitExtended * trackerhit ); + TrackExtended( plcio::Track * track ); + ~TrackExtended(); + + plcio::Track * getTrack(); + const float * getSeedDirection(); + const float * getSeedPosition(); + ClusterExtendedVec & getClusterVec(); + ClusterExtended * getSuperCluster(); + TrackerHitExtendedVec & getTrackerHitExtendedVec(); + void addCluster(ClusterExtended * cluster); + void setSuperCluster(ClusterExtended * superCluster); + void setSeedDirection( float * seedDirection ); + void setSeedPosition( float * seedPosition); + void addTrackerHitExtended( TrackerHitExtended * trackerhit); + void ClearTrackerHitExtendedVec(); + + void setX0(float x0); + void setY0(float y0); + void setR0(float r0); + void setD0(float d0); + void setZ0(float z0); + void setBz(float bz); + void setPhi0(float phi0); + void setPhi(float phi); + void setOmega(float omega); + void setTanLambda(float tanLambda); + + + void setStart(float * xStart); + void setEnd(float * xEnd); + + + float getX0(); + float getY0(); + float getD0(); + float getZ0(); + float getOmega(); + float getTanLambda(); + float getPhi(); + float getR0(); + float getBz(); + float getPhi0(); + + float * getStart(); + float * getEnd(); + + void setGroupTracks( GroupTracks * group ); + GroupTracks * getGroupTracks(); + + float getChi2(); + void setChi2(float chi2); + + int getNDF(); + void setNDF(int ndf); + + float * getCovMatrix(); + void setCovMatrix( float * cov); + + private: + + ClusterExtended *_superCluster; + ClusterExtendedVec _clusterVec; + GroupTracks * _group; + Track * _track; + float _seedDirection[3]; + float _seedPosition[3]; + TrackerHitExtendedVec _trackerHitVector; + + float _x0; + float _y0; + float _r0; + float _bz; + float _phi0; + + float _xStart[3]; + float _xEnd[3]; + + float _d0; // d0 in canonical parameterisation + float _z0; // z0 in canonical parameterisation + float _omega; // omega in canonical parameterisation + float _tanLambda; // tanlambda in canonical parameterisation + float _phi; // phi in canonical parameterisation + + float _chi2; // chi2 of the fit + + float _cov[15]; // covariance matrix + + int _ndf; // NDF + + int _type; // Track type; + +}; + +typedef std::vector<TrackExtended*> TrackExtendedVec; + +#endif diff --git a/Utilities/DataHelper/src/TrackerHitExtended.cc b/Utilities/DataHelper/src/TrackerHitExtended.cc new file mode 100644 index 00000000..d2169bf3 --- /dev/null +++ b/Utilities/DataHelper/src/TrackerHitExtended.cc @@ -0,0 +1,133 @@ +#include "DataHelper/TrackExtended.h" +#include "DataHelper/TrackerHitExtended.h" + +TrackerHitExtended::TrackerHitExtended(const edm4hep::TrackerHit& trackerhit): + _trackerHit(trackerhit){ + _trackAR = NULL; + _trackVecAR.clear(); + _usedInFit = false; +} + +TrackerHitExtended::~TrackerHitExtended() { + +} + +void TrackerHitExtended::setTrackExtended(TrackExtended * trackAR) { + _trackAR = trackAR; + _trackVecAR.clear(); + _trackVecAR.push_back(trackAR); + +} + +void TrackerHitExtended::addTrackExtended(TrackExtended * trackAR) { + _trackVecAR.push_back(trackAR); +} + +TrackExtendedVec & TrackerHitExtended::getTrackExtendedVec() { + return _trackVecAR; +} + + +void TrackerHitExtended::setTrackerHitTo(TrackerHitExtended * hitTo) { + _hitTo = hitTo; +} + +void TrackerHitExtended::setTrackerHitFrom(TrackerHitExtended * hitFrom) { + _hitFrom = hitFrom; + +} + +void TrackerHitExtended::setGenericDistance(float genericDistance) { + _genericDistance = genericDistance; +} + +//void TrackerHitExtended::setTrackerHit(edm4hep::TrackerHit * hit) { +// _trackerHit = hit; +//} + +void TrackerHitExtended::setYresTo(float yresTo) { + _yresTo = yresTo; +} + +void TrackerHitExtended::setYresFrom(float yresFrom) { + _yresFrom = yresFrom; +} + +void TrackerHitExtended::setDirVec(float * dirVec) { + _dirVec[0] = dirVec[0]; + _dirVec[1] = dirVec[1]; + _dirVec[2] = dirVec[2]; +} + +void TrackerHitExtended::setType(int ityp) { + _type = ityp; +} + +void TrackerHitExtended::setDet(int idet) { + _idet = idet; +} + +edm4hep::TrackerHit * TrackerHitExtended::getTrackerHit() { + return &_trackerHit; +} + +TrackExtended * TrackerHitExtended::getTrackExtended() { + return _trackAR; +} + +TrackerHitExtended * TrackerHitExtended::getTrackerHitFrom() { + return _hitFrom; +} +TrackerHitExtended * TrackerHitExtended::getTrackerHitTo() { + return _hitTo; +} +float TrackerHitExtended::getGenericDistance() { + return _genericDistance; +} +float TrackerHitExtended::getYresTo() { + return _yresTo; +} +float TrackerHitExtended::getYresFrom() { + return _yresFrom; +} + +float * TrackerHitExtended::getDirVec() { + return _dirVec; +} + +void TrackerHitExtended::clearTrackVec() { + _trackVecAR.clear(); +} + +float TrackerHitExtended::getResolutionRPhi() { + return _rphiReso; +} + +float TrackerHitExtended::getResolutionZ() { + return _zReso; +} + +void TrackerHitExtended::setResolutionRPhi(float rphiReso) { + _rphiReso = rphiReso; +} + +void TrackerHitExtended::setResolutionZ(float zReso) { + _zReso = zReso; +} + +int TrackerHitExtended::getType() { + return _type; +} + +int TrackerHitExtended::getDet() { + return _idet; +} + +void TrackerHitExtended::setUsedInFit(bool usedInFit) { + _usedInFit = usedInFit; +} + +bool TrackerHitExtended::getUsedInFit() { + return _usedInFit; +} + diff --git a/Utilities/DataHelper/src/TrackerHitExtended.cc~ b/Utilities/DataHelper/src/TrackerHitExtended.cc~ new file mode 100644 index 00000000..551b1923 --- /dev/null +++ b/Utilities/DataHelper/src/TrackerHitExtended.cc~ @@ -0,0 +1,133 @@ +#include "TrackExtended.h" +#include "TrackerHitExtended.h" + +TrackerHitExtended::TrackerHitExtended(const edm4hep::TrackerHit& trackerhit): + _trackerHit(trackerhit){ + _trackAR = NULL; + _trackVecAR.clear(); + _usedInFit = false; +} + +TrackerHitExtended::~TrackerHitExtended() { + +} + +void TrackerHitExtended::setTrackExtended(TrackExtended * trackAR) { + _trackAR = trackAR; + _trackVecAR.clear(); + _trackVecAR.push_back(trackAR); + +} + +void TrackerHitExtended::addTrackExtended(TrackExtended * trackAR) { + _trackVecAR.push_back(trackAR); +} + +TrackExtendedVec & TrackerHitExtended::getTrackExtendedVec() { + return _trackVecAR; +} + + +void TrackerHitExtended::setTrackerHitTo(TrackerHitExtended * hitTo) { + _hitTo = hitTo; +} + +void TrackerHitExtended::setTrackerHitFrom(TrackerHitExtended * hitFrom) { + _hitFrom = hitFrom; + +} + +void TrackerHitExtended::setGenericDistance(float genericDistance) { + _genericDistance = genericDistance; +} + +//void TrackerHitExtended::setTrackerHit(edm4hep::TrackerHit * hit) { +// _trackerHit = hit; +//} + +void TrackerHitExtended::setYresTo(float yresTo) { + _yresTo = yresTo; +} + +void TrackerHitExtended::setYresFrom(float yresFrom) { + _yresFrom = yresFrom; +} + +void TrackerHitExtended::setDirVec(float * dirVec) { + _dirVec[0] = dirVec[0]; + _dirVec[1] = dirVec[1]; + _dirVec[2] = dirVec[2]; +} + +void TrackerHitExtended::setType(int ityp) { + _type = ityp; +} + +void TrackerHitExtended::setDet(int idet) { + _idet = idet; +} + +edm4hep::TrackerHit * TrackerHitExtended::getTrackerHit() { + return &_trackerHit; +} + +TrackExtended * TrackerHitExtended::getTrackExtended() { + return _trackAR; +} + +TrackerHitExtended * TrackerHitExtended::getTrackerHitFrom() { + return _hitFrom; +} +TrackerHitExtended * TrackerHitExtended::getTrackerHitTo() { + return _hitTo; +} +float TrackerHitExtended::getGenericDistance() { + return _genericDistance; +} +float TrackerHitExtended::getYresTo() { + return _yresTo; +} +float TrackerHitExtended::getYresFrom() { + return _yresFrom; +} + +float * TrackerHitExtended::getDirVec() { + return _dirVec; +} + +void TrackerHitExtended::clearTrackVec() { + _trackVecAR.clear(); +} + +float TrackerHitExtended::getResolutionRPhi() { + return _rphiReso; +} + +float TrackerHitExtended::getResolutionZ() { + return _zReso; +} + +void TrackerHitExtended::setResolutionRPhi(float rphiReso) { + _rphiReso = rphiReso; +} + +void TrackerHitExtended::setResolutionZ(float zReso) { + _zReso = zReso; +} + +int TrackerHitExtended::getType() { + return _type; +} + +int TrackerHitExtended::getDet() { + return _idet; +} + +void TrackerHitExtended::setUsedInFit(bool usedInFit) { + _usedInFit = usedInFit; +} + +bool TrackerHitExtended::getUsedInFit() { + return _usedInFit; +} + diff --git a/Utilities/DataHelper/src/TrackerHitExtended.h~ b/Utilities/DataHelper/src/TrackerHitExtended.h~ new file mode 100644 index 00000000..e2a85261 --- /dev/null +++ b/Utilities/DataHelper/src/TrackerHitExtended.h~ @@ -0,0 +1,84 @@ +#ifndef TRACKERHITAR_H +#define TRACKERHITAR_H 1 + +//#include "lcio.h" +//#include "EVENT/LCIO.h" +#include "edm4hep/TrackerHit.h" +//#include "TrackExtended.h" +#include <vector> + + +//using namespace edm4hep; + +class TrackExtended; +typedef std::vector<TrackExtended*> TrackExtendedVec; + +/** + * Class extending native LCIO class TrackerHit. <br> + * Class TrackerHitExtended is used in TrackwiseClustering <br> + * and Wolf processors. <br> + * @author A. Raspereza (DESY)<br> + * @version $Id: TrackerHitExtended.h,v 1.7 2007-09-05 09:39:49 rasp Exp $<br> + */ +class TrackerHitExtended { + + public: + + TrackerHitExtended(const edm4hep::TrackerHit& trackerhit); + ~TrackerHitExtended(); + void setTrackExtended(TrackExtended * trackAR); + void addTrackExtended(TrackExtended * trackAR); + void setTrackerHitTo(TrackerHitExtended * hitTo); + void setTrackerHitFrom(TrackerHitExtended * hitFrom); + void setGenericDistance(float genericDistance); + void setTrackerHit(edm4hep::TrackerHit hit); + void setYresTo(float yresTo); + void setYresFrom(float yresFrom); + void setDirVec(float * dirVec); + void clearTrackVec(); + void setResolutionRPhi(float rphiReso); + void setResolutionZ(float zReso); + void setType(int type); + void setDet(int idet); + void setUsedInFit(bool usedInFit); + + edm4hep::TrackerHit* getTrackerHit(); + TrackExtended * getTrackExtended(); + TrackExtendedVec & getTrackExtendedVec(); + TrackerHitExtended * getTrackerHitFrom(); + TrackerHitExtended * getTrackerHitTo(); + float getGenericDistance(); + float getYresTo(); + float getYresFrom(); + float * getDirVec(); + float getResolutionRPhi(); + float getResolutionZ(); + int getType(); + int getDet(); + bool getUsedInFit(); + + private: + + edm4hep::TrackerHit _trackerHit; + TrackExtended * _trackAR; + TrackerHitExtended * _hitTo; + TrackerHitExtended * _hitFrom; + TrackExtendedVec _trackVecAR; + + float _rphiReso; + float _zReso; + float _yresTo; + float _yresFrom; + float _genericDistance; + float _dirVec[3]; + + int _type; + int _idet; + + bool _usedInFit; + +}; + +typedef std::vector<TrackerHitExtended*> TrackerHitExtendedVec; + +#endif -- GitLab