Merge pull request #141 from fucd/track

REC: compatible tool for tracker hit

Detector/CMakeLists.txt

@@ -6,3 +6,4 @@ add_subdirectory(DetEcalMatrix)
 add_subdirectory(DetInterface)
 add_subdirectory(DetSegmentation)
 add_subdirectory(GeomSvc)
+add_subdirectory(Identifier)

Detector/Identifier/CMakeLists.txt

+################################################################################
+# Package: Identifier
+################################################################################
+
+gaudi_add_header_only_library(Identifier)
+
+install(TARGETS Identifier
+  EXPORT CEPCSWTargets
+  RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" COMPONENT bin
+  LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib
+  COMPONENT dev)

Detector/Identifier/include/Identifier/CEPCConf.h

+/* CEPC conference description ID */
+#ifndef CEPCConf_H
+#define CEPCConf_H
+
+#include <string>
+
+namespace CEPCConf{
+  struct DetID{ // compatible for old codes from Marlin, will convert from compact file, default initial values here
+    static const int NOTUSED = 0;
+    static const int VXD     = 1;
+    static const int SIT     = 2;
+    static const int FTD     = 3;
+    static const int TPC     = 4;
+    static const int SET     = 5;
+    static const int ETD     = 6;
+    static const int DC      = 7; 
+    
+    static const int ECAL        = 20;
+    static const int ECAL_PLUG   = 21;
+    static const int HCAL        = 22;
+    static const int HCAL_RING   = 23;
+    static const int LCAL        = 24;
+    static const int BCAL        = 25;
+    static const int LHCAL       = 26;
+    static const int YOKE        = 27;
+    static const int COIL        = 28;
+    static const int ECAL_ENDCAP = 29;
+    static const int HCAL_ENDCAP = 30;
+    static const int YOKE_ENDCAP = 31;
+    
+    static const int bwd    = -1;
+    static const int barrel =  0;
+    static const int fwd    = +1;
+  };
+  
+  struct TrkHitTypeBit{
+    static const int ONE_DIMENSIONAL      = 29;
+    static const int COMPOSITE_SPACEPOINT = 30;
+    static const int PLANAR               = 3; // 3 is compatible with old tracking codes, 31 or 28 is better in future to modify uniformly
+  };
+  
+  struct TrkHitQualityBit{
+    static const int USED_IN_FIT          = 30;
+    static const int USED_IN_TRACK        = 29;
+    static const int DOUBLE_HIT_CANDIDATE = 28;
+    static const int GOOD                 = 27;
+  };
+}
+#endif

Utilities/DataHelper/CMakeLists.txt

@@ -15,10 +15,12 @@ gaudi_add_library(DataHelperLib
                           src/TrackerHitExtended.cc
                           src/TrackExtended.cc
                           src/TrackHitPair.cc
+			  src/TrackerHitHelper.cpp
                   LINK EDM4HEP::edm4hep
                        EDM4HEP::edm4hepDict
                        ${GSL_LIBRARIES}
                        ${CLHEP_LIBRARIES}
+		       Identifier
 )
 
 install(TARGETS DataHelperLib

Utilities/DataHelper/include/DataHelper/TrackerHitHelper.h

+#ifndef TrackerHitHelper_H
+#define TrackerHitHelper_H
+
+#include "edm4hep/TrackerHit.h"
+#include <array>
+
+namespace CEPC{
+  std::array<float, 6> GetCovMatrix(edm4hep::TrackerHit& hit, bool useSpacePointerBuilderMethod = false);
+  float                GetResolutionRPhi(edm4hep::TrackerHit& hit);
+  float                GetResolutionZ(edm4hep::TrackerHit& hit);
+}
+
+#endif

Utilities/DataHelper/src/TrackerHitHelper.cpp

+#include "DataHelper/TrackerHitHelper.h"
+#include "Identifier/CEPCConf.h"
+
+#include "TMatrixF.h"
+#include "TMatrixFSym.h"
+#include "CLHEP/Matrix/SymMatrix.h"
+#include "CLHEP/Matrix/Matrix.h"
+#include "CLHEP/Vector/ThreeVector.h"
+#include "CLHEP/Vector/Rotation.h"
+#include <bitset>
+
+std::array<float,6> CEPC::GetCovMatrix(edm4hep::TrackerHit& hit, bool useSpacePointBuilderMethod){
+  if(hit.isAvailable()){
+    int type = hit.getType();
+    if(std::bitset<32>(type)[CEPCConf::TrkHitTypeBit::COMPOSITE_SPACEPOINT]){
+      return hit.getCovMatrix();
+    }
+    else if(std::bitset<32>(type)[CEPCConf::TrkHitTypeBit::PLANAR]){
+      std::array<float,6> cov{0.,0.,0.,0.,0.,0.};
+      float thetaU = hit.getCovMatrix(0);
+      float phiU   = hit.getCovMatrix(1);
+      float dU     = hit.getCovMatrix(2);
+      float thetaV = hit.getCovMatrix(3);
+      float phiV   = hit.getCovMatrix(4);
+      float dV     = hit.getCovMatrix(5);
+      
+      if(!useSpacePointBuilderMethod){
+	TMatrixF diffs(2,3);
+	TMatrixF diffsT(3,2);
+	diffs(0,0) = sin(thetaU)*cos(phiU);
+	diffs(0,1) = sin(thetaU)*sin(phiU);
+	diffs(0,2) = cos(thetaU);
+	diffs(1,0) = sin(thetaV)*cos(phiV);
+	diffs(1,1) = sin(thetaV)*sin(phiV);
+	diffs(1,2) = cos(thetaV);
+	
+	diffsT.Transpose(diffs);
+	
+	TMatrixF covMatrixUV(2,2);
+	covMatrixUV(0,0) = dU*dU;
+	covMatrixUV(0,1) = 0;
+	covMatrixUV(1,0) = 0;
+	covMatrixUV(1,1) = dV*dV;
+	
+	TMatrixF covMatrixXYZ(3,3);
+	covMatrixXYZ = diffsT*covMatrixUV*diffs;
+	cov[0] = covMatrixXYZ(0,0);
+	cov[1] = covMatrixXYZ(1,0);
+	cov[2] = covMatrixXYZ(1,1);
+	cov[3] = covMatrixXYZ(2,0);
+	cov[4] = covMatrixXYZ(2,1);
+	cov[5] = covMatrixXYZ(2,2);
+      }
+      else{ // Method used in SpacePointBuilder, results are almost same with above
+	CLHEP::Hep3Vector u_sensor(sin(thetaU)*cos(phiU), sin(thetaU)*sin(phiU), cos(thetaU));
+	CLHEP::Hep3Vector v_sensor(sin(thetaV)*cos(phiV), sin(thetaV)*sin(phiV), cos(thetaV));
+	CLHEP::Hep3Vector w_sensor = u_sensor.cross(v_sensor);
+	CLHEP::HepRotation rot_sensor(u_sensor, v_sensor, w_sensor);
+	CLHEP::HepMatrix rot_sensor_matrix;
+	rot_sensor_matrix = rot_sensor;
+	CLHEP::HepSymMatrix cov_plane(3,0);
+	cov_plane(1,1) = dU*dU;
+	cov_plane(2,2) = dV*dV;
+	CLHEP::HepSymMatrix cov_xyz= cov_plane.similarity(rot_sensor_matrix);
+	cov[0] = cov_xyz[0][0];
+	cov[1] = cov_xyz[1][0];
+	cov[2] = cov_xyz[1][1];
+	cov[3] = cov_xyz[2][0];
+	cov[4] = cov_xyz[2][1];
+	cov[5] = cov_xyz[2][2];
+      }
+      return cov;
+    }
+    else{
+      std::cout << "Warning: not SpacePoint and Planar, return original cov matrix preliminaryly." << std::endl;
+      return hit.getCovMatrix();
+    }
+  }
+  std::array<float,6> cov{0.,0.,0.,0.,0.,0.};
+  return cov;
+}
+
+float CEPC::GetResolutionRPhi(edm4hep::TrackerHit& hit){
+  if(hit.isAvailable()){
+    int type = hit.getType();
+    if(std::bitset<32>(type)[CEPCConf::TrkHitTypeBit::COMPOSITE_SPACEPOINT]){
+      return sqrt(hit.getCovMatrix(0)+hit.getCovMatrix(2));
+    }
+    else if(std::bitset<32>(type)[CEPCConf::TrkHitTypeBit::PLANAR]){
+      return hit.getCovMatrix(2);
+    }
+    else{
+      std::cout << "Warning: not SpacePoint and Planar, return value from original cov matrix preliminaryly." << std::endl;
+      return sqrt(hit.getCovMatrix(0)+hit.getCovMatrix(2));
+    }
+  }
+  return 0.;
+}
+
+float CEPC::GetResolutionZ(edm4hep::TrackerHit& hit){
+  if(hit.isAvailable()){
+    int type = hit.getType();
+    if(std::bitset<32>(type)[CEPCConf::TrkHitTypeBit::COMPOSITE_SPACEPOINT]){
+      return sqrt(hit.getCovMatrix(5));
+    }
+    else if(std::bitset<32>(type)[CEPCConf::TrkHitTypeBit::PLANAR]){
+      return hit.getCovMatrix(5);
+    }
+    else{
+      std::cout << "Warning: not SpacePoint and Planar, return value from original cov matrix preliminaryly." << std::endl;
+      return sqrt(hit.getCovMatrix(5));
+    }
+  }
+  return 0.;
+}