From 28fdb29331963bfa789fa2d1b13d64147b795c50 Mon Sep 17 00:00:00 2001 From: Chengdong Fu Date: Mon, 29 Jan 2024 15:11:20 +0800 Subject: [PATCH 1/3] 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 973a55628..730f61e34 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: ( "<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 cov; diff --git a/Digitisers/SimpleDigi/src/PlanarDigiAlg.h b/Digitisers/SimpleDigi/src/PlanarDigiAlg.h index 0e53e707a..981a171c2 100644 --- a/Digitisers/SimpleDigi/src/PlanarDigiAlg.h +++ b/Digitisers/SimpleDigi/src/PlanarDigiAlg.h @@ -84,6 +84,9 @@ class PlanarDigiAlg : public GaudiAlgorithm Gaudi::Property _usePlanarTag{ this, "UsePlanarTag", true }; Gaudi::Property _eThreshold{ this, "EnergyThreshold", 0 }; Gaudi::Property _maxPull{ this, "PullCutToRetry", 1000. }; + Gaudi::Property _parameterize{ this, "ParameterizeResolution", false}; + Gaudi::Property _parU{ this, "ParametersU", {0} }; + Gaudi::Property _parV{ this, "ParametersV", {0} }; // Input collections DataHandle _headerCol{"EventHeaderCol", Gaudi::DataHandle::Reader, this}; From eadd8b47515abaaf546aea0d01d3a920aeebcaee Mon Sep 17 00:00:00 2001 From: Chengdong Fu Date: Mon, 29 Jan 2024 15:12:26 +0800 Subject: [PATCH 2/3] update for new LCG --- Detector/DetCRD/scripts/CRD-Sim.py | 2 ++ Detector/DetCRD/scripts/CRD_VXD_MOST2-sim.py | 2 ++ Detector/DetCRD/scripts/CRD_o1_v02-SimRec.py | 6 ++++++ Detector/DetCRD/scripts/Standalone-Sim-RotCrystal.py | 2 ++ 4 files changed, 12 insertions(+) diff --git a/Detector/DetCRD/scripts/CRD-Sim.py b/Detector/DetCRD/scripts/CRD-Sim.py index 0e734b785..e658c1277 100644 --- a/Detector/DetCRD/scripts/CRD-Sim.py +++ b/Detector/DetCRD/scripts/CRD-Sim.py @@ -1,4 +1,6 @@ #!/usr/bin/env python +import os + from Gaudi.Configuration import * from Configurables import k4DataSvc diff --git a/Detector/DetCRD/scripts/CRD_VXD_MOST2-sim.py b/Detector/DetCRD/scripts/CRD_VXD_MOST2-sim.py index 5d595621a..cb3ceeb76 100644 --- a/Detector/DetCRD/scripts/CRD_VXD_MOST2-sim.py +++ b/Detector/DetCRD/scripts/CRD_VXD_MOST2-sim.py @@ -1,4 +1,6 @@ #!/usr/bin/env python +import os + from Gaudi.Configuration import * from Configurables import k4DataSvc diff --git a/Detector/DetCRD/scripts/CRD_o1_v02-SimRec.py b/Detector/DetCRD/scripts/CRD_o1_v02-SimRec.py index 253ea21a7..05844c9f8 100644 --- a/Detector/DetCRD/scripts/CRD_o1_v02-SimRec.py +++ b/Detector/DetCRD/scripts/CRD_o1_v02-SimRec.py @@ -1,6 +1,10 @@ #!/usr/bin/env python +import os + from Gaudi.Configuration import * +NTupleSvc().Output = ["MyTuples DATAFILE='test.root' OPT='NEW' TYP='ROOT'"] + from Configurables import k4DataSvc dsvc = k4DataSvc("EventDataSvc") @@ -160,6 +164,8 @@ from Configurables import DCHDigiAlg digiDC = DCHDigiAlg("DCHDigi") digiDC.DigiDCHitCollection = dchitname +# TODO: DCHDigiAlg need fixed, only WriteAna = True can run +digiDC.WriteAna = True #digiDC.OutputLevel = DEBUG # two strip tracker hits -> one space point diff --git a/Detector/DetCRD/scripts/Standalone-Sim-RotCrystal.py b/Detector/DetCRD/scripts/Standalone-Sim-RotCrystal.py index 1386affd7..c2d88fdeb 100644 --- a/Detector/DetCRD/scripts/Standalone-Sim-RotCrystal.py +++ b/Detector/DetCRD/scripts/Standalone-Sim-RotCrystal.py @@ -1,4 +1,6 @@ #!/usr/bin/env python +import os + from Gaudi.Configuration import * from Configurables import k4DataSvc From 90be1cd8baa61c6ad0a1df2f5e496acf479da908 Mon Sep 17 00:00:00 2001 From: Chengdong Fu Date: Mon, 29 Jan 2024 15:13:29 +0800 Subject: [PATCH 3/3] parameterized resolution in comment --- Detector/DetCRD/scripts/CRD_o1_v01-SimRec.py | 22 ++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/Detector/DetCRD/scripts/CRD_o1_v01-SimRec.py b/Detector/DetCRD/scripts/CRD_o1_v01-SimRec.py index 19836013b..52a6209da 100644 --- a/Detector/DetCRD/scripts/CRD_o1_v01-SimRec.py +++ b/Detector/DetCRD/scripts/CRD_o1_v01-SimRec.py @@ -1,6 +1,10 @@ #!/usr/bin/env python +import os + from Gaudi.Configuration import * +NTupleSvc().Output = ["MyTuples DATAFILE='test.root' OPT='NEW' TYP='ROOT'"] + from Configurables import k4DataSvc dsvc = k4DataSvc("EventDataSvc") @@ -125,6 +129,10 @@ digiVXD.ResolutionU = [0.0028, 0.006, 0.004, 0.004, 0.004, 0.004] digiVXD.ResolutionV = [0.0028, 0.006, 0.004, 0.004, 0.004, 0.004] digiVXD.UsePlanarTag = True +# if Parameterized spatial resolution, set ParameterizeResolution to True +digiVXD.ParameterizeResolution = False +digiVXD.ParametersU = [5.60959e-03, 5.74913e-03, 7.03433e-03, 1.99516, -663.952, 3.752e-03, 0, -0.0704734, 0.0454867e-03, 1.07359] +digiVXD.ParametersV = [5.60959e-03, 5.74913e-03, 7.03433e-03, 1.99516, -663.952, 3.752e-03, 0, -0.0704734, 0.0454867e-03, 1.07359] #digiVXD.OutputLevel = DEBUG digiSIT = PlanarDigiAlg("SITDigi") @@ -135,6 +143,10 @@ digiSIT.ResolutionU = [0.0072] digiSIT.ResolutionV = [0.086] digiSIT.UsePlanarTag = True +# if Parameterized spatial resolution, set ParameterizeResolution to True +digiSIT.ParameterizeResolution = False +digiSIT.ParametersU = [2.29655e-03, 0.965899e-03, 0.584699e-03, 17.0856, 84.566, 12.4695e-03, -0.0643059, 0.168662, 1.87998e-03, 0.514452] +digiSIT.ParametersV = [1.44629e-02, 2.20108e-03, 1.03044e-02, 4.39195e+00, 3.29641e+00, 1.55167e+18, -5.41954e+01, 5.72986e+00, -6.80699e-03, 5.04095e-01] #digiSIT.OutputLevel = DEBUG digiSET = PlanarDigiAlg("SETDigi") @@ -145,6 +157,10 @@ digiSET.ResolutionU = [0.0072] digiSET.ResolutionV = [0.086] digiSET.UsePlanarTag = True +# if Parameterized spatial resolution, set ParameterizeResolution to True +digiSET.ParameterizeResolution = False +digiSET.ParametersU = [2.29655e-03, 0.965899e-03, 0.584699e-03, 17.0856, 84.566, 12.4695e-03, -0.0643059, 0.168662, 1.87998e-03, 0.514452] +digiSET.ParametersV = [1.44629e-02, 2.20108e-03, 1.03044e-02, 4.39195e+00, 3.29641e+00, 1.55167e+18, -5.41954e+01, 5.72986e+00, -6.80699e-03, 5.04095e-01] #digiSET.OutputLevel = DEBUG digiFTD = PlanarDigiAlg("FTDDigi") @@ -155,11 +171,17 @@ digiFTD.ResolutionU = [0.003, 0.003, 0.0072, 0.0072, 0.0072, 0.0072, 0.0072] digiFTD.ResolutionV = [0.003, 0.003, 0.0072, 0.0072, 0.0072, 0.0072, 0.0072] digiFTD.UsePlanarTag = True +# if Parameterized spatial resolution, set ParameterizeResolution to True +digiFTD.ParameterizeResolution = False +digiFTD.ParametersU = [2.29655e-03, 0.965899e-03, 0.584699e-03, 17.0856, 84.566, 12.4695e-03, -0.0643059, 0.168662, 1.87998e-03, 0.514452] +digiFTD.ParametersV = [1.44629e-02, 2.20108e-03, 1.03044e-02, 4.39195e+00, 3.29641e+00, 1.55167e+18, -5.41954e+01, 5.72986e+00, -6.80699e-03, 5.04095e-01] #digiFTD.OutputLevel = DEBUG from Configurables import DCHDigiAlg digiDC = DCHDigiAlg("DCHDigi") digiDC.DigiDCHitCollection = dchitname +# TODO: DCHDigiAlg need fixed, only WriteAna = True can run +digiDC.WriteAna = True #digiDC.OutputLevel = DEBUG # two strip tracker hits -> one space point