33 #include <TMatrixFfwd.h>
35 #include <TMatrixTUtils.h>
38 #include <boost/graph/adjacency_list.hpp>
39 #include <boost/graph/connected_components.hpp>
50 using namespace boost;
58 inline constexpr
T square(
const T&
x ) {
return x*
x; }
91 , m_clusterlist(nullptr)
92 , m_clusterhitassoc(nullptr)
93 , m_makeZClustering(
true)
109 cout <<
PHWHERE <<
"DST Node missing, doing nothing." << endl;
123 auto trkrclusters = findNode::getClass<TrkrClusterContainer>(dstNode,
"TRKR_CLUSTER");
132 dstNode->addNode(DetNode);
138 DetNode->
addNode(TrkrClusterContainerNode);
141 auto clusterhitassoc = findNode::getClass<TrkrClusterHitAssoc>(topNode,
"TRKR_CLUSTERHITASSOC");
150 dstNode->addNode(DetNode);
165 cout <<
"====================== MvtxClusterizer::InitRun() =====================" << endl;
166 cout <<
" Z-dimension Clustering = " << boolalpha <<
m_makeZClustering << noboolalpha << endl;
167 cout <<
"===========================================================================" << endl;
176 m_hits = findNode::getClass<TrkrHitSetContainer>(topNode,
"TRKR_HITSET");
179 cout <<
PHWHERE <<
"ERROR: Can't find node TRKR_HITSET" << endl;
184 m_clusterlist = findNode::getClass<TrkrClusterContainer>(topNode,
"TRKR_CLUSTER");
187 cout <<
PHWHERE <<
" ERROR: Can't find TRKR_CLUSTER." << endl;
192 m_clusterhitassoc = findNode::getClass<TrkrClusterHitAssoc>(topNode,
"TRKR_CLUSTERHITASSOC");
195 cout <<
PHWHERE <<
" ERROR: Can't find TRKR_CLUSTERHITASSOC" << endl;
210 cout <<
"Entering MvtxClusterizer::ClusterMvtx " << endl;
213 if (!geom_container)
return;
223 hitsetitr != hitsetrange.second;
228 if(
Verbosity() > 1) cout <<
"MvtxClusterizer found hitsetkey " << hitsetitr->first << endl;
234 std::vector <std::pair< TrkrDefs::hitkey, TrkrHit*> > hitvec;
238 hitr != hitrangei.second;
241 hitvec.push_back(make_pair(hitr->first, hitr->second));
243 if (
Verbosity() > 2) cout <<
"hitvec.size(): " << hitvec.size() << endl;
246 typedef adjacency_list<vecS, vecS, undirectedS> Graph;
250 for (
unsigned int i = 0; i < hitvec.size(); i++)
252 for (
unsigned int j = 0; j < hitvec.size(); j++)
261 vector<int> component(num_vertices(G));
264 connected_components(G, &component[0]);
268 set<int> cluster_ids;
270 multimap<int, std::pair<TrkrDefs::hitkey, TrkrHit*> > clusters;
271 for (
unsigned int i = 0; i < component.size(); i++)
273 cluster_ids.insert(component[i]);
274 clusters.insert(make_pair(component[i], hitvec[i]));
278 for (set<int>::iterator clusiter = cluster_ids.begin(); clusiter != cluster_ids.end(); ++clusiter)
280 int clusid = *clusiter;
281 auto clusrange = clusters.equal_range(clusid);
283 if (
Verbosity() > 2) cout <<
"Filling cluster id " << clusid <<
" of " << std::distance(cluster_ids.begin(),clusiter )<< endl;
288 auto clus = std::make_unique<TrkrClusterv3>();
289 clus->setClusKey(ckey);
298 const unsigned int nhits = std::distance( clusrange.first, clusrange.second );
300 double locclusx = NAN;
301 double locclusz = NAN;
309 for (
auto mapiter = clusrange.first; mapiter != clusrange.second; ++mapiter)
318 auto local_coords = layergeom->get_local_coords_from_pixel(row,col);
325 local_coords.SetY( 1
e-4 );
328 locxsum += local_coords.X();
329 loczsum += local_coords.Z();
336 locclusx = locxsum / nhits;
337 locclusz = loczsum / nhits;
341 const double pitch = layergeom->get_pixel_x();
342 const double length = layergeom->get_pixel_z();
343 const double phisize = phibins.size() * pitch;
344 const double zsize = zbins.size() *
length;
346 static const double invsqrt12 = 1./std::sqrt(12);
355 double phierror = pitch * invsqrt12;
357 static constexpr std::array<double, 7> scalefactors_phi = {{ 0.36, 0.6,0.37,0.49,0.4,0.37,0.33 }};
358 if(phibins.size() == 1 && zbins.size() == 1) phierror*=scalefactors_phi[0];
359 else if(phibins.size() == 2 && zbins.size() == 1) phierror*=scalefactors_phi[1];
360 else if(phibins.size() == 1 && zbins.size() == 2) phierror*=scalefactors_phi[2];
361 else if( phibins.size() == 2 && zbins.size() == 2 ) phierror*=scalefactors_phi[0];
362 else if( phibins.size() == 2 && zbins.size() == 3 ) phierror*=scalefactors_phi[1];
363 else if( phibins.size() == 3 && zbins.size() == 2 ) phierror*=scalefactors_phi[2];
364 else if( phibins.size() == 3 && zbins.size() == 3 ) phierror*=scalefactors_phi[3];
372 static constexpr std::array<double, 4> scalefactors_z = {{ 0.47, 0.48, 0.71, 0.55 }};
373 double zerror = length*invsqrt12;
374 if( zbins.size() == 2 && phibins.size() == 2 ) zerror*=scalefactors_z[0];
375 else if( zbins.size() == 2 && phibins.size() == 3 ) zerror*=scalefactors_z[1];
376 else if( zbins.size() == 3 && phibins.size() == 2 ) zerror*=scalefactors_z[2];
377 else if( zbins.size() == 3 && phibins.size() == 3 ) zerror*=scalefactors_z[3];
380 cout <<
" MvtxClusterizer: layer " << layer <<
" rad " << layergeom->get_radius() <<
" phibins " << phibins.size() <<
" pitch " << pitch <<
" phisize " << phisize
381 <<
" zbins " << zbins.size() <<
" length " << length <<
" zsize " << zsize << endl;
383 clus->setLocalX(locclusx);
384 clus->setLocalY(locclusz);
386 clus->setActsLocalError(0,0,
square(phierror));
387 clus->setActsLocalError(0,1,0.);
388 clus->setActsLocalError(1,0,0.);
389 clus->setActsLocalError(1,1,
square(zerror));
393 clus->setSubSurfKey(0);
416 TrkrClusterContainer *clusterlist = findNode::getClass<TrkrClusterContainer>(topNode,
"TRKR_CLUSTER");
417 if (!clusterlist)
return;
419 cout <<
"================= Aftyer MvtxClusterizer::process_event() ====================" << endl;
421 cout <<
" There are " << clusterlist->
size() <<
" clusters recorded: " << endl;
425 cout <<
"===========================================================================" << endl;