Public Member Functions | Public Attributes | Private Member Functions | Private Attributes | List of all members
shower::TCShowerAlg Class Reference

#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::Vertex >> const &vertexlist, 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::PFParticle > const &hit_fm, art::FindManyP< recob::Cluster > const &hitcls_fm, art::FindManyP< recob::Track > const &trkpfp_fm, art::FindManyP< anab::Calorimetry > const &fmcal)
 

Public Attributes

TVector3 shwDir
 
TVector3 dcosVtxErr
 
TVector3 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 35 of file TCShowerAlg.h.

Constructor & Destructor Documentation

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

Definition at line 45 of file TCShowerAlg.cxx.

46  : fCalorimetryAlg(pset.get<fhicl::ParameterSet>("CalorimetryAlg"))
47  , fProjectionMatchingAlg(pset.get<fhicl::ParameterSet>("ProjectionMatchingAlg"))
48  {}
pma::ProjectionMatchingAlg fProjectionMatchingAlg
Definition: TCShowerAlg.h:67
calo::CalorimetryAlg fCalorimetryAlg
Definition: TCShowerAlg.h:66

Member Function Documentation

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

Definition at line 455 of file TCShowerAlg.cxx.

457  {
458 
459  for (size_t i = 0; i < showerhits.size(); ++i) {
460  if (hit.key() == showerhits[i].key()) return false;
461  }
462 
463  return true;
464 
465  } // addShowerHit
key_type key() const noexcept
Definition: Ptr.h:216
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 363 of file TCShowerAlg.cxx.

372  {
373  int pull = 0;
374  return goodHit(clockData,
375  detProp,
376  hit,
377  maxDist,
378  minDistVert,
379  trk_wire1,
380  trk_tick1,
381  trk_wire2,
382  trk_tick2,
383  pull);
384 
385  } // goodHit
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 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 393 of file TCShowerAlg.cxx.

403  {
405 
406  double wirePitch = geom->WirePitch(hit->WireID());
407  double tickToDist = detProp.DriftVelocity(detProp.Efield(), detProp.Temperature());
408  tickToDist *= 1.e-3 * sampling_rate(clockData); // 1e-3 is conversion of 1/us to 1/ns
409  double UnitsPerTick = tickToDist / wirePitch;
410 
411  double x0 = hit->WireID().Wire;
412  double y0 = hit->PeakTime() * UnitsPerTick;
413 
414  double x1 = trk_wire1.at(hit->WireID());
415  double y1 = trk_tick1.at(hit->WireID()) * UnitsPerTick;
416 
417  double x2 = trk_wire2.at(hit->WireID());
418  double y2 = trk_tick2.at(hit->WireID()) * UnitsPerTick;
419 
420  double distToVert = std::hypot(x0 - x1, y0 - y1);
421  if (distToVert < minDistVert) return -1;
422 
423  if (x2 > x1) {
424  if (distToVert < 50)
425  maxDist /= 4;
426  else if (distToVert < 100)
427  maxDist /= 2; // trying to exclude photon showers
428  }
429 
430  // exclude cluster if it's "behind" the vertex
431  double a = std::hypot(x2 - x1, y2 - y1);
432  double b = std::hypot(x0 - x1, y0 - y1);
433  double c = std::hypot(x0 - x2, y0 - y2);
434  double costheta = -(pow(c, 2) - pow(a, 2) - pow(b, 2)) / (2 * a * b);
435  if (costheta < 0) return -1;
436 
437  double dist =
438  std::abs((y2 - y1) * x0 - (x2 - x1) * y0 + x2 * y1 - y2 * x1) / std::hypot(y2 - y1, x2 - x1);
439 
440  if (dist < maxDist) {
441  if (((y2 - y1) * x0 - (x2 - x1) * y0 + x2 * y1 - y2 * x1) > 0)
442  pull = 1;
443  else
444  pull = -1;
445  return 1;
446  }
447 
448  return 0;
449 
450  } // goodHit
geo::WireID WireID() const
Definition: Hit.h:233
constexpr T pow(T x)
Definition: pow.h:72
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:580
geo::Length_t WirePitch(geo::PlaneID const &planeid) const
Returns the distance between two consecutive wires.
std::enable_if_t< std::is_arithmetic_v< T >, T > hypot(T x, T y)
Definition: hypot.h:60
T abs(T value)
const double a
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:218
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
static bool * b
Definition: config.cpp:1043
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
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::Vertex >> const &  vertexlist,
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::PFParticle > const &  hit_fm,
art::FindManyP< recob::Cluster > const &  hitcls_fm,
art::FindManyP< recob::Track > const &  trkpfp_fm,
art::FindManyP< anab::Calorimetry > const &  fmcal 
)

Definition at line 51 of file TCShowerAlg.cxx.

64  {
65  totalEnergy.resize(2);
66  totalEnergyErr.resize(2);
67  dEdx.resize(2);
68  dEdxErr.resize(2);
69 
71 
72  std::vector<pfpStuff> allpfps;
73 
74  // put together pfparticle information
75  for (size_t i = 0; i < pfplist.size(); ++i) {
76  pfpStuff thispfp;
77  thispfp.hits.clear();
78  thispfp.clsIDs.clear();
79  thispfp.pfp = pfplist[i];
80 
81  std::vector<art::Ptr<recob::Vertex>> thisvtxlist = vtxpfp_fm.at(pfplist[i].key());
82  if (thisvtxlist.size()) thispfp.vtx = thisvtxlist[0];
83 
84  std::vector<art::Ptr<recob::Track>> thistrklist = trkpfp_fm.at(pfplist[i].key());
85  if (thistrklist.size()) thispfp.trk = thistrklist[0];
86 
87  std::vector<art::Ptr<recob::Cluster>> thisclusterlist = clspfp_fm.at(pfplist[i].key());
88  std::vector<int> clustersize;
89 
90  for (size_t j = 0; j < thisclusterlist.size(); ++j) {
91 
92  thispfp.clsIDs.push_back(thisclusterlist[j]->ID());
93 
94  std::vector<art::Ptr<recob::Hit>> thishitlist = cls_fm.at(thisclusterlist[j].key());
95  clustersize.push_back((int)thishitlist.size());
96 
97  for (size_t k = 0; k < thishitlist.size(); ++k) {
98  thispfp.hits.push_back(thishitlist[k]);
99  } // loop through hits
100 
101  } // loop through clusters
102 
103  if (clustersize.size() == 3) {
104  if (!thispfp.vtx) continue;
105  if (!thispfp.trk) continue;
106 
107  allpfps.push_back(thispfp);
108 
109  double tick = detProp.ConvertXToTicks(thispfp.vtx->position().X(), geo::PlaneID(0, 0, 2));
110  int wire = geom->WireCoordinate(
111  thispfp.vtx->position().Y(), thispfp.vtx->position().Z(), geo::PlaneID(0, 0, 2));
112 
113  std::cout << "pfp " << thispfp.pfp->Self() + 1 << " cluster sizes " << clustersize[0] << ":"
114  << clustersize[1] << ":" << clustersize[2] << " vertex " << thispfp.vtx->ID()
115  << " " << tick << ":" << wire << " z " << thispfp.vtx->position().Z()
116  << std::endl;
117 
118  } // add pfp to list
119 
120  } // loop through pfparticles
121 
122  std::sort(allpfps.begin(), allpfps.end(), comparePFP);
123  std::reverse(allpfps.begin(), allpfps.end());
124 
125  std::cout << "sorted pfps: ";
126  for (size_t i = 0; i < allpfps.size(); ++i)
127  std::cout << allpfps[i].pfp->Self() + 1 << " ";
128  std::cout << std::endl;
129 
130  bool showerCandidate = false;
131 
132  for (size_t i = 0; i < allpfps.size(); ++i) {
133 
134  showerHits.clear();
135 
136  art::Ptr<recob::Vertex> pfpvtx = allpfps[i].vtx;
137  art::Ptr<recob::Track> pfptrk = allpfps[i].trk;
138  std::vector<art::Ptr<recob::Hit>> pfphits = allpfps[i].hits;
139  std::vector<art::Ptr<recob::Cluster>> pfpcls = clspfp_fm.at(allpfps[i].pfp.key());
140 
141  std::cout << "pfp " << allpfps[i].pfp->Self() + 1 << " hits " << pfphits.size() << std::endl;
142 
143  TVector3 vtx;
144  vtx[0] = pfpvtx->position().X();
145  vtx[1] = pfpvtx->position().Y();
146  vtx[2] = pfpvtx->position().Z();
147 
148  if (pfptrk->Vertex().Z() < pfptrk->End().Z()) {
149  shwvtx = pfptrk->Vertex<TVector3>();
150  shwDir = pfptrk->VertexDirection<TVector3>();
151  }
152  else {
153  shwvtx = pfptrk->End<TVector3>();
154  shwDir = -pfptrk->EndDirection<TVector3>();
155  }
156 
157  int tolerance = 60; // how many shower like cluster you need to define a shower
158  double pullTolerance = 0.7; // hits should be evenly distributed around the track
159  double maxDist = 20; // how far a shower like cluster can be from the track
160  double minDistVert = 15; // exclude tracks near the vertex
161 
162  // TODO: count number of shower-like hits hits "behind" vertex
163  // if high, pick an earlier cluster and refit hits
164  // this is going to take some restructuring
165 
166  if (pfphits.size() < 30) continue;
167  if (pfphits.size() > 500) continue;
168  // adjust tolerances for short tracks
169  if (pfphits.size() < 90) {
170  tolerance = 50;
171  pullTolerance = 0.9;
172  }
173  if (pfphits.size() > 400)
174  tolerance = 200;
175  else if (pfphits.size() > 100)
176  tolerance = 100; // RSF added used to be 100, 100
177 
178  // add pfp hits to shower
179  for (size_t ii = 0; ii < pfphits.size(); ++ii) {
180  if (addShowerHit(pfphits[ii], showerHits)) showerHits.push_back(pfphits[ii]);
181  } // loop over pfphits
182 
183  int nShowerHits = 0;
184  double showerHitPull = 0;
185 
186  TVector3 pfpStart = shwvtx;
187  TVector3 pfpPt2 = shwvtx + shwDir; // a second point along the track
188 
189  // track vertex
190  std::map<geo::PlaneID, double> trk_tick1;
191  std::map<geo::PlaneID, double> trk_wire1;
192 
193  // second track point
194  std::map<geo::PlaneID, double> trk_tick2;
195  std::map<geo::PlaneID, double> trk_wire2;
196 
197  for (auto iPlane = geom->begin_plane_id(); iPlane != geom->end_plane_id(); ++iPlane) {
198  trk_tick1[*iPlane] = detProp.ConvertXToTicks(pfpStart[0], *iPlane);
199  trk_wire1[*iPlane] = geom->WireCoordinate(pfpStart[1], pfpStart[2], *iPlane);
200  trk_tick2[*iPlane] = detProp.ConvertXToTicks(pfpPt2[0], *iPlane);
201  trk_wire2[*iPlane] = geom->WireCoordinate(pfpPt2[1], pfpPt2[2], *iPlane);
202  }
203 
204  for (size_t j = 0; j < clusterlist.size(); ++j) {
205  std::vector<art::Ptr<recob::Hit>> cls_hitlist = cls_fm.at(clusterlist[j].key());
206 
207  if (clusterlist[j]->ID() > 0 && cls_hitlist.size() > 10) continue;
208  if (cls_hitlist.size() > 50) continue;
209 
210  bool isGoodCluster = false; // true if the hit belongs to a cluster that
211  // should be added to the shower
212 
213  bool skipit = false; // skip clusters already in the pfp
214  for (size_t k = 0; k < allpfps[i].clsIDs.size(); ++k) {
215  if (allpfps[i].clsIDs[k] == clusterlist[j]->ID()) skipit = true;
216  }
217  if (skipit) continue;
218 
219  for (size_t jj = 0; jj < cls_hitlist.size(); ++jj) {
220  int isGoodHit = goodHit(clockData,
221  detProp,
222  cls_hitlist[jj],
223  maxDist,
224  minDistVert,
225  trk_wire1,
226  trk_tick1,
227  trk_wire2,
228  trk_tick2);
229 
230  if (isGoodHit == -1) {
231  isGoodCluster = false;
232  break;
233  }
234  else if (isGoodHit == 1) {
235  isGoodCluster = true;
236  }
237 
238  } // loop over hits in cluster
239 
240  // add hits to shower
241  if (isGoodCluster) {
242  for (size_t jj = 0; jj < cls_hitlist.size(); ++jj) {
243  nShowerHits++;
244 
245  int showerHitPullAdd = 0;
246  goodHit(clockData,
247  detProp,
248  cls_hitlist[jj],
249  maxDist,
250  minDistVert,
251  trk_wire1,
252  trk_tick1,
253  trk_wire2,
254  trk_tick2,
255  showerHitPullAdd);
256  showerHitPull += showerHitPullAdd;
257 
258  if (addShowerHit(cls_hitlist[jj], showerHits)) showerHits.push_back(cls_hitlist[jj]);
259 
260  } // loop over hits in cluster
261  } // cluster contains hit close to track
262 
263  } // loop through cluserlist
264 
265  showerHitPull /= nShowerHits;
266 
267  std::cout << "shower hits " << showerHits.size() << " " << nShowerHits << " shower pull "
268  << showerHitPull << std::endl;
269 
270  if (nShowerHits > tolerance && std::abs(showerHitPull) < pullTolerance) {
271  showerCandidate = true;
272  std::cout << "SHOWER CANDIDATE" << std::endl;
273  // loop over hits to find those that aren't associated with any clusters
274  if (nShowerHits > 400) maxDist *= 2; // TODO: optimize this threshold
275  for (size_t k = 0; k < hitlist.size(); ++k) {
276  std::vector<art::Ptr<recob::Cluster>> hit_clslist = hitcls_fm.at(hitlist[k].key());
277  if (hit_clslist.size()) continue;
278 
279  int isGoodHit = goodHit(clockData,
280  detProp,
281  hitlist[k],
282  maxDist * 2,
283  minDistVert * 2,
284  trk_wire1,
285  trk_tick1,
286  trk_wire2,
287  trk_tick2);
288  if (isGoodHit == 1 && addShowerHit(hitlist[k], showerHits))
289  showerHits.push_back(hitlist[k]);
290  } // loop over hits
291 
292  // loop over clusters to see if any fall within the shower
293  for (size_t k = 0; k < clusterlist.size(); ++k) {
294  std::vector<art::Ptr<recob::Hit>> cls_hitlist = cls_fm.at(clusterlist[k].key());
295  if (clusterlist[k]->ID() > 0 && cls_hitlist.size() > 50) continue;
296 
297  double thisDist = maxDist;
298  double thisMin = minDistVert;
299 
300  if (cls_hitlist.size() < 10) {
301  thisDist *= 4;
302  thisMin *= 4;
303  }
304  else if (cls_hitlist.size() < 30)
305  thisDist *= 2;
306 
307  int nhits = 0;
308  int ngoodhits = 0;
309 
310  // are the cluster hits close?
311  for (size_t kk = 0; kk < cls_hitlist.size(); ++kk) {
312  nhits++;
313  int isGoodHit = goodHit(clockData,
314  detProp,
315  cls_hitlist[kk],
316  thisDist,
317  thisMin,
318  trk_wire1,
319  trk_tick1,
320  trk_wire2,
321  trk_tick2);
322  if (isGoodHit == -1) {
323  ngoodhits = 0;
324  break;
325  }
326  else if (isGoodHit == 1) {
327  ngoodhits++;
328  }
329  } // loop over cluster hits
330 
331  double fracGood = (double)ngoodhits / nhits;
332 
333  bool isGoodTrack = fracGood > 0.4;
334 
335  if (isGoodTrack) {
336  for (size_t kk = 0; kk < cls_hitlist.size(); ++kk) {
337  if (addShowerHit(cls_hitlist[kk], showerHits)) showerHits.push_back(cls_hitlist[kk]);
338  } // loop over hits to add them to showe
339  }
340 
341  } // loop through clusterlist
342 
343  } // decide if shower
344 
345  if (showerCandidate) {
346  std::cout << "THIS IS THE SHOWER PFP: " << allpfps[i].pfp->Self() + 1 << std::endl;
347  break;
348  }
349 
350  } // loop through allpfps
351 
352  if (showerCandidate) return 1;
353 
354  return 0;
355  } // makeShowers
geo::Length_t WireCoordinate(double YPos, double ZPos, geo::PlaneID const &planeid) const
Returns the index of the nearest wire to the specified position.
std::vector< double > totalEnergyErr
Definition: TCShowerAlg.h:43
bool addShowerHit(art::Ptr< recob::Hit > hit, std::vector< art::Ptr< recob::Hit >> showerhits) const
unsigned int ID
std::vector< double > dEdx
Definition: TCShowerAlg.h:44
auto const tolerance
struct vector vector
Vector_t VertexDirection() const
Definition: Track.h:132
T abs(T value)
def key(type, name=None)
Definition: graph.py:13
Point_t const & Vertex() const
Definition: Track.h:124
tick_as<> tick
Tick number, represented by std::ptrdiff_t.
Definition: electronics.h:75
plane_id_iterator end_plane_id() const
Returns an iterator pointing after the last plane ID in the detector.
static unsigned int reverse(QString &chars, unsigned char *level, unsigned int a, unsigned int b)
Definition: qstring.cpp:11649
Vector_t EndDirection() const
Definition: Track.h:133
std::vector< double > totalEnergy
Definition: TCShowerAlg.h:42
Point_t const & End() const
Definition: Track.h:125
detail::Node< FrameID, bool > PlaneID
Definition: CRTID.h:125
plane_id_iterator begin_plane_id() const
Returns an iterator pointing to the first plane ID in the detector.
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
std::vector< art::Ptr< recob::Hit > > showerHits
Definition: TCShowerAlg.h:47
std::vector< double > dEdxErr
Definition: TCShowerAlg.h:45
const Point_t & position() const
Return vertex 3D position.
Definition: Vertex.h:60
QTextStream & endl(QTextStream &s)

Member Data Documentation

int shower::TCShowerAlg::bestplane

Definition at line 46 of file TCShowerAlg.h.

TVector3 shower::TCShowerAlg::dcosVtxErr

Definition at line 39 of file TCShowerAlg.h.

std::vector<double> shower::TCShowerAlg::dEdx

Definition at line 44 of file TCShowerAlg.h.

std::vector<double> shower::TCShowerAlg::dEdxErr

Definition at line 45 of file TCShowerAlg.h.

calo::CalorimetryAlg shower::TCShowerAlg::fCalorimetryAlg
private

Definition at line 66 of file TCShowerAlg.h.

pma::ProjectionMatchingAlg shower::TCShowerAlg::fProjectionMatchingAlg
private

Definition at line 67 of file TCShowerAlg.h.

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

Definition at line 47 of file TCShowerAlg.h.

TVector3 shower::TCShowerAlg::shwDir

Definition at line 38 of file TCShowerAlg.h.

TVector3 shower::TCShowerAlg::shwvtx

Definition at line 40 of file TCShowerAlg.h.

std::vector<double> shower::TCShowerAlg::totalEnergy

Definition at line 42 of file TCShowerAlg.h.

std::vector<double> shower::TCShowerAlg::totalEnergyErr

Definition at line 43 of file TCShowerAlg.h.

TVector3 shower::TCShowerAlg::xyzErr

Definition at line 41 of file TCShowerAlg.h.


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