15 #if __cplusplus < 201402L
16 #include <boost/make_unique.hpp>
59 m_seedFinderCfg.seedFilter = std::make_unique<Acts::SeedFilter<SpacePoint>>(
64 m_file =
new TFile(
"seedingOutfile.root",
"recreate");
84 _iteration_map = findNode::getClass<TrkrClusterIterationMapv1>(topNode,
"CLUSTER_ITERATION_MAP");
86 std::cerr <<
PHWHERE <<
"Cluster Iteration Map missing, aborting." << std::endl;
91 auto eventTimer = std::make_unique<PHTimer>(
"eventTimer");
93 eventTimer->restart();
96 std::cout <<
"Processing PHActsSiliconSeeding event "
99 std::vector<const SpacePoint*> spVec;
103 auto seederTime = eventTimer->get_accumulated_time();
104 eventTimer->restart();
109 auto circleFitTime = eventTimer->get_accumulated_time();
116 std::cout <<
"Finished PHActsSiliconSeeding process_event"
121 std::cout <<
"PHActsSiliconSeeding Acts seed time "
122 << seederTime << std::endl;
123 std::cout <<
"PHActsSiliconSeeding circle fit time "
124 << circleFitTime << std::endl;
125 std::cout <<
"PHActsSiliconSeeding total event time "
126 << circleFitTime + seederTime << std::endl;
144 <<
" bad initial circle fits" << std::endl;
146 <<
" bad second circle fits" << std::endl;
158 auto covConverter = [=](
const SpacePoint&
sp, float, float, float)
167 std::unique_ptr<Acts::SpacePointGrid<SpacePoint>> grid =
168 Acts::SpacePointGridCreator::createGrid<SpacePoint>(
m_gridCfg);
183 auto groupIt = spGroup.
begin();
184 auto endGroup = spGroup.end();
186 for(; !(groupIt == endGroup); ++groupIt)
199 int numGoodSeeds = 0;
202 for(
auto& seeds : seedVector)
205 for(
auto&
seed : seeds)
208 std::cout <<
"Seed " << numSeeds <<
" has "
209 <<
seed.sp().size() <<
" measurements "
214 std::vector<TrkrCluster*> clusters;
215 std::vector<Acts::Vector3D> globalPositions;
217 for(
auto& spacePoint :
seed.sp())
219 auto cluskey = spacePoint->m_clusKey;
227 std::cout <<
"Adding cluster with x,y "
228 << spacePoint->x() <<
", " << spacePoint->y()
229 <<
" mm in detector "
238 auto fitTimer = std::make_unique<PHTimer>(
"trackfitTimer");
248 auto circlefittime = fitTimer->get_accumulated_time();
261 clusters, globalPositions);
263 auto svtxtracktime = fitTimer->get_accumulated_time();
266 std::cout <<
"Circle fit time " << circlefittime <<
" and svtx time "
267 << svtxtracktime << std::endl;
279 std::cout <<
"Total number of seeds found in "
280 << seedVector.size() <<
" volume regions gives "
281 << numSeeds <<
" Acts seeds " << std::endl;
295 std::vector<TrkrCluster*>& clusters,
296 std::vector<Acts::Vector3D>& clusGlobPos)
298 auto fitTimer = std::make_unique<PHTimer>(
"trackfitTimer");
305 auto possibleStubs = fitTimer->get_accumulated_time();
308 int numSeedsPerActsSeed = 0;
316 double trackCharge =
charge;
317 double trackPhi = atan2(py,px);
318 double trackEta = atanh(pz / sqrt(px * px + py * py + pz * pz));
329 for(
const auto& [stub, stubClusterPairs] : stubs)
331 std::vector<TrkrCluster*> stubClusters = stubClusterPairs.first;
332 std::vector<Acts::Vector3D> stubClusterPositions = stubClusterPairs.second;
336 numSeedsPerActsSeed++;
338 auto svtxTrack = std::make_unique<SvtxTrack_v2>();
342 for(
const auto clus : stubClusters)
344 const auto cluskey = clus->getClusKey();
345 svtxTrack->insert_cluster_key(
cluskey);
352 { std::cout <<
"Cluskey adding : " <<
cluskey << std::endl; }
361 float pt = 0.3 * 1.4 * R / 100.;
363 trackPx = pt * cos(trackPhi);
364 trackPy = pt * sin(trackPhi);
365 trackPz = pt * sinh(trackEta);
376 std::cout <<
"Setting silicon seed with track id "
379 << trackX <<
", " << trackY <<
", " << trackZ
380 << std::endl <<
" and (px,py,pz) " << trackPx
381 <<
", " << trackPy <<
", " << trackPz << std::endl
382 <<
" with charge " << trackCharge << std::endl;
385 svtxTrack->set_x(trackX);
386 svtxTrack->set_y(trackY);
387 svtxTrack->set_z(trackZ);
388 svtxTrack->set_px(trackPx);
389 svtxTrack->set_py(trackPy);
390 svtxTrack->set_pz(trackPz);
391 svtxTrack->set_charge(trackCharge);
398 auto makeTracks = fitTimer->get_accumulated_time();
407 std::cout <<
"make stub time " << possibleStubs <<
" create track time "
408 << makeTracks << std::endl;
413 std::cout <<
"Found " << numSeedsPerActsSeed <<
" seeds for one Acts seed"
418 std::map<const unsigned int, std::pair<std::vector<TrkrCluster*>,
419 std::vector<Acts::Vector3D>>>
421 std::map<
const unsigned int,
422 std::pair<std::vector<TrkrCluster*>,
423 std::vector<Acts::Vector3D>>> allSeeds)
425 std::map<const unsigned int, std::pair<std::vector<TrkrCluster*>,
426 std::vector<Acts::Vector3D>>> returnStub;
434 std::vector<Acts::Vector3D> mvtxClusters;
435 std::vector<TrkrCluster*> mvtxTrkrClusters;
436 double R=NAN, X0=NAN, Y0=NAN;
439 for(
const auto& [key, clusterPairs] : allSeeds)
441 auto clusterVec = clusterPairs.first;
442 auto clusterPosVec = clusterPairs.second;
445 for(
int i=0; i<clusterVec.size(); i++)
449 mvtxClusters.push_back(clusterPosVec.at(i));
450 mvtxTrkrClusters.push_back(clusterVec.at(i));
460 if(mvtxClusters.size() == 0)
461 {
return returnStub; }
463 for(
const auto& [key, clusterPairs] : allSeeds)
465 auto clusterVec = clusterPairs.first;
466 auto clusterPosVec = clusterPairs.second;
468 for(
int i=0; i<clusterVec.size(); i++)
475 double residual = sqrt( pow( globalPos(0) - X0, 2) +
476 pow( globalPos(1) - Y0, 2)) -
R;
479 std::cout <<
"Residual for cluster " << clusterVec.at(i)->getClusKey()
480 <<
" and position " << globalPos(0) <<
", "
481 << globalPos(1) <<
" is " << residual << std::endl;
484 double r = sqrt( pow(globalPos(0), 2) +
485 pow(globalPos(1), 2));
489 firstLayerBestResidual = residual;
490 firstlayerkey = clusterVec.at(i)->getClusKey();
491 firstlayerpos = clusterPosVec.at(i);
495 secondLayerBestResidual = residual;
496 secondlayerkey = clusterVec.at(i)->getClusKey();
497 secondlayerpos = clusterPosVec.at(i);
503 std::vector<TrkrCluster*> bestClusters = mvtxTrkrClusters;
504 std::vector<Acts::Vector3D> bestClusterPos = mvtxClusters;
508 bestClusterPos.push_back(firstlayerpos);
513 bestClusterPos.push_back(secondlayerpos);
518 std::cout <<
"Best cluster set is " << std::endl;
519 for(
auto& cluster : bestClusters)
520 std::cout << cluster->getClusKey() << std::endl;
524 returnStub.insert(std::make_pair(0, std::make_pair(bestClusters,bestClusterPos)));
529 std::map<const unsigned int, std::pair<std::vector<TrkrCluster*>, std::vector<Acts::Vector3D>>>
531 std::vector<Acts::Vector3D>& clusGlobPos)
534 std::vector<TrkrCluster*> mvtxClusters;
535 std::vector<TrkrCluster*> inttFirstLayerClusters;
536 std::vector<TrkrCluster*> inttSecondLayerClusters;
537 std::vector<Acts::Vector3D> mvtxClusPos;
538 std::vector<Acts::Vector3D> inttFirstLayerClusPos;
539 std::vector<Acts::Vector3D> inttSecondLayerClusPos;
541 std::map<const unsigned int, std::pair<std::vector<TrkrCluster*>,std::vector<Acts::Vector3D>>> stubs;
542 unsigned int combo = 0;
544 for(
int i = 0; i < allClusters.size(); i++)
546 const auto cluskey = allClusters.at(i)->getClusKey();
547 const auto globalPos = clusGlobPos.at(i);
551 mvtxClusters.push_back(allClusters.at(i));
552 mvtxClusPos.push_back(globalPos);
556 const double r = sqrt(pow(globalPos(0), 2) +
557 pow(globalPos(1), 2));
560 inttFirstLayerClusters.push_back(allClusters.at(i));
561 inttFirstLayerClusPos.push_back(globalPos);
565 inttSecondLayerClusters.push_back(allClusters.at(i));
566 inttSecondLayerClusPos.push_back(globalPos);
573 if(inttFirstLayerClusters.size() == 0 and
574 inttSecondLayerClusters.size() == 0)
576 stubs.insert(std::make_pair(combo, std::make_pair(mvtxClusters, mvtxClusPos)));
581 else if(inttFirstLayerClusters.size() == 1 and
582 inttSecondLayerClusters.size() == 0)
584 std::vector<TrkrCluster*> dumVec = mvtxClusters;
585 dumVec.push_back(inttFirstLayerClusters.at(0));
586 std::vector<Acts::Vector3D> dumclusVec = mvtxClusPos;
587 dumclusVec.push_back(inttFirstLayerClusPos.at(0));
588 stubs.insert(std::make_pair(combo, std::make_pair(dumVec,dumclusVec)));
590 else if(inttFirstLayerClusters.size() == 0 and
591 inttSecondLayerClusters.size() == 1)
593 std::vector<TrkrCluster*> dumVec = mvtxClusters;
594 dumVec.push_back(inttSecondLayerClusters.at(0));
595 std::vector<Acts::Vector3D> dumclusVec = mvtxClusPos;
596 dumclusVec.push_back(inttSecondLayerClusPos.at(0));
597 stubs.insert(std::make_pair(combo, std::make_pair(dumVec,dumclusVec)));
605 for(
int i=0; i<inttFirstLayerClusters.size(); i++)
607 for(
int j=0; j<inttSecondLayerClusters.size(); j++)
610 std::vector<TrkrCluster*> dumVec = mvtxClusters;
611 dumVec.push_back(inttFirstLayerClusters.at(i));
612 dumVec.push_back(inttSecondLayerClusters.at(j));
613 std::vector<Acts::Vector3D> dumclusVec = mvtxClusPos;
614 dumclusVec.push_back(inttFirstLayerClusPos.at(i));
615 dumclusVec.push_back(inttSecondLayerClusPos.at(j));
616 stubs.insert(std::make_pair(combo, std::make_pair(dumVec,dumclusVec)));
626 std::vector<Acts::Vector3D>& clusGlobPos,
627 double&
x,
double&
y,
double&
z,
628 double& px,
double& py,
double& pz)
631 for(
const auto clus : clusters)
632 std::cout <<
"Evaluating cluster : " << clus->getClusKey()
643 std::cout <<
"Circle R, X0, Y0 : " << R <<
", " << X0
644 <<
", " << Y0 << std::endl;
656 std::cout <<
"x,y circle fit : " << x <<
", "
673 double phi = atan2(-1 * (X0-x), Y0-y);
682 std::cout <<
"track seed phi : " << phi << std::endl;
691 double theta = atan(1./m);
698 std::cout <<
"Track seed theta: " << theta << std::endl;
702 float pt = 0.3 * 1.4 * R / 100.;
703 float eta = -log(tan(theta/2.));
704 float p = pt * cosh(eta);
708 px = p * sin(theta) * cos(phi);
709 py = p * sin(theta) * sin(phi);
714 std::cout <<
"Momentum vector estimate: (" << px <<
" , "
715 << py <<
", " << pz <<
") " << std::endl;
718 auto fitTimer = std::make_unique<PHTimer>(
"inttMatchTimer");
723 auto additionalClusters =
findInttMatches(clusGlobPos, R, X0, Y0, z, m);
727 for(
auto&
cluskey : additionalClusters)
731 auto addClusters = fitTimer->get_accumulated_time();
735 std::cout <<
"find intt clusters time " << addClusters << std::endl;
743 std::vector<Acts::Vector3D>& clusters,
758 for(
const auto& glob : clusters)
760 h_hits->Fill(glob(0), glob(1));
762 sqrt(pow(glob(0),2) + pow(glob(1),2)));
765 sqrt(pow(glob(0),2) + pow(glob(1),2)));
779 R, X0, Y0, xplus, yplus,
783 if(std::isnan(xplus))
787 std::cout <<
"Circle intersection calc failed, skipping"
790 <<
" and circ rad " << R <<
" with center " << X0
791 <<
", " << Y0 << std::endl;
798 const auto lastclusglob = clusters.at(clusters.size() - 1);
799 const double lastClusPhi = atan2(lastclusglob(1), lastclusglob(0));
800 const double plusPhi = atan2(yplus, xplus);
801 const double minusPhi = atan2(yminus, xminus);
803 if(fabs(lastClusPhi - plusPhi) < fabs(lastClusPhi - minusPhi))
805 xProj[
layer] = xplus;
806 yProj[
layer] = yplus;
810 xProj[
layer] = xminus;
811 yProj[
layer] = yminus;
818 h_zprojHits->Fill(zProj[layer], sqrt(pow(xProj[layer],2) +
819 pow(yProj[layer],2)));
824 std::cout <<
"Projected point is : " << xProj[
layer] <<
", "
825 << yProj[
layer] <<
", " << zProj[
layer] << std::endl;
833 std::vector<Acts::Vector3D>& clusters,
834 const double xProj[],
835 const double yProj[],
836 const double zProj[])
838 std::vector<TrkrDefs::cluskey> matchedClusters;
841 for(
int inttlayer = 0; inttlayer <
m_nInttLayers; inttlayer++)
844 const double projR = sqrt(pow(xProj[inttlayer], 2) +
845 pow(yProj[inttlayer], 2));
846 const double projPhi = atan2(yProj[inttlayer], xProj[inttlayer]);
847 const double projRphi = projR * projPhi;
849 for (
auto hitsetitr = hitsetrange.first;
850 hitsetitr != hitsetrange.second;
855 double ladderLocation[3] = {0.,0.,0.};
863 const double ladderphi = atan2(ladderLocation[1], ladderLocation[0]) + layerGeom->get_strip_phi_tilt();
864 const auto stripZSpacing = layerGeom->get_strip_z_spacing();
866 float dphi = ladderphi - projPhi;
868 { dphi -= 2. *
M_PI; }
869 else if (dphi < -1 *
M_PI)
870 { dphi += 2. *
M_PI; }
878 TVector3 projectionLocal(0,0,0);
879 TVector3 projectionGlobal(xProj[inttlayer],yProj[inttlayer],zProj[inttlayer]);
880 projectionLocal = layerGeom->get_local_from_world_coords(ladderzindex,
885 for(
auto clusIter = range.first; clusIter != range.second; ++clusIter )
887 const auto cluskey = clusIter->first;
888 const auto cluster = clusIter->second;
892 if(fabs(projectionLocal[1] - cluster->getLocalX()) <
m_rPhiSearchWin and
893 fabs(projectionLocal[2] - cluster->getLocalY()) < stripZSpacing / 2.)
896 matchedClusters.push_back(
cluskey);
900 clusters.push_back(globalPos);
903 float inttClusRphi = sqrt(pow(globalPos(0),2)+pow(globalPos(1),2)) * atan2(globalPos(1),globalPos(0));
906 h_nInttProj->Fill(projectionLocal[1] - cluster->getLocalX(),
907 projectionLocal[2] - cluster->getLocalY());
908 h_hits->Fill(globalPos(0), globalPos(1));
910 sqrt(pow(globalPos(0),2)+pow(globalPos(1),2)));
912 h_resids->Fill(zProj[inttlayer] - globalPos(2),
913 projRphi - inttClusRphi);
917 std::cout <<
"Matched INTT cluster with cluskey " <<
cluskey
918 <<
" and position " << globalPos.transpose()
919 << std::endl <<
" with projections rphi "
920 << projRphi <<
" and inttclus rphi " << inttClusRphi
921 <<
" and proj z " << zProj[inttlayer] <<
" and inttclus z "
922 << globalPos(2) <<
" in layer " << inttlayer
924 <<
" in rphi and strip z spacing " << stripZSpacing
932 return matchedClusters;
935 const double circRadius,
954 double D = layerRadius*layerRadius - circRadius*circRadius + circX0*circX0 + circY0*circY0;
955 double a = 1.0 + (circX0*circX0) / (circY0*circY0);
956 double b = - D * circX0/( circY0*circY0);
957 double c = D*D / (4.0*circY0*circY0) - layerRadius*layerRadius;
959 xplus = (-b + sqrt(b*b - 4.0* a * c) ) / (2.0 * a);
960 xminus = (-b - sqrt(b*b - 4.0* a * c) ) / (2.0 * a);
966 yplus = - (2*circX0*xplus -
D) / (2.0*circY0);
967 yminus = -(2*circX0*xminus -
D) / (2.0*circY0);
971 const double Y0,
double&
x,
987 double miny = (sqrt(pow(X0, 2) * pow(R, 2) * pow(Y0, 2) + pow(R, 2)
988 * pow(Y0,4)) + pow(X0,2) * Y0 + pow(Y0, 3))
989 / (pow(X0, 2) + pow(Y0, 2));
991 double miny2 = (-sqrt(pow(X0, 2) * pow(R, 2) * pow(Y0, 2) + pow(R, 2)
992 * pow(Y0,4)) + pow(X0,2) * Y0 + pow(Y0, 3))
993 / (pow(X0, 2) + pow(Y0, 2));
995 double minx = sqrt(pow(R, 2) - pow(miny - Y0, 2)) + X0;
996 double minx2 = -sqrt(pow(R, 2) - pow(miny2 - Y0, 2)) + X0;
999 std::cout <<
"minx1 and x2 : " << minx <<
", " << minx2 << std::endl
1000 <<
"miny1 and y2 : " << miny <<
", " << miny2 << std::endl;
1003 if(fabs(minx) < fabs(minx2))
1008 if(fabs(miny) < fabs(miny2))
1015 std::cout <<
"Minimum x and y positions " << x <<
", "
1022 const double circPhi)
1035 double trackPhi = 0;
1037 for(
const auto&
pos : globalPos)
1039 double clusPhi = atan2(
pos(1),
pos(0));
1043 if(fabs(fabs(clusPhi) -
M_PI) < 0.2)
1045 trackPhi += clusPhi;
1048 trackPhi /= globalPos.size();
1052 trackPhi -= 2. *
M_PI;
1056 for(
int i=0; i<4; i++)
1058 if(trackPhi > quadrants[i] && trackPhi <= quadrants[i+1])
1066 std::cout <<
"quadrant was not set... shouldn't be possible"
1069 if(quadrant == 1 or quadrant == 2)
1071 if(circPhi > trackPhi)
1083 if(normCircPhi > normTrackPhi)
1090 std::cout <<
"Track seed charge determined to be "
1091 << charge <<
" in quadrant " << quadrant << std::endl;
1098 double &
A,
double &
B)
1102 double xsum = 0,x2sum = 0,ysum = 0,xysum = 0;
1103 for(
const auto&
pos : globPos)
1106 double r = sqrt(pow(
pos(0),2) + pow(
pos(1), 2));
1110 x2sum=x2sum+pow(r,2);
1115 double r = sqrt(pow(
pos(0),2) + pow(
pos(1), 2));
1116 std::cout <<
" r " << r <<
" z " <<
pos(2)
1122 A = (globPos.size()*xysum-xsum*ysum) / (globPos.size()*x2sum-xsum*xsum);
1125 B = (x2sum*ysum-xsum*xysum) / (x2sum*globPos.size()-xsum*xsum);
1131 std::cout <<
"intt z_fit layer " << i <<
" is "
1140 double&
R,
double& X0,
double& Y0)
1157 int iter, IterMAX=99;
1159 double Mz, Mxy, Mxx, Myy, Mxz, Myz, Mzz, Cov_xy, Var_z;
1160 double A0, A1, A2,
A22, A3,
A33;
1162 double DET, Xcenter, Ycenter;
1169 for(
auto& globalPos : globalPositions)
1171 meanX += globalPos(0);
1172 meanY += globalPos(1);
1179 Mxx=Myy=Mxy=Mxz=Myz=Mzz=0.;
1181 for(
auto&
pos : globalPositions)
1184 double Xi =
pos(0) - meanX;
1185 double Yi =
pos(1) - meanY;
1186 double Zi = Xi * Xi + Yi * Yi;
1204 Cov_xy = Mxx * Myy - Mxy * Mxy;
1205 Var_z = Mzz - Mz * Mz;
1207 A2 = -3 * Mz * Mz - Mzz;
1208 A1 = Var_z * Mz + 4 * Cov_xy * Mz - Mxz * Mxz - Myz * Myz;
1209 A0 = Mxz * (Mxz * Myy - Myz * Mxy) +
1210 Myz * (Myz * Mxx - Mxz * Mxy) - Var_z * Cov_xy;
1214 for (x=0., y=A0, iter=0; iter<IterMAX; iter++)
1216 double Dy = A1 + x * (A22 + A33 *
x);
1217 double xnew = x - y /
Dy;
1219 double ynew = A0 + xnew * (A1 + xnew * (A2 + xnew * A3));
1220 if (fabs(ynew)>=fabs(y))
break;
1226 DET = x*x - x*Mz + Cov_xy;
1227 Xcenter = (Mxz*(Myy -
x) - Myz*Mxy)/DET/2;
1228 Ycenter = (Myz*(Mxx -
x) - Mxz*Mxy)/DET/2;
1232 X0 = Xcenter + meanX;
1233 Y0 = Ycenter + meanY;
1234 R = sqrt(Xcenter*Xcenter + Ycenter*Ycenter + Mz);
1248 auto cov = sl.covariance();
1249 float x = globalPos.x();
1250 float y = globalPos.y();
1251 float z = globalPos.z();
1252 float r = std::sqrt(x * x + y * y);
1253 float varianceRphi = cov(0,0);
1254 float varianceZ = cov(1,1);
1257 sl.referenceSurface().geometryId(), varianceRphi, varianceZ});
1260 std::cout <<
"Space point has "
1261 << x <<
", " << y <<
", " << z
1262 <<
" with variances " << cov(0,0)
1264 <<
" and cluster key "
1265 << cluskey <<
" and geo id "
1266 << sl.referenceSurface().geometryId() << std::endl;
1274 std::vector<const SpacePoint*> spVec;
1275 unsigned int numSiliconHits = 0;
1278 for (
auto hitsetitr = hitsetrange.first;
1279 hitsetitr != hitsetrange.second;
1283 for(
auto clusIter = range.first; clusIter != range.second; ++clusIter )
1285 const auto cluskey = clusIter->first;
1286 const auto cluster = clusIter->second;
1316 spVec.push_back(
sp);
1325 std::cout <<
"Total number of silicon hits to seed find with is "
1326 << numSiliconHits << std::endl;
1336 auto iter = surfMap.find(hitsetkey);
1337 if(iter != surfMap.end())
1339 return iter->second;
1395 m_surfMaps = findNode::getClass<ActsSurfaceMaps>(topNode,
"ActsSurfaceMaps");
1398 std::cout <<
PHWHERE <<
"Acts surface maps not on node tree, can't continue."
1403 m_geomContainerIntt = findNode::getClass<PHG4CylinderGeomContainer>(topNode,
"CYLINDERGEOM_INTT");
1406 std::cout <<
PHWHERE <<
"CYLINDERGEOM_INTT node not found on node tree"
1412 m_tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode,
"ActsTrackingGeometry");
1415 std::cout <<
PHWHERE <<
"No ActsTrackingGeometry on node tree. Bailing."
1422 m_clusterMap = findNode::getClass<TrkrClusterContainer>(topNode,
1423 "TRKR_CLUSTER_TRUTH");
1425 m_clusterMap = findNode::getClass<TrkrClusterContainer>(topNode,
1430 std::cout <<
PHWHERE <<
"No cluster container on the node tree. Bailing."
1434 m_hitsets = findNode::getClass<TrkrHitSetContainer>(topNode,
"TRKR_HITSET");
1437 std::cout <<
PHWHERE <<
"No hitset container on node tree. Bailing."
1452 std::cerr <<
"DST node is missing, quitting" << std::endl;
1453 throw std::runtime_error(
"Failed to find DST node in PHActsSiliconSeeding::createNodes");
1501 h_nInttProj =
new TH2F(
"nInttProj",
";l_{0}^{proj}-l_{0}^{clus} [cm]; l_{1}^{proj}-l_{1}^{clus} [cm]",
1502 10000,-10,10,10000,-50,50);
1503 h_nMvtxHits =
new TH1I(
"nMvtxHits",
";N_{MVTX}",6,0,6);
1504 h_nInttHits =
new TH1I(
"nInttHits",
";N_{INTT}",80,0,80);
1505 h_nHits =
new TH2I(
"nHits",
";N_{MVTX};N_{INTT}",10,0,10,
1507 h_nActsSeeds =
new TH1I(
"nActsSeeds",
";N_{Seeds}",400,0,400);
1508 h_nSeeds =
new TH1I(
"nActsGoodSeeds",
";N_{Seeds}",400,0,400);
1509 h_nTotSeeds =
new TH1I(
"nTotSeeds",
";N_{Seeds}",500,0,500);
1510 h_nInputMeas =
new TH1I(
"nInputMeas",
";N_{Meas}",2000,0,2000);
1515 h_hits =
new TH2F(
"hits",
";x [cm]; y [cm]",1000,-20,20,
1517 h_zhits =
new TH2F(
"zhits",
";z [cm]; r [cm]",1000,-30,30,
1519 h_projHits =
new TH2F(
"projhits",
";x [cm]; y [cm]",1000,-20,20,
1521 h_zprojHits =
new TH2F(
"zprojhits",
";z [cm]; r [cm]",1000,-30,30,
1523 h_resids =
new TH2F(
"resids",
";z_{resid} [cm]; rphi_{resid} [cm]",
1529 double returnPhi =
phi;
1531 returnPhi += 2 *
M_PI;