analysis Namespace Reference

Classes
class	CalorimetryAnalysis

Functions
void	TrkDirectionAtXYZ (const recob::Track trk, const double x, const double y, const double z, float out[3])
void	True2RecoMappingXYZ (float t, float x, float y, float z, float out[3])
float	XYZtoPlanecoordinate (const float x, const float y, const float z, const int plane)
float	distance3d (const float &x1, const float &y1, const float &z1, const float &x2, const float &y2, const float &z2)
float	ContainedLength (const TVector3 &v0, const TVector3 &v1, const std::vector< geoalgo::AABox > &boxes)
void	FillHits (const detinfo::DetectorClocksData &clockData, const std::vector< art::Ptr< recob::Hit >> &hits, std::vector< float > &charge_u, std::vector< float > &charge_v, std::vector< float > &charge_y, std::vector< unsigned > &wire_u, std::vector< unsigned > &wire_v, std::vector< unsigned > &wire_y, std::vector< unsigned > &channel_u, std::vector< unsigned > &channel_v, std::vector< unsigned > &channel_y, std::vector< int > &multiplicity_u, std::vector< int > &multiplicity_v, std::vector< int > &multiplicity_y, std::vector< float > &width_u, std::vector< float > &width_v, std::vector< float > &width_y, std::vector< float > &time_u, std::vector< float > &time_v, std::vector< float > &time_y, std::vector< unsigned > &nhit_u, std::vector< unsigned > &nhit_v, std::vector< unsigned > &nhit_y, bool use_integral=true)
std::vector< unsigned >	FinddEdxHits (const anab::Calorimetry &calo, const std::vector< art::Ptr< recob::Hit >> &hits)
float	calcPitch (const geo::GeometryCore *geo, const geo::PlaneID &plane, TVector3 location, TVector3 direction, const spacecharge::SpaceCharge *SCEService=NULL)

Function Documentation

float analysis::calcPitch	(	const geo::GeometryCore *	geo,
		const geo::PlaneID &	plane,
		TVector3	location,
		TVector3	direction,
		const spacecharge::SpaceCharge *	SCEService = NULL
	)

Definition at line 1543 of file CalorimetryAnalysis_module.cc.

                                                                                                                                                           {
  float angletovert = geo->WireAngleToVertical(geo->View(plane), plane) - 0.5*::util::pi<>();

  // apply space charge change to dir
  if (SCEService) {
    geo::Point_t location_p(location);
    geo::Vector_t offset_v = SCEService->GetPosOffsets(location_p);
    TVector3 offset = TVector3(offset_v.X(), offset_v.Y(), offset_v.Z());
    // fix the x-sign
    if (location.X() < 0.) offset.SetX(-offset.X());

    // location a ~wire's distance along the track 
    TVector3 location_dx = location + geo->WirePitch(plane) * direction;

    // Apply space charge to the offset location
    geo::Point_t location_dx_p(location_dx);
    geo::Vector_t offset_dx_v = SCEService->GetPosOffsets(location_dx_p);
    TVector3 offset_dx(offset_dx_v.X(), offset_dx_v.Y(), offset_dx_v.Z());
    // fix the x-sign
    if (location_dx.X() < 0.) offset_dx.SetX(-offset_dx.X());

    // this is the actual direction
    direction = (offset_dx + location_dx - offset - location).Unit();
  }

  // get the pitch
  float cosgamma = abs(cos(angletovert) * direction.Z() + sin(angletovert) * direction.Y());
  float pitch = geo->WirePitch(plane) / cosgamma;

  // map pitch back onto particle trajectory
  if (SCEService) {
    geo::Point_t location_p(location);
    geo::Vector_t offset_v = SCEService->GetPosOffsets(location_p);
    TVector3 offset(offset_v.X(), offset_v.Y(), offset_v.Z());
    // fix the x-sign
    if (location.X() < 0.) offset.SetX(-offset.X());

    TVector3 location_sce = location + offset;

    TVector3 location_sce_dx = location_sce + direction * pitch;

    // map this back onto the particle trajectory
    geo::Point_t location_sce_dx_p(location_sce_dx);
    geo::Vector_t back_v = SCEService->GetCalPosOffsets(location_sce_dx_p, plane.TPC);
    TVector3 back(back_v.X(), back_v.Y(), back_v.Z());
    TVector3 location_dx = location_sce_dx + back;

    pitch = (location - location_dx).Mag();
  }

  return pitch;
}

float analysis::ContainedLength	(	const TVector3 &	v0,
		const TVector3 &	v1,
		const std::vector< geoalgo::AABox > &	boxes
	)

Definition at line 2509 of file CalorimetryAnalysis_module.cc.

                                                             {
  static const geoalgo::GeoAlgo algo;

  // if points are the same, return 0
  if ((v0 - v1).Mag() < 1e-6) return 0;

  // construct individual points
  geoalgo::Point_t p0(v0);
  geoalgo::Point_t p1(v1);

  // construct line segment
  geoalgo::LineSegment line(p0, p1);

  double length = 0;

  // total contained length is sum of lengths in all boxes
  // assuming they are non-overlapping
  for (auto const &box: boxes) {
    int n_contained = box.Contain(p0) + box.Contain(p1);
    // both points contained -- length is total length (also can break out of loop)
    if (n_contained == 2) {
      length = (v1 - v0).Mag();
      break;
    }
    // one contained -- have to find intersection point (which must exist)
    if (n_contained == 1) {
      auto intersections = algo.Intersection(line, box);
      // Because of floating point errors, it can sometimes happen
      // that there is 1 contained point but no "Intersections"
      // if one of the points is right on the edge
      if (intersections.size() == 0) {
        // determine which point is on the edge
        double tol = 1e-5;
        bool p0_edge = algo.SqDist(p0, box) < tol;
        bool p1_edge = algo.SqDist(p1, box) < tol;
        assert(p0_edge || p1_edge);
        // contained one is on edge -- can treat both as not contained
        //
        // In this case, no length
        if ((p0_edge && box.Contain(p0)) || (box.Contain(p1) && p1_edge))
          continue;
        // un-contaned one is on edge -- treat both as contained
        else if ((p0_edge && box.Contain(p1)) || (box.Contain(p0) && p1_edge)) {
          length = (v1 - v0).Mag();
          break;
        }
        else {
          assert(false); // bad
        }
      }
      // floating point errors can also falsely cause 2 intersection points
      //
      // in this case, one of the intersections must be very close to the 
      // "contained" point, so the total contained length will be about
      // the same as the distance between the two intersection points
      else if (intersections.size() == 2) {
        length += (intersections.at(0).ToTLorentzVector().Vect() - intersections.at(1).ToTLorentzVector().Vect()).Mag();
        continue;
      }
      // "Correct"/ideal case -- 1 intersection point
      else if (intersections.size() == 1) {
        // get TVector at intersection point
        TVector3 int_tv(intersections.at(0).ToTLorentzVector().Vect());
        length += ( box.Contain(p0) ? (v0 - int_tv).Mag() : (v1 - int_tv).Mag() ); 
      }
      else assert(false); // bad
    }
    // none contained -- either must have zero or two intersections
    if (n_contained == 0) {
      auto intersections = algo.Intersection(line, box);
      if (!(intersections.size() == 0 || intersections.size() == 2)) {
        // more floating point error fixes...
        //
        // figure out which points are near the edge
        double tol = 1e-5;
        bool p0_edge = algo.SqDist(p0, box) < tol;
        bool p1_edge = algo.SqDist(p1, box) < tol;
        // and which points are near the intersection
        TVector3 vint = intersections.at(0).ToTLorentzVector().Vect();

        bool p0_int = (v0 - vint).Mag() < tol;
        bool p1_int = (v1 - vint).Mag() < tol;
        // exactly one of them should produce the intersection
        assert((p0_int && p0_edge) != (p1_int && p1_edge));
        // void variables when assert-ions are turned off
        (void) p0_int; (void) p1_int;

        // both close to edge -- full length is contained
        if (p0_edge && p1_edge) {
          length += (v0 - v1).Mag();
        }
        // otherwise -- one of them is not on an edge, no length is contained
        else {}
      }
      // assert(intersections.size() == 0 || intersections.size() == 2);
      else if (intersections.size() == 2) {
        TVector3 start(intersections.at(0).ToTLorentzVector().Vect());
        TVector3 end(intersections.at(1).ToTLorentzVector().Vect());
        length += (start - end).Mag();
      }
    }
  }

  return length;
}//ContainedLength

float analysis::distance3d	(	const float &	x1,
		const float &	y1,
		const float &	z1,
		const float &	x2,
		const float &	y2,
		const float &	z2
	)

Definition at line 2500 of file CalorimetryAnalysis_module.cc.

  {
    return sqrt((x1-x2)*(x1-x2) +
                (y1-y2)*(y1-y2) +
                (z1-z2)*(z1-z2));

  }

void analysis::FillHits	(	const detinfo::DetectorClocksData &	clockData,
		const std::vector< art::Ptr< recob::Hit >> &	hits,
		std::vector< float > &	charge_u,
		std::vector< float > &	charge_v,
		std::vector< float > &	charge_y,
		std::vector< unsigned > &	wire_u,
		std::vector< unsigned > &	wire_v,
		std::vector< unsigned > &	wire_y,
		std::vector< unsigned > &	channel_u,
		std::vector< unsigned > &	channel_v,
		std::vector< unsigned > &	channel_y,
		std::vector< int > &	multiplicity_u,
		std::vector< int > &	multiplicity_v,
		std::vector< int > &	multiplicity_y,
		std::vector< float > &	width_u,
		std::vector< float > &	width_v,
		std::vector< float > &	width_y,
		std::vector< float > &	time_u,
		std::vector< float > &	time_v,
		std::vector< float > &	time_y,
		std::vector< unsigned > &	nhit_u,
		std::vector< unsigned > &	nhit_v,
		std::vector< unsigned > &	nhit_y,
		bool	use_integral = true
	)

Definition at line 2364 of file CalorimetryAnalysis_module.cc.

                                 {

  struct HitInfo {
    float charge;
    int multiplicity;
    float width;
    int start;
    int end;
    float time;
    int nhit;
    int wire;
    HitInfo(): charge(0), multiplicity(0), width(0), start(0), end(0), time(0.), nhit(0) {}
  };

  std::map<unsigned, HitInfo> hit_map_u;
  std::map<unsigned, HitInfo> hit_map_v;
  std::map<unsigned, HitInfo> hit_map_y;

  for (const art::Ptr<recob::Hit> hit: hits) {
    if (hit->WireID().Plane == 0) {
      // if using summed ADC method to count charge, don't double-count hits in the same snippet
      if (!use_integral && hit_map_u[hit->Channel()].start == hit->StartTick() && hit_map_u[hit->Channel()].end == hit->EndTick()) {
        continue;
      }
      hit_map_u[hit->Channel()].charge += (use_integral) ? hit->Integral() : hit->SummedADC();
      hit_map_u[hit->Channel()].multiplicity = std::max(hit_map_u[hit->Channel()].multiplicity, (int) hit->Multiplicity());
      hit_map_u[hit->Channel()].start = hit_map_u[hit->Channel()].nhit ? std::min(hit_map_u[hit->Channel()].start, hit->StartTick()) : hit->StartTick();
      hit_map_u[hit->Channel()].end = hit_map_u[hit->Channel()].nhit ? std::max(hit_map_u[hit->Channel()].end, hit->EndTick()) : hit->EndTick();
      hit_map_u[hit->Channel()].width = (hit_map_u[hit->Channel()].end - hit_map_u[hit->Channel()].start);
      hit_map_u[hit->Channel()].time = (hit->PeakTime() + hit_map_u[hit->Channel()].time * hit_map_u[hit->Channel()].nhit) / (hit_map_u[hit->Channel()].nhit + 1);
      hit_map_u[hit->Channel()].nhit ++;
      hit_map_u[hit->Channel()].wire = hit->WireID().Wire;
    }
    else if (hit->WireID().Plane == 1) {
      // if using summed ADC method to count charge, don't double-count hits in the same snippet
      if (!use_integral && hit_map_v[hit->Channel()].start == hit->StartTick() && hit_map_v[hit->Channel()].end == hit->EndTick()) {
        continue;
      }
      hit_map_v[hit->Channel()].charge += (use_integral) ? hit->Integral() : hit->SummedADC();
      hit_map_v[hit->Channel()].multiplicity = std::max(hit_map_v[hit->Channel()].multiplicity, (int) hit->Multiplicity());
      hit_map_v[hit->Channel()].start = hit_map_v[hit->Channel()].nhit ? std::min(hit_map_v[hit->Channel()].start, hit->StartTick()) : hit->StartTick();
      hit_map_v[hit->Channel()].end = hit_map_v[hit->Channel()].nhit ? std::max(hit_map_v[hit->Channel()].end, hit->EndTick()) : hit->EndTick();
      hit_map_v[hit->Channel()].width = (hit_map_v[hit->Channel()].end - hit_map_v[hit->Channel()].start);
      hit_map_v[hit->Channel()].time = (hit->PeakTime() + hit_map_v[hit->Channel()].time * hit_map_v[hit->Channel()].nhit) / (hit_map_v[hit->Channel()].nhit + 1);
      hit_map_v[hit->Channel()].nhit ++;
      hit_map_v[hit->Channel()].wire = hit->WireID().Wire;
    }
    else {
      // if using summed ADC method to count charge, don't double-count hits in the same snippet
      if (!use_integral && hit_map_y[hit->Channel()].start == hit->StartTick() && hit_map_y[hit->Channel()].end == hit->EndTick()) {
        continue;
      }
      hit_map_y[hit->Channel()].charge += (use_integral) ? hit->Integral() : hit->SummedADC();
      hit_map_y[hit->Channel()].multiplicity = std::max(hit_map_y[hit->Channel()].multiplicity, (int) hit->Multiplicity());
      hit_map_y[hit->Channel()].start = hit_map_y[hit->Channel()].nhit ? std::min(hit_map_y[hit->Channel()].start, hit->StartTick()) : hit->StartTick();
      hit_map_y[hit->Channel()].end = hit_map_y[hit->Channel()].nhit ? std::max(hit_map_y[hit->Channel()].end, hit->EndTick()) : hit->EndTick();
      hit_map_y[hit->Channel()].width = (hit_map_y[hit->Channel()].end - hit_map_y[hit->Channel()].start);
      hit_map_y[hit->Channel()].time = (hit->PeakTime() + hit_map_y[hit->Channel()].time * hit_map_y[hit->Channel()].nhit) / (hit_map_y[hit->Channel()].nhit + 1);
      hit_map_y[hit->Channel()].nhit ++;
      hit_map_y[hit->Channel()].wire = hit->WireID().Wire;
    }
  }

  for (auto const &pair: hit_map_u) {
    channel_u.push_back(pair.first);
    charge_u.push_back(pair.second.charge);
    multiplicity_u.push_back(pair.second.multiplicity);
    width_u.push_back(pair.second.width);
    time_u.push_back(pair.second.time);
    nhit_u.push_back(pair.second.nhit);
    wire_u.push_back(pair.second.wire);
  }
  for (auto const &pair: hit_map_v) {
    channel_v.push_back(pair.first);
    charge_v.push_back(pair.second.charge);
    multiplicity_v.push_back(pair.second.multiplicity);
    width_v.push_back(pair.second.width);
    time_v.push_back(pair.second.time);
    nhit_v.push_back(pair.second.nhit);
    wire_v.push_back(pair.second.wire);
  }
  for (auto const &pair: hit_map_y) {
    channel_y.push_back(pair.first);
    charge_y.push_back(pair.second.charge);
    multiplicity_y.push_back(pair.second.multiplicity);
    width_y.push_back(pair.second.width);
    time_y.push_back(pair.second.time);
    nhit_y.push_back(pair.second.nhit);
    wire_y.push_back(pair.second.wire);
  }

}

std::vector<unsigned> analysis::FinddEdxHits	(	const anab::Calorimetry &	calo,
		const std::vector< art::Ptr< recob::Hit >> &	hits
	)

Definition at line 1522 of file CalorimetryAnalysis_module.cc.

                                                                                                       {
  std::vector<unsigned> ret;

  // NOTE: TpIndices are __actually__ the art::Ptr keys of the Hit's
  // Don't ask me why
  for (size_t ind: calo.TpIndices()) {
    bool found = false;
    for (const art::Ptr<recob::Hit> &hit: hits) {
      if (hit.key() == ind) {
        ret.push_back(hit->Channel());
        found = true;
        break;
      }
    }
    if (!found) {
      throw cet::exception("dEdx w/out Hit!");
    }
  }
  return ret;
}

void analysis::TrkDirectionAtXYZ	(	const recob::Track	trk,
		const double	x,
		const double	y,
		const double	z,
		float	out[3]
	)

Definition at line 2327 of file CalorimetryAnalysis_module.cc.

  {
    float min_dist = 1000;
    size_t i_min = -1;
    for(size_t i=0; i < trk.NumberTrajectoryPoints(); i++)
    {
      if (trk.HasValidPoint(i))
      { // check this point is valid
        auto point_i = trk.LocationAtPoint(i);
        float distance = distance3d((double)point_i.X(), (double)point_i.Y(), (double)point_i.Z(),
                  x, y, z);
        if (distance < min_dist)
        {
          min_dist = distance;
          i_min = i;
        }
      }// if point is valid
    }// for all track points

    if (i_min == (size_t)-1) return;

    auto direction = trk.DirectionAtPoint(i_min);
    out[0] = (float)direction.X();
    out[1] = (float)direction.Y();
    out[2] = (float)direction.Z();

    float norm;
    norm = out[0]*out[0] + out[1]*out[1] + out[2]*out[2];
    if (fabs(norm -1) > 0.001)
    {
      std::cout << "i_min = " << i_min << std::endl;
      std::cout << "minimum distance = " << min_dist << std::endl;
      std::cout << "out[0], out[1], out[2] = " << out[0] << " , " << out[1] << " , " << out[2] << std::endl;
      std::cout << "norm = " << norm << std::endl;
    }
  }

void analysis::True2RecoMappingXYZ	(	float	t,
		float	x,
		float	y,
		float	z,
		float	out[3]
	)

Definition at line 2481 of file CalorimetryAnalysis_module.cc.

                                                                          {
  (void) t;
  out[0] = x;
  out[1] = y;
  out[2] = z;
}

float analysis::XYZtoPlanecoordinate	(	const float	x,
		const float	y,
		const float	z,
		const int	plane
	)

Definition at line 2488 of file CalorimetryAnalysis_module.cc.

  {
    auto const* geom = ::lar::providerFrom<geo::Geometry>();
    geo::Point_t p(x,y,z);
    geo::TPCID tpc = geom->FindTPCAtPosition(p);
    if (!tpc) return -1e6;
    geo::PlaneID planeID(tpc, plane);
    double _wire2cm = geom->WirePitch(planeID);

    return geom->WireCoordinate(p, planeID) * _wire2cm;
  }