Update/cleanup and 2024

law.cfg

@@ -22,6 +22,7 @@ topsf.tasks.inference_v2.post_fit_shapes
 topsf.tasks.inference_v2.impacts
 topsf.tasks.inference_v2.plot_impacts
 topsf.tasks.inference_v2.plot_shapes
+topsf.tasks.corrections
 
 [logging]
 
@@ -35,7 +36,7 @@ columnflow.columnar_util-perf: INFO
 
 default_analysis: topsf.config.run3.analysis_sf.analysis_sf
 default_config: run3_sf_2022_postEE_nano_v12
-default_dataset: tt_fh_powheg
+default_dataset: tt_sl_powheg
 run3_analysis: topsf.config.run3.analysis_sf.analysis_sf
 run3_config: run3_sf_2022_preEE_nano_v12
 
@@ -44,6 +45,7 @@ default_keep_reduced_events: True
 production_modules: columnflow.production.{categories,normalization,mc_weight,pileup,processes,seeds}, columnflow.production.cms.{btag,electron,mc_weight,muon,pdf,pileup,scale,seeds}, topsf.production.{default,gen_top}
 calibration_modules: columnflow.calibration.cms.{jets,met}, topsf.calibration.{default,skip_jec}
 selection_modules: columnflow.selection.cms.{json_filter,met_filters}, topsf.selection.{default,categories,jet,bjet,fatjet,lepton,wp}
+hist_production_modules: topsf.weights.default
 ml_modules: columnflow.ml
 inference_modules: columnflow.inference, topsf.inference.{default,uhh2}
 
@@ -58,7 +60,7 @@ skip_ensure_proxy: False
 
 # some remote workflow parameter defaults
 htcondor_flavor: $CF_HTCONDOR_FLAVOR
-htcondor_share_software: False
+htcondor_share_software: True
 slurm_flavor: $CF_SLURM_FLAVOR
 slurm_partition: $CF_SLURM_PARTITION
 
@@ -81,10 +83,11 @@ log_array_function_runtime: False
 [outputs]
 
 # list of all used file systems
-wlcg_file_systems: wlcg_fs, wlcg_fs_cernbox, wlcg_fs_infn_redirector, wlcg_fs_global_redirector
+wlcg_file_systems: wlcg_fs, wlcg_fs_desy, wlcg_fs_cernbox, wlcg_fs_desy_store, wlcg_fs_infn_redirector, wlcg_fs_global_redirector
 
 # list of file systems used by columnflow.tasks.external.GetDatasetLFNs.iter_nano_files to
 # look for the correct fs per nano input file (in that order)
+# lfn_sources: local_desy_dcache, wlcg_fs_desy_store, wlcg_fs_infn_redirector, wlcg_fs_global_redirector
 lfn_sources: local_desy_dcache, wlcg_fs_infn_redirector, wlcg_fs_global_redirector
 
 # output locations per task family
@@ -122,6 +125,7 @@ task_cf.CreateDatacards: local
 
 job_file_dir: $CF_JOB_BASE
 job_file_dir_cleanup: False
+job_file_dir_mkdtemp: sub_{{task_id}}_XXX
 
 
 [local_fs]
@@ -139,6 +143,7 @@ base: /pnfs/desy.de/cms/tier2
 # set this to your desired location, e.g.:
 # base: root://eosuser.cern.ch/eos/user/$CF_CERN_USER_FIRSTCHAR/$CF_CERN_USER/$CF_STORE_NAME
 base: &::wlcg_fs_desy::base
+base_mkdir_rec: &::wlcg_fs_desy::webdav_base
 create_file_dir: True
 use_cache: $CF_WLCG_USE_CACHE
 cache_root: $CF_WLCG_CACHE_ROOT
@@ -152,6 +157,23 @@ xrootd_base: root://dcache-cms-xrootd.desy.de:1094/pnfs/desy.de/cms/tier2/store/
 gsiftp_base: gsiftp://dcache-cms-gridftp.desy.de/pnfs/desy.de/cms/tier2/store/user/$CF_CERN_USER/$CF_STORE_NAME
 webdav_base: davs://dcache-cms-webdav-wan.desy.de:2880/pnfs/desy.de/cms/tier2/store/user/$CF_CERN_USER/$CF_STORE_NAME
 base: &::webdav_base
+base_filecopy: &::webdav_base
+base_stat: &::webdav_base
+
+
+[wlcg_fs_desy_store]
+
+xrootd_base: root://dcache-cms-xrootd.desy.de:1094/pnfs/desy.de/cms/tier2
+gsiftp_base: gsiftp://dcache-door-cms04.desy.de:2811/pnfs/desy.de/cms/tier2
+webdav_base: davs://dcache-cms-webdav-wan.desy.de:2880/pnfs/desy.de/cms/tier2
+base: &::webdav_base
+use_cache: $CF_WLCG_USE_CACHE
+cache_root: $CF_WLCG_CACHE_ROOT
+cache_cleanup: $CF_WLCG_CACHE_CLEANUP
+cache_max_size: 15GB
+cache_global_lock: True
+cache_mtime_patience: -1
+rucio_report_access: T2_DE_DESY
 
 
 [wlcg_fs_cernbox]
@@ -170,6 +192,7 @@ cache_cleanup: $CF_WLCG_CACHE_CLEANUP
 cache_max_size: 15GB
 cache_global_lock: True
 cache_mtime_patience: -1
+rucio_report_access: True
 
 
 [wlcg_fs_global_redirector]
@@ -181,6 +204,7 @@ cache_cleanup: $CF_WLCG_CACHE_CLEANUP
 cache_max_size: 15GB
 cache_global_lock: True
 cache_mtime_patience: -1
+rucio_report_access: True
 
 
 [versions]
@@ -213,7 +237,7 @@ cache_mtime_patience: -1
 
 # for quick changes of the scheduler
 local_scheduler: false
-scheduler_host: naf-cms12.desy.de
+scheduler_host: naf-cms14.desy.de
 scheduler_port: 8082
 
 [luigi_resources]

topsf/calibration/default.py

@@ -4,82 +4,218 @@
 Calibration methods.
 """
 import functools
+import law
 
 from columnflow.calibration import Calibrator, calibrator
-from columnflow.calibration.cms.jets import jets_ak4, jets_ak8
+from columnflow.calibration.cms.jets import jec_ak4, jer_ak4, jec_ak8, jer_ak8
+from columnflow.calibration.cms.met import met_phi
 from columnflow.production.cms.mc_weight import mc_weight
 from columnflow.production.cms.jet import msoftdrop
 from columnflow.production.cms.seeds import deterministic_seeds
 from columnflow.util import maybe_import
 from columnflow.columnar_util import set_ak_column
+from columnflow.production.cms.jet import jet_id, fatjet_id
 
 from topsf.calibration.jets import (
     jet_lepton_cleaner,
     jec_subjets,
-    # jer_subjets
+    jer_subjets
 )
+from topsf.util import has_tag, record_calls
+
+logger = law.logger.get_logger(__name__)
 
 np = maybe_import("numpy")
 ak = maybe_import("awkward")
 
 set_ak_column_f32 = functools.partial(set_ak_column, value_type=np.float32)
 
+jec_ak4_Puppi = jec_ak4.derive(
+    "jec_ak4_Puppi",
+    cls_dict={
+        "met_name": "PuppiMET",
+        "raw_met_name": "RawPuppiMET",
+    }
+)
+jer_ak4_Puppi = jer_ak4.derive(
+    "jer_ak4_Puppi",
+    cls_dict={
+        "met_name": "PuppiMET",
+        "raw_met_name": "RawPuppiMET",
+    }
+)
+
+jec_ak8_Puppi = jec_ak8.derive(
+    "jec_ak8_Puppi",
+    cls_dict={
+        "propagate_met": False,
+        "met_name": "DO_NOT_USE",
+        "raw_met_name": "DO_NOT_USE",
+    },
+)
+
+jer_ak8_Puppi = jer_ak8.derive(
+    "jer_ak8_Puppi",
+    cls_dict={
+        "propagate_met": False,
+        "met_name": "DO_NOT_USE",
+        "raw_met_name": "DO_NOT_USE",
+    }
+)
+
+jec_subjets_Puppi = jec_subjets.derive(
+    "jec_subjets_Puppi",
+    cls_dict={
+        "propagate_met": False,
+        "met_name": "DO_NOT_USE",
+        "raw_met_name": "DO_NOT_USE",
+    }
+)
+
+jer_subjets_Puppi = jer_subjets.derive(
+    "jer_subjets_Puppi",
+    cls_dict={
+        "propagate_met": False,
+        "met_name": "DO_NOT_USE",
+        "raw_met_name": "DO_NOT_USE",
+    }
+)
+
 
 @calibrator(
     uses={
         mc_weight,
         deterministic_seeds,
         jet_lepton_cleaner,
-        jets_ak4,
-        jets_ak8,
-        jec_subjets,
-        # jer_subjets,
         msoftdrop,
+        "Muon.pt", "Muon.tunepRelPt",
     },
     produces={
         mc_weight,
         deterministic_seeds,
         jet_lepton_cleaner,
-        jets_ak4,
-        jets_ak8,
-        jec_subjets,
-        # jer_subjets,
         msoftdrop,
+        "Muon.pt", "Muon.rawPt",
     },
 )
 def default(self: Calibrator, events: ak.Array, **kwargs) -> ak.Array:
-    events = self[jet_lepton_cleaner](events, **kwargs)  # set Jet.pt to raw Pt -> run before jets (JER) calibrator
-    if self.dataset_inst.is_mc:
-        events = self[mc_weight](events, **kwargs)
+    run_list = []
+    with record_calls(self, run_list):
+        # highPt muons: use tuneP pT
+        events = set_ak_column_f32(events, "Muon.rawPt", events.Muon.pt)
+        events = set_ak_column_f32(events, "Muon.pt", events.Muon.tunepRelPt * events.Muon.pt)
+        logger.info_once(
+            "Finished recalculating muon pt with tuneP for highPt muons. Stored original pt in Muon.rawPt."
+        )
+        if self.dataset_inst.is_mc:
+            events = self[mc_weight](events, **kwargs)
+        events = self[deterministic_seeds](events, **kwargs)
+
+        events = self[jet_lepton_cleaner](events, **kwargs)  # set Jet.pt to raw Pt -> run before jets (JER) calibrator
+        # run JEC calibrators for AK4 and AK8 jets
+        # fake subjet area column by setting it to an array with the same structure as the subjet pt column containing 0.5
+        # (needed to be able to use same code as for top-level AK4/AK8 jets, as the producer formally requires an `area`
+        # column, despite not actually using it)
+        events = set_ak_column_f32(events, "SubJet.area", 0.5 * ak.ones_like(events.SubJet.pt))
+        if self.config_inst.x.year == 2024:
+            events = self[jec_ak4_Puppi](events, **kwargs)
+            events = self[jec_ak8_Puppi](events, **kwargs)
+            events = self[jec_subjets_Puppi](events, **kwargs)
+            if self.dataset_inst.is_mc:
+                events = self[jer_ak4_Puppi](events, **kwargs)
+                events = self[jer_ak8_Puppi](events, **kwargs)
+                # events = self[jer_subjets_Puppi](events, **kwargs)
+        else:
+            events = self[jec_ak4](events, **kwargs)
+            events = self[jec_ak8](events, **kwargs)
+            events = self[jec_subjets](events, **kwargs)
+            if self.dataset_inst.is_mc:
+                events = self[jer_ak4](events, **kwargs)
+                events = self[jer_ak8](events, **kwargs)
+                # events = self[jer_subjets](events, **kwargs)
+
+        events = self[msoftdrop](events, **kwargs)
+        if self.config_inst.x.year in {2022, 2023}:
+            events = self[met_phi](events, **kwargs)
+        elif self.config_inst.x.year == 2024:
+            logger.warning_once("met_phi calibrator not run for 2024 config, as it is not yet available.")
+        else:
+            raise ValueError(f"Unsupported year {self.config_inst.x.year} in default calibrator")
+
+        if not has_tag("skip_jet_ids", self.config_inst, self.dataset_inst, operator=any):
+            logger.debug("Recalulating (fat)jet IDs.")
+            events = self[jet_id](events, **kwargs)
+            events = self[fatjet_id](events, **kwargs)
+
+    logger.info_once(
+        "Finished default calibration steps:\n" +
+        "\n".join(run_list)
+    )
 
-    events = self[jets_ak4](events, **kwargs)  # call_force ?
-    events = self[jets_ak8](events, **kwargs)  # call_force ?
-
-    events = self[deterministic_seeds](events, **kwargs)
-
-    # fake subjet area column by setting it to an array with the same structure as the subjet pt column containing 0.5
-    # (needed to be able to use same code as for top-level AK4/AK8 jets, as the producer formally requires an `area`
-    # column, despite not actually using it)
-    events = set_ak_column_f32(events, "SubJet.area", 0.5 * ak.ones_like(events.SubJet.pt))
+    return events
 
-    events = self[jec_subjets](events, **kwargs)
-    # if self.dataset_inst.is_mc:
-    #     events = self[jer_subjets](events, **kwargs)
-    events = self[msoftdrop](events, **kwargs)
 
-    return events
+@default.init
+def default_init(self: Calibrator) -> None:
+    # add met_phi only for 2022/23 configs
+    if self.config_inst.x.year in {2022, 2023}:
+        self.uses |= {
+            met_phi,
+            jec_ak4,
+            jec_ak8,
+            jer_ak4,
+            jer_ak8,
+            jec_subjets,
+            jer_subjets,
+        }
+        self.produces |= {
+            met_phi,
+            jec_ak4,
+            jec_ak8,
+            jer_ak4,
+            jer_ak8,
+            jec_subjets,
+            jer_subjets,
+        }
+    elif self.config_inst.x.year == 2024:
+        self.uses |= {
+            jec_ak4_Puppi,
+            jec_ak8_Puppi,
+            jer_ak4_Puppi,
+            jer_ak8_Puppi,
+            jec_subjets_Puppi,
+            jer_subjets_Puppi,
+        }
+        self.produces |= {
+            jec_ak4_Puppi,
+            jec_ak8_Puppi,
+            jer_ak4_Puppi,
+            jer_ak8_Puppi,
+            jec_subjets_Puppi,
+            jer_subjets_Puppi,
+        }
+
+    if not has_tag("skip_jet_ids", self.config_inst, self.dataset_inst, operator=any):
+        self.uses |= {
+            jet_id,
+            fatjet_id,
+        }
+        self.produces |= {
+            jet_id,
+            fatjet_id,
+        }
 
 
 @calibrator(
-    uses={mc_weight, deterministic_seeds, jets_ak4, jets_ak8},
-    produces={mc_weight, deterministic_seeds, jets_ak4, jets_ak8},
+    uses={mc_weight, deterministic_seeds},
+    produces={mc_weight, deterministic_seeds},
 )
 def no_jet_cleaning(self: Calibrator, events: ak.Array, **kwargs) -> ak.Array:
     if self.dataset_inst.is_mc:
         events = self[mc_weight](events, **kwargs)
 
-    events = self[jets_ak4](events, **kwargs)
-    events = self[jets_ak8](events, **kwargs)  # call_force ?
+    # events = self[jets_ak4](events, **kwargs)
+    # events = self[jets_ak8](events, **kwargs)  # call_force ?
     events = self[deterministic_seeds](events, **kwargs)
 
     return events

topsf/calibration/jets.py

@@ -44,8 +44,10 @@ def jet_lepton_cleaner(self: Calibrator, events: ak.Array, **kwargs) -> ak.Array
 
     # revert JEC for jet pt and jet mass,
     # set correction factor to 0
-    events = set_ak_column(events, "Jet.pt", events.Jet.pt * (1 - events.Jet.rawFactor))
-    events = set_ak_column(events, "Jet.mass", events.Jet.mass * (1 - events.Jet.rawFactor))
+    raw_pt = events.Jet.pt * (1 - events.Jet.rawFactor)
+    raw_mass = events.Jet.mass * (1 - events.Jet.rawFactor)
+    events = set_ak_column(events, "Jet.pt", raw_pt)
+    events = set_ak_column(events, "Jet.mass", raw_mass)
     events = set_ak_column(events, "Jet.rawFactor", 0)
 
     # build jet lorentz vectors
@@ -158,6 +160,10 @@ def jet_lepton_cleaner(self: Calibrator, events: ak.Array, **kwargs) -> ak.Array
         value = ak.fill_none(ak.nan_to_none(getattr(jet_lv, var)), 0.0)
         events = set_ak_column(events, f"Jet.{var}", value)
 
+    # do not set to updated raw factor as out dated JECs are already reverted
+    # raw_factor = ak.nan_to_num(1 - raw_pt / events.Jet.pt, nan=0.0)
+    # events = set_ak_column(events, "Jet.rawFactor", raw_factor)
+
     return events