11 #include <svd/modules/svdReconstruction/SVDCoGTimeEstimatorModule.h>
27 setDescription(
"From SVDShaperDigit to SVDRecoDigit. Strip charge is evaluated as the max of the 6 samples; hit time is evaluated as a corrected Centre of Gravity (CoG) time.");
28 setPropertyFlags(c_ParallelProcessingCertified);
30 addParam(
"SVDEventInfo", m_svdEventInfoName,
31 "SVDEventInfo name",
string(
""));
32 addParam(
"ShaperDigits", m_storeShaperDigitsName,
33 "ShaperDigits collection name",
string(
""));
34 addParam(
"RecoDigits", m_storeRecoDigitsName,
35 "RecoDigits collection name",
string(
""));
36 addParam(
"StripPeakTimeCorrection", m_corrPeakTime,
37 "Correct for the different peaking times of the strips, obtained from local run calibration",
true);
38 addParam(
"CalibrationWithEventT0", m_calEventT0,
39 "Use the timing informations of the EventT0 in order to calibrate the CoG.",
52 m_storeTrueHits.isOptional();
53 m_storeMCParticles.isOptional();
56 m_storeShaper.isRequired(m_storeShaperDigitsName);
58 if (!m_storeSVDEvtInfo.isOptional(m_svdEventInfoName)) m_svdEventInfoName =
"SVDEventInfoSim";
59 m_storeSVDEvtInfo.isRequired(m_svdEventInfoName);
64 RelationArray relRecoDigitShaperDigits(m_storeReco, m_storeShaper);
66 RelationArray relRecoDigitTrueHits(m_storeReco, m_storeTrueHits);
67 RelationArray relRecoDigitMCParticles(m_storeReco, m_storeMCParticles);
68 RelationArray relShaperDigitTrueHits(m_storeShaper, m_storeTrueHits);
69 RelationArray relShaperDigitMCParticles(m_storeShaper, m_storeMCParticles);
78 m_storeTrueHitsName = m_storeTrueHits.getName();
79 m_storeMCParticlesName = m_storeMCParticles.getName();
81 m_relRecoDigitShaperDigitName = relRecoDigitShaperDigits.
getName();
82 m_relRecoDigitTrueHitName = relRecoDigitTrueHits.
getName();
83 m_relRecoDigitMCParticleName = relRecoDigitMCParticles.
getName();
84 m_relShaperDigitTrueHitName = relShaperDigitTrueHits.
getName();
85 m_relShaperDigitMCParticleName = relShaperDigitMCParticles.
getName();
87 B2DEBUG(1,
" 1. COLLECTIONS:");
88 B2DEBUG(1,
" --> MCParticles: " << m_storeMCParticlesName);
89 B2DEBUG(1,
" --> Digits: " << m_storeShaperDigitsName);
90 B2DEBUG(1,
" --> RecoDigits: " << m_storeRecoDigitsName);
91 B2DEBUG(1,
" --> TrueHits: " << m_storeTrueHitsName);
92 B2DEBUG(1,
" --> DigitMCRel: " << m_relShaperDigitMCParticleName);
93 B2DEBUG(1,
" --> RecoDigitMCRel: " << m_relRecoDigitMCParticleName);
94 B2DEBUG(1,
" --> RecoDigitDigitRel: " << m_relRecoDigitShaperDigitName);
95 B2DEBUG(1,
" --> DigitTrueRel: " << m_relShaperDigitTrueHitName);
96 B2DEBUG(1,
" --> RecoDigitTrueRel: " << m_relRecoDigitTrueHitName);
108 std::vector<float> probabilities = {0.5};
111 if (!m_storeShaper || !m_storeShaper.getEntries() || !m_storeSVDEvtInfo.isValid())
return;
113 SVDModeByte modeByte = m_storeSVDEvtInfo->getModeByte();
114 size_t nDigits = m_storeShaper.getEntries();
116 RelationArray relRecoDigitShaperDigit(m_storeReco, m_storeShaper,
117 m_relRecoDigitShaperDigitName);
118 if (relRecoDigitShaperDigit) relRecoDigitShaperDigit.
clear();
121 RelationArray relShaperDigitMCParticle(m_storeShaper, m_storeMCParticles, m_relShaperDigitMCParticleName);
122 RelationArray relShaperDigitTrueHit(m_storeShaper, m_storeTrueHits, m_relShaperDigitTrueHitName);
124 RelationArray relRecoDigitMCParticle(m_storeReco, m_storeMCParticles,
125 m_relRecoDigitMCParticleName);
126 if (relRecoDigitMCParticle) relRecoDigitMCParticle.
clear();
129 RelationArray relRecoDigitTrueHit(m_storeReco, m_storeTrueHits,
130 m_relRecoDigitTrueHitName);
131 if (relRecoDigitTrueHit) relRecoDigitTrueHit.
clear();
134 createRelationLookup(relShaperDigitMCParticle, m_mcRelation, nDigits);
135 createRelationLookup(relShaperDigitTrueHit, m_trueRelation, nDigits);
140 m_NumberOfAPVSamples = m_storeSVDEvtInfo->getNSamples();
142 B2DEBUG(1,
"number of APV samples = " << m_NumberOfAPVSamples);
146 m_StopCreationReco =
false;
149 if (m_StopCreationReco)
152 samples_vec = shaper.getSamples();
155 VxdID thisSensorID = shaper.getSensorID();
156 bool thisSide = shaper.isUStrip();
157 int thisCellID = shaper.getCellID();
160 m_weightedMeanTime = CalculateWeightedMeanPeakTime(samples_vec);
161 if (m_StopCreationReco)
163 m_amplitude = CalculateAmplitude(samples_vec);
164 m_amplitudeError = CalculateAmplitudeError(thisSensorID, thisSide, thisCellID);
167 m_weightedMeanTimeError = CalculateWeightedMeanPeakTimeError(samples_vec);
169 m_chi2 = CalculateChi2();
172 if (m_amplitude > 255)
173 B2DEBUG(10,
"ERROR: m_amplitude = " << m_amplitude <<
", should be <= 255");
177 m_amplitude = m_PulseShapeCal.getChargeFromADC(thisSensorID, thisSide, thisCellID, m_amplitude);
178 m_amplitudeError = m_PulseShapeCal.getChargeFromADC(thisSensorID, thisSide, thisCellID, m_amplitudeError);
181 m_weightedMeanTime -= m_PulseShapeCal.getPeakTime(thisSensorID, thisSide, thisCellID);
185 m_weightedMeanTime = m_TimeCal.getCorrectedTime(thisSensorID, thisSide, thisCellID, m_weightedMeanTime, triggerBin);
186 m_weightedMeanTimeError = m_TimeCal.getCorrectedTimeError(thisSensorID, thisSide, thisCellID, m_weightedMeanTime,
187 m_weightedMeanTimeError, triggerBin);
191 if (m_amplitude > 100000) {
192 B2DEBUG(42,
"Charge = " << m_amplitude);
193 B2DEBUG(42,
"corresponding ADC = " << m_PulseShapeCal.getADCFromCharge(thisSensorID, thisSide, thisCellID, m_amplitude));
194 B2DEBUG(42,
"thisLayerNumber = " << thisSensorID.
getLayerNumber());
197 B2DEBUG(42,
"thisSide = " << thisSide);
198 B2DEBUG(42,
"thisCellID = " << thisCellID);
199 B2DEBUG(42,
"-----");
203 m_storeReco.appendNew(
SVDRecoDigit(shaper.getSensorID(), shaper.isUStrip(), shaper.getCellID(), m_amplitude, m_amplitudeError,
204 m_weightedMeanTime, m_weightedMeanTimeError, probabilities, m_chi2, modeByte));
207 int recoDigitIndex = m_storeReco.getEntries() - 1;
208 vector<pair<unsigned int, float> > digit_weights;
209 digit_weights.reserve(1);
210 digit_weights.emplace_back(shaper.getArrayIndex(), 1.0);
211 relRecoDigitShaperDigit.
add(recoDigitIndex, digit_weights.begin(), digit_weights.end());
214 map<unsigned int, float> mc_relations;
215 map<unsigned int, float> truehit_relations;
218 fillRelationMap(m_mcRelation, mc_relations, shaper.getArrayIndex());
219 fillRelationMap(m_trueRelation, truehit_relations, shaper.getArrayIndex());
222 if (!mc_relations.empty()) {
223 relRecoDigitMCParticle.
add(recoDigitIndex, mc_relations.begin(), mc_relations.end());
225 if (!truehit_relations.empty()) {
226 relRecoDigitTrueHit.
add(recoDigitIndex, truehit_relations.begin(), truehit_relations.end());
245 float averagetime = 0;
246 float sumAmplitudes = 0;
248 for (
int k = 0; k < m_NumberOfAPVSamples; k ++) {
249 averagetime += k * samples[k];
250 sumAmplitudes += samples[k];
252 if (sumAmplitudes != 0) {
253 averagetime /= (sumAmplitudes);
254 averagetime *= m_DeltaT;
257 m_StopCreationReco =
true;
258 B2WARNING(
"Trying to divide by 0 (ZERO)! Sum of amplitudes is NULL! Skipping this SVDShaperDigit!");
268 for (
int k = 0; k < m_NumberOfAPVSamples; k ++) {
269 if (samples[k] > amplitude)
270 amplitude = samples[k];
287 for (
int k = 0; k < m_NumberOfAPVSamples; k ++) {
289 tmpResSq += TMath::Power(k * m_DeltaT - m_weightedMeanTime, 2);
292 return m_amplitudeError / Atot * TMath::Sqrt(tmpResSq);
299 stripnoise = m_NoiseCal.getNoise(ThisSensorID, ThisSide, ThisCellID);
316 if (!relation)
return;
318 lookup.resize(digits);
319 for (
const auto& element : relation) {
320 lookup[element.getFromIndex()] = &element;
325 std::map<unsigned int, float>& relation,
unsigned int index)
328 if (!lookup.empty() && lookup[index]) {
330 const unsigned int size = element.getSize();
332 for (
unsigned int i = 0; i < size; ++i) {
335 if (element.getWeight(i) < 0)
continue;
336 relation[element.getToIndex(i)] += element.getWeight(i);