From 1e94cfc21205fa97ab6b7f0b4d4424f4058b5f72 Mon Sep 17 00:00:00 2001 From: Dave Casper <dcasper@uci.edu> Date: Thu, 7 Jul 2022 00:05:03 -0700 Subject: [PATCH] Added information to tree [skip ci] --- PhysicsAnalysis/NeutrinoSearch/CMakeLists.txt | 2 +- .../python/FlukaSearchConfig.py | 140 +++++++ ...noSearchConfig.py => GenieSearchConfig.py} | 4 +- .../NeutrinoSearch/src/NeutrinoSearchAlg.cxx | 362 +++++++++++++++++- .../NeutrinoSearch/src/NeutrinoSearchAlg.h | 55 ++- 5 files changed, 533 insertions(+), 30 deletions(-) create mode 100644 PhysicsAnalysis/NeutrinoSearch/python/FlukaSearchConfig.py rename PhysicsAnalysis/NeutrinoSearch/python/{NeutrinoSearchConfig.py => GenieSearchConfig.py} (96%) diff --git a/PhysicsAnalysis/NeutrinoSearch/CMakeLists.txt b/PhysicsAnalysis/NeutrinoSearch/CMakeLists.txt index 33d12b4c..cff2ab85 100644 --- a/PhysicsAnalysis/NeutrinoSearch/CMakeLists.txt +++ b/PhysicsAnalysis/NeutrinoSearch/CMakeLists.txt @@ -5,7 +5,7 @@ atlas_add_component( 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 + LINK_LIBRARIES AthenaBaseComps StoreGateLib xAODFaserWaveform ScintIdentifier FaserCaloIdentifier GeneratorObjects FaserActsGeometryLib TrackerSimEvent TrackerSimData TrackerIdentifier TrackerReadoutGeometry TrkTrack GeoPrimitives TrackerRIO_OnTrack ) atlas_install_python_modules(python/*.py) diff --git a/PhysicsAnalysis/NeutrinoSearch/python/FlukaSearchConfig.py b/PhysicsAnalysis/NeutrinoSearch/python/FlukaSearchConfig.py new file mode 100644 index 00000000..4e196a37 --- /dev/null +++ b/PhysicsAnalysis/NeutrinoSearch/python/FlukaSearchConfig.py @@ -0,0 +1,140 @@ +""" + 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='FlukaSearch.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/fluka/210001/rec/r0008/FaserMC-MDC_Fluka_unit30_Nm_71m_m3750_v3-210001-00000-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210001/rec/r0008/FaserMC-MDC_Fluka_unit30_Nm_71m_m3750_v3-210001-00001-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210001/rec/r0008/FaserMC-MDC_Fluka_unit30_Nm_71m_m3750_v3-210001-00002-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210001/rec/r0008/FaserMC-MDC_Fluka_unit30_Nm_71m_m3750_v3-210001-00003-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210001/rec/r0008/FaserMC-MDC_Fluka_unit30_Nm_71m_m3750_v3-210001-00004-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210001/rec/r0008/FaserMC-MDC_Fluka_unit30_Nm_71m_m3750_v3-210001-00005-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210001/rec/r0008/FaserMC-MDC_Fluka_unit30_Nm_71m_m3750_v3-210001-00006-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210001/rec/r0008/FaserMC-MDC_Fluka_unit30_Nm_71m_m3750_v3-210001-00007-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210001/rec/r0008/FaserMC-MDC_Fluka_unit30_Nm_71m_m3750_v3-210001-00008-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210001/rec/r0008/FaserMC-MDC_Fluka_unit30_Nm_71m_m3750_v3-210001-00009-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210001/rec/r0008/FaserMC-MDC_Fluka_unit30_Nm_71m_m3750_v3-210001-00010-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210001/rec/r0008/FaserMC-MDC_Fluka_unit30_Nm_71m_m3750_v3-210001-00011-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210001/rec/r0008/FaserMC-MDC_Fluka_unit30_Nm_71m_m3750_v3-210001-00012-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210001/rec/r0008/FaserMC-MDC_Fluka_unit30_Nm_71m_m3750_v3-210001-00013-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210001/rec/r0008/FaserMC-MDC_Fluka_unit30_Nm_71m_m3750_v3-210001-00014-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210001/rec/r0008/FaserMC-MDC_Fluka_unit30_Nm_71m_m3750_v3-210001-00015-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210001/rec/r0008/FaserMC-MDC_Fluka_unit30_Nm_71m_m3750_v3-210001-00016-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210001/rec/r0008/FaserMC-MDC_Fluka_unit30_Nm_71m_m3750_v3-210001-00017-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210001/rec/r0008/FaserMC-MDC_Fluka_unit30_Nm_71m_m3750_v3-210001-00018-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210001/rec/r0008/FaserMC-MDC_Fluka_unit30_Nm_71m_m3750_v3-210001-00019-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210001/rec/r0008/FaserMC-MDC_Fluka_unit30_Nm_71m_m3750_v3-210001-00020-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210001/rec/r0008/FaserMC-MDC_Fluka_unit30_Nm_71m_m3750_v3-210001-00021-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210001/rec/r0008/FaserMC-MDC_Fluka_unit30_Nm_71m_m3750_v3-210001-00022-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210001/rec/r0008/FaserMC-MDC_Fluka_unit30_Nm_71m_m3750_v3-210001-00023-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210001/rec/r0008/FaserMC-MDC_Fluka_unit30_Nm_71m_m3750_v3-210001-00024-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210001/rec/r0008/FaserMC-MDC_Fluka_unit30_Nm_71m_m3750_v3-210001-00025-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210001/rec/r0008/FaserMC-MDC_Fluka_unit30_Nm_71m_m3750_v3-210001-00026-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210001/rec/r0008/FaserMC-MDC_Fluka_unit30_Nm_71m_m3750_v3-210001-00027-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210002/rec/r0008/FaserMC-MDC_Fluka_unit30_Pm_71m_m3750_v3-210002-00000-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210002/rec/r0008/FaserMC-MDC_Fluka_unit30_Pm_71m_m3750_v3-210002-00001-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210002/rec/r0008/FaserMC-MDC_Fluka_unit30_Pm_71m_m3750_v3-210002-00002-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210002/rec/r0008/FaserMC-MDC_Fluka_unit30_Pm_71m_m3750_v3-210002-00003-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210002/rec/r0008/FaserMC-MDC_Fluka_unit30_Pm_71m_m3750_v3-210002-00004-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210002/rec/r0008/FaserMC-MDC_Fluka_unit30_Pm_71m_m3750_v3-210002-00005-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210002/rec/r0008/FaserMC-MDC_Fluka_unit30_Pm_71m_m3750_v3-210002-00006-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210002/rec/r0008/FaserMC-MDC_Fluka_unit30_Pm_71m_m3750_v3-210002-00007-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210002/rec/r0008/FaserMC-MDC_Fluka_unit30_Pm_71m_m3750_v3-210002-00008-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210002/rec/r0008/FaserMC-MDC_Fluka_unit30_Pm_71m_m3750_v3-210002-00009-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210002/rec/r0008/FaserMC-MDC_Fluka_unit30_Pm_71m_m3750_v3-210002-00010-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210002/rec/r0008/FaserMC-MDC_Fluka_unit30_Pm_71m_m3750_v3-210002-00011-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210002/rec/r0008/FaserMC-MDC_Fluka_unit30_Pm_71m_m3750_v3-210002-00012-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210002/rec/r0008/FaserMC-MDC_Fluka_unit30_Pm_71m_m3750_v3-210002-00013-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210002/rec/r0008/FaserMC-MDC_Fluka_unit30_Pm_71m_m3750_v3-210002-00014-s0005-r0008-xAOD.root', + '/run/media/dcasper/Data/faser/fluka/210002/rec/r0008/FaserMC-MDC_Fluka_unit30_Pm_71m_m3750_v3-210002-00015-s0005-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, UseFlukaWeights=True)) + + # # 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/python/NeutrinoSearchConfig.py b/PhysicsAnalysis/NeutrinoSearch/python/GenieSearchConfig.py similarity index 96% rename from PhysicsAnalysis/NeutrinoSearch/python/NeutrinoSearchConfig.py rename to PhysicsAnalysis/NeutrinoSearch/python/GenieSearchConfig.py index 9f52bf53..a80de32b 100644 --- a/PhysicsAnalysis/NeutrinoSearch/python/NeutrinoSearchConfig.py +++ b/PhysicsAnalysis/NeutrinoSearch/python/GenieSearchConfig.py @@ -24,7 +24,7 @@ def NeutrinoSearchAlgCfg(flags, **kwargs): acc.addEventAlgo(NeutrinoSearchAlg) thistSvc = CompFactory.THistSvc() - thistSvc.Output += ["HIST2 DATAFILE='NeutrinoSearch.root' OPT='RECREATE'"] + thistSvc.Output += ["HIST2 DATAFILE='GenieSearch.root' OPT='RECREATE'"] acc.addService(thistSvc) return acc @@ -74,7 +74,7 @@ if __name__ == "__main__": acc.merge(PoolReadCfg(ConfigFlags)) # algorithm - acc.merge(NeutrinoSearchAlgCfg(ConfigFlags)) + acc.merge(NeutrinoSearchAlgCfg(ConfigFlags, UseGenieWeights=True)) # # Hack to avoid problem with our use of MC databases when isMC = False # replicaSvc = acc.getService("DBReplicaSvc") diff --git a/PhysicsAnalysis/NeutrinoSearch/src/NeutrinoSearchAlg.cxx b/PhysicsAnalysis/NeutrinoSearch/src/NeutrinoSearchAlg.cxx index cf68a84e..8f4525c4 100644 --- a/PhysicsAnalysis/NeutrinoSearch/src/NeutrinoSearchAlg.cxx +++ b/PhysicsAnalysis/NeutrinoSearch/src/NeutrinoSearchAlg.cxx @@ -1,4 +1,5 @@ #include "NeutrinoSearchAlg.h" +#include "TrkTrack/Track.h" #include "TrackerRIO_OnTrack/FaserSCT_ClusterOnTrack.h" #include "TrackerIdentifier/FaserSCT_ID.h" #include "ScintIdentifier/VetoNuID.h" @@ -10,6 +11,8 @@ #include "Identifier/Identifier.h" #include "TrackerReadoutGeometry/SCT_DetectorManager.h" #include "TrackerReadoutGeometry/SiDetectorElement.h" +#include "TrackerPrepRawData/FaserSCT_Cluster.h" +#include "xAODTruth/TruthParticle.h" #include <cmath> @@ -23,12 +26,16 @@ NeutrinoSearchAlg::NeutrinoSearchAlg(const std::string &name, StatusCode NeutrinoSearchAlg::initialize() { + ATH_CHECK(m_truthEventContainer.initialize()); + ATH_CHECK(m_truthParticleContainer.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(m_clusterContainer.initialize()); + ATH_CHECK(m_simDataCollection.initialize()); ATH_CHECK(detStore()->retrieve(m_sctHelper, "FaserSCT_ID")); ATH_CHECK(detStore()->retrieve(m_vetoNuHelper, "VetoNuID")); @@ -40,20 +47,37 @@ StatusCode NeutrinoSearchAlg::initialize() ATH_CHECK(detStore()->retrieve(m_detMgr, "SCT")); ATH_CHECK(m_extrapolationTool.retrieve()); + if (m_useFlukaWeights) + { + m_baseEventCrossSection = (m_flukaCrossSection * kfemtoBarnsPerMilliBarn)/m_flukaCollisions; + } + else if (m_useGenieWeights) + { + m_baseEventCrossSection = 1.0/m_genieLuminosity; + } + else + { + m_baseEventCrossSection = 1.0; + } + 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("VetoUpstream", &m_vetoUpstream, "vetoUpstream/D"); m_tree->Branch("VetoPmt10", &m_veto10, "veto10/D"); m_tree->Branch("VetoPmt11", &m_veto11, "veto11/D"); + m_tree->Branch("VetoDownstream", &m_vetoDownstream, "vetoDownstream/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("TriggerTotal", &m_triggerTotal, "triggerTotal/D"); m_tree->Branch("PreshowerPmt0", &m_preshower0, "preshower0/D"); m_tree->Branch("PreshowerPmt1", &m_preshower1, "preshower1/D"); @@ -62,14 +86,39 @@ StatusCode NeutrinoSearchAlg::initialize() 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("EcalTotal", &m_ecalTotal, "ecalTotal/D"); + m_tree->Branch("Clust0", &m_station0Clusters, "clust0/I"); + m_tree->Branch("Clust1", &m_station1Clusters, "clust0/I"); + m_tree->Branch("Clust2", &m_station2Clusters, "clust0/I"); + m_tree->Branch("Clust3", &m_station3Clusters, "clust0/I"); + + m_tree->Branch("x", &m_x, "x/D"); + m_tree->Branch("y", &m_y, "y/D"); + m_tree->Branch("z", &m_z, "z/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("charge", &m_charge, "charge/I"); m_tree->Branch("chi2", &m_chi2, "chi2/D"); + m_tree->Branch("ndof", &m_ndof, "ndof/I"); + m_tree->Branch("longTracks", &m_longTracks, "longTracks/I"); + m_tree->Branch("pTruthLepton", &m_truthLeptonMomentum, "pTruthLepton/D"); + m_tree->Branch("truthBarcode", &m_truthBarcode, "truthBarcode/I"); + m_tree->Branch("truthPdg", &m_truthPdg, "truthPdg/I"); + m_tree->Branch("CrossSection", &m_crossSection, "crossSection/D"); + ATH_CHECK(histSvc()->regTree("/HIST2/tree", m_tree)); + if (m_enforceBlinding) + { + ATH_MSG_INFO("Blinding will be enforced for real data."); + } + else + { + ATH_MSG_INFO("Blinding will NOT be enforced for real data."); + } + return StatusCode::SUCCESS; } @@ -85,34 +134,89 @@ StatusCode NeutrinoSearchAlg::execute(const EventContext &ctx) const m_veto01 = 0; m_veto10 = 0; m_veto11 = 0; + m_vetoUpstream = 0; + m_vetoDownstream = 0; m_trigger00 = 0; m_trigger01 = 0; m_trigger10 = 0; m_trigger11 = 0; + m_triggerTotal = 0; m_preshower0 = 0; m_preshower1 = 0; m_ecal00 = 0; m_ecal01 = 0; m_ecal10 = 0; m_ecal11 = 0; + m_ecalTotal = 0; + m_station0Clusters = 0; + m_station1Clusters = 0; + m_station2Clusters = 0; + m_station3Clusters = 0; + m_crossSection = 0; + m_chi2 = 0; + m_ndof = 0; + m_px = 0; + m_py = 0; + m_pz = 0; + m_charge = 0; + m_x = 0; + m_y = 0; + m_z = 0; + m_longTracks = 0; + m_truthLeptonMomentum = 0; + m_truthBarcode = 0; + m_truthPdg = 0; + + bool realData = true; + + // Work out effective cross section for MC + SG::ReadHandle<xAOD::TruthEventContainer> truthEventContainer { m_truthEventContainer, ctx }; + if (truthEventContainer.isValid() && truthEventContainer->size() > 0) + { + realData = false; + if (m_useFlukaWeights) + { + double flukaWeight = truthEventContainer->at(0)->weights()[0]; + ATH_MSG_ALWAYS("Found fluka weight = " << flukaWeight); + m_crossSection = m_baseEventCrossSection * flukaWeight; + } + else if (m_useGenieWeights) + { + m_crossSection = m_baseEventCrossSection; + } + else + { + ATH_MSG_WARNING("Monte carlo event with no weighting scheme specified. Setting crossSection (weight) to " << m_baseEventCrossSection << " fb."); + m_crossSection = m_baseEventCrossSection; + } + } - SG::ReadHandle<xAOD::WaveformHitContainer> vetoNuContainer { m_vetoNuContainer, ctx }; - ATH_CHECK(vetoNuContainer.isValid()); - - SG::ReadHandle<xAOD::WaveformHitContainer> vetoContainer { m_vetoContainer, ctx }; - ATH_CHECK(vetoContainer.isValid()); + // Find the primary lepton (if any) + SG::ReadHandle<xAOD::TruthParticleContainer> truthParticleContainer { m_truthParticleContainer, ctx }; + if (truthParticleContainer.isValid() && truthParticleContainer->size() > 0) + { + for (auto particle : *truthParticleContainer) + { + if ( particle->absPdgId() == 11 || particle->absPdgId() == 13 || particle->absPdgId() == 15 ) + { + if (particle->status() == 1 && (particle->nParents() == 0 || particle->nParents() == 2) ) + m_truthLeptonMomentum = particle->p4().P(); - SG::ReadHandle<xAOD::WaveformHitContainer> triggerContainer { m_triggerContainer, ctx }; - ATH_CHECK(triggerContainer.isValid()); + break; + } + } + } - SG::ReadHandle<xAOD::WaveformHitContainer> preshowerContainer { m_preshowerContainer, ctx }; - ATH_CHECK(preshowerContainer.isValid()); + SG::ReadHandle<xAOD::WaveformHitContainer> vetoNuContainer { m_vetoNuContainer, ctx }; + ATH_CHECK(vetoNuContainer.isValid()); - SG::ReadHandle<xAOD::WaveformHitContainer> ecalContainer { m_ecalContainer, ctx }; - ATH_CHECK(ecalContainer.isValid()); + // If real data, check for blinding before we do anything else + bool blinded = realData; for (auto hit : *vetoNuContainer) { + if (!waveformHitOK(hit)) continue; + blinded = false; auto id = hit->identify(); if (m_vetoNuHelper->plate(id) == 0) { @@ -129,8 +233,23 @@ StatusCode NeutrinoSearchAlg::execute(const EventContext &ctx) const } } + if (m_enforceBlinding && blinded) return StatusCode::SUCCESS; + + 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 : *vetoContainer) { + if (!waveformHitOK(hit)) continue; auto id = hit->identify(); auto station = m_vetoHelper->station(id); auto plate = m_vetoHelper->plate(id); @@ -139,10 +258,12 @@ StatusCode NeutrinoSearchAlg::execute(const EventContext &ctx) const if (plate == 0) { m_veto00 += hit->integral(); + m_vetoUpstream += hit->integral(); } else if (plate == 1) { m_veto01 += hit->integral(); + m_vetoUpstream += hit->integral(); } else { @@ -154,10 +275,12 @@ StatusCode NeutrinoSearchAlg::execute(const EventContext &ctx) const if (plate == 0) { m_veto10 += hit->integral(); + m_vetoDownstream += hit->integral(); } else if (plate == 1) { m_veto11 += hit->integral(); + m_vetoDownstream += hit->integral(); } else { @@ -173,6 +296,7 @@ StatusCode NeutrinoSearchAlg::execute(const EventContext &ctx) const for (auto hit : *triggerContainer) { + if (!waveformHitOK(hit)) continue; auto id = hit->identify(); auto plate = m_triggerHelper->plate(id); auto pmt = m_triggerHelper->pmt(id); @@ -181,10 +305,12 @@ StatusCode NeutrinoSearchAlg::execute(const EventContext &ctx) const if (pmt == 0) { m_trigger00 += hit->integral(); + m_triggerTotal += hit->integral(); } else if (pmt == 1) { m_trigger01 += hit->integral(); + m_triggerTotal += hit->integral(); } else { @@ -196,10 +322,12 @@ StatusCode NeutrinoSearchAlg::execute(const EventContext &ctx) const if (pmt == 0) { m_trigger10 += hit->integral(); + m_triggerTotal += hit->integral(); } else if (pmt == 1) { m_trigger11 += hit->integral(); + m_triggerTotal += hit->integral(); } else { @@ -215,6 +343,7 @@ StatusCode NeutrinoSearchAlg::execute(const EventContext &ctx) const for (auto hit : *preshowerContainer) { + if (!waveformHitOK(hit)) continue; auto id = hit->identify(); if (m_preshowerHelper->plate(id) == 0) { @@ -233,6 +362,7 @@ StatusCode NeutrinoSearchAlg::execute(const EventContext &ctx) const for (auto hit : *ecalContainer) { + if (!waveformHitOK(hit)) continue; auto id = hit->identify(); auto row = m_ecalHelper->row(id); auto mod = m_ecalHelper->module(id); @@ -241,10 +371,12 @@ StatusCode NeutrinoSearchAlg::execute(const EventContext &ctx) const if (mod == 0) { m_ecal00 += hit->integral(); + m_ecalTotal += hit->integral(); } else if (mod == 1) { m_ecal01 += hit->integral(); + m_ecalTotal += hit->integral(); } else { @@ -256,10 +388,12 @@ StatusCode NeutrinoSearchAlg::execute(const EventContext &ctx) const if (mod == 0) { m_ecal10 += hit->integral(); + m_ecalTotal += hit->integral(); } else if (mod == 1) { m_ecal11 += hit->integral(); + m_ecalTotal += hit->integral(); } else { @@ -273,25 +407,207 @@ StatusCode NeutrinoSearchAlg::execute(const EventContext &ctx) const } } + SG::ReadHandle<Tracker::FaserSCT_ClusterContainer> clusterContainer { m_clusterContainer, ctx }; + ATH_CHECK(clusterContainer.isValid()); + + for (auto collection : *clusterContainer) + { + Identifier id = collection->identify(); + int station = m_sctHelper->station(id); + int clusters = (int) collection->size(); + switch (station) + { + case 0: + m_station0Clusters += clusters; + break; + case 1: + m_station1Clusters += clusters; + break; + case 2: + m_station2Clusters += clusters; + break; + case 3: + m_station3Clusters += clusters; + break; + default: + ATH_MSG_FATAL("Unknown tracker station number " << station); + break; + } + } + + SG::ReadHandle<TrackCollection> trackCollection {m_trackCollection, ctx}; ATH_CHECK(trackCollection.isValid()); - const Trk::Track* trackCandidate {nullptr}; + const Trk::TrackParameters* candidateParameters {nullptr}; + const Trk::Track* candidateTrack {nullptr}; - for (auto track : *trackCollection) + for (const Trk::Track* track : *trackCollection) { if (track == nullptr) continue; + std::set<int> stationMap; + std::set<std::pair<int, int>> layerMap; + + // Check for hit in the three downstream stations + for (auto measurement : *(track->measurementsOnTrack())) + { + const Tracker::FaserSCT_ClusterOnTrack* cluster = dynamic_cast<const Tracker::FaserSCT_ClusterOnTrack*>(measurement); + if (cluster != nullptr) + { + Identifier id = cluster->identify(); + int station = m_sctHelper->station(id); + int layer = m_sctHelper->layer(id); + stationMap.emplace(station); + layerMap.emplace(station, layer); + } + } + if (stationMap.count(1) == 0 || stationMap.count(2) == 0 || stationMap.count(3) == 0) continue; + + int nLayers = std::count_if(layerMap.begin(), layerMap.end(), [](std::pair<int,int> p){return p.first != 0;}); + if (nLayers < m_minTrackerLayers) continue; + m_longTracks++; + const Trk::TrackParameters* upstreamParameters {nullptr}; + for (auto params : *(track->trackParameters())) + { + if (upstreamParameters == nullptr || params->position().z() < upstreamParameters->position().z()) upstreamParameters = params; + } + if (candidateParameters == nullptr || upstreamParameters->momentum().mag() > candidateParameters->momentum().mag()) + { + candidateParameters = upstreamParameters; + candidateTrack = track; + m_chi2 = track->fitQuality()->chiSquared(); + m_ndof = track->fitQuality()->numberDoF(); + } } - 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) + SG::ReadHandle<TrackerSimDataCollection> simDataCollection {m_simDataCollection, ctx}; +// ATH_MSG_INFO("SimData valid? " << simDataCollection.isValid()); + if (candidateTrack != nullptr && simDataCollection.isValid()) { - m_tree->Fill(); + std::map<int, float> truthMap; + for (auto measurement : *(candidateTrack->measurementsOnTrack())) + { + const Tracker::FaserSCT_ClusterOnTrack* cluster = dynamic_cast<const Tracker::FaserSCT_ClusterOnTrack*>(measurement); + if (cluster != nullptr) + { + // ATH_MSG_INFO("ClusterOnTrack is OK"); + cluster->dump(msg()); + +// Hack to work around issue with cluster navigation + + auto idRDO = cluster->identify(); + + if (simDataCollection->count(idRDO) > 0) + { + // ATH_MSG_INFO("rdo entry found"); + const auto& simdata = simDataCollection->find(idRDO)->second; + const auto& deposits = simdata.getdeposits(); + //loop through deposits and record contributions + HepMcParticleLink primary{}; + for( const auto& depositPair : deposits) + { + // ATH_MSG_INFO("Deposit found"); + float eDep = depositPair.second; + int barcode = depositPair.first->barcode(); + // if( depositPair.second > highestDep) + // { + // highestDep = depositPair.second; + // barcode = depositPair.first->barcode(); + // primary = depositPair.first; + // depositPair.first->print(std::cout); + // ATH_MSG_INFO("pdg id "<<depositPair.first->pdg_id()); + // } + if (truthMap.count(barcode) > 0) + { + truthMap[barcode] += eDep; + } + else + { + truthMap[barcode] = eDep; + } + } + } + + + + + // // const Tracker::FaserSCT_Cluster* origCluster = dynamic_cast<const Tracker::FaserSCT_Cluster*>(cluster->prepRawData()); + // auto origCluster = cluster->prepRawData(); + // if (origCluster != nullptr) + // { + // ATH_MSG_INFO("Orig Cluster is OK"); + // auto rdoList = origCluster->rdoList(); + // for (auto idRDO : rdoList) + // { + // ATH_MSG_INFO("rdoList not empty"); + // if (simDataCollection->count(idRDO) > 0) + // { + // ATH_MSG_INFO("rdo entry found"); + // const auto& simdata = simDataCollection->find(idRDO)->second; + // const auto& deposits = simdata.getdeposits(); + // //loop through deposits and record contributions + // HepMcParticleLink primary{}; + // for( const auto& depositPair : deposits) + // { + // ATH_MSG_INFO("Deposit found"); + // float eDep = depositPair.second; + // int barcode = depositPair.first->barcode(); + // // if( depositPair.second > highestDep) + // // { + // // highestDep = depositPair.second; + // // barcode = depositPair.first->barcode(); + // // primary = depositPair.first; + // // depositPair.first->print(std::cout); + // // ATH_MSG_INFO("pdg id "<<depositPair.first->pdg_id()); + // // } + // if (truthMap.count(barcode) > 0) + // { + // truthMap[barcode] += eDep; + // } + // else + // { + // truthMap[barcode] = eDep; + // } + // } + // } + // } + // } + } + } + std::vector<std::pair<int, float>> truth(truthMap.begin(), truthMap.end()); + std::sort(truth.begin(), truth.end(), [](auto v1, auto v2) { return v1.second > v2.second; }); + if (truth.size()>0) ATH_MSG_ALWAYS("Selected track truth info:"); + for (auto v : truth) + { + auto truthParticle = (*(std::find_if(truthParticleContainer->cbegin(), truthParticleContainer->cend(), [v](const xAOD::TruthParticle* p){ return p->barcode() == v.first; }))); + if (m_truthPdg == 0) m_truthPdg = truthParticle->pdgId(); + if (m_truthBarcode == 0) m_truthBarcode = v.first; + ATH_MSG_ALWAYS("truth info: barcode = " << v.first << " ( " << truthParticle->p4().P()/1000 << " GeV/c, Id code = " << truthParticle->pdgId() << ") -> deposited energy: " << v.second/1000); + } + } + + if (candidateParameters != nullptr) + { + m_x = candidateParameters->position().x(); + m_y = candidateParameters->position().y(); + m_z = candidateParameters->position().z(); + m_px = candidateParameters->momentum().x(); + m_py = candidateParameters->momentum().y(); + m_pz = candidateParameters->momentum().z(); + m_charge = (int) candidateParameters->charge(); } + // Here we apply the signal selection + // Very simple/unrealistic to start + if (m_vetoUpstream == 0 || m_vetoDownstream == 0 || + m_triggerTotal == 0 || + m_preshower0 == 0 || m_preshower1 == 0 || + m_ecalTotal == 0 || + candidateParameters == nullptr) + return StatusCode::SUCCESS; + + m_tree->Fill(); + return StatusCode::SUCCESS; } @@ -301,4 +617,8 @@ StatusCode NeutrinoSearchAlg::finalize() return StatusCode::SUCCESS; } - +bool NeutrinoSearchAlg::waveformHitOK(const xAOD::WaveformHit* hit) const +{ + if (hit->hit_status() == 0 || hit->hit_status() == (1<<xAOD::WaveformStatus::WAVE_OVERFLOW)) return true; + return false; +} \ No newline at end of file diff --git a/PhysicsAnalysis/NeutrinoSearch/src/NeutrinoSearchAlg.h b/PhysicsAnalysis/NeutrinoSearch/src/NeutrinoSearchAlg.h index 6f47a394..46672712 100644 --- a/PhysicsAnalysis/NeutrinoSearch/src/NeutrinoSearchAlg.h +++ b/PhysicsAnalysis/NeutrinoSearch/src/NeutrinoSearchAlg.h @@ -5,6 +5,11 @@ #include "AthenaBaseComps/AthHistogramming.h" #include "TrkTrack/TrackCollection.h" #include "xAODFaserWaveform/WaveformHitContainer.h" +#include "xAODFaserWaveform/WaveformHit.h" +#include "xAODTruth/TruthEventContainer.h" +#include "xAODTruth/TruthParticleContainer.h" +#include "TrackerPrepRawData/FaserSCT_ClusterContainer.h" +#include "TrackerSimData/TrackerSimDataCollection.h" #include "FaserActsGeometryInterfaces/IFaserActsExtrapolationTool.h" @@ -27,17 +32,24 @@ public: virtual StatusCode initialize() override; virtual StatusCode execute(const EventContext &ctx) const override; virtual StatusCode finalize() override; + bool waveformHitOK(const xAOD::WaveformHit* hit) const; const ServiceHandle <ITHistSvc> &histSvc() const; private: ServiceHandle <ITHistSvc> m_histSvc; + + SG::ReadHandleKey<xAOD::TruthEventContainer> m_truthEventContainer { this, "EventContainer", "TruthEvents", "Truth event container name." }; + SG::ReadHandleKey<xAOD::TruthParticleContainer> m_truthParticleContainer { this, "ParticleContainer", "TruthParticles", "Truth particle container name." }; + SG::ReadHandleKey<TrackerSimDataCollection> m_simDataCollection {this, "TrackerSimDataCollection", "SCT_SDO_Map"}; + 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" }; + SG::ReadHandleKey<Tracker::FaserSCT_ClusterContainer> m_clusterContainer { this, "ClusterContainer", "SCT_ClusterContainer", "Tracker cluster container name" }; ToolHandle<IFaserActsExtrapolationTool> m_extrapolationTool { this, "ExtrapolationTool", "FaserActsExtrapolationTool" }; const TrackerDD::SCT_DetectorManager* m_detMgr {nullptr}; @@ -50,11 +62,24 @@ private: 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" }; +// 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" }; + IntegerProperty m_minTrackerLayers { this, "MinTrackerLayers", 7, "Minimum number of layers with hits on track" }; + + BooleanProperty m_useFlukaWeights { this, "UseFlukaWeights", false, "Flag to weight events according to value stored in HepMC::GenEvent" }; + BooleanProperty m_useGenieWeights { this, "UseGenieWeights", false, "Flag to weight events according to Genie luminosity" }; + IntegerProperty m_flukaCollisions { this, "FlukaCollisions", 137130000, "Number of proton-proton collisions in FLUKA sample." }; + DoubleProperty m_flukaCrossSection { this, "FlukaCrossSection", 80.0, "Fluka p-p inelastic cross-section in millibarns." }; + DoubleProperty m_genieLuminosity { this, "GenieLuminosity", 150.0, "Genie luminosity in inverse fb." }; + +// BooleanProperty m_enforceBlinding { this, "EnforceBlinding", true, "Ignore data events with no VetoNu signals." }; + const bool m_enforceBlinding {true}; + + double m_baseEventCrossSection {1.0}; + const double kfemtoBarnsPerMilliBarn {1.0e12}; mutable TTree* m_tree; mutable unsigned int m_event_number; @@ -62,24 +87,42 @@ private: mutable double m_vetoNu1; mutable double m_veto00; mutable double m_veto01; + mutable double m_vetoUpstream; mutable double m_veto10; mutable double m_veto11; + mutable double m_vetoDownstream; mutable double m_trigger00; mutable double m_trigger01; mutable double m_trigger10; mutable double m_trigger11; + mutable double m_triggerTotal; 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_ecalTotal; + mutable int m_station0Clusters; + mutable int m_station1Clusters; + mutable int m_station2Clusters; + mutable int m_station3Clusters; + mutable double m_x; + mutable double m_y; + mutable double m_z; mutable double m_px; mutable double m_py; mutable double m_pz; - mutable double m_pe; + mutable int m_charge; mutable double m_chi2; + mutable int m_ndof; + mutable int m_longTracks; + mutable double m_truthLeptonMomentum; + mutable int m_truthBarcode; + mutable int m_truthPdg; + mutable double m_crossSection; + }; inline const ServiceHandle <ITHistSvc> &NeutrinoSearchAlg::histSvc() const { -- GitLab