#include "TCShowerAlg.h"

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

int	makeShowers (detinfo::DetectorClocksData const &dataClock, detinfo::DetectorPropertiesData const &detProp, std::vector< art::Ptr< recob::PFParticle >> const &pfplist, std::vector< art::Ptr< recob::Cluster >> const &clusterlist, std::vector< art::Ptr< recob::Hit >> const &hitlist, art::FindManyP< recob::Hit > const &cls_fm, art::FindManyP< recob::Cluster > const &clspfp_fm, art::FindManyP< recob::Vertex > const &vtxpfp_fm, art::FindManyP< recob::Cluster > const &hitcls_fm, art::FindManyP< recob::Track > const &trkpfp_fm)

Public Attributes
geo::Vector_t	shwDir

TVector3	dcosVtxErr

geo::Point_t	shwvtx

TVector3	xyzErr

std::vector< double >	totalEnergy

std::vector< double >	totalEnergyErr

std::vector< double >	dEdx

std::vector< double >	dEdxErr

int	bestplane

std::vector< art::Ptr< recob::Hit > >	showerHits

Private Member Functions
int	goodHit (detinfo::DetectorClocksData const &dataClock, detinfo::DetectorPropertiesData const &detProp, art::Ptr< recob::Hit > const &, double maxDist, double minDistVert, std::map< geo::PlaneID, double > const &trk_wire1, std::map< geo::PlaneID, double > const &trk_tick1, std::map< geo::PlaneID, double > const &trk_wire2, std::map< geo::PlaneID, double > const &trk_tick2) const

int	goodHit (detinfo::DetectorClocksData const &dataClock, detinfo::DetectorPropertiesData const &detProp, art::Ptr< recob::Hit > const &, double maxDist, double minDistVert, std::map< geo::PlaneID, double > const &trk_wire1, std::map< geo::PlaneID, double > const &trk_tick1, std::map< geo::PlaneID, double > const &trk_wire2, std::map< geo::PlaneID, double > const &trk_tick2, int &pull) const

bool	addShowerHit (art::Ptr< recob::Hit > hit, std::vector< art::Ptr< recob::Hit >> showerhits) const

Private Attributes
calo::CalorimetryAlg	fCalorimetryAlg

pma::ProjectionMatchingAlg	fProjectionMatchingAlg

Detailed Description

Definition at line 38 of file TCShowerAlg.h.

Constructor & Destructor Documentation

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

Definition at line 44 of file TCShowerAlg.cxx.

     : fCalorimetryAlg(pset.get<fhicl::ParameterSet>("CalorimetryAlg"))
     , fProjectionMatchingAlg(pset.get<fhicl::ParameterSet>("ProjectionMatchingAlg"))
   {}

Member Function Documentation

bool shower::TCShowerAlg::addShowerHit	(	art::Ptr< recob::Hit >	hit,
		std::vector< art::Ptr< recob::Hit >>	showerhits
	)		const

private

Definition at line 446 of file TCShowerAlg.cxx.

References art::Ptr< T >::key().

Referenced by makeShowers().

   {
 
     for (size_t i = 0; i < showerhits.size(); ++i) {
       if (hit.key() == showerhits[i].key()) return false;
     }
 
     return true;
 
   } // addShowerHit

int shower::TCShowerAlg::goodHit	(	detinfo::DetectorClocksData const &	dataClock,
		detinfo::DetectorPropertiesData const &	detProp,
		art::Ptr< recob::Hit > const &	hit,
		double	maxDist,
		double	minDistVert,
		std::map< geo::PlaneID, double > const &	trk_wire1,
		std::map< geo::PlaneID, double > const &	trk_tick1,
		std::map< geo::PlaneID, double > const &	trk_wire2,
		std::map< geo::PlaneID, double > const &	trk_tick2
	)		const

private

Definition at line 356 of file TCShowerAlg.cxx.

Referenced by makeShowers().

   {
     int pull = 0;
     return goodHit(clockData,
                    detProp,
                    hit,
                    maxDist,
                    minDistVert,
                    trk_wire1,
                    trk_tick1,
                    trk_wire2,
                    trk_tick2,
                    pull);
 
   } // goodHit

int shower::TCShowerAlg::goodHit	(	detinfo::DetectorClocksData const &	dataClock,
		detinfo::DetectorPropertiesData const &	detProp,
		art::Ptr< recob::Hit > const &	hit,
		double	maxDist,
		double	minDistVert,
		std::map< geo::PlaneID, double > const &	trk_wire1,
		std::map< geo::PlaneID, double > const &	trk_tick1,
		std::map< geo::PlaneID, double > const &	trk_wire2,
		std::map< geo::PlaneID, double > const &	trk_tick2,
		int &	pull
	)		const

private

Definition at line 385 of file TCShowerAlg.cxx.

References util::abs(), larg4::dist(), detinfo::DetectorPropertiesData::DriftVelocity(), detinfo::DetectorPropertiesData::Efield(), recob::Hit::PeakTime(), detinfo::sampling_rate(), detinfo::DetectorPropertiesData::Temperature(), geo::WireID::Wire, recob::Hit::WireID(), geo::GeometryCore::WirePitch(), x1, x2, y1, and y2.

   {
     art::ServiceHandle<geo::Geometry const> geom;
 
     double wirePitch = geom->WirePitch(hit->WireID());
     double tickToDist = detProp.DriftVelocity(detProp.Efield(), detProp.Temperature());
     tickToDist *= 1.e-3 * sampling_rate(clockData); // 1e-3 is conversion of 1/us to 1/ns
     double UnitsPerTick = tickToDist / wirePitch;
 
     double x0 = hit->WireID().Wire;
     double y0 = hit->PeakTime() * UnitsPerTick;
 
     double x1 = trk_wire1.at(hit->WireID());
     double y1 = trk_tick1.at(hit->WireID()) * UnitsPerTick;
 
     double x2 = trk_wire2.at(hit->WireID());
     double y2 = trk_tick2.at(hit->WireID()) * UnitsPerTick;
 
     double distToVert = std::hypot(x0 - x1, y0 - y1);
     if (distToVert < minDistVert) return -1;
 
     if (x2 > x1) {
       if (distToVert < 50)
         maxDist /= 4;
       else if (distToVert < 100)
         maxDist /= 2; // trying to exclude photon showers
     }
 
     // exclude cluster if it's "behind" the vertex
     double a = std::hypot(x2 - x1, y2 - y1);
     double b = std::hypot(x0 - x1, y0 - y1);
     double c = std::hypot(x0 - x2, y0 - y2);
     double costheta = -(pow(c, 2) - pow(a, 2) - pow(b, 2)) / (2 * a * b);
     if (costheta < 0) return -1;
 
     double dist =
       std::abs((y2 - y1) * x0 - (x2 - x1) * y0 + x2 * y1 - y2 * x1) / std::hypot(y2 - y1, x2 - x1);
 
     if (dist < maxDist) {
       if (((y2 - y1) * x0 - (x2 - x1) * y0 + x2 * y1 - y2 * x1) > 0)
         pull = 1;
       else
         pull = -1;
       return 1;
     }
 
     return 0;
 
   } // goodHit

int shower::TCShowerAlg::makeShowers	(	detinfo::DetectorClocksData const &	dataClock,
		detinfo::DetectorPropertiesData const &	detProp,
		std::vector< art::Ptr< recob::PFParticle >> const &	pfplist,
		std::vector< art::Ptr< recob::Cluster >> const &	clusterlist,
		std::vector< art::Ptr< recob::Hit >> const &	hitlist,
		art::FindManyP< recob::Hit > const &	cls_fm,
		art::FindManyP< recob::Cluster > const &	clspfp_fm,
		art::FindManyP< recob::Vertex > const &	vtxpfp_fm,
		art::FindManyP< recob::Cluster > const &	hitcls_fm,
		art::FindManyP< recob::Track > const &	trkpfp_fm
	)

Definition at line 49 of file TCShowerAlg.cxx.

References util::abs(), addShowerHit(), detinfo::DetectorPropertiesData::ConvertXToTicks(), dEdx, dEdxErr, recob::Track::End(), recob::Track::EndDirection(), goodHit(), geo::GeometryCore::Iterate(), recob::Vertex::position(), showerHits, shwDir, shwvtx, totalEnergy, totalEnergyErr, lar::dump::vector(), recob::Track::Vertex(), recob::Track::VertexDirection(), and geo::GeometryCore::WireCoordinate().

Referenced by shower::TCShower::getShowersWithoutSlices(), and shower::TCShower::getShowersWithSlices().

   {
     totalEnergy.resize(2);
     totalEnergyErr.resize(2);
     dEdx.resize(2);
     dEdxErr.resize(2);
 
     art::ServiceHandle<geo::Geometry const> geom;
 
     std::vector<pfpStuff> allpfps;
 
     // put together pfparticle information
     for (size_t i = 0; i < pfplist.size(); ++i) {
       pfpStuff thispfp;
       thispfp.hits.clear();
       thispfp.clsIDs.clear();
       thispfp.pfp = pfplist[i];
 
       std::vector<art::Ptr<recob::Vertex>> thisvtxlist = vtxpfp_fm.at(pfplist[i].key());
       if (thisvtxlist.size()) thispfp.vtx = thisvtxlist[0];
 
       std::vector<art::Ptr<recob::Track>> thistrklist = trkpfp_fm.at(pfplist[i].key());
       if (thistrklist.size()) thispfp.trk = thistrklist[0];
 
       std::vector<art::Ptr<recob::Cluster>> thisclusterlist = clspfp_fm.at(pfplist[i].key());
       std::vector<int> clustersize;
 
       for (size_t j = 0; j < thisclusterlist.size(); ++j) {
 
         thispfp.clsIDs.push_back(thisclusterlist[j]->ID());
 
         std::vector<art::Ptr<recob::Hit>> thishitlist = cls_fm.at(thisclusterlist[j].key());
         clustersize.push_back((int)thishitlist.size());
 
         for (size_t k = 0; k < thishitlist.size(); ++k) {
           thispfp.hits.push_back(thishitlist[k]);
         } // loop through hits
 
       } // loop through clusters
 
       if (clustersize.size() == 3) {
         if (!thispfp.vtx) continue;
         if (!thispfp.trk) continue;
 
         allpfps.push_back(thispfp);
 
         double tick = detProp.ConvertXToTicks(thispfp.vtx->position().X(), geo::PlaneID(0, 0, 2));
         int wire = geom->WireCoordinate(thispfp.vtx->position(), geo::PlaneID(0, 0, 2));
 
         std::cout << "pfp " << thispfp.pfp->Self() + 1 << " cluster sizes " << clustersize[0] << ":"
                   << clustersize[1] << ":" << clustersize[2] << " vertex " << thispfp.vtx->ID()
                   << " " << tick << ":" << wire << " z " << thispfp.vtx->position().Z()
                   << std::endl;
 
       } // add pfp to list
 
     } // loop through pfparticles
 
     std::sort(allpfps.begin(), allpfps.end(), comparePFP);
     std::reverse(allpfps.begin(), allpfps.end());
 
     std::cout << "sorted pfps: ";
     for (size_t i = 0; i < allpfps.size(); ++i)
       std::cout << allpfps[i].pfp->Self() + 1 << " ";
     std::cout << std::endl;
 
     bool showerCandidate = false;
 
     for (size_t i = 0; i < allpfps.size(); ++i) {
 
       showerHits.clear();
 
       art::Ptr<recob::Vertex> pfpvtx = allpfps[i].vtx;
       art::Ptr<recob::Track> pfptrk = allpfps[i].trk;
       std::vector<art::Ptr<recob::Hit>> pfphits = allpfps[i].hits;
       std::vector<art::Ptr<recob::Cluster>> pfpcls = clspfp_fm.at(allpfps[i].pfp.key());
 
       std::cout << "pfp " << allpfps[i].pfp->Self() + 1 << " hits " << pfphits.size() << std::endl;
 
       TVector3 vtx;
       vtx[0] = pfpvtx->position().X();
       vtx[1] = pfpvtx->position().Y();
       vtx[2] = pfpvtx->position().Z();
 
       if (pfptrk->Vertex().Z() < pfptrk->End().Z()) {
         shwvtx = pfptrk->Vertex();
         shwDir = pfptrk->VertexDirection();
       }
       else {
         shwvtx = pfptrk->End();
         shwDir = -pfptrk->EndDirection();
       }
 
       int tolerance = 60;         // how many shower like cluster you need to define a shower
       double pullTolerance = 0.7; // hits should be evenly distributed around the track
       double maxDist = 20;        // how far a shower like cluster can be from the track
       double minDistVert = 15;    // exclude tracks near the vertex
 
       // TODO: count number of shower-like hits hits "behind" vertex
       // if high, pick an earlier cluster and refit hits
       // this is going to take some restructuring
 
       if (pfphits.size() < 30) continue;
       if (pfphits.size() > 500) continue;
       // adjust tolerances for short tracks
       if (pfphits.size() < 90) {
         tolerance = 50;
         pullTolerance = 0.9;
       }
       if (pfphits.size() > 400)
         tolerance = 200;
       else if (pfphits.size() > 100)
         tolerance = 100; // RSF added used to be 100, 100
 
       // add pfp hits to shower
       for (size_t ii = 0; ii < pfphits.size(); ++ii) {
         if (addShowerHit(pfphits[ii], showerHits)) showerHits.push_back(pfphits[ii]);
       } // loop over pfphits
 
       int nShowerHits = 0;
       double showerHitPull = 0;
 
       geo::Point_t pfpStart = shwvtx;
       geo::Point_t pfpPt2 = shwvtx + shwDir; // a second point along the track
 
       // track vertex
       std::map<geo::PlaneID, double> trk_tick1;
       std::map<geo::PlaneID, double> trk_wire1;
 
       // second track point
       std::map<geo::PlaneID, double> trk_tick2;
       std::map<geo::PlaneID, double> trk_wire2;
 
       for (auto const& planeid : geom->Iterate<geo::PlaneID>()) {
         trk_tick1[planeid] = detProp.ConvertXToTicks(pfpStart.X(), planeid);
         trk_wire1[planeid] = geom->WireCoordinate(pfpStart, planeid);
         trk_tick2[planeid] = detProp.ConvertXToTicks(pfpPt2.X(), planeid);
         trk_wire2[planeid] = geom->WireCoordinate(pfpPt2, planeid);
       }
 
       for (size_t j = 0; j < clusterlist.size(); ++j) {
         std::vector<art::Ptr<recob::Hit>> cls_hitlist = cls_fm.at(clusterlist[j].key());
 
         if (clusterlist[j]->ID() > 0 && cls_hitlist.size() > 10) continue;
         if (cls_hitlist.size() > 50) continue;
 
         bool isGoodCluster = false; // true if the hit belongs to a cluster that
                                     // should be added to the shower
 
         bool skipit = false; // skip clusters already in the pfp
         for (size_t k = 0; k < allpfps[i].clsIDs.size(); ++k) {
           if (allpfps[i].clsIDs[k] == clusterlist[j]->ID()) skipit = true;
         }
         if (skipit) continue;
 
         for (size_t jj = 0; jj < cls_hitlist.size(); ++jj) {
           int isGoodHit = goodHit(clockData,
                                   detProp,
                                   cls_hitlist[jj],
                                   maxDist,
                                   minDistVert,
                                   trk_wire1,
                                   trk_tick1,
                                   trk_wire2,
                                   trk_tick2);
 
           if (isGoodHit == -1) {
             isGoodCluster = false;
             break;
           }
           else if (isGoodHit == 1) {
             isGoodCluster = true;
           }
 
         } // loop over hits in cluster
 
         // add hits to shower
         if (isGoodCluster) {
           for (size_t jj = 0; jj < cls_hitlist.size(); ++jj) {
             nShowerHits++;
 
             int showerHitPullAdd = 0;
             goodHit(clockData,
                     detProp,
                     cls_hitlist[jj],
                     maxDist,
                     minDistVert,
                     trk_wire1,
                     trk_tick1,
                     trk_wire2,
                     trk_tick2,
                     showerHitPullAdd);
             showerHitPull += showerHitPullAdd;
 
             if (addShowerHit(cls_hitlist[jj], showerHits)) showerHits.push_back(cls_hitlist[jj]);
 
           } // loop over hits in cluster
         }   // cluster contains hit close to track
 
       } // loop through cluserlist
 
       showerHitPull /= nShowerHits;
 
       std::cout << "shower hits " << showerHits.size() << " " << nShowerHits << " shower pull "
                 << showerHitPull << std::endl;
 
       if (nShowerHits > tolerance && std::abs(showerHitPull) < pullTolerance) {
         showerCandidate = true;
         std::cout << "SHOWER CANDIDATE" << std::endl;
         // loop over hits to find those that aren't associated with any clusters
         if (nShowerHits > 400) maxDist *= 2; // TODO: optimize this threshold
         for (size_t k = 0; k < hitlist.size(); ++k) {
           std::vector<art::Ptr<recob::Cluster>> hit_clslist = hitcls_fm.at(hitlist[k].key());
           if (hit_clslist.size()) continue;
 
           int isGoodHit = goodHit(clockData,
                                   detProp,
                                   hitlist[k],
                                   maxDist * 2,
                                   minDistVert * 2,
                                   trk_wire1,
                                   trk_tick1,
                                   trk_wire2,
                                   trk_tick2);
           if (isGoodHit == 1 && addShowerHit(hitlist[k], showerHits))
             showerHits.push_back(hitlist[k]);
         } // loop over hits
 
         // loop over clusters to see if any fall within the shower
         for (size_t k = 0; k < clusterlist.size(); ++k) {
           std::vector<art::Ptr<recob::Hit>> cls_hitlist = cls_fm.at(clusterlist[k].key());
           if (clusterlist[k]->ID() > 0 && cls_hitlist.size() > 50) continue;
 
           double thisDist = maxDist;
           double thisMin = minDistVert;
 
           if (cls_hitlist.size() < 10) {
             thisDist *= 4;
             thisMin *= 4;
           }
           else if (cls_hitlist.size() < 30)
             thisDist *= 2;
 
           int nhits = 0;
           int ngoodhits = 0;
 
           // are the cluster hits close?
           for (size_t kk = 0; kk < cls_hitlist.size(); ++kk) {
             nhits++;
             int isGoodHit = goodHit(clockData,
                                     detProp,
                                     cls_hitlist[kk],
                                     thisDist,
                                     thisMin,
                                     trk_wire1,
                                     trk_tick1,
                                     trk_wire2,
                                     trk_tick2);
             if (isGoodHit == -1) {
               ngoodhits = 0;
               break;
             }
             else if (isGoodHit == 1) {
               ngoodhits++;
             }
           } // loop over cluster hits
 
           double fracGood = (double)ngoodhits / nhits;
 
           bool isGoodTrack = fracGood > 0.4;
 
           if (isGoodTrack) {
             for (size_t kk = 0; kk < cls_hitlist.size(); ++kk) {
               if (addShowerHit(cls_hitlist[kk], showerHits)) showerHits.push_back(cls_hitlist[kk]);
             } // loop over hits to add them to showe
           }
 
         } // loop through clusterlist
 
       } // decide if shower
 
       if (showerCandidate) {
         std::cout << "THIS IS THE SHOWER PFP: " << allpfps[i].pfp->Self() + 1 << std::endl;
         break;
       }
 
     } // loop through allpfps
 
     if (showerCandidate) return 1;
 
     return 0;
   } // makeShowers