TrackCalorimetryAlg.cxx

Go to the documentation of this file.

#include "TrackCalorimetryAlg.h"
#include "fhiclcpp/ParameterSet.h"

#include "larcorealg/Geometry/WireReadoutGeom.h"
#include "larcoreobj/SimpleTypesAndConstants/geo_types.h"
#include "lardata/ArtDataHelper/TrackUtils.h" // lar::util::TrackPitchInView()
#include "lardataalg/DetectorInfo/DetectorPropertiesData.h"
#include "lardataobj/AnalysisBase/Calorimetry.h"
#include "lardataobj/RecoBase/Hit.h"
#include "lardataobj/RecoBase/Track.h"
#include "lardataobj/RecoBase/TrackingTypes.h"

calo::TrackCalorimetryAlg::TrackCalorimetryAlg(fhicl::ParameterSet const& p)
  : caloAlg(p.get<fhicl::ParameterSet>("CalorimetryAlg"))
{
  fNHitsToDetermineStart = p.get<unsigned int>("NHitsToDetermineStart", 3);
}

void calo::TrackCalorimetryAlg::ExtractCalorimetry(
  detinfo::DetectorClocksData const& clock_data,
  detinfo::DetectorPropertiesData const& det_prop,
  std::vector<recob::Track> const& trackVector,
  std::vector<recob::Hit> const& hitVector,
  std::vector<std::vector<size_t>> const& hit_indices_per_track,
  std::vector<anab::Calorimetry>& caloVector,
  std::vector<size_t>& assnTrackCaloVector,
  geo::WireReadoutGeom const& wireReadoutGeom)
{
  //loop over the track list
  for (size_t i_track = 0; i_track < trackVector.size(); i_track++) {

    recob::Track const& track = trackVector[i_track];
    std::vector<float> path_length_fraction_vec(CreatePathLengthFractionVector(track));

    //sort hits into each plane
    std::vector<std::vector<size_t>> hit_indices_per_plane(wireReadoutGeom.Nplanes());
    for (auto const& i_hit : hit_indices_per_track[i_track])
      hit_indices_per_plane[hitVector[i_hit].WireID().Plane].push_back(i_hit);

    //loop over the planes
    for (unsigned int i_plane = 0; i_plane < wireReadoutGeom.Nplanes(); ++i_plane) {
      geo::PlaneID const planeID{0, 0, i_plane};

      //project down the track into wire/tick space for this plane
      std::vector<std::pair<geo::WireID, float>> traj_points_in_plane(
        track.NumberTrajectoryPoints());
      for (size_t i_trjpt = 0; i_trjpt < track.NumberTrajectoryPoints(); i_trjpt++) {
        double x_pos = track.LocationAtPoint(i_trjpt).X();
        float tick = det_prop.ConvertXToTicks(x_pos, (int)i_plane, 0, 0);
        traj_points_in_plane[i_trjpt] = std::make_pair(
          wireReadoutGeom.Plane(planeID).NearestWireID(track.LocationAtPoint(i_trjpt)), tick);
      }

      HitPropertiesMultiset_t HitPropertiesMultiset;
      //now loop through hits
      for (auto const& i_hit : hit_indices_per_plane[i_plane])
        AnalyzeHit(clock_data,
                   det_prop,
                   hitVector[i_hit],
                   track,
                   traj_points_in_plane,
                   path_length_fraction_vec,
                   HitPropertiesMultiset,
                   wireReadoutGeom);

      MakeCalorimetryObject(
        HitPropertiesMultiset, track, i_track, caloVector, assnTrackCaloVector, planeID);

    } //end loop over planes

  } //end loop over tracks

} //end ExtractCalorimetry

class dist_projected {
public:
  dist_projected(recob::Hit const& h, geo::WireReadoutGeom const& wireReadoutGeom)
    : hit(h), wireReadoutGeom(wireReadoutGeom)
  {}
  bool operator()(std::pair<geo::WireID, float> const& i, std::pair<geo::WireID, float> const& j)
  {
    float dw_i = ((int)(i.first.Wire) - (int)(hit.WireID().Wire)) *
                 wireReadoutGeom.Plane(i.first.asPlaneID()).WirePitch();
    float dw_j = ((int)(j.first.Wire) - (int)(hit.WireID().Wire)) *
                 wireReadoutGeom.Plane(j.first.asPlaneID()).WirePitch();
    float dt_i = i.second - hit.PeakTime();
    float dt_j = j.second - hit.PeakTime();
    return std::hypot(dw_i, dt_i) < std::hypot(dw_j, dt_j);
  }

private:
  recob::Hit const& hit;
  geo::WireReadoutGeom const& wireReadoutGeom;
};

std::vector<float> calo::TrackCalorimetryAlg::CreatePathLengthFractionVector(
  recob::Track const& track)
{
  std::vector<float> trk_path_length_frac_vec(track.NumberTrajectoryPoints());

  float cumulative_path_length = 0;
  const float total_path_length = track.Length();
  for (size_t i_trj = 1; i_trj < track.NumberTrajectoryPoints(); i_trj++) {
    cumulative_path_length += (track.LocationAtPoint(i_trj) - track.LocationAtPoint(i_trj - 1)).R();
    trk_path_length_frac_vec[i_trj] = cumulative_path_length / total_path_length;
  }

  return trk_path_length_frac_vec;
}

void calo::TrackCalorimetryAlg::AnalyzeHit(
  detinfo::DetectorClocksData const& clock_data,
  detinfo::DetectorPropertiesData const& det_prop,
  recob::Hit const& hit,
  recob::Track const& track,
  std::vector<std::pair<geo::WireID, float>> const& traj_points_in_plane,
  std::vector<float> const& path_length_fraction_vec,
  HitPropertiesMultiset_t& HitPropertiesMultiset,
  geo::WireReadoutGeom const& wireReadoutGeom)
{
  size_t traj_iter = std::distance(traj_points_in_plane.begin(),
                                   std::min_element(traj_points_in_plane.begin(),
                                                    traj_points_in_plane.end(),
                                                    dist_projected(hit, wireReadoutGeom)));
  float pitch = lar::util::TrackPitchInView(
    track, wireReadoutGeom.Plane(hit.WireID().asPlaneID()).View(), traj_iter);

  HitPropertiesMultiset.emplace(hit.Integral(),
                                hit.Integral() / pitch,
                                caloAlg.dEdx_AREA(clock_data, det_prop, hit, pitch),
                                pitch,
                                track.LocationAtPoint<TVector3>(traj_iter),
                                path_length_fraction_vec[traj_iter]);
}

bool calo::TrackCalorimetryAlg::IsInvertedTrack(HitPropertiesMultiset_t const& hpm)
{

  if (hpm.size() <= fNHitsToDetermineStart) return false;

  float charge_start = 0, charge_end = 0;
  unsigned int counter = 0;
  for (HitPropertiesMultiset_t::iterator it_hpm = hpm.begin(); it_hpm != hpm.end(); it_hpm++) {
    charge_start += (*it_hpm).charge;
    counter++;
    if (counter == fNHitsToDetermineStart) break;
  }

  counter = 0;
  for (HitPropertiesMultiset_t::reverse_iterator it_hpm = hpm.rbegin(); it_hpm != hpm.rend();
       it_hpm++) {
    charge_end += (*it_hpm).charge;
    counter++;
    if (counter == fNHitsToDetermineStart) break;
  }

  return charge_start > charge_end;
}

void calo::TrackCalorimetryAlg::MakeCalorimetryObject(HitPropertiesMultiset_t const& hpm,
                                                      recob::Track const& track,
                                                      size_t const& i_track,
                                                      std::vector<anab::Calorimetry>& caloVector,
                                                      std::vector<size_t>& assnTrackCaloVector,
                                                      geo::PlaneID const& planeID)
{
  size_t n_hits = hpm.size();
  std::vector<float> dEdxVector, dQdxVector, resRangeVector, deadWireVector, pitchVector;
  std::vector<TVector3> XYZVector;

  dEdxVector.reserve(n_hits);
  dQdxVector.reserve(n_hits);
  resRangeVector.reserve(n_hits);
  deadWireVector.reserve(n_hits);
  pitchVector.reserve(n_hits);
  XYZVector.reserve(n_hits);

  float kinetic_energy = 0, track_length = track.Length();
  if (IsInvertedTrack(hpm)) {
    for (auto it_hpm = hpm.begin(); it_hpm != hpm.end(); ++it_hpm) {
      dEdxVector.push_back(it_hpm->dEdx);
      dQdxVector.push_back(it_hpm->dQdx);
      resRangeVector.push_back(it_hpm->path_fraction * track_length);
      pitchVector.push_back(it_hpm->pitch);
      XYZVector.push_back(it_hpm->xyz);
      kinetic_energy += dEdxVector.back() * pitchVector.back();
    }
  }
  else {
    for (auto it_hpm = hpm.rbegin(); it_hpm != hpm.rend(); ++it_hpm) {
      dEdxVector.push_back(it_hpm->dEdx);
      dQdxVector.push_back(it_hpm->dQdx);
      resRangeVector.push_back((1 - it_hpm->path_fraction) * track_length);
      pitchVector.push_back(it_hpm->pitch);
      XYZVector.push_back(it_hpm->xyz);
      kinetic_energy += dEdxVector.back() * pitchVector.back();
    }
  }

  caloVector.emplace_back(kinetic_energy,
                          dEdxVector,
                          dQdxVector,
                          resRangeVector,
                          deadWireVector,
                          track_length,
                          pitchVector,
                          recob::tracking::convertCollToPoint(XYZVector),
                          planeID);
  assnTrackCaloVector.emplace_back(i_track);
}

void calo::TrackCalorimetryAlg::PrintHitPropertiesMultiset(HitPropertiesMultiset_t const& hpm)
{

  for (auto const& hit : hpm)
    hit.Print();

  std::cout << "Inverted? " << IsInvertedTrack(hpm) << std::endl;
}