#include "LArPandoraShowerAlg.h"

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

void	OrderShowerHits (detinfo::DetectorPropertiesData const &detProp, std::vector< art::Ptr< recob::Hit >> &hits, geo::Point_t const &ShowerPosition, geo::Vector_t const &ShowerDirection) const

void	OrderShowerSpacePointsPerpendicular (std::vector< art::Ptr< recob::SpacePoint >> &showersps, geo::Point_t const &vertex, geo::Vector_t const &direction) const

void	OrderShowerSpacePoints (std::vector< art::Ptr< recob::SpacePoint >> &showersps, geo::Point_t const &vertex, geo::Vector_t const &direction) const

void	OrderShowerSpacePoints (std::vector< art::Ptr< recob::SpacePoint >> &showersps, geo::Point_t const &vertex) const

geo::Point_t	ShowerCentre (std::vector< art::Ptr< recob::SpacePoint >> const &showersps) const

geo::Point_t	ShowerCentre (detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, std::vector< art::Ptr< recob::SpacePoint >> const &showersps, art::FindManyP< recob::Hit > const &fmh, float &totalCharge) const

geo::Point_t	ShowerCentre (detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, std::vector< art::Ptr< recob::SpacePoint >> const &showerspcs, art::FindManyP< recob::Hit > const &fmh) const

double	DistanceBetweenSpacePoints (art::Ptr< recob::SpacePoint > const &sp_a, art::Ptr< recob::SpacePoint > const &sp_b) const

double	SpacePointCharge (art::Ptr< recob::SpacePoint > const &sp, art::FindManyP< recob::Hit > const &fmh) const

double	SpacePointTime (art::Ptr< recob::SpacePoint > const &sp, art::FindManyP< recob::Hit > const &fmh) const

TVector2	HitCoordinates (detinfo::DetectorPropertiesData const &detProp, art::Ptr< recob::Hit > const &hit) const

double	SpacePointProjection (art::Ptr< recob::SpacePoint > const &sp, geo::Point_t const &vertex, geo::Vector_t const &direction) const

double	SpacePointPerpendicular (art::Ptr< recob::SpacePoint > const &sp, geo::Point_t const &vertex, geo::Vector_t const &direction) const

double	SpacePointPerpendicular (art::Ptr< recob::SpacePoint > const &sp, geo::Point_t const &vertex, geo::Vector_t const &direction, double proj) const

double	RMSShowerGradient (std::vector< art::Ptr< recob::SpacePoint >> &sps, const geo::Point_t &ShowerCentre, const geo::Vector_t &Direction, const unsigned int nSegments) const

double	CalculateRMS (const std::vector< float > &perps) const

double	SCECorrectPitch (double const &pitch, geo::Point_t const &pos, geo::Vector_t const &dir, unsigned int const &TPC) const

double	SCECorrectEField (double const &EField, geo::Point_t const &pos) const

std::map< art::Ptr< recob::Hit >, std::vector< art::Ptr< recob::Hit > > >	OrganizeHits (const std::vector< art::Ptr< recob::Hit >> &hits) const

void	DebugEVD (art::Ptr< recob::PFParticle > const &pfparticle, art::Event const &Event, const reco::shower::ShowerElementHolder &ShowerEleHolder, std::string const &evd_disp_name_append="") const

Private Attributes
bool	fUseCollectionOnly

art::InputTag	fPFParticleLabel

bool	fSCEXFlip

spacecharge::SpaceCharge const *	fSCE

geo::WireReadoutGeom const *	fChannelMap = &art::ServiceHandle<geo::WireReadout>()->Get()

art::ServiceHandle< art::TFileService >	tfs

const std::string	fInitialTrackInputLabel

const std::string	fShowerStartPositionInputLabel

const std::string	fShowerDirectionInputLabel

const std::string	fInitialTrackSpacePointsInputLabel

Detailed Description

Definition at line 40 of file LArPandoraShowerAlg.h.

Constructor & Destructor Documentation

shower::LArPandoraShowerAlg::LArPandoraShowerAlg ( const fhicl::ParameterSet & pset )

explicit

Definition at line 26 of file LArPandoraShowerAlg.cxx.

   : fUseCollectionOnly(pset.get<bool>("UseCollectionOnly"))
   , fPFParticleLabel(pset.get<art::InputTag>("PFParticleLabel"))
   , fSCEXFlip(pset.get<bool>("SCEXFlip"))
   , fSCE(lar::providerFrom<spacecharge::SpaceChargeService>())
   , fInitialTrackInputLabel(pset.get<std::string>("InitialTrackInputLabel"))
   , fShowerStartPositionInputLabel(pset.get<std::string>("ShowerStartPositionInputLabel"))
   , fShowerDirectionInputLabel(pset.get<std::string>("ShowerDirectionInputLabel"))
   , fInitialTrackSpacePointsInputLabel(pset.get<std::string>("InitialTrackSpacePointsInputLabel"))
 {}

Member Function Documentation

double shower::LArPandoraShowerAlg::CalculateRMS ( const std::vector< float > & perps ) const

Definition at line 466 of file LArPandoraShowerAlg.cxx.

References sum.

Referenced by RMSShowerGradient().

 {
 
   if (perps.size() < 2) return std::numeric_limits<double>::lowest();
 
   float sum = 0;
   for (const auto& perp : perps) {
     sum += perp * perp;
   }
 
   return std::sqrt(sum / (perps.size() - 1));
 }

void shower::LArPandoraShowerAlg::DebugEVD	(	art::Ptr< recob::PFParticle > const &	pfparticle,
		art::Event const &	Event,
		const reco::shower::ShowerElementHolder &	ShowerEleHolder,
		std::string const &	evd_disp_name_append = `""`
	)		const

Definition at line 601 of file LArPandoraShowerAlg.cxx.

References util::abs(), reco::shower::ShowerElementHolder::CheckElement(), art::Event::event(), art::fill_ptr_vector(), fInitialTrackInputLabel, fInitialTrackSpacePointsInputLabel, recob::Track::FlagsAtPoint(), fPFParticleLabel, fShowerDirectionInputLabel, fShowerStartPositionInputLabel, reco::shower::ShowerElementHolder::GetElement(), art::ProductRetriever::getValidHandle(), art::Ptr< T >::key(), recob::Track::LocationAtPoint(), recob::TrajectoryPointFlagTraits::NoPoint, recob::Track::NumberTrajectoryPoints(), proj, art::Event::run(), recob::PFParticle::Self(), SpacePointProjection(), art::Event::subRun(), tfs, x, x_max, x_min, y, and z.

 {
 
   std::cout << "Making Debug Event Display" << std::endl;
 
   //Function for drawing reco showers to check direction and initial track selection
 
   // Get run info to make unique canvas names
   int run = Event.run();
   int subRun = Event.subRun();
   int event = Event.event();
   int PFPID = pfparticle->Self();
 
   // Create the canvas
   TString canvasName = Form("canvas_%i_%i_%i_%i", run, subRun, event, PFPID);
   if (evd_disp_name_append.length() > 0) canvasName += "_" + evd_disp_name_append;
   TCanvas* canvas = tfs->make<TCanvas>(canvasName, canvasName);
 
   // Initialise variables
   double x = 0;
   double y = 0;
   double z = 0;
 
   double x_min = std::numeric_limits<double>::max(), x_max = -std::numeric_limits<double>::max();
   double y_min = std::numeric_limits<double>::max(), y_max = -std::numeric_limits<double>::max();
   double z_min = std::numeric_limits<double>::max(), z_max = -std::numeric_limits<double>::max();
 
   // Get a bunch of associations (again)
   // N.B. this is a horribly inefficient way of doing things but as this is only
   // going to be used to debug I don't care, I would rather have generality in this case
 
   auto const pfpHandle = Event.getValidHandle<std::vector<recob::PFParticle>>(fPFParticleLabel);
 
   // Get the spacepoint - PFParticle assn
   art::FindManyP<recob::SpacePoint> fmspp(pfpHandle, Event, fPFParticleLabel);
   if (!fmspp.isValid()) {
     throw cet::exception("LArPandoraShowerAlg") << "Trying to get the spacepoint and failed. Somet\
       hing is not configured correctly. Stopping ";
   }
 
   // Get the SpacePoints
   std::vector<art::Ptr<recob::SpacePoint>> const& spacePoints = fmspp.at(pfparticle.key());
 
   //We cannot progress with no spacepoints.
   if (spacePoints.empty()) { return; }
 
   // Get info from shower property holder
   geo::Point_t showerStartPosition = {-999, -999, -999};
   geo::Vector_t showerDirection = {-999, -999, -999};
   std::vector<art::Ptr<recob::SpacePoint>> trackSpacePoints;
 
   //######################
   //### Start Position ###
   //######################
   double startXYZ[3] = {-999, -999, -999};
   if (!ShowerEleHolder.CheckElement(fShowerStartPositionInputLabel)) {
     mf::LogError("LArPandoraShowerAlg") << "Start position not set, returning " << std::endl;
     // return;
   }
   else {
     ShowerEleHolder.GetElement(fShowerStartPositionInputLabel, showerStartPosition);
     // Create 3D point at vertex, chosed to be origin for ease of use of display
     startXYZ[0] = showerStartPosition.X();
     startXYZ[1] = showerStartPosition.Y();
     startXYZ[2] = showerStartPosition.Z();
   }
   auto startPoly = std::make_unique<TPolyMarker3D>(1, startXYZ);
 
   //########################
   //### Shower Direction ###
   //########################
 
   double xDirPoints[2] = {-999, -999};
   double yDirPoints[2] = {-999, -999};
   double zDirPoints[2] = {-999, -999};
 
   //initialise counter point
   int point = 0;
 
   // Make 3D points for each spacepoint in the shower
   auto allPoly = std::make_unique<TPolyMarker3D>(spacePoints.size());
 
   if (!ShowerEleHolder.CheckElement(fShowerDirectionInputLabel) &&
       !ShowerEleHolder.CheckElement("ShowerStartPosition")) {
     mf::LogError("LArPandoraShowerAlg") << "Direction not set, returning " << std::endl;
   }
   else {
 
     // Get the min and max projections along the direction to know how long to draw
     ShowerEleHolder.GetElement(fShowerDirectionInputLabel, showerDirection);
 
     // the direction line
     double minProj = std::numeric_limits<double>::max();
     double maxProj = -std::numeric_limits<double>::max();
 
     //initialise counter point
     int point = 0;
 
     for (auto spacePoint : spacePoints) {
       auto const pos = spacePoint->position();
       x = pos.X();
       y = pos.Y();
       z = pos.Z();
       allPoly->SetPoint(point, x, y, z);
       ++point;
 
       x_min = std::min(x, x_min);
       x_max = std::max(x, x_max);
       y_min = std::min(y, y_min);
       y_max = std::max(y, y_max);
       z_min = std::min(z, z_min);
       z_max = std::max(z, z_max);
 
       // Calculate the projection of (point-startpoint) along the direction
       double proj = shower::LArPandoraShowerAlg::SpacePointProjection(
         spacePoint, showerStartPosition, showerDirection);
       maxProj = std::max(proj, maxProj);
       minProj = std::min(proj, minProj);
     } // loop over spacepoints
 
     xDirPoints[0] = (showerStartPosition.X() + minProj * showerDirection.X());
     xDirPoints[1] = (showerStartPosition.X() + maxProj * showerDirection.X());
 
     yDirPoints[0] = (showerStartPosition.Y() + minProj * showerDirection.Y());
     yDirPoints[1] = (showerStartPosition.Y() + maxProj * showerDirection.Y());
 
     zDirPoints[0] = (showerStartPosition.Z() + minProj * showerDirection.Z());
     zDirPoints[1] = (showerStartPosition.Z() + maxProj * showerDirection.Z());
   }
 
   auto dirPoly = std::make_unique<TPolyLine3D>(2, xDirPoints, yDirPoints, zDirPoints);
 
   //#########################
   //### Initial Track SPs ###
   //#########################
 
   auto trackPoly = std::make_unique<TPolyMarker3D>(trackSpacePoints.size());
   if (!ShowerEleHolder.CheckElement(fInitialTrackSpacePointsInputLabel)) {
     mf::LogError("LArPandoraShowerAlg") << "TrackSpacePoints not set, returning " << std::endl;
     //    return;
   }
   else {
     ShowerEleHolder.GetElement(fInitialTrackSpacePointsInputLabel, trackSpacePoints);
     point = 0; // re-initialise counter
     for (auto spacePoint : trackSpacePoints) {
       auto const pos = spacePoint->position();
       x = pos.X();
       y = pos.Y();
       z = pos.Z();
       trackPoly->SetPoint(point, x, y, z);
       ++point;
       x_min = std::min(x, x_min);
       x_max = std::max(x, x_max);
       y_min = std::min(y, y_min);
       y_max = std::max(y, y_max);
       z_min = std::min(z, z_min);
       z_max = std::max(z, z_max);
     } // loop over track spacepoints
   }
 
   //#########################
   //### Other PFParticles ###
   //#########################
 
   //  we want to draw all of the PFParticles in the event
   //Get the PFParticles
   std::vector<art::Ptr<recob::PFParticle>> pfps;
   art::fill_ptr_vector(pfps, pfpHandle);
 
   // initialse counters
   // Split into tracks and showers to make it clearer what pandora is doing
   int pfpTrackCounter = 0;
   int pfpShowerCounter = 0;
 
   // initial loop over pfps to find nuber of spacepoints for tracks and showers
   for (auto const& pfp : pfps) {
     std::vector<art::Ptr<recob::SpacePoint>> const& sps = fmspp.at(pfp.key());
     // If running pandora cheating it will call photons pdg 22
     int pdg = abs(pfp->PdgCode()); // Track or shower
     if (pdg == 11 || pdg == 22) { pfpShowerCounter += sps.size(); }
     else {
       pfpTrackCounter += sps.size();
     }
   }
 
   auto pfpPolyTrack = std::make_unique<TPolyMarker3D>(pfpTrackCounter);
   auto pfpPolyShower = std::make_unique<TPolyMarker3D>(pfpShowerCounter);
 
   // initialise counters
   int trackPoints = 0;
   int showerPoints = 0;
 
   for (auto const& pfp : pfps) {
     std::vector<art::Ptr<recob::SpacePoint>> const& sps = fmspp.at(pfp.key());
     int pdg = abs(pfp->PdgCode()); // Track or shower
     for (auto const& sp : sps) {
       auto const pos = sp->position();
       x = pos.X();
       y = pos.Y();
       z = pos.Z();
       x_min = std::min(x, x_min);
       x_max = std::max(x, x_max);
       y_min = std::min(y, y_min);
       y_max = std::max(y, y_max);
       z_min = std::min(z, z_min);
       z_max = std::max(z, z_max);
 
       // If running pandora cheating it will call photons pdg 22
       if (pdg == 11 || pdg == 22) {
         pfpPolyShower->SetPoint(showerPoints, x, y, z);
         ++showerPoints;
       }
       else {
         pfpPolyTrack->SetPoint(trackPoints, x, y, z);
         ++trackPoints;
       }
     } // loop over sps
 
     //pfpPolyTrack->Draw();
 
   } // if (fDrawAllPFPs)
 
   //#################################
   //### Initial Track Traj Points ###
   //#################################
 
   auto TrackTrajPoly = std::make_unique<TPolyMarker3D>(TPolyMarker3D(1));
   auto TrackInitTrajPoly = std::make_unique<TPolyMarker3D>(TPolyMarker3D(1));
 
   if (ShowerEleHolder.CheckElement(fInitialTrackInputLabel)) {
 
     //Get the track
     recob::Track InitialTrack;
     ShowerEleHolder.GetElement(fInitialTrackInputLabel, InitialTrack);
 
     if (InitialTrack.NumberTrajectoryPoints() != 0) {
 
       point = 0;
       // Make 3D points for each trajectory point in the track stub
       for (unsigned int traj = 0; traj < InitialTrack.NumberTrajectoryPoints(); ++traj) {
 
         //ignore bogus info.
         auto flags = InitialTrack.FlagsAtPoint(traj);
         if (flags.isSet(recob::TrajectoryPointFlagTraits::NoPoint)) { continue; }
 
         geo::Point_t TrajPositionPoint = InitialTrack.LocationAtPoint(traj);
         x = TrajPositionPoint.X();
         y = TrajPositionPoint.Y();
         z = TrajPositionPoint.Z();
         TrackTrajPoly->SetPoint(point, x, y, z);
         ++point;
       } // loop over trajectory points
 
       geo::Point_t TrajInitPositionPoint = InitialTrack.LocationAtPoint(0);
       x = TrajInitPositionPoint.X();
       y = TrajInitPositionPoint.Y();
       z = TrajInitPositionPoint.Z();
       TrackInitTrajPoly->SetPoint(TrackInitTrajPoly->GetN(), x, y, z);
     }
   }
 
   gStyle->SetOptStat(0);
   TH3F axes("axes", "", 1, x_min, x_max, 1, y_min, y_max, 1, z_min, z_max);
   axes.SetDirectory(0);
   axes.GetXaxis()->SetTitle("X");
   axes.GetYaxis()->SetTitle("Y");
   axes.GetZaxis()->SetTitle("Z");
   axes.Draw();
 
   // Draw all of the things
   pfpPolyShower->SetMarkerStyle(20);
   pfpPolyShower->SetMarkerColor(4);
   pfpPolyShower->Draw();
   pfpPolyTrack->SetMarkerStyle(20);
   pfpPolyTrack->SetMarkerColor(6);
   pfpPolyTrack->Draw();
   allPoly->SetMarkerStyle(20);
   allPoly->Draw();
   trackPoly->SetMarkerStyle(20);
   trackPoly->SetMarkerColor(2);
   trackPoly->Draw();
   startPoly->SetMarkerStyle(21);
   startPoly->SetMarkerSize(0.5);
   startPoly->SetMarkerColor(3);
   startPoly->Draw();
   dirPoly->SetLineWidth(1);
   dirPoly->SetLineColor(6);
   dirPoly->Draw();
   TrackTrajPoly->SetMarkerStyle(22);
   TrackTrajPoly->SetMarkerColor(7);
   TrackTrajPoly->Draw();
   TrackInitTrajPoly->SetMarkerStyle(22);
   TrackInitTrajPoly->SetMarkerColor(4);
   TrackInitTrajPoly->Draw();
 
   // Save the canvas. Don't usually need this when using TFileService but this in the alg
   // not a module and didn't work without this so im going with it.
   canvas->Write();
 }

double shower::LArPandoraShowerAlg::DistanceBetweenSpacePoints	(	art::Ptr< recob::SpacePoint > const &	sp_a,
		art::Ptr< recob::SpacePoint > const &	sp_b
	)		const

Definition at line 310 of file LArPandoraShowerAlg.cxx.

References recob::SpacePoint::position().

Referenced by ShowerRecoTools::ShowerIncrementalTrackHitFinder::PruneFrontOfSPSPool(), and ShowerRecoTools::ShowerIncrementalTrackHitFinder::PruneTrack().

 {
   return (sp_a->position() - sp_b->position()).R();
 }

TVector2 shower::LArPandoraShowerAlg::HitCoordinates	(	detinfo::DetectorPropertiesData const &	detProp,
		art::Ptr< recob::Hit > const &	hit
	)		const

Definition at line 351 of file LArPandoraShowerAlg.cxx.

References geo::WireID::asPlaneID(), detinfo::DetectorPropertiesData::ConvertTicksToX(), fChannelMap, recob::Hit::PeakTime(), geo::WireReadoutGeom::Plane(), geo::WireID::Wire, recob::Hit::WireID(), and geo::PlaneGeo::WirePitch().

Referenced by ShowerRecoTools::Shower2DLinearRegressionTrackHitFinder::FindInitialTrackHits(), and OrderShowerHits().

 {
   //Get the pitch
   const geo::WireID WireID = hit->WireID();
   const geo::PlaneID planeid = WireID.asPlaneID();
   double pitch = fChannelMap->Plane(planeid).WirePitch();
 
   return TVector2(WireID.Wire * pitch, detProp.ConvertTicksToX(hit->PeakTime(), planeid));
 }

void shower::LArPandoraShowerAlg::OrderShowerHits	(	detinfo::DetectorPropertiesData const &	detProp,
		std::vector< art::Ptr< recob::Hit >> &	hits,
		geo::Point_t const &	ShowerPosition,
		geo::Vector_t const &	ShowerDirection
	)		const

Definition at line 40 of file LArPandoraShowerAlg.cxx.

References util::abs(), geo::WireID::asPlaneID(), fChannelMap, HitCoordinates(), hits(), geo::WireReadoutGeom::Plane(), recob::Hit::WireID(), and geo::PlaneGeo::WirePitch().

Referenced by ShowerRecoTools::Shower2DLinearRegressionTrackHitFinder::CalculateElement().

 {
   std::map<double, art::Ptr<recob::Hit>> OrderedHits;
   art::Ptr<recob::Hit> startHit = hits.front();
 
   //Get the wireID
   const geo::WireID startWireID = startHit->WireID();
 
   //Get the plane
   const geo::PlaneID planeid = startWireID.asPlaneID();
 
   //Get the pitch
   auto const& plane = fChannelMap->Plane(planeid);
   double pitch = plane.WirePitch();
 
   TVector2 Shower2DStartPosition = {plane.WireCoordinate(ShowerStartPosition) * pitch,
                                     ShowerStartPosition.X()};
 
   //Vector of the plane
   auto const PlaneDirection = plane.GetIncreasingWireDirection();
 
   //get the shower 2D direction
   TVector2 Shower2DDirection = {ShowerDirection.Dot(PlaneDirection), ShowerDirection.X()};
 
   Shower2DDirection = Shower2DDirection.Unit();
 
   for (auto const& hit : hits) {
 
     //Get the wireID
     const geo::WireID WireID = hit->WireID();
 
     if (WireID.asPlaneID() != planeid) { break; }
 
     //Get the hit Vector.
     TVector2 hitcoord = HitCoordinates(detProp, hit);
 
     //Order the hits based on the projection
     TVector2 pos = hitcoord - Shower2DStartPosition;
     OrderedHits[pos * Shower2DDirection] = hit;
   }
 
   //Transform the shower.
   std::vector<art::Ptr<recob::Hit>> showerHits;
   std::transform(OrderedHits.begin(),
                  OrderedHits.end(),
                  std::back_inserter(showerHits),
                  [](std::pair<double, art::Ptr<recob::Hit>> const& hit) { return hit.second; });
 
   //Sometimes get the order wrong. Depends on direction compared to the plane Correct for it here:
   art::Ptr<recob::Hit> frontHit = showerHits.front();
   art::Ptr<recob::Hit> backHit = showerHits.back();
 
   //Get the hit Vector.
   TVector2 fronthitcoord = HitCoordinates(detProp, frontHit);
   TVector2 frontpos = fronthitcoord - Shower2DStartPosition;
 
   //Get the hit Vector.
   TVector2 backhitcoord = HitCoordinates(detProp, backHit);
   TVector2 backpos = backhitcoord - Shower2DStartPosition;
 
   double frontproj = frontpos * Shower2DDirection;
   double backproj = backpos * Shower2DDirection;
   if (std::abs(backproj) < std::abs(frontproj)) {
     std::reverse(showerHits.begin(), showerHits.end());
   }
 
   hits = showerHits;
   return;
 }

void shower::LArPandoraShowerAlg::OrderShowerSpacePoints	(	std::vector< art::Ptr< recob::SpacePoint >> &	showersps,
		geo::Point_t const &	vertex,
		geo::Vector_t const &	direction
	)		const

Definition at line 142 of file LArPandoraShowerAlg.cxx.

References SpacePointProjection().

Referenced by ShowerRecoTools::ShowerLengthPercentile::CalculateElement(), ShowerRecoTools::Shower3DCylinderTrackHitFinder::CalculateElement(), ShowerRecoTools::ShowerPFPVertexStartPosition::CalculateElement(), ShowerRecoTools::ShowerIncrementalTrackHitFinder::CalculateElement(), RMSShowerGradient(), and ShowerRecoTools::ShowerIncrementalTrackHitFinder::RunTestOfIncrementalSpacePointFinder().

 {
 
   std::map<double, art::Ptr<recob::SpacePoint>> OrderedSpacePoints;
 
   //Loop over the spacepoints and get the pojected distance from the vertex.
   for (auto const& sp : showersps) {
 
     // Get the projection of the space point along the direction
     double len = SpacePointProjection(sp, vertex, direction);
 
     //Add to the list
     OrderedSpacePoints[len] = sp;
   }
 
   //Return an ordered list.
   showersps.clear();
   for (auto const& sp : OrderedSpacePoints) {
     showersps.push_back(sp.second);
   }
 }

void shower::LArPandoraShowerAlg::OrderShowerSpacePoints	(	std::vector< art::Ptr< recob::SpacePoint >> &	showersps,
		geo::Point_t const &	vertex
	)		const

Definition at line 167 of file LArPandoraShowerAlg.cxx.

 {
   std::map<double, art::Ptr<recob::SpacePoint>> OrderedSpacePoints;
 
   //Loop over the spacepoints and get the pojected distance from the vertex.
   for (auto const& sp : showersps) {
 
     //Get the distance away from the start
     double mag = (sp->position() - vertex).R();
 
     //Add to the list
     OrderedSpacePoints[mag] = sp;
   }
 
   //Return an ordered list.
   showersps.clear();
   for (auto const& sp : OrderedSpacePoints) {
     showersps.push_back(sp.second);
   }
 }

void shower::LArPandoraShowerAlg::OrderShowerSpacePointsPerpendicular	(	std::vector< art::Ptr< recob::SpacePoint >> &	showersps,
		geo::Point_t const &	vertex,
		geo::Vector_t const &	direction
	)		const

Definition at line 115 of file LArPandoraShowerAlg.cxx.

References SpacePointPerpendicular().

Referenced by ShowerRecoTools::ShowerLengthPercentile::CalculateElement().

 {
 
   std::map<double, art::Ptr<recob::SpacePoint>> OrderedSpacePoints;
 
   //Loop over the spacepoints and get the pojected distance from the vertex.
   for (auto const& sp : showersps) {
 
     // Get the perpendicular distance
     double perp = SpacePointPerpendicular(sp, vertex, direction);
 
     //Add to the list
     OrderedSpacePoints[perp] = sp;
   }
 
   //Return an ordered list.
   showersps.clear();
   for (auto const& sp : OrderedSpacePoints) {
     showersps.push_back(sp.second);
   }
 }

std::map< art::Ptr< recob::Hit >, std::vector< art::Ptr< recob::Hit > > > shower::LArPandoraShowerAlg::OrganizeHits ( const std::vector< art::Ptr< recob::Hit >> & hits ) const

Definition at line 528 of file LArPandoraShowerAlg.cxx.

References recob::Hit::EndTick(), hits(), recob::Hit::Integral(), util::details::operator==(), util::quantities::concepts::operator>(), recob::Hit::StartTick(), geo::WireID::Wire, and recob::Hit::WireID().

Referenced by ShowerRecoTools::ShowerUnidirectiondEdx::CalculateElement(), and ShowerRecoTools::ShowerTrajPointdEdx::CalculateElement().

 {
   // In this case, we need to only accept one hit in each snippet
   // Snippets are counted by the Start, End, and Wire. If all these are the same for a hit, then they are on the same snippet.
   //
   // If there are multiple valid hits on the same snippet, we need a way to pick the best one.
   // (TODO: find a good way). The current method is to take the one with the highest charge integral.
   typedef std::pair<unsigned, std::vector<unsigned>> OrganizedHits;
   struct HitIdentifier {
     int startTick;
     int endTick;
     int wire;
     float integral;
 
     // construct
     explicit HitIdentifier(const recob::Hit& hit)
       : startTick(hit.StartTick())
       , endTick(hit.EndTick())
       , wire(hit.WireID().Wire)
       , integral(hit.Integral())
     {}
 
     // Defines whether two hits are on the same snippet
     inline bool operator==(const HitIdentifier& rhs) const
     {
       return startTick == rhs.startTick && endTick == rhs.endTick && wire == rhs.wire;
     }
 
     // Defines which hit to pick between two both on the same snippet
     inline bool operator>(const HitIdentifier& rhs) const { return integral > rhs.integral; }
   };
 
   unsigned nplanes = 6; // TODO FIXME
   std::vector<std::vector<OrganizedHits>> hits_org(nplanes);
   std::vector<std::vector<HitIdentifier>> hit_idents(nplanes);
   for (unsigned i = 0; i < hits.size(); i++) {
     HitIdentifier this_ident(*hits[i]);
 
     // check if we have found a hit on this snippet before
     bool found_snippet = false;
     for (unsigned j = 0; j < hits_org[hits[i]->WireID().Plane].size(); j++) {
       if (this_ident == hit_idents[hits[i]->WireID().Plane][j]) {
         found_snippet = true;
         if (this_ident > hit_idents[hits[i]->WireID().Plane][j]) {
           hits_org[hits[i]->WireID().Plane][j].second.push_back(
             hits_org[hits[i]->WireID().Plane][j].first);
           hits_org[hits[i]->WireID().Plane][j].first = i;
           hit_idents[hits[i]->WireID().Plane][j] = this_ident;
         }
         else {
           hits_org[hits[i]->WireID().Plane][j].second.push_back(i);
         }
         break;
       }
     }
     if (!found_snippet) {
       hits_org[hits[i]->WireID().Plane].push_back({i, {}});
       hit_idents[hits[i]->WireID().Plane].push_back(this_ident);
     }
   }
   std::map<art::Ptr<recob::Hit>, std::vector<art::Ptr<recob::Hit>>> ret;
   for (auto const& planeHits : hits_org) {
     for (const OrganizedHits& hit_org : planeHits) {
       std::vector<art::Ptr<recob::Hit>> secondaryHits{};
       for (const unsigned secondaryID : hit_org.second) {
         secondaryHits.push_back(hits[secondaryID]);
       }
       ret[hits[hit_org.first]] = std::move(secondaryHits);
     }
   }
   return ret;
 }

double shower::LArPandoraShowerAlg::RMSShowerGradient	(	std::vector< art::Ptr< recob::SpacePoint >> &	sps,
		const geo::Point_t &	ShowerCentre,
		const geo::Vector_t &	Direction,
		const unsigned int	nSegments
	)		const

Definition at line 399 of file LArPandoraShowerAlg.cxx.

References util::abs(), CalculateRMS(), util::counter(), OrderShowerSpacePoints(), ShowerCentre(), SpacePointPerpendicular(), and SpacePointProjection().

Referenced by ShowerRecoTools::ShowerDirectionCheater::CalculateElement(), and ShowerRecoTools::ShowerPCADirection::CalculateElement().

 {
   if (nSegments == 0)
     throw cet::exception("LArPandoraShowerAlg")
       << "Unable to calculate RMS Shower Gradient with 0 segments" << std::endl;
 
   if (sps.size() < 3) return 0;
 
   //Order the spacepoints
   this->OrderShowerSpacePoints(sps, ShowerCentre, Direction);
 
   //Get the length of the shower.
   const double minProj = this->SpacePointProjection(sps[0], ShowerCentre, Direction);
   const double maxProj = this->SpacePointProjection(sps[sps.size() - 1], ShowerCentre, Direction);
 
   const double length = (maxProj - minProj);
   const double segmentsize = length / nSegments;
 
   if (segmentsize < std::numeric_limits<double>::epsilon()) return 0;
 
   std::map<int, std::vector<float>> len_segment_map;
 
   //Split the the spacepoints into segments.
   for (auto const& sp : sps) {
 
     //Get the the projected length
     const double len = this->SpacePointProjection(sp, ShowerCentre, Direction);
 
     //Get the length to the projection
     const double len_perp = this->SpacePointPerpendicular(sp, ShowerCentre, Direction, len);
 
     int sg_len = std::round(len / segmentsize);
     len_segment_map[sg_len].push_back(len_perp);
   }
 
   int counter = 0;
   float sumx = 0.f;
   float sumy = 0.f;
   float sumx2 = 0.f;
   float sumxy = 0.f;
 
   //Get the rms of the segments and caclulate the gradient.
   for (auto const& segment : len_segment_map) {
     if (segment.second.size() < 2) continue;
     float RMS = this->CalculateRMS(segment.second);
 
     // Check if the calculation failed
     if (RMS < 0) continue;
 
     //Calculate the gradient using regression
     sumx += segment.first;
     sumy += RMS;
     sumx2 += segment.first * segment.first;
     sumxy += RMS * segment.first;
     ++counter;
   }
 
   const float denom = counter * sumx2 - sumx * sumx;
 
   return std::abs(denom) < std::numeric_limits<float>::epsilon() ?
            0 :
            (counter * sumxy - sumx * sumy) / denom;
 }

double shower::LArPandoraShowerAlg::SCECorrectEField	(	double const &	EField,
		geo::Point_t const &	pos
	)		const

Definition at line 508 of file LArPandoraShowerAlg.cxx.

References spacecharge::SpaceCharge::EnableSimEfieldSCE(), fSCE, and spacecharge::SpaceCharge::GetEfieldOffsets().

Referenced by ShowerRecoTools::ShowerTrajPointdEdx::CalculateElement().

 {
 
   // Check the space charge service is properly configured
   if (!fSCE || !fSCE->EnableSimEfieldSCE()) {
     throw cet::exception("LArPandoraShowerALG")
       << "Trying to correct SCE EField when service is not configured" << std::endl;
   }
   // Gets relative E field Distortions
   geo::Vector_t EFieldOffsets = fSCE->GetEfieldOffsets(pos);
   // Add 1 in X direction as this is the direction of the drift field
   EFieldOffsets += geo::Vector_t{1, 0, 0};
   // Convert to Absolute E Field from relative
   EFieldOffsets *= EField;
   // We only care about the magnitude for recombination
   return EFieldOffsets.r();
 }

double shower::LArPandoraShowerAlg::SCECorrectPitch	(	double const &	pitch,
		geo::Point_t const &	pos,
		geo::Vector_t const &	dir,
		unsigned int const &	TPC
	)		const

Definition at line 479 of file LArPandoraShowerAlg.cxx.

References dir, spacecharge::SpaceCharge::EnableCalSpatialSCE(), fSCE, fSCEXFlip, spacecharge::SpaceCharge::GetCalPosOffsets(), and spacecharge::SpaceCharge::GetPosOffsets().

Referenced by ShowerRecoTools::ShowerTrajPointdEdx::CalculateElement().

 {
 
   if (!fSCE || !fSCE->EnableCalSpatialSCE()) {
     throw cet::exception("LArPandoraShowerALG")
       << "Trying to correct SCE pitch when service is not configured" << std::endl;
   }
   // As the input pos is sce corrected already, find uncorrected pos
   const geo::Point_t uncorrectedPos = pos + fSCE->GetPosOffsets(pos);
   //Get the size of the correction at pos
   const geo::Vector_t posOffset = fSCE->GetCalPosOffsets(uncorrectedPos, TPC);
 
   //Get the position of next hit
   const geo::Point_t nextPos = uncorrectedPos + pitch * dir;
   //Get the offsets at the next pos
   const geo::Vector_t nextPosOffset = fSCE->GetCalPosOffsets(nextPos, TPC);
 
   //Calculate the corrected pitch
   const int xFlip(fSCEXFlip ? -1 : 1);
   geo::Vector_t pitchVec{pitch * dir.X() + xFlip * (nextPosOffset.X() - posOffset.X()),
                          pitch * dir.Y() + (nextPosOffset.Y() - posOffset.Y()),
                          pitch * dir.Z() + (nextPosOffset.Z() - posOffset.Z())};
 
   return pitchVec.r();
 }

geo::Point_t shower::LArPandoraShowerAlg::ShowerCentre ( std::vector< art::Ptr< recob::SpacePoint >> const & showersps ) const

Definition at line 190 of file LArPandoraShowerAlg.cxx.

References x, y, and z.

Referenced by ShowerRecoTools::ShowerPCAPropergationStartPosition::CalculateElement(), ShowerRecoTools::ShowerPFPVertexStartPosition::CalculateElement(), ShowerRecoTools::ShowerDirectionCheater::CalculateElement(), ShowerRecoTools::ShowerPCADirection::CalculateShowerPCA(), ShowerRecoTools::ShowerIncrementalTrackHitFinder::FitSegmentAndCalculateResidual(), RMSShowerGradient(), ShowerCentre(), and ShowerRecoTools::ShowerTrackPCADirection::ShowerPCAVector().

 {
   if (showersps.empty()) return {};
 
   double x{};
   double y{};
   double z{};
   for (auto const& sp : showersps) {
     auto const pos = sp->position();
     x += pos.X();
     y += pos.Y();
     z += pos.Z();
   }
   geo::Point_t centre_position{x, y, z};
   centre_position *= (1. / showersps.size());
   return centre_position;
 }

geo::Point_t shower::LArPandoraShowerAlg::ShowerCentre	(	detinfo::DetectorClocksData const &	clockData,
		detinfo::DetectorPropertiesData const &	detProp,
		std::vector< art::Ptr< recob::SpacePoint >> const &	showersps,
		art::FindManyP< recob::Hit > const &	fmh,
		float &	totalCharge
	)		const

Definition at line 221 of file LArPandoraShowerAlg.cxx.

References detinfo::DetectorPropertiesData::ElectronLifetime(), fUseCollectionOnly, hits(), geo::kCollection, pmtana::mean(), n, detinfo::sampling_rate(), x, y, and z.

 {
   geo::Point_t chargePoint{};
 
   //Loop over the spacepoints and get the charge weighted center.
   for (auto const& sp : showersps) {
 
     //Get the associated hits
     std::vector<art::Ptr<recob::Hit>> const& hits = fmh.at(sp.key());
 
     //Average the charge unless sepcified.
     float charge = 0;
     float charge2 = 0;
     for (auto const& hit : hits) {
 
       if (fUseCollectionOnly) {
         if (hit->SignalType() == geo::kCollection) {
           charge = hit->Integral();
           //Correct for the lifetime: Need to do other detproperites
           charge *= std::exp((sampling_rate(clockData) * hit->PeakTime()) /
                              (detProp.ElectronLifetime() * 1e3));
           break;
         }
       }
       else {
 
         //Correct for the lifetime FIX: Need  to do other detproperties somehow
         double Q = hit->Integral() * std::exp((sampling_rate(clockData) * hit->PeakTime()) /
                                               (detProp.ElectronLifetime() * 1e3));
 
         charge += Q;
         charge2 += Q * Q;
       }
     }
 
     if (!fUseCollectionOnly) {
       //Calculate the unbiased standard deviation and mean.
       float mean = charge / ((float)hits.size());
 
       float rms = 1;
 
       if (hits.size() > 1) {
         rms = std::sqrt((charge2 - charge * charge) / ((float)(hits.size() - 1)));
       }
 
       charge = 0;
       int n = 0;
       for (auto const& hit : hits) {
         double lifetimecorrection = std::exp((sampling_rate(clockData) * hit->PeakTime()) /
                                              (detProp.ElectronLifetime() * 1e3));
         if (hit->Integral() * lifetimecorrection > (mean - 2 * rms) &&
             hit->Integral() * lifetimecorrection < (mean + 2 * rms)) {
           charge += hit->Integral() * lifetimecorrection;
           ++n;
         }
       }
 
       if (n == 0) {
         mf::LogWarning("LArPandoraShowerAlg") << "no points used to make the charge value. \n";
       }
 
       charge /= n;
     }
 
     //Get the position of the spacepoint
     auto const pos = sp->position();
     double const x = pos.X();
     double const y = pos.Y();
     double const z = pos.Z();
 
     chargePoint.SetXYZ(
       chargePoint.X() + charge * x, chargePoint.Y() + charge * y, chargePoint.Z() + charge * z);
     totalCharge += charge;
 
     if (charge == 0) {
       mf::LogWarning("LArPandoraShowerAlg") << "Averaged charge, within 2 sigma, for a spacepoint "
                                                "is zero, Maybe this not a good method. \n";
     }
   }
 
   double intotalcharge = 1 / totalCharge;
   return chargePoint * intotalcharge;
 }

geo::Point_t shower::LArPandoraShowerAlg::ShowerCentre	(	detinfo::DetectorClocksData const &	clockData,
		detinfo::DetectorPropertiesData const &	detProp,
		std::vector< art::Ptr< recob::SpacePoint >> const &	showerspcs,
		art::FindManyP< recob::Hit > const &	fmh
	)		const

Definition at line 209 of file LArPandoraShowerAlg.cxx.

References ShowerCentre().

 {
   float totalCharge = 0;
   return shower::LArPandoraShowerAlg::ShowerCentre(
     clockData, detProp, showerspcs, fmh, totalCharge);
 }

double shower::LArPandoraShowerAlg::SpacePointCharge	(	art::Ptr< recob::SpacePoint > const &	sp,
		art::FindManyP< recob::Hit > const &	fmh
	)		const

Definition at line 318 of file LArPandoraShowerAlg.cxx.

References hits(), and art::Ptr< T >::key().

Referenced by ShowerRecoTools::ShowerPCADirection::CalculateShowerPCA(), ShowerRecoTools::ShowerTrackPCADirection::ShowerPCAVector(), and ShowerRecoTools::ShowerIncrementalTrackHitFinder::ShowerPCAVector().

 {
   double Charge = 0;
 
   //Average over the charge even though there is only one
   std::vector<art::Ptr<recob::Hit>> const& hits = fmh.at(sp.key());
   for (auto const& hit : hits) {
     Charge += hit->Integral();
   }
 
   Charge /= (float)hits.size();
 
   return Charge;
 }

double shower::LArPandoraShowerAlg::SpacePointPerpendicular	(	art::Ptr< recob::SpacePoint > const &	sp,
		geo::Point_t const &	vertex,
		geo::Vector_t const &	direction
	)		const

Definition at line 373 of file LArPandoraShowerAlg.cxx.

References proj, and SpacePointProjection().

Referenced by ShowerRecoTools::ShowerLengthPercentile::CalculateElement(), ShowerRecoTools::Shower3DCylinderTrackHitFinder::FindTrackSpacePoints(), OrderShowerSpacePointsPerpendicular(), and RMSShowerGradient().

 {
   // Get the projection of the spacepoint
   double proj = shower::LArPandoraShowerAlg::SpacePointProjection(sp, vertex, direction);
 
   return shower::LArPandoraShowerAlg::SpacePointPerpendicular(sp, vertex, direction, proj);
 }

double shower::LArPandoraShowerAlg::SpacePointPerpendicular	(	art::Ptr< recob::SpacePoint > const &	sp,
		geo::Point_t const &	vertex,
		geo::Vector_t const &	direction,
		double	proj
	)		const

Definition at line 383 of file LArPandoraShowerAlg.cxx.

References recob::SpacePoint::position().

 {
   // Get the position of the spacepoint
   auto pos = sp->position() - vertex;
 
   // Take away the projection * distance to find the perpendicular vector
   pos = pos - proj * direction;
 
   // Get the the projected length
   return pos.R();
 }

double shower::LArPandoraShowerAlg::SpacePointProjection	(	art::Ptr< recob::SpacePoint > const &	sp,
		geo::Point_t const &	vertex,
		geo::Vector_t const &	direction
	)		const

Definition at line 362 of file LArPandoraShowerAlg.cxx.

References recob::SpacePoint::position().

Referenced by ShowerRecoTools::ShowerLengthPercentile::CalculateElement(), ShowerRecoTools::ShowerIncrementalTrackHitFinder::CalculateElement(), shower::LArPandoraShowerCheatingAlg::CheatDebugEVD(), DebugEVD(), ShowerRecoTools::Shower3DCylinderTrackHitFinder::FindTrackSpacePoints(), OrderShowerSpacePoints(), RMSShowerGradient(), and SpacePointPerpendicular().

 {
   // Get the position of the spacepoint
   auto const pos = sp->position() - vertex;
 
   // Get the the projected length
   return pos.Dot(direction);
 }

double shower::LArPandoraShowerAlg::SpacePointTime	(	art::Ptr< recob::SpacePoint > const &	sp,
		art::FindManyP< recob::Hit > const &	fmh
	)		const

Definition at line 335 of file LArPandoraShowerAlg.cxx.

References hits(), and art::Ptr< T >::key().

Referenced by ShowerRecoTools::ShowerPCADirection::CalculateShowerPCA(), ShowerRecoTools::ShowerTrackPCADirection::ShowerPCAVector(), and ShowerRecoTools::ShowerIncrementalTrackHitFinder::ShowerPCAVector().

 {
   double Time = 0;
 
   //Average over the hits
   std::vector<art::Ptr<recob::Hit>> const& hits = fmh.at(sp.key());
   for (auto const& hit : hits) {
     Time += hit->PeakTime();
   }
 
   Time /= (float)hits.size();
   return Time;
 }

Member Data Documentation

geo::WireReadoutGeom const* shower::LArPandoraShowerAlg::fChannelMap = &art::ServiceHandle<geo::WireReadout>()->Get()

private

Definition at line 127 of file LArPandoraShowerAlg.h.

Referenced by HitCoordinates(), and OrderShowerHits().

const std::string shower::LArPandoraShowerAlg::fInitialTrackInputLabel

private

Definition at line 130 of file LArPandoraShowerAlg.h.

Referenced by DebugEVD().

const std::string shower::LArPandoraShowerAlg::fInitialTrackSpacePointsInputLabel

private

Definition at line 133 of file LArPandoraShowerAlg.h.

Referenced by DebugEVD().

art::InputTag shower::LArPandoraShowerAlg::fPFParticleLabel

private

Definition at line 123 of file LArPandoraShowerAlg.h.

Referenced by DebugEVD().

spacecharge::SpaceCharge const* shower::LArPandoraShowerAlg::fSCE

private

Definition at line 126 of file LArPandoraShowerAlg.h.

Referenced by SCECorrectEField(), and SCECorrectPitch().

bool shower::LArPandoraShowerAlg::fSCEXFlip

private

Definition at line 124 of file LArPandoraShowerAlg.h.

Referenced by SCECorrectPitch().

const std::string shower::LArPandoraShowerAlg::fShowerDirectionInputLabel

private

Definition at line 132 of file LArPandoraShowerAlg.h.

Referenced by DebugEVD().

const std::string shower::LArPandoraShowerAlg::fShowerStartPositionInputLabel

private

Definition at line 131 of file LArPandoraShowerAlg.h.

Referenced by DebugEVD().

bool shower::LArPandoraShowerAlg::fUseCollectionOnly

private

Definition at line 122 of file LArPandoraShowerAlg.h.

Referenced by ShowerCentre().

art::ServiceHandle<art::TFileService> shower::LArPandoraShowerAlg::tfs

private

Definition at line 128 of file LArPandoraShowerAlg.h.

Referenced by DebugEVD().

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

larpandora/v10_00_14/source/larpandora/LArPandoraEventBuilding/LArPandoraShower/Algs/LArPandoraShowerAlg.h
larpandora/v10_00_14/source/larpandora/LArPandoraEventBuilding/LArPandoraShower/Algs/LArPandoraShowerAlg.cxx

Public Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation