diff --git a/TileCalorimeter/TileMonitoring/CMakeLists.txt b/TileCalorimeter/TileMonitoring/CMakeLists.txt index 517fe6d3d5d498e69d1760adbc76e93337b543fb..eb7cffbbde6a5e4c6319a44b961faa60b87334e1 100644 --- a/TileCalorimeter/TileMonitoring/CMakeLists.txt +++ b/TileCalorimeter/TileMonitoring/CMakeLists.txt @@ -56,4 +56,5 @@ atlas_add_component( TileMonitoring # Install files from the package: atlas_install_headers( TileMonitoring ) +atlas_install_python_modules( python/*.py ) atlas_install_joboptions( share/*.py ) diff --git a/TileCalorimeter/TileMonitoring/python/TileJetMonitorAlgorithm.py b/TileCalorimeter/TileMonitoring/python/TileJetMonitorAlgorithm.py new file mode 100644 index 0000000000000000000000000000000000000000..267556693d48bf909b997bb662b81e6de972e710 --- /dev/null +++ b/TileCalorimeter/TileMonitoring/python/TileJetMonitorAlgorithm.py @@ -0,0 +1,268 @@ +# +# Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration +# + +''' +@file TileJetMonitorAlgorithm.py +@author +@date +@brief Python configuration of TileJetMonitorAlgorithm algorithm for the Run III +''' + + +def TileJetMonitoringConfig(inputFlags, **kwargs): + + ''' Function to configure TileJetMonitorAlgorithm algorithm in the monitoring system.''' + + + from AthenaCommon import CfgMgr + + # Define one top-level monitoring algorithm. The new configuration + # framework uses a component accumulator. + from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator + result = ComponentAccumulator() + + from TileGeoModel.TileGMConfig import TileGMCfg + result.merge(TileGMCfg(inputFlags)) + + from LArGeoAlgsNV.LArGMConfig import LArGMCfg + result.merge(LArGMCfg(inputFlags)) + + from TileConditions.TileConditionsConfig import tileCondCfg + result.merge(tileCondCfg(inputFlags)) + + # The following class will make a sequence, configure algorithms, and link + # them to GenericMonitoringTools + from AthenaMonitoring import AthMonitorCfgHelper + helper = AthMonitorCfgHelper(inputFlags,'TileMonitoring') + + # Adding an TileJetMonitorAlgorithm algorithm to the helper + from TileMonitoring.TileMonitoringConf import TileJetMonitorAlgorithm + tileJetMonAlg = helper.addAlgorithm(TileJetMonitorAlgorithm, 'TileJetMonAlg') + + from AthenaCommon.Constants import DEBUG, INFO, VERBOSE + # tileJetMonAlg.OutputLevel = VERBOSE + tileJetMonAlg.TriggerChain = '' + + for k, v in kwargs.items(): + setattr(tileJetMonAlg, k, v) + + DoEnergyProfiles = kwargs.get('DoEnergyProfiles', tileJetMonAlg.getDefaultProperty('DoEnergyProfiles')) + Do1DHistograms = kwargs.get('Do1DHistograms', tileJetMonAlg.getDefaultProperty('Do1DHistograms')) + DoEnergyDiffHistograms = kwargs.get('DoEnergyDiffHistograms', tileJetMonAlg.getDefaultProperty('DoEnergyDiffHistograms')) + + + if not inputFlags.DQ.DataType == 'heavyioncollision': + + tileJetMonAlg.JVT = CfgMgr.JetVertexTaggerTool() + + jetCleaningTool = CfgMgr.JetCleaningTool() + jetCleaningTool.CutLevel = "LooseBad" + jetCleaningTool.DoUgly = False + + tileJetMonAlg.JetCleaningTool = jetCleaningTool + result.addPublicTool(jetCleaningTool) + + jetPtMin = 20000 + jetTrackingEtaLimit = 2.4 + eventCleaningTool = CfgMgr.ECUtils__EventCleaningTool() + eventCleaningTool.JetCleaningTool = jetCleaningTool + eventCleaningTool.PtCut = jetPtMin + eventCleaningTool.EtaCut = jetTrackingEtaLimit + eventCleaningTool.JvtDecorator = "passJvt" + eventCleaningTool.OrDecorator = "passOR" + eventCleaningTool.CleaningLevel = jetCleaningTool.CutLevel + + tileJetMonAlg.EventCleaningTool = eventCleaningTool + tileJetMonAlg.JetTrackingEtaLimit = jetTrackingEtaLimit + tileJetMonAlg.JetPtMin = jetPtMin + + tileJetMonAlg.DoEventCleaning = True + tileJetMonAlg.DoJetCleaning = True + + else: + + tileJetMonAlg.DoEventCleaning = False + tileJetMonAlg.DoJetCleaning = False + + + # 1) Configure histogram with TileJetMonAlg algorithm execution time + executeTimeGroup = helper.addGroup(tileJetMonAlg, 'TileJetMonExecuteTime', 'Tile/') + executeTimeGroup.defineHistogram('TIME_execute', path = 'Jet', type='TH1F', + title = 'Time for execute TileJetMonAlg algorithm;time [#mus]', + xbins = 300, xmin = 0, xmax = 300000) + + + + from TileMonitoring.TileMonitoringCfgHelper import addValueVsModuleAndChannelMaps, getPartitionName + runNumber = inputFlags.Input.RunNumber[0] + + + # 2) Configure 2D histograms (profiles/maps) with Tile channel time vs module and channel per partion (DQ summary) + channelTimeDQGroup = helper.addGroup(tileJetMonAlg, 'TileJetChanTimeDQ', 'Tile/Jet/') + addValueVsModuleAndChannelMaps(channelTimeDQGroup, name = 'tileJetChanTime', title = 'Average time with jets', + path = 'DQ', type = 'TProfile2D', value='time', run = str(runNumber)) + + + gains = ['LG', 'HG'] + partitions = ['LBA', 'LBC', 'EBA', 'EBC'] + + # 3a) Configure 1D histograms with Tile channel time per partition + channelTimeGroup = helper.addGroup(tileJetMonAlg, 'TileJetChanTime', 'Tile/Jet/ChanTime/') + for partition in partitions: + title = 'Partition ' + partition + ': Tile Channel Time;time [ns];N' + name = 'channelTime' + partition + path = partition + channelTimeGroup.defineHistogram(name, title = title, path = path, type = 'TH1F', + xbins = 600, xmin = -30.0, xmax = 30.0) + + # 3b) Configure 1D histograms with Tile channel time per partition for extended barrels without scintillators + for partition in ['EBA', 'EBC']: + title = 'Partition ' + partition + ': Tile Channel Time (without scintillators);time [ns];N' + name = 'channelTime' + partition + '_NoScint' + path = partition + channelTimeGroup.defineHistogram(name, title = title, path = path, type = 'TH1F', + xbins = 600, xmin = -30.0, xmax = 30.0) + + + # Energy upper limits of the cell-time histograms + energiesHG = [500, 1000, 2000, 4000, 6000, 8000, 10000, 13000, 16000, 20000] + energiesLG = [25000, 30000, 40000, 50000, 65000, 80000] + energiesALL = {'LG' : energiesLG, 'HG' : energiesHG} + tileJetMonAlg.CellEnergyUpperLimitsHG = energiesHG + tileJetMonAlg.CellEnergyUpperLimitsLG = energiesLG + + # 4) Configure histograms with Tile cell time in energy slices per partition and gain + cellTimeGroup = helper.addGroup(tileJetMonAlg, 'TileJetCellTime', 'Tile/Jet/CellTime/') + for partition in partitions: + for gain in gains: + index = 0 + energies = energiesALL[gain] + for index in xrange(0, len(energies) + 1): + toEnergy = energies[index] if index < len(energies) else None + fromEnergy = energies[index - 1] if index > 0 else None + name = 'Cell_time_' + partition + '_' + gain + '_slice_' + str(index) + title = 'Partition ' + partition + ': ' + gain + ' Tile Cell time in energy range' + if not toEnergy: + title += ' > ' + str(fromEnergy) + ' MeV; time [ns]' + elif not fromEnergy: + title += ' < ' + str(toEnergy) + ' MeV; time [ns]' + else: + title += ' [' + str(fromEnergy) + ' .. ' + str(toEnergy) + ') MeV; time [ns]' + cellTimeGroup.defineHistogram(name, title = title, path = partition, type = 'TH1F', + xbins = 600, xmin = -30.0, xmax = 30.0) + + + + if DoEnergyProfiles: + + # 5) Configure 1D histograms (profiles) with Tile cell energy profile in energy slices per partition and gain + cellEnergyProfileGroup = helper.addGroup(tileJetMonAlg, 'TileJetCellEnergyProfile', 'Tile/Jet/CellTime/') + for partition in partitions: + for gain in gains: + name = 'index_' + partition + '_' + gain + name += ',energy_' + partition + '_' + gain + name += ';Cell_ene_' + partition + '_' + gain + '_prof' + title = 'Partition ' + partition + ': ' + gain + ' Tile Cell energy profile;Slice;Energy [MeV]' + xmax = len(energiesALL[gain]) + 0.5 + nbins = len(energiesALL[gain]) + 1 + cellEnergyProfileGroup.defineHistogram(name, title = title, path = partition, type = 'TProfile', + xbins = nbins, xmin = -0.5, xmax = xmax) + + + else: + + # 6) Configure 1D histograms with Tile cell energy in energy slices per partition, gain and slice + cellEnergyGroup = helper.addGroup(tileJetMonAlg, 'TileJetCellEnergy', 'Tile/Jet/CellTime/') + for partition in partitions: + for gain in gains: + energies = energiesALL[gain] + for index in xrange(0, len(energies) + 1): + toEnergy = energies[index] if index < len(energies) else 2 * energies[index - 1] + fromEnergy = energies[index - 1] if index > 0 else -1000 + name = 'Cell_ene_' + partition + '_' + gain + '_slice_' + str(index) + title = 'Partition ' + partition + ': ' + gain + ' Tile Cell Energy' + title += ' in energy range [' + str(fromEnergy) + ' .. ' + str(toEnergy) + ') MeV;Energy [MeV]' + cellEnergyGroup.defineHistogram(name, title = title, path = partition, type = 'TH1F', + xbins = 100, xmin = fromEnergy, xmax = toEnergy) + + + + from TileCalibBlobObjs.Classes import TileCalibUtils as Tile + + if Do1DHistograms: + + # 7) Configure 1D histograms with Tile channel time per channel + channelTime1DGroup = helper.addGroup(tileJetMonAlg, 'TileJetChanTime1D', 'Tile/Jet/ChanTime/') + + for ros in xrange(1, Tile.MAX_ROS): + for module in xrange(0, Tile.MAX_DRAWER): + for channel in xrange(0, Tile.MAX_CHAN): + moduleName = Tile.getDrawerString(ros, module) + title = 'Time in ' + moduleName + ' channel ' + str(channel) + ';time [ns];N' + name = moduleName + '_ch_' + str(channel) + '_1d' + path = getPartitionName(ros) + '/' + moduleName + channelTime1DGroup.defineHistogram(name, title = title, path = path, type = 'TH1F', + xbins = 600, xmin = -30.0, xmax = 30.0) + + + + if DoEnergyDiffHistograms: + + # 7) Configure 1D histograms with Tile cell relative energy difference between two channels per even channel + energyDiffGroup = helper.addGroup(tileJetMonAlg, 'TileJetEnergyDiff', 'Tile/Jet/EnergyDiff/') + + for ros in xrange(1, Tile.MAX_ROS): + for module in xrange(0, Tile.MAX_DRAWER): + for channel in xrange(0, Tile.MAX_CHAN): + if not channel % 2: + for gain in gains: + moduleName = Tile.getDrawerString(ros, module) + title = 'Tile Cell Energy difference in ' + moduleName + ' channel ' + str(channel) + ' ' + gain + title += ';#frac{ene1 - ene2}{ene1 + ene2}' + name = moduleName + '_enediff_' + gain + '_ch1_' + str(channel) + path = getPartitionName(ros) + '/' + moduleName + energyDiffGroup.defineHistogram(name, title = title, path = path, type = 'TH1F', + xbins = 100, xmin = -1.0, xmax = 1.0) + + + + accumalator, sequence = helper.result() + result.merge(accumalator) + return result, sequence + +if __name__=='__main__': + + # Setup the Run III behavior + from AthenaCommon.Configurable import Configurable + Configurable.configurableRun3Behavior = 1 + + # Setup logs + from AthenaCommon.Logging import log + from AthenaCommon.Constants import DEBUG,INFO + log.setLevel(INFO) + + # Set the Athena configuration flags + from AthenaConfiguration.AllConfigFlags import ConfigFlags + + # from AthenaConfiguration.TestDefaults import defaultTestFiles + # ConfigFlags.Input.Files = defaultTestFiles.ESD + + ConfigFlags.Input.Files = ['/cvmfs/atlas-nightlies.cern.ch/repo/data/data-art/Tier0ChainTests/q431/21.0/myESD.pool.root'] + ConfigFlags.Input.isMC = False + + ConfigFlags.Output.HISTFileName = 'TileJetMonitorOutput.root' + ConfigFlags.lock() + + # Initialize configuration object, add accumulator, merge, and run. + from AthenaConfiguration.MainServicesConfig import MainServicesSerialCfg + from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg + cfg = MainServicesSerialCfg() + cfg.merge(PoolReadCfg(ConfigFlags)) + + tileJetMonitorAccumulator, tileMonitoringSequence = TileJetMonitoringConfig(ConfigFlags, + Do1DHistograms = True, + DoEnergyDiffHistograms = True) + cfg.merge(tileJetMonitorAccumulator) + + cfg.run() diff --git a/TileCalorimeter/TileMonitoring/python/TileMonitoringCfgHelper.py b/TileCalorimeter/TileMonitoring/python/TileMonitoringCfgHelper.py new file mode 100644 index 0000000000000000000000000000000000000000..b2fb690ed036261d0f3ab52c18842921fbb4a84c --- /dev/null +++ b/TileCalorimeter/TileMonitoring/python/TileMonitoringCfgHelper.py @@ -0,0 +1,91 @@ +# +# Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration +# + +''' +@file TileMonitoringHelper.py +@author +@date +@brief Helper functions for Run 3 Tile monitoring algorithm configuration +''' + + +_cellNameEB = ['E3', 'E4', 'D4', 'D4', 'C10', 'C10', 'A12', 'A12', 'B11', 'B11', 'A13', 'A13', + 'E1', 'E2', 'B12', 'B12', 'D5', 'D5', 'E3*', 'E4*', 'A14', 'A14', 'B13', 'B13', + '', '', '', '', '', '', 'B14', 'A15', 'A15', '', '', 'B14', + 'B15', 'D6', 'D6', 'B15', 'A16', 'A16', '', '', '', '', '', ''] + + +_cellNameLB = ['D0', 'A1', 'B1', 'B1', 'A1', 'A2', 'B2', 'B2', 'A2', 'A3', 'A3', 'B3', + 'B3', 'D1', 'D1', 'A4', 'B4', 'B4', 'A4', 'A5', 'A5', 'B5', 'B5', 'A6', + 'D2', 'D2', 'A6', 'B6', 'B6', 'A7', '', '', 'A7', 'B7', 'B7', 'A8', + 'A9', 'A9', 'A8', 'B8', 'B8', 'D3', 'B9', '', 'D3', 'A10', 'A10', 'B9'] + + +_partitionName = {0: 'AUX', 1 : 'LBA', 2 : 'LBC', 3 : 'EBA', 4 : 'EBC'} + + +def getCellName(partition, channel): + ''' + This function returns name of Tile cell for given partition and channel. + + Arguments: + partition -- Tile partition name (LBA, LBC, EBA, EBC) + channel -- Tile channel number ([0..47]) + ''' + return _cellNameLB[channel] if partition.startswith('L') else _cellNameEB[channel] + + +def getPartitionName(ros): + ''' + This function returns name of Tile partition for given ROS. + + Arguments: + ros -- Tile ROS ([0..5]) + ''' + return _partitionName[ros] + + +def addValueVsModuleAndChannelMaps(group, name, title, path, subDirectory = False, type = 'TH2D', value = '', trigger = '', run = ''): + ''' + This function configures 2D histograms (maps) with Tile monitored value vs module and channel per partion. + + Arguments: + group -- Group (technically, a GenericMonitoringTool instance) + name -- Name of histogram (actual name is constructed dynamicaly like: name + partition + trigger) + title -- Title of histogram (actual title is constructed dynamicaly like: run + trigger + partion + title) + path -- Path in file for histogram (relative to the path of given group) + subDirectory -- Put the configured histograms into sub directory named like partion (True, False) + type -- Type of histogram (TH2D, TProfile2D) + value -- Name of monitored value (needed for TProfile2D) + trigger -- Name of trigger (given it will be put into title and name of histogram) + run -- Run number (given it will be put into the title) + ''' + + from TileCalibBlobObjs.Classes import TileCalibUtils as Tile + + for ros in range(1, Tile.MAX_ROS): + partition = getPartitionName(ros) + labels = [] + for module in range(1, Tile.MAX_DRAWER + 1): # modules start from 1 + label = partition + '0' + str(module) if module < 10 else partition + str(module) + labels.append(label) + + for channel in range(0, Tile.MAX_CHAN): + cellName = getCellName(partition, channel) + label = cellName + '_' + 'ch' + str(channel) if cellName else 'ch' + str(channel) + labels.append(label) + + fullName = 'module' + partition + ',channel' + partition + if 'Profile' in type: fullName += (',' + value + partition) + fullName += ';' + name + partition + trigger + + fullPath = path + '/' + partition if subDirectory else path + + fullTitle = 'Partition ' + partition + ': ' + title + if trigger: fullTitle = 'Trigger ' + trigger + ' ' + fullTitle + if run: fullTitle = 'Run ' + run + ' ' + fullTitle + + group.defineHistogram( fullName, path = fullPath, type = type, title = fullTitle, + xbins = 64, xmin = 0.5, xmax = 64.5, ybins = 48, ymin = -0.5, ymax = 47.5, + labels = labels ) diff --git a/TileCalorimeter/TileMonitoring/src/TileJetMonitorAlgorithm.cxx b/TileCalorimeter/TileMonitoring/src/TileJetMonitorAlgorithm.cxx new file mode 100644 index 0000000000000000000000000000000000000000..cacfea01c2b879421e4d258e2ca4c4b1efbd7a65 --- /dev/null +++ b/TileCalorimeter/TileMonitoring/src/TileJetMonitorAlgorithm.cxx @@ -0,0 +1,460 @@ +/* + Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration +*/ + +#include "TileJetMonitorAlgorithm.h" + +#include "TileEvent/TileCell.h" +#include "TileIdentifier/TileHWID.h" +#include "TileIdentifier/TileFragHash.h" +#include "TileCalibBlobObjs/TileCalibUtils.h" +#include "TileConditions/TileCablingSvc.h" +#include "TileConditions/TileCablingService.h" + +#include "CaloIdentifier/TileID.h" +#include "CaloEvent/CaloCell.h" +//#include "CaloEvent/CaloCluster.h" + + +#include "StoreGate/ReadHandle.h" +#include <xAODCore/ShallowCopy.h> +#include "xAODCaloEvent/CaloCluster.h" +//#include "JetUtils/JetCaloQualityUtils.h" +//#include "JetUtils/JetCellAccessor.h" + + +TileJetMonitorAlgorithm::TileJetMonitorAlgorithm( const std::string& name, ISvcLocator* pSvcLocator ) + : AthMonitorAlgorithm(name, pSvcLocator) + , m_tileID{nullptr} + , m_tileHWID{nullptr} + , m_cabling{nullptr} +{} + + +TileJetMonitorAlgorithm::~TileJetMonitorAlgorithm() {} + + +StatusCode TileJetMonitorAlgorithm::initialize() { + + ATH_CHECK( AthMonitorAlgorithm::initialize() ); + + ATH_MSG_INFO("in initialize()"); + + ATH_CHECK( detStore()->retrieve(m_tileID) ); + ATH_CHECK( detStore()->retrieve(m_tileHWID) ); + + //=== get TileCablingSvc + ServiceHandle<TileCablingSvc> cablingSvc("TileCablingSvc", name()); + ATH_CHECK( cablingSvc.retrieve() ); + + //=== cache pointers to cabling helpers + m_cabling = cablingSvc->cablingService(); + + if (!m_cabling) { + ATH_MSG_ERROR( "Pointer to TileCablingService is zero: " << m_cabling); + return StatusCode::FAILURE; + } + + ATH_MSG_INFO("value of m_doJetCleaning: " << m_doJetCleaning); + + //=== get TileBadChanTool + ATH_MSG_DEBUG("::Retrieving tile bad channel tool"); + ATH_CHECK(m_tileBadChanTool.retrieve()); + ATH_MSG_DEBUG("::Retrieved tile bad channel tool"); + + if (m_doJetCleaning) { + ATH_MSG_DEBUG("::initializing JVT updater"); + ATH_CHECK(m_jvt.retrieve()); + ATH_MSG_DEBUG("::initialized JVT updater"); + + ATH_MSG_DEBUG("::initializing JetCleaningTool"); + ATH_CHECK(m_jetCleaningTool.retrieve()); + ATH_CHECK(m_eventCleaningTool.retrieve()); + ATH_MSG_DEBUG("::initialized JetCleaningTool"); + } else { + m_jvt.disable(); + m_jetCleaningTool.disable(); + m_eventCleaningTool.disable(); + } + + ATH_CHECK( m_jetContainerKey.initialize() ); + + return StatusCode::SUCCESS; +} + + +StatusCode TileJetMonitorAlgorithm::fillHistograms( const EventContext& ctx ) const { + + // In case you want to measure the execution time + auto timer = Monitored::Timer("TIME_execute"); + + if (!isGoodEvent(ctx)) { + ATH_MSG_DEBUG("::fillHistograms(), event skipped "); + return StatusCode::SUCCESS; + } else { + ATH_MSG_DEBUG("::fillHistograms(), event accepted "); + } + + SG::ReadHandle<xAOD::JetContainer> jetContainer(m_jetContainerKey); + if (!jetContainer.isValid()) { + ATH_MSG_WARNING("Can't retrieve Jet Container: " << m_jetContainerKey.key()); + return StatusCode::SUCCESS; + } + + uint32_t lumiBlock = GetEventInfo(ctx)->lumiBlock(); + + ATH_MSG_VERBOSE("::fillHistograms(), lumiblock " << lumiBlock); + + std::set<Identifier> usedCells; // cells already used in the given event + + int iJet = 0; + for (const xAOD::Jet* jet : *jetContainer) { + if ((jet->pt() > m_jetPtMin) && (fabs(jet->eta()) < m_jetEtaMax)) { + if (isGoodJet(*jet)) { + ATH_MSG_DEBUG("::fillHistograms, jet " << iJet + << ", eta " << jet->eta() + << ", phi " << jet->phi() + << ", constituents " << jet->numConstituents()); + CHECK(fillTimeHistograms(*jet, lumiBlock, usedCells)); + } else { + ATH_MSG_DEBUG("::fillHistogram, BAD jet " << iJet + << ", eta " << jet->eta() + << ", phi " << jet->phi() + << ", constituents " << jet->numConstituents()); + } + } + ++iJet; + } + + + fill("TileJetMonExecuteTime", timer); + + ATH_MSG_VERBOSE("::fillHistograms(), end-of-loops "); + + return StatusCode::SUCCESS; +} + + +/*---------------------------------------------------------*/ +StatusCode TileJetMonitorAlgorithm::fillTimeHistograms(const xAOD::Jet& jet, uint32_t lumiBlock, std::set<Identifier>& usedCells) const { +/*---------------------------------------------------------*/ + + ATH_MSG_VERBOSE( "in fillTimeHistograms()" ); + + if( jet.numConstituents() == 0 || !jet.getConstituents().isValid()) return StatusCode::SUCCESS; + + int cellIndex(-1); + + ToolHandle<GenericMonitoringTool> tileJetChannTimeDQTool = getGroup("TileJetChanTimeDQ"); + + std::array<std::string, 2> gainName{"LG", "HG"}; + std::array<std::string, 5> partitionName{"AUX", "LBA", "LBC", "EBA", "EBC"}; + + for (const xAOD::JetConstituent* jet_constituent : jet.getConstituents()) { + if( jet_constituent->type() == xAOD::Type::CaloCluster ){ + const xAOD::CaloCluster* calo_cluster = static_cast<const xAOD::CaloCluster*>(jet_constituent->rawConstituent()); + if (calo_cluster && calo_cluster->getCellLinks()) { + for (const CaloCell* cell : *calo_cluster) { + cellIndex++; + if (cell->caloDDE()->getSubCalo() == CaloCell_ID::TILE) { // a Tile Cell + ATH_MSG_DEBUG("Cell " << cellIndex << " IS TILECAL !!"); + const TileCell *tilecell = (TileCell*) cell; + Identifier id = tilecell->ID(); + if (usedCells.find(id) == usedCells.end()) { + usedCells.insert(id); + } else { + continue; + } + // int section= m_tileID->section(id); + // int module = m_tileID->module(id); // ranges 0..63 + auto module = Monitored::Scalar<int>("module", m_tileID->module(id)); + int sample = m_tileID->sample(id); // ranges 0..3 (A, BC, D, E) + int ros1 = 0; + int ros2 = 0; + int chan1 = -1; + int chan2 = -1; + uint32_t bad1 = 0; + uint32_t bad2 = 0; + int gain1 = tilecell->gain1(); + int gain2 = tilecell->gain2(); + unsigned int qbit1 = tilecell->qbit1(); + unsigned int qbit2 = tilecell->qbit2(); + + const CaloDetDescrElement * caloDDE = tilecell->caloDDE(); + IdentifierHash hash1 = caloDDE->onl1(); + if (hash1 != TileHWID::NOT_VALID_HASH) { + HWIdentifier adc_id = m_tileHWID->adc_id(hash1, gain1); + ros1 = m_tileHWID->ros(adc_id); + chan1 = m_tileHWID->channel(adc_id); + bad1 = m_tileBadChanTool->encodeStatus(m_tileBadChanTool->getAdcStatus(adc_id)); + } + + // How is it here with partition? D0 spans two partitions.... + // It should be ok to treat it in this way: + IdentifierHash hash2 = caloDDE->onl2(); + if (hash2 != TileHWID::NOT_VALID_HASH) { + HWIdentifier adc_id = m_tileHWID->adc_id(hash2, gain2); + ros2 = m_tileHWID->ros(adc_id); + chan2 = m_tileHWID->channel(adc_id); + bad2 = m_tileBadChanTool->encodeStatus(m_tileBadChanTool->getAdcStatus(adc_id)); + } + + bool is_good1 = isGoodChannel(ros1, module, chan1, bad1, qbit1, id); + bool is_good2 = isGoodChannel(ros2, module, chan2, bad2, qbit2, id); + float ene1 = is_good1 ? tilecell->ene1() : -1; + float ene2 = is_good2 ? tilecell->ene2() : -1; + + ATH_MSG_DEBUG(".... " << TileCalibUtils::getDrawerString(ros1, module) + << ", ch1 " << chan1 + << ", ch2 " << chan2 + << ", qbit " << qbit1 << "/" << qbit2 + << ", is_bad " << bad1 << "/" << bad2 + << ", isGood " << is_good1 + << "/" << is_good2 + << ", ene " << tilecell->energy()); + /* + Now really fill the histograms time vs lumiblock and 1dim time + */ + + // first channel + if (is_good1 && (ene1 > m_energyChanMin) && (ene1 < m_energyChanMax)) { + if (m_do1DHistograms) { + std::string name = TileCalibUtils::getDrawerString(ros1, module) + "_ch_" + std::to_string(chan1) + "_1d"; + auto channelTime = Monitored::Scalar<float>(name, tilecell->time1()); + fill("TileJetChanTime1D", channelTime); + } + + if (m_do2DHistograms) { + ATH_MSG_WARNING("This histograms are not implemented yet!"); + } + + // info for DQ histograms + auto moduleDQ = Monitored::Scalar<int>("module" + partitionName[ros1], module + 1); + auto channelDQ = Monitored::Scalar<int>("channel" + partitionName[ros1], chan1); + auto timeDQ = Monitored::Scalar<float>("time" + partitionName[ros1], tilecell->time1()); + Monitored::fill(tileJetChannTimeDQTool, moduleDQ, channelDQ, timeDQ); + + // general histograms, only require non-affected channels + if (bad1 < 2) { + ATH_MSG_DEBUG( "Filling in time1 for " << TileCalibUtils::getDrawerString(ros1, module) + << ", ch " << chan1 + << ", ene " << ene1 + << ", LB " << lumiBlock + << ", time: " << tilecell->time1()); + + std::string name("channelTime" + partitionName[ros1]); + auto channelTime = Monitored::Scalar<float>(name, tilecell->time1()); + fill("TileJetChanTime", channelTime); + + if ((ros1 >= 3) && (sample < TileID::SAMP_E)) { + std::string nameNoScint("channelTime" + partitionName[ros1] + "_NoScint"); + auto channelTimeNoScint = Monitored::Scalar<float>(nameNoScint, tilecell->time1()); + fill("TileJetChanTime", channelTimeNoScint); + } + } + } + + // second channel + if (is_good2 && (ene2 > m_energyChanMin) && (ene2 < m_energyChanMax)) { + if (m_do1DHistograms) { + std::string name = TileCalibUtils::getDrawerString(ros2, module) + "_ch_" + std::to_string(chan2) + "_1d"; + auto channelTime = Monitored::Scalar<float>(name, tilecell->time2()); + fill("TileJetChanTime1D", channelTime); + } + + if (m_do2DHistograms) { + ATH_MSG_WARNING("This histograms are not implemented yet!"); + } + + // info for DQ histograms + auto moduleDQ = Monitored::Scalar<int>("module" + partitionName[ros2], module + 1); + auto channelDQ = Monitored::Scalar<int>("channel" + partitionName[ros2], chan2); + auto timeDQ = Monitored::Scalar<float>("time" + partitionName[ros2], tilecell->time2()); + Monitored::fill(tileJetChannTimeDQTool, moduleDQ, channelDQ, timeDQ); + + // general histograms, only require non-affected channels + if (bad2 < 2) { + ATH_MSG_DEBUG( "Filling in time2 for " << TileCalibUtils::getDrawerString(ros2, module) + << ", ch " << chan2 + << ", ene " << ene2 + << ", LB " << lumiBlock + << ", time: " << tilecell->time2() + <<" (qbit2 " << qbit2 << ", ch1 " << chan1 << ", ene1 " << ene1 << ", bad1 " << bad1 << ", qbit1 " << qbit1 << ")" ); + + std::string name("channelTime" + partitionName[ros2]); + auto channelTime = Monitored::Scalar<float>(name, tilecell->time2()); + fill("TileJetChanTime", channelTime); + + if ((ros2 > 2) && (sample < TileID::SAMP_E)) { + std::string nameNoScint("channelTime" + partitionName[ros2] + "_NoScint"); + auto channelTimeNoScint = Monitored::Scalar<float>(nameNoScint, tilecell->time2()); + fill("TileJetChanTime", channelTimeNoScint); + } + } + } + + /* + Now filling the cell-based histograms, + HG-HG and LG-LG combinations only + */ + if ((is_good1) && (is_good2) && (gain1 == gain2)) { + + if (m_doEnergyDiffHistograms && (tilecell->energy() > m_energyDiffThreshold)) { + // EneDiff histograms + + int evenChannnel = (chan1 % 2 == 0) ? chan1 : chan2; + std::string name = TileCalibUtils::getDrawerString(ros1, module) + "_enediff_" + + gainName[gain1] + "_ch1_" + std::to_string(evenChannnel); + auto energyDifference = Monitored::Scalar<float>(name, tilecell->eneDiff() / tilecell->energy()); + fill("TileJetEnergyDiff", energyDifference); + } + + if ((bad1 < 2) && (bad2 < 2)) { + + // cell-time histograms, only overall, require not affected channels + if ((ros1 < 3) || (sample != TileID::SAMP_E)) { + int index = findIndex(gain1, tilecell->energy()); + ATH_MSG_DEBUG( "Filling in cell-time for " << TileCalibUtils::getDrawerString(ros1, module) + << ", ch1 " << chan1 + << ", ch2 " << chan2 + << ", ene " << tilecell->energy() + << ", index " << index + << ", time: " << tilecell->time()); + + std::string name("Cell_time_" + partitionName[ros1] + "_" + gainName[gain1] + "_slice_" + std::to_string(index)); + auto cellTime = Monitored::Scalar<float>(name, tilecell->time()); + fill("TileJetCellTime", cellTime); + + if (m_doEnergyProfiles) { + std::string indexName("index_" + partitionName[ros1] + "_" + gainName[gain1]); + auto energyIndex = Monitored::Scalar<float>(indexName, index); + + std::string energyName("energy_" + partitionName[ros1] + "_" + gainName[gain1]); + auto cellEnergy = Monitored::Scalar<float>(energyName, tilecell->energy()); + + fill("TileJetCellEnergyProfile", energyIndex, cellEnergy); + } else { + std::string name("Cell_ene_" + partitionName[ros1] + "_" + gainName[gain1] + "_slice_" + std::to_string(index)); + auto cellEnergy = Monitored::Scalar<float>(name, tilecell->energy()); + fill("TileJetCellEnergy", cellEnergy); + } + + } + } + } + } else { + ATH_MSG_DEBUG("Cell " << cellIndex << " is NOT Tilecal"); + } + } + } + } + } + + return StatusCode::SUCCESS; +} + + + +bool TileJetMonitorAlgorithm::isGoodChannel(int ros, int module, int channel, uint32_t bad, unsigned int qbit, Identifier id) const { + + if ((ros < 1) || (ros >= (int) TileCalibUtils::MAX_ROS)) return false; // invalid partition + + if ((module < 0) || (module >= (int) TileCalibUtils::MAX_DRAWER)) return false; // invalid module number + + if (m_cabling->channel2hole(ros, channel) <= 0) + return false; // non-existing PMT (empty hole) + + if (((qbit & TileCell::MASK_BADCH) != 0) || // masked as bad + ((qbit & TileCell::MASK_TIME) != TileCell::MASK_TIME) || // flagged + ((qbit & TileCell::MASK_ALGO) == TileFragHash::OptFilterDsp)) // in DSP + + return false; + /* + bad is the status in the DB (see http://alxr.usatlas.bnl.gov/lxr-stb6/source/atlas/TileCalorimeter/TileConditions/src/TileBadChanTool.cxx#390). + Meaning: + 0 = good, 1 = noisy, 2 = affected, 3 = bad, 4 = otherwise + */ + if (bad > 2) return false; + + /* + Now check for special C10, merged E1, E4' + C10 spec is ok only if channel = 5 (i.e. pmt=6). The other is pmt=5 + E1 merged and E4' should be dropped if channel = 12 (i.e. pmt=13) + */ + return ((( channel != 4) && (channel != 12)) || m_cabling->TileGap_connected(id)); +} + + +bool TileJetMonitorAlgorithm::isGoodEvent(const EventContext& ctx) const { + /* check for errors in LAr and Tile, see https://twiki.cern.ch/twiki/bin/viewauth/Atlas/DataPreparationCheckListForPhysicsAnalysis + */ + if (! m_doEventCleaning) return true; + + ATH_MSG_DEBUG("::isGoodEvent()...."); + + SG::ReadHandle<xAOD::EventInfo> eventInfo = GetEventInfo(ctx); + + if (eventInfo->errorState(xAOD::EventInfo::LAr) == xAOD::EventInfo::Error) return false; + if (eventInfo->errorState(xAOD::EventInfo::Tile) == xAOD::EventInfo::Error) return false; + + /* see https://twiki.cern.ch/twiki/bin/view/AtlasProtected/HowToCleanJets2017 + */ + if (! m_doJetCleaning) return true; + + const xAOD::JetContainer* jetContainer = nullptr; + if (! evtStore()->retrieve(jetContainer, "AntiKt4EMTopoJets").isSuccess()){ + ATH_MSG_INFO("Cannot retrieve AntiKt4EMTopoJets. However, returning true."); + return true; + } + + auto jetsSC = xAOD::shallowCopyContainer(*jetContainer); + std::unique_ptr< xAOD::JetContainer > jetsCopy(jetsSC.first); + std::unique_ptr< xAOD::ShallowAuxContainer > jetsCopyAux(jetsSC.second); + int iJet = 0; + for (auto jet : *jetsCopy) { + ATH_MSG_DEBUG("Jet " << iJet << ", pT " << jet->pt()/1000.0 << " GeV, eta " + << jet->eta()); + jet->auxdata<char>(m_jvtDecorator) = passesJvt(*jet); + jet->auxdata<char>(m_orDecorator) = true; + ATH_MSG_DEBUG("... done with jet " << iJet); + ++iJet; + } + + bool accept = m_eventCleaningTool->acceptEvent(&*jetsCopy); + + return accept; +} + + +bool TileJetMonitorAlgorithm::passesJvt(const xAOD::Jet& jet) const { + + if (jet.pt() > m_jetPtMin + && jet.pt() < m_jetPtMax + && fabs(jet.getAttribute<float>("DetectorEta")) < m_jetTrackingEtaLimit + && m_jvt->updateJvt(jet) < m_jvtThreshold) + return false; + return true; + +} + +bool TileJetMonitorAlgorithm::isGoodJet(const xAOD::Jet& jet) const { + + if (! m_doJetCleaning) return true; + if (jet.pt() < m_jetPtMin) return false; + if (! passesJvt(jet)) return false; + if (! m_jetCleaningTool->keep(jet)) return false; + return true; + +} + +unsigned int TileJetMonitorAlgorithm::findIndex(const int gain, const float energy) const { + + if (gain == 1) { + return (std::upper_bound(m_cellEnergyUpperLimitsHG.begin(), m_cellEnergyUpperLimitsHG.end(), energy) + - m_cellEnergyUpperLimitsHG.begin()); + } else { + return (std::upper_bound(m_cellEnergyUpperLimitsLG.begin(), m_cellEnergyUpperLimitsLG.end(), energy) + - m_cellEnergyUpperLimitsLG.begin()); + } + +} diff --git a/TileCalorimeter/TileMonitoring/src/TileJetMonitorAlgorithm.h b/TileCalorimeter/TileMonitoring/src/TileJetMonitorAlgorithm.h new file mode 100644 index 0000000000000000000000000000000000000000..96c4ec49c263ad79be016634ba592c6e27de21a1 --- /dev/null +++ b/TileCalorimeter/TileMonitoring/src/TileJetMonitorAlgorithm.h @@ -0,0 +1,93 @@ +/* + Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef TILEMONITORING_TILEJETMONITORALGORITHM +#define TILEMONITORING_TILEJETMONITORALGORITHM + +#include "AthenaMonitoring/AthMonitorAlgorithm.h" +#include "AthenaMonitoring/Monitored.h" + +#include "TileConditions/ITileBadChanTool.h" + +#include "StoreGate/ReadHandleKey.h" +#include "xAODJet/JetContainer.h" + +// JVT +#include "JetInterface/IJetUpdateJvt.h" +// Jet cleaning +#include "JetInterface/IJetSelector.h" +// Event cleaning +#include "JetSelectorTools/IEventCleaningTool.h" + + +class TileCablingService; +class TileID; +class TileHWID; + +/** @class TileJetMonitorAlgorithm + * @brief Class for Tile Jet based monitoring + */ + +class TileJetMonitorAlgorithm : public AthMonitorAlgorithm { + + public: + + TileJetMonitorAlgorithm( const std::string& name, ISvcLocator* pSvcLocator ); + virtual ~TileJetMonitorAlgorithm(); + virtual StatusCode initialize() override; + virtual StatusCode fillHistograms( const EventContext& ctx ) const override; + + private: + + StatusCode fillTimeHistograms(const xAOD::Jet& jet, uint32_t lumiBlock, std::set<Identifier>& usedCells) const; + unsigned int findIndex(const int gain, const float energy) const; + bool isGoodChannel(int part, int module, int channel, uint32_t bad, unsigned int qbit, Identifier id) const; + bool passesJvt(const xAOD::Jet& jet) const; + bool isGoodJet(const xAOD::Jet& jet) const; + bool isGoodEvent(const EventContext& ctx) const; + + Gaudi::Property<float> m_jetPtMin{this, "JetPtMin", 20000., "Threshold in MeV"}; + Gaudi::Property<float> m_jetPtMax{this, "JetPtMax", 120000, ""}; + Gaudi::Property<float> m_jetEtaMax{this, "JetEtaMax", 1.6, ""}; + Gaudi::Property<float> m_energyChanMin{this, "ChannelEnergyMin", 2000, ""}; + Gaudi::Property<float> m_energyChanMax{this, "ChannelEnergyMax", 4000, ""}; + Gaudi::Property<bool> m_do1DHistograms{this, "Do1DHistograms", false, ""}; + Gaudi::Property<bool> m_do2DHistograms{this, "Do2DHistograms", false, ""}; + Gaudi::Property<bool> m_doEnergyDiffHistograms{this, "DoEnergyDiffHistograms", false, ""}; + Gaudi::Property<float> m_energyDiffThreshold{this, "EnergyDiffThreshold", 2000, ""}; + Gaudi::Property<bool> m_doEnergyProfiles{this, "DoEnergyProfiles", true, ""}; + Gaudi::Property<bool> m_doEventCleaning{this, "DoEventCleaning", true, ""}; + Gaudi::Property<bool> m_doJetCleaning{this, "DoJetCleaning", false, ""}; + Gaudi::Property<float> m_jetTrackingEtaLimit{this, "JetTrackingEtaLimit", 2.4, ""}; + Gaudi::Property<float> m_jvtThreshold{this, "JvtThreshold", 0.59, ""}; + Gaudi::Property<std::string> m_jvtDecorator{this, "JvtDecorator", "passJvt", ""}; + Gaudi::Property<std::string> m_orDecorator{this, "OrDecorator", "passOR", ""}; + Gaudi::Property<std::vector<float>> m_cellEnergyUpperLimitsHG{this, + "CellEnergyUpperLimitsHG", {}, "Energy upper limits of the HG cell-time histograms"}; + Gaudi::Property<std::vector<float>> m_cellEnergyUpperLimitsLG{this, + "CellEnergyUpperLimitsLG", {}, "Energy upper limits of the LG cell-time histograms"}; + + ToolHandle<ITileBadChanTool> m_tileBadChanTool{this, + "TileBadChanTool", "TileBadChanTool", "Tile bad channel tool"}; + + // JVT + ToolHandle<IJetUpdateJvt> m_jvt{this, "JVT", "", ""}; + + // event/jet cleaning + ToolHandle<IJetSelector> m_jetCleaningTool{this, "JetCleaningTool", "", ""}; + ToolHandle<ECUtils::IEventCleaningTool> m_eventCleaningTool{this, "EventCleaningTool", "", ""}; + + + SG::ReadHandleKey<xAOD::JetContainer> m_jetContainerKey{this, + "JetContainer", "AntiKt4EMTopoJets", "Jet container for monitoring"}; + + + const TileID* m_tileID; + const TileHWID* m_tileHWID; + + const TileCablingService* m_cabling; //!< TileCabling instance + +}; + +#endif // TILEMONITORING_TILEJETMONITORALGORITHM diff --git a/TileCalorimeter/TileMonitoring/src/components/TileMonitoring_entries.cxx b/TileCalorimeter/TileMonitoring/src/components/TileMonitoring_entries.cxx index b9951b659034397bc07a997b10f391440bde59f4..a60e851c537bd182ad055c4fa5d7435b8cfa6a3f 100644 --- a/TileCalorimeter/TileMonitoring/src/components/TileMonitoring_entries.cxx +++ b/TileCalorimeter/TileMonitoring/src/components/TileMonitoring_entries.cxx @@ -22,6 +22,7 @@ #include "TileMonitoring/TileTBBeamMonTool.h" #include "TileMonitoring/TileTBMonTool.h" #include "TileMonitoring/TileTBCellMonTool.h" +#include "../TileJetMonitorAlgorithm.h" DECLARE_COMPONENT( TileFatherMonTool ) DECLARE_COMPONENT( TilePaterMonTool ) @@ -47,4 +48,5 @@ DECLARE_COMPONENT( TileRawChannelNoiseMonTool ) DECLARE_COMPONENT( TileTBBeamMonTool ) DECLARE_COMPONENT( TileTBMonTool ) DECLARE_COMPONENT( TileTBCellMonTool ) +DECLARE_COMPONENT( TileJetMonitorAlgorithm )