From fe87f25966b8ccdcf9de407a33fbb20eb765cbd9 Mon Sep 17 00:00:00 2001
From: Chengdong Fu <fucd@ihep.ac.cn>
Date: Fri, 19 Feb 2021 15:57:59 +0800
Subject: [PATCH] add optional control for cov matrix
---
Digitisers/SimpleDigi/src/PlanarDigiAlg.cpp | 36 ++++---
Digitisers/SimpleDigi/src/PlanarDigiAlg.h | 3 +
.../include/DataHelper/TrackerHitHelper.h | 1 +
Utilities/DataHelper/src/TrackerHitHelper.cpp | 100 +++++++++---------
4 files changed, 77 insertions(+), 63 deletions(-)
diff --git a/Digitisers/SimpleDigi/src/PlanarDigiAlg.cpp b/Digitisers/SimpleDigi/src/PlanarDigiAlg.cpp
index 482082a0..4a5b1f31 100644
--- a/Digitisers/SimpleDigi/src/PlanarDigiAlg.cpp
+++ b/Digitisers/SimpleDigi/src/PlanarDigiAlg.cpp
@@ -1,5 +1,6 @@
/* -*- Mode: C++; tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 2 -*- */
#include "PlanarDigiAlg.h"
+#include "DataHelper/TrackerHitHelper.h"
#include "GearSvc/IGearSvc.h"
#include "EventSeeder/IEventSeeder.h"
#include "TrackSystemSvc/ITrackSystemSvc.h"
@@ -16,7 +17,7 @@
#include "UTIL/CellIDEncoder.h"
#include <UTIL/Operators.h>
*/
-
+#include "Identifier/CEPCConf.h"
#include "UTIL/ILDConf.h"
// STUFF needed for GEAR
@@ -174,7 +175,7 @@ StatusCode PlanarDigiAlg::execute()
int module = encoder[lcio::ILDCellID0::module];
int sensor = encoder[lcio::ILDCellID0::sensor];
- debug() << "Hit = "<< i << " has celId " << celId << endmsg;
+ debug() << "Hit = " << i << " has celId " << celId << endmsg;
debug() << "side = " << side << endmsg;
debug() << "layerNumber = " << layer << endmsg;
debug() << "moduleNumber = " << module << endmsg;
@@ -307,15 +308,16 @@ StatusCode PlanarDigiAlg::execute()
debug() << " U[0] = "<< u_direction[0] << " U[1] = "<< u_direction[1]
<< " V[0] = "<< v_direction[0] << " V[1] = "<< v_direction[1]
<< endmsg ;
- // fucd
- std::array<float, 6> cov;
- cov[0] = u_direction[0];
- cov[1] = u_direction[1];
- cov[2] = resU;
- cov[3] = v_direction[0];
- cov[4] = v_direction[1];
- cov[5] = resV;
- trkHit.setCovMatrix(cov);
+ // fucd: next TODO: cov[0] = resU*reU, cov[2] = resV*resV, cov[5] = 0
+ if(_usePlanarTag){
+ std::array<float, 6> cov;
+ cov[0] = u_direction[0];
+ cov[1] = u_direction[1];
+ cov[2] = resU;
+ cov[3] = v_direction[0];
+ cov[4] = v_direction[1];
+ cov[5] = resV;
+ trkHit.setCovMatrix(cov);
/* zoujh: TODO - generate TrackerHitPlane with podio
trkHit->setU( u_direction ) ;
trkHit->setV( v_direction ) ;
@@ -325,11 +327,17 @@ StatusCode PlanarDigiAlg::execute()
if( _isStrip ) trkHit->setdV( 0 ); // no error in v direction for strip hits as there is no meesurement information in v direction
else trkHit->setdV( resV ) ;
*/
- trkHit.setType(8);
- if( _isStrip || (resU!=0&&resV==0) ){
- trkHit.setType( UTIL::set_bit( trkHit.getType() , UTIL::ILDTrkHitTypeBit::ONE_DIMENSIONAL ) ) ;
+ std::bitset<32> type;
+ type.set(CEPCConf::TrkHitTypeBit::PLANAR);
+ trkHit.setType((int)type.to_ulong());
+ }
+ else{
+ trkHit.setCovMatrix(CEPC::ConvertToCovXYZ(resU, u_direction[0], u_direction[1], resV, v_direction[0], v_direction[1]));
}
+ if( _isStrip || (resU!=0&&resV==0) ){
+ trkHit.setType( UTIL::set_bit( trkHit.getType() , UTIL::ILDTrkHitTypeBit::ONE_DIMENSIONAL ) ) ;
+ }
trkHit.setEDep( SimTHit.getEDep() );
// make the relation
diff --git a/Digitisers/SimpleDigi/src/PlanarDigiAlg.h b/Digitisers/SimpleDigi/src/PlanarDigiAlg.h
index 668333f2..10ddd9dd 100644
--- a/Digitisers/SimpleDigi/src/PlanarDigiAlg.h
+++ b/Digitisers/SimpleDigi/src/PlanarDigiAlg.h
@@ -79,6 +79,9 @@ protected:
Gaudi::Property<FloatVec> _resV{ this, "ResolutionV", {0.0040} };
// whether hits are 1D strip hits
Gaudi::Property<bool> _isStrip{ this, "IsStrip", false };
+ // whether use Planar tag for type and cov, if true, CEPCConf::TrkHitTypeBit::PLANAR bit is set as true
+ // cov[0]=thetaU, cov[1]=phiU, cov[2]=resU, cov[0]=thetaV, cov[1]=phiV, cov[2]=resV
+ Gaudi::Property<bool> _usePlanarTag{ this, "UsePlanarTag", true };
// Input collections
DataHandle<edm4hep::EventHeaderCollection> _headerCol{"EventHeaderCol", Gaudi::DataHandle::Reader, this};
diff --git a/Utilities/DataHelper/include/DataHelper/TrackerHitHelper.h b/Utilities/DataHelper/include/DataHelper/TrackerHitHelper.h
index 60057ada..fcd9dd14 100644
--- a/Utilities/DataHelper/include/DataHelper/TrackerHitHelper.h
+++ b/Utilities/DataHelper/include/DataHelper/TrackerHitHelper.h
@@ -8,6 +8,7 @@ namespace CEPC{
std::array<float, 6> GetCovMatrix(edm4hep::TrackerHit& hit, bool useSpacePointerBuilderMethod = false);
float GetResolutionRPhi(edm4hep::TrackerHit& hit);
float GetResolutionZ(edm4hep::TrackerHit& hit);
+ std::array<float, 6> ConvertToCovXYZ(float dU, float thetaU, float phiU, float dV, float thetaV, float phiV, bool useSpacePointBuilderMethod = false);
}
#endif
diff --git a/Utilities/DataHelper/src/TrackerHitHelper.cpp b/Utilities/DataHelper/src/TrackerHitHelper.cpp
index 58f5278d..e12eb886 100644
--- a/Utilities/DataHelper/src/TrackerHitHelper.cpp
+++ b/Utilities/DataHelper/src/TrackerHitHelper.cpp
@@ -2,7 +2,6 @@
#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"
@@ -16,60 +15,13 @@ std::array<float,6> CEPC::GetCovMatrix(edm4hep::TrackerHit& hit, bool useSpacePo
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;
+ return ConvertToCovXYZ(dU, thetaU, phiU, dV, thetaV, phiV, useSpacePointBuilderMethod);
}
else{
std::cout << "Warning: not SpacePoint and Planar, return original cov matrix preliminaryly." << std::endl;
@@ -113,3 +65,53 @@ float CEPC::GetResolutionZ(edm4hep::TrackerHit& hit){
}
return 0.;
}
+
+std::array<float, 6> CEPC::ConvertToCovXYZ(float dU, float thetaU, float phiU, float dV, float thetaV, float phiV, bool useSpacePointBuilderMethod){
+ std::array<float,6> cov{0.,0.,0.,0.,0.,0.};
+ 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;
+}
--
GitLab