From 28fdb29331963bfa789fa2d1b13d64147b795c50 Mon Sep 17 00:00:00 2001
From: Chengdong Fu <fucd@ihep.ac.cn>
Date: Mon, 29 Jan 2024 15:11:20 +0800
Subject: [PATCH] import parameterized spatial resolution model
---
Digitisers/SimpleDigi/src/PlanarDigiAlg.cpp | 84 +++++++++++++++------
Digitisers/SimpleDigi/src/PlanarDigiAlg.h | 3 +
2 files changed, 63 insertions(+), 24 deletions(-)
diff --git a/Digitisers/SimpleDigi/src/PlanarDigiAlg.cpp b/Digitisers/SimpleDigi/src/PlanarDigiAlg.cpp
index 973a5562..730f61e3 100644
--- a/Digitisers/SimpleDigi/src/PlanarDigiAlg.cpp
+++ b/Digitisers/SimpleDigi/src/PlanarDigiAlg.cpp
@@ -54,13 +54,20 @@ PlanarDigiAlg::PlanarDigiAlg(const std::string& name, ISvcLocator* svcLoc)
StatusCode PlanarDigiAlg::initialize()
{
+ if (_parameterize) {
+ if (_parU.size()!=10 || _parV.size()!=10) {
+ fatal() << "parameters number must be 10! now " << _parU.size() << " for U and " << _parV.size() << " for V" << endmsg;
+ return StatusCode::FAILURE;
+ }
+ }
+ else {
+ if (_resU.size() != _resV.size()) {
+ fatal() << "Inconsistent number of resolutions given for U and V coordinate: "
+ << "ResolutionU :" << _resU.size() << " != ResolutionV : " << _resV.size()
+ << endmsg;
- if( _resU.size() != _resV.size() ) {
- fatal() << "Inconsistent number of resolutions given for U and V coordinate: "
- << "ResolutionU :" << _resU.size() << " != ResolutionV : " << _resV.size()
- << endmsg;
-
- return StatusCode::FAILURE;
+ return StatusCode::FAILURE;
+ }
}
// get the GEAR manager
@@ -167,6 +174,7 @@ StatusCode PlanarDigiAlg::execute()
int i = 0;
for( auto SimTHit : *STHcol ) {
if (SimTHit.getEDep()<=_eThreshold) continue;
+ debug() << "MCParticle id " << SimTHit.getMCParticle().id() << endmsg;
const int celId = SimTHit.getCellID() ;
@@ -206,28 +214,60 @@ StatusCode PlanarDigiAlg::execute()
// << endmsg;
// continue;
// }
- if( (_resU.size() > 1 && layer > _resU.size()-1) || (_resV.size() > 1 && layer > _resV.size()-1) ) {
- fatal() << "layer exceeds resolution vector, please check input parameters ResolutionU and ResolutionV" << endmsg;
- return StatusCode::FAILURE;
- }
+ gear::MeasurementSurface const* ms = _GEAR->getMeasurementSurfaceStore().GetMeasurementSurface( encoder.lowWord() );
+ gear::CartesianCoordinateSystem* cartesian = dynamic_cast< gear::CartesianCoordinateSystem* >( ms->getCoordinateSystem() );
+ CLHEP::Hep3Vector uVec = cartesian->getLocalXAxis();
+ CLHEP::Hep3Vector vVec = cartesian->getLocalYAxis();
- float resU = ( _resU.size() > 1 ? _resU.value().at( layer ) : _resU.value().at(0) ) ;
- float resV = ( _resV.size() > 1 ? _resV.value().at( layer ) : _resV.value().at(0) ) ;
+ float u_direction[2] ;
+ u_direction[0] = uVec.theta();
+ u_direction[1] = uVec.phi();
+
+ float v_direction[2] ;
+ v_direction[0] = vVec.theta();
+ v_direction[1] = vVec.phi();
+
+ debug() << " U[0] = "<< u_direction[0] << " U[1] = "<< u_direction[1]
+ << " V[0] = "<< v_direction[0] << " V[1] = "<< v_direction[1]
+ << endmsg ;
+
+ float resU(0), resV(0);
+
+ if (!_parameterize) {
+ if( (_resU.size() > 1 && layer > _resU.size()-1) || (_resV.size() > 1 && layer > _resV.size()-1) ) {
+ fatal() << "layer exceeds resolution vector, please check input parameters ResolutionU and ResolutionV" << endmsg;
+ return StatusCode::FAILURE;
+ }
+
+ resU = ( _resU.size() > 1 ? _resU.value().at( layer ) : _resU.value().at(0) ) ;
+ resV = ( _resV.size() > 1 ? _resV.value().at( layer ) : _resV.value().at(0) ) ;
+ }
+ else { // Riccardo's parameterized model
+ auto mom = SimTHit.getMomentum();
+ CLHEP::Hep3Vector momVec(mom[0], mom[1], mom[2]);
+ const double alpha = uVec.azimAngle(momVec, vVec);
+ const double cotanAlpha = 1./tan(alpha);
+ // TODO: title angle (PI/2), magnetic field (3)
+ const double tanLorentzAngle = (side==0) ? 0. : 0.053 * 3 * cos(M_PI/2.);
+ const double x = fabs(-cotanAlpha - tanLorentzAngle);
+ resU = _parU[0] + _parU[1] * x + _parU[2] * exp(-_parU[9] * x) * cos(_parU[3] * x + _parU[4])
+ + _parU[5] * exp(-0.5 * pow(((x - _parU[6]) / _parU[7]), 2)) + _parU[8] * pow(x, 0.5);
+
+ const double beta = vVec.azimAngle(momVec, uVec);
+ const double cotanBeta = 1./tan(beta);
+ const double y = fabs(-cotanBeta);
+ resV = _parV[0] + _parV[1] * y + _parV[2] * exp(-_parV[9] * y) * cos(_parV[3] * y + _parV[4])
+ + _parV[5] * exp(-0.5 * pow(((y - _parV[6]) / _parV[7]), 2)) + _parV[8] * pow(y, 0.5);
+ }
debug() << " --- will smear hit with resU = " << resU << " and resV = " << resV << endmsg;
auto& pos = SimTHit.getPosition() ;
- //edm4hep::Vector3d smearedPos;
-
- //GearSurfaces::MeasurementSurface* ms = _GEAR->getMeasurementSurfaceStore().GetMeasurementSurface( SimTHit->getCellID0() );
-
- gear::MeasurementSurface const* ms = _GEAR->getMeasurementSurfaceStore().GetMeasurementSurface( encoder.lowWord() );;
CLHEP::Hep3Vector globalPoint(pos[0],pos[1],pos[2]);
CLHEP::Hep3Vector localPoint = ms->getCoordinateSystem()->getLocalPoint(globalPoint);
CLHEP::Hep3Vector localPointSmeared = localPoint;
-
debug() << std::setprecision(8) << "Position of hit before smearing global: ( "<<pos[0]<<" "<<pos[1]<<" "<<pos[2]<< " ) "
<< "local: ( " << localPoint.x() << " " << localPoint.y() << " " << localPoint.z() << " )" << endmsg;
@@ -260,7 +300,6 @@ StatusCode PlanarDigiAlg::execute()
//check if hit is in boundaries
if ( ms->isLocalInBoundary( localPointSmeared ) && fabs(dx)<=_maxPull*resU && fabs(dy)<=_maxPull*resV ) {
- //if ( ms->isLocalInBoundary( localPointSmeared ) ) {
accept_hit = true;
break;
}
@@ -285,10 +324,6 @@ StatusCode PlanarDigiAlg::execute()
<< localPointSmeared.x() << " " << localPointSmeared.y() << " " << localPointSmeared.z() << " )"
<< endmsg;
- //smearedPos[0] = globalPointSmeared.x();
- //smearedPos[1] = globalPointSmeared.y();
- //smearedPos[2] = globalPointSmeared.z();
-
//make the TrackerHitPlaneImpl
auto trkHit = trkhitVec->create();
@@ -296,7 +331,7 @@ StatusCode PlanarDigiAlg::execute()
edm4hep::Vector3d smearedPos(globalPointSmeared.x(), globalPointSmeared.y(), globalPointSmeared.z());
trkHit.setPosition( smearedPos ) ;
-
+ /*
gear::CartesianCoordinateSystem* cartesian = dynamic_cast< gear::CartesianCoordinateSystem* >( ms->getCoordinateSystem() );
CLHEP::Hep3Vector uVec = cartesian->getLocalXAxis();
CLHEP::Hep3Vector vVec = cartesian->getLocalYAxis();
@@ -312,6 +347,7 @@ 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: next TODO: cov[0] = resU*reU, cov[2] = resV*resV, cov[5] = 0
if(_usePlanarTag){
std::array<float, 6> cov;
diff --git a/Digitisers/SimpleDigi/src/PlanarDigiAlg.h b/Digitisers/SimpleDigi/src/PlanarDigiAlg.h
index 0e53e707..981a171c 100644
--- a/Digitisers/SimpleDigi/src/PlanarDigiAlg.h
+++ b/Digitisers/SimpleDigi/src/PlanarDigiAlg.h
@@ -84,6 +84,9 @@ protected:
Gaudi::Property<bool> _usePlanarTag{ this, "UsePlanarTag", true };
Gaudi::Property<float> _eThreshold{ this, "EnergyThreshold", 0 };
Gaudi::Property<float> _maxPull{ this, "PullCutToRetry", 1000. };
+ Gaudi::Property<bool> _parameterize{ this, "ParameterizeResolution", false};
+ Gaudi::Property<FloatVec> _parU{ this, "ParametersU", {0} };
+ Gaudi::Property<FloatVec> _parV{ this, "ParametersV", {0} };
// Input collections
DataHandle<edm4hep::EventHeaderCollection> _headerCol{"EventHeaderCol", Gaudi::DataHandle::Reader, this};
--
GitLab