70 art::ServiceHandle<geo::Geometry const> geom;
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];
84 std::vector<art::Ptr<recob::Track>> thistrklist = trkpfp_fm.at(pfplist[i].key());
85 if (thistrklist.size()) thispfp.trk = thistrklist[0];
87 std::vector<art::Ptr<recob::Cluster>> thisclusterlist = clspfp_fm.at(pfplist[i].key());
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);
110 int wire = geom->WireCoordinate(
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);
123 std::reverse(allpfps.begin(), allpfps.end());
126 for (
size_t i = 0; i < allpfps.size(); ++i)
127 std::cout << allpfps[i].pfp->Self() + 1 <<
" ";
128 std::cout << std::endl;
130 bool showerCandidate =
false;
132 for (
size_t i = 0; i < allpfps.size(); ++i) {
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());
141 std::cout <<
"pfp " << allpfps[i].pfp->Self() + 1 <<
" hits " << pfphits.size() << std::endl;
144 vtx[0] = pfpvtx->position().X();
145 vtx[1] = pfpvtx->position().Y();
146 vtx[2] = pfpvtx->position().Z();
148 if (pfptrk->Vertex().Z() < pfptrk->End().Z()) {
149 shwvtx = pfptrk->Vertex<TVector3>();
150 shwDir = pfptrk->VertexDirection<TVector3>();
153 shwvtx = pfptrk->End<TVector3>();
154 shwDir = -pfptrk->EndDirection<TVector3>();
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;
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);
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 "
268 << showerHitPull << std::endl;
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;
std::vector< double > totalEnergyErr
bool addShowerHit(art::Ptr< recob::Hit > hit, std::vector< art::Ptr< recob::Hit >> showerhits) const
std::vector< double > dEdx
IDparameter< geo::PlaneID > PlaneID
Member type of validated geo::PlaneID parameter.
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
tick_as<> tick
Tick number, represented by std::ptrdiff_t.
std::vector< double > totalEnergy
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
BEGIN_PROLOG could also be cout