IDEA DCH full digitization chain by muhammadsahil5030 · Pull Request #78 · key4hep/k4RecTracker

muhammadsahil5030 · 2026-04-07T15:40:54Z

BEGINRELEASENOTES

Description

This PR introduces a full digitization chain for the IDEA Drift Chamber (DCH).

The implementation follows the standard detector signal formation sequence:
ionization → drift → amplification → signal formation → electronics → digitization.

New Gaudi algorithm:
- DCHdigi_v03
  - cluster generation based on Bethe–Bloch
  - drift time from 2D LUT (x–t relation)
  - Polya-distributed avalanche gain
  - waveform formation (single-electron pulse model)
  - signal propagation along wire (delay + attenuation)
  - electronics response (transfer function)
  - ADC digitization
Drift time handling:
- DCHXT2DLUT
  - reads xt_mean and xt_error from ROOT file
  - provides interpolated drift time + smearing
Noise simulation:
- DCHFFTNoise
  - FFT-based colored noise generation
  - compatible with waveform sampling
EDM4hep extensions:
- SenseWireHit
- waveform storage via RawTimeSeriesCollection
Steering:
- configurable via runDCHdigiV3.py

Outputs

SenseWireHitCollection
SenseWireHitSimTrackerHitLinkCollection
Digitized waveforms (left/right) as RawTimeSeriesCollection

ENDRELEASENOTES

jmcarcell · 2026-04-07T18:53:37Z

DCHdigi/test/test_DCHdigi/runDCHdigiV3.py

@@ -0,0 +1,93 @@
+#


There should be a test that runs this steering file, otherwise we don't know if it runs.

Updated the debug output file name to match the new naming convention.

tmadlener

This doesn't look like it has been clang-formatted yet (and CI will likely fail also due to that.

I have started to go through the code and there are a few things that could be improved. I didn't look at all of it, but some of the comments might apply on a broader scope.

tmadlener · 2026-04-08T08:48:49Z

DCHdigi/include/BetheBloch.h

+namespace BB
+{
+//Atomic number of material
+double Z=2;
+
+//Relative atomic mass
+double A=4.0026;
+
+//Atomic number of particle
+double z_particle=1;
+
+//Electron rest mass MeV
+double me = 0.511;
+
+//Electron charge in eV
+double e = 1;
+
+double K = 0.307075;
+
+//double I=41.8*1e-6;
+//variables for delta
+double c=11.1393;
+double x0=2.2017;
+double x1=3.6122;
+double a=0.13443;
+double k1=5.8347;
+const double k2=2.0*TMath::Log(10);
+
+//double getDelta(double p, double m){
+//	double delta;
+//	double x= TMath::Log10(p/m);
+//
+//	cout<<x<<"\t"<<p<<endl;
+//	if(x>=x1) delta=k2*x-c;
+//	if(x<x1&&x>=x0) delta=k2*x-c+a*pow((x1-x),k1);
+//	if(x<x0) delta=0.0;
+//


Are these used anywhere else than below in bethe_bloch? In that case I would simply move them there. Also these could / should probably marked constepr (where possible) or at least const

tmadlener · 2026-04-08T08:55:46Z

DCHdigi/include/BetheBloch.h

+//	return delta;
+//}
+
+double bethe_bloch( double *xp,double *par){


It doesn't look like this has to go through a ROOT interface at any point. Hence, we should be able to pass the xp and par as ,std::arrays. Checking the only call site of this, it might even just be double xp and std::array<double, 2> par. Alternatively already explicitly use m and I as arguments to not hide behind a way too generic interface.

tmadlener · 2026-04-08T08:57:08Z

DCHdigi/include/DCHdigi_v03.h

+	  extension::SenseWireHitCollection, 
+	  extension::SenseWireHitSimTrackerHitLinkCollection,


Should this still use the extension types, or rather the ones that went to EDM4hep?

tmadlener · 2026-04-08T08:57:30Z

DCHdigi/include/DCHdigi_v03.h

+class TF1;
+class TFile;
+
+/// constant to convert from mm (EDM4hep) to DD4hep (cm)


Stray comment?

tmadlener · 2026-04-08T09:01:33Z

DCHdigi/include/DCHdigi_v03.h

+  TVector3 Convert_EDM4hepVector_to_TVector3(const edm4hep::Vector3d& v, double scale) const {
+    return {v[0] * scale, v[1] * scale, v[2] * scale};
+  };
+  edm4hep::Vector3d Convert_TVector3_to_EDM4hepVector(const TVector3& v, double scale) const {
+    return {v.x() * scale, v.y() * scale, v.z() * scale};
+  };


At least for the scaling we have utilities in EDM4hep: https://github.com/key4hep/EDM4hep/blob/6d24c90b72c2434cacf3824a9b05b758d7e3911f/utils/include/edm4hep/utils/vector_utils.h#L293-L296

I.e. it should be possible to do:

return edm4hep::Vector3d{v.x(), v.y(), v.z()} * scale;

tmadlener · 2026-04-08T09:14:52Z

DCHdigi/include/DCHUtilities.h

+// This function is used to calculate the dNdx using BB.h file:
+// when we have a beam-test TGraph gNcldx(bg)->(clusters/cm),
+inline double DCHdigi_v03::get_dNcldx_per_cm(double betagamma, const edm4hep::MCParticle& mc) const
+{
+  // 1) reconstruct momentum from β
+  const double m_GeV = mc.getMass();
+  const double p_GeV = betagamma * m_GeV;
+
+  // 2) call your header: bethe_bloch(xp, par)
+  double xp[1]  = { p_GeV };                           // GeV (header multiplies by 1e3 internally)
+  double par[2] = { m_GeV * 1.0e3,                     // mass in MeV
+                    m_MeanExcEnergy_eV.value() * 1e-6  // I in MeV
+                  };
+  const double dEdx_MeVcm2_per_g = BB::bethe_bloch(xp, par);
+  const double dEdx_MeV_per_cm = dEdx_MeVcm2_per_g * m_GasDensity_g_cm3.value();
+  const double lambda_per_cm = (dEdx_MeV_per_cm * 1.0e6) / m_W_eff_eV.value();
+
+  return (lambda_per_cm > 0.0) ? lambda_per_cm : 0.0;
+}


This should probably live in the cpp file of the algorithm to which it belongs.

tmadlener · 2026-04-08T09:25:14Z

DCHdigi/include/DCHUtilities.h

+
+// --- Experimental cluster size probabilities w(n) for He-iSobutane (90:10) ---
+// Based on Fischle et al., NIM A301 (1991)
+inline const std::vector<double>& w_cluster_He_iC4H10()


This can be done as constexpr when using a std::array<double, 35>. That way we don't have to keep track of init and we don't need a static vector. Plus we can guarantee that all the (admittedly probably negligible) work is done at compile time.

See e.g. https://godbolt.org/z/GsPTsG39E

tmadlener · 2026-04-08T09:29:39Z

DCHdigi/include/DCHXT2DLUT.h

+class DCHXT2DLUT {
+public:
+  DCHXT2DLUT() = default;
+  ~DCHXT2DLUT();
+
+  // Load xt_mean and xt_error from ROOT file (par.root)
+  bool load(const std::string& rootFile,
+            const std::string& meanName  = "xt_mean",
+            const std::string& sigmaName = "xt_error");
+
+  bool isLoaded() const { return m_loaded; }


The way this is used you could get rid of load and in turn isLoaded. In the algorithm you construct it and immediately afterwards call load. You never call load again. So you could just make the constructor do the loading in which case you can guarantee that either you have a valid object or you throw an exception in case you don't.

m_file should probably be a unique_ptr in which case you could save some manual management.

tmadlener · 2026-04-08T09:30:38Z

DCHdigi/src/DCHdigi_v03.cpp

+  // advance internal state to minimize possibility of creating correlations
+  rng_engine.discard(10);
+  for (int i = 0; i < 10; ++i)
+    myRandom.Rndm();


You will probably have to do more than 10 iterations here if you really want to guarantee no correlations. I don't know exactly how many, but I seem to remember O(1000).

tmadlener · 2026-04-08T09:32:43Z

DCHdigi/src/DCHdigi_v03.cpp

+  if (!m_uidSvc)
+    ThrowException("Unable to get UniqueIDGenSvc");
+
+  if (0 > m_z_resolution.value())
+    ThrowException("Z resolution input value can not be negative!");
+
+  if (0 > m_xy_resolution.value())
+    ThrowException("Radial (XY) resolution input value can not be negative!");


Many of these can probably be done with K4_GAUDI_CHECK. Which would also make it possible to get rid of ThrowException (for which I don't see the added value, since throw std::runtime_error("...."); is only a small number of additional characters, but withotu hiding what you actually do.

muhammadsahil5030 added 8 commits April 7, 2026 16:51

Add files via upload

16e0f7a

Add files via upload

bbabcab

Add files via upload

5b9530b

Add files via upload

b2e98ec

Update runDCHdigiV3.py

8ee20b9

Update runDCHdigiV3.py

03c2b61

Merge branch 'key4hep:main' into main

babbdd9

Update runDCHdigiV3.py

6b0b249

jmcarcell reviewed Apr 7, 2026

View reviewed changes

muhammadsahil5030 added 3 commits April 7, 2026 23:00

Update execution command in runDCHdigiV3.py

b871cab

Update test script for DCH v3 simulation and digitization

1b45a16

Change debug output file name in check_DCHdigi_output.py

c61604c

Updated the debug output file name to match the new naming convention.

tmadlener reviewed Apr 8, 2026

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

IDEA DCH full digitization chain#78

IDEA DCH full digitization chain#78
muhammadsahil5030 wants to merge 11 commits intokey4hep:mainfrom
muhammadsahil5030:main

muhammadsahil5030 commented Apr 7, 2026

Uh oh!

jmcarcell Apr 7, 2026

Uh oh!

tmadlener left a comment

Uh oh!

tmadlener Apr 8, 2026

Uh oh!

tmadlener Apr 8, 2026

Uh oh!

tmadlener Apr 8, 2026

Uh oh!

tmadlener Apr 8, 2026

Uh oh!

tmadlener Apr 8, 2026

Uh oh!

tmadlener Apr 8, 2026

Uh oh!

tmadlener Apr 8, 2026

Uh oh!

tmadlener Apr 8, 2026

Uh oh!

tmadlener Apr 8, 2026

Uh oh!

tmadlener Apr 8, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

		extension::SenseWireHitCollection,
		extension::SenseWireHitSimTrackerHitLinkCollection,

Conversation

muhammadsahil5030 commented Apr 7, 2026

Description

Outputs

Uh oh!

Choose a reason for hiding this comment

Uh oh!

tmadlener left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants