trkf::CosmicTracker Class Reference

Inheritance diagram for trkf::CosmicTracker:

Public Member Functions
	CosmicTracker (fhicl::ParameterSet const &pset)

Private Member Functions
void	produce (art::Event &evt)

Private Attributes
cluster::ClusterMatchTQ	fClusterMatch
trkf::CosmicTrackerAlg	fCTAlg
std::string	fClusterModuleLabel
	label for input cluster collection More...
std::string	fSortDir
	sort space points More...
bool	fStitchTracks
	Stitch tracks from different TPCs. More...
double	fDisCut
	Distance cut for track merging. More...
double	fAngCut
	Angle cut for track merging. More...
bool	fTrajOnly
	Only use trajectory points from TrackTrajectoryAlg for debugging. More...

Detailed Description

Definition at line 216 of file CosmicTracker_module.cc.

Constructor & Destructor Documentation

trkf::CosmicTracker::CosmicTracker ( fhicl::ParameterSet const &  pset )

explicit

Definition at line 243 of file CosmicTracker_module.cc.

    : EDProducer{pset}
    , fClusterMatch(pset.get<fhicl::ParameterSet>("ClusterMatch"))
    , fCTAlg(pset.get<fhicl::ParameterSet>("CTAlg"))
  {
    fClusterModuleLabel = pset.get<std::string>("ClusterModuleLabel");
    fSortDir = pset.get<std::string>("SortDirection", "+z");
    fStitchTracks = pset.get<bool>("StitchTracks");
    fDisCut = pset.get<double>("DisCut");
    fAngCut = pset.get<double>("AngCut");
    fTrajOnly = pset.get<bool>("TrajOnly");

    produces<std::vector<recob::Track>>();
    produces<std::vector<recob::SpacePoint>>();
    produces<art::Assns<recob::Track, recob::Cluster>>();
    produces<art::Assns<recob::Track, recob::SpacePoint>>();
    produces<art::Assns<recob::SpacePoint, recob::Hit>>();
    produces<art::Assns<recob::Track, recob::Hit>>();
  }

Member Function Documentation

void trkf::CosmicTracker::produce ( art::Event &  evt )

private

Definition at line 265 of file CosmicTracker_module.cc.

  {

    // get services
    art::ServiceHandle<geo::Geometry const> geom;

    std::unique_ptr<std::vector<recob::Track>> tcol(new std::vector<recob::Track>);
    std::unique_ptr<std::vector<recob::SpacePoint>> spcol(new std::vector<recob::SpacePoint>);
    std::unique_ptr<art::Assns<recob::Track, recob::SpacePoint>> tspassn(
      new art::Assns<recob::Track, recob::SpacePoint>);
    std::unique_ptr<art::Assns<recob::Track, recob::Cluster>> tcassn(
      new art::Assns<recob::Track, recob::Cluster>);
    std::unique_ptr<art::Assns<recob::Track, recob::Hit>> thassn(
      new art::Assns<recob::Track, recob::Hit>);
    std::unique_ptr<art::Assns<recob::SpacePoint, recob::Hit>> shassn(
      new art::Assns<recob::SpacePoint, recob::Hit>);

    //double timetick = detprop->SamplingRate()*1e-3;    //time sample in us
    //double presamplings = detprop->TriggerOffset(); // presamplings in ticks
    //double plane_pitch = geom->PlanePitch(0,1);   //wire plane pitch in cm
    //double wire_pitch = geom->WirePitch();    //wire pitch in cm
    //double Efield_drift = larprop->Efield(0);  // Electric Field in the drift region in kV/cm
    //double Temperature = detprop->Temperature();  // LAr Temperature in K

    //double driftvelocity = larprop->DriftVelocity(Efield_drift,Temperature);    //drift velocity in the drift region (cm/us)
    //double timepitch = driftvelocity*timetick;                         //time sample (cm)

    // get input Cluster object(s).
    art::Handle<std::vector<recob::Cluster>> clusterListHandle;
    std::vector<art::Ptr<recob::Cluster>> clusterlist;
    if (evt.getByLabel(fClusterModuleLabel, clusterListHandle))
      art::fill_ptr_vector(clusterlist, clusterListHandle);

    art::FindManyP<recob::Hit> fm(clusterListHandle, evt, fClusterModuleLabel);

    // find matched clusters
    auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
    auto const detProp =
      art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt, clockData);
    auto const matchedclusters = fClusterMatch.MatchedClusters(clockData, detProp, clusterlist, fm);

    // get track space points
    std::vector<std::vector<trkPoint>> trkpts(matchedclusters.size());
    for (size_t itrk = 0; itrk < matchedclusters.size(); ++itrk) { //loop over tracks

      std::vector<art::Ptr<recob::Hit>> hitlist;
      for (size_t iclu = 0; iclu < matchedclusters[itrk].size(); ++iclu) { //loop over clusters

        std::vector<art::Ptr<recob::Hit>> hits = fm.at(matchedclusters[itrk][iclu]);
        for (size_t ihit = 0; ihit < hits.size(); ++ihit) {
          hitlist.push_back(hits[ihit]);
        }
      }
      // reconstruct space points and directions
      fCTAlg.SPTReco(clockData, detProp, hitlist);
      if (!fTrajOnly) {
        if (!fCTAlg.trkPos.size()) continue;
        for (size_t i = 0; i < hitlist.size(); ++i) {
          trkPoint trkpt;
          trkpt.pos = fCTAlg.trkPos[i];
          trkpt.dir = fCTAlg.trkDir[i];
          trkpt.hit = hitlist[i];
          trkpts[itrk].push_back(trkpt);
        }
        if (fSortDir == "+x") std::sort(trkpts[itrk].begin(), trkpts[itrk].end(), sp_sort_x0);
        if (fSortDir == "-x") std::sort(trkpts[itrk].begin(), trkpts[itrk].end(), sp_sort_x1);
        if (fSortDir == "+y") std::sort(trkpts[itrk].begin(), trkpts[itrk].end(), sp_sort_y0);
        if (fSortDir == "-y") std::sort(trkpts[itrk].begin(), trkpts[itrk].end(), sp_sort_y1);
        if (fSortDir == "+z") std::sort(trkpts[itrk].begin(), trkpts[itrk].end(), sp_sort_z0);
        if (fSortDir == "-z") std::sort(trkpts[itrk].begin(), trkpts[itrk].end(), sp_sort_z1);
      }

      if (fTrajOnly) { //debug only
        if (!fCTAlg.trajPos.size()) continue;
        size_t spStart = spcol->size();
        std::vector<recob::SpacePoint> spacepoints;
        for (size_t ipt = 0; ipt < fCTAlg.trajPos.size(); ++ipt) {
          art::PtrVector<recob::Hit> sp_hits;
          double hitcoord[3];
          hitcoord[0] = fCTAlg.trajPos[ipt].X();
          hitcoord[1] = fCTAlg.trajPos[ipt].Y();
          hitcoord[2] = fCTAlg.trajPos[ipt].Z();
          double err[6] = {util::kBogusD};
          recob::SpacePoint mysp(hitcoord,
                                 err,
                                 util::kBogusD,
                                 spStart + spacepoints.size()); //3d point at end of track
          spacepoints.push_back(mysp);
          spcol->push_back(mysp);
          util::CreateAssn(evt, *spcol, fCTAlg.trajHit[ipt], *shassn);
        } // ihit
        size_t spEnd = spcol->size();
        if (fSortDir == "+x") {
          std::sort(spacepoints.begin(), spacepoints.end(), spt_sort_x0);
          std::sort(spcol->begin() + spStart, spcol->begin() + spEnd, spt_sort_x0);
        }
        if (fSortDir == "-x") {
          std::sort(spacepoints.begin(), spacepoints.end(), spt_sort_x1);
          std::sort(spcol->begin() + spStart, spcol->begin() + spEnd, spt_sort_x1);
        }
        if (fSortDir == "+y") {
          std::sort(spacepoints.begin(), spacepoints.end(), spt_sort_y0);
          std::sort(spcol->begin() + spStart, spcol->begin() + spEnd, spt_sort_y0);
        }
        if (fSortDir == "-y") {
          std::sort(spacepoints.begin(), spacepoints.end(), spt_sort_y1);
          std::sort(spcol->begin() + spStart, spcol->begin() + spEnd, spt_sort_y1);
        }
        if (fSortDir == "+z") {
          std::sort(spacepoints.begin(), spacepoints.end(), spt_sort_z0);
          std::sort(spcol->begin() + spStart, spcol->begin() + spEnd, spt_sort_z0);
        }
        if (fSortDir == "-z") {
          std::sort(spacepoints.begin(), spacepoints.end(), spt_sort_z1);
          std::sort(spcol->begin() + spStart, spcol->begin() + spEnd, spt_sort_z1);
        }

        if (spacepoints.size() > 0) {

          // make a vector of the trajectory points along the track
          std::vector<TVector3> xyz(spacepoints.size());
          for (size_t s = 0; s < spacepoints.size(); ++s) {
            xyz[s] = TVector3(spacepoints[s].XYZ());
          }
          //Calculate track direction cosines
          TVector3 startpointVec, endpointVec, DirCos;
          startpointVec = xyz[0];
          endpointVec = xyz.back();
          DirCos = endpointVec - startpointVec;
          //SetMag casues a crash if the magnitude of the vector is zero
          try {
            DirCos.SetMag(1.0); //normalize vector
          }
          catch (...) {
            std::cout << "The Spacepoint is infinitely small" << std::endl;
            continue;
          }
          std::vector<TVector3> dircos(spacepoints.size(), DirCos);

          tcol->push_back(
            recob::Track(recob::TrackTrajectory(recob::tracking::convertCollToPoint(xyz),
                                                recob::tracking::convertCollToVector(dircos),
                                                recob::Track::Flags_t(xyz.size()),
                                                false),
                         0,
                         -1.,
                         0,
                         recob::tracking::SMatrixSym55(),
                         recob::tracking::SMatrixSym55(),
                         tcol->size()));

          // make associations between the track and space points
          util::CreateAssn(evt, *tcol, *spcol, *tspassn, spStart, spEnd);

          // and the hits and track
          std::vector<art::Ptr<recob::Hit>> trkhits;
          for (size_t ihit = 0; ihit < hitlist.size(); ++ihit) {
            trkhits.push_back(hitlist[ihit]);
          }
          util::CreateAssn(evt, *tcol, trkhits, *thassn);
        }
      }
    } //itrk

    // by default creates a spacepoint and direction for each hit
    if (!fTrajOnly) {
      std::vector<std::vector<unsigned int>> trkidx;
      //stitch tracks from different TPCs
      if (fStitchTracks) {
        for (size_t itrk1 = 0; itrk1 < trkpts.size(); ++itrk1) {
          for (size_t itrk2 = itrk1 + 1; itrk2 < trkpts.size(); ++itrk2) {
            if (MatchTrack(trkpts[itrk1], trkpts[itrk2], fDisCut, fAngCut)) {
              int found1 = -1;
              int found2 = -1;
              for (size_t i = 0; i < trkidx.size(); ++i) {
                for (size_t j = 0; j < trkidx[i].size(); ++j) {
                  if (trkidx[i][j] == itrk1) found1 = i;
                  if (trkidx[i][j] == itrk2) found2 = i;
                }
              }
              if (found1 == -1 && found2 == -1) {
                std::vector<unsigned int> tmp;
                tmp.push_back(itrk1);
                tmp.push_back(itrk2);
                trkidx.push_back(tmp);
              }
              else if (found1 == -1 && found2 != -1) {
                trkidx[found2].push_back(itrk1);
              }
              else if (found1 != -1 && found2 == -1) {
                trkidx[found1].push_back(itrk2);
              }
              else if (found1 != found2) { //merge two vectors
                trkidx[found1].insert(
                  trkidx[found1].end(), trkidx[found2].begin(), trkidx[found2].end());
                trkidx.erase(trkidx.begin() + found2);
              }
            } //found match
          }   //itrk2
        }     //itrk1
        for (size_t itrk = 0; itrk < trkpts.size(); ++itrk) {
          bool found = false;
          for (size_t i = 0; i < trkidx.size(); ++i) {
            for (size_t j = 0; j < trkidx[i].size(); ++j) {
              if (trkidx[i][j] == itrk) found = true;
            }
          }
          if (!found) {
            std::vector<unsigned int> tmp;
            tmp.push_back(itrk);
            trkidx.push_back(tmp);
          }
        }
      } //stitch
      else {
        trkidx.resize(trkpts.size());
        for (size_t i = 0; i < trkpts.size(); ++i) {
          trkidx[i].push_back(i);
        }
      }

      //make recob::track and associations
      for (size_t i = 0; i < trkidx.size(); ++i) {
        //all track points
        std::vector<trkPoint> finaltrkpts;
        //all the clusters associated with the current track
        std::vector<art::Ptr<recob::Cluster>> clustersPerTrack;
        //all hits
        std::vector<art::Ptr<recob::Hit>> hitlist;
        for (size_t j = 0; j < trkidx[i].size(); ++j) {
          for (size_t k = 0; k < trkpts[trkidx[i][j]].size(); ++k) {
            finaltrkpts.push_back(trkpts[trkidx[i][j]][k]);
            hitlist.push_back(trkpts[trkidx[i][j]][k].hit);
          } //k
          for (size_t iclu = 0; iclu < matchedclusters[trkidx[i][j]].size(); ++iclu) {
            art::Ptr<recob::Cluster> cluster(clusterListHandle,
                                             matchedclusters[trkidx[i][j]][iclu]);
            clustersPerTrack.push_back(cluster);
          }

        } //j
        if (fStitchTracks) {
          if (fSortDir == "+x") std::sort(finaltrkpts.begin(), finaltrkpts.end(), sp_sort_x0);
          if (fSortDir == "-x") std::sort(finaltrkpts.begin(), finaltrkpts.end(), sp_sort_x1);
          if (fSortDir == "+y") std::sort(finaltrkpts.begin(), finaltrkpts.end(), sp_sort_y0);
          if (fSortDir == "-y") std::sort(finaltrkpts.begin(), finaltrkpts.end(), sp_sort_y1);
          if (fSortDir == "+z") std::sort(finaltrkpts.begin(), finaltrkpts.end(), sp_sort_z0);
          if (fSortDir == "-z") std::sort(finaltrkpts.begin(), finaltrkpts.end(), sp_sort_z1);
        }
        size_t spStart = spcol->size();
        std::vector<recob::SpacePoint> spacepoints;
        for (size_t ipt = 0; ipt < finaltrkpts.size(); ++ipt) {
          art::PtrVector<recob::Hit> sp_hits;
          sp_hits.push_back(finaltrkpts[ipt].hit);
          double hitcoord[3];
          hitcoord[0] = finaltrkpts[ipt].pos.X();
          hitcoord[1] = finaltrkpts[ipt].pos.Y();
          hitcoord[2] = finaltrkpts[ipt].pos.Z();
          double err[6] = {util::kBogusD};
          recob::SpacePoint mysp(hitcoord,
                                 err,
                                 util::kBogusD,
                                 spStart + spacepoints.size()); //3d point at end of track
          spacepoints.push_back(mysp);
          spcol->push_back(mysp);
          util::CreateAssn(evt, *spcol, sp_hits, *shassn);
        } // ipt
        size_t spEnd = spcol->size();
        if (spacepoints.size() > 0) {
          // make a vector of the trajectory points along the track
          std::vector<TVector3> xyz(spacepoints.size());
          std::vector<TVector3> dircos(spacepoints.size());
          for (size_t s = 0; s < spacepoints.size(); ++s) {
            xyz[s] = TVector3(spacepoints[s].XYZ());
            dircos[s] = finaltrkpts[s].dir;
            //flip direction if needed.
            if (spacepoints.size() > 1) {
              if (s == 0) {
                TVector3 xyz1 = TVector3(spacepoints[s + 1].XYZ());
                TVector3 dir = xyz1 - xyz[s];
                if (dir.Angle(dircos[s]) > 0.8 * TMath::Pi()) { dircos[s] = -dircos[s]; }
              }
              else {
                TVector3 dir = xyz[s] - xyz[s - 1];
                if (dir.Angle(dircos[s]) > 0.8 * TMath::Pi()) { dircos[s] = -dircos[s]; }
              }
            }
          }
          tcol->push_back(
            recob::Track(recob::TrackTrajectory(recob::tracking::convertCollToPoint(xyz),
                                                recob::tracking::convertCollToVector(dircos),
                                                recob::Track::Flags_t(xyz.size()),
                                                false),
                         0,
                         -1.,
                         0,
                         recob::tracking::SMatrixSym55(),
                         recob::tracking::SMatrixSym55(),
                         tcol->size()));

          // make associations between the track and space points
          util::CreateAssn(evt, *tcol, *spcol, *tspassn, spStart, spEnd);

          // now the track and clusters
          util::CreateAssn(evt, *tcol, clustersPerTrack, *tcassn);

          // and the hits and track
          std::vector<art::Ptr<recob::Hit>> trkhits;
          for (size_t ihit = 0; ihit < hitlist.size(); ++ihit) {
            trkhits.push_back(hitlist[ihit]);
          }
          util::CreateAssn(evt, *tcol, trkhits, *thassn);
        }

      } //i
    }
    mf::LogVerbatim("Summary") << std::setfill('-') << std::setw(175) << "-" << std::setfill(' ');
    mf::LogVerbatim("Summary") << "CosmicTracker Summary:";
    for (unsigned int i = 0; i < tcol->size(); ++i)
      mf::LogVerbatim("Summary") << tcol->at(i);
    mf::LogVerbatim("Summary") << "CosmicTracker Summary End:";

    evt.put(std::move(tcol));
    evt.put(std::move(spcol));
    evt.put(std::move(tspassn));
    evt.put(std::move(tcassn));
    evt.put(std::move(thassn));
    evt.put(std::move(shassn));

    return;
  }

Member Data Documentation

double trkf::CosmicTracker::fAngCut

private

Angle cut for track merging.

Definition at line 232 of file CosmicTracker_module.cc.

cluster::ClusterMatchTQ trkf::CosmicTracker::fClusterMatch

private

Definition at line 223 of file CosmicTracker_module.cc.

std::string trkf::CosmicTracker::fClusterModuleLabel

private

label for input cluster collection

Definition at line 226 of file CosmicTracker_module.cc.

trkf::CosmicTrackerAlg trkf::CosmicTracker::fCTAlg

private

Definition at line 224 of file CosmicTracker_module.cc.

double trkf::CosmicTracker::fDisCut

private

Distance cut for track merging.

Definition at line 231 of file CosmicTracker_module.cc.

std::string trkf::CosmicTracker::fSortDir

private

sort space points

Definition at line 228 of file CosmicTracker_module.cc.

bool trkf::CosmicTracker::fStitchTracks

private

Stitch tracks from different TPCs.

Definition at line 230 of file CosmicTracker_module.cc.

bool trkf::CosmicTracker::fTrajOnly

private

Only use trajectory points from TrackTrajectoryAlg for debugging.

Definition at line 234 of file CosmicTracker_module.cc.

---

The documentation for this class was generated from the following file:

• CosmicTracker_module.cc