12 #include <top/modules/TOPXTalkChargeShareSetter/TOPXTalkChargeShareSetterModule.h>
15 #include <framework/datastore/StoreArray.h>
18 #include <top/dataobjects/TOPRawDigit.h>
19 #include <top/dataobjects/TOPRawWaveform.h>
20 #include <top/dataobjects/TOPDigit.h>
40 setDescription(
"Crosstalk & chargeshare flag setter");
41 setPropertyFlags(c_ParallelProcessingCertified);
44 addParam(
"sumChargeShare", m_sumChargeShare,
45 "sum up charge of PrimaryChargeShare and SecondaryChargeShare", (
bool)
false);
46 addParam(
"timeCut", m_timeCut,
47 "cut range of hittiming for chargeshare flag [ns]", (
float) 1);
48 addParam(
"preValleyDepthLoose", m_preValleyDepthLoose,
49 "loose threshold for depth of pre valley [ADC counts], for corss talk identification", 20);
50 addParam(
"preValleyDepthTight", m_preValleyDepthTight,
51 "tight threshold for depth of pre valley [ADC counts], identified as cross talk with loose threshold for the second peak amplitude",
53 addParam(
"2ndPeakAmplitudeLoose", m_2ndPeakAmplitudeLoose,
54 "loose threshold for amplitude of the second peak [ADC counts] for cross talk identification. Defined as ADC count difference between the valley just after the main peak and the second peak, Used when the \"preValleyDepthTight\" was satisfied",
56 addParam(
"2ndPeakAmplitudeRatioTight", m_2ndPeakAmplitudeRatioTight,
57 "tight threshold for amplitude ratio of the second peak to the main peak height [ADC counts]", 0.2);
58 addParam(
"nSampleBefore", m_nSampleBefore,
59 "the number of samples by which the pre-valley should exist from the CFD timing, used for cross talk identification."
60 "Set negative value not to check existence of the pre-valley.", 5);
61 addParam(
"nSampleAfter", m_nSampleAfter,
62 "the number of samples by which the second peak should exist from the CFD timing, used for cross talk identification."
63 "Set negative value not to check existence of the second peak.", 10);
67 TOPXTalkChargeShareSetterModule::~TOPXTalkChargeShareSetterModule() {}
69 void TOPXTalkChargeShareSetterModule::initialize()
74 if (m_preValleyDepthLoose < 0 || m_preValleyDepthTight < 0 || m_nSampleBefore < 0) m_checkPreValleyForXTalkId =
false;
75 if (m_2ndPeakAmplitudeLoose < 0 || m_2ndPeakAmplitudeRatioTight < 0 || m_nSampleAfter < 0) m_checkPostValleyForXTalkId =
false;
78 void TOPXTalkChargeShareSetterModule::beginRun()
82 void TOPXTalkChargeShareSetterModule::event()
88 std::map<int, TOPDigit*> hitInfoMap;
89 for (
auto& digit : digits) {
91 if (digit.getHitQuality() == TOPDigit::c_Junk)
continue;
93 const auto* rawDigit = digit.getRelated<
TOPRawDigit>();
94 if (!rawDigit)
continue;
97 if (waveform and isCrossTalk(waveform->getWaveform(), TMath::FloorNint(digit.getRawTime()), digit.getPulseHeight()))
98 digit.setHitQuality(TOPDigit::c_CrossTalk);
101 int slotId = digit.getModuleID();
102 unsigned int channelId = digit.getChannel();
103 double rawTime = digit.getRawTime();
104 for (
auto& digit2 : digits) {
105 if (digit2.getChannel() != channelId || digit2.getModuleID() != slotId)
continue;
106 if (digit2.getHitQuality() == TOPDigit::c_CrossTalk
107 and TMath::Abs(rawTime - digit2.getRawTime() - m_nCrossTalkRingingSamples / 2.) < m_nCrossTalkRingingSamples / 2.)
108 digit.setHitQuality(TOPDigit::c_CrossTalk);
112 if (digit.getHitQuality() == TOPDigit::c_CrossTalk)
continue;
114 int globalPixelId = (digit.getModuleID() - 1) * 512 + digit.getPixelID() - 1;
115 while (hitInfoMap.count(globalPixelId) > 0)globalPixelId += 10000;
116 hitInfoMap[globalPixelId] = &digit;
122 for (
auto& digit : digits) {
124 if (digit.getHitQuality() == TOPDigit::c_Junk
125 || digit.getHitQuality() == TOPDigit::c_CrossTalk)
continue;
127 short pixelId = digit.getPixelID();
128 short slotId = digit.getModuleID();
129 short pmtId = digit.getPMTNumber();
130 double hitTime = digit.getTime();
131 double charge = digit.getIntegral();
132 bool isPrimaryChargeShare =
true;
133 std::vector<TOPDigit*> vSecondaryCandidates;
135 int adjacentPixelIds[] = { pixelId - 1 - c_NPixelsPerRow, pixelId - c_NPixelsPerRow, pixelId + 1 - c_NPixelsPerRow, pixelId + 1,
136 pixelId + 1 + c_NPixelsPerRow, pixelId + c_NPixelsPerRow, pixelId - 1 + c_NPixelsPerRow, pixelId - 1
139 for (
const auto& adjacentPixelId : adjacentPixelIds) {
140 if (adjacentPixelId > 0 && adjacentPixelId <= 512) {
141 int globalPixelId = (slotId - 1) * 512 + adjacentPixelId - 1;
142 while (hitInfoMap.count(globalPixelId) > 0) {
144 if (pmtId != hitInfoMap[globalPixelId]->getPMTNumber()
145 or TMath::Abs(hitTime - hitInfoMap[globalPixelId]->getTime()) > m_timeCut) {
146 globalPixelId += 10000;
continue;
149 float adjacentIntegral = hitInfoMap[globalPixelId]->getIntegral();
151 if (charge > adjacentIntegral
152 or (charge == adjacentIntegral && pixelId > hitInfoMap[globalPixelId]->getPixelID())) {
153 vSecondaryCandidates.push_back(hitInfoMap[globalPixelId]);
154 hitInfoMap[globalPixelId]->setSecondaryChargeShare();
156 isPrimaryChargeShare =
false;
159 globalPixelId += 10000;
161 if (!isPrimaryChargeShare)
break;
165 if (isPrimaryChargeShare && vSecondaryCandidates.size() > 0) {
166 digit.setPrimaryChargeShare();
168 if (m_sumChargeShare) {
169 for (
auto& vSecondaryCandidate : vSecondaryCandidates) {
170 digit.setPulseHeight(digit.getPulseHeight() + vSecondaryCandidate->getPulseHeight());
171 vSecondaryCandidate->setPulseHeight(0.0);
173 digit.setIntegral(digit.getIntegral() + vSecondaryCandidate->getIntegral());
174 vSecondaryCandidate->setIntegral(0.0);
181 void TOPXTalkChargeShareSetterModule::endRun()
185 void TOPXTalkChargeShareSetterModule::terminate()
190 bool TOPXTalkChargeShareSetterModule::isCrossTalk(std::vector<short> wfm,
int iRawTime,
short height)
193 int nWfmSampling = wfm.size();
194 for (
int iWin = 0 ; iWin < 16 ; iWin++) {
196 int jRawTime = iRawTime - TMath::FloorNint(iRawTime / (TOPRawDigit::c_WindowSize)) * (TOPRawDigit::c_WindowSize) +
197 (TOPRawDigit::c_WindowSize) / 4 * iWin;
198 if (jRawTime > 0 && jRawTime < nWfmSampling - 1 && wfm[jRawTime] < height / 2. && wfm[jRawTime + 1] > height / 2.) {
199 bool preValleyExist =
false;
200 short preValleyDepth = -1;
201 if (!m_checkPreValleyForXTalkId) preValleyExist =
true;
203 for (
int iSample = jRawTime ; iSample - 1 > 0 ; iSample--) {
204 if (jRawTime - iSample > m_nSampleBefore)
return false;
205 else if (wfm[iSample] - wfm[iSample - 1] >= 0)
continue;
207 preValleyDepth = (-1) * wfm[iSample];
208 if (preValleyDepth < m_preValleyDepthLoose)
return false;
210 if (!m_checkPostValleyForXTalkId)
return true;
211 preValleyExist =
true;
217 if (!preValleyExist)
return false;
221 short valley_adc = 9999;
222 for (
int jSample = jRawTime ; jSample < nWfmSampling - 1 ; jSample++) {
223 if (jSample - jRawTime > m_nSampleAfter)
return false;
224 if ((wfm[jSample + 1] - wfm[jSample])*sign > 0)
continue;
226 if (sign < 0 && valley_adc > height + 1)
227 valley_adc = wfm[jSample];
228 if (sign > 0 && valley_adc < height) {
229 if (wfm[jSample] - valley_adc > (height * m_2ndPeakAmplitudeRatioTight)
230 || (preValleyDepth > m_preValleyDepthTight
231 && (wfm[jSample] - valley_adc) > m_2ndPeakAmplitudeLoose))
235 sign = (sign > 0 ? -1 : 1);