Merge branch 'ReconSimulatedData-BrahmaSutra-Dec10-2024' into 'master'

Sequence to reconstruct and monitor simulated data from RICH in DD4HEP

See merge request !413

Rich/Panoptes/options/Rich_MC_DD4HEP_Recon_PID_Example.py

+###############################################################################
+# (c) Copyright 2019 CERN for the benefit of the LHCb Collaboration           #
+#                                                                             #
+# This software is distributed under the terms of the GNU General Public      #
+# Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING".   #
+#                                                                             #
+# In applying this licence, CERN does not waive the privileges and immunities #
+# granted to it by virtue of its status as an Intergovernmental Organization  #
+# or submit itself to any jurisdiction.                                       #
+###############################################################################
+
+import glob, os
+from GaudiConf import IOHelper
+from Gaudi.Configuration import *
+
+from ROOT import gROOT
+#from Gaudi.Configuration import VERBOSE, DEBUG, INFO
+from GaudiKernel.SystemOfUnits import GeV
+
+from Moore import options
+from DDDB.CheckDD4Hep import UseDD4Hep
+from Configurables import DDDBConf, LHCbApp
+
+options.input_type = 'ROOT'
+
+LHCbApp().CondDBtag = "sim10/run3-ideal"
+DDDBConf().GeometryVersion = "run3/trunk"
+
+options.geometry_version = "run3/trunk"
+options.conditions_version = "sim10/run3-ideal"
+
+#options.input_type = "XDIGI"
+
+#options.input_raw_format = 4.3
+
+options.simulation = True
+
+LHCbApp().Simulation = True
+LHCbApp().DataType = "Upgrade"
+
+#########################################################
+
+InputArea = "/afs/cern.ch/work/s/seaso/DataFiles/BooleOutput/RichUpgrade/Aug2024/DD4HEP/prod/Run3-trunk/DigiFiles-TimeWindow-25/"
+
+Input_Local_File0 = "Rich-Boole-TimeGateCut-25000-Prod-set-0-SetC-Run3-trunk-Sim10-Ideal-Aug-2024-DD4HEP-Extended.digi"
+Input_Local_File1 = "Rich-Boole-TimeGateCut-25000-Prod-set-1-2-SetC-Run3-trunk-Sim10-Ideal-Aug-2024-DD4HEP-Extended.digi"
+Input_Local_File2 = "Rich-Boole-TimeGateCut-25000-Prod-set-3-4-SetC-Run3-trunk-Sim10-Ideal-Aug-2024-DD4HEP-Extended.digi"
+
+Input_Local_File3 = "Rich-Boole-TimeGateCut-25000-Prod-set-5-SetC-Run3-trunk-Sim10-Ideal-Aug-2024-DD4HEP-Extended.digi"
+Input_Local_File4 = "Rich-Boole-TimeGateCut-25000-Prod-set-6-SetC-Run3-trunk-Sim10-Ideal-Aug-2024-DD4HEP-Extended.digi"
+Input_Local_File5 = "Rich-Boole-TimeGateCut-25000-Prod-set-7-8-SetC-Run3-trunk-Sim10-Ideal-Aug-2024-DD4HEP-Extended.digi"
+Input_Local_File6 = "Rich-Boole-TimeGateCut-25000-Prod-set-9-10-SetC-Run3-trunk-Sim10-Ideal-Aug-2024-DD4HEP-Extended.digi"
+
+Input_Local_File7 = "Rich-Boole-TimeGateCut-25000-Prod-set-11-12-SetC-Run3-trunk-Sim10-Ideal-Nov-2024-DD4HEP-Extended.digi"
+Input_Local_File8 = "Rich-Boole-TimeGateCut-25000-Prod-set-13-15-16-SetC-Run3-trunk-Sim10-Ideal-Nov-2024-DD4HEP-Extended.digi"
+
+#print("InputFiles ",Input_Local_File0," ",Input_Local_File1," ",Input_Local_File2," ",Input_Local_File3,"  ",Input_Local_File4,"  ",Input_Local_File5,"  ",Input_Local_File6)
+#print("InputFiles ",Input_Local_File7," ",Input_Local_File8)
+
+options.input_files = [
+    InputArea + Input_Local_File0, InputArea + Input_Local_File1,
+    InputArea + Input_Local_File2, InputArea + Input_Local_File3,
+    InputArea + Input_Local_File4, InputArea + Input_Local_File5,
+    InputArea + Input_Local_File6, InputArea + Input_Local_File7,
+    InputArea + Input_Local_File8
+]
+
+#########################################################
+
+options.evt_max = -1
+
+options.print_freq = 10
+
+############################################################
+# Configuring the detectors for reconstruction.
+
+# Define test name
+# myName = os.getenv("QMTTEST_NAME", "UNDEFINED")
+
+# Due to issues with hit sorting in FT disable Mat corrections
+# See for more details https://gitlab.cern.ch/lhcb/Panoptes/-/issues/7
+from PyConf.Algorithms import PrStoreSciFiHits
+PrStoreSciFiHits.global_bind(ApplyMatContractionCalibration=False)
+# To be removed once the issue is resolved properly.
+
+from PyConf.Algorithms import (PrForwardTrackingVelo, PrMatchNN,
+                               PrHybridSeeding)
+
+if UseDD4Hep:
+    from Configurables import LHCb__Det__LbDD4hep__DD4hepSvc as DD4hepSvc
+    dd4hep = DD4hepSvc(DetectorList=[
+        "/world", "Magnet", "VP", "FT", "UT", "Rich1", "Rich2", "Ecal", "Hcal"
+    ])
+
+    #if not noUT: dd4hep.DetectorList += ["UT"]
+
+################################################################
+#Activate tracking
+
+importOptions("$PANOPTESROOT/options/Tracking.py")
+
+import os, multiprocessing
+from Moore import options, run_reconstruction
+from Panoptes.BrahmaSutraConfiguration import BrahmaSutraConf
+
+from ROOT import gROOT
+from Gaudi.Configuration import VERBOSE, DEBUG, INFO
+from GaudiKernel.SystemOfUnits import GeV
+
+options.output_level = INFO
+#options.output_level = DEBUG
+
+# Silence ROOT
+gROOT.ProcessLine("gErrorIgnoreLevel = kWarning;")
+
+#isQMTTest = 'QMTTEST_NAME' in os.environ
+#myName = os.getenv("QMTTEST_NAME", "RichMCRecon")
+myOutputArea = "/afs/cern.ch/work/s/seaso/DataFiles/PanoptesOutput/Dec10-2024/prod/"
+myLocalName = "Rich_MC_Recon_Example_DataSet_0_16"
+myLocalNtupleName = "Rich_MC_ReconNtuple_Example_DataSet_0_16"
+
+options.histo_file = myOutputArea + myLocalName + ".root"
+options.ntuple_file = myOutputArea + myLocalNtupleName + ".root"
+
+# Override previously set defaults
+options.n_threads = 1
+
+# Is UT active in this task ?
+#noUT = ("noUT" in myName)
+
+# Override run number for MC samples
+forcedRunN = (999999 if options.simulation else 0)
+
+# Momentum Cuts
+
+minPRads = (2 * GeV, 2 * GeV, 2 * GeV)
+min_p = min(minPRads)
+
+# Enable additional monitoring algorithms ?
+#do_moni = isQMTTest
+do_moni = True
+
+# Output level for photon counting alg
+#out_level = (INFO if isQMTTest else VERBOSE)
+out_level = INFO
+options.output_level = INFO
+# Silence ROOT
+gROOT.ProcessLine("gErrorIgnoreLevel = kWarning;")
+
+from RecoConf.hlt2_tracking import make_hlt2_tracks as track_maker
+
+with PrForwardTrackingVelo.bind(MinP=min_p),\
+     PrMatchNN.bind(MinP=min_p),\
+     track_maker.bind(use_pr_kf=True, light_reco=True, fast_reco=False),\
+     BrahmaSutraConf.bind(do_data_monitoring=do_moni,ParticleGun=False):
+
+    run_reconstruction(options, BrahmaSutraConf)

Rich/Panoptes/python/Panoptes/BrahmaSutraConfiguration.py

+###############################################################################
+# (c) Copyright 2021 CERN for the benefit of the LHCb Collaboration           #
+#                                                                             #
+# This software is distributed under the terms of the GNU General Public      #
+# Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING".   #
+#                                                                             #
+# In applying this licence, CERN does not waive the privileges and immunities #
+# granted to it by virtue of its status as an Intergovernmental Organization  #
+# or submit itself to any jurisdiction.                                       #
+###############################################################################
+from GaudiKernel.SystemOfUnits import GeV
+
+from PyConf.tonic import configurable
+from PyConf.application import make_odin
+from PyConf import configurable
+
+from RecoConf.data_from_file import mc_unpacker
+
+from RecoConf.mc_checking import (
+    make_links_lhcbids_mcparticles_tracking_system,
+    make_links_tracks_mcparticles, make_links_lhcbids_mcparticles_VP_FT)
+
+from RecoConf.rich_reconstruction import (
+    default_rich_reco_options, make_rich_pids, get_detector_bool_opts)
+from RecoConf.standalone import reco_prefilters
+from RecoConf.rich_data_monitoring import (
+    default_rich_monitoring_options,
+    make_rich_pixel_monitors,
+    make_rich_track_monitors,
+)
+
+from RecoConf.rich_mc_checking import (
+    default_rich_checking_options, make_rich_checkers, default_rich_checkers)
+
+from PyConf.Algorithms import (
+    TracksToSelection, Rich__Future__MC__TrackToMCParticleRelations as
+    TkToMCPRels, Rich__Future__Rec__MC__Moni__RichPIDTupleCreatorAlg as
+    RichPIDTuple)
+
+from PyConf.Tools import TrackSelector
+
+from Moore.config import Reconstruction
+from MooreOnlineConf.utils import update_and_reset
+
+from DDDB.CheckDD4Hep import UseDD4Hep
+
+from GaudiKernel.SystemOfUnits import GeV
+
+
+@configurable
+def BrahmaSutraConf(offlineMode=True,
+                    noUT=False,
+                    do_data_monitoring=True,
+                    ParticleGun=False):
+    """
+    The tracking and RICH reco for the RICH reconstruction from MC data
+    """
+    # Make the tracks
+    if not noUT:
+        print("Configuring tracking with UT")
+        from RecoConf.hlt2_tracking import make_hlt2_tracks
+        hlt2_tracks = make_hlt2_tracks()
+    else:
+        print("Configuring tracking without UT")
+        from RecoConf.hlt2_tracking import make_hlt2_tracks_without_UT
+        hlt2_tracks = make_hlt2_tracks_without_UT()
+
+    #print("BrahmaSutraConf: Available Track Types", hlt2_tracks.keys())
+
+
+# The track type we use
+    track_version = "v1"
+
+    tkType = "Long"
+
+    # Get the fitted long tracks to use
+    # In ParticleGun events only 'Best' tracks available.
+
+    if (not ParticleGun):
+        tks = hlt2_tracks["Best" + tkType]
+    else:
+        tks = hlt2_tracks["Best"]
+
+    # CK Theta res range to use
+    #CKResHistoRange = (0.0250, 0.0055, 0.0035)
+    CKResHistoRange = (0.0055, 0.0035)
+
+    # Default RICH reco options
+    reco_opts = default_rich_reco_options()
+    reco_opts["PhotonSelection"] = "Nominal"
+
+    # Track selection
+
+    tkSel = TracksToSelection(InputLocation=tks[track_version])
+
+    # RICH photon reco
+    conf = make_rich_pids(
+        track_name=tkType,
+        input_tracks=tkSel.OutputLocation,
+        options=reco_opts)
+
+    # The detector and radiator options. Not used for now.
+    #rad_opts = get_radiator_bool_opts(reco_opts, tkType)
+    det_opts = get_detector_bool_opts(reco_opts, tkType)
+
+    algs_to_run = []
+
+    # Online histogram handling. Not used for MC data
+    if not offlineMode:
+        # Should be first in the list if required
+        algs_to_run += [update_and_reset()]
+
+    if do_data_monitoring:
+        # Data monitoring options
+        moni_opts = default_rich_monitoring_options()
+        moni_set = "Standard"
+        # Force same CK theta range as counting monitor
+        moni_opts["CKResHistoRange"] = CKResHistoRange
+        moni_algs = make_rich_track_monitors(
+            conf=conf,
+            reco_opts=reco_opts,
+            moni_opts=moni_opts,
+            moni_set=moni_set,
+        )
+        moni_algs.update(
+            make_rich_pixel_monitors(
+                conf=conf,
+                reco_opts=reco_opts,
+                moni_opts=moni_opts,
+                moni_set=moni_set,
+            ))
+        algs_to_run += [v for _, v in sorted(moni_algs.items())]
+
+        check_opts = default_rich_checking_options()
+
+        # Force same CK theta range as calibration monitor
+        check_opts["CKResHistoRange"] = CKResHistoRange
+        conf["OriginalV1Tracks"] = tks[track_version]
+
+        check_algs = make_rich_checkers(
+            conf=conf,
+            reco_opts=reco_opts,
+            check_opts=check_opts,
+            moni_set=moni_set)
+
+        algs_to_run += [v for _, v in sorted(check_algs.items())]
+
+        if (not ParticleGun):
+            make_pid_tuple = make_rich_PID_Ntuple(
+                conf=conf,
+                reco_opts=reco_opts,
+                check_opts=check_opts,
+                moni_set=moni_set)
+
+            algs_to_run += [v for _, v in sorted(make_pid_tuple.items())]
+
+    return Reconstruction('rich_MC_reconstruction', algs_to_run,
+                          reco_prefilters())
+
+
+@configurable
+def make_rich_PID_Ntuple(conf, reco_opts, check_opts, moni_set="Standard"):
+    """
+    Activates PID Ntuple from RICH. Also returns dict of the activated algorithm
+
+    Args:
+        conf       (dict): Reconstruction configuration (data) to run monitoring on
+        reco_opts  (dict): Reconstruction options
+        check_opts (dict): MC checking options
+        moni_set (string): Monitor set to activate
+
+    Returns:
+        dict of activated monitoring algorithms
+    """
+
+    # The dict of configured monitors to return
+    results = {}
+
+    # Momentum selections for performance plots
+    momentumCuts = check_opts["MomentumRanges"]
+
+    # The track name for this configuration
+    track_name = conf["TrackName"]
+
+    # The detector and radiator options
+    det_opts = get_detector_bool_opts(reco_opts, track_name)
+    #rad_opts = get_radiator_bool_opts(reco_opts, track_name)
+
+    # get the list of checkers to activate
+    checkers = default_rich_checkers(moni_set)
+
+    # MC Info
+    mcps = mc_unpacker("MCParticles")
+    richSummaries = mc_unpacker("MCRichDigitSummaries")
+
+    # Track linker stuff
+    if check_opts["UseUT"]:
+        links_to_lhcbids = make_links_lhcbids_mcparticles_tracking_system()
+    else:
+        links_to_lhcbids = make_links_lhcbids_mcparticles_VP_FT()
+
+    links_to_tracks = make_links_tracks_mcparticles(
+        InputTracks=conf["OriginalV1Tracks"], LinksToLHCbIDs=links_to_lhcbids)
+
+    # Make the RICH Track -> MCParticle relations table
+    tkMCPRels = TkToMCPRels(
+        name='RichTkToMCPRelations_{hash}',
+        TracksLocation=conf["OriginalV1Tracks"],
+        MCParticlesLinkLocation=links_to_tracks,
+        MCParticlesLocation=mcps)
+
+    # Momentum cuts for plotting etc (by radiator)
+    #pCuts = {
+    #    "MinP": (0.5 * GeV, 0.5 * GeV, 0.5 * GeV),
+    #    "MaxP": (120.0 * GeV, 120.0 * GeV, 120.0 * GeV)
+    #}
+    pCutsPID = [2.0 * GeV, 100 * GeV]
+
+    # PID performance
+    key = "PIDPerformance"
+    if key in checkers:
+        #print ("BrahmaSutra: Activating PID Ntuple Creation")
+        # Make a PID monitor for this selection
+        results[key] = RichPIDTuple(
+            name="RichPIDTupleCreatorAlg",
+            TrackSelector=TrackSelector(
+                MinPCut=pCutsPID[0], MaxPCut=pCutsPID[1]),
+            TracksLocation=conf["InputTracks"],
+            RichPIDsLocation=conf["RichPIDs"],
+            TrackToMCParticlesRelations=tkMCPRels.TrackToMCParticlesRelations)
+    return results