72 std::vector<pfpStuff> allpfps;
75 for (
size_t i = 0; i < pfplist.size(); ++i) {
78 thispfp.clsIDs.clear();
79 thispfp.pfp = pfplist[i];
81 std::vector<art::Ptr<recob::Vertex>> thisvtxlist = vtxpfp_fm.at(pfplist[i].
key());
82 if (thisvtxlist.size()) thispfp.vtx = thisvtxlist[0];
85 if (thistrklist.size()) thispfp.trk = thistrklist[0];
88 std::vector<int> clustersize;
90 for (
size_t j = 0; j < thisclusterlist.size(); ++j) {
92 thispfp.clsIDs.push_back(thisclusterlist[j]->
ID());
94 std::vector<art::Ptr<recob::Hit>> thishitlist = cls_fm.at(thisclusterlist[j].
key());
95 clustersize.push_back((
int)thishitlist.size());
97 for (
size_t k = 0;
k < thishitlist.size(); ++
k) {
98 thispfp.hits.push_back(thishitlist[
k]);
103 if (clustersize.size() == 3) {
104 if (!thispfp.vtx)
continue;
105 if (!thispfp.trk)
continue;
107 allpfps.push_back(thispfp);
109 double tick = detProp.ConvertXToTicks(thispfp.vtx->position().X(),
geo::PlaneID(0, 0, 2));
111 thispfp.vtx->position().Y(), thispfp.vtx->position().Z(),
geo::PlaneID(0, 0, 2));
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()
122 std::sort(allpfps.begin(), allpfps.end(), comparePFP);
125 std::cout <<
"sorted pfps: ";
126 for (
size_t i = 0; i < allpfps.size(); ++i)
127 std::cout << allpfps[i].pfp->Self() + 1 <<
" ";
130 bool showerCandidate =
false;
132 for (
size_t i = 0; i < allpfps.size(); ++i) {
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());
141 std::cout <<
"pfp " << allpfps[i].pfp->Self() + 1 <<
" hits " << pfphits.size() <<
std::endl;
148 if (pfptrk->
Vertex().Z() < pfptrk->
End().Z()) {
158 double pullTolerance = 0.7;
160 double minDistVert = 15;
166 if (pfphits.size() < 30)
continue;
167 if (pfphits.size() > 500)
continue;
169 if (pfphits.size() < 90) {
173 if (pfphits.size() > 400)
175 else if (pfphits.size() > 100)
179 for (
size_t ii = 0; ii < pfphits.size(); ++ii) {
184 double showerHitPull = 0;
186 TVector3 pfpStart =
shwvtx;
190 std::map<geo::PlaneID, double> trk_tick1;
191 std::map<geo::PlaneID, double> trk_wire1;
194 std::map<geo::PlaneID, double> trk_tick2;
195 std::map<geo::PlaneID, double> trk_wire2;
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);
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());
207 if (clusterlist[j]->
ID() > 0 && cls_hitlist.size() > 10)
continue;
208 if (cls_hitlist.size() > 50)
continue;
210 bool isGoodCluster =
false;
214 for (
size_t k = 0;
k < allpfps[i].clsIDs.size(); ++
k) {
215 if (allpfps[i].clsIDs[
k] == clusterlist[j]->
ID()) skipit =
true;
217 if (skipit)
continue;
219 for (
size_t jj = 0; jj < cls_hitlist.size(); ++jj) {
220 int isGoodHit =
goodHit(clockData,
230 if (isGoodHit == -1) {
231 isGoodCluster =
false;
234 else if (isGoodHit == 1) {
235 isGoodCluster =
true;
242 for (
size_t jj = 0; jj < cls_hitlist.size(); ++jj) {
245 int showerHitPullAdd = 0;
256 showerHitPull += showerHitPullAdd;
265 showerHitPull /= nShowerHits;
267 std::cout <<
"shower hits " <<
showerHits.size() <<
" " << nShowerHits <<
" shower pull " 270 if (nShowerHits > tolerance &&
std::abs(showerHitPull) < pullTolerance) {
271 showerCandidate =
true;
272 std::cout <<
"SHOWER CANDIDATE" <<
std::endl;
274 if (nShowerHits > 400) maxDist *= 2;
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;
279 int isGoodHit =
goodHit(clockData,
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;
297 double thisDist = maxDist;
298 double thisMin = minDistVert;
300 if (cls_hitlist.size() < 10) {
304 else if (cls_hitlist.size() < 30)
311 for (
size_t kk = 0; kk < cls_hitlist.size(); ++kk) {
313 int isGoodHit =
goodHit(clockData,
322 if (isGoodHit == -1) {
326 else if (isGoodHit == 1) {
331 double fracGood = (double)ngoodhits / nhits;
333 bool isGoodTrack = fracGood > 0.4;
336 for (
size_t kk = 0; kk < cls_hitlist.size(); ++kk) {
345 if (showerCandidate) {
346 std::cout <<
"THIS IS THE SHOWER PFP: " << allpfps[i].pfp->Self() + 1 <<
std::endl;
352 if (showerCandidate)
return 1;
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
bool addShowerHit(art::Ptr< recob::Hit > hit, std::vector< art::Ptr< recob::Hit >> showerhits) const
std::vector< double > dEdx
Vector_t VertexDirection() const
Point_t const & Vertex() const
tick_as<> tick
Tick number, represented by std::ptrdiff_t.
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)
Vector_t EndDirection() const
std::vector< double > totalEnergy
Point_t const & End() const
detail::Node< FrameID, bool > PlaneID
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
std::vector< double > dEdxErr
const Point_t & position() const
Return vertex 3D position.
QTextStream & endl(QTextStream &s)