SBNRecoUtils Namespace Reference

Functions
std::vector< geo::BoxBoundedGeo >	ActiveVolumes (const geo::GeometryCore *geometry)
std::vector< std::vector < geo::BoxBoundedGeo > >	TPCVolumes (const geo::GeometryCore *geometry)
geo::BoxBoundedGeo	DetectorVolume (const geo::GeometryCore *geometry)
double	MaxLength (const geo::GeometryCore *geometry)
int	TrueParticleID (const core::ProviderManager &manager, const art::Ptr< recob::Hit > hit, bool rollup_unsaved_ids=1)
int	TrueParticleIDFromTotalTrueEnergy (const core::ProviderManager &manager, const std::vector< art::Ptr< recob::Hit > > &hits, bool rollup_unsaved_ids=1)
int	TrueParticleIDFromTotalRecoCharge (const core::ProviderManager &manager, const std::vector< art::Ptr< recob::Hit > > &hits, bool rollup_unsaved_ids=1)
int	TrueParticleIDFromTotalRecoHits (const core::ProviderManager &manager, const std::vector< art::Ptr< recob::Hit > > &hits, bool rollup_unsaved_ids=1)
bool	IsInsideTPC (const core::ProviderManager &manager, TVector3 position, double distance_buffer)
double	CalculateTrackLength (const core::ProviderManager &manager, const art::Ptr< recob::Track > track)
double	TrackCompletion (const core::ProviderManager &manager, int mcparticle_id, const std::vector< art::Ptr< recob::Hit >> &reco_track_hits)

Function Documentation

std::vector< geo::BoxBoundedGeo > SBNRecoUtils::ActiveVolumes

(

const geo::GeometryCore *

geometry

)

Definition at line 3 of file GeoUtil.cc.

                                                                                       {
  std::vector<std::vector<geo::BoxBoundedGeo>> tpc_volumes = TPCVolumes(geometry);;
  std::vector<geo::BoxBoundedGeo> active_volumes;

  // make each cryostat volume a box enclosing all tpc volumes
  for (const std::vector<geo::BoxBoundedGeo> &tpcs: tpc_volumes) {
    double XMin = std::min_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MinX() < rhs.MinX(); })->MinX();
    double YMin = std::min_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MinY() < rhs.MinY(); })->MinY();
    double ZMin = std::min_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MinZ() < rhs.MinZ(); })->MinZ();
    
    double XMax = std::max_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MaxX() < rhs.MaxX(); })->MaxX();
    double YMax = std::max_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MaxY() < rhs.MaxY(); })->MaxY();
    double ZMax = std::max_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MaxZ() < rhs.MaxZ(); })->MaxZ();
    
    active_volumes.emplace_back(XMin, XMax, YMin, YMax, ZMin, ZMax);
  } 
  return active_volumes;
}

double SBNRecoUtils::CalculateTrackLength	(	const core::ProviderManager &	manager,
		const art::Ptr< recob::Track >	track
	)

Definition at line 163 of file sbnana/sbnanalysis/ana/SBNOscReco/RecoUtils/RecoUtils.cc.

                                                                                                           {
  double length = 0;
  if (track->NumberTrajectoryPoints()==1) return length; //Nothing to calculate if there is only one point

  for (size_t i_tp = 0; i_tp < track->NumberTrajectoryPoints()-1; i_tp++){ //Loop from the first to 2nd to last point
    TVector3 this_point(track->TrajectoryPoint(i_tp).position.X(),track->TrajectoryPoint(i_tp).position.Y(),track->TrajectoryPoint(i_tp).position.Z());
    if (!SBNRecoUtils::IsInsideTPC(manager, this_point,0)){
      std::cout<<"SBNRecoUtils::CalculateTrackLength - Current trajectory point not in the TPC volume.  Skip over this point in the track length calculation"<<std::endl;
      continue;
    }
    TVector3 next_point(track->TrajectoryPoint(i_tp+1).position.X(),track->TrajectoryPoint(i_tp+1).position.Y(),track->TrajectoryPoint(i_tp+1).position.Z());
    if (!SBNRecoUtils::IsInsideTPC(manager, next_point,0)){
      std::cout<<"SBNRecoUtils::CalculateTrackLength - Next trajectory point not in the TPC volume.  Skip over this point in the track length calculation"<<std::endl;
      continue;
    }

    length+=(next_point-this_point).Mag();
  }
  return length;
}

geo::BoxBoundedGeo SBNRecoUtils::DetectorVolume

(

const geo::GeometryCore *

geometry

)

Definition at line 35 of file GeoUtil.cc.

                                                                             {
  std::vector<geo::BoxBoundedGeo> active_volumes = ActiveVolumes(geometry);

  double XMin = std::min_element(active_volumes.begin(), active_volumes.end(), [](auto &lhs, auto &rhs) { return lhs.MinX() < rhs.MinX(); })->MinX();
  double YMin = std::min_element(active_volumes.begin(), active_volumes.end(), [](auto &lhs, auto &rhs) { return lhs.MinY() < rhs.MinY(); })->MinY();
  double ZMin = std::min_element(active_volumes.begin(), active_volumes.end(), [](auto &lhs, auto &rhs) { return lhs.MinZ() < rhs.MinZ(); })->MinZ();
  
  double XMax = std::max_element(active_volumes.begin(), active_volumes.end(), [](auto &lhs, auto &rhs) { return lhs.MaxX() < rhs.MaxX(); })->MaxX();
  double YMax = std::max_element(active_volumes.begin(), active_volumes.end(), [](auto &lhs, auto &rhs) { return lhs.MaxY() < rhs.MaxY(); })->MaxY();
  double ZMax = std::max_element(active_volumes.begin(), active_volumes.end(), [](auto &lhs, auto &rhs) { return lhs.MaxZ() < rhs.MaxZ(); })->MaxZ();

  return geo::BoxBoundedGeo(XMin, XMax, YMin, YMax, ZMin, ZMax);
}

bool SBNRecoUtils::IsInsideTPC	(	const core::ProviderManager &	manager,
		TVector3	position,
		double	distance_buffer
	)

Definition at line 113 of file sbnana/sbnanalysis/ana/SBNOscReco/RecoUtils/RecoUtils.cc.

                                                                                                           {
  double vtx[3] = {position.X(), position.Y(), position.Z()};
  bool inside = false;
  const geo::GeometryCore* geom = manager.GetGeometryProvider();
  geo::TPCID idtpc = geom->FindTPCAtPosition(vtx);

  if (geom->HasTPC(idtpc))
  {
    const geo::TPCGeo& tpcgeo = geom->GetElement(idtpc);
    double minx = tpcgeo.MinX(); double maxx = tpcgeo.MaxX();
    double miny = tpcgeo.MinY(); double maxy = tpcgeo.MaxY();
    double minz = tpcgeo.MinZ(); double maxz = tpcgeo.MaxZ();

    for (size_t c = 0; c < geom->Ncryostats(); c++)
    {
      const geo::CryostatGeo& cryostat = geom->Cryostat(c);
      for (size_t t = 0; t < cryostat.NTPC(); t++)
      {
        const geo::TPCGeo& tpcg = cryostat.TPC(t);
        if (tpcg.MinX() < minx) minx = tpcg.MinX();
        if (tpcg.MaxX() > maxx) maxx = tpcg.MaxX();
        if (tpcg.MinY() < miny) miny = tpcg.MinY();
        if (tpcg.MaxY() > maxy) maxy = tpcg.MaxY();
        if (tpcg.MinZ() < minz) minz = tpcg.MinZ();
        if (tpcg.MaxZ() > maxz) maxz = tpcg.MaxZ();
      }
    }

    //x
    double dista = fabs(minx - position.X());
    double distb = fabs(position.X() - maxx);
    if ((position.X() > minx) && (position.X() < maxx) &&
        (dista > distance_buffer) && (distb > distance_buffer)) inside = true;
    //y
    dista = fabs(maxy - position.Y());
    distb = fabs(position.Y() - miny);
    if (inside && (position.Y() > miny) && (position.Y() < maxy) &&
        (dista > distance_buffer) && (distb > distance_buffer)) inside = true;
    else inside = false;
    //z
    dista = fabs(maxz - position.Z());
    distb = fabs(position.Z() - minz);
    if (inside && (position.Z() > minz) && (position.Z() < maxz) &&
        (dista > distance_buffer) && (distb > distance_buffer)) inside = true;
    else inside = false;
  }

  return inside;
}

double SBNRecoUtils::MaxLength

(

const geo::GeometryCore *

geometry

)

Definition at line 22 of file GeoUtil.cc.

                                                              {
  std::vector<geo::BoxBoundedGeo> active_volumes = ActiveVolumes(geometry);
  double max_length = -1;
  for (const geo::BoxBoundedGeo &vol: active_volumes) {
    double X = vol.MaxX() - vol.MinX();
    double Y = vol.MaxY() - vol.MinY();
    double Z = vol.MaxZ() - vol.MinZ();
    double this_length = sqrt(X*X + Y*Y + Z*Z);
    if (this_length > max_length) max_length = this_length;
  }
  return max_length;
}

std::vector< std::vector< geo::BoxBoundedGeo > > SBNRecoUtils::TPCVolumes

(

const geo::GeometryCore *

geometry

)

Definition at line 49 of file GeoUtil.cc.

                                                                                               {
  std::vector<std::vector<geo::BoxBoundedGeo>> tpc_volumes;
  for (auto const &cryo: geometry->IterateCryostats()) {
    geo::GeometryCore::TPC_iterator iTPC = geometry->begin_TPC(cryo.ID()),
                                    tend = geometry->end_TPC(cryo.ID());
    std::vector<geo::BoxBoundedGeo> this_tpc_volumes;
    while (iTPC != tend) {
      geo::TPCGeo const& TPC = *iTPC;
      this_tpc_volumes.push_back(TPC.ActiveBoundingBox());
      iTPC++;
    }
     tpc_volumes.push_back(std::move(this_tpc_volumes));
  }
  return tpc_volumes;
}

double SBNRecoUtils::TrackCompletion	(	const core::ProviderManager &	manager,
		int	mcparticle_id,
		const std::vector< art::Ptr< recob::Hit >> &	reco_track_hits
	)

Definition at line 184 of file sbnana/sbnanalysis/ana/SBNOscReco/RecoUtils/RecoUtils.cc.

                                                                                                                                              {
  // get handle to back tracker
  cheat::BackTracker *bt = manager.GetBackTrackerProvider();

  // get all the IDE's of the truth track
  const std::vector<const sim::IDE*> mcparticle_ides = bt->TrackIdToSimIDEs_Ps(mcparticle_id);
  // sum it up
  double mcparticle_energy = 0.;
  for (auto const &ide: mcparticle_ides) {
    mcparticle_energy += ide->energy;
  }

  // get all the hits of the reco track that match the truth track
  const std::vector<art::Ptr<recob::Hit>> matched_reco_track_hits = bt->TrackIdToHits_Ps(mcparticle_id, reco_track_hits);

  // for each of the hits get the energy coming from the track
  double matched_reco_energy = 0.;
  for (auto const &matched_reco_track_hit: matched_reco_track_hits) {
    std::vector<sim::IDE> this_hit_IDEs = bt->HitToAvgSimIDEs(*matched_reco_track_hit);
    for (auto const &ide: this_hit_IDEs) {
      if (ide.trackID == mcparticle_id) {
        matched_reco_energy += ide.energy;
      }
    }
  }

  return matched_reco_energy / mcparticle_energy;
}

int SBNRecoUtils::TrueParticleID	(	const core::ProviderManager &	manager,
		const art::Ptr< recob::Hit >	hit,
		bool	rollup_unsaved_ids = 1
	)

Definition at line 3 of file sbnana/sbnanalysis/ana/SBNOscReco/RecoUtils/RecoUtils.cc.

                                                                                                                        {
  std::map<int,double> id_to_energy_map;
  const cheat::BackTracker *bt_serv = manager.GetBackTrackerProvider();
  std::vector<sim::TrackIDE> track_ides = bt_serv->HitToTrackIDEs(hit);
  for (unsigned int idIt = 0; idIt < track_ides.size(); ++idIt) {
    int id = track_ides.at(idIt).trackID;
    if (rollup_unsaved_ids) id = std::abs(id);
    double energy = track_ides.at(idIt).energy;
    id_to_energy_map[id]+=energy;
  }
  //Now loop over the map to find the maximum contributor
  double likely_particle_contrib_energy = -99999;
  int likely_track_id = 0;
  for (std::map<int,double>::iterator mapIt = id_to_energy_map.begin(); mapIt != id_to_energy_map.end(); mapIt++){
    double particle_contrib_energy = mapIt->second;
    if (particle_contrib_energy > likely_particle_contrib_energy){
      likely_particle_contrib_energy = particle_contrib_energy;
      likely_track_id = mapIt->first;
    }
  }
  return likely_track_id;
}

int SBNRecoUtils::TrueParticleIDFromTotalRecoCharge	(	const core::ProviderManager &	manager,
		const std::vector< art::Ptr< recob::Hit > > &	hits,
		bool	rollup_unsaved_ids = 1
	)

Definition at line 58 of file sbnana/sbnanalysis/ana/SBNOscReco/RecoUtils/RecoUtils.cc.

                                                                                                                                                         {
  // Make a map of the tracks which are associated with this object and the charge each contributes
  std::map<int,double> trackMap;
  for (std::vector<art::Ptr<recob::Hit> >::const_iterator hitIt = hits.begin(); hitIt != hits.end(); ++hitIt) {
    art::Ptr<recob::Hit> hit = *hitIt;
    int trackID = TrueParticleID(manager, hit, rollup_unsaved_ids);
    trackMap[trackID] += hit->Integral();
  }

  // Pick the track with the highest charge as the 'true track'
  double highestCharge = 0;
  int objectTrack = -99999;
  for (std::map<int,double>::iterator trackIt = trackMap.begin(); trackIt != trackMap.end(); ++trackIt) {
    if (trackIt->second > highestCharge) {
      highestCharge = trackIt->second;
      objectTrack  = trackIt->first;
    }
  }
  return objectTrack;
}

int SBNRecoUtils::TrueParticleIDFromTotalRecoHits	(	const core::ProviderManager &	manager,
		const std::vector< art::Ptr< recob::Hit > > &	hits,
		bool	rollup_unsaved_ids = 1
	)

Definition at line 81 of file sbnana/sbnanalysis/ana/SBNOscReco/RecoUtils/RecoUtils.cc.

                                                                                                                                                       {
  // Make a map of the tracks which are associated with this object and the number of hits they are the primary contributor to
  std::map<int,int> trackMap;
  for (std::vector<art::Ptr<recob::Hit> >::const_iterator hitIt = hits.begin(); hitIt != hits.end(); ++hitIt) {
    art::Ptr<recob::Hit> hit = *hitIt;
    int trackID = TrueParticleID(manager, hit, rollup_unsaved_ids);
    trackMap[trackID]++;
  }

  // Pick the track which is the primary contributor to the most hits as the 'true track'
  int objectTrack = -99999;
  int highestCount = -1;
  int NHighestCounts = 0;
  for (std::map<int,int>::iterator trackIt = trackMap.begin(); trackIt != trackMap.end(); ++trackIt) {
    if (trackIt->second > highestCount) {
      highestCount = trackIt->second;
      objectTrack  = trackIt->first;
      NHighestCounts = 1;
    }
    else if (trackIt->second == highestCount){
      NHighestCounts++;
    }
  }
  if (NHighestCounts > 1){
    std::cout<<"SBNRecoUtils::TrueParticleIDFromTotalRecoHits - There are " << NHighestCounts << " particles which tie for highest number of contributing hits (" << highestCount<<" hits).  Using SBNRecoUtils::TrueParticleIDFromTotalTrueEnergy instead."<<std::endl;
    objectTrack = SBNRecoUtils::TrueParticleIDFromTotalTrueEnergy(manager, hits,rollup_unsaved_ids);
  }
  return objectTrack;
}

int SBNRecoUtils::TrueParticleIDFromTotalTrueEnergy	(	const core::ProviderManager &	manager,
		const std::vector< art::Ptr< recob::Hit > > &	hits,
		bool	rollup_unsaved_ids = 1
	)

Definition at line 28 of file sbnana/sbnanalysis/ana/SBNOscReco/RecoUtils/RecoUtils.cc.

                                                                                                                                                         {
  const cheat::BackTracker *bt_serv = manager.GetBackTrackerProvider();
  std::map<int,double> trackIDToEDepMap;
  for (std::vector<art::Ptr<recob::Hit> >::const_iterator hitIt = hits.begin(); hitIt != hits.end(); ++hitIt) {
    art::Ptr<recob::Hit> hit = *hitIt;
    std::vector<sim::TrackIDE> trackIDs = bt_serv->HitToTrackIDEs(hit);
    for (unsigned int idIt = 0; idIt < trackIDs.size(); ++idIt) {
      int id = trackIDs[idIt].trackID;
      if (rollup_unsaved_ids) id = std::abs(id);
      trackIDToEDepMap[id] += trackIDs[idIt].energy;
    }
  }

  //Loop over the map and find the track which contributes the highest energy to the hit vector
  double maxenergy = -1;
  int objectTrack = -99999;
  for (std::map<int,double>::iterator mapIt = trackIDToEDepMap.begin(); mapIt != trackIDToEDepMap.end(); mapIt++){
    double energy = mapIt->second;
    double trackid = mapIt->first;
    if (energy > maxenergy){
      maxenergy = energy;
      objectTrack = trackid;
    }
  }

  return objectTrack;
}