fix: store systematic in array option

.github/workflows/BTA_workflow.yml

@@ -107,7 +107,7 @@ jobs:
          python runner.py --wf BTA_addAllTracks --json metadata/test_bta_run3.json --executor iterative --overwrite --campaign Summer23 --year 2023
     - name: BTA addPFMuons workflow test(BTA for SF)
       run: |
-         python runner.py --wf BTA_addPFMuons --json metadata/test_bta_run3.json --executor iterative --overwrite --campaign Summer23 --year 2023
+         python runner.py --wf BTA_addPFMuons --json metadata/test_bta_run3.json --executor iterative --overwrite --campaign Summer23 --year 2023 
 
     - name: BTA_ttbar workflow test
       run: |

.github/workflows/QCD_workflow.yml

@@ -117,7 +117,7 @@ jobs:
         elif [[ $string == *"ci:weight_only"* ]]; then
             opts=$(echo "$opts" | sed 's/--isSyst all/--isSyst weight_only/g') 
         fi
-        python runner.py --workflow QCD_sf --json metadata/test_bta_run3.json  --executor iterative  --campaign Summer23 --year 2023 $opts
+        python runner.py --workflow QCD_sf --json metadata/test_bta_run3.json  --executor iterative  --campaign Summer23 --year 2023 $opts  --isSyst all 
 
     - name: QCD soft mu workflows with correctionlib
       run: |
@@ -133,5 +133,5 @@ jobs:
             opts=$(echo "$opts" | sed 's/--isSyst all/--isSyst weight_only/g') 
         fi
 
-        python runner.py --workflow QCD_smu_sf --json metadata/test_bta_run3.json --executor iterative  --campaign Summer23 --year 2023 --overwrite $opts
+        python runner.py --workflow QCD_smu_sf --json metadata/test_bta_run3.json --executor iterative  --campaign Summer23 --year 2023 --overwrite $opts --isSyst all 
 

.github/workflows/ctag_DY_workflow.yml

@@ -115,7 +115,7 @@ jobs:
         elif [[ $string == *"ci:weight_only"* ]]; then
             opts=$(echo "$opts" | sed 's/--isSyst all/--isSyst weight_only/g') 
         fi
-        python runner.py --workflow ctag_DY_sf --json metadata/test_bta_run3.json --limit 1 --executor iterative --campaign Summer23 --year 2023  $opts
+        python runner.py --workflow ctag_DY_sf --json metadata/test_bta_run3.json --limit 1 --executor iterative --campaign Summer23 --year 2023  $opts  --isSyst all --isArray
 
 
     - name: ctag electron DY workflows
@@ -131,4 +131,4 @@ jobs:
         elif [[ $string == *"ci:weight_only"* ]]; then
             opts=$(echo "$opts" | sed 's/--isSyst all/--isSyst weight_only/g') 
         fi
-        python runner.py --workflow ectag_DY_sf --json metadata/test_bta_run3.json --limit 1 --executor iterative --campaign Summer23 --year 2023 $opts 
+        python runner.py --workflow ectag_DY_sf --json metadata/test_bta_run3.json --limit 1 --executor iterative --campaign Summer23 --year 2023 $opts   --isSyst all  --isArray

.github/workflows/ctag_Wc_workflow.yml

@@ -117,7 +117,7 @@ jobs:
         elif [[ $string == *"ci:weight_only"* ]]; then
             opts=$(echo "$opts" | sed 's/--isSyst all/--isSyst weight_only/g') 
         fi
-        python runner.py --workflow ctag_Wc_sf --json metadata/test_bta_run3.json  --executor iterative --campaign Summer23 --year 2023 $opts
+        python runner.py --workflow ctag_Wc_sf --json metadata/test_bta_run3.json  --executor iterative --campaign Summer23 --year 2023 $opts  --isSyst all  --isArray 
 
     - name: ctag electron W+c workflows
       run: |
@@ -132,7 +132,7 @@ jobs:
         elif [[ $string == *"ci:weight_only"* ]]; then
             opts=$(echo "$opts" | sed 's/--isSyst all/--isSyst weight_only/g') 
         fi
-        python runner.py --workflow ectag_Wc_sf --json metadata/test_bta_run3.json --executor iterative  --overwrite  --campaign Summer23 --year 2023 $opts
+        python runner.py --workflow ectag_Wc_sf --json metadata/test_bta_run3.json --executor iterative  --overwrite  --campaign Summer23 --year 2023 $opts  --isSyst all  --isArray 
 
     - name: ctag muon W+c fixed WP workflows with correctionlib
       run: |
@@ -147,7 +147,7 @@ jobs:
         elif [[ $string == *"ci:weight_only"* ]]; then
             opts=$(echo "$opts" | sed 's/--isSyst all/--isSyst weight_only/g') 
         fi
-        python runner.py --workflow ctag_Wc_WP_sf --json metadata/test_bta_run3.json  --executor iterative  --overwrite  --campaign Summer23 --year 2023 $opts
+        python runner.py --workflow ctag_Wc_WP_sf --json metadata/test_bta_run3.json  --executor iterative  --overwrite  --campaign Summer23 --year 2023 $opts  --isSyst all 
 
     - name: ctag electron W+c fixed WP workflows with correctionlib
       run: |
@@ -162,4 +162,4 @@ jobs:
         elif [[ $string == *"ci:weight_only"* ]]; then
             opts=$(echo "$opts" | sed 's/--isSyst all/--isSyst weight_only/g') 
         fi
-        python runner.py --workflow ectag_Wc_WP_sf --json metadata/test_bta_run3.json --executor iterative  --overwrite  --campaign Summer23 --year 2023 $opts
+        python runner.py --workflow ectag_Wc_WP_sf --json metadata/test_bta_run3.json --executor iterative  --overwrite  --campaign Summer23 --year 2023 $opts --isSyst all 

docs/developer.md

@@ -129,7 +129,7 @@ if not self.noHist:
     )
 # Output arrays - store the pruned objects in the output arrays
 if self.isArray:
-    array_writer(self, pruned_ev, events, systematics[0], dataset, isRealData)
+    array_writer(self, pruned_ev, events, weights, systematics,dataset, isRealData)
   ```
 
 </p>

scripts/fetch.py

@@ -17,7 +17,6 @@
     "--input",
     default=None,
     type=str,
-    # required=True,
     help="List of samples in DAS (default: %(default)s)",
 )
 parser.add_argument(
@@ -106,7 +105,7 @@
     "--campaign",
     help="campaign name (same as the campaign in runner.py)",
     default=None,
-    require=True,
+    required=True,
     type=str,
 )
 parser.add_argument("--year", help="year", default=None, type=str)

src/BTVNanoCommissioning/utils/array_writer.py

@@ -8,6 +8,7 @@ def array_writer(
     processor_class,  # the NanoProcessor class ("self")
     pruned_event,  # the event with specific calculated variables stored
     nano_event,  # entire NanoAOD/PFNano event with many variables
+    weights,  # weight for the event
     systname,  # name of systematic shift
     dataset,  # dataset name
     isRealData,  # boolean
@@ -49,6 +50,16 @@ def array_writer(
     othersMC=["Pileup_nTrueInt", "Pileup_nPU"],  # other fields, for MC only
     empty=False,
 ):
+    if not isRealData and systname != ["nominal"]:
+        pruned_event["weight"] = weights.weight()
+        for ind_wei in weights.weightStatistics.keys():
+            pruned_event[f"{ind_wei}_weight"] = weights.partial_weight(
+                include=[ind_wei]
+            )
+        if len(systname) > 1:
+            for syst in systname[1:]:
+                pruned_event[f"weight_syst_{syst}"] = weights.weight(modifier=syst)
+
     if empty:
         print("WARNING: No events selected. Writing blank file.")
         out_branch = []
@@ -81,7 +92,7 @@ def array_writer(
 
     # Write to root files
     print("Branches to write:", out_branch)
-    outdir = f"{processor_class.name}/{systname}/{dataset}/"
+    outdir = f"{processor_class.name}/{systname[0]}/{dataset}/"
     os.system(f"mkdir -p {outdir}")
 
     with uproot.recreate(

src/BTVNanoCommissioning/utils/correction.py

@@ -613,7 +613,7 @@ def JME_shifts(
             met["orig_pt"], met["orig_phi"] = nocorrmet["pt"], nocorrmet["phi"]
 
             ## JEC variations
-            if systematic != False:
+            if not isRealData and systematic != False:
                 if systematic != "JERC_split":
                     jesuncmap = correct_map["JME"][f"{jecname}_Total_AK4PFPuppi"]
                     jesunc = ak.unflatten(jesuncmap.evaluate(j.eta, j.pt), nj)
@@ -2079,6 +2079,7 @@ def common_shifts(self, events):
     return shifts
 
 
+# common weights
 def weight_manager(pruned_ev, SF_map, isSyst):
     """
     Example for Scaling Factors (SFs):
@@ -2098,7 +2099,7 @@ def weight_manager(pruned_ev, SF_map, isSyst):
         add_pdf_weight(weights, pruned_ev.LHEPdfWeight, isSyst)
     if "LHEScaleWeight" in pruned_ev.fields:
         add_scalevar_weight(weights, pruned_ev.LHEScaleWeight, isSyst)
-    if "TTTo" in pruned_ev.metadata["dataset"]:
+    if "TT" in pruned_ev.metadata["dataset"]:
         weights.add(
             "ttbar_weight",
             top_pT_reweighting(pruned_ev.GenPart),
@@ -2108,8 +2109,18 @@ def weight_manager(pruned_ev, SF_map, isSyst):
             )
             * 2.0
             + ak.ones_like(top_pT_reweighting(pruned_ev.GenPart)),
-            ak.ones_like(top_pT_reweighting(pruned_ev.GenPart)),
         )
+        if isSyst != False:
+            weights.add(
+                "ttbar_weight",
+                top_pT_reweighting(pruned_ev.GenPart),
+                (
+                    top_pT_reweighting(pruned_ev.GenPart)
+                    - ak.ones_like(top_pT_reweighting(pruned_ev.GenPart))
+                )
+                * 2.0,
+                ak.ones_like(top_pT_reweighting(pruned_ev.GenPart)),
+            )
 
     if "hadronFlavour" in pruned_ev.Jet.fields:
 

src/BTVNanoCommissioning/utils/histogrammer.py

@@ -76,7 +76,7 @@ def histogrammer(events, workflow, year="2022", campaign="Summer22"):
         # _hist_dict["nJetSVs"] = Hist.Hist(syst_axis, n_axis, Hist.storage.Weight())
 
     elif "QCD_smu" == workflow:
-        obj_list = ["jet0", "hl", "soft_l"]
+        obj_list = ["mujet", "hl", "soft_l"]
         # _hist_dict["dr_SVjet0"] = Hist.Hist(
         #     syst_axis, flav_axis, dr_SV_axis, Hist.storage.Weight()
         # )
@@ -705,17 +705,30 @@ def histo_writter(pruned_ev, output, weights, systematics, isSyst, SF_map):
         "DeepJetC",
     ]  # exclude b-tag SFs for btag inputs
     # define Jet flavor
-    nj = 1 if type(pruned_ev.SelJet.pt[0]) == float else len(pruned_ev.SelJet.pt[0])
+
+    # Reduce the jet to the correct dimension in the plot
+    nj = 4 if "jet4" in output.keys() else 2 if "jet2" in output.keys() else 1
+    pruned_ev.SelJet = pruned_ev.SelJet if nj == 1 else pruned_ev.SelJet[:, :nj]
+    if "var" in str(ak.type(pruned_ev.SelJet.pt)) and nj == 1:
+        pruned_ev.SelJet = pruned_ev.SelJet[:, 0]
     if "hadronFlavour" in pruned_ev.SelJet.fields:
         isRealData = False
         genflavor = ak.values_astype(
-            pruned_ev.SelJet.hadronFlavour + (1 * (pruned_ev.SelJet.partonFlavour == 0) & (pruned_ev.SelJet.hadronFlavour == 0)),
+            pruned_ev.SelJet.hadronFlavour
+            + 1
+            * (
+                (pruned_ev.SelJet.partonFlavour == 0)
+                & (pruned_ev.SelJet.hadronFlavour == 0)
+            ),
             int,
         )
         if "MuonJet" in pruned_ev.fields:
             smflav = ak.values_astype(
-                1 * (pruned_ev.MuonJet.partonFlavour == 0)
-                & (pruned_ev.MuonJet.hadronFlavour == 0)
+                1
+                * (
+                    (pruned_ev.MuonJet.partonFlavour == 0)
+                    & (pruned_ev.MuonJet.hadronFlavour == 0)
+                )
                 + pruned_ev.MuonJet.hadronFlavour,
                 int,
             )
@@ -924,30 +937,22 @@ def histo_writter(pruned_ev, output, weights, systematics, isSyst, SF_map):
                         weight=weights.partial_weight(exclude=exclude_btv),
                     )
 
-        # pT ratio
-        if "hl" in pruned_ev.fields:
-            output["hl_ptratio"].fill(
-                syst,
-                genflavor[:, 0],
-                ratio=pruned_ev.hl.pt / pruned_ev.SelJet[:, 0].pt,
-                weight=weight,
-            )
-        if "sl" in pruned_ev.fields:
-            output["sl_ptratio"].fill(
-                syst,
-                genflavor[:, 0],
-                ratio=pruned_ev.sl.pt / pruned_ev.SelJet[:, 0].pt,
-                weight=weight,
-            )
         if "dr_poslnegl" in output.keys():
+            # DY histograms
             output["dr_poslnegl"].fill(
                 syst, pruned_ev.posl.delta_r(pruned_ev.negl), weight=weight
             )
             output["dr_posljet"].fill(
-                syst, genflavor, pruned_ev.posl.delta_r(pruned_ev.SelJet), weight=weight
+                syst,
+                genflavor,
+                pruned_ev.posl.delta_r(pruned_ev.SelJet),
+                weight=weight,
             )
             output["dr_negljet"].fill(
-                syst, genflavor, pruned_ev.negl.delta_r(pruned_ev.SelJet), weight=weight
+                syst,
+                genflavor,
+                pruned_ev.negl.delta_r(pruned_ev.SelJet),
+                weight=weight,
             )
         # Muon enriched jet histograms
         if "MuonJet" in pruned_ev.fields:
@@ -993,6 +998,21 @@ def histo_writter(pruned_ev, output, weights, systematics, isSyst, SF_map):
                 dr=pruned_ev.MuonJet.delta_r(pruned_ev.SoftMuon),
                 weight=weight,
             )
+            if "hl" in pruned_ev.fields:
+                output["hl_ptratio"].fill(
+                    syst,
+                    smflav,
+                    ratio=pruned_ev.hl.pt / pruned_ev.MuonJet.pt,
+                    weight=weight,
+                )
+            if "sl" in pruned_ev.fields:
+                output["sl_ptratio"].fill(
+                    syst,
+                    smflav,
+                    ratio=pruned_ev.sl.pt / pruned_ev.MuonJet.pt,
+                    weight=weight,
+                )
+
             if "SelMuon" in pruned_ev.fields and "hl" not in pruned_ev.fields:
                 output["dr_lmujethmu"].fill(
                     syst,

src/BTVNanoCommissioning/workflows/QCD_soft_mu_validation.py

@@ -131,7 +131,8 @@ def process_shift(self, events, shift_name):
                     self,
                     events[event_level],
                     events,
-                    "nominal",
+                    None,
+                    ["nominal"],
                     dataset,
                     isRealData,
                     empty=True,
@@ -258,20 +259,17 @@ def process_shift(self, events, shift_name):
             systematics = ["nominal"] + list(weights.variations)
         else:
             systematics = [shift_name]
-        if not isRealData:
-            pruned_ev["weight"] = weights.weight()
-            for ind_wei in weights.weightStatistics.keys():
-                pruned_ev[f"{ind_wei}_weight"] = weights.partial_weight(
-                    include=[ind_wei]
-                )
+
         # Configure histograms
         if not self.noHist:
             output = histo_writter(
                 pruned_ev, output, weights, systematics, self.isSyst, self.SF_map
             )
         # Output arrays
         if self.isArray:
-            array_writer(self, pruned_ev, events, systematics[0], dataset, isRealData)
+            array_writer(
+                self, pruned_ev, events, weights, systematics, dataset, isRealData
+            )
         return {dataset: output}
 
     def postprocess(self, accumulator):

src/BTVNanoCommissioning/workflows/QCD_validation.py

@@ -96,7 +96,6 @@ def process_shift(self, events, shift_name):
             valid_events = (ak.count(pruned_ev.Jet.pt, axis=1) > 0) & (
                 ak.count(pruned_ev.JetSVs.pt, axis=1) > 0
             )
-            print("valid_events", valid_events, len(valid_events), len(events))
             filtered_events = pruned_ev[valid_events]
 
             # Pad pruned_ev.JetSVs.jetIdx to match the length of pruned_ev.Jet
@@ -190,20 +189,17 @@ def process_shift(self, events, shift_name):
             systematics = ["nominal"] + list(weights.variations)
         else:
             systematics = [shift_name]
-        if not isRealData:
-            pruned_ev["weight"] = weights.weight()
-            for ind_wei in weights.weightStatistics.keys():
-                pruned_ev[f"{ind_wei}_weight"] = weights.partial_weight(
-                    include=[ind_wei]
-                )
+
         # Configure histograms
         if not self.noHist:
             output = histo_writter(
                 pruned_ev, output, weights, systematics, self.isSyst, self.SF_map
             )
         # Output arrays
         if self.isArray:
-            array_writer(self, pruned_ev, events, systematics[0], dataset, isRealData)
+            array_writer(
+                self, pruned_ev, events, weights, systematics, dataset, isRealData
+            )
 
         return {dataset: output}
 

src/BTVNanoCommissioning/workflows/ctag_DY_valid_sf.py

@@ -198,7 +198,8 @@ def process_shift(self, events, shift_name):
                     self,
                     events[event_level],
                     events,
-                    "nominal",
+                    None,
+                    ["nominal"],
                     dataset,
                     isRealData,
                     empty=True,
@@ -222,7 +223,7 @@ def process_shift(self, events, shift_name):
         )
         # Keep the structure of events and pruned the object size
         pruned_ev = events[event_level]
-        pruned_ev["SelJet"] = event_jet[event_level][:, 0]
+        pruned_ev["SelJet"] = event_jet[event_level]
         if isMu:
             pruned_ev["MuonPlus"] = sposmu
             pruned_ev["MuonMinus"] = snegmu
@@ -260,20 +261,17 @@ def process_shift(self, events, shift_name):
             systematics = ["nominal"] + list(weights.variations)
         else:
             systematics = [shift_name]
-        if not isRealData:
-            pruned_ev["weight"] = weights.weight()
-            for ind_wei in weights.weightStatistics.keys():
-                pruned_ev[f"{ind_wei}_weight"] = weights.partial_weight(
-                    include=[ind_wei]
-                )
+
         # Configure histograms
         if not self.noHist:
             output = histo_writter(
                 pruned_ev, output, weights, systematics, self.isSyst, self.SF_map
             )
         # Output arrays
         if self.isArray:
-            array_writer(self, pruned_ev, events, systematics[0], dataset, isRealData)
+            array_writer(
+                self, pruned_ev, events, weights, systematics, dataset, isRealData
+            )
 
         return {dataset: output}