PndApolloniusTriplet.h

Go to the documentation of this file.

#pragma once

#define _USE_MATH_DEFINES

#include <math.h>
#include <tuple>

#include "PndSttGeometryMap.h"
#include "FairHit.h"
#include "PndSttHit.h"
#include "FairLogger.h"
#include "PndHoughApollonius.h"
#include "PndTrackEvaluatorDetStt.h"
#include "PndSttCA.h"
// general
using std::cout;
using std::endl;

namespace PndApollonius {

struct Thresholds {
  double fDistanceThresholdSTTSkewed = 0.8;
  double fDistanceThresholdSTTCombinedSkewed = 1.;

  double fDistanceThresholdSTTNarrow = 0.3;
  double fDistanceThresholdSTTFar = 0.5;

  double fDistanceThresholdMVDNarrow = 0.5;
  double fDistanceThresholdMVDMid1 = 0.7;
  double fDistanceThresholdMVDMid2 = 1.;
  double fDistanceThresholdMVDFar = 2.;
  double fDistanceThresholdGEM = 1.;

  double fSTTArcLengthCut = 26.;
  double fMVDArcLengthCutNarrow = 5.;
  double fMVDArcLengthCutMid = 7.;
  double fMVDArcLengthCutFar = 10.;
};

struct Triplet {
  std::array<PndSttHit *, 3> fTripletHits;
  friend std::ostream &operator<<(std::ostream &output, const Triplet &t)

  {
    output << "Triplet: ";
    for (auto hit : t.fTripletHits)
      output << hit->GetTubeID() << "(" << hit->GetX() << "," << hit->GetY() << ") /";
    output << std::endl;

    return output;
  }
};

struct TripletSolution {

  enum detID { MVDpixel, MVDstrip, STT, GEM, NOTDEFINED };
  std::vector<std::string> detID_Names{"MVD_Pixel", "MVD_Strip", "STT", "GEM", "NOTDEFINED"};

  TripletSolution(){};

  TripletSolution(TVector3 track) : fTrack(track){};

  void AddHits(detID detector, std::vector<FairHit *> data)
  {
    fHits[detector] = data;
    fAllHits.insert(fAllHits.end(), data.begin(), data.end());
  }

  void AddHits(std::vector<PndSttHit *> sttHits)
  {
    fHits[STT].insert(fHits[STT].end(), sttHits.begin(), sttHits.end());
    fAllHits.insert(fAllHits.end(), sttHits.begin(), sttHits.end());
  }

  void AddHit(detID detector, FairHit *data)
  {
    if (fHits.count(detector) == 0) {
      fHits[detector];
    }
    if (std::find(fHits[detector].begin(), fHits[detector].end(), data) == fHits[detector].end()) { // add only if hit is not in dataset
      fHits[detector].push_back(data);
      fAllHits.push_back(data);
    }
  }

  int GetNHits() { return fHits[MVDpixel].size() + fHits[MVDstrip].size() + fHits[STT].size() + fHits[GEM].size(); }

  TVector2 HitOnTrack(FairHit *hit) { return TVector2(hit->GetX() - fTrack.X(), hit->GetY() - fTrack.Y()); }

  void IsClockwise();

  void SortHits(TVector2 &firstHit, detID detector);

  void SortStt(PndSttGeometryMap *fGeometryMap);

  FairHit *GetFirstHit(PndSttGeometryMap *fGeometryMap);

  void SortAllHits(PndSttGeometryMap *fGeometryMap);

  bool IsReasonableHit(FairHit *CurrentHit, FairHit *PriviousHit);

  friend std::ostream &operator<<(std::ostream &output, const TripletSolution &t)
  {
    output << "Track x/y/r " << t.fTrack.X() << "/" << t.fTrack.Y() << "/" << t.fTrack.Z() << " ";
    if (t.fTriplet.fTripletHits[0] != nullptr) {
      output << " trip: " << t.fTriplet.fTripletHits[0]->GetTubeID() << "/" << t.fTriplet.fTripletHits[1]->GetTubeID() << "/" << t.fTriplet.fTripletHits[2]->GetTubeID() << " ";
    }
    if (t.fHits.count(TripletSolution::detID::MVDpixel) > 0) {
      output << t.detID_Names.at(TripletSolution::detID::MVDpixel) << " : ";
      for (auto hit : t.fHits.at(TripletSolution::detID::MVDpixel)) {
        output << hit->GetEntryNr() << "/";
      }
      output << " ";
    }
    if (t.fHits.count(TripletSolution::detID::MVDstrip) > 0) {
      output << t.detID_Names.at(TripletSolution::detID::MVDstrip) << " : ";
      for (auto hit : t.fHits.at(TripletSolution::detID::MVDstrip)) {
        output << hit->GetEntryNr() << "/";
      }
      output << " ";
    }
    if (t.fHits.count(TripletSolution::detID::STT) > 0) {
      output << t.detID_Names.at(TripletSolution::detID::STT) << " : ";
      for (auto hit : t.fHits.at(TripletSolution::detID::STT)) {
        PndSttHit *myHit = static_cast<PndSttHit *>(hit);
        output << myHit->GetTubeID() << "/";
      }
      output << " ";
    }
    if (t.fHits.count(TripletSolution::detID::GEM) > 0) {
      output << t.detID_Names.at(TripletSolution::detID::GEM) << " : ";
      for (auto hit : t.fHits.at(TripletSolution::detID::GEM)) {
        output << hit->GetEntryNr() << "/";
      }
      output << " ";
    }

    output << "MSD: " << t.fMeanSquare;
    return output;
  }

  Triplet fTriplet{nullptr, nullptr, nullptr}; //< triplet used to create solution
  TVector3 fTrack;                             //< circle coordinates of track in x,y,r in [cm]
  std::map<detID, std::vector<FairHit *>> fHits;
  std::vector<FairHit *> fAllHits;
  double fMeanSquare = 1000000.;
  bool fClockwise = true;
};

struct TripletValues {
  std::vector<PndSttHit *> fFirstRow;
  std::vector<PndSttHit *> fMidRow;
  std::vector<PndSttHit *> fLastRow;

  TripletValues &operator+=(const TripletValues &right)
  {
    fFirstRow.insert(fFirstRow.end(), right.fFirstRow.begin(), right.fFirstRow.end());
    fMidRow.insert(fMidRow.end(), right.fMidRow.begin(), right.fMidRow.end());
    fLastRow.insert(fLastRow.end(), right.fLastRow.begin(), right.fLastRow.end());

    return *this;
  }

  friend std::ostream &operator<<(std::ostream &output, const TripletValues &t)
  {
    output << "FirstRow: ";
    std::for_each(t.fFirstRow.begin(), t.fFirstRow.end(), [](PndSttHit *hit) { std::cout << hit->GetTubeID() << "(" << hit->GetX() << "," << hit->GetY() << ") /"; });
    output << std::endl;
    output << "MidRow: ";
    std::for_each(t.fMidRow.begin(), t.fMidRow.end(), [](PndSttHit *hit) { std::cout << hit->GetTubeID() << "(" << hit->GetX() << "," << hit->GetY() << ") /"; });
    output << std::endl;
    output << "LastRow: ";
    std::for_each(t.fLastRow.begin(), t.fLastRow.end(), [](PndSttHit *hit) { std::cout << hit->GetTubeID() << "(" << hit->GetX() << "," << hit->GetY() << ") /"; });
    output << std::endl;

    return output;
  }
};

struct ReductionMaps {
  std::map<int, PndSttHit *> fMapCAToMidTube;
  std::map<int, PndSttHit *> fMapCAToOuterTube;
  std::map<PndSttHit *, int> fMapInnerTubeToCA;
  std::map<PndSttHit *, int> fMapMidTubeToCA;

  void Reset()
  {
    fMapCAToMidTube.clear();
    fMapCAToOuterTube.clear();
    fMapInnerTubeToCA.clear();
    fMapMidTubeToCA.clear();
  }
};

struct ApolloniusTripletFunctions {
  std::map<int, std::vector<PndSttHit *>> GetAllTubesByRow(std::vector<PndSttHit *> hits, PndSttGeometryMap *fGeometryMap, int &fAllHitsCounter);

  std::map<int, std::vector<std::vector<PndSttHit *>>>
  GetTubeStructure(std::vector<PndSttHit *> hits, PndSttGeometryMap *fGeometryMap, int &fAllHitsCounter, bool &fIsStrongCurling);

  ReductionMaps CreateMaps(TripletValues &triplets, std::map<FairLink, int> &fMapHitstoCATracklet);

  std::vector<Triplet> ReduceCombinatorics(TripletValues &triplets, std::map<FairLink, int> &fMapHitstoCATracklet, TClonesArray *sttHits, PndSttCA *fCATrackFinder,
                                           PndSttGeometryMap *fGeometryMap, double &fMinDistance, bool &fWithCombiReduction);

  bool IsTripletUsed(std::vector<TripletSolution> &solutions, Triplet &triplet);

  std::vector<TripletSolution> GenerateTripletTracks(Triplet triplet, std::vector<PndSttHit *> &sttHits, PndSttGeometryMap *fGeometryMap);

  std::vector<double> calcIntersectionPointCircleLine(std::vector<double> circle, double m, double b, double Ax);

  TripletSolution FindHitsCloseToCircle(std::vector<PndSttHit *> &sttHits, TVector3 &circle, double &fDistanceThresholdSTTFar, PndSttGeometryMap *fGeometryMap);

  double MeanSquareDistance(TripletSolution &solution);

  double DistanceCirclePoint(TVector3 &circle, FairHit *hit);

  double DistanceCircleSttHit(TVector3 &circle, PndSttHit *sttHit);

  double DistanceCircleSttHit(TVector3 &circle, PndSttHit *sttHit, double &sqaredDistance);

  double DistanceLineSttHit(TVector3 &circle, PndSttHit *sttHit, double &squaredDistance) { return sqrt(squaredDistance); };

  double DistanceLineSttHit(TVector3 &circle, PndSttHit *sttHit);

  double SquaredDistanceLineSttHit(TVector3 &circle, PndSttHit *sttHit);

  double DistanceLinePoint(TVector3 &circle, FairHit *hit);

  void CheckContinuitySolutions(std::vector<TripletSolution> &solutions, PndSttGeometryMap *fGeometryMap);

  bool IsContinuous(TripletSolution &solution, PndSttGeometryMap *fGeometryMap);

  void AddOtherDetectors(vector<TripletSolution> &solutions, PndSttStrawMap *fStrawMap, std::map<TString, TClonesArray *> &fBranchMap, TString sttname);

  std::vector<TripletSolution> FindBestSolutions(std::vector<TripletSolution> &solutions);

  std::vector<TripletSolution> CheckCombinedSolutions(std::vector<TripletSolution> &solutions, int nExpectedTracks);

  std::vector<TripletSolution> CheckSolutions(std::vector<TripletSolution> &solutions, std::vector<std::vector<int>> combinations);

  std::vector<int> GetUniqueTubeIDs(std::vector<TripletSolution> &solutions, std::vector<int> combinations);

  std::vector<std::vector<int>> GetKOutOfN(int k, int n);

  std::vector<TripletSolution> CombineIdenticalSolutions(std::vector<TripletSolution> &solutions);

  bool ContainsTriplet(Triplet &triplet, TripletSolution &solution);

  std::vector<PndSttHit *> TubeReduction(std::vector<std::pair<int, int>> &Rows, std::map<int, std::vector<std::vector<PndSttHit *>>> &tubeStructure,
                                         std::map<FairLink, int> &fMapHitstoCATracklet, int position);
  std::map<FairLink, int> GetHitsToCAMap(PndSttCA *fCATrackFinder);

  double CalcRatioSameHits(TripletSolution &sol, std::vector<TripletSolution> &AlreadyFoundSolutions);
  std::vector<TripletSolution> CombineIdenticalSolutionsFinal(std::vector<TripletSolution> &solutions, Double_t ratioOfSameHits = 0.6);

  PndTrack FromTripletSolutionToPndTrack(PndApollonius::TripletSolution &sol, Double_t B);

  PndTrack FromTripletSolutionToPndTrack(PndApollonius::TripletSolution &sol, Double_t B, PndTrackCand &cand);
};

} // namespace PndApollonius