WIP: QG scale factor workflows

.gitignore

@@ -20,3 +20,4 @@ dask-worker-space/
 .job_wrapper_failure
 jobs_*
 .success
+arrays*/

condor_lxplus/execute.sh

@@ -1,8 +1,11 @@
 #!/bin/bash -xe
+export EOS_MGM_URL=root://eosuser.cern.ch
 
 JOBID=$1
 COMMDIR=$2
 
+echo "Running job $JOBID"
+
 export HOME=`pwd`
 if [ -d /afs/cern.ch/user/${USER:0:1}/$USER ]; then
   export HOME=/afs/cern.ch/user/${USER:0:1}/$USER  # crucial on lxplus condor but cannot set on cmsconnect
@@ -18,13 +21,17 @@ voms-proxy-info
 export PATH="$4:$PATH" 
 
 # Build the sample json given the job id
+echo "Building sample json for job $JOBID"
 python -c "import json; json.dump(json.load(open('$WORKDIR/split_samples.json'))['$JOBID'], open('$WORKDIR/sample.json', 'w'), indent=4)"
 
+echo "Sample json:"
+cat $WORKDIR/sample.json
 declare -A ARGS
 for key in workflow output samplejson year campaign isSyst isArray noHist overwrite voms chunk skipbadfiles outputDir remoteRepo; do
     ARGS[$key]=$(jq -r ".$key" $WORKDIR/arguments.json)
 done
 
+echo "Arguments:"
 # Unparse arguments and send to runner.py
 OPTS="--wf ${ARGS[workflow]} --year ${ARGS[year]} --campaign ${ARGS[campaign]} --chunk ${ARGS[chunk]}"
 if [ "${ARGS[voms]}" != "null" ]; then
@@ -48,5 +55,14 @@ OPTS="$OPTS --executor iterative --overwrite --outputdir $3"
 echo "Now launching: python runner.py $OPTS"
 python runner.py $OPTS
 
+OUTPUT=$?
+
+if [ $OUTPUT -ne 0 ]; then
+    echo "Job $JOBID failed with exit code $OUTPUT"
+    exit $OUTPUT
+fi
+
+echo "Job $JOBID completed successfully"
 touch $WORKDIR/.success
 
+exit $OUTPUT

condor_lxplus/submitter.py

@@ -59,6 +59,11 @@ def get_condor_submitter_parser(parser):
         default=None,
         help="If specified, access BTVNanoCommsioning from a remote tarball (downloaded via https), instead of from a transferred sandbox",
     )
+    parser.add_argument(
+        "--max_retries",
+        default=3,
+        help="Maximum number of retries for failed jobs. Failed jobs are forced to another machine.",
+    )
     return parser
 
 
@@ -260,6 +265,15 @@ def get_main_parser():
     with open(os.path.join(job_dir, "jobnum_list.txt"), "w") as f:
         f.write("\n".join([str(i) for i in range(counter)]))
 
+    if args.max_retries > 0:
+        retry_str = f"""
+on_exit_remove      = (ExitBySignal == False) && (ExitCode == 0)
+max_retries         = {args.max_retries}
+requirements        = Machine =!= LastRemoteHost
+        """
+    else:
+        retry_str = ""
+
     ## store the jdl file
     jdl_template = """Universe   = vanilla
 Executable = {executable}
@@ -281,6 +295,8 @@ def get_main_parser():
 transfer_input_files    = {transfer_input_files}
 transfer_output_files   = .success
 
+{retry_str}
+
 Queue JOBNUM from {jobnum_file}
 """.format(
         executable=f"{base_dir}/condor_lxplus/execute.sh",
@@ -289,6 +305,7 @@ def get_main_parser():
         envpath=envpath,
         log_dir=f"{base_dir}/{job_dir}/log",
         transfer_input_files=f"{base_dir}/{job_dir}/arguments.json,{base_dir}/{job_dir}/split_samples.json,{base_dir}/{job_dir}/jobnum_list.txt",
+        retry_str=retry_str,
         jobnum_file=f"{base_dir}/{job_dir}/jobnum_list.txt",
     )
     with open(os.path.join(job_dir, "submit.jdl"), "w") as f:

runner.py

@@ -298,7 +298,8 @@ def debug_parser(parser):
         coffeaoutput = (
             f'{histoutdir}/hists_{args.workflow}_{(sample_json).rstrip(".json")}.coffea'
         )
-    os.system(f"mkdir -p {histoutdir}")
+    if not args.noHist:
+        os.system(f"mkdir -p {histoutdir}")
     # load dataset
     with open(args.samplejson) as f:
         sample_dict = json.load(f)

src/BTVNanoCommissioning/utils/array_writer.py

@@ -50,14 +50,16 @@ def array_writer(
     othersMC=["Pileup_nTrueInt", "Pileup_nPU"],  # other fields, for MC only
     empty=False,
 ):
-    if not isRealData and systname != ["nominal"]:
+    if not isRealData and not empty:
         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:]:
+            for syst in systname:
+                if syst == "nominal":
+                    continue
                 pruned_event[f"weight_syst_{syst}"] = weights.weight(modifier=syst)
 
     if empty:

src/BTVNanoCommissioning/utils/histogrammer.py

@@ -425,6 +425,8 @@ def histogrammer(events, workflow, year="2022", campaign="Summer22"):
             _hist_dict[f"dr_{i}jet"] = Hist.Hist(
                 syst_axis, flav_axis, dr_axis, Hist.storage.Weight()
             )
+    elif "qgtag" in workflow:
+        pass
 
     ### Common kinematic variables histogram creation
     if "Wc_sf" not in workflow:

src/BTVNanoCommissioning/utils/sample.py

@@ -285,4 +285,20 @@
             "TTto2L2Nu_TuneCP5Up_13p6TeV_powheg-pythia8",
         ],
     },
+    "qgtag_dijet": {
+        "data": ["ZeroBias", "JetMET0", "JetMET1"],
+        "MC": [
+            "QCD-4Jets_HT-40to70_TuneCP5_13p6TeV_madgraphMLM-pythia8",
+            "QCD-4Jets_HT-70to100_TuneCP5_13p6TeV_madgraphMLM-pythia8",
+            "QCD-4Jets_HT-100to200_TuneCP5_13p6TeV_madgraphMLM-pythia8",
+            "QCD-4Jets_HT-200to400_TuneCP5_13p6TeV_madgraphMLM-pythia8",
+            "QCD-4Jets_HT-400to600_TuneCP5_13p6TeV_madgraphMLM-pythia8",
+            "QCD-4Jets_HT-600to800_TuneCP5_13p6TeV_madgraphMLM-pythia8",
+            "QCD-4Jets_HT-800to1000_TuneCP5_13p6TeV_madgraphMLM-pythia8",
+            "QCD-4Jets_HT-1000to1200_TuneCP5_13p6TeV_madgraphMLM-pythia8",
+            "QCD-4Jets_HT-1200to1500_TuneCP5_13p6TeV_madgraphMLM-pythia8",
+            "QCD-4Jets_HT-1500to2000_TuneCP5_13p6TeV_madgraphMLM-pythia8",
+            "QCD-4Jets_HT-2000_TuneCP5_13p6TeV_madgraphMLM-pythia8",
+        ],
+    },
 }

src/BTVNanoCommissioning/utils/selection.py

@@ -37,7 +37,7 @@ def jet_id(events, campaign, max_eta=2.5, min_pt=20):
                 (events.Jet.jetId >= 2) & (events.Jet.neHEF < 0.99),
                 ak.where(
                     (abs(events.Jet.eta) > 3.0),
-                    (events.Jet.jetId & (1 << 1)) & (events.Jet.neEmEF < 0.4),
+                    (events.Jet.jetId >= 2) & (events.Jet.neEmEF < 0.4),
                     ak.zeros_like(events.Jet.pt, dtype=bool),
                 ),
             ),
@@ -86,7 +86,7 @@ def jet_id(events, campaign, max_eta=2.5, min_pt=20):
             & ((events.Jet.pt > 50) | (events.Jet.puId >= 7))
         )
     else:
-        jetmask = (events.Jet.pt > min_pt) & (abs(events.Jet.eta) <= max_eta) & (jetid)
+        jetmask = (events.Jet.pt > min_pt) & (abs(events.Jet.eta) <= max_eta)# & (jetid)
     return jetmask
 
 

src/BTVNanoCommissioning/workflows/__init__.py

@@ -43,6 +43,11 @@
     NanoProcessor as BTA_ttbar_processor,
 )  # ttbar -kinFit
 
+## QG - Dijet producer
+from BTVNanoCommissioning.workflows.qgtag_dijet_producer import (
+    NanoProcessor as QGtagDijetProcessor,
+)
+
 # from BTVNanoCommissioning.workflows.example import (
 #     NanoProcessor as ExampleProcessor,
 # )
@@ -98,6 +103,12 @@
 # DY
 workflows["ctag_DY_sf"] = partial(CTAGDYValidSFProcessor, selectionModifier="DYM")
 workflows["ectag_DY_sf"] = partial(CTAGDYValidSFProcessor, selectionModifier="DYE")
+workflows["qgtag_DY_sf"] = partial(CTAGDYValidSFProcessor, selectionModifier="QG")
+
+# QG
+workflows["qgtag_dijet"] = partial(QGtagDijetProcessor, selectionModifier="DiPFJetAve")
+workflows["qgtag_dijet_zb"] = partial(QGtagDijetProcessor, selectionModifier="ZB")
+workflows["qgtag_dijet_pfjet"] = partial(QGtagDijetProcessor, selectionModifier="PFJet")
 
 # Tutorial
 # workflows["example"] = ExampleProcessor

src/BTVNanoCommissioning/workflows/ctag_DY_valid_sf.py

@@ -68,7 +68,7 @@ def process_shift(self, events, shift_name):
 
         isMu = False
         isEle = False
-        if "DYM" in self.selMod:
+        if "DYM" in self.selMod or "QG" in self.selMod:
             triggers = ["Mu17_TrkIsoVVL_Mu8_TrkIsoVVL_DZ_Mass8"]
             isMu = True
         elif "DYE" in self.selMod:
@@ -77,7 +77,7 @@ def process_shift(self, events, shift_name):
         else:
             raise ValueError(self.selMod, "is not a valid selection modifier.")
 
-        histname = {"DYM": "ctag_DY_sf", "DYE": "ectag_DY_sf"}
+        histname = {"DYM": "ctag_DY_sf", "DYE": "ectag_DY_sf", "QG": "qgtag_DY_sf"}
         output = {} if self.noHist else histogrammer(events, histname[self.selMod])
 
         if isRealData:
@@ -128,17 +128,6 @@ def process_shift(self, events, shift_name):
             axis=-1,
         )
 
-        pl_iso = ak.all(
-            events.Jet.metric_table(pos_dilep) > 0.4, axis=2, mask_identity=True
-        )
-        nl_iso = ak.all(
-            events.Jet.metric_table(neg_dilep) > 0.4, axis=2, mask_identity=True
-        )
-        jet_sel = ak.fill_none(
-            jet_id(events, self._campaign) & pl_iso & nl_iso,
-            False,
-            axis=-1,
-        )
 
         pos_dilep = ak.pad_none(pos_dilep, 1, axis=1)
         neg_dilep = ak.pad_none(neg_dilep, 1, axis=1)
@@ -155,15 +144,21 @@ def process_shift(self, events, shift_name):
         nl_iso = ak.all(
             events.Jet.metric_table(neg_dilep[:, 0]) > 0.4, axis=2, mask_identity=True
         )
-        event_jet = events.Jet[
-            ak.fill_none(
-                jet_id(events, self._campaign) & pl_iso & nl_iso,
-                False,
-                axis=-1,
-            )
-        ]
+
+        if "QG" in self.selMod:
+            jetmask = jet_id(events, self._campaign, max_eta=5.13)
+        else:
+            jetmask = jet_id(events, self._campaign)
+
+        jet_sel = ak.fill_none(
+            pl_iso & nl_iso & jetmask,
+            False,
+            axis=-1,
+        )
+
+        event_jet = events.Jet[jet_sel]
+
         req_jets = ak.count(event_jet.pt, axis=1) >= 1
-        # event_jet = ak.pad_none(event_jet, 1, axis=1)
 
         # store jet index for PFCands, create mask on the jet index
         jetindx = ak.mask(
@@ -173,8 +168,25 @@ def process_shift(self, events, shift_name):
         jetindx = ak.pad_none(jetindx, 1)
         jetindx = jetindx[:, 0]
 
+        selection = (
+            req_lumi & req_trig & req_dilep & req_dilepmass & req_jets & req_metfilter
+        )
+
+        if "QG" in self.selMod:
+            temp_jet = ak.pad_none(event_jet, 1, axis=1)
+
+            req_lead_jet = ak.fill_none(
+                (
+                    np.abs(temp_jet[:, 0].delta_phi(pos_dilep[:, 0] + neg_dilep[:, 0]))
+                    > 2.7
+                ),
+                False,
+                axis=-1,
+            )
+            selection = selection & req_lead_jet
+
         event_level = ak.fill_none(
-            req_lumi & req_trig & req_dilep & req_dilepmass & req_jets & req_metfilter,
+            selection,
             False,
         )
         if len(events[event_level]) == 0:
@@ -208,7 +220,11 @@ 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]
+        if self.selMod == "QG":
+            pruned_ev["SelJet"] = event_jet[event_level][:, 0]
+        else:
+            pruned_ev["SelJet"] = event_jet[event_level]
+
         if isMu:
             pruned_ev["MuonPlus"] = sposmu
             pruned_ev["MuonMinus"] = snegmu
@@ -233,7 +249,7 @@ def process_shift(self, events, shift_name):
         pruned_ev["dr_mu1jet"] = sposmu.delta_r(sel_jet)
         pruned_ev["dr_mu2jet"] = snegmu.delta_r(sel_jet)
         # Find the PFCands associate with selected jets. Search from jetindex->JetPFCands->PFCand
-        if "PFCands" in events.fields:
+        if "PFCands" in events.fields and "QG" not in self.selMod:
             pruned_ev["PFCands"] = PFCand_link(events, event_level, jetindx)
 
         ####################
@@ -254,8 +270,40 @@ def process_shift(self, events, shift_name):
             )
         # Output arrays
         if self.isArray:
+            if "QG" in self.selMod:
+                othersData = [
+                    "SV_*",
+                    "PV_npvs",
+                    "PV_npvsGood",
+                    "Rho_*",
+                    "SoftMuon_dxySig",
+                    "Muon_sip3d",
+                    "run",
+                    "luminosityBlock",
+                ]
+                for trigger in triggers:
+                    othersData.append(f"HLT_{trigger}")
+            else:
+                othersData = [
+                    "PFCands_*",
+                    "MuonJet_*",
+                    "SV_*",
+                    "PV_npvs",
+                    "PV_npvsGood",
+                    "Rho_*",
+                    "SoftMuon_dxySig",
+                    "Muon_sip3d",
+                ]
+
             array_writer(
-                self, pruned_ev, events, weights, systematics, dataset, isRealData
+                self,
+                pruned_ev,
+                events,
+                weights,
+                systematics,
+                dataset,
+                isRealData,
+                othersData=othersData,
             )
 
         return {dataset: output}