MmFastDigiTool.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
*/
#include "MmFastDigiTool.h"
#include "xAODMuonViews/ChamberViewer.h"
#include "CLHEP/Random/RandGaussZiggurat.h"
#include "CLHEP/Random/RandFlat.h"
 
namespace {
    constexpr double percentage(unsigned int numerator, unsigned int denom) {
        return 100. * numerator / std::max(denom, 1u);
    }
}
namespace MuonR4 {
    
    MmFastDigiTool::MmFastDigiTool(const std::string& type, const std::string& name, const IInterface* pIID):
        MuonDigitizationTool{type,name, pIID} {}
 
    StatusCode MmFastDigiTool::initialize() {
        ATH_CHECK(MuonDigitizationTool::initialize());
        ATH_CHECK(m_writeKey.initialize());
        ATH_CHECK(m_uncertCalibKey.initialize());
        ATH_CHECK(m_effiDataKey.initialize(!m_effiDataKey.empty()));
        return StatusCode::SUCCESS;
    }
    StatusCode MmFastDigiTool::finalize() {
        std::stringstream statstr{};
        unsigned allHits{0};
        for (unsigned int g = 0; g < m_allHits.size(); ++g) {
            allHits += m_allHits[g];
            statstr<<" *** Layer "<<(g+1)<<" "<<percentage(m_acceptedHits[g],m_allHits[g])
                   <<"% of "<<m_allHits[g]<<std::endl; 
        }
        if(!allHits) return StatusCode::SUCCESS;
        ATH_MSG_INFO("Tried to convert "<<allHits<<" hits. Successes rate per layer  "<<std::endl<<statstr.str());
        return StatusCode::SUCCESS;
    }
    StatusCode MmFastDigiTool::digitize(const EventContext& ctx,
                                        const TimedHits& hitsToDigit,
                                        xAOD::MuonSimHitContainer* sdoContainer) const {
        const MmIdHelper& idHelper{m_idHelperSvc->mmIdHelper()};
        // Prepare the temporary cache
        DigiCache digitCache{};
        const Muon::DigitEffiData* efficiencyMap{nullptr};
        ATH_CHECK(retrieveConditions(ctx, m_effiDataKey, efficiencyMap));
        const NswErrorCalibData* errorCalibDB{nullptr};
        ATH_CHECK(retrieveConditions(ctx, m_uncertCalibKey, errorCalibDB));
 
        CLHEP::HepRandomEngine* rndEngine = getRandomEngine(ctx);
        xAOD::ChamberViewer viewer{hitsToDigit, m_idHelperSvc.get()};
        do {
            DeadTimeMap deadTimes{};
            for (const TimedHit& simHit : viewer) {
                if (m_digitizeMuonOnly && std::abs(simHit->pdgId()) != 13){
                    continue;
                }
                const Identifier hitId{simHit->identify()};
 
                const MuonGMR4::MmReadoutElement* readOutEle = m_detMgr->getMmReadoutElement(hitId);
                const Amg::Vector2D locPos{xAOD::toEigen(simHit->localPosition()).block<2,1>(0,0)};
 
                const unsigned int hitGapInNsw = (idHelper.multilayer(hitId) -1) * 4 + idHelper.gasGap(hitId) -1;
                ++m_allHits[hitGapInNsw];
 
                const MuonGMR4::StripDesign& design{readOutEle->stripLayer(hitId).design()};
                bool isValid{false};
 
                int channelNumber = design.stripNumber(locPos);
                if(channelNumber<0){
                    if (!design.insideTrapezoid(locPos)) {
                        ATH_MSG_WARNING("Hit "<<m_idHelperSvc->toString(hitId)<<" "<<Amg::toString(locPos)
                                      <<" is outside bounds "<<std::endl<<design<<" rejecting it");
                    }
                    continue;
                }
                const Identifier clusId = idHelper.channelID(hitId, 
                                                             idHelper.multilayer(hitId), 
                                                             idHelper.gasGap(hitId), 
                                                             channelNumber, isValid);
 
                if(!isValid) {
                    ATH_MSG_WARNING("Invalid strip identifier for layer " << m_idHelperSvc->toStringGasGap(hitId) 
                                 << " channel " << channelNumber << " locPos " << Amg::toString(locPos));
                    continue;
                }
 
                if(efficiencyMap && efficiencyMap->getEfficiency(clusId) < CLHEP::RandFlat::shoot(rndEngine, 0., 1.)){
                    continue;
                }
                if (!passDeadTime(clusId, hitTime(simHit), m_deadTime, deadTimes)) {
                    continue;
                }
                NswErrorCalibData::Input errorCalibInput{};
                errorCalibInput.stripId = clusId;
                errorCalibInput.locTheta = M_PI- simHit->localDirection().theta();
                errorCalibInput.clusterAuthor=66; // cluster time projection method
                const double uncert = errorCalibDB->clusterUncertainty(errorCalibInput);
                ATH_MSG_VERBOSE("mm hit has theta " << errorCalibInput.locTheta / Gaudi::Units::deg << " and uncertainty " << uncert);
 
 
                constexpr int dummyResponseTime = 100;
                constexpr float dummyDepositedCharge = 66666;
 
                const double newLocalX = CLHEP::RandGaussZiggurat::shoot(rndEngine, locPos.x(), uncert);
                const int newChannel = design.stripNumber(newLocalX * Amg::Vector2D::UnitX());
                if (newChannel < 0) {                
                    continue;
                }
 
                const Identifier digitId = idHelper.channelID(hitId, 
                                                              idHelper.multilayer(hitId), 
                                                              idHelper.gasGap(hitId), 
                                                              newChannel, isValid);
 
                if(!isValid) {
                    ATH_MSG_WARNING("Invalid strip identifier for layer " << m_idHelperSvc->toStringGasGap(hitId) 
                                    << " channel " << newChannel << " locPos " << Amg::toString(locPos));
                    continue;
                }
 
                const double pull = (newLocalX - (*design.center(newChannel)).x()) / (design.stripPitch() * std::cos(design.stereoAngle()));
                const double w1 = CLHEP::RandFlat::shoot(rndEngine, 0., 0.5 *(1. - pull)); 
                const double w2 = 1 - pull -2*w1;
                const double w3 = pull + w1;
                MmDigitCollection* outColl = fetchCollection(hitId, digitCache);
 
                const Identifier digitIdB = idHelper.channelID(hitId, 
                                                               idHelper.multilayer(hitId), 
                                                               idHelper.gasGap(hitId), 
                                                               newChannel - 1, isValid);
 
                if (isValid) {
                    outColl->push_back(std::make_unique<MmDigit>(digitIdB, dummyResponseTime, w1 * dummyDepositedCharge));
                }
                outColl->push_back(std::make_unique<MmDigit>(digitId,dummyResponseTime, w2 * dummyDepositedCharge));
 
                const Identifier digitIdA = idHelper.channelID(hitId, 
                                                               idHelper.multilayer(hitId), 
                                                               idHelper.gasGap(hitId), 
                                                               newChannel + 1, isValid);
                if (isValid) {
                    outColl->push_back(std::make_unique<MmDigit>(digitIdA, dummyResponseTime, w3 * dummyDepositedCharge));
                }
 
                xAOD::MuonSimHit* sdoHit = addSDO(simHit, sdoContainer);
                sdoHit->setIdentifier(clusId);
                ++m_acceptedHits[hitGapInNsw];
            }
        } while(viewer.next());
        ATH_CHECK(writeDigitContainer(ctx, m_writeKey, std::move(digitCache), idHelper.module_hash_max()));
        return StatusCode::SUCCESS;
    }   
}