Belle II Software development
SVDNNClusterizerModule.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/SVDNNClusterizerModule.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/SVDRecoDigit.h>
21#include <svd/dataobjects/SVDCluster.h>
22#include <mva/dataobjects/DatabaseRepresentationOfWeightfile.h>
23
24#include <svd/reconstruction/NNWaveFitTool.h>
25
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(SVDNNClusterizer);
40
41//-----------------------------------------------------------------
42// Implementation
43//-----------------------------------------------------------------
44
46{
47 B2DEBUG(200, "SVDNNClusterizerModule ctor");
48 //Set module properties
49 setDescription("Clusterize SVDRecoDigits and reconstruct hits");
51
52 // 1. Collections.
54 "RecoDigits 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
62 // 2. Calibration and time fitter sources
63 addParam("TimeFitterName", m_timeFitterName,
64 "Name of time fitter data file", string("SVDTimeNet_6samples"));
65 addParam("CalibratePeak", m_calibratePeak, "Use calibrattion (vs. default) for peak widths and positions", bool(false));
66
67 // 3. Clustering
68 // FIXME: Idiotic names of parameters kept for compatibility with the old clusterizer.
69 addParam("NoiseSN", m_cutAdjacent,
70 "SN for digits to be considered for clustering", m_cutAdjacent);
71 addParam("SeedSN", m_cutSeed,
72 "SN for digits to be considered as seed", m_cutSeed);
73 addParam("ClusterSN", m_cutCluster,
74 "Minimum SN for clusters", m_cutCluster);
75 addParam("HeadTailSize", m_sizeHeadTail,
76 "Cluster size at which to switch to Analog head tail algorithm", m_sizeHeadTail);
77
78}
79
81{
82 //Register collections
87
88 storeClusters.registerInDataStore();
89 storeRecoDigits.isRequired();
90 storeTrueHits.isOptional();
91 storeMCParticles.isOptional();
92
93 RelationArray relClusterRecoDigits(storeClusters, storeRecoDigits);
94 RelationArray relClusterTrueHits(storeClusters, storeTrueHits);
95 RelationArray relClusterMCParticles(storeClusters, storeMCParticles);
96 RelationArray relRecoDigitTrueHits(storeRecoDigits, storeTrueHits);
97 RelationArray relRecoDigitMCParticles(storeRecoDigits, storeMCParticles);
98
99 relClusterRecoDigits.registerInDataStore();
100 //Relations to simulation objects only if the ancestor relations exist
101 if (relRecoDigitTrueHits.isOptional())
102 relClusterTrueHits.registerInDataStore();
103 if (relRecoDigitMCParticles.isOptional())
104 relClusterMCParticles.registerInDataStore();
105
106 //Store names to speed up creation later
107 m_storeClustersName = storeClusters.getName();
108 m_storeRecoDigitsName = storeRecoDigits.getName();
109 m_storeTrueHitsName = storeTrueHits.getName();
110 m_storeMCParticlesName = storeMCParticles.getName();
111
112 m_relClusterRecoDigitName = relClusterRecoDigits.getName();
113 m_relClusterTrueHitName = relClusterTrueHits.getName();
114 m_relClusterMCParticleName = relClusterMCParticles.getName();
115 m_relRecoDigitTrueHitName = relRecoDigitTrueHits.getName();
116 m_relRecoDigitMCParticleName = relRecoDigitMCParticles.getName();
117
118 B2INFO(" 1. COLLECTIONS:");
119 B2INFO(" --> MCParticles: " << m_storeMCParticlesName);
120 B2INFO(" --> Digits: " << m_storeRecoDigitsName);
121 B2INFO(" --> Clusters: " << m_storeClustersName);
122 B2INFO(" --> TrueHits: " << m_storeTrueHitsName);
123 B2INFO(" --> DigitMCRel: " << m_relRecoDigitMCParticleName);
124 B2INFO(" --> ClusterMCRel: " << m_relClusterMCParticleName);
125 B2INFO(" --> ClusterDigitRel: " << m_relClusterRecoDigitName);
126 B2INFO(" --> DigitTrueRel: " << m_relRecoDigitTrueHitName);
127 B2INFO(" --> ClusterTrueRel: " << m_relClusterTrueHitName);
128 B2INFO(" 2. CALIBRATION DATA:");
129 B2INFO(" --> Time NN: " << m_timeFitterName);
130 B2INFO(" 4. CLUSTERING:");
131 B2INFO(" --> Neighbour cut: " << m_cutAdjacent);
132 B2INFO(" --> Seed cut: " << m_cutSeed);
133 B2INFO(" --> Cluster charge cut: " << m_cutCluster);
134 B2INFO(" --> HT for clusters >: " << m_sizeHeadTail);
135
136 // Properly initialize the NN time fitter
137 // FIXME: Should be moved to beginRun
138 // FIXME: No support for 3/6 sample switching within a run/event
140 m_fitter.setNetwrok(dbXml->m_data);
141}
142
144 RelationLookup& lookup, size_t digits)
145{
146 lookup.clear();
147 //If we don't have a relation we don't build a lookuptable
148 if (!relation) return;
149 //Resize to number of digits and set all values
150 lookup.resize(digits);
151 for (const auto& element : relation) {
152 lookup[element.getFromIndex()] = &element;
153 }
154}
155
157 std::map<unsigned int, float>& relation, unsigned int index)
158{
159 //If the lookup table is not empty and the element is set
160 if (!lookup.empty() && lookup[index]) {
161 const RelationElement& element = *lookup[index];
162 const unsigned int size = element.getSize();
163 //Add all Relations to the map
164 for (unsigned int i = 0; i < size; ++i) {
165 //negative weights are from ignored particles, we don't like them and
166 //thus ignore them :D
167 if (element.getWeight(i) < 0) continue;
168 relation[element.getToIndex(i)] += element.getWeight(i);
169 }
170 }
171}
172
174{
175 const StoreArray<SVDRecoDigit> storeRecoDigits(m_storeRecoDigitsName);
176 // If no recodigits, nothing to do
177 if (!storeRecoDigits || !storeRecoDigits.getEntries()) return;
178
179 size_t nDigits = storeRecoDigits.getEntries();
180 B2DEBUG(90, "Initial size of RecoDigits array: " << nDigits);
181
182 const StoreArray<MCParticle> storeMCParticles(m_storeMCParticlesName);
183 const StoreArray<SVDTrueHit> storeTrueHits(m_storeTrueHitsName);
184
185 RelationArray relRecoDigitMCParticle(storeRecoDigits, storeMCParticles, m_relRecoDigitMCParticleName);
186 RelationArray relRecoDigitTrueHit(storeRecoDigits, storeTrueHits, m_relRecoDigitTrueHitName);
187
189 storeClusters.clear();
190
191 RelationArray relClusterMCParticle(storeClusters, storeMCParticles,
193 if (relClusterMCParticle) relClusterMCParticle.clear();
194
195 RelationArray relClusterRecoDigit(storeClusters, storeRecoDigits,
197 if (relClusterRecoDigit) relClusterRecoDigit.clear();
198
199 RelationArray relClusterTrueHit(storeClusters, storeTrueHits,
201 if (relClusterTrueHit) relClusterTrueHit.clear();
202
203 //Build lookup tables for relations
204 createRelationLookup(relRecoDigitMCParticle, m_mcRelation, nDigits);
205 createRelationLookup(relRecoDigitTrueHit, m_trueRelation, nDigits);
206
207 // Create fit tool object
209
210 // I. Group digits by sensor/side.
211 vector<pair<unsigned short, unsigned short> > sensorDigits;
212 VxdID lastSensorID(0);
213 size_t firstSensorDigit = 0;
214 for (size_t iDigit = 0; iDigit < nDigits; ++iDigit) {
215 const SVDRecoDigit& adigit = *storeRecoDigits[iDigit];
216 VxdID sensorID = adigit.getSensorID();
217 sensorID.setSegmentNumber(adigit.isUStrip() ? 1 : 0);
218 if (sensorID != lastSensorID) { // we have a new sensor side
219 sensorDigits.push_back(make_pair(firstSensorDigit, iDigit));
220 firstSensorDigit = iDigit;
221 lastSensorID = sensorID;
222 }
223 }
224 // save last VxdID
225 sensorDigits.push_back(make_pair(firstSensorDigit, nDigits));
226
227 // ICYCLE OVER SENSORS
228 for (auto id_indices : sensorDigits) {
229 // Retrieve parameters from sensorDigits
230 unsigned int firstDigit = id_indices.first;
231 unsigned int lastDigit = id_indices.second;
232 // Get VXDID from the first digit
233 const SVDRecoDigit& sampleRecoDigit = *storeRecoDigits[firstDigit];
234 VxdID sensorID = sampleRecoDigit.getSensorID();
235 bool isU = sampleRecoDigit.isUStrip();
236
237 // Retrieve sensor parameters from GeoCache
238 const SensorInfo& sensorInfo = dynamic_cast<const SensorInfo&>(VXD::GeoCache::getInstance().getSensorInfo(sensorID));
239
240 // 4. Cycle through digits and form clusters on the way.
241
242 // Identify contiguous groups of strips
243 // We maintain two pointers: firstDigit and lastDigit. We also keep track of last seen
244 // strip number.
245 // We start with firstDigit = lastDigit = 0
246 // If there is a gap in strip numbers or the current digit/strip is invalid, we save the
247 // current group (fd, cd), unless empty, and set fd = cd.
248 // If current strip is valid, we include it: ld = cd+1.
249
250 vector<pair<size_t, size_t> > stripGroups;
251 size_t firstClusterDigit = firstDigit;
252 size_t lastClusterDigit = firstDigit;
253 short lastStrip = -2;
254
255 B2DEBUG(300, "Clustering digits " << firstDigit << " to " << lastDigit);
256 for (size_t iDigit = firstDigit; iDigit < lastDigit; ++iDigit) {
257
258 const SVDRecoDigit& recoDigit = *storeRecoDigits[iDigit];
259 unsigned short currentStrip = recoDigit.getCellID();
260 B2DEBUG(300, "Digit " << iDigit << ", strip: " << currentStrip << ", lastStrip: " << lastStrip);
261 B2DEBUG(300, "First CD: " << firstClusterDigit << " Last CD: " << lastClusterDigit);
262
263 // See if we have a gap in strip numbers.
264 bool consecutive = ((currentStrip - lastStrip) == 1);
265 lastStrip = currentStrip;
266
267 B2DEBUG(300, (consecutive ? "consecutive" : "gap"));
268
269 // If we have a gap, we save if there is something to save
270 if (!consecutive && (firstClusterDigit < lastClusterDigit)) {
271 B2DEBUG(300, "Saving (" << firstClusterDigit << ", " << lastClusterDigit << ")");
272 stripGroups.emplace_back(firstClusterDigit, lastClusterDigit);
273 }
274
275 // Update counters
276 if (consecutive) {// consecutive : include in cluster
277 lastClusterDigit = iDigit + 1;
278 } else { // valid after a gap = new cluster
279 firstClusterDigit = iDigit;
280 lastClusterDigit = iDigit + 1;
281 }
282 } // for digits
283 // save last if any
284 if (firstClusterDigit < lastClusterDigit) {
285 B2DEBUG(300, "Saving (" << firstClusterDigit << ", " << lastDigit << ")");
286 stripGroups.emplace_back(firstClusterDigit, lastDigit);
287 }
288
289 ostringstream os;
290 os << "StripGroups: " << endl;
291 for (auto item : stripGroups) {
292 os << "(" << item.first << ", " << item.second << "), ";
293 }
294 B2DEBUG(300, os.str());
295
296 // 5. Now we loop through stripGroups and form clusters.
297 // We are still keeping digit indices, so we are keeping hold on the relations as wwell.
298 double pitch = isU ? sensorInfo.getUPitch() : sensorInfo.getVPitch();
299
300 // FIXME: Pretty stupid. There should be ClusterCandidate type to hold all this, as
301 // well as normed samples. It could calculate a lot of things by itself and make the
302 // code cleaner.
303 vector<unsigned short> stripNumbers;
304 vector<float> stripPositions;
305 vector<float> stripNoises;
306 vector<float> stripGains;
307 vector<float> timeShifts;
308 vector<float> waveWidths;
309 vector<apvSamples> storedNormedSamples;
310 vector<SVDRecoDigit::OutputProbArray> storedPDFs;
311
312 // Now reconstruct clusters
313 for (auto clusterBounds : stripGroups) {
314
315 unsigned short clusterSize = clusterBounds.second - clusterBounds.first;
316 assert(clusterSize > 0);
317
318 stripNumbers.clear();
319 stripPositions.clear();
320 stripNoises.clear();
321 stripGains.clear();
322 timeShifts.clear();
323 waveWidths.clear();
324
325 for (size_t iDigit = clusterBounds.first; iDigit < clusterBounds.second; ++iDigit) {
326 // Now we take the normed samples and calculate the common
327 // cluster probability distribution and, from it, the time and its error, and
328 // the probability of the hit being from beam background.
329 const SVDRecoDigit& recoDigit = *storeRecoDigits[iDigit];
330
331 unsigned short stripNo = recoDigit.getCellID();
332 stripNumbers.push_back(stripNo);
333 // Is the calibrations interface sensible?
334 double stripNoiseADU = m_noiseCal.getNoise(sensorID, isU, stripNo);
335 stripNoises.push_back(
336 m_pulseShapeCal.getChargeFromADC(sensorID, isU, stripNo, stripNoiseADU)
337 );
338 // Some calibrations magic.
339 // FIXME: Only use calibration on real data. Until simulations correspond to
340 // default calibrtion, we cannot use it.
341 double peakWidth = 270;
342 double timeShift = isU ? 4.0 : 0.0;
343 if (m_calibratePeak) {
344 peakWidth = 1.988 * m_pulseShapeCal.getWidth(sensorID, isU, stripNo);
345 timeShift = m_pulseShapeCal.getPeakTime(sensorID, isU, stripNo)
346 - 0.25 * peakWidth;
347 }
348 waveWidths.push_back(peakWidth);
349 timeShifts.push_back(timeShift);
350 stripPositions.push_back(
351 isU ? sensorInfo.getUCellPosition(stripNo) : sensorInfo.getVCellPosition(stripNo)
352 );
353 // Recover samples from ShaperDigits and store them.
354 apvSamples normedSamples;
355 const SVDShaperDigit* shaperDigit = recoDigit.getRelatedTo<SVDShaperDigit>();
356 if (!shaperDigit) // PANIC, this should not happen.
357 B2FATAL("Missing SVDRecoDigits->SVDShaperDigits relation. This should not happen.");
358 auto samples = shaperDigit->getSamples();
359 transform(samples.begin(), samples.end(), normedSamples.begin(),
360 bind(divides<float>(), _1, stripNoiseADU));
361
362 // These are from ShaperDigit, we need to zeroSuppress again,
363 // just in case.
364 zeroSuppress(normedSamples, m_cutAdjacent);
365 storedNormedSamples.emplace_back(normedSamples);
366 storedPDFs.emplace_back(recoDigit.getProbabilities());
367 }
368 // The formulas don't take into account differences in strip noises. So we take RMS
369 // strip noise as cluster noise.
370 // clusterNoise is noise in electrons.
371 float clusterNoise = sqrt(
372 1.0 / clusterSize
373 * inner_product(stripNoises.begin(), stripNoises.end(), stripNoises.begin(), 0.0)
374 );
375 B2DEBUG(200, "RMS cluster noise: " << clusterNoise);
376
377 // This will hold component pdfs. We may want to rememeber them to study
378 // homogeneity of cluster times.
379 shared_ptr<nnFitterBinData> pStrip;
380 // This will aggregate the components pdfs to get cluster time pdf
381 nnFitterBinData pCluster(m_fitter.getBinCenters().size());
382 fill(pCluster.begin(), pCluster.end(), double(1.0));
383
384 for (size_t iClusterStrip = 0; iClusterStrip < clusterSize; ++iClusterStrip) {
385 size_t iDigit = clusterBounds.first + iClusterStrip;
386 ostringstream os1;
387 os1 << "Input to NNFitter: iDigit = " << iDigit << endl << "Samples: ";
388 copy(storedNormedSamples[iClusterStrip].begin(), storedNormedSamples[iClusterStrip].end(),
389 ostream_iterator<double>(os1, " "));
390 os1 << endl;
391 os1 << "PDF from RecoDigit: " << endl;
392 copy(storedPDFs[iClusterStrip].begin(), storedPDFs[iClusterStrip].end(), ostream_iterator<double>(os1, " "));
393 os1 << endl;
394 fitTool.multiply(pCluster, storedPDFs[iClusterStrip]);
395 os1 << "Accummulated: " << endl;
396 copy(pCluster.begin(), pCluster.end(), ostream_iterator<double>(os1, " "));
397 B2DEBUG(200, os1.str());
398 }
399 // We can now calculate clsuter time and its error.
400 double clusterTime, clusterTimeErr;
401 tie(clusterTime, clusterTimeErr) = fitTool.getTimeShift(pCluster);
402 B2DEBUG(200, "Time: " << clusterTime << " +/- " << clusterTimeErr);
403 // Now we have the cluster time pdf, so we can calculate amplitudes.
404 // In the next cycle thrugh cluster's digits, we calculate ampltidues and their
405 // errors.
406 vector<double> stripAmplitudes(stripNoises.size());
407 vector<double> stripAmplitudeErrors(stripNoises.size());
408 double clusterChi2 = 0.0;
409 for (size_t iClusterStrip = 0; iClusterStrip < clusterSize; ++iClusterStrip) {
410 size_t iDigit = clusterBounds.first + iClusterStrip;
411 double snAmp, snAmpError, chi2;
412 tie(snAmp, snAmpError, chi2) =
413 fitTool.getAmplitudeChi2(storedNormedSamples[iClusterStrip], clusterTime, waveWidths[iClusterStrip]);
414 // We have to de-normalize: Multiply amplitudes in SN units by noises in electrons.
415 stripAmplitudes[iClusterStrip] = stripNoises[iClusterStrip] * snAmp;
416 stripAmplitudeErrors[iClusterStrip] = stripNoises[iClusterStrip] * snAmpError;
417 clusterChi2 += chi2;
418 B2DEBUG(200, "Digit " << iDigit << " Noise: " << stripNoises[iClusterStrip]
419 << " Amplitude: " << stripAmplitudes[iClusterStrip]
420 << " +/- " << stripAmplitudeErrors[iClusterStrip]
421 << " Chi2: " << chi2
422 );
423 }
424 // Cluster charge, S/N, chi2
425 float clusterCharge = accumulate(stripAmplitudes.begin(), stripAmplitudes.end(), 0.0);
426 float clusterChargeError = sqrt(
427 inner_product(stripAmplitudeErrors.begin(), stripAmplitudeErrors.end(),
428 stripAmplitudeErrors.begin(), 0.0)
429 );
430 float clusterSN = (clusterChargeError > 0) ? clusterCharge / clusterChargeError : clusterCharge;
431 // FIXME: Return chi2 and ndf from GetApmplitude, this is ugly.
432 clusterChi2 /= clusterSize;
433 // Cluster seed
434 size_t seedIndex = distance(stripAmplitudes.begin(), max_element(
435 stripAmplitudes.begin(), stripAmplitudes.end()));
436 float clusterSeedCharge = stripAmplitudes[seedIndex];
437 B2DEBUG(200, "Cluster parameters:");
438 B2DEBUG(200, "Charge: " << clusterCharge << " +/- " << clusterChargeError);
439 B2DEBUG(200, "Seed: " << clusterSeedCharge << " +/- " << stripAmplitudeErrors[seedIndex]);
440 B2DEBUG(200, "S/N: " << clusterSN);
441 B2DEBUG(200, "chi2: " << clusterChi2);
442
443 // Hit position
444 float clusterPosition, clusterPositionError;
445
446 // Now censoring:
447 // NB: We don't censor strip charges again. The requirement of 3 consecutive samples
448 // above threshold is a good enough safeguard.
449 if ((clusterCharge < clusterNoise * m_cutCluster) || (clusterSeedCharge < clusterNoise * m_cutSeed))
450 continue;
451
452 if (clusterSize < m_sizeHeadTail) { // centre of gravity
453 clusterPosition = 1.0 / clusterCharge * inner_product(
454 stripAmplitudes.begin(), stripAmplitudes.end(), stripPositions.begin(), 0.0
455 );
456 // Strip charge equal to zero-suppression threshold
457 float phantomCharge = m_cutAdjacent * clusterNoise;
458 if (clusterSize == 1) {
459 clusterPositionError = pitch * phantomCharge / (clusterCharge + phantomCharge);
460 } else {
461 clusterPositionError = pitch * phantomCharge / clusterCharge;
462 }
463 } else { // Head-tail
464 float leftStripCharge = stripAmplitudes.front();
465 float leftPos = stripPositions.front();
466 float rightStripCharge = stripAmplitudes.back();
467 float rightPos = stripPositions.back();
468 float centreCharge = (clusterCharge - leftStripCharge - rightStripCharge) / (clusterSize - 2);
469 leftStripCharge = (leftStripCharge < centreCharge) ? leftStripCharge : centreCharge;
470 rightStripCharge = (rightStripCharge < centreCharge) ? rightStripCharge : centreCharge;
471 clusterPosition = 0.5 * (leftPos + rightPos)
472 + 0.5 * (rightStripCharge - leftStripCharge) / centreCharge * pitch;
473 float sn = centreCharge / m_cutAdjacent / clusterNoise;
474 // Rough estimates of Landau noise
475 float landauHead = leftStripCharge / centreCharge;
476 double landauTail = rightStripCharge / centreCharge;
477 clusterPositionError = 0.5 * pitch * sqrt(1.0 / sn / sn
478 + 0.5 * landauHead * landauHead + 0.5 * landauTail * landauTail);
479 }
480 //Lorentz shift correction
481 clusterPosition -= sensorInfo.getLorentzShift(isU, clusterPosition);
482
483 // Finally, we save the cluster and its relations.
484 map<unsigned int, float> mc_relations;
485 map<unsigned int, float> truehit_relations;
486 vector<pair<unsigned int, float> > digit_weights;
487 digit_weights.reserve(clusterSize);
488
489 for (size_t iDigit = clusterBounds.first; iDigit < clusterBounds.second; ++iDigit) {
490 // Store relations to MCParticles and SVDTrueHits
491 fillRelationMap(m_mcRelation, mc_relations, iDigit);
492 fillRelationMap(m_trueRelation, truehit_relations, iDigit);
493 //Add digit to the Cluster->Digit relation list
494 digit_weights.emplace_back(iDigit, stripAmplitudes[iDigit - clusterBounds.first]);
495 }
496
497 //Store the cluster into Datastore ...
498 int clsIndex = storeClusters.getEntries();
499 VxdID clusterSensorID = sensorID;
500 clusterSensorID.setSegmentNumber(0);
501 storeClusters.appendNew(
502 SVDCluster(sensorID, isU, clusterPosition, clusterPositionError, clusterTime,
503 clusterTimeErr, clusterCharge, clusterSeedCharge, clusterSize, clusterSN, clusterChi2)
504 );
505
506 //Create relations to the cluster
507 if (!mc_relations.empty()) {
508 relClusterMCParticle.add(clsIndex, mc_relations.begin(), mc_relations.end());
509 }
510 if (!truehit_relations.empty()) {
511 relClusterTrueHit.add(clsIndex, truehit_relations.begin(), truehit_relations.end());
512 }
513 relClusterRecoDigit.add(clsIndex, digit_weights.begin(), digit_weights.end());
514 } // CYCLE OVER CLUSTERS for stripGroups
515
516 } // CYCLE OVER SENSORS for items in sensorDigits
517
518 B2DEBUG(100, "Number of clusters: " << storeClusters.getEntries());
519
520} // event()
521
522
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.
TO * getRelatedTo(const std::string &name="", const std::string &namedRelation="") const
Get the object to which this object has a relation.
The SVD Cluster class This class stores all information about reconstructed SVD clusters.
Definition: SVDCluster.h:29
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 RecoDigit class.
Definition: SVDRecoDigit.h:43
OutputProbArray getProbabilities() const
Get signal time pdf.
Definition: SVDRecoDigit.h:144
VxdID getSensorID() const
Get the sensor ID.
Definition: SVDRecoDigit.h:98
short int getCellID() const
Get strip ID.
Definition: SVDRecoDigit.h:114
bool isUStrip() const
Get strip direction.
Definition: SVDRecoDigit.h:109
The SVD ShaperDigit class.
APVFloatSamples getSamples() const
Get array of samples.
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 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::string m_storeRecoDigitsName
Name of the collection to use for the SVDRecoDigits.
virtual void initialize() override
Initialize the module.
virtual void event() override
do the clustering
std::vector< const RelationElement * > RelationLookup
Container for a RelationArray Lookup table.
std::string m_relRecoDigitTrueHitName
Name of the relation between SVDRecoDigits and SVDTrueHits.
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_relRecoDigitMCParticleName
Name of the relation between SVDRecoDigits and MCParticles.
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 SVDRecoDigit index to a map.
std::string m_storeMCParticlesName
Name of the collection to use for the MCParticles.
RelationLookup m_trueRelation
Lookup table for SVDRecoDigit->SVDTrueHit relation.
SVDPulseShapeCalibrations m_pulseShapeCal
Calibrations: pusle shape and gain.
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_relClusterMCParticleName
Name of the relation between SVDClusters and MCParticles.
SVDNNClusterizerModule()
Constructor defining the parameters.
RelationLookup m_mcRelation
Lookup table for SVDRecoDigit->MCParticle relation.
std::string m_relClusterRecoDigitName
Name of the relation between SVDClusters and SVDRecoDigits.
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 ROOT::Math::XYZVector & getLorentzShift(double uCoord, double vCoord) const
Calculate Lorentz shift along a given coordinate in a magnetic field at a given position.
Definition: SensorInfo.cc:104
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
double getVCellPosition(int vID) const
Return the position of a specific strip/pixel in v direction.
double getUPitch(double v=0) const
Return the pitch of the sensor.
double getUCellPosition(int uID, int vID=-1) const
Return the position of a specific strip/pixel in u direction.
double getVPitch(double v=0) const
Return the pitch of the sensor.
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
Abstract base class for different kinds of events.
STL namespace.