From 909958029d67028530243379e8d8fa8dd042be7f Mon Sep 17 00:00:00 2001 From: Brieuc Francois Date: Thu, 21 Mar 2024 14:32:49 +0100 Subject: [PATCH] Mit tutorial continued (#143) * Add script to plot recoil mass out of full sim * Push work in progress for FCCeeGeneralOverview.md * Add script to plot calo energy resolution * Update plot_calo_energy_resolution.py * Add ALLEEGRO and IDEA steering file --- .../FCCeeGeneralOverview.md | 197 +++++++- .../plot_calo_energy_resolution.py | 43 ++ .../FCCeeGeneralOverview/plot_recoil_mass.py | 26 + .../FCCeeGeneralOverview/recoil_mass.png | Bin 0 -> 9655 bytes .../FCCeeGeneralOverview/run_ALLEGRO_RECO.py | 448 ++++++++++++++++++ .../FCCeeGeneralOverview/run_IDEA_DIGI.py | 67 +++ 6 files changed, 765 insertions(+), 16 deletions(-) create mode 100644 full-detector-simulations/FCCeeGeneralOverview/plot_calo_energy_resolution.py create mode 100644 full-detector-simulations/FCCeeGeneralOverview/plot_recoil_mass.py create mode 100644 full-detector-simulations/FCCeeGeneralOverview/recoil_mass.png create mode 100644 full-detector-simulations/FCCeeGeneralOverview/run_ALLEGRO_RECO.py create mode 100644 full-detector-simulations/FCCeeGeneralOverview/run_IDEA_DIGI.py diff --git a/full-detector-simulations/FCCeeGeneralOverview/FCCeeGeneralOverview.md b/full-detector-simulations/FCCeeGeneralOverview/FCCeeGeneralOverview.md index ead4a927..f0a52d16 100644 --- a/full-detector-simulations/FCCeeGeneralOverview/FCCeeGeneralOverview.md +++ b/full-detector-simulations/FCCeeGeneralOverview/FCCeeGeneralOverview.md @@ -3,51 +3,211 @@ Welcome to this general overview of the FCC-ee Full Simulation. -This tutorial aims at showing you how to run the state of the art full simulation of the various detector concepts currently under study for FCC-ee: CLD, ALLEGRO and IDEA. The DD4hep geometry descriptions of these detectors are hosted in the [k4geo](https://github.com/key4hep/k4geo/tree/main/FCCee) GitHub repository and made centrally available with the Key4hep stack under the `$K4GEO` environment variable. +This tutorial aims at showing you how to run the state of the art full simulation of the various detector concepts currently under study for FCC-ee: CLD, ALLEGRO and IDEA. The [DD4hep](https://dd4hep.web.cern.ch/dd4hep/usermanuals/DD4hepManual/DD4hepManual.pdf) geometry descriptions of these detectors are hosted in the [k4geo](https://github.com/key4hep/k4geo/tree/main/FCCee) GitHub repository and made centrally available with the Key4hep stack under the `$K4GEO` environment variable. This tutorial should work on any `Alma9` machine with `cvmfs` access. So, let's start playing with Full Sim! +## Setting-up the environment +```bash +# connect to an Alma9 machine with cvmfs access +ssh -X username@submit-test.mit.edu +# set-up the Key4hep environment +source /cvmfs/sw-nightlies.hsf.org/key4hep/setup.sh +# create a repository with the tutorial +mkdir FCC_Full_Sim_Tutorial +cd FCC_Full_Sim_Tutorial +git clone https://github.com/HEP-FCC/fcc-tutorials +cd fcc-tutorials/full-detector-simulations/FCCeeGeneralOverview/ +``` + ## Towards Full Sim physics analyses with CLD The CLD detector has a complete geometry description and reconstruction chain. It is thus a very good candidate to start full sim physics analyses. To illustrate that, we will process some physics events through its Geant4 simulation and reconstruction, look at automatically generated diagnostic plots and produce ourselves a higher level quantity plot. ### Running CLD simulation -Let's first run the CLD Geant4 simulation, through ddsim, for some e+e- → - - +Let's first run the CLD Geant4 simulation, through ddsim, for some $e^{+}e^{-} \to Z(\mu\mu)H(\text{inclusive})$ events, taken from the generation used for Delphes simulation. + ```bash -source /cvmfs/sw-nightlies.hsf.org/key4hep/setup.sh - +# retrieve Z(mumu)H(X) MC generator events +wget http://fccsw.web.cern.ch/fccsw/tutorials/MIT2024/wzp6_ee_mumuH_ecm240_GEN.stdhep.gz +gunzip wzp6_ee_mumuH_ecm240_GEN.stdhep.gz +# run the Geant4 simulation +ddsim -I wzp6_ee_mumuH_ecm240_GEN.stdhep -N 10 -O wzp6_ee_mumuH_ecm240_CLD_SIM.root --compactFile $K4GEO/FCCee/CLD/compact/CLD_o2_v05/CLD_o2_v05.xml +# NB: we run only on 10 events (-N 10) here for the sake of time ``` +This will produce an output file in edm4hep dataformat with Geant4 SimHits that can then be fed to the reconstruction step. Note that ddsim can also digest other MC output format like `hepevt`, `hepmc`, `pairs` (GuineaPig output), ..., and of course also has particle gun**s** as we will see later. More information can be obtained with `ddsim -h`. + ### Running CLD reconstruction -Inspect the ROOT file content with +Let's now apply the CLD reconstruction (from ILCSoft through the Gaudi wrappers and data format converters) on top of the SIM step: +``` +cd ../../../ +git clone https://github.com/key4hep/CLDConfig.git +cd CLDConfig/CLDConfig +k4run CLDReconstruction.py --inputFiles ../../fcc-tutorials/full-detector-simulations/FCCeeGeneralOverview/wzp6_ee_mumuH_ecm240_CLD_SIM.root --outputBasename ../../fcc-tutorials/full-detector-simulations/FCCeeGeneralOverview/wzp6_ee_mumuH_ecm240_CLD_RECO +# Change the EvtMax variable if you want to run on more events (-1 means all events) +# And do not forget to modify the geoservice.detectors variable if you do not use the central detector +cd ../../fcc-tutorials/full-detector-simulations/FCCeeGeneralOverview/ +``` + +This created an edm4hep ROOT file with a bunch of new RECO level collections, including `edm4hep::ReconstructedParticle` from Particle Flow (PandoraPFA). You can inspect the ROOT file content with + +``` +podio-dump wzp6_ee_mumuH_ecm240_CLD_RECO_edm4hep.root +``` +A detailed documentation on the collection content still has to be written. + +NB: this step also produces a file named `wzp6_ee_mumuH_ecm240_CLD_RECO_aida.root` where you can find a lot of debugging distributions such as the pulls of the track fit. + +### Plotting the Higgs recoil mass + +Let's now use the reconstructed sample to produce some physics quantities. As an example, we will plot the Higgs recoil mass using the Python bindings of edm4hep. The following very simple script already does a decent job: + + +``` +from podio import root_io +import ROOT +ROOT.gROOT.SetBatch(True) + +input_file_path = "wzp6_ee_mumuH_ecm240_CLD_RECO.root" +podio_reader = root_io.Reader(input_file_path) + +th1_recoil = ROOT.TH1F("Recoil Mass", "Recoil Mass", 100, 110, 160) + +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: + 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: + th1_recoil.Fill((p_cm - p_mumu).M()) + +canvas_recoil = ROOT.TCanvas("Recoil Mass", "Recoil Mass") +th1_recoil.Draw() +canvas_recoil.Print("recoil_mass.png") +``` + +Let's run it on a sample with slightly more stat: + ``` -podio-dump wzp6_ee_mumuH_ecm240_edm4hep_recoed_edm4hep.root +wget https://fccsw.web.cern.ch/fccsw/tutorials/MIT2024/wzp6_ee_mumuH_ecm240_CLD_RECO_moreStat.root +python plot_recoil_mass.py +display recoil_mass.png ``` -### Analyzing the output +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. -Let's now use the reconstructed sample to produce some physics quantities. +## Towards detector optimization with Full Sim: ALLEGRO ECAL -TOBEWRITTEN +In this section we will learn how to run with a modified version of a detector (needed for optimization), taking the ALLEGRO ECAL as an example. -If we had to analyze a large amount of events we would of course use FCCAnalyses or plain C++. +### Modifying the sub-detector content -For a realistic analysis you would also need to access much more information than what is shown here. To understand what member can be called on a given edm4hep object, you can use this [Doxygen page](https://edm4hep.web.cern.ch/classedm4hep_1_1_reconstructed_particle-members.html) +So far we have been using the 'central' version of the detector, as directly provided by Key4hep. The first thing we have to do to run a modified version is to clone the repository where the DD4hep detector geometries are hosted and update the environment variable `$K4GEO` so that it points to your local folder: + +``` +# go in a 'clean' place i.e. outside of the git repository in which you were before +cd ../../../ +# load the Key4hep environment if not yet done so +source /cvmfs/sw-nightlies.hsf.org/key4hep/setup.sh +git clone https://github.com/key4hep/k4geo +export K4GEO=$PWD/k4geo/ +``` +NB: a typical DD4hep detector implementation has a C++ detector builder where the shapes, volumes and placedVolumes are defined plus an xml "compact" file that configures the detector (e.g. setting the dimensions) and from which the C++ builder retrieves all the free parameters. When modifying the compact file nothing has to be recompiled, when touching the C++ detector builder, `k4geo` has to be recompiled and some more environment variables have to be set in order to use your local detector builder. This can be done as follows: +``` +# No need to run this for the tutorial, it is just here for completeness +cd k4geo +# modify a detector builder, e.g. detector/calorimeter/ECalBarrel_NobleLiquid_InclinedTrapezoids_o1_v02_geo.cpp +mkdir build install +cd build +cmake .. -DCMAKE_INSTALL_PREFIX=../install +make install -j 8 +cd .. +k4_local_repo # this updates environment variables and PATHs such as LD_LIBRARY_PATH +cd .. +``` +In the xml compact file, the detector builder to use is specified by the `type` keyword in the `detector` beacon (e.g. ``) and it should match the detector type defined in the C++ builder with the instruction `DECLARE_DETELEMENT(MYCOOLDETECTOR)`. -## Optimizing a detector with Full Sim +Now lets first modify the sub-detector content of ALLEGRO. Since we will deal with the calorimeter, let's remove the drift chamber to run faster. For this, you just need to open the main detector "compact file" (`xml` file with detector parameters) with your favorite text editor and remove the import of the drift chamber (line 39), for instance: + +``` +# copy paste won't work here +vim $K4GEO/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/ALLEGRO_o1_v02.xml +:39 +dd +:wq +``` + +### Runnning a first ALLEGRO ECAL simulation + +Let's now run a first particle gun simulation: + +``` +cd fcc-tutorials/full-detector-simulations/FCCeeGeneralOverview/ +ddsim --enableGun --gun.distribution uniform --gun.energy "10*GeV" --gun.particle e- --numberOfEvents 100 --outputFile electron_gun_10GeV_ALLEGRO_SIM.root --random.enableEventSeed --random.seed 42 --compactFile $K4GEO/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/ALLEGRO_o1_v02.xml +``` + +And apply the ALLEGRO reconstruction, including an MVA based calibration: +``` +wget https://fccsw.web.cern.ch/fccsw/tutorials/MIT2024/lgbm_calibration-CaloClusters.onnx # needed for the MVA calibration +k4run run_ALLEGRO_RECO.py --EventDataSvc.input="electron_gun_10GeV_ALLEGRO_SIM.root" --out.filename="electron_gun_10GeV_ALLEGRO_RECO.root" +``` + +Now let's plot the energy resolution for raw clusters and MVA calibrated clusters: + +``` +python plot_calo_energy_resolution.py electron_gun_10GeV_ALLEGRO_RECO.root +# display with your favorite png renderer +``` + +Look at both distributions to see how the MVA calibration modifies the cluster energy. NB: we are in the middle of a transition for the Geant4 interface and we are here using the new version which is not yet fully validated, this explains the drop in energy. + +### Changing Liquid Argon to Liquid Krypton + +In a sampling calorimeter, the ratio between the energy deposited in the dead absorbers and sensitive media (sampling fraction) is an important parameter for energy resolution: the more energy in the sensitive media the better. Liquid Krypton being denser than Liquid Argon, it is an appealing choice for the detector design. Though it is more expensive and difficult to procure in real life, testing this option in Full Sim is extremely cheap: + +``` +vim $K4GEO/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/ECalBarrel_thetamodulemerged.xml +# change "LAr" to "LKr" in line 114 and 121 +``` + +And to accomodate this change, we have to update the sampling fraction in the steering file: + +``` +vim run_ALLEGRO_RECO.py +# comment out line 105 uncomment line 106 +``` +The derivation of new sampling fractions upon geometry change is done by ECAL experts and is out of scope for this tutorial. + +Let's now re-run the simulation with the modified detector: +``` +ddsim --enableGun --gun.distribution uniform --gun.energy "10*GeV" --gun.particle e- --numberOfEvents 100 --outputFile electron_gun_10GeV_ALLEGRO_LKr_SIM.root --random.enableEventSeed --random.seed 42 --compactFile $K4GEO/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/ALLEGRO_o1_v02.xml +k4run run_ALLEGRO_RECO.py --EventDataSvc.input="electron_gun_10GeV_ALLEGRO_LKr_SIM.root" --out.filename="electron_gun_10GeV_ALLEGRO_LKr_RECO.root" +python plot_calo_energy_resolution.py electron_gun_10GeV_ALLEGRO_LKr_RECO.root +``` + +See how the resolution changed. + + +## Towards IDEA tracking with the detailed Drift Chamber + +``` +ddsim --enableGun --gun.distribution uniform --gun.energy "10*GeV" --gun.particle e- --numberOfEvents 100 --outputFile electron_gun_10GeV_IDEA_SIM.root --random.enableEventSeed --random.seed 42 --compactFile $K4GEO/FCCee/IDEA/compact/IDEA_o1_v02/IDEA_o1_v02.xml +``` -## Running IDEA simulation with detailed Drift Chamber ## Detector visualization @@ -75,10 +235,15 @@ Load the file from your computer with the triple-dot button of this webpage: htt From there, you can navigate the geometry hierarchy to see how volumes are nested, display all or parts of the detector (right click on a volume > Draw > all), play with the camera settings (right click on an empty space in the display window > Show Controls, Clipping > Enable X and Enable Y), etc. See the following picture for illustration: -```{image} https://fccsw.web.cern.ch/fccsw/tutorials/MIT2024/pictures/JSROOT_Screenshot.png +```{figure} https://fccsw.web.cern.ch/fccsw/tutorials/MIT2024/pictures/JSROOT_Screenshot.png :align: center ``` This tool is useful but not perfect and will not meet all the needs (especially if you want to overlay a physics event). To go further, other solutions are described in the dedicated [Visualization](https://hep-fcc.github.io/fcc-tutorials/master/full-detector-simulations/Visualization/Visualization.html) tutorial. +## Additional Resources + - [Key4hep tutorial](https://key4hep.github.io/key4hep-doc/setup-and-getting-started/README.html) + - [FCC Full Sim webpage](https://fcc-ee-detector-full-sim.docs.cern.ch/) + - Bi-weekly FCC Full Sim [working meeting](https://indico.cern.ch/category/16938/) + diff --git a/full-detector-simulations/FCCeeGeneralOverview/plot_calo_energy_resolution.py b/full-detector-simulations/FCCeeGeneralOverview/plot_calo_energy_resolution.py new file mode 100644 index 00000000..d4cf12f5 --- /dev/null +++ b/full-detector-simulations/FCCeeGeneralOverview/plot_calo_energy_resolution.py @@ -0,0 +1,43 @@ +import sys +import ROOT + +# prevent ROOT to display anything +ROOT.gROOT.SetBatch(ROOT.kTRUE) + +f = ROOT.TFile(sys.argv[1]) +events = f.Get("events") + +# raw clusters +c = ROOT.TCanvas("c_clusterEnergyResolution", "") +h = ROOT.TH1F("h_clusterEnergyResolution", ";ECal Barrel Cluster Energy [GeV]; Number of Clusters", 100, 5, 15) +events.Draw("CaloClusters.energy >> h_clusterEnergyResolution") +fit_range_min = h.GetXaxis().GetBinCenter(h.GetMaximumBin()) - 3 * h.GetRMS() +fit_range_max = h.GetXaxis().GetBinCenter(h.GetMaximumBin()) + 3 * h.GetRMS() +fit_result = h.Fit("gaus", "SQ", "", fit_range_min, fit_range_max) +mean = str(round(fit_result.Get().Parameter(1), 2)) +resolution = str(round(fit_result.Get().Parameter(2), 2)) +legend = ROOT.TLegend(0.47, 0.65, 0.8, 0.8) +legend.SetBorderSize(0) +legend.SetFillStyle(0) +legend.AddEntry(ROOT.nullptr, "#color[2]{#mu: %s GeV}"%mean, "") +legend.AddEntry(ROOT.nullptr, "#color[2]{#sigma = %s GeV}"%resolution, "") +legend.Draw() +c.Print(sys.argv[1].replace(".root", "_clusterEnergyResolution.png")) + +# MVA calibrated clusters +c = ROOT.TCanvas("c_calibratedClusterEnergyResolution", "") +h = ROOT.TH1F("h_calibratedClusterEnergyResolution", ";ECal Barrel Calibrated Cluster Energy [GeV]; Number of Clusters", 100, 5, 15) +events.Draw("CalibratedCaloClusters.energy >> h_calibratedClusterEnergyResolution") +fit_range_min = h.GetXaxis().GetBinCenter(h.GetMaximumBin()) - 3 * h.GetRMS() +fit_range_max = h.GetXaxis().GetBinCenter(h.GetMaximumBin()) + 3 * h.GetRMS() +fit_result = h.Fit("gaus", "SQ", "", fit_range_min, fit_range_max) +mean = str(round(fit_result.Get().Parameter(1), 2)) +resolution = str(round(fit_result.Get().Parameter(2), 2)) +legend = ROOT.TLegend(0.47, 0.65, 0.8, 0.8) +legend.SetBorderSize(0) +legend.SetFillStyle(0) +legend.AddEntry(ROOT.nullptr, "#color[2]{#mu: %s GeV}"%mean, "") +legend.AddEntry(ROOT.nullptr, "#color[2]{#sigma = %s GeV}"%resolution, "") +legend.Draw() +c.Print(sys.argv[1].replace(".root", "_calibratedClusterEnergyResolution.png")) + diff --git a/full-detector-simulations/FCCeeGeneralOverview/plot_recoil_mass.py b/full-detector-simulations/FCCeeGeneralOverview/plot_recoil_mass.py new file mode 100644 index 00000000..bdb93620 --- /dev/null +++ b/full-detector-simulations/FCCeeGeneralOverview/plot_recoil_mass.py @@ -0,0 +1,26 @@ +from podio import root_io +import ROOT +ROOT.gROOT.SetBatch(True) + +input_file_path = "./wzp6_ee_mumuH_ecm240_CLD_RECO_moreStat.root" +podio_reader = root_io.Reader(input_file_path) + +th1_recoil = ROOT.TH1F("Recoil Mass", "Recoil Mass", 100, 110, 160) + +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: + 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: + th1_recoil.Fill((p_cm - p_mumu).M()) + +canvas_recoil = ROOT.TCanvas("Recoil Mass", "Recoil Mass") +th1_recoil.Draw() +canvas_recoil.Print("recoil_mass.png") + + diff --git a/full-detector-simulations/FCCeeGeneralOverview/recoil_mass.png b/full-detector-simulations/FCCeeGeneralOverview/recoil_mass.png new file mode 100644 index 0000000000000000000000000000000000000000..29057686b643f45bbc874d95e8a5764c300e929c GIT binary patch literal 9655 zcmch7cU%+O+U|g$6zLIAy6i0|p$SL_m8wV+K~OpZqJW|IqDT>tB1#Q43KAqBozNl? z5otoe&^rl`4FTyjcW|F`uD^4>@4NS(d-F#!nVGD4*Lt4!dCSTRBLf{KhI0%605ILq z)w~S=)O7$rHAF`Zp1Ih3x&{0=X|JcF34VgFU5$n?01yOjXkNSH4__WvPZ7JDEf3o? zX#{fo5pB)BZ{)||?zH8;9{SXwj2D1YHTaDD1}KaXrE}f!+K;N$a;JWn)*kV0b}a<^9C3~=`_cmCos{7K2yoR10Ji!H z0N?>V18_Ap@Arcc@h}F3iZts$7h*E_zKBRvQ<|SYYLP$M!?&Q!O2{GPps<5k0m@cw z%W@+zF?0iAsD9{OeJ^Bl@a{f&b*ZE*^kA($(GoT_KtrpL>^BiW?38Tb<%<_?*&ipX z?arl0Ht)(+zWs&l0tfTD~;>zyeU7N#s$*>o^tWRy{aE9tHm4UdoBr1cYjueqm zV$T1a4p{lU6AWnI2`?WMBnTnlBt#g9LY-F6s!TxFu#o@hnu#ru$BbB`7c z!;}t=BpX++Z(NQ~XirvPgucoHh+HQig?t=kOFU=Y;HS zEb7S@x}BNa`67R`X+zmuWd|*eQ;xo<>8uIo3)~*}vk5!gOrzi%$j!=VY!JCAttWK5 zVX_HlG4)ovG!G;E!YRT9kL@wRkyq9>RmZ1A)OR~_)QSGB$+eTA%Hww6avhRv5-Aj% z`XMeS@mbh@ZP@)%GxdXM7R0cv_mxbJ{*aj|WpblT$YzmYZlTWMJ~+>uRRx;pOPJzPgpHagH(riL4Zhgid%4+;9 zt(}Nv*hyQBZf_%%a{)WleEOd3~R^oG$pr3Oct4^o&cji(IZyUwr z1Be6)Wv`5KDfYZU$k8sbWDw~#I?Xo|7jzDy_>Cqv_~#8{2n(;|ZY!hYX`;dbXwR@> z{913=(fTR@1fpTx=a^v42R}R2_lJAWGmR}mVVCbgrrzeL1pI7UuNhrV5ivZzJWws* zGBum1o^wUzbVRG_t#-YMhHW7{nYDmPFjmz{?~=CA0m$ydRV+(ff&cD&+9ZhZiC$`v zLM^Hu@*=LNWv5M2+or(=r2L}mv-Q0VoMBu6_tv$){_=~tl~QPR08SYs2tL{Y66!nK zv^}xE{_>dgSQhnz?kAxf`fU5)Xct5=mT&V`XW}G7LDvDV+n+|9y}vV?=pZ#KR@QhKKFj(Grob5T3PO=t`hguKuT?|SM&mWixkY1O0l{_dib&)3 zL{MM8RHFu(a#V>q3xB|Vq5KJ7q%5=aczotg8<$P|*;mCk`Nt{yla$a8#jJhDguVfP znHnKtn9~>{T&NPyABy*ZXiKUFe_d6R?w@&=W*xhjOm18^KoW4zJL+~pnzbY8rKv8f zPJ@qlu{#=ZmNs%aQx`|Ype*o8^Pz5G(`iKegtDk+B9>9IF{z>VJ(b{ftdy8Jj9x-l z`mm0QFbmRgMA2Gkz~^w6nA-xb5~BE32k0d0poQ>I9#T5{S>dzMI#DOM?CWM8Ro}!3 zMsMZP8fb%ZS8>{j2B@H+1qAhuzy0d~I(-_G{O(XKz48QJhZO=S-2e0Q0Tzxr^}lZ< zpj9KVDZ2I9enyu!=Kz1QSaVS9r_R#LZWB0^mWp zkAYOe!W+FiO^65|D;Ab7-9NPtoXvCY96NQllT7f)MiLrT)Z2J&ENp|HgidHLjL7P~ zu#(HqPZ4O|TQH3(+7eJlVROxlfGO8E-=7wI8#}G(mU(IXj)virgrAz#SWy|HCZPnyiMJG*l#FIk3?STyqvi= zTBASNhxn1GEU#~p>KLetlXV;Qwel;XdNtUcjYI4;8|u?xpZ$rmezBpVtaExUA@`Lu z17I1o=rvDIHQz2upe%lVUI8Ub&sQ|w_ybq|N;5t^^-<;UMuUg+q2NZJt1_Dv4H2D+h~F4?SY|y z151zv68-P9^T_*Sp0}`{zsu@ZX}oq+w91<1XOy`k&+?*9>g9!yNxToA3Ac99Oiy1c z8`B@fh7X8qUThGc&fmDcFi`0v?`H}Ly|S-1SnW3Y1;5v)aXTa|@vda&Wovim9$iV* zdvV2Xquvxd$X7n~EFJEf^cYj+?v1b695ka!Jk)9JK*L^X1uTPho~_mMdMVbE3CPu4bN>4zBllv$UGq*Z-__n#c9T zsz0Y$J+ycFn%z^9=*}{Yg5>0I>Vv1(sKP5#&T!5?#;=bSzoVLZ5IiSwp+j8i_?B>o zj|tQOOkbP<=rY0kDBqR6xi2}^`*88?pNm^PTVv`0KT$lE^AC-IDNTEq)_qmqs!Qt& zE;yz5xV|x{TvN03$cj)+dhBd{{>hK+pN(%8fE-RiPgZ;mV*djDp{RzeQuZC$2AR=L z*En>?zH(Ue@ik)zj_%D?u@b}aQ?*9Q>g%PICSk@=!wH+9z%&a24j%kX~WcF;;ZYOWFmV7K=l zS48(Nj4HATHb1JBFGvjQ!ClPC>WeAvcuPWW2J&?a8UN)k+^7Aj;`$hIA@Z1HDnUOr zQ5Cm@X&mZoQ0w;0BUu+&7M;#jo~2}LLEq-Q?tA@pq0!%#2TL9Ctcf2$jln)^43w?q z-f8>dx8WduV}skXdjFwMv22*<>!g0@x9BvG6LvN9Y~hw?{Bo zsWDnJ`0`MUzvr42-ffgQi#6(<{*wo?9G*MANaWJ~k?0fHF`pDi>+Pq6Cf5p9#wL^a z*M)E6JRkWRp0%gIzGggh-_y~4r=M@##vIbt|FH)L12pf6RT-d(qZT zpe^s@tDoXb+bV@{-%la1vgm@m18+#DSyygSzwW4=S-`w9k|D#>^hCN!z6lxrRd_@hT z!SybeEZtnWOw!YK2)KC?>vP5-q(k8ir`y;%d%Ihfq8FIxoKz0eI!|*;M+A(ST|dqM z{1(C6JR$DeA?gPP-2FdJ<6;HTQrbngJA;)PGII%q+rj7F6wDRR3pR_ikN0QY>@9DZ zQ9io6u&i>cfxkbw0%{`D+0z^~eL`~YJWL(UKlKdFKjRQQg5swXMkMpx4%Dr_io97g z{-c;yDSKpcE7ZD#7{pp29O@tK&!tR;lRgOKWSgbo8mw;Nu?AYkZ+%M1OBg>jv^4<# zvLN;ZCyO9FufzH_ETbjlYKlzO6a|X3$2mXR*u!s5ewDnwdBEL^B1M%#J!tJ4o&Vy^ zE)fv89YpAo7X927YYiWwT#_&vVSeo0TXIg65w@t^^d)365ltVZd_LidnZd4YksVWD zy{U#5PE~q*`J!o9eeIQzr-2>UgrelNZ^m>chi=s`s{(a`juy=OBSmsHliMHZ=U`%A zSv6P-+qv=@GPWjb+^q1$&*B`Z3(AB%Z;#gCWRq5W5A5BDc7}a=Ct= zJGc=*nw+)}lN@(~3QvM@OoLkvS5d*cd=yDGMtqPOVZEDfSld*!m^}Z$_7kl9`igE~ zi%EU`>u4KzS*b*~Q>8^2zPRNkf*RBI;=*(8xTJk{6<4Piv*!<8-jf1zWX|zFY_cW? zKSro+!@XbhS8Bc7U6n|$r-@VzF~D6+?mO`8voZqiqM)5S-{cSfNR)i;()9HGMZYUG zWcb~H7nww7^{4?_fe@zR$_xGdt7-_XsQ!V(*@0K1jusQjGiD&XGid%o@nytfG9qeY zU{gwkpAT2;b6C*8|LxY;aDZIY6Pdej1Uv#Cd@S#Hs`atT?{jJF%&_l1;Yb+sXDEd1 zs$$UHg&pR=f5(jRZD#Wy$=thFeM!_(=o%H_%P*Jng7eb@UC+U?uq*TjmTf+D%G_Lp z5rAVriMa{@?FA$%pstXU79f0Z$*eNFeGj_8e;>-6nCMjA)H~G`(s~_92J*N`P9%QOHL#1WMG6R&=Hy`Lk zSt>^BN$KA;|EAip;J;!&QH-|yr-2gl2!C{z@f}LB+Ng=g3JU{5M=PYUC{rcvXALa} zZd^O<&+Yfqju~o?!t9mHyS{*$dmdMM7Kyy6({Fte`(txxb_%smX)|g2qKl0{L}P2*8c_Gk z@P6aco&#S648a;`q((gUytPd(DZ-oVA6yi*yl-hgpReJ1<*haIbs( zbzTys6Ba;@xxW|~9>;@7s`1Xu;q$E(sMt3~VeJPlsGp-}lh@%B5E;^WNBp*hZ57KcD7VpwbGeSsC<6^=De&vhxqPV|Kc6IkH}CJ`e>4d zRoam_T1HZ;NlmDUe#m&XT2!-4e`Tp{B%E2PT+#8kLs=w+B6W&vR>-pH z9y_PkFvU9(YQg!q#IO(RZroT~*+;uWHY1Ccud>19+19s^QAksA-}~(Zc7H@i@(zBJ z5P&5--Sw?0Ro;6#-RVwRconuc>=x?>v1`Hl3z4B<8&L;0+(y5ak;i?#Lu&D>^ z-TT|5aRtjIwmQ?UO>Q+9gwZ#(Ml`kpnyAk0)d0n>8{F_##WlTw-#dc#kO=6z8kuJo3LXSU1{Ud#=&@^&d@>qPyj~SZWr@k)E8Gjk(SBt=Dtk zFF#4s60E!ajN0z`n&y%g&_x}3G22vqKgPE-nfdY=(8vK0EHeUfInnK{@*Cz#Aq zERkkb%D`|}*C@Px&A9PxZt|!2fAmz2TPeWGA|@cbD;K0wvbgKU=$pje@0-8(Tfhzs zMvlp;i>cabXn`G2cuISc_9*})XyH!)Nq>q0fWSv^E6{b$Dg^av?ftX|VC!;CXyjzUi~yLLry4Aw-5bVJ z9Wz#NHOo$wi|(%-i>;`D)Ug61Q&THEA|hmxg+<$4R~3TI3Yp9;F9%>&J&n)3e1u2n zx1d<~9+x&B)7sF(RD~M>R8we6j5Uj?>1IgFuY_X%UT`zWO}-1Uurr z6`l!i#r>v5l<+gbUqx%p6@0*bGO(%}+XhA$7h`U#AF6ljkx28Hjs+T^E;SE!_oNhP zNN*qj^cehJ8FI>}uki<|LS35~{6<(>U2lZ%Ti6=0(s8#@XX#RnW&K_u{##A?uj2Ce zgQ?7|!19X_fVUG=QT=>i%G_J^zXJC^;XXbC)7MiN1!XtE@#(cs)+_)F3MQZ7&}&q{ zS9Zw19i)dCJenKMAC>}mC8ihn=vw+}oLeQwkqrMQj*n)?nzMypBNl0Y>TA|Ai zl8*{fG}a{^065tp1#W_uj`AOh54qQZSzzsGGucECyuWhTTmDPHOktkB% zr(-;B$(snqgqfVsz{!G4sc8gYT0OiF`}d!E%%0l$X_%U_(Scm?5NKinQIL327Hy}8 z&At*jg56>T%9HY7V%2H@)H%L2`>k3GYc@01fB~@Gy_OH{oceIeDK)p0lf`jH{vmJF2WSIz6JXJa;idbn7m8_Pc4~|s zvu0uBiDCdm$}xF@tN`XCD8|3_4t|ICzo7uirR-GJbt=q9sFdgsgyDa5>;FU^U-i5a z%;~4hIk_Q!fvWK<^8P2_zs__XOMsldiK!{x2Vqh_3t&uXK(7O_)qJ`cXk}`uWD3y& zGb8XqH4W%VV9DwZ*2&FeVOcPWlw9($svUR7c>q(f&sxFONhvMXDdKN@J1gQB0Kn4N zNVM}Znjejd%PW_WzhgSMeR%AFdRHc!z>{LJ_F$FQ-R-oF(a+y1SET_g&Eu{cLmVyK zpA7h9pS&{hD9A67;tp2n-Lj&Z7j!!WwS~B==m3!`xI8J-R7aY~joj}(<1klYp!^_@ z6?+Ef9fcx+(KJgR_x55k7PR{elhc^%71E{>r`f$yGasv>dKR#q<;@Eyph38I3+8Q| zA-cE=(XsX`Ws)D3-#-%c7O z+1X)@sHSMZ^`@Z)6gaI|{=Gf^pAGYWp0UKS_kg{ga0gH(%0cnwjr)Hk5nMcw%)(+M z(j_yyy#b|cQV)7Uu-lQ?c`#q8ar1NV4sNnebG#@4Uj_~{y2Ac|J5#eDtM1; zX2-Wdznv-^A?DPeVr&e_>Mleg7uW2>u8-Lb4M~mqef|h(t6bjdZf{V>G@XeniE}!@4ztI7m@*sAuD(0nTA(ORkLAKhR%bVrxt^Jd&;P1V-sHOtD zl|(`G&#{A&<36@0T7C3$K><)m+hum*^^s2>2Y{mEAG5M7(KJ}$kl!n4p63NutXR`%Sp0V6Tp$DGVeiL8EAp_3QV4@iRgbvPfX;g zNF-8Q=NCmtXqh*30IHRi20XXdEa3F4Ky;+Ws??0M8JH;qgg>^>$^Rv);O%E6QjH*5 zc%K<+z|tYEHLwXq3iN;LHd+D2byq-zl>h*L27*K)Nr9OsqLxqfNqhV(pu(7@Kf#0A zkI1KI86(-Tl&GD+{7D;$B8}8E*&LSsI28_>v*6C{D-DWs1++VwCfp?Ri4nKcvt}|b zug%zU0qw3(!v-i6%$z^B6juOd;e8lT2LhWCvapn*AnVV{GcYiYH-xPB+fdiBx6L5) zGipS~l+OIp8#GEB+{Q|rdlu2#H( z%5?I+l;CIOSNJz5$2mZgwGIRf-E~T!5dTUFj|0J=Mcgmd&ij5`HIKhEKPNPDW@d8P zm#2Nd84t#HeBD%kd37Gx2SXdB&ae98WGkpy9%$)?GcqzDXgT+r%SeIe1d+&CcRGb1 z3H^ZORlG*OkJV|TbY_>a;<|?RFb)8gOyCE(a<2^UXu;ncf!$i+-_-{&?=pRiV#s@0 z-CtN~0EM{S1Y^nP-!1A9Fsbh=}#n*kV3Fj&(~l-9sJp z`_&-{1u&IP)CV&>2{$R>t}t~v4bJe#f})sf+8;3&Ze!jYXXTbAao62biXyx3QiXGv zA7^`i$KEMuA%c@%ZDdV4hhxeB+-A(Q6A%5+myoyy2cl)t*?_P81px75Wx&;5YAJbe z5A9|3+3;Ns$W%@|H-LJY2W;mj0(I;sWvBs*dk|Yit1~BcsdYp|F$GG{(@8W5kAEd~ zu(TGgO4+AkjG$VCD<5?XUE}$B865n$)J~^o6#uK9nUL1{jukRx{GM*xl8KE$A|?-z z@qEnmXx3%wGYxHBx^~PWU$ERM{;Xd$2q$e%m)}f=%~cJy>2LH5@LZ2DnHKM;^R03m z5Zh^OeQKQ!&iz~yF*fLnJV0-Y+jK{$xA+YmCl;*U{vkFZ{0IJt#Kujb&W9=fJ7})N zK3fD^%wrY?u8$Y*7(>=2iCU1`isZR+2EwV5LB6H@pgZ?!)#kz@izH^#ysGy;mBY97 zTtpVtFRg57sK<5g=1X={M^IztLtU;wUe~ih!X43L_1s)+59|2sS(67pYMXt=JSwA4e^C9FqOzq2&r8CI3~%1pix^x|@4auJz9DSgWpcQ3I7MN> zw%jky8FyLv;&AXz;ME@q);(GDj4XQ`_q#;;-yOWOfp;=8x7MaUG%9QbfZoRR`~FUtQhpcbB|-z*lNP|i!S})%k9Z9 z3ivApzZZb*M~N~ma0j@Fac(fX=ZW(A)SKY?0(2&h{s;Eq> zXVYRWxVOgxxq3e+uBShHgPXij;U@Ww97fqU&T!D3UA+Ln`+O}VK4{q$aSw4npatrd zK+0e`EX|E#KGI3wwsVdHi;kviLiV4Ir#)Sv$H?752dsMsT+1i#oAZ>2#L6>{*u|IDab}!qnK#6B13rPXb;xY_g9>oy0 zI_oRpUinx$I zOU>qA4nC?Lb~Ui<$u3h6cRqj?z`ts7 zBX0X}b!v5~d&!>gt$v*_9JYXq@M_4}mmIkjs)#GL>ugt%B=mR|&Wc;AD2W|qR$?Es2z2CeGuOQZ?9^)dg{tzVzYOS+Eae>o`D)p^<)lW_aXid zE9|g!S5mCsL|i?%6ohYzHobO=#I4S=0^dETZ~}d&j~~*roaP0wck)zafgWSAROkw( zh0W5JC^@IBpEEu-t}AhY=jPMpn3}q|Z-j)Urh~IFo!nkKrAUlZgc4OGo$x+d#$b>~ zqf7#v{8fwQ;iH?bqgM{EGyhG2F$pK#!TX~%E>+`WC}LHTOk+A2lklYDQ~_+ERY!I$ z+fTIn_K7nf(gU{0XBT&JJmiDEXsbqLm%TanW6Qg=d)nThGhiYv;H&)RPG9eOAR%h1 zdoImuqttJQ0>WvjNId9Ne59r_yQxy1mJdJcq!)lphGfa&X_H!K% zV##pz%O!Cj4Wm_7jhDmfz0W%Pt43?LDIYv-mkTcKTiD*QOh}vp2UVRPqhL7Fye?VQ zu3Pt>U<|;Y%DK$vPTWWfh3KsWa#w3E3qZWIyI&5fd^XQo_Eo;r56%Q01nubCDWLL%c|U)=RGI4*n7+%MXXIvUwrel zoC^LQ^Uosgs}a^gLm9zN>6)epbawok;!T7z*=022g0qkG_UrX{E}8AMjw!!j%UvjxUzn-DK