diff --git a/Detector/DetCRD/scripts/CRD-Sim.py b/Detector/DetCRD/scripts/CRD-Sim.py
index 0e734b785f0c6442b320bc75e43c3d10472bc8a5..e658c12772b0d2a8e304055eafae7ad1e5c7a56f 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 5d595621af36d9284d9ac4e1f78adf85c0838f5d..cb3ceeb761e542e1cf6e43c2a469db62a0feac3f 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_v01-SimRec.py b/Detector/DetCRD/scripts/CRD_o1_v01-SimRec.py
index 19836013be8fcc3b1dcbcd6993473de64b1a0ba5..52a6209da65917bb7ced0acacc15d198fee15973 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.TrackerHitCollection = vxdhitname
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.TrackerHitAssociationCollection = "SITTrackerHitAssociation"
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.TrackerHitAssociationCollection = "SETTrackerHitAssociation"
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.TrackerHitAssociationCollection = "FTDTrackerHitAssociation"
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
diff --git a/Detector/DetCRD/scripts/CRD_o1_v02-SimRec.py b/Detector/DetCRD/scripts/CRD_o1_v02-SimRec.py
index 253ea21a71125f20f9216094dd2071e2fdd950e4..05844c9f8d5e2a2e94ae13349a12dda57d0cfcc2 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 @@ digiFTD.UsePlanarTag = True
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 1386affd7623cda211e7d1d8f983fa749b26b9f3..c2d88fdebb76f8b3a70d029979009f2bc2212365 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
diff --git a/Digitisers/SimpleDigi/src/PlanarDigiAlg.cpp b/Digitisers/SimpleDigi/src/PlanarDigiAlg.cpp
index 973a556282a16dccc8933c32cb10307926ec1725..730f61e34df936eab5d54bb4e84fff1ca6866a05 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 0e53e707abcff7188e4c3031261c7163107d5abd..981a171c2a8f6503242f82894db75ff0e8b5735a 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};