257 std::vector<SiPixelCluster>
output;
259 unsigned int meanExp = nSplitted;
261 output.push_back(aCluster);
265 std::vector<float> clx(meanExp);
266 std::vector<float> cly(meanExp);
267 std::vector<float> cls(meanExp);
268 std::vector<float> oldclx(meanExp);
269 std::vector<float> oldcly(meanExp);
270 std::vector<SiPixelCluster::Pixel> originalpixels = aCluster.
pixels();
271 std::vector<std::pair<int, SiPixelCluster::Pixel>> pixels;
272 for (
unsigned int j = 0;
j < originalpixels.size();
j++) {
276 int perDiv = originalpixels[
j].adc / sub;
278 std::cout <<
"Splitting " <<
j <<
" in [ " << pixels.size() <<
" , " << pixels.size() + sub
279 <<
" ], expected numb of clusters: " << meanExp <<
" original pixel (x,y) " << originalpixels[
j].x
280 <<
" " << originalpixels[
j].y <<
" sub " << sub << std::endl;
281 for (
int k = 0;
k < sub;
k++) {
283 perDiv = originalpixels[
j].adc - perDiv *
k;
287 std::vector<int> clusterForPixel(pixels.size());
289 for (
unsigned int j = 0;
j < meanExp;
j++) {
292 clx[
j] = originalpixels[0].x +
j;
293 cly[
j] = originalpixels[0].y +
j;
297 int remainingSteps = 100;
298 while (!stop && remainingSteps > 0) {
301 std::vector<std::vector<float>> distanceMapX(originalpixels.size(), std::vector<float>(meanExp));
302 std::vector<std::vector<float>> distanceMapY(originalpixels.size(), std::vector<float>(meanExp));
303 std::vector<std::vector<float>> distanceMap(originalpixels.size(), std::vector<float>(meanExp));
304 for (
unsigned int j = 0;
j < originalpixels.size();
j++) {
306 std::cout <<
"Original Pixel pos " <<
j <<
" " << pixels[
j].second.x <<
" " << pixels[
j].second.y << std::endl;
307 for (
unsigned int i = 0;
i < meanExp;
i++) {
308 distanceMapX[
j][
i] = 1.f * originalpixels[
j].x - clx[
i];
309 distanceMapY[
j][
i] = 1.f * originalpixels[
j].y - cly[
i];
316 dist += (2.f * distanceMapX[
j][
i] /
sizeX) * (2.
f * distanceMapX[
j][i] /
sizeX);
320 dist += 1.f * (
std::abs(distanceMapY[
j][i]) -
sizeY / 2.f + 1.f) *
323 dist += 1.f * (2.f * distanceMapY[
j][
i] /
sizeY) * (2.
f * distanceMapY[
j][i] /
sizeY);
325 distanceMap[
j][
i] =
sqrt(dist);
327 std::cout <<
"Cluster " << i <<
" Original Pixel " <<
j <<
" distances: " << distanceMapX[
j][
i] <<
" " 328 << distanceMapY[
j][
i] <<
" " << distanceMap[
j][
i] << std::endl;
334 std::multimap<float, int> scores;
341 std::vector<float> weightOfPixel(pixels.size());
342 for (std::multimap<float, int>::iterator it = scores.begin(); it != scores.end(); it++) {
343 int pixel_index = it->second;
345 std::cout <<
"Original Pixel " << pixel_index <<
" with score " << it->first << std::endl;
347 int subpixel_counter = 0;
348 for (
auto subpixel = pixels.begin(); subpixel != pixels.end(); ++subpixel, ++subpixel_counter) {
349 if (subpixel->first > pixel_index) {
351 }
else if (subpixel->first != pixel_index) {
356 for (
unsigned int subcluster_index = 0; subcluster_index < meanExp; subcluster_index++) {
357 float nsig = (cls[subcluster_index] - expectedADC) /
359 float clQest = 1.f / (1.f +
std::exp(nsig)) + 1
e-6
f;
360 float clDest = 1.f / (distanceMap[pixel_index][subcluster_index] + 0.05f);
363 std::cout <<
" Q: " << clQest <<
" D: " << clDest <<
" " << distanceMap[pixel_index][subcluster_index]
365 float est = clQest * clDest;
367 cl = subcluster_index;
371 cls[
cl] += subpixel->second.adc;
372 clusterForPixel[subpixel_counter] =
cl;
373 weightOfPixel[subpixel_counter] = maxEst;
375 std::cout <<
"Pixel weight j cl " << weightOfPixel[subpixel_counter] <<
" " << subpixel_counter <<
" " << cl
382 for (
unsigned int subcluster_index = 0; subcluster_index < meanExp; subcluster_index++) {
383 if (
std::abs(clx[subcluster_index] - oldclx[subcluster_index]) > 0.01f)
385 if (
std::abs(cly[subcluster_index] - oldcly[subcluster_index]) > 0.01f)
387 oldclx[subcluster_index] = clx[subcluster_index];
388 oldcly[subcluster_index] = cly[subcluster_index];
389 clx[subcluster_index] = 0;
390 cly[subcluster_index] = 0;
391 cls[subcluster_index] = 1
e-99;
393 for (
unsigned int pixel_index = 0; pixel_index < pixels.size(); pixel_index++) {
394 if (clusterForPixel[pixel_index] < 0)
397 std::cout <<
"j " << pixel_index <<
" x " << pixels[pixel_index].second.x <<
" * y " 398 << pixels[pixel_index].second.y <<
" * ADC " << pixels[pixel_index].second.adc <<
" * W " 399 << weightOfPixel[pixel_index] << std::endl;
400 clx[clusterForPixel[pixel_index]] += pixels[pixel_index].second.x * pixels[pixel_index].second.adc;
401 cly[clusterForPixel[pixel_index]] += pixels[pixel_index].second.y * pixels[pixel_index].second.adc;
402 cls[clusterForPixel[pixel_index]] += pixels[pixel_index].second.adc;
404 for (
unsigned int subcluster_index = 0; subcluster_index < meanExp; subcluster_index++) {
405 if (cls[subcluster_index] != 0) {
406 clx[subcluster_index] /= cls[subcluster_index];
407 cly[subcluster_index] /= cls[subcluster_index];
410 std::cout <<
"Center for cluster " << subcluster_index <<
" x,y " << clx[subcluster_index] <<
" " 411 << cly[subcluster_index] << std::endl;
412 cls[subcluster_index] = 0;
416 std::cout <<
"maxstep " << remainingSteps << std::endl;
418 std::vector<std::vector<SiPixelCluster::Pixel>> pixelsForCl(meanExp);
419 for (
int cl = 0; cl < (
int)meanExp; cl++) {
420 for (
unsigned int j = 0;
j < pixels.size();
j++) {
421 if (clusterForPixel[
j] == cl and pixels[
j].
second.adc != 0) {
425 for (
unsigned int k =
j + 1; k < pixels.size(); k++) {
426 if (pixels[k].
second.adc != 0 and pixels[k].second.x == pixels[
j].second.x and
427 pixels[k].second.y == pixels[
j].second.y and clusterForPixel[k] == cl) {
429 std::cout <<
"Resetting all sub-pixel for location " << pixels[
k].second.x <<
", " << pixels[
k].second.y
430 <<
" at index " << k <<
" associated to cl " << clusterForPixel[
k] << std::endl;
431 pixels[
j].second.adc += pixels[
k].second.adc;
432 pixels[
k].second.adc = 0;
435 for (
unsigned int p = 0;
p < pixels.size(); ++
p)
437 std::cout <<
"index, x, y, ADC: " <<
p <<
", " << pixels[
p].second.x <<
", " << pixels[
p].second.y <<
", " 438 << pixels[
p].second.adc <<
" associated to cl " << clusterForPixel[
p] << std::endl
439 <<
"Adding pixel " << pixels[
j].second.x <<
", " << pixels[
j].second.y <<
" to cluster " << cl
441 pixelsForCl[
cl].push_back(pixels[
j].
second);
450 for (
int cl = 0; cl < (
int)meanExp; cl++) {
452 std::cout <<
"Pixels of cl " << cl <<
" ";
453 for (
unsigned int j = 0;
j < pixelsForCl[
cl].size();
j++) {
456 std::cout << pixelsForCl[
cl][
j].x <<
"," << pixelsForCl[
cl][
j].y <<
"|";
458 output.emplace_back(newpix, pixelsForCl[cl][
j].
adc);
460 output.back().add(newpix, pixelsForCl[cl][
j].
adc);
465 if (!pixelsForCl[cl].
empty()) {
std::multimap< float, int > secondDistScore(const std::vector< std::vector< float >> &distanceMap)
U second(std::pair< T, U > const &p)
Abs< T >::type abs(const T &t)
constexpr int adc(sample_type sample)
get the ADC sample (12 bits)
const std::vector< Pixel > pixels() const