Merge branch 'swp31326' into 'master'

Manual sweep of !31326 (NSW Calib Tool)

See merge request !31980

MuonSpectrometer/MuonCalib/NSWCalib/NSWCalibTools/NSWCalibTools/INSWCalibTool.h

@@ -6,12 +6,29 @@
 
 #include "GaudiKernel/IAlgTool.h"
 #include "GeoPrimitives/GeoPrimitives.h"
+#include "Identifier/Identifier.h"
 
 #include <cmath>
 #include <vector>
 
 static const InterfaceID IID_INSWCalibTool("Muon::INSWCalibTool",1,0);
 
+namespace NSWCalib { 
+
+  struct CalibratedStrip {
+    double charge;
+    double time;
+    double resTime;
+    double distDrift;
+    double resTransDistDrift;
+    double resLongDistDrift;
+    double dx;
+    Identifier identifier;
+  };
+
+}
+
+
 namespace Muon {
 
   class MM_RawData;
@@ -23,8 +40,8 @@ namespace Muon {
     static const InterfaceID& interfaceID()  {return IID_INSWCalibTool;}
 
   public:  // interface methods
-
-    virtual StatusCode calibrate(const Muon::MM_RawData* mmRawData, const Amg::Vector3D& globalPos, double& dist_drift, double& distRes_drift, double& calib_charge) = 0;
+    
+    virtual StatusCode calibrate(const Muon::MM_RawData* mmRawData, const Amg::Vector3D& globalPos, NSWCalib::CalibratedStrip& calibStrip) const = 0;
 
   };
   

MuonSpectrometer/MuonCalib/NSWCalib/NSWCalibTools/src/NSWCalibTool.cxx

@@ -4,9 +4,12 @@
 
 #include "NSWCalibTool.h"
 #include "GaudiKernel/SystemOfUnits.h"
+#include "GaudiKernel/PhysicalConstants.h"
 
 namespace {
   static constexpr double const& toRad = M_PI/180;
+  static constexpr double const& pitchErr = 0.425 * 0.425 / 12;
+  static constexpr double const& reciprocalSpeedOfLight = 1. / Gaudi::Units::c_light; // mm/ns
 }
 
 Muon::NSWCalibTool::NSWCalibTool(const std::string& t,
@@ -18,9 +21,12 @@ Muon::NSWCalibTool::NSWCalibTool(const std::string& t,
   declareInterface<INSWCalibTool>(this);
 
   declareProperty("MagFieldSvc",   m_magFieldSvc, "Magnetic Field Service");
-  declareProperty("DriftVelocity", m_vDrift = 47., "Drift Velocity");
+  declareProperty("DriftVelocity", m_vDrift = 0.047, "Drift Velocity");
   declareProperty("TimeResolution", m_timeRes = 25., "Time Resolution");
-
+  declareProperty("longDiff",m_longDiff=0.019); //mm/mm
+  declareProperty("transDiff",m_transDiff=0.036); //mm/mm
+  declareProperty("ionUncertainty",m_ionUncertainty=4.0); //ns
+  declareProperty("timeOffset", m_timeOffset = -100); //ns                       
 }
 
 
@@ -42,29 +48,27 @@ StatusCode Muon::NSWCalibTool::initialize()
 
 }
 
-StatusCode Muon::NSWCalibTool::calibrate( const Muon::MM_RawData* mmRawData, const Amg::Vector3D& globalPos, double& dist_drift, double& distRes_drift, double& calib_charge)
+StatusCode Muon::NSWCalibTool::calibrate( const Muon::MM_RawData* mmRawData, const Amg::Vector3D& globalPos, NSWCalib::CalibratedStrip& calibStrip) const
 {
-  double charge, time, vDriftCorrected;
-  charge = mmRawData->charge();
-  time = mmRawData->time();
+  double  vDriftCorrected;
+  calibStrip.charge = mmRawData->charge();
+  calibStrip.time = mmRawData->time() - globalPos.norm() * reciprocalSpeedOfLight + m_timeOffset;
+
+  calibStrip.identifier = mmRawData->identify();
 
   Amg::Vector3D magneticField;
-  calib_charge = charge;
 
   m_magFieldSvc->getField(&globalPos,&magneticField);
-  m_mmBFieldX = magneticField.x();
-  m_mmBFieldY = magneticField.y();
-  m_mmBFieldZ = magneticField.z();
   
   double lorentzAngle = m_lorentzAngleFunction->Eval((magneticField.y() > 0. ? 1. : -1.) * std::fabs (magneticField.y()) ) * toRad; // in radians;
 
   vDriftCorrected = m_vDrift * std::cos(lorentzAngle);
-  dist_drift = vDriftCorrected * time;
-  dist_drift =  dist_drift * Gaudi::Units::perThousand;
-  distRes_drift = vDriftCorrected * m_timeRes;
-  distRes_drift = distRes_drift * Gaudi::Units::perThousand;
+  calibStrip.distDrift = vDriftCorrected * calibStrip.time;
+  calibStrip.resTransDistDrift = pitchErr + std::pow(m_transDiff * calibStrip.distDrift, 2); 
+  calibStrip.resLongDistDrift = std::pow(m_ionUncertainty * vDriftCorrected, 2) + std::pow(m_longDiff * calibStrip.distDrift, 2);
+  calibStrip.dx = std::sin(lorentzAngle) * calibStrip.time * m_vDrift;
 
-  return StatusCode::SUCCESS;  
+  return StatusCode::SUCCESS;
 }
 
 

MuonSpectrometer/MuonCalib/NSWCalib/NSWCalibTools/src/NSWCalibTool.h

@@ -33,7 +33,8 @@ namespace Muon {
 
     virtual ~NSWCalibTool() = default;
     
-    StatusCode calibrate( const Muon::MM_RawData* mmRawData, const Amg::Vector3D& globalPos, double& dist_drift, double& distRes_drift, double& calib_charge);
+    StatusCode calibrate( const Muon::MM_RawData* mmRawData, const Amg::Vector3D& globalPos, NSWCalib::CalibratedStrip& calibStrip) const;
+
     virtual StatusCode initialize();
     virtual StatusCode finalize();
 
@@ -45,13 +46,13 @@ namespace Muon {
    
 
     TF1* m_lorentzAngleFunction;
- 
-    float m_mmBFieldX;
-    float m_mmBFieldY;
-    float m_mmBFieldZ;
-    
+     
     float m_vDrift;
     float m_timeRes;
+    float m_longDiff;
+    float m_transDiff;
+    float m_ionUncertainty;
+    double m_timeOffset;
   };
   
 }

MuonSpectrometer/MuonCnv/MuonMM_CnvTools/src/MmRdoToPrepDataToolCore.cxx

@@ -137,11 +137,8 @@ StatusCode Muon::MmRdoToPrepDataToolCore::processCollection(const MM_RawDataColl
       ATH_MSG_WARNING("Could not get the global strip position for MM");
       continue;
     }
-    double dist_drift, distRes_drift, charge_calib;
-    ATH_CHECK (m_calibTool->calibrate(rdo, globalPos, dist_drift, distRes_drift, charge_calib));
-
-    const int time = rdo->time();
-    const int charge = charge_calib;
+    NSWCalib::CalibratedStrip calibStrip;
+    ATH_CHECK (m_calibTool->calibrate(rdo, globalPos, calibStrip));
 
     const Amg::Vector3D globalDir(globalPos.x(), globalPos.y(), globalPos.z());
     Trk::LocalDirection localDir;
@@ -149,7 +146,7 @@ StatusCode Muon::MmRdoToPrepDataToolCore::processCollection(const MM_RawDataColl
     Amg::Vector2D lpos;
     psurf.globalToLocal(globalPos,globalPos,lpos);
     psurf.globalToLocalDirection(globalDir, localDir);
-    float inAngle_XZ = fabs( localDir.angleXZ() / CLHEP::degree);
+    float inAngle_XZ = std::abs( localDir.angleXZ() / CLHEP::degree);
    
     ATH_MSG_DEBUG(" Surface centre x " << psurf.center().x() << " y " << psurf.center().y() << " z " << psurf.center().z() );
     ATH_MSG_DEBUG(" localPos x " << localPos.x() << " localPos y " << localPos.y() << " lpos recalculated 0 " << lpos[0] << " lpos y " << lpos[1]);
@@ -159,26 +156,28 @@ StatusCode Muon::MmRdoToPrepDataToolCore::processCollection(const MM_RawDataColl
 
     // get for now a temporary error matrix -> to be fixed
     double resolution = 0.07;
-    if (fabs(inAngle_XZ)>3) resolution = ( -.001/3.*fabs(inAngle_XZ) ) + .28/3.;
+    if (std::abs(inAngle_XZ)>3) resolution = ( -.001/3.*std::abs(inAngle_XZ) ) + .28/3.;
     double errX = 0;
     const MuonGM::MuonChannelDesign* design = detEl->getDesign(layid);
     if( !design ){
       ATH_MSG_WARNING("Failed to get design for " << m_idHelperSvc->toString(layid) );
     }else{
-      errX = fabs(design->inputPitch)/sqrt(12);
+      errX = std::abs(design->inputPitch)/std::sqrt(12);
       ATH_MSG_DEBUG(" strips inputPitch " << design->inputPitch << " error " << errX);
     }
 // add strip width to error
-    resolution = sqrt(resolution*resolution+errX*errX);
+    resolution = std::sqrt(resolution*resolution+errX*errX);
 
-    Amg::MatrixX* cov = new Amg::MatrixX(1,1);
+    Amg::MatrixX* cov = new Amg::MatrixX(2,2);
     cov->setIdentity();
-    (*cov)(0,0) = resolution*resolution;  
+    (*cov)(0,0) = calibStrip.resTransDistDrift;  
+    (*cov)(1,1) = calibStrip.resLongDistDrift;
+    localPos.x() += calibStrip.dx;
 
     if(!merge) {
-      prdColl->push_back(new MMPrepData(prdId,hash,localPos,rdoList,cov,detEl,time,charge));
+      prdColl->push_back(new MMPrepData(prdId, hash, localPos, rdoList, cov, detEl, calibStrip.time, calibStrip.charge, calibStrip.distDrift));
     } else {
-       MMPrepData mpd = MMPrepData(prdId,hash,localPos,rdoList,cov,detEl,time,charge);
+       MMPrepData mpd = MMPrepData(prdId, hash, localPos, rdoList, cov, detEl, calibStrip.time, calibStrip.charge, calibStrip.distDrift);
        // set the hash of the MMPrepData such that it contains the correct value in case it gets used in SimpleMMClusterBuilderTool::getClusters
        mpd.setHashAndIndex(hash,0);
        MMprds.push_back(mpd);