Belle II Software development
SVDClusterizerDirectModule.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9#include <svd/modules/svdReconstruction/SVDClusterizerDirectModule.h>
10
11#include <framework/datastore/StoreArray.h>
12#include <framework/logging/Logger.h>
13
14#include <vxd/geometry/GeoCache.h>
15#include <svd/geometry/SensorInfo.h>
16
17#include <mdst/dataobjects/MCParticle.h>
18#include <svd/dataobjects/SVDTrueHit.h>
19#include <svd/dataobjects/SVDShaperDigit.h>
20#include <svd/dataobjects/SVDCluster.h>
21#include <mva/dataobjects/DatabaseRepresentationOfWeightfile.h>
22
23#include <svd/reconstruction/NNWaveFitTool.h>
24
25#include <unordered_map>
26#include <algorithm>
27#include <numeric>
28#include <functional>
29#include <cassert>
30
31using namespace std;
32using namespace std::placeholders;
33using namespace Belle2;
34using namespace Belle2::SVD;
35
36//-----------------------------------------------------------------
37// Register the Module
38//-----------------------------------------------------------------
39REG_MODULE(SVDClusterizerDirect);
40
41//-----------------------------------------------------------------
42// Implementation
43//-----------------------------------------------------------------
44
46{
47 B2DEBUG(200, "SVDClusterizerDirectModule ctor");
48 //Set module properties
49 setDescription("Clusterize SVDShaperDigits and reconstruct hits");
51
52 // 1. Collections.
54 "6-digits collection name", string(""));
55 addParam("Clusters", m_storeClustersName,
56 "Cluster collection name", string(""));
57 addParam("TrueHits", m_storeTrueHitsName,
58 "TrueHit collection name", string(""));
59 addParam("MCParticles", m_storeMCParticlesName,
60 "MCParticles collection name", string(""));
61 addParam("SVDEventInfo", m_svdEventInfoName,
62 "SVDEventInfo name", string(""));
63
64 // 2. Calibration and time fitter sources
65 addParam("TimeFitterName", m_timeFitterName,
66 "Name of time fitter data file", string("SVDTimeNet_6samples"));
67 addParam("CalibratePeak", m_calibratePeak, "Use calibrattion (vs. default) for peak widths and positions", bool(false));
68
69 // 3. Clustering
70 // FIXME: Idiotic names of parameters kept for compatibility with the old clusterizer.
71 addParam("NoiseSN", m_cutAdjacent,
72 "SN for digits to be considered for clustering", m_cutAdjacent);
73 addParam("SeedSN", m_cutSeed,
74 "SN for digits to be considered as seed", m_cutSeed);
75 addParam("ClusterSN", m_cutCluster,
76 "Minimum SN for clusters", m_cutCluster);
77 addParam("HeadTailSize", m_sizeHeadTail,
78 "Cluster size at which to switch to Analog head tail algorithm", m_sizeHeadTail);
79
80}
81
83{
84 //Register collections
89
90 storeClusters.registerInDataStore();
91 storeShaperDigits.isRequired();
92 storeTrueHits.isOptional();
93 storeMCParticles.isOptional();
94 m_storeSVDEvtInfo.isRequired();
95
96 if (!m_storeSVDEvtInfo.isOptional(m_svdEventInfoName)) m_svdEventInfoName = "SVDEventInfoSim";
98
99 RelationArray relClusterShaperDigits(storeClusters, storeShaperDigits);
100 RelationArray relClusterTrueHits(storeClusters, storeTrueHits);
101 RelationArray relClusterMCParticles(storeClusters, storeMCParticles);
102 RelationArray relShaperDigitTrueHits(storeShaperDigits, storeTrueHits);
103 RelationArray relShaperDigitMCParticles(storeShaperDigits, storeMCParticles);
104
105 relClusterShaperDigits.registerInDataStore();
106 //Relations to simulation objects only if the ancestor relations exist
107 if (relShaperDigitTrueHits.isOptional())
108 relClusterTrueHits.registerInDataStore();
109 if (relShaperDigitMCParticles.isOptional())
110 relClusterMCParticles.registerInDataStore();
111
112 //Store names to speed up creation later
113 m_storeClustersName = storeClusters.getName();
114 m_storeShaperDigitsName = storeShaperDigits.getName();
115 m_storeTrueHitsName = storeTrueHits.getName();
116 m_storeMCParticlesName = storeMCParticles.getName();
117
118 m_relClusterShaperDigitName = relClusterShaperDigits.getName();
119 m_relClusterTrueHitName = relClusterTrueHits.getName();
120 m_relClusterMCParticleName = relClusterMCParticles.getName();
121 m_relShaperDigitTrueHitName = relShaperDigitTrueHits.getName();
122 m_relShaperDigitMCParticleName = relShaperDigitMCParticles.getName();
123
124 B2INFO(" 1. COLLECTIONS:");
125 B2INFO(" --> MCParticles: " << m_storeMCParticlesName);
126 B2INFO(" --> Digits: " << m_storeShaperDigitsName);
127 B2INFO(" --> Clusters: " << m_storeClustersName);
128 B2INFO(" --> TrueHits: " << m_storeTrueHitsName);
129 B2INFO(" --> DigitMCRel: " << m_relShaperDigitMCParticleName);
130 B2INFO(" --> ClusterMCRel: " << m_relClusterMCParticleName);
131 B2INFO(" --> ClusterDigitRel: " << m_relClusterShaperDigitName);
132 B2INFO(" --> DigitTrueRel: " << m_relShaperDigitTrueHitName);
133 B2INFO(" --> ClusterTrueRel: " << m_relClusterTrueHitName);
134 B2INFO(" 2. CALIBRATION DATA:");
135 B2INFO(" --> Time NN: " << m_timeFitterName);
136 B2INFO(" 4. CLUSTERING:");
137 B2INFO(" --> Neighbour cut: " << m_cutAdjacent);
138 B2INFO(" --> Seed cut: " << m_cutSeed);
139 B2INFO(" --> Cluster charge cut: " << m_cutCluster);
140 B2INFO(" --> HT for clusters >: " << m_sizeHeadTail);
141
142 // Properly initialize the NN time fitter
143 // FIXME: Should be moved to beginRun
144 // FIXME: No support for 3/6 sample switching within a run/event
146 m_fitter.setNetwrok(dbXml->m_data);
147}
148
150 RelationLookup& lookup, size_t digits)
151{
152 lookup.clear();
153 //If we don't have a relation we don't build a lookuptable
154 if (!relation) return;
155 //Resize to number of digits and set all values
156 lookup.resize(digits);
157 for (const auto& element : relation) {
158 lookup[element.getFromIndex()] = &element;
159 }
160}
161
163 std::map<unsigned int, float>& relation, unsigned int index)
164{
165 //If the lookup table is not empty and the element is set
166 if (!lookup.empty() && lookup[index]) {
167 const RelationElement& element = *lookup[index];
168 const unsigned int size = element.getSize();
169 //Add all Relations to the map
170 for (unsigned int i = 0; i < size; ++i) {
171 //negative weights are from ignored particles, we don't like them and
172 //thus ignore them :D
173 if (element.getWeight(i) < 0) continue;
174 relation[element.getToIndex(i)] += element.getWeight(i);
175 }
176 }
177}
178
180{
181
183 // If no digits or no SVDEventInfo, nothing to do
184 if (!storeShaperDigits || !storeShaperDigits.getEntries() || !m_storeSVDEvtInfo.isValid()) return;
185
186 SVDModeByte modeByte = m_storeSVDEvtInfo->getModeByte();
187
188 size_t nDigits = storeShaperDigits.getEntries();
189 B2DEBUG(90, "Initial size of StoreDigits array: " << nDigits);
190
191 const StoreArray<MCParticle> storeMCParticles(m_storeMCParticlesName);
192 const StoreArray<SVDTrueHit> storeTrueHits(m_storeTrueHitsName);
193
194 RelationArray relShaperDigitMCParticle(storeShaperDigits, storeMCParticles, m_relShaperDigitMCParticleName);
195 RelationArray relShaperDigitTrueHit(storeShaperDigits, storeTrueHits, m_relShaperDigitTrueHitName);
196
198 storeClusters.clear();
199
200 RelationArray relClusterMCParticle(storeClusters, storeMCParticles,
202 if (relClusterMCParticle) relClusterMCParticle.clear();
203
204 RelationArray relClusterShaperDigit(storeClusters, storeShaperDigits,
206 if (relClusterShaperDigit) relClusterShaperDigit.clear();
207
208 RelationArray relClusterTrueHit(storeClusters, storeTrueHits,
210 if (relClusterTrueHit) relClusterTrueHit.clear();
211
212 //Build lookup tables for relations
213 createRelationLookup(relShaperDigitMCParticle, m_mcRelation, nDigits);
214 createRelationLookup(relShaperDigitTrueHit, m_trueRelation, nDigits);
215
216 // Create fit tool object
218
219 // I. Group digits by sensor/side.
220 vector<pair<unsigned short, unsigned short> > sensorDigits;
221 VxdID lastSensorID(0);
222 size_t firstSensorDigit = 0;
223 for (size_t iDigit = 0; iDigit < nDigits; ++iDigit) {
224 const SVDShaperDigit& adigit = *storeShaperDigits[iDigit];
225 VxdID sensorID = adigit.getSensorID();
226 sensorID.setSegmentNumber(adigit.isUStrip() ? 1 : 0);
227 if (sensorID != lastSensorID) { // we have a new sensor side
228 sensorDigits.push_back(make_pair(firstSensorDigit, iDigit));
229 firstSensorDigit = iDigit;
230 lastSensorID = sensorID;
231 }
232 }
233 // save last VxdID
234 sensorDigits.push_back(make_pair(firstSensorDigit, nDigits));
235
236 // ICYCLE OVER SENSORS
237 for (auto id_indices : sensorDigits) {
238 // Retrieve parameters from sensorDigits
239 unsigned int firstDigit = id_indices.first;
240 unsigned int lastDigit = id_indices.second;
241 // Get VXDID from the first digit
242 const SVDShaperDigit& sampleDigit = *storeShaperDigits[firstDigit];
243 VxdID sensorID = sampleDigit.getSensorID();
244 bool isU = sampleDigit.isUStrip();
245
246 // Retrieve sensor parameters from GeoCache
247 const SensorInfo& info = dynamic_cast<const SensorInfo&>(VXD::GeoCache::getInstance().getSensorInfo(sensorID));
248 double pitch = isU ? info.getUPitch() : info.getVPitch();
249
250 // 4. Cycle through digits and form clusters on the way.
251
252 // Identify contiguous groups of strips
253 // We maintain two pointers: firstDigit and lastDigit. We also keep track of last seen
254 // strip number.
255 // We start with firstDigit = lastDigit = 0
256 // If there is a gap in strip numbers or the current digit/strip is invalid, we save the
257 // current group (fd, cd), unless empty, and set fd = cd.
258 // If current strip is valid, we include it: ld = cd+1.
259
260 vector<pair<size_t, size_t> > stripGroups;
261 unordered_map<size_t, apvSamples> storedNormedSamples;
262 size_t firstClusterDigit = firstDigit;
263 size_t lastClusterDigit = firstDigit;
264 short lastStrip = -2;
265
266 B2DEBUG(300, "Clustering digits " << firstDigit << " to " << lastDigit);
267 for (size_t iDigit = firstDigit; iDigit < lastDigit; ++iDigit) {
268
269 const SVDShaperDigit& digit = *storeShaperDigits[iDigit];
270 unsigned short currentStrip = digit.getCellID();
271 B2DEBUG(300, "Digit " << iDigit << ", strip: " << currentStrip << ", lastStrip: " << lastStrip);
272 B2DEBUG(300, "First CD: " << firstClusterDigit << " Last CD: " << lastClusterDigit);
273 // Get calibration
274 float stripNoiseADU = m_noiseCal.getNoise(sensorID, isU, currentStrip);
275
276 // If the strip is not masked away, save normalized samples (sample/stripNoise)
277 apvSamples normedSamples;
278
279 auto samples = digit.getSamples();
280 transform(samples.begin(), samples.end(), normedSamples.begin(),
281 bind(divides<float>(), _1, stripNoiseADU));
282 bool validDigit = pass3Samples(normedSamples, m_cutAdjacent);
283
284 // If this is a valid digit, store normed samples.
285 // It's not just that we don't have to normalize again - we don't need other data than strip
286 // number and normed samples from the digits!
287 zeroSuppress(normedSamples, m_cutAdjacent);
288 storedNormedSamples.insert(make_pair(iDigit, normedSamples));
289
290
291 // See if we have a gap in strip numbers.
292 bool consecutive = ((currentStrip - lastStrip) == 1);
293 lastStrip = currentStrip;
294
295 B2DEBUG(300, (validDigit ? "Valid " : "Invalid ") << (consecutive ? "consecutive" : "gap"));
296
297 // If this is not a valid digit or we have a gap, we save if there is something to save
298 if ((!validDigit || !consecutive) && (firstClusterDigit < lastClusterDigit)) {
299 B2DEBUG(300, "Saving (" << firstClusterDigit << ", " << lastClusterDigit << ")");
300 stripGroups.emplace_back(firstClusterDigit, lastClusterDigit);
301 }
302
303 // Update counters
304 if (validDigit) {
305 if (consecutive) {// consecutive and valid : include in cluster
306 lastClusterDigit = iDigit + 1;
307 } else { // valid after a gap = new cluster
308 firstClusterDigit = iDigit;
309 lastClusterDigit = iDigit + 1;
310 }
311 } else { // invalid digit : skip
312 firstClusterDigit = iDigit + 1;
313 lastClusterDigit = iDigit + 1;
314 }
315 } // for digits
316 // save last if any
317 if (firstClusterDigit < lastClusterDigit) {
318 B2DEBUG(300, "Saving (" << firstClusterDigit << ", " << lastDigit << ")");
319 stripGroups.emplace_back(firstClusterDigit, lastDigit);
320 }
321
322 ostringstream os;
323 os << sensorID << "NormedSamples: " << endl;
324 for (auto item : storedNormedSamples) {
325 os << item.first << ": ";
326 copy(item.second.begin(), item.second.end(), ostream_iterator<double>(os, " "));
327 os << endl;
328 }
329 os << "StripGroups: " << endl;
330 for (auto item : stripGroups) {
331 os << "(" << item.first << ", " << item.second << "), ";
332 }
333 B2DEBUG(300, os.str());
334
335 // 5. Now we loop through stripGroups and form clusters.
336 // We are still keeping digit indices, so we are keeping hold on the relations as wwell.
337
338 // FIXME: Pretty stupid. There should be ClusterCandidate type to hold all this, as
339 // well as normed samples. It could calculate a lot of things by itself and make the
340 // code cleaner.
341 vector<unsigned short> stripNumbers;
342 vector<float> stripPositions;
343 vector<float> stripNoises;
344 vector<float> timeShifts;
345 vector<float> waveWidths;
346
347 for (auto clusterBounds : stripGroups) {
348
349 unsigned short clusterSize = clusterBounds.second - clusterBounds.first;
350 assert(clusterSize > 0);
351
352 stripNumbers.clear();
353 stripPositions.clear();
354 stripNoises.clear();
355 timeShifts.clear();
356 waveWidths.clear();
357
358 for (size_t iDigit = clusterBounds.first; iDigit < clusterBounds.second; ++iDigit) {
359 // Now we take the normed samples and calculate the common
360 // cluster probability distribution and, from it, the time and its error, and
361 // the probability of the hit being from beam background.
362 const SVDShaperDigit& digit = *storeShaperDigits[iDigit];
363
364 unsigned short stripNo = digit.getCellID();
365 stripNumbers.push_back(stripNo);
366 // Is the calibrations interface sensible?
367 double stripNoiseADU = m_noiseCal.getNoise(sensorID, isU, stripNo);
368 stripNoises.push_back(
369 m_pulseShapeCal.getChargeFromADC(sensorID, isU, stripNo, stripNoiseADU)
370 );
371 // Some calibrations magic.
372 // FIXME: Only use calibration on real data. Until simulations correspond to
373 // default calibrtion, we cannot use it.
374 double peakWidth = 270;
375 double timeShift = isU ? 2.5 : -2.2;
376 if (m_calibratePeak) {
377 peakWidth = 1.988 * m_pulseShapeCal.getWidth(sensorID, isU, stripNo);
378 timeShift = m_pulseShapeCal.getPeakTime(sensorID, isU, stripNo)
379 - 0.25 * peakWidth;
380 }
381 // Correct with trigger bin information
382 const double triggerBinSep = 4 * 1.96516; //in ns
383 double apvPhase = triggerBinSep * (0.5 + static_cast<int>(modeByte.getTriggerBin()));
384 timeShift = timeShift + apvPhase;
385 waveWidths.push_back(peakWidth);
386 timeShifts.push_back(timeShift);
387 stripPositions.push_back(
388 isU ? info.getUCellPosition(stripNo) : info.getVCellPosition(stripNo)
389 );
390 stripPositions.push_back(
391 isU ? info.getUCellPosition(stripNo) : info.getVCellPosition(stripNo)
392 );
393 }
394 // The formulas don't take into account differences in strip noises. So we take RMS
395 // strip noise as cluster noise.
396 // cluster noise is in e-!
397 float clusterNoise = sqrt(
398 1.0 / clusterSize
399 * inner_product(stripNoises.begin(), stripNoises.end(), stripNoises.begin(), 0.0)
400 );
401 B2DEBUG(200, "RMS cluster noise: " << clusterNoise);
402
403 // This will hold component pdfs. We may want to rememeber them to study
404 // homogeneity of cluster times.
405 shared_ptr<nnFitterBinData> pStrip;
406 // This will aggregate the components pdfs to get cluster time pdf
407 nnFitterBinData pCluster(m_fitter.getBinCenters().size());
408 fill(pCluster.begin(), pCluster.end(), double(1.0));
409 for (size_t iClusterStrip = 0; iClusterStrip < clusterSize; ++iClusterStrip) {
410 size_t iDigit = clusterBounds.first + iClusterStrip;
411 ostringstream os1;
412 os1 << "Input to NNFitter: iDigit = " << iDigit << endl << "Samples: ";
413 copy(storedNormedSamples[iDigit].begin(), storedNormedSamples[iDigit].end(),
414 ostream_iterator<double>(os1, " "));
415 os1 << endl;
416 pStrip = m_fitter.getFit(storedNormedSamples[iDigit], waveWidths[iClusterStrip]);
417 os1 << "Output from NNWaveFitter: " << endl;
418 copy(pStrip->begin(), pStrip->end(), ostream_iterator<double>(os1, " "));
419 os1 << endl;
420 // Apply strip time shift to pdf
421 fitTool.shiftInTime(*pStrip, -timeShifts[iClusterStrip]);
422 fitTool.multiply(pCluster, *pStrip);
423 os1 << "Accummulated: " << endl;
424 copy(pCluster.begin(), pCluster.end(), ostream_iterator<double>(os1, " "));
425 B2DEBUG(200, os1.str());
426 }
427 // We can now calculate clsuter time and its error.
428 double clusterTime, clusterTimeErr;
429 tie(clusterTime, clusterTimeErr) = fitTool.getTimeShift(pCluster);
430 B2DEBUG(200, "Time: " << clusterTime << " +/- " << clusterTimeErr);
431 // Now we have the cluster time pdf, so we can calculate amplitudes.
432 // In the next cycle thrugh cluster's digits, we calculate ampltidues and their
433 // errors.
434 vector<double> stripAmplitudes(stripNoises.size());
435 vector<double> stripAmplitudeErrors(stripNoises.size());
436 double clusterChi2 = 0.0;
437 for (size_t iClusterStrip = 0; iClusterStrip < clusterSize; ++iClusterStrip) {
438 size_t iDigit = clusterBounds.first + iClusterStrip;
439 double snAmp, snAmpError, chi2;
440 tie(snAmp, snAmpError, chi2) =
441 fitTool.getAmplitudeChi2(storedNormedSamples[iDigit], clusterTime, waveWidths[iClusterStrip]);
442 // We have to re-normalize: Multiply amplitudes in SN units by noises in electrons.
443 stripAmplitudes[iClusterStrip] =
444 stripNoises[iClusterStrip] * snAmp;
445 stripAmplitudeErrors[iClusterStrip] =
446 stripNoises[iClusterStrip] * snAmpError;
447 clusterChi2 += chi2;
448 B2DEBUG(200, "Digit " << iDigit << " Noise: " << stripNoises[iClusterStrip]
449 << " Amplitude: " << stripAmplitudes[iClusterStrip]
450 << " +/- " << stripAmplitudeErrors[iClusterStrip]
451 << " Chi2: " << chi2
452 );
453 }
454 // Cluster charge and S/N
455 float clusterCharge = accumulate(stripAmplitudes.begin(), stripAmplitudes.end(), 0.0);
456 float clusterChargeError = sqrt(
457 inner_product(stripAmplitudeErrors.begin(), stripAmplitudeErrors.end(),
458 stripAmplitudeErrors.begin(), 0.0)
459 );
460 float clusterSN = (clusterChargeError > 0) ? clusterCharge / clusterChargeError : clusterCharge;
461 clusterChi2 /= clusterSize;
462 // Cluster seed
463 size_t seedIndex = distance(stripAmplitudes.begin(), max_element(
464 stripAmplitudes.begin(), stripAmplitudes.end()));
465 //unsigned short clusterSeedStrip = stripNumbers[seedIndex];
466 float clusterSeedCharge = stripAmplitudes[seedIndex];
467 B2DEBUG(200, "Cluster parameters:");
468 B2DEBUG(200, "Charge: " << clusterCharge << " +/- " << clusterChargeError);
469 B2DEBUG(200, "Seed: " << clusterSeedCharge << " +/- " << stripAmplitudeErrors[seedIndex]);
470 B2DEBUG(200, "S/N: " << clusterSN);
471
472
473 // Hit position
474 float clusterPosition, clusterPositionError;
475
476 // Now censoring:
477 // NB: We don't censor strip charges again. The requirement of 3 consecutive samples
478 // above threshold is a good enough safeguard.
479 if ((clusterCharge < clusterNoise * m_cutCluster) || (clusterSeedCharge < clusterNoise * m_cutSeed))
480 continue;
481
482 if (clusterSize < m_sizeHeadTail) { // centre of gravity
483 clusterPosition = 1.0 / clusterCharge * inner_product(
484 stripAmplitudes.begin(), stripAmplitudes.end(), stripPositions.begin(), 0.0
485 );
486 // Strip charge equal to zero-suppression threshold
487 float phantomCharge = m_cutAdjacent * clusterNoise;
488 if (clusterSize == 1) {
489 clusterPositionError = pitch * phantomCharge / (clusterCharge + phantomCharge);
490 } else {
491 clusterPositionError = pitch * phantomCharge / clusterCharge;
492 }
493 } else { // Head-tail
494 float leftStripCharge = stripAmplitudes.front();
495 float leftPos = stripPositions.front();
496 float rightStripCharge = stripAmplitudes.back();
497 float rightPos = stripPositions.back();
498 float centreCharge = (clusterCharge - leftStripCharge - rightStripCharge) / (clusterSize - 2);
499 leftStripCharge = (leftStripCharge < centreCharge) ? leftStripCharge : centreCharge;
500 rightStripCharge = (rightStripCharge < centreCharge) ? rightStripCharge : centreCharge;
501 clusterPosition = 0.5 * (leftPos + rightPos)
502 + 0.5 * (rightStripCharge - leftStripCharge) / centreCharge * pitch;
503 float sn = centreCharge / m_cutAdjacent / clusterNoise;
504 // Rough estimates of Landau noise
505 float landauHead = leftStripCharge / centreCharge;
506 double landauTail = rightStripCharge / centreCharge;
507 clusterPositionError = 0.5 * pitch * sqrt(1.0 / sn / sn
508 + 0.5 * landauHead * landauHead + 0.5 * landauTail * landauTail);
509 }
510 //Lorentz shift correction
511 clusterPosition -= info.getLorentzShift(isU, clusterPosition);
512
513 // Finally, we save the cluster and its relations.
514 map<unsigned int, float> mc_relations;
515 map<unsigned int, float> truehit_relations;
516 vector<pair<unsigned int, float> > digit_weights;
517 digit_weights.reserve(clusterSize);
518
519 for (size_t iDigit = clusterBounds.first; iDigit < clusterBounds.second; ++iDigit) {
520 // Store relations to MCParticles and SVDTrueHits
521 fillRelationMap(m_mcRelation, mc_relations, iDigit);
522 fillRelationMap(m_trueRelation, truehit_relations, iDigit);
523 //Add digit to the Cluster->Digit relation list
524 digit_weights.emplace_back(iDigit, stripAmplitudes[iDigit - clusterBounds.first]);
525 }
526
527 //Store the cluster into Datastore ...
528 int clsIndex = storeClusters.getEntries();
529 VxdID clusterSensorID = sensorID;
530 clusterSensorID.setSegmentNumber(0);
531 storeClusters.appendNew(
532 SVDCluster(sensorID, isU, clusterPosition, clusterPositionError, clusterTime,
533 clusterTimeErr, clusterCharge, clusterSeedCharge, clusterSize, clusterSN, clusterChi2)
534 );
535
536 //Create relations to the cluster
537 if (!mc_relations.empty()) {
538 relClusterMCParticle.add(clsIndex, mc_relations.begin(), mc_relations.end());
539 }
540 if (!truehit_relations.empty()) {
541 relClusterTrueHit.add(clsIndex, truehit_relations.begin(), truehit_relations.end());
542 }
543 relClusterShaperDigit.add(clsIndex, digit_weights.begin(), digit_weights.end());
544 } // CYCLE OVER CLUSTERS for stripGroups
545
546 } // CYCLE OVER SENSORS for items in sensorDigits
547
548 B2DEBUG(100, "Number of clusters: " << storeClusters.getEntries());
549
550} // event()
551
552
Class for accessing objects in the database.
Definition: DBObjPtr.h:21
Base class for Modules.
Definition: Module.h:72
void setDescription(const std::string &description)
Sets the description of the module.
Definition: Module.cc:214
void setPropertyFlags(unsigned int propertyFlags)
Sets the flags for the module properties.
Definition: Module.cc:208
@ c_ParallelProcessingCertified
This module can be run in parallel processing mode safely (All I/O must be done through the data stor...
Definition: Module.h:80
Low-level class to create/modify relations between StoreArrays.
Definition: RelationArray.h:62
void add(index_type from, index_type to, weight_type weight=1.0)
Add a new element to the relation.
void clear() override
Clear all elements from the relation.
Class to store a single element of a relation.
The SVD Cluster class This class stores all information about reconstructed SVD clusters.
Definition: SVDCluster.h:29
Class to store SVD mode information.
Definition: SVDModeByte.h:69
baseType getTriggerBin() const
Get the triggerBin id.
Definition: SVDModeByte.h:140
float getNoise(const VxdID &sensorID, const bool &isU, const unsigned short &strip) const
This is the method for getting the noise.
double getChargeFromADC(const Belle2::VxdID &sensorID, const bool &isU, const unsigned short &strip, const double &pulseADC) const
Return the charge (number of electrons/holes) collected on a specific strip, given the number of ADC ...
float getPeakTime(const VxdID &sensorID, const bool &isU, const unsigned short &strip) const
Return the peaking time of the strip.
float getWidth(const VxdID &sensorID, const bool &isU, const unsigned short &strip) const
Return the width of the pulse shape for a given strip.
The SVD ShaperDigit class.
VxdID getSensorID() const
Get the sensor ID.
APVFloatSamples getSamples() const
Get array of samples.
short int getCellID() const
Get strip ID.
bool isUStrip() const
Get strip direction.
The class holds arrays of bins and bin centers, and a wave generator object containing information on...
Definition: NNWaveFitTool.h:91
std::tuple< double, double, double > getAmplitudeChi2(const apvSamples &samples, double timeShift, double tau)
Return std::tuple with signal amplitude, its error, and chi2 of the fit.
std::tuple< double, double > getTimeShift(const nnFitterBinData &p)
Return std::tuple containing time shift and its error.
void shiftInTime(nnFitterBinData &p, double timeShift)
Shift the probability array in time.
void multiply(nnFitterBinData &p, const nnFitterBinData &p1)
Multiply probabilities.
void setNetwrok(const std::string &xmlData)
Set proper network definition file.
const NNWaveFitTool & getFitTool() const
Get a handle to a NNWaveFit object.
Definition: NNWaveFitter.h:98
const nnFitterBinData & getBinCenters() const
Get bin times of the network output.
Definition: NNWaveFitter.h:104
std::shared_ptr< nnFitterBinData > getFit(const apvSamples &samples, double tau)
Fitting method Send data and get rseult structure.
std::string m_relShaperDigitMCParticleName
Name of the relation between SVDShaperDigits and MCParticles.
virtual void initialize() override
Initialize the module.
std::string m_storeShaperDigitsName
Name of the collection to use for the SVDShaperDigits.
virtual void event() override
do the clustering
std::vector< const RelationElement * > RelationLookup
Container for a RelationArray Lookup table.
double m_cutCluster
Cluster cut in units of m_elNoise.
std::string m_storeTrueHitsName
Name of the collection to use for the SVDTrueHits.
std::string m_timeFitterName
Name of the time fitter (db label)
void fillRelationMap(const RelationLookup &lookup, std::map< unsigned int, float > &relation, unsigned int index)
Add the relation from a given SVDShaperDigit index to a map.
std::string m_storeMCParticlesName
Name of the collection to use for the MCParticles.
std::string m_relShaperDigitTrueHitName
Name of the relation between SVDShaperDigits and SVDTrueHits.
RelationLookup m_trueRelation
Lookup table for SVDShaperDigit->SVDTrueHit relation.
SVDPulseShapeCalibrations m_pulseShapeCal
Calibrations: pusle shape and gain.
std::string m_svdEventInfoName
Name of the SVDEventInfo object.
SVDNoiseCalibrations m_noiseCal
Calibrations: noise.
int m_sizeHeadTail
Size of the cluster at which we switch from Center of Gravity to Analog Head Tail.
void createRelationLookup(const RelationArray &relation, RelationLookup &lookup, size_t digits)
Create lookup maps for relations We do not use the RelationIndex as we know much more about the relat...
std::string m_storeClustersName
Name of the collection to use for the SVDClusters.
double m_cutSeed
Seed cut in units of m_elNoise.
std::string m_relClusterShaperDigitName
Name of the relation between SVDClusters and SVDShaperDigits.
std::string m_relClusterMCParticleName
Name of the relation between SVDClusters and MCParticles.
RelationLookup m_mcRelation
Lookup table for SVDShaperDigit->MCParticle relation.
SVDClusterizerDirectModule()
Constructor defining the parameters.
StoreObjPtr< SVDEventInfo > m_storeSVDEvtInfo
Storage for SVDEventInfo object.
bool m_calibratePeak
Use peak widths and peak time calibrations? Unitl this is also simulated, set to true only for testbe...
double m_cutAdjacent
Noise (cluster member) cut in units of m_elNoise.
std::string m_relClusterTrueHitName
Name of the relation between SVDClusters and SVDTrueHits.
Specific implementation of SensorInfo for SVD Sensors which provides additional sensor specific infor...
Definition: SensorInfo.h:25
const std::string & getName() const
Return name under which the object is saved in the DataStore.
bool isRequired(const std::string &name="")
Ensure this array/object has been registered previously.
bool isOptional(const std::string &name="")
Tell the DataStore about an optional input.
bool registerInDataStore(DataStore::EStoreFlags storeFlags=DataStore::c_WriteOut)
Register the object/array in the DataStore.
Accessor to arrays stored in the data store.
Definition: StoreArray.h:113
T * appendNew()
Construct a new T object at the end of the array.
Definition: StoreArray.h:246
int getEntries() const
Get the number of objects in the array.
Definition: StoreArray.h:216
void clear() override
Delete all entries in this array.
Definition: StoreArray.h:207
const SensorInfoBase & getSensorInfo(Belle2::VxdID id) const
Return a referecne to the SensorInfo of a given SensorID.
Definition: GeoCache.cc:67
static GeoCache & getInstance()
Return a reference to the singleton instance.
Definition: GeoCache.cc:214
Class to uniquely identify a any structure of the PXD and SVD.
Definition: VxdID.h:33
void setSegmentNumber(baseType segment)
Set the sensor segment.
Definition: VxdID.h:113
void addParam(const std::string &name, T &paramVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
Definition: Module.h:560
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:650
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28
Namespace to encapsulate code needed for simulation and reconstrucion of the SVD.
Definition: GeoSVDCreator.h:23
std::array< apvSampleBaseType, nAPVSamples > apvSamples
vector od apvSample BaseType objects
std::vector< double > nnFitterBinData
Vector of values defined for bins, such as bin times or bin probabilities.
Definition: NNWaveFitTool.h:30
void zeroSuppress(T &a, double thr)
pass zero suppression
bool pass3Samples(const T &a, double thr)
pass 3-samples
Abstract base class for different kinds of events.
STL namespace.