MCFTDigitMonitor.cpp

// Include files 

// local
#include "MCFTDigitMonitor.h"
#include "SiPMResponse.h"

#include <boost/range/irange.hpp>

//-----------------------------------------------------------------------------
// Implementation file for class : MCFTDigitMonitor
//
// 2012-07-05 : Eric Cogneras
//-----------------------------------------------------------------------------

// Declaration of the Algorithm Factory
DECLARE_ALGORITHM_FACTORY( MCFTDigitMonitor )

MCFTDigitMonitor::MCFTDigitMonitor(const std::string& name,ISvcLocator* pSvcLocator) :
       Consumer(name, pSvcLocator, 
       KeyValue{"DigitLocation", LHCb::MCFTDigitLocation::Default}) {}

//=============================================================================
// Initialization
//=============================================================================
StatusCode MCFTDigitMonitor::initialize() {

    //StatusCode sc = GaudiHistoAlg::initialize(); // must be executed first
    //if ( sc.isFailure() ) return sc;  // error printed already by GaudiAlgorithm

    /// Retrieve and initialize DeFT
    m_deFT = getDet<DeFTDetector>( DeFTDetectorLocation::Default );
    if ( m_deFT == nullptr )
        return Error("Could not initialize DeFTDetector.", StatusCode::FAILURE);
    if( m_deFT->version() < 61 )
        return Error("This version requires FTDet v6.1 or higher", StatusCode::FAILURE);

    // Plot the electronics response function
    SiPMResponse* sipmResponse = tool<SiPMResponse>("SiPMResponse","SiPMResponse");
    for( auto time : boost::irange(-60, 120) ) {
        plot(time, "ElectronicsResponse",
                "Electronics response function; #it{t} [ns]; Relative gain",
                -60., 120., 180, sipmResponse->response(time) );
    }

    return StatusCode::SUCCESS;
}


//=============================================================================
// Main execution
//=============================================================================
void MCFTDigitMonitor::operator()(const LHCb::MCFTDigits& mcDigitsCont) const {

    // Count signal, noise and spillover digits
    int nSignal = 0, nNoise = 0, nSpillover = 0;

    // Loop over MCFTDigits
    for( auto mcDigit : mcDigitsCont ) {

        // Check if digit is from pure spillover or from pure noise
        bool isSpillover = true;
        std::set<const LHCb::MCHit*> mcHits;
        for (const auto& deposit : mcDigit->deposits() ) {
            const auto& mcHit = deposit->mcHit();
            if( mcHit != nullptr ) {
                mcHits.insert(mcHit);
                if( mcHit->parent()->registry()->identifier() == "/Event/MC/FT/Hits" )
                    isSpillover = false;
            }
        }
        if ( mcHits.empty() ) isSpillover = false;

        std::string hitType = "Signal/";
        if( mcHits.empty() ) {
            nNoise++;
            hitType = "Noise/";
        } else if( isSpillover ) {
            nSpillover++;
            hitType = "Spillover/";
        } else {
            nSignal++;
        }
        // plot number of hits
        plot( mcHits.size(), hitType+"nMCHitsPerDigit",
                "Number of MCHits per digit; MCHits/digit; Number of digits" ,
                -0.5 , 20.5, 21);

        // plot digit PEs
        plot(mcDigit->photoElectrons(), hitType+"nPhotoElectrons",
                "PhotoElectrons per digit; PEs; Number of digits" ,
                0., 30., 100);

        // plot digit adc counts
        plot(mcDigit->adcCount(), hitType+"ADCCounts",
                "Charge in 2bit ADC; Charge [ADC]; Number of digits" ,
                -0.5, 3.5, 4);

        // deposits distribution in (x,y)
        if( !mcHits.empty() ) {
            plot2D((*mcHits.begin())->midPoint().x(), (*mcHits.begin())->midPoint().y(),
                    hitType+"XYDistribution", "Digits XY Distribution; #it{x} [mm]; #it{y} [mm]" ,
                    -3600. , 3600. , -2700. , 2700. , 100, 100);
            // average photoElectron counts per hits
            plot(mcDigit->photoElectrons()/float(mcHits.size()), hitType+"PEsPerMCHit",
                    "PE counts / MCHit; PE count / nHits;Number of digits" ,
                    0. , 30., 100);
        }

        // plot occupancy
        const LHCb::FTChannelID chanID = mcDigit -> channelID();
        plot((float)chanID.station(), hitType+"DigitsPerStation",
                "Digits per station; Station; Digits" , 0.5, 3.5, 3);
        plot((float)chanID.module(), hitType+"DigitsPerModule",
                "Digits per module; Module; Digits" , -0.5, 5.5, 6);
        plot((float)chanID.sipmInModule(), hitType+"DigitsPerSiPM",
                "Digits per SiPM; SiPMID; Digits", -0.5, 15.5, 16);
        plot((float)chanID.channel(), hitType+"DigitsPerChannel",
                "Digits per channel; Channel; Digits", -0.5, 127.5, 128);

        const DeFTModule* module = m_deFT->findModule( chanID );
        if ( module != nullptr ) {
            int pseudoChannel = module->pseudoChannel( chanID );
            plot(pseudoChannel, hitType+"DigitsPerPseudoChannel",
                    "Digits per pseudo channel;Pseudo channel;Digits/(64 channels)",
                    0., 12288., 192);
        }

        // Count the number of photons for this digit
        int nPhotonsPerDigit = 0;
        for( const auto& deposit : mcDigit -> deposits() )
            nPhotonsPerDigit += deposit->nPhotons();
        plot((float)nPhotonsPerDigit, hitType+"nPhotonsPerDigit",
                "Photons per digit;Photons;Number of digits", 0.5, 60.5, 60);

        // photons/digit vs the final photoElectron charge
        plot2D((float)nPhotonsPerDigit, mcDigit->photoElectrons(), hitType+"PEvsPhotons",
                "PEvsPhotons; Photons; PE" , 0.5, 40.5 ,0.5, 40.5, 40, 40);
    }

    // number of clusters
    plot(mcDigitsCont.size(), "nDigits",
            "Number of digits; Digits/event; Events", 0.,75e3, 100);
    plot(nSignal, "nSignalDigits",
            "Number of signal digits; Digits/event; Events", 0.,75e3, 100);
    plot(nSpillover, "nSpilloverDigits",
            "Number of spillover digits; Digits/event; Events", 0.,75e3, 100);
    plot(nNoise, "nNoiseDigits",
            "Number of noise digits; Digits/event; Events", 0.,75e3, 100);

    return;
}

//=============================================================================