11 #include <pxd/modules/pxdReconstruction/PXDClusterizerModule.h>
13 #include <framework/datastore/DataStore.h>
14 #include <framework/datastore/StoreArray.h>
15 #include <framework/datastore/RelationArray.h>
16 #include <framework/logging/Logger.h>
18 #include <vxd/geometry/GeoCache.h>
20 #include <mdst/dataobjects/MCParticle.h>
21 #include <pxd/dataobjects/PXDDigit.h>
22 #include <pxd/dataobjects/PXDCluster.h>
23 #include <pxd/dataobjects/PXDTrueHit.h>
24 #include <pxd/geometry/SensorInfo.h>
26 #include <pxd/reconstruction/PXDClusterPositionEstimator.h>
41 PXDClusterizerModule::PXDClusterizerModule() :
Module()
42 , m_elNoise(0.7), m_cutSeed(5.0), m_cutAdjacent(3.0), m_cutCluster(8.0)
43 , m_cutAdjacentSignal(0), m_sizeHeadTail(3), m_clusterCacheSize(0)
50 "Noise added by the electronics, set in ADU",
m_elNoise);
52 "SN for digits to be considered for clustering",
m_cutAdjacent);
57 "Maximum desired number of sensor rows", 0);
59 "Minimum cluster size to switch to Analog head tail algorithm for cluster center",
80 storeDigits.isRequired();
81 storeTrueHits.isOptional();
82 storeMCParticles.isOptional();
85 RelationArray relClusterTrueHits(storeClusters, storeTrueHits);
86 RelationArray relClusterMCParticles(storeClusters, storeMCParticles);
87 RelationArray relDigitMCParticles(storeDigits, storeMCParticles);
113 "PXDClusterizer Parameters (in default system units, *=cannot be set directly):");
114 B2DEBUG(20,
" --> ElectronicNoise: " <<
m_elNoise);
116 B2DEBUG(20,
" --> SeedSN: " <<
m_cutSeed);
148 storeClusters.
clear();
150 RelationArray relClusterMCParticle(storeClusters, storeMCParticles,
152 if (relClusterMCParticle) relClusterMCParticle.
clear();
156 if (relClusterDigit) relClusterDigit.
clear();
158 RelationArray relClusterTrueHit(storeClusters, storeTrueHits,
160 if (relClusterTrueHit) relClusterTrueHit.
clear();
179 for (
int i = 0; i < nPixels; i++) {
180 const PXDDigit*
const storeDigit = storeDigits[i];
181 Pixel px(storeDigit, i);
188 }
else if (px <= lastPixel) {
190 B2FATAL(
"Pixels are not sorted correctly, please include the "
191 "PXDDigitSorter module before running the Clusterizer or fix "
192 "the input to be ordered by v,u in ascending order");
220 }
catch (std::out_of_range& e) {
221 B2WARNING(
"PXD clustering: Ignoring pixel " << px.
getU() <<
"," << px.
getV() <<
": " << e.what());
231 if (!relation)
return;
233 lookup.resize(digits);
235 lookup[element.getFromIndex()] = &element;
243 if (!lookup.empty() && lookup[index]) {
245 const unsigned int size = element.getSize();
247 for (
unsigned int i = 0; i < size; ++i) {
250 if (element.getWeight(i) < 0)
continue;
251 relation[element.getToIndex(i)] += element.getWeight(i);
266 RelationArray relClusterMCParticle(storeClusters, storeMCParticles,
270 RelationArray relClusterTrueHit(storeClusters, storeTrueHits,
277 map<unsigned int, float> mc_relations;
278 map<unsigned int, float> truehit_relations;
279 vector<pair<unsigned int, float> > digit_weights;
288 mc_relations.clear();
289 truehit_relations.clear();
290 digit_weights.clear();
291 digit_weights.reserve(cls.size());
293 const Pixel& seed = cls.getSeed();
297 projU.
add(px.getU(), info.getUCellPosition(px.getU()), px.getCharge());
298 projV.
add(px.getV(), info.getVCellPosition(px.getV()), px.getCharge());
305 digit_weights.emplace_back(px.getIndex(), px.getCharge());
310 const double pitchU = info.getUPitch();
311 const double pitchV = info.getVPitch(projV.
getPos());
315 double posUU = cls.getCharge() * pitchU * pitchU / 12.0;
316 double posVV = cls.getCharge() * pitchV * pitchV / 12.0;
318 for (
const Pixel& px : cls.pixels()) {
319 const double du = info.getUCellPosition(px.getU()) - projU.
getPos();
320 const double dv = info.getVCellPosition(px.getV()) - projV.
getPos();
321 posUU += px.getCharge() * du * du;
322 posVV += px.getCharge() * dv * dv;
323 posUV += px.getCharge() * du * dv;
325 rho = posUV / sqrt(posUU * posVV);
333 TVector3 lorentzShift = info.getLorentzShift(projU.
getPos(), projV.
getPos());
336 B2DEBUG(20,
"Lorentz shift: " << lorentzShift.X() <<
" " << lorentzShift.Y());
343 set<Pixel> pixelSet(cls.pixels().begin(), cls.pixels().end());
346 vector<float> sectorEtaValues = {0, 0, 0, 0};
356 vector<int> sectorShapeIndices = { -1, -1, -1, -1};
369 rho, cls.getCharge(), seed.getCharge(),
371 sectorEtaValues, sectorShapeIndices
375 if (!mc_relations.empty()) relClusterMCParticle.
add(clsIndex, mc_relations.begin(), mc_relations.end());
376 if (!truehit_relations.empty()) relClusterTrueHit.
add(clsIndex, truehit_relations.begin(), truehit_relations.end());
377 relClusterDigit.
add(clsIndex, digit_weights.begin(), digit_weights.end());
384 const ClusterProjection& secondary,
double minPitch,
double centerPitch,
double maxPitch)
392 + (maxCharge * maxPitch - minCharge * minPitch) / centerCharge));
395 const double landauHead = minCharge / centerCharge * minPitch;
396 const double landauTail = maxCharge / centerCharge * maxPitch;
397 primary.
setError(0.5 * sqrt(1.0 / snHead / snHead + 1.0 / snTail / snTail
398 + 0.5 * landauHead * landauHead + 0.5 * landauTail * landauTail));
399 }
else if (primary.
getSize() <= 2) {
404 primary.
setError(2.0 * centerPitch / sn);