11 #include <svd/modules/svdReconstruction/SVDSimpleClusterizerModule.h>
13 #include <framework/datastore/DataStore.h>
14 #include <framework/datastore/RelationArray.h>
15 #include <framework/datastore/RelationIndex.h>
16 #include <framework/logging/Logger.h>
18 #include <svd/geometry/SensorInfo.h>
19 #include <svd/dataobjects/SVDEventInfo.h>
36 m_cutSeed(5.0), m_cutAdjacent(3.0), m_sizeHeadTail(3), m_cutCluster(0), m_useDB(true)
39 setDescription(
"Clusterize SVDRecoDigits fitted by the Center of Gravity estimator");
40 setPropertyFlags(c_ParallelProcessingCertified);
43 addParam(
"RecoDigits", m_storeRecoDigitsName,
44 "SVDRecoDigits collection name",
string(
""));
45 addParam(
"Clusters", m_storeClustersName,
46 "SVDCluster collection name",
string(
""));
47 addParam(
"SVDTrueHits", m_storeTrueHitsName,
48 "TrueHit collection name",
string(
""));
49 addParam(
"MCParticles", m_storeMCParticlesName,
50 "MCParticles collection name",
string(
""));
51 addParam(
"ShaperDigits", m_storeShaperDigitsName,
52 "SVDShaperDigits collection name",
56 addParam(
"AdjacentSN", m_cutAdjacent,
57 "SN for digits to be considered for clustering", m_cutAdjacent);
58 addParam(
"SeedSN", m_cutSeed,
59 "SN for digits to be considered as seed", m_cutSeed);
60 addParam(
"HeadTailSize", m_sizeHeadTail,
61 "Cluster size at which to switch to Analog head tail algorithm", m_sizeHeadTail);
62 addParam(
"ClusterSN", m_cutCluster,
63 "minimum value of the SNR of the cluster", m_cutCluster);
64 addParam(
"timeAlgorithm", m_timeAlgorithm,
65 " int to choose time algorithm: 0 = 6-sample CoG (default for 6-sample acquisition mode), 1 = 3-sample CoG (default for 3-sample acquisition mode), 2 = 3-sample ELS",
67 addParam(
"Calibrate3SampleWithEventT0", m_calibrate3SampleWithEventT0,
68 " if true returns the calibrated time instead of the raw time for 3-sample time algorithms",
69 m_calibrate3SampleWithEventT0);
70 addParam(
"useDB", m_useDB,
71 "if false use clustering module parameters", m_useDB);
76 void SVDSimpleClusterizerModule::initialize()
79 m_storeClusters.registerInDataStore(m_storeClustersName, DataStore::c_ErrorIfAlreadyRegistered);
80 m_storeDigits.isRequired(m_storeRecoDigitsName);
81 m_storeTrueHits.isOptional(m_storeTrueHitsName);
82 m_storeMCParticles.isOptional(m_storeMCParticlesName);
84 RelationArray relClusterDigits(m_storeClusters, m_storeDigits);
85 RelationArray relClusterTrueHits(m_storeClusters, m_storeTrueHits);
86 RelationArray relClusterMCParticles(m_storeClusters, m_storeMCParticles);
87 RelationArray relDigitTrueHits(m_storeDigits, m_storeTrueHits);
88 RelationArray relDigitMCParticles(m_storeDigits, m_storeMCParticles);
98 m_storeClustersName = m_storeClusters.getName();
99 m_storeRecoDigitsName = m_storeDigits.getName();
100 m_storeTrueHitsName = m_storeTrueHits.getName();
101 m_storeMCParticlesName = m_storeMCParticles.getName();
103 m_relClusterRecoDigitName = relClusterDigits.
getName();
104 m_relClusterTrueHitName = relClusterTrueHits.
getName();
105 m_relClusterMCParticleName = relClusterMCParticles.
getName();
106 m_relRecoDigitTrueHitName = relDigitTrueHits.
getName();
107 m_relRecoDigitMCParticleName = relDigitMCParticles.
getName();
110 B2DEBUG(1,
"SVDSimpleClusterizer Parameters (in default system unit, *=cannot be set directly):");
112 B2DEBUG(1,
" 1. COLLECTIONS:");
113 B2DEBUG(1,
" --> MCParticles: " << DataStore::arrayName<MCParticle>(m_storeMCParticlesName));
114 B2DEBUG(1,
" --> SVDRecoDigits: " << DataStore::arrayName<SVDRecoDigit>(m_storeRecoDigitsName));
115 B2DEBUG(1,
" --> SVDClusters: " << DataStore::arrayName<SVDCluster>(m_storeClustersName));
116 B2DEBUG(1,
" --> SVDTrueHits: " << DataStore::arrayName<SVDTrueHit>(m_storeTrueHitsName));
117 B2DEBUG(1,
" --> DigitMCRel: " << m_relRecoDigitMCParticleName);
118 B2DEBUG(1,
" --> ClusterMCRel: " << m_relClusterMCParticleName);
119 B2DEBUG(1,
" --> ClusterDigitRel: " << m_relClusterRecoDigitName);
120 B2DEBUG(1,
" --> DigitTrueRel: " << m_relRecoDigitTrueHitName);
121 B2DEBUG(1,
" --> ClusterTrueRel: " << m_relClusterTrueHitName);
122 B2DEBUG(1,
" 2. CLUSTERING:");
123 B2DEBUG(1,
" --> Neighbour cut: " << m_cutAdjacent);
124 B2DEBUG(1,
" --> Seed cut: " << m_cutSeed);
125 B2DEBUG(1,
" --> Size HeadTail: " << m_sizeHeadTail);
130 void SVDSimpleClusterizerModule::event()
132 int nDigits = m_storeDigits.getEntries();
136 m_storeClusters.clear();
138 RelationArray relClusterMCParticle(m_storeClusters, m_storeMCParticles,
139 m_relClusterMCParticleName);
140 if (relClusterMCParticle) relClusterMCParticle.
clear();
142 RelationArray relClusterDigit(m_storeClusters, m_storeDigits,
143 m_relClusterRecoDigitName);
144 if (relClusterDigit) relClusterDigit.
clear();
146 RelationArray relClusterTrueHit(m_storeClusters, m_storeTrueHits,
147 m_relClusterTrueHitName);
148 if (relClusterTrueHit) relClusterTrueHit.
clear();
152 m_cutSeed = m_ClusterCal.getMinSeedSNR(m_storeDigits[0]->getSensorID(), m_storeDigits[0]->isUStrip());
153 m_cutAdjacent = m_ClusterCal.getMinAdjSNR(m_storeDigits[0]->getSensorID(), m_storeDigits[0]->isUStrip());
154 m_cutCluster = m_ClusterCal.getMinClusterSNR(m_storeDigits[0]->getSensorID(), m_storeDigits[0]->isUStrip());
158 m_sizeHeadTail, m_cutSeed, m_cutAdjacent, m_cutCluster, m_timeAlgorithm, m_storeShaperDigitsName, m_storeRecoDigitsName);
162 while (i < nDigits) {
165 VxdID thisSensorID = m_storeDigits[i]->getSensorID();
166 bool thisSide = m_storeDigits[i]->isUStrip();
167 int thisCellID = m_storeDigits[i]->getCellID();
170 m_cutSeed = m_ClusterCal.getMinSeedSNR(thisSensorID, thisSide);
171 m_cutAdjacent = m_ClusterCal.getMinAdjSNR(thisSensorID, thisSide);
172 m_cutCluster = m_ClusterCal.getMinClusterSNR(thisSensorID, thisSide);
176 float thisNoise = m_NoiseCal.getNoiseInElectrons(thisSensorID, thisSide, thisCellID);
177 float thisCharge = m_storeDigits[i]->getCharge();
178 B2DEBUG(10,
"Noise = " << thisNoise <<
" e-, Charge = " << thisCharge);
180 if ((
float)thisCharge / thisNoise < m_cutAdjacent) {
188 aStrip.
charge = thisCharge;
189 aStrip.
cellID = thisCellID;
190 aStrip.
noise = thisNoise;
192 aStrip.
time = m_storeDigits[i]->getTime();
193 aStrip.
timeError = m_storeDigits[i]->getTimeError();
196 if (! clusterCandidate.
add(thisSensorID, thisSide, aStrip)) {
199 if (clusterCandidate.
size() > 0) {
202 writeClusters(clusterCandidate);
207 clusterCandidate =
SimpleClusterCandidate(thisSensorID, thisSide, m_sizeHeadTail, m_cutSeed, m_cutAdjacent, m_cutCluster,
209 m_storeShaperDigitsName, m_storeRecoDigitsName);
212 if (! clusterCandidate.
add(thisSensorID, thisSide, aStrip))
213 B2WARNING(
"this state is forbidden!!");
220 if (clusterCandidate.
size() > 0) {
223 writeClusters(clusterCandidate);
226 B2DEBUG(1,
"Number of clusters: " << m_storeClusters.getEntries());
233 RelationArray relClusterDigit(m_storeClusters, m_storeDigits, m_relClusterRecoDigitName);
235 RelationArray relClusterMCParticle(m_storeClusters, m_storeMCParticles, m_relClusterMCParticleName);
236 RelationArray relClusterTrueHit(m_storeClusters, m_storeTrueHits, m_relClusterTrueHitName);
242 VxdID sensorID = cluster.getSensorID();
243 bool isU = cluster.isUSide();
244 float seedCharge = cluster.getSeedCharge();
245 float charge = cluster.getCharge();
246 float size = cluster.size();
247 float SNR = cluster.getSNR();
248 float position = cluster.getPosition();
249 float positionError = m_ClusterCal.getCorrectedClusterPositionError(sensorID, isU, size, cluster.getPositionError());
251 float time = cluster.getTime();
252 float timeError = cluster.getTimeError();
253 int firstFrame = cluster.getFirstFrame();
257 std::string m_svdEventInfoName =
"SVDEventInfo";
259 m_svdEventInfoName =
"SVDEventInfoSim";
261 if (!eventinfo) B2ERROR(
"No SVDEventInfo!");
268 float caltime = time;
269 if (m_timeAlgorithm == 1 and m_calibrate3SampleWithEventT0)
270 caltime = m_3CoGTimeCal.getCorrectedTime(sensorID, isU, -1, time, -1);
271 else if (m_timeAlgorithm == 2 and m_calibrate3SampleWithEventT0)
272 caltime = m_3ELSTimeCal.getCorrectedTime(sensorID, isU, -1, time, -1);
278 time = eventinfo->getTimeInFTSWReference(caltime, firstFrame);
282 sensorID, isU, position, positionError, time, timeError, charge, seedCharge, size, SNR, -1, firstFrame
286 int clsIndex = m_storeClusters.getEntries() - 1;
288 map<int, float> mc_relations;
289 map<int, float> truehit_relations;
291 vector<pair<int, float> > digit_weights;
292 digit_weights.reserve(size);
294 std::vector<stripInCluster> strips = cluster.getStripsInCluster();
296 for (
auto strip : strips) {
299 if (relDigitMCParticle) {
301 for (relMC_type& mcRel : relDigitMCParticle.
getElementsFrom(m_storeDigits[strip.recoDigitIndex])) {
304 if (mcRel.weight < 0)
continue;
305 mc_relations[mcRel.indexTo] += mcRel.weight;
309 if (relDigitTrueHit) {
311 for (relTrueHit_type& trueRel : relDigitTrueHit.
getElementsFrom(m_storeDigits[strip.recoDigitIndex])) {
314 if (trueRel.weight < 0)
continue;
315 truehit_relations[trueRel.indexTo] += trueRel.weight;
319 digit_weights.push_back(make_pair(strip.recoDigitIndex, strip.charge));
324 if (!mc_relations.empty()) {
325 relClusterMCParticle.
add(clsIndex, mc_relations.begin(), mc_relations.end());
327 if (!truehit_relations.empty()) {
328 relClusterTrueHit.
add(clsIndex, truehit_relations.begin(), truehit_relations.end());
331 relClusterDigit.
add(clsIndex, digit_weights.begin(), digit_weights.end());