Inclination angle effects in MCH digitizer (#13038)

* first prototype for inclination and impact on resolution

Author: mwinn2

Detectors/MUON/MCH/Simulation/include/MCHSimulation/Hit.h

@@ -25,19 +25,25 @@ class Hit : public ::o2::BasicXYZEHit<float>
 
  public:
   Hit(int trackId = 0, short detElemId = 0, math_utils::Point3D<float> entrancePoint = {}, const math_utils::Point3D<float> exitPoint = {},
-      float eloss = 0.0, float length = 0.0, float tof = 0.0)
-    : ::o2::BasicXYZEHit<float>(entrancePoint.x(), entrancePoint.y(), entrancePoint.z(), tof, eloss, trackId, detElemId), mLength{length}, mExitPoint(exitPoint)
+      float eloss = 0.0, float length = 0.0, float entranceTof = 0.0, float eTot = 0.0)
+    : ::o2::BasicXYZEHit<float>(entrancePoint.x(), entrancePoint.y(), entrancePoint.z(), entranceTof, eloss, trackId, detElemId), mLength{length}, mExitPoint(exitPoint), mEtot{eTot}
   {
   }
 
   math_utils::Point3D<float> entrancePoint() const { return GetPos(); }
   math_utils::Point3D<float> exitPoint() const { return mExitPoint; }
 
+  // time in s
+  float entranceTof() const { return GetTime(); }
+  // particle energy in GeV
+  float eTot() const { return mEtot; }
+
   short detElemId() const { return GetDetectorID(); }
 
  private:
   float mLength = {};
   math_utils::Point3D<float> mExitPoint = {};
+  float mEtot = {};
   ClassDefNV(Hit, 1);
 };
 

Detectors/MUON/MCH/Simulation/include/MCHSimulation/Response.h

@@ -32,11 +32,12 @@ class Response
  public:
   Response(Station station);
   ~Response() = default;
-
   float getChargeSpread() const { return mChargeSpread; }
   float getPitch() const { return mPitch; }
   float getSigmaIntegration() const { return mSigmaIntegration; }
   bool isAboveThreshold(float charge) const { return charge > mChargeThreshold; }
+  bool isAngleEffect() const { return mAngleEffect; }
+  bool isMagnetEffect() const { return mMagnetEffect; }
 
   /** Converts energy deposition into a charge.
    *
@@ -65,6 +66,20 @@ class Response
   /// compute the number of samples corresponding to the charge in ADC units
   uint32_t nSamples(float charge) const;
 
+  /// compute deteriation of y-resolution due to track inclination and B-field
+  float inclandbfield(float thetawire, float betagamma, float bx) const;
+
+ private:
+  MathiesonOriginal mMathieson{}; ///< Mathieson function
+  float mPitch = 0.f;             ///< anode-cathode pitch (cm)
+  float mChargeSlope = 0.f;       ///< charge slope used in E to charge conversion
+  float mChargeSpread = 0.f;      ///< width of the charge distribution (cm)
+  float mSigmaIntegration = 0.f;  ///< number of sigmas used for charge distribution
+  float mChargeCorr = 0.f;        ///< amplitude of charge correlation between cathodes
+  float mChargeThreshold = 0.f;   ///< minimum fraction of charge considered
+  bool mAngleEffect = true;       ///< switch for angle effect influencing charge deposition
+  bool mMagnetEffect = true;      ///< switch for magnetic field influencing charge deposition
+
   /// Ratio of particle mean eloss with respect MIP's Khalil Boudjemline, sep 2003, PhD.Thesis and Particle Data Book
   float eLossRatio(float logbetagamma) const;
   /// ToDo: check Aliroot formula vs PDG, if really log_10 and not ln or bug in Aliroot
@@ -74,22 +89,11 @@ class Response
 
   /// Angle effect: Normalisation form theta=10 degres to theta between 0 and 10 (Khalil BOUDJEMLINE sep 2003 Ph.D Thesis)
   /// Angle with respect to the wires assuming that chambers are perpendicular to the z axis.
-  float angleEffectNorma(float elossratio) const;
+  float angleEffectNorma(float angle) const;
 
   /// Magnetic field effect: Normalisation form theta=16 degres (eq. 10 degrees B=0) to theta between -20 and 20 (Lamia Benhabib jun 2006 )
   /// Angle with respect to the wires assuming that chambers are perpendicular to the z axis.
   float magAngleEffectNorma(float angle, float bfield) const;
-
- private:
-  MathiesonOriginal mMathieson{}; ///< Mathieson function
-  float mPitch = 0.f;             ///< anode-cathode pitch (cm)
-  float mChargeSlope = 0.f;       ///< charge slope used in E to charge conversion
-  float mChargeSpread = 0.f;      ///< width of the charge distribution (cm)
-  float mSigmaIntegration = 0.f;  ///< number of sigmas used for charge distribution
-  float mChargeCorr = 0.f;        ///< amplitude of charge correlation between cathodes
-  float mChargeThreshold = 0.f;   ///< minimum fraction of charge considered
-  bool mAngleEffect = true;       ///< switch for angle effect influencing charge deposition
-  bool mMagnetEffect = true;      ///< switch for magnetic field influencing charge deposition
 };
 } // namespace mch
 } // namespace o2

Detectors/MUON/MCH/Simulation/src/DEDigitizer.cxx

@@ -17,6 +17,8 @@
 
 #include "DetectorsRaw/HBFUtils.h"
 #include "MCHSimulation/DigitizerParam.h"
+#include <TGeoGlobalMagField.h>
+#include "Field/MagneticField.h"
 
 /// Convert collision time to ROF time (ROF duration = 4 BC)
 std::pair<o2::InteractionRecord, uint8_t> time2ROFtime(const o2::InteractionRecord& time)
@@ -25,9 +27,11 @@ std::pair<o2::InteractionRecord, uint8_t> time2ROFtime(const o2::InteractionReco
   return std::make_pair(o2::InteractionRecord(time.bc - bc, time.orbit), bc);
 }
 
+//_________________________________________________________________________________________________
+
 namespace o2::mch
 {
-
+//_________________________________________________________________________________________________
 DEDigitizer::DEDigitizer(int deId, math_utils::Transform3D transformation, std::mt19937& random)
   : mDeId{deId},
     mResponse{deId < 300 ? Station::Type1 : Station::Type2345},
@@ -78,6 +82,26 @@ void DEDigitizer::processHit(const Hit& hit, const InteractionRecord& collisionT
   auto localX = mResponse.getAnod(lpos.X());
   auto localY = lpos.Y();
 
+  // calculate angle between track and wire assuming wire perpendicular to z-axis
+  // take global coordinates to avoid issues with rotations of detection elements
+  // neglect rotation of chambers w.r.t. beam
+  auto thetawire = atan((exitPoint.Y() - entrancePoint.Y()) / (entrancePoint.Z() - exitPoint.Z()));
+
+  // local b-field
+  double b[3] = {0., 0., 0.};
+  double x[3] = {entrancePoint.X(), entrancePoint.Y(), entrancePoint.Z()};
+  if (TGeoGlobalMagField::Instance()->GetField()) {
+    TGeoGlobalMagField::Instance()->Field(x, b);
+  } else {
+    LOG(fatal) << "no b field in MCH DEDigitizer";
+  }
+
+  // calculate track betagamma
+  // assume beta = 1 and mass of charged muon to calculate track betagamma from particle energy
+  auto betagamma = hit.eTot() / 0.1056583745;
+  auto yAngleEffect = mResponse.inclandbfield(thetawire, betagamma, b[0]);
+  localY += yAngleEffect;
+
   // borders of charge integration area
   auto dxy = mResponse.getSigmaIntegration() * mResponse.getChargeSpread();
   auto xMin = localX - dxy;

Detectors/MUON/MCH/Simulation/src/Response.cxx

@@ -47,7 +47,6 @@ Response::Response(Station station)
   mChargeCorr = ResponseParam::Instance().chargeCorrelation;
   mChargeThreshold = ResponseParam::Instance().chargeThreshold;
 }
-
 //_____________________________________________________________________
 float Response::etocharge(float edepos) const
 {
@@ -90,29 +89,66 @@ uint32_t Response::nSamples(float charge) const
   return std::round(std::pow(charge / signalParam[1], 1. / signalParam[2]) + signalParam[0]);
 }
 //_____________________________________________________________________
+float Response::inclandbfield(float thetawire, float betagamma, float bx) const
+{
+  float yAngleEffect = 0;
+  // auxiliary variables for b-field and inclination angle effect
+  float eLossParticleElossMip = 0.0;
+  float sigmaEffect10degrees = 0.0;
+  float sigmaEffectThetadegrees = 0.0;
+
+  if (isAngleEffect()) {
+    if (!isMagnetEffect()) {
+      thetawire = abs(thetawire);
+      if ((betagamma > 3.2) && (thetawire * TMath::RadToDeg() <= 15.)) {
+        betagamma = log(betagamma); // check if ln or log10
+        eLossParticleElossMip = eLossRatio(betagamma);
+        sigmaEffect10degrees = angleEffect10(eLossParticleElossMip);
+        sigmaEffectThetadegrees = sigmaEffect10degrees / angleEffectNorma(thetawire * TMath::RadToDeg());
+        if (o2::mch::Station() == o2::mch::Station::Type1) {
+          sigmaEffectThetadegrees /= 1.09833 + 0.017 * (thetawire * TMath::RadToDeg());
+        }
+        yAngleEffect = 0.0001 * gRandom->Gaus(0, sigmaEffectThetadegrees); // error due to the angle effect in cm
+      }
+    } else {
+      if ((betagamma > 3.2) && (abs(thetawire * TMath::RadToDeg()) <= 15.)) {
+        betagamma = log(betagamma);
+        eLossParticleElossMip = eLossRatio(betagamma);
+        sigmaEffect10degrees = angleEffect10(eLossParticleElossMip);
+        sigmaEffectThetadegrees = sigmaEffect10degrees / magAngleEffectNorma(thetawire * TMath::RadToDeg(), bx / 10.); // check b-field unit in aliroot and O2
+        if (o2::mch::Station() == o2::mch::Station::Type1) {
+          sigmaEffectThetadegrees /= 1.09833 + 0.017 * (thetawire * TMath::RadToDeg());
+        }
+        yAngleEffect = 0.0001 * gRandom->Gaus(0, sigmaEffectThetadegrees);
+      }
+    }
+  }
+  return yAngleEffect;
+}
+//_____________________________________________________________________
 float Response::eLossRatio(float logbetagamma) const
 {
   // Ratio of particle mean eloss with respect MIP's Khalil Boudjemline, sep 2003, PhD.Thesis and Particle Data Book
   /// copied from aliroot AliMUONv1.cxx
   float eLossRatioParam[5] = {1.02138, -9.54149e-02, +7.83433e-02, -9.98208e-03, +3.83279e-04};
-  return eLossRatioParam[0] + eLossRatioParam[1] * logbetagamma + eLossRatioParam[2] * std::pow(logbetagamma, 2) + eLossRatioParam[3] * std::pow(logbetagamma, 3) + eLossRatioParam[4] * std::pow(logbetagamma, 4);
+  return eLossRatioParam[0] + eLossRatioParam[1] * logbetagamma + eLossRatioParam[2] * logbetagamma * logbetagamma + eLossRatioParam[3] * logbetagamma * logbetagamma * logbetagamma + eLossRatioParam[4] * logbetagamma * logbetagamma * logbetagamma * logbetagamma;
 }
 //_____________________________________________________________________
 float Response::angleEffect10(float elossratio) const
 {
   /// Angle effect in tracking chambers at theta =10 degres as a function of ElossRatio (Khalil BOUDJEMLINE sep 2003 Ph.D Thesis) (in micrometers)
   /// copied from aliroot AliMUONv1.cxx
   float angleEffectParam[3] = {1.90691e+02, -6.62258e+01, 1.28247e+01};
-  return angleEffectParam[0] + angleEffectParam[1] * elossratio + angleEffectParam[2] * std::pow(elossratio, 2);
+  return angleEffectParam[0] + angleEffectParam[1] * elossratio + angleEffectParam[2] * elossratio * elossratio;
 }
 //_____________________________________________________________________
-float Response::angleEffectNorma(float elossratio) const
+float Response::angleEffectNorma(float angle) const
 {
   /// Angle effect: Normalisation form theta=10 degres to theta between 0 and 10 (Khalil BOUDJEMLINE sep 2003 Ph.D Thesis)
   /// Angle with respect to the wires assuming that chambers are perpendicular to the z axis.
   /// copied from aliroot AliMUONv1.cxx
   float angleEffectParam[4] = {4.148, -6.809e-01, 5.151e-02, -1.490e-03};
-  return angleEffectParam[0] + angleEffectParam[1] * elossratio + angleEffectParam[2] * std::pow(elossratio, 2) + angleEffectParam[3] * std::pow(elossratio, 3);
+  return angleEffectParam[0] + angleEffectParam[1] * angle + angleEffectParam[2] * angle * angle + angleEffectParam[3] * angle * angle * angle;
 }
 //_____________________________________________________________________
 float Response::magAngleEffectNorma(float angle, float bfield) const
@@ -122,5 +158,5 @@ float Response::magAngleEffectNorma(float angle, float bfield) const
   /// copied from aliroot AliMUONv1.cxx
   float angleEffectParam[7] = {8.6995, 25.4022, 13.8822, 2.4717, 1.1551, -0.0624, 0.0012};
   float aux = std::abs(angle - angleEffectParam[0] * bfield);
-  return 121.24 / ((angleEffectParam[1] + angleEffectParam[2] * std::abs(bfield)) + angleEffectParam[3] * aux + angleEffectParam[4] * std::pow(aux, 2) + angleEffectParam[5] * std::pow(aux, 3) + angleEffectParam[6] * std::pow(aux, 4));
+  return 121.24 / ((angleEffectParam[1] + angleEffectParam[2] * std::abs(bfield)) + angleEffectParam[3] * aux + angleEffectParam[4] * aux * aux + angleEffectParam[5] * aux * aux * aux + angleEffectParam[6] * aux * aux * aux * aux);
 }

Detectors/MUON/MCH/Simulation/src/Stepper.cxx

@@ -59,6 +59,7 @@ void Stepper::process(const TVirtualMC& vmc)
   if (t.isEntering()) {
     float x, y, z;
     vmc.TrackPosition(x, y, z);
+    mTof = (float)vmc.TrackTime();
     mEntrancePoint.SetXYZ(x, y, z);
     resetStep();
   }
@@ -68,9 +69,11 @@ void Stepper::process(const TVirtualMC& vmc)
 
   if (t.isExiting() || t.isStopped()) {
     float x, y, z;
+    float etot;
     vmc.TrackPosition(x, y, z);
+    etot = (float)vmc.Etot();
     mHits->emplace_back(stack->GetCurrentTrackNumber(), detElemId, mEntrancePoint,
-                        math_utils::Point3D<float>{x, y, z}, mTrackEloss, mTrackLength);
+                        math_utils::Point3D<float>{x, y, z}, mTrackEloss, mTrackLength, mTof, etot);
     resetStep();
   }
 }

Detectors/MUON/MCH/Simulation/src/Stepper.h

@@ -43,6 +43,7 @@ class Stepper
  private:
   float mTrackEloss{0.0};
   float mTrackLength{0.0};
+  float mTof{0.0};
   std::vector<o2::mch::Hit>* mHits{nullptr};
   math_utils::Point3D<float> mEntrancePoint;
 };

Detectors/MUON/MCH/Simulation/test/testDigitizer.cxx

@@ -28,7 +28,9 @@
 #include "MCHSimulation/Hit.h"
 #include "SimulationDataFormat/MCCompLabel.h"
 #include "SimulationDataFormat/MCTruthContainer.h"
+#include "Field/MagneticField.h"
 #include "TGeoManager.h"
+#include "TGeoGlobalMagField.h"
 #include "boost/format.hpp"
 #include <boost/test/data/test_case.hpp>
 #include <boost/property_tree/ptree.hpp>
@@ -50,6 +52,20 @@ struct GEOMETRY {
   }
 };
 
+void initField(float l3Current, float dipoleCurrent)
+{
+  /// Create the magnetic field map if not already done
+  if (TGeoGlobalMagField::Instance()->GetField()) {
+    return;
+  }
+
+  auto field =
+    o2::field::MagneticField::createFieldMap(l3Current, dipoleCurrent, o2::field::MagneticField::kConvLHC, false, 3500.,
+                                             "A-A", "$(O2_ROOT)/share/Common/maps/mfchebKGI_sym.root");
+  TGeoGlobalMagField::Instance()->SetField(field);
+  TGeoGlobalMagField::Instance()->Lock();
+}
+
 namespace
 {
 short detElemId1 = 101;
@@ -100,6 +116,8 @@ BOOST_AUTO_TEST_CASE(DigitizerTest)
   o2::conf::ConfigurableParam::setValue("MCHDigitizer", "seed", 123);
   o2::mch::Digitizer digitizer(transformation);
 
+  initField(-30000.0, -6000.0);
+
   int trackId1 = 0;
   int trackId2 = 1;
   int trackId3 = 2;
@@ -108,13 +126,14 @@ BOOST_AUTO_TEST_CASE(DigitizerTest)
   float eloss = 1.e-6;
   float length = 0.f;
   float tof = 0.f;
+  float energ = 10.0;
 
   std::vector<o2::mch::Hit> hits1(2);
-  hits1.at(0) = o2::mch::Hit(trackId1, detElemId1, entrancePoint1, exitPoint1, eloss, length, tof);
-  hits1.at(1) = o2::mch::Hit(trackId2, detElemId2, entrancePoint2, exitPoint2, eloss, length, tof);
+  hits1.at(0) = o2::mch::Hit(trackId1, detElemId1, entrancePoint1, exitPoint1, eloss, length, tof, energ);
+  hits1.at(1) = o2::mch::Hit(trackId2, detElemId2, entrancePoint2, exitPoint2, eloss, length, tof, energ);
 
   std::vector<o2::mch::Hit> hits2(1);
-  hits2.at(0) = o2::mch::Hit(trackId3, detElemId3, entrancePoint3, exitPoint3, eloss, length, tof);
+  hits2.at(0) = o2::mch::Hit(trackId3, detElemId3, entrancePoint3, exitPoint3, eloss, length, tof, energ);
 
   digitizer.processHits(hits1, collisionTime1, 0, 0);
   digitizer.processHits(hits2, collisionTime2, 0, 0);

Steer/DigitizerWorkflow/src/MCHDigitizerSpec.cxx

@@ -45,7 +45,7 @@ namespace mch
 class MCHDPLDigitizerTask : public o2::base::BaseDPLDigitizer
 {
  public:
-  MCHDPLDigitizerTask() : o2::base::BaseDPLDigitizer(o2::base::InitServices::GEOM) {}
+  MCHDPLDigitizerTask() : o2::base::BaseDPLDigitizer(o2::base::InitServices::FIELD | o2::base::InitServices::GEOM) {}
 
   void initDigitizerTask(framework::InitContext& ic) override
   {