Digi: import parameterized model for silicon spatial resolution

Detector/DetCRD/scripts/CRD-Sim.py

@@ -1,4 +1,6 @@
 #!/usr/bin/env python
+import os
+
 from Gaudi.Configuration import *
 
 from Configurables import k4DataSvc

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

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

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

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

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,18 +324,14 @@ 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();
 
     trkHit.setCellID( encoder.lowWord() );
 
     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;

Digitisers/SimpleDigi/src/PlanarDigiAlg.h

@@ -84,6 +84,9 @@ class PlanarDigiAlg : public GaudiAlgorithm
   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};