256 std::vector<SiPixelCluster>
output;
258 unsigned int meanExp = nSplitted;
260 output.push_back(aCluster);
264 std::vector<float> clx(meanExp);
265 std::vector<float> cly(meanExp);
266 std::vector<float> cls(meanExp);
267 std::vector<float> oldclx(meanExp);
268 std::vector<float> oldcly(meanExp);
269 std::vector<SiPixelCluster::Pixel> originalpixels = aCluster.
pixels();
270 std::vector<std::pair<int, SiPixelCluster::Pixel>> pixels;
271 for (
unsigned int j = 0;
j < originalpixels.size();
j++) {
275 int perDiv = originalpixels[
j].adc / sub;
277 std::cout <<
"Splitting " <<
j <<
" in [ " << pixels.size() <<
" , " << pixels.size() + sub
278 <<
" ], expected numb of clusters: " << meanExp <<
" original pixel (x,y) " << originalpixels[
j].x
279 <<
" " << originalpixels[
j].y <<
" sub " << sub << std::endl;
280 for (
int k = 0;
k < sub;
k++) {
282 perDiv = originalpixels[
j].adc - perDiv *
k;
286 std::vector<int> clusterForPixel(pixels.size());
288 for (
unsigned int j = 0;
j < meanExp;
j++) {
291 clx[
j] = originalpixels[0].x +
j;
292 cly[
j] = originalpixels[0].y +
j;
296 int remainingSteps = 100;
297 while (!stop && remainingSteps > 0) {
300 std::vector<std::vector<float>> distanceMapX(originalpixels.size(), std::vector<float>(meanExp));
301 std::vector<std::vector<float>> distanceMapY(originalpixels.size(), std::vector<float>(meanExp));
302 std::vector<std::vector<float>> distanceMap(originalpixels.size(), std::vector<float>(meanExp));
303 for (
unsigned int j = 0;
j < originalpixels.size();
j++) {
305 std::cout <<
"Original Pixel pos " <<
j <<
" " << pixels[
j].second.x <<
" " << pixels[
j].second.y << std::endl;
306 for (
unsigned int i = 0;
i < meanExp;
i++) {
307 distanceMapX[
j][
i] = 1.f * originalpixels[
j].x - clx[
i];
308 distanceMapY[
j][
i] = 1.f * originalpixels[
j].y - cly[
i];
315 dist += (2.f * distanceMapX[
j][
i] /
sizeX) * (2.
f * distanceMapX[
j][
i] /
sizeX);
322 dist += 1.f * (2.f * distanceMapY[
j][
i] /
sizeY) * (2.
f * distanceMapY[
j][
i] /
sizeY);
324 distanceMap[
j][
i] =
sqrt(dist);
326 std::cout <<
"Cluster " <<
i <<
" Original Pixel " <<
j <<
" distances: " << distanceMapX[
j][
i] <<
" " 327 << distanceMapY[
j][
i] <<
" " << distanceMap[
j][
i] << std::endl;
333 std::multimap<float, int> scores;
340 std::vector<float> weightOfPixel(pixels.size());
341 for (std::multimap<float, int>::iterator it = scores.begin(); it != scores.end(); it++) {
342 int pixel_index = it->second;
344 std::cout <<
"Original Pixel " << pixel_index <<
" with score " << it->first << std::endl;
346 int subpixel_counter = 0;
347 for (
auto subpixel = pixels.begin(); subpixel != pixels.end(); ++subpixel, ++subpixel_counter) {
348 if (subpixel->first > pixel_index) {
350 }
else if (subpixel->first != pixel_index) {
355 for (
unsigned int subcluster_index = 0; subcluster_index < meanExp; subcluster_index++) {
356 float nsig = (cls[subcluster_index] - expectedADC) /
358 float clQest = 1.f / (1.f +
std::exp(nsig)) + 1
e-6
f;
359 float clDest = 1.f / (distanceMap[pixel_index][subcluster_index] + 0.05f);
362 std::cout <<
" Q: " << clQest <<
" D: " << clDest <<
" " << distanceMap[pixel_index][subcluster_index]
364 float est = clQest * clDest;
366 cl = subcluster_index;
370 cls[
cl] += subpixel->second.adc;
371 clusterForPixel[subpixel_counter] =
cl;
372 weightOfPixel[subpixel_counter] = maxEst;
374 std::cout <<
"Pixel weight j cl " << weightOfPixel[subpixel_counter] <<
" " << subpixel_counter <<
" " <<
cl 381 for (
unsigned int subcluster_index = 0; subcluster_index < meanExp; subcluster_index++) {
382 if (
std::abs(clx[subcluster_index] - oldclx[subcluster_index]) > 0.01f)
384 if (
std::abs(cly[subcluster_index] - oldcly[subcluster_index]) > 0.01f)
386 oldclx[subcluster_index] = clx[subcluster_index];
387 oldcly[subcluster_index] = cly[subcluster_index];
388 clx[subcluster_index] = 0;
389 cly[subcluster_index] = 0;
390 cls[subcluster_index] = 1
e-99;
392 for (
unsigned int pixel_index = 0; pixel_index < pixels.size(); pixel_index++) {
393 if (clusterForPixel[pixel_index] < 0)
396 std::cout <<
"j " << pixel_index <<
" x " << pixels[pixel_index].second.x <<
" * y " 397 << pixels[pixel_index].second.y <<
" * ADC " << pixels[pixel_index].second.adc <<
" * W " 398 << weightOfPixel[pixel_index] << std::endl;
399 clx[clusterForPixel[pixel_index]] += pixels[pixel_index].second.x * pixels[pixel_index].second.adc;
400 cly[clusterForPixel[pixel_index]] += pixels[pixel_index].second.y * pixels[pixel_index].second.adc;
401 cls[clusterForPixel[pixel_index]] += pixels[pixel_index].second.adc;
403 for (
unsigned int subcluster_index = 0; subcluster_index < meanExp; subcluster_index++) {
404 if (cls[subcluster_index] != 0) {
405 clx[subcluster_index] /= cls[subcluster_index];
406 cly[subcluster_index] /= cls[subcluster_index];
409 std::cout <<
"Center for cluster " << subcluster_index <<
" x,y " << clx[subcluster_index] <<
" " 410 << cly[subcluster_index] << std::endl;
411 cls[subcluster_index] = 0;
415 std::cout <<
"maxstep " << remainingSteps << std::endl;
417 std::vector<std::vector<SiPixelCluster::Pixel>> pixelsForCl(meanExp);
418 for (
int cl = 0;
cl < (
int)meanExp;
cl++) {
419 for (
unsigned int j = 0;
j < pixels.size();
j++) {
420 if (clusterForPixel[
j] ==
cl and pixels[
j].
second.adc != 0) {
424 for (
unsigned int k =
j + 1;
k < pixels.size();
k++) {
425 if (pixels[
k].
second.adc != 0 and pixels[
k].second.x == pixels[
j].second.x and
426 pixels[
k].second.y == pixels[
j].second.y and clusterForPixel[
k] ==
cl) {
428 std::cout <<
"Resetting all sub-pixel for location " << pixels[
k].second.x <<
", " << pixels[
k].second.y
429 <<
" at index " <<
k <<
" associated to cl " << clusterForPixel[
k] << std::endl;
430 pixels[
j].second.adc += pixels[
k].second.adc;
431 pixels[
k].second.adc = 0;
434 for (
unsigned int p = 0;
p < pixels.size(); ++
p)
436 std::cout <<
"index, x, y, ADC: " <<
p <<
", " << pixels[
p].second.x <<
", " << pixels[
p].second.y <<
", " 437 << pixels[
p].second.adc <<
" associated to cl " << clusterForPixel[
p] << std::endl
438 <<
"Adding pixel " << pixels[
j].second.x <<
", " << pixels[
j].second.y <<
" to cluster " <<
cl 440 pixelsForCl[
cl].push_back(pixels[
j].
second);
449 for (
int cl = 0;
cl < (
int)meanExp;
cl++) {
452 for (
unsigned int j = 0;
j < pixelsForCl[
cl].size();
j++) {
455 std::cout << pixelsForCl[
cl][
j].x <<
"," << pixelsForCl[
cl][
j].y <<
"|";
464 if (!pixelsForCl[
cl].
empty()) {
std::multimap< float, int > secondDistScore(const std::vector< std::vector< float >> &distanceMap)
U second(std::pair< T, U > const &p)
const std::vector< Pixel > pixels() const
Abs< T >::type abs(const T &t)
uint16_t *__restrict__ uint16_t const *__restrict__ adc