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 00000000..29057686 Binary files /dev/null and b/full-detector-simulations/FCCeeGeneralOverview/recoil_mass.png differ diff --git a/full-detector-simulations/FCCeeGeneralOverview/run_ALLEGRO_RECO.py b/full-detector-simulations/FCCeeGeneralOverview/run_ALLEGRO_RECO.py new file mode 100644 index 00000000..436762ad --- /dev/null +++ b/full-detector-simulations/FCCeeGeneralOverview/run_ALLEGRO_RECO.py @@ -0,0 +1,448 @@ +from Configurables import ApplicationMgr +from Configurables import EventCounter +from Configurables import AuditorSvc, ChronoAuditor +from Configurables import k4DataSvc, PodioInput +from Configurables import PodioOutput +from Configurables import CaloTowerToolFCCee +from Configurables import CreateCaloClustersSlidingWindowFCCee +from Configurables import CorrectCaloClusters +from Configurables import CalibrateCaloClusters +from Configurables import CreateEmptyCaloCellsCollection +from Configurables import CreateCaloCellPositionsFCCee +from Configurables import CellPositionsECalBarrelModuleThetaSegTool +from Configurables import RedoSegmentation +from Configurables import CreateCaloCells +from Configurables import CalibrateCaloHitsTool +from Configurables import CalibrateInLayersTool +from Configurables import SimG4Alg +from Configurables import SimG4PrimariesFromEdmTool +from Configurables import SimG4ConstantMagneticFieldTool +from Configurables import SimG4Svc +from Configurables import SimG4FullSimActions +from Configurables import SimG4SaveParticleHistory +from Configurables import GeoSvc +from Configurables import HepMCToEDMConverter +from Configurables import RewriteBitfield +from Gaudi.Configuration import INFO, VERBOSE, DEBUG +# from Gaudi.Configuration import * + +import os + +from GaudiKernel.SystemOfUnits import GeV, tesla, mm +from GaudiKernel.PhysicalConstants import pi, halfpi, twopi +from math import cos, sin, tan + +# Loading the output of the SIM step +evtsvc = k4DataSvc('EventDataSvc') +evtsvc.input = "electron_gun_10GeV_ALLEGRO_SIM.root" + +input_reader = PodioInput('InputReader') + +podioevent = k4DataSvc("EventDataSvc") + +Nevts = -1 # -1 means all events +pdgCode = 11 # from what you simulated with ddsim +addNoise = False +dumpGDML = False +runHCal = False +# for big productions, save significant space removing hits and cells +# however, hits and cluster cells might be wanted for small productions for detailed event displays +# also, cluster cells are needed for the MVA training +saveHits = False +saveCells = False +saveClusterCells = True + +# cluster energy corrections +# simple parametrisations of up/downstream losses +applyUpDownstreamCorrections = False +# BDT regression from total cluster energy and fraction of energy in each layer (after correction for sampling fraction) +applyMVAClusterEnergyCalibration = True + +# Input for simulations (momentum is expected in GeV!) +# Parameters for the particle gun simulations, dummy if use_pythia is set +# to True +# theta from 80 to 100 degrees corresponds to -0.17 < eta < 0.17 +# reminder: cell granularity in theta = 0.5625 degrees +# (in strips: 0.5625/4=0.14) + +# Detector geometry +geoservice = GeoSvc("GeoSvc") +# if K4GEO is empty, this should use relative path to working directory +path_to_detector = os.environ.get("K4GEO", "") +print(path_to_detector) +detectors_to_use = [ + 'FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/ALLEGRO_o1_v02.xml' +] +# prefix all xmls with path_to_detector +geoservice.detectors = [ + os.path.join(path_to_detector, _det) for _det in detectors_to_use +] +geoservice.OutputLevel = INFO + +# from Configurables import GeoToGdmlDumpSvc +if dumpGDML: + from Configurables import GeoToGdmlDumpSvc + gdmldumpservice = GeoToGdmlDumpSvc("GeoToGdmlDumpSvc") + +# Detector readouts +# ECAL +ecalBarrelReadoutName = "ECalBarrelModuleThetaMerged" +ecalBarrelReadoutName2 = "ECalBarrelModuleThetaMerged2" +ecalEndcapReadoutName = "ECalEndcapPhiEta" +# HCAL +if runHCal: + hcalBarrelReadoutName = "HCalBarrelReadout" + hcalBarrelReadoutName2 = "BarHCal_Readout_phitheta" + hcalEndcapReadoutName = "HCalEndcapReadout" +else: + hcalBarrelReadoutName = "" + hcalBarrelReadoutName2 = "" + hcalEndcapReadoutName = "" + +# Digitization (Merging hits into cells, EM scale calibration) +# EM scale calibration (sampling fraction) +calibEcalBarrel = CalibrateInLayersTool("CalibrateECalBarrel", + samplingFraction=[0.3839086333240162] * 1 + [0.1356521216954895] * 1 + [0.143942924309219] * 1 + [0.1487773592062343] * 1 + [0.15275058492860488] * 1 + [0.15646237236776503] * 1 + [0.15932192403562903] * 1 + [0.16303503690217488] * 1 + [0.16553100078190436] * 1 + [0.16715213916224647] * 1 + [0.1714089975995641] * 1 + [0.16746055755998251] * 1, # ddsim LAr + #samplingFraction=[0.44473486865369044] * 1 + [0.2040814884295127] * 1 + [0.21240641052591336] * 1 + [0.2183384713006271] * 1 + [0.22316850922812564] * 1 + [0.22791094940058731] * 1 + [0.23210241481573746] * 1 + [0.23665653192329097] * 1 + [0.24008008036913617] * 1 + [0.243785972750989] * 1 + [0.24535176801370764] * 1 + [0.23949137241404994] * 1, # ddsim LKr + readoutName=ecalBarrelReadoutName, + layerFieldName="layer") + +calibEcalEndcap = CalibrateCaloHitsTool( + "CalibrateECalEndcap", invSamplingFraction="4.27") +if runHCal: + calibHcells = CalibrateCaloHitsTool( + "CalibrateHCal", invSamplingFraction="31.4") + calibHcalEndcap = CalibrateCaloHitsTool( + "CalibrateHCalEndcap", invSamplingFraction="31.7") + +# Create cells in ECal barrel +# 1. step - merge hits into cells with theta and module segmentation +# (module is a 'physical' cell i.e. lead + LAr + PCB + LAr +lead) +# 2. step - rewrite the cellId using the merged theta-module segmentation +# (merging several modules and severla theta readout cells). +# Add noise at this step if you derived the noise already assuming merged cells + +# Step 1: merge hits into cells according to initial segmentation +ecalBarrelCellsName = "ECalBarrelCells" +createEcalBarrelCells = CreateCaloCells("CreateECalBarrelCells", + doCellCalibration=True, + calibTool=calibEcalBarrel, + addCellNoise=False, + filterCellNoise=False, + addPosition=True, + OutputLevel=INFO, + hits=ecalBarrelReadoutName, + cells=ecalBarrelCellsName) + +# Step 2a: compute new cellID of cells based on new readout +# (merged module-theta segmentation with variable merging vs layer) +resegmentEcalBarrel = RedoSegmentation("ReSegmentationEcal", + # old bitfield (readout) + oldReadoutName=ecalBarrelReadoutName, + # specify which fields are going to be altered (deleted/rewritten) + oldSegmentationIds=["module", "theta"], + # new bitfield (readout), with new segmentation (merged modules and theta cells) + newReadoutName=ecalBarrelReadoutName2, + OutputLevel=INFO, + debugPrint=200, + inhits=ecalBarrelCellsName, + outhits="ECalBarrelCellsMerged") + +# Step 2b: merge new cells with same cellID together +# do not apply cell calibration again since cells were already +# calibrated in Step 1 +ecalBarrelCellsName2 = "ECalBarrelCells2" +createEcalBarrelCells2 = CreateCaloCells("CreateECalBarrelCells2", + doCellCalibration=False, + addCellNoise=False, + filterCellNoise=False, + OutputLevel=INFO, + hits="ECalBarrelCellsMerged", + cells=ecalBarrelCellsName2) + +# Add to Ecal barrel cells the position information +# (good for physics, all coordinates set properly) + +cellPositionEcalBarrelTool = CellPositionsECalBarrelModuleThetaSegTool( + "CellPositionsECalBarrel", + readoutName=ecalBarrelReadoutName, + OutputLevel=INFO +) +ecalBarrelPositionedCellsName = "ECalBarrelPositionedCells" +createEcalBarrelPositionedCells = CreateCaloCellPositionsFCCee( + "CreateECalBarrelPositionedCells", + OutputLevel=INFO +) +createEcalBarrelPositionedCells.positionsTool = cellPositionEcalBarrelTool +createEcalBarrelPositionedCells.hits.Path = ecalBarrelCellsName +createEcalBarrelPositionedCells.positionedHits.Path = ecalBarrelPositionedCellsName + +cellPositionEcalBarrelTool2 = CellPositionsECalBarrelModuleThetaSegTool( + "CellPositionsECalBarrel2", + readoutName=ecalBarrelReadoutName2, + OutputLevel=INFO +) +createEcalBarrelPositionedCells2 = CreateCaloCellPositionsFCCee( + "CreateECalBarrelPositionedCells2", + OutputLevel=INFO +) +createEcalBarrelPositionedCells2.positionsTool = cellPositionEcalBarrelTool2 +createEcalBarrelPositionedCells2.hits.Path = ecalBarrelCellsName2 +createEcalBarrelPositionedCells2.positionedHits.Path = "ECalBarrelPositionedCells2" + + +# Create cells in ECal endcap +createEcalEndcapCells = CreateCaloCells("CreateEcalEndcapCaloCells", + doCellCalibration=True, + calibTool=calibEcalEndcap, + addCellNoise=False, + filterCellNoise=False, + OutputLevel=INFO) +createEcalEndcapCells.hits.Path = ecalEndcapReadoutName +createEcalEndcapCells.cells.Path = "ECalEndcapCells" + +if runHCal: + # Create cells in HCal + # 1 - merge hits into cells with the default readout + hcalBarrelCellsName = "HCalBarrelCells" + createHcalBarrelCells = CreateCaloCells("CreateHCalBarrelCells", + doCellCalibration=True, + calibTool=calibHcells, + addCellNoise=False, + filterCellNoise=False, + addPosition=True, + hits=hcalBarrelHitsName, + cells=hcalBarrelCellsName, + OutputLevel=INFO) + + # 2 - attach positions to the cells + from Configurables import CellPositionsHCalBarrelPhiThetaSegTool + cellPositionHcalBarrelTool = CellPositionsHCalBarrelPhiThetaSegTool( + "CellPositionsHCalBarrel", + readoutName=hcalBarrelReadoutName, + OutputLevel=INFO + ) + hcalBarrelPositionedCellsName = "HCalBarrelPositionedCells" + createHcalBarrelPositionedCells = CreateCaloCellPositionsFCCee( + "CreateHcalBarrelPositionedCells", + OutputLevel=INFO + ) + createHcalBarrelPositionedCells.positionsTool = cellPositionHcalBarrelTool + createHcalBarrelPositionedCells.hits.Path = hcalBarrelCellsName + createHcalBarrelPositionedCells.positionedHits.Path = hcalBarrelPositionedCellsName + + # 3 - compute new cellID of cells based on new readout - removing row information + hcalBarrelCellsName2 = "HCalBarrelCells2" + rewriteHCalBarrel = RewriteBitfield("RewriteHCalBarrel", + # old bitfield (readout) + oldReadoutName=hcalBarrelReadoutName, + # specify which fields are going to be deleted + removeIds=["row"], + # new bitfield (readout), with new segmentation + newReadoutName=hcalBarrelReadoutName2, + debugPrint=10, + OutputLevel=INFO) + # clusters are needed, with deposit position and cellID in bits + rewriteHCalBarrel.inhits.Path = hcalBarrelCellsName + rewriteHCalBarrel.outhits.Path = hcalBarrelCellsName2 + + # 4 - attach positions to the new cells + from Configurables import CellPositionsHCalBarrelPhiThetaSegTool + hcalBarrelPositionedCellsName2 = "HCalBarrelPositionedCells2" + cellPositionHcalBarrelTool2 = CellPositionsHCalBarrelPhiThetaSegTool( + "CellPositionsHCalBarrel2", + readoutName=hcalBarrelReadoutName2, + OutputLevel=INFO + ) + createHcalBarrelPositionedCells2 = CreateCaloCellPositionsFCCee( + "CreateHCalBarrelPositionedCells2", + OutputLevel=INFO + ) + createHcalBarrelPositionedCells2.positionsTool = cellPositionHcalBarrelTool2 + createHcalBarrelPositionedCells2.hits.Path = hcalBarrelCellsName2 + createHcalBarrelPositionedCells2.positionedHits.Path = hcalBarrelPositionedCellsName2 + + # createHcalEndcapCells = CreateCaloCells("CreateHcalEndcapCaloCells", + # doCellCalibration=True, + # calibTool=calibHcalEndcap, + # addCellNoise=False, + # filterCellNoise=False, + # OutputLevel=INFO) + # createHcalEndcapCells.hits.Path="HCalEndcapHits" + # createHcalEndcapCells.cells.Path="HCalEndcapCells" + +else: + hcalBarrelCellsName = "emptyCaloCells" + hcalBarrelPositionedCellsName = "emptyCaloCells" + hcalBarrelCellsName2 = "emptyCaloCells" + hcalBarrelPositionedCellsName2 = "emptyCaloCells" + cellPositionHcalBarrelTool = None + cellPositionHcalBarrelTool2 = None + +# Empty cells for parts of calorimeter not implemented yet +createemptycells = CreateEmptyCaloCellsCollection("CreateEmptyCaloCells") +createemptycells.cells.Path = "emptyCaloCells" + +# Produce sliding window clusters (ECAL only) +towers = CaloTowerToolFCCee("towers", + deltaThetaTower=4 * 0.009817477 / 4, deltaPhiTower=2 * 2 * pi / 1536., + ecalBarrelReadoutName=ecalBarrelReadoutName, + ecalEndcapReadoutName=ecalEndcapReadoutName, + ecalFwdReadoutName="", + hcalBarrelReadoutName="", + hcalExtBarrelReadoutName="", + hcalEndcapReadoutName="", + hcalFwdReadoutName="", + OutputLevel=INFO) +towers.ecalBarrelCells.Path = ecalBarrelPositionedCellsName +towers.ecalEndcapCells.Path = "ECalEndcapCells" +towers.ecalFwdCells.Path = "emptyCaloCells" + +towers.hcalBarrelCells.Path = "emptyCaloCells" +towers.hcalExtBarrelCells.Path = "emptyCaloCells" +towers.hcalEndcapCells.Path = "emptyCaloCells" +towers.hcalFwdCells.Path = "emptyCaloCells" + +# Cluster variables +windT = 9 +windP = 17 +posT = 5 +posP = 11 +dupT = 7 +dupP = 13 +finT = 9 +finP = 17 +# Minimal energy to create a cluster in GeV (FCC-ee detectors have to reconstruct low energy particles) +threshold = 0.040 + +createClusters = CreateCaloClustersSlidingWindowFCCee("CreateClusters", + towerTool=towers, + nThetaWindow=windT, nPhiWindow=windP, + nThetaPosition=posT, nPhiPosition=posP, + nThetaDuplicates=dupT, nPhiDuplicates=dupP, + nThetaFinal=finT, nPhiFinal=finP, + energyThreshold=threshold, + energySharingCorrection=False, + attachCells=True, + OutputLevel=INFO + ) +createClusters.clusters.Path = "CaloClusters" +createClusters.clusterCells.Path = "CaloClusterCells" + +correctCaloClusters = CorrectCaloClusters("correctCaloClusters", + inClusters=createClusters.clusters.Path, + outClusters="Corrected" + createClusters.clusters.Path, + numLayers=[12], + firstLayerIDs=[0], + lastLayerIDs=[11], + readoutNames=[ecalBarrelReadoutName], + # do not split the following line or it will break scripts that update the values of the corrections + upstreamParameters = [[0.03900891447361534, -4.322941016402328, -139.1811369546787, 0.498342628339746, -3.3545078429754813, -13.99996971344221]], + upstreamFormulas=[ + ['[0]+[1]/(x-[2])', '[0]+[1]/(x-[2])']], + # do not split the following line or it will break scripts that update the values of the corrections + downstreamParameters = [[-0.0027661744480442195, 0.006059143775380306, 0.9788596364251927, -1.4951749409378743, -0.08491999337012696, 16.017621428757778]], + downstreamFormulas=[ + ['[0]+[1]*x', '[0]+[1]/sqrt(x)', '[0]+[1]/x']], + OutputLevel=INFO + ) + +calibrateCaloClusters = CalibrateCaloClusters("calibrateCaloClusters", + inClusters=createClusters.clusters.Path, + outClusters="Calibrated" + createClusters.clusters.Path, + systemIDs=[4], + numLayers=[12], + firstLayerIDs=[0], + readoutNames=[ + ecalBarrelReadoutName], + layerFieldNames=["layer"], + calibrationFile="./lgbm_calibration-CaloClusters.onnx", + OutputLevel=INFO + ) + +cellPositionHcalBarrelNoSegTool = None +cellPositionHcalExtBarrelTool = None + +# Output +out = PodioOutput("out", + OutputLevel=INFO) + +if runHCal: + out.outputCommands = ["keep *", "drop emptyCaloCells"] +else: + out.outputCommands = ["keep *", "drop HCal*", "drop emptyCaloCells"] + +if not saveCells: + out.outputCommands.append("drop ECal*Cells*") +if not saveClusterCells: + out.outputCommands.append("drop *ClusterCells*") +if not saveHits: + out.outputCommands.append("drop ECal*Hits*") + +out.filename = "electron_gun_10GeV_ALLEGRO_RECO.root" + + +# CPU information +chra = ChronoAuditor() +audsvc = AuditorSvc() +audsvc.Auditors = [chra] +createEcalBarrelCells.AuditExecute = True +createEcalBarrelPositionedCells.AuditExecute = True +if runHCal: + createHcalBarrelCells.AuditExecute = True +out.AuditExecute = True + +event_counter = EventCounter('event_counter') +event_counter.Frequency = 10 + +ExtSvc = [geoservice, podioevent, audsvc] +if dumpGDML: + ExtSvc += [gdmldumpservice] + +TopAlg = [ + event_counter, + input_reader, + createEcalBarrelCells, + createEcalBarrelPositionedCells, + resegmentEcalBarrel, + createEcalBarrelCells2, + createEcalBarrelPositionedCells2, + createEcalEndcapCells +] + +if runHCal: + TopAlg += [ + createHcalBarrelCells, + createHcalBarrelPositionedCells, + rewriteHCalBarrel, + createHcalBarrelPositionedCells2, + # createHcalEndcapCells + ] + +TopAlg += [ + createemptycells, + createClusters, +] + +if applyUpDownstreamCorrections: + TopAlg += [ + correctCaloClusters, + ] + +if applyMVAClusterEnergyCalibration: + TopAlg += [ + calibrateCaloClusters, + ] + calibrateCaloClusters.AuditExecute = True + +TopAlg += [ + out +] + +ApplicationMgr( + TopAlg=TopAlg, + EvtSel='NONE', + EvtMax=Nevts, + ExtSvc=ExtSvc, + StopOnSignal=True, +) diff --git a/full-detector-simulations/FCCeeGeneralOverview/run_IDEA_DIGI.py b/full-detector-simulations/FCCeeGeneralOverview/run_IDEA_DIGI.py new file mode 100644 index 00000000..296f86fb --- /dev/null +++ b/full-detector-simulations/FCCeeGeneralOverview/run_IDEA_DIGI.py @@ -0,0 +1,67 @@ +import os + +from Gaudi.Configuration import INFO, DEBUG + +# input +from Configurables import k4DataSvc, PodioInput +evtsvc = k4DataSvc('EventDataSvc') +evtsvc.input = "electron_gun_10GeV_IDEA_SIM.root" + +# Detector geometry (needed for the digitizer because we need to go in local coordinates) +from Configurables import GeoSvc +geoservice = GeoSvc("GeoSvc") +path_to_detector = os.environ.get("K4GEO", "") +print(path_to_detector) +detectors_to_use=[ + 'FCCee/IDEA/compact/IDEA_o1_v02/IDEA_o1_v02.xml' + ] +# prefix all xmls with path_to_detector +geoservice.detectors = [os.path.join(path_to_detector, _det) for _det in detectors_to_use] +geoservice.OutputLevel = INFO + +# Digitize drift chamber sim hits +from Configurables import DCHsimpleDigitizerExtendedEdm +dch_digitizer = DCHsimpleDigitizerExtendedEdm("DCHsimpleDigitizerExtendedEdm", + inputSimHits = "", + outputDigiHits = savetrackertool.SimTrackHits.Path.replace("sim", "digi"), + outputSimDigiAssociation = "DC_simDigiAssociation", + readoutName = "CDCHHits", + xyResolution = 0.1, # mm + zResolution = 1, # mm + debugMode = False, + OutputLevel = INFO +) + + +# Output +from Configurables import PodioOutput +out = PodioOutput("out", + OutputLevel=INFO) +out.outputCommands = ["keep *"] +out.filename = "output_simplifiedDriftChamber_MagneticField_"+str(magneticField)+"_pMin_"+str(momentum*1000)+"_MeV"+"_ThetaMinMax_"+str(thetaMin)+"_"+str(thetaMax)+"_pdgId_"+str(pdgCode)+"_stepLength_default.root" + +#CPU information +from Configurables import AuditorSvc, ChronoAuditor +chra = ChronoAuditor() +audsvc = AuditorSvc() +audsvc.Auditors = [chra] +genAlg.AuditExecute = True +hepmc_converter.AuditExecute = True +out.AuditExecute = True + +from Configurables import EventCounter +event_counter = EventCounter('event_counter') +event_counter.Frequency = 1 + +from Configurables import ApplicationMgr +ApplicationMgr( + TopAlg = [ + event_counter, + dch_digitizer, + out + ], + EvtSel = 'NONE', + EvtMax = 100, + ExtSvc = [geoservice, podioevent, , audsvc], + StopOnSignal = True, + )