30 #include <TMatrixFfwd.h>
32 #include <TMatrixTUtils.h>
48 template<
class T>
inline constexpr
T square(
const T&
x ) {
return x*
x; }
50 typedef std::pair<unsigned short, unsigned short>
iphiz;
51 typedef std::pair<unsigned short, iphiz>
ihit;
59 unsigned int layer = 0;
61 unsigned int sector = 0;
65 unsigned short phioffset = 0;
66 unsigned short zbins = 0;
67 unsigned short zoffset = 0;
68 unsigned short maxHalfSizeZ = 0;
69 unsigned short maxHalfSizePhi = 0;
70 double m_drift_velocity_scale = 1.0;
73 std::map<TrkrDefs::cluskey, TrkrCluster *> *clusterlist =
nullptr;
74 std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey> *clusterhitassoc =
nullptr;
77 pthread_mutex_t mythreadlock;
79 void remove_hit(
double adc,
int phibin,
int zbin, std::multimap<unsigned short, ihit> &all_hit_map, std::vector<std::vector<unsigned short>> &adcval)
81 typedef std::multimap<unsigned short, ihit>::iterator hit_iterator;
82 std::pair<hit_iterator, hit_iterator> iterpair = all_hit_map.equal_range(adc);
83 hit_iterator
it = iterpair.first;
84 for (; it != iterpair.second; ++
it) {
85 if (it->second.second.first == phibin && it->second.second.second == zbin) {
86 all_hit_map.erase(it);
90 adcval[phibin][zbin] = 0;
93 void remove_hits(std::vector<ihit> &ihit_list, std::multimap<unsigned short, ihit> &all_hit_map,std::vector<std::vector<unsigned short>> &adcval)
95 for(
auto iter = ihit_list.begin(); iter != ihit_list.end();++iter){
96 unsigned short adc = iter->first;
97 unsigned short phibin = iter->second.first;
98 unsigned short zbin = iter->second.second;
99 remove_hit(adc,phibin,zbin,all_hit_map,adcval);
103 void find_z_range(
int phibin,
int zbin,
const thread_data& my_data,
const std::vector<std::vector<unsigned short>> &adcval,
int& zdown,
int& zup){
105 const int FitRangeZ = (
int) my_data.maxHalfSizeZ;
106 const int NZBinsMax = (
int) my_data.zbins;
109 for(
int iz=0; iz< FitRangeZ; iz++){
112 if(cz <= 0 || cz >= NZBinsMax){
117 if(adcval[phibin][cz] <= 0) {
122 if(adcval[phibin][cz]+adcval[phibin][cz+1] <
123 adcval[phibin][cz+2]+adcval[phibin][cz+3]){
130 for(
int iz=0; iz< FitRangeZ; iz++){
132 if(cz <= 0 || cz >= NZBinsMax){
136 if(adcval[phibin][cz] <= 0) {
141 if(adcval[phibin][cz]+adcval[phibin][cz-1] <
142 adcval[phibin][cz-2]+adcval[phibin][cz-3]){
151 void find_phi_range(
int phibin,
int zbin,
const thread_data& my_data,
const std::vector<std::vector<unsigned short>> &adcval,
int& phidown,
int& phiup)
154 int FitRangePHI = (
int) my_data.maxHalfSizePhi;
155 int NPhiBinsMax = (
int) my_data.phibins;
158 for(
int iphi=0; iphi< FitRangePHI; iphi++){
159 int cphi = phibin + iphi;
160 if(cphi < 0 || cphi >= NPhiBinsMax){
166 if(adcval[cphi][zbin] <= 0) {
171 if(cphi<NPhiBinsMax-4){
172 if(adcval[cphi][zbin]+adcval[cphi+1][zbin] <
173 adcval[cphi+2][zbin]+adcval[cphi+3][zbin]){
181 for(
int iphi=0; iphi< FitRangePHI; iphi++){
182 int cphi = phibin - iphi;
183 if(cphi < 0 || cphi >= NPhiBinsMax){
188 if(adcval[cphi][zbin] <= 0) {
194 if(adcval[cphi][zbin]+adcval[cphi-1][zbin] <
195 adcval[cphi-2][zbin]+adcval[cphi-3][zbin]){
204 void get_cluster(
int phibin,
int zbin,
const thread_data& my_data,
const std::vector<std::vector<unsigned short>> &adcval, std::vector<ihit> &ihit_list)
210 find_z_range(phibin, zbin, my_data, adcval, zdown, zup);
213 for(
int iz=zbin - zdown ; iz<= zbin + zup; iz++){
216 find_phi_range(phibin, iz, my_data, adcval, phidown, phiup);
217 for (
int iphi = phibin - phidown; iphi <= (phibin + phiup); iphi++){
218 iphiz iCoord(std::make_pair(iphi,iz));
219 ihit thisHit(adcval[iphi][iz],iCoord);
220 ihit_list.push_back(thisHit);
233 std::map<unsigned int, std::vector<Surface>>::iterator mapIter;
239 <<
"Error: hitsetkey not found in clusterSurfaceMap, hitsetkey = "
240 << hitsetkey << std::endl;
244 double world_phi = atan2(world[1], world[0]);
245 double world_z = world[2];
247 std::vector<Surface> surf_vec = mapIter->second;
248 unsigned int surf_index = 999;
252 double fraction = (world_phi +
M_PI) / (2.0 *
M_PI);
253 double rounded_nsurf =
round( (
double) (surf_vec.size()/2) * fraction - 0.5);
254 unsigned int nsurf = (
unsigned int) rounded_nsurf;
256 nsurf += surf_vec.size()/2;
261 Surface this_surf = surf_vec[nsurf];
262 auto vec3d = this_surf->center(tGeometry->
geoContext);
263 std::vector<double> surf_center = {vec3d(0) / 10.0, vec3d(1) / 10.0, vec3d(2) / 10.0};
264 double surf_z = surf_center[2];
265 double surf_phi = atan2(surf_center[1], surf_center[0]);
267 if( (world_phi > surf_phi - surfStepPhi / 2.0 && world_phi < surf_phi + surfStepPhi / 2.0 ) &&
268 (world_z > surf_z - surfStepZ / 2.0 && world_z < surf_z + surfStepZ / 2.0) )
276 <<
"Error: TPC surface index not defined, skipping cluster!"
278 std::cout <<
" coordinates: " << world[0] <<
" " << world[1] <<
" " << world[2]
279 <<
" radius " << sqrt(world[0]*world[0]+world[1]*world[1]) << std::endl;
280 std::cout <<
" world_phi " << world_phi <<
" world_z " << world_z << std::endl;
281 std::cout <<
" surf coords: " << surf_center[0] <<
" " << surf_center[1] <<
" " << surf_center[2] << std::endl;
282 std::cout <<
" surf_phi " << surf_phi <<
" surf_z " << surf_z << std::endl;
283 std::cout <<
" number of surfaces " << surf_vec.size() << std::endl;
287 return surf_vec[surf_index];
291 void calc_cluster_parameter(
const std::vector<ihit> &ihit_list,
int iclus,
const thread_data& my_data )
296 const double z_min = -105.5;
297 const double z_max = 105.5;
301 double phi_sum = 0.0;
302 double adc_sum = 0.0;
304 double phi2_sum = 0.0;
306 double radius = my_data.layergeom->get_radius();
309 int phibinlo = 666666;
312 int clus_size = ihit_list.size();
314 if(clus_size == 1)
return;
316 std::vector<TrkrDefs::hitkey> hitkeyvec;
317 for(
auto iter = ihit_list.begin(); iter != ihit_list.end();++iter){
318 double adc = iter->first;
320 if (adc <= 0)
continue;
322 int iphi = iter->second.first + my_data.phioffset;
323 int iz = iter->second.second + my_data.zoffset;
324 if(iphi > phibinhi) phibinhi = iphi;
325 if(iphi < phibinlo) phibinlo = iphi;
326 if(iz > zbinhi) zbinhi = iz;
327 if(iz < zbinlo) zbinlo = iz;
330 double phi_center = my_data.layergeom->get_phicenter(iphi);
331 phi_sum += phi_center * adc;
332 phi2_sum +=
square(phi_center)*adc;
335 double z = my_data.layergeom->get_zcenter(iz);
340 z*my_data.m_drift_velocity_scale +
341 (z>0 ? z_max:
z_min)*(1.-my_data.m_drift_velocity_scale);
351 hitkeyvec.push_back(hitkey);
353 if (adc_sum < 10)
return;
356 double clusphi = phi_sum / adc_sum;
357 double clusz = z_sum / adc_sum;
359 const double phi_cov = phi2_sum/adc_sum -
square(clusphi);
360 const double z_cov = z2_sum/adc_sum -
square(clusz);
366 auto clus = std::make_unique<TrkrClusterv3>();
367 clus->setClusKey(ckey);
370 const double phi_err_square = (phibinhi == phibinlo) ?
371 square(radius*my_data.layergeom->get_phistep())/12:
372 square(radius)*phi_cov/(adc_sum*0.14);
374 const double z_err_square = (zbinhi == zbinlo) ?
375 square(my_data.layergeom->get_zstep())/12:
376 z_cov/(adc_sum*0.14);
385 clusz -= (clusz<0) ? my_data.par0_neg:my_data.par0_pos;
389 clus->setAdc(adc_sum);
391 float clusx = radius * cos(clusphi);
392 float clusy = radius * sin(clusphi);
399 ERR[1][1] = phi_err_square;
403 ERR[2][2] = z_err_square;
424 clus->setSubSurfKey(subsurfkey);
430 double clusRadius = sqrt(clusx * clusx + clusy * clusy);
431 double rClusPhi = clusRadius * clusphi;
432 double surfRadius = sqrt(center(0)*center(0) + center(1)*center(1));
433 double surfPhiCenter = atan2(center[1], center[0]);
434 double surfRphiCenter = surfPhiCenter * surfRadius;
435 double surfZCenter = center[2];
437 auto local = surface->globalToLocal(my_data.tGeometry->geoContext,
450 localPos(0) = rClusPhi - surfRphiCenter;
451 localPos(1) = clusz - surfZCenter;
454 clus->setLocalX(localPos(0));
455 clus->setLocalY(localPos(1));
456 clus->setActsLocalError(0,0, ERR[1][1]);
457 clus->setActsLocalError(1,0, ERR[2][1]);
458 clus->setActsLocalError(0,1, ERR[1][2]);
459 clus->setActsLocalError(1,1, ERR[2][2]);
464 if(my_data.clusterlist)
466 const auto [iter, inserted] = my_data.clusterlist->insert(std::make_pair(ckey, clus.get()));
472 if( inserted ) clus.release();
477 <<
"Error: duplicated cluster key: " << ckey <<
" - new cluster not inserted in map"
483 if(my_data.do_assoc && my_data.clusterhitassoc){
484 for (
unsigned int i = 0; i < hitkeyvec.size(); i++){
485 my_data.clusterhitassoc->insert(std::make_pair(ckey, hitkeyvec[i]));
490 void *ProcessSector(
void *threadarg) {
492 auto my_data =
static_cast<thread_data*
>(threadarg);
493 const auto& pedestal = my_data->pedestal;
494 const auto&
phibins = my_data->phibins;
495 const auto& phioffset = my_data->phioffset;
496 const auto& zbins = my_data->zbins ;
497 const auto& zoffset = my_data->zoffset ;
503 std::vector<std::vector<unsigned short>> adcval(
phibins, std::vector<unsigned short>(zbins, 0));
504 std::multimap<unsigned short, ihit> all_hit_map;
505 std::vector<ihit> hit_vect;
508 hitr != hitrangei.second;
513 float_t fadc = (hitr->second->getAdc()) - pedestal;
516 unsigned short adc = 0;
518 adc = (
unsigned short) fadc;
521 if(phibin >=
phibins)
continue;
523 if(zbin >= zbins)
continue;
526 iphiz iCoord(std::make_pair(phibin,zbin));
527 ihit thisHit(adc,iCoord);
529 all_hit_map.insert(std::make_pair(adc, thisHit));
532 adcval[phibin][zbin] = (
unsigned short) adc;
537 while(all_hit_map.size()>0){
539 auto iter = all_hit_map.rbegin();
540 if(iter == all_hit_map.rend()){
543 ihit hiHit = iter->second;
544 int iphi = hiHit.second.first;
545 int iz = hiHit.second.second;
550 std::vector<ihit> ihit_list;
551 get_cluster(iphi, iz, *my_data, adcval, ihit_list);
558 calc_cluster_parameter(ihit_list,nclus++, *my_data );
559 remove_hits(ihit_list,all_hit_map, adcval);
561 pthread_exit(
nullptr);
571 bool reject_it =
false;
575 int PhiBinsSector = PhiBins/12;
578 double PhiBinSize = 2.0* radius *
M_PI / (double) PhiBins;
581 int sector_lo = sector * PhiBinsSector;
582 int sector_hi = sector_lo + PhiBinsSector - 1;
586 if(phibin < sector_lo + sector_fiducial_bins || phibin > sector_hi - sector_fiducial_bins)
609 std::cout <<
PHWHERE <<
"DST Node missing, doing nothing." << std::endl;
614 auto trkrclusters = findNode::getClass<TrkrClusterContainer>(dstNode,
"TRKR_CLUSTER");
623 dstNode->addNode(DetNode);
629 DetNode->
addNode(TrkrClusterContainerNode);
632 auto clusterhitassoc = findNode::getClass<TrkrClusterHitAssoc>(topNode,
"TRKR_CLUSTERHITASSOC");
633 if (!clusterhitassoc)
641 dstNode->addNode(DetNode);
657 std::cout <<
"TpcClusterizer::Process_Event" << std::endl;
663 std::cout <<
PHWHERE <<
"DST Node missing, doing nothing." << std::endl;
668 m_hits = findNode::getClass<TrkrHitSetContainer>(topNode,
"TRKR_HITSET");
671 std::cout <<
PHWHERE <<
"ERROR: Can't find node TRKR_HITSET" << std::endl;
676 m_clusterlist = findNode::getClass<TrkrClusterContainer>(topNode,
"TRKR_CLUSTER");
679 std::cout <<
PHWHERE <<
" ERROR: Can't find TRKR_CLUSTER." << std::endl;
684 m_clusterhitassoc = findNode::getClass<TrkrClusterHitAssoc>(topNode,
"TRKR_CLUSTERHITASSOC");
687 std::cout <<
PHWHERE <<
" ERROR: Can't find TRKR_CLUSTERHITASSOC" << std::endl;
692 findNode::getClass<PHG4CylinderCellGeomContainer>(topNode,
"CYLINDERCELLGEOM_SVTX");
695 std::cout <<
PHWHERE <<
"ERROR: Can't find node CYLINDERCELLGEOM_SVTX" << std::endl;
699 m_tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode,
700 "ActsTrackingGeometry");
704 <<
"ActsTrackingGeometry not found on node tree. Exiting"
709 m_surfMaps = findNode::getClass<ActsSurfaceMaps>(topNode,
714 <<
"ActsSurfaceMaps not found on node tree. Exiting"
724 const int num_hitsets = std::distance(hitsetrange.first,hitsetrange.second);
734 std::vector<thread_pair_t> threads;
735 threads.reserve( num_hitsets );
738 pthread_attr_init(&attr);
739 pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
741 if (pthread_mutex_init(&mythreadlock,
nullptr) != 0)
743 printf(
"\n mutex init failed\n");
748 hitsetitr != hitsetrange.second;
751 const auto hitsetid = hitsetitr->first;
759 thread_pair_t& thread_pair = threads.emplace_back();
761 thread_pair.data.layergeom = layergeom;
762 thread_pair.data.hitset = hitset;
763 thread_pair.data.layer =
layer;
764 thread_pair.data.pedestal =
pedestal;
765 thread_pair.data.sector = sector;
766 thread_pair.data.side = side;
775 thread_pair.data.par0_neg =
par0_neg;
776 thread_pair.data.par0_pos =
par0_pos;
778 unsigned short NPhiBins = (
unsigned short) layergeom->
get_phibins();
779 unsigned short NPhiBinsSector = NPhiBins/12;
780 unsigned short NZBins = (
unsigned short)layergeom->
get_zbins();
781 unsigned short NZBinsSide = NZBins/2;
782 unsigned short NZBinsMin = 0;
783 unsigned short PhiOffset = NPhiBinsSector * sector;
789 NZBinsMin = NZBins / 2;
792 unsigned short ZOffset = NZBinsMin;
794 thread_pair.data.phibins = NPhiBinsSector;
795 thread_pair.data.phioffset = PhiOffset;
796 thread_pair.data.zbins = NZBinsSide;
797 thread_pair.data.zoffset = ZOffset ;
799 int rc = pthread_create(&thread_pair.thread, &attr, ProcessSector, (
void *)&thread_pair.data);
801 std::cout <<
"Error:unable to create thread," << rc << std::endl;
805 pthread_attr_destroy(&attr);
808 for(
const auto& thread_pair:threads )
810 int rc2 = pthread_join(thread_pair.thread,
nullptr);
812 { std::cout <<
"Error:unable to join," << rc2 << std::endl; }
816 std::cout <<
"TPC Clusterizer found " <<
m_clusterlist->
size() <<
" Clusters " << std::endl;