Merge branch 'debug_tan_theta_y' into 'master'

debug tan(theta_y) peak

See merge request faser/calypso!389

PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx

@@ -335,6 +335,18 @@ StatusCode NtupleDumperAlg::initialize()
   m_tree->Branch("t_st3_x", &m_t_st_x[3]);
   m_tree->Branch("t_st3_y", &m_t_st_y[3]);
   m_tree->Branch("t_st3_z", &m_t_st_z[3]);
+  m_tree->Branch("t_st0_px", &m_t_st_px[0]);
+  m_tree->Branch("t_st0_py", &m_t_st_py[0]);
+  m_tree->Branch("t_st0_pz", &m_t_st_pz[0]);
+  m_tree->Branch("t_st1_px", &m_t_st_px[1]);
+  m_tree->Branch("t_st1_py", &m_t_st_py[1]);
+  m_tree->Branch("t_st1_pz", &m_t_st_pz[1]);
+  m_tree->Branch("t_st2_px", &m_t_st_px[2]);
+  m_tree->Branch("t_st2_py", &m_t_st_py[2]);
+  m_tree->Branch("t_st2_pz", &m_t_st_pz[2]);
+  m_tree->Branch("t_st3_px", &m_t_st_px[3]);
+  m_tree->Branch("t_st3_py", &m_t_st_py[3]);
+  m_tree->Branch("t_st3_pz", &m_t_st_pz[3]);
   m_tree->Branch("isFiducial", &m_isFiducial);
 
   m_tree->Branch("truthParticleBarcode", &m_truthParticleBarcode);
@@ -1054,6 +1066,14 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
           m_t_st_y[station].push_back(positions[station].y());
           m_t_st_z[station].push_back(positions[station].z());
         }
+
+        auto momenta = m_fiducialParticleTool->getTruthMomenta(truthParticle->barcode());
+        for (int station = 0; station < 4; ++station) {
+          m_t_st_px[station].push_back(momenta[station].x());
+          m_t_st_py[station].push_back(momenta[station].y());
+          m_t_st_pz[station].push_back(momenta[station].z());
+        }
+
         if (truthParticle->hasProdVtx()) {
           m_t_prodVtx_x.push_back(truthParticle->prodVtx()->x());
           m_t_prodVtx_y.push_back(truthParticle->prodVtx()->y());
@@ -1550,6 +1570,11 @@ NtupleDumperAlg::clearTree() const
     m_t_st_y[station].clear();
     m_t_st_z[station].clear();
   }
+  for (int station = 0; station < 4; ++station) {
+    m_t_st_px[station].clear();
+    m_t_st_py[station].clear();
+    m_t_st_pz[station].clear();
+  }
   m_truthParticleBarcode.clear();
   m_truthParticleMatchedTracks.clear();
   m_truthParticleIsFiducial.clear();
@@ -1634,6 +1659,11 @@ void NtupleDumperAlg::clearTrackTruth() const {
     m_t_st_y[station].push_back(999999);
     m_t_st_z[station].push_back(999999);
   }
+  for (int station = 0; station < 4; ++station) {
+    m_t_st_px[station].push_back(999999);
+    m_t_st_py[station].push_back(999999);
+    m_t_st_pz[station].push_back(999999);
+  }
   m_isFiducial.push_back(false);
 }
 

PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h

@@ -317,6 +317,9 @@ private:
   mutable std::array<std::vector<double>, 4> m_t_st_x; // vector of the x components of the simulated hits of the truth particle for each station
   mutable std::array<std::vector<double>, 4> m_t_st_y; // vector of the y components of the simulated hits of the truth particle for each station
   mutable std::array<std::vector<double>, 4> m_t_st_z; // vector of the z components of the simulated hits of the truth particle for each station
+  mutable std::array<std::vector<double>, 4> m_t_st_px; // x components of the true momentum at each station
+  mutable std::array<std::vector<double>, 4> m_t_st_py; // y components of the true momentum at each station
+  mutable std::array<std::vector<double>, 4> m_t_st_pz; // z components of the true momentum at each station
   mutable std::vector<bool> m_isFiducial; // track is fiducial if there are simulated hits for stations 1 - 3 and the distance from the center is smaller than 100 mm
 
   mutable std::vector<int> m_truthParticleBarcode; // vector of barcodes of all truth particles with a momentum larger 50 GeV

Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/IFiducialParticleTool.h

@@ -23,6 +23,12 @@ public:
    */
   virtual std::array<HepGeom::Point3D<double>, 4>
   getTruthPositions(int barcode) const = 0;
+
+  /** Return average truth momentum in each station.
+   * @param barcode of a xAOD::TruthParticle
+   */
+  virtual std::array<HepGeom::Point3D<double>, 4>
+  getTruthMomenta(int barcode) const = 0;
 };
 
 #endif // FASERACTSKALMANFILTER_IFIDUCIALEVENTSELECTIONTOOL_H

Tracking/Acts/FaserActsKalmanFilter/src/FiducialParticleTool.cxx

@@ -6,8 +6,7 @@
 
 constexpr float NaN = std::numeric_limits<double>::quiet_NaN();
 
-FiducialParticleTool::FiducialParticleTool(const std::string &type,
-                                           const std::string &name,
+FiducialParticleTool::FiducialParticleTool(const std::string &type, const std::string &name,
                                            const IInterface *parent)
     : base_class(type, name, parent) {}
 
@@ -32,23 +31,32 @@ bool FiducialParticleTool::isFiducial(int barcode) const {
   return isFiducial;
 }
 
-HepGeom::Point3D<double>
-FiducialParticleTool::getGlobalPosition(const FaserSiHit &hit) const {
-  Identifier waferId =
-      m_sctHelper->wafer_id(hit.getStation(), hit.getPlane(), hit.getRow(),
-                            hit.getModule(), hit.getSensor());
+HepGeom::Point3D<double> FiducialParticleTool::getGlobalPosition(const FaserSiHit &hit) const {
+  Identifier waferId = m_sctHelper->wafer_id(hit.getStation(), hit.getPlane(), hit.getRow(),
+                                             hit.getModule(), hit.getSensor());
   auto localStartPos = hit.localStartPosition();
   auto localEndPos = hit.localEndPosition();
   HepGeom::Point3D<double> localPos = 0.5 * (localEndPos + localStartPos);
-  const TrackerDD::SiDetectorElement *element =
-      m_detMgr->getDetectorElement(waferId);
-  auto globalPosition =
-      Amg::EigenTransformToCLHEP(element->transformHit()) * localPos;
+  const TrackerDD::SiDetectorElement *element = m_detMgr->getDetectorElement(waferId);
+  auto globalPosition = Amg::EigenTransformToCLHEP(element->transformHit()) * localPos;
   return globalPosition;
 }
 
-std::array<HepGeom::Point3D<double>, 4>
-FiducialParticleTool::getTruthPositions(int barcode) const {
+HepGeom::Point3D<double> FiducialParticleTool::getMomentum(const FaserSiHit &hit) const {
+  Identifier waferId = m_sctHelper->wafer_id(hit.getStation(), hit.getPlane(), hit.getRow(),
+                                             hit.getModule(), hit.getSensor());
+  auto localStartPos = hit.localStartPosition();
+  auto localEndPos = hit.localEndPosition();
+  const TrackerDD::SiDetectorElement *element = m_detMgr->getDetectorElement(waferId);
+  auto globalStartPosition = Amg::EigenTransformToCLHEP(element->transformHit()) * localStartPos;
+  auto globalEndPosition = Amg::EigenTransformToCLHEP(element->transformHit()) * localEndPos;
+  double p = hit.particleLink()->momentum().rho();
+  HepGeom::Point3D<double> globalDirection = globalEndPosition - globalStartPosition;
+  HepGeom::Point3D<double> momentum = globalDirection / globalDirection.mag() * p;
+  return momentum;
+}
+
+std::array<HepGeom::Point3D<double>, 4> FiducialParticleTool::getTruthPositions(int barcode) const {
   // initialize positions as NaN
   std::array<HepGeom::Point3D<double>, 4> positions{};
   for (auto &station : positions) {
@@ -89,3 +97,46 @@ FiducialParticleTool::getTruthPositions(int barcode) const {
   }
   return positions;
 }
+
+std::array<HepGeom::Point3D<double>, 4> FiducialParticleTool::getTruthMomenta(int barcode) const {
+  // initialize positions as NaN
+  std::array<HepGeom::Point3D<double>, 4> momenta{};
+  for (auto &station : momenta) {
+    station.setX(NaN);
+    station.setY(NaN);
+    station.setZ(NaN);
+  }
+
+  // get simulated hits
+  SG::ReadHandle<FaserSiHitCollection> siHitCollection(m_siHitCollectionKey);
+  if (!siHitCollection.isValid()) {
+    ATH_MSG_WARNING("FaserSiHitCollection not valid.");
+    return momenta;
+  }
+
+  // create map with truth momenta in each station
+  std::array<std::vector<HepGeom::Point3D<double>>, 4> hitMap{};
+  for (const FaserSiHit &hit : *siHitCollection) {
+    if (hit.trackNumber() == barcode) {
+      HepGeom::Point3D<double> momentum = getMomentum(hit);
+      hitMap[hit.getStation()].push_back(momentum);
+    }
+  }
+
+  // calculate average momentum in each station
+  for (int station = 0; station < 4; ++station) {
+    std::vector<HepGeom::Point3D<double>> &hits{hitMap[station]};
+    if (hits.empty()) {
+      continue;
+    } else {
+      // calculate average momentum of all FaserSiHits in a station
+      auto const count = static_cast<double>(hits.size());
+      HepGeom::Point3D<double> sums{};
+      for (const HepGeom::Point3D<double> &hit : hits) {
+        sums += hit;
+      }
+      momenta[station] = sums / count;
+    }
+  }
+  return momenta;
+}

Tracking/Acts/FaserActsKalmanFilter/src/FiducialParticleTool.h

@@ -23,6 +23,10 @@ public:
 
   std::array<HepGeom::Point3D<double>, 4>
   getTruthPositions(int barcode) const override;
+  std::array<HepGeom::Point3D<double>, 4>
+  getTruthMomenta(int barcode) const override;
+  HepGeom::Point3D<double>
+  getMomentum(const FaserSiHit &hit) const;
 
 private:
   HepGeom::Point3D<double> getGlobalPosition(const FaserSiHit &hit) const;