9#include <svd/modules/svdReconstruction/SVDClusterizerDirectModule.h>
11#include <framework/datastore/StoreArray.h>
12#include <framework/logging/Logger.h>
14#include <vxd/geometry/GeoCache.h>
15#include <svd/geometry/SensorInfo.h>
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>
23#include <svd/reconstruction/NNWaveFitTool.h>
25#include <unordered_map>
32using namespace std::placeholders;
47 B2DEBUG(200,
"SVDClusterizerDirectModule ctor");
54 "6-digits collection name",
string(
""));
56 "Cluster collection name",
string(
""));
58 "TrueHit collection name",
string(
""));
60 "MCParticles collection name",
string(
""));
62 "SVDEventInfo name",
string(
""));
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));
72 "SN for digits to be considered for clustering",
m_cutAdjacent);
74 "SN for digits to be considered as seed",
m_cutSeed);
78 "Cluster size at which to switch to Analog head tail algorithm",
m_sizeHeadTail);
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);
124 B2INFO(
" 1. COLLECTIONS:");
134 B2INFO(
" 2. CALIBRATION DATA:");
136 B2INFO(
" 4. CLUSTERING:");
154 if (!relation)
return;
156 lookup.resize(digits);
157 for (
const auto& element : relation) {
158 lookup[element.getFromIndex()] = &element;
163 std::map<unsigned int, float>& relation,
unsigned int index)
166 if (!lookup.empty() && lookup[index]) {
168 const unsigned int size = element.getSize();
170 for (
unsigned int i = 0; i < size; ++i) {
173 if (element.getWeight(i) < 0)
continue;
174 relation[element.getToIndex(i)] += element.getWeight(i);
188 size_t nDigits = storeShaperDigits.
getEntries();
189 B2DEBUG(90,
"Initial size of StoreDigits array: " << nDigits);
198 storeClusters.
clear();
200 RelationArray relClusterMCParticle(storeClusters, storeMCParticles,
202 if (relClusterMCParticle) relClusterMCParticle.
clear();
204 RelationArray relClusterShaperDigit(storeClusters, storeShaperDigits,
206 if (relClusterShaperDigit) relClusterShaperDigit.
clear();
208 RelationArray relClusterTrueHit(storeClusters, storeTrueHits,
210 if (relClusterTrueHit) relClusterTrueHit.
clear();
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) {
227 if (sensorID != lastSensorID) {
228 sensorDigits.push_back(make_pair(firstSensorDigit, iDigit));
229 firstSensorDigit = iDigit;
230 lastSensorID = sensorID;
234 sensorDigits.push_back(make_pair(firstSensorDigit, nDigits));
237 for (
auto id_indices : sensorDigits) {
239 unsigned int firstDigit = id_indices.first;
240 unsigned int lastDigit = id_indices.second;
242 const SVDShaperDigit& sampleDigit = *storeShaperDigits[firstDigit];
248 double pitch = isU ? info.getUPitch() : info.getVPitch();
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;
266 B2DEBUG(300,
"Clustering digits " << firstDigit <<
" to " << lastDigit);
267 for (
size_t iDigit = firstDigit; iDigit < lastDigit; ++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);
280 transform(samples.begin(), samples.end(), normedSamples.begin(),
281 bind(divides<float>(), _1, stripNoiseADU));
288 storedNormedSamples.insert(make_pair(iDigit, normedSamples));
292 bool consecutive = ((currentStrip - lastStrip) == 1);
293 lastStrip = currentStrip;
295 B2DEBUG(300, (validDigit ?
"Valid " :
"Invalid ") << (consecutive ?
"consecutive" :
"gap"));
298 if ((!validDigit || !consecutive) && (firstClusterDigit < lastClusterDigit)) {
299 B2DEBUG(300,
"Saving (" << firstClusterDigit <<
", " << lastClusterDigit <<
")");
300 stripGroups.emplace_back(firstClusterDigit, lastClusterDigit);
306 lastClusterDigit = iDigit + 1;
308 firstClusterDigit = iDigit;
309 lastClusterDigit = iDigit + 1;
312 firstClusterDigit = iDigit + 1;
313 lastClusterDigit = iDigit + 1;
317 if (firstClusterDigit < lastClusterDigit) {
318 B2DEBUG(300,
"Saving (" << firstClusterDigit <<
", " << lastDigit <<
")");
319 stripGroups.emplace_back(firstClusterDigit, lastDigit);
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,
" "));
329 os <<
"StripGroups: " << endl;
330 for (
auto item : stripGroups) {
331 os <<
"(" << item.first <<
", " << item.second <<
"), ";
333 B2DEBUG(300, os.str());
341 vector<unsigned short> stripNumbers;
342 vector<float> stripPositions;
343 vector<float> stripNoises;
344 vector<float> timeShifts;
345 vector<float> waveWidths;
347 for (
auto clusterBounds : stripGroups) {
349 unsigned short clusterSize = clusterBounds.second - clusterBounds.first;
350 assert(clusterSize > 0);
352 stripNumbers.clear();
353 stripPositions.clear();
358 for (
size_t iDigit = clusterBounds.first; iDigit < clusterBounds.second; ++iDigit) {
364 unsigned short stripNo = digit.
getCellID();
365 stripNumbers.push_back(stripNo);
368 stripNoises.push_back(
374 double peakWidth = 270;
375 double timeShift = isU ? 2.5 : -2.2;
382 const double triggerBinSep = 4 * 1.96516;
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)
390 stripPositions.push_back(
391 isU ? info.getUCellPosition(stripNo) : info.getVCellPosition(stripNo)
397 float clusterNoise =
sqrt(
399 * inner_product(stripNoises.begin(), stripNoises.end(), stripNoises.begin(), 0.0)
401 B2DEBUG(200,
"RMS cluster noise: " << clusterNoise);
405 shared_ptr<nnFitterBinData> pStrip;
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;
412 os1 <<
"Input to NNFitter: iDigit = " << iDigit << endl <<
"Samples: ";
413 copy(storedNormedSamples[iDigit].begin(), storedNormedSamples[iDigit].end(),
414 ostream_iterator<double>(os1,
" "));
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,
" "));
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());
428 double clusterTime, clusterTimeErr;
429 tie(clusterTime, clusterTimeErr) = fitTool.
getTimeShift(pCluster);
430 B2DEBUG(200,
"Time: " << clusterTime <<
" +/- " << clusterTimeErr);
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]);
443 stripAmplitudes[iClusterStrip] =
444 stripNoises[iClusterStrip] * snAmp;
445 stripAmplitudeErrors[iClusterStrip] =
446 stripNoises[iClusterStrip] * snAmpError;
448 B2DEBUG(200,
"Digit " << iDigit <<
" Noise: " << stripNoises[iClusterStrip]
449 <<
" Amplitude: " << stripAmplitudes[iClusterStrip]
450 <<
" +/- " << stripAmplitudeErrors[iClusterStrip]
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)
460 float clusterSN = (clusterChargeError > 0) ? clusterCharge / clusterChargeError : clusterCharge;
461 clusterChi2 /= clusterSize;
463 size_t seedIndex = distance(stripAmplitudes.begin(), max_element(
464 stripAmplitudes.begin(), stripAmplitudes.end()));
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);
474 float clusterPosition, clusterPositionError;
483 clusterPosition = 1.0 / clusterCharge * inner_product(
484 stripAmplitudes.begin(), stripAmplitudes.end(), stripPositions.begin(), 0.0
488 if (clusterSize == 1) {
489 clusterPositionError = pitch * phantomCharge / (clusterCharge + phantomCharge);
491 clusterPositionError = pitch * phantomCharge / clusterCharge;
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;
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);
511 clusterPosition -= info.getLorentzShift(isU, clusterPosition);
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);
519 for (
size_t iDigit = clusterBounds.first; iDigit < clusterBounds.second; ++iDigit) {
524 digit_weights.emplace_back(iDigit, stripAmplitudes[iDigit - clusterBounds.first]);
529 VxdID clusterSensorID = sensorID;
532 SVDCluster(sensorID, isU, clusterPosition, clusterPositionError, clusterTime,
533 clusterTimeErr, clusterCharge, clusterSeedCharge, clusterSize, clusterSN, clusterChi2)
537 if (!mc_relations.empty()) {
538 relClusterMCParticle.
add(clsIndex, mc_relations.begin(), mc_relations.end());
540 if (!truehit_relations.empty()) {
541 relClusterTrueHit.
add(clsIndex, truehit_relations.begin(), truehit_relations.end());
543 relClusterShaperDigit.
add(clsIndex, digit_weights.begin(), digit_weights.end());
548 B2DEBUG(100,
"Number of clusters: " << storeClusters.
getEntries());
Class for accessing objects in the database.
void setDescription(const std::string &description)
Sets the description of the module.
void setPropertyFlags(unsigned int propertyFlags)
Sets the flags for the module properties.
@ c_ParallelProcessingCertified
This module can be run in parallel processing mode safely (All I/O must be done through the data stor...
Low-level class to create/modify relations between StoreArrays.
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.
Class to store SVD mode information.
baseType getTriggerBin() const
Get the triggerBin id.
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.
void setNetwrok(const std::string &xmlData)
Set proper network definition file.
const NNWaveFitTool & getFitTool() const
Get a handle to a NNWaveFit object.
const nnFitterBinData & getBinCenters() const
Get bin times of the network output.
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.
NNWaveFitter m_fitter
Time fitter.
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...
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.
T * appendNew()
Construct a new T object at the end of the array.
int getEntries() const
Get the number of objects in the array.
void clear() override
Delete all entries in this array.
const SensorInfoBase & getSensorInfo(Belle2::VxdID id) const
Return a referecne to the SensorInfo of a given SensorID.
static GeoCache & getInstance()
Return a reference to the singleton instance.
Class to uniquely identify a any structure of the PXD and SVD.
void setSegmentNumber(baseType segment)
Set the sensor segment.
void addParam(const std::string &name, T ¶mVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
double sqrt(double a)
sqrt for double
Namespace to encapsulate code needed for simulation and reconstrucion of the SVD.
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.
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.