diff --git a/Calorimeter/CaloDigiAlgs/CMakeLists.txt b/Calorimeter/CaloDigiAlgs/CMakeLists.txt index 5a4de58a7469a2acb320de9c758ea1aef98ebcac..89d20eae8e67a52daa6a23b60825e06f43610ec3 100644 --- a/Calorimeter/CaloDigiAlgs/CMakeLists.txt +++ b/Calorimeter/CaloDigiAlgs/CMakeLists.txt @@ -9,7 +9,9 @@ atlas_subdir( CaloDigiAlgs ) atlas_add_component( CaloDigiAlgs src/*.cxx src/*.h src/components/*.cxx - LINK_LIBRARIES AthenaBaseComps Identifier FaserCaloIdentifier StoreGateLib WaveRawEvent FaserCaloSimEvent WaveDigiToolsLib) + LINK_LIBRARIES AthenaBaseComps Identifier FaserCaloIdentifier + WaveformConditionsToolsLib StoreGateLib WaveRawEvent + FaserCaloSimEvent WaveDigiToolsLib) atlas_install_python_modules( python/*.py ) diff --git a/Calorimeter/CaloDigiAlgs/python/CaloDigiAlgsConfig.py b/Calorimeter/CaloDigiAlgs/python/CaloDigiAlgsConfig.py index ce45e6257885956f5edb6284780ff9756b51b9fc..793f73401b54c91724bdb6283fa4c6099811628c 100644 --- a/Calorimeter/CaloDigiAlgs/python/CaloDigiAlgsConfig.py +++ b/Calorimeter/CaloDigiAlgs/python/CaloDigiAlgsConfig.py @@ -6,6 +6,7 @@ from AthenaConfiguration.ComponentFactory import CompFactory from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg +from WaveformConditionsTools.WaveformCableMappingConfig import WaveformCableMappingCfg # One stop shopping for normal FASER data def CaloWaveformDigitizationCfg(flags): @@ -17,6 +18,8 @@ def CaloWaveformDigitizationCfg(flags): acc.merge(CaloWaveformDigiCfg(flags, "CaloWaveformDigiAlg")) acc.merge(CaloWaveformDigitizationOutputCfg(flags)) + acc.merge(WaveformCableMappingCfg(flags)) + return acc # Return configured digitization algorithm from SIM hits @@ -54,6 +57,8 @@ def CaloWaveformDigitizationOutputCfg(flags, **kwargs): ] acc.merge(OutputStreamCfg(flags, "RDO")) ostream = acc.getEventAlgo("OutputStreamRDO") - ostream.TakeItemsFromInput = True # Copies all data from input file to output + # ostream.TakeItemsFromInput = True # Copies all data from input file to output + # ostream.TakeItemsFromInput = False + # Try turning this off ostream.ItemList += ItemList return acc diff --git a/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx b/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx index c54e1e37344d963f28a12f9c9f8b5d86781f93c4..3e10288eaa2e68f82f1c72be38fde04cb73c66fd 100644 --- a/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx +++ b/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx @@ -20,6 +20,7 @@ CaloWaveformDigiAlg::initialize() { // Initalize tools ATH_CHECK( m_digiTool.retrieve() ); + ATH_CHECK( m_mappingTool.retrieve() ); // Set key to read waveform from ATH_CHECK( m_caloHitContainerKey.initialize() ); @@ -104,6 +105,22 @@ CaloWaveformDigiAlg::execute(const EventContext& ctx) const { } + // This is a bit of a hack to make sure all waveforms have + // at least baseline entries. Should really add this to the + // logic above + for (const auto& w : waveforms) { + if (w.second.size() > 0) continue; + + // Waveform was empty, fill with baseline + int channel = m_mappingTool->getChannelMapping(w.first); + ATH_MSG_DEBUG("Writing baseline into empty waveform in channel "<< channel); + int i = m_digiTool->nsamples(); + while(i--) { // Use while to avoid unused variable warning with for + int baseline = m_digiTool->generate_baseline(m_base_mean, m_base_rms); + waveforms[w.first].push_back(baseline); + } + } + //m_chrono->chronoStop("Digit"); //m_chrono->chronoStart("Write"); @@ -113,6 +130,7 @@ CaloWaveformDigiAlg::execute(const EventContext& ctx) const { RawWaveform* wfm = new RawWaveform(); wfm->setWaveform(0, w.second); wfm->setIdentifier(w.first); + wfm->setChannel(m_mappingTool->getChannelMapping(w.first)); wfm->setSamples(nsamples); waveformContainerHandle->push_back(wfm); } diff --git a/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.h b/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.h index 8eeb0a44eba0dde0ede044036b57406cfcf25efc..25111bc9655ae99d97e7ea56a98524102ce151f3 100644 --- a/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.h +++ b/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.h @@ -10,6 +10,7 @@ // Tool classes #include "WaveDigiTools/IWaveformDigitisationTool.h" +#include "WaveformConditionsTools/IWaveformCableMappingTool.h" // Handles #include "StoreGate/ReadHandleKey.h" @@ -79,6 +80,11 @@ class CaloWaveformDigiAlg : public AthReentrantAlgorithm { ToolHandle<IWaveformDigitisationTool> m_digiTool {this, "WaveformDigitisationTool", "WaveformDigitisationTool"}; + /** + * @name Mapping tool + */ + ToolHandle<IWaveformCableMappingTool> m_mappingTool + {this, "WaveformCableMappingTool", "WaveformCableMappingTool"}; /** * @name Input HITS using SG::ReadHandleKey diff --git a/Calorimeter/FaserCaloSimEvent/FaserCaloSimEvent/CaloHit.h b/Calorimeter/FaserCaloSimEvent/FaserCaloSimEvent/CaloHit.h index b1e955f845d8f82ea8dcff82a6bbd5883f96b4b1..0b6fbe42ad0dd24cb4adbcb8072a65055d6ad2be 100644 --- a/Calorimeter/FaserCaloSimEvent/FaserCaloSimEvent/CaloHit.h +++ b/Calorimeter/FaserCaloSimEvent/FaserCaloSimEvent/CaloHit.h @@ -16,6 +16,8 @@ #include "CLHEP/Geometry/Point3D.h" #include "GeneratorObjects/HepMcParticleLink.h" +class Identifier; + class CaloHit { /////////////////////////////////////////////////////////////////// @@ -80,8 +82,13 @@ public: // Const methods: /////////////////////////////////////////////////////////////////// + // This returns the HitId, used in Geant. This is not the detector Identifier unsigned int identify() const; + // This is the detector readout Identifier (pmt Identifier) + // provided by ScintHitIdHelper::getIdentifier + Identifier getIdentifier() const; + // local start position of the energy deposit: HepGeom::Point3D<double> localStartPosition() const; diff --git a/Calorimeter/FaserCaloSimEvent/src/CaloHit.cxx b/Calorimeter/FaserCaloSimEvent/src/CaloHit.cxx index c8e18d7fc737268e388bf469be037ea43655a6d3..86e895117d03eb537127c68d74efef7b08bf8788 100644 --- a/Calorimeter/FaserCaloSimEvent/src/CaloHit.cxx +++ b/Calorimeter/FaserCaloSimEvent/src/CaloHit.cxx @@ -142,6 +142,10 @@ int CaloHit::getModule() const { return CaloHitIdHelper::GetHelper()->getModule(m_ID); } +Identifier CaloHit::getIdentifier() const { + return CaloHitIdHelper::GetHelper()->getIdentifier(m_ID); +} + void CaloHit::print() const { std::cout << "*** Calorimeter Hit" << std::endl; std::cout << " Station Number " << getRow() << std::endl; diff --git a/Control/CalypsoExample/Digitization/scripts/faserMDC_digi.py b/Control/CalypsoExample/Digitization/scripts/faserMDC_digi.py new file mode 100755 index 0000000000000000000000000000000000000000..be31e19994c9e93391f12228f80ca1096c222d51 --- /dev/null +++ b/Control/CalypsoExample/Digitization/scripts/faserMDC_digi.py @@ -0,0 +1,144 @@ +#!/usr/bin/env python +# +# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +# Run with: +# ./faser_digi.py filepath runtype +# +# filepath - fully qualified path, including url if needed, to the input HITS file +# example: "root://eospublic.cern.ch//eos/experiment/faser/sim/GeniePilot/HITS/1/faser.150fbInv.1.001.HITS.pool.root" +# +# For MDC, we assume this is TI12 geometry +# +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("-t", "--tag", default="", + help="Specify digi tag (to append to output filename)") +parser.add_argument("-n", "--nevents", type=int, default=-1, + help="Specify number of events to process (default: all)") +parser.add_argument("-v", "--verbose", action='store_true', + help="Turn on DEBUG output") + +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 + +else: + runtype = "TI12MC" + print(f"Assuming {runtype} geometry for MDC") + +print(f"Starting digitization of {filepath.name} with type {runtype}") +if args.nevents > 0: + print(f"Reconstructing {args.nevents} events by command-line option") + +# Start digitization + +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 = True # Needed to bypass autoconfig +ConfigFlags.IOVDb.DatabaseInstance = "OFLP200" # Use MC conditions for now + +ConfigFlags.Input.ProjectName = "mc20" +ConfigFlags.GeoModel.Align.Dynamic = False +ConfigFlags.Beam.NumberOfCollisions = 0. +ConfigFlags.Digitization.TruthOutput = True + +# TI12 old geometry +if runtype == "TI12OldMC": + ConfigFlags.GeoModel.FaserVersion = "FASER-01" + ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-01" + +# Testbeam setup +elif runtype == "TestBeamMC" : + ConfigFlags.GeoModel.FaserVersion = "FASER-TB00" + ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-TB00" + +# New TI12 geometry (ugh) +elif runtype == "TI12MC": + ConfigFlags.GeoModel.FaserVersion = "FASERNU-03" + 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 +# Remove any filetype modifier +if filestem[-5:] == "-HITS": + filestem = filestem[:-5] + +if len(args.tag) > 0: + filestem += f"-{args.tag}" + +ConfigFlags.Output.RDOFileName = f"{filestem}-RDO.root" + +# +# Play around with this? +# ConfigFlags.Concurrency.NumThreads = 2 +# ConfigFlags.Concurrency.NumConcurrentEvents = 2 +ConfigFlags.lock() + +# +# Configure components +from CalypsoConfiguration.MainServicesConfig import MainServicesCfg +from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg +from AthenaPoolCnvSvc.PoolWriteConfig import PoolWriteCfg + +acc = MainServicesCfg(ConfigFlags) +acc.merge(PoolReadCfg(ConfigFlags)) +acc.merge(PoolWriteCfg(ConfigFlags)) + +# +# Needed, or move to MainServicesCfg? +from FaserGeoModel.FaserGeoModelConfig import FaserGeometryCfg +acc.merge(FaserGeometryCfg(ConfigFlags)) + +# Set up algorithms +from FaserSCT_Digitization.FaserSCT_DigitizationConfigNew import FaserSCT_DigitizationCfg +acc.merge(FaserSCT_DigitizationCfg(ConfigFlags)) + +from CaloDigiAlgs.CaloDigiAlgsConfig import CaloWaveformDigitizationCfg +acc.merge(CaloWaveformDigitizationCfg(ConfigFlags)) + +from ScintDigiAlgs.ScintDigiAlgsConfig import ScintWaveformDigitizationCfg +acc.merge(ScintWaveformDigitizationCfg(ConfigFlags)) + +# Configure verbosity +# ConfigFlags.dump() +if args.verbose: + acc.foreach_component("*").OutputLevel = VERBOSE + +else: + acc.foreach_component("*").OutputLevel = INFO + +acc.foreach_component("*ClassID*").OutputLevel = INFO + +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/Digitization/scripts/submit_faserMDC_digi.sh b/Control/CalypsoExample/Digitization/scripts/submit_faserMDC_digi.sh new file mode 100755 index 0000000000000000000000000000000000000000..2070b43f3979c92823d2cdd78e0077cb7343c4d9 --- /dev/null +++ b/Control/CalypsoExample/Digitization/scripts/submit_faserMDC_digi.sh @@ -0,0 +1,128 @@ +#!/bin/bash +# Used with a condor file to submit to vanilla universe +# +# Usage: +# submit_faserMDC_digi.sh filepath [release_directory] [working_directory] +# +# filepath - 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) +# +# 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 gen/g???? it will be passed to the job +# +#---------------------------------------- +# +# Parse command-line options +file_path=${1} +release_directory=${2} +working_directory=${3} +# +# Set defaults if arguments aren't provided +if [ -z "$file_path" ] +then + echo "No file specified!" + echo "Usage: submit_faserMDC_digi.sh file [release dir] [output dir]" + exit 1 +fi +# +if [ -z "$release_directory" ] +then + release_directory=`pwd` +fi +# +if [ -z "$working_directory" ] +then + working_directory=`pwd` +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 short 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_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=calypso/asetup.faser Athena,22.0.49 +source build/x86*/setup.sh +# +# +# Try to find a release tag +cd calypso +recotag=`git describe` +if [[ "$recotag" == "reco/r"???? ]]; then + tag=`echo "$recotag" | cut -c 6-11` + echo "Found reco tag: $tag" +fi +if [[ "$recotag" == "digi/d"???? ]]; then + tag=`echo "$recotag" | cut -c 6-11` + echo "Found reco tag: $tag" +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 "$tag" ]]; then + faserMDC_digi.py "$file_path" +else + faserMDC_digi.py "--tag=$tag" "$file_path" +fi +# +# Print out ending time +date diff --git a/Control/CalypsoExample/Generation/CMakeLists.txt b/Control/CalypsoExample/Generation/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..4ffb523fe90a093c43dc84e121ae7792078e8047 --- /dev/null +++ b/Control/CalypsoExample/Generation/CMakeLists.txt @@ -0,0 +1,12 @@ +################################################################################ +# Package: Generation +################################################################################ + +# Declare the package name: +atlas_subdir( Generation ) + +# Install files from the package: +atlas_install_python_modules( python/*.py ) +atlas_install_scripts( scripts/*.sh scripts/*.py ) + + diff --git a/Control/CalypsoExample/Generation/data/mdc/FaserMC-MDC_PG_emi_100GeV-101104.json b/Control/CalypsoExample/Generation/data/mdc/FaserMC-MDC_PG_emi_100GeV-101104.json new file mode 100644 index 0000000000000000000000000000000000000000..434153301a8580ef241b7186c448003f67016f01 --- /dev/null +++ b/Control/CalypsoExample/Generation/data/mdc/FaserMC-MDC_PG_emi_100GeV-101104.json @@ -0,0 +1,13 @@ +{ + "file_length": 500, + "mass": 0.511, + "maxE": 100.0, + "minE": 100.0, + "pid": 11, + "radius": 100.0, + "run": 101104, + "sampler": "const", + "segment": 0, + "short": "MDC_PG_emi_100GeV", + "zpos": -1000.0 +} diff --git a/Control/CalypsoExample/Generation/data/mdc/FaserMC-MDC_PG_emi_logE-101102.json b/Control/CalypsoExample/Generation/data/mdc/FaserMC-MDC_PG_emi_logE-101102.json new file mode 100644 index 0000000000000000000000000000000000000000..56410bd94a49de9dab87b529b3eb5db1f2976591 --- /dev/null +++ b/Control/CalypsoExample/Generation/data/mdc/FaserMC-MDC_PG_emi_logE-101102.json @@ -0,0 +1,13 @@ +{ + "file_length": 200, + "mass": 0.511, + "maxE": 1000.0, + "minE": 10.0, + "pid": 11, + "radius": 100.0, + "run": 101102, + "sampler": "log", + "segment": 0, + "short": "MDC_PG_emi_logE", + "zpos": -1000.0 +} diff --git a/Control/CalypsoExample/Generation/data/mdc/FaserMC-MDC_PG_epl_100GeV-101103.json b/Control/CalypsoExample/Generation/data/mdc/FaserMC-MDC_PG_epl_100GeV-101103.json new file mode 100644 index 0000000000000000000000000000000000000000..4b540f260b55a33fe744c32bea39b505f9fbc397 --- /dev/null +++ b/Control/CalypsoExample/Generation/data/mdc/FaserMC-MDC_PG_epl_100GeV-101103.json @@ -0,0 +1,13 @@ +{ + "file_length": 500, + "mass": 0.511, + "maxE": 100.0, + "minE": 100.0, + "pid": -11, + "radius": 100.0, + "run": 101103, + "sampler": "const", + "segment": 0, + "short": "MDC_PG_epl_100GeV", + "zpos": -1000.0 +} diff --git a/Control/CalypsoExample/Generation/data/mdc/FaserMC-MDC_PG_epl_logE-101101.json b/Control/CalypsoExample/Generation/data/mdc/FaserMC-MDC_PG_epl_logE-101101.json new file mode 100644 index 0000000000000000000000000000000000000000..95fd9fa706c2cd385859df1d00beefaa9e802571 --- /dev/null +++ b/Control/CalypsoExample/Generation/data/mdc/FaserMC-MDC_PG_epl_logE-101101.json @@ -0,0 +1,13 @@ +{ + "file_length": 200, + "mass": 0.511, + "maxE": 1000.0, + "minE": 10.0, + "pid": -11, + "radius": 100.0, + "run": 101101, + "sampler": "log", + "segment": 0, + "short": "MDC_PG_epl_logE", + "zpos": -1000.0 +} diff --git a/Control/CalypsoExample/Generation/data/mdc/FaserMC-MDC_PG_gam_100GeV-102201.json b/Control/CalypsoExample/Generation/data/mdc/FaserMC-MDC_PG_gam_100GeV-102201.json new file mode 100644 index 0000000000000000000000000000000000000000..04ddbb05ec2e43348693ea0a16203c921f0ac2de --- /dev/null +++ b/Control/CalypsoExample/Generation/data/mdc/FaserMC-MDC_PG_gam_100GeV-102201.json @@ -0,0 +1,13 @@ +{ + "file_length": 500, + "mass": 0.0, + "maxE": 100.0, + "minE": 100.0, + "pid": 22, + "radius": 100.0, + "run": 102201, + "sampler": "const", + "segment": 0, + "short": "MDC_PG_gam_100GeV", + "zpos": -1000.0 +} diff --git a/Control/CalypsoExample/Generation/data/mdc/FaserMC-MDC_PG_mumi_logE-101302.json b/Control/CalypsoExample/Generation/data/mdc/FaserMC-MDC_PG_mumi_logE-101302.json new file mode 100644 index 0000000000000000000000000000000000000000..1387658b5efee880662a42560523f1ebff81e68f --- /dev/null +++ b/Control/CalypsoExample/Generation/data/mdc/FaserMC-MDC_PG_mumi_logE-101302.json @@ -0,0 +1,13 @@ +{ + "file_length": 2000, + "mass": 105.66, + "maxE": 2000.0, + "minE": 10.0, + "pid": 13, + "radius": 100.0, + "run": 101302, + "sampler": "log", + "segment": 0, + "short": "MDC_PG_mumi_logE", + "zpos": null +} diff --git a/Control/CalypsoExample/Generation/data/mdc/FaserMC-MDC_PG_mupl_logE-101301.json b/Control/CalypsoExample/Generation/data/mdc/FaserMC-MDC_PG_mupl_logE-101301.json new file mode 100644 index 0000000000000000000000000000000000000000..dda3d48ab90dbb3e963ac81996e26ad0864b2a7b --- /dev/null +++ b/Control/CalypsoExample/Generation/data/mdc/FaserMC-MDC_PG_mupl_logE-101301.json @@ -0,0 +1,13 @@ +{ + "file_length": 2000, + "mass": 105.66, + "maxE": 2000.0, + "minE": 10.0, + "pid": -13, + "radius": 100.0, + "run": 101301, + "sampler": "log", + "segment": 0, + "short": "MDC_PG_mupl_logE", + "zpos": null +} diff --git a/Control/CalypsoExample/Generation/data/mdc/FaserMC-MDC_PG_pimi_100GeV-121102.json b/Control/CalypsoExample/Generation/data/mdc/FaserMC-MDC_PG_pimi_100GeV-121102.json new file mode 100644 index 0000000000000000000000000000000000000000..cd8abc7eba3cbfbf56570fe320b8222b829bb190 --- /dev/null +++ b/Control/CalypsoExample/Generation/data/mdc/FaserMC-MDC_PG_pimi_100GeV-121102.json @@ -0,0 +1,13 @@ +{ + "file_length": 500, + "mass": 139.6, + "maxE": 100.0, + "minE": 100.0, + "pid": -211, + "radius": 100.0, + "run": 121102, + "sampler": "const", + "segment": 0, + "short": "MDC_PG_pimi_100GeV", + "zpos": -1000.0 +} diff --git a/Control/CalypsoExample/Generation/data/mdc/FaserMC-MDC_PG_pipl_100GeV-121101.json b/Control/CalypsoExample/Generation/data/mdc/FaserMC-MDC_PG_pipl_100GeV-121101.json new file mode 100644 index 0000000000000000000000000000000000000000..9874ba205e3142a184bffd54d793c407b621b03f --- /dev/null +++ b/Control/CalypsoExample/Generation/data/mdc/FaserMC-MDC_PG_pipl_100GeV-121101.json @@ -0,0 +1,13 @@ +{ + "file_length": 500, + "mass": 139.6, + "maxE": 100.0, + "minE": 100.0, + "pid": 211, + "radius": 100.0, + "run": 121101, + "sampler": "const", + "segment": 0, + "short": "MDC_PG_pipl_100GeV", + "zpos": -1000.0 +} diff --git a/Control/CalypsoExample/Generation/python/__init__.py b/Control/CalypsoExample/Generation/python/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..f078213647f851a7dfbd2d5ca63157209ff91b22 --- /dev/null +++ b/Control/CalypsoExample/Generation/python/__init__.py @@ -0,0 +1,2 @@ +# Generation + diff --git a/Control/CalypsoExample/Generation/python/faserMDC_pgparser.py b/Control/CalypsoExample/Generation/python/faserMDC_pgparser.py new file mode 100644 index 0000000000000000000000000000000000000000..7c41638ea2c603bebafea3ba3a72d62ea2cb2eed --- /dev/null +++ b/Control/CalypsoExample/Generation/python/faserMDC_pgparser.py @@ -0,0 +1,92 @@ +# +# Copyright (C) 2022 CERN for the benefit of the ATLAS collaboration +# Copyright (C) 2022 CERN for the benefit of the FASER collaboration +# +# Parser function for MDC particle gun samples +# +def faserMDC_pgparser(): + + import sys + import json + import argparse + + parser = argparse.ArgumentParser(description="Run FASER ParticleGun Simulation") + + parser.add_argument("--conf", action='append', + help="Specify configuration file with default values") + parser.add_argument("--run", default=123456, type=int, + help="Run number to generate") + parser.add_argument("--segment", default=00000, type=int, + help="Segment number to generate") + parser.add_argument("--file_length", default=1000, type=int, + help="Total events per file segement") + + parser.add_argument("--short", default="MDCPG_logE", + help="Short description for filename") + parser.add_argument("--tag", default=None, + help="Generator tag (g0000)") + + parser.add_argument("--pid", default=13, type=int, + help="Specify PDG ID of particle (note plus/minus different)") + parser.add_argument("--mass", default=105.66, type=float, + help="Specify particle mass (in MeV)") + parser.add_argument("--radius", default=100., type=float, + help="Specify radius (in mm)") + parser.add_argument("--angle", default=0.015, type=float, + help="Specify angular width (in Rad)") + parser.add_argument("--zpos", default=None, type=float, + help="Specify z position of particles (in mm) (helpful to avoid FASERnu)") + + parser.add_argument("--sampler", default="log", + help="Specify energy sampling (log, lin, const)") + parser.add_argument("--minE", default=10., type=float, + help="Minimum energy in GeV (for log or lin sampler)") + parser.add_argument("--maxE", default=1000., type=float, + help="Maximum energy (or constant) in GeV") + + parser.add_argument("--nevts", default=-1, type=int, + help="Number of events to generate (for debugging)") + parser.add_argument("--dump", action='store_true', + help="Write out full configuration") + parser.add_argument("--noexec", action='store_true', + help="Exit after parsing configuration (no execution)") + + pg_args = parser.parse_args() + + # Get defaults + if pg_args.conf is not None: + for conf_fname in pg_args.conf: + with open(conf_fname, 'r') as f: + parser.set_defaults(**json.load(f)) + + # Reload arguments to override config file with command line + pg_args = parser.parse_args() + + # Print out configuration if requested + if pg_args.dump: + tmp_args = vars(pg_args).copy() + del tmp_args['dump'] # Don't dump the dump value + del tmp_args['conf'] # Don't dump the conf file name either + del tmp_args['nevts'] # Debugging, not part of configuration + del tmp_args['noexec'] # Debugging, not part of configuration + print("Configuration:") + print(json.dumps(tmp_args, indent=4, sort_keys=False)) + + if pg_args.noexec: + sys.exit(0) + + # + # Add some derived quantities + # + + # Create the file name also (could add gentag here) + pg_args.outfile = f"FaserMC-{pg_args.short}-{pg_args.run:06}-{pg_args.segment:05}" + + if pg_args.tag: + pg_args.outfile += f"-{pg_args.tag}" + + pg_args.outfile += "-HITS.root" + + return pg_args + +# All done diff --git a/Control/CalypsoExample/Generation/scripts/faserMDC_particlegun.py b/Control/CalypsoExample/Generation/scripts/faserMDC_particlegun.py new file mode 100755 index 0000000000000000000000000000000000000000..fe827c0fdea36cd271a579c29a2f7b47776a7bf5 --- /dev/null +++ b/Control/CalypsoExample/Generation/scripts/faserMDC_particlegun.py @@ -0,0 +1,198 @@ +#!/usr/bin/env python +""" +Produce particle gun samples +Derived from G4FaserAlgConfigNew + +Usage: +faserMDC_particlegun.py --conf=<config_file> + +Copyright (C) 2002-2021 CERN for the benefit of the ATLAS and FASER collaborations +""" + +if __name__ == '__main__': + + import sys + import time + a = time.time() +# +# Parse command-line options +# + from Generation.faserMDC_pgparser import faserMDC_pgparser + args = faserMDC_pgparser() +# +# Figure out events to run and skip +# + nskipped = args.segment*args.file_length + if args.nevts > 0: + nevents = args.nevts + else: + nevents = args.file_length +# +# Print out what we are doing +# + print(f"Generating {nevents} in file {args.outfile}") +# +# Set up logging and config behaviour +# + from AthenaCommon.Logging import log + from AthenaCommon.Constants import DEBUG, VERBOSE + from AthenaCommon.Configurable import Configurable + log.setLevel(DEBUG) + Configurable.configurableRun3Behavior = 1 +# +# Import and set config flags +# + from CalypsoConfiguration.AllConfigFlags import ConfigFlags + ConfigFlags.Exec.MaxEvents = nevents + ConfigFlags.Exec.SkipEvents = nskipped + from AthenaConfiguration.Enums import ProductionStep + ConfigFlags.Common.ProductionStep = ProductionStep.Simulation +# +# All these must be specified to avoid auto-configuration +# + ConfigFlags.Input.RunNumber = [args.run] #Isn't updating - todo: investigate + ConfigFlags.Input.OverrideRunNumber = True + ConfigFlags.Input.LumiBlockNumber = [(args.segment+1)] + ConfigFlags.Input.isMC = True +# +# Output file name +# + ConfigFlags.Output.HITSFileName = args.outfile +# +# Sim ConfigFlags +# + ConfigFlags.Sim.Layout = "FASER" + ConfigFlags.Sim.PhysicsList = "FTFP_BERT" + ConfigFlags.Sim.ReleaseGeoModel = False + ConfigFlags.Sim.IncludeParentsInG4Event = True # Controls whether BeamTruthEvent is written to output HITS file + ConfigFlags.addFlag("Sim.Gun",{"Generator" : "SingleParticle"}) # Property bag for particle gun keyword:argument pairs + ConfigFlags.addFlag("Sim.Beam.xangle", 0) # Potential beam crossing angles + ConfigFlags.addFlag("Sim.Beam.yangle", 0) + ConfigFlags.addFlag("Sim.Beam.xshift", 0) # Potential beam shift + ConfigFlags.addFlag("Sim.Beam.yshift", 0) + + ConfigFlags.GeoModel.FaserVersion = "FASERNU-03" # Geometry set-up + ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-02" # Conditions set-up + ConfigFlags.addFlag("Input.InitialTimeStamp", 0) # To avoid autoconfig + ConfigFlags.GeoModel.Align.Dynamic = False + +# +# Preset particle gun parameters +# + from math import atan + from AthenaCommon.SystemOfUnits import GeV, TeV, cm, m + from AthenaCommon.PhysicalConstants import pi + + import ParticleGun as PG + ConfigFlags.Sim.Gun = { + "Generator" : "SingleParticle", + "pid" : args.pid, "mass" : args.mass, + "theta" : PG.GaussianSampler(0, args.angle, oneside = True), + "phi" : [0, 2*pi], "radius" : args.radius, + "randomSeed" : args.outfile } + + if args.zpos: + ConfigFlags.Sim.Gun["z"] = args.zpos + + # Determine energy sampling + if args.sampler == "lin": + ConfigFlags.Sim.Gun["energy"] = PG.UniformSampler(args.minE*GeV, args.maxE*GeV) + elif args.sampler == "log": + ConfigFlags.Sim.Gun["energy"] = PG.LogSampler(args.minE*GeV, args.maxE*GeV) + elif args.sampler == "const": + ConfigFlags.Sim.Gun["energy"] = PG.ConstSampler(args.maxE*GeV) + else: + print(f"Sampler {args.sampler} not known!") + sys.exit(1) + + # import sys + # ConfigFlags.fillFromArgs(sys.argv[1:]) + +# +# MDC geometry configuration +# + detectors = ['Veto', 'VetoNu', 'Preshower', 'FaserSCT', 'Ecal', 'Trigger', 'Dipole', 'Emulsion'] +# +# Setup detector flags +# + from CalypsoConfiguration.DetectorConfigFlags import setupDetectorsFromList + setupDetectorsFromList(ConfigFlags, detectors, toggle_geometry=True) +# +# Finalize flags +# + ConfigFlags.lock() +# +# Initialize a new component accumulator +# + from CalypsoConfiguration.MainServicesConfig import MainServicesCfg + cfg = MainServicesCfg(ConfigFlags) +# +# Configure the particle gun as requested, or using defaults +# + +# +# Particle gun generators - the generator, energy, angle, particle type, position, etc can be modified by passing keyword arguments +# + from FaserParticleGun.FaserParticleGunConfig import FaserParticleGunCfg + cfg.merge(FaserParticleGunCfg(ConfigFlags)) + from McEventSelector.McEventSelectorConfig import McEventSelectorCfg + cfg.merge(McEventSelectorCfg(ConfigFlags)) + +# +# Output file +# + from AthenaPoolCnvSvc.PoolWriteConfig import PoolWriteCfg + cfg.merge(PoolWriteCfg(ConfigFlags)) + +# +# Shift LOS +# + + if (ConfigFlags.Sim.Beam.xangle or ConfigFlags.Sim.Beam.yangle or + ConfigFlags.Sim.Beam.xshift or ConfigFlags.Sim.Beam.yshift): + + MCEventKey = "BeamTruthEventShifted" + import McParticleEvent.Pythonizations + from GeneratorUtils.ShiftLOSConfig import ShiftLOSCfg + cfg.merge(ShiftLOSCfg(ConfigFlags, OutputMCEventKey = MCEventKey, + xcross = ConfigFlags.Sim.Beam.xangle, + ycross = ConfigFlags.Sim.Beam.yangle, + xshift = ConfigFlags.Sim.Beam.xshift, + yshift = ConfigFlags.Sim.Beam.yshift)) + + else: + MCEventKey = "BeamTruthEvent" + +# +# Add the G4FaserAlg +# + from G4FaserAlg.G4FaserAlgConfigNew import G4FaserAlgCfg + cfg.merge(G4FaserAlgCfg(ConfigFlags, InputTruthCollection = MCEventKey)) +# +# Dump config +# + from AthenaConfiguration.ComponentFactory import CompFactory + cfg.addEventAlgo(CompFactory.JobOptsDumperAlg(FileName="G4FaserTestConfig.txt")) + cfg.getService("StoreGateSvc").Dump = True + cfg.getService("ConditionStore").Dump = True + cfg.printConfig(withDetails=True, summariseProps = False) # gags on ParticleGun if summariseProps = True? + + ConfigFlags.dump() + #f = open("test.pkl","wb") + #cfg.store(f) + #f.close() +# +# Execute and finish +# + + #cfg.foreach_component("*").OutputLevel = "INFO" # Use warning for production + + sc = cfg.run() + + b = time.time() + log.info("Run G4FaserAlg in " + str(b-a) + " seconds") +# +# Success should be 0 +# + sys.exit(not sc.isSuccess()) + diff --git a/Control/CalypsoExample/Generation/scripts/faser_particlegun.py b/Control/CalypsoExample/Generation/scripts/faser_particlegun.py new file mode 100755 index 0000000000000000000000000000000000000000..ee38005df79c3f6aab97b5c800f3949521110a02 --- /dev/null +++ b/Control/CalypsoExample/Generation/scripts/faser_particlegun.py @@ -0,0 +1,178 @@ +#!/usr/bin/env python +""" +Produce particle gun samples +Derived from G4FaserAlgConfigNew +This is a general low-level script, +although this could be useful for testing. + +Usage: +faser_particlegun.py <options> + +Control using command-line options: +--evtMax=1000 +--skipEvt=1000 + +Output.HITSFileName=<name> +Sim.Gun='{"pid" : 11, "z": -1500.}' + +Copyright (C) 2002-2021 CERN for the benefit of the ATLAS and FASER collaborations +""" + +if __name__ == '__main__': + + import time + a = time.time() +# +# Set up logging and config behaviour +# + from AthenaCommon.Logging import log + from AthenaCommon.Constants import DEBUG, VERBOSE + from AthenaCommon.Configurable import Configurable + log.setLevel(DEBUG) + Configurable.configurableRun3Behavior = 1 +# +# Import and set config flags +# + from CalypsoConfiguration.AllConfigFlags import ConfigFlags + ConfigFlags.Exec.MaxEvents = 10 # can be overridden from command line with --evtMax=<number> + ConfigFlags.Exec.SkipEvents = 0 # can be overridden from command line with --skipEvt=<number> + from AthenaConfiguration.Enums import ProductionStep + ConfigFlags.Common.ProductionStep = ProductionStep.Simulation +# +# All these must be specified to avoid auto-configuration +# + ConfigFlags.Input.RunNumber = [12345] #Isn't updating - todo: investigate + ConfigFlags.Input.OverrideRunNumber = True + ConfigFlags.Input.LumiBlockNumber = [1] + ConfigFlags.Input.isMC = True +# +# Output file name +# + ConfigFlags.Output.HITSFileName = "my.HITS.pool.root" # can be overridden from command line with Output.HITSFileName=<name> +# +# Sim ConfigFlags +# + ConfigFlags.Sim.Layout = "FASER" + ConfigFlags.Sim.PhysicsList = "FTFP_BERT" + ConfigFlags.Sim.ReleaseGeoModel = False + ConfigFlags.Sim.IncludeParentsInG4Event = True # Controls whether BeamTruthEvent is written to output HITS file + ConfigFlags.addFlag("Sim.Gun",{"Generator" : "SingleParticle"}) # Property bag for particle gun keyword:argument pairs + ConfigFlags.addFlag("Sim.Beam.xangle", 0) # Potential beam crossing angles + ConfigFlags.addFlag("Sim.Beam.yangle", 0) + + ConfigFlags.GeoModel.FaserVersion = "FASERNU-02" # Geometry set-up + ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-02" # Conditions set-up + ConfigFlags.addFlag("Input.InitialTimeStamp", 0) # To avoid autoconfig + ConfigFlags.GeoModel.Align.Dynamic = False + +# +# Preset particle gun parameters +# + from math import atan + from AthenaCommon.SystemOfUnits import GeV, TeV, cm, m + from AthenaCommon.PhysicalConstants import pi + + # 11 - electron, 13 - muon, 22 - photon + import ParticleGun as PG + ConfigFlags.Sim.Gun = { + "Generator" : "SingleParticle", "pid" : 13, + "energy" : PG.LogSampler(10*GeV, 1*TeV), + "theta" : PG.GaussianSampler(0, atan((10*cm)/(7*m)), oneside = True), + "phi" : [0, 2*pi], "mass" : 0.511, "radius" : -10*cm, #"z": -2.0*m, + "randomSeed" : 12345} + +# +# Command-line overrides +# +# These have the format: Sim.Gun='{"pid" : 11}' +# Filename: Output.HITSFileName="test.muon.001.root" +# Also number of events: --evtMax 100 +# Starting at z = -1.5m (-1500.) will miss the veto (for electrons) + + import sys + ConfigFlags.fillFromArgs(sys.argv[1:]) + +# +# By being a little clever, we can steer the geometry setup from the command line using GeoModel.FaserVersion +# +# MDC configuration +# + detectors = ['Veto', 'Preshower', 'FaserSCT', 'Ecal', 'Trigger', 'Dipole', 'Emulsion'] +# +# Setup detector flags +# + from CalypsoConfiguration.DetectorConfigFlags import setupDetectorsFromList + setupDetectorsFromList(ConfigFlags, detectors, toggle_geometry=True) +# +# Finalize flags +# + ConfigFlags.lock() +# +# Initialize a new component accumulator +# + from CalypsoConfiguration.MainServicesConfig import MainServicesCfg + cfg = MainServicesCfg(ConfigFlags) +# +# Configure the particle gun as requested, or using defaults +# + +# +# Particle gun generators - the generator, energy, angle, particle type, position, etc can be modified by passing keyword arguments +# + from FaserParticleGun.FaserParticleGunConfig import FaserParticleGunCfg + cfg.merge(FaserParticleGunCfg(ConfigFlags)) + from McEventSelector.McEventSelectorConfig import McEventSelectorCfg + cfg.merge(McEventSelectorCfg(ConfigFlags)) + +# +# Output file +# + from AthenaPoolCnvSvc.PoolWriteConfig import PoolWriteCfg + cfg.merge(PoolWriteCfg(ConfigFlags)) + +# +# Shift LOS +# + + if ConfigFlags.Sim.Beam.xangle or ConfigFlags.Sim.Beam.yangle: + MCEventKey = "BeamTruthEventShifted" + import McParticleEvent.Pythonizations + from GeneratorUtils.ShiftLOSConfig import ShiftLOSCfg + cfg.merge(ShiftLOSCfg(ConfigFlags, OutputMCEventKey = MCEventKey, + xcross = ConfigFlags.Sim.Beam.xangle, ycross = ConfigFlags.Sim.Beam.yangle)) + else: + MCEventKey = "BeamTruthEvent" + +# +# Add the G4FaserAlg +# + from G4FaserAlg.G4FaserAlgConfigNew import G4FaserAlgCfg + cfg.merge(G4FaserAlgCfg(ConfigFlags, InputTruthCollection = MCEventKey)) +# +# Dump config +# + from AthenaConfiguration.ComponentFactory import CompFactory + cfg.addEventAlgo(CompFactory.JobOptsDumperAlg(FileName="G4FaserTestConfig.txt")) + cfg.getService("StoreGateSvc").Dump = True + cfg.getService("ConditionStore").Dump = True + cfg.printConfig(withDetails=True, summariseProps = False) # gags on ParticleGun if summariseProps = True? + + ConfigFlags.dump() + #f = open("test.pkl","wb") + #cfg.store(f) + #f.close() +# +# Execute and finish +# + + #cfg.foreach_component("*").OutputLevel = "INFO" # Use warning for production + + sc = cfg.run() + + b = time.time() + log.info("Run G4FaserAlg in " + str(b-a) + " seconds") +# +# Success should be 0 +# + sys.exit(not sc.isSuccess()) + diff --git a/Control/CalypsoExample/Generation/scripts/submit_faserMDC_particlegun.sh b/Control/CalypsoExample/Generation/scripts/submit_faserMDC_particlegun.sh new file mode 100755 index 0000000000000000000000000000000000000000..1aa00c891ed46604d5ba294e6ea180c2ba00163a --- /dev/null +++ b/Control/CalypsoExample/Generation/scripts/submit_faserMDC_particlegun.sh @@ -0,0 +1,134 @@ +#!/bin/bash +# Used with a condor file to submit to vanilla universe +# +# Usage: +# submit_faserMDC_particlegun.sh config_file segment [release_directory] [working_directory] +# +# config_file - full file name (with path) +# segment - segment number (file segment) +# release_directory - optional path to release install directory (default pwd) +# working_directory - optional path to output directory location (default pwd) +# +# 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 gen/g???? it will be passed to the job +# +#---------------------------------------- +# +# Parse command-line options +config_path=${1} +segment=${2} +release_directory=${3} +working_directory=${4} +# +# Set defaults if arguments aren't provided +if [ -z "$config_path" ] +then + echo "No config_path specified!" + echo "Usage: submit_faserMDC_particlegun.sh config_file segment [release dir] [output dir]" + exit 1 +fi +# +if [ -z "$segment" ] +then + segment=0 +fi +# +if [ -z "$release_directory" ] +then + release_directory=`pwd` +fi +# +if [ -z "$working_directory" ] +then + working_directory=`pwd` +fi +# +# Apply padding to segment number +printf -v seg_str "%05d" $segment +# +starting_directory=`pwd` +# +# Now extract the file stem +# +# First, get the filename +config_file=$(basename "$config_path") +# +# Now split based on '.' to get stem +defaultIFS=$IFS +IFS='.' +read config_file_stem ext <<< "$config_file" +# Set the IFS delimeter back or else echo doesn't work... +IFS=$defaultIFS +# +# Make output directory if needed +output_directory="$working_directory/${config_file_stem}" +mkdir -p "$output_directory" +# +# This magic redirects everything in this script to our log file +exec >& "$output_directory/${config_file_stem}-${seg_str}.log" +echo `date` - $HOSTNAME +echo "File: $config_file" +echo "Segment: $seg_str" +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" +# This doesn't seem to work, as we need the --input argument +#asetup +#source build/x8*/setup.sh +# +# Do this by hand +asetup --input=calypso/asetup.faser Athena,22.0.49 +source build/x86*/setup.sh +# +# +# Try to find a release tag +cd calypso +gentag=`git describe` +if [[ "$gentag" == "gen/g"???? ]]; then + gtag=`echo "$gentag" | cut -c 5-10` + echo "Found gen tag: $gtag" +fi +if [[ "$gentag" == "digi/d"???? ]]; then + tag=`echo "$recotag" | cut -c 6-11` + echo "Found reco tag: $tag" +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 "${config_file_stem}-${seg_str}" ]]; then + echo "Directory ${config_file_stem}-${seg_str} already exists" +else + mkdir "${config_file_stem}-${seg_str}" +fi +cd "${config_file_stem}-${seg_str}" +# +# Run job +if [[ -z "$gtag" ]]; then + faserMDC_particlegun.py "--conf=$config_path" "--segment=$seg_str" +else + faserMDC_particlegun.py "--conf=$config_path" "--segment=$seg_str" "--tag=$gtag" +fi +# +# Print out ending time +date diff --git a/Control/CalypsoExample/Reconstruction/scripts/faserMDC_reco.py b/Control/CalypsoExample/Reconstruction/scripts/faserMDC_reco.py new file mode 100755 index 0000000000000000000000000000000000000000..30b40afce32b6ce8e124a9d2fcee3b32f0f7ae6e --- /dev/null +++ b/Control/CalypsoExample/Reconstruction/scripts/faserMDC_reco.py @@ -0,0 +1,207 @@ +#!/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 - optionally specify the data type (TI12Data, TI12Data02, or TestBeamData). +# In a normal file system location, this will be extracted from the directory name, +# but runtype will override this assignment. +# TI12Data02 is needed for the IFT geometry. +# MDC will assume this geometry. +# +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)") +parser.add_argument("-v", "--verbose", action='store_true', + help="Turn on DEBUG output") + + +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 + +# Assume based on MDC reco +else: + runtype = "TI12Data02" + +# Assume this is MC +args.isMC = True + +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 = args.isMC # 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" or runtype == "TestBeam2021": + ConfigFlags.GeoModel.FaserVersion = "FASER-TB00" + ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-TB00" + +# New TI12 geometry (ugh) +elif runtype == "TI12Data02": + ConfigFlags.GeoModel.FaserVersion = "FASER-02" + 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 + +if args.isMC: + from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg + acc.merge(PoolReadCfg(ConfigFlags)) +else: + 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)) + +# Not ready for primetime +# from CaloRecAlgs.CaloRecAlgsConfig import CalorimeterReconstructionCfg +# acc.merge(CalorimeterReconstructionCfg(ConfigFlags)) + +# Tracker clusters +from TrackerPrepRawDataFormation.TrackerPrepRawDataFormationConfig import FaserSCT_ClusterizationCfg +acc.merge(FaserSCT_ClusterizationCfg(ConfigFlags)) + +# SpacePoints +from TrackerSpacePointFormation.TrackerSpacePointFormationConfig import TrackerSpacePointFinderCfg +acc.merge(TrackerSpacePointFinderCfg(ConfigFlags)) + +print("Configuring TrackerSegmentFit (new)") +# Try Dave's new fitter +from TrackerSegmentFit.TrackerSegmentFitConfig import SegmentFitAlgCfg +acc.merge(SegmentFitAlgCfg(ConfigFlags)) + +# +# Configure output +from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg +itemList = [ "xAOD::EventInfo#*" + , "xAOD::EventAuxInfo#*" + , "xAOD::FaserTriggerData#*" + , "xAOD::FaserTriggerDataAux#*" + , "FaserSCT_RDO_Container#*" + , "Tracker::FaserSCT_ClusterContainer#*" + , "FaserSCT_SpacePointContainer#*" + #, "FaserSCT_SpacePointOverlapCollection#*" + , "TrackCollection#*" +] +acc.merge(OutputStreamCfg(ConfigFlags, "xAOD", itemList)) + +# Waveform reconstruction output +from WaveRecAlgs.WaveRecAlgsConfig import WaveformReconstructionOutputCfg +acc.merge(WaveformReconstructionOutputCfg(ConfigFlags)) + +# Calorimeter reconstruction output +# from CaloRecAlgs.CaloRecAlgsConfig import CalorimeterReconstructionOutputCfg +# acc.merge(CalorimeterReconstructionOutputCfg(ConfigFlags)) + +# Check what we have +print( "Writing out xAOD objects:" ) +print( acc.getEventAlgo("OutputStreamxAOD").ItemList ) + +# Hack to avoid problem with our use of MC databases when isMC = False +if not args.isMC: + replicaSvc = acc.getService("DBReplicaSvc") + replicaSvc.COOLSQLiteVetoPattern = "" + replicaSvc.UseCOOLSQLite = True + replicaSvc.UseCOOLFrontier = False + replicaSvc.UseGeomSQLite = True + +# Configure verbosity +# ConfigFlags.dump() +if args.verbose: + acc.foreach_component("*").OutputLevel = 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 + +else: + acc.foreach_component("*").OutputLevel = INFO + +acc.foreach_component("*ClassID*").OutputLevel = INFO + +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_faserMDC_reco.sh b/Control/CalypsoExample/Reconstruction/scripts/submit_faserMDC_reco.sh new file mode 100755 index 0000000000000000000000000000000000000000..60c39ca370330fcad5432523d33740722b754a8f --- /dev/null +++ b/Control/CalypsoExample/Reconstruction/scripts/submit_faserMDC_reco.sh @@ -0,0 +1,134 @@ +#!/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 desc 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_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=calypso/asetup.faser Athena,22.0.49 +source build/x86*/setup.sh +# +# +# Try to find a release tag +cd calypso +recotag=`git describe` +if [[ "$recotag" == "reco/r"???? ]]; then + tag=`echo "$recotag" | cut -c 6-11` + echo "Found reco tag: $tag" +fi +if [[ "$recotag" == "reco/p"???? ]]; then + tag=`echo "$recotag" | cut -c 6-11` + echo "Found reco tag: $tag" +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 + faserMDC_reco.py "--nevents=$nevents" "$file_path" +else + faserMDC_reco.py "--nevents=$nevents" "--reco=$tag" "$file_path" +fi +# +# Print out ending time +date diff --git a/Control/CalypsoExample/Reconstruction/scripts/submit_faser_reco.sh b/Control/CalypsoExample/Reconstruction/scripts/submit_faser_reco.sh index 64cca5567690c0680a4966b7e64f1d8ea869b1ca..277d6eaa79a3bc510cb413aa742b292837e28403 100755 --- a/Control/CalypsoExample/Reconstruction/scripts/submit_faser_reco.sh +++ b/Control/CalypsoExample/Reconstruction/scripts/submit_faser_reco.sh @@ -97,8 +97,12 @@ source build/x8*/setup.sh cd calypso recotag=`git describe` if [[ "$recotag" == "reco/r"???? ]]; then - rtag=`echo "$recotag" | cut -c 6-11` - echo "Found reco tag: $rtag" + tag=`echo "$recotag" | cut -c 6-11` + echo "Found reco tag: $tag" +fi +if [[ "$recotag" == "reco/p"???? ]]; then + tag=`echo "$recotag" | cut -c 6-11` + echo "Found reco tag: $tag" fi # # Move to the run directory @@ -123,7 +127,7 @@ cd "$file_stem" if [[ -z "$rtag" ]]; then faser_reco.py "--nevents=$nevents" "$file_path" else - faser_reco.py "--nevents=$nevents" "--reco=$rtag" "$file_path" + faser_reco.py "--nevents=$nevents" "--reco=$tag" "$file_path" fi # # Print out ending time diff --git a/Scintillator/ScintDigiAlgs/CMakeLists.txt b/Scintillator/ScintDigiAlgs/CMakeLists.txt index 3ebe3d691decc51518eb6f64ecae03efd76aaf73..5c1874e2502241f919efb21ade0e7c26de7369ca 100644 --- a/Scintillator/ScintDigiAlgs/CMakeLists.txt +++ b/Scintillator/ScintDigiAlgs/CMakeLists.txt @@ -9,7 +9,9 @@ atlas_subdir( ScintDigiAlgs ) atlas_add_component( ScintDigiAlgs src/*.cxx src/*.h src/components/*.cxx - LINK_LIBRARIES AthenaBaseComps Identifier ScintIdentifier StoreGateLib WaveRawEvent ScintSimEvent WaveDigiToolsLib) + LINK_LIBRARIES AthenaBaseComps Identifier ScintIdentifier + WaveformConditionsToolsLib StoreGateLib WaveRawEvent + ScintSimEvent WaveDigiToolsLib) atlas_install_python_modules( python/*.py ) diff --git a/Scintillator/ScintDigiAlgs/python/ScintDigiAlgsConfig.py b/Scintillator/ScintDigiAlgs/python/ScintDigiAlgsConfig.py index 292a19c0890672bb3a1530b080d936c8716a2fff..ddf3074618aff121e704b2214309a14952436837 100644 --- a/Scintillator/ScintDigiAlgs/python/ScintDigiAlgsConfig.py +++ b/Scintillator/ScintDigiAlgs/python/ScintDigiAlgsConfig.py @@ -5,6 +5,8 @@ from AthenaConfiguration.ComponentFactory import CompFactory from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg +from WaveformConditionsTools.WaveformCableMappingConfig import WaveformCableMappingCfg + # Crystallball parameter dictionary used in simulated digitized wave reconstruction. # Crystalball function Parameters estimated from Deion's slides uploaded at # https://indico.cern.ch/event/1099652/contributions/4626975/attachments/2352595/4013927/Faser-Physics-run3933-plots.pdf (20/01/2022) @@ -13,12 +15,14 @@ dict_CB_param = {} dict_CB_param["Trigger"]=dict(CB_alpha=-0.38, CB_n=25, CB_mean=815, CB_sigma=7.7, CB_norm = 500 ) dict_CB_param["Timing"]=dict(CB_alpha=-0.32, CB_n=65, CB_mean=846, CB_sigma=5.3, CB_norm = 500) # copy from Preshower; Timing was not in TestBeam dict_CB_param["Veto"]=dict(CB_alpha=-0.38, CB_n=25, CB_mean=815, CB_sigma=7.7, CB_norm = 1000) # copy from Trigger; Veto was not in TestBeam, but in sim "Veto" is the TestBeam Trigger component +dict_CB_param["VetoNu"]=dict(CB_alpha=-0.38, CB_n=25, CB_mean=815, CB_sigma=7.7, CB_norm = 1000) # copy from Trigger; Veto was not in TestBeam, but in sim "Veto" is the TestBeam Trigger component dict_CB_param["Preshower"]=dict(CB_alpha=-0.32, CB_n=65, CB_mean=846, CB_sigma=5.3, CB_norm = 500) dict_baseline_params = { "Trigger" : {"mean" : 15000, "rms" : 3}, "Timing" : {"mean" : 15000, "rms" : 3}, "Veto" : {"mean" : 15000, "rms" : 3}, + "VetoNu" : {"mean" : 15000, "rms" : 3}, "Preshower" : {"mean" : 15000, "rms" : 3}, } @@ -33,8 +37,10 @@ def ScintWaveformDigitizationCfg(flags): if "TB" not in flags.GeoModel.FaserVersion: acc.merge(ScintWaveformDigiCfg(flags, "TimingWaveformDigiAlg", "Trigger")) acc.merge(ScintWaveformDigiCfg(flags, "VetoWaveformDigiAlg", "Veto")) + acc.merge(ScintWaveformDigiCfg(flags, "VetoNuWaveformDigiAlg", "VetoNu")) acc.merge(ScintWaveformDigiCfg(flags, "PreshowerWaveformDigiAlg", "Preshower")) acc.merge(ScintWaveformDigitizationOutputCfg(flags)) + acc.merge(WaveformCableMappingCfg(flags)) return acc # Return configured digitization algorithm from SIM hits @@ -72,6 +78,7 @@ def ScintWaveformDigitizationOutputCfg(flags, **kwargs): ] acc.merge(OutputStreamCfg(flags, "RDO")) ostream = acc.getEventAlgo("OutputStreamRDO") - ostream.TakeItemsFromInput = True # Copies all data from input file to output + # ostream.TakeItemsFromInput = True # Copies all data from input file to output + # ostream.TakeItemsFromInput = False ostream.ItemList += ItemList return acc diff --git a/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx b/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx index d5235a92ccdacbc2631e4dafc66a3f9dea48dc9d..7eb12286567c42591314592ece065136a754d43d 100644 --- a/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx +++ b/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx @@ -17,6 +17,7 @@ ScintWaveformDigiAlg::initialize() { // Initalize tools ATH_CHECK( m_digiTool.retrieve() ); + ATH_CHECK( m_mappingTool.retrieve() ); // Set key to read waveform from ATH_CHECK( m_scintHitContainerKey.initialize() ); @@ -26,6 +27,7 @@ ScintWaveformDigiAlg::initialize() { // Set up helpers ATH_CHECK(detStore()->retrieve(m_vetoID, "VetoID")); + ATH_CHECK(detStore()->retrieve(m_vetoNuID, "VetoNuID")); ATH_CHECK(detStore()->retrieve(m_triggerID, "TriggerID")); ATH_CHECK(detStore()->retrieve(m_preshowerID, "PreshowerID")); @@ -80,6 +82,8 @@ ScintWaveformDigiAlg::execute(const EventContext& ctx) const { auto first = *scintHitHandle->begin(); if (first.isVeto()) { waveforms = m_digiTool->create_waveform_map(m_vetoID); + } else if (first.isVetoNu()) { + waveforms = m_digiTool->create_waveform_map(m_vetoNuID); } else if (first.isTrigger()) { waveforms = m_digiTool->create_waveform_map(m_triggerID); } else if (first.isPreshower()) { @@ -87,12 +91,34 @@ ScintWaveformDigiAlg::execute(const EventContext& ctx) const { } // Loop over time samples + // Should really do sums first, as repeating these in the loop is a waste for (const auto& tk : m_timekernel) { - std::map<unsigned int, float> counts; + std::map<Identifier, float> counts; // Convolve hit energy with evaluated kernel and sum for each hit id (i.e. channel) + Identifier id; for (const auto& hit : *scintHitHandle) { - counts[hit.identify()] += tk * hit.energyLoss(); + + // Special handling for trigger scintillator + if (first.isTrigger()) { + // Hits are created per plate + // Need to split between 2 PMTs + + // Identifier from hit is plate ID + Identifier plate_id = m_triggerID->plate_id(hit.getIdentifier()); + + // Need to do something better here, but for now just fill both + id = m_triggerID->pmt_id(plate_id, 0); // ID for PMT 0 + counts[id] += tk * hit.energyLoss(); + id = m_triggerID->pmt_id(plate_id, 1); // ID for PMT 1 + counts[id] += tk * hit.energyLoss(); + + } else { + // All others have only 1 PMT + // Use detector id (not hit id) throughout + id = hit.getIdentifier(); + counts[id] += tk * hit.energyLoss(); + } } // Subtract count from basleine and add result to correct waveform vector @@ -107,11 +133,27 @@ ScintWaveformDigiAlg::execute(const EventContext& ctx) const { } // Convert hit id to Identifier and store - Identifier id = ScintHitIdHelper::GetHelper()->getIdentifier(c.first); + id = c.first; waveforms[id].push_back(value); } } + // This is a bit of a hack to make sure all waveforms have + // at least baseline entries. Should really add this to the + // logic above + for (const auto& w : waveforms) { + if (w.second.size() > 0) continue; + + // Waveform was empty, fill with baseline + int channel = m_mappingTool->getChannelMapping(w.first); + ATH_MSG_DEBUG("Writing baseline into empty waveform in channel "<< channel); + int i = m_digiTool->nsamples(); + while(i--) { // Use while to avoid unused variable warning with for + int baseline = m_digiTool->generate_baseline(m_base_mean, m_base_rms); + waveforms[w.first].push_back(baseline); + } + } + // Loop over wavefrom vectors to make and store raw waveform unsigned int nsamples = m_digiTool->nsamples(); for (const auto& w : waveforms) { @@ -119,6 +161,14 @@ ScintWaveformDigiAlg::execute(const EventContext& ctx) const { wfm->setWaveform(0, w.second); wfm->setIdentifier(w.first); wfm->setSamples(nsamples); + // Only save this waveform if channel mapping is valid + // This will avoid making a waveform for the missing veto counter + int channel = m_mappingTool->getChannelMapping(w.first); + if (channel < 0) { + ATH_MSG_DEBUG("Skipping waveform with unknown channel!"); + continue; + } + wfm->setChannel(channel); waveformContainerHandle->push_back(wfm); } diff --git a/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.h b/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.h index de03370698620986377c695852af07c7d8c79b18..528325621b279bdc3f570c6e5ade93fe549b6c12 100644 --- a/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.h +++ b/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.h @@ -10,6 +10,7 @@ // Tool classes #include "WaveDigiTools/IWaveformDigitisationTool.h" +#include "WaveformConditionsTools/IWaveformCableMappingTool.h" // Handles #include "StoreGate/ReadHandleKey.h" @@ -21,6 +22,7 @@ //Helpers #include "ScintIdentifier/VetoID.h" +#include "ScintIdentifier/VetoNuID.h" #include "ScintIdentifier/TriggerID.h" #include "ScintIdentifier/PreshowerID.h" #include "ScintSimEvent/ScintHitIdHelper.h" @@ -83,6 +85,7 @@ class ScintWaveformDigiAlg : public AthReentrantAlgorithm { /// Detector ID helpers const VetoID* m_vetoID{nullptr}; + const VetoNuID* m_vetoNuID{nullptr}; const TriggerID* m_triggerID{nullptr}; const PreshowerID* m_preshowerID{nullptr}; @@ -92,6 +95,11 @@ class ScintWaveformDigiAlg : public AthReentrantAlgorithm { ToolHandle<IWaveformDigitisationTool> m_digiTool {this, "WaveformDigitisationTool", "WaveformDigitisationTool"}; + /** + * @name Mapping tool + */ + ToolHandle<IWaveformCableMappingTool> m_mappingTool + {this, "WaveformCableMappingTool", "WaveformCableMappingTool"}; /** * @name Input HITS using SG::ReadHandleKey diff --git a/Scintillator/ScintSimEvent/ScintSimEvent/ScintHit.h b/Scintillator/ScintSimEvent/ScintSimEvent/ScintHit.h index a2e02b596f4485d330a24ae6a5190832c9f0fa24..864217ffb2628ff0873ae4823345b6f5e7e10621 100644 --- a/Scintillator/ScintSimEvent/ScintSimEvent/ScintHit.h +++ b/Scintillator/ScintSimEvent/ScintSimEvent/ScintHit.h @@ -16,6 +16,8 @@ #include "CLHEP/Geometry/Point3D.h" #include "GeneratorObjects/HepMcParticleLink.h" +class Identifier; + class ScintHit { /////////////////////////////////////////////////////////////////// @@ -80,8 +82,13 @@ public: // Const methods: /////////////////////////////////////////////////////////////////// + // This is the HitId, used in Geant. This is not the detector Identifier unsigned int identify() const; + // This is the detector readout Identifier (pmt Identifier) + // provided by ScintHitIdHelper::getIdentifier + Identifier getIdentifier() const; + // local start position of the energy deposit: HepGeom::Point3D<double> localStartPosition() const; diff --git a/Scintillator/ScintSimEvent/src/ScintHit.cxx b/Scintillator/ScintSimEvent/src/ScintHit.cxx index 3f5506ed8af50ea471fbba91e8d5c674ab803648..8c2201a133f75676c44c284c61061556048d9600 100644 --- a/Scintillator/ScintSimEvent/src/ScintHit.cxx +++ b/Scintillator/ScintSimEvent/src/ScintHit.cxx @@ -142,6 +142,10 @@ bool ScintHit::isVetoNu() const { return ScintHitIdHelper::GetHelper()->isVetoNu(m_ID); } +Identifier ScintHit::getIdentifier() const { + return ScintHitIdHelper::GetHelper()->getIdentifier(m_ID); +} + HepGeom::Point3D<double> ScintHit::localStartPosition() const { return HepGeom::Point3D<double>((double) m_stX, (double) m_stY, (double) m_stZ); diff --git a/Tracker/TrackerDigitization/FaserSCT_Digitization/python/FaserSCT_DigitizationConfigNew.py b/Tracker/TrackerDigitization/FaserSCT_Digitization/python/FaserSCT_DigitizationConfigNew.py index c9f267da7d9d2ba0c6bc973c3fa949b65c217ecc..6ae2b277c89ed09edee634d4fdc6808faec2fa28 100644 --- a/Tracker/TrackerDigitization/FaserSCT_Digitization/python/FaserSCT_DigitizationConfigNew.py +++ b/Tracker/TrackerDigitization/FaserSCT_Digitization/python/FaserSCT_DigitizationConfigNew.py @@ -242,7 +242,7 @@ def FaserSCT_OutputCfg(flags): acc.merge(TruthDigitizationOutputCfg(flags)) acc.merge(OutputStreamCfg(flags, "RDO")) ostream = acc.getEventAlgo("OutputStreamRDO") - ostream.TakeItemsFromInput = True + #ostream.TakeItemsFromInput = True # Don't write hits to RDO by default ostream.ItemList += ItemList return acc diff --git a/Waveform/WaveRecAlgs/python/WaveRecAlgsConfig.py b/Waveform/WaveRecAlgs/python/WaveRecAlgsConfig.py index 7d3219eb0da3e1026fae15754bb9f52bd1618c65..414773bf9dc2e9b8140b33d81a80220c4288e4bf 100644 --- a/Waveform/WaveRecAlgs/python/WaveRecAlgsConfig.py +++ b/Waveform/WaveRecAlgs/python/WaveRecAlgsConfig.py @@ -12,30 +12,18 @@ WaveformReconstructionTool = CompFactory.WaveformReconstructionTool ClockReconstructionTool = CompFactory.ClockReconstructionTool # One stop shopping for normal FASER data -def WaveformReconstructionCfg(flags, naive = False): +def WaveformReconstructionCfg(flags): """ Return all algorithms and tools for Waveform reconstruction """ acc = ComponentAccumulator() if not flags.Input.isMC: acc.merge(WaveformClockRecCfg(flags, "ClockRecAlg")) - if flags.Input.isMC and naive: - if "TB" not in flags.GeoModel.FaserVersion: - acc.merge(PseudoScintHitToWaveformRecCfg(flags, "PseudoTriggerHitWaveformRecAlg", "Trigger")) - - acc.merge(PseudoScintHitToWaveformRecCfg(flags, "PseudoVetoHitToWaveformRecAlg", "Veto")) - acc.merge(PseudoScintHitToWaveformRecCfg(flags, "PseudoPresehowerHitWaveformRecAlg", "Preshower")) - acc.merge(PseudoCaloHitToWaveformRecCfg(flags, "PseudoCaloHitWaveformRecAlg")) - return acc - acc.merge(WaveformHitRecCfg(flags, "TriggerWaveformRecAlg", "Trigger")) acc.merge(WaveformHitRecCfg(flags, "VetoWaveformRecAlg", "Veto")) - if flags.Input.isMC: - print("Turning off VetoNu reco in MC!") - else: - acc.merge(WaveformHitRecCfg(flags, "VetoNuWaveformRecAlg", "VetoNu")) acc.merge(WaveformHitRecCfg(flags, "PreshowerWaveformRecAlg", "Preshower")) acc.merge(WaveformHitRecCfg(flags, "CaloWaveformRecAlg", "Calo")) + acc.merge(WaveformHitRecCfg(flags, "VetoNuWaveformRecAlg", "VetoNu")) acc.merge(WaveformTimingCfg(flags)) diff --git a/Waveform/WaveRecTools/CMakeLists.txt b/Waveform/WaveRecTools/CMakeLists.txt index f7f9672886ef12ab2fd9441447bf0f0a0f583418..7c42c848a64280325967fc058b7d9e8ac31170e3 100644 --- a/Waveform/WaveRecTools/CMakeLists.txt +++ b/Waveform/WaveRecTools/CMakeLists.txt @@ -14,7 +14,7 @@ atlas_add_library( WaveRecToolsLib PUBLIC_HEADERS WaveRecTools PRIVATE_INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} LINK_LIBRARIES AthenaBaseComps AthenaKernel GeoPrimitives - WaveformConditionsToolsLib WaveRawEvent xAODFaserWaveform + WaveformConditionsToolsLib WaveRawEvent xAODFaserWaveform Identifier PRIVATE_LINK_LIBRARIES ${ROOT_LIBRARIES} ) diff --git a/Waveform/WaveRecTools/WaveRecTools/IPseudoSimToWaveformRecTool.h b/Waveform/WaveRecTools/WaveRecTools/IPseudoSimToWaveformRecTool.h index 6ba6e4eacdb8288a5f3667b02c3823aeaeb54d7e..fda076b2ae7a2a22bc5743b32eded9b92a006da3 100644 --- a/Waveform/WaveRecTools/WaveRecTools/IPseudoSimToWaveformRecTool.h +++ b/Waveform/WaveRecTools/WaveRecTools/IPseudoSimToWaveformRecTool.h @@ -23,6 +23,8 @@ #include "xAODFaserWaveform/WaveformHitContainer.h" #include "xAODFaserWaveform/WaveformHit.h" +class Identifier; + ///Interface for Pseudo waveform rec tools class IPseudoSimToWaveformRecTool : virtual public IAlgTool { @@ -31,9 +33,11 @@ public: // InterfaceID DeclareInterfaceID(IPseudoSimToWaveformRecTool, 1, 0); + /* IPseudoSimToWaveformRecTool(): m_msgSvc ( "MessageSvc", "ITrkEventCnvTool" ) {} + */ virtual ~IPseudoSimToWaveformRecTool() = default; @@ -42,8 +46,9 @@ public: StatusCode reconstruct(const CONT* hitCollection, xAOD::WaveformHitContainer* waveContainer) const; -private: - ServiceHandle<IMessageSvc> m_msgSvc; + // Make the actual hits (separate this so it doesn't need to be in templated function + virtual StatusCode make_hits(const std::map<Identifier, float>& sums, + xAOD::WaveformHitContainer* waveContainer) const = 0; }; diff --git a/Waveform/WaveRecTools/WaveRecTools/IPseudoSimToWaveformRecTool.icc b/Waveform/WaveRecTools/WaveRecTools/IPseudoSimToWaveformRecTool.icc index 9f85c6b020b3efc9d81949cbf914b2a424bda19b..29a689bdc2e3c7c665f08df1462ba6afb237fd84 100644 --- a/Waveform/WaveRecTools/WaveRecTools/IPseudoSimToWaveformRecTool.icc +++ b/Waveform/WaveRecTools/WaveRecTools/IPseudoSimToWaveformRecTool.icc @@ -2,35 +2,12 @@ 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; + std::map<Identifier, 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); + idSums[hit.getIdentifier()] += hit.energyLoss(); } - return StatusCode::SUCCESS; + return make_hits(idSums, container); } diff --git a/Waveform/WaveRecTools/src/PseudoSimToWaveformRecTool.cxx b/Waveform/WaveRecTools/src/PseudoSimToWaveformRecTool.cxx index 6d894aeb64125274507105bd71ce37c7ad2bed31..93ae1397e3c591b2f54b3765f0138042878e25bc 100644 --- a/Waveform/WaveRecTools/src/PseudoSimToWaveformRecTool.cxx +++ b/Waveform/WaveRecTools/src/PseudoSimToWaveformRecTool.cxx @@ -10,18 +10,9 @@ #include "PseudoSimToWaveformRecTool.h" +#include "Identifier/Identifier.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) @@ -32,7 +23,43 @@ PseudoSimToWaveformRecTool::PseudoSimToWaveformRecTool(const std::string& type, StatusCode PseudoSimToWaveformRecTool::initialize() { ATH_MSG_INFO( name() << "::initalize()" ); + + ATH_CHECK( m_mappingTool.retrieve() ); + return StatusCode::SUCCESS; } +StatusCode +PseudoSimToWaveformRecTool::make_hits(const std::map<Identifier, float>& sums, + xAOD::WaveformHitContainer* container) const { + +// Check the container + if (!container) { + ATH_MSG_ERROR("HitCollection passed to reconstruct() is null!"); + return StatusCode::FAILURE; + } + + // Get the nominal trigger time (in ns) from config + // For now, hack up the timing to match the configuration + //float trigger_time = m_timingTool->nominalTriggerTime(); + //float offset; + + for (const auto& id : sums) { + xAOD::WaveformHit* hit = new xAOD::WaveformHit(); + container->push_back(hit); + + hit->set_identifier(id.first); + hit->set_channel(m_mappingTool->getChannelMapping(id.first)); + hit->set_peak(0); + hit->set_mean(0); + hit->set_width(0); + hit->set_integral(id.second); + hit->set_localtime(0); + hit->set_trigger_time(0); + hit->set_raw_peak(0); + hit->set_raw_integral(0); + } + + return StatusCode::SUCCESS; +} diff --git a/Waveform/WaveRecTools/src/PseudoSimToWaveformRecTool.h b/Waveform/WaveRecTools/src/PseudoSimToWaveformRecTool.h index f3376f6ec4aa93becef3da8fc52c19d003d6dced..e117479737d3dad28c1a8dfb7065ab1742d13884 100644 --- a/Waveform/WaveRecTools/src/PseudoSimToWaveformRecTool.h +++ b/Waveform/WaveRecTools/src/PseudoSimToWaveformRecTool.h @@ -13,10 +13,13 @@ #include "AthenaBaseComps/AthAlgTool.h" #include "WaveRecTools/IPseudoSimToWaveformRecTool.h" +// Tool classes +#include "WaveformConditionsTools/IWaveformCableMappingTool.h" + //Gaudi #include "GaudiKernel/ToolHandle.h" -//STL +class Identifier; class PseudoSimToWaveformRecTool: public extends<AthAlgTool, IPseudoSimToWaveformRecTool> { public: @@ -28,8 +31,12 @@ class PseudoSimToWaveformRecTool: public extends<AthAlgTool, IPseudoSimToWavefor /// Retrieve the necessary services in initialize StatusCode initialize(); + virtual StatusCode make_hits(const std::map<Identifier, float>& sums, + xAOD::WaveformHitContainer* waveContainer) const; + private: - // None + ToolHandle<IWaveformCableMappingTool> m_mappingTool + {this, "WaveformCableMappingTool", "WaveformCableMappingTool"}; }; diff --git a/Waveform/WaveRecTools/src/WaveformReconstructionTool.cxx b/Waveform/WaveRecTools/src/WaveformReconstructionTool.cxx index 25cdc0a3928b08ef0956698d9f2667ce89286e75..8882027bfc13a7b28a7043f2b9feb70ed8e25877 100644 --- a/Waveform/WaveRecTools/src/WaveformReconstructionTool.cxx +++ b/Waveform/WaveRecTools/src/WaveformReconstructionTool.cxx @@ -60,7 +60,7 @@ WaveformReconstructionTool::reconstructPrimary( // Set digitizer channel and identifier newhit->set_channel(wave.channel()); - newhit->set_id(wave.identify32().get_compact()); + newhit->set_identifier(wave.identify()); // Make sure we have ADC counts if (wave.adc_counts().size() == 0) { @@ -260,7 +260,7 @@ WaveformReconstructionTool::reconstructSecondary( // Fill values newhit->set_channel(wave.channel()); - newhit->set_id(wave.identify32().get_compact()); + newhit->set_identifier(wave.identify()); newhit->set_status_bit(xAOD::WaveformStatus::SECONDARY); newhit->set_baseline_mean(baseline.mean); newhit->set_baseline_rms(baseline.rms); @@ -737,7 +737,7 @@ WaveformReconstructionTool::fitGaussian(const WaveformFitResult& raw, const std: // Define fit function and preset range TF1 gfunc("gfunc", "gaus"); gfunc.SetParameters(raw.peak, raw.mean, raw.sigma); - gfunc.SetParError(0, std::sqrt(raw.peak)); + gfunc.SetParError(0, raw.peak/10.); gfunc.SetParError(1, raw.sigma); gfunc.SetParError(2, raw.sigma / 5.); gfunc.SetParLimits(2, 0., 20.); // Constrain width @@ -790,7 +790,7 @@ WaveformReconstructionTool::fitCBall(const WaveformFitResult& gfit, // Define fit function and preset values TF1 cbfunc("cbfunc", "crystalball"); cbfunc.SetParameter(0, cbfit.peak); // Peak height - cbfunc.SetParError(0, std::sqrt(cbfit.peak)); + cbfunc.SetParError(0, cbfit.peak/10.); cbfunc.SetParameter(1, cbfit.mean); // Mean cbfunc.SetParError(1, cbfit.sigma); cbfunc.SetParameter(2, cbfit.sigma); // Width diff --git a/Waveform/WaveformConditions/WaveformConditionsTools/src/WaveformCableMappingTool.cxx b/Waveform/WaveformConditions/WaveformConditionsTools/src/WaveformCableMappingTool.cxx index 530857354d1e1e82f5c522a72e8b6e0c7d6b5269..3922d0605339fce204cff01c9dbd69d3c6f21a97 100644 --- a/Waveform/WaveformConditions/WaveformConditionsTools/src/WaveformCableMappingTool.cxx +++ b/Waveform/WaveformConditions/WaveformConditionsTools/src/WaveformCableMappingTool.cxx @@ -143,7 +143,7 @@ WaveformCableMappingTool::getCableMapping(void) const { int WaveformCableMappingTool::getChannelMapping(const EventContext& ctx, const Identifier id) const { // Print where you are - ATH_MSG_DEBUG("in getChannelMapping()"); + ATH_MSG_DEBUG("in getChannelMapping("<< id <<")"); int channel = -1; // Read Cond Handle @@ -211,6 +211,10 @@ WaveformCableMappingTool::getChannelMapping(const EventContext& ctx, const Ident identifier = -1; } + // ATH_MSG_DEBUG("Comparing " << det_type << " ID: " + // << identifier.get_identifier() + // << " to requested " << id); + // Is this the identifier we are looking for? if (id != identifier) continue; @@ -222,9 +226,11 @@ WaveformCableMappingTool::getChannelMapping(const EventContext& ctx, const Ident } // End of loop over attributes - if (channel < 0) + if (channel < 0) { ATH_MSG_WARNING("No channel found for identifier " << id << "!"); + } + return channel; } diff --git a/xAOD/xAODFaserWaveform/xAODFaserWaveform/versions/WaveformHit_v1.h b/xAOD/xAODFaserWaveform/xAODFaserWaveform/versions/WaveformHit_v1.h index 8e497633d1d673e6bb6220552874e5dbea1ce0e1..e3615dbedd60ab1863d1d749cd7595dc961cacb2 100644 --- a/xAOD/xAODFaserWaveform/xAODFaserWaveform/versions/WaveformHit_v1.h +++ b/xAOD/xAODFaserWaveform/xAODFaserWaveform/versions/WaveformHit_v1.h @@ -48,12 +48,17 @@ namespace xAOD { void set_channel(unsigned int value); /// 32-bit version of channel identifier + /// From Identifier::get_identifier32()::get_compact() unsigned int id() const; void set_id(unsigned int value); + /// Interface for proper identifier Identifier identify() const { return Identifier(this->id()); } + void set_identifier(const Identifier& id) { + set_id(id.get_identifier32().get_compact()); + } /// All values are in units of ns and mV