ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHActsSiliconSeeding.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHActsSiliconSeeding.cc
1 #include "PHActsSiliconSeeding.h"
2 
4 
7 #include <phool/getClass.h>
8 #include <phool/phool.h>
9 #include <phool/PHDataNode.h>
10 #include <phool/PHNode.h>
11 #include <phool/PHNodeIterator.h>
12 #include <phool/PHObject.h>
13 #include <phool/PHTimer.h>
14 
15 #if __cplusplus < 201402L
16 #include <boost/make_unique.hpp>
17 #endif
18 
19 #include <intt/CylinderGeomIntt.h>
20 #include <intt/InttDefs.h>
21 
24 
29 #include <trackbase/TrkrCluster.h>
31 #include <trackbase/TrkrHitSet.h>
33 #include <trackbase/TrkrDefs.h>
35 
39 #include <Acts/Seeding/Seed.hpp>
41 
43  : SubsysReco(name)
44 {}
45 
47 {
50 
52  std::make_shared<Acts::BinFinder<SpacePoint>>(Acts::BinFinder<SpacePoint>());
54  std::make_shared<Acts::BinFinder<SpacePoint>>(Acts::BinFinder<SpacePoint>());
55 
58 
59  m_seedFinderCfg.seedFilter = std::make_unique<Acts::SeedFilter<SpacePoint>>(
61 
62  if(m_seedAnalysis)
63  {
64  m_file = new TFile("seedingOutfile.root","recreate");
66  }
67 
69 }
71 {
74 
77 
79 }
80 
82 {
83  if(_n_iteration>0){
84  _iteration_map = findNode::getClass<TrkrClusterIterationMapv1>(topNode, "CLUSTER_ITERATION_MAP");
85  if (!_iteration_map){
86  std::cerr << PHWHERE << "Cluster Iteration Map missing, aborting." << std::endl;
88  }
89  }
90 
91  auto eventTimer = std::make_unique<PHTimer>("eventTimer");
92  eventTimer->stop();
93  eventTimer->restart();
94 
95  if(Verbosity() > 0)
96  std::cout << "Processing PHActsSiliconSeeding event "
97  << m_event << std::endl;
98 
99  std::vector<const SpacePoint*> spVec;
100  auto seedVector = runSeeder(spVec);
101 
102  eventTimer->stop();
103  auto seederTime = eventTimer->get_accumulated_time();
104  eventTimer->restart();
105 
106  makeSvtxTracks(seedVector);
107 
108  eventTimer->stop();
109  auto circleFitTime = eventTimer->get_accumulated_time();
110 
111  for(auto sp : spVec)
112  { delete sp; }
113  spVec.clear();
114 
115  if(Verbosity() > 0)
116  std::cout << "Finished PHActsSiliconSeeding process_event"
117  << std::endl;
118 
119  if(Verbosity() > 0)
120  {
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;
127  }
128 
129  m_event++;
130 
132 }
133 
135 {
136  if(m_seedAnalysis)
137  {
138  writeHistograms();
139  }
140 
141  if(Verbosity() > 1)
142  {
143  std::cout << "There were " << m_nBadInitialFits
144  << " bad initial circle fits" << std::endl;
145  std::cout << "There were " << m_nBadUpdates
146  << " bad second circle fits" << std::endl;
147  }
148 
150 }
151 
152 GridSeeds PHActsSiliconSeeding::runSeeder(std::vector<const SpacePoint*>& spVec)
153 {
154 
156 
158  auto covConverter = [=](const SpacePoint& sp, float, float, float)
159  -> Acts::Vector2D { return {sp.m_varianceRphi, sp.m_varianceZ};
160  };
161 
162  spVec = getMvtxSpacePoints();
163 
164  if(m_seedAnalysis)
165  { h_nInputMeas->Fill(spVec.size()); }
166 
167  std::unique_ptr<Acts::SpacePointGrid<SpacePoint>> grid =
168  Acts::SpacePointGridCreator::createGrid<SpacePoint>(m_gridCfg);
169 
170  auto spGroup = Acts::BinnedSPGroup<SpacePoint>(spVec.begin(),
171  spVec.end(),
172  covConverter,
175  std::move(grid),
177 
182  GridSeeds seedVector;
183  auto groupIt = spGroup.begin();
184  auto endGroup = spGroup.end();
185 
186  for(; !(groupIt == endGroup); ++groupIt)
187  {
188  seedVector.push_back(seedFinder.createSeedsForGroup(groupIt.bottom(),
189  groupIt.middle(),
190  groupIt.top()));
191  }
192 
193  return seedVector;
194 }
195 
197 {
198  int numSeeds = 0;
199  int numGoodSeeds = 0;
200 
202  for(auto& seeds : seedVector)
203  {
205  for(auto& seed : seeds)
206  {
207  if(Verbosity() > 1)
208  std::cout << "Seed " << numSeeds << " has "
209  << seed.sp().size() << " measurements "
210  << std::endl;
211 
212  numSeeds++;
213 
214  std::vector<TrkrCluster*> clusters;
215  std::vector<Acts::Vector3D> globalPositions;
216  ActsTransformations transformer;
217  for(auto& spacePoint : seed.sp())
218  {
219  auto cluskey = spacePoint->m_clusKey;
220  clusters.push_back(m_clusterMap->findCluster(cluskey));
221 
222  globalPositions.push_back(transformer.getGlobalPosition(
225 
226  if(Verbosity() > 1) {
227  std::cout << "Adding cluster with x,y "
228  << spacePoint->x() <<", " << spacePoint->y()
229  << " mm in detector "
231  << std::endl;
232  }
233  }
234 
235  double x = NAN, y = NAN, z = seed.z() / Acts::UnitConstants::cm;
236  double px, py, pz;
237 
238  auto fitTimer = std::make_unique<PHTimer>("trackfitTimer");
239  fitTimer->stop();
240  fitTimer->restart();
241 
244  int charge = circleFitSeed(clusters, globalPositions,
245  x, y, z,
246  px, py, pz);
247  fitTimer->stop();
248  auto circlefittime = fitTimer->get_accumulated_time();
249  fitTimer->restart();
251  if(std::isnan(x))
252  {
254  continue;
255  }
256 
257  numGoodSeeds++;
258 
260  px, py, pz, charge,
261  clusters, globalPositions);
262  fitTimer->stop();
263  auto svtxtracktime = fitTimer->get_accumulated_time();
264  if(Verbosity() > 0)
265  {
266  std::cout << "Circle fit time " << circlefittime << " and svtx time "
267  << svtxtracktime << std::endl;
268  }
269  }
270  }
271 
272  if(m_seedAnalysis)
273  {
274  h_nSeeds->Fill(numGoodSeeds);
275  h_nActsSeeds->Fill(numSeeds);
276  }
277  if(Verbosity() > 1)
278  {
279  std::cout << "Total number of seeds found in "
280  << seedVector.size() << " volume regions gives "
281  << numSeeds << " Acts seeds " << std::endl;
282  }
283 
284  return;
285 
286 }
287 
289  const double y,
290  const double z,
291  const double px,
292  const double py,
293  const double pz,
294  const int charge,
295  std::vector<TrkrCluster*>& clusters,
296  std::vector<Acts::Vector3D>& clusGlobPos)
297 {
298  auto fitTimer = std::make_unique<PHTimer>("trackfitTimer");
299  fitTimer->stop();
300  fitTimer->restart();
301 
302  auto stubs = makePossibleStubs(clusters, clusGlobPos);
303 
304  fitTimer->stop();
305  auto possibleStubs = fitTimer->get_accumulated_time();
306  fitTimer->restart();
307 
308  int numSeedsPerActsSeed = 0;
309 
310  double trackX = x;
311  double trackY = y;
312  double trackZ = z;
313  double trackPx = px;
314  double trackPy = py;
315  double trackPz = pz;
316  double trackCharge = charge;
317  double trackPhi = atan2(py,px);
318  double trackEta = atanh(pz / sqrt(px * px + py * py + pz * pz));
319 
323 
324  if(stubs.size() > 1 and m_cleanSeeds)
325  {
326  stubs = identifyBestSeed(stubs);
327  }
328 
329  for(const auto& [stub, stubClusterPairs] : stubs)
330  {
331  std::vector<TrkrCluster*> stubClusters = stubClusterPairs.first;
332  std::vector<Acts::Vector3D> stubClusterPositions = stubClusterPairs.second;
333 
334  int nMvtx = 0;
335  int nIntt = 0;
336  numSeedsPerActsSeed++;
337 
338  auto svtxTrack = std::make_unique<SvtxTrack_v2>();
339 
340  svtxTrack->set_id(m_trackMap->size());
341 
342  for(const auto clus : stubClusters)
343  {
344  const auto cluskey = clus->getClusKey();
345  svtxTrack->insert_cluster_key(cluskey);
347  { nMvtx++; }
349  { nIntt++; }
350 
351  if(Verbosity() > 2)
352  { std::cout << "Cluskey adding : " << cluskey << std::endl; }
353  }
354 
356  double R, X0, Y0;
357  circleFitByTaubin(stubClusterPositions, R, X0, Y0);
358 
361  float pt = 0.3 * 1.4 * R / 100.;
362 
363  trackPx = pt * cos(trackPhi);
364  trackPy = pt * sin(trackPhi);
365  trackPz = pt * sinh(trackEta);
366 
368  if(m_seedAnalysis)
369  {
370  h_nInttHits->Fill(nIntt);
371  h_nMvtxHits->Fill(nMvtx);
372  h_nHits->Fill(nMvtx, nIntt);
373  }
374 
375  if(Verbosity() > 1) {
376  std::cout << "Setting silicon seed with track id "
377  << m_trackMap->size()
378  << " and (x,y,z) = "
379  << trackX << ", " << trackY << ", " << trackZ
380  << std::endl << " and (px,py,pz) " << trackPx
381  << ", " << trackPy << ", " << trackPz << std::endl
382  << " with charge " << trackCharge << std::endl;
383  }
384 
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);
392 
393  m_trackMap->insert(svtxTrack.get());
394 
395  }
396 
397  fitTimer->stop();
398  auto makeTracks = fitTimer->get_accumulated_time();
399  fitTimer->restart();
400 
401 
402  if(m_seedAnalysis)
403  { h_nTotSeeds->Fill(numSeedsPerActsSeed); }
404 
405  if(Verbosity() > 0)
406  {
407  std::cout << "make stub time " << possibleStubs << " create track time "
408  << makeTracks << std::endl;
409  }
410 
411  if(Verbosity() > 1)
412  {
413  std::cout << "Found " << numSeedsPerActsSeed << " seeds for one Acts seed"
414  << std::endl;
415  }
416 }
417 
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)
424 {
425  std::map<const unsigned int, std::pair<std::vector<TrkrCluster*>,
426  std::vector<Acts::Vector3D>>> returnStub;
427 
428  double firstLayerBestResidual = std::numeric_limits<double>::max();
429  double secondLayerBestResidual = std::numeric_limits<double>::max();
432  Acts::Vector3D firstlayerpos, secondlayerpos;
433 
434  std::vector<Acts::Vector3D> mvtxClusters;
435  std::vector<TrkrCluster*> mvtxTrkrClusters;
436  double R=NAN, X0=NAN, Y0=NAN;
437 
439  for(const auto& [key, clusterPairs] : allSeeds)
440  {
441  auto clusterVec = clusterPairs.first;
442  auto clusterPosVec = clusterPairs.second;
443 
445  for(int i=0; i<clusterVec.size(); i++)
446  {
447  if(TrkrDefs::getTrkrId(clusterVec.at(i)->getClusKey()) == TrkrDefs::mvtxId)
448  {
449  mvtxClusters.push_back(clusterPosVec.at(i));
450  mvtxTrkrClusters.push_back(clusterVec.at(i));
451  }
452  }
453 
454  circleFitByTaubin(mvtxClusters, R, X0, Y0);
455 
456  break;
457  }
458 
460  if(mvtxClusters.size() == 0)
461  { return returnStub; }
462 
463  for(const auto& [key, clusterPairs] : allSeeds)
464  {
465  auto clusterVec = clusterPairs.first;
466  auto clusterPosVec = clusterPairs.second;
467 
468  for(int i=0; i<clusterVec.size(); i++)
469  {
470  if(TrkrDefs::getTrkrId(clusterVec.at(i)->getClusKey()) != TrkrDefs::inttId)
471  { continue; }
472 
473  Acts::Vector3D globalPos = clusterPosVec.at(i);
474 
475  double residual = sqrt( pow( globalPos(0) - X0, 2) +
476  pow( globalPos(1) - Y0, 2)) - R;
477 
478  if(Verbosity() > 2) {
479  std::cout << "Residual for cluster " << clusterVec.at(i)->getClusKey()
480  << " and position " << globalPos(0) << ", "
481  << globalPos(1) << " is " << residual << std::endl;
482  }
483 
484  double r = sqrt( pow(globalPos(0), 2) +
485  pow(globalPos(1), 2));
486 
487  if(r < m_nInttLayerRadii[2] and residual < firstLayerBestResidual)
488  {
489  firstLayerBestResidual = residual;
490  firstlayerkey = clusterVec.at(i)->getClusKey();
491  firstlayerpos = clusterPosVec.at(i);
492  }
493  if( r > m_nInttLayerRadii[2] and residual < secondLayerBestResidual)
494  {
495  secondLayerBestResidual = residual;
496  secondlayerkey = clusterVec.at(i)->getClusKey();
497  secondlayerpos = clusterPosVec.at(i);
498  }
499  }
500  }
501 
502  // Form the set of clusters that was identified as the smallest residuals
503  std::vector<TrkrCluster*> bestClusters = mvtxTrkrClusters;
504  std::vector<Acts::Vector3D> bestClusterPos = mvtxClusters;
506  {
507  bestClusters.push_back(m_clusterMap->findCluster(firstlayerkey));
508  bestClusterPos.push_back(firstlayerpos);
509  }
510  if(secondlayerkey < std::numeric_limits<unsigned long long>::max())
511  {
512  bestClusters.push_back(m_clusterMap->findCluster(secondlayerkey));
513  bestClusterPos.push_back(secondlayerpos);
514  }
515 
516  if(Verbosity() > 2)
517  {
518  std::cout << "Best cluster set is " << std::endl;
519  for(auto& cluster : bestClusters)
520  std::cout << cluster->getClusKey() << std::endl;
521  }
522 
524  returnStub.insert(std::make_pair(0, std::make_pair(bestClusters,bestClusterPos)));
525 
526  return returnStub;
527 }
528 
529 std::map<const unsigned int, std::pair<std::vector<TrkrCluster*>, std::vector<Acts::Vector3D>>>
530 PHActsSiliconSeeding::makePossibleStubs(std::vector<TrkrCluster*>& allClusters,
531  std::vector<Acts::Vector3D>& clusGlobPos)
532 {
533 
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;
540 
541  std::map<const unsigned int, std::pair<std::vector<TrkrCluster*>,std::vector<Acts::Vector3D>>> stubs;
542  unsigned int combo = 0;
543 
544  for(int i = 0; i < allClusters.size(); i++)
545  {
546  const auto cluskey = allClusters.at(i)->getClusKey();
547  const auto globalPos = clusGlobPos.at(i);
548 
550  {
551  mvtxClusters.push_back(allClusters.at(i));
552  mvtxClusPos.push_back(globalPos);
553  }
554  else
555  {
556  const double r = sqrt(pow(globalPos(0), 2) +
557  pow(globalPos(1), 2));
558  if(r < 8.)
559  {
560  inttFirstLayerClusters.push_back(allClusters.at(i));
561  inttFirstLayerClusPos.push_back(globalPos);
562  }
563  else
564  {
565  inttSecondLayerClusters.push_back(allClusters.at(i));
566  inttSecondLayerClusPos.push_back(globalPos);
567  }
568  }
569  }
570 
573  if(inttFirstLayerClusters.size() == 0 and
574  inttSecondLayerClusters.size() == 0)
575  {
576  stubs.insert(std::make_pair(combo, std::make_pair(mvtxClusters, mvtxClusPos)));
577  }
578 
581  else if(inttFirstLayerClusters.size() == 1 and
582  inttSecondLayerClusters.size() == 0)
583  {
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)));
589  }
590  else if(inttFirstLayerClusters.size() == 0 and
591  inttSecondLayerClusters.size() == 1)
592  {
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)));
598  }
599  else
600  {
605  for(int i=0; i<inttFirstLayerClusters.size(); i++)
606  {
607  for(int j=0; j<inttSecondLayerClusters.size(); j++)
608  {
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)));
617  combo++;
618  }
619  }
620  }
621 
622  return stubs;
623 }
624 
625 int PHActsSiliconSeeding::circleFitSeed(std::vector<TrkrCluster*>& clusters,
626  std::vector<Acts::Vector3D>& clusGlobPos,
627  double& x, double& y, double& z,
628  double& px, double& py, double& pz)
629 {
630  if(Verbosity() > 2) {
631  for(const auto clus : clusters)
632  std::cout << "Evaluating cluster : " << clus->getClusKey()
633  << std::endl;
634  }
635 
636 
639  double R, X0, Y0;
640  circleFitByTaubin(clusGlobPos,R, X0, Y0);
641 
642  if(Verbosity() > 2)
643  std::cout << "Circle R, X0, Y0 : " << R << ", " << X0
644  << ", " << Y0 << std::endl;
645 
646  findRoot(R, X0, Y0, x, y);
647 
652 
653  if(fabs(x) > m_maxSeedPCA or fabs(y) > m_maxSeedPCA)
654  {
655  if(Verbosity() > 1)
656  std::cout << "x,y circle fit : " << x << ", "
657  << y << std::endl;
658 
659  x = NAN;
660  y = NAN;
662  return 1;
663  }
664 
665  int charge = getCharge(clusGlobPos, atan2(Y0,X0));
666 
672 
673  double phi = atan2(-1 * (X0-x), Y0-y);
674  if(charge > 0)
675  {
676  phi += M_PI;
677  if(phi > M_PI)
678  phi -= 2. * M_PI;
679  }
680 
681  if(Verbosity() > 2)
682  std::cout << "track seed phi : " << phi << std::endl;
683 
684  double m, B;
685 
687 
688  lineFit(clusGlobPos, m, B);
689  z = B;
690 
691  double theta = atan(1./m);
692 
694  if(theta < 0)
695  theta += M_PI;
696 
697  if(Verbosity() > 2)
698  std::cout << "Track seed theta: " << theta << std::endl;
699 
702  float pt = 0.3 * 1.4 * R / 100.;
703  float eta = -log(tan(theta/2.));
704  float p = pt * cosh(eta);
705 
708  px = p * sin(theta) * cos(phi);
709  py = p * sin(theta) * sin(phi);
710  pz = p * cos(theta);
711 
712  if(Verbosity() > 2)
713  {
714  std::cout << "Momentum vector estimate: (" << px <<" , "
715  << py << ", " << pz << ") " << std::endl;
716  }
717 
718  auto fitTimer = std::make_unique<PHTimer>("inttMatchTimer");
719  fitTimer->stop();
720  fitTimer->restart();
721 
723  auto additionalClusters = findInttMatches(clusGlobPos, R, X0, Y0, z, m);
724 
727  for(auto& cluskey : additionalClusters)
728  { clusters.push_back(m_clusterMap->findCluster(cluskey)); }
729 
730  fitTimer->stop();
731  auto addClusters = fitTimer->get_accumulated_time();
732 
733  if(Verbosity() > 0)
734  {
735  std::cout << "find intt clusters time " << addClusters << std::endl;
736  }
737 
738  return charge;
739 
740 }
741 
742 std::vector<TrkrDefs::cluskey> PHActsSiliconSeeding::findInttMatches(
743  std::vector<Acts::Vector3D>& clusters,
744  const double R,
745  const double X0,
746  const double Y0,
747  const double B,
748  const double m)
749 {
750 
751  double xProj[m_nInttLayers];
752  double yProj[m_nInttLayers];
753  double zProj[m_nInttLayers];
754 
756  if(m_seedAnalysis)
757  {
758  for(const auto& glob : clusters)
759  {
760  h_hits->Fill(glob(0), glob(1));
761  h_zhits->Fill(glob(2),
762  sqrt(pow(glob(0),2) + pow(glob(1),2)));
763  h_projHits->Fill(glob(0), glob(1));
764  h_zprojHits->Fill(glob(2),
765  sqrt(pow(glob(0),2) + pow(glob(1),2)));
766  }
767  }
768 
769 
770 
772  for(int layer = 0; layer < m_nInttLayers; ++layer)
773  {
774  double xplus = 0;
775  double yplus = 0;
776  double xminus = 0;
777  double yminus = 0;
779  R, X0, Y0, xplus, yplus,
780  xminus, yminus);
781 
783  if(std::isnan(xplus))
784  {
785  if(Verbosity() > 2)
786  {
787  std::cout << "Circle intersection calc failed, skipping"
788  << std::endl;
789  std::cout << "layer radius " << m_nInttLayerRadii[layer]
790  << " and circ rad " << R << " with center " << X0
791  << ", " << Y0 << std::endl;
792  }
793  continue;
794  }
795 
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);
802 
803  if(fabs(lastClusPhi - plusPhi) < fabs(lastClusPhi - minusPhi))
804  {
805  xProj[layer] = xplus;
806  yProj[layer] = yplus;
807  }
808  else
809  {
810  xProj[layer] = xminus;
811  yProj[layer] = yminus;
812  }
813 
814  zProj[layer] = m * m_nInttLayerRadii[layer] + B;
815 
816  if(m_seedAnalysis) {
817  h_projHits->Fill(xProj[layer], yProj[layer]);
818  h_zprojHits->Fill(zProj[layer], sqrt(pow(xProj[layer],2) +
819  pow(yProj[layer],2)));
820  }
821 
822  if(Verbosity() > 2)
823  {
824  std::cout << "Projected point is : " << xProj[layer] << ", "
825  << yProj[layer] << ", " << zProj[layer] << std::endl;
826  }
827  }
828 
829  return matchInttClusters(clusters, xProj, yProj, zProj);
830 }
831 
832 std::vector<TrkrDefs::cluskey> PHActsSiliconSeeding::matchInttClusters(
833  std::vector<Acts::Vector3D>& clusters,
834  const double xProj[],
835  const double yProj[],
836  const double zProj[])
837 {
838  std::vector<TrkrDefs::cluskey> matchedClusters;
840 
841  for(int inttlayer = 0; inttlayer < m_nInttLayers; inttlayer++)
842  {
843  auto hitsetrange = m_hitsets->getHitSets(TrkrDefs::TrkrId::inttId, inttlayer+3);
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;
848 
849  for (auto hitsetitr = hitsetrange.first;
850  hitsetitr != hitsetrange.second;
851  ++hitsetitr)
852  {
853  const int ladderzindex = InttDefs::getLadderZId(hitsetitr->first);
854  const int ladderphiindex = InttDefs::getLadderPhiId(hitsetitr->first);
855  double ladderLocation[3] = {0.,0.,0.};
856 
857  // Add three to skip the mvtx layers for comparison
858  // to projections
859  auto layerGeom = dynamic_cast<CylinderGeomIntt*>
860  (m_geomContainerIntt->GetLayerGeom(inttlayer+3));
861 
862  layerGeom->find_segment_center(ladderzindex, ladderphiindex, ladderLocation);
863  const double ladderphi = atan2(ladderLocation[1], ladderLocation[0]) + layerGeom->get_strip_phi_tilt();
864  const auto stripZSpacing = layerGeom->get_strip_z_spacing();
865 
866  float dphi = ladderphi - projPhi;
867  if(dphi > M_PI)
868  { dphi -= 2. * M_PI; }
869  else if (dphi < -1 * M_PI)
870  { dphi += 2. * M_PI; }
871 
875  if(fabs(dphi) > 0.2)
876  { continue; }
877 
878  TVector3 projectionLocal(0,0,0);
879  TVector3 projectionGlobal(xProj[inttlayer],yProj[inttlayer],zProj[inttlayer]);
880  projectionLocal = layerGeom->get_local_from_world_coords(ladderzindex,
881  ladderphiindex,
882  projectionGlobal);
883 
884  auto range = m_clusterMap->getClusters(hitsetitr->first);
885  for(auto clusIter = range.first; clusIter != range.second; ++clusIter )
886  {
887  const auto cluskey = clusIter->first;
888  const auto cluster = clusIter->second;
889 
892  if(fabs(projectionLocal[1] - cluster->getLocalX()) < m_rPhiSearchWin and
893  fabs(projectionLocal[2] - cluster->getLocalY()) < stripZSpacing / 2.)
894  {
895 
896  matchedClusters.push_back(cluskey);
898  const auto globalPos = transform.getGlobalPosition(cluster, m_surfMaps,
899  m_tGeometry);
900  clusters.push_back(globalPos);
901 
903  float inttClusRphi = sqrt(pow(globalPos(0),2)+pow(globalPos(1),2)) * atan2(globalPos(1),globalPos(0));
904  if(m_seedAnalysis)
905  {
906  h_nInttProj->Fill(projectionLocal[1] - cluster->getLocalX(),
907  projectionLocal[2] - cluster->getLocalY());
908  h_hits->Fill(globalPos(0), globalPos(1));
909  h_zhits->Fill(globalPos(2),
910  sqrt(pow(globalPos(0),2)+pow(globalPos(1),2)));
911 
912  h_resids->Fill(zProj[inttlayer] - globalPos(2),
913  projRphi - inttClusRphi);
914  }
915  if(Verbosity() > 4)
916  {
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
923  << " with search windows " << m_rPhiSearchWin
924  << " in rphi and strip z spacing " << stripZSpacing
925  << std::endl;
926  }
927  }
928  }
929  }
930  }
931 
932  return matchedClusters;
933 }
935  const double circRadius,
936  const double circX0,
937  const double circY0,
938  double& xplus,
939  double& yplus,
940  double& xminus,
941  double& yminus)
942 {
953 
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;
958 
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);
961 
962  // both values of x are valid
963  // but for each of those values, there are two possible y values on circle 1
964  // but only one of those falls on the radical line:
965 
966  yplus = - (2*circX0*xplus - D) / (2.0*circY0);
967  yminus = -(2*circX0*xminus - D) / (2.0*circY0);
968 }
969 
970 void PHActsSiliconSeeding::findRoot(const double R, const double X0,
971  const double Y0, double& x,
972  double& y)
973 {
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));
990 
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));
994 
995  double minx = sqrt(pow(R, 2) - pow(miny - Y0, 2)) + X0;
996  double minx2 = -sqrt(pow(R, 2) - pow(miny2 - Y0, 2)) + X0;
997 
998  if(Verbosity() > 1)
999  std::cout << "minx1 and x2 : " << minx << ", " << minx2 << std::endl
1000  << "miny1 and y2 : " << miny << ", " << miny2 << std::endl;
1001 
1003  if(fabs(minx) < fabs(minx2))
1004  x = minx;
1005  else
1006  x = minx2;
1007 
1008  if(fabs(miny) < fabs(miny2))
1009  y = miny;
1010  else
1011  y = miny2;
1012 
1013  if(Verbosity() > 1)
1014  {
1015  std::cout << "Minimum x and y positions " << x << ", "
1016  << y << std::endl;
1017  }
1018 
1019 }
1020 
1021 int PHActsSiliconSeeding::getCharge(const std::vector<Acts::Vector3D>& globalPos,
1022  const double circPhi)
1023 {
1024 
1031  int charge = 0;
1032 
1035  double trackPhi = 0;
1036 
1037  for(const auto& pos : globalPos)
1038  {
1039  double clusPhi = atan2(pos(1), pos(0));
1040 
1043  if(fabs(fabs(clusPhi) - M_PI) < 0.2)
1044  clusPhi = normPhi2Pi(clusPhi);
1045  trackPhi += clusPhi;
1046  }
1047 
1048  trackPhi /= globalPos.size();
1049 
1051  if(trackPhi > M_PI)
1052  trackPhi -= 2. * M_PI;
1053 
1054  float quadrants[5] = {-M_PI,-M_PI / 2., 0, M_PI/2., M_PI};
1055  int quadrant = -1;
1056  for(int i=0; i<4; i++)
1057  {
1058  if(trackPhi > quadrants[i] && trackPhi <= quadrants[i+1])
1059  {
1060  quadrant = i;
1061  break;
1062  }
1063  }
1064 
1065  if(quadrant == -1)
1066  std::cout << "quadrant was not set... shouldn't be possible"
1067  << std::endl;
1068 
1069  if(quadrant == 1 or quadrant == 2)
1070  {
1071  if(circPhi > trackPhi)
1072  charge = -1;
1073  else
1074  charge = 1;
1075  }
1076  else
1077  {
1080  double normTrackPhi = normPhi2Pi(trackPhi);
1081  double normCircPhi = normPhi2Pi(circPhi);
1082 
1083  if(normCircPhi > normTrackPhi)
1084  charge = -1;
1085  else
1086  charge = 1;
1087  }
1088 
1089  if(Verbosity() > 1)
1090  std::cout << "Track seed charge determined to be "
1091  << charge << " in quadrant " << quadrant << std::endl;
1092 
1093  return charge;
1094 
1095 }
1096 
1097 void PHActsSiliconSeeding::lineFit(const std::vector<Acts::Vector3D>& globPos,
1098  double &A, double &B)
1099 {
1100  // copied from: https://www.bragitoff.com
1101  // we want to fit z vs radius
1102  double xsum = 0,x2sum = 0,ysum = 0,xysum = 0;
1103  for(const auto& pos : globPos)
1104  {
1105  double z = pos(2);
1106  double r = sqrt(pow(pos(0),2) + pow(pos(1), 2));
1107 
1108  xsum=xsum+r; // calculate sigma(xi)
1109  ysum=ysum+z; // calculate sigma(yi)
1110  x2sum=x2sum+pow(r,2); // calculate sigma(x^2i)
1111  xysum=xysum+r*z; // calculate sigma(xi*yi)
1112 
1113  if(Verbosity() > 4)
1114  {
1115  double r = sqrt(pow(pos(0),2) + pow(pos(1), 2));
1116  std::cout << " r " << r << " z " << pos(2)
1117  << std::endl;
1118  }
1119  }
1120 
1122  A = (globPos.size()*xysum-xsum*ysum) / (globPos.size()*x2sum-xsum*xsum);
1123 
1125  B = (x2sum*ysum-xsum*xysum) / (x2sum*globPos.size()-xsum*xsum);
1126 
1127  if(Verbosity() > 4)
1128  {
1129  for(int i =0; i <m_nInttLayers; i++)
1130  {
1131  std::cout << "intt z_fit layer " << i << " is "
1132  << A * m_nInttLayerRadii[i] + B << std::endl;
1133  }
1134  }
1135 
1136  return;
1137 }
1138 
1139 void PHActsSiliconSeeding::circleFitByTaubin(const std::vector<Acts::Vector3D>& globalPositions,
1140  double& R, double& X0, double& Y0)
1141 {
1157  int iter, IterMAX=99;
1158 
1159  double Mz, Mxy, Mxx, Myy, Mxz, Myz, Mzz, Cov_xy, Var_z;
1160  double A0, A1, A2, A22, A3, A33;
1161  double x, y;
1162  double DET, Xcenter, Ycenter;
1163 
1164  // Compute x- and y- sample means
1165  double meanX = 0;
1166  double meanY = 0;
1167  double weight = 0;
1168 
1169  for(auto& globalPos : globalPositions)
1170  {
1171  meanX += globalPos(0);
1172  meanY += globalPos(1);
1173  weight++;
1174  }
1175 
1176  meanX /= weight;
1177  meanY /= weight;
1178 
1179  Mxx=Myy=Mxy=Mxz=Myz=Mzz=0.;
1180 
1181  for(auto& pos : globalPositions)
1182  {
1183 
1184  double Xi = pos(0) - meanX;
1185  double Yi = pos(1) - meanY;
1186  double Zi = Xi * Xi + Yi * Yi;
1187 
1188  Mxy += Xi*Yi;
1189  Mxx += Xi*Xi;
1190  Myy += Yi*Yi;
1191  Mxz += Xi*Zi;
1192  Myz += Yi*Zi;
1193  Mzz += Zi*Zi;
1194  }
1195 
1196  Mxx /= weight;
1197  Myy /= weight;
1198  Mxy /= weight;
1199  Mxz /= weight;
1200  Myz /= weight;
1201  Mzz /= weight;
1202 
1203  Mz = Mxx + Myy;
1204  Cov_xy = Mxx * Myy - Mxy * Mxy;
1205  Var_z = Mzz - Mz * Mz;
1206  A3 = 4 * 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;
1211  A22 = A2 + A2;
1212  A33 = A3 + A3 + A3;
1213 
1214  for (x=0., y=A0, iter=0; iter<IterMAX; iter++) // usually, 4-6 iterations are enough
1215  {
1216  double Dy = A1 + x * (A22 + A33 * x);
1217  double xnew = x - y / Dy;
1218  if ((xnew == x)||(!std::isfinite(xnew))) break;
1219  double ynew = A0 + xnew * (A1 + xnew * (A2 + xnew * A3));
1220  if (fabs(ynew)>=fabs(y)) break;
1221  x = xnew; y = ynew;
1222  }
1223 
1224  // computing parameters of the fitting circle
1225 
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;
1229 
1230  // assembling the output
1231 
1232  X0 = Xcenter + meanX;
1233  Y0 = Ycenter + meanY;
1234  R = sqrt(Xcenter*Xcenter + Ycenter*Ycenter + Mz);
1235 
1236 }
1237 
1239  const SourceLink& sl)
1240 {
1241  Acts::Vector2D localPos(sl.location()(0), sl.location()(1));
1242  Acts::Vector3D globalPos(0,0,0);
1243  Acts::Vector3D mom(1,1,1);
1244 
1245  globalPos = sl.referenceSurface().localToGlobal(m_tGeometry->geoContext,
1246  localPos, mom);
1247 
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);
1255 
1256  SpacePointPtr spPtr(new SpacePoint{cluskey, x, y, z, r,
1257  sl.referenceSurface().geometryId(), varianceRphi, varianceZ});
1258 
1259  if(Verbosity() > 2)
1260  std::cout << "Space point has "
1261  << x << ", " << y << ", " << z
1262  << " with variances " << cov(0,0)
1263  << ", " << cov(1,1)
1264  << " and cluster key "
1265  << cluskey << " and geo id "
1266  << sl.referenceSurface().geometryId() << std::endl;
1267 
1268  return spPtr;
1269 
1270 }
1271 
1272 std::vector<const SpacePoint*> PHActsSiliconSeeding::getMvtxSpacePoints()
1273 {
1274  std::vector<const SpacePoint*> spVec;
1275  unsigned int numSiliconHits = 0;
1276 
1277  auto hitsetrange = m_hitsets->getHitSets(TrkrDefs::TrkrId::mvtxId);
1278  for (auto hitsetitr = hitsetrange.first;
1279  hitsetitr != hitsetrange.second;
1280  ++hitsetitr)
1281  {
1282  auto range = m_clusterMap->getClusters(hitsetitr->first);
1283  for( auto clusIter = range.first; clusIter != range.second; ++clusIter )
1284  {
1285  const auto cluskey = clusIter->first;
1286  const auto cluster = clusIter->second;
1287  if(_iteration_map != NULL && _n_iteration >0 ){
1288  // std::cout << "map exists entries: " << _iteration_map->size() << std::endl;
1290  continue; // skip hits used in a previous iteration
1291  }
1292  }
1293 
1295  const auto surface = getSurface(hitsetkey);
1296  if(!surface)
1297  continue;
1298 
1299  Acts::BoundVector loc = Acts::BoundVector::Zero();
1300  loc[Acts::eBoundLoc0] = cluster->getLocalX() * Acts::UnitConstants::cm;
1301  loc[Acts::eBoundLoc1] = cluster->getLocalY() * Acts::UnitConstants::cm;
1302 
1303  Acts::BoundMatrix cov = Acts::BoundMatrix::Zero();
1305  cluster->getActsLocalError(0,0) * Acts::UnitConstants::cm2;
1307  cluster->getActsLocalError(0,1) * Acts::UnitConstants::cm2;
1309  cluster->getActsLocalError(1,0) * Acts::UnitConstants::cm2;
1311  cluster->getActsLocalError(1,1) * Acts::UnitConstants::cm2;
1312 
1313  SourceLink sl(cluskey, surface, loc, cov);
1314 
1315  auto sp = makeSpacePoint(cluskey, sl).release();
1316  spVec.push_back(sp);
1317  numSiliconHits++;
1318  }
1319  }
1320 
1321  if(m_seedAnalysis)
1322  { h_nInputMvtxMeas->Fill(numSiliconHits); }
1323 
1324  if(Verbosity() > 1)
1325  std::cout << "Total number of silicon hits to seed find with is "
1326  << numSiliconHits << std::endl;
1327 
1328 
1329  return spVec;
1330 }
1332 {
1335  auto surfMap = m_surfMaps->siliconSurfaceMap;
1336  auto iter = surfMap.find(hitsetkey);
1337  if(iter != surfMap.end())
1338  {
1339  return iter->second;
1340  }
1341 
1343  return nullptr;
1344 
1345 }
1347 {
1349 
1350  config.bFieldInZ = m_bField;
1351  config.minPt = m_minSeedPt / m_gridFactor;
1352  config.rMax = m_rMax;
1353  config.zMax = m_zMax;
1354  config.zMin = m_zMin;
1355  config.deltaRMax = m_deltaRMax;
1356  config.cotThetaMax = m_cotThetaMax;
1357 
1358  return config;
1359 }
1360 
1362 {
1364 
1366  config.rMax = m_rMax;
1367  config.rMin = m_rMin;
1368  config.zMin = m_zMin;
1369  config.zMax = m_zMax;
1370 
1372  config.deltaRMin = 1.;
1373  config.deltaRMax = m_deltaRMax;
1374 
1376  config.collisionRegionMin = -300.;
1377  config.collisionRegionMax = 300.;
1378  config.sigmaScattering = 5.;
1380  config.cotThetaMax = m_cotThetaMax;
1381  config.minPt = m_minSeedPt;
1382  config.bFieldInZ = m_bField;
1383 
1385  config.radLengthPerSeed = 0.05;
1386 
1388  config.impactMax = 20;
1389 
1390  return config;
1391 }
1392 
1394 {
1395  m_surfMaps = findNode::getClass<ActsSurfaceMaps>(topNode, "ActsSurfaceMaps");
1396  if(!m_surfMaps)
1397  {
1398  std::cout << PHWHERE << "Acts surface maps not on node tree, can't continue."
1399  << std::endl;
1401  }
1402 
1403  m_geomContainerIntt = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_INTT");
1404  if(!m_geomContainerIntt)
1405  {
1406  std::cout << PHWHERE << "CYLINDERGEOM_INTT node not found on node tree"
1407  << std::endl;
1409  }
1410 
1411 
1412  m_tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode, "ActsTrackingGeometry");
1413  if(!m_tGeometry)
1414  {
1415  std::cout << PHWHERE << "No ActsTrackingGeometry on node tree. Bailing."
1416  << std::endl;
1418  }
1419 
1420 
1421  if(m_useTruthClusters)
1422  m_clusterMap = findNode::getClass<TrkrClusterContainer>(topNode,
1423  "TRKR_CLUSTER_TRUTH");
1424  else
1425  m_clusterMap = findNode::getClass<TrkrClusterContainer>(topNode,
1426  "TRKR_CLUSTER");
1427 
1428  if(!m_clusterMap)
1429  {
1430  std::cout << PHWHERE << "No cluster container on the node tree. Bailing."
1431  << std::endl;
1433  }
1434  m_hitsets = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
1435  if(!m_hitsets)
1436  {
1437  std::cout << PHWHERE << "No hitset container on node tree. Bailing."
1438  << std::endl;
1440  }
1442 }
1444 {
1445 
1446  PHNodeIterator iter(topNode);
1447 
1448  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
1449 
1450  if (!dstNode)
1451  {
1452  std::cerr << "DST node is missing, quitting" << std::endl;
1453  throw std::runtime_error("Failed to find DST node in PHActsSiliconSeeding::createNodes");
1454  }
1455 
1456  PHCompositeNode *svtxNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "SVTX"));
1457 
1458  if (!svtxNode)
1459  {
1460  svtxNode = new PHCompositeNode("SVTX");
1461  dstNode->addNode(svtxNode);
1462  }
1463 
1464  m_trackMap = findNode::getClass<SvtxTrackMap>(topNode,_track_map_name);
1465  if(!m_trackMap)
1466  {
1468  PHIODataNode<PHObject> *trackNode =
1470  svtxNode->addNode(trackNode);
1471 
1472  }
1473 
1475 }
1476 
1478 {
1479  m_file->cd();
1480  h_nInttProj->Write();
1481  h_nMvtxHits->Write();
1482  h_nSeeds->Write();
1483  h_nActsSeeds->Write();
1484  h_nTotSeeds->Write();
1485  h_nInttHits->Write();
1486  h_nInputMeas->Write();
1487  h_nHits->Write();
1488  h_nInputMvtxMeas->Write();
1489  h_nInputInttMeas->Write();
1490  h_hits->Write();
1491  h_zhits->Write();
1492  h_projHits->Write();
1493  h_zprojHits->Write();
1494  h_resids->Write();
1495  m_file->Write();
1496  m_file->Close();
1497 }
1498 
1500 {
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,
1506  80,0,80);
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);
1511  h_nInputMvtxMeas = new TH1I("nInputMvtxMeas",";N_{meas}^{mvtx}",
1512  150,0,150);
1513  h_nInputInttMeas = new TH1I("nInputInttMeas",";N_{meas}^{intt}",
1514  150,0,150);
1515  h_hits = new TH2F("hits",";x [cm]; y [cm]",1000,-20,20,
1516  1000,-20,20);
1517  h_zhits = new TH2F("zhits",";z [cm]; r [cm]",1000,-30,30,
1518  1000,-30,30);
1519  h_projHits = new TH2F("projhits",";x [cm]; y [cm]",1000,-20,20,
1520  1000,-20,20);
1521  h_zprojHits = new TH2F("zprojhits",";z [cm]; r [cm]",1000,-30,30,
1522  1000,-30,30);
1523  h_resids = new TH2F("resids",";z_{resid} [cm]; rphi_{resid} [cm]",
1524  100,-1,1,100,-1,1);
1525 }
1526 
1528 {
1529  double returnPhi = phi;
1530  if(returnPhi < 0)
1531  returnPhi += 2 * M_PI;
1532  return returnPhi;
1533 }
1534 
1535 
1537 {
1538  if(!spacing)
1539  {
1540  m_gridFactor = 1.;
1541  m_rMax = 50.;
1542  m_cotThetaMax = 1.335647;
1543  m_maxSeedPCA = 0.1;
1544  }
1545 
1546 }