From 4125224a4a0ad176a5ea4095d8de62c93a2b2d66 Mon Sep 17 00:00:00 2001
From: Chengdong Fu <fucd@ihep.ac.cn>
Date: Wed, 3 Feb 2021 10:52:49 +0800
Subject: [PATCH] compatible tool for tracker hit

---
 Utilities/DataHelper/CMakeLists.txt           |   1 +
 .../include/DataHelper/TrackerHitHelper.h     |  13 ++
 .../DataHelper/include/Identifier/CEPCConf.h  |  49 ++++++++
 Utilities/DataHelper/src/TrackerHitHelper.cpp | 116 ++++++++++++++++++
 4 files changed, 179 insertions(+)
 create mode 100644 Utilities/DataHelper/include/DataHelper/TrackerHitHelper.h
 create mode 100644 Utilities/DataHelper/include/Identifier/CEPCConf.h
 create mode 100644 Utilities/DataHelper/src/TrackerHitHelper.cpp

diff --git a/Utilities/DataHelper/CMakeLists.txt b/Utilities/DataHelper/CMakeLists.txt
index 0583b24b..00c42ddd 100644
--- a/Utilities/DataHelper/CMakeLists.txt
+++ b/Utilities/DataHelper/CMakeLists.txt
@@ -15,6 +15,7 @@ gaudi_add_library(DataHelperLib
                           src/TrackerHitExtended.cc
                           src/TrackExtended.cc
                           src/TrackHitPair.cc
+			  src/TrackerHitHelper.cpp
                   LINK EDM4HEP::edm4hep
                        EDM4HEP::edm4hepDict
                        ${GSL_LIBRARIES}
diff --git a/Utilities/DataHelper/include/DataHelper/TrackerHitHelper.h b/Utilities/DataHelper/include/DataHelper/TrackerHitHelper.h
new file mode 100644
index 00000000..c83fb7a4
--- /dev/null
+++ b/Utilities/DataHelper/include/DataHelper/TrackerHitHelper.h
@@ -0,0 +1,13 @@
+#ifndef TrackerHitHelper_H
+#define TrackerHitHelper_H
+
+#include "edm4hep/TrackerHit.h"
+#include <array>
+
+namespace CEPC{
+  std::array<float, 6> GetCovMatrix(edm4hep::TrackerHit& hit);
+  float                GetResolutionRPhi(edm4hep::TrackerHit& hit);
+  float                GetResolutionZ(edm4hep::TrackerHit& hit);
+}
+
+#endif
diff --git a/Utilities/DataHelper/include/Identifier/CEPCConf.h b/Utilities/DataHelper/include/Identifier/CEPCConf.h
new file mode 100644
index 00000000..2dceac7e
--- /dev/null
+++ b/Utilities/DataHelper/include/Identifier/CEPCConf.h
@@ -0,0 +1,49 @@
+/* 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
diff --git a/Utilities/DataHelper/src/TrackerHitHelper.cpp b/Utilities/DataHelper/src/TrackerHitHelper.cpp
new file mode 100644
index 00000000..640be272
--- /dev/null
+++ b/Utilities/DataHelper/src/TrackerHitHelper.cpp
@@ -0,0 +1,116 @@
+#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){
+  if(hit.isAvailable()){
+    int type = hit.getType();
+    std::cout << CEPCConf::TrkHitTypeBit::COMPOSITE_SPACEPOINT << " " << CEPCConf::TrkHitTypeBit::PLANAR << " " << std::endl;
+    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);
+    
+#ifndef MethodUsedInSpacePointBuilder
+      std::cout << "======================" << std::endl;
+      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];
+#endif
+      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.;
+}
-- 
GitLab