MdtSensitiveDetector.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
*/
 
#include "MuonSensitiveDetectorsR4/Utils.h"
#include "MdtSensitiveDetector.h"
 
 
#include <MCTruth/TrackHelper.h>
#include <G4Geantino.hh>
#include <G4ChargedGeantino.hh>
 
#include <limits>
#include <iostream>
#include <GeoPrimitives/CLHEPtoEigenConverter.h>
#include <GeoModelHelpers/throwExcept.h>
#include <xAODMuonSimHit/MuonSimHitAuxContainer.h>
#include <GaudiKernel/SystemOfUnits.h>
#include <StoreGate/ReadHandle.h>
 
using namespace MuonGMR4;
using namespace CxxUtils;
using namespace ActsTrk;
namespace MuonG4R4{
 
MdtSensitiveDetector::MdtSensitiveDetector(const std::string& name, 
                                           const std::string& output_key,
                                           const std::string& trfStore_key,
                                           const MuonGMR4::MuonDetectorManager* detMgr):
    G4VSensitiveDetector{name},
    AthMessaging{name},
    m_writeHandle{output_key},
    m_trfCacheKey{trfStore_key},
    m_detMgr{detMgr} {
    m_trfCacheKey.initialize().ignore();
}
// Implemenation of memebr functions
void MdtSensitiveDetector::Initialize(G4HCofThisEvent*) {
  if (m_writeHandle.isValid()) {
      ATH_MSG_VERBOSE("Simulation hit container "<<m_writeHandle.fullKey()<<" is already written");
      return;
  }
  if (!m_writeHandle.recordNonConst(std::make_unique<xAOD::MuonSimHitContainer>(),
                                    std::make_unique<xAOD::MuonSimHitAuxContainer>()).isSuccess()) {
      THROW_EXCEPTION(" Failed to record "<<m_writeHandle.fullKey());     
  }
  ATH_MSG_DEBUG("Output container "<<m_writeHandle.fullKey()<<" has been successfully created");
}
 
const MuonGMR4::MdtReadoutElement* MdtSensitiveDetector::getReadoutElement(const G4TouchableHistory* touchHist) const {
   const std::string stationVolume = touchHist->GetVolume(2)->GetName();
   using namespace MuonGMR4;
   const std::vector<std::string> volumeTokens = tokenize(stationVolume, "_");
   ATH_MSG_VERBOSE("Name of the station volume is "<<stationVolume);
   if (volumeTokens.size() != 5) {
      THROW_EXCEPTION(" Cannot deduce the station name from "<<stationVolume);
   }
   const std::string stName = volumeTokens[0].substr(0,3);
   const int stationEta = atoi(volumeTokens[2]);
   const int stationPhi = atoi(volumeTokens[3]);
   const int multiLayer = atoi(volumeTokens[4]);
   const MdtIdHelper& idHelper{m_detMgr->idHelperSvc()->mdtIdHelper()};
   const Identifier detElId = idHelper.multilayerID(idHelper.elementID(stName,stationEta, stationPhi), multiLayer);
   const MdtReadoutElement* readOutEle = m_detMgr->getMdtReadoutElement(detElId);
   if (!readOutEle) {
      THROW_EXCEPTION(" Failed to retrieve a valid detector element from "
                     <<m_detMgr->idHelperSvc()->toStringDetEl(detElId)<<" "<<stationVolume);
   }
   return readOutEle;
}
 
G4bool MdtSensitiveDetector::ProcessHits(G4Step* aStep,G4TouchableHistory* /*ROHist*/) {
    G4Track* currentTrack = aStep->GetTrack();
 
    // MDTs sensitive to charged particle only
    if (currentTrack->GetDefinition()->GetPDGCharge() == 0.0) {
      if (currentTrack->GetDefinition()!= G4Geantino::GeantinoDefinition()) return true;
      else if (currentTrack->GetDefinition() != G4ChargedGeantino::ChargedGeantinoDefinition()) return true;
    }
 
    
    constexpr double velCutOff = 10.*Gaudi::Units::micrometer / Gaudi::Units::second;
    if (currentTrack->GetVelocity() < velCutOff) return true;
 
    const G4TouchableHistory* touchHist = static_cast<const G4TouchableHistory*>(currentTrack->GetTouchable());
    const MdtReadoutElement* reEle{getReadoutElement(touchHist)};
 
    ActsGeometryContext gctx{};
 
    SG::ReadHandle<DetectorAlignStore> trfStoreHandle{m_trfCacheKey};
    if (!trfStoreHandle.isValid()) {
      ATH_MSG_FATAL("Failed to retrieve "<<m_trfCacheKey.fullKey()<<".");
      return false;
    }
    gctx.setStore(std::make_unique<DetectorAlignStore>(*trfStoreHandle));
    
    const Identifier HitID = getIdentifier(gctx, reEle, touchHist);
    if (!HitID.is_valid()) {
        ATH_MSG_VERBOSE("No valid hit found");
        return true;
    }
 
    const Amg::Transform3D globalToLocal{reEle->globalToLocalTrans(gctx, reEle->measurementHash(HitID))};
 
    // transform pre and post step positions to local positions
    const Amg::Vector3D trackPosition{Amg::Hep3VectorToEigen(currentTrack->GetPosition())};
    const Amg::Vector3D trackDirection{Amg::Hep3VectorToEigen(currentTrack->GetMomentumDirection())};
 
    const Amg::Vector3D trackLocPos{globalToLocal * trackPosition};  
    const Amg::Vector3D trackLocDir{globalToLocal.linear()* trackDirection};
  
    const double lambda = Amg::intersect<3>(Amg::Vector3D::Zero(), Amg::Vector3D::UnitZ(),
                                       trackLocPos, trackLocDir).value_or(0);
  
    const Amg::Vector3D driftHit{trackLocPos + lambda * trackLocDir};
 
    const double globalTime{currentTrack->GetGlobalTime() +  lambda / currentTrack->GetVelocity()};
  
    TrackHelper trHelp{currentTrack};
 
    ATH_MSG_VERBOSE(" Dumping of hit "<<m_detMgr->idHelperSvc()->toString(HitID)
                  <<", barcode: "<<trHelp.GenerateParticleLink().barcode()
                  <<", "<<(*currentTrack)
                  <<", driftCircle: "<<Amg::toString(driftHit, 2)
                  <<", direction "<<Amg::toString(trackLocDir, 2)
                  <<" to SimHit container ahead. ");
 
    xAOD::MuonSimHit* hit = new xAOD::MuonSimHit();
    m_writeHandle->push_back(hit);  
    hit->setIdentifier(HitID); 
    hit->setLocalPosition(xAOD::toStorage(driftHit));  
    hit->setLocalDirection(xAOD::toStorage(trackLocDir));
    hit->setMass(currentTrack->GetDefinition()->GetPDGMass());
    hit->setGlobalTime(globalTime);
    hit->setPdgId(currentTrack->GetDefinition()->GetPDGEncoding());
    hit->setEnergyDeposit(aStep->GetTotalEnergyDeposit());
    hit->setKineticEnergy(currentTrack->GetKineticEnergy());
    hit->setGenParticleLink(trHelp.GenerateParticleLink());
 
  return true;
}
Identifier MdtSensitiveDetector::getIdentifier(const ActsGeometryContext& gctx,
                                               const MuonGMR4::MdtReadoutElement* readOutEle,
                                               const G4TouchableHistory* touchHist) const {
   const Amg::Transform3D localToGlobal{getTransform(touchHist, 0)};
   Amg::Vector3D refTubePos = (readOutEle->globalToLocalTrans(gctx, readOutEle->measurementHash(1,1)) * localToGlobal).translation();
   ATH_MSG_VERBOSE("Position of the tube wire w.r.t. the first tube in the multi layer "<<Amg::toString(refTubePos, 2));
   static const double layerPitch = 1./ std::sin(60*Gaudi::Units::deg);
   const int layer = std::round(refTubePos.x() * layerPitch / readOutEle->tubePitch()) +1;
   if (layer <= 0) {
      THROW_EXCEPTION("Tube hit in nirvana -- It seems that the tube position "
                      <<Amg::toString(refTubePos, 2)<<", perp: "<<refTubePos.perp()
                      <<" is outside of the volume envelope "
                      <<m_detMgr->idHelperSvc()->toStringDetEl(readOutEle->identify())<<". ");      
   }
   refTubePos  = (readOutEle->globalToLocalTrans(gctx, readOutEle->measurementHash(layer,1)) * localToGlobal).translation();
   const double tubePitches = refTubePos.y() / readOutEle->tubePitch();
   unsigned int tube = std::round(tubePitches) + 1;
   tube = std::max(1u, std::min(readOutEle->numTubesInLay(), tube));
  const Amg::Transform3D closureCheck{readOutEle->globalToLocalTrans(gctx, 
                                      readOutEle->measurementHash(layer, tube))*localToGlobal};
    if (!Amg::isIdentity(closureCheck)) {
        ATH_MSG_WARNING("Correction needed "<<layer<<","<<tube<<" "<<Amg::toString(closureCheck));
        if (closureCheck.translation().y() > 0) ++tube;
        else --tube;
    }
   
   const IdentifierHash tubeHash = readOutEle->measurementHash(layer, tube);
   const Identifier tubeId = readOutEle->measurementId(tubeHash);
   {
        const Amg::Transform3D closureCheck{readOutEle->globalToLocalTrans(gctx, tubeHash)*localToGlobal};
        if (!Amg::isIdentity(closureCheck)) {
            THROW_EXCEPTION("Tube hit in Nirvana --  It seems that the tube position "
                             <<Amg::toString(refTubePos, 2)<<", perp: "<<refTubePos.perp()
                             <<" is outside of the volume envelope "
                             <<m_detMgr->idHelperSvc()->toStringDetEl(readOutEle->identify())<<". "
                             <<"Layer: "<<layer<<", tube: "<<tube<<" "
                             <<Amg::toString(closureCheck)
                             <<"tube volume : "<<touchHist->GetVolume(0)->GetName()
                             <<" mdt chamber: "<<touchHist->GetVolume(2)->GetName());    
        }
   }
   ATH_MSG_VERBOSE("Tube & layer number candidate "<<tube<<", "<<layer<<" back and forth transformation "
                     <<Amg::toString(closureCheck));
   return tubeId;  
}
 
}