293 int threshold = param.at(0);
294 int max_gap = param.at(1);
295 int range = param.at(2);
296 int min_length = param.at(3);
297 int muon_length = param.at(4);
301 const int h = 3500;
const int w = 2000;
302 constexpr
int accu_h = h + w + 1 ;
const int accu_w = 180;
303 const int x_c = (w/2);
const int y_c = (h/2);
306 int accu[accu_h][accu_w] = {{0}};
307 int (*adata)[accu_w];
308 adata = &accu[(accu_h-1)/2];
311 vector<vector<int>> data = coords;
312 vector<vector<int>> deaccu;
313 vector<vector<int>> outlines;
314 vector<vector<int>> outhit_idx;
317 int count = coords.size();
318 for ( ; count>0; count--){
319 int idx = rndgen.Uniform(count);
320 int max_val = threshold-1;
321 if ((coords.at(idx)).
empty())
323 int x = coords[idx][0],
y = coords[idx][1], rho = 0, theta = 0;
327 for (
int j=0; j<accu_w; j++){
328 int r = int(round((x-x_c)*cos(j*M_PI/accu_w) + (y-y_c)*sin(j*M_PI/accu_w)));
329 int val = ++(adata[
r][j]);
333 theta = j*180/accu_w;
336 if (max_val < threshold){
337 (coords.at(idx)).clear();
341 vector<vector<int>> endpoint(2, vector<int>(4));
342 vector<int> lines_idx;
344 for (
int k=0;
k<2;
k++){
346 while (
gap < max_gap){
348 if ( (idx+i) == int(data.size()) || (idx+i) <0)
350 if ((data.at(idx+i)).
empty())
352 int x1 = data[idx+i][0], y1 = int(data[idx+i][1]), wire_idx = data[idx+i][2];
354 if (endpoint[
k][0]!= 0){
355 last_x = endpoint[
k][0];
356 diffx =
abs(last_x - x1);
361 int y_val = int(round((rho - (x1 - x_c)*cos(theta*M_PI/180.0))/sin(theta*M_PI/180.0) + y_c));
362 if (
abs(y_val-y1) <= range){
364 endpoint[
k] = {x1, y1, wire_idx, idx+i};
365 (coords.at(idx+i)).clear();
366 (data.at(idx+i)).clear();
368 lines_idx.push_back(wire_idx);
377 for (
int n = (deaccu.size()-1);
n>=0;
n--){
378 int x1 = deaccu[
n][0], y1 = int(deaccu[
n][1]);
379 int y_val = int(round((rho - (x1 - x_c)*cos(theta*M_PI/180.0))/sin(theta*M_PI/180.0) + y_c));
380 if (y1 >= (y_val-range) && y1 <= (y_val+range)){
381 for (
int m=0;
m<accu_w;
m++){
382 int r = int(round((x1-x_c)*cos(
m*M_PI/accu_w) + (y1-y_c)*sin(
m*M_PI/accu_w)));
385 deaccu.erase(deaccu.begin() +
n);
389 int x0_end = endpoint[0][0], y0_end = endpoint[0][1], x1_end = endpoint[1][0], y1_end = endpoint[1][1];
390 int wire0_end = endpoint[0][2], wire1_end = endpoint[1][2];
391 int idx0_end = endpoint[0][3], idx1_end = endpoint[1][3];
392 if ((x0_end==0 && y0_end==0) || (x1_end==0 && y1_end==0))
394 vector<int> outline = {x0_end, y0_end, x1_end, y1_end, wire0_end, wire1_end, idx0_end, idx1_end, rho, theta};
396 outlines.push_back(outline);
398 outhit_idx.push_back(lines_idx);
403 for (
size_t i=0; i<outlines.size(); i++){
405 for (
size_t j=i+1; j<outlines.size() && same ==
false; j++){
406 int xi_coords[2] = {outlines[i][0], outlines[i][2]};
int xj_coords[2] = {outlines[j][0], outlines[j][2]};
407 int yi_coords[2] = {outlines[i][1], outlines[i][3]};
int yj_coords[2] = {outlines[j][1], outlines[j][3]};
408 int rhoi = outlines[i][8], rhoj = outlines[j][8];
409 int thetai = outlines[i][9], thetaj = outlines[j][9];
414 for (
int k=0;
k<2 && same ==
false;
k++){
415 for (
int l=0; l<2 && same ==
false; l++){
417 if ((xi_coords[
k] < (xj_coords[l] + var)) && (xi_coords[
k] > (xj_coords[l] - var)))
419 if ((yi_coords[
k] < (yj_coords[l] + var)) && (yi_coords[
k] > (yj_coords[l] - var)))
421 if ((rhoi < (rhoj + rho_var)) && (rhoi > (rhoj - rho_var)))
423 if ((thetai < (thetaj + theta_var)) && (thetai > (thetaj - theta_var)))
428 outlines[j][2] = outlines[i][0];
429 outlines[j][3] = outlines[i][1];
432 outlines[j][0] = outlines[i][0];
433 outlines[j][1] = outlines[i][1];
438 outlines[j][2] = outlines[i][2];
439 outlines[j][3] = outlines[i][3];
442 outlines[j][0] = outlines[i][2];
443 outlines[j][1] = outlines[i][3];
448 (outlines.at(i)).clear();
450 (outhit_idx.at(j)).insert( (outhit_idx.at(j)).
end(), (outhit_idx.at(i)).
begin(), (outhit_idx.at(i)).
end());
451 (outhit_idx.at(i)).clear();
459 for (
size_t i=0; i < outlines.size(); i++){
460 if ((outlines.at(i)).
empty())
462 int x0_end = outlines[i][0], y0_end = outlines[i][1], x1_end = outlines[i][2], y1_end = outlines[i][3];
464 int y_diff =
abs(y0_end-y1_end);
465 if (y_diff > muon_length){
466 lines.push_back(outlines.at(i));
468 hit_idx.push_back(outhit_idx.at(i));
472 float length =
Distance(x0_end,y0_end,x1_end,y1_end);
473 if (length > min_length){
474 lines.push_back(outlines.at(i));
476 hit_idx.push_back(outhit_idx.at(i));
481 data.clear(); deaccu.clear(); outlines.clear(); outhit_idx.clear();
process_name opflash particleana ie x
tuple m
now if test mode generate materials, CRT shell, world, gdml header else just generate CRT shell for u...
auto counter(T begin, T end)
Returns an object to iterate values from begin to end in a range-for loop.
process_name opflash particleana ie ie y
float Distance(int x1, int y1, int x2, int y2)
auto end(FixedBins< T, C > const &) noexcept
auto begin(FixedBins< T, C > const &) noexcept
std::size_t count(Cont const &cont)
bool empty(FixedBins< T, C > const &) noexcept
1.8.5