27 #include <TMatrixFfwd.h>
29 #include <TMatrixTUtils.h>
31 #include <boost/graph/adjacency_list.hpp>
32 #include <boost/graph/connected_components.hpp>
40 using namespace boost;
48 inline constexpr
T square(
const T&
x ) {
return x*
x; }
53 if (get_z_clustering(layer))
80 , m_clusterlist(nullptr)
81 , m_clusterhitassoc(nullptr)
82 , _fraction_of_mip(0.5)
83 , _thresholds_by_layer()
84 , _make_z_clustering()
111 cout <<
PHWHERE <<
"DST Node missing, doing nothing." << endl;
117 auto trkrclusters = findNode::getClass<TrkrClusterContainer>(dstNode,
"TRKR_CLUSTER");
126 dstNode->addNode(DetNode);
132 DetNode->
addNode(TrkrClusterContainerNode);
135 auto clusterhitassoc = findNode::getClass<TrkrClusterHitAssoc>(topNode,
"TRKR_CLUSTERHITASSOC");
144 dstNode->addNode(DetNode);
165 cout <<
"====================== InttClusterizer::InitRun() =====================" << endl;
166 cout <<
" Fraction of expected thickness MIP energy = " <<
_fraction_of_mip << endl;
171 cout <<
" Cluster Threshold in Layer #" << iter->first <<
" = " << 1.0e6 * iter->second <<
" keV" << endl;
177 cout <<
" Z-dimension Clustering in Layer #" << iter->first <<
" = " << boolalpha << iter->second << noboolalpha << endl;
183 cout <<
" Energy weighting clusters in Layer #" << iter->first <<
" = " << boolalpha << iter->second << noboolalpha << endl;
185 cout <<
"===========================================================================" << endl;
195 m_hits = findNode::getClass<TrkrHitSetContainer>(topNode,
"TRKR_HITSET");
198 cout <<
PHWHERE <<
"ERROR: Can't find node TRKR_HITSET" << endl;
203 m_clusterlist = findNode::getClass<TrkrClusterContainer>(topNode,
"TRKR_CLUSTER");
206 cout <<
PHWHERE <<
" ERROR: Can't find TRKR_CLUSTER." << endl;
211 m_clusterhitassoc = findNode::getClass<TrkrClusterHitAssoc>(topNode,
"TRKR_CLUSTERHITASSOC");
214 cout <<
PHWHERE <<
" ERROR: Can't find TRKR_CLUSTERHITASSOC" << endl;
234 if (!geom_container)
return;
238 layeriter != layerrange.second;
242 int layer = layeriter->second->get_layer();
267 cout <<
"Entering InttClusterizer::ClusterLadderCells " << endl;
275 if (!geom_container)
return;
285 hitsetitr != hitsetrange.second;
291 if(
Verbosity() > 1) cout <<
"InttClusterizer found hitsetkey " << hitsetitr->first << endl;
305 std::vector <std::pair< TrkrDefs::hitkey, TrkrHit*> > hitvec;
308 hitr != hitrangei.second;
311 hitvec.push_back(make_pair(hitr->first, hitr->second));
314 cout <<
"hitvec.size(): " << hitvec.size() << endl;
316 typedef adjacency_list<vecS, vecS, undirectedS> Graph;
320 for (
unsigned int i = 0; i < hitvec.size(); i++)
322 for (
unsigned int j = i + 1; j < hitvec.size(); j++)
335 vector<int> component(num_vertices(G));
338 connected_components(G, &component[0]);
342 set<int> cluster_ids;
344 multimap<int, std::pair<TrkrDefs::hitkey, TrkrHit*> > clusters;
345 for (
unsigned int i = 0; i < component.size(); i++)
347 cluster_ids.insert(component[i]);
348 clusters.insert(make_pair(component[i], hitvec[i]));
352 for (set<int>::iterator clusiter = cluster_ids.begin(); clusiter != cluster_ids.end(); ++clusiter)
354 int clusid = *clusiter;
357 pair<multimap<int, std::pair<TrkrDefs::hitkey, TrkrHit*>>::iterator,
358 multimap<int, std::pair<TrkrDefs::hitkey, TrkrHit*>>::iterator> clusrange = clusters.equal_range(clusid);
359 multimap<int, std::pair<TrkrDefs::hitkey, TrkrHit*>>::iterator mapiter = clusrange.first;
363 auto clus = std::make_unique<TrkrClusterv3>();
364 clus->setClusKey(ckey);
367 cout <<
"Filling cluster with key " << ckey << endl;
374 double xlocalsum = 0.0;
375 double ylocalsum = 0.0;
376 double zlocalsum = 0.0;
377 unsigned int clus_adc = 0.0;
380 for (mapiter = clusrange.first; mapiter != clusrange.second; ++mapiter)
390 unsigned int hit_adc = (mapiter->second).
second->getAdc();
393 double local_hit_location[3] = {0., 0., 0.};
400 xlocalsum += local_hit_location[0] * (double) hit_adc;
401 ylocalsum += local_hit_location[1] * (double) hit_adc;
402 zlocalsum += local_hit_location[2] * (double) hit_adc;
406 xlocalsum += local_hit_location[0];
407 ylocalsum += local_hit_location[1];
408 zlocalsum += local_hit_location[2];
417 if (
Verbosity() > 2) cout <<
" nhits = " << nhits << endl;
420 cout <<
" From geometry object: hit x " << local_hit_location[0] <<
" hit y " << local_hit_location[1] <<
" hit z " << local_hit_location[2] << endl;
421 cout <<
" nhits " << nhits <<
" clusx = " << xlocalsum / nhits <<
" clusy " << ylocalsum / nhits <<
" clusz " << zlocalsum / nhits <<
" hit_adc " << hit_adc << endl;
426 static const float invsqrt12 = 1./sqrt(12);
435 float phierror = pitch * invsqrt12;
437 static constexpr std::array<double, 3> scalefactors_phi = {{ 0.85, 0.4, 0.33 }};
438 if( phibins.size() == 1 && layer < 5) phierror*=scalefactors_phi[0];
439 else if( phibins.size() == 2 && layer < 5) phierror*=scalefactors_phi[1];
440 else if( phibins.size() == 2 && layer > 4) phierror*=scalefactors_phi[2];
442 const float zerror = length * invsqrt12;
444 double cluslocaly = NAN;
445 double cluslocalz = NAN;
449 cluslocaly = ylocalsum / (double) clus_adc;
450 cluslocalz = zlocalsum / (double) clus_adc;
454 cluslocaly = ylocalsum / nhits;
455 cluslocalz = zlocalsum / nhits;
459 clus->setAdc(clus_adc);
463 clus->setLocalX(cluslocaly);
464 clus->setLocalY(cluslocalz);
467 clus->setSubSurfKey(0);
468 clus->setActsLocalError(0,0,
square(phierror));
469 clus->setActsLocalError(0,1, 0.);
470 clus->setActsLocalError(1,0, 0.);
471 clus->setActsLocalError(1,1,
square(zerror));
481 cout <<
"After InttClusterizer, cluster-hit associations are:" << endl;
492 TrkrClusterContainer *clusterlist = findNode::getClass<TrkrClusterContainer>(topNode,
"TRKR_CLUSTER");
493 if (!clusterlist)
return;
495 cout <<
"================= After InttClusterizer::process_event() ====================" << endl;
497 cout <<
" There are " << clusterlist->
size() <<
" clusters recorded: " << endl;
501 cout <<
"===========================================================================" << endl;