270 std::vector<SiPixelCluster>
output;
272 unsigned int meanExp = nSplitted;
274 output.push_back(aCluster);
278 std::vector<float> clx(meanExp);
279 std::vector<float> cly(meanExp);
280 std::vector<float> cls(meanExp);
281 std::vector<float> oldclx(meanExp);
282 std::vector<float> oldcly(meanExp);
283 std::vector<SiPixelCluster::Pixel> originalpixels = aCluster.
pixels();
284 std::vector<std::pair<int, SiPixelCluster::Pixel> > pixels;
285 for (
unsigned int j = 0; j < originalpixels.size(); j++) {
288 if (sub < 1) sub = 1;
289 int perDiv = originalpixels[j].adc / sub;
291 std::cout <<
"Splitting " << j <<
" in [ " << pixels.size() <<
" , " 292 << pixels.size() + sub <<
" ], expected numb of clusters: " 293 << meanExp <<
" original pixel (x,y) " 294 << originalpixels[j].x <<
" " << originalpixels[j].y
295 <<
" sub " << sub << std::endl;
296 for (
int k = 0;
k < sub;
k++) {
297 if (
k == sub - 1) perDiv = originalpixels[j].adc - perDiv *
k;
299 originalpixels[j].
y, perDiv)));
302 std::vector<int> clusterForPixel(pixels.size());
304 for (
unsigned int j = 0; j < meanExp; j++) {
307 clx[j] = originalpixels[0].x + j;
308 cly[j] = originalpixels[0].y + j;
312 int remainingSteps = 100;
313 while (!stop && remainingSteps > 0) {
316 std::vector<std::vector<float> > distanceMapX(originalpixels.size(),
317 std::vector<float>(meanExp));
318 std::vector<std::vector<float> > distanceMapY(originalpixels.size(),
319 std::vector<float>(meanExp));
320 std::vector<std::vector<float> > distanceMap(originalpixels.size(),
321 std::vector<float>(meanExp));
322 for (
unsigned int j = 0; j < originalpixels.size(); j++) {
324 std::cout <<
"Original Pixel pos " << j <<
" " << pixels[j].second.x <<
" " 325 << pixels[j].second.y << std::endl;
326 for (
unsigned int i = 0;
i < meanExp;
i++) {
327 distanceMapX[j][
i] = 1.f * originalpixels[j].x - clx[
i];
328 distanceMapY[j][
i] = 1.f * originalpixels[j].y - cly[
i];
336 (2.f*distanceMapX[j][
i]/
sizeX)*(2.
f*distanceMapX[j][i]/
sizeX);
340 dist += 1.f * (
std::abs(distanceMapY[j][i]) -
sizeY/2.f + 1.f) *
343 dist += 1.f * (2.f*distanceMapY[j][
i]/
sizeY) *
344 (2.
f*distanceMapY[j][i]/
sizeY);
346 distanceMap[j][
i] =
sqrt(dist);
348 std::cout <<
"Cluster " << i <<
" Original Pixel " << j
349 <<
" distances: " << distanceMapX[j][
i] <<
" " 350 << distanceMapY[j][
i] <<
" " << distanceMap[j][
i]
357 std::multimap<float, int> scores;
364 std::vector<float> weightOfPixel(pixels.size());
365 for (std::multimap<float, int>::iterator it = scores.begin();
366 it != scores.end(); it++) {
367 int pixel_index = it->second;
369 std::cout <<
"Original Pixel " << pixel_index <<
" with score " << it->first << std::endl;
371 int subpixel_counter = 0;
372 for (
auto subpixel = pixels.begin(); subpixel != pixels.end();
373 ++subpixel, ++subpixel_counter) {
374 if (subpixel->first > pixel_index) {
376 }
else if (subpixel->first != pixel_index) {
381 for (
unsigned int subcluster_index = 0;
382 subcluster_index < meanExp; subcluster_index++) {
384 (cls[subcluster_index] - expectedADC) /
387 float clQest = 1.f / (1.f +
std::exp(nsig)) + 1
e-6
f;
388 float clDest = 1.f / (distanceMap[pixel_index][subcluster_index] + 0.05f);
391 std::cout <<
" Q: " << clQest <<
" D: " << clDest <<
" " 392 << distanceMap[pixel_index][subcluster_index] << std::endl;
393 float est = clQest * clDest;
395 cl = subcluster_index;
399 cls[
cl] += subpixel->second.adc;
400 clusterForPixel[subpixel_counter] =
cl;
401 weightOfPixel[subpixel_counter] = maxEst;
403 std::cout <<
"Pixel weight j cl " << weightOfPixel[subpixel_counter]
404 <<
" " << subpixel_counter
405 <<
" " << cl << std::endl;
411 for (
unsigned int subcluster_index = 0;
412 subcluster_index < meanExp; subcluster_index++) {
413 if (
std::abs(clx[subcluster_index] - oldclx[subcluster_index]) > 0.01f)
415 if (
std::abs(cly[subcluster_index] - oldcly[subcluster_index]) > 0.01f)
417 oldclx[subcluster_index] = clx[subcluster_index];
418 oldcly[subcluster_index] = cly[subcluster_index];
419 clx[subcluster_index] = 0;
420 cly[subcluster_index] = 0;
421 cls[subcluster_index] = 1
e-99;
423 for (
unsigned int pixel_index = 0;
424 pixel_index < pixels.size(); pixel_index++) {
425 if (clusterForPixel[pixel_index] < 0)
continue;
427 std::cout <<
"j " << pixel_index <<
" x " << pixels[pixel_index].second.x <<
" * y " 428 << pixels[pixel_index].second.y <<
" * ADC " 429 << pixels[pixel_index].second.adc <<
" * W " 430 << weightOfPixel[pixel_index] << std::endl;
431 clx[clusterForPixel[pixel_index]] += pixels[pixel_index].second.x * pixels[pixel_index].second.adc;
432 cly[clusterForPixel[pixel_index]] += pixels[pixel_index].second.y * pixels[pixel_index].second.adc;
433 cls[clusterForPixel[pixel_index]] += pixels[pixel_index].second.adc;
435 for (
unsigned int subcluster_index = 0;
436 subcluster_index < meanExp; subcluster_index++) {
437 if (cls[subcluster_index] != 0) {
438 clx[subcluster_index] /= cls[subcluster_index];
439 cly[subcluster_index] /= cls[subcluster_index];
442 std::cout <<
"Center for cluster " << subcluster_index <<
" x,y " 443 << clx[subcluster_index] <<
" " 444 << cly[subcluster_index] << std::endl;
445 cls[subcluster_index] = 0;
450 std::vector<std::vector<SiPixelCluster::Pixel> > pixelsForCl(meanExp);
451 for (
int cl = 0; cl < (
int)meanExp; cl++) {
452 for (
unsigned int j = 0; j < pixels.size(); j++) {
453 if (clusterForPixel[j] == cl and
454 pixels[j].
second.adc != 0) {
458 for (
unsigned int k = j + 1; k < pixels.size(); k++) {
459 if (pixels[k].
second.adc != 0
460 and pixels[k].second.x == pixels[j].second.x
461 and pixels[k].second.y == pixels[j].second.y
462 and clusterForPixel[k] == cl) {
464 std::cout <<
"Resetting all sub-pixel for location " 465 << pixels[
k].second.x <<
", " << pixels[
k].second.y
466 <<
" at index " << k <<
" associated to cl " 467 << clusterForPixel[
k] << std::endl;
468 pixels[j].second.adc += pixels[
k].second.adc;
469 pixels[
k].second.adc = 0;
472 for (
unsigned int p = 0;
p < pixels.size(); ++
p)
474 std::cout <<
"index, x, y, ADC: " <<
p <<
", " 475 << pixels[
p].second.x <<
", " << pixels[
p].second.y
476 <<
", " << pixels[
p].second.adc
477 <<
" associated to cl " << clusterForPixel[
p] << std::endl
478 <<
"Adding pixel " << pixels[j].second.x <<
", " << pixels[j].second.y
479 <<
" to cluster " << cl << std::endl;
480 pixelsForCl[
cl].push_back(pixels[j].
second);
489 for (
int cl = 0; cl < (
int)meanExp; cl++) {
491 for (
unsigned int j = 0; j < pixelsForCl[
cl].size(); j++) {
493 pixelsForCl[cl][j].
y);
495 std::cout << pixelsForCl[
cl][j].x <<
"," << pixelsForCl[
cl][j].y <<
"|";
497 output.emplace_back(newpix, pixelsForCl[cl][j].
adc);
499 output.back().add(newpix, pixelsForCl[cl][j].
adc);
503 if (!pixelsForCl[cl].
empty()) {
int adc(sample_type sample)
get the ADC sample (12 bits)
U second(std::pair< T, U > const &p)
Abs< T >::type abs(const T &t)
std::multimap< float, int > secondDistScore(const std::vector< std::vector< float > > &distanceMap)
const std::vector< Pixel > pixels() const