diff --git a/CMakeLists.txt b/CMakeLists.txt index 2d070546..cbd5b919 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -55,6 +55,7 @@ endfunction() add_subdirectory(DCHdigi) add_subdirectory(ARCdigi) add_subdirectory(VTXdigi) +add_subdirectory(VTXdigi_Modular) add_subdirectory(VTXdigiDetailed) set(CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake ${CMAKE_MODULE_PATH}) add_subdirectory(Tracking) diff --git a/Makefile b/Makefile index aba113ff..b934efce 100644 --- a/Makefile +++ b/Makefile @@ -11,6 +11,25 @@ make: printf "#!/bin/bash\nif [ -n \"\$$KEY4HEP_STACK\" ];\nthen\n echo '----> Info: Key4hep stack already set up. Skipping...'\nelse\n source /cvmfs/sw-nightlies.hsf.org/key4hep/setup.sh\nfi\nexport LD_LIBRARY_PATH=${CURDIR}/install/lib:${CURDIR}/install/lib64:\$$LD_LIBRARY_PATH\nexport PYTHONPATH=${CURDIR}/install/python:\$$PYTHONPATH\n" > ${CURDIR}/setup.sh ; \ chmod +x ${CURDIR}/setup.sh +.PHONY: install +install: + cd build ; \ + make install -j8 ; \ + cd .. ; \ + export LD_LIBRARY_PATH=${CURDIR}/install/lib:${CURDIR}/install/lib64:$$LD_LIBRARY_PATH ; \ + export PYTHONPATH=${CURDIR}/install/python:$$PYTHONPATH ; \ + printf "#!/bin/bash\nif [ -n \"\$$KEY4HEP_STACK\" ];\nthen\n echo '----> Info: Key4hep stack already set up. Skipping...'\nelse\n source /cvmfs/sw-nightlies.hsf.org/key4hep/setup.sh\nfi\nexport LD_LIBRARY_PATH=${CURDIR}/install/lib:${CURDIR}/install/lib64:\$$LD_LIBRARY_PATH\nexport PYTHONPATH=${CURDIR}/install/python:\$$PYTHONPATH\n" > ${CURDIR}/setup.sh ; \ + chmod +x ${CURDIR}/setup.sh + +debug: + cd build ; \ + make install -j8 -d ; \ + cd .. ; \ + export LD_LIBRARY_PATH=${CURDIR}/install/lib:${CURDIR}/install/lib64:$$LD_LIBRARY_PATH ; \ + export PYTHONPATH=${CURDIR}/install/python:$$PYTHONPATH ; \ + printf "#!/bin/bash\nif [ -n \"\$$KEY4HEP_STACK\" ];\nthen\n echo '----> Info: Key4hep stack already set up. Skipping...'\nelse\n source /cvmfs/sw-nightlies.hsf.org/key4hep/setup.sh\nfi\nexport LD_LIBRARY_PATH=${CURDIR}/install/lib:${CURDIR}/install/lib64:\$$LD_LIBRARY_PATH\nexport PYTHONPATH=${CURDIR}/install/python:\$$PYTHONPATH\n" > ${CURDIR}/setup.sh ; \ + chmod +x ${CURDIR}/setup.sh + .PHONY: clean clean: @ (rm -r build install && rm setup.sh) || true diff --git a/README.md b/README.md index 6a36922d..d1b8265a 100644 --- a/README.md +++ b/README.md @@ -74,6 +74,7 @@ make get_data * `ARCdigi`: ARC digitization (for now, this step produces 'reco' collection) * `VTXdigi`: vertex detector digitization (for now, this step produces 'reco' collection) * `VTXdigiDetailed`: vertex detector and silicon sensors tracker detector digitization with detailed charge readout (for now, this step produces 'reco' collection) +* `VTXdigi_Modular`: silicon pixel digitizer, shares charge deposition across neighboring pixels. Supports different charge-sharing implementations, mainly lookup table-based implementation for now. Envisioned to be merged with the drift-bases implementation from `VTXdigiDetailed` in the future. (for now, this step produces 'reco' collection) * `Tracking`: tracking algorithms orchestrating [GenFit](https://github.com/GenFit/GenFit) ## Execute Examples diff --git a/VTXdigi_Modular/CMakeLists.txt b/VTXdigi_Modular/CMakeLists.txt new file mode 100644 index 00000000..e7de86b4 --- /dev/null +++ b/VTXdigi_Modular/CMakeLists.txt @@ -0,0 +1,56 @@ +set(PackageName VTXdigi_Modular) + +project(${PackageName}) + +file(GLOB sources + ${PROJECT_SOURCE_DIR}/src/*.cpp +) + +file(GLOB headers + ${PROJECT_SOURCE_DIR}/include/*.h +) + +gaudi_add_module(${PackageName} + SOURCES ${sources} + LINK + Gaudi::GaudiKernel + EDM4HEP::edm4hep + k4FWCore::k4FWCore + k4FWCore::k4Interface + DD4hep::DDRec +) + +target_include_directories(${PackageName} PUBLIC + $ + $ +) + +set_target_properties(${PackageName} PROPERTIES PUBLIC_HEADER "${headers}") + +file(GLOB scripts + ${PROJECT_SOURCE_DIR}/test/*.py +) + +file(COPY ${scripts} DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/test) + +install(TARGETS ${PackageName} + EXPORT ${CMAKE_PROJECT_NAME}Targets + RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" COMPONENT bin + LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib + PUBLIC_HEADER DESTINATION "${CMAKE_INSTALL_INCLUDEDIR}/@{CMAKE_PROJECT_NAME}" COMPONENT dev +) + +install(FILES ${scripts} DESTINATION test) + +set(test_name "test_runVTXdigi_Modular") +set(test_script "${CMAKE_CURRENT_SOURCE_DIR}/test/test.sh") +add_test(NAME ${test_name} + COMMAND bash ${test_script} + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} +) +set_tests_properties(${test_name} PROPERTIES + ENVIRONMENT "SOURCE_DIR_TEST=${CMAKE_CURRENT_SOURCE_DIR}/test;K4RECTRACKER=${CMAKE_INSTALL_PREFIX}/share/k4RecTracker;PATH=${CMAKE_INSTALL_PREFIX}/bin:$ENV{PATH};PYTHONPATH=${CMAKE_INSTALL_PREFIX}/python:${CMAKE_CURRENT_SOURCE_DIR}/test:$ENV{PYTHONPATH};LD_LIBRARY_PATH=${CMAKE_INSTALL_PREFIX}/lib:${CMAKE_INSTALL_PREFIX}/lib64:$ENV{LD_LIBRARY_PATH}" +) +if(COMMAND set_test_env) + set_test_env(${test_name}) +endif() diff --git a/VTXdigi_Modular/README.md b/VTXdigi_Modular/README.md new file mode 100644 index 00000000..6d52b9e8 --- /dev/null +++ b/VTXdigi_Modular/README.md @@ -0,0 +1,175 @@ +# VTXdidi_Modular: Silicon pixel digitizer +Developed by Jona Dilg (jona.dilg@cern.ch), Armin Ilg. 2026. + +## Description +The vertex digitizer creates digiHits from simHits by `TrackerHitPlane` by simulating the deposition and subsequent collection of charges in the sensor. The exact implementation of charge collection is user-definable. + +### Output +Produces EDM4hep `TrackerHitPlane` hits (referred to as digiHits), either per pixel hit or per cluster. For each `TrackerHitPlane`, creates a `TrackerHitSimTrackerHitLink` to each `SimTrackerHit` that contributed charge to any involved pixel. + +### Charge Collector Implementations +The algorithm allows a choice of different charge collectors, defined by the `ChargeCollectionMethod` Gaudi property. + +Each charge collector implements a method to distribute a simHit's deposited charge to the pixels around that simHit. + +#### `SinglePixel` +Simply fills the total charge deposited by a simHit into the pixel that the simHit position lies in. This will always produce a cluster size of 1 pix/clst and the binary sensor resolution expectation of pitch/sqrt(12) (at vertical incidence angles). + +#### `Debug` +Takes the pixel that the simHit position lies in as the central pixel. Deposited charge is distributed to pixels in an L-shape: +- 50% charge to the central pixel (index `[i_u,i_v]`) +- 30% charge to the pixel to the right of the central pixel (index `[i_u+1,i_v]`) +- 10% charge each to the two pixels above the central pixel (indices `[i_u,i_v+1]` and `[i_u,i_v+2]`) +This is intended to be used to test the correct `u`/`v` orientation of the pixel + +#### `LookupTable` +Shares charge among a hits surrounding pixels according to a lookup table. These lookup tables are generated from detailed simulations performed in sensor R&D. + +- **Requires** Gaudi property `LookupTableFile`. + +- A validation of this algorithm on test beam data is in progress. + +- An exemplary, Gaussian based LUT is placed in the examples folder. + +- Further information: + - This talk in the DRD3 WG4 (Simulation) Meeting contains more information on the implementation and preliminary results: https://indico.cern.ch/event/1658032/contributions/6968509/attachments/3239139/5776904/2026-03-16_DRD3-WG4.pdf + - A previous version of the talk was also held in the FCC Full Sim Working Meeting: https://indico.cern.ch/event/1613709/contributions/6814309/attachments/3190168/5677333/2025-12-10_FCC-FullSym-WorkingMeeting-3.pdf + + +## Gaudi Properties +- `SubDetectorName` - Name of the subdetector (eg. `Vertex`) +- `SubDetectorChildName` - Name of the subdetector child (eg. `VertexBarrel`), if applicable. If undefined, the subdetector itself is assumed to contain layers as children. +- `Clusterize` - Enable clustering of pixel hits. Defaults to `True` + - `True` - Clusterize pixel hits. Creates an EDM4hep `TrackerHitPlane` for each cluster. The hit position is the center of gravity of the cluster. The timestamp is the earliest timestamp among simHits that contributed to this cluster (possibly smeared by a gaussian according to `TimeSmearing`). + - `False` - Generate an EDM4hep `TrackerHitPlane` for each pixel hit, without clustering. The hit position is the pixel center. The timestamp is the earliest timestamp among the contributing simHits (possibly smeared). +- `ChargeCollectionMethod` - Select the charge collection implementation. Defaults to `SinglePixel` (for now) +- `Layers` - Select which layers of the (sub)detector to be digitized (eg. `[0, 1, 2]`). Defaults to digitizing all layers of the subdetector. +- `Threshold` - Sets a pixel threshold. After charges in an event have been collected in pixels, all pixels with `collected charges < threshold` are ignored. In terms of e-. Defaults to 0. +- `ChargeSmearing` - Sets the sigma of the Gaussian smearing that is applied to each pixels collected charge. In terms of e-. Defaults to 0, where no smearing is applied. Note that Pixels that do not collect any charge in an event are not smeared (eg. this can not be used to generate random pixel noise, only to smear the resolution on charge that is actually being collected), for performance. +- `TimeSmearing` - Sets the sigma of the Gaussian smearing that is applied to each pixel hit's timestamp before clustering / creating per-pixel digiHits. **UNIT???**. Defaults to 0. +- `ClusterPositionUncertainty` - Each digiHit has an estimate of it's spatial resolution uncertainty. To first order, this is simply the sensors spatial resolution in `u` and `v`. Yet, the spatial resolution changes with cluster shape and charge spread across the cluster (and with particle angle). Estimating a digiHits spatial resolution is important for the tracking algorithm, where a good estimation of the spatial resolution of each cluster allows for more precise refitting. Three different methods for estimating the spatial resolution are implemented (defaults to `[]`): + - Setting `[sigma_u, sigma_v]` simply applies the same spatial resolution to every digiHit. This is the most basic option. + - Setting `[]` does a basic charge-weighted uncertainty estimation, acting as an upper limit on the spatial resolution of a given digiHit. The charge-weighted uncertainty estimation is based on the assumption that single-pixel clusters have a resolution of pitch/sqrt(12), while larger clusters improve the resolution. This is, in turn, based on the assumption that the hit-distribution inside single-pixel clusters is completely flat, which does not hold for sensors with any significant charge sharing (such as small-pitch TPSCo 65nm CIS). With strong charge sharing, only hits close to the center of a pixel produce single-pixel clusters, while pixels closer to the edge produce multi-pixel clusters. This biases the resolution of single-pixel clusters, improving it. Because of the assumption, this method calculates an upper limit on a clusters spatial resolution. Importantly, this upper limit does not necessarily represent the dependence on the cluster size well, and can introduce a bias with better spatial resolution in multi-pixel clusters that might be unphysical. + - Setting `[sigma_u1, sigma_u2, sigma_u3, sigma_v1, sigma_v2, sigma_v3]` applies different residuals based on the cluster lengths in `u` and `v`. Clusters with a length above 3 will be assigned the `sigma_i3` as well. As of now, this is the most realistic estimate. This is only a good assumption if the particle incidence angle is close to vertical, as the cluster length is strongly dependent on the incidence angle. **Needs to be tested for shallow angles.** + +- `DebugHistograms` - Set to `True` to enable producing debugging histograms with the `Gaudi__Histograming__Sink__Root` service. Can decrease performance significantly, but is useful for debugging. Defaults to `False`. Include the following in the options file to produce the histograms file: + ``` + from Configurables import Gaudi__Histograming__Sink__Root as RootHistoSink + rootHistSvc = RootHistoSink("RootHistoSink") + rootHistSvc.FileName = "OutputFile.root" + # add rootHistSvc to the ExtSvc ;ist in ApplicationMgr + ``` +- `InfoPrintInterval` - Set the interval at which the event number is printed to the `INFO` debug level at the start of the event. In terms of events. Defaults to 100. + +- `LookupTableFile` - Lookup table file, necessary for the `LookupTable` charge collector. +- `LookupTableSegmentStepLength` - In the `LookupTable` charge collector implementation, a particles path through the sensor is sampled in smalls steps, sharing charges among surrounding pixels for each step. This sets that step length. Should be smaller than the LUT binning. In mm. Defaults to 0.0005 mm. +- `LookupTableIgnorePitch` - Ignore the sensor thickness and pixel pitch values stored in the LUT file. Useful for slightly stretching/shrinking the LUT to fit curved sensors where the sensor length is not an integer multiple of the pixel pitch. Defaults to `False`, where an error is thrown if the values from the detector geometry do not match the LUT file. If set to `True`, only a warning is printed is in case of a mismatch. + +--- +## Steering File Changes +The default steering files have very coarse (but fast) Geant4 settings that are not sufficient to use the full capability of precise digitisation. The following settings are necessary to achieve precise simulations: +```python +TODO: + +``` + +*Make sure to apply these changes properly. Simply pasting them into the top o a default steering file will overwrite these with the existing settings that come later. Best is to use a minimal steering file and rely on the ddsim defaults.* + +--- +## Example Options File +This is an example of the necessary changes to an options file. Can be used with the IDEA option 1, version 03 (o1v03) options file [run_digi_reco.py](https://github.com/HEP-FCC/FCC-config/blob/main/FCCee/FullSim/IDEA/IDEA_o1_v03/run_digi_reco.py). + +**Note** that in IDEA o1v03, the inner 3 VTX layers `[0,1,2]` have ARCADIA 20um pitch sensors, while the outer 2 VTX layers `[3,4]` have ATLASpix sensors. The digitizer can only handle one sensor type, so we digitize only the inner VTX. We would need another instance of the digitizer to digitize the outer VTX. + +```python +from Configurables import VTXdigi_Modular + +innervtxb_digitizer = VTXdigi_Modular( + "VTXBdigitizer", + + SimTrackHitCollectionName = ["VertexBarrelCollection"], + HeaderName = ["EventHeader"], + SimTrkHitRelationsCollection = ["VTXBSimDigiLinks"], + TrackerHitCollectionName = ["VTXBDigis"], + SubDetectorName = "Vertex", + SubDetectorChildName = "VertexBarrel", + + Layers = [0,1,2], + ChargeCollectionMethod = "LookupTable", + LookupTableFile = "PATH/TO/LOOKUPTABLE.init", + ClusterPositionUncertainty = [], + Threshold = 100, # in e + ChargeSmearing = 0, # in e + TimeSmearing = 0, # UNIT??? + + DebugHistograms = True, + OutputLevel = INFO) +``` +To create the debug histograms, you need the corresponding service: +```python +from Configurables import Gaudi__Histograming__Sink__Root as RootHistoSink + +rootHistSvc = RootHistoSink( + "RootHistoSink", + + FileName = "PATH/TO/WRITE/OUTPUTFILE.root" +) +``` +And finally, you need to add these services to the application manager: +```python +application_mgr = ApplicationMgr( + ... + TopAlg = [ + ... + innervtxb_digitizer, + ... + ], + ExtSvc = [ + ... + rootHistSvc, + ... + ], + ... +) +``` + + +--- +## Maintainers notes +- **Reference frames** + - The global detector reference frame uses `x,y,z`, `z` points along the beams and `y` points upwards + - The local sensor frame is `u,v,w` where `u,v` span the sensor plane and `w` is the sensors normal vector. For barrel sensors, `v` is typically parallel to `z` +- **Vectors** can be given as + 1) `dd4hep::rec::Vector3D` - fully featured vector, overloads operators `*+-` etc + 2) `edm4hep::Vector3d` - natively used by edm4hep (where simHit, digiHit are from) + -> generally use `dd4hep::rec::Vector3D`, convert via `ConvertVector()` where `edm4hep::Vector3d` is needed +- **Hit naming** scheme + - `simTrackerHit` refers to a `edm4hep::SimTrackerHit` + - `simHit` is always a `VTXdigi_tools::SimHitWrapper` (which contains a edm4hep::SimTrackerHit and some pre-computed info like surface and layer number) + - `digiHit` refers to `edm4hep::TrackerHitPlane` (the output of this digitizer) +- There are some crude tool tests for the functions in `VTXdigi_tools`. These are not exhaustive, and were only used early on in refactoring the code. Can be executed by calling `VTXdigi_tools::ToolTest()` in `VTXdigi_Modular::initialize()`. Not recommended. +### Units +- All **lengths** are in mm (edm4hep already does this) + BUT: dd4hep uses cm internally, so convert when passing values to/from dd4hep via `dd4hep::mm = 0.1` + Conversion: [A] = cm, [B] = mm + - `A = dd4hep::mm * B` + - `B = 1/dd4hep::mm * A` +- `Energies` are given in keV, but deposited charge is always handled in terms of electron-hole pairs (3.65 eV per eh-pair). This is much more common in sensor R&D. + +### Preliminary UML diagram +Universal Modelling Language (UML) diagram of the VTXdigi_Modular code structure. + + +### Futher improvements +The code has a number of `TODO:` and `FIXME:` marked. Larger items are: +- *(priority)* **Title** -- Description +- *(medium)* **Hit position uncertainty estimation** -- So far, a clusters spatial resolution estimation is done either via the cluster length in *u* and *v*, or via basic charge-weighted uncertainty estimation. Both of these options are sub-optimal and produce results much worse than what is possible. Optimally, all cluster shape and charge sharing information would be taken into account in calculating a clusters position and position resolution. + - An ML model trained on truth information might be the easiest way to go, here. + - Alternatively, an eta-function based approach can be used to correct the hit position beyond simple center-of-gravity by applying an empirically-determined function that accounts for non-linear charge sharing effects. The numeric eta function can be determined from the LUT by projecting it onto the *uv*-plane and measuring the charge-sharing ratio between the central and neighboring pixels. With this, the charge-weighted uncertainty estimation can be improved to account for charge sharing even in single-pixel clusters. *(Talk to reconstruction experts about how this is and might be done for a better understanding of what the digitiser should do)* + - Another option is to assume all hits come as straight lines from the IP. This is used to reconstruct each hits incidence angle to the sensor. A lookup table describing the resolution in dependence on the incidence angle is supplied, from which each hits position uncertainty is calculated. These lookup tables are generated from ddsim simulations, where MIPs (eg. 10 GeV muons) are shot at a single sensor plane at different angles. +- *(medium)* **Transfer propagation-based digitiser** -- Transfer the propagation-based digitizer currently implemented in `VTXDigi_Detailed` to a `ChargeCollector` implementation. Making this possible is the reason this specific architecture was chosen. Not urgent, but will make maintaining the repo a bit simpler. +- *(medium)* **Include some LUTs** -- We hope to include some LUTs for the Standard, n-Blanket and n-Gap layout of the TPSCo 65nm CIS chips in the repo. This is likely not a problem with NDAs, but has to be checked back with seniors. +- *(low)* **N-bit charge information** -- Currently charge information is completely analogue (collected charge per pixel is stored as float). Enable purely digital (N=1) or semi-digital (N>1) readout where the charge information is stored in N bits. This is much closer to how charge information is handled in real systems. +- *(low)* **Timing information in LUT** -- Currently, LUTs do not store any information on the signal shape charges in a given voxel will produce on the electrodes. Timestamps are simply taken from the simHits and smeared. Yet, charges from different voxels can take drastically different times (TPSCo 65nm CIS: O(ns)) to collect. Thus, the timestamps and amplitude of pixel hits might be changed. To simulate this, each LUT entry would need to encode not only the amount of shared charge, but also a parametrisation of the pixel response function. This is computationally expensive, and not implemented in Allpix Squared either. An exact algorithm on how to implement this time encoding is not yet clear. +- *(low)* **Gaussian-smearing based digitizer** -- Implementing such an algorithm here would simplify repository. Because the charge collector acts on pixels (and not on positions like the Gaussian-smearing), this is non-trivial and would require a check in the beginning of the event loop. +- diff --git a/VTXdigi_Modular/UML.svg b/VTXdigi_Modular/UML.svg new file mode 100644 index 00000000..87b2b722 --- /dev/null +++ b/VTXdigi_Modular/UML.svg @@ -0,0 +1,3449 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/VTXdigi_Modular/include/IChargeCollector.h b/VTXdigi_Modular/include/IChargeCollector.h new file mode 100644 index 00000000..9a84d6fe --- /dev/null +++ b/VTXdigi_Modular/include/IChargeCollector.h @@ -0,0 +1,28 @@ +// VTXdigi_Modular/include/IChargeCollector.h +#pragma once + +#include "TGeoMatrix.h" + +struct VTXdigi_Modular; + +namespace VTXdigi_tools { + + class SimHitWrapper; // forward-declare things in include/VTXdigi_tools.h + class HitMap; + +class IChargeCollector { +public: + virtual ~IChargeCollector() = default; + virtual void FillHit(const SimHitWrapper& simHit, HitMap& hitMap, const TGeoHMatrix& trafoMatrix) const = 0; + float GetChargeCollectionDepthCenter() const { return m_chargeCollectionDepthCenter; } + +protected: + explicit IChargeCollector(const VTXdigi_Modular& digitizer) : m_digitizer(digitizer) {} + + const VTXdigi_Modular& m_digitizer; + float m_chargeCollectionDepthCenter=0; // Defines the vertical center of the charge collection region in the sensitive volume. Needed for correct digiHit positions (and residual plots) with LUTs where the charge collection varies along their depth (like TPSCo 65nm CIS). in mm, wrt. the sensor local w coordinate (w=0 at center of sensitive volume) +}; + +std::unique_ptr CreateChargeCollector(const VTXdigi_Modular& digitizer, const std::string& algorithm); + +} // namespace VTXdigi_tools \ No newline at end of file diff --git a/VTXdigi_Modular/include/VTXdigi_Modular.h b/VTXdigi_Modular/include/VTXdigi_Modular.h new file mode 100644 index 00000000..89f2e30f --- /dev/null +++ b/VTXdigi_Modular/include/VTXdigi_Modular.h @@ -0,0 +1,352 @@ +// VTXdigi_Modular/include/VTXdigi_Modular.h +#pragma once + +#include "IChargeCollector.h" +#include "VTXdigi_tools.h" + +// GAUDI +#include "GAUDI_VERSION.h" +#include "Gaudi/Property.h" +#include "Gaudi/Accumulators/Histogram.h" +#include "GaudiKernel/IRndmGenSvc.h" + +// K4FWCORE +#include "k4Interface/IGeoSvc.h" +#include "k4FWCore/Transformer.h" +#include "k4Interface/IUniqueIDGenSvc.h" + +// EDM4HEP +#include "edm4hep/SimTrackerHitCollection.h" +#include "edm4hep/EventHeaderCollection.h" +#include "edm4hep/TrackerHitPlaneCollection.h" +#include "edm4hep/TrackerHitSimTrackerHitLinkCollection.h" + +// DD4HEP +#include "DDRec/SurfaceManager.h" + +/** @class VTXdigi_Modular + * + * Creates trackerHits from simHits. Produces clusters from simHits, outputs either the cluster centre or all hits in the cluster as digitized hits. + * + * @author Jona Dilg, Armin Ilg + * @date 2025-09 + */ + +/* -- Forward declarations -- */ +namespace VTXdigi_tools { + class IChargeCollector; +} + +struct VTXdigi_Modular final : k4FWCore::MultiTransformer (const edm4hep::SimTrackerHitCollection&, const edm4hep::EventHeaderCollection&)> { + + VTXdigi_Modular(const std::string& name, ISvcLocator* svcLoc); + + StatusCode initialize() override; + StatusCode finalize() override; + + std::tuple operator() (const edm4hep::SimTrackerHitCollection& simHits, const edm4hep::EventHeaderCollection& headers) const override; + + /** @brief Increment histograms. To be called from the charge collector, once per simHit */ + void FillHistograms_fromChargeCollector_perSimHit(const int layer, const dd4hep::rec::Vector3D& pathTravel, const float pathLength_Geant4, const dd4hep::rec::Vector3D& truthPos_local, const TGeoHMatrix& trafoMatrix, const bool createdInGenerator) const; + + /* -- Accessors for charge collector -- */ + + inline std::array ActiveVolumeDimensions() const { return {m_sensorLength.first, m_sensorLength.second, m_sensorActiveThickness}; } + + inline std::pair PixelPitch() const { return m_pixelPitch; } + + inline std::pair PixelCount() const { return m_pixelCount; } + + inline float Threshold() const { return m_threshold; } + + inline bool LUT_ignorePitch() const { return m_LUT_ignorePitch.value(); } + + inline bool LUT_shiftTruthPos() const { return m_LUT_shiftTruthPos.value(); } + + /** @brief Draw a random number for charge smearing + * FIXME: this is not thread safe, but I don't know how this is done in Gaudi (while retaining thread safety & reproducibility with a given seed). + */ + inline float DrawChargeSmearingNumber() const { return static_cast(m_rndm_charge()); } + + inline std::string LutFileName() const { return m_LUT_FileName; } + inline float LutStepLength() const { return m_LUT_stepLength; } + +private: + + /* ---- Initialization & finalization functions ---- */ + + void InitServicesAndGeometry(); + void InitLayersAndSensors(); + void InitHistograms(); + + void PrintCountersSummary() const; + + /* ---- Core algorithm functions ---- */ + + /** @brief Check the event before starting digitization (eg. if there even are any simHits) */ + bool CheckEventSetup(const edm4hep::SimTrackerHitCollection& simHits, const edm4hep::EventHeaderCollection& headers) const; + + /** @brief Check if a simHit is in a layer that is configured for digitization */ + bool CheckSimhitLayer(const edm4hep::SimTrackerHit& simHit) const; + + /** @brief Smear collected charge on each pixel, and remove pixels that do not pass the threshold. */ + void ApplyNoiseAndThreshold(VTXdigi_tools::HitMap& hitMap) const; + + /** @brief Clusterize all pixelHits on a hitMap + * @note Uses either next-neightbor clustering or no clustering, depending on the value of the Gaudi property "Clusterize" */ + std::vector Clusterize(const VTXdigi_tools::HitMap& hitMap) const; + + /** @brief Create a digiHit per cluster */ + void CreateDigiHits(edm4hep::TrackerHitPlaneCollection& digiHits, edm4hep::TrackerHitSimTrackerHitLinkCollection& digiHitLinks, const dd4hep::DDSegmentation::CellID& cellID, const TGeoHMatrix& trafoMatrix, const std::vector& clusters) const; + + void FillHistograms_perSimHit(const VTXdigi_tools::SimHitWrapper& hit) const; + void FillHistograms_perPixel(const dd4hep::DDSegmentation::CellID& cellID, const VTXdigi_tools::Pixel& pix, const std::pair clusterPos_local) const; + void FillHistograms_perDigiHit(const std::unordered_set& simHits, const edm4hep::TrackerHitPlane& digiHit, const TGeoHMatrix& trafoMatrix, const int clusterSize, const int clusterSize_u, const int clusterSize_v) const; + void FillHistograms_perSensor(const std::vector& simHits, const edm4hep::TrackerHitPlaneCollection& digiHits, const TGeoHMatrix& trafoMatrix, const dd4hep::DDSegmentation::CellID& cellID) const; + + /* -- Properties -- */ + + const std::string m_undefinedString = "UNDEFINED"; + + Gaudi::Property m_subDetName{this, "SubDetectorName", m_undefinedString, "Name of the subdetector (eg. \"Vertex\")"}; + Gaudi::Property m_subDetChildName{this, "SubDetectorChildName", m_undefinedString, "Name of the subdetector child (eg. \"VertexBarrel\"), if applicable. If undefined, the subdetector itself is assumed to contain layers as children."}; + + Gaudi::Property m_geoServiceName{this, "GeoSvcName", "GeoSvc", "The name of the GeoSvc instance"}; + Gaudi::Property m_encodingStringVariable{this, "EncodingStringParameterName", "GlobalTrackerReadoutID", "The name of the DD4hep constant that contains the Encoding string for tracking detectors"}; + Gaudi::Property m_clusterize{this, "Clusterize", true, "Whether to clusterize hits or not. If false, each pixel hit is output as a separate trackerHit. If true, the cluster centre and total charge are output."}; + + + /* -- Properties mainlyrelated to the main event loop -- */ + + Gaudi::Property> m_layers{this, "Layers", {}, "Which layers to run on (0-indexed). If empty, all layers are run."}; + + /* -- Properties and members related to the various charge collection algorithms-- */ + + Gaudi::Property m_chargeCollectionMethod{this, "ChargeCollectionMethod", "Drift", "Method used for charge collection: \"Fast\", \"Drift\", \"LookupTable\", etc."}; + Gaudi::Property m_threshold{this, "Threshold", 0.0f, "Pixel threshold for firing (in e-)."}; + Gaudi::Property> m_positionUncertainty{this, "ClusterPositionUncertainty", {}, "Sensor spatial resolution in u and v direction (in mm). Used for the position uncertainty in digiHits"}; + Gaudi::Property m_smearing_charge{this, "ChargeSmearing", 0.0f, "Gaussian smearing to be applied to a pixels collected charge (in e-). Applied after charge collection but before thresholding. If 0, no noise is applied. Quadratically add pixel noise and threshold smearing if necessary."}; + Gaudi::Property m_smearing_time{this, "TimeSmearing", 0.0f, "Gaussian smearing to be applied to a pixels time (in ns). Applied to the digiHits time stamp. If 0, no time smearing is applied."}; + + Gaudi::Property m_debugHistograms{this, "DebugHistograms", false, "Flag to create and fill debug histograms. Not recommended for multithreading, might lead to crashes. Default is false."}; + Gaudi::Property m_infoPrintInterval{this, "InfoPrintInterval", 100, "Interval for printing information during processing."}; + + /* LUT */ + Gaudi::Property m_LUT_FileName{this, "LookupTableFile", "", "File to load the lookup table from. Must be given if ChargeCollectionMethod is set to \"LookupTable\"."}; + Gaudi::Property m_LUT_stepLength{this, "LookupTableSegmentStepLength", 0.0005f, "Length of the segments that a particle path through the sensor is split into. The deposited charge is distributed evenly over the segments, and each segments charge is distributed according to the in-pixel bin the segment center falls into. In mm. Defaults to 0.0005 mm."}; + Gaudi::Property m_LUT_ignorePitch{this, "LookupTableIgnorePitch", false, "Ignore the sensor thickness and pixel pitch values stored in the LUT file. Useful for slightly stretching/shrinking the LUT to fit curved sensors where the sensor length is not an integer multiple of the pixel pitch. If empty, the LUT file values are used."}; + Gaudi::Property m_LUT_shiftTruthPos{this, "LookupTableShiftTruthPosition", false, "Internally shift the truth position of the simHit. Only affects the output histograms, does not affect any collection. If turned to false, angled particle trajectories will bias the residual plots in case of LUT tables with uneven charge collection across the sensor thickness."}; + + /* -- Services, geometry variables -- */ + + SmartIF m_randomService; + SmartIF m_geoService; + std::unique_ptr m_cellIdDecoder; + const dd4hep::Detector* m_detector = nullptr; + const dd4hep::rec::SurfaceMap* m_surfaceMap; // map from cellID (unsigned long, without segmentation bits) to simSurface (dd4hep::rec::ISurface*) + dd4hep::VolumeManager m_volumeManager; // volume manager to get the physical cell sensitive volume + dd4hep::DetElement m_subDetector; // subdetector DetElement. contains layers as children + + /* -- Member variables -- */ + + std::unique_ptr m_chargeCollector = nullptr; + + std::pair m_pixelCount = {0, 0}; + std::pair m_pixelPitch = {0.0f, 0.0f}; + float m_sensorActiveThickness = 0.0f; // also in mm + float m_inactiveMaterialAbove = 0.0f; // in mm, inactive material above the active volume in sensor normal direction + float m_inactiveMaterialBelow = 0.0f; + + std::pair m_sensorLength = {0.0f, 0.0f}; + TGeoRotation m_sensorNormalRotation = TGeoRotation("sensorNormalRotation"); // rotation to rotate the sensor local coordinate system. Initialised to unit matrix. + + Rndm::Numbers m_rndm_charge; // TODO: Is this multithreading safe? + Rndm::Numbers m_rndm_time; + + /* -- Counters -- */ + + mutable Gaudi::Accumulators::Counter<> m_counter_eventsRead{this, "Events read"}; + mutable Gaudi::Accumulators::Counter<> m_counter_eventsRejected_noSimHits{this, " - Events rejected (no simHits)"}; + mutable Gaudi::Accumulators::Counter<> m_counter_eventsAccepted{this, " = Events accepted"}; + + mutable Gaudi::Accumulators::Counter<> m_counter_simHitsRead{this, "SimTrackerHits read"}; + mutable Gaudi::Accumulators::Counter<> m_counter_simHitsRejected_LayerNotToBeDigitized{this, " - SimTrackerHits rejected (layer not to be digitized)"}; + mutable Gaudi::Accumulators::Counter<> m_counter_simHitsAccepted{this, " = SimTrackerHits accepted"}; + mutable Gaudi::Accumulators::Counter<> m_counter_accSimHitsCreatedInGenerator{this, "( accepted SimTrackerHits from particles created in generator)"}; + mutable Gaudi::Accumulators::Counter<> m_counter_accSimHitsDirectMcParticleLink{this, "( accepted SimTrackerHits with direct MC particle link)"}; + + + + mutable Gaudi::Accumulators::Counter<> m_counter_digiHitsCreated{this, "Digi hits created"}; + + /* -- Histograms -- */ + + enum { + hist1d_simHitE, + hist1d_simHitCharge, + hist1d_simHitMomentum_keV, + hist1d_simHitMomentum_MeV, + hist1d_simHitMomentum_GeV, + hist1d_simHit_x, + hist1d_simHit_y, + hist1d_simHit_z, + hist1d_simHit_z_createdInGenerator, + hist1d_simHit_z_createdInSimulation, + hist1d_simHit_Vertex_x, + hist1d_simHit_Vertex_y, + hist1d_simHit_Vertex_z, + hist1d_simhit_MomentumDirection_x, + hist1d_simhit_MomentumDirection_y, + hist1d_simhit_MomentumDirection_z, + hist1d_simHit_InitialMomentumDirection_x, + hist1d_simHit_InitialMomentumDirection_y, + hist1d_simHit_InitialMomentumDirection_z, + hist1d_digiHitCharge, + hist1d_simHitPDG, + hist1d_digiHitsPerSimHit, + hist1d_clusterSize, + hist1d_clusterSize_u, + hist1d_clusterSize_v, + hist1d_clusterSize_createdInGenerator, + hist1d_clusterSize_createdInSimulation, + hist1d_residual_u, + hist1d_residual_u_maxEParticleOnSensor, + hist1d_residual_u_noParents, + hist1d_residual_u_directMcParticleLink, + hist1d_residual_u_mcOriginFarFromSimHit, + hist1d_residual_u_directMcParticleLink_mcOriginFarFromSimHit, + hist1d_residual_u_singlePixelCluster, + hist1d_residual_u_multiPixelCluster, + hist1d_residual_u_createdInGenerator, + hist1d_residual_u_createdInSimulation, + hist1d_residual_v, + hist1d_residual_v_maxEParticleOnSensor, + hist1d_residual_v_noParents, + hist1d_residual_v_directMcParticleLink, + hist1d_residual_v_mcOriginFarFromSimHit, + hist1d_residual_v_directMcParticleLink_mcOriginFarFromSimHit, + hist1d_residual_v_singlePixelCluster, + hist1d_residual_v_multiPixelCluster, + hist1d_residual_v_createdInGenerator, + hist1d_residual_v_createdInSimulation, + hist1d_residual_w, + hist1d_residual_w_noParents, + hist1d_residual_r, + hist1d_clusterPosUncertainty_u, + hist1d_clusterPosUncertainty_u_noParents, + hist1d_clusterPosUncertainty_v, + hist1d_clusterPosUncertainty_v_noParents, + hist1d_pixelDistToClusterCenter_u, + hist1d_pixelDistToClusterCenter_v, + hist1d_pathTravel_u, + hist1d_pathTravel_v, + hist1d_pathTravel_r, + hist1d_simHitTimeStamp, + hist1d_digiHitTimeStamp, + hist1d_simHitsPerDigiHit, + hist1d_highestEnergyParticleOnSensor_energy, + hist1dArrayLen + }; + mutable std::unordered_map< + int, // layer number + std::array< + std::unique_ptr< + Gaudi::Accumulators::StaticHistogram< + 1, + Gaudi::Accumulators::atomicity::full, + float + > + >, + hist1dArrayLen + > + > m_hist1d; + + enum{ + histProfile1d_digiHitCharge_vs_global_z, + histProfile1d_clusterSize_vs_global_z, + histProfile1d_clusterSize_u_vs_global_z, + histProfile1d_clusterSize_v_vs_global_z, + histProfile1d_residual_u_vs_global_z, + histProfile1d_residual_v_vs_global_z, + histProfile1d_residual_r_vs_global_z, + histProfile1d_pathTravel_u_vs_global_z, + histProfile1d_pathTravel_v_vs_global_z, + histProfile1d_pathTravel_r_vs_global_z, + histProfile1dArrayLen }; + mutable std::unordered_map< + int, // layer number + std::array< + std::unique_ptr< + Gaudi::Accumulators::StaticProfileHistogram< + 1, + Gaudi::Accumulators::atomicity::full, + float + > + >, + histProfile1dArrayLen + > + > m_histProfile1d; + + enum { + hist2d_digiHitCharge_vs_global_z, + hist2d_hitMap_simHits, + hist2d_hitMap_simHits_createdInGenerator, + hist2d_hitMap_simHits_createdInSimulation, + hist2d_hitMap_simHits_eDepAboveThreshold, + hist2d_hitMap_pixelHits, + hist2d_clusterSize_vs_global_z, + hist2d_clusterSize_u_vs_global_z, + hist2d_clusterSize_v_vs_global_z, + hist2d_clusterSize_vs_global_z_createdInGenerator, + hist2d_clusterSize_vs_global_z_createdInSim, + hist2d_residual_u_vs_global_z, + hist2d_residual_v_vs_global_z, + hist2d_residual_r_vs_global_z, + hist2d_residual_vs_clusterPosUncertainty_u, + hist2d_residual_vs_clusterPosUncertainty_v, + hist2d_pathTravel_u_vs_global_z, + hist2d_pathTravel_u_vs_global_z_createdInGenerator, + hist2d_pathTravel_u_vs_global_z_createdInSimulation, + hist2d_pathTravel_v_vs_global_z, + hist2d_pathTravel_v_vs_global_z_createdInGenerator, + hist2d_pathTravel_v_vs_global_z_createdInSimulation, + hist2d_pathTravel_w_vs_global_z, + hist2d_pathTravel_w_vs_global_z_createdInGenerator, + hist2d_pathTravel_w_vs_global_z_createdInSimulation, + hist2d_pathTravel_r_vs_global_z, + hist2d_simHit_xy, + hist2d_simHit_xz, + hist2d_simHit_yz, + hist2dArrayLen + }; + mutable std::unordered_map< + int, // layer number + std::array< + std::unique_ptr< + Gaudi::Accumulators::StaticHistogram< + 2, + Gaudi::Accumulators::atomicity::full, + float + > + >, + hist2dArrayLen + > + > m_hist2d; + + enum { + hist1dglobal_pathTravel_r, + hist1dglobal_pathTravel_r_Geant4, + hist1dglobal_pathTravel_r_ratio, + hist1dglobalArrayLen + }; + std::array< + std::unique_ptr< + Gaudi::Accumulators::StaticHistogram< + 1, + Gaudi::Accumulators::atomicity::full, + float + > + >, + hist1dglobalArrayLen + > m_hist1dglobal; +}; // class VTXdigi_Modular diff --git a/VTXdigi_Modular/include/VTXdigi_tools.h b/VTXdigi_Modular/include/VTXdigi_tools.h new file mode 100644 index 00000000..afde3001 --- /dev/null +++ b/VTXdigi_Modular/include/VTXdigi_tools.h @@ -0,0 +1,228 @@ +// VTXdigi_Modular/include/VTXdigi_tools.h +#pragma once + +#include "GaudiKernel/GaudiException.h" +#include "GaudiKernel/RndmGenerators.h" + +#include "DDRec/Surface.h" + +#include "edm4hep/SimTrackerHit.h" + +#include +#include +#include +#include + +namespace VTXdigi_tools { + +constexpr float kChargePerkeV = 273.97f; // in electrons, for silicon (1 eh-pair ~ 3.65 eV) + +/* -- Tool tests -- */ +bool ToolTest(); + +/* -- SimHitWrapper -- */ + +/** @brief Check if a MCParticle was created in the generator */ +bool CreatedInGenerator(const edm4hep::MCParticle& mcParticle); +/** @brief Check if a simTrackerHit was created in the generator */ +inline bool CreatedInGenerator(const edm4hep::SimTrackerHit& simTrackerHit) { return CreatedInGenerator(simTrackerHit.getParticle()); }; + +/** @brief Class to contain all information about a simTrackerHit that is needed for the digitization. + * @note this is where the simTrackerHit is actually stored, everything else (pixelHit / cluster) will store pointers to this. */ +class SimHitWrapper { + edm4hep::SimTrackerHit m_simTrackerHit; + dd4hep::DDSegmentation::CellID m_cellID; // without segmentation bits + float m_charge; + int m_layerNumber; + mutable dd4hep::rec::Vector3D m_truthPos; // simHit truth position, local coordinates. Mutable because it might be adjusted in const ChargeCollector::FillHit() to account for charge collection biases (see ChargeCollector_impl.h ChargeCollector_LUT::MoveTruthPosition() for more info). + +public: + SimHitWrapper(edm4hep::SimTrackerHit simTrackerHit, const std::unique_ptr& cellIdDecoder); + SimHitWrapper(const SimHitWrapper& other) = default; + SimHitWrapper(SimHitWrapper&& other) = default; + SimHitWrapper() = default; + + /** @brief Set the truth position of the simHit in local coordinates */ + inline void SetTruthPos(const dd4hep::rec::Vector3D& pos) const { m_truthPos = pos; } // only used for histogramming after filling the hits, so not really a problem that this is mutable + + friend void swap(SimHitWrapper& a, SimHitWrapper& b) noexcept; + inline const edm4hep::SimTrackerHit* hitPtr() const { return &m_simTrackerHit; } + + /** @brief Access the truth position of the simHit in local coordinates + * @note might have been adjusted by ChargeCollector::FillHit() to correct for charge collection effects (eg. not completely depleted sensors where charges are only collected close to the upper sensor surface as in TPSCo 65nm CIS). */ + inline const dd4hep::rec::Vector3D truthPos() const { return m_truthPos; } + + inline dd4hep::DDSegmentation::CellID cellID() const { return m_cellID; } + inline float charge() const { return m_charge; } + inline int layer() const { return m_layerNumber; } + + /** @brief Check if the simHit particle was created in the generator (via MCParticle, which might not be the "real" particle, if SaveAllParticles is turned off in ddsim) */ + inline bool CreatedInGenerator() const { return ::VTXdigi_tools::CreatedInGenerator(m_simTrackerHit); } + + /** @brief Check if the linked MC particle is the particle that produced this hit. If false, the hit was produced by a secondary particle, which had no MC particle */ + inline bool hasDirectMcParticleLink() const { return !m_simTrackerHit.isProducedBySecondary(); } +}; + +void swap(SimHitWrapper& a, SimHitWrapper& b) noexcept; + +/** @brief Find simHits from different MCParticles. If two simHits originate from the same MCParticle (eg. have shared parents), return only the one from the MCParticle further up the family tree. */ +// std::unordered_set FindSimHitsWithIndividualParents(const std::vector& simHits); + +/* -- Pixel -- */ + +/** @brief A pixel in the hit map. Can have multiple contributing simHits */ +struct Pixel { + float charge; + std::unordered_set simHits; + std::pair index; // This info is saved in (a) the map key, and (b) here inside the Pixel object. This is inefficient. But it makes the code a bit nicer not having to pass the index around separately. + + Pixel(std::pair pix) : charge(0.f), index(pix) { + simHits.reserve(2); // avoid too many reallocations, will rarely see more than 2 simTrackerHits contributing to the same pixel + } + Pixel() : charge(0.f), index({-1, -1}) { + simHits.reserve(2); // avoid too many reallocations, will rarely see more than 2 simTrackerHits contributing to the same pixel + } +}; + +/* -- Cluster -- */ + +/** @brief Holds clusters: charge, pointers to all contributing pixels, pointers to all contributing simHits */ +struct Cluster { + std::vector pixels; // raw pointer into HitMap::Pixels. Valid as long as the HitMap is not modified after clustering. See HitMap::ComputeClusters() for details. + std::unordered_set simHits; + float charge = 0.f; + + inline int GetSize() const { return pixels.size(); }; + int GetSize(const int axis) const; // axis = 0 for u, 1 for v. + + /** @brief Compute the center position of a cluster via charge-weighed center of gravity in terms of pixel indices */ + std::pair ComputePos() const; + + /** @brief Compute the uncertainty of the cluster position via charge-weighed center of gravity in terms of pixel indices */ + std::pair ComputePosUncertainty_ChargeWeighted() const; + std::pair ComputePosUncertainty_ChargeWeighted(const std::pair& clusterPos) const; +}; + +/** @brief Get the indices of all direct neighbors of a pixel */ +std::array, 4> GetDirectNeighbors(const std::pair& i_uv); + +/* -- HitMap -- */ + +struct Hash_PairInt { + size_t operator()(const std::pair& i_uv) const noexcept { + return (static_cast(i_uv.first) << 32) ^ static_cast(i_uv.second); + } +}; + +using PixelMap = std::unordered_map, Pixel, Hash_PairInt>; + +/** @brief HitMap of all pixel hits on a sensor + * @note uses a std::unordered_map to only store pixels that have charge, which is more memory efficient for large pixel counts and low occupancy. */ +class HitMap { + /* I tried implementing a vector that contains every pixel, but for large pixel counts this is very memory-inefficient. Instead, this class uses a std::map to only store pixels that have charge. This is more memory efficient for sparse hits, with O(1) simHit/sensor/event this is at least a factor 100 faster that the vector approach. Maybe not the case to ttbar run with O(0.1%) pixel occupancy. */ + + PixelMap m_pixels; // hit pixels, stored by value + std::pair m_pixCount; // size of the sensor in pixels + +public: + HitMap(std::pair pixCount); + + /** @brief Add charge and a simHit to a pixel */ + void FillCharge(std::pair i_uv, float charge, const SimHitWrapper& simHitWrapper); + + void ApplyChargeSmearing(const Rndm::Numbers& rndm_charge); + void ApplyThreshold(const float threshold); + + /** @brief Get one pixel's collected charge */ + float GetCharge(std::pair i_uv) const; + + /** @brief Get the total charge across all pixels */ + float GetTotalCharge() const; + + /** @brief Return a const reference to the internal pixel map */ + inline const PixelMap& Hits() const { return m_pixels; }; + + /** @brief Return the number of pixels with charge */ + inline int GetTotalPixelsWithCharge() const { return m_pixels.size(); }; + + /** @brief Clusterize hit pixels in the HitMap, using direct neighbors + * @note Returns Cluster objects whose `pixels` members hold raw pointers into m_pixels. This is safe because ComputeClusters() is const and called only after all FillCharge() insertions are complete -> the map will not rehash while the returned Clusters are alive. Do not call FillCharge() on this HitMap after calling ComputeClusters(). */ + std::vector ComputeClusters() const; + + /** @brief Return a vector of clusters, where each cluster is a single pixel */ + std::vector ComputeClusters_singePixels() const; + + inline void Reset() { m_pixels.clear(); }; + +private: + /** @brief Returns true if the pixel is out of bounds */ + inline bool _OutOfBounds(std::pair i_uv) const; +}; // class HitMap + +/* -- helpers -- */ + +/** @brief Convert a dd4hep::rec::Vector3D to edm4hep::Vector3d and vice-versa */ +dd4hep::rec::Vector3D ConvertVector(edm4hep::Vector3d vec); +/** @copydoc dd4hep::rec::Vector3D ConvertVector(edm4hep::Vector3d vec) */ +dd4hep::rec::Vector3D ConvertVector(edm4hep::Vector3f vec); +/** @copydoc dd4hep::rec::Vector3D ConvertVector(edm4hep::Vector3d vec) */ +edm4hep::Vector3d ConvertVector(dd4hep::rec::Vector3D vec); + +/** @brief Compute the transformation matrix from global detector to local sensor frame for a given sensor cellID */ +TGeoHMatrix ComputeSensorTrafoMatrix(const dd4hep::DDSegmentation::CellID& cellID, const dd4hep::VolumeManager& volumeManager, const TGeoRotation& sensorNormalRotation); + +/** @brief Transform a position from global detector coordinates to local sensor coordinates, using the sensor transformation matrix */ +dd4hep::rec::Vector3D GlobalToLocal(const dd4hep::rec::Vector3D& global, const TGeoHMatrix& M); +/** @brief Transform a position from local sensor coordinates to global detector coordinates, using the sensor transformation matrix */ +dd4hep::rec::Vector3D LocalToGlobal(const dd4hep::rec::Vector3D& local, const TGeoHMatrix& M); + +dd4hep::DDSegmentation::CellID GetCellID_short(const edm4hep::SimTrackerHit& simTrackerHit); +dd4hep::DDSegmentation::CellID GetCellID_short(const dd4hep::DDSegmentation::CellID& cellID); + +int GetLayer(const dd4hep::DDSegmentation::CellID& cellID, const std::unique_ptr& cellIdDecoder); +int GetLayer(const edm4hep::SimTrackerHit& simTrackerHit, const std::unique_ptr& cellIdDecoder); + +/* -- Binning tools -- */ + +/** @brief Given a histogram definition (x0, binWidth, nBins) and a value x, compute the bin index i in which x falls. + * @return Int, -1 if x is out of range. + * @note Bins are 0-indexed (vs ROOT's 1-indexing) */ +int ComputeBinIndex(float x, float binX0, float binWidth, int binN); + +/** @brief Given a binning definition (x0, binWidth, nBins), compute the center position of a given bin index i in the histogram + * @return Float, the center position of the bin + * @note Bins are 0-indexed (vs ROOT's 1-indexing) */ +float ComputeBinCenter(int i, float binX0, float binWidth); +/** @brief Given a binning definition (x0, x1, nBins), compute the center position of a given bin index i in the histogram + * @return Float, the center position of the bin + * @note Bins are 0-indexed (vs ROOT's 1-indexing) */ +float ComputeBinCenter(int i, float binX0, float binX1, int binN); + +/** @brief Compute the pixel indices (i_u, i_v) for a given (local) position inside the sensor */ +std::pair ComputePixelIndices(const dd4hep::rec::Vector3D& pos, const std::pair pixelPitch, const std::pair pixelCount); + +/** @brief Compute the in-pixel indices (j_u, j_v, j_w) for a given (local) position inside the pixel and layer index + * @note Assumption: each layer has only 1 type if sensor */ +std::array ComputeInPixelIndices(const dd4hep::rec::Vector3D& pos, const std::array& binCount, const std::pair& pixelPitch, const std::array& activeVolumeDimensions); + +/** @brief Compute the center position of a given pixel (i_u,i_v) in sensor-local coordinates (u,v,w) + * + * @note The w coordinate is set to depletedRegionDepthCenter. 0 for center, +25 for sensor surface, +20 for TPSCo 65nm maps. +*/ +dd4hep::rec::Vector3D ComputePosFromPixIndex_local(const std::pair pixelIndex, const std::pair sensorLength, const std::pair pixelPitch, float depletedRegionDepthCenter); +/** @brief Compute the center position of a given pixel (i_u,i_v) in sensor-local coordinates (u,v,0) */ +dd4hep::rec::Vector3D ComputePosFromPixIndex_local(const std::pair pixelIndex, const std::pair sensorLength, const std::pair pixelPitch); + +/** @brief Compute the position of a given (pixel-)index (i_u,i_v) in sensor-local coordinates (u,v,0) . + * @note index 0 indicates the center of the pixel, index -0.5 the lower edge and +0.5 the upper edge. + * @note Does not check if the position is within the sensor bounds! + * @note The w coordinate is set to depletedRegionDepthCenter. 0 for center, +25 for sensor surface, +20 for TPSCo 65nm maps. + */ +dd4hep::rec::Vector3D ComputePosFromPixIndex_local(const std::pair index, const std::pair sensorLength, const std::pair pixelPitch, float depletedRegionDepthCenter); +/** @brief Compute the position of a given (pixel-)index (i_u,i_v) in sensor-local coordinates (u,v,0) . + * @note index 0 indicates the center of the pixel, index -0.5 the lower edge and +0.5 the upper edge. + * @note Does not check if the position is within the sensor bounds! + */ +dd4hep::rec::Vector3D ComputePosFromPixIndex_local(const std::pair index, const std::pair sensorLength, const std::pair pixelPitch); + +} // namespace VTXdigi_tools \ No newline at end of file diff --git a/VTXdigi_Modular/src/ChargeCollector_impl.cpp b/VTXdigi_Modular/src/ChargeCollector_impl.cpp new file mode 100644 index 00000000..660d6ff7 --- /dev/null +++ b/VTXdigi_Modular/src/ChargeCollector_impl.cpp @@ -0,0 +1,542 @@ +// VTXdigi_Modular/src/ChargeCollector_impl.cpp + +#include "../src/ChargeCollector_impl.h" + +namespace VTXdigi_tools { +using ::VTXdigi_Modular; // "unqualified name introduction from global namespace" (just so I remember what to call this in C++ speak) +using ::endmsg; // makes the Copilot autocomplete work better + +std::unique_ptr CreateChargeCollector(const VTXdigi_Modular& digitizer, const std::string& algorithm) { + std::unique_ptr chargeCollector; + + if (algorithm == "LookupTable") { + chargeCollector = std::make_unique(digitizer); + } else if (algorithm == "Drift") { + throw std::runtime_error("ChargeCollector_Drift not implemented yet."); + } else if (algorithm == "Fast") { + throw std::runtime_error("ChargeCollector_Fast not implemented yet."); + } else if (algorithm == "SinglePixel") { + chargeCollector = std::make_unique(digitizer); + } else if (algorithm == "Debug") { + chargeCollector = std::make_unique(digitizer); + } + else { + throw std::runtime_error("Unknown ChargeCollector type: " + algorithm); + } + + digitizer.info() << " - Created charge collector with algorithm " << algorithm << ". Charge collection depth center at " << chargeCollector->GetChargeCollectionDepthCenter() << " mm." << endmsg; + return chargeCollector; + // if (algorithm == "LookupTable") { + // return std::make_unique(digitizer); + // } else if (algorithm == "Drift") { + // // return std::make_unique(digitizer); + // throw std::runtime_error("ChargeCollector_Drift not implemented yet."); + // } else if (algorithm == "Fast") { + // // return std::make_unique(); + // throw std::runtime_error("ChargeCollector_Fast not implemented yet."); + // } else if (algorithm == "SinglePixel") { + // return std::make_unique(digitizer); + // } else if (algorithm == "Debug") { + // return std::make_unique(digitizer); + // } + // throw std::runtime_error("Unknown ChargeCollector type: " + algorithm); +} + +bool ConstructPath(Path& path, const SimHitWrapper& simHit, const TGeoHMatrix& trafoMatrix, const VTXdigi_Modular& digitizer) { + const float eps = 1e-6f; // reasonable for number O(0.01) (like sensor thickness in mm) with float precision + + path.simPos = simHit.truthPos(); + + double momentum_global[3] = { + static_cast(simHit.hitPtr()->getMomentum().x), + static_cast(simHit.hitPtr()->getMomentum().y), + static_cast(simHit.hitPtr()->getMomentum().z) + }; + double momentum_local[3]; + trafoMatrix.MasterToLocalVect(momentum_global, momentum_local); + + /* Step 1 - travel vector */ + const double scaleFactor_travel = digitizer.ActiveVolumeDimensions().at(2) / std::abs(momentum_local[2]); + path.travel = scaleFactor_travel * dd4hep::rec::Vector3D(momentum_local[0], momentum_local[1], momentum_local[2]); + + /* Step 2 - entry point */ + if (std::abs(path.simPos.z()) > digitizer.ActiveVolumeDimensions().at(2)/2.f + eps) { + digitizer.warning() << "SimHit position is outside the sensor volume (local w = " << path.simPos.z() << " mm, sensor thickness = " << digitizer.ActiveVolumeDimensions().at(2) << " mm). This should never happen. Forcing it to w=0." << endmsg; + path.simPos.z() = 0.f; // ensures no divide by zero etc + } + float shiftDist_w; + if (path.travel.z() >= 0.f) { + shiftDist_w = path.simPos.z() + 0.5f * digitizer.ActiveVolumeDimensions().at(2); + } + else { + shiftDist_w = path.simPos.z() - 0.5f * digitizer.ActiveVolumeDimensions().at(2); + } + const float scaleFactor_entry = shiftDist_w / path.travel.z(); + path.entry = path.simPos - scaleFactor_entry * path.travel; + + /* Step 3 - clip path to sensor edges (in u/v) */ + std::pair t = std::make_pair(0.f, 1.f); // parametrize path as entry + t*travel; t in [0,1] + t = ComputePathClippingFactors(t, path.entry.x(), path.travel.x(), digitizer.ActiveVolumeDimensions().at(0)); + t = ComputePathClippingFactors(t, path.entry.y(), path.travel.y(), digitizer.ActiveVolumeDimensions().at(1)); + if (t.first != 0.f || t.second != 1.f) { + if (0.f <= t.first && t.first < t.second && t.second <= 1.f) { + /* valid clipping */ + digitizer.debug() << " - Clipping SimHitPath with t [" << t.first << ", " << t.second << "]. PathLength changed to " << static_cast((t.second - t.first) * path.travel.r()*1000) << " um from " << static_cast(path.travel.r()*1000) << " um" << endmsg; + + path.entry = path.entry + t.first * path.travel; + path.travel = (t.second - t.first) * path.travel; + } + else [[unlikely]] { + /* invalid clipping, shouldn't happen */ + digitizer.warning() << "VTXdigi_tools::Path::Path() - invalid clipping factors t = [" << t.first << ", " << t.second << "]. Path might lie completely outside the sensor." << endmsg; + digitizer.debug() << " -> entry (" << path.entry.x() << ", " << path.entry.y() << ", " << path.entry.z() << ") mm, exit (" << path.entry.x() + path.travel.x() << ", " << path.entry.y() + path.travel.y() << ", " << path.entry.z() + path.travel.z() << ") mm, sensor dim. (+-" << digitizer.ActiveVolumeDimensions().at(0)/2 << ", +-" << digitizer.ActiveVolumeDimensions().at(1)/2 << ") mm" << endmsg; + digitizer.debug() << " -> Path length " << static_cast(path.travel.r()*1000) << " um, in G4 " << static_cast(simHit.hitPtr()->getPathLength()*1000) << " um" << endmsg; + return false; + } + } + + /* Step 4 -check that path is not much longer than the length it had in Geant4 */ + path.lengthG4 = simHit.hitPtr()->getPathLength(); + if (path.travel.r() > kPathLengthTolerance * path.lengthG4) { + digitizer.debug() << " - Shortening path length from " << static_cast(path.travel.r()*1000) << " um to " << static_cast(path.lengthG4*1000) << " um (the respective path length in Geant4)." << endmsg; + + /* make sure the path stays centred around the simTrackerHit position */ + const float t_simPos = ( (path.simPos - path.entry).dot(path.travel) ) / (path.travel.r() * path.travel.r()); + + const float t_length_halved = 0.5f * path.lengthG4 / path.travel.r(); // length of the new path in terms of t [0,1] on old path, halved + const float t_center = std::max(t_length_halved, std::min(t_simPos, 1.f - t_length_halved)); // center of new path clamped to [t_length_half, 1 - t_length_half] while not exceeding [0,1] + + const float t_min = t_center - t_length_halved; + const float t_max = t_center + t_length_halved; + + path.entry = path.entry + t_min * path.travel; + path.travel = (t_max - t_min) * path.travel; + } + path.length = path.travel.r(); + + + digitizer.debug() << " - Constructed path, length " << path.travel.r()*1000 << " um (G4-length " << path.lengthG4*1000 << " um), entry (" << path.entry.x() << ", " << path.entry.y() << ", " << path.entry.z() << ") mm, exit (" << path.entry.x() + path.travel.x() << ", " << path.entry.y() + path.travel.y() << ", " << path.entry.z() + path.travel.z() << ") mm, " << endmsg; + return true; // indicate valid path constructed +} + +std::pair ComputePathClippingFactors(std::pair t, const float entry_ax, const float travel_ax, const float sensorLength_ax) { + /* only need the components that are parallel to the axis (u/v) that we are clipping */ + const bool positiveDir = travel_ax >= 0.f; // false -> path points in negative direction along this axis + + const float minPos = std::min(entry_ax, entry_ax + travel_ax); + if (minPos < -0.5f * sensorLength_ax) { + /* path extends out of sensor in negative direction*/ + + const float t_clip = (-minPos - 0.5f * sensorLength_ax) / std::abs(travel_ax); + if (positiveDir){ + t.first = std::max(t.first, t_clip); + } else { + t.second = std::min(t.second, 1-t_clip); + } + } + + const float maxPos = std::max(entry_ax, entry_ax + travel_ax); + if (maxPos > 0.5f * sensorLength_ax) { + const float t_clip = (maxPos - 0.5f * sensorLength_ax) / std::abs(travel_ax); + + if (positiveDir) { + t.second = std::min(t.second, 1-t_clip); + } else { + t.first = std::max(t.first, t_clip); + } + } + + return t; +} + + + +/* -- LUT approach -- */ + +LookupTable::LookupTable(const std::string& lutFileName, const VTXdigi_Modular& digitizer) { + digitizer.debug() << " - Constructing LUT from file \"" << lutFileName << "\"." << endmsg; + + /* parse the LUT in Allpix Squared format + * See https://indico.cern.ch/event/1489052/contributions/6475539/attachments/3063712/5418424/Allpix_workshop_Lemoine.pdf (slide 10) for more info on fields in the LUT file */ + + if (lutFileName.empty()) + throw std::runtime_error("VTXdigi_tools::LookupTable::LookupTable(): LUT file name is empty. A LUT file must be given to load the lookup table."); + + const int headerLines = 5; + + digitizer.debug() << " - Opening LUT file \"" << lutFileName << "\"." << endmsg; + std::ifstream lutFile(lutFileName); + if (!lutFile.is_open()) + throw std::runtime_error("VTXdigi_tools::LookupTable::LookupTable(): Could not open LUT file \"" + lutFileName + "\"."); + + std::string line; + int lineNumber = 0; + + /* Parse pixel-pitch, thickness, in-pixel bin count from header (all in 5th line) */ + for (; lineNumber < 5; ++lineNumber) + std::getline(lutFile, line); + std::istringstream headerStringStream(line); + std::string headerEntry; + std::vector headerLineEntries; + + while (std::getline(headerStringStream, headerEntry, ' ')) { + if (!headerEntry.empty()) + headerLineEntries.push_back(headerEntry); + } + + if (headerLineEntries.size() != 11) + throw std::runtime_error("VTXdigi_tools::LookupTable::LookupTable(): Invalid number of entries in LUT file in 5th header line: found " + std::to_string(headerLineEntries.size()) + " entries, expected 11."); + + for (int j=0; j<3; j++) { + m_binCount.at(j) = std::stoi(headerLineEntries.at(7+j)); + } + digitizer.debug() << " - found in-pixel bin count of (" << m_binCount.at(0) << ", " << m_binCount.at(1) << ", " << m_binCount.at(2) << ") from LUT file header." << endmsg; + + /* -> compare the values we just parsed to the values retrieved from the detector geometry */ + const float eps = 1e-7f; // reasonable for number O(0.01) (like sensor thickness in mm) with float precision + + const float sensorThickness = std::stof(headerLineEntries.at(0)) / 1000.f; // convert from um to mm + if (std::abs(sensorThickness - digitizer.ActiveVolumeDimensions().at(2)) > eps) { + if (!digitizer.LUT_ignorePitch()) { + throw std::runtime_error("VTXdigi_tools::LookupTable::LookupTable(): Sensor thickness mismatch between LUT file and detector geometry: LUT file specifies " + std::to_string(sensorThickness) + " mm, but geometry has " + std::to_string(digitizer.ActiveVolumeDimensions().at(2)) + " mm active volume thickness."); + } + else { + digitizer.warning() << "Sensor thickness mismatch between LUT file and detector geometry. LUT file: " << sensorThickness << "mm, geometry (active volume thickness): " << digitizer.ActiveVolumeDimensions().at(2) << "mm. Ignored because LookupTableIgnorePitch is set to true." << endmsg; + } + } + + const std::pair pitch = {std::stof(headerLineEntries.at(1)) / 1000.f, std::stof(headerLineEntries.at(2)) / 1000.f}; + if (std::abs(pitch.first - digitizer.PixelPitch().first) > eps || std::abs(pitch.second - digitizer.PixelPitch().second) > eps) { + if (!digitizer.LUT_ignorePitch()) + throw std::runtime_error("VTXdigi_tools::LookupTable::LookupTable(): Pixel pitch mismatch between LUT file and detector geometry: LUT file specifies (" + std::to_string(pitch.first) + ", " + std::to_string(pitch.second) + ") mm, but geometry has (" + std::to_string(digitizer.PixelPitch().first) + ", " + std::to_string(digitizer.PixelPitch().second) + ") mm."); + else + digitizer.warning() << "Pixel pitch mismatch between LUT file and detector geometry. LUT file: " << pitch.first << "mm, geometry: " << digitizer.PixelPitch().first << "mm. Ignored because LookupTableIgnorePitch is set to true." << endmsg; + } + + digitizer.debug() << " - Found matching pixel pitch and sensor thickness in LUT file." << endmsg; + + /* Get matrix size (5x5, 7x7, ...) from the length of the first line after the header */ + if (std::getline(lutFile, line)) { + m_matrixSize = static_cast(std::sqrt(std::count(line.begin(), line.end(), ' ') - 2)); // not very robust, but works for valid Allpix2 files. first 3 entries are bin indices + } + else { + throw std::runtime_error("VTXdigi_tools::LookupTable::LookupTable(): Could not read first line after header in LUT file: " + digitizer.LutFileName()); + } + if (m_matrixSize < 3 || m_matrixSize % 2 == 0) + throw std::runtime_error("VTXdigi_tools::LookupTable::LookupTable(): Matrix size must be an odd integer >= 3, but is " + std::to_string(m_matrixSize) + "."); + m_matrixSize_half = (m_matrixSize - 1) / 2; + digitizer.debug() << " - Inferred matrix size of " << m_matrixSize << " from first line." << endmsg; + + /* Set up the matrix vector */ + m_matrices.resize(m_binCount.at(0) * m_binCount.at(1) * m_binCount.at(2) * m_matrixSize * m_matrixSize, 0.f); + + /* set up mapping from Allpix2 LUT format + * (row-major, starts on bottom left) + * to the format expected by the LookupTable class + * (row-major, starts on top-left) */ + std::unordered_map indexMapping; // i: index in local format; indexMapping[i]: index in Allpix2 format + for (int i_u = 0; i_u < m_matrixSize; i_u++) { + for (int i_v = 0; i_v < m_matrixSize; i_v++) { + const int i_AP2 = i_u + (m_matrixSize - 1 - i_v) * m_matrixSize; + int i_VTXdigi = i_u + i_v * m_matrixSize; + indexMapping[i_VTXdigi] = i_AP2; + } + } + + /* Now we can finally start parsing the matrix values */ + lutFile.clear(); + lutFile.seekg(0, std::ios::beg); + for (int i=0; i matricesEntrySum_perWBin; + matricesEntrySum_perWBin.resize(m_binCount.at(2), 0.f); + float matricesEntrySum = 0.f; + + lineNumber = headerLines + 1; + while (std::getline(lutFile, line)) { + if (line.empty() || line[0] == '#') + throw std::runtime_error("VTXdigi_tools::LookupTable::LookupTable(): Empty or comment line found in LUT file at line " + std::to_string(lineNumber+1) + ". All lines (past the 5 header lines) must contain valid matrix data."); + + std::istringstream stringStream(line); + std::vector lineEntries; + std::string entryString; + + /* read the line & do sanity checks */ + while (std::getline(stringStream, entryString, ' ')) { + if (!entryString.empty()) + lineEntries.push_back(entryString); + } + + if (static_cast(lineEntries.size()) != 3 + m_matrixSize*m_matrixSize) + throw std::runtime_error("VTXdigi_tools::LookupTable::LookupTable(): Invalid number of entries in LUT file at line " + std::to_string(lineNumber+1) + ": found " + std::to_string(lineEntries.size()) + " entries, but expected " + std::to_string(3 + m_matrixSize*m_matrixSize) + " (3 for bin indices, " + std::to_string(m_matrixSize*m_matrixSize) + " for matrix values)."); + + /* First 3 entries are in-pixel binning indices */ + Index_inPix j_uvw({std::stoi(lineEntries[0])-1, std::stoi(lineEntries[1])-1, std::stoi(lineEntries[2])-1});// Allpix2 input is 1-indexed. Insane, I know. + + if (j_uvw.at(0) < 0 || j_uvw.at(0) >= m_binCount[0] || + j_uvw.at(1) < 0 || j_uvw.at(1) >= m_binCount[1] || + j_uvw.at(2) < 0 || j_uvw.at(2) >= m_binCount[2]) { + throw std::runtime_error("Invalid in-pixel bin indices in LUT file at line " + std::to_string(lineNumber+1) + ": got (" + std::to_string(j_uvw.at(0)) + ", " + std::to_string(j_uvw.at(1)) + ", " + std::to_string(j_uvw.at(2)) + "), but expected ranges are [0, " + std::to_string(m_binCount[0]-1) + "], [0, " + std::to_string(m_binCount[1]-1) + "], [0, " + std::to_string(m_binCount[2]-1) + "]."); + } + + /* Parse matrix values & set it */ + std::vector matrixEntries(m_matrixSize*m_matrixSize, 0.); + float matrixEntrySum = 0.f; + for (int i = 0; i < m_matrixSize*m_matrixSize; i++) { + float entry = std::stof(lineEntries[3 + indexMapping[i]]); // NaN check done on sum + if (entry < kLutEntryMinimum) + entry = 0.f; // avoid very small entries for performace + matrixEntries[i] = entry; + matrixEntrySum += entry; + } + digitizer.verbose() << " - Parsed matrix for in-pixel bin (" << j_uvw.at(0) << ", " << j_uvw.at(1) << ", " << j_uvw.at(2) << "), entry sum " << std::to_string(matrixEntrySum) << ", setting it now..." << endmsg; + if (std::isnan(matrixEntrySum)) + throw std::runtime_error("VTXdigi_tools::LookupTable::LookupTable(): Charge sharing matrix for in-pixel bin (" + std::to_string(j_uvw.at(0)) + "," + std::to_string(j_uvw.at(1)) + "," + std::to_string(j_uvw.at(2)) + ") contains NaN values (sum of entries is NaN)."); + + matricesEntrySum += matrixEntrySum; + matricesEntrySum_perWBin[j_uvw.at(2)] += matrixEntrySum; + SetMatrix(j_uvw, matrixEntries); + + lineNumber++; + } // loop over lines containing a matrix each + + if (lineNumber - (headerLines+1) != m_binCount[0] * m_binCount[1] * m_binCount[2]) + throw std::runtime_error("Invalid number of matrices loaded from file: expected " + std::to_string(m_binCount[0] * m_binCount[1] * m_binCount[2]) + " matrices (inferred from bin count in header) but found " + std::to_string(lineNumber - headerLines) + " lines."); + + const float matrixNumber = static_cast(lineNumber - headerLines); + + // From which w-level are charges collected? + float collectedFromW = 0.f; + for (int j_w = 0; j_w < m_binCount.at(2); j_w++) { + const float binHeight = sensorThickness / m_binCount.at(2); + const float w_bin_center = j_w*binHeight + 0.5f*binHeight - 0.5f*sensorThickness; + collectedFromW += matricesEntrySum_perWBin.at(j_w) * w_bin_center; + } + collectedFromW /= matricesEntrySum; + + if (std::abs(collectedFromW) > 0.5f * sensorThickness) throw std::runtime_error("Invalid charge collection depth inferred from LUT file: " + std::to_string(collectedFromW) + " mm. This is outside the sensor volume, which extends from " + std::to_string(-0.5f*sensorThickness) + " mm to " + std::to_string(0.5f*sensorThickness) + " mm."); + if (std::abs(collectedFromW) >= 1.e-5f) { + m_chargeCollectionDepthCenter = collectedFromW; // the member is initalised to 0.f + } + + digitizer.info() << " - Loaded lookup table from file. Matrices parsed: " << (lineNumber - headerLines) << ". Charge collected from sensitive volume: " << matricesEntrySum/matrixNumber*100 << " percent (rest is lost, eg recombination). The center of the collection region is at " << m_chargeCollectionDepthCenter << endmsg; +} + +void LookupTable::SetMatrix(const Index_inPix& j_uvw, const std::vector& weights) { + if (static_cast(weights.size()) != m_matrixSize*m_matrixSize) + throw std::runtime_error("VTXdigi_tools::LookupTable::SetMatrix: weights size (" + std::to_string(weights.size()) + ") does not match matrix size (" + std::to_string(m_matrixSize*m_matrixSize) + ")"); + + /* check if matrix is valid */ + float sum = 0.f; + for (int row = 0; row < m_matrixSize; ++row) { + for (int col = 0; col < m_matrixSize; ++col) { + sum += weights.at(row*m_matrixSize + col); + } + } + if (std::isnan(sum)) + throw std::runtime_error("VTXdigi_tools::LookupTable::SetMatrix: Charge sharing matrix for in-pixel bin (" + std::to_string(j_uvw.at(0)) + "," + std::to_string(j_uvw.at(1)) + "," + std::to_string(j_uvw.at(2)) + ") contains NaN values."); + if (sum < 0 || sum > 1.f + 1.e-5f) + throw std::runtime_error("VTXdigi_tools::LookupTable::SetMatrix: Charge sharing matrix for in-pixel bin (" + std::to_string(j_uvw.at(0)) + "," + std::to_string(j_uvw.at(1)) + "," + std::to_string(j_uvw.at(2)) + ") has a weight sum of " + std::to_string(sum) + ", but needs to lie in [0,1]."); + + for (int row = 0; row < m_matrixSize; ++row) { + for (int col = 0; col < m_matrixSize; ++col) { + /* weights are given in row-major order, starting at top left. + * We store charge sharing matrices in col-major order, starting at bottom left (lowest bin index) */ + m_matrices.at(FindIndex(j_uvw, col, row)) = weights.at((m_matrixSize-1-row)*m_matrixSize + col); + } + } +} + +void LookupTable::SetAllMatrices(const std::vector& weights) { + for (int j_u = 0; j_u < m_binCount.at(0); ++j_u) { + for (int j_v = 0; j_v < m_binCount.at(1); ++j_v) { + for (int j_w = 0; j_w < m_binCount.at(2); ++j_w) { + SetMatrix({j_u, j_v, j_w}, weights); + } + } + } +} + +int LookupTable::FindIndex (const Index_inPix& j, const int col, const int row) const { + if (j[0] < 0 || j[0] >= m_binCount[0] + || j[1] < 0 || j[1] >= m_binCount[1] + || j[2] < 0 || j[2] >= m_binCount[2] ) [[unlikely]] { + throw std::runtime_error("VTXdigi_tools::LookupTable::FindIndex: in-pix bin out of range"); + } + if (col < 0 || col >= m_matrixSize || row < 0 || row >= m_matrixSize) [[unlikely]] { + throw std::runtime_error("VTXdigi_tools::LookupTable::FindIndex: col or row out of range"); + } + + int index_matrix = j[0] + m_binCount[0] * (j[1] + m_binCount[1] * j[2]); + int index_element = col * m_matrixSize + row; + return index_matrix * m_matrixSize * m_matrixSize + index_element; +} + +ChargeCollector_LUT::ChargeCollector_LUT(const VTXdigi_Modular& digitizer) : IChargeCollector(digitizer), + m_LUT(digitizer.LutFileName(), digitizer), + m_stepLength(digitizer.LutStepLength()), + m_shiftTruthPos(digitizer.LUT_shiftTruthPos()) { + + /* LUT is constructed in place (from file) */ + m_chargeCollectionDepthCenter = m_LUT.GetChargeCollectionDepthCenter(); + + m_digitizer.info() << " - ChargeCollector_LUT constructed successfully." << endmsg; +} + +void ChargeCollector_LUT::FillHit(const SimHitWrapper& simHit, HitMap& hitMap, const TGeoHMatrix& trafoMatrix) const { + /* TODO: implement voxel-traversal instead of numerically splitting the charge. + I would start from the simple algo here: https://www.redblobgames.com/grids/line-drawing.html */ + + Path path; + if (!ConstructPath(path, simHit, trafoMatrix, m_digitizer)) [[unlikely]] + return; + + if (m_shiftTruthPos) { + MoveTruthPosition(simHit, path); // shifts the sim hit position to the depth in the sensor where most charge is collected, to get usesful residual plots. + } + + const int stepCount = std::max(1, static_cast(std::ceil(path.length / m_stepLength))); + const float segmentCharge = simHit.charge() / stepCount; + + Index_segment nextSeg, seg = ComputeSegmentIndices(0, stepCount, path); + int segmentsInBin = 1; + + /* TODO: currently, DistributeSegmentCharge is THE bottleneck. Optimise this by collecting all entries from each vector in a local size x size matrix, and copy that to m_LUT once all segments in a pixel have been filled. */ + + for (int i_next = 1; i_next < stepCount; ++i_next) { + nextSeg = ComputeSegmentIndices(i_next, stepCount, path); + + /* collect segments in same bin */ + if (seg == nextSeg) { + // m_digitizer.verbose() << " - Segment lies in the same pixel and in-pixel bin as previous segment, continuing." << endmsg; + ++segmentsInBin; + continue; + } + else { + DistributeSegmentCharge(hitMap, seg, segmentCharge, segmentsInBin, simHit); + seg = nextSeg; + segmentsInBin = 1; + } + } // loop over segments + + m_digitizer.FillHistograms_fromChargeCollector_perSimHit(simHit.layer(), path.travel, path.lengthG4, simHit.truthPos(), trafoMatrix, simHit.CreatedInGenerator()); // fill histograms once per sim hit, with info from the path (eg. travel vector, which contains info on the angle of incidence) + + DistributeSegmentCharge(hitMap, seg, segmentCharge, segmentsInBin, simHit); +} + +Index_segment ChargeCollector_LUT::ComputeSegmentIndices(const int step, const int stepCount, const Path& path) const { + if (step < 0 || step >= stepCount) + throw std::runtime_error("VTXdigi_tools::ChargeCollector_LUT::ComputeSegmentIndices: step (=" + std::to_string(step) + ") out of range [0, " + std::to_string(stepCount-1) + "].");\ + + Index_segment seg; + float t = (step + 0.5f) / stepCount; // +0.5 to get center of segment + const dd4hep::rec::Vector3D pos = path.entry + t * path.travel; + + seg.i = ComputePixelIndices(pos, m_digitizer.PixelPitch(), m_digitizer.PixelCount()); + + seg.j = ComputeInPixelIndices(pos, m_LUT.GetBinCount(), m_digitizer.PixelPitch(), m_digitizer.ActiveVolumeDimensions()); + + return seg; +} + +void ChargeCollector_LUT::DistributeSegmentCharge(HitMap& hitMap, const Index_segment& i_seg, const float segmentCharge, const int segmentsInBin, const SimHitWrapper& simHit) const { + + /* cache things, this is the hottest loop */ + const int lutSize = m_LUT.GetSize(); + const int i_u_origin = i_seg.i.first - m_LUT.GetSizeHalf(); // pix index of leftmost pixel in LUT matrix + const int i_v_origin = i_seg.i.second - m_LUT.GetSizeHalf(); + const int pixelCount_u = static_cast(m_digitizer.PixelCount().first); + const int pixelCount_v = static_cast(m_digitizer.PixelCount().second); + const float charge = segmentCharge * segmentsInBin; + + const int col_min = std::max(0, -i_u_origin); + const int col_max = std::min(lutSize, pixelCount_u - i_u_origin); + + const int row_min = std::max(0, -i_v_origin); + const int row_max = std::min(lutSize, pixelCount_v - i_v_origin); + + for (int col = col_min; col < col_max; ++col) { + const int i_u = i_u_origin + col; // convert from col in [0, matrixSize) to pixel offset in [-matrixSize_half, matrixSize_half]. Note size_half = (size-1)/2 + + for (int row = row_min; row < row_max; ++row) { + const int i_v = i_v_origin + row; + + const float chargeToAdd = m_LUT.GetWeight(i_seg.j, col, row) * charge; + hitMap.FillCharge({i_u, i_v}, chargeToAdd, simHit); + } + } +} + +void ChargeCollector_LUT::MoveTruthPosition(const SimHitWrapper& simHit, const Path& path) const { + if (m_chargeCollectionDepthCenter == 0.f) + return; + + /* shift the sim hit position along the path to the depth that is closest to the target w (ie. target depth) */ + + float t = (m_chargeCollectionDepthCenter - path.entry.z()) / path.travel.z(); // how far along the path do we need to go to get to the target depth? + // t in [0,1] means it's within the path, otherwise it's outside of the path and we will shift to the closest end (entry or exit) + + if (t < 0.f) + t = 0.f; + else if (t > 1.f) + t = 1.f; + + simHit.SetTruthPos(path.entry + t * path.travel); +} + +/* -- Single pixel approach -- */ + +ChargeCollector_SinglePixel::ChargeCollector_SinglePixel(const VTXdigi_Modular& digitizer) : IChargeCollector(digitizer) { + m_digitizer.debug() << "ChargeCollector_SinglePixel constructed." << endmsg; +} + +void ChargeCollector_SinglePixel::FillHit(const SimHitWrapper& simHit, HitMap& hitMap, const TGeoHMatrix& trafoMatrix) const { + (void) trafoMatrix; // Not used in this implementation of ChargeCollector, but we need to keep it as argument to conform to the interface. Silences the unused parameter warning. + + const std::pair i_uv = ComputePixelIndices(simHit.truthPos(), m_digitizer.PixelPitch(), m_digitizer.PixelCount()); + + if (!(i_uv.first == -1 || i_uv.second == -1 || i_uv.first >= static_cast(m_digitizer.PixelCount().first) || i_uv.second >= static_cast(m_digitizer.PixelCount().second))) + hitMap.FillCharge(i_uv, simHit.charge(), simHit); +} + +/* -- Debug approach -- */ + +ChargeCollector_Debug::ChargeCollector_Debug(const VTXdigi_Modular& digitizer) : IChargeCollector(digitizer) { + m_digitizer.debug() << "ChargeCollector_Debug constructed." << endmsg; +} + +void ChargeCollector_Debug::FillHit(const SimHitWrapper& simHit, HitMap& hitMap, const TGeoHMatrix& trafoMatrix) const { + (void) trafoMatrix; // Not used in this implementation of ChargeCollector, but we need to keep it as argument to conform to the interface. Silences the unused parameter warning. + + const dd4hep::rec::Vector3D pos_local = simHit.truthPos(); + const float charge = simHit.charge(); + const std::pair i_uv = ComputePixelIndices(pos_local, m_digitizer.PixelPitch(), m_digitizer.PixelCount()); + + m_digitizer.verbose() << " - SimHit at local position (" << pos_local.x() << ", " << pos_local.y() << ", " << pos_local.z() << ")" << endmsg; + m_digitizer.verbose() << " - and pixel indices (" << i_uv.first << ", " << i_uv.second << ")" << endmsg; + if (i_uv.first == -1 || i_uv.second == -1 || i_uv.first >= static_cast(m_digitizer.PixelCount().first) || i_uv.second >= static_cast(m_digitizer.PixelCount().second)) { + m_digitizer.warning() << "simHit local position (" << pos_local.x() << ", " << pos_local.y() << ", " << pos_local.z() << ") is out of sensor bounds U: [" << -m_digitizer.ActiveVolumeDimensions().at(0)/2 << ", " << m_digitizer.ActiveVolumeDimensions().at(0)/2 << "], V: [" << -m_digitizer.ActiveVolumeDimensions().at(1)/2 << ", " << m_digitizer.ActiveVolumeDimensions().at(1)/2 << "]. This simHit will be skipped." << endmsg; + } + else { + m_digitizer.verbose() << " - Filling charge " << simHit.charge() << " e." << endmsg; + + hitMap.FillCharge(i_uv, 0.5*charge, simHit); + if (i_uv.first + 1 < static_cast(m_digitizer.PixelCount().first)) + hitMap.FillCharge({i_uv.first + 1, i_uv.second}, 0.3*charge, simHit); + if (i_uv.second + 1 < static_cast(m_digitizer.PixelCount().second)) + hitMap.FillCharge({i_uv.first, i_uv.second + 1}, 0.1*charge, simHit); + if (i_uv.second + 2 < static_cast(m_digitizer.PixelCount().second)) + hitMap.FillCharge({i_uv.first, i_uv.second + 2}, 0.1*charge, simHit); + + m_digitizer.verbose() << " - Total charge collected in hitMap: " << hitMap.GetTotalCharge() << " e." << endmsg; + } +} + +// /* -- Drift approach -- */ + +} // namespace VTXdigi_tools + + diff --git a/VTXdigi_Modular/src/ChargeCollector_impl.h b/VTXdigi_Modular/src/ChargeCollector_impl.h new file mode 100644 index 00000000..310cef5d --- /dev/null +++ b/VTXdigi_Modular/src/ChargeCollector_impl.h @@ -0,0 +1,145 @@ +// VTXdigi_Modular/src/ChargeCollector_impl.h +#pragma once + +#include "../include/VTXdigi_Modular.h" + +namespace VTXdigi_tools { +class SimHitWrapper; // forward-declaration for include/VTXdigi_tools.h +class HitMap; // forward-declaration for include/VTXdigi_tools.h + +using Index_pix = std::pair; +using Index_inPix = std::array; + +constexpr float kPathLengthTolerance = 1.05f; // tolerance factor for how much longer the computed path can be compared to the Geant4 path length. If the computed path is longer than the Geant4 path, either the linear path approximation breaks down, or the particle begins or ends inside the sensor volume +constexpr float kLutEntryMinimum = 1.e-3f; // LUT entries below this value are set to zero, to minimise unnecessary computations in hot loop. result is quite sensitive to this, so choose carefully. 1e-5 seems to be a good compromise between accuracy and performance for the TPSCo 65nm CIS LUT + +/** @brief holds pixel indices i and in-pixel bin indices j (eg. for a segment) */ +struct Index_segment { + Index_pix i; + Index_inPix j; + + inline bool operator==(const Index_segment& other) const { + return (i == other.i) && (j == other.j); + } +}; + +/** @brief Computes & then holds position & information about a simHits path through the sensor */ +struct Path { + Path() = default; + + dd4hep::rec::Vector3D entry; // in local coordinates, in mm + dd4hep::rec::Vector3D travel; // in local coordinates, in mm + dd4hep::rec::Vector3D simPos; // in local coordinates, in mm + + float length; // in mm + float lengthG4; + int nSegments; + float charge; +}; + +/** @brief Compute the factors by which to clip a path along a given axis (to clip it to the sensor volume) */ +std::pair ComputePathClippingFactors(std::pair t, const float entry_ax, const float travel_ax, const float sensorLength_ax); + +/** @brief Construct path information from a simHit and the sensor's transformation matrix + * @note returns true if path is valid, false otherwise. false means the path would not intersect the sensor volume */ +bool ConstructPath(Path& path, const SimHitWrapper& simHit, const TGeoHMatrix& trafoMatrix, const VTXdigi_Modular& digitizer); + + +/* -- Charge collector algorithm: LUT-based -- */ + +class LookupTable { + /* TODO: use sparse storage (instead of storing ALL lut entries, even though ~80% are empty (for TPSCo 65nm CIS)). Make sure to keep each matrix contiguous, because we iterate over those. (maybe even move to iterator over matrix instead of getting each weight individually). This will also improve performance in HitMap::FillCharge() because ~80% less operations have to be performed. */ + + Index_inPix m_binCount; + int m_matrixSize; + int m_matrixSize_half; // =(matrixSize-1)/2 (for better performance in hot loop) + std::vector m_matrices; // charge sharing matrices, one per in-pixel bin (flattened for better performance) + float m_chargeCollectionDepthCenter; + +public: + + /** @brief Construct lookup table from a file + * @note Checks that parameters from the LUT file match those in the digitiser */ + LookupTable(const std::string& lutFileName, const VTXdigi_Modular& digitizer); // load weights from file + LookupTable() = default; + + float GetChargeCollectionDepthCenter() const { return m_chargeCollectionDepthCenter; } + + /** @brief Set the charge sharing matrix for a specific in-pixel bin + * @param weights A flat vector containing the matrix weights in row-major order (ie. row-by-row) (length must be matrixSize*matrixSize) */ + void SetMatrix(const Index_inPix& j_uvw, const std::vector& weights); + + /** @brief Set all charge sharing matrices to the same values */ + void SetAllMatrices(const std::vector& weights); + + /** @brief Access matrix as const reference */ + const std::vector& GetMatrix(const Index_inPix& j_uvw) const; + + /** @brief Access a specific weight in a charge sharing matrix */ + inline float GetWeight(const Index_inPix& j_uvw, const int col, const int row) const { + return m_matrices[FindIndex(j_uvw, col, row)]; + }; + + inline int GetSize() const { return m_matrixSize; } + inline int GetSizeHalf() const { return m_matrixSize_half; } + + inline int GetBinCount(int i) const { return m_binCount.at(i); } + inline Index_inPix GetBinCount() const { return m_binCount; } +private: + + /** @brief Convert 3D in-pixel bin indices and a matrix row/column to a flat index for m_matrices */ + int FindIndex (const Index_inPix& j, const int col, const int row) const; +}; // class LookupTable + +class ChargeCollector_LUT : public IChargeCollector { + + LookupTable m_LUT; + const float m_stepLength; // in mm + const bool m_shiftTruthPos; // if true, the truth position in the simHitWrapper is shifted to the depth in the sensor where most charge is collected, to get useful residual plots. Mirrors VTXdigi_Modular::m_LUT_shiftTruthPosition Gaudi property. + +public: + + explicit ChargeCollector_LUT(const VTXdigi_Modular& digitizer); + void FillHit(const SimHitWrapper& simHit, HitMap& hitMap, const TGeoHMatrix& trafoMatrix) const override; + +private: + + /** @brief Compute the pixel- and in-pixel bin indices for a given segment in a path */ + Index_segment ComputeSegmentIndices(const int step, const int stepCount, const Path& path) const; + + /** @brief Fill charge into the hitmap, according to a in-pixel bin's charge sharing weights */ + void DistributeSegmentCharge(HitMap& hitMap, const Index_segment& i_seg, const float charge, const int segmentsInBin, const SimHitWrapper& simHit) const; + + /** @brief Move the truth position in the simHitWrapper along the computed path to the depleted region center (as defined by the digitizer Gaudi property) + * @note This is necessary to get meaningful residual plots for sensors where charge is not collected across the whole sensor thickness. In TPSCo 65nm CIS the depletion region only extends ~10 um down into the sensor, so charges are only collected close to the upper sensor surface. For tracks passing through the sensor, the simHit position is defined at the sensor center. So for angled tracks, the residual (which is truth - reco) would be dominated by the track angle instead of the charge collection effects we want to study. By moving the truth position to the depleted region center, we get residuals that are dominated by charge collection effects. + * @note Only used for plotting the residuals, does not change the output collections in any way. + * @note If the path does not reach the depleted region center, the truth position is moved to the path end closest to the depleted region center */ + void MoveTruthPosition(const SimHitWrapper& simHit, const Path& path) const; +}; + + +/* -- Charge collector algorithms for debugging -- */ + +class ChargeCollector_SinglePixel : public IChargeCollector { +public: + explicit ChargeCollector_SinglePixel(const VTXdigi_Modular& digitizer); + void FillHit(const SimHitWrapper& simHit, HitMap& hitMap, const TGeoHMatrix& trafoMatrix) const override; +}; + +class ChargeCollector_Debug : public IChargeCollector { +public: + explicit ChargeCollector_Debug(const VTXdigi_Modular& digitizer); + void FillHit(const SimHitWrapper& simHit, HitMap& hitMap, const TGeoHMatrix& trafoMatrix) const override; +}; + + + +/* -- Charge collector algorithm: Propagation-based -- */ + +class ChargeCollector_Drift : public IChargeCollector { +public: + explicit ChargeCollector_Drift(const VTXdigi_Modular& digitizer); + void FillHit(const SimHitWrapper& simHit, HitMap& hitMap, const TGeoHMatrix& trafoMatrix) const override; +}; + +} // namespace VTXdigi_tools \ No newline at end of file diff --git a/VTXdigi_Modular/src/VTXdigi_Modular.cpp b/VTXdigi_Modular/src/VTXdigi_Modular.cpp new file mode 100644 index 00000000..da7fcc1e --- /dev/null +++ b/VTXdigi_Modular/src/VTXdigi_Modular.cpp @@ -0,0 +1,1751 @@ +// VTXdigi_Modular/src/VTXdigi_Modular.cpp +#include "VTXdigi_Modular.h" + +DECLARE_COMPONENT(VTXdigi_Modular) + +VTXdigi_Modular::VTXdigi_Modular(const std::string& name, ISvcLocator* svcLoc) + : MultiTransformer(name, svcLoc, + {KeyValues("SimTrackHitCollectionName", {"UNDEFINED_SimTrackHitCollectionName"}), + KeyValues("HeaderName", {"UNDEFINED_HeaderName"}),}, + {KeyValues("TrackerHitCollectionName", {"UNDEFINED_TrackerHitCollectionName"}), + KeyValues("SimTrkHitRelationsCollection", {"UNDEFINED_SimTrkHitRelationsCollection"})}) { + info() << "Constructed successfully" << endmsg; +} + +StatusCode VTXdigi_Modular::initialize() { + info() << "INITIALIZING VTXdigi_Modular..." << endmsg; + + info() << "OutputLevel set to " << msgSvc()->outputLevel(name()) << endmsg; + // TODO: implement if-clause for debug/verbose messages in hot loops (to avoid constructing the message string when the message won't be printed) -> this will improve performance significantly + + // if (!VTXdigi_tools::ToolTest()) + // return StatusCode::FAILURE; + + InitServicesAndGeometry(); + + InitLayersAndSensors(); + + if (m_debugHistograms) + InitHistograms(); + + // This needs to come in after the properties, geometry and services have all been initialized + verbose() << "Initializing charge collection method: " << m_chargeCollectionMethod.value() << endmsg; + m_chargeCollector = VTXdigi_tools::CreateChargeCollector(*this, m_chargeCollectionMethod); + + info() << " - Initialized successfully." << endmsg; + return StatusCode::SUCCESS; +} + +StatusCode VTXdigi_Modular::finalize() { + info() << "FINALIZING VTXdigi_Modular..." << endmsg; + + PrintCountersSummary(); + + debug() << " - finalized successfully." << endmsg; + return StatusCode::SUCCESS; +} + + +/* ---- Event loop ---- */ + +std::tuple VTXdigi_Modular::operator() + (const edm4hep::SimTrackerHitCollection& simTrackerHits, const edm4hep::EventHeaderCollection& headers) const { + if (!CheckEventSetup(simTrackerHits, headers)) { + return std::make_tuple(edm4hep::TrackerHitPlaneCollection(), edm4hep::TrackerHitSimTrackerHitLinkCollection()); + } + + /* TODO: Implement fast digitization, where simHit pos is simply smeared by a Gaussian. + * this makes a lot of the init etc unnecessary, so maybe keep it in a separate class? We cannot simply add a ChargeCollector implementation to do this, because ChargeCollectors act on pixels, but fast digitization acts on the position itself (ofc we could emulate this in clusters, but thats incredibly inefficient and stinks) */ + + std::unordered_map> sensorSimHits; // map from sensor (shortened cellID) to hits + for (const edm4hep::SimTrackerHit& simTrackerHit : simTrackerHits) { + + + if (CheckSimhitLayer(simTrackerHit)){ + const dd4hep::DDSegmentation::CellID cellID = VTXdigi_tools::GetCellID_short(simTrackerHit); + sensorSimHits[cellID].emplace_back(simTrackerHit, m_cellIdDecoder); // simTrackerHits are copied here. Pointers to these are passed around (eg. in hit/pixel/cluster objects). + if (VTXdigi_tools::CreatedInGenerator(simTrackerHit)) { + ++m_counter_accSimHitsCreatedInGenerator; + } + if (!simTrackerHit.isProducedBySecondary()) { // this is the same as SimHitWrapper::hasDirectMcParticleLink() + ++m_counter_accSimHitsDirectMcParticleLink; + } + } + } + debug() << " - Found " << simTrackerHits.size() << " simTrackerHits on " << sensorSimHits.size() << " individual sensors. Digitising sensor-wise..." << endmsg; + + auto digiHits = edm4hep::TrackerHitPlaneCollection(); + auto digiHitLinks = edm4hep::TrackerHitSimTrackerHitLinkCollection(); + std::vector hitsSensor; + + /* loop over sensors */ + for (const auto& [cellID, simHits] : sensorSimHits) { + debug() << " - Processing sensor with cellID " << cellID << " (layer " << simHits.back().layer() << "). Has " << simHits.size() << " simHits." << endmsg; + + const TGeoHMatrix trafoMatrix = VTXdigi_tools::ComputeSensorTrafoMatrix(cellID, m_volumeManager, m_sensorNormalRotation); // transformation from global detector to local sensor frame + + /* Fill hit map with charge from all simHits on this sensor, using the selected charge collection method */ + VTXdigi_tools::HitMap hitMap(m_pixelCount); + + for (const VTXdigi_tools::SimHitWrapper& simHit : simHits) { + const dd4hep::rec::Vector3D pos_global = VTXdigi_tools::ConvertVector(simHit.hitPtr()->getPosition()); + simHit.SetTruthPos(VTXdigi_tools::GlobalToLocal(pos_global, trafoMatrix)); // do this only now to not compute the trafo matrix twice. TruthPos might be shifted by the charge collection algorithm later. + debug() << " - Processing simHit, charge dep. " << simHit.charge() << " e at (" << simHit.truthPos().x() << ", " << simHit.truthPos().y() << ", " << simHit.truthPos().z() << ") local, (" << pos_global.x() << ", " << pos_global.y() << ", " << pos_global.z() << ") global" << endmsg; + + m_chargeCollector->FillHit(simHit, hitMap, trafoMatrix); // uses the selected charge collection method + + if (m_debugHistograms.value()) + FillHistograms_perSimHit(simHit); + } + + + if (m_smearing_charge.value() > 0.f) + hitMap.ApplyChargeSmearing(m_rndm_charge); + if (m_threshold.value() > 0.f) + hitMap.ApplyThreshold(m_threshold.value()); + + std::vector clusters = Clusterize(hitMap); + + CreateDigiHits(digiHits, digiHitLinks, cellID, trafoMatrix, clusters); + if (m_debugHistograms.value()) + FillHistograms_perSensor(simHits, digiHits, trafoMatrix, cellID); + } /* loop over sensors */ + + debug() << " - Finished digitization. Created " << digiHits.size() << " digiHits from " << simTrackerHits.size() << " simTrackerHits." << endmsg; + return std::make_tuple(std::move(digiHits), std::move(digiHitLinks)); +} // operator() + +/* ---- Initialization & finalization functions ---- */ + +void VTXdigi_Modular::InitServicesAndGeometry() { + /* Sets the members: + * - m_uidService + * - m_geoService + * - m_cellIdDecoder + * - m_detector + * - m_surfaceMap + * - m_volumeManager + * - m_subDetector + */ + if (m_threshold.value() < 0.f) + throw GaudiException("Threshold " + std::to_string(m_threshold.value()) + " e- is negative.", "VTXdigi_Modular::InitServicesAndGeometry()", StatusCode::FAILURE); + if (m_smearing_charge.value() < 0.f) + throw GaudiException("Charge smearing sigma " + std::to_string(m_smearing_charge.value()) + " e- is negative.", "VTXdigi_Modular::InitServicesAndGeometry()", StatusCode::FAILURE); + + if (m_threshold.value() <= 5 * m_smearing_charge.value()) + warning() << "Threshold " << m_threshold.value() << " e- is less than 5 times the charge smearing sigma " << m_smearing_charge.value() << " e-. This digitiser does only apply smearing to pixels that collect any charge from simHits, so it cannot simulate random firing pixels. (doing this by drawing a noise for every pixel in the detector for every event would be INCREDIBLY slow. A work-around to simulate random pixels firing might be implemented in the future)." << endmsg; + if (m_smearing_time.value() < 0.f) + throw GaudiException("Time smearing sigma " + std::to_string(m_smearing_time.value()) + " ns is negative.", "VTXdigi_Modular::InitServicesAndGeometry()", StatusCode::FAILURE); + + if (m_positionUncertainty.value().empty()) + info() << "Cluster position uncertainty not set. Doing simple charge-weighted uncertainty estimation." << endmsg; + else if (m_positionUncertainty.value().size() == 2) + info() << "Cluster position uncertainty set to (" << m_positionUncertainty.value().at(0) << " mm, " << m_positionUncertainty.value().at(1) << " mm) in u and v direction." << endmsg; + else if (m_positionUncertainty.value().size() == 6) + info() << "Cluster position uncertainty set to (" << m_positionUncertainty.value().at(0) << ", " << m_positionUncertainty.value().at(1) << ", " << m_positionUncertainty.value().at(2) << ") mm and (" << m_positionUncertainty.value().at(3) << ", " << m_positionUncertainty.value().at(4) << ", " << m_positionUncertainty.value().at(5) << ") mm for cluster lengths of (1, 2, 3), in u and v direction, respectively." << endmsg; + else + throw GaudiException("Property ClusterPositionUncertainty must be either empty (for charge-weighted estimation), have exactly 2 values (for fixed uncertainty in u and v), or six values (for cluster-length based estimation).", "VTXdigi_Modular::InitServicesAndGeometry()", StatusCode::FAILURE); + + m_randomService = service("RndmGenSvc", false); + if (!m_randomService) + throw GaudiException("Unable to get RndmGenSvc.", "VTXdigi_Modular::InitServicesAndGeometry()", StatusCode::FAILURE); + + if (m_rndm_charge.initialize(m_randomService, Rndm::Gauss(0., m_smearing_charge.value())).isFailure()) + throw GaudiException("Unable to initialize random number generator for charge smearing.", "VTXdigi_Modular::InitServicesAndGeometry()", StatusCode::FAILURE); + if (m_rndm_time.initialize(m_randomService, Rndm::Gauss(0., m_smearing_time.value())).isFailure()) + throw GaudiException("Unable to initialize random number generator for time smearing.", "VTXdigi_Modular::InitServicesAndGeometry()", StatusCode::FAILURE); + + if (m_subDetName.value() == m_undefinedString) + throw GaudiException("Property SubDetectorName is not set!", "VTXdigi_Modular::InitServicesAndGeometry()", StatusCode::FAILURE); + + m_geoService = serviceLocator()->service(m_geoServiceName); + if (!m_geoService) + throw GaudiException("Unable to retrieve the GeoSvc", "VTXdigi_Modular::InitServicesAndGeometry()", StatusCode::FAILURE); + + /* TODO: get cellID without using geoService, a la Jessy (if this is advantageous) */ + std::string cellIDstr = m_geoService->constantAsString(m_encodingStringVariable.value()); + m_cellIdDecoder = std::make_unique(cellIDstr); + if (!m_cellIdDecoder) + throw GaudiException("Unable to retrieve the cellID decoder", "VTXdigi_Modular::InitServicesAndGeometry()", StatusCode::FAILURE); + + m_detector = m_geoService->getDetector(); + if (!m_detector) + throw GaudiException("Unable to retrieve the DD4hep detector from GeoSvc", "VTXdigi_Modular::InitServicesAndGeometry()", StatusCode::FAILURE); + + const dd4hep::rec::SurfaceManager* surfaceManager = m_detector->extension(); + if (!surfaceManager) + throw GaudiException("Unable to retrieve the SurfaceManager from the DD4hep detector", "VTXdigi_Modular::InitServicesAndGeometry()", StatusCode::FAILURE); + + m_surfaceMap = surfaceManager->map(m_subDetName.value()); + if (!m_surfaceMap) + throw GaudiException("Unable to retrieve the simSurface map for subdetector " + m_subDetName.value(), "VTXdigi_Modular::InitServicesAndGeometry()", StatusCode::FAILURE); + + m_volumeManager = m_detector->volumeManager(); + if (!m_volumeManager.isValid()) + throw GaudiException("Unable to retrieve the VolumeManager from the DD4hep detector", "VTXdigi_Modular::InitServicesAndGeometry()", StatusCode::FAILURE); + + + { /* DD4hep has a transformation from the global detector coordinates to each sensors local system. The definition of the local system might change. + * We have to make sure that a sensor always sits in the u-v plane in the local system, eg. we might have to swap the axes of the local system. + * We accomplish this with a transformation matrix. This is valid for all sensors in the subdetector */ + + /* Set the sensor local transformation matrix, st. the sensors normal vector is parallel to w-axis (by simply looking at the first sensor in the map) */ + auto surfaceMapIter = m_surfaceMap->begin(); + if (surfaceMapIter == m_surfaceMap->end()) + throw GaudiException("Surface map for subdetector " + m_subDetName.value() + " is empty.", "VTXdigi_Modular::InitServicesAndGeometry()", StatusCode::FAILURE); + const unsigned long cellID = surfaceMapIter->first; + const dd4hep::rec::ISurface* surface = surfaceMapIter->second; + + TGeoHMatrix sensorTrafoMatrix = m_volumeManager.lookupDetElement(cellID).nominal().worldTransformation(); + double tempVec[3]; + + const double epsilon = 1.0e-6; // reasonable for comparing to 1 + + /* first: rotate sensor normal onto W-axis*/ + sensorTrafoMatrix.MasterToLocalVect(surface->normal().unit(), tempVec); + dd4hep::rec::Vector3D n_local(tempVec[0], tempVec[1], tempVec[2]); + + if (std::abs(n_local.x() - 1.0) < epsilon) { + debug() << " - Local sensor normal vector is (1,0,0). Defining rotation matrix to rotate it to (0,0,1)." << endmsg; + m_sensorNormalRotation = TGeoRotation("rot",90.,90.,0.); + } + else if (std::abs(n_local.y() - 1.0) < epsilon) { + debug() << " - Local sensor normal vector is (0,1,0). Defining rotation matrix to rotate it to (0,0,1)." << endmsg; + m_sensorNormalRotation = TGeoRotation("rot",0.,-90.,0.); + } + else if (std::abs(n_local.z() - 1.0) < epsilon) { + debug() << " - Local sensor normal vector is already parallel to (0,0,1)." << endmsg; + // no rotation needed + } + else { + error() << "Sensor local normal vector is not aligned with any local axis!" << endmsg; + } + + /* TODO: this only makes sure that the sensor normal vector is perpendicular to the surface. + * It does NOT make sure that the local x and y are not swapped and have the correct polarity. + * I tried coming up with the correct transformation matrix based on the axis rotations, but tbh I am quite stumped by how the rotation acts on the vectors, and how I can fix it. */ + TGeoHMatrix M = VTXdigi_tools::ComputeSensorTrafoMatrix(cellID, m_volumeManager, m_sensorNormalRotation); + M.MasterToLocalVect(surface->u().unit(), tempVec); + if (std::abs(tempVec[0]-1.) > epsilon || std::abs(tempVec[1]) > epsilon || std::abs(tempVec[2]) > epsilon) + warning() << "After rotation, local sensor U direction (" << tempVec[0] << "," << tempVec[1] << "," << tempVec[2] << ") is not parallel to (1,0,0) axis!" << endmsg; + M.MasterToLocalVect(surface->v().unit(), tempVec); + if (std::abs(tempVec[0]) > epsilon || std::abs(tempVec[1]-1.) > epsilon || std::abs(tempVec[2]) > epsilon) + warning() << "After rotation, local sensor V direction (" << tempVec[0] << "," << tempVec[1] << "," << tempVec[2] << ") is not parallel to (0,1,0) axis!" << endmsg; + M.MasterToLocalVect(surface->normal().unit(), tempVec); + if (std::abs(tempVec[0]) > epsilon || std::abs(tempVec[1]) > epsilon || std::abs(tempVec[2]-1.) > epsilon) + warning() << "After rotation, local sensor normal direction (" << tempVec[0] << "," << tempVec[1] << "," << tempVec[2] << ") is not parallel to (0,0,1) axis!" << endmsg; + } + + /* IDEA / Allegro: + * - subDet is child of detector, eg. Vertex + * - subDetChild is child of subDet, eg. VertexBarrel + * - subDetChildChild are layers + * If this is not the case, layers might be direct children of subDet: */ + + debug() << " - Retrieving subdetector " << m_subDetName.value() << " from DD4hep detector..." << endmsg; + verbose() << " - The detector has the following subDetectors: " << endmsg; + for (const auto& [subDetName, subDet] : m_detector->detectors()) { + verbose() << " - " << subDetName << endmsg; + } + + if (m_subDetChildName.value() != m_undefinedString) { + /* IDEA/Allegro setup */ + const dd4hep::DetElement subDet = m_detector->detector(m_subDetName.value()); + m_subDetector = subDet.child(m_subDetChildName.value()); + } + else { + /* layers are direct children of subDet */ + m_subDetector = m_detector->detector(m_subDetName.value()); + } + if (!m_subDetector) + throw GaudiException("Unable to retrieve the subdetector DetElement " + m_subDetName.value(), "VTXdigi_Modular::InitServicesAndGeometry()", StatusCode::FAILURE); + + debug() << " - Retrieved all necessary services and dd4hep detector elements." << endmsg; +} + +void VTXdigi_Modular::InitLayersAndSensors() { + /* Get the detector & sensor geometry information from the DD4hep detector + * + * Sets the members: + * - m_pixelPitch + * - m_sensorActiveThickness + * - m_layers + * - */ + verbose() << " - Retrieving layers, subdetector geometry, and sensor size and pitch..." << endmsg; + + /* Find layers that we want to digitize */ + verbose() << " - Determining relevant layers... " << endmsg; + std::vector availableLayers; + for (const auto& [layerName, layerObj] : m_subDetector.children()) { + dd4hep::VolumeID layerVolumeID = layerObj.volumeID(); + const int layer = m_cellIdDecoder->get(layerVolumeID, "layer"); + availableLayers.push_back(layer); + } + if (availableLayers.empty()) + throw GaudiException("No layers found in subdetector " + m_subDetName.value(), "VTXdigi_Modular::InitLayersAndSensors()", StatusCode::FAILURE); + + if (m_layers.value().empty()) { + /* digitize all layers in m_subDetector */ + m_layers.value() = availableLayers; + } else { + /* If layers-to-digitize is specified: check that all requested layers exist */ + + for (const auto layer : m_layers.value()) { + if (std::find(availableLayers.begin(), availableLayers.end(), layer) == availableLayers.end()) { + throw GaudiException("Requested layer " + std::to_string(layer) + " to digitize, but this layer does not exist in subdetector " + m_subDetName.value(), "VTXdigi_Modular::InitLayersAndSensors()", StatusCode::FAILURE); + } + } + debug() << " - Digitizing only specific layers, as defined in Gaudi property by the user. All requested layers were found." << endmsg; + } + info() << " - Digitizing " << m_layers.value().size() << " layers: " << m_layers.value() << " in subdetector " << m_subDetName.value() << "." << endmsg; + + /* Get pixel pitch from the segmentation (from the readout that matches our simHitCollection) + * I don't know of a better way to do this (that also works for the IDEA detector model) */ + std::string simHitCollectionName; + if (this->getProperty("SimTrackHitCollectionName", simHitCollectionName).isFailure()) + throw GaudiException("Could not retrieve SimTrackHitCollectionName property while checking geometry consistency.", "VTXdigi_Modular::InitLayersAndSensors()", StatusCode::FAILURE); + + dd4hep::Detector::HandleMap readoutHandleMap = m_detector->readouts(); + int readoutCount = 0; + std::string matchedReadoutKey; + /* loop over readouts, see if one matches our simHitCollection */ + for (const auto& [readoutKey, readoutHandle] : readoutHandleMap) { + if (simHitCollectionName.find(readoutKey) != std::string::npos) { + ++readoutCount; + if (readoutCount != 1) + continue; + matchedReadoutKey = readoutKey; + + const dd4hep::Segmentation& segmentation = m_detector->readout(readoutKey).segmentation(); + if (!segmentation.isValid()) + throw GaudiException("Segmentation for readout " + readoutKey + " is not valid while checking geometry consistency.", "VTXdigi_Modular::InitLayersAndSensors()", StatusCode::FAILURE); + const auto cellDimensions = segmentation.cellDimensions(0); // this assumes all cells have the same dimensions (ie. only one sensor type in this readout) + m_pixelPitch.first = cellDimensions.at(0) * 10; // convert cm to mm + m_pixelPitch.second = cellDimensions.at(1) * 10; + /* TODO: get pixel count from segmentation*/ + } + else { + verbose() << " - Readout \"" << readoutKey << "\" does NOT MATCH SimTrackHitCollectionName \"" << simHitCollectionName << "\". skipping." << endmsg; + } + } // end loop over readouts + if (readoutCount == 0) + throw GaudiException("Could not find any readout matching SimTrackHitCollectionName " + simHitCollectionName + " in detector while checking geometry consistency.", "VTXdigi_Modular::InitLayersAndSensors()", StatusCode::FAILURE); + else if (readoutCount > 1) + warning() << "Found multiple (" << readoutCount << ") readouts matching SimTrackHitCollectionName \"" << simHitCollectionName << "\" in detector while checking geometry consistency. Used the first one found. Enable verbose messages for more info." << endmsg; + + + /* TODO: Check that local sensor coordinates (u,v,w) are correctly defined wrt. to the global coordinates (already done in VTXdigi_Allpix2 master branch, but VERY clunky). This requires deep understanding of coordinate systems. I think there is a easy way to do this. I have not figured it out yet */ + + + /* Loop over all sensors in the subDetector and check that they have the same dimensions + * This takes << 1s for the IDEA VTX */ + int moduleNumber = 0, sensorNumber = 0; + bool membersDefined = false; + for (const auto& [layerKey, layerObj] : m_subDetector.children()) { + dd4hep::VolumeID layerVolumeID = layerObj.volumeID(); + int layer = m_cellIdDecoder->get(layerVolumeID, "layer"); + /* If list of layers is given: skip layers that are not on the list*/ + if (!m_layers.value().empty()) { + if (std::find(m_layers.value().begin(), m_layers.value().end(), layer) == m_layers.value().end()) { + verbose() << " - Skipping layer " << layerKey << " (layer " << layer << ", volumeID " << layerVolumeID << " with " << layerObj.children().size() << " modules) as it is not in the LayersToDigitize list." << endmsg; + continue; + } + } + verbose() << " - Found layer \"" << layerKey << "\" (layer " << layer << ", volumeID " << layerVolumeID << ", " << layerObj.children().size() << " modules) for subDetector \"" << m_subDetName.value() << "\"." << endmsg; + + for (const auto& [moduleKey, moduleObj] : layerObj.children()) { + ++moduleNumber; + for (const auto& [sensorKey, sensorObj] : moduleObj.children()) { + ++sensorNumber; + + // FIRST: find thickness of the active sensor volume + // using the active volume solid of this sensor + dd4hep::VolumeID sensorVolumeID = sensorObj.volumeID(); + dd4hep::Volume sensorVolume = sensorObj.volume(); + std::array solidDimensions{}; // dimensions of active volume. indices must not correspond to local sensor (u,v,w) axes, but might be swapped + + // sensorVolume.solid() returns a dd4hep::Solid_type object, generalised as dd4hep::Solid. This can either be a box or a trapezoid (for sensors in IDEA / ALLEGRO) + try { + dd4hep::Box sensorBox = sensorVolume.solid(); // directions 0,1,2 do not necessarily correspond to the local sensor u,v,w directions, but might be swapped. + solidDimensions[0] = sensorBox.x() * 2 * 10; + solidDimensions[1] = sensorBox.y() * 2 * 10; + solidDimensions[2] = sensorBox.z() * 2 * 10; + + verbose() << " - Sensor \"" << sensorKey << "\" (layer " << layerKey << ", module " << moduleKey << ", volumeID " << sensorVolumeID << ") has a dd4hep::Box solid." << endmsg; + } + catch (...) { + try { + dd4hep::Trd1 sensorTrd1 = sensorVolume.solid(); + solidDimensions[0] = (sensorTrd1.dX1() + sensorTrd1.dX2()) / 2 * 2 * 10; // average of upper and lower base of trapezoid. Convert half-length in cm to full length in mm + solidDimensions[1] = sensorTrd1.dY() * 2 * 10; // Convert half-length in cm to full length in mm + solidDimensions[2] = sensorTrd1.dZ() * 2 * 10; + // Note: there is some weirdness in dX1() and dX2() with Trd1. I did not dig into this. Assume that these might be a bit funky. + verbose() << " - Sensor \"" << sensorKey << "\" (layer " << layerKey << ", module " << moduleKey << ", volumeID " << sensorVolumeID << ") has a dd4hep::Trd1 solid." << endmsg; + } + catch (...) { + throw GaudiException("Unknown sensor solid type found (neither dd4hep::Box nor dd4hep::Trd1).", "VTXdigi_Modular::InitLayersAndSensors()", StatusCode::FAILURE); + } + } + + const uint thicknessIndex = std::distance(solidDimensions.begin(), std::min_element(solidDimensions.begin(), solidDimensions.end())); + const float solidThickness = solidDimensions.at(thicknessIndex); + const float solidLength_0 = solidDimensions.at((thicknessIndex+1) % 3); + const float solidLength_1 = solidDimensions.at((thicknessIndex+2) % 3); + + // SECOND: find lengths of the sensor (these will match the segmentation), and amount of inactive material above and below + // using the sensitive surface of this sensor (which lies in the middle of the active volume) + const auto surfaceIt = m_surfaceMap->find(sensorObj.volumeID()); + if (surfaceIt == m_surfaceMap->end()) { + throw GaudiException("Could not find surface for sensor " + sensorKey + " (volumeID " + std::to_string(sensorVolumeID) + ") in layer " + std::to_string(layer) + " of subDetector " + m_subDetName.value() + " while checking geometry consistency.", "VTXdigi_Modular::InitLayersAndSensors()", StatusCode::FAILURE); + } + dd4hep::rec::ISurface* surface = surfaceIt->second; + if (!surface) { + throw GaudiException("Surface pointer for sensor " + sensorKey + " (volumeID " + std::to_string(sensorVolumeID) + ") in layer " + std::to_string(layer) + " of subDetector " + m_subDetName.value() + " is null while checking geometry consistency.", "VTXdigi_Modular::InitLayersAndSensors()", StatusCode::FAILURE); + } + + const float surfaceLength_u = surface->length_along_u() * 10; // convert cm to mm + const float surfaceLength_v = surface->length_along_v() * 10; + const float surfaceThickness_above = surface->outerThickness() * 10; // sensor thickness measured from w=0 upwards, including inactive material above the active volume. + const float surfaceThickness_below = surface->innerThickness() * 10; // same, but below + //Note: the sensor local coordinate system (u,v,w) is centered on the active volume, so inactive material upper/lower might be assymetric + + // THIRD: consistency checks on sensor dimensions + if (solidThickness > (surfaceThickness_above + surfaceThickness_below)) { + throw GaudiException("Solid sensor thickness " + std::to_string(solidThickness) + " mm is larger than total sensor thickness " + std::to_string(surfaceThickness_above) + " + " + std::to_string(surfaceThickness_below) + " = " + std::to_string(surfaceThickness_above + surfaceThickness_below) + " mm (including inactive material) in sensor " + sensorKey + " (volumeID " + std::to_string(sensorVolumeID) + ") in layer " + std::to_string(layer) + " of subDetector " + m_subDetName.value() + ". This indicates an inconsistency in the geometry description.", "VTXdigi_Modular::InitLayersAndSensors()", StatusCode::FAILURE); + } + + float tolerance = 0.01f; // 10 um should be small enough to find geometry inconsistencies, but does not trip on trapezoid weirdness in IDEA ultra light in O(2 um) + if ( + !(( std::abs(solidLength_0-surfaceLength_u) < tolerance) && (std::abs(solidLength_1-surfaceLength_v) < tolerance)) && + !((std::abs(solidLength_1-surfaceLength_u) < tolerance) && (std::abs(solidLength_0-surfaceLength_v) < tolerance)) + ) + { + throw GaudiException("Sensor solid (active volume) dimensions " + std::to_string(solidLength_0) + "mm and " + std::to_string(solidLength_1) + "mm could not be matched to sensitive surface dimensions (" + std::to_string(surfaceLength_u) + " x " + std::to_string(surfaceLength_v) + ") mm (including inactive material) within a tolerance of " + std::to_string(tolerance) + " mm. (sensor " + sensorKey + ", volumeID " + std::to_string(sensorVolumeID) + ", layer " + std::to_string(layer) + ", subDetector " + m_subDetName.value() + "). This indicates an inconsistency in the geometry description.", "VTXdigi_Modular::InitLayersAndSensors()", StatusCode::FAILURE); + } + + // FOURTH: apply the parameters we found + if (!membersDefined) { + // For the first sensor we find, set the digitizer class members + m_sensorActiveThickness = solidThickness; + m_inactiveMaterialAbove = surfaceThickness_above - solidThickness/2.f; + m_inactiveMaterialBelow = surfaceThickness_below - solidThickness/2.f; + m_sensorLength.first = surfaceLength_u; + m_sensorLength.second = surfaceLength_v; + + float pixelCountU = m_sensorLength.first / m_pixelPitch.first; + float pixelCountV = m_sensorLength.second / m_pixelPitch.second; + if (abs(pixelCountU - std::round(pixelCountU)) > 0.0001 || abs(pixelCountV - std::round(pixelCountV)) > 0.0001) + throw GaudiException("Sensor side length (" + std::to_string(m_sensorLength.first) + " x " + std::to_string(m_sensorLength.second) + ") mm and pixel pitch (" + std::to_string(m_pixelPitch.first) + " x " + std::to_string(m_pixelPitch.second) + ") mm result in a non-integer pixel count (" + std::to_string(pixelCountU) + " x " + std::to_string(pixelCountV) + ") in subDetector " + m_subDetName.value() + ".", "VTXdigi_Modular::InitLayersAndSensors()", StatusCode::FAILURE); + m_pixelCount.first = std::round(pixelCountU); + m_pixelCount.second = std::round(pixelCountV); + membersDefined = true; + debug() << " - Setting expected dimensions according to first found sensor. Key: " << sensorKey << ", volumeID: " << sensorVolumeID << " (sensor " << sensorNumber << " in layer " << layer << ")" << endmsg; + debug() << " - Solid: thickness " << solidThickness << "mm, lengths " << solidLength_0 << " / " << solidLength_1 << "mm (arbitray ordering)" << endmsg; + debug() << " - Sensitive surface: thickness " << surfaceThickness_above << "mm / " << surfaceThickness_below << "mm (above/below 0), lengths " << surfaceLength_u << "mm / " << surfaceLength_v << "mm (u/v)" << endmsg; + } + else { + // For all other sensors: check for consistency with first sensor + if (std::abs(surfaceLength_u - m_sensorLength.first) > 0.001 || std::abs(surfaceLength_v - m_sensorLength.second) > 0.001 || std::abs(solidThickness - m_sensorActiveThickness) > 0.001) + throw GaudiException("Sensor dimension mismatch found in sensor " + sensorKey + " (volumeID " + std::to_string(sensorVolumeID) + ") in layer " + std::to_string(layer) + " of subDetector " + m_subDetName.value() + ": expected dimensions of (" + std::to_string(m_sensorLength.first) + " x " + std::to_string(m_sensorLength.second) + " x " + std::to_string(m_sensorActiveThickness) + ") mm3, but found (" + std::to_string(surfaceLength_u) + " x " + std::to_string(surfaceLength_v) + " x " + std::to_string(solidThickness) + ") mm3. This algorithm expects exactly one type of sensor per subDetector. Use different instances of the algorithm if different layers consist of different sensors.", "VTXdigi_Modular::InitLayersAndSensors()", StatusCode::FAILURE); + } + } // loop over sensors + } // loop over modules + } // loop over layers + + info() << " - Retrieved sensor parameters: area (" << m_sensorLength.first << " x " << m_sensorLength.second << ") mm2, thickness " << m_sensorActiveThickness << " mm with (" << m_inactiveMaterialAbove << "/" << m_inactiveMaterialBelow << ") mm of inactive material above/below, pixel pitch (" << m_pixelPitch.first << " x " << m_pixelPitch.second << ") mm, pixel count (" << m_pixelCount.first << " x " << m_pixelCount.second << "). All " << sensorNumber << " sensors in the relevant layers share these parameters." << endmsg; +} // InitLayersAndSensors() + +void VTXdigi_Modular::InitHistograms() { + /* Define axes globally to make adjusting them easier + * TODO: Make some of these adjustable via Gaudi Parameters? Might not be necessary.*/ + Gaudi::Accumulators::Axis axis_xy{4000, -20, 20}; + Gaudi::Accumulators::Axis axis_z{4000, -200, 200}; + Gaudi::Accumulators::Axis axis_cosTheta{100, 0, 1}; + Gaudi::Accumulators::Axis axis_theta{4*180, 0, 180}; + Gaudi::Accumulators::Axis axis_phi{4*180, -180, 180}; + + Gaudi::Accumulators::Axis axis_MomentumFraction{400, -1.0001f, 1.0001f}; + + Gaudi::Accumulators::Axis axis_moduleID{2000, -0.5f, 1999.5f}; + Gaudi::Accumulators::Axis axis_clusterSize{60, 0.5f, 60.5f}; + Gaudi::Accumulators::Axis axis_E{1000, 0, m_sensorActiveThickness*2000.f}; + Gaudi::Accumulators::Axis axis_charge{1000, 0, m_sensorActiveThickness*500000.f}; + Gaudi::Accumulators::Axis axis_particleE{1000, 0, 10.f}; // GeV + Gaudi::Accumulators::Axis axis_momentum_keV{10000, 0.f, 1000.f}; + Gaudi::Accumulators::Axis axis_momentum_MeV{10000, 0.f, 1000.f}; + Gaudi::Accumulators::Axis axis_momentum_GeV{10000, 0.f, 1000.f}; + Gaudi::Accumulators::Axis axis_time{1000, 0.f, 1000.f}; + + Gaudi::Accumulators::Axis axis_residual{1600, -200.f, 200.f}; + Gaudi::Accumulators::Axis axis_residual_abs{800, 0.f, 200.f}; + Gaudi::Accumulators::Axis axis_residual_pixels{2020, -50.5f, 50.5f}; // 20 in-pix bins + Gaudi::Accumulators::Axis axis_pdg{1401, -700.5f, 700.5f}; + + Gaudi::Accumulators::Axis axis_pixels_u{ + static_cast(m_pixelCount.first), + -0.5f, + static_cast(m_pixelCount.first+0.5)}; + Gaudi::Accumulators::Axis axis_pixels_v{ + static_cast(m_pixelCount.second), + -0.5f, + static_cast(m_pixelCount.second+0.5)}; + + Gaudi::Accumulators::Axis axis_pathLength{500, 0.f, m_sensorActiveThickness*1000.f*10.f}; + Gaudi::Accumulators::Axis axis_pathTravel{1000, -m_sensorActiveThickness*1000.f*10.f, m_sensorActiveThickness*1000.f*10.f}; + + /* Fill histograms per layer */ + for (int layer : m_layers.value()) { + if (layer == 0) { + axis_z = Gaudi::Accumulators::Axis{100, -96.5, 96.5}; // Want to cover layer 0 of IDEA vertex det perfectly to avoid binning-edge-effects, so we use a the correct length of 185 mm + } + + std::array< std::unique_ptr< Gaudi::Accumulators::StaticHistogram< 1, Gaudi::Accumulators::atomicity::full, float > >, hist1dArrayLen > hist1d; + + hist1d.at(hist1d_simHitE).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_energyDep", + "SimHit deposited energy - Layer " + std::to_string(layer) + ";Energy [keV];Entries", + axis_E + } + ); + hist1d.at(hist1d_simHitCharge).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_chargeDep", + "SimHit deposited charge - Layer " + std::to_string(layer) + ";Charge [e-];Entries", + axis_charge + } + ); + hist1d.at(hist1d_simHitMomentum_keV).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_momentum_keV", + "SimHit particle momentum at hit position - Layer " + std::to_string(layer) + ";Momentum [keV/c];Entries", + axis_momentum_keV + } + ); + hist1d.at(hist1d_simHitMomentum_MeV).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_momentum_MeV", + "SimHit particle momentum at hit position - Layer " + std::to_string(layer) + ";Momentum [MeV/c];Entries", + axis_momentum_MeV + } + ); + hist1d.at(hist1d_simHitMomentum_GeV).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_momentum_GeV", + "SimHit particle momentum at hit position - Layer " + std::to_string(layer) + ";Momentum [GeV/c];Entries", + axis_momentum_GeV + } + ); + + hist1d.at(hist1d_simHit_x).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_global_x", + "Global x-position of simHits - Layer " + std::to_string(layer) + ";x [mm];Entries", + axis_xy + } + ); + hist1d.at(hist1d_simHit_y).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_global_y", + "Global y-position of simHits - Layer " + std::to_string(layer) + ";y [mm];Entries", + axis_xy + } + ); + hist1d.at(hist1d_simHit_z).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_global_z", + "Global z-position of simHits - Layer " + std::to_string(layer) + ";z [mm];Entries", + axis_z + } + ); + hist1d.at(hist1d_simHit_z_createdInGenerator).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_global_z_createdInGenerator", + "Global z-position of simHits from particles created in the generator - Layer " + std::to_string(layer) + ";z [mm];Entries", + axis_z + } + ); + hist1d.at(hist1d_simHit_z_createdInSimulation).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_global_z_createdInSimulation", + "Global z-position of simHits from particles created in simulation - Layer " + std::to_string(layer) + ";z [mm];Entries", + axis_z + } + ); + + hist1d.at(hist1d_simHit_Vertex_x).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_Vertex_x", + "X-position of the production vertex of simHits - Layer " + std::to_string(layer) + ";x [mm];Entries", + axis_xy + } + ); + hist1d.at(hist1d_simHit_Vertex_y).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_Vertex_y", + "Y-position of the production vertex of simHits - Layer " + std::to_string(layer) + ";y [mm];Entries", + axis_xy + } + ); + hist1d.at(hist1d_simHit_Vertex_z).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_Vertex_z", + "Z-position of the production vertex of simHits - Layer " + std::to_string(layer) + ";z [mm];Entries", + axis_z + } + ); + + + hist1d.at(hist1d_simhit_MomentumDirection_x).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_MomentumDirection_x", + "X-component of the normalised momentum of simHits (at simHit position) - Layer " + std::to_string(layer) + ";Fraction of momentum in x direction [a.u.];Entries", + axis_MomentumFraction + } + ); + hist1d.at(hist1d_simhit_MomentumDirection_y).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_MomentumDirection_y", + "Y-component of the normalised momentum of simHits (at simHit position) - Layer " + std::to_string(layer) + ";Fraction of momentum in y direction [a.u.];Entries", + axis_MomentumFraction + } + ); + hist1d.at(hist1d_simhit_MomentumDirection_z).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_MomentumDirection_z", + "Z-component of the normalised momentum of simHits (at simHit position) - Layer " + std::to_string(layer) + ";Fraction of momentum in z direction [a.u.];Entries", + axis_MomentumFraction + } + ); + + hist1d.at(hist1d_simHit_InitialMomentumDirection_x).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_InitialMomentumDirection_x", + "X-component of the initial momentum of simHits (at their creation vertex) - Layer " + std::to_string(layer) + ";Fraction of momentum in x direction [a.u.];Entries", + axis_MomentumFraction + } + ); + hist1d.at(hist1d_simHit_InitialMomentumDirection_y).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_InitialMomentumDirection_y", + "Y-component of the initial momentum of simHits (at their creation vertex) - Layer " + std::to_string(layer) + ";Fraction of momentum in y direction [a.u.];Entries", + axis_MomentumFraction + } + ); + hist1d.at(hist1d_simHit_InitialMomentumDirection_z).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_InitialMomentumDirection_z", + "Z-component of the initial momentum of simHits (at their creation vertex) - Layer " + std::to_string(layer) + ";Fraction of momentum in z direction [a.u.];Entries", + axis_MomentumFraction + } + ); + + + + + + hist1d.at(hist1d_simHitPDG).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_PDG", + "SimHit particle PDG code - Layer " + std::to_string(layer) + ";PDG;Entries", + axis_pdg + } + ); + + + hist1d.at(hist1d_digiHitCharge).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_chargeDep", + "DigiHit charge - Layer " + std::to_string(layer) + ";Charge [e-];Entries", + axis_charge + } + ); + hist1d.at(hist1d_digiHitsPerSimHit).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_digiHits_per_simHit", + "Raw number of digiHits created per simHit - Layer " + std::to_string(layer) + ";DigiHits per SimHit;Entries", + axis_clusterSize // not technically correct, but works + } + ); + + hist1d.at(hist1d_clusterSize).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_clusterSize", + "Number of pixels per cluster (after clustering algorithm) - Layer " + std::to_string(layer) + ";Cluster size [pixels];Entries", + axis_clusterSize + } + ); + hist1d.at(hist1d_clusterSize_u).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_clusterSize_u", + "Cluster length in u direction (after clustering algorithm) - Layer " + std::to_string(layer) + ";Cluster size [pixels];Entries", + axis_clusterSize + } + ); + hist1d.at(hist1d_clusterSize_v).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_clusterSize_v", + "Cluster length in v direction (after clustering algorithm) - Layer " + std::to_string(layer) + ";Cluster size [pixels];Entries", + axis_clusterSize + } + ); + hist1d.at(hist1d_clusterSize_createdInGenerator).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_clusterSize_createdInGenerator", + "Number of pixels per cluster (for simHits created in generator) - Layer " + std::to_string(layer) + ";Cluster size [pixels];Entries", + axis_clusterSize + } + ); + hist1d.at(hist1d_clusterSize_createdInSimulation).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_clusterSize_createdInSimulation", + "Number of pixels per cluster (for simHits created in simulation) - Layer " + std::to_string(layer) + ";Cluster size [pixels];Entries", + axis_clusterSize + } + ); + + hist1d.at(hist1d_residual_u).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_residual_u", + "Residual (u_digiHit - u_simHit) in local u direction - Layer " + std::to_string(layer) + ";Residual u [um];Entries", + axis_residual + } + ); + hist1d.at(hist1d_residual_u_maxEParticleOnSensor).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_residual_u_maxEParticleOnSensor", + "Residual (u_digiHit - u_simHit) in local u direction, to the simHit with the highest energy MCParticle on the sensor - Layer " + std::to_string(layer) + ";Residual u [um];Entries", + axis_residual + } + ); + hist1d.at(hist1d_residual_u_noParents).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_residual_u_noParents", + "Residual (u_digiHit - u_simHit) in local u direction, for simHits without any parents - Layer " + std::to_string(layer) + ";Residual u [um];Entries", + axis_residual + } + ); + hist1d.at(hist1d_residual_u_directMcParticleLink).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_residual_u_directMcParticleLink", + "Residual (u_digiHit - u_simHit) in local u direction, for simHits with a direct MC particle link (vs via an intermediate MCParticle, that was removed to save space) - Layer " + std::to_string(layer) + ";Residual u [um];Entries", + axis_residual + } + ); + hist1d.at(hist1d_residual_u_mcOriginFarFromSimHit).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_residual_u_mcOriginFarFromSimHit", + "Residual (u_digiHit - u_simHit) in local u direction, for simHits from MCParticles that were created more than 500 um from the simHit - Layer " + std::to_string(layer) + ";Residual u [um];Entries", + axis_residual + } + ); + hist1d.at(hist1d_residual_u_directMcParticleLink_mcOriginFarFromSimHit).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_residual_u_directMcParticleLink_mcOriginFarFromSimHit", + "Residual (u_digiHit - u_simHit) in local u direction, for simHits with a direct MC particle link that were created more than 500 um from the simHit - Layer " + std::to_string(layer) + ";Residual u [um];Entries", + axis_residual + } + ); + hist1d.at(hist1d_residual_u_singlePixelCluster).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_residual_u_singlePixelCluster", + "Residual (u_digiHit - u_simHit) in local u direction for clusters with exactly one pixel - Layer " + std::to_string(layer) + ";Residual u [um];Entries", + axis_residual + } + ); + hist1d.at(hist1d_residual_u_multiPixelCluster).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_residual_u_multiPixelCluster", + "Residual (u_digiHit - u_simHit) in local u direction for clusters with more than one pixel - Layer " + std::to_string(layer) + ";Residual u [um];Entries", + axis_residual + } + ); + hist1d.at(hist1d_residual_u_createdInGenerator).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_residual_u_createdInGenerator", + "Residual (u_digiHit - u_simHit) in local u direction, from particles created in the generator - Layer " + std::to_string(layer) + ";Residual u [um];Entries", + axis_residual + } + ); + hist1d.at(hist1d_residual_u_createdInSimulation).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_residual_u_createdInSimulation", + "Residual (u_digiHit - u_simHit) in local u direction, from particles created in simulation - Layer " + std::to_string(layer) + ";Residual u [um];Entries", + axis_residual + } + ); + hist1d.at(hist1d_residual_v).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_residual_v", + "Residual (v_digiHit - v_simHit) in local v direction - Layer " + std::to_string(layer) + ";Residual v [um];Entries", + axis_residual + } + ); + hist1d.at(hist1d_residual_v_maxEParticleOnSensor).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_residual_v_maxEParticleOnSensor", + "Residual (v_digiHit - v_simHit) in local v direction, to the simHit with the highest energy MCParticle on the sensor - Layer " + std::to_string(layer) + ";Residual v [um];Entries", + axis_residual + } + ); + hist1d.at(hist1d_residual_v_noParents).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_residual_v_noParents", + "Residual (v_digiHit - v_simHit) in local v direction, for simHits without any parents - Layer " + std::to_string(layer) + ";Residual v [um];Entries", + axis_residual + } + ); + hist1d.at(hist1d_residual_v_directMcParticleLink).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_residual_v_directMcParticleLink", + "Residual (v_digiHit - v_simHit) in local v direction, for simHits with a direct MC particle link (vs via an intermediate MCParticle, that was removed to save space) - Layer " + std::to_string(layer) + ";Residual v [um];Entries", + axis_residual + } + ); + hist1d.at(hist1d_residual_v_mcOriginFarFromSimHit).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_residual_v_mcOriginFarFromSimHit", + "Residual (v_digiHit - v_simHit) in local v direction, for simHits with a direct MC particle link and MCParticles that were created more than 500 um from the simHit - Layer " + std::to_string(layer) + ";Residual v [um];Entries", + axis_residual + } + ); + hist1d.at(hist1d_residual_v_directMcParticleLink_mcOriginFarFromSimHit).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_residual_v_directMcParticleLink_mcOriginFarFromSimHit", + "Residual (v_digiHit - v_simHit) in local v direction, for simHits with a direct MC particle link and MCParticles that were created more than 500 um from the simHit - Layer " + std::to_string(layer) + ";Residual v [um];Entries", + axis_residual + } + ); + hist1d.at(hist1d_residual_v_singlePixelCluster).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_residual_v_singlePixelCluster", + "Residual (v_digiHit - v_simHit) in local v direction for clusters with exactly one pixel - Layer " + std::to_string(layer) + ";Residual v [um];Entries", + axis_residual + } + ); + hist1d.at(hist1d_residual_v_multiPixelCluster).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_residual_v_multiPixelCluster", + "Residual (v_digiHit - v_simHit) in local v direction for clusters with more than one pixel - Layer " + std::to_string(layer) + ";Residual v [um];Entries", + axis_residual + } + ); + hist1d.at(hist1d_residual_v_createdInGenerator).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_residual_v_createdInGenerator", + "Residual (v_digiHit - v_simHit) in local v direction, from particles created in the generator - Layer " + std::to_string(layer) + ";Residual v [um];Entries", + axis_residual + } + ); + hist1d.at(hist1d_residual_v_createdInSimulation).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_residual_v_createdInSimulation", + "Residual (v_digiHit - v_simHit) in local v direction, from particles created in simulation - Layer " + std::to_string(layer) + ";Residual v [um];Entries", + axis_residual + } + ); + hist1d.at(hist1d_residual_w).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_residual_w", + "Residual (w_digiHit - w_simHit) in local w direction (ie. sensor normal direction) - Layer " + std::to_string(layer) + ";Residual w [um];Entries", + axis_residual + } + ); + hist1d.at(hist1d_residual_w_noParents).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_residual_w_noParents", + "Residual (w_digiHit - w_simHit) in local w direction, for simHits without any parents - Layer " + std::to_string(layer) + ";Residual w [um];Entries", + axis_residual + } + ); + hist1d.at(hist1d_residual_r).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_residual_r", + "Residual magnitude (|r_digiHit - r_simHit|) - Layer " + std::to_string(layer) + ";Residual r [um];Entries", + axis_residual_abs + } + ); + + hist1d.at(hist1d_clusterPosUncertainty_u).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_clusterPosUncertainty_u", + "DigiHit position uncertainty in u direction - Layer " + std::to_string(layer) + ";Position uncertainty u [um];Entries", + axis_residual + } + ); + hist1d.at(hist1d_clusterPosUncertainty_u_noParents).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_clusterPosUncertainty_u_noParents", + "DigiHit position uncertainty in u direction, for simHits without any parents - Layer " + std::to_string(layer) + ";Position uncertainty u [um];Entries", + axis_residual + } + ); + hist1d.at(hist1d_clusterPosUncertainty_v).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_clusterPosUncertainty_v", + "DigiHit position uncertainty in v direction - Layer " + std::to_string(layer) + ";Position uncertainty v [um];Entries", + axis_residual + } + ); + hist1d.at(hist1d_clusterPosUncertainty_v_noParents).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_clusterPosUncertainty_v_noParents", + "DigiHit position uncertainty in v direction, for simHits without any parents - Layer " + std::to_string(layer) + ";Position uncertainty v [um];Entries", + axis_residual + } + ); + + hist1d.at(hist1d_pixelDistToClusterCenter_u).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_pixelDistanceToClusterCenter_u", + "Distance of pixel to center of cluster in u direction - Layer " + std::to_string(layer) + ";Distance to cluster center u [pix];Entries", + axis_residual_pixels + } + ); + hist1d.at(hist1d_pixelDistToClusterCenter_v).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_pixelDistanceToClusterCenter_v", + "Distance of pixel to center of cluster in v direction - Layer " + std::to_string(layer) + ";Distance to cluster center v [pix];Entries", + axis_residual_pixels + } + ); + + hist1d.at(hist1d_pathTravel_u).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_pathTravel_u", + "Path travel length inside the sensor active volume in u direction - Layer " + std::to_string(layer) + ";Path length in u [um];Entries", + axis_pathTravel + } + ); + hist1d.at(hist1d_pathTravel_v).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_pathTravel_v", + "Path travel length inside the sensor active volume in v direction - Layer " + std::to_string(layer) + ";Path length in v [um];Entries", + axis_pathTravel + } + ); + hist1d.at(hist1d_pathTravel_r).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_pathTravel_r", + "Path travel length inside the sensor active volume - Layer " + std::to_string(layer) + ";Path length [um];Entries", + axis_pathLength + } + ); + + hist1d.at(hist1d_simHitTimeStamp).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_timeStamp", + "SimHit time stamp - Layer " + std::to_string(layer) + ";Time [?s];Entries", + axis_time + } + ); + hist1d.at(hist1d_digiHitTimeStamp).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_timeStamp", + "Cluster time stamp - Layer " + std::to_string(layer) + ";Time [?s];Entries", + axis_time + } + ); + + + hist1d.at(hist1d_simHitsPerDigiHit).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_numberOfSimHits", + "Number of simHits per digiHit - Layer " + std::to_string(layer) + ";Number of SimHits;Entries", + axis_clusterSize // not technically correct, but works + } + ); + + hist1d.at(hist1d_highestEnergyParticleOnSensor_energy).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_highestEnergyParticleOnSensor_energy", + "Energy of the particle with the highest energy on the sensor - Layer " + std::to_string(layer) + ";Energy [GeV];Entries", + axis_particleE + } + ); + + + + m_hist1d.emplace(layer, std::move(hist1d)); + + /* -- 1d-profile-hist -- */ + + std::array< std::unique_ptr< Gaudi::Accumulators::StaticProfileHistogram<1,Gaudi::Accumulators::atomicity::full,float>>, histProfile1dArrayLen> histProfile1d; + + histProfile1d.at(histProfile1d_digiHitCharge_vs_global_z).reset( + new Gaudi::Accumulators::StaticProfileHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_chargeDep_vs_global_z", + "DigiHit charge - Layer " + std::to_string(layer) + ";SimHit global z position [mm];Charge [e-]", + axis_z + } + ); + + histProfile1d.at(histProfile1d_clusterSize_vs_global_z).reset( + new Gaudi::Accumulators::StaticProfileHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_clusterSize_vs_global_z", + "Cluster size - Layer " + std::to_string(layer) + ";SimHit global z position [mm];Pixels per cluster", + axis_z + } + ); + histProfile1d.at(histProfile1d_clusterSize_u_vs_global_z).reset( + new Gaudi::Accumulators::StaticProfileHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_clusterSize_u_vs_global_z", + "Cluster length along u - Layer " + std::to_string(layer) + ";SimHit global z position [mm];Cluster length in u [pix]", + axis_z + } + ); + histProfile1d.at(histProfile1d_clusterSize_v_vs_global_z).reset( + new Gaudi::Accumulators::StaticProfileHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_clusterSize_v_vs_global_z", + "Cluster length along v - Layer " + std::to_string(layer) + ";SimHit global z position [mm];Cluster length in v [pix]", + axis_z + } + ); + + histProfile1d.at(histProfile1d_residual_u_vs_global_z).reset( + new Gaudi::Accumulators::StaticProfileHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_residual_u_vs_global_z", + "Residual u (u_digiHit - u_simHit) vs. global z position - Layer " + std::to_string(layer) + ";SimHit global z position [mm];Residual u [um]", + axis_z + } + ); + histProfile1d.at(histProfile1d_residual_v_vs_global_z).reset( + new Gaudi::Accumulators::StaticProfileHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_residual_v_vs_global_z", + "Residual v (v_digiHit - v_simHit) vs. global z position - Layer " + std::to_string(layer) + ";SimHit global z position [mm];Residual v [um]", + axis_z + } + ); + histProfile1d.at(histProfile1d_residual_r_vs_global_z).reset( + new Gaudi::Accumulators::StaticProfileHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_residual_r_vs_global_z", + "Residual r (|r_digiHit - r_simHit|) vs. global z position - Layer " + std::to_string(layer) + ";SimHit global z position [mm];Residual r [um]", + axis_z + } + ); + + histProfile1d.at(histProfile1d_pathTravel_u_vs_global_z).reset( + new Gaudi::Accumulators::StaticProfileHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_pathTravel_u_vs_global_z", + "Path travel length inside the sensor active volume in u direction vs. global z position - Layer " + std::to_string(layer) + ";SimHit global z position [mm];Path length in u [um]", + axis_z + } + ); + histProfile1d.at(histProfile1d_pathTravel_v_vs_global_z).reset( + new Gaudi::Accumulators::StaticProfileHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_pathTravel_v_vs_global_z", + "Path travel length inside the sensor active volume in v direction vs. global z position - Layer " + std::to_string(layer) + ";SimHit global z position [mm];Path length in v [um]", + axis_z + } + ); + histProfile1d.at(histProfile1d_pathTravel_r_vs_global_z).reset( + new Gaudi::Accumulators::StaticProfileHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_pathTravel_r_vs_global_z", + "Path travel length inside the sensor active volume vs. global z position - Layer " + std::to_string(layer) + ";SimHit global z position [mm];Path length [um]", + axis_z + } + ); + + m_histProfile1d.emplace(layer, std::move(histProfile1d)); + + /* -- 2d-hist -- */ + + std::array< std::unique_ptr< Gaudi::Accumulators::StaticHistogram< 2, Gaudi::Accumulators::atomicity::full, float > >, hist2dArrayLen> hist2d; + + hist2d.at(hist2d_digiHitCharge_vs_global_z).reset( + new Gaudi::Accumulators::StaticHistogram<2, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_chargeDep_vs_global_z_2D", + "DigiHit charge - Layer " + std::to_string(layer) + ";Global z position [mm];Charge [e-]", + axis_z, + axis_charge + } + ); + + hist2d.at(hist2d_hitMap_simHits).reset( + new Gaudi::Accumulators::StaticHistogram<2, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_hitMap", + "SimHit hit map (pixel in which this simHit position lies) - Layer " + std::to_string(layer) + ";Pixel u;Pixel v;Entries", + axis_pixels_u, + axis_pixels_v + } + ); + hist2d.at(hist2d_hitMap_simHits_createdInGenerator).reset( + new Gaudi::Accumulators::StaticHistogram<2, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_hitMap_createdInGenerator", + "SimHit hit map (pixel in which this simHit position lies) for particles created in the generator - Layer " + std::to_string(layer) + ";Pixel u;Pixel v;Entries", + axis_pixels_u, + axis_pixels_v + } + ); + hist2d.at(hist2d_hitMap_simHits_createdInSimulation).reset( + new Gaudi::Accumulators::StaticHistogram<2, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_hitMap_createdInSimulation", + "SimHit hit map (pixel in which this simHit position lies) for particles created in the simulation - Layer " + std::to_string(layer) + ";Pixel u;Pixel v;Entries", + axis_pixels_u, + axis_pixels_v + } + ); + hist2d.at(hist2d_hitMap_simHits_eDepAboveThreshold).reset( + new Gaudi::Accumulators::StaticHistogram<2, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_hitMap_eDepAboveThreshold", + "SimHit hit map (pixel in which this simHit position lies) for simHits with total energy deposition above threshold - Layer " + std::to_string(layer) + ";Pixel u;Pixel v;Entries", + axis_pixels_u, + axis_pixels_v + } + ); + hist2d.at(hist2d_hitMap_pixelHits).reset( + new Gaudi::Accumulators::StaticHistogram<2, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_hitMap", + "Hit map (pixels that collect charge over the threshold) - Layer " + std::to_string(layer) + ";Pixel u;Pixel v;Entries", + axis_pixels_u, + axis_pixels_v + } + ); + hist2d.at(hist2d_clusterSize_vs_global_z).reset( + new Gaudi::Accumulators::StaticHistogram<2, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_clusterSize_vs_global_z_2D", + "Cluster size vs. global hit z - Layer " + std::to_string(layer) + ";Global z position [mm];Pixels per cluster;Entries", + axis_z, + axis_clusterSize + } + ); + hist2d.at(hist2d_clusterSize_u_vs_global_z).reset( + new Gaudi::Accumulators::StaticHistogram<2, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_clusterSize_u_vs_global_z_2D", + "Cluster length along u vs. global hit z - Layer " + std::to_string(layer) + ";Global z position [mm];Cluster length in u [pix];Entries", + axis_z, + axis_clusterSize + } + ); + hist2d.at(hist2d_clusterSize_v_vs_global_z).reset( + new Gaudi::Accumulators::StaticHistogram<2, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_clusterSize_v_vs_global_z_2D", + "Cluster length along v vs. global hit z - Layer " + std::to_string(layer) + ";Global z position [mm];Cluster length in v [pix];Entries", + axis_z, + axis_clusterSize + } + ); + hist2d.at(hist2d_clusterSize_vs_global_z_createdInGenerator).reset( + new Gaudi::Accumulators::StaticHistogram<2, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_clusterSize_vs_global_z_createdInGenerator_2D", + "Cluster size vs. global hit z (simHit particle created in generator) - Layer " + std::to_string(layer) + ";Global z position [mm];Pixels per cluster;Entries", + axis_z, + axis_clusterSize + } + ); + hist2d.at(hist2d_clusterSize_vs_global_z_createdInSim).reset( + new Gaudi::Accumulators::StaticHistogram<2, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_clusterSize_vs_global_z_createdInSimulation_2D", + "Cluster size vs. global hit z (simHit particle created in simulation) - Layer " + std::to_string(layer) + ";Global z position [mm];Pixels per cluster;Entries", + axis_z, + axis_clusterSize + } + ); + + hist2d.at(hist2d_residual_u_vs_global_z).reset( + new Gaudi::Accumulators::StaticHistogram<2, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_residual_u_vs_global_z_2D", + "Residual u (u_digiHit - u_simHit) vs. global z position - Layer " + std::to_string(layer) + ";Global z position [mm];Residual u [um];Entries", + axis_z, + axis_residual + } + ); + hist2d.at(hist2d_residual_v_vs_global_z).reset( + new Gaudi::Accumulators::StaticHistogram<2, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_residual_v_vs_global_z_2D", + "Residual v (v_digiHit - v_simHit) vs. global z position - Layer "+ std::to_string(layer) + ";Global z position [mm];Residual v [um];Entries", + axis_z, + axis_residual + } + ); + hist2d.at(hist2d_residual_r_vs_global_z).reset( + new Gaudi::Accumulators::StaticHistogram<2, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_residual_r_vs_global_z_2D", + "Residual r (|r_digiHit - r_simHit|) vs. global z position - Layer " + std::to_string(layer) + ";Global z position [mm];Residual r [um];Entries", + axis_z, + axis_residual_abs + } + ); + + hist2d.at(hist2d_residual_vs_clusterPosUncertainty_u).reset( + new Gaudi::Accumulators::StaticHistogram<2, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_residual_vs_clusterPosUncertainty_u_2D", + "Residual r (|r_digiHit - r_simHit|) vs. cluster position uncertainty in u direction - Layer " + std::to_string(layer) + ";Cluster position uncertainty u [um];Residual u [um];Entries", + axis_residual_abs, + axis_residual_abs + } + ); + hist2d.at(hist2d_residual_vs_clusterPosUncertainty_v).reset( + new Gaudi::Accumulators::StaticHistogram<2, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/digiHit_residual_vs_clusterPosUncertainty_v_2D", + "Residual r (|r_digiHit - r_simHit|) vs. cluster position uncertainty in v direction - Layer " + std::to_string(layer) + ";Cluster position uncertainty v [um];Residual v [um];Entries", + axis_residual_abs, + axis_residual_abs + } + ); + + hist2d.at(hist2d_pathTravel_u_vs_global_z).reset( + new Gaudi::Accumulators::StaticHistogram<2, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_pathTravel_u_vs_global_z_2D", + "Path travel length inside the sensor active volume in u direction vs. global z position - Layer " + std::to_string(layer) + ";Global z position [mm];Path length in u [um];Entries", + axis_z, + axis_pathTravel + } + ); + hist2d.at(hist2d_pathTravel_u_vs_global_z_createdInGenerator).reset( + new Gaudi::Accumulators::StaticHistogram<2, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_pathTravel_u_vs_global_z_createdInGenerator_2D", + "Path travel length inside the sensor active volume in u direction vs. global z position, for particles created in the generator - Layer " + std::to_string(layer) + ";Global z position [mm];Path length in u [um];Entries", + axis_z, + axis_pathTravel + } + ); + hist2d.at(hist2d_pathTravel_u_vs_global_z_createdInSimulation).reset( + new Gaudi::Accumulators::StaticHistogram<2, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_pathTravel_u_vs_global_z_createdInSimulation_2D", + "Path travel length inside the sensor active volume in u direction vs. global z position, for particles created in the simulation - Layer " + std::to_string(layer) + ";Global z position [mm];Path length in u [um];Entries", + axis_z, + axis_pathTravel + } + ); + hist2d.at(hist2d_pathTravel_v_vs_global_z).reset( + new Gaudi::Accumulators::StaticHistogram<2, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_pathTravel_v_vs_global_z_2D", + "Path travel length inside the sensor active volume in v direction vs. global z position - Layer " + std::to_string(layer) + ";Global z position [mm];Path length in v [um];Entries", + axis_z, + axis_pathTravel + } + ); + hist2d.at(hist2d_pathTravel_v_vs_global_z_createdInGenerator).reset( + new Gaudi::Accumulators::StaticHistogram<2, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_pathTravel_v_vs_global_z_createdInGenerator_2D", + "Path travel length inside the sensor active volume in v direction vs. global z position, for particles created in the generator - Layer " + std::to_string(layer) + ";Global z position [mm];Path length in v [um];Entries", + axis_z, + axis_pathTravel + } + ); + hist2d.at(hist2d_pathTravel_v_vs_global_z_createdInSimulation).reset( + new Gaudi::Accumulators::StaticHistogram<2, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_pathTravel_v_vs_global_z_createdInSimulation_2D", + "Path travel length inside the sensor active volume in v direction vs. global z position, for particles created in the simulation - Layer " + std::to_string(layer) + ";Global z position [mm];Path length in v [um];Entries", + axis_z, + axis_pathTravel + } + ); + hist2d.at(hist2d_pathTravel_w_vs_global_z).reset( + new Gaudi::Accumulators::StaticHistogram<2, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_pathTravel_w_vs_global_z_2D", + "Path travel length inside the sensor active volume in w direction vs. global z position - Layer " + std::to_string(layer) + ";Global z position [mm];Path length in w [um];Entries", + axis_z, + axis_pathTravel + } + ); + hist2d.at(hist2d_pathTravel_w_vs_global_z_createdInGenerator).reset( + new Gaudi::Accumulators::StaticHistogram<2, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_pathTravel_w_vs_global_z_createdInGenerator_2D", + "Path travel length inside the sensor active volume in w direction vs. global z position, for particles created in the generator - Layer " + std::to_string(layer) + ";Global z position [mm];Path length in w [um];Entries", + axis_z, + axis_pathTravel + } + ); + hist2d.at(hist2d_pathTravel_w_vs_global_z_createdInSimulation).reset( + new Gaudi::Accumulators::StaticHistogram<2, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_pathTravel_w_vs_global_z_createdInSimulation_2D", + "Path travel length inside the sensor active volume in w direction vs. global z position, for particles created in the simulation - Layer " + std::to_string(layer) + ";Global z position [mm];Path length in w [um];Entries", + axis_z, + axis_pathTravel + } + ); + hist2d.at(hist2d_pathTravel_r_vs_global_z).reset( + new Gaudi::Accumulators::StaticHistogram<2, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_pathTravel_r_vs_global_z_2D", + "Path travel length inside the sensor active volume vs. global z position - Layer " + std::to_string(layer) + ";Global z position [mm];Path length [um];Entries", + axis_z, + axis_pathLength + } + ); + + hist2d.at(hist2d_simHit_xy).reset( + new Gaudi::Accumulators::StaticHistogram<2, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_xy_2D", + "SimHit x vs y position (global) - Layer " + std::to_string(layer) + ";X position [mm];Y position [mm];Entries", + axis_z, + axis_z + } + ); + hist2d.at(hist2d_simHit_xz).reset( + new Gaudi::Accumulators::StaticHistogram<2, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_xz_2D", + "SimHit x vs z position (global) - Layer " + std::to_string(layer) + ";X position [mm];Z position [mm];Entries", + axis_z, + axis_z + } + ); + hist2d.at(hist2d_simHit_yz).reset( + new Gaudi::Accumulators::StaticHistogram<2, Gaudi::Accumulators::atomicity::full, float> {this, + "Layer" + std::to_string(layer) + "/simHit_yz_2D", + "SimHit y vs z position (global) - Layer " + std::to_string(layer) + ";Y position [mm];Z position [mm];Entries", + axis_z, + axis_z + } + ); + + m_hist2d.emplace(layer, std::move(hist2d)); + + } /* loop over layers */ + + /* global (eg covering all layers) histograms */ + + m_hist1dglobal.at(hist1dglobal_pathTravel_r).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Global/simHit_pathLength", + "Path travel length in sensor active volume, as computed by VTXdigi_tools reconstruction;Path length [um];Entries", + axis_pathLength + } + ); + m_hist1dglobal.at(hist1dglobal_pathTravel_r_Geant4).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Global/simHit_pathLength_Geant4", + "Path length in sensor active volume, as given by Geant4;Path length [um];Entries", + axis_pathLength + } + ); + m_hist1dglobal.at(hist1dglobal_pathTravel_r_ratio).reset( + new Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, float> {this, + "Global/simHit_pathLength_ratio", + "Path length in sensor active volume divided by path length given by Geant4;Path length / path length Geant4;Entries", + {500, 0.f, 2.f} + } + ); + +} + +/* ---- Eventloop functions ---- */ + +bool VTXdigi_Modular::CheckEventSetup(const edm4hep::SimTrackerHitCollection& simTrackerHits, const edm4hep::EventHeaderCollection& headers) const { + if (m_counter_eventsRead.value() % m_infoPrintInterval.value() == 0) + info() << "PROCESSING event [run " << headers.at(0).getRunNumber() << ", event " << headers.at(0).getEventNumber() << ", found " << simTrackerHits.size() << " simHits]. " << m_counter_eventsRead.value() << " events so far." << endmsg; + /* events are not necessarily numbered sequentially... */ + ++m_counter_eventsRead; + + if (simTrackerHits.size()==0) { + debug() << " - No SimTrackerHits in collection, returning empty output collections" << endmsg; + ++m_counter_eventsRejected_noSimHits; + return false; + } + + ++m_counter_eventsAccepted; + return true; +} + +bool VTXdigi_Modular::CheckSimhitLayer(const edm4hep::SimTrackerHit& simTrackerHit) const { + ++m_counter_simHitsRead; + const int layer = m_cellIdDecoder->get(simTrackerHit.getCellID(), "layer"); + if (m_layers.value().size()>0) { + if (std::find(m_layers.value().begin(), m_layers.value().end(), layer) == m_layers.value().end()) { + ++m_counter_simHitsRejected_LayerNotToBeDigitized; + return false; + } + } + + ++m_counter_simHitsAccepted; + return true; +} + +std::vector VTXdigi_Modular::Clusterize(const VTXdigi_tools::HitMap& hitMap) const { + if (!hitMap.Hits().size()) { + debug() << " - No pixels with charge found on this sensor." << endmsg; + return {}; + } + + if (m_clusterize.value()) { + debug() << " - Clusterizing " << hitMap.Hits().size() << " hits with a total charge of " << hitMap.GetTotalCharge() << " e." << endmsg; + return hitMap.ComputeClusters(); + } + else { + return hitMap.ComputeClusters_singePixels(); + } +} + +void VTXdigi_Modular::CreateDigiHits(edm4hep::TrackerHitPlaneCollection& digiHits, edm4hep::TrackerHitSimTrackerHitLinkCollection& digiHitLinks, const dd4hep::DDSegmentation::CellID& cellID, const TGeoHMatrix& trafoMatrix, const std::vector& clusters) const { + + // Called for each cluster on a sensor, so cellID and direction vectors are same among these clusters + const dd4hep::rec::Vector3D direction_u_3d = VTXdigi_tools::LocalToGlobal(dd4hep::rec::Vector3D(1, 0, 0), trafoMatrix); + const edm4hep::Vector2f direction_u = edm4hep::Vector2f(direction_u_3d.theta(), direction_u_3d.phi()); // edm4hep expects direction in to given as (theta, phi) + + const dd4hep::rec::Vector3D direction_v_3d = VTXdigi_tools::LocalToGlobal(dd4hep::rec::Vector3D(0, 1, 0), trafoMatrix); + const edm4hep::Vector2f direction_v = edm4hep::Vector2f(direction_v_3d.theta(), direction_v_3d.phi()); // edm4hep expects direction in to given as (theta, phi) + + for (auto& cluster : clusters) { + if (cluster.simHits.empty()) { + error() << "Cluster with no contributing simHits found." << endmsg; + } + if (cluster.charge <= 0) { + error() << "Cluster with non-positive charge found." << endmsg; + } + + edm4hep::MutableTrackerHitPlane digiHit = digiHits.create(); + ++m_counter_digiHitsCreated; + /* TODO: do we want cellID with segmentation? (instead of using the short cellID that only encodes the sensor and not the pixel) */ + digiHit.setCellID(cellID); + digiHit.setEDep(cluster.charge / VTXdigi_tools::kChargePerkeV); + + // position + const std::pair clusterPos_index = cluster.ComputePos(); + const dd4hep::rec::Vector3D clusterPos_local = VTXdigi_tools::ComputePosFromPixIndex_local(clusterPos_index, m_sensorLength, m_pixelPitch, m_chargeCollector->GetChargeCollectionDepthCenter()); + const dd4hep::rec::Vector3D clusterPos_global = VTXdigi_tools::LocalToGlobal(clusterPos_local, trafoMatrix); + debug() << " - Found cluster with " << cluster.pixels.size() << " pixels, charge " << cluster.charge << ", center at (" << clusterPos_index.first << ", " << clusterPos_index.second << "). Has " << cluster.simHits.size() << " contributing simHits." << endmsg; + digiHit.setPosition(VTXdigi_tools::ConvertVector(clusterPos_global)); + + // pos uncertainty + digiHit.setU(direction_u); + digiHit.setV(direction_v); + if (m_positionUncertainty.value().empty()) { + std::pair clusterPos_unc = cluster.ComputePosUncertainty_ChargeWeighted(clusterPos_index); + digiHit.setDu(clusterPos_unc.first * m_pixelPitch.first); + digiHit.setDv(clusterPos_unc.second * m_pixelPitch.second); + } + else if (m_positionUncertainty.value().size() == 2) { + digiHit.setDu(m_positionUncertainty.value().at(0)); + digiHit.setDv(m_positionUncertainty.value().at(1)); + } + else if (m_positionUncertainty.value().size() == 6) { + int clusterSize_u = cluster.GetSize(0); + if (clusterSize_u == 1) + digiHit.setDu(m_positionUncertainty.value().at(0)); + else if (clusterSize_u == 2) + digiHit.setDu(m_positionUncertainty.value().at(1)); + else + digiHit.setDu(m_positionUncertainty.value().at(2)); + + int clusterSize_v = cluster.GetSize(1); + if (clusterSize_v == 1) + digiHit.setDv(m_positionUncertainty.value().at(3)); + else if (clusterSize_v == 2) + digiHit.setDv(m_positionUncertainty.value().at(4)); + else + digiHit.setDv(m_positionUncertainty.value().at(5)); + } + debug() << " - Set digiHit position uncertainty to (" << digiHit.getDu() << ", " << digiHit.getDv() << ") mm." << endmsg; + + // collect timestamp + const VTXdigi_tools::Pixel* seedPixel = cluster.pixels.front(); + for (const VTXdigi_tools::Pixel* pix : cluster.pixels) { + if (pix->charge > seedPixel->charge) + seedPixel = pix; // use timestamp of pixel with the highest charge + } + float timeStamp = std::numeric_limits::max(); + for (const auto& simHit : seedPixel->simHits) { + timeStamp = std::min(timeStamp, simHit->hitPtr()->getTime()); // use earliest timestamp among simHits contributing to that pixel + } + if (m_smearing_time > 0.f) + timeStamp += m_rndm_time; + digiHit.setTime(timeStamp); + + // Create links to simHits + for (const auto& simHit : cluster.simHits) { + auto link = digiHitLinks.create(); + link.setFrom(digiHit); + link.setTo(*(simHit->hitPtr())); + } + + if (m_debugHistograms.value()) { + FillHistograms_perDigiHit(cluster.simHits, digiHit, trafoMatrix, cluster.pixels.size(), cluster.GetSize(0), cluster.GetSize(1)); + + for (const VTXdigi_tools::Pixel* pix : cluster.pixels) + FillHistograms_perPixel(cellID, *pix, clusterPos_index); + } + } /* loop over clusters */ +} + +/* ---- Histogramming functions ---- */ + +void VTXdigi_Modular::FillHistograms_perSimHit(const VTXdigi_tools::SimHitWrapper& simHit) const { + /* executed once for each simHit, no cuts applied before */ + + const int layer = simHit.layer(); + const TGeoHMatrix trafoMatrix = VTXdigi_tools::ComputeSensorTrafoMatrix(simHit.cellID(), m_volumeManager, m_sensorNormalRotation); + + const dd4hep::rec::Vector3D simHitPos_global = VTXdigi_tools::ConvertVector(simHit.hitPtr()->getPosition()); + const dd4hep::rec::Vector3D simHitPos_local = VTXdigi_tools::GlobalToLocal(simHitPos_global, trafoMatrix); + const dd4hep::rec::Vector3D simHitProdVertex_global = VTXdigi_tools::ConvertVector(simHit.hitPtr()->getParticle().getVertex()); + const std::pair i_uv = VTXdigi_tools::ComputePixelIndices(simHitPos_local, m_pixelPitch, m_pixelCount); + + const dd4hep::rec::Vector3D simHitMomentum = VTXdigi_tools::ConvertVector(simHit.hitPtr()->getMomentum()); // in GeV + const dd4hep::rec::Vector3D simHitMomentumInitial = VTXdigi_tools::ConvertVector(simHit.hitPtr()->getParticle().getMomentum()); // in GeV + + ++(*m_hist1d.at(layer).at(hist1d_simHitE))[simHit.hitPtr()->getEDep() * (dd4hep::GeV / dd4hep::keV)]; + ++(*m_hist1d.at(layer).at(hist1d_simHitCharge))[simHit.charge()]; + ++(*m_hist1d.at(layer).at(hist1d_simHitMomentum_keV))[simHitMomentum.r() * 1.E6]; + ++(*m_hist1d.at(layer).at(hist1d_simHitMomentum_MeV))[simHitMomentum.r() * 1.E3]; + ++(*m_hist1d.at(layer).at(hist1d_simHitMomentum_GeV))[simHitMomentum.r()]; + + ++(*m_hist1d.at(layer).at(hist1d_simHitTimeStamp))[simHit.hitPtr()->getTime()]; + + ++(*m_hist1d.at(layer).at(hist1d_simHit_x))[simHitPos_global.x()]; + ++(*m_hist1d.at(layer).at(hist1d_simHit_y))[simHitPos_global.y()]; + ++(*m_hist1d.at(layer).at(hist1d_simHit_z))[simHitPos_global.z()]; + + ++(*m_hist1d.at(layer).at(hist1d_simHit_Vertex_x))[simHitProdVertex_global.x()]; + ++(*m_hist1d.at(layer).at(hist1d_simHit_Vertex_y))[simHitProdVertex_global.y()]; + ++(*m_hist1d.at(layer).at(hist1d_simHit_Vertex_z))[simHitProdVertex_global.z()]; + + ++(*m_hist1d.at(layer).at(hist1d_simhit_MomentumDirection_x))[simHitMomentum.x() / simHitMomentum.r()]; + ++(*m_hist1d.at(layer).at(hist1d_simhit_MomentumDirection_y))[simHitMomentum.y() / simHitMomentum.r()]; + ++(*m_hist1d.at(layer).at(hist1d_simhit_MomentumDirection_z))[simHitMomentum.z() / simHitMomentum.r()]; + + ++(*m_hist1d.at(layer).at(hist1d_simHit_InitialMomentumDirection_x))[simHitMomentumInitial.x() / simHitMomentumInitial.r()]; + ++(*m_hist1d.at(layer).at(hist1d_simHit_InitialMomentumDirection_y))[simHitMomentumInitial.y() / simHitMomentumInitial.r()]; + ++(*m_hist1d.at(layer).at(hist1d_simHit_InitialMomentumDirection_z))[simHitMomentumInitial.z() / simHitMomentumInitial.r()]; + + if ( simHit.CreatedInGenerator() ) { + ++(*m_hist1d.at(layer).at(hist1d_simHit_z_createdInGenerator))[simHitPos_global.z()]; + ++(*m_hist2d.at(layer).at(hist2d_hitMap_simHits_createdInGenerator))[{i_uv.first, i_uv.second}]; + } + else { + ++(*m_hist1d.at(layer).at(hist1d_simHit_z_createdInSimulation))[simHitPos_global.z()]; + ++(*m_hist2d.at(layer).at(hist2d_hitMap_simHits_createdInSimulation))[{i_uv.first, i_uv.second}]; + } + + if (simHit.charge() >= m_threshold.value()) { + ++(*m_hist2d.at(layer).at(hist2d_hitMap_simHits_eDepAboveThreshold))[{i_uv.first, i_uv.second}]; + } + + ++(*m_hist1d.at(layer).at(hist1d_simHitPDG))[simHit.hitPtr()->getParticle().getPDG()]; + + ++(*m_hist2d.at(layer).at(hist2d_hitMap_simHits))[{i_uv.first, i_uv.second}]; + + ++(*m_hist2d.at(layer).at(hist2d_simHit_xy))[{simHitPos_global.x(), simHitPos_global.y()}]; + ++(*m_hist2d.at(layer).at(hist2d_simHit_xz))[{simHitPos_global.x(), simHitPos_global.z()}]; + ++(*m_hist2d.at(layer).at(hist2d_simHit_yz))[{simHitPos_global.y(), simHitPos_global.z()}]; +} + +void VTXdigi_Modular::FillHistograms_perPixel(const dd4hep::DDSegmentation::CellID& cellID, const VTXdigi_tools::Pixel& pix, const std::pair clusterPos_local) const { + /* executed once for each pixel */ + + const int layer = VTXdigi_tools::GetLayer(cellID, m_cellIdDecoder); + const std::pair i_uv = pix.index; + + ++(*m_hist2d.at(layer).at(hist2d_hitMap_pixelHits))[{i_uv.first, i_uv.second}]; + + if (clusterPos_local.first != 0 && clusterPos_local.second != 0) { + /* cluster at (0,0) means that hits weren't clustered */ + const float distToClusterCenter_u = static_cast(pix.index.first) - clusterPos_local.first; // in pix + const float distToClusterCenter_v = static_cast(pix.index.second) - clusterPos_local.second; + + ++(*m_hist1d.at(layer).at(hist1d_pixelDistToClusterCenter_u))[distToClusterCenter_u]; // convert to um + ++(*m_hist1d.at(layer).at(hist1d_pixelDistToClusterCenter_v))[distToClusterCenter_v]; + } +} + +void VTXdigi_Modular::FillHistograms_perDigiHit(const std::unordered_set& simHits, const edm4hep::TrackerHitPlane& digiHit, const TGeoHMatrix& trafoMatrix, const int clusterSize, const int clusterSize_u, const int clusterSize_v) const { + /* executed once for each digiHit */ + const int layer = VTXdigi_tools::GetLayer(digiHit.getCellID(), m_cellIdDecoder); + const dd4hep::rec::Vector3D pos_global = VTXdigi_tools::ConvertVector(digiHit.getPosition()); + const dd4hep::rec::Vector3D pos_local = VTXdigi_tools::GlobalToLocal(pos_global, trafoMatrix); + + ++(*m_hist1d.at(layer).at(hist1d_digiHitCharge))[digiHit.getEDep() * VTXdigi_tools::kChargePerkeV]; + + ++(*m_hist1d.at(layer).at(hist1d_clusterSize))[clusterSize]; + ++(*m_hist1d.at(layer).at(hist1d_clusterSize_u))[clusterSize_u]; + ++(*m_hist1d.at(layer).at(hist1d_clusterSize_v))[clusterSize_v]; + + (*m_histProfile1d.at(layer).at(histProfile1d_clusterSize_vs_global_z))[pos_global.z()] += clusterSize; + (*m_histProfile1d.at(layer).at(histProfile1d_clusterSize_u_vs_global_z))[pos_global.z()] += clusterSize_u; + (*m_histProfile1d.at(layer).at(histProfile1d_clusterSize_v_vs_global_z))[pos_global.z()] += clusterSize_v; + + ++(*m_hist2d.at(layer).at(hist2d_clusterSize_vs_global_z))[{pos_global.z(), clusterSize}]; + ++(*m_hist2d.at(layer).at(hist2d_clusterSize_u_vs_global_z))[{pos_global.z(), clusterSize_u}]; + ++(*m_hist2d.at(layer).at(hist2d_clusterSize_v_vs_global_z))[{pos_global.z(), clusterSize_v}]; + + ++(*m_hist1d.at(layer).at(hist1d_clusterPosUncertainty_u))[ digiHit.getDu() * 1000.f ]; // convert to um + ++(*m_hist1d.at(layer).at(hist1d_clusterPosUncertainty_v))[ digiHit.getDv() * 1000.f ]; + + ++(*m_hist1d.at(layer).at(hist1d_digiHitTimeStamp))[ digiHit.getTime() ]; + + ++(*m_hist1d.at(layer).at(hist1d_simHitsPerDigiHit))[ simHits.size() ]; + + // per simhit that is contributing to this digiHit (cluster) + for (const auto& simHit : simHits) { + const edm4hep::MCParticle mcParticle = simHit->hitPtr()->getParticle(); + + const dd4hep::rec::Vector3D simHitPos_local = simHit->truthPos(); + const dd4hep::rec::Vector3D simHitPos_global = VTXdigi_tools::LocalToGlobal(simHitPos_local, trafoMatrix); + const dd4hep::rec::Vector3D residual_local = pos_local - simHitPos_local; // residual = observed - predicted + + const float hit_z = simHitPos_global.z(); + + (*m_histProfile1d.at(layer).at(histProfile1d_residual_u_vs_global_z))[hit_z] += residual_local.x()*1000.f; // convert from mm to um + (*m_histProfile1d.at(layer).at(histProfile1d_residual_v_vs_global_z))[hit_z] += residual_local.y()*1000.f; + (*m_histProfile1d.at(layer).at(histProfile1d_residual_r_vs_global_z))[hit_z] += residual_local.r()*1000.f; + + ++(*m_hist2d.at(layer).at(hist2d_residual_u_vs_global_z))[{hit_z, residual_local.x()*1000.f}]; + ++(*m_hist2d.at(layer).at(hist2d_residual_v_vs_global_z))[{hit_z, residual_local.y()*1000.f}]; + ++(*m_hist2d.at(layer).at(hist2d_residual_r_vs_global_z))[{hit_z, residual_local.r()*1000.f}]; + + ++(*m_hist1d.at(layer).at(hist1d_residual_u))[ residual_local.x()*1000.f ]; + ++(*m_hist1d.at(layer).at(hist1d_residual_v))[ residual_local.y()*1000.f ]; + ++(*m_hist1d.at(layer).at(hist1d_residual_w))[ residual_local.z()*1000.f ]; + ++(*m_hist1d.at(layer).at(hist1d_residual_r))[ residual_local.r()*1000.f ]; + + ++(*m_hist2d.at(layer).at(hist2d_residual_vs_clusterPosUncertainty_u))[{std::abs(digiHit.getDu()) * 1000.f, std::abs(residual_local.x())*1000.f}]; // convert to um + ++(*m_hist2d.at(layer).at(hist2d_residual_vs_clusterPosUncertainty_v))[{std::abs(digiHit.getDv()) * 1000.f, std::abs(residual_local.y())*1000.f}]; + + if (clusterSize == 1) { + ++(*m_hist1d.at(layer).at(hist1d_residual_u_singlePixelCluster))[ residual_local.x()*1000.f ]; + ++(*m_hist1d.at(layer).at(hist1d_residual_v_singlePixelCluster))[ residual_local.y()*1000.f ]; + } + else if (clusterSize > 1) { + ++(*m_hist1d.at(layer).at(hist1d_residual_u_multiPixelCluster))[ residual_local.x()*1000.f ]; + ++(*m_hist1d.at(layer).at(hist1d_residual_v_multiPixelCluster))[ residual_local.y()*1000.f ]; + } + + if ( simHit->CreatedInGenerator() ) { + ++(*m_hist1d.at(layer).at(hist1d_clusterSize_createdInGenerator))[clusterSize]; + ++(*m_hist2d.at(layer).at(hist2d_clusterSize_vs_global_z_createdInGenerator))[{hit_z, clusterSize}]; + + ++(*m_hist1d.at(layer).at(hist1d_residual_u_createdInGenerator))[ residual_local.x()*1000.f ]; + ++(*m_hist1d.at(layer).at(hist1d_residual_v_createdInGenerator))[ residual_local.y()*1000.f ]; + } + else { + ++(*m_hist1d.at(layer).at(hist1d_clusterSize_createdInSimulation))[clusterSize]; + ++(*m_hist2d.at(layer).at(hist2d_clusterSize_vs_global_z_createdInSim))[{hit_z, clusterSize}]; + + ++(*m_hist1d.at(layer).at(hist1d_residual_u_createdInSimulation))[ residual_local.x()*1000.f ]; + ++(*m_hist1d.at(layer).at(hist1d_residual_v_createdInSimulation))[ residual_local.y()*1000.f ]; + } + + if (mcParticle.parents_size() == 0) { + ++(*m_hist1d.at(layer).at(hist1d_residual_u_noParents))[ residual_local.x()*1000.f ]; + ++(*m_hist1d.at(layer).at(hist1d_residual_v_noParents))[ residual_local.y()*1000.f ]; + ++(*m_hist1d.at(layer).at(hist1d_residual_w_noParents))[ residual_local.z()*1000.f ]; + + ++(*m_hist1d.at(layer).at(hist1d_clusterPosUncertainty_u_noParents))[ digiHit.getDu() * 1000.f ]; // convert to um + ++(*m_hist1d.at(layer).at(hist1d_clusterPosUncertainty_v_noParents))[ digiHit.getDv() * 1000.f ]; + } + + const bool hasDirectMcParticleLink = simHit->hasDirectMcParticleLink(); + + const dd4hep::rec::Vector3D mcParticleProdVertex_global = VTXdigi_tools::ConvertVector(mcParticle.getVertex()); + const dd4hep::rec::Vector3D mcParticleProdVertex_difference = pos_global - mcParticleProdVertex_global; + const float mcParticleDistance = mcParticleProdVertex_difference.r(); + + const bool mcOriginFarFromSimHit = mcParticleDistance > 0.5f; // in mm + + if (hasDirectMcParticleLink) { + ++(*m_hist1d.at(layer).at(hist1d_residual_u_directMcParticleLink))[ residual_local.x()*1000.f ]; + ++(*m_hist1d.at(layer).at(hist1d_residual_v_directMcParticleLink))[ residual_local.y()*1000.f ]; + } + if (mcOriginFarFromSimHit) { + ++(*m_hist1d.at(layer).at(hist1d_residual_u_mcOriginFarFromSimHit))[ residual_local.x()*1000.f ]; + ++(*m_hist1d.at(layer).at(hist1d_residual_v_mcOriginFarFromSimHit))[ residual_local.y()*1000.f ]; + } + if (hasDirectMcParticleLink && mcOriginFarFromSimHit) { + ++(*m_hist1d.at(layer).at(hist1d_residual_u_directMcParticleLink_mcOriginFarFromSimHit))[ residual_local.x()*1000.f ]; + ++(*m_hist1d.at(layer).at(hist1d_residual_v_directMcParticleLink_mcOriginFarFromSimHit))[ residual_local.y()*1000.f ]; + } + } // loop over contributing simHits +} + +void VTXdigi_Modular::FillHistograms_perSensor(const std::vector& simHits, const edm4hep::TrackerHitPlaneCollection& digiHits, const TGeoHMatrix& trafoMatrix, const dd4hep::DDSegmentation::CellID& cellID) const { + /* executed once for each sensor, after all clusters have been created */ + const int layer = VTXdigi_tools::GetLayer(cellID, m_cellIdDecoder); + + if (simHits.empty()) { + debug() << " - No simHits found on this sensor." << endmsg; + return; + } + + // among simHits, find the one with the MCParticle with the highest energy + // (needed for comparing to Allpix Squared residuals, when simulating single sensor with particle gun) + size_t maxE_index = 0; + float maxE = simHits.at(0).hitPtr()->getParticle().getEnergy(); // assume this is in GeV. Documentation is not clear + + for (size_t i = 1; i < simHits.size(); ++i) { + const VTXdigi_tools::SimHitWrapper& simHit = simHits.at(i); + const edm4hep::MCParticle mcParticle = simHit.hitPtr()->getParticle(); + if (mcParticle.getEnergy() > maxE && simHit.hasDirectMcParticleLink()) { + // we want to find the primary particle (assuming we only shot a single particle in the simulation) -> neede to reject simHits with indirect MCParticle link + maxE_index = i; + maxE = mcParticle.getEnergy(); + } + } // loop to find simHit with highest MCParticle energy. + + const VTXdigi_tools::SimHitWrapper& simHit = simHits.at(maxE_index); + const dd4hep::rec::Vector3D simHitPos_local = simHit.truthPos(); + + ++(*m_hist1d.at(layer).at(hist1d_highestEnergyParticleOnSensor_energy))[ maxE ]; // in GeV + + for (const auto& digiHit : digiHits) { + const dd4hep::rec::Vector3D pos_global = VTXdigi_tools::ConvertVector(digiHit.getPosition()); + const dd4hep::rec::Vector3D pos_local = VTXdigi_tools::GlobalToLocal(pos_global, trafoMatrix); + + const dd4hep::rec::Vector3D residual_local = pos_local - simHitPos_local; // residual = observed - predicted + + ++(*m_hist1d.at(layer).at(hist1d_residual_u_maxEParticleOnSensor))[ residual_local.x()*1000.f ]; + ++(*m_hist1d.at(layer).at(hist1d_residual_v_maxEParticleOnSensor))[ residual_local.y()*1000.f ]; + } +} + +void VTXdigi_Modular::FillHistograms_fromChargeCollector_perSimHit(const int layer, const dd4hep::rec::Vector3D& pathTravel, const float pathLength_Geant4, const dd4hep::rec::Vector3D& truthPos_local, const TGeoHMatrix& trafoMatrix, const bool createdInGenerator) const { + if (!m_debugHistograms.value()) return; + + const float pathLength = pathTravel.r(); // in mm + const float factor_um_per_mm = 1000.f; + const dd4hep::rec::Vector3D truthPos_global = VTXdigi_tools::LocalToGlobal(truthPos_local, trafoMatrix); + + ++(*m_hist1dglobal.at(hist1dglobal_pathTravel_r))[pathLength*factor_um_per_mm]; // convert from mm to um + ++(*m_hist1dglobal.at(hist1dglobal_pathTravel_r_Geant4))[pathLength_Geant4*factor_um_per_mm]; + if (pathLength_Geant4 != 0.f) { + ++(*m_hist1dglobal.at(hist1dglobal_pathTravel_r_ratio))[pathLength / pathLength_Geant4]; + } + + ++(*m_hist1d.at(layer).at(hist1d_pathTravel_u))[pathTravel.x()*factor_um_per_mm]; + ++(*m_hist1d.at(layer).at(hist1d_pathTravel_v))[pathTravel.y()*factor_um_per_mm]; + ++(*m_hist1d.at(layer).at(hist1d_pathTravel_r))[pathLength*factor_um_per_mm]; + + (*m_histProfile1d.at(layer).at(histProfile1d_pathTravel_u_vs_global_z))[truthPos_global.z()] += pathTravel.x()*factor_um_per_mm; + (*m_histProfile1d.at(layer).at(histProfile1d_pathTravel_v_vs_global_z))[truthPos_global.z()] += pathTravel.y()*factor_um_per_mm; + (*m_histProfile1d.at(layer).at(histProfile1d_pathTravel_r_vs_global_z))[truthPos_global.z()] += pathLength*factor_um_per_mm; + + ++(*m_hist2d.at(layer).at(hist2d_pathTravel_u_vs_global_z))[{truthPos_global.z(), pathTravel.x()*factor_um_per_mm}]; + ++(*m_hist2d.at(layer).at(hist2d_pathTravel_v_vs_global_z))[{truthPos_global.z(), pathTravel.y()*factor_um_per_mm}]; + ++(*m_hist2d.at(layer).at(hist2d_pathTravel_w_vs_global_z))[{truthPos_global.z(), pathTravel.z()*factor_um_per_mm}]; + ++(*m_hist2d.at(layer).at(hist2d_pathTravel_r_vs_global_z))[{truthPos_global.z(), pathLength*factor_um_per_mm}]; + + if (createdInGenerator) { + ++(*m_hist2d.at(layer).at(hist2d_pathTravel_u_vs_global_z_createdInGenerator))[{truthPos_global.z(), pathTravel.x()*factor_um_per_mm}]; + ++(*m_hist2d.at(layer).at(hist2d_pathTravel_v_vs_global_z_createdInGenerator))[{truthPos_global.z(), pathTravel.y()*factor_um_per_mm}]; + ++(*m_hist2d.at(layer).at(hist2d_pathTravel_w_vs_global_z_createdInGenerator))[{truthPos_global.z(), pathTravel.z()*factor_um_per_mm}]; + } + else { + ++(*m_hist2d.at(layer).at(hist2d_pathTravel_u_vs_global_z_createdInSimulation))[{truthPos_global.z(), pathTravel.x()*factor_um_per_mm}]; + ++(*m_hist2d.at(layer).at(hist2d_pathTravel_v_vs_global_z_createdInSimulation))[{truthPos_global.z(), pathTravel.y()*factor_um_per_mm}]; + ++(*m_hist2d.at(layer).at(hist2d_pathTravel_w_vs_global_z_createdInSimulation))[{truthPos_global.z(), pathTravel.z()*factor_um_per_mm}]; + } +} + +void VTXdigi_Modular::PrintCountersSummary() const { + const int colWidths[] = {65, 10}; + info() << " Counters summary: " << endmsg; + info() << " | " << std::setw(colWidths[0]) << std::left << "Events read" + << " | " << std::setw(colWidths[1]) << std::right << m_counter_eventsRead.value() << " |" << endmsg; + info() << " | " << std::setw(colWidths[0]) << std::left << "Events rejected (no simHits)" + << " | " << std::setw(colWidths[1]) << std::right << m_counter_eventsRejected_noSimHits.value() << " |" << endmsg; + info() << " | " << std::setw(colWidths[0]) << std::left << "Events accepted" + << " | " << std::setw(colWidths[1]) << std::right << m_counter_eventsAccepted.value() << " |" << endmsg; + + info() << " | " << std::setw(colWidths[0]) << std::left << "SimTrackerHits read" + << " | " << std::setw(colWidths[1]) << std::right << m_counter_simHitsRead.value() << " |" << endmsg; + info() << " | " << std::setw(colWidths[0]) << std::left << "SimTrackerHits rejected (layer ignored)" + << " | " << std::setw(colWidths[1]) << std::right << m_counter_simHitsRejected_LayerNotToBeDigitized.value() << " |" << endmsg; + info() << " | " << std::setw(colWidths[0]) << std::left << "SimTrackerHits accepted" + << " | " << std::setw(colWidths[1]) << std::right << m_counter_simHitsAccepted.value() << " |" << endmsg; + info() << " | " << std::setw(colWidths[0]) << std::left << "( accepted SimTrackerHits from particles created in generator )" + << " | " << std::setw(colWidths[1]) << std::right << m_counter_accSimHitsCreatedInGenerator.value() << " |" << endmsg; + info() << " | " << std::setw(colWidths[0]) << std::left << "( accepted SimTrackerHits with direct link to MCParticle )" + << " | " << std::setw(colWidths[1]) << std::right << m_counter_accSimHitsDirectMcParticleLink.value() << " |" << endmsg; + + info() << " | " << std::setw(colWidths[0]) << std::left << "Digi hits created" + << " | " << std::setw(colWidths[1]) << std::right << m_counter_digiHitsCreated.value() << " |" << endmsg; +} \ No newline at end of file diff --git a/VTXdigi_Modular/src/VTXdigi_tools.cpp b/VTXdigi_Modular/src/VTXdigi_tools.cpp new file mode 100644 index 00000000..dc3e652f --- /dev/null +++ b/VTXdigi_Modular/src/VTXdigi_tools.cpp @@ -0,0 +1,742 @@ +// VTXdigi_Modular/src/VTXdigi_tools.cpp +#include "VTXdigi_tools.h" + +namespace VTXdigi_tools { + +SimHitWrapper::SimHitWrapper(edm4hep::SimTrackerHit simTrackerHit, const std::unique_ptr& cellIdDecoder) : m_simTrackerHit(simTrackerHit) { + m_cellID = GetCellID_short(m_simTrackerHit); + + m_charge = static_cast(m_simTrackerHit.getEDep() * (dd4hep::GeV / dd4hep::keV) * kChargePerkeV); // convert energy deposit (in keV) to number of electrons + + m_layerNumber = GetLayer(m_cellID, cellIdDecoder); +} + +void swap(SimHitWrapper& a, SimHitWrapper& b) noexcept { + std::swap(a.m_simTrackerHit, b.m_simTrackerHit); + std::swap(a.m_cellID, b.m_cellID); + std::swap(a.m_charge, b.m_charge); + std::swap(a.m_layerNumber, b.m_layerNumber); +} // swap(Hit&, Hit&) + +// std::unordered_set FindSimHitsWithIndividualParents(const std::unordered_set simHits) { +// // TODO: implement this function, which is needed to find simHits from different MCParticles. If two simHits originate from the same MCParticle, return only the one from the MCParticle further up the family tree. + +// std::unordered_set individualSimHits; + +// for (const auto* simHit : simHits) { +// // if no simHit in individualSimHits shares a parent with simHit, add simHit to individualSimHits + +// for (const auto* otherSimHit : individualSimHits) { +// const edm4hep::MCParticle& mcP = simHit->hitPtr()->getParticle(); +// const edm4hep::MCParticle& otherMcP = otherSimHit->hitPtr()->getParticle(); + +// // first: check obvious cases +// if (mcP == otherMcP) { +// // no need to do anything, simHits are from the same MCParticle +// break; +// } + +// // second: check if simHit and otherSimHit are direct children of each other? +// if (mcP.parents_size() == 1) { +// const edm4hep::MCParticle& mcP_parent = mcP.getParents()[0]; +// if (mcP_parent == otherMcP) { +// // no need to do anything, other is created by the parent of simHit +// break; +// } +// } +// if (otherMcP.parents_size() == 1) { +// const edm4hep::MCParticle& otherMcP_parent = otherMcP.getParents()[0]; +// if (otherMcP_parent == mcP) { +// // replace otherSimHit with simHit, because simHit is from a parent of the otherSimHit + +// break; +// } +// } + + + +// if (simHit->hitPtr()->getParticle().getParents()[0] == otherSimHit->hitPtr()->getParticle().getParents()[0]) { +// hasSharedParent = true; +// break; +// } +// } +// if (!hasSharedParent) { +// individualSimHits.insert(simHit); +// } +// } + +// } + + +// SimulatorStatus bits (see https://edm4hep.web.cern.ch/classedm4hep_1_1_mutable_m_c_particle.html) +// 29 : "Backscatter", +// 30 : "CreatedInSimulation", +// 26 : "DecayedInCalorimeter", +// 27 : "DecayedInTracker", +// 22 : "HandledInFastSim", +// 25 : "LeftWorld", +// 23 : "Overlay", +// 24 : "Stopped", +// 28 : "VertexIsNotEndpointOfParent", + +/* -- helpers -- */ + +/** @brief Check if a MCParticle was created in the generator */ +bool CreatedInGenerator(const edm4hep::MCParticle& mcParticle) { + int32_t simulatorStatus = mcParticle.getSimulatorStatus(); + int mask_bit = edm4hep::MCParticle::BITCreatedInSimulation; // should be bit 30 + const int32_t mask = 1 << mask_bit; + + if ((simulatorStatus & mask) == 0) { + return true; // bit is not set -> created in generator + } + return false; // bit is set -> created in simulation +} + +dd4hep::rec::Vector3D ConvertVector(edm4hep::Vector3d vec) { + return dd4hep::rec::Vector3D(vec.x, vec.y, vec.z); +} +dd4hep::rec::Vector3D ConvertVector(edm4hep::Vector3f vec) { + return dd4hep::rec::Vector3D(static_cast(vec.x), static_cast(vec.y), static_cast(vec.z)); +} +edm4hep::Vector3d ConvertVector(dd4hep::rec::Vector3D vec) { + return edm4hep::Vector3d(vec.x(), vec.y(), vec.z()); +} + +TGeoHMatrix ComputeSensorTrafoMatrix(const dd4hep::DDSegmentation::CellID& cellID, const dd4hep::VolumeManager& volumeManager, const TGeoRotation& sensorNormalRotation) { + TGeoHMatrix M = volumeManager.lookupDetElement(cellID).nominal().worldTransformation(); + + /* rotate the local coordinate system st. sensor U is (1,0,0), V is (0,1,0) and normal vector is (0,0,1) */ + M.Multiply(sensorNormalRotation); + + /* rotation is unitless, but need to convert translation from cm to mm (dd4hep::mm = 0.1) */ + double* transl = M.GetTranslation(); + transl[0] = transl[0] / dd4hep::mm; + transl[1] = transl[1] / dd4hep::mm; + transl[2] = transl[2] / dd4hep::mm; + M.SetTranslation(transl); + + return M; +} + +dd4hep::rec::Vector3D GlobalToLocal(const dd4hep::rec::Vector3D& global, const TGeoHMatrix& M) { + double local[3]; + M.MasterToLocal(global, local); + return dd4hep::rec::Vector3D(local[0], local[1], local[2]); +} + +dd4hep::rec::Vector3D LocalToGlobal(const dd4hep::rec::Vector3D& local, const TGeoHMatrix& M) { + double global[3]; + M.LocalToMaster(local, global); + return dd4hep::rec::Vector3D(global[0], global[1], global[2]); +} + +dd4hep::DDSegmentation::CellID GetCellID_short(const edm4hep::SimTrackerHit& simTrackerHit) { + /* Mask removes segmentation bits, now cellID is unique for each sensor */ + std::uint64_t m_mask = (static_cast(1) << 32) - 1; + return simTrackerHit.getCellID() & m_mask; +} +dd4hep::DDSegmentation::CellID GetCellID_short(const dd4hep::DDSegmentation::CellID& cellID) { + /* Mask removes segmentation bits, now cellID is unique for each sensor */ + std::uint64_t m_mask = (static_cast(1) << 32) - 1; + return cellID & m_mask; +} + +int GetLayer(const dd4hep::DDSegmentation::CellID& cellID, const std::unique_ptr& cellIdDecoder) { + return static_cast(cellIdDecoder->get(cellID, "layer")); +} +int GetLayer(const edm4hep::SimTrackerHit& simTrackerHit, const std::unique_ptr& cellIdDecoder) { + return GetLayer(GetCellID_short(simTrackerHit), cellIdDecoder); +} + +/* -- Binning things -- */ + +int ComputeBinIndex(float x, float binX0, float binWidth, int binN) { + /** Get the bin index for a given x value + * binX0 is the lower edge of the first bin + * binWidth is the width of the bins + * binN is the number of bins + * return -1 if x is out of range + */ + + if (binN <= 0) throw std::runtime_error("VTXdigi_tools::ComputeBinIndex(): binN must be positive"); + if (binWidth <= 0.0) throw std::runtime_error("VTXdigi_tools::ComputeBinIndex(): binWidth must be positive"); + + float relativePos = (x - binX0) / binWidth; // shift to [0, binN] + if (relativePos < 0.0f || relativePos > static_cast(binN)) + return -1; + if (relativePos == static_cast(binN)) + return binN - 1; // include upper edge in last bin (makes sense for pixels) + return static_cast(relativePos); +} // ComputeBinIndex() + +float ComputeBinCenter(int i, float binX0, float binWidth) { + if (i < 0) throw std::runtime_error("VTXdigi_tools::ComputeBinCenter(): bin index must be non-negative"); + if (binWidth <= 0.0) throw std::runtime_error("VTXdigi_tools::ComputeBinCenter(): binWidth must be positive"); + + return binX0 + (static_cast(i) + 0.5f) * binWidth; // add 0.5*binWidth to shift from lower edge to center +} // ComputeBinCenter() +float ComputeBinCenter(int i, float binX0, float binX1, int binN) { + if (binN <= 0) throw std::runtime_error("VTXdigi_tools::ComputeBinCenter(): binN must be positive"); + if (binX1 <= binX0) throw std::runtime_error("VTXdigi_tools::ComputeBinCenter(): binX1 must be greater than binX0"); + if (i < 0 || i >= binN) throw std::runtime_error("VTXdigi_tools::ComputeBinCenter(): bin index out of bounds"); + + const float binWidth = (binX1 - binX0) / static_cast(binN); + return ComputeBinCenter(i, binX0, binWidth); +} // ComputeBinCenter() + +std::pair ComputePixelIndices(const dd4hep::rec::Vector3D& pos, const std::pair pixelPitch, const std::pair pixelCount) { + const float length_u_half = 0.5 * pixelPitch.first * pixelCount.first; + int i_u = ComputeBinIndex( + pos.x(), + -length_u_half, + pixelPitch.first, + pixelCount.first); + + const float length_v_half = 0.5 * pixelPitch.second * pixelCount.second; + int i_v = ComputeBinIndex( + pos.y(), + -length_v_half, + pixelPitch.second, + pixelCount.second); + + return {i_u, i_v}; +} // ComputePixelIndices() + +std::array ComputeInPixelIndices(const dd4hep::rec::Vector3D& pos, const std::array& binCount, const std::pair& pixelPitch, const std::array& activeVolumeDimensions) { + std::array indices; + + const float posShifted_u = pos.x() + 0.5 * activeVolumeDimensions[0]; // shift to [0, length_u] + if (posShifted_u < 0.0 || posShifted_u > activeVolumeDimensions[0]) { + indices[0] = -1; // out of bounds + } + else { + float posInPixel_u = std::fmod(posShifted_u, pixelPitch.first); + if (posInPixel_u < 0.0) posInPixel_u += pixelPitch.first; // ensure positive remainder + indices[0] = ComputeBinIndex(posInPixel_u, 0.0, pixelPitch.first / binCount[0], binCount[0]); + } + + const float posShifted_v = pos.y() + 0.5 * activeVolumeDimensions[1]; + if (posShifted_v < 0.0 || posShifted_v > activeVolumeDimensions[1]) { + indices[1] = -1; // out of bounds + } + else { + float posInPixel_v = std::fmod(posShifted_v, pixelPitch.second); + if (posInPixel_v < 0.0) posInPixel_v += pixelPitch.second; + indices[1] = ComputeBinIndex(posInPixel_v, 0.0, pixelPitch.second / binCount[1], binCount[1]); + } + + // vertical (w) binning: shift to [0, thickness] + const float posShifted_w = pos.z() + 0.5 * activeVolumeDimensions[2]; + indices[2] = ComputeBinIndex(posShifted_w, 0.0, activeVolumeDimensions[2] / binCount[2], binCount[2]); // no fmod, so out-of-bounds is caught + + return indices; +} // ComputeInPixelIndices() + +dd4hep::rec::Vector3D ComputePosFromPixIndex_local(const std::pair pixelIndex, const std::pair sensorLength, const std::pair pixelPitch, float depletedRegionDepthCenter) { + /* returns the position of the center of pixel i_u, i_v in the local sensor frame */ + + float u = (static_cast(pixelIndex.first) + 0.5f) * pixelPitch.first - 0.5f * sensorLength.first; // in mm + float v = (static_cast(pixelIndex.second) + 0.5f) * pixelPitch.second - 0.5f * sensorLength.second; + float w = depletedRegionDepthCenter; + + return dd4hep::rec::Vector3D(u, v, w); +} + +dd4hep::rec::Vector3D ComputePosFromPixIndex_local(const std::pair pixelIndex, const std::pair sensorLength, const std::pair pixelPitch) { + return ComputePosFromPixIndex_local(pixelIndex, sensorLength, pixelPitch, 0.f); +} + +dd4hep::rec::Vector3D ComputePosFromPixIndex_local(const std::pair index, const std::pair sensorLength, const std::pair pixelPitch, float depletedRegionDepthCenter) { + /* returns the position of the center of pixel i_u, i_v in the local sensor frame */ + + float u = (index.first + 0.5f) * pixelPitch.first - 0.5f * sensorLength.first; // in mm. Add 0.5*pixelPitch to shift from pixel edge to center, since index 0 is defined as the center of the pixel. + float v = (index.second + 0.5f) * pixelPitch.second - 0.5f * sensorLength.second; + float w = depletedRegionDepthCenter; + + return dd4hep::rec::Vector3D(u, v, w); +} + +dd4hep::rec::Vector3D ComputePosFromPixIndex_local(const std::pair index, const std::pair sensorLength, const std::pair pixelPitch) { + return ComputePosFromPixIndex_local(index, sensorLength, pixelPitch, 0.f); +} + +/* -- HitMap -- */ + +HitMap::HitMap(std::pair pixelCount) : m_pixCount(pixelCount) { + const int inverseOccupancy = 2000; // assume occupancy, 5e-4 is quite conservative for Z-run + m_pixels.reserve(pixelCount.first * pixelCount.second / inverseOccupancy); // avoid too many reallocations +} + +void HitMap::FillCharge(std::pair i_uv, float charge, const SimHitWrapper& simHitWrapper) { + if (charge < 1.e-6f) + return; // skip very small charge additions for performance (this is NECESSARY to skip in-pix bins with weight ~0) + if (_OutOfBounds(i_uv)) [[unlikely]] + throw std::runtime_error("HitMap::FillCharge: pixel i_u or i_v ( " + std::to_string(i_uv.first) + ", " + std::to_string(i_uv.second) + ") out of range"); + + auto [iter, inserted] = m_pixels.try_emplace(i_uv, Pixel(i_uv)); + iter->second.charge += charge; + iter->second.simHits.insert(&simHitWrapper); +} + +void HitMap::ApplyChargeSmearing(const Rndm::Numbers& rndm_charge) { + auto hitIter = m_pixels.begin(); + while (hitIter != m_pixels.end()) { + hitIter->second.charge = std::max(hitIter->second.charge + static_cast(rndm_charge()), 0.f); // don't allow negative charge after smearing + ++hitIter; + } +} + +void HitMap::ApplyThreshold(const float threshold) { + auto hitIter = m_pixels.begin(); + while (hitIter != m_pixels.end()) { + if (hitIter->second.charge < threshold) + hitIter = m_pixels.erase(hitIter); // erase returns the iterator to the next element, so this is safe to do while iterating + else + ++hitIter; + } +} + +float HitMap::GetCharge(std::pair i_uv) const { + if (_OutOfBounds(i_uv)) [[unlikely]] { + throw std::runtime_error("HitMap::GetCharge: pixel i_u or i_v ( " + std::to_string(i_uv.first) + ", " + std::to_string(i_uv.second) + ") out of range"); + } + auto it = m_pixels.find(i_uv); + if (it == m_pixels.end()) + return 0.f; // if pixel not found, charge is 0 + return it->second.charge; +} + +float HitMap::GetTotalCharge() const { + float totalCharge = 0.f; + for (const auto& [i_uv, pixHit] : m_pixels) { + totalCharge += pixHit.charge; + } + return totalCharge; +} + +inline bool HitMap::_OutOfBounds(std::pair i_uv) const { + return ( + i_uv.first < 0 + || i_uv.first >= static_cast(m_pixCount.first) + || i_uv.second < 0 + || i_uv.second >= static_cast(m_pixCount.second) + ); +} + +/* -- Clusterization -- */ + +std::pair Cluster::ComputePos() const { + std::pair pos{0.f, 0.f}; + for (const Pixel* pix : pixels) { + pos.first += pix->index.first * pix->charge; + pos.second += pix->index.second * pix->charge; + } + pos.first /= charge; + pos.second /= charge; + return pos; +} + +int Cluster::GetSize(const int axis) const { + int min = std::numeric_limits::max(); + int max = std::numeric_limits::min(); + + if (axis == 0) { // u + for (const Pixel* pix : pixels) { + min = std::min(min, pix->index.first); + max = std::max(max, pix->index.first); + } + } + else if (axis == 1) { // v + for (const Pixel* pix : pixels) { + min = std::min(min, pix->index.second); + max = std::max(max, pix->index.second); + } + } + else { + throw std::runtime_error("Cluster::GetClusterSize: axis must be 0 (u) or 1 (v), got " + std::to_string(axis)); + } + + return max - min + 1; // +1 because of counting: if min=max, cluster size is 1, not 0 +} + +std::pair Cluster::ComputePosUncertainty_ChargeWeighted(const std::pair& clusterPos) const { + float sig2_u=0.f, sig2_v=0.f; + for (const Pixel* pix : pixels) { + float du = (pix->index.first - clusterPos.first); + float dv = (pix->index.second - clusterPos.second); + sig2_u += pix->charge * du * du; + sig2_v += pix->charge * dv * dv; + } + sig2_u /= charge; + sig2_v /= charge; + return {std::sqrt(sig2_u), std::sqrt(sig2_v)}; +} + +std::pair Cluster::ComputePosUncertainty_ChargeWeighted() const { + return ComputePosUncertainty_ChargeWeighted(ComputePos()); +} + + +std::array, 4> GetDirectNeighbors(const std::pair& i_uv) { + return {{ + {i_uv.first - 1, i_uv.second}, // left + {i_uv.first + 1, i_uv.second}, // right + {i_uv.first, i_uv.second - 1}, // down + {i_uv.first, i_uv.second + 1} // up + }}; +} + +std::vector HitMap::ComputeClusters_singePixels() const { + std::vector clusters; + + for (const auto& p : m_pixels) { + const Pixel* pixel = &(p.second); // cluster stores pointers to pixels (pixels are stored in HitMap as values) + + clusters.emplace_back(); + clusters.back().pixels.push_back(pixel); + clusters.back().charge = pixel->charge; + for (const SimHitWrapper* simHitWrapper : pixel->simHits) { + clusters.back().simHits.insert(simHitWrapper); + } + } // loop over pixelHits + + return clusters; +} + +std::vector HitMap::ComputeClusters() const { + /* Breadth First Search (BFS) implementation for clustering */ + + std::vector clusters; + std::unordered_set, Hash_PairInt> visited; + + for (const auto& p : m_pixels) { + const std::pair seed_uv = p.first; + if (visited.contains(seed_uv)) + continue; + + clusters.emplace_back(); // create new cluster + clusters.back().pixels.reserve(10); // 10 should include >90% of clusters. i guess. + + std::queue> queue; + queue.push(seed_uv); + visited.insert(seed_uv); + + while (!queue.empty()) { + const std::pair current_uv = queue.front(); + queue.pop(); + + /* Add pixl to cluster */ + const Pixel* pixel = &(m_pixels.at(current_uv)); // get pixel pointer from map + clusters.back().pixels.push_back(pixel); + clusters.back().charge += pixel->charge; + for (const SimHitWrapper* simHitWrapper : pixel->simHits) { + clusters.back().simHits.insert(simHitWrapper); + } + /* Add all neighboring pixels to queue */ + for (const auto& neighbor_uv : GetDirectNeighbors(current_uv)) { + if (!m_pixels.contains(neighbor_uv)) + continue; + if (visited.contains(neighbor_uv)) + continue; + queue.push(neighbor_uv); + visited.insert(neighbor_uv); + } + + } // loop over queue + } // loop over cluster-seeds + return clusters; +} + +/* -- Tool tests -- */ + +bool ToolTest() { + std::cout << " | Running VTXdigi tool tests" << std::endl; + bool passed = true; + + std::cout << " | VTXdigi_tools::ComputeBinIndex()"; + { + bool passedInternal = true; + const float binX0 = -1.0f; + const float binWidth = 0.5f; + const int binN = 6; + + const float inputs[5] = { -1.0f, 0.0f, 2.0f, -1.1f, 2.1f }; + const int expectedOutputs[5] = { 0, 2, 5, -1, -1 }; + + for (size_t i = 0; i < 5; ++i) { + int result = ComputeBinIndex(inputs[i], binX0, binWidth, binN); + if (result != expectedOutputs[i]) { + std::cout << " - FAILED " << std::endl << " | -> Expected bin index " << expectedOutputs[i] << " for x=" << inputs[i] << ", got " << result << std::endl; + passedInternal = false; + } + } + if (passedInternal) + std::cout << " - PASSED" << std::endl; + passed = passed && passedInternal; + } + + std::cout << " | VTXdigi_tools::ComputePixelIndices()"; + { + bool passedInternal = true; + + const std::pair pixelPitch = { 1.0f, 2.0f }; + const std::pair pixelCount = { 10, 10 }; + + std::array inputs = {{dd4hep::rec::Vector3D(0.f, 0.f, 0.f)}}; + inputs = { + dd4hep::rec::Vector3D( -4.5f, -9.0f, 0.0f ), + dd4hep::rec::Vector3D(0.0f, 0.0f, 0.0f), + dd4hep::rec::Vector3D(4.4f, 8.9f, 0.0f), + dd4hep::rec::Vector3D(-5.0f, -10.0f, 0.0f), // edge tests + dd4hep::rec::Vector3D(5.0f, 10.0f, 0.0f), + dd4hep::rec::Vector3D(-5.1f, 0.0f, 0.0f), // out of bounds tests + dd4hep::rec::Vector3D(5.1f, 0.0f, 0.0f), + dd4hep::rec::Vector3D(0.0f, -10.1f, 0.0f), + dd4hep::rec::Vector3D(0.0f, 10.1f, 0.0f), + dd4hep::rec::Vector3D(5.1f, 10.1f, 0.0f), + }; + std::array, inputs.size()> expectedOutputs = {{{0, 0}}}; + expectedOutputs= {{ + {0,0}, + {5,5}, + {9,9}, + {0,0}, + {9,9}, + {-1,5}, + {-1,5}, + {5,-1}, + {5,-1}, + {-1,-1} + }}; + + for (size_t i = 0; i < inputs.size(); ++i) { + std::pair result = ComputePixelIndices(inputs.at(i), pixelPitch, pixelCount); + if (result != expectedOutputs[i]) { + std::cout << " - FAILED " << std::endl << " | -> Expected pixel indices (" << expectedOutputs[i].first << ", " << expectedOutputs[i].second << ") for (u,v,w)=(" << inputs.at(i).x() << ", " << inputs.at(i).y() << ", " << inputs.at(i).z() << "), got (" << result.first << ", " << result.second << ")" << std::endl; + passedInternal = false; + } + } + + if (passedInternal) + std::cout << " - PASSED" << std::endl; + passed = passed && passedInternal; + } + + std::cout << " | VTXdigi_tools::ComputePosFromPixIndex_local(std::pair)"; + { + bool passedInternal = true; + + const std::pair sensorLength = { 10.0, 20.0 }; // 10 x 10 pixels + const std::pair pixelPitch = { 1.0, 2.0 }; + + std::array, 6> inputs = {{{0, 0}}}; + inputs = {{ + {0, 0}, + {4, 4}, + {5, 5}, + {9, 9}, + {-1, -1}, // out of bounds test + {10, 10}, + }}; + std::array expectedOutputs = {{dd4hep::rec::Vector3D(0.f, 0.f, 0.f)}}; + expectedOutputs = {{ + dd4hep::rec::Vector3D( -4.5, -9.0, 0.0 ), + dd4hep::rec::Vector3D( -0.5, -1.0, 0.0 ), + dd4hep::rec::Vector3D( 0.5, 1.0, 0.0 ), + dd4hep::rec::Vector3D( 4.5, 9.0, 0.0 ), + dd4hep::rec::Vector3D( -5.5, -11.0, 0.0 ), + dd4hep::rec::Vector3D( 5.5, 11.0, 0.0 ), + }}; + + for (size_t i = 0; i < inputs.size(); ++i) { + dd4hep::rec::Vector3D result = ComputePosFromPixIndex_local(inputs.at(i), sensorLength, pixelPitch); + if (std::abs(result.x() - expectedOutputs[i].x()) > 1e-6 || std::abs(result.y() - expectedOutputs[i].y()) > 1e-6) { + if (passedInternal) + std::cout << " - FAILED " << std::endl; + std::cout << " | -> Expected local position (" << expectedOutputs[i].x() << ", " << expectedOutputs[i].y() << ") for index (i_u,i_v)=(" << inputs.at(i).first << ", " << inputs.at(i).second << "), got (" << result.x() << ", " << result.y() << ")" << std::endl; + passedInternal = false; + } + } + + if (passedInternal) + std::cout << " - PASSED" << std::endl; + passed = passed && passedInternal; + } + + std::cout << " | VTXdigi_tools::ComputePosFromPixIndex_local(std::pair)"; + { + bool passedInternal = true; + + const std::pair pixelPitch = { 1.0, 2.0 }; + const std::pair sensorLength = { 10.0, 20.0 }; // 10 x 10 pixels + + std::array, 9> inputs = {{{0., 0.}}}; + inputs = {{ + {0., 0.}, + {4., 4.}, + {4.5, 4.5}, + {5., 5.}, + {9.0, 9.0}, + {-0.5, -0.5}, + {9.5, 9.5}, + {-1., -1.}, // expect out of bounds to not be treated at all (done elsewhere) + {10., 10.}, + }}; + std::array expectedOutputs = {{dd4hep::rec::Vector3D(0.f, 0.f, 0.f)}}; + expectedOutputs = {{ + dd4hep::rec::Vector3D( -4.5, -9.0, 0.0 ), + dd4hep::rec::Vector3D( -0.5, -1.0, 0.0 ), + dd4hep::rec::Vector3D( 0., 0., 0. ), + dd4hep::rec::Vector3D( 0.5, 1.0, 0.0 ), + dd4hep::rec::Vector3D( 4.5, 9.0, 0.0 ), + dd4hep::rec::Vector3D( -5.0, -10.0, 0.0 ), + dd4hep::rec::Vector3D( 5.0, 10.0, 0.0 ), + dd4hep::rec::Vector3D( -5.5, -11.0, 0.0 ), + dd4hep::rec::Vector3D( 5.5, 11.0, 0.0 ), + }}; + + for (size_t i = 0; i < inputs.size(); ++i) { + dd4hep::rec::Vector3D result = ComputePosFromPixIndex_local(inputs.at(i), sensorLength, pixelPitch); + if (std::abs(result.x() - expectedOutputs[i].x()) > 1e-6 || std::abs(result.y() - expectedOutputs[i].y()) > 1e-6) { + if (passedInternal) + std::cout << " - FAILED " << std::endl; + std::cout << " | -> Expected local position (" << expectedOutputs[i].x() << ", " << expectedOutputs[i].y() << ") for index (i_u,i_v)=(" << inputs.at(i).first << ", " << inputs.at(i).second << "), got (" << result.x() << ", " << result.y() << ")" << std::endl; + passedInternal = false; + } + } + + if (passedInternal) + std::cout << " - PASSED" << std::endl; + passed = passed && passedInternal; + } + + std::cout << " | VTXdigi_tools::ComputeInPixelIndices()"; + { + bool passedInternal = true; + + const std::array binCount = { 4, 4, 4 }; + const std::pair pixelPitch = { 1.0f, 2.0f }; // bin-width is 0.25 + + const std::array activeVolumeDimensions = { 10.0f, 20.0f, 1.0f }; + + std::array inputs = {{dd4hep::rec::Vector3D(0.f, 0.f, 0.f)}}; + inputs = { + dd4hep::rec::Vector3D(0.1, 0.1, -0.4), + dd4hep::rec::Vector3D(0.1, 1.9, 0.1), + dd4hep::rec::Vector3D(-0.1, -0.1, -0.4), + dd4hep::rec::Vector3D(-0.1, -0.1, 0.1), + dd4hep::rec::Vector3D(0., 0., 0.), // edges + dd4hep::rec::Vector3D(1., 2.,0.5), + dd4hep::rec::Vector3D(-5.1, -10.1, 0.), // out-of-bounds + dd4hep::rec::Vector3D(5.1, 10.1, 0.), + dd4hep::rec::Vector3D(0.1, 0.1, -0.6), + dd4hep::rec::Vector3D(0.1, 0.1, 0.6), + }; + std::array, inputs.size()> expectedOutputs = {{{0, 0, 0}}}; + expectedOutputs = {{ + {0, 0, 0}, + {0, 3, 2}, + {3, 3, 0}, + {3, 3, 2}, + {0, 0, 2}, // edges + {0, 0, 3}, + {-1, -1, 2}, // out-of-bounds + {-1, -1, 2}, + {0, 0, -1}, + {0, 0, -1}, + }}; + + for (size_t i = 0; i < inputs.size(); ++i) { + std::array result = ComputeInPixelIndices(inputs.at(i), binCount, pixelPitch, activeVolumeDimensions); + if (result != expectedOutputs[i]) { + std::cout << " - FAILED " << std::endl << " | -> Expected in-pixel indices (" << expectedOutputs[i][0] << ", " << expectedOutputs[i][1] << ", " << expectedOutputs[i][2] << ") for (u,v,w)=(" << inputs.at(i).x() << ", " << inputs.at(i).y() << ", " << inputs.at(i).z() << "), got (" << result[0] << ", " << result[1] << ", " << result[2] << ")" << std::endl; + passedInternal = false; + } + } + + if (passedInternal) + std::cout << " - PASSED" << std::endl; + passed = passed && passedInternal; + } + + std::cout << " | VTXdigi_tools::HitMap()"; + { + bool passedInternal = true; + + SimHitWrapper simHitWrapper; // not used, just need something to pass to FillCharge() + + HitMap hitMap({10,10}); + + hitMap.FillCharge({0,0},1, simHitWrapper); + hitMap.FillCharge({1,1},2, simHitWrapper); + + if (hitMap.GetCharge({0,0}) != 1) { + std::cout << " - FAILED " << std::endl << " | -> Expected charge 1 at (0,0), got " << hitMap.GetCharge({0,0}) << std::endl; + passedInternal = false; + } + if (hitMap.GetCharge({1,1}) != 2) { + std::cout << " - FAILED " << std::endl << " | -> Expected charge 2 at (1,1), got " << hitMap.GetCharge({1,1}) << std::endl; + passedInternal = false; + } + if (hitMap.GetTotalCharge() != 3) { + std::cout << " - FAILED " << std::endl << " | -> Expected total charge 3, got " << hitMap.GetTotalCharge() << std::endl; + passedInternal = false; + } + + hitMap.FillCharge({0,0},100, simHitWrapper); // test adding charge to existing pixel + + if (hitMap.GetCharge({0,0}) != 101) { + std::cout << " - FAILED " << std::endl << " | -> Expected charge 101 at (0,0), got " << hitMap.GetCharge({0,0}) << std::endl; + passedInternal = false; + } + if (hitMap.GetTotalCharge() != 103) { + std::cout << " - FAILED " << std::endl << " | -> Expected total charge 103, got " << hitMap.GetTotalCharge() << std::endl; + passedInternal = false; + } + + try { + hitMap.FillCharge({-1,0},1, simHitWrapper); + std::cout << " - FAILED " << std::endl << " | -> Expected out-of-bounds exception for pixel (-1,0)" << std::endl; + passedInternal = false; + } + catch (const std::runtime_error& e) {} + try { + hitMap.FillCharge({0,-1},1, simHitWrapper); + std::cout << " - FAILED " << std::endl << " | -> Expected out-of-bounds exception for pixel (0,-1)" << std::endl; + passedInternal = false; + } + catch (const std::runtime_error& e) {} + try { + hitMap.FillCharge({10,0},1, simHitWrapper); + std::cout << " - FAILED " << std::endl << " | -> Expected out-of-bounds exception for pixel (10,0)" << std::endl; + passedInternal = false; + } + catch (const std::runtime_error& e) {} + try { + hitMap.FillCharge({0,10},1, simHitWrapper); + std::cout << " - FAILED " << std::endl << " | -> Expected out-of-bounds exception for pixel (0,10)" << std::endl; + passedInternal = false; + } + catch (const std::runtime_error& e) {} + + hitMap.Reset(); + + if (hitMap.GetTotalCharge() != 0) { + std::cout << " - FAILED " << std::endl << " | -> Expected total charge 0 after Reset(), got " << hitMap.GetTotalCharge() << std::endl; + passedInternal = false; + } + + if (passedInternal) + std::cout << " - PASSED" << std::endl; + passed = passed && passedInternal; + } + + if (passed) + std::cout << " | All tests passed." << std::endl; + return passed; +} + +} // namespace VTXdigi_tools + + diff --git a/VTXdigi_Modular/test/steering_digi.py b/VTXdigi_Modular/test/steering_digi.py new file mode 100644 index 00000000..e69de29b diff --git a/VTXdigi_Modular/test/steering_sim.py b/VTXdigi_Modular/test/steering_sim.py new file mode 100644 index 00000000..ccbbc297 --- /dev/null +++ b/VTXdigi_Modular/test/steering_sim.py @@ -0,0 +1,637 @@ +from DDSim.DD4hepSimulation import DD4hepSimulation +from g4units import mm, GeV, MeV + +################################### +# user options +simulateCalo = False # set to False to skip the calo SD action +################################### + +SIM = DD4hepSimulation() +SIM.printLevel = 3 +SIM.runType = "batch" + +## Lorentz boost for the crossing angle, in radian! +SIM.crossingAngleBoost = 0.015 + +SIM.enableDetailedShowerMode = False +SIM.enableG4GPS = False +SIM.enableG4Gun = False +SIM.enableGun = False + +## Physics list to use in simulation +simulateThinSilicon = True # set to False to not use settings for simulation of thin silicon sensor in barrel + +if not simulateThinSilicon: + SIM.physicsList = "FTFP_BERT" +else: + SIM.physicsList = "FTFP_BERT_EMZ" + +simulateThinSiliconDetectors = ["VertexBarrel"] + +## FourVector of translation for the Smearing of the Vertex position: x y z t +SIM.vertexOffset = [0.0, 0.0, 0.0, 0.0] +## FourVector of the Sigma for the Smearing of the Vertex position: x y z t +SIM.vertexSigma = [0.0, 0.0, 0.0, 0.0] + +SIM.action.calo = "Geant4VoidSensitiveAction" +SIM.action.mapActions["DRcalo"] = "Geant4VoidSensitiveAction" + +## set the default event action +SIM.action.event = [] + +## Set the drift chamber action +SIM.action.mapActions['DCH_v2'] = "Geant4TrackerAction" + +## set the default run action +SIM.action.run = [] + +## set the default stack action +SIM.action.stack = [] + +## set the default step action +SIM.action.step = [] + +## set the default track action +SIM.action.track = [] + +## set the default tracker action +SIM.action.tracker = ( + "Geant4TrackerWeightedAction", + {"HitPositionCombination": 2, "CollectSingleDeposits": False}, +) + +## set a different tracker action for silicon detectors that are thin and require single deposits for accurate delta electron simulation +for region in simulateThinSiliconDetectors: + SIM.action.mapActions[region] = ( + "Geant4TrackerWeightedAction", + {"HitPositionCombination": 2, "CollectSingleDeposits": True}, + ) + +## List of patterns matching sensitive detectors of type Tracker. +SIM.action.trackerSDTypes = ["tracker"] + + +################################################################################ +## Configuration for the magnetic field (stepper) +################################################################################ +SIM.field.delta_chord = 0.25 +SIM.field.delta_intersection = 0.001 +SIM.field.delta_one_step = 0.01 +SIM.field.eps_max = 0.001 +SIM.field.eps_min = 5e-05 +SIM.field.equation = "Mag_UsualEqRhs" +SIM.field.largest_step = 10000.0 +SIM.field.min_chord_step = 0.01 +SIM.field.stepper = "ClassicalRK4" + + +################################################################################ +## Configuration for sensitive detector filters +## +## Set the default filter for 'tracker' +## >>> SIM.filter.tracker = "edep1kev" +## Use no filter for 'calorimeter' by default +## >>> SIM.filter.calo = "" +## +## Assign a filter to a sensitive detector via pattern matching +## >>> SIM.filter.mapDetFilter['FTD'] = "edep1kev" +## +## Or more than one filter: +## >>> SIM.filter.mapDetFilter['FTD'] = ["edep1kev", "geantino"] +## +## Don't use the default filter or anything else: +## >>> SIM.filter.mapDetFilter['TPC'] = None ## or "" or [] +## +## Create a custom filter. The dictionary is used to instantiate the filter later on +## >>> SIM.filter.filters['edep3kev'] = dict(name="EnergyDepositMinimumCut/3keV", parameter={"Cut": 3.0*keV} ) +## +## +################################################################################ + +## +## default filter for calorimeter sensitive detectors; +## this is applied if no other filter is used for a calorimeter +## +# note: do not turn on the calo filter, otherwise all optical photons will be killed! +SIM.filter.calo = "" + +## list of filter objects: map between name and parameter dictionary +SIM.filter.filters = { + "geantino": {"name": "GeantinoRejectFilter/GeantinoRejector", "parameter": {}}, + "edep1kev": {"name": "EnergyDepositMinimumCut", "parameter": {"Cut": 0.001}}, + "edep0": {"name": "EnergyDepositMinimumCut/Cut0", "parameter": {"Cut": 0.0}}, +} + +## a map between patterns and filter objects, using patterns to attach filters to sensitive detector +SIM.filter.mapDetFilter["DCH_v2"] = "edep0" +SIM.filter.mapDetFilter["Muon-System"] = "edep0" + +## default filter for tracking sensitive detectors; this is applied if no other filter is used for a tracker +SIM.filter.tracker = "edep1kev" + + +################################################################################ +## Configuration for the Detector Construction. +################################################################################ +SIM.geometry.dumpGDML = "" +SIM.geometry.dumpHierarchy = 0 + +## Print Debug information about Elements +SIM.geometry.enableDebugElements = False + +## Print Debug information about Materials +SIM.geometry.enableDebugMaterials = False + +## Print Debug information about Placements +SIM.geometry.enableDebugPlacements = False + +## Print Debug information about Reflections +SIM.geometry.enableDebugReflections = False + +## Print Debug information about Regions +SIM.geometry.enableDebugRegions = False + +## Print Debug information about Shapes +SIM.geometry.enableDebugShapes = False + +## Print Debug information about Surfaces +SIM.geometry.enableDebugSurfaces = False + +## Print Debug information about Volumes +SIM.geometry.enableDebugVolumes = False + +## Print information about placements +SIM.geometry.enablePrintPlacements = False + +## Print information about Sensitives +SIM.geometry.enablePrintSensitives = False + +################################################################################ +## Configuration for the GuineaPig InputFiles +################################################################################ + +## Set the number of pair particles to simulate per event. +## Only used if inputFile ends with ".pairs" +## If "-1" all particles will be simulated in a single event +## +SIM.guineapig.particlesPerEvent = "-1" + + +################################################################################ +## Configuration for the DDG4 ParticleGun +################################################################################ + +## direction of the particle gun, 3 vector +SIM.gun.direction = (1.0, 0.1, 0.1) + +## choose the distribution of the random direction for theta +## +## Options for random distributions: +## +## 'uniform' is the default distribution, flat in theta +## 'cos(theta)' is flat in cos(theta) +## 'eta', or 'pseudorapidity' is flat in pseudorapity +## 'ffbar' is distributed according to 1+cos^2(theta) +## +## Setting a distribution will set isotrop = True +## +SIM.gun.distribution = None + +## Total energy (including mass) for the particle gun. +## +## If not None, it will overwrite the setting of momentumMin and momentumMax +SIM.gun.energy = 10.0 * GeV + +## Maximal pseudorapidity for random distibution (overrides thetaMin) +SIM.gun.etaMax = None + +## Minimal pseudorapidity for random distibution (overrides thetaMax) +SIM.gun.etaMin = None + +## isotropic distribution for the particle gun +## +## use the options phiMin, phiMax, thetaMin, and thetaMax to limit the range of randomly distributed directions +## if one of these options is not None the random distribution will be set to True and cannot be turned off! +## +SIM.gun.isotrop = False + +## Maximal momentum when using distribution (default = 0.0) +# SIM.gun.momentumMax = 10000.0 + +## Minimal momentum when using distribution (default = 0.0) +# SIM.gun.momentumMin = 0.0 +SIM.gun.multiplicity = 1 +SIM.gun.particle = "e-" + +## Maximal azimuthal angle for random distribution +SIM.gun.phiMax = None + +## Minimal azimuthal angle for random distribution +SIM.gun.phiMin = None + +## position of the particle gun, 3 vector +SIM.gun.position = (0.0, 0.0, 0.0) + +## Maximal polar angle for random distribution +SIM.gun.thetaMax = None + +## Minimal polar angle for random distribution +SIM.gun.thetaMin = None + + +################################################################################ +## Configuration for the hepmc3 InputFiles +################################################################################ + +## Set the name of the attribute contraining color flow information index 0. +SIM.hepmc3.Flow1 = "flow1" + +## Set the name of the attribute contraining color flow information index 1. +SIM.hepmc3.Flow2 = "flow2" + +## Set to false if the input should be opened with the hepmc2 ascii reader. +## +## If ``True`` a '.hepmc' file will be opened with the HEPMC3 Reader Factory. +## +## Defaults to true if DD4hep was build with HEPMC3 support. +## +SIM.hepmc3.useHepMC3 = True + + +################################################################################ +## Configuration for Input Files. +################################################################################ + +## Set one or more functions to configure input steps. +## +## The functions must take a ``DD4hepSimulation`` object as their only argument and return the created generatorAction +## ``gen`` (for example). +## +## For example one can add this to the ddsim steering file: +## +## def exampleUserPlugin(dd4hepSimulation): +## '''Example code for user created plugin. +## +## :param DD4hepSimulation dd4hepSimulation: The DD4hepSimulation instance, so all parameters can be accessed +## :return: GeneratorAction +## ''' +## from DDG4 import GeneratorAction, Kernel +## # Geant4InputAction is the type of plugin, Cry1 just an identifier +## gen = GeneratorAction(Kernel(), 'Geant4InputAction/Cry1' , True) +## # CRYEventReader is the actual plugin, steeringFile its constructor parameter +## gen.Input = 'CRYEventReader|' + 'steeringFile' +## # we can give a dictionary of Parameters that has to be interpreted by the setParameters function of the plugin +## gen.Parameters = {'DataFilePath': '/path/to/files/data'} +## gen.enableUI() +## return gen +## +## SIM.inputConfig.userInputPlugin = exampleUserPlugin +## +## Repeat function definition and assignment to add multiple input steps +## +## +SIM.inputConfig.userInputPlugin = [] + + +################################################################################ +## Configuration for the generator-level InputFiles +################################################################################ + +## Set the name of the collection containing the MCParticle input. +## Default is "MCParticle". +## +SIM.lcio.mcParticleCollectionName = "MCParticle" + + +################################################################################ +## Configuration for the LCIO output file settings +################################################################################ + +## The event number offset to write in slcio output file. E.g setting it to 42 will start counting events from 42 instead of 0 +SIM.meta.eventNumberOffset = 0 + +## Event parameters to write in every event. Use C/F/I ids to specify parameter type. E.g parameterName/F=0.42 to set a float parameter +SIM.meta.eventParameters = [] + +## The run number offset to write in slcio output file. E.g setting it to 42 will start counting runs from 42 instead of 0 +SIM.meta.runNumberOffset = 0 + + +################################################################################ +## Configuration for the output levels of DDG4 components +################################################################################ + +## Output level for geometry. +SIM.output.geometry = 2 + +## Output level for input sources +SIM.output.inputStage = 3 + +## Output level for Geant4 kernel +SIM.output.kernel = 3 + +## Output level for ParticleHandler +SIM.output.part = 3 + +## Output level for Random Number Generator setup +SIM.output.random = 6 + + +################################################################################ +## Configuration for Output Files. +################################################################################ + +## Use the DD4HEP output plugin regardless of outputfilename. +SIM.outputConfig.forceDD4HEP = False + +## Use the EDM4HEP output plugin regardless of outputfilename. +SIM.outputConfig.forceEDM4HEP = False + +## Use the LCIO output plugin regardless of outputfilename. +SIM.outputConfig.forceLCIO = False + +## Set a function to configure the outputFile. +## +## The function must take a ``DD4hepSimulation`` object as its only argument and return ``None``. +## +## For example one can add this to the ddsim steering file: +## +## def exampleUserPlugin(dd4hepSimulation): +## '''Example code for user created plugin. +## +## :param DD4hepSimulation dd4hepSimulation: The DD4hepSimulation instance, so all parameters can be accessed +## :return: None +## ''' +## from DDG4 import EventAction, Kernel +## dd = dd4hepSimulation # just shorter variable name +## evt_root = EventAction(Kernel(), 'Geant4Output2ROOT/' + dd.outputFile, True) +## evt_root.HandleMCTruth = True or False +## evt_root.Control = True +## output = dd.outputFile +## if not dd.outputFile.endswith(dd.outputConfig.myExtension): +## output = dd.outputFile + dd.outputConfig.myExtension +## evt_root.Output = output +## evt_root.enableUI() +## Kernel().eventAction().add(evt_root) +## return None +## +## SIM.outputConfig.userOutputPlugin = exampleUserPlugin +## # arbitrary options can be created and set via the steering file or command line +## SIM.outputConfig.myExtension = '.csv' +## + + +def Geant4Output2EDM4hep_DRC_plugin(dd4hepSimulation): + from DDG4 import EventAction, Kernel + + evt_root = EventAction( + Kernel(), "Geant4Output2EDM4hep_DRC/" + dd4hepSimulation.outputFile, True + ) + evt_root.Control = True + output = dd4hepSimulation.outputFile + evt_root.Output = output + evt_root.enableUI() + Kernel().eventAction().add(evt_root) + return None + + +SIM.outputConfig.userOutputPlugin = Geant4Output2EDM4hep_DRC_plugin + +################################################################################ +## Configuration for the Particle Handler/ MCTruth treatment +################################################################################ + +## Enable lots of printout on simulated hits and MC-truth information +SIM.part.enableDetailedHitsAndParticleInfo = False + +## Keep all created particles +SIM.part.keepAllParticles = False + +## Minimal distance between particle vertex and endpoint of parent after +## which the vertexIsNotEndpointOfParent flag is set +## +SIM.part.minDistToParentVertex = 2.2e-14 + +## MinimalKineticEnergy to store particles created in the tracking region +SIM.part.minimalKineticEnergy = 1.0 + +## Printout at End of Tracking +SIM.part.printEndTracking = False + +## Printout at Start of Tracking +SIM.part.printStartTracking = False + +## List of processes to save, on command line give as whitespace separated string in quotation marks +SIM.part.saveProcesses = ["Decay"] + +## Optionally enable an extended Particle Handler +SIM.part.userParticleHandler = "Geant4TCUserParticleHandler" + + +################################################################################ +## Configuration for the PhysicsList and Monte Carlo particle selection. +## +## To load arbitrary plugins, add a function to be executed. +## +## The function must take the DDG4.Kernel() object as the only argument. +## +## For example, add a function definition and the call to a steering file:: +## +## def setupCerenkov(kernel): +## from DDG4 import PhysicsList +## seq = kernel.physicsList() +## cerenkov = PhysicsList(kernel, 'Geant4CerenkovPhysics/CerenkovPhys') +## cerenkov.MaxNumPhotonsPerStep = 10 +## cerenkov.MaxBetaChangePerStep = 10.0 +## cerenkov.TrackSecondariesFirst = True +## cerenkov.VerboseLevel = 2 +## cerenkov.enableUI() +## seq.adopt(cerenkov) +## ph = PhysicsList(kernel, 'Geant4OpticalPhotonPhysics/OpticalGammaPhys') +## ph.addParticleConstructor('G4OpticalPhoton') +## ph.VerboseLevel = 2 +## ph.enableUI() +## seq.adopt(ph) +## return None +## +## SIM.physics.setupUserPhysics(setupCerenkov) +## +## # End of example +## +################################################################################ + +## Set of Generator Statuses that are used to mark unstable particles that should decay inside of Geant4. +## +SIM.physics.alternativeDecayStatuses = set() + +## If true, add decay processes for all particles. +## +## Only enable when creating a physics list not based on an existing Geant4 list! +## +SIM.physics.decays = False + +## The name of the Geant4 Physics list. +SIM.physics.list = "FTFP_BERT" + +## location of particle.tbl file containing extra particles and their lifetime information +## +## For example in $DD4HEP/examples/DDG4/examples/particle.tbl +## +SIM.physics.pdgfile = None + +## The global geant4 rangecut for secondary production +## +## Default is 0.7 mm as is the case in geant4 10 +## +## To disable this plugin and be absolutely sure to use the Geant4 default range cut use "None" +## +## Set printlevel to DEBUG to see a printout of all range cuts, +## but this only works if range cut is not "None" +## +SIM.physics.rangecut = None + +## Set of PDG IDs that will not be passed from the input record to Geant4. +## +## Quarks, gluons and W's Z's etc should not be treated by Geant4 +## +SIM.physics.rejectPDGs = { + 1, + 2, + 3, + 4, + 5, + 6, + 3201, + 3203, + 4101, + 4103, + 21, + 23, + 24, + 5401, + 25, + 2203, + 5403, + 3101, + 3103, + 4403, + 2101, + 5301, + 2103, + 5303, + 4301, + 1103, + 4303, + 5201, + 5203, + 3303, + 4201, + 4203, + 5101, + 5103, + 5503, +} + +## Set of PDG IDs for particles that should not be passed to Geant4 if their properTime is 0. +## +## The properTime of 0 indicates a documentation to add FSR to a lepton for example. +## +SIM.physics.zeroTimePDGs = {17, 11, 13, 15} + + +def setupOpticalPhysics(kernel): + from DDG4 import PhysicsList + + seq = kernel.physicsList() + cerenkov = PhysicsList(kernel, "Geant4CerenkovPhysics/CerenkovPhys") + cerenkov.TrackSecondariesFirst = True + cerenkov.VerboseLevel = 1 + cerenkov.enableUI() + seq.adopt(cerenkov) + + opt = PhysicsList(kernel, "Geant4OpticalPhotonPhysics/OpticalGammaPhys") + opt.addParticleConstructor("G4OpticalPhoton") + opt.VerboseLevel = 1 + # set BoundaryInvokeSD to true when using DRC wafer as the SD + # opt.BoundaryInvokeSD = True + opt.enableUI() + seq.adopt(opt) + + return None + +if simulateCalo: + SIM.physics.setupUserPhysics(setupOpticalPhysics) + +def setupDRCFastSim(kernel): + from DDG4 import DetectorConstruction, Geant4, PhysicsList + + geant4 = Geant4(kernel) + seq = geant4.detectorConstruction() + # Create a model for fast simulation + model = DetectorConstruction(kernel, str("Geant4DRCFiberModel/ShowerModel")) + # Mandatory model parameters + model.RegionName = "FastSimOpFiberRegion" + model.Enable = True + model.ApplicableParticles = ["opticalphoton"] + model.enableUI() + seq.adopt(model) + # Now build the physics list: + phys = kernel.physicsList() + ph = PhysicsList(kernel, str("Geant4FastPhysics/FastPhysicsList")) + ph.EnabledParticles = ["opticalphoton"] + ph.BeVerbose = True + ph.enableUI() + phys.adopt(ph) + phys.dump() + + +# turn-on fastsim if the skipScint option of the SD is set to false +# SIM.physics.setupUserPhysics(setupDRCFastSim) + +################################################################################ +## Properties for the random number generator +################################################################################ + +## If True, calculate random seed for each event basedon eventID and runID +## Allows reproducibility even whenSkippingEvents +SIM.random.enableEventSeed = False +SIM.random.file = None +SIM.random.luxury = 1 +SIM.random.replace_gRandom = True +SIM.random.seed = None +SIM.random.type = None + + +################################################################################ +## Configuration for setting commands to run during different phases. +## +## In this section, one can configure commands that should be run during the different phases of the Geant4 execution. +## +## 1. Configuration +## 2. Initialization +## 3. Pre Run +## 4. Post Run +## 5. Terminate / Finalization +## +## For example, one can add +## +## >>> SIM.ui.commandsConfigure = ['/physics_lists/em/SyncRadiation true'] +## +## Further details should be taken from the Geant4 documentation. +## +################################################################################ + +## List of UI commands to run during the 'Configure' phase. +SIM.ui.commandsConfigure = [] + +## List of UI commands to run during the 'Initialize' phase. +SIM.ui.commandsInitialize = [] + +## List of UI commands to run during the 'PostRun' phase. +SIM.ui.commandsPostRun = [] + +## List of UI commands to run during the 'PreRun' phase. +SIM.ui.commandsPreRun = [] + +## List of UI commands to run during the 'Terminate' phase. +SIM.ui.commandsTerminate = [] \ No newline at end of file diff --git a/VTXdigi_Modular/test/test.sh b/VTXdigi_Modular/test/test.sh new file mode 100644 index 00000000..4e2c1b33 --- /dev/null +++ b/VTXdigi_Modular/test/test.sh @@ -0,0 +1,102 @@ +#!/bin/bash +############################################################################### +# Test script +# +# Requirements: +# - key4hep stack must be sourced +# - Environment variables K4GEO and CLDCONFIG must be set +############################################################################### + +Dir_TestSrc="${SOURCE_DIR_TEST:-$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)}" + +export Dir_TestSrc +export PYTHONPATH="${Dir_TestSrc}:$PYTHONPATH" + +# Check required environment variables +if [ -z "${KEY4HEP_STACK:-}" ] || [ -z "${K4GEO:-}" ] || [ -z "${CLDCONFIG:-}" ]; then + echo "Key4hep stack not found in this environment. Setting up latest nightly..." + source /cvmfs/sw-nightlies.hsf.org/key4hep/setup.sh +fi + +# Configuration +NEvts=10 +Momentum=10 +Particle="mu-" +File_SimOutput="Step1_test_vtxdigi_edm4hep.root" +File_DigiOutput="Test_vtxdigi_output" +File_Log="VTXdigi_Modular.log" + +echo "========================================================================" +echo " VTXDigi_Modular execution test" +echo "========================================================================" +echo " Events: ${NEvts}" +echo " Particle: ${Particle} @ ${Momentum} GeV (isotropic particle gun)" +echo " IDEA detector, inner VTX" +echo "========================================================================" + +# Step 1: Simulation with ddsim +echo " [1/3] Running simulation with ddsim..." + +File_SimSteering="${Dir_TestSrc}/steering_sim.py" +File_Detector=${K4GEO}/FCCee/IDEA/compact/IDEA_o1_v03/IDEA_o1_v03.xml + +ddsim --steeringFile ${File_SimSteering} \ + --compactFile ${File_Detector} \ + --enableGun \ + --gun.distribution uniform \ + --gun.energy "${Momentum}*GeV" \ + --gun.particle ${Particle} \ + --numberOfEvents ${NEvts} \ + --outputFile ${File_SimOutput} \ + # > File_Log 2>&1 + +if [ $? -ne 0 ]; then + echo "ERROR: Simulation failed" + echo "Logs in ${File_Log}" + exit 1 +fi +echo " ✓ Simulation completed: ${File_SimOutput}" + +# Step 2: Create digitization steering file +echo "" +echo "[2/3] Running digitization with VTXdigitizerDetailed..." + +File_DigiSteering="${Dir_TestSrc}/steering_digi.py" + +k4run ${File_DigiSteering} \ +-n ${N_EVENTS} \ +--inputFiles ${File_SimOutput} \ +--outputBasename ${} + +if [ $? -ne 0 ]; then + echo "ERROR: Digitization failed" + echo "Logs in ${File_Log}" + exit 1 +fi +echo " ✓ Digitization completed: ${File_DigiOutput}_REC.edm4hep.root" + +# Step 3: Verify output +echo "" +echo "[3/3] Verifying output..." +if [ ! -f "${File_DigiOutput}_REC.edm4hep.root" ]; then + echo "ERROR: Output file not created" + echo "Logs in ${File_Log}" + exit 1 +fi + +# Check collections in output file +echo "" +echo "Test if collections are in output file:" +podio-dump ${File_DigiOutput}_REC.edm4hep.root -c -d 2>/dev/null | grep -E "(VertexBarrel|VertexEndcap|InnerTrackerBarrel|InnerTrackerEndcap|OuterTrackerBarrel|OuterTrackerEndcap)" || true + +# Cleanup +echo "" +echo "Cleaning up temporary files..." +rm -rf ${File_SimOutput} cdb.log ${File_DigiOutput}_aida.root __pycache__ + +echo "" +echo "========================================================================" +echo " ✓ Test PASSED" +echo " Output file: ${File_DigiOutput}_REC.edm4hep.root" +echo " Log file: ${File_Log}" +echo "========================================================================"