diff --git a/Control/CalypsoExample/Reconstruction/CMakeLists.txt b/Control/CalypsoExample/Reconstruction/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..b5644c450e793fb6fdcb32f502b37c4d7a54231d --- /dev/null +++ b/Control/CalypsoExample/Reconstruction/CMakeLists.txt @@ -0,0 +1,23 @@ +################################################################################ +# Package: Reconstruction +################################################################################ + +# Declare the package name: +atlas_subdir( Reconstruction ) + +# Component(s) in the package: +#atlas_add_component( GeoModelTest +# src/GeoModelTestAlg.cxx +# src/components/GeoModelTest_entries.cxx +# INCLUDE_DIRS ${GEOMODEL_INCLUDE_DIRS} +# LINK_LIBRARIES ${GEOMODEL_LIBRARIES} AthenaBaseComps GeoModelFaserUtilities ScintReadoutGeometry TrackerReadoutGeometry MagFieldInterfaces MagFieldElements MagFieldConditions ) + +#atlas_add_test( ReconstructionTest +# SCRIPT python ${CMAKE_CURRENT_SOURCE_DIR}/python/GeoModelTestConfig.py +# PROPERTIES WORKING_DIRECTORY ${CMAKE_BINARY_DIR} +# PROPERTIES TIMEOUT 300 ) + +# Install files from the package: +#atlas_install_joboptions( share/*.py ) +#atlas_install_python_modules( python/*.py ) +atlas_install_scripts( scripts/*.sh scripts/*.py ) diff --git a/Control/CalypsoExample/Reconstruction/README.md b/Control/CalypsoExample/Reconstruction/README.md new file mode 100644 index 0000000000000000000000000000000000000000..dfca9d50a0658b07db71361c1a37329dd0b68c17 --- /dev/null +++ b/Control/CalypsoExample/Reconstruction/README.md @@ -0,0 +1,44 @@ +This package stores the production reconstruction scripts. + +The following scripts run production reconstruction jobs: +* `faser_reco.py` - this is the python script for the main reconstruction job. + +* `submit_faser_reco.sh` - bash script to set up the environment and run the faser_reco job. This can be called from condor or other batch scheduling systems. + + +Production reco is intended to be run from a specific git tag. + +To see the available tags, use +``` +git tag -l "reco/*" +``` +and to check out a tag simply use the tag name in the checkout command +``` +git checkout reco/r0001 +``` + +To make a production reco tag, the easiest is to just tag a commit +``` +git tag -a reco/r0001 9fceb02 -m "tag message" +``` + +This needs to be pushed to the upstream master +``` +git push origin reco/r0001 +git push upstream reco/r0001 +``` +If you need to delete a tag from the repository +``` +git push origin --delete reco/r0001 +``` + +To checkout a branch based on a tag (to be able to make fixes): +``` +git checkout -b mybranch reco/r0001 +``` + +To find the updated files between two tags, and the changes in those files +``` +git diff tag1 tag2 --stat +git diff tag1 tag2 -- filename +``` diff --git a/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py b/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py new file mode 100755 index 0000000000000000000000000000000000000000..a5c03abb9b11173dcde1b373d5e66d6963dd7255 --- /dev/null +++ b/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py @@ -0,0 +1,175 @@ +#!/usr/bin/env python +# +# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +# Run with: +# ./faser_reco.py filepath [runtype] +# +# filepath - fully qualified path, including url if needed, to the input raw data file +# example: "root://hepatl30//atlas/local/torrence/faser/commissioning/TestBeamData/Run-004150/Faser-Physics-004150-00000.raw" +# +# runtype - optional flag to specify the data type (TI12Data or TestBeamData). +# In a normal file system location, this will be extracted from the directory name, +# but runtype will override this assignment. +# +import sys +import argparse + +parser = argparse.ArgumentParser(description="Run FASER reconstruction") + +parser.add_argument("file_path", + help="Fully qualified path of the raw input file") +parser.add_argument("run_type", nargs="?", default="", + help="Specify run type (if it can't be parsed from path)") +parser.add_argument("-r", "--reco", default="", + help="Specify reco tag (to append to output filename)") +parser.add_argument("-n", "--nevents", type=int, default=-1, + help="Specify number of events to process (default: all)") + +args = parser.parse_args() + +from pathlib import Path + +filepath=Path(args.file_path) + +# runtype has been provided +if len(args.run_type) > 0: + runtype=args.run_type + +# Extract runtype from path +# Should be directory above run +# i.e.: TestBeamData/Run-004150/Faser-Physics-004150-00000.raw" +else: + if len(filepath.parts) < 3: + print("Can't determine run type from path - specify on command line instead") + sys.exit(-1) + + runtype = filepath.parts[-3] + +print(f"Starting reconstruction of {filepath.name} with type {runtype}") +if args.nevents > 0: + print(f"Reconstructing {args.nevents} events by command-line option") + +# Start reconstruction + +from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator +from AthenaConfiguration.ComponentFactory import CompFactory +from AthenaCommon.Constants import VERBOSE, INFO + +from AthenaCommon.Configurable import Configurable +from CalypsoConfiguration.AllConfigFlags import ConfigFlags + +Configurable.configurableRun3Behavior = True + +# Flags for this job +ConfigFlags.Input.isMC = False # Needed to bypass autoconfig +ConfigFlags.IOVDb.DatabaseInstance = "OFLP200" # Use MC conditions for now + +ConfigFlags.Input.ProjectName = "data20" +ConfigFlags.GeoModel.Align.Dynamic = False + +# TI12 Cosmics geometry +if runtype == "TI12Data": + ConfigFlags.GeoModel.FaserVersion = "FASER-01" + ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-01" + +# Testbeam setup +elif runtype == "TestBeamData": + ConfigFlags.GeoModel.FaserVersion = "FASER-TB00" + ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-02" + +else: + print("Invalid run type found:", runtype) + print("Specify correct type or update list") + sys.exit(-1) + + +# Must use original input string here, as pathlib mangles double // in path names +ConfigFlags.Input.Files = [ args.file_path ] + +filestem = filepath.stem +if len(args.reco) > 0: + filestem += f"-{args.reco}" + +ConfigFlags.addFlag("Output.xAODFileName", f"{filestem}-xAOD.root") +ConfigFlags.Output.ESDFileName = f"{filestem}-ESD.root" + +# +# Play around with this? +# ConfigFlags.Concurrency.NumThreads = 2 +# ConfigFlags.Concurrency.NumConcurrentEvents = 2 +ConfigFlags.lock() + +# +# Configure components +from CalypsoConfiguration.MainServicesConfig import MainServicesCfg +from AthenaPoolCnvSvc.PoolWriteConfig import PoolWriteCfg + +acc = MainServicesCfg(ConfigFlags) +acc.merge(PoolWriteCfg(ConfigFlags)) + +# +# Set up RAW data access +from FaserByteStreamCnvSvc.FaserByteStreamCnvSvcConfig import FaserByteStreamCnvSvcCfg +acc.merge(FaserByteStreamCnvSvcCfg(ConfigFlags)) + +# +# Needed, or move to MainServicesCfg? +from FaserGeoModel.FaserGeoModelConfig import FaserGeometryCfg +acc.merge(FaserGeometryCfg(ConfigFlags)) + +# Set up algorithms +from WaveRecAlgs.WaveRecAlgsConfig import WaveformReconstructionCfg +acc.merge(WaveformReconstructionCfg(ConfigFlags)) + +# Tracker clusters +from TrackerPrepRawDataFormation.TrackerPrepRawDataFormationConfig import FaserSCT_ClusterizationCfg +acc.merge(FaserSCT_ClusterizationCfg(ConfigFlags)) + +# ... try SpacePoints +#from TrackerSpacePointFormation.TrackerSpacePointFormationConfig import TrackerSpacePointFinderCfg +#acc.merge(TrackerSpacePointFinderCfg(ConfigFlags)) + +# Try Dave's fitter +from TrackerClusterFit.TrackerClusterFitConfig import ClusterFitAlgCfg +acc.merge(ClusterFitAlgCfg(ConfigFlags)) + +# +# Configure output +from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg +itemList = [ "xAOD::EventInfo#*" + , "xAOD::EventAuxInfo#*" + , "xAOD::FaserTriggerData#*" + , "xAOD::FaserTriggerDataAux#*" + , "FaserSCT_RDO_Container#*" + , "Tracker::FaserSCT_ClusterContainer#*" + #, "Tracker::SCT_SpacePointContainer#*" + #, "Tracker::SCT_SpacePointOverlapCollection#*" + , "TrackCollection#*" +] +acc.merge(OutputStreamCfg(ConfigFlags, "xAOD", itemList)) + +# Waveform reconstruction +from WaveRecAlgs.WaveRecAlgsConfig import WaveformReconstructionOutputCfg +acc.merge(WaveformReconstructionOutputCfg(ConfigFlags)) + +# Check what we have +print( "Writing out xAOD objects:" ) +print( acc.getEventAlgo("OutputStreamxAOD").ItemList ) + +# Configure verbosity +# ConfigFlags.dump() +# logging.getLogger('forcomps').setLevel(VERBOSE) +# acc.foreach_component("*").OutputLevel = VERBOSE +acc.foreach_component("*").OutputLevel = INFO +acc.foreach_component("*ClassID*").OutputLevel = INFO +# log.setLevel(VERBOSE) + +#acc.getService("FaserByteStreamInputSvc").DumpFlag = True +#acc.getService("FaserEventSelector").OutputLevel = VERBOSE +#acc.getService("FaserByteStreamInputSvc").OutputLevel = VERBOSE +#acc.getService("FaserByteStreamCnvSvc").OutputLevel = VERBOSE +#acc.getService("FaserByteStreamAddressProviderSvc").OutputLevel = VERBOSE +acc.getService("MessageSvc").Format = "% F%40W%S%7W%R%T %0W%M" + +# Execute and finish +sys.exit(int(acc.run(maxEvents=args.nevents).isFailure())) diff --git a/Control/CalypsoExample/Reconstruction/scripts/submit_faser_reco.sh b/Control/CalypsoExample/Reconstruction/scripts/submit_faser_reco.sh new file mode 100755 index 0000000000000000000000000000000000000000..6bb4b895732639182709dc0d237458704ddeb4fe --- /dev/null +++ b/Control/CalypsoExample/Reconstruction/scripts/submit_faser_reco.sh @@ -0,0 +1,130 @@ +#!/bin/bash +# Used with a condor file to submit to vanilla universe +# +# Usage: +# submit_faser_reco.sh file_path [release_directory] [working_directory] [nevents] +# +# file_path - full file name (with path) +# release_directory - optional path to release install directory (default pwd) +# working_directory - optional path to output directory location (default pwd) +# nevents - optional number of events to process (default: -1 - all) +# +# The release directory must already be set up +# (so an unqualified asetup can set up the release properly) +# +# Script will use git describe to find the release tag. +# If this matches reco/r???? it will be passed to the reco job +# +#---------------------------------------- +# +# Parse command-line options +file_path=${1} +release_directory=${2} +working_directory=${3} +nevents=${4} +# +# Set defaults if arguments aren't provided +if [ -z "$file_path" ] +then + echo "No file_path specified!" + exit 1 +fi +# +if [ -z "$release_directory" ] +then + release_directory=`pwd` +fi +# +if [ -z "$working_directory" ] +then + working_directory=`pwd` +fi +# +if [ -z "$nevents" ] +then + nevents="-1" +fi +# +starting_directory=`pwd` +# +# Now extract the run number and file stem +# +# First, get the filename +file_name=$(basename "$file_path") +# +# Now split based on '.' to get stem +defaultIFS=$IFS +IFS='.' +read file_stem ext <<< "$file_name" +# +# Finally extract the run number +IFS='-' +# Read the split words into an array based on delimiter +read faser type run_number segment <<< "$file_stem" +# +# Set the IFS delimeter back or else echo doesn't work... +IFS=$defaultIFS +# +# Make output directory if needed +output_directory="$working_directory/Run-$run_number" +mkdir -p "$output_directory" +# +# This magic redirects everything in this script to our log file +exec >& "$output_directory/$file_stem.log" +echo `date` - $HOSTNAME +echo "File: $file_name" +echo "Release: $release_directory" +echo "Output: $output_directory" +echo "Starting: $starting_directory" +# +# Set up the release (do this automatically)? +export ATLAS_LOCAL_ROOT_BASE=/cvmfs/atlas.cern.ch/repo/ATLASLocalRootBase +source ${ATLAS_LOCAL_ROOT_BASE}/user/atlasLocalSetup.sh +# +# Try automatic +# Always go back to the starting directory in case paths are relative +cd "$starting_directory" +cd "$release_directory" +asetup +source build/x8*/setup.sh +# +# Do this by hand +# asetup --input="$release_directory/calypso/asetup.faser" Athena,22.0.40 +# source "$release_directory/build/x8*/setup.sh" +# +# +# Try to find a release tag +cd calypso +recotag=`git describe` +if [[ "$recotag" == "reco/r"???? ]]; then + rtag=`echo "$recotag" | cut -c 6-11` + echo "Found reco tag: $rtag" +fi +# +# Move to the run directory +cd "$starting_directory" +cd "$output_directory" +# +# Remove any previous directory if it exists +#if [[ -e "$file_stem" ]]; then +# echo "Remove previous directory $file_stem" +# rm -rf "$file_stem" +#fi +# +# Make run directory +if [[ -e "$file_stem" ]]; then + echo "Directory $file_stem already exists" +else + mkdir "$file_stem" +fi +cd "$file_stem" +# +# Run job +if [[ -z "$rtag" ]]; then + faser_reco.py "--nevents=$nevents" "$file_path" +else + faser_reco.py "--nevents=$nevents" "--reco=$rtag" "$file_path" +fi +# +# Print out ending time +date diff --git a/Event/FaserByteStreamCnvSvc/python/FaserByteStreamCnvSvcConfig.py b/Event/FaserByteStreamCnvSvc/python/FaserByteStreamCnvSvcConfig.py index 633017999803b3dda6eca426bc3acdf4689c92bb..e2315ae1c6f5ff50b973e6c08377a3cee9da9ef4 100644 --- a/Event/FaserByteStreamCnvSvc/python/FaserByteStreamCnvSvcConfig.py +++ b/Event/FaserByteStreamCnvSvc/python/FaserByteStreamCnvSvcConfig.py @@ -46,6 +46,10 @@ def FaserByteStreamCnvSvcCfg(configFlags, **kwargs): from FaserByteStreamCnvSvcBase.FaserByteStreamCnvSvcBaseConfig import FaserByteStreamCnvSvcBaseCfg result.merge(FaserByteStreamCnvSvcBaseCfg(configFlags)) +# Configure WaveformByteStream converter + from WaveByteStream.WaveByteStreamConfig import WaveByteStreamCfg + result.merge(WaveByteStreamCfg(configFlags)) + # Load ByteStreamEventStorageInputSvc bsInputSvc = CompFactory.FaserByteStreamInputSvc result.addService(bsInputSvc(name = "FaserByteStreamInputSvc")) diff --git a/Waveform/WaveEventCnv/WaveByteStream/CMakeLists.txt b/Waveform/WaveEventCnv/WaveByteStream/CMakeLists.txt index 6fc3d0b95846b2e395cccbbe6938336def278e4e..bdf35c03e86b0fbac4807da3f4ace952b81e2dc1 100644 --- a/Waveform/WaveEventCnv/WaveByteStream/CMakeLists.txt +++ b/Waveform/WaveEventCnv/WaveByteStream/CMakeLists.txt @@ -13,7 +13,7 @@ atlas_subdir( WaveByteStream ) atlas_add_component( WaveByteStream src/*.cxx src/components/*.cxx - LINK_LIBRARIES AthenaKernel GaudiKernel FaserByteStreamCnvSvcBaseLib FaserEventStorageLib WaveRawEvent + LINK_LIBRARIES Identifier ScintIdentifier FaserCaloIdentifier AthenaKernel GaudiKernel FaserByteStreamCnvSvcBaseLib FaserEventStorageLib WaveRawEvent PRIVATE_LINK_LIBRARIES AthenaBaseComps ) atlas_install_python_modules( python/*.py ) diff --git a/Waveform/WaveEventCnv/WaveByteStream/python/WaveByteStreamConfig.py b/Waveform/WaveEventCnv/WaveByteStream/python/WaveByteStreamConfig.py index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..cc13137d8e59eef7801d0ffa8fa60230721b372f 100644 --- a/Waveform/WaveEventCnv/WaveByteStream/python/WaveByteStreamConfig.py +++ b/Waveform/WaveEventCnv/WaveByteStream/python/WaveByteStreamConfig.py @@ -0,0 +1,43 @@ +# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + +from AthenaConfiguration.ComponentFactory import CompFactory +from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator + +def WaveByteStreamCfg(configFlags, **kwargs): + + acc = ComponentAccumulator() + + # Decide which waveform mapping to instantiate from Geo tag + + if configFlags.Input.isMC: + # Nothing to do for MC + print("WaveByteStreamConfig.WaveByteStreamCfg - Called for isMC True, nothing to do...") + return acc + + print("WaveByteStreamConfig.WaveByteStreamCfg - Found FaserVersion: ", configFlags.GeoModel.FaserVersion,) + + # Channels are ordered front->back, bottom->top, right->left + if configFlags.GeoModel.FaserVersion == "FASER-TB00": + print(" - setting up testbeam detector") + + waveform_tool = CompFactory.RawWaveformDecoderTool("RawWaveformDecoderTool") + waveform_tool.CaloChannels = [5, 3, 1, 4, 2, 0] + waveform_tool.VetoChannels = [] + waveform_tool.TriggerChannels = [8, 9] + waveform_tool.PreshowerChannels = [6, 7] + acc.addPublicTool(waveform_tool) + + elif configFlags.GeoModel.FaserVersion == "FASER-01": + print(" - setting up TI12 detector") + + waveform_tool = CompFactory.RawWaveformDecoderTool("RawWaveformDecoderTool") + waveform_tool.CaloChannels = [1, 0, 3, 2] + waveform_tool.VetoChannels = [4, 5, 6, 7] + waveform_tool.TriggerChannels = [9, 8, 11, 10] + waveform_tool.PreshowerChannels = [12, 13] + acc.addPublicTool(waveform_tool) + + else: + print(" - unknown version: user must set up Waveform channel mapping by hand!") + + return acc diff --git a/Waveform/WaveEventCnv/WaveByteStream/src/RawWaveformDecoderTool.cxx b/Waveform/WaveEventCnv/WaveByteStream/src/RawWaveformDecoderTool.cxx index 2081244c434b40e789bf268284a7b56b13a89929..ec7a6641e3df94a8ac0b3b989e5c99303e3b62d7 100644 --- a/Waveform/WaveEventCnv/WaveByteStream/src/RawWaveformDecoderTool.cxx +++ b/Waveform/WaveEventCnv/WaveByteStream/src/RawWaveformDecoderTool.cxx @@ -20,31 +20,16 @@ RawWaveformDecoderTool::RawWaveformDecoderTool(const std::string& type, { declareInterface<RawWaveformDecoderTool>(this); + // These must be configured by job options declareProperty("CaloChannels", m_caloChannels); - m_caloChannels.push_back(0); - m_caloChannels.push_back(1); - m_caloChannels.push_back(2); - m_caloChannels.push_back(3); - declareProperty("VetoChannels", m_vetoChannels); - m_vetoChannels.push_back(4); - m_vetoChannels.push_back(5); - m_vetoChannels.push_back(6); - m_vetoChannels.push_back(7); - declareProperty("TriggerChannels", m_triggerChannels); - m_triggerChannels.push_back(8); - m_triggerChannels.push_back(9); - m_triggerChannels.push_back(10); - m_triggerChannels.push_back(11); - declareProperty("PreshowerChannels", m_preshowerChannels); - m_preshowerChannels.push_back(12); - m_preshowerChannels.push_back(13); + // Test channels is provided for conveniene, not normally used declareProperty("TestChannels", m_testChannels); - m_testChannels.push_back(14); + // Clock should always be in channel 15 declareProperty("ClockChannels", m_clockChannels); m_clockChannels.push_back(15); @@ -58,6 +43,82 @@ StatusCode RawWaveformDecoderTool::initialize() { ATH_MSG_DEBUG("RawWaveformDecoderTool::initialize()"); + + // Set up helpers + ATH_CHECK(detStore()->retrieve(m_ecalID, "EcalID")); + ATH_CHECK(detStore()->retrieve(m_vetoID, "VetoID")); + ATH_CHECK(detStore()->retrieve(m_triggerID, "TriggerID")); + ATH_CHECK(detStore()->retrieve(m_preshowerID, "PreshowerID")); + + // Take each channel list and create identifiers + + // Loop through digitizer channel lists creating Identifiers + m_identifierMap.clear(); + + // First, calorimeter. Can either be 4 or 6 channels + // Bottom row first from left to right, then top row + int index=0; + int module=0; + int row=0; + int pmt=0; + + int max_modules = m_caloChannels.size() / 2; + for (auto const& chan : m_caloChannels) { + row = index / max_modules; + module = index % max_modules; + index++; + // Only store in map if digitizer channel is valid + if (chan < 0) continue; + m_identifierMap[chan] = m_ecalID->pmt_id(row, module, pmt); + ATH_MSG_DEBUG("Mapped digitizer channel " << chan << " to calo ID: " << m_identifierMap[chan]); + } + + // Next, veto detector. Have station and plate here. + int station=0; + int plate=0; + pmt=0; + index=0; + + int max_stations=m_vetoChannels.size() / 2; + for (auto const& chan : m_vetoChannels) { + station = index / max_stations; + plate = index % max_stations; + index++; + // Only store in map if digitizer channel is valid + if (chan < 0) continue; + m_identifierMap[chan] = m_vetoID->pmt_id(station, plate, pmt); + ATH_MSG_DEBUG("Mapped digitizer channel " << chan << " to veto ID: " << m_identifierMap[chan]); + } + + // Next, trigger detector. Have pmt and plate. + pmt=0; + station=0; + index=0; + int max_plates=m_triggerChannels.size() / 2; + for (auto const& chan : m_triggerChannels) { + plate = index / max_plates; + pmt = index % max_plates; + index++; + // Only store in map if digitizer channel is valid + if (chan < 0) continue; + m_identifierMap[chan] = m_triggerID->pmt_id(station, plate, pmt); + ATH_MSG_DEBUG("Mapped dititizer channel " << chan << " to trigger ID: " << m_identifierMap[chan]); + } + + // Finally, preshower detector. + pmt=0; + station=0; + plate=0; + index=0; + for (auto const& chan : m_preshowerChannels) { + plate = index; + index++; + // Only store in map if digitizer channel is valid + if (chan < 0) continue; + m_identifierMap[chan] = m_preshowerID->pmt_id(station, plate, pmt); + ATH_MSG_DEBUG("Mapped digitizer channel " << chan << " to preshower ID: " << m_identifierMap[chan]); + } + return StatusCode::SUCCESS; } @@ -111,7 +172,7 @@ RawWaveformDecoderTool::convert(const DAQFormats::EventFull* re, ATH_MSG_DEBUG("Found valid digitizer fragment"); } - std::vector<unsigned int>* channelList; + std::vector<int>* channelList; if (key == std::string("CaloWaveforms")) { channelList = &m_caloChannels; @@ -130,7 +191,7 @@ RawWaveformDecoderTool::convert(const DAQFormats::EventFull* re, return StatusCode::FAILURE; } - for (unsigned int channel: *channelList) { + for (int channel: *channelList) { ATH_MSG_DEBUG("Converting channel "+std::to_string(channel)+" for "+key); // Check if this has data @@ -161,6 +222,11 @@ RawWaveformDecoderTool::convert(const DAQFormats::EventFull* re, << *frag); } + // Set ID if one exists (clock, for instance, doesn't have an identifier) + if (m_identifierMap.count(channel) == 1) { + wfm->setIdentifier(m_identifierMap[channel]); + } + container->push_back(wfm); // Sanity check diff --git a/Waveform/WaveEventCnv/WaveByteStream/src/RawWaveformDecoderTool.h b/Waveform/WaveEventCnv/WaveByteStream/src/RawWaveformDecoderTool.h index 252698c5fb975370945883d112a72c21ef7798a2..e2852f13772113d8e52a62c9655a18365cb80f61 100644 --- a/Waveform/WaveEventCnv/WaveByteStream/src/RawWaveformDecoderTool.h +++ b/Waveform/WaveEventCnv/WaveByteStream/src/RawWaveformDecoderTool.h @@ -13,6 +13,12 @@ #include "EventFormats/DAQFormats.hpp" #include "WaveRawEvent/RawWaveformContainer.h" +#include "Identifier/Identifier.h" +#include "FaserCaloIdentifier/EcalID.h" +#include "ScintIdentifier/VetoID.h" +#include "ScintIdentifier/TriggerID.h" +#include "ScintIdentifier/PreshowerID.h" + // This class provides conversion between bytestream and Waveform objects class RawWaveformDecoderTool : public AthAlgTool { @@ -32,13 +38,48 @@ class RawWaveformDecoderTool : public AthAlgTool { private: // List of channels to include in each container - std::vector<unsigned int> m_caloChannels; - std::vector<unsigned int> m_vetoChannels; - std::vector<unsigned int> m_triggerChannels; - std::vector<unsigned int> m_preshowerChannels; - std::vector<unsigned int> m_testChannels; - std::vector<unsigned int> m_clockChannels; + // List order must correspond to offline channel order + // All L/R designations refer to looking at the detector from + // the beam direction. + // + // In general, the ordering is layer (longitudinal), row (vertical), module (horizontal) + // Layers increase with longitudianl position downstream + // Rows increase from bottom to top + // Modules increase from right to left + // + // For all lists, use invalid channel (-1) to indicate detectors + // missing in sequence (i.e. 3 of 4 veto counters) + // + // TI12 detector: + // Calorimeter order: + // bottom right, bottom left, top right, top left + // Veto + // front to back. + // Trigger + // bottom right PMT, bottom left PMT, top right PMT, top left PMT + // Preshower + // front to back + // + // 2021 Testbeam detector: + // Calo order: + // bottom right, bottom center, bottom left, top R, top C, top L + // All others are just in order front to back + + std::vector<int> m_caloChannels; + std::vector<int> m_vetoChannels; + std::vector<int> m_triggerChannels; + std::vector<int> m_preshowerChannels; + std::vector<int> m_testChannels; + std::vector<int> m_clockChannels; + + // Identifiers keyed by digitizer channel + std::map<unsigned int, Identifier> m_identifierMap; + // ID helpers + const EcalID* m_ecalID{nullptr}; + const VetoID* m_vetoID{nullptr}; + const TriggerID* m_triggerID{nullptr}; + const PreshowerID* m_preshowerID{nullptr}; }; diff --git a/Waveform/WaveEventCnv/WaveEventAthenaPool/WaveEventAthenaPool/RawWaveform_p0.h b/Waveform/WaveEventCnv/WaveEventAthenaPool/WaveEventAthenaPool/RawWaveform_p0.h index e4721fefeec5ef7be0f945bca38a790c81151ae6..9eee456fbdae20d6b3619e7b99bb6ee640dc1ef1 100644 --- a/Waveform/WaveEventCnv/WaveEventAthenaPool/WaveEventAthenaPool/RawWaveform_p0.h +++ b/Waveform/WaveEventCnv/WaveEventAthenaPool/WaveEventAthenaPool/RawWaveform_p0.h @@ -9,6 +9,8 @@ #include <iostream> #include <iomanip> +#include "Identifier/Identifier32.h" + class RawWaveform_p0 { public: RawWaveform_p0(); @@ -17,7 +19,7 @@ class RawWaveform_p0 { friend class RawWaveformCnv_p0; private: - unsigned int m_ID; + Identifier32::value_type m_ID; unsigned int m_channel; std::vector<unsigned int> m_adc_counts; diff --git a/Waveform/WaveEventCnv/WaveEventAthenaPool/src/RawWaveformCnv_p0.cxx b/Waveform/WaveEventCnv/WaveEventAthenaPool/src/RawWaveformCnv_p0.cxx index d298c0ff390393b6b8822bcc43167eee7c89ee9d..dedc71cbac77ccbdfa48dff63013cc6c16ab97f8 100644 --- a/Waveform/WaveEventCnv/WaveEventAthenaPool/src/RawWaveformCnv_p0.cxx +++ b/Waveform/WaveEventCnv/WaveEventAthenaPool/src/RawWaveformCnv_p0.cxx @@ -9,7 +9,8 @@ RawWaveformCnv_p0::persToTrans(const RawWaveform_p0* persObj, RawWaveform* trans // Just fill available data here // Rest of it patched up in RawWaveformContainerCnv_p0 - transObj->setIdentifier(persObj->m_ID); + // Persistent object stores ID value, so instantiate Identifier here + transObj->setIdentifier(Identifier32(persObj->m_ID)); transObj->setChannel(persObj->m_channel); transObj->setCounts(persObj->m_adc_counts); @@ -23,7 +24,9 @@ RawWaveformCnv_p0::transToPers(const RawWaveform* transObj, RawWaveform_p0* pers // log << MSG::DEBUG << (*transObj) << endmsg; // log << MSG::DEBUG << "Persistent waveform (before):" << endmsg; // persObj->print(); - persObj->m_ID = transObj->identify(); + + // Save actual ID number + persObj->m_ID = transObj->identify32().get_compact(); persObj->m_channel = transObj->channel(); // Use copy here diff --git a/Waveform/WaveRawEvent/CMakeLists.txt b/Waveform/WaveRawEvent/CMakeLists.txt index 48473bf8813e27d85bd73a1c915ad16f34eabbc0..cf5334c539d23320189bf8d10a8466b9422cf440 100644 --- a/Waveform/WaveRawEvent/CMakeLists.txt +++ b/Waveform/WaveRawEvent/CMakeLists.txt @@ -16,12 +16,12 @@ atlas_add_library( WaveRawEvent INCLUDE_DIRS ${CLHEP_INCLUDE_DIRS} PRIVATE_INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} DEFINITIONS ${CLHEP_DEFINITIONS} - LINK_LIBRARIES ${CLHEP_LIBRARIES} AthAllocators AthenaKernel CxxUtils StoreGateLib SGtests + LINK_LIBRARIES ${CLHEP_LIBRARIES} AthAllocators AthenaKernel Identifier CxxUtils StoreGateLib SGtests PRIVATE_LINK_LIBRARIES ${ROOT_LIBRARIES} ScintIdentifier EventFormats ) atlas_add_dictionary( WaveRawEventDict WaveRawEvent/WaveRawEventDict.h WaveRawEvent/selection.xml INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} ${CLHEP_INCLUDE_DIRS} ${GEANT4_INCLUDE_DIRS} - LINK_LIBRARIES ${ROOT_LIBRARIES} ${CLHEP_LIBRARIES} AthAllocators CxxUtils StoreGateLib ScintIdentifier WaveRawEvent ) + LINK_LIBRARIES ${ROOT_LIBRARIES} ${CLHEP_LIBRARIES} AthAllocators Identifier CxxUtils StoreGateLib ScintIdentifier WaveRawEvent ) diff --git a/Waveform/WaveRawEvent/WaveRawEvent/RawWaveform.h b/Waveform/WaveRawEvent/WaveRawEvent/RawWaveform.h index 6ea6d9037afe587ebe21d52eb9280b23e7a397c8..3770e05513a4d6f81ca16aa08785ea0c36a459cb 100644 --- a/Waveform/WaveRawEvent/WaveRawEvent/RawWaveform.h +++ b/Waveform/WaveRawEvent/WaveRawEvent/RawWaveform.h @@ -20,10 +20,13 @@ #include <iomanip> #include "AthenaKernel/CLASS_DEF.h" +#include "Identifier/Identifier.h" +#include "Identifier/Identifier32.h" +#include "Identifier/Identifiable.h" class DigitizerDataFragment; -class RawWaveform { +class RawWaveform : public Identifiable { /////////////////////////////////////////////////////////////////// // Public methods: @@ -60,7 +63,9 @@ public: unsigned int channel() const; const std::vector<unsigned int>& adc_counts() const; - unsigned int identify() const; + // Return channel identifier + Identifier identify() const; + Identifier32 identify32() const; // some print-out: void print() const; @@ -69,7 +74,7 @@ public: {return m_ID < rhs.m_ID;} // Set functions - void setIdentifier (unsigned int id); + void setIdentifier (Identifier id); // Identifier class void setHeader (const DigitizerDataFragment* frag); void setWaveform (unsigned int channel, const std::vector<uint16_t> waveform); @@ -99,7 +104,7 @@ private: unsigned int m_channel; std::vector<unsigned int> m_adc_counts; - unsigned int m_ID; + Identifier32 m_ID; }; @@ -134,8 +139,11 @@ RawWaveform::n_samples() const { return m_samples; } inline const std::vector<unsigned int>& RawWaveform::adc_counts() const { return m_adc_counts; } -inline unsigned int -RawWaveform::identify() const { return m_ID; } +inline Identifier +RawWaveform::identify() const { return Identifier(m_ID); } + +inline Identifier32 +RawWaveform::identify32() const { return m_ID; } std::ostream &operator<<(std::ostream &out, const RawWaveform &wfm); diff --git a/Waveform/WaveRawEvent/src/RawWaveform.cxx b/Waveform/WaveRawEvent/src/RawWaveform.cxx index 13700a735c3653867a7e7ab5df8037e3566664ae..d6ef63493018955a79d42eb2c4f7fe2eebc40de2 100644 --- a/Waveform/WaveRawEvent/src/RawWaveform.cxx +++ b/Waveform/WaveRawEvent/src/RawWaveform.cxx @@ -25,8 +25,8 @@ RawWaveform::RawWaveform( ) : RawWaveform::~RawWaveform() {} void -RawWaveform::setIdentifier(unsigned int id) { - m_ID = id; +RawWaveform::setIdentifier(Identifier id) { + m_ID = id.get_identifier32(); } void diff --git a/Waveform/WaveRecAlgs/CMakeLists.txt b/Waveform/WaveRecAlgs/CMakeLists.txt index c72671bf705ace3ad1935034564982a3b60a94a2..a982b17cf28379f5c6a5d4208bc6134efc042994 100644 --- a/Waveform/WaveRecAlgs/CMakeLists.txt +++ b/Waveform/WaveRecAlgs/CMakeLists.txt @@ -9,7 +9,7 @@ atlas_subdir( WaveRecAlgs ) atlas_add_component( WaveRecAlgs src/*.cxx src/*.h src/components/*.cxx - LINK_LIBRARIES AthenaBaseComps StoreGateLib WaveRawEvent xAODFaserWaveform WaveRecToolsLib ) + LINK_LIBRARIES AthenaBaseComps StoreGateLib WaveRawEvent xAODFaserWaveform WaveRecToolsLib ScintSimEvent FaserCaloSimEvent) atlas_install_python_modules( python/*.py ) diff --git a/Waveform/WaveRecAlgs/python/WaveRecAlgsConfig.py b/Waveform/WaveRecAlgs/python/WaveRecAlgsConfig.py index 0c2d5ffc6c5d7cfff88cc74cf8d4b35f5f635ef8..39a99ed0c39d01e08c5f555e4e3efd98e0bc70a1 100644 --- a/Waveform/WaveRecAlgs/python/WaveRecAlgsConfig.py +++ b/Waveform/WaveRecAlgs/python/WaveRecAlgsConfig.py @@ -11,10 +11,18 @@ WaveformReconstructionTool = CompFactory.WaveformReconstructionTool ClockReconstructionTool = CompFactory.ClockReconstructionTool # One stop shopping for normal FASER data -def WaveformReconstructionCfg(flags): +def WaveformReconstructionCfg(flags, naive = True): """ Return all algorithms and tools for Waveform reconstruction """ acc = ComponentAccumulator() + if flags.Input.isMC and naive: + if not "TB" in flags.GeoModel.FaserVersion: + acc.merge(PseudoScintHitToWaveformRecCfg(flags, "PseudoTimingHitWaveformRecAlg", "Trigger")) + acc.merge(PseudoScintHitToWaveformRecCfg(flags, "PseudoVetoHitToWaveformRecAlg", "Veto")) + acc.merge(PseudoScintHitToWaveformRecCfg(flags, "PseudoPresehowerHitWaveformRecAlg", "Preshower")) + acc.merge(PseudoCaloHitToWaveformRecCfg(flags, "PseudoCaloHitWaveformRecAlg")) + return acc + acc.merge(WaveformClockRecCfg(flags, "ClockRecAlg")) acc.merge(WaveformHitRecCfg(flags, "VetoWaveformRecAlg", "Veto")) @@ -68,3 +76,41 @@ def WaveformReconstructionOutputCfg(flags, **kwargs): # ostream = acc.getEventAlgo("OutputStreamRDO") # ostream.TakeItemsFromInput = True # Don't know what this does return acc + +# Return configured "reconstruction" algorithm from scint SIM hits +# Specify data source (Veto, Trigger, Preshower) +def PseudoScintHitToWaveformRecCfg(flags, name="PseudoScintHitToWaveformRecAlg", source="", **kwargs): + + acc = ComponentAccumulator() + + tool = CompFactory.PseudoSimToWaveformRecTool(name=source+"PseudoSimToWaveformRecTool", **kwargs) + + kwargs.setdefault("ScintHitContainerKey", source+"Hits") + kwargs.setdefault("WaveformHitContainerKey", source+"WaveformHits") + kwargs.setdefault("PseudoSimToWaveformRecTool", tool) + + + recoAlg = CompFactory.PseudoScintHitToWaveformRecAlg(name, **kwargs) + recoAlg.PseudoSimToWaveformRecTool = tool + + acc.addEventAlgo(recoAlg) + + return acc + +# Return configured "reconstruction" algorithm from calo SIM hits +def PseudoCaloHitToWaveformRecCfg(flags, name="PseudoCaloHitToWaveformRecAlg", **kwargs): + + acc = ComponentAccumulator() + + tool = CompFactory.PseudoSimToWaveformRecTool(name="CaloPseudoSimToWaveformRecTool", **kwargs) + + kwargs.setdefault("CaloHitContainerKey", "EcalHits") + kwargs.setdefault("WaveformHitContainerKey", "CaloWaveformHits") + kwargs.setdefault("PseudoSimToWaveformRecTool", tool) + + recoAlg = CompFactory.PseudoCaloHitToWaveformRecAlg(name, **kwargs) + recoAlg.PseudoSimToWaveformRecTool = tool + + acc.addEventAlgo(recoAlg) + + return acc diff --git a/Waveform/WaveRecAlgs/share/PseudoSimToWaveformRecExample_jobOptions.py b/Waveform/WaveRecAlgs/share/PseudoSimToWaveformRecExample_jobOptions.py new file mode 100644 index 0000000000000000000000000000000000000000..cbcb5836850f8f1958d5e55e22eca9f6afe1f0e4 --- /dev/null +++ b/Waveform/WaveRecAlgs/share/PseudoSimToWaveformRecExample_jobOptions.py @@ -0,0 +1,66 @@ +#!/usr/bin/env python + +import sys + +if __name__ == "__main__": + + fileroot = "SimToRec" + naive = True + + from AthenaCommon.Logging import log, logging + from AthenaCommon.Constants import DEBUG, VERBOSE, INFO + from CalypsoConfiguration.AllConfigFlags import ConfigFlags + from AthenaCommon.Configurable import Configurable + + Configurable.configurableRun3Behavior = True + + log.setLevel(INFO) + + ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-02" # Always needed; must match FaserVersion + ConfigFlags.IOVDb.DatabaseInstance = "OFLP200" # Use MC conditions for now + ConfigFlags.Input.ProjectName = "mc21" # Needed to bypass autoconfig + ConfigFlags.Input.isMC = True # Needed to bypass autoconfig + ConfigFlags.GeoModel.FaserVersion = "FASER-TB00" # FASER geometry + ConfigFlags.Common.isOnline = False + ConfigFlags.GeoModel.Align.Dynamic = False + + ConfigFlags.Input.Files = [ + "my.HITS.pool.root" + #"/bundle/data/FASER/LC_output/BatchOutput/TestBeam/TB.Elec.8.r5.e100.SIM.root" + ] + + + ConfigFlags.addFlag("Output.xAODFileName", f"{fileroot}.xAOD.root") + ConfigFlags.Output.ESDFileName = f"{fileroot}.ESD.root" + + ConfigFlags.lock() + + from CalypsoConfiguration.MainServicesConfig import MainServicesCfg + + acc = MainServicesCfg(ConfigFlags) + + from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg + from AthenaPoolCnvSvc.PoolWriteConfig import PoolWriteCfg + + acc.merge(PoolReadCfg(ConfigFlags)) + acc.merge(PoolWriteCfg(ConfigFlags)) + + from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg + itemList = [ + "xAOD::EventInfo#*", + "xAOD::WaveformHitContainer#*", + "xAOD::WaveformHitAuxContainer#*", + ] + + acc.merge(OutputStreamCfg(ConfigFlags, "xAOD", itemList, disableEventTag=True)) + + from WaveRecAlgs.WaveRecAlgsConfig import WaveformReconstructionCfg + acc.merge(WaveformReconstructionCfg(ConfigFlags, naive)) + + #acc.foreach_component("*").OutputLevel = VERBOSE + + # Execute and finish + sc = acc.run(maxEvents=1000) + + # Success should be 0 + sys.exit(not sc.isSuccess()) diff --git a/Waveform/WaveRecAlgs/src/PseudoCaloHitToWaveformRecAlg.cxx b/Waveform/WaveRecAlgs/src/PseudoCaloHitToWaveformRecAlg.cxx new file mode 100644 index 0000000000000000000000000000000000000000..5ead65452727c36307d8b4c038812ca9c3683f47 --- /dev/null +++ b/Waveform/WaveRecAlgs/src/PseudoCaloHitToWaveformRecAlg.cxx @@ -0,0 +1,68 @@ +#include "PseudoCaloHitToWaveformRecAlg.h" + +PseudoCaloHitToWaveformRecAlg::PseudoCaloHitToWaveformRecAlg(const std::string& name, + ISvcLocator* pSvcLocator) + : AthReentrantAlgorithm(name, pSvcLocator) { + +} + +StatusCode +PseudoCaloHitToWaveformRecAlg::initialize() { + ATH_MSG_INFO(name() << "::initalize()" ); + + // Initalize tools + ATH_CHECK( m_recoTool.retrieve() ); + + // Set key to read waveform from + ATH_CHECK( m_caloHitContainerKey.initialize() ); + + // Set key to write container + ATH_CHECK( m_waveformHitContainerKey.initialize() ); + + return StatusCode::SUCCESS; +} + +StatusCode +PseudoCaloHitToWaveformRecAlg::finalize() { + ATH_MSG_INFO(name() << "::finalize()"); + + return StatusCode::SUCCESS; +} + +StatusCode +PseudoCaloHitToWaveformRecAlg::execute(const EventContext& ctx) const { + ATH_MSG_DEBUG("Executing"); + + ATH_MSG_DEBUG("Run: " << ctx.eventID().run_number() + << " Event: " << ctx.eventID().event_number()); + + // Find the input HIT collection + SG::ReadHandle<CaloHitCollection> caloHitHandle(m_caloHitContainerKey, ctx); + + ATH_CHECK( caloHitHandle.isValid() ); + ATH_MSG_DEBUG("Found ReadHandle for CaloHitCollection " << m_caloHitContainerKey); + + if (caloHitHandle->size() == 0) { + ATH_MSG_DEBUG("CaloHitCollection found with zero length!"); + return StatusCode::SUCCESS; + } + + + // Find the output waveform container + SG::WriteHandle<xAOD::WaveformHitContainer> waveformHitContainerHandle(m_waveformHitContainerKey, ctx); + ATH_CHECK( waveformHitContainerHandle.record( std::make_unique<xAOD::WaveformHitContainer>(), + std::make_unique<xAOD::WaveformHitAuxContainer>() ) ); + + ATH_MSG_DEBUG("WaveformsHitContainer '" << waveformHitContainerHandle.name() << "' initialized"); + + + // "Reconstruct" the hits + CHECK( m_recoTool->reconstruct<CaloHitCollection>(caloHitHandle.ptr(), + waveformHitContainerHandle.ptr()) ); + + + ATH_MSG_DEBUG("WaveformsHitContainer '" << waveformHitContainerHandle.name() << "' filled with "<< waveformHitContainerHandle->size() <<" items"); + + return StatusCode::SUCCESS; +} + diff --git a/Waveform/WaveRecAlgs/src/PseudoCaloHitToWaveformRecAlg.h b/Waveform/WaveRecAlgs/src/PseudoCaloHitToWaveformRecAlg.h new file mode 100644 index 0000000000000000000000000000000000000000..df8b4f2005aca4cb75e679580df72bfab037dd2a --- /dev/null +++ b/Waveform/WaveRecAlgs/src/PseudoCaloHitToWaveformRecAlg.h @@ -0,0 +1,75 @@ +#ifndef WAVERECALGS_PSEUDOCALOHITTOWAVEFORMRECALG_H +#define WAVERECALGS_PSEUDOCALOHITTOWAVEFORMRECALG_H + +// Base class +#include "AthenaBaseComps/AthReentrantAlgorithm.h" + +// Data classes +#include "FaserCaloSimEvent/CaloHitCollection.h" + +#include "xAODFaserWaveform/WaveformHit.h" +#include "xAODFaserWaveform/WaveformHitContainer.h" +#include "xAODFaserWaveform/WaveformHitAuxContainer.h" + +// Tool classes +#include "WaveRecTools/IPseudoSimToWaveformRecTool.h" + +// Handles +#include "StoreGate/ReadHandleKey.h" +#include "StoreGate/WriteHandleKey.h" + +// Gaudi +#include "GaudiKernel/ServiceHandle.h" +#include "GaudiKernel/ToolHandle.h" + +// STL +#include <string> + +class PseudoCaloHitToWaveformRecAlg : public AthReentrantAlgorithm { + + public: + // Constructor + PseudoCaloHitToWaveformRecAlg(const std::string& name, ISvcLocator* pSvcLocator); + virtual ~PseudoCaloHitToWaveformRecAlg() = default; + + /** @name Usual algorithm methods */ + //@{ + virtual StatusCode initialize() override; + virtual StatusCode execute(const EventContext& ctx) const override; + virtual StatusCode finalize() override; + //@} + + private: + + /** @name Disallow default instantiation, copy, assignment */ + //@{ + PseudoCaloHitToWaveformRecAlg() = delete; + PseudoCaloHitToWaveformRecAlg(const PseudoCaloHitToWaveformRecAlg&) = delete; + PseudoCaloHitToWaveformRecAlg &operator=(const PseudoCaloHitToWaveformRecAlg&) = delete; + //@} + + /** + * @name Reconstruction tool + */ + ToolHandle<IPseudoSimToWaveformRecTool> m_recoTool + {this, "PseudoSimToWaveformRecTool", "PseudoSimToWaveformRecTool"}; + + /** + * @name Input simulated HITS using SG::ReadHandleKey + */ + //@{ + SG::ReadHandleKey<CaloHitCollection> m_caloHitContainerKey + {this, "CaloHitContainerKey", ""}; + //@} + + /** + * @name Output data using SG::WriteHandleKey + */ + //@{ + SG::WriteHandleKey<xAOD::WaveformHitContainer> m_waveformHitContainerKey + {this, "WaveformHitContainerKey", ""}; + //@} + +}; + +#endif // WAVERECALGS_PSEUDOCALOHITTOWAVEFORMRECALG_H diff --git a/Waveform/WaveRecAlgs/src/PseudoScintHitToWaveformRecAlg.cxx b/Waveform/WaveRecAlgs/src/PseudoScintHitToWaveformRecAlg.cxx new file mode 100644 index 0000000000000000000000000000000000000000..2b896e63c8a67b845072b74489ffe5a39b162751 --- /dev/null +++ b/Waveform/WaveRecAlgs/src/PseudoScintHitToWaveformRecAlg.cxx @@ -0,0 +1,68 @@ +#include "PseudoScintHitToWaveformRecAlg.h" + +PseudoScintHitToWaveformRecAlg::PseudoScintHitToWaveformRecAlg(const std::string& name, + ISvcLocator* pSvcLocator) + : AthReentrantAlgorithm(name, pSvcLocator) { + +} + +StatusCode +PseudoScintHitToWaveformRecAlg::initialize() { + ATH_MSG_INFO(name() << "::initalize()" ); + + // Initalize tools + ATH_CHECK( m_recoTool.retrieve() ); + + // Set key to read waveform from + ATH_CHECK( m_scintHitContainerKey.initialize() ); + + // Set key to write container + ATH_CHECK( m_waveformHitContainerKey.initialize() ); + + return StatusCode::SUCCESS; +} + +StatusCode +PseudoScintHitToWaveformRecAlg::finalize() { + ATH_MSG_INFO(name() << "::finalize()"); + + return StatusCode::SUCCESS; +} + +StatusCode +PseudoScintHitToWaveformRecAlg::execute(const EventContext& ctx) const { + ATH_MSG_DEBUG("Executing"); + + ATH_MSG_DEBUG("Run: " << ctx.eventID().run_number() + << " Event: " << ctx.eventID().event_number()); + + // Find the input HIT collection + SG::ReadHandle<ScintHitCollection> scintHitHandle(m_scintHitContainerKey, ctx); + + ATH_CHECK( scintHitHandle.isValid() ); + ATH_MSG_DEBUG("Found ReadHandle for ScintHitCollection " << m_scintHitContainerKey); + + if (scintHitHandle->size() == 0) { + ATH_MSG_DEBUG("ScintHitCollection found with zero length!"); + return StatusCode::SUCCESS; + } + + + // Find the output waveform container + SG::WriteHandle<xAOD::WaveformHitContainer> waveformHitContainerHandle(m_waveformHitContainerKey, ctx); + ATH_CHECK( waveformHitContainerHandle.record( std::make_unique<xAOD::WaveformHitContainer>(), + std::make_unique<xAOD::WaveformHitAuxContainer>() ) ); + + ATH_MSG_DEBUG("WaveformsHitContainer '" << waveformHitContainerHandle.name() << "' initialized"); + + + // "Reconstruct" the hits + CHECK( m_recoTool->reconstruct<ScintHitCollection>(scintHitHandle.ptr(), + waveformHitContainerHandle.ptr()) ); + + + ATH_MSG_DEBUG("WaveformsHitContainer '" << waveformHitContainerHandle.name() << "' filled with "<< waveformHitContainerHandle->size() <<" items"); + + return StatusCode::SUCCESS; +} + diff --git a/Waveform/WaveRecAlgs/src/PseudoScintHitToWaveformRecAlg.h b/Waveform/WaveRecAlgs/src/PseudoScintHitToWaveformRecAlg.h new file mode 100644 index 0000000000000000000000000000000000000000..741c291c92913fa7e16732e1891e6ce90a016c98 --- /dev/null +++ b/Waveform/WaveRecAlgs/src/PseudoScintHitToWaveformRecAlg.h @@ -0,0 +1,75 @@ +#ifndef WAVERECALGS_PSEUDOSCINTHITTOWAVEFORMRECALG_H +#define WAVERECALGS_PSEUDOSCINTHITTOWAVEFORMRECALG_H + +// Base class +#include "AthenaBaseComps/AthReentrantAlgorithm.h" + +// Data classes +#include "ScintSimEvent/ScintHitCollection.h" + +#include "xAODFaserWaveform/WaveformHit.h" +#include "xAODFaserWaveform/WaveformHitContainer.h" +#include "xAODFaserWaveform/WaveformHitAuxContainer.h" + +// Tool classes +#include "WaveRecTools/IPseudoSimToWaveformRecTool.h" + +// Handles +#include "StoreGate/ReadHandleKey.h" +#include "StoreGate/WriteHandleKey.h" + +// Gaudi +#include "GaudiKernel/ServiceHandle.h" +#include "GaudiKernel/ToolHandle.h" + +// STL +#include <string> + +class PseudoScintHitToWaveformRecAlg : public AthReentrantAlgorithm { + + public: + // Constructor + PseudoScintHitToWaveformRecAlg(const std::string& name, ISvcLocator* pSvcLocator); + virtual ~PseudoScintHitToWaveformRecAlg() = default; + + /** @name Usual algorithm methods */ + //@{ + virtual StatusCode initialize() override; + virtual StatusCode execute(const EventContext& ctx) const override; + virtual StatusCode finalize() override; + //@} + + private: + + /** @name Disallow default instantiation, copy, assignment */ + //@{ + PseudoScintHitToWaveformRecAlg() = delete; + PseudoScintHitToWaveformRecAlg(const PseudoScintHitToWaveformRecAlg&) = delete; + PseudoScintHitToWaveformRecAlg &operator=(const PseudoScintHitToWaveformRecAlg&) = delete; + //@} + + /** + * @name Reconstruction tool + */ + ToolHandle<IPseudoSimToWaveformRecTool> m_recoTool + {this, "PseudoSimToWaveformRecTool", "PseudoSimToWaveformRecTool"}; + + /** + * @name Input simulated HITS using SG::ReadHandleKey + */ + //@{ + SG::ReadHandleKey<ScintHitCollection> m_scintHitContainerKey + {this, "ScintHitContainerKey", ""}; + //@} + + /** + * @name Output data using SG::WriteHandleKey + */ + //@{ + SG::WriteHandleKey<xAOD::WaveformHitContainer> m_waveformHitContainerKey + {this, "WaveformHitContainerKey", ""}; + //@} + +}; + +#endif // WAVERECALGS_PSEUDOSCINTHITTOWAVEFORMRECALG_H diff --git a/Waveform/WaveRecAlgs/src/components/WaveRecAlgs_entries.cxx b/Waveform/WaveRecAlgs/src/components/WaveRecAlgs_entries.cxx index 21f97520a8e2da1093b5b4f929143cff32fd83a9..46d7f8980ad87690bc2a6eed5fd2de15a04a30dd 100644 --- a/Waveform/WaveRecAlgs/src/components/WaveRecAlgs_entries.cxx +++ b/Waveform/WaveRecAlgs/src/components/WaveRecAlgs_entries.cxx @@ -1,5 +1,9 @@ #include "../RawWaveformRecAlg.h" #include "../WaveClockRecAlg.h" +#include "../PseudoScintHitToWaveformRecAlg.h" +#include "../PseudoCaloHitToWaveformRecAlg.h" DECLARE_COMPONENT( RawWaveformRecAlg ) DECLARE_COMPONENT( WaveClockRecAlg ) +DECLARE_COMPONENT( PseudoScintHitToWaveformRecAlg ) +DECLARE_COMPONENT( PseudoCaloHitToWaveformRecAlg ) diff --git a/Waveform/WaveRecTools/WaveRecTools/IPseudoSimToWaveformRecTool.h b/Waveform/WaveRecTools/WaveRecTools/IPseudoSimToWaveformRecTool.h new file mode 100644 index 0000000000000000000000000000000000000000..6ba6e4eacdb8288a5f3667b02c3823aeaeb54d7e --- /dev/null +++ b/Waveform/WaveRecTools/WaveRecTools/IPseudoSimToWaveformRecTool.h @@ -0,0 +1,53 @@ +/* + Copyright (C) 2021 CERN for the benefit of the FASER collaboration +*/ + +/** + * @file IPseudoSimToWaveformRecTool.h + * Header file for the IPseudoSimToWaveformRecTool class + * @author Carl Gwilliam, 2021 + */ + + +#ifndef WAVERECTOOLS_IPSEUDOSIMTOWAVEFORMRECTOOL_H +#define WAVERECTOOLS_IPSEUDOSIMTOWAVEFORMRECTOOL_H + +// Base class +#include "GaudiKernel/IAlgTool.h" +#include "GaudiKernel/ToolHandle.h" + +#include "GaudiKernel/ServiceHandle.h" +#include "GaudiKernel/IMessageSvc.h" +#include "GaudiKernel/MsgStream.h" + +#include "xAODFaserWaveform/WaveformHitContainer.h" +#include "xAODFaserWaveform/WaveformHit.h" + +///Interface for Pseudo waveform rec tools +class IPseudoSimToWaveformRecTool : virtual public IAlgTool +{ +public: + + // InterfaceID + DeclareInterfaceID(IPseudoSimToWaveformRecTool, 1, 0); + + IPseudoSimToWaveformRecTool(): + m_msgSvc ( "MessageSvc", "ITrkEventCnvTool" ) + {} + + virtual ~IPseudoSimToWaveformRecTool() = default; + + // Convert sim hit energies to pseudo-waveform + template<class CONT> + StatusCode reconstruct(const CONT* hitCollection, + xAOD::WaveformHitContainer* waveContainer) const; + +private: + ServiceHandle<IMessageSvc> m_msgSvc; + +}; + +#include "WaveRecTools/IPseudoSimToWaveformRecTool.icc" + + +#endif // SCINTRECTOOLS_IPSEUDOSIMTOWAVEFORMRECTOOL_H diff --git a/Waveform/WaveRecTools/WaveRecTools/IPseudoSimToWaveformRecTool.icc b/Waveform/WaveRecTools/WaveRecTools/IPseudoSimToWaveformRecTool.icc new file mode 100644 index 0000000000000000000000000000000000000000..9f85c6b020b3efc9d81949cbf914b2a424bda19b --- /dev/null +++ b/Waveform/WaveRecTools/WaveRecTools/IPseudoSimToWaveformRecTool.icc @@ -0,0 +1,36 @@ +template<class CONT> +StatusCode IPseudoSimToWaveformRecTool::reconstruct(const CONT* hitCollection, + xAOD::WaveformHitContainer* container) const { + + + // Check the container + if (!container) { + MsgStream log(&(*m_msgSvc), name()); + log << MSG::ERROR << "HitCollection passed to reconstruct() is null!" << endmsg; + return StatusCode::FAILURE; + } + + std::map<int, float> idSums; + + // Sum hits in each "channel" + for (const auto& hit : *hitCollection) { + idSums[hit.identify()] += hit.energyLoss(); + } + + for (const auto& id : idSums) { + xAOD::WaveformHit* hit = new xAOD::WaveformHit(); + container->push_back(hit); + + hit->set_id(id.first); + hit->set_channel(0); + hit->set_peak(0); + hit->set_mean(0); + hit->set_width(0); + hit->set_integral(id.second); + hit->set_localtime(0); + hit->set_raw_peak(0); + hit->set_raw_integral(0); + } + + return StatusCode::SUCCESS; +} diff --git a/Waveform/WaveRecTools/src/PseudoSimToWaveformRecTool.cxx b/Waveform/WaveRecTools/src/PseudoSimToWaveformRecTool.cxx new file mode 100644 index 0000000000000000000000000000000000000000..6d894aeb64125274507105bd71ce37c7ad2bed31 --- /dev/null +++ b/Waveform/WaveRecTools/src/PseudoSimToWaveformRecTool.cxx @@ -0,0 +1,38 @@ +/* + Copyright (C) 2021 CERN for the benefit of the FASER collaboration +*/ + +/** + * @file PseudoSimToWaveformRec.cxx + * Implementation file for the PseudoSimToWaveformRec class + * @ author E. Torrence, 2021 + **/ + +#include "PseudoSimToWaveformRecTool.h" + +#include "xAODFaserWaveform/WaveformHit.h" + +#include "TH1F.h" +#include "TF1.h" +#include "TFitResult.h" +#include "TFitResultPtr.h" +#include "TGraph.h" + +#include <vector> +#include <tuple> +#include <math.h> + +// Constructor +PseudoSimToWaveformRecTool::PseudoSimToWaveformRecTool(const std::string& type, const std::string& name, const IInterface* parent) : + base_class(type, name, parent) +{ +} + +// Initialization +StatusCode +PseudoSimToWaveformRecTool::initialize() { + ATH_MSG_INFO( name() << "::initalize()" ); + return StatusCode::SUCCESS; +} + + diff --git a/Waveform/WaveRecTools/src/PseudoSimToWaveformRecTool.h b/Waveform/WaveRecTools/src/PseudoSimToWaveformRecTool.h new file mode 100644 index 0000000000000000000000000000000000000000..f3376f6ec4aa93becef3da8fc52c19d003d6dced --- /dev/null +++ b/Waveform/WaveRecTools/src/PseudoSimToWaveformRecTool.h @@ -0,0 +1,36 @@ +/* + Copyright (C) 2021 CERN for the benefit of the FASER collaboration +*/ + +/** @file PseduoSimToWaveformRecTool.h + * Header file for PseudoSimToWaveformRecTool.h + * + */ +#ifndef WAVERECTOOLS_PSEUDOSIMTOWAVEFORMRECTOOL_H +#define WAVERECTOOLS_PSEUDOSIMTOWAVEFORMRECTOOL_H + +//Athena +#include "AthenaBaseComps/AthAlgTool.h" +#include "WaveRecTools/IPseudoSimToWaveformRecTool.h" + +//Gaudi +#include "GaudiKernel/ToolHandle.h" + +//STL + +class PseudoSimToWaveformRecTool: public extends<AthAlgTool, IPseudoSimToWaveformRecTool> { + public: + + /// Normal constructor for an AlgTool; 'properties' are also declared here + PseudoSimToWaveformRecTool(const std::string& type, + const std::string& name, const IInterface* parent); + + /// Retrieve the necessary services in initialize + StatusCode initialize(); + + private: + // None + +}; + +#endif // WAVERECTOOLS_PSEUDOSIMTOWAVEFORMRECTOOL_H diff --git a/Waveform/WaveRecTools/src/WaveformReconstructionTool.cxx b/Waveform/WaveRecTools/src/WaveformReconstructionTool.cxx index 934f20803d18a1786ace8dedf18062203db8ba56..d34a654f505243f5333de060c8e7da942ab01093 100644 --- a/Waveform/WaveRecTools/src/WaveformReconstructionTool.cxx +++ b/Waveform/WaveRecTools/src/WaveformReconstructionTool.cxx @@ -59,7 +59,10 @@ WaveformReconstructionTool::reconstruct(const RawWaveform& raw_wave, // We always want to create at least one hit, so create it here xAOD::WaveformHit* hit = new xAOD::WaveformHit(); container->push_back(hit); + + // Set digitizer channel and identifier hit->set_channel(raw_wave.channel()); + hit->set_id(raw_wave.identify32().get_compact()); // // Make some sanity checks on the waveform data @@ -80,9 +83,6 @@ WaveformReconstructionTool::reconstruct(const RawWaveform& raw_wave, return StatusCode::SUCCESS; } - // Set channel - hit->set_channel(raw_wave.channel()); - // // Find baseline WaveformBaselineData baseline; diff --git a/Waveform/WaveRecTools/src/components/WaveRecTools_entries.cxx b/Waveform/WaveRecTools/src/components/WaveRecTools_entries.cxx index a96a42c9b5e4a0007558650314ef6e1144281877..d801b8a121d49472a72c875389d778ccfb175b60 100644 --- a/Waveform/WaveRecTools/src/components/WaveRecTools_entries.cxx +++ b/Waveform/WaveRecTools/src/components/WaveRecTools_entries.cxx @@ -1,5 +1,7 @@ #include "../ClockReconstructionTool.h" #include "../WaveformReconstructionTool.h" +#include "../PseudoSimToWaveformRecTool.h" DECLARE_COMPONENT( ClockReconstructionTool ) DECLARE_COMPONENT( WaveformReconstructionTool ) +DECLARE_COMPONENT( PseudoSimToWaveformRecTool ) diff --git a/xAOD/xAODFaserWaveform/CMakeLists.txt b/xAOD/xAODFaserWaveform/CMakeLists.txt index 1185251139a20ea833a672371173ad00c2987ccd..b5b6abf4398952c40739dab07f360785a9aecf0b 100644 --- a/xAOD/xAODFaserWaveform/CMakeLists.txt +++ b/xAOD/xAODFaserWaveform/CMakeLists.txt @@ -10,7 +10,7 @@ find_package( xAODUtilities ) atlas_add_library( xAODFaserWaveform xAODFaserWaveform/*.h xAODFaserWaveform/versions/*.h Root/*.cxx PUBLIC_HEADERS xAODFaserWaveform - LINK_LIBRARIES xAODCore ) + LINK_LIBRARIES Identifier xAODCore ) atlas_add_xaod_smart_pointer_dicts( INPUT xAODFaserWaveform/selection.xml @@ -21,6 +21,6 @@ atlas_add_xaod_smart_pointer_dicts( atlas_add_dictionary( xAODFaserWaveformDict xAODFaserWaveform/xAODFaserWaveformDict.h ${_selectionFile} - LINK_LIBRARIES xAODCore xAODFaserWaveform + LINK_LIBRARIES Identifier xAODCore xAODFaserWaveform EXTRA_FILES Root/dict/*.cxx ) diff --git a/xAOD/xAODFaserWaveform/Root/WaveformHitAuxContainer_v1.cxx b/xAOD/xAODFaserWaveform/Root/WaveformHitAuxContainer_v1.cxx index 76aed56247567f13e39fbcbf7471950dc93b3c23..2cd4ffcf6859a6af5341683af8333e1ba57be1ed 100644 --- a/xAOD/xAODFaserWaveform/Root/WaveformHitAuxContainer_v1.cxx +++ b/xAOD/xAODFaserWaveform/Root/WaveformHitAuxContainer_v1.cxx @@ -11,6 +11,7 @@ namespace xAOD { : AuxContainerBase() { AUX_VARIABLE(channel); + AUX_VARIABLE(id); AUX_VARIABLE(localtime); AUX_VARIABLE(peak); AUX_VARIABLE(width); diff --git a/xAOD/xAODFaserWaveform/Root/WaveformHit_v1.cxx b/xAOD/xAODFaserWaveform/Root/WaveformHit_v1.cxx index 4d4ce278978e68024637c117bc1e9456c2e22730..8b26eebffaa1018193523608e0494ef5fcc29d44 100644 --- a/xAOD/xAODFaserWaveform/Root/WaveformHit_v1.cxx +++ b/xAOD/xAODFaserWaveform/Root/WaveformHit_v1.cxx @@ -15,6 +15,8 @@ namespace xAOD { AUXSTORE_PRIMITIVE_SETTER_AND_GETTER( WaveformHit_v1, unsigned int, channel, set_channel ) + AUXSTORE_PRIMITIVE_SETTER_AND_GETTER( WaveformHit_v1, unsigned int, id, set_id ) + AUXSTORE_PRIMITIVE_SETTER_AND_GETTER( WaveformHit_v1, float, localtime, set_localtime ) AUXSTORE_PRIMITIVE_SETTER_AND_GETTER( WaveformHit_v1, float, peak, set_peak ) diff --git a/xAOD/xAODFaserWaveform/xAODFaserWaveform/versions/WaveformHitAuxContainer_v1.h b/xAOD/xAODFaserWaveform/xAODFaserWaveform/versions/WaveformHitAuxContainer_v1.h index 0ea61b39a567a61882a390aa35277363771e970f..c272fcd2c8d8ce8323b819ea8ca02ceeb02b45f9 100644 --- a/xAOD/xAODFaserWaveform/xAODFaserWaveform/versions/WaveformHitAuxContainer_v1.h +++ b/xAOD/xAODFaserWaveform/xAODFaserWaveform/versions/WaveformHitAuxContainer_v1.h @@ -29,6 +29,7 @@ namespace xAOD { /// @name Basic variables ///@ { std::vector<unsigned int> channel; + std::vector<unsigned int> id; std::vector<float> localtime; std::vector<float> peak; std::vector<float> width; diff --git a/xAOD/xAODFaserWaveform/xAODFaserWaveform/versions/WaveformHit_v1.h b/xAOD/xAODFaserWaveform/xAODFaserWaveform/versions/WaveformHit_v1.h index 8b69c776da83546cde3a787a25bcf2b3a830cf7f..9ea7b6dc573d53a533df44d02ff3edbf92c8a9ce 100644 --- a/xAOD/xAODFaserWaveform/xAODFaserWaveform/versions/WaveformHit_v1.h +++ b/xAOD/xAODFaserWaveform/xAODFaserWaveform/versions/WaveformHit_v1.h @@ -16,6 +16,7 @@ extern "C" { // Core EDM include(s): #include "AthContainers/AuxElement.h" +#include "Identifier/Identifier.h" namespace xAOD { @@ -42,10 +43,18 @@ namespace xAOD { /// @name Access WaveformHit elements /// @{ - /// Waveform channel + /// Waveform digitizer channel unsigned int channel() const; void set_channel(unsigned int value); + /// 32-bit version of channel identifier + unsigned int id() const; + void set_id(unsigned int value); + + Identifier identify() const { + return Identifier(this->id()); + } + /// Best results float localtime() const; void set_localtime(float value);