MIT tutorial first draft ready (#144)

HEP-FCC · Mar 21, 2024 · 156ac16 · 156ac16
1 parent 9099580
commit 156ac16
Show file tree

Hide file tree

Showing 2 changed files with 95 additions and 22 deletions.
diff --git a/full-detector-simulations/FCCeeGeneralOverview/FCCeeGeneralOverview.md b/full-detector-simulations/FCCeeGeneralOverview/FCCeeGeneralOverview.md
@@ -15,7 +15,7 @@ So, let's start playing with Full Sim!
 ssh -X [email protected]
 # set-up the Key4hep environment
 source /cvmfs/sw-nightlies.hsf.org/key4hep/setup.sh
-# create a repository with the tutorial
+# create a repository for the tutorial
 mkdir FCC_Full_Sim_Tutorial
 cd FCC_Full_Sim_Tutorial
 git clone https://github.com/HEP-FCC/fcc-tutorials
@@ -56,7 +56,7 @@ k4run CLDReconstruction.py --inputFiles ../../fcc-tutorials/full-detector-simula
 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
+This created an edm4hep ROOT file with a bunch of new DIGI/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
@@ -69,7 +69,7 @@ NB: this step also produces a file named `wzp6_ee_mumuH_ecm240_CLD_RECO_aida.roo
 
 ### 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:
+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 [PODIO](https://github.com/key4hep/key4hep-tutorials/blob/main/edm4hep_analysis/edm4hep_api_intro.md#reading-edm4hep-files). The following very simple script already does a decent job:
 
 <!-- Explain a bit the script -->
 ```
@@ -169,10 +169,11 @@ Now let's plot the energy resolution for raw clusters and MVA calibrated cluster
 
 ```
 python plot_calo_energy_resolution.py electron_gun_10GeV_ALLEGRO_RECO.root
-# display with your favorite png renderer
+display electron_gun_10GeV_ALLEGRO_RECO_clusterEnergyResolution.png
+
 ```
 
-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.
+Look at both distributions to see how the MVA calibration improves the performance (both in terms of response and in terms of resolution). NB: we are in the middle of a transition for the Geant4 interface and we are here using the new workflow (ddsim) which is not yet fully validated. For instance, the low efficiency is under investigation.
 
 ### Changing Liquid Argon to Liquid Krypton
 
@@ -187,7 +188,7 @@ And to accomodate this change, we have to update the sampling fraction in the st
 
 ```
 vim run_ALLEGRO_RECO.py
-# comment out line 105 uncomment line 106
+# comment out line 105 and 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.
 
@@ -196,15 +197,38 @@ 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
+display electron_gun_10GeV_ALLEGRO_LKr_RECO_calibratedClusterEnergyResolution.png
 ```
 
-See how the resolution changed.
+You can see that the energy resolution did not improve. This can be due to multiple reasons: the derivation of the sampling fractions was done with little statistics for the LKr option, the MVA calibration was not re-trained, we plot only 100 events, ... But at least now you know how to change some of the detector free parameters :-)
 
 
 ## Towards IDEA tracking with the detailed Drift Chamber
 
+The IDEA tracking system Full Sim description is getting complete (we are just missing the Silicon Wrapper which will arrive soon) and we now have to design tracking algorithms. In this exercise we will produce a dataset containing Vertex and Drift Chamber digitized hits and which can therefore be used to devlop tracking algorithms.
+
+Let's run the IDEA simulation and digitization: 
+
 ```
 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
+ k4run run_IDEA_DIGI.py --EventDataSvc.input="electron_gun_10GeV_IDEA_SIM.root" --out.filename="electron_gun_10GeV_IDEA_DIGI.root"
+```
+
+The following collections will be useful for tracking:
+```
+CDCHDigis                       extension::DriftChamberDigi
+VTXDDigis                       edm4hep::TrackerHit3D
+VTXIBDigis                      edm4hep::TrackerHit3D
+VTXOBDigis                      edm4hep::TrackerHit3D
+```
+Note that the drift chamber digitized hit has two positions (on the left or on the right of the wire) because the only accessible information from a drift chamber is the distance to the wire which is degenerated. This left/right ambiguity is alleviated in the tracking algorithm.
+
+You can now play with the digitized hits:
+```
+root electron_gun_10GeV_IDEA_DIGI.root
+events->Draw("leftHitSimHitDeltaDistToWire") // this is the difference of the distance to the wire between the left digi and the sim hit (only produced in debug mode)
+events->Draw("CDCHDigis.rightPosition.x:CDCHDigis.rightPosition.y:CDCHDigis.rightPosition.z") // you can see that there was no magnetic field in this simulation and you can see the secondaries created inside the tracking volume
+events->Draw("CDCHHits.eDep/CDCHHits.pathLength", "CDCHHits.eDep/CDCHHits.pathLength < 4e-7") // this shows the de/dx from the simHits
 ```
 
 
@@ -246,4 +270,3 @@ This tool is useful but not perfect and will not meet all the needs (especially
  - [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/run_IDEA_DIGI.py b/full-detector-simulations/FCCeeGeneralOverview/run_IDEA_DIGI.py
@@ -1,13 +1,15 @@
 import os
 
-from Gaudi.Configuration import INFO, DEBUG
+from Gaudi.Configuration import *
 
-# input
+# Loading the input SIM file
 from Configurables import k4DataSvc, PodioInput
 evtsvc = k4DataSvc('EventDataSvc')
 evtsvc.input = "electron_gun_10GeV_IDEA_SIM.root"
+inp = PodioInput('InputReader')
 
-# Detector geometry (needed for the digitizer because we need to go in local coordinates)
+################## Simulation setup
+# Detector geometry
 from Configurables import GeoSvc
 geoservice = GeoSvc("GeoSvc")
 path_to_detector = os.environ.get("K4GEO", "")
@@ -19,34 +21,78 @@
 geoservice.detectors = [os.path.join(path_to_detector, _det) for _det in detectors_to_use]
 geoservice.OutputLevel = INFO
 
-# Digitize drift chamber sim hits
+# digitize vertex hits
+from Configurables import VTXdigitizer
+import math
+innerVertexResolution_x = 0.003 # [mm], assume 5 µm resolution for ARCADIA sensor
+innerVertexResolution_y = 0.003 # [mm], assume 5 µm resolution for ARCADIA sensor
+innerVertexResolution_t = 1000 # [ns]
+outerVertexResolution_x = 0.050/math.sqrt(12) # [mm], assume ATLASPix3 sensor with 50 µm pitch
+outerVertexResolution_y = 0.150/math.sqrt(12) # [mm], assume ATLASPix3 sensor with 150 µm pitch
+outerVertexResolution_t = 1000 # [ns]
+
+vtxib_digitizer = VTXdigitizer("VTXIBdigitizer",
+    inputSimHits = "VTXIBCollection",
+    outputDigiHits = "VTXIBDigis",
+    detectorName = "Vertex",
+    readoutName = "VTXIBCollection",
+    xResolution = innerVertexResolution_x, # mm, r-phi direction
+    yResolution = innerVertexResolution_y, # mm, z direction
+    tResolution = innerVertexResolution_t,
+    forceHitsOntoSurface = False,
+    OutputLevel = INFO
+)
+
+vtxob_digitizer = VTXdigitizer("VTXOBdigitizer",
+    inputSimHits = "VTXOBCollection",
+    outputDigiHits = "VTXOBDigis",
+    detectorName = "Vertex",
+    readoutName = "VTXOBCollection",
+    xResolution = outerVertexResolution_x, # mm, r-phi direction
+    yResolution = outerVertexResolution_y, # mm, z direction
+    tResolution = outerVertexResolution_t, # ns
+    forceHitsOntoSurface = False,
+    OutputLevel = INFO
+)
+
+vtxd_digitizer  = VTXdigitizer("VTXDdigitizer",
+    inputSimHits = "VTXDCollection",
+    outputDigiHits = "VTXDDigis",
+    detectorName = "Vertex",
+    readoutName = "VTXDCollection",
+    xResolution = outerVertexResolution_x, # mm, r direction
+    yResolution = outerVertexResolution_y, # mm, phi direction
+    tResolution = outerVertexResolution_t, # ns
+    forceHitsOntoSurface = False,
+    OutputLevel = INFO
+)
+
+# digitize drift chamber hits
 from Configurables import DCHsimpleDigitizerExtendedEdm
 dch_digitizer = DCHsimpleDigitizerExtendedEdm("DCHsimpleDigitizerExtendedEdm",
-    inputSimHits = "",
-    outputDigiHits = savetrackertool.SimTrackHits.Path.replace("sim", "digi"),
+    inputSimHits = "CDCHHits",
+    outputDigiHits = "CDCHDigis",
     outputSimDigiAssociation = "DC_simDigiAssociation",
     readoutName = "CDCHHits",
     xyResolution = 0.1, # mm
     zResolution = 1, # mm
-    debugMode = False,
+    debugMode = True,
     OutputLevel = INFO
 )
 
-
-# Output
+################ 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"
+
+out.filename = "electron_gun_10GeV_IDEA_DIGI.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
@@ -56,12 +102,16 @@
 from Configurables import ApplicationMgr
 ApplicationMgr(
     TopAlg = [
+			  inp,
               event_counter,
+              vtxib_digitizer,
+              vtxob_digitizer,
+              vtxd_digitizer,
               dch_digitizer,
               out
               ],
     EvtSel = 'NONE',
-    EvtMax   = 100,
-    ExtSvc = [geoservice, podioevent, , audsvc],
+    EvtMax   = -1,
+    ExtSvc = [geoservice, evtsvc, audsvc],
     StopOnSignal = True,
  )