|  | 
| 21 | 21 | #include "PulseDigi.h" | 
| 22 | 22 | 
 | 
| 23 | 23 | class HGCROCRawSample { | 
| 24 |  | -	public: | 
| 25 |  | -		struct RawEntry { | 
| 26 |  | -			double adc{0}; | 
| 27 |  | -			double toa{0}; | 
| 28 |  | -			double tot{0}; | 
| 29 |  | -			bool totProgress{false}; | 
| 30 |  | -			bool totComplete{false}; | 
| 31 |  | -		}; | 
| 32 |  | - | 
| 33 |  | -		HGCROCRawSample(std::size_t n_samp) { rawEntries.resize(n_samp); } | 
| 34 |  | - | 
| 35 |  | -		void setAmplitude(std::size_t idx, double amp) { rawEntries[idx].adc = amp; } | 
| 36 |  | -		void setTOA(std::size_t idx, double toa) { rawEntries[idx].toa = toa; } | 
| 37 |  | -		void setTOT(std::size_t idx, double cross_t) { rawEntries[idx].tot = cross_t - rawEntries[idx].toa; } | 
| 38 |  | -		void setTotProgress(std::size_t idx) { rawEntries[idx].totProgress = true; } | 
| 39 |  | -		void setTotComplete(std::size_t idx) { rawEntries[idx].totComplete = true; } | 
| 40 |  | - | 
| 41 |  | -		const std::vector<RawEntry>& getEntries() const { return rawEntries; } | 
| 42 |  | - | 
| 43 |  | -	private: | 
| 44 |  | -		std::vector<RawEntry> rawEntries; | 
|  | 24 | +public: | 
|  | 25 | +  struct RawEntry { | 
|  | 26 | +    double adc{0}; | 
|  | 27 | +    double toa{0}; | 
|  | 28 | +    double tot{0}; | 
|  | 29 | +    bool totProgress{false}; | 
|  | 30 | +    bool totComplete{false}; | 
|  | 31 | +  }; | 
|  | 32 | + | 
|  | 33 | +  HGCROCRawSample(std::size_t n_samp) { rawEntries.resize(n_samp); } | 
|  | 34 | + | 
|  | 35 | +  void setAmplitude(std::size_t idx, double amp) { rawEntries[idx].adc = amp; } | 
|  | 36 | +  void setTOA(std::size_t idx, double toa) { rawEntries[idx].toa = toa; } | 
|  | 37 | +  void setTOT(std::size_t idx, double cross_t) { | 
|  | 38 | +    rawEntries[idx].tot = cross_t - rawEntries[idx].toa; | 
|  | 39 | +  } | 
|  | 40 | +  void setTotProgress(std::size_t idx) { rawEntries[idx].totProgress = true; } | 
|  | 41 | +  void setTotComplete(std::size_t idx) { rawEntries[idx].totComplete = true; } | 
|  | 42 | + | 
|  | 43 | +  const std::vector<RawEntry>& getEntries() const { return rawEntries; } | 
|  | 44 | + | 
|  | 45 | +private: | 
|  | 46 | +  std::vector<RawEntry> rawEntries; | 
| 45 | 47 | }; | 
| 46 | 48 | 
 | 
| 47 | 49 | namespace eicrecon { | 
| 48 | 50 | 
 | 
| 49 |  | -	void PulseDigi::init() {} | 
| 50 |  | - | 
| 51 |  | -	void PulseDigi::process(const PulseDigi::Input& input, const PulseDigi::Output& output) const { | 
| 52 |  | -		const auto [in_pulses] = input; | 
| 53 |  | -		auto [out_digi_hits]   = output; | 
| 54 |  | - | 
| 55 |  | -		for (const auto& pulse : *in_pulses) { | 
| 56 |  | -			double pulse_t     = pulse.getTime(); | 
| 57 |  | -			double pulse_dt    = pulse.getInterval(); | 
| 58 |  | -			std::size_t n_amps = pulse.getAmplitude().size(); | 
| 59 |  | - | 
| 60 |  | -			// Estimate the number of samples. | 
| 61 |  | -			const std::size_t timeIdx_begin = | 
| 62 |  | -				static_cast<std::size_t>(std::floor(pulse_t / m_cfg.time_window)); | 
| 63 |  | -			const std::size_t timeIdx_end = static_cast<std::size_t>( | 
| 64 |  | -					std::floor((pulse_t + (n_amps - 1) * pulse_dt) / m_cfg.time_window)); | 
| 65 |  | - | 
| 66 |  | -			HGCROCRawSample raw_sample(timeIdx_end - timeIdx_begin + 1); | 
| 67 |  | -	 | 
| 68 |  | -			// For ADC, amplitude is measured with a fixed phase. | 
| 69 |  | -			// This was reproduced by sample_tick and adc_counter as follows. | 
| 70 |  | -			// Amplitude is measured whenever adc_counter reaches sample_tick. | 
| 71 |  | -			int sample_tick = std::llround(m_cfg.time_window / pulse_dt); | 
| 72 |  | -			int adc_counter = | 
| 73 |  | -				std::llround((timeIdx_begin * m_cfg.time_window + m_cfg.adc_phase - pulse_t) / pulse_dt); | 
| 74 |  | - | 
| 75 |  | -			bool tot_progress  = false; | 
| 76 |  | -			bool tot_complete  = false; | 
| 77 |  | -			std::size_t toaIdx = 0; | 
| 78 |  | - | 
| 79 |  | -			for (std::size_t i = 0; i < n_amps; i++) { | 
| 80 |  | -				double t = pulse_t + i * pulse_dt; | 
| 81 |  | -				const std::size_t sampleIdx = | 
| 82 |  | -					static_cast<std::size_t>(std::floor(t / m_cfg.time_window)) - timeIdx_begin; | 
| 83 |  | - | 
| 84 |  | -				adc_counter++; | 
| 85 |  | - | 
| 86 |  | -				// Measure amplitudes for ADC | 
| 87 |  | -				if (adc_counter == sample_tick) { | 
| 88 |  | -					raw_sample.setAmplitude(sampleIdx, pulse.getAmplitude()[i]); | 
| 89 |  | -					adc_counter = 0; | 
| 90 |  | -					if(tot_progress) raw_sample.setTotProgress(sampleIdx); | 
| 91 |  | -				} | 
| 92 |  | - | 
| 93 |  | -				// Measure up-crossing time for TOA | 
| 94 |  | -				if (!tot_progress && pulse.getAmplitude()[i] > m_cfg.toa_thres) { | 
| 95 |  | -					toaIdx = sampleIdx; | 
| 96 |  | -					raw_sample.setTOA(sampleIdx, | 
| 97 |  | -							get_crossing_time(m_cfg.toa_thres, pulse_dt, t, pulse.getAmplitude()[i], | 
| 98 |  | -								pulse.getAmplitude()[i - 1])); | 
| 99 |  | -					tot_progress = true; | 
| 100 |  | -					tot_complete = false; | 
| 101 |  | -					raw_sample.setTotProgress(sampleIdx); | 
| 102 |  | -				} | 
| 103 |  | - | 
| 104 |  | -				// Measure down-crossing time for TOT | 
| 105 |  | -				if (tot_progress && !tot_complete && pulse.getAmplitude()[i] < m_cfg.tot_thres) {  | 
| 106 |  | -					raw_sample.setTOT(toaIdx, | 
| 107 |  | -							get_crossing_time(m_cfg.tot_thres, pulse_dt, t, pulse.getAmplitude()[i], | 
| 108 |  | -								pulse.getAmplitude()[i - 1])); | 
| 109 |  | -					tot_progress = false; | 
| 110 |  | -					tot_complete = true; | 
| 111 |  | -					raw_sample.setTotComplete(sampleIdx); | 
| 112 |  | -				} | 
| 113 |  | -			} | 
| 114 |  | - | 
| 115 |  | -			// Fill HGCROCSamples and RawHGCROCHit | 
| 116 |  | -			auto out_digi_hit = out_digi_hits->create(); | 
| 117 |  | -			out_digi_hit.setCellID(pulse.getCellID()); | 
| 118 |  | - | 
| 119 |  | -			const auto& entries = raw_sample.getEntries(); | 
| 120 |  | - | 
| 121 |  | -			for (const auto& entry : entries) { | 
| 122 |  | -				edm4eic::HGCROCSample sample;	 | 
| 123 |  | -				auto adc = std::max(std::llround(entry.adc / m_cfg.dyRangeADC * m_cfg.capADC), 0LL); | 
| 124 |  | -				sample.ADC = adc > m_cfg.capADC ? m_cfg.capADC : adc; | 
| 125 |  | -				auto toa = std::max(std::llround(entry.toa / m_cfg.dyRangeTOA * m_cfg.capTOA), 0LL); | 
| 126 |  | -                                sample.timeOfArrival = toa > m_cfg.capTOA ? m_cfg.capTOA : toa; | 
| 127 |  | -				auto tot = std::max(std::llround(entry.tot / m_cfg.dyRangeTOT * m_cfg.capTOT), 0LL); | 
| 128 |  | -                                sample.timeOverThreshold = tot > m_cfg.capTOT ? m_cfg.capTOT : tot; | 
| 129 |  | -				sample.TOTInProgress = entry.totProgress; | 
| 130 |  | -				sample.TOTComplete = entry.totComplete; | 
| 131 |  | -				out_digi_hit.addToSamples(sample); | 
| 132 |  | -			} | 
| 133 |  | -		} | 
| 134 |  | -	} // PulseDigi:process | 
| 135 |  | - | 
| 136 |  | -	double PulseDigi::get_crossing_time(double thres, double dt, double t, double amp1, | 
| 137 |  | -			double amp2) const { | 
| 138 |  | -		double numerator   = (thres - amp2) * dt; | 
| 139 |  | -		double denominator = amp2 - amp1; | 
| 140 |  | -		double added       = t; | 
| 141 |  | -		return (numerator / denominator) + added; | 
| 142 |  | -	} | 
|  | 51 | +void PulseDigi::init() {} | 
|  | 52 | + | 
|  | 53 | +void PulseDigi::process(const PulseDigi::Input& input, const PulseDigi::Output& output) const { | 
|  | 54 | +  const auto [in_pulses] = input; | 
|  | 55 | +  auto [out_digi_hits]   = output; | 
|  | 56 | + | 
|  | 57 | +  for (const auto& pulse : *in_pulses) { | 
|  | 58 | +    double pulse_t     = pulse.getTime(); | 
|  | 59 | +    double pulse_dt    = pulse.getInterval(); | 
|  | 60 | +    std::size_t n_amps = pulse.getAmplitude().size(); | 
|  | 61 | + | 
|  | 62 | +    // Estimate the number of samples. | 
|  | 63 | +    const std::size_t timeIdx_begin = | 
|  | 64 | +        static_cast<std::size_t>(std::floor(pulse_t / m_cfg.time_window)); | 
|  | 65 | +    const std::size_t timeIdx_end = static_cast<std::size_t>( | 
|  | 66 | +        std::floor((pulse_t + (n_amps - 1) * pulse_dt) / m_cfg.time_window)); | 
|  | 67 | + | 
|  | 68 | +    HGCROCRawSample raw_sample(timeIdx_end - timeIdx_begin + 1); | 
|  | 69 | + | 
|  | 70 | +    // For ADC, amplitude is measured with a fixed phase. | 
|  | 71 | +    // This was reproduced by sample_tick and adc_counter as follows. | 
|  | 72 | +    // Amplitude is measured whenever adc_counter reaches sample_tick. | 
|  | 73 | +    int sample_tick = std::llround(m_cfg.time_window / pulse_dt); | 
|  | 74 | +    int adc_counter = | 
|  | 75 | +        std::llround((timeIdx_begin * m_cfg.time_window + m_cfg.adc_phase - pulse_t) / pulse_dt); | 
|  | 76 | + | 
|  | 77 | +    bool tot_progress  = false; | 
|  | 78 | +    bool tot_complete  = false; | 
|  | 79 | +    std::size_t toaIdx = 0; | 
|  | 80 | + | 
|  | 81 | +    for (std::size_t i = 0; i < n_amps; i++) { | 
|  | 82 | +      double t = pulse_t + i * pulse_dt; | 
|  | 83 | +      const std::size_t sampleIdx = | 
|  | 84 | +          static_cast<std::size_t>(std::floor(t / m_cfg.time_window)) - timeIdx_begin; | 
|  | 85 | + | 
|  | 86 | +      adc_counter++; | 
|  | 87 | + | 
|  | 88 | +      // Measure amplitudes for ADC | 
|  | 89 | +      if (adc_counter == sample_tick) { | 
|  | 90 | +        raw_sample.setAmplitude(sampleIdx, pulse.getAmplitude()[i]); | 
|  | 91 | +        adc_counter = 0; | 
|  | 92 | +        if (tot_progress) | 
|  | 93 | +          raw_sample.setTotProgress(sampleIdx); | 
|  | 94 | +      } | 
|  | 95 | + | 
|  | 96 | +      // Measure up-crossing time for TOA | 
|  | 97 | +      if (!tot_progress && pulse.getAmplitude()[i] > m_cfg.toa_thres) { | 
|  | 98 | +        toaIdx = sampleIdx; | 
|  | 99 | +        raw_sample.setTOA(sampleIdx, | 
|  | 100 | +                          get_crossing_time(m_cfg.toa_thres, pulse_dt, t, pulse.getAmplitude()[i], | 
|  | 101 | +                                            pulse.getAmplitude()[i - 1])); | 
|  | 102 | +        tot_progress = true; | 
|  | 103 | +        tot_complete = false; | 
|  | 104 | +        raw_sample.setTotProgress(sampleIdx); | 
|  | 105 | +      } | 
|  | 106 | + | 
|  | 107 | +      // Measure down-crossing time for TOT | 
|  | 108 | +      if (tot_progress && !tot_complete && pulse.getAmplitude()[i] < m_cfg.tot_thres) { | 
|  | 109 | +        raw_sample.setTOT(toaIdx, | 
|  | 110 | +                          get_crossing_time(m_cfg.tot_thres, pulse_dt, t, pulse.getAmplitude()[i], | 
|  | 111 | +                                            pulse.getAmplitude()[i - 1])); | 
|  | 112 | +        tot_progress = false; | 
|  | 113 | +        tot_complete = true; | 
|  | 114 | +        raw_sample.setTotComplete(sampleIdx); | 
|  | 115 | +      } | 
|  | 116 | +    } | 
|  | 117 | + | 
|  | 118 | +    // Fill HGCROCSamples and RawHGCROCHit | 
|  | 119 | +    auto out_digi_hit = out_digi_hits->create(); | 
|  | 120 | +    out_digi_hit.setCellID(pulse.getCellID()); | 
|  | 121 | + | 
|  | 122 | +    const auto& entries = raw_sample.getEntries(); | 
|  | 123 | + | 
|  | 124 | +    for (const auto& entry : entries) { | 
|  | 125 | +      edm4eic::HGCROCSample sample; | 
|  | 126 | +      auto adc   = std::max(std::llround(entry.adc / m_cfg.dyRangeADC * m_cfg.capADC), 0LL); | 
|  | 127 | +      sample.ADC = adc > m_cfg.capADC ? m_cfg.capADC : adc; | 
|  | 128 | +      auto toa   = std::max(std::llround(entry.toa / m_cfg.dyRangeTOA * m_cfg.capTOA), 0LL); | 
|  | 129 | +      sample.timeOfArrival = toa > m_cfg.capTOA ? m_cfg.capTOA : toa; | 
|  | 130 | +      auto tot = std::max(std::llround(entry.tot / m_cfg.dyRangeTOT * m_cfg.capTOT), 0LL); | 
|  | 131 | +      sample.timeOverThreshold = tot > m_cfg.capTOT ? m_cfg.capTOT : tot; | 
|  | 132 | +      sample.TOTInProgress     = entry.totProgress; | 
|  | 133 | +      sample.TOTComplete       = entry.totComplete; | 
|  | 134 | +      out_digi_hit.addToSamples(sample); | 
|  | 135 | +    } | 
|  | 136 | +  } | 
|  | 137 | +} // PulseDigi:process | 
|  | 138 | + | 
|  | 139 | +double PulseDigi::get_crossing_time(double thres, double dt, double t, double amp1, | 
|  | 140 | +                                    double amp2) const { | 
|  | 141 | +  double numerator   = (thres - amp2) * dt; | 
|  | 142 | +  double denominator = amp2 - amp1; | 
|  | 143 | +  double added       = t; | 
|  | 144 | +  return (numerator / denominator) + added; | 
|  | 145 | +} | 
| 143 | 146 | } // namespace eicrecon | 
0 commit comments