Realistic BTOF digitization

[ssedd1123]

Briefly, what does this PR introduce?

It's an attempt to make LGAD sensor response for Barrel TOF more realistic by doing the following,

1. Spread charge deposition from strips that are hit directly to nearby strips.
2. Perform digitization by converting charge deposited to electric pulse.
3. Digitize the pulse by converting it to TDC and ADC value. ADC propto pulse height and TDC propto time when voltage crosses certain threshold.
4. Convert the TDC and ADC value back to time and Edep by linear fit.
5. Return the BTOF hit point as "TOFBarrelRecHit". Position is estimated by either weighted sum of ADC or just center of the cell with max ADC. It's weighted average by default, but the behavior can be changed by parameters.

Noise, edge correction and time talk correction will be developed in the future.

What kind of change does this PR introduce?

• [x ] Bug fix (issue Enhanced Realism in TOF Reconstruction & Digitization #1656 )
• New feature (issue #__)
• Documentation update
• [x ] Other: Implementation of digitization and clusterization for LGADs in Barrel TOF

Please check if this PR fulfills the following:

• Tests for the changes have been added
• Documentation has been added / updated
• [x ] Changes have been communicated to collaborators

Does this PR introduce breaking changes? What changes might users need to make to their code?

Does this PR change default behavior?

It replaces the content of "TOFBarrelRecHit" with results of this new clusterization. It was originally just geant point + random smearing. I have communicated to the reconstruction group and they advise that I submit this pull request so they can study the effect of having this new digitization/clusterization algorithm. They will decide if I should save the data in a new branch or replace the origin content.

[ssedd1123]

My takeaway from the meeting is that there were the following concerns raised:

1. The 25ns window starting at t=0 works currently for single event simulations where each event contains a single bunch crossing starting at t=0, but will not work for full frame simulations. The could would need to be adapted to work for that case, and we have some old input files to test against (ref. Simon's slide 2).
2. The current code assumes that both 100 MHz EIC and 40MHz HGCROC clocks are in phase=0 at t=0. This is too optimistic, as this case occurs only once in a 5 bunch crossings. This could be addressed by storing waveforms like it's done now, but also storing strobe counters for both clocks. Those are randomly sampled in the simulation part and affect points at which ADC counts are sampled. The reconstruction would have to invert that operation.
3. The current algorithms named TOFPulseGeneration and TOFPulseDigitization contain an implementation of HGCROC/EICROC simulator. If there is a piece of code that characterizes that device (namely the part that converts currents/voltages to ADC), it should be factored out into a common implementation for re-use in other detectors.

@ssedd1123 @simonge @johnlajoie @wdconinc Please comment to correct if I'm missing anything. We also need to identify what is a minimal viable product to ship this PR and what can go into further iterations.

My understanding is the same as yours. One more potentially problematic issue is that the PulseGeneration class now only passes a digitized version of the puluse. It would be ideal if I can pass the original pulse and let PulseDigitization class handle the digitization. I can pass Landau parameters, but then we are restricted to using Landau pulse only.

In my opinion, the optimum way is to pass the pulse as TF1 or just a functor. The voltage will be calculated from the TF1/functor in PulseDigitization when needed. However, that would require me to pass a generic TObject. @veprbl @wdconinc Do you know if it can be done? If not, passing a digitized pulse is not the end of the world either. We can interpolate between time bins as long as the bins are not too wide.

For now, I am just going to save an extra bit of information stating which 40 MHz cycle each pulse belongs to. I will keep the pulse from TOFPulseGenerator digitized the way it is. That is the simplest and cleanest way to extend the pulse classes to record all events in the same signal frame. Once I figure out a way to pass TF1 between analysis classes, I will try to revise these pulse classes further.

[simonge]

My takeaway from the meeting is that there were the following concerns raised:

1. The 25ns window starting at t=0 works currently for single event simulations where each event contains a single bunch crossing starting at t=0, but will not work for full frame simulations. The could would need to be adapted to work for that case, and we have some old input files to test against (ref. Simon's slide 2).
2. The current code assumes that both 100 MHz EIC and 40MHz HGCROC clocks are in phase=0 at t=0. This is too optimistic, as this case occurs only once in a 5 bunch crossings. This could be addressed by storing waveforms like it's done now, but also storing strobe counters for both clocks. Those are randomly sampled in the simulation part and affect points at which ADC counts are sampled. The reconstruction would have to invert that operation.
3. The current algorithms named TOFPulseGeneration and TOFPulseDigitization contain an implementation of HGCROC/EICROC simulator. If there is a piece of code that characterizes that device (namely the part that converts currents/voltages to ADC), it should be factored out into a common implementation for re-use in other detectors.

@ssedd1123 @simonge @johnlajoie @wdconinc Please comment to correct if I'm missing anything. We also need to identify what is a minimal viable product to ship this PR and what can go into further iterations.

We got a bit carried away into the details of the multiple clocks, synchronization, bunch crossing shapes and even the information provided by the TOF detector for the physics and vertex reconstruction. Fundamentally all I want to see from a minimum viable product is something which addresses your point 1. Ideally the pulse generation wouldn't care which 40MHz clock cycle it occurred in and that everything to do with clocks gets handled by the digitization.

I'm not sure your point 2 was discussed or maybe I missed the subtlety of the point until now. Ideally this would always be handled in the MC event sample but having an option to remove potential biases from standard samples would be useful too. If a pulse isn't restricted to the 25ns window and the time is just digitized to 20ps I don't think there should be any bias based on the strobe anyway. If the 140/120ps counters are introduced in the digitization code the precision of the measurement varies over the 25ns period and maybe the deadtime too.

In my opinion, the optimum way is to pass the pulse as TF1 or just a functor. The voltage will be calculated from the TF1/functor in PulseDigitization when needed. However, that would require me to pass a generic TObject. @veprbl @wdconinc Do you know if it can be done? If not, passing a digitized pulse is not the end of the world either. We can interpolate between time bins as long as the bins are not too wide.

I think as long as the binning of the digitized pulse is equal to or smaller than the precision of the signal, passing a digitized pulse is fine, the ASIC will be sampling the signal anyway. If the pulse shape itself moves away from a purely analytical function at some point to replicate some of the random processes/noise in the detector/sensor, having a workflow with digitized pulses makes more sense to me.

[veprbl]

My takeaway from the meeting is that there were the following concerns raised:

1. The 25ns window starting at t=0 works currently for single event simulations where each event contains a single bunch crossing starting at t=0, but will not work for full frame simulations. The could would need to be adapted to work for that case, and we have some old input files to test against (ref. Simon's slide 2).
2. The current code assumes that both 100 MHz EIC and 40MHz HGCROC clocks are in phase=0 at t=0. This is too optimistic, as this case occurs only once in a 5 bunch crossings. This could be addressed by storing waveforms like it's done now, but also storing strobe counters for both clocks. Those are randomly sampled in the simulation part and affect points at which ADC counts are sampled. The reconstruction would have to invert that operation.
3. The current algorithms named TOFPulseGeneration and TOFPulseDigitization contain an implementation of HGCROC/EICROC simulator. If there is a piece of code that characterizes that device (namely the part that converts currents/voltages to ADC), it should be factored out into a common implementation for re-use in other detectors.

@ssedd1123 @simonge @johnlajoie @wdconinc Please comment to correct if I'm missing anything. We also need to identify what is a minimal viable product to ship this PR and what can go into further iterations.

My understanding is the same as yours. One more potentially problematic issue is that the PulseGeneration class now only passes a digitized version of the puluse. It would be ideal if I can pass the original pulse and let PulseDigitization class handle the digitization. I can pass Landau parameters, but then we are restricted to using Landau pulse only.

In my opinion, the optimum way is to pass the pulse as TF1 or just a functor. The voltage will be calculated from the TF1/functor in PulseDigitization when needed. However, that would require me to pass a generic TObject. @veprbl @wdconinc Do you know if it can be done? If not, passing a digitized pulse is not the end of the world either. We can interpolate between time bins as long as the bins are not too wide.

For now, I am just going to save an extra bit of information stating which 40 MHz cycle each pulse belongs to. I will keep the pulse from TOFPulseGenerator digitized the way it is. That is the simplest and cleanest way to extend the pulse classes to record all events in the same signal frame. Once I figure out a way to pass TF1 between analysis classes, I will try to revise these pulse classes further.

I was thinking about the same topic. A significant portion of the discussion is about how the signals are sampled, and I get why you would want to pass an analytical signal shape into your digitization routine. However, while the framework technically allows to pass arbitrary C++ types between factories, we chose to only pass immutable PODIO structures. The types for those have to be defined in advance and we require a careful consideration and discussion for each item that we add to our data model https://github.com/eic/EDM4eic/blob/main/edm4eic.yaml. So it would be possible to make a type for a function (e.g. represented as a string + vector of floating point parameters). However, I doubt about universality of such approach, e.g. signals may come as a result of MC simulations and not have a nice analytical expression. Would it make sense to just reuse the waveform type and pass some extra oversampled values instead?

[ssedd1123]

@veprbl @simonge @wdconinc

I believe my pull request is now ready. Here's an overview of the changes I made:

1. I renamed all "TOFPulse" classes to "LGADPulse" to make them more general.
2. Dedicated configuration classes have been introduced for BTOFChargeSharing, LGADPulseGeneration, and LGADPulseDigitization to ensure consistency.
3. The time series generated by LGADPulseGeneration now includes clock cycle information. Specifically, pulse.getTime() returns the time (in ns) at which the cycle begins. For example, pulse.getTime() == 0 for hits arriving at 0 < t < 25 ns, pulse.getTime() == 25ns for hits arriving at 25 < t < 50 ns, and so on.
4. LGADPulseDigitization has been updated to handle this additional information. Hits arriving at 0 < t < 25 ns return TDC values between 0–1023, while hits arriving at 25 < t < 50 ns return TDC values between 1024–2047, and so on.
   Both LGADPulseGeneration and LGADPulseDigitization now support handling negative time.

I believe these modifications are viable, but there are some edge cases that I should let you know. If a pulse is generated near the edges of bin cycles, and the voltage exceeds the threshold to trigger a TDC in the first EICROC cycle but peaks in the second EICROC cycle, the code will return the correct TDC value but an incorrect ADC value because the peak is being cut off. This leads to information loss. However, I don't see this as an issue because it mirrors a real-world problem that the hardware designers need to deal with. Since the electronics are still in the design phase, there is no definitive answer on how they will behave in this scenario. I propose leaving this issue open and revisiting it once we now how the electronics will behave in the real world.

Another known issue is the absence of dead time in the current implementation. Some of my colleagues are trying to find out how much dead time there will be for each cell, but, as mentioned above, the specifics are uncertain and subject to change. Dead time can be added to the implementation later once the behavior of the electronics is better understood.

[ruse-traveler]

Very nice! I'm quite happy with how this has evolved, and I think the pulse generation-digitization workflow is in great shape! I have a few minor comments that I've left on those files.

But something that I think warrants a little more discussion is the charge sharing: I think it's actually pretty close to being generalized to work for more than just the BTOF, and have included some links to examples of how you might go about this in a comment.

That being said, I'd be willing to defer this change to a follow-up PR so long as we have a clear plan (i.e. an issue we can track and a point-of-contact for who's going to be working on it).

[ruse-traveler]

Like @simonge suggested in #1635 (comment), I think this might be where you could generalize this algorithm to work for more than just the BTOF. Correct me if I'm wrong, but at the end of the day the check you're doing here is (schematically):

bool isNeighbor = (
  (sensorID_hit == sensorID_test) &&
  (
    (std::abs(xID_hit - xID_test) == 1) ||
    (std::abs(yID_hit - yID_test) == 1)
  )
);
if (isNeighbor) answer.push_back( CellID_test );

If you could rewrite the check on CellID here in terms of a function like above, then you can build on some already-developed tools to generalize this to work for any detector. The difference between geometries would be encoded in the expression the user sets in the plugin (or at the commandline).

For example, the CalorimeterIslandCluster::is_neighbor method uses the EvaluatorSvc to evaluate a provided expression (set by CalorimeterIslandClusterConfig::adjacencyMatrix) that returns true or false for if 2 calorimeter hits are neighboring based on their Cell IDs. You can see how the expression gets turned into a std::function here, and then some examples of the expressions provided can be seen in almost all of the calorimeter plugins (the EEEMCal is a very straightforward example).

[ssedd1123]

Thank you for your suggestions. I am impressed. I did not know you can compile string in runtime. Is that an API to a JIT compiler?

Regardless of how it works, I will try it out and let you know of the results.