TPC: add 83mKr calibration generator with Geant4 ionisation-fluctuation support

- TPCDetParam.UseGeant4Edep: use Geant4's own energy deposit for ionisation
  instead of Bethe-Bloch/NA49, with a configurable SpecialCutsGeV threshold
- Adds the Kr-83m decay generator (GeneratorKrDecay) as a compiled Generator
  subclass in Detectors/TPC/simulation, with a thin krGenerator.C macro for
  use as an o2-sim external generator
- Adds plotCluster.C for viewing the resulting calibration spectrum

Author: Ankur Yadav

Detectors/TPC/base/include/TPCBase/ParameterDetector.h

@@ -33,6 +33,8 @@ struct ParameterDetector : public o2::conf::ConfigurableParamHelper<ParameterDet
   TimeBin TmaxTriggered = 550;   ///< Maximum time bin in case of triggered readout mode
   float DriftTimeOffset = 7.3;   ///< drift time offset in time bins (we observe ~2.4\mus before October 2023 and ~1.45 \mus after)
   bool ExcludeFCGap = true;      ///< exclude electrons created in the gap between the IFC vessel and OFC vessel and FC strips
+  bool UseGeant4Edep = false;    ///< use Geant4 energy deposit directly for ionisation (for Kr-83m calibration runs)
+  float SpecialCutsGeV = 1e-6f;  ///< CUTELE/CUTGAM/DCUTE/BCUTE applied to TPC drift-gas media by SetSpecialPhysicsCuts() when UseGeant4Edep is set (Kr-83m calibration isolation test); default 1e-6 GeV (1 keV) matches the prior hardcoded value
 
   O2ParamDef(ParameterDetector, "TPCDetParam");
 };

Detectors/TPC/simulation/CMakeLists.txt

@@ -18,11 +18,13 @@ o2_add_library(TPCSimulation
                        src/DigitTime.cxx
                        src/ElectronTransport.cxx
                        src/GEMAmplification.cxx
+                       src/GeneratorKrDecay.cxx
                        src/Point.cxx
                        src/SAMPAProcessing.cxx
                        src/IDCSim.cxx
                PUBLIC_LINK_LIBRARIES O2::DetectorsBase O2::SimulationDataFormat
                                      O2::TPCBase O2::TPCSpaceCharge O2::TPCSpaceChargeIO O2::TPCCalibration
+                                     O2::Generators
                                      ROOT::Physics)
 
 o2_target_root_dictionary(TPCSimulation
@@ -34,6 +36,7 @@ o2_target_root_dictionary(TPCSimulation
                                   include/TPCSimulation/DigitTime.h
                                   include/TPCSimulation/ElectronTransport.h
                                   include/TPCSimulation/GEMAmplification.h
+                                  include/TPCSimulation/GeneratorKrDecay.h
                                   include/TPCSimulation/Point.h
                                   include/TPCSimulation/SAMPAProcessing.h
                                   include/TPCSimulation/IDCSim.h)
@@ -54,9 +57,18 @@ if(BUILD_TESTING)
                                                O2::DataFormatsTPC
                          LABELS tpc)
 
+  o2_add_test_root_macro(macro/krGenerator.C
+                         PUBLIC_LINK_LIBRARIES O2::TPCSimulation
+                         LABELS tpc)
+
   o2_add_test_root_macro(macro/toyCluster.C
                          PUBLIC_LINK_LIBRARIES O2::TPCBase
                                                O2::DataFormatsTPC
                          LABELS tpc)
 
+  o2_add_test_root_macro(macro/plotCluster.C
+                         PUBLIC_LINK_LIBRARIES O2::TPCReconstruction
+                                               O2::DataFormatsTPC
+                         LABELS tpc)
+
 endif()

Detectors/TPC/simulation/README.md

@@ -12,7 +12,7 @@ For the digitization as conducted by the [TPC digitizer](include/TPCSimulation/D
 
 * The energy loss of each individual GEANT hit is converted into a number of electrons by dividing by the effective ionization potential W_i. Each of these electrons is in the following treated individually.
 * The electron is projected onto the readout plane, taking into account diffusion, i.e. smearing its position by a 3D gaussian function ([ElectronTransport](include/TPCSimulation/ElectronTransport.h)). Then, the position is transformed into the local coordinate system of the Readout Chamber (ROC).
-* Having arrived at the amplification stage, the electrons undergo amplification in the GEM stack ([GEMAmplification](include/TPCSimulation/GEMAmplification.h)), taking into account fluctuations of the gain. These fluctuations follow a Polya distribution. For performance considerations, two different versions of the amplification are available (one effective single-stage amplification and a successive simulation of the collection, amplification, and extraction processes in the individual GEMs.
+* Having arrived at the amplification stage, the electrons undergo amplification in the GEM stack ([GEMAmplification](include/TPCSimulation/GEMAmplification.h)), taking into account fluctuations of the gain. These fluctuations follow a Polya distribution. For performance considerations, two different versions of the amplification are available (one effective single-stage amplification and a successive simulation of the collection, amplification, and extraction processes in the individual GEMs).
 * Capacitive coupling of the amplification structure to the readout anode leads to another contribution to the signal, the so-called Common Mode effect. Since the bottom electrode of GEM 4 is unsegmented, capacitive coupling occurs within a full ROC ([CommonMode](include/TPCSimulation/CommonMode.h))
 * The charge signal is then folded with the transfer function of the front-end cards ([SAMPAProcessing](include/TPCSimulation/SAMPAProcessing.h)) and written to the intermediate storage container structure ([DigitContainer](include/TPCSimulation/DigitContainer.h)/[DigitTime](include/TPCSimulation/DigitTime.h)/[DigitGlobalPad](include/TPCSimulation/DigitGlobalPad.h)), which is described below.
 
@@ -24,7 +24,7 @@ The digitization is conducted for each TPC sector individually in order to ensur
 The input can be created by running the simulation `o2-sim`, which produces the file `o2sim.root` with the hits stored in separated branches for all sectors.
 It should be noted that due to diffusion and the space-charge distortions, charge leakage between sectors can occur. In order to avoid the unnecessary processing of individual hits, several measures are taken
 
-* for a given sector the hits within an additional safety margin of +/- 10 degree are processed. For this reason, For this reason, the hits are not stored for a given sector, but shifted by 10 degrees. Hence only two branches need to be loaded for the digitization of a given sector.
+* for a given sector the hits within an additional safety margin of +/- 10 degree are processed. For this reason, the hits are not stored for a given sector, but shifted by 10 degrees. Hence only two branches need to be loaded for the digitization of a given sector.
 * Individual hits are only processed when they are within the processed by 3 sigma of the expected width from diffusion
 
 Hits passing that requirement are further processed by the [TPC digitizer](include/TPCSimulation/Digitizer.h) and undergo the above described physics processes.
@@ -36,7 +36,48 @@ The [DigitContainer](include/TPCSimulation/DigitContainer.h) is a circular buffe
 The [DigitTime](include/TPCSimulation/DigitTime.h) is then a flat contained of [DigitGlobalPad](include/TPCSimulation/DigitGlobalPad.h), where the latter correspond to one pad on the pad plane. Accordingly, the buffering of the actual ADC values is conducted using this object.
 Similarly, the MC labels are passed throughout the chain, and finally sorted by the number of occurrences, i.e. the track with the largest contribution to the digit is mentioned first etc.
 
-Correlations among digits from different events can only occur within the integration time of the detector (plus additional 50% contigiency), and therefore the digits are written to disk when the processed event is more than 750 time bins after. This means, that saturation effects are applied to the ADC values stored in the [DigitGlobalPad](include/TPCSimulation/DigitGlobalPad.h) and the relevant information is transformed in a [Digit](../../../DataFormats/Detectors/TPC/include/DataFormatsTPC/Digit.h) which is written to disk.
+Correlations among digits from different events can only occur within the integration time of the detector (plus additional 50% contingency), and therefore the digits are written to disk when the processed event is more than 750 time bins after. This means, that saturation effects are applied to the ADC values stored in the [DigitGlobalPad](include/TPCSimulation/DigitGlobalPad.h) and the relevant information is transformed in a [Digit](../../../DataFormats/Detectors/TPC/include/DataFormatsTPC/Digit.h) which is written to disk.
 
 ### Output data
-The digitizer workflow produces the file `tpcdigits.root` by default, data is stored in separated branches for all sectors.
+The digitizer workflow produces the file `tpcdigits.root` by default, data is stored in separated branches for all sectors.
+
+# 83mKr calibration generator
+
+[GeneratorKrDecay](include/TPCSimulation/GeneratorKrDecay.h) is a `Generator` producing
+83mKr decay vertices uniformly distributed in the TPC drift volume, for gain-map and
+energy-resolution calibration simulation with `o2-sim`.
+
+Kr-83m decays via two sequential internal transitions. For each generated vertex, one of
+eight decay channels (combinations of the two transitions' internal-conversion/gamma
+modes) is sampled according to branching fractions derived from the transition energies
+and internal conversion coefficients (ICC), then the corresponding conversion
+electrons/Auger electrons/fluorescence photons ([KrDecayTable](include/TPCSimulation/GeneratorKrDecay.h))
+are emitted as primary tracks.
+
+Transition energies and ICC values are read at runtime from the installed Geant4
+photon-evaporation data file, `$G4LEVELGAMMADATA/z36.a83` (the environment variable is
+set automatically by Geant4's own setup, inherited by any O2 session). If the variable
+is unset or the file cannot be parsed, hardcoded fallback values from
+`PhotonEvaporation5.7/z36.a83` are used instead and a warning is printed.
+
+The number of decays generated per event defaults to 1000 and can be overridden via the
+`KR_N_PER_EVENT` environment variable.
+
+## Usage
+
+`GeneratorKrDecay` is used as an `o2-sim` external generator, via the thin macro
+[macro/krGenerator.C](macro/krGenerator.C). `fileName` must resolve to that macro; `$O2path`
+below is the path to your O2 source checkout:
+
+```shell
+export KR_N_PER_EVENT=5000  # decays per event; defaults to 1000 if unset
+o2-sim -g external -m TPC -n 100 \
+  --configKeyValues "GeneratorExternal.fileName=$O2path/Detectors/TPC/simulation/macro/krGenerator.C;GeneratorExternal.funcName=krGenerator();TPCDetParam.UseGeant4Edep=1"
+```
+
+To change Geant4-related parameters (e.g. `StepFunction`, fluorescence, Auger cascade),
+add `G4.configMacroFile=<path to a Geant4 UI macro>` to `--configKeyValues`.
+
+`TPCDetParam.UseGeant4Edep=1` is required: it switches `Detector::ProcessHits()` to use Geant4's
+own energy deposit directly instead of the default Bethe-Bloch/NA49 sampling, which is what
+correctly resolves the Kr-83m decay channels into their discrete energy peaks.

Detectors/TPC/simulation/include/TPCSimulation/Detector.h

@@ -146,6 +146,8 @@ class Detector : public o2::base::DetImpl<Detector>
   void PostTrack() override { ; }
   void PreTrack() override { ; }
 
+  void SetSpecialPhysicsCuts() override;
+
   void SetGeoFileName(const TString file) { mGeoFileName = file; }
   const TString& GetGeoFileName() const { return mGeoFileName; }
 

Detectors/TPC/simulation/include/TPCSimulation/GeneratorKrDecay.h

@@ -0,0 +1,105 @@
+// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
+// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+// All rights not expressly granted are reserved.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+
+/// \file GeneratorKrDecay.h
+/// \brief Generator for 83mKr decays, for TPC gain-map calibration simulation
+/// \author Ankur Yadav <ankur.yadav@cern.ch>
+
+#ifndef ALICEO2_TPC_GeneratorKrDecay_H_
+#define ALICEO2_TPC_GeneratorKrDecay_H_
+
+#include "Generators/Generator.h"
+#include <array>
+#include <vector>
+
+namespace o2
+{
+namespace tpc
+{
+
+// O2 status encoding from MCGenProperties.h
+// bits 0-8: hepmc(9), bits 9-18: gen(10), bits 19-28: reserved(10), bits 29-31: sentinel=5
+inline int krO2EncodedStatus(int hepmc, int gen = 0)
+{
+  return (5 << 29) | ((gen & 0x3FF) << 9) | (hepmc & 0x1FF);
+}
+
+namespace KrGenConfig
+{
+inline double rInner = 83.5;
+inline double rOuter = 246.5; // TPC outermost pad row outer edge ~247 cm; stay inside
+inline double halfZ = 249.7;
+inline int nPerEvent = 1000; // override at runtime via KR_N_PER_EVENT env var
+} // namespace KrGenConfig
+
+struct KrProduct {
+  int pdg;
+  double eKin;
+};
+
+/// Table of 83mKr internal-conversion/gamma decay channels and their
+/// branching fractions, derived at runtime from Geant4's level/gamma data
+/// (falls back to hardcoded PhotonEvaporation5.7/z36.a83 values if
+/// unavailable). Used by GeneratorKrDecay to sample one decay channel per
+/// generated vertex.
+class KrDecayTable
+{
+ public:
+  struct Channel {
+    double fraction;
+    int nProducts;
+    KrProduct products[6]; // max 5 used; 6 for safety
+  };
+  // Eight physically motivated channels (T1 mode x T2 mode):
+  //   T1: ICC_total=2035 -> 99.951% IC (75.163% outer-shell, 24.788% K-shell), 0.049% gamma
+  //       K-shell: 65.2% K-fluorescence (Kalpha), 34.8% K-Auger
+  //   T2: ICC_total=17.09 -> 94.472% IC, 5.528% gamma
+  //   Source: G4 PhotonEvaporation5.7/z36.a83, RadioactiveDecay5.6/z36.a83
+  static const int kNChannels = 8;
+  Channel channels[kNChannels];
+  double cumulative[kNChannels];
+  KrDecayTable();
+  const Channel& sample() const;
+};
+
+} // namespace tpc
+} // namespace o2
+
+namespace o2
+{
+namespace eventgen
+{
+
+/// FairGenerator producing 83mKr decay vertices uniformly distributed in the
+/// TPC drift volume, for gain-map/energy-resolution calibration simulation.
+/// Each vertex emits the conversion electrons/photons of one randomly
+/// sampled o2::tpc::KrDecayTable::Channel.
+class GeneratorKrDecay : public Generator
+{
+ public:
+  GeneratorKrDecay();
+  ~GeneratorKrDecay() override;
+  Bool_t Init() override;
+  Bool_t generateEvent() override;
+  Bool_t importParticles() override;
+
+ private:
+  o2::tpc::KrDecayTable* mTable;
+  std::vector<std::array<double, 3>> mVertices;
+  // No ClassDefOverride: the base Generator class's dictionary is sufficient
+  // for a runtime-only generator that is never streamed via ROOT I/O
+  // (matches o2::eventgen::GeneratorGeantinos and other Generator subclasses).
+};
+
+} // namespace eventgen
+} // namespace o2
+
+#endif // ALICEO2_TPC_GeneratorKrDecay_H_

Detectors/TPC/simulation/macro/krGenerator.C

@@ -0,0 +1,27 @@
+// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
+// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+// All rights not expressly granted are reserved.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+
+/// \file krGenerator
+/// \brief This macro instantiates the compiled 83mKr TPC calibration
+///        generator (o2::eventgen::GeneratorKrDecay), for use with
+///        o2-sim -g external --extGenFile krGenerator.C --extGenFunc krGenerator
+/// \author Ankur Yadav <ankur.yadav@cern.ch>
+
+#if !defined(__CLING__) || defined(__ROOTCLING__)
+#include "FairGenerator.h"
+#include "TPCSimulation/GeneratorKrDecay.h"
+#endif
+
+FairGenerator* krGenerator()
+{
+  auto gen = new o2::eventgen::GeneratorKrDecay();
+  return gen;
+}