From 758573cb4b47f3a8f2457966dbc9018dc48ed9d4 Mon Sep 17 00:00:00 2001 From: BrieucF Date: Fri, 5 Apr 2024 10:19:48 +0200 Subject: [PATCH] [Full Sim] Add a C++ higgs recoil mass plotter --- .../FCCeeGeneralOverview.md | 57 +++++++++++++++++-- .../FCCeeGeneralOverview/plot_recoil_mass.C | 31 ++++++++++ .../FCCeeGeneralOverview/plot_recoil_mass.py | 6 +- 3 files changed, 87 insertions(+), 7 deletions(-) create mode 100644 full-detector-simulations/FCCeeGeneralOverview/plot_recoil_mass.C diff --git a/full-detector-simulations/FCCeeGeneralOverview/FCCeeGeneralOverview.md b/full-detector-simulations/FCCeeGeneralOverview/FCCeeGeneralOverview.md index 5829a485..0b5d402c 100644 --- a/full-detector-simulations/FCCeeGeneralOverview/FCCeeGeneralOverview.md +++ b/full-detector-simulations/FCCeeGeneralOverview/FCCeeGeneralOverview.md @@ -13,8 +13,8 @@ So, let's start playing with Full Sim! ```bash # connect to a machine with cvmfs access and running an OS supported by Key4hep (Alma9 here) ssh -X username@submit-test.mit.edu # or ssh -X username@lxplus.cern.ch -# set-up the Key4hep environment -# make sure you are in bash or zsh shell +# set-up the Key4hep environment, using the nightlies since we need the latest and greatest version of the packages +# (make sure you are in bash or zsh shell) source /cvmfs/sw-nightlies.hsf.org/key4hep/setup.sh # create a repository for the tutorial mkdir FCC_Full_Sim_Tutorial @@ -75,11 +75,12 @@ Let's now use the reconstructed sample to produce some physics quantities. As an ``` +import sys from podio import root_io import ROOT ROOT.gROOT.SetBatch(True) -input_file_path = "wzp6_ee_mumuH_ecm240_CLD_RECO.root" +input_file_path = sys.argv[1] podio_reader = root_io.Reader(input_file_path) th1_recoil = ROOT.TH1F("Recoil Mass", "Recoil Mass", 100, 110, 160) @@ -101,14 +102,62 @@ th1_recoil.Draw() canvas_recoil.Print("recoil_mass.png") ``` -Let's run it on a sample with slightly more stat: +The equivalent C++ ROOT macro looks like this: +``` +#include "podio/ROOTFrameReader.h" +#include "podio/Frame.h" + +#include "edm4hep/ReconstructedParticleCollection.h" + +int plot_recoil_mass(std::string input_file_path) { + auto reader = podio::ROOTFrameReader(); + reader.openFile(input_file_path); + + TH1* th1_recoil = new TH1F("Recoil Mass", "Recoil Mass", 100, 110., 160.); + + ROOT::Math::LorentzVector> p_cm(0., 0., 0., 240.); + for (size_t i = 0; i < reader.getEntries("events"); ++i) { + auto event = podio::Frame(reader.readNextEntry("events")); + auto& pfos = event.get("TightSelectedPandoraPFOs"); + int n_good_muons = 0; + ROOT::Math::LorentzVector> p_mumu; + for (const auto& pfo : pfos) { + if (std::abs(pfo.getPDG()) == 13 and pfo.getEnergy() > 20.) { + n_good_muons++; + p_mumu += ROOT::Math::LorentzVector>(pfo.getMomentum().x, pfo.getMomentum().y, pfo.getMomentum().z, pfo.getEnergy()); + } + } + if(n_good_muons == 2) + th1_recoil->Fill((p_cm - p_mumu).M()); + } + TCanvas* canvas_recoil = new TCanvas("Recoil Mass", "Recoil Mass"); + th1_recoil->Draw(); + canvas_recoil->Print("recoil_mass.png"); + return 0; +} +``` + +Let's get a similar sample but with more stat: ``` cd ../../fcc-tutorials/full-detector-simulations/FCCeeGeneralOverview/ wget https://fccsw.web.cern.ch/fccsw/tutorials/MIT2024/wzp6_ee_mumuH_ecm240_CLD_RECO_moreStat.root +``` + +and produce the simple Higgs recoil mass plot in Python: +``` python plot_recoil_mass.py wzp6_ee_mumuH_ecm240_CLD_RECO_moreStat.root display recoil_mass.png ``` +or in C++ with the ROOT interpreter: +``` +root -b +.L plot_recoil_mass.C +plot_recoil_mass("wzp6_ee_mumuH_ecm240_CLD_RECO_moreStat.root") +.q +display recoil_mass.png +``` + This illustrates how easy it is already to do physics with Full Sim. Of course, if we had to do a realistic analysis, we would run on more events, properly select muons from the Z, include backgrounds, ..., and we would therefore use FCCAnalyses or plain C++ but it is not the topic of this tutorial. If you want to go further, the following [Doxygen page](https://edm4hep.web.cern.ch/classedm4hep_1_1_reconstructed_particle-members.html) will help you in understanding what members can be called on a given edm4hep object. ## Towards detector optimization with Full Sim: ALLEGRO ECAL diff --git a/full-detector-simulations/FCCeeGeneralOverview/plot_recoil_mass.C b/full-detector-simulations/FCCeeGeneralOverview/plot_recoil_mass.C new file mode 100644 index 00000000..a0e170c0 --- /dev/null +++ b/full-detector-simulations/FCCeeGeneralOverview/plot_recoil_mass.C @@ -0,0 +1,31 @@ +#include "podio/ROOTFrameReader.h" +#include "podio/Frame.h" + +#include "edm4hep/ReconstructedParticleCollection.h" + +int plot_recoil_mass(std::string input_file_path) { + auto reader = podio::ROOTFrameReader(); + reader.openFile(input_file_path); + + TH1* th1_recoil = new TH1F("Recoil Mass", "Recoil Mass", 100, 110., 160.); + + ROOT::Math::LorentzVector> p_cm(0., 0., 0., 240.); + for (size_t i = 0; i < reader.getEntries("events"); ++i) { + auto event = podio::Frame(reader.readNextEntry("events")); + auto& pfos = event.get("TightSelectedPandoraPFOs"); + int n_good_muons = 0; + ROOT::Math::LorentzVector> p_mumu; + for (const auto& pfo : pfos) { + if (std::abs(pfo.getPDG()) == 13 and pfo.getEnergy() > 20.) { + n_good_muons++; + p_mumu += ROOT::Math::LorentzVector>(pfo.getMomentum().x, pfo.getMomentum().y, pfo.getMomentum().z, pfo.getEnergy()); + } + } + if(n_good_muons == 2) + th1_recoil->Fill((p_cm - p_mumu).M()); + } + TCanvas* canvas_recoil = new TCanvas("Recoil Mass", "Recoil Mass"); + th1_recoil->Draw(); + canvas_recoil->Print("recoil_mass.png"); + return 0; +} diff --git a/full-detector-simulations/FCCeeGeneralOverview/plot_recoil_mass.py b/full-detector-simulations/FCCeeGeneralOverview/plot_recoil_mass.py index 782f771f..e1f6eaa2 100644 --- a/full-detector-simulations/FCCeeGeneralOverview/plot_recoil_mass.py +++ b/full-detector-simulations/FCCeeGeneralOverview/plot_recoil_mass.py @@ -6,15 +6,15 @@ input_file_path = sys.argv[1] podio_reader = root_io.Reader(input_file_path) -th1_recoil = ROOT.TH1F("Recoil Mass", "Recoil Mass", 100, 110, 160) +th1_recoil = ROOT.TH1F("Recoil Mass", "Recoil Mass", 100, 110., 160.) -p_cm = ROOT.Math.LorentzVector('ROOT::Math::PxPyPzE4D')(0, 0, 0, 240) +p_cm = ROOT.Math.LorentzVector('ROOT::Math::PxPyPzE4D')(0., 0., 0., 240.) for event in podio_reader.get("events"): pfos = event.get("TightSelectedPandoraPFOs") n_good_muons = 0 p_mumu = ROOT.Math.LorentzVector('ROOT::Math::PxPyPzE4D')() for pfo in pfos: - if abs(pfo.getPDG()) == 13 and pfo.getEnergy() > 20: + if abs(pfo.getPDG()) == 13 and pfo.getEnergy() > 20.: n_good_muons += 1 p_mumu += ROOT.Math.LorentzVector('ROOT::Math::PxPyPzE4D')(pfo.getMomentum().x, pfo.getMomentum().y, pfo.getMomentum().z, pfo.getEnergy()) if n_good_muons == 2: