28 #include <nlohmann/json.hpp>
34 nearZero(
const double &
v){
44 static inline std::string
to_string(
const std::vector<T>&
v)
46 std::ostringstream oss;
50 std::copy(
v.begin(),
v.end() - 1, std::ostream_iterator<T>(oss,
", "));
64 declareInterface<IFPGATrackSimRoadFinderTool>(
this);
71 const std::vector<Gaudi::Details::PropertyBase*> props = this->getProperties();
72 for( Gaudi::Details::PropertyBase* prop : props ) {
73 if (prop->ownerTypeName()==this->type()) {
74 ATH_MSG_DEBUG(
"Property:\t" << prop->name() <<
"\t : \t" << prop->toString());
88 if (nLogicalLayers == 0){
89 ATH_MSG_ERROR(
"Number of logical layers is zero in FPGATrackSimGenScanTool::initialize");
90 return StatusCode::FAILURE;
118 return StatusCode::FAILURE;
127 return StatusCode::SUCCESS;
134 std::vector<std::shared_ptr<const FPGATrackSimRoad>> &roads)
163 std::vector<HitPairSet> pairsets;
181 if (pairsets.size() >1) {
185 s +=
"(" + std::to_string(hit->
layer) +
"," + std::to_string(hit->
hitptr->getR()) +
"), ";
196 for (std::unique_ptr<FPGATrackSimRoad>&
r :
m_roads) {
197 const std::vector<std::vector<std::shared_ptr<const FPGATrackSimHit>>>& theseHits =
r->getAllHits();
201 std::vector<std::vector<std::shared_ptr<const FPGATrackSimHit>>>
vec(5);
202 for (
size_t ihit = 0; ihit < toUse.size(); ++ihit) {
203 unsigned int layer = toUse[ihit];
204 if (
layer >= theseHits.size() || theseHits[
layer].empty()) {
205 ATH_MSG_ERROR(
"Hit index out of range in keepHitsStrategy: layer=" <<
layer <<
", hits.size()=" << theseHits.size());
206 return StatusCode::FAILURE;
208 vec[ihit].push_back(theseHits[
layer][0]);
210 r->setHits(std::move(
vec));
215 for (
auto &
r :
m_roads) roads.push_back(std::move(
r));
223 return StatusCode::SUCCESS;
228 std::vector<HitPairSet> &output_pairsets)
233 std::vector<std::vector<const StoredHit *>> hitsByLayer(
m_binnedhits->getNLayers());
253 if (passedPairFilter)
256 std::vector<HitPairSet> pairsets;
266 output_pairsets.push_back(pairset);
275 output_pairsets.push_back(std::move(
pairs));
280 output_pairsets.push_back(std::move(filteredpairs));
283 return StatusCode::SUCCESS;
290 const std::vector<const StoredHit *>& newhits)
294 std::vector<bool> pairset_used(inputstate.
pairsets.size(),
false);
296 for (
auto &newhit : newhits) {
299 std::set<const StoredHit *> vetoList;
302 for (
unsigned ps_idx = 0; ps_idx < inputstate.
pairsets.size(); ++ps_idx) {
303 auto &pairset = inputstate.
pairsets[ps_idx];
308 pairset_used[ps_idx]=
true;
309 outputstate.
pairsets.push_back(newset);
311 for (
auto vetohit : pairset.hitlist) {
312 vetoList.insert(vetohit);
319 if (vetoList.count(prevhit) == 0) {
324 outputstate.
pairsets.push_back(std::move(newset));
331 if (lyridx <= allowed_missed_hits) {
340 for (
unsigned ps_idx = 0; ps_idx < inputstate.
pairsets.size(); ++ps_idx) {
341 auto &pairset = inputstate.
pairsets[ps_idx];
342 if ((!pairset_used[ps_idx])&&(lyridx < (pairset.hitlist.size() + allowed_missed_hits))) {
343 outputstate.
pairsets.push_back(pairset);
359 if (lyridx > allowed_missed_hits) {
372 std::vector<HitPairSet> &output_pairsets)
377 std::vector<std::vector<const StoredHit *>> hitsByLayer(
m_binnedhits->getNLayers());
385 for (
unsigned lyridx = 0; lyridx <
m_binnedhits->getNLayers(); lyridx++) {
395 return StatusCode::SUCCESS;
401 std::vector<std::vector<const StoredHit *>>& hitsByLayer)
406 hitsByLayer.at(hit.
layer).push_back(&hit);
409 return StatusCode::SUCCESS;
422 std::vector<const StoredHit *>
const * lastlyr = 0;
423 std::vector<const StoredHit *>
const * lastlastlyr = 0;
430 for (
const StoredHit *
const &ptr2 : *lastlyr) {
435 for (
const StoredHit *
const &ptr2 : *lastlastlyr) {
441 lastlastlyr = lastlyr;
442 lastlyr = &hitsByLayer[lyr];
446 return StatusCode::SUCCESS;
454 int lyr =
std::min(pair.first->layer,pair.second->layer);
472 return StatusCode::SUCCESS;
478 std::vector<HitPairSet> &pairsets,
488 if ((std::abs(
int(pair.second->layer) -
int(pair.first->layer)) > 1)
489 && (pairset.hasLayer(
std::min(pair.first->layer,pair.second->layer) + 1)))
503 int size = pairset.addPair(pair);
505 ATH_MSG_VERBOSE(
"addPair " << pairsets.size() <<
" " << pairset.pairList.size() <<
" " <<
size);
516 pairsets.push_back(std::move(newpairset));
520 return StatusCode::SUCCESS;
535 cutvar(std::string
name,
double val,
double cut, std::vector<TH1D *>& histset) :
537 bool passed() {
return std::abs(m_val) < m_cut; }
538 void fill(
unsigned cat) { m_histset[cat]->Fill(m_val); }
542 std::vector<TH1D *> &m_histset;
546 std::vector<cutvar> allcuts;
547 allcuts.push_back(cutvar(
"MatchPhi", pairset.
MatchPhi(pair),
550 allcuts.push_back(cutvar(
"MatchEta", pairset.
MatchEta(pair),
553 allcuts.push_back(cutvar(
"DeltaDeltaPhi", pairset.
DeltaDeltaPhi(pair),
556 allcuts.push_back(cutvar(
"DeltaDeltaEta", pairset.
DeltaDeltaEta(pair),
559 allcuts.push_back(cutvar(
"PhiCurvature", pairset.
PhiCurvature(pair),
562 allcuts.push_back(cutvar(
"EtaCurvature", pairset.
EtaCurvature(pair),
566 allcuts.push_back(cutvar(
569 allcuts.push_back(cutvar(
573 allcuts.push_back(cutvar(
576 allcuts.push_back(cutvar(
581 unsigned passedCuts = std::count_if(allcuts.begin(), allcuts.end(),
582 [](cutvar&
cut) { return cut.passed(); });
583 bool passedAll = (passedCuts == allcuts.size());
586 bool passedAllButOne = (passedCuts == allcuts.size() - 1);
587 for (cutvar&
cut: allcuts) {
591 (passedAll || (passedAllButOne && !
cut.passed())));
597 for (cutvar &
cut : allcuts)
599 s +=
cut.m_name +
" : (" +
cut.passed() +
", " +
cut.m_val +
"), ";
603 << *pairset.
lastpair().second <<
"\n "
604 << *pair.first <<
"\n "
615 std::vector<std::vector<std::shared_ptr<const FPGATrackSimHit>>>
616 sorted_hits(
m_binnedhits->getNLayers(),std::vector<std::shared_ptr<const FPGATrackSimHit>>());
619 hitLayers |= 1 << hit->layer;
620 sorted_hits[hit->layer].push_back(hit->hitptr);
629 m_roads.emplace_back(std::make_unique<FPGATrackSimRoad>());
634 r->setHits(std::move(sorted_hits));
642 r->setHitLayers(hitLayers);
646 r->setFitParams(fittedpars);
665 if (pairList.size() > 0)
668 LastPhiCurvature = PhiCurvature(pair);
669 LastEtaCurvature = EtaCurvature(pair);
672 pairList.push_back(pair);
674 hitLayers |= (0x1 << pair.first->layer);
675 hitLayers |= (0x1 << pair.second->layer);
676 if (!hasHit(pair.first))
677 hitlist.push_back(pair.first);
678 if (!hasHit(pair.second))
679 hitlist.push_back(pair.second);
680 return this->pairList.size();
686 if ((lastpair().
first->hitptr==pair.first->hitptr)||
687 (lastpair().
first->hitptr==pair.second->hitptr)||
688 (lastpair().
second->hitptr==pair.first->hitptr)||
689 (lastpair().
second->hitptr==pair.second->hitptr))
692 double dr = (lastpair().second->hitptr->getR() - pair.first->hitptr->getR()) / 2.0;
693 double lastpairextrap = lastpair().second->phiShift - lastpair().dPhi() / lastpair().dR() *
dr;
694 double newpairextrap = pair.first->phiShift + pair.
dPhi() / pair.
dR() *
dr;
695 return lastpairextrap - newpairextrap;
700 if ((lastpair().
first->hitptr==pair.first->hitptr)||
701 (lastpair().
first->hitptr==pair.second->hitptr)||
702 (lastpair().
second->hitptr==pair.first->hitptr)||
703 (lastpair().
second->hitptr==pair.second->hitptr))
706 double dr = (lastpair().second->hitptr->getR() - pair.first->hitptr->getR()) / 2.0;
707 double lastpairextrap = lastpair().second->etaShift - lastpair().dEta() / lastpair().dR() *
dr;
708 double newpairextrap = pair.first->etaShift + pair.
dEta() / pair.
dR() *
dr;
709 return lastpairextrap - newpairextrap;
714 return (pair.
dPhi() * lastpair().dR() - lastpair().
dPhi() * pair.
dR()) / (lastpair().dR() * pair.
dR());
719 return (pair.
dEta() * lastpair().dR() - lastpair().
dEta() * pair.
dR()) / (lastpair().dR() * pair.
dR());
724 return 2 * DeltaDeltaPhi(pair) / (lastpair().dR() + pair.
dR());
728 return 2 * DeltaDeltaEta(pair) / (lastpair().dR() + pair.
dR());
732 return PhiCurvature(pair) - LastPhiCurvature;
736 return EtaCurvature(pair) - LastEtaCurvature;
741 return lastpair().PhiInExtrap(r_in) + 0.5 * PhiCurvature(pair) * (r_in -
r) * (r_in -
r);
746 return pair.
PhiOutExtrap(r_out) + 0.5 * PhiCurvature(pair) * (r_out -
r) * (r_out -
r);
754 double N =
hits.size();
758 double sum_PhiR2 = 0;
771 double r = hit->hitptr->getR();
775 double dphi = hit->phiShift;
776 double dphi2 = dphi*dphi;
777 double deta = hit->etaShift;
778 double deta2 = deta*deta;
783 sum_PhiR2 += dphi*r2;
799 double r6_t0 = (-sum_R2 * sum_R4 + sum_R3*sum_R3);
800 double r5_t0 = (sum_R*sum_R4 - sum_R2*sum_R3);
801 double r4_t0 = (-sum_R * sum_R3 + sum_R2*sum_R2);
802 double r4_t1 = (-
N*sum_R4 + sum_R2*sum_R2);
803 double r3_t0 = (
N*sum_R3 - sum_R*sum_R2);
804 double r2_t0 = (-
N*sum_R2 + sum_R*sum_R);
807 const double denom_phi =
N * r6_t0 + sum_R * r5_t0 + sum_R2 * r4_t0;
808 if (nearZero(denom_phi)){
809 ATH_MSG_ERROR(
"Divide by zero (phi) trapped in FPGATrackSimGenScanTool::fitRoad");
814 std::vector<double> phivars(3);
815 phivars[0] = (sum_Phi*r6_t0 + sum_PhiR*r5_t0 + sum_PhiR2*r4_t0)/denom_phi;
816 phivars[1] = (sum_Phi*r5_t0 + sum_PhiR*r4_t1 + sum_PhiR2*r3_t0)/denom_phi;
817 phivars[2] = (sum_Phi*r4_t0 + sum_PhiR*r3_t0 + sum_PhiR2*r2_t0)/denom_phi;
821 const double denom_eta =
N*sum_R2 - sum_R*sum_R;
822 if (nearZero(denom_eta)){
823 ATH_MSG_ERROR(
"Divide by zero (eta) trapped in FPGATrackSimGenScanTool::fitRoad");
827 std::vector<double> etavars(2);
828 etavars[0] = (-sum_R*sum_EtaR + sum_R2*sum_Eta)/denom_eta;
829 etavars[1] = (
N*sum_EtaR - sum_R*sum_Eta)/denom_eta;
834 parshift[0] = etavars[0] + etavars[1]*
m_rin;
835 parshift[1]= etavars[0] + etavars[1]*
m_rout;
836 parshift[2]= -(phivars[0]/
m_rin + phivars[1] + phivars[2]*
m_rin) ;
837 parshift[3]= -(phivars[0]/
m_rout + phivars[1] + phivars[2]*
m_rout) ;
839 parshift[4]= -phivars[2]/4.0*
y*
y ;
842 const auto& lastStep =
m_binnedhits->getBinTool().lastStep();
844 parset[
par]+=parshift[
par];
845 inBin = inBin && (std::abs(parshift[
par]) < lastStep->binWidth(
par)/2.0);
848 double ec0 = etavars[0];
849 double ec1 = etavars[1];
850 double eta_chi2 = sum_Eta2 - 2*ec0*sum_Eta - 2*ec1*sum_EtaR +
N*ec0*ec0 + 2*ec0*ec1*sum_R + ec1*ec1*sum_R2;
852 double pc0 = phivars[0];
853 double pc1 = phivars[1];
854 double pc2 = phivars[2];
856 double phi_chi2 = sum_Phi2 - 2*pc0*sum_Phi - 2*pc1*sum_PhiR - 2*pc2*sum_PhiR2 +
N*pc0*pc0 + 2*pc0*pc1*sum_R + 2*pc0*pc2*sum_R2 + 2*pc1*pc2*sum_R3 + pc1*pc1*sum_R2 + pc2*pc2*sum_R4;
860 double r = hit->hitptr->getR();
861 ATH_MSG_VERBOSE(
"Fitted r= " <<
r <<
" phishift " << hit->phiShift <<
" =?= " << pc0+pc1*
r+pc2*
r*
r <<
" etashift " << hit->etaShift <<
" =?= " << ec0+ec1*
r);
865 ATH_MSG_DEBUG(
"Fitted parset inBin="<< inBin <<
" nhits=" <<
hits.size() <<
" " << parset <<
" chi2 " << eta_chi2 <<
"," << phi_chi2);
868 trackpars =
m_binnedhits->getBinTool().binDesc()->parSetToTrackPars(parset);
880 std::vector<unsigned> toUse;
884 if (hitmask == 0x1f) {
887 else if (hitmask == 0x1e) {
890 else if (hitmask == 0x1d) {
893 else if (hitmask == 0x1b) {
896 else if (hitmask == 0x17) {
899 else if (hitmask == 0x0f) {
906 if (hitmask == 0x1f) {
909 else if (hitmask == 0x1e) {
912 else if (hitmask == 0x1d) {
915 else if (hitmask == 0x1b) {
918 else if (hitmask == 0x17) {
921 else if (hitmask == 0x0f) {
928 if (hitmask == 0x1f) {
931 else if (hitmask == 0x1e) {
934 else if (hitmask == 0x1d) {
937 else if (hitmask == 0x1b) {
940 else if (hitmask == 0x17) {
943 else if (hitmask == 0x0f) {
950 if (hitmask == 0x1f) {
953 else if (hitmask == 0x1e) {
956 else if (hitmask == 0x1d) {
959 else if (hitmask == 0x1b) {
962 else if (hitmask == 0x17) {
965 else if (hitmask == 0x0f) {