diff --git a/PhysicsAnalysis/NeutrinoSearch/CMakeLists.txt b/PhysicsAnalysis/NeutrinoSearch/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..33d12b4c20a0c4b7bed9a23dbd5b27f18157ff8a --- /dev/null +++ b/PhysicsAnalysis/NeutrinoSearch/CMakeLists.txt @@ -0,0 +1,12 @@ +atlas_subdir(NeutrinoSearch) + +atlas_add_component( + NeutrinoSearch + src/NeutrinoSearchAlg.h + src/NeutrinoSearchAlg.cxx + src/component/NeutrinoSearch_entries.cxx + LINK_LIBRARIES AthenaBaseComps StoreGateLib xAODFaserWaveform ScintIdentifier FaserCaloIdentifier GeneratorObjects FaserActsGeometryLib TrackerSimEvent TrackerIdentifier TrackerReadoutGeometry TrkTrack GeoPrimitives TrackerRIO_OnTrack +) + +atlas_install_python_modules(python/*.py) +# atlas_install_scripts(test/*.py) diff --git a/PhysicsAnalysis/NeutrinoSearch/python/NeutrinoSearchConfig.py b/PhysicsAnalysis/NeutrinoSearch/python/NeutrinoSearchConfig.py new file mode 100644 index 0000000000000000000000000000000000000000..9f52bf53cae7dadf9a6487522764123c33cef78b --- /dev/null +++ b/PhysicsAnalysis/NeutrinoSearch/python/NeutrinoSearchConfig.py @@ -0,0 +1,104 @@ +""" + Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration +""" + +from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator +from AthenaConfiguration.ComponentFactory import CompFactory +from MagFieldServices.MagFieldServicesConfig import MagneticFieldSvcCfg + +def NeutrinoSearchAlgCfg(flags, **kwargs): + # Initialize GeoModel + from FaserGeoModel.FaserGeoModelConfig import FaserGeometryCfg + acc = FaserGeometryCfg(flags) + + acc.merge(MagneticFieldSvcCfg(flags)) + # acc.merge(FaserActsTrackingGeometrySvcCfg(flags)) + # acc.merge(FaserActsAlignmentCondAlgCfg(flags)) + + actsExtrapolationTool = CompFactory.FaserActsExtrapolationTool("FaserActsExtrapolationTool") + actsExtrapolationTool.MaxSteps = 1000 + actsExtrapolationTool.TrackingGeometryTool = CompFactory.FaserActsTrackingGeometryTool("TrackingGeometryTool") + + NeutrinoSearchAlg = CompFactory.NeutrinoSearchAlg("NeutrinoSearchAlg",**kwargs) + NeutrinoSearchAlg.ExtrapolationTool = actsExtrapolationTool + acc.addEventAlgo(NeutrinoSearchAlg) + + thistSvc = CompFactory.THistSvc() + thistSvc.Output += ["HIST2 DATAFILE='NeutrinoSearch.root' OPT='RECREATE'"] + acc.addService(thistSvc) + + return acc + +if __name__ == "__main__": + + import sys + from AthenaCommon.Logging import log, logging + from AthenaCommon.Constants import DEBUG, VERBOSE, INFO + from AthenaCommon.Configurable import Configurable + from CalypsoConfiguration.AllConfigFlags import ConfigFlags + from AthenaConfiguration.TestDefaults import defaultTestFiles + from CalypsoConfiguration.MainServicesConfig import MainServicesCfg + from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg + # from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg + + # Set up logging and new style config + log.setLevel(DEBUG) + Configurable.configurableRun3Behavior = True + + # Configure + ConfigFlags.Input.Files = [ + '/run/media/dcasper/Data/faser/genie/rec/FaserMC-MDC_Genie_all_150invfb_v1-200001-00000-00007-s0006-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/genie/rec/FaserMC-MDC_Genie_all_150invfb_v1-200001-00008-00015-s0006-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/genie/rec/FaserMC-MDC_Genie_all_150invfb_v1-200001-00016-00023-s0006-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/genie/rec/FaserMC-MDC_Genie_all_150invfb_v1-200001-00024-00031-s0006-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/genie/rec/FaserMC-MDC_Genie_all_150invfb_v1-200001-00032-00039-s0006-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/genie/rec/FaserMC-MDC_Genie_all_150invfb_v1-200001-00040-00047-s0006-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/genie/rec/FaserMC-MDC_Genie_all_150invfb_v1-200001-00048-00055-s0006-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/genie/rec/FaserMC-MDC_Genie_all_150invfb_v1-200001-00056-00063-s0006-r0008-xAOD.root' + ] + ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-02" # Always needed; must match FaserVersionS + ConfigFlags.IOVDb.DatabaseInstance = "OFLP200" # Use MC conditions for now + ConfigFlags.Input.ProjectName = "data21" # Needed to bypass autoconfig + ConfigFlags.Input.isMC = True # Needed to bypass autoconfig + ConfigFlags.GeoModel.FaserVersion = "FASERNU-03" # FASER geometry + ConfigFlags.Common.isOnline = False + ConfigFlags.GeoModel.Align.Dynamic = False + ConfigFlags.Beam.NumberOfCollisions = 0. + + ConfigFlags.Detector.GeometryFaserSCT = True + + ConfigFlags.lock() + + # Core components + acc = MainServicesCfg(ConfigFlags) + acc.merge(PoolReadCfg(ConfigFlags)) + + # algorithm + acc.merge(NeutrinoSearchAlgCfg(ConfigFlags)) + + # # Hack to avoid problem with our use of MC databases when isMC = False + # replicaSvc = acc.getService("DBReplicaSvc") + # replicaSvc.COOLSQLiteVetoPattern = "" + # replicaSvc.UseCOOLSQLite = True + # replicaSvc.UseCOOLFrontier = False + # replicaSvc.UseGeomSQLite = True + + # Timing + #acc.merge(MergeRecoTimingObjCfg(ConfigFlags)) + + # Dump config + # logging.getLogger('forcomps').setLevel(VERBOSE) + # acc.foreach_component("*").OutputLevel = VERBOSE + # acc.foreach_component("*ClassID*").OutputLevel = INFO + # acc.getCondAlgo("FaserSCT_AlignCondAlg").OutputLevel = VERBOSE + # acc.getCondAlgo("FaserSCT_DetectorElementCondAlg").OutputLevel = VERBOSE + # acc.getService("StoreGateSvc").Dump = True + # acc.getService("ConditionStore").Dump = True + # acc.printConfig(withDetails=True) + # ConfigFlags.dump() + + # Execute and finish + sc = acc.run(maxEvents=-1) + + # Success should be 0 + sys.exit(not sc.isSuccess()) diff --git a/PhysicsAnalysis/NeutrinoSearch/src/NeutrinoSearchAlg.cxx b/PhysicsAnalysis/NeutrinoSearch/src/NeutrinoSearchAlg.cxx new file mode 100644 index 0000000000000000000000000000000000000000..cf68a84ec2b0d1786ae0f848b61ab9263a1e1a80 --- /dev/null +++ b/PhysicsAnalysis/NeutrinoSearch/src/NeutrinoSearchAlg.cxx @@ -0,0 +1,304 @@ +#include "NeutrinoSearchAlg.h" +#include "TrackerRIO_OnTrack/FaserSCT_ClusterOnTrack.h" +#include "TrackerIdentifier/FaserSCT_ID.h" +#include "ScintIdentifier/VetoNuID.h" +#include "ScintIdentifier/VetoID.h" +#include "ScintIdentifier/TriggerID.h" +#include "ScintIdentifier/PreshowerID.h" +#include "FaserCaloIdentifier/EcalID.h" +#include "TrackerPrepRawData/FaserSCT_Cluster.h" +#include "Identifier/Identifier.h" +#include "TrackerReadoutGeometry/SCT_DetectorManager.h" +#include "TrackerReadoutGeometry/SiDetectorElement.h" +#include <cmath> + + + +NeutrinoSearchAlg::NeutrinoSearchAlg(const std::string &name, + ISvcLocator *pSvcLocator) + : AthReentrantAlgorithm(name, pSvcLocator), + AthHistogramming(name), + m_histSvc("THistSvc/THistSvc", name) {} + + +StatusCode NeutrinoSearchAlg::initialize() +{ + ATH_CHECK(m_trackCollection.initialize()); + ATH_CHECK(m_vetoNuContainer.initialize()); + ATH_CHECK(m_vetoContainer.initialize()); + ATH_CHECK(m_triggerContainer.initialize()); + ATH_CHECK(m_preshowerContainer.initialize()); + ATH_CHECK(m_ecalContainer.initialize()); + + ATH_CHECK(detStore()->retrieve(m_sctHelper, "FaserSCT_ID")); + ATH_CHECK(detStore()->retrieve(m_vetoNuHelper, "VetoNuID")); + ATH_CHECK(detStore()->retrieve(m_vetoHelper, "VetoID")); + ATH_CHECK(detStore()->retrieve(m_triggerHelper, "TriggerID")); + ATH_CHECK(detStore()->retrieve(m_preshowerHelper, "PreshowerID")); + ATH_CHECK(detStore()->retrieve(m_ecalHelper, "EcalID")); + + ATH_CHECK(detStore()->retrieve(m_detMgr, "SCT")); + ATH_CHECK(m_extrapolationTool.retrieve()); + + m_tree = new TTree("tree", "tree"); + m_tree->Branch("event_number", &m_event_number, "event_number/I"); + m_tree->Branch("VetoNuPmt0", &m_vetoNu0, "vetoNu0/D"); + m_tree->Branch("VetoNuPmt1", &m_vetoNu1, "vetoNu1/D"); + + m_tree->Branch("VetoPmt00", &m_veto00, "veto00/D"); + m_tree->Branch("VetoPmt01", &m_veto01, "veto01/D"); + m_tree->Branch("VetoPmt10", &m_veto10, "veto10/D"); + m_tree->Branch("VetoPmt11", &m_veto11, "veto11/D"); + + m_tree->Branch("TriggerPmt00", &m_trigger00, "trigger00/D"); + m_tree->Branch("TriggerPmt01", &m_trigger01, "trigger01/D"); + m_tree->Branch("TriggerPmt10", &m_trigger10, "trigger10/D"); + m_tree->Branch("TriggerPmt11", &m_trigger11, "trigger11/D"); + + m_tree->Branch("PreshowerPmt0", &m_preshower0, "preshower0/D"); + m_tree->Branch("PreshowerPmt1", &m_preshower1, "preshower1/D"); + + m_tree->Branch("EcalPmt00", &m_ecal00, "ecal00/D"); + m_tree->Branch("EcalPmt01", &m_ecal01, "ecal01/D"); + m_tree->Branch("EcalPmt10", &m_ecal10, "ecal10/D"); + m_tree->Branch("EcalPmt11", &m_ecal11, "ecal11/D"); + + m_tree->Branch("px", &m_px, "px/D"); + m_tree->Branch("py", &m_py, "py/D"); + m_tree->Branch("pz", &m_pz, "pz/D"); + m_tree->Branch("pe", &m_pe, "pe/D"); + m_tree->Branch("chi2", &m_chi2, "chi2/D"); + ATH_CHECK(histSvc()->regTree("/HIST2/tree", m_tree)); + + return StatusCode::SUCCESS; +} + + +StatusCode NeutrinoSearchAlg::execute(const EventContext &ctx) const +{ + m_event_number = ctx.eventID().event_number(); + + // Waveform data + m_vetoNu0 = 0; + m_vetoNu1 = 0; + m_veto00 = 0; + m_veto01 = 0; + m_veto10 = 0; + m_veto11 = 0; + m_trigger00 = 0; + m_trigger01 = 0; + m_trigger10 = 0; + m_trigger11 = 0; + m_preshower0 = 0; + m_preshower1 = 0; + m_ecal00 = 0; + m_ecal01 = 0; + m_ecal10 = 0; + m_ecal11 = 0; + + SG::ReadHandle<xAOD::WaveformHitContainer> vetoNuContainer { m_vetoNuContainer, ctx }; + ATH_CHECK(vetoNuContainer.isValid()); + + SG::ReadHandle<xAOD::WaveformHitContainer> vetoContainer { m_vetoContainer, ctx }; + ATH_CHECK(vetoContainer.isValid()); + + SG::ReadHandle<xAOD::WaveformHitContainer> triggerContainer { m_triggerContainer, ctx }; + ATH_CHECK(triggerContainer.isValid()); + + SG::ReadHandle<xAOD::WaveformHitContainer> preshowerContainer { m_preshowerContainer, ctx }; + ATH_CHECK(preshowerContainer.isValid()); + + SG::ReadHandle<xAOD::WaveformHitContainer> ecalContainer { m_ecalContainer, ctx }; + ATH_CHECK(ecalContainer.isValid()); + + for (auto hit : *vetoNuContainer) + { + auto id = hit->identify(); + if (m_vetoNuHelper->plate(id) == 0) + { + m_vetoNu0 += hit->integral(); + } + else if (m_vetoNuHelper->plate(id) == 1) + { + m_vetoNu1 += hit->integral(); + } + else + { + ATH_MSG_FATAL("Invalid VetoNu plate number: " << m_vetoNuHelper->plate(id)); + return StatusCode::FAILURE; + } + } + + for (auto hit : *vetoContainer) + { + auto id = hit->identify(); + auto station = m_vetoHelper->station(id); + auto plate = m_vetoHelper->plate(id); + if (station == 0) + { + if (plate == 0) + { + m_veto00 += hit->integral(); + } + else if (plate == 1) + { + m_veto01 += hit->integral(); + } + else + { + ATH_MSG_FATAL("Invalid Veto plate number: " << plate); + } + } + else if (station == 1) + { + if (plate == 0) + { + m_veto10 += hit->integral(); + } + else if (plate == 1) + { + m_veto11 += hit->integral(); + } + else + { + ATH_MSG_FATAL("Invalid Veto plate number: " << plate); + } + } + else + { + ATH_MSG_FATAL("Invalid Veto station number: " << station); + return StatusCode::FAILURE; + } + } + + for (auto hit : *triggerContainer) + { + auto id = hit->identify(); + auto plate = m_triggerHelper->plate(id); + auto pmt = m_triggerHelper->pmt(id); + if (plate == 0) + { + if (pmt == 0) + { + m_trigger00 += hit->integral(); + } + else if (pmt == 1) + { + m_trigger01 += hit->integral(); + } + else + { + ATH_MSG_FATAL("Invalid Trigger pmt number: " << pmt); + } + } + else if (plate == 1) + { + if (pmt == 0) + { + m_trigger10 += hit->integral(); + } + else if (pmt == 1) + { + m_trigger11 += hit->integral(); + } + else + { + ATH_MSG_FATAL("Invalid Trigger pmt number: " << pmt); + } + } + else + { + ATH_MSG_FATAL("Invalid Trigger plate number: " << plate); + return StatusCode::FAILURE; + } + } + + for (auto hit : *preshowerContainer) + { + auto id = hit->identify(); + if (m_preshowerHelper->plate(id) == 0) + { + m_preshower0 += hit->integral(); + } + else if (m_preshowerHelper->plate(id) == 1) + { + m_preshower1 += hit->integral(); + } + else + { + ATH_MSG_FATAL("Invalid Preshower plate number: " << m_preshowerHelper->plate(id)); + return StatusCode::FAILURE; + } + } + + for (auto hit : *ecalContainer) + { + auto id = hit->identify(); + auto row = m_ecalHelper->row(id); + auto mod = m_ecalHelper->module(id); + if (row == 0) + { + if (mod == 0) + { + m_ecal00 += hit->integral(); + } + else if (mod == 1) + { + m_ecal01 += hit->integral(); + } + else + { + ATH_MSG_FATAL("Invalid Ecal module number: " << mod); + } + } + else if (row == 1) + { + if (mod == 0) + { + m_ecal10 += hit->integral(); + } + else if (mod == 1) + { + m_ecal11 += hit->integral(); + } + else + { + ATH_MSG_FATAL("Invalid Ecal module number: " << mod); + } + } + else + { + ATH_MSG_FATAL("Invalid Ecal row number: " << row); + return StatusCode::FAILURE; + } + } + + SG::ReadHandle<TrackCollection> trackCollection {m_trackCollection, ctx}; + ATH_CHECK(trackCollection.isValid()); + + const Trk::Track* trackCandidate {nullptr}; + + for (auto track : *trackCollection) + { + if (track == nullptr) continue; + } + + if (m_vetoNu0 > m_vetoNuThreshold || m_vetoNu1 > m_vetoNuThreshold || + m_veto00 > m_vetoThreshold || m_veto01 > m_vetoThreshold || m_veto10 > m_vetoThreshold || m_veto11 > m_vetoThreshold || + m_trigger00 > m_triggerThreshold || m_trigger01 > m_triggerThreshold || m_trigger10 > m_triggerThreshold || m_trigger11 > m_triggerThreshold || + m_preshower0 > m_preshowerThreshold || m_preshower1 > m_preshowerThreshold || + m_ecal00 > m_ecalThreshold || m_ecal01 > m_ecalThreshold || m_ecal10 > m_ecalThreshold || m_ecal11 > m_ecalThreshold) + { + m_tree->Fill(); + } + + return StatusCode::SUCCESS; +} + + +StatusCode NeutrinoSearchAlg::finalize() +{ + return StatusCode::SUCCESS; +} + + diff --git a/PhysicsAnalysis/NeutrinoSearch/src/NeutrinoSearchAlg.h b/PhysicsAnalysis/NeutrinoSearch/src/NeutrinoSearchAlg.h new file mode 100644 index 0000000000000000000000000000000000000000..6f47a39494ed51da7f9ab656a5cecc9bd7a0028a --- /dev/null +++ b/PhysicsAnalysis/NeutrinoSearch/src/NeutrinoSearchAlg.h @@ -0,0 +1,89 @@ +#ifndef NEUTRINOSEARCH_NEUTRINOSEARCHALG_H +#define NEUTRINOSEARCH_NEUTRINOSEARCHALG_H + +#include "AthenaBaseComps/AthReentrantAlgorithm.h" +#include "AthenaBaseComps/AthHistogramming.h" +#include "TrkTrack/TrackCollection.h" +#include "xAODFaserWaveform/WaveformHitContainer.h" +#include "FaserActsGeometryInterfaces/IFaserActsExtrapolationTool.h" + + +class TTree; +class FaserSCT_ID; +class VetoNuID; +class VetoID; +class TriggerID; +class PreshowerID; +class EcalID; +namespace TrackerDD +{ + class SCT_DetectorManager; +} + +class NeutrinoSearchAlg : public AthReentrantAlgorithm, AthHistogramming { +public: + NeutrinoSearchAlg(const std::string &name, ISvcLocator *pSvcLocator); + virtual ~NeutrinoSearchAlg() = default; + virtual StatusCode initialize() override; + virtual StatusCode execute(const EventContext &ctx) const override; + virtual StatusCode finalize() override; + const ServiceHandle <ITHistSvc> &histSvc() const; + +private: + + ServiceHandle <ITHistSvc> m_histSvc; + SG::ReadHandleKey<TrackCollection> m_trackCollection { this, "TrackCollection", "CKFTrackCollection", "Input track collection name" }; + SG::ReadHandleKey<xAOD::WaveformHitContainer> m_vetoNuContainer { this, "VetoNuContainer", "VetoNuWaveformHits", "VetoNu hit container name" }; + SG::ReadHandleKey<xAOD::WaveformHitContainer> m_vetoContainer { this, "VetoContainer", "VetoWaveformHits", "Veto hit container name" }; + SG::ReadHandleKey<xAOD::WaveformHitContainer> m_triggerContainer { this, "TriggerContainer", "TriggerWaveformHits", "Trigger hit container name" }; + SG::ReadHandleKey<xAOD::WaveformHitContainer> m_preshowerContainer { this, "PreshowerContainer", "PreshowerWaveformHits", "Preshower hit container name" }; + SG::ReadHandleKey<xAOD::WaveformHitContainer> m_ecalContainer { this, "EcalContainer", "CaloWaveformHits", "Ecal hit container name" }; + + ToolHandle<IFaserActsExtrapolationTool> m_extrapolationTool { this, "ExtrapolationTool", "FaserActsExtrapolationTool" }; + const TrackerDD::SCT_DetectorManager* m_detMgr {nullptr}; + + const FaserSCT_ID* m_sctHelper; + const VetoNuID* m_vetoNuHelper; + const VetoID* m_vetoHelper; + const TriggerID* m_triggerHelper; + const PreshowerID* m_preshowerHelper; + const EcalID* m_ecalHelper; + + // TODO: use realistic thresholds for MIP + DoubleProperty m_vetoNuThreshold { this, "VetoNuThreshold", 0, "Threshold for VetoNu pmts" }; + DoubleProperty m_vetoThreshold { this, "VetoThreshold", 0, "Threshold for Veto pmts" }; + DoubleProperty m_triggerThreshold { this, "TriggerThreshold", 0, "Threshold for Trigger pmts" }; + DoubleProperty m_preshowerThreshold { this, "PreshowerThreshold", 0, "Threshold for Preshower pmts" }; + DoubleProperty m_ecalThreshold { this, "EcalThreshold", 0, "Threshold for Ecal pmts" }; + + mutable TTree* m_tree; + mutable unsigned int m_event_number; + mutable double m_vetoNu0; + mutable double m_vetoNu1; + mutable double m_veto00; + mutable double m_veto01; + mutable double m_veto10; + mutable double m_veto11; + mutable double m_trigger00; + mutable double m_trigger01; + mutable double m_trigger10; + mutable double m_trigger11; + mutable double m_preshower0; + mutable double m_preshower1; + mutable double m_ecal00; + mutable double m_ecal01; + mutable double m_ecal10; + mutable double m_ecal11; + + mutable double m_px; + mutable double m_py; + mutable double m_pz; + mutable double m_pe; + mutable double m_chi2; +}; + +inline const ServiceHandle <ITHistSvc> &NeutrinoSearchAlg::histSvc() const { + return m_histSvc; +} + +#endif // NEUTRINOSEARCH_NEUTRINOSEARCHALG_H diff --git a/PhysicsAnalysis/NeutrinoSearch/src/component/NeutrinoSearch_entries.cxx b/PhysicsAnalysis/NeutrinoSearch/src/component/NeutrinoSearch_entries.cxx new file mode 100644 index 0000000000000000000000000000000000000000..d47072b60eb09ca5e38a86763f85ff0c5bc36e94 --- /dev/null +++ b/PhysicsAnalysis/NeutrinoSearch/src/component/NeutrinoSearch_entries.cxx @@ -0,0 +1,3 @@ +#include "../NeutrinoSearchAlg.h" + +DECLARE_COMPONENT(NeutrinoSearchAlg) \ No newline at end of file