Skip to content

Commit

Permalink
fix met in notebook (#347)
Browse files Browse the repository at this point in the history
* fix met in notebook
  • Loading branch information
jpata authored Oct 5, 2024
1 parent b77be93 commit 3200518
Show file tree
Hide file tree
Showing 6 changed files with 165 additions and 47 deletions.
6 changes: 3 additions & 3 deletions mlpf/data_cms/postprocessing2.py
Original file line number Diff line number Diff line change
Expand Up @@ -165,9 +165,9 @@ def get_charge(pid):

def compute_gen_met(g):
genpart = [elem for elem in g.nodes if elem[0] == "cp"]
px = np.sum([g.nodes[elem]["pt"] * np.cos(g.nodes[elem]["phi"]) for elem in genpart])
py = np.sum([g.nodes[elem]["pt"] * np.sin(g.nodes[elem]["phi"]) for elem in genpart])
met = np.sqrt(px**2 + py**2)
sum_px = np.sum([g.nodes[elem]["pt"] * np.cos(g.nodes[elem]["phi"]) for elem in genpart])
sum_py = np.sum([g.nodes[elem]["pt"] * np.sin(g.nodes[elem]["phi"]) for elem in genpart])
met = np.sqrt(sum_px**2 + sum_py**2)
return met


Expand Down
20 changes: 10 additions & 10 deletions mlpf/data_cms/prepare_args.py
Original file line number Diff line number Diff line change
Expand Up @@ -15,20 +15,20 @@
# ("SinglePiMinusFlatPt0p7To1000_cfi", 800000, 801010, "genjob_pu55to75.sh", outdir + "/pu55to75"),

# ("TTbar_14TeV_TuneCUETP8M1_cfi", 702000, 720000, "genjob_nopu.sh", outdir + "/nopu"),
# ("MultiParticlePFGun50_cfi", 800000, 805000, "genjob_nopu.sh", outdir + "/nopu"),
# ("VBF_TuneCP5_14TeV_pythia8_cfi", 900000, 920010, "genjob_nopu.sh", outdir + "/nopu"),
# ("QCDForPF_14TeV_TuneCUETP8M1_cfi", 1000000,1020010, "genjob_nopu.sh", outdir + "/nopu"),
# ("ZTT_All_hadronic_14TeV_TuneCUETP8M1_cfi", 1100000,1100010, "genjob_nopu.sh", outdir + "/nopu"),

# ("SingleElectronFlatPt1To1000_pythia8_cfi", 900000, 901010, "genjob_nopu.sh", outdir + "/nopu"),
# ("SingleGammaFlatPt1To1000_pythia8_cfi", 1000000,1001010, "genjob_nopu.sh", outdir + "/nopu"),
# ("SingleMuFlatPt1To1000_pythia8_cfi", 1100000,1101010, "genjob_nopu.sh", outdir + "/nopu"),
# ("SingleNeutronFlatPt0p7To1000_cfi", 1200000,1200110, "genjob_nopu.sh", outdir + "/nopu"),
# ("SinglePi0Pt1To1000_pythia8_cfi", 1300000,1301010, "genjob_nopu.sh", outdir + "/nopu"),
# ("SinglePiMinusFlatPt0p7To1000_cfi", 1400000,1401010, "genjob_nopu.sh", outdir + "/nopu"),
# ("SingleProtonMinusFlatPt0p7To1000_cfi", 1500000,1500010, "genjob_nopu.sh", outdir + "/nopu"),
# ("SingleTauFlatPt1To1000_cfi", 1600000,1600110, "genjob_nopu.sh", outdir + "/nopu"),
# ("SingleK0FlatPt1To1000_pythia8_cfi", 1700000,1701010, "genjob_nopu.sh", outdir + "/nopu"),
# ("MultiParticlePFGun50_cfi", 800000, 801000, "genjob_nopu.sh", outdir + "/nopu"),
# ("SingleElectronFlatPt1To1000_pythia8_cfi", 900000, 905010, "genjob_nopu.sh", outdir + "/nopu"),
# ("SingleGammaFlatPt1To1000_pythia8_cfi", 1000000,1005010, "genjob_nopu.sh", outdir + "/nopu"),
# ("SingleMuFlatPt1To1000_pythia8_cfi", 1100000,1105010, "genjob_nopu.sh", outdir + "/nopu"),
# ("SingleNeutronFlatPt0p7To1000_cfi", 1200000,1205010, "genjob_nopu.sh", outdir + "/nopu"),
# ("SinglePi0Pt1To1000_pythia8_cfi", 1300000,1305010, "genjob_nopu.sh", outdir + "/nopu"),
# ("SinglePiMinusFlatPt0p7To1000_cfi", 1400000,1405010, "genjob_nopu.sh", outdir + "/nopu"),
# ("SingleProtonMinusFlatPt0p7To1000_cfi", 1500000,1505010, "genjob_nopu.sh", outdir + "/nopu"),
# ("SingleTauFlatPt1To1000_cfi", 1600000,1605010, "genjob_nopu.sh", outdir + "/nopu"),
# ("SingleK0FlatPt1To1000_pythia8_cfi", 1700000,1705010, "genjob_nopu.sh", outdir + "/nopu"),
]

if __name__ == "__main__":
Expand Down
2 changes: 1 addition & 1 deletion notebooks/cms/cms-validate-postprocessing.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -479,7 +479,7 @@
"def met(pt, phi):\n",
" px = pt * np.cos(phi)\n",
" py = pt * np.sin(phi)\n",
" pt = np.sqrt(awkward.sum(px**2 + py**2, axis=1))\n",
" pt = np.sqrt(awkward.sum(px, axis=1)**2 + awkward.sum(py, axis=1)**2)\n",
" return pt"
]
},
Expand Down
172 changes: 145 additions & 27 deletions notebooks/cms/cms-validate-root.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,16 @@
"from plot_utils import pid_to_text"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "99f9bec3-cfe9-4395-85d4-67bda0f501f9",
"metadata": {},
"outputs": [],
"source": [
"# uproot.open(\"/local/joosep/mlpf/cms/20240823_simcluster/nopu/QCDForPF_14TeV_TuneCUETP8M1_cfi/root/pfntuple_1000000.root\")[\"pfana/pftree\"].keys()"
]
},
{
"cell_type": "code",
"execution_count": null,
Expand All @@ -46,7 +56,7 @@
"outputs": [],
"source": [
"def load_tree(ttree):\n",
" particles_pythia = ttree.arrays([\"gen_pt\", \"gen_eta\", \"gen_phi\", \"gen_energy\", \"gen_pdgid\", \"gen_status\"])\n",
" particles_pythia = ttree.arrays([\"gen_pt\", \"gen_eta\", \"gen_phi\", \"gen_energy\", \"gen_pdgid\", \"gen_status\", \"gen_daughters\"])\n",
" particles_cp = ttree.arrays([\"caloparticle_pt\", \"caloparticle_eta\", \"caloparticle_phi\", \"caloparticle_energy\", \"caloparticle_pid\"])\n",
" genjet = ttree.arrays([\"genjet_pt\", \"genjet_eta\", \"genjet_phi\", \"genjet_energy\"])\n",
" genmet = ttree.arrays([\"genmet_pt\"])\n",
Expand All @@ -60,10 +70,10 @@
"metadata": {},
"outputs": [],
"source": [
"#https://jpata.web.cern.ch/jpata/mlpf/cms/20240823_simcluster/nopu/TTbar_14TeV_TuneCUETP8M1_cfi/root\n",
"#Download from https://jpata.web.cern.ch/jpata/mlpf/cms/20240823_simcluster/nopu/\n",
"\n",
"tts = [\n",
" load_tree(uproot.open(fn)[\"pfana/pftree\"]) for fn in glob.glob(\"/local/joosep/mlpf/cms/20240823_simcluster/nopu/TTbar_14TeV_TuneCUETP8M1_cfi/root/pfntuple_7000*.root\")\n",
" load_tree(uproot.open(fn)[\"pfana/pftree\"]) for fn in glob.glob(\"/local/joosep/mlpf/cms/20240823_simcluster/nopu/QCDForPF_14TeV_TuneCUETP8M1_cfi/root/pfntuple_100000*.root\")\n",
"]\n",
"tts = awkward.concatenate(tts, axis=0)"
]
Expand Down Expand Up @@ -101,7 +111,7 @@
"metadata": {},
"outputs": [],
"source": [
"b = np.logspace(-2,3,100)\n",
"b = np.logspace(-2,4,100)\n",
"plt.figure(figsize=(5,5))\n",
"\n",
"abs_pid = np.abs(particles_pythia[\"gen_pdgid\"])\n",
Expand Down Expand Up @@ -144,17 +154,36 @@
"source": [
"#Sum pT in event\n",
"b = np.linspace(0,1000,100)\n",
"plt.figure(figsize=(5,5))\n",
"fig, axs = plt.subplots(1,2, figsize=(10,5))\n",
"\n",
"\n",
"plt.sca(axs[0])\n",
"\n",
"plt.title(\"sum $p_T$ in event\")\n",
"plt.hist2d(\n",
" awkward.to_numpy(awkward.sum(particles_pythia[mask_pythia_nonu][\"gen_pt\"], axis=1)),\n",
" awkward.to_numpy(awkward.sum(particles_cp[mask_cp][\"caloparticle_pt\"], axis=1)),\n",
" bins=(b,b), cmap=\"hot_r\", #norm=matplotlib.colors.LogNorm()\n",
")\n",
"plt.plot([0,2000],[0,2000], color=\"black\")\n",
"plt.plot([0,2000],[0,2000], color=\"black\", lw=0.5)\n",
"plt.xlabel(\"Pythia, no nu\")\n",
"plt.ylabel(\"CaloParticle\")\n",
"\n",
"plt.sca(axs[1])\n",
"\n",
"#it seems like caloparticles and p\n",
"pythia_ptcut = (particles_pythia[\"gen_pt\"]>0.25)\n",
"\n",
"plt.title(\"$p_T>0.25$ GeV\")\n",
"plt.hist2d(\n",
" awkward.to_numpy(awkward.sum(particles_pythia[mask_pythia_nonu & pythia_ptcut][\"gen_pt\"], axis=1)),\n",
" awkward.to_numpy(awkward.sum(particles_cp[mask_cp][\"caloparticle_pt\"], axis=1)),\n",
" bins=(b,b), cmap=\"hot_r\", #norm=matplotlib.colors.LogNorm()\n",
")\n",
"plt.plot([0,2000],[0,2000], color=\"black\", lw=0.5)\n",
"plt.xlabel(\"Pythia, no nu\")\n",
"plt.ylabel(\"CaloParticle\")"
"plt.ylabel(\"CaloParticle\")\n",
"plt.tight_layout()\n",
"plt.suptitle(\"sum $p_T$\")"
]
},
{
Expand All @@ -174,6 +203,7 @@
"source": [
"jets_coll = {}\n",
"jets_coll[\"genjet\"] = genjets\n",
"jet_constituents = {}\n",
"\n",
"vec = vector.awk(\n",
" awkward.zip(\n",
Expand All @@ -188,6 +218,22 @@
"jetdef = fastjet.JetDefinition(fastjet.antikt_algorithm, 0.4)\n",
"cluster = fastjet.ClusterSequence(vec.to_xyzt(), jetdef)\n",
"jets_coll[\"pythia_nonu\"] = cluster.inclusive_jets(min_pt=3)\n",
"jet_constituents[\"pythia_nonu\"] = cluster.constituent_index(min_pt=3)\n",
"\n",
"vec = vector.awk(\n",
" awkward.zip(\n",
" { \n",
" \"pt\": particles_pythia[mask_pythia_nonu & pythia_ptcut][\"gen_pt\"],\n",
" \"eta\": particles_pythia[mask_pythia_nonu & pythia_ptcut][\"gen_eta\"],\n",
" \"phi\": particles_pythia[mask_pythia_nonu & pythia_ptcut][\"gen_phi\"],\n",
" \"energy\": particles_pythia[mask_pythia_nonu & pythia_ptcut][\"gen_energy\"],\n",
" }\n",
" )\n",
")\n",
"jetdef = fastjet.JetDefinition(fastjet.antikt_algorithm, 0.4)\n",
"cluster = fastjet.ClusterSequence(vec.to_xyzt(), jetdef)\n",
"jets_coll[\"pythia_nonu_ptcut\"] = cluster.inclusive_jets(min_pt=3)\n",
"jet_constituents[\"pythia_nonu_ptcut\"] = cluster.constituent_index(min_pt=3)\n",
"\n",
"vec = vector.awk(\n",
" awkward.zip(\n",
Expand Down Expand Up @@ -216,11 +262,13 @@
"jetdef = fastjet.JetDefinition(fastjet.antikt_algorithm, 0.4)\n",
"cluster = fastjet.ClusterSequence(vec.to_xyzt(), jetdef)\n",
"jets_coll[\"cp\"] = cluster.inclusive_jets(min_pt=3)\n",
"jet_constituents[\"cp\"] = cluster.constituent_index(min_pt=3)\n",
"\n",
"genjet_to_pythia = jet_utils.match_two_jet_collections(jets_coll, \"genjet\", \"pythia\", 0.1)\n",
"genjet_to_cp = jet_utils.match_two_jet_collections(jets_coll, \"genjet\", \"cp\", 0.1)\n",
"pythia_to_cp = jet_utils.match_two_jet_collections(jets_coll, \"pythia\", \"cp\", 0.1)\n",
"pythia_nonu_to_cp = jet_utils.match_two_jet_collections(jets_coll, \"pythia_nonu\", \"cp\", 0.1)"
"pythia_nonu_to_cp = jet_utils.match_two_jet_collections(jets_coll, \"pythia_nonu\", \"cp\", 0.1)\n",
"pythia_nonu_ptcut_to_cp = jet_utils.match_two_jet_collections(jets_coll, \"pythia_nonu_ptcut\", \"cp\", 0.1)"
]
},
{
Expand All @@ -230,11 +278,12 @@
"metadata": {},
"outputs": [],
"source": [
"plt.figure(figsize=(5,5))\n",
"b = np.logspace(0,3,100)\n",
"b = np.logspace(0,4,100)\n",
"plt.hist(awkward.flatten(jets_coll[\"genjet\"].pt), bins=b, histtype=\"step\", lw=1, label=\"genjet\")\n",
"plt.hist(awkward.flatten(jets_coll[\"pythia\"].pt), bins=b, histtype=\"step\", lw=1, label=\"Pythia\")\n",
"plt.hist(awkward.flatten(jets_coll[\"cp\"].pt), bins=b, histtype=\"step\", lw=1, label=\"CaloParticle\")\n",
"plt.hist(awkward.flatten(jets_coll[\"pythia\"].pt), bins=b, histtype=\"step\", lw=1, label=\"Pythia\")\n",
"plt.hist(awkward.flatten(jets_coll[\"pythia_nonu\"].pt), bins=b, histtype=\"step\", lw=1, label=\"Pythia, NoNu\")\n",
"plt.hist(awkward.flatten(jets_coll[\"pythia_nonu_ptcut\"].pt), bins=b, histtype=\"step\", lw=1, label=\"Pythia, NoNu, pt>0.25 GeV\")\n",
"plt.xscale(\"log\")\n",
"plt.xlabel(\"jet $p_T$ [GeV]\")\n",
"plt.ylabel(\"number of jets\")\n",
Expand Down Expand Up @@ -282,20 +331,70 @@
"\n",
"plt.hist(\n",
" awkward.flatten(\n",
" (jets_coll[\"cp\"][pythia_to_cp[\"cp\"]].pt / jets_coll[\"pythia\"][pythia_to_cp[\"pythia\"]].pt)\n",
" ), bins=b, histtype=\"step\", lw=1, label=\"CaloParticle vs Pythia\"\n",
");\n",
"\n",
"plt.hist(\n",
" awkward.flatten(\n",
" (jets_coll[\"cp\"][pythia_nonu_to_cp[\"cp\"]].pt / jets_coll[\"pythia_nonu\"][pythia_nonu_to_cp[\"pythia_nonu\"]].pt)\n",
" ), bins=b, histtype=\"step\", lw=1, label=\"CaloParticle vs. Pythia\"\n",
" ), bins=b, histtype=\"step\", lw=1, label=\"CaloParticle vs. Pythia NoNu\"\n",
");\n",
"\n",
"plt.hist(\n",
" awkward.flatten(\n",
" (jets_coll[\"cp\"][pythia_to_cp[\"cp\"]].pt / jets_coll[\"pythia\"][pythia_to_cp[\"pythia\"]].pt)\n",
" ), bins=b, histtype=\"step\", lw=1, label=\"CaloParticle vs Pythia NoNu\"\n",
" (jets_coll[\"cp\"][pythia_nonu_ptcut_to_cp[\"cp\"]].pt / jets_coll[\"pythia_nonu_ptcut\"][pythia_nonu_ptcut_to_cp[\"pythia_nonu_ptcut\"]].pt)\n",
" ), bins=b, histtype=\"step\", lw=1, label=\"CaloParticle vs. Pythia NoNu, pt>0.25GeV\"\n",
");\n",
"\n",
"\n",
"plt.xscale(\"log\")\n",
"plt.yscale(\"log\")\n",
"plt.xlabel(\"jet $p_T$ / Pythia jet $p_T$\")\n",
"plt.legend(loc=1, fontsize=12)\n",
"plt.ylim(1,1e6)\n",
"plt.axvline(1.0, color=\"black\", ls=\"--\", lw=0.5)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4358b911-705f-4192-b696-74d6ac0afdf5",
"metadata": {},
"outputs": [],
"source": [
"plt.figure(figsize=(5,5))\n",
"b = np.linspace(0.9,1.1,600)\n",
"\n",
"pt_pythia = jets_coll[\"pythia\"][pythia_to_cp[\"pythia\"]].pt\n",
"pt_cp = jets_coll[\"cp\"][pythia_to_cp[\"cp\"]].pt\n",
"plt.hist(\n",
" awkward.flatten(\n",
" (pt_cp / pt_pythia)[pt_pythia>50]\n",
" ), bins=b, histtype=\"step\", lw=1, label=\"CaloParticle vs Pythia\"\n",
");\n",
"\n",
"pt_pythia = jets_coll[\"pythia_nonu\"][pythia_nonu_to_cp[\"pythia_nonu\"]].pt\n",
"pt_cp = jets_coll[\"cp\"][pythia_nonu_to_cp[\"cp\"]].pt\n",
"plt.hist(\n",
" awkward.flatten(\n",
" (pt_cp / pt_pythia)[pt_pythia>50]\n",
" ), bins=b, histtype=\"step\", lw=1, label=\"CaloParticle vs. Pythia NoNu\"\n",
");\n",
"\n",
"pt_pythia = jets_coll[\"pythia_nonu_ptcut\"][pythia_nonu_ptcut_to_cp[\"pythia_nonu_ptcut\"]].pt\n",
"pt_cp = jets_coll[\"cp\"][pythia_nonu_ptcut_to_cp[\"cp\"]].pt\n",
"plt.hist(\n",
" awkward.flatten(\n",
" (pt_cp / pt_pythia)[pt_pythia>50]\n",
" ), bins=b, histtype=\"step\", lw=1, label=\"CaloParticle vs. Pythia NoNu, pt>0.25GeV\"\n",
");\n",
"\n",
"#plt.xscale(\"log\")\n",
"plt.yscale(\"log\")\n",
"plt.xlabel(\"jet $p_T$ / Pythia jet $p_T$\")\n",
"plt.legend(loc=1, fontsize=12)\n",
"plt.ylim(1,1e5)\n",
"plt.axvline(1.0, color=\"black\", ls=\"--\", lw=0.5)"
]
},
Expand All @@ -317,10 +416,11 @@
"def met(pt, phi):\n",
" px = pt * np.cos(phi)\n",
" py = pt * np.sin(phi)\n",
" pt = np.sqrt(awkward.sum(px**2 + py**2, axis=1))\n",
" pt = np.sqrt(awkward.sum(px, axis=1)**2 + awkward.sum(py, axis=1)**2)\n",
" return pt\n",
"\n",
"met_pythia = met(particles_pythia[mask_pythia_nonu][\"gen_pt\"], particles_pythia[mask_pythia_nonu][\"gen_phi\"])\n",
"met_pythia_nonu = met(particles_pythia[mask_pythia_nonu][\"gen_pt\"], particles_pythia[mask_pythia_nonu][\"gen_phi\"])\n",
"met_pythia_nonu_ptcut = met(particles_pythia[mask_pythia_nonu & pythia_ptcut][\"gen_pt\"], particles_pythia[mask_pythia_nonu & pythia_ptcut][\"gen_phi\"])\n",
"met_cp = met(particles_cp[mask_cp][\"caloparticle_pt\"], particles_cp[mask_cp][\"caloparticle_phi\"])"
]
},
Expand All @@ -331,29 +431,47 @@
"metadata": {},
"outputs": [],
"source": [
"b = np.logspace(0,3,100)\n",
"b = np.linspace(0,100,21)\n",
"plt.figure(figsize=(5,5))\n",
"plt.hist(genmet, bins=b, histtype=\"step\", label=\"GenMET\")\n",
"plt.hist(met_pythia, bins=b, histtype=\"step\", label=\"Pythia, no nu\")\n",
"plt.hist(genmet, bins=b, histtype=\"step\", label=\"genMetTrue\")\n",
"plt.hist(met_pythia_nonu, bins=b, histtype=\"step\", label=\"Pythia, no nu\")\n",
"plt.hist(met_cp, bins=b, histtype=\"step\", label=\"CaloParticle\")\n",
"plt.xscale(\"log\")\n",
"plt.legend(loc=2, fontsize=12)\n",
"#plt.xscale(\"log\")\n",
"plt.yscale(\"log\")\n",
"plt.legend(loc=1, fontsize=12)\n",
"plt.xlabel(\"MET [GeV]\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c104e1ce-7ff7-4107-8ff7-1e7b63df4d0e",
"id": "459d4032-58b6-4464-b396-7e800306fcdd",
"metadata": {},
"outputs": [],
"source": [
"plt.figure(figsize=(5,5))\n",
"plt.scatter(\n",
" met_pythia, met_cp\n",
")\n",
"plt.plot([0,500], [0,500], color=\"black\", ls=\"--\", marker=\".\")"
"fig, axs = plt.subplots(1,2, figsize=(10,5))\n",
"\n",
"plt.sca(axs[0])\n",
"plt.scatter(genmet, met_pythia_nonu)\n",
"plt.xlabel(\"genMetTrue\")\n",
"plt.ylabel(\"Pythia MET\")\n",
"\n",
"plt.sca(axs[1])\n",
"plt.scatter(met_pythia_nonu, met_cp)\n",
"plt.xlabel(\"Pythia MET\")\n",
"plt.ylabel(\"CaloParticle MET\")\n",
"\n",
"plt.tight_layout()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "00f745fa-0c10-4169-9e68-187350bea53f",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
Expand Down
6 changes: 3 additions & 3 deletions scripts/clic/postprocessing.py
Original file line number Diff line number Diff line change
Expand Up @@ -792,9 +792,9 @@ def get_p4(part, prefix="MCParticles"):


def compute_met(p4):
px = awkward.sum(p4.px, axis=1)
py = awkward.sum(p4.py, axis=1)
met = np.sqrt(px**2 + py**2)
sum_px = awkward.sum(p4.px, axis=1)
sum_py = awkward.sum(p4.py, axis=1)
met = np.sqrt(sum_px**2 + sum_py**2)
return met


Expand Down
6 changes: 3 additions & 3 deletions scripts/fccee_cld/postprocessing.py
Original file line number Diff line number Diff line change
Expand Up @@ -706,9 +706,9 @@ def get_p4(part, prefix="MCParticles"):

def compute_met(part, prefix="MCParticles"):
p4 = get_p4(part, prefix)
px = awkward.sum(p4.px, axis=1)
py = awkward.sum(p4.py, axis=1)
met = np.sqrt(px**2 + py**2)
sum_px = awkward.sum(p4.px, axis=1)
sum_py = awkward.sum(p4.py, axis=1)
met = np.sqrt(sum_px**2 + sum_py**2)
return met


Expand Down

0 comments on commit 3200518

Please sign in to comment.