SensitiveDetector.h

/**************************************************************************
 * basf2 (Belle II Analysis Software Framework)                           *
 * Author: The Belle II Collaboration                                     *
 *                                                                        *
 * See git log for contributors and copyright holders.                    *
 * This file is licensed under LGPL-3.0, see LICENSE.md.                  *
 **************************************************************************/
 
#pragma once
 
#include <vxd/simulation/SensitiveDetectorBase.h>
#include <framework/datastore/StoreArray.h>
#include <framework/datastore/RelationArray.h>
#include <array>
 
#ifdef VXD_SENSITIVEDETECTOR_DEBUG
#include <vxd/simulation/SensitiveDetectorDebugHelper.h>
#endif
 
namespace Belle2 {
  namespace VXD {

    template <class SimHitClass, class TrueHitClass> class SensitiveDetector: public SensitiveDetectorBase {
    public:
      explicit SensitiveDetector(VXD::SensorInfoBase* sensorInfo);
 
    private:
      int saveTrueHit(const SensorTraversal& traversal) override;

      int saveSimHit(const SensorTraversal& traversal,
                     const SensorTraversal::range& points) override;

      void saveRelations(const SensorTraversal& traversal, int trueHitIndex,
                         std::vector<std::pair<unsigned int, float>> simHitIndices) override;

      std::array<float, 3> vecToFloat(const G4ThreeVector& vec)
      {
        return std::array<float, 3> {{(float)vec.x(), (float)vec.y(), (float)vec.z()}};
      }

      StoreArray<SimHitClass> m_simhits;
      StoreArray<TrueHitClass> m_truehits;
      StoreArray<MCParticle> m_mcparticles;
      RelationArray m_relMCSimHits {m_mcparticles, m_simhits};
      RelationArray m_relMCTrueHits {m_mcparticles, m_truehits};
      RelationArray m_relTrueSimHits {m_truehits, m_simhits};
    };
 
    template <class SimHitClass, class TrueHitClass>
    SensitiveDetector<SimHitClass, TrueHitClass>::SensitiveDetector(VXD::SensorInfoBase* sensorInfo):
      VXD::SensitiveDetectorBase(sensorInfo)
    {
      //Register output collections.
      //m_mcparticles.isRequired();
      m_simhits.registerInDataStore();
      m_truehits.registerInDataStore();
      m_relMCSimHits.registerInDataStore();
      m_relMCTrueHits.registerInDataStore();
      m_relTrueSimHits.registerInDataStore();
      registerMCParticleRelation(m_relMCSimHits, RelationArray::c_negativeWeight);
      registerMCParticleRelation(m_relMCTrueHits, RelationArray::c_deleteElement);
    }
 
    template <class SimHitClass, class TrueHitClass>
    int SensitiveDetector<SimHitClass, TrueHitClass>::saveTrueHit(const SensorTraversal& traversal)
    {
      //We have the full sensor traversal information, we only lack the midpoint so lets get it
      StepInformation midPoint = findMidPoint(traversal);
 
      //By now everything is done: just convert the position and momenta in float[3] and create a truehit
      auto posEntry = vecToFloat(traversal.front().position);
      auto posMidPoint = vecToFloat(midPoint.position);
      auto posExit = vecToFloat(traversal.back().position);
      auto momEntry = vecToFloat(traversal.front().momentum);
      auto momMidPoint = vecToFloat(midPoint.momentum);
      auto momExit = vecToFloat(traversal.back().momentum);
      //And we should convert the electrons back in energy ...
      const double energyDep = traversal.getElectrons() * Const::ehEnergy;
 
      //create the truehit
      int trueHitIndex = m_truehits.getEntries();
      m_truehits.appendNew(getSensorID(), posEntry.data(), posMidPoint.data(), posExit.data(),
                           momEntry.data(), momMidPoint.data(), momExit.data(), energyDep, midPoint.time);
#ifdef VXD_SENSITIVEDETECTOR_DEBUG
      SensitiveDetectorDebugHelper::getInstance().addTrueHit(m_truehits[trueHitIndex]);
#endif
      return trueHitIndex;
    }
 
    template <class SimHitClass, class TrueHitClass>
    int SensitiveDetector<SimHitClass, TrueHitClass>::saveSimHit(const SensorTraversal& traversal,
        const SensorTraversal::range& points)
    {
      //Convert position to floats here
      auto posIn = vecToFloat(points.first->position);
      auto posOut = vecToFloat(points.second->position);
      auto electronProfile = simplifyEnergyDeposit(points);
 
      //Create the simhit
      int simHitIndex = m_simhits.getEntries();
      SimHitClass* simhit = m_simhits.appendNew(getSensorID(), traversal.getPDGCode(), points.first->time,
                                                posIn.data(), posOut.data());
      simhit->setEnergyDeposit(electronProfile);
#ifdef VXD_SENSITIVEDETECTOR_DEBUG
      int startPoint = std::distance(traversal.begin(), (SensorTraversal::const_iterator)points.first);
      int endPoint = std::distance(traversal.begin(), (SensorTraversal::const_iterator)points.second);
      SensitiveDetectorDebugHelper::getInstance().addSimHit(simhit, startPoint, endPoint);
#endif
      return simHitIndex;
    }
 
    template <class SimHitClass, class TrueHitClass>
    void SensitiveDetector<SimHitClass, TrueHitClass>::saveRelations(const SensorTraversal& traversal,
        int trueHitIndex, std::vector<std::pair<unsigned int, float>> simHitIndices)
    {
      m_relMCSimHits.add(traversal.getTrackID(), simHitIndices.begin(), simHitIndices.end());
      //If there is no truehit there are obviously no relations ;)
      if (trueHitIndex >= 0) {
        m_relMCTrueHits.add(traversal.getTrackID(), trueHitIndex, traversal.getElectrons());
        m_relTrueSimHits.add(trueHitIndex, simHitIndices.begin(), simHitIndices.end());
      }
    }
 
  } //VXD Namespace
} //Belle2 namespace