18 std::vector<art::Ptr<recob::Hit>>
hits;
19 std::vector<int> clsIDs;
23 bool comparePFP(
const pfpStuff& l,
const pfpStuff&
r)
28 double const lz = l.hits.size();
29 double const rz = r.hits.size();
32 constexpr
int hitthres = 80;
34 if (lz > hitthres && rz <= hitthres)
return false;
36 if (lz <= hitthres && rz > hitthres)
return true;
45 : fCalorimetryAlg(pset.
get<
fhicl::ParameterSet>(
"CalorimetryAlg"))
46 , fProjectionMatchingAlg(pset.
get<
fhicl::ParameterSet>(
"ProjectionMatchingAlg"))
67 std::vector<pfpStuff> allpfps;
70 for (
size_t i = 0; i < pfplist.size(); ++i) {
73 thispfp.clsIDs.clear();
74 thispfp.pfp = pfplist[i];
76 std::vector<art::Ptr<recob::Vertex>> thisvtxlist = vtxpfp_fm.at(pfplist[i].key());
77 if (thisvtxlist.size()) thispfp.vtx = thisvtxlist[0];
80 if (thistrklist.size()) thispfp.trk = thistrklist[0];
83 std::vector<int> clustersize;
85 for (
size_t j = 0; j < thisclusterlist.size(); ++j) {
87 thispfp.clsIDs.push_back(thisclusterlist[j]->ID());
89 std::vector<art::Ptr<recob::Hit>> thishitlist = cls_fm.at(thisclusterlist[j].key());
90 clustersize.push_back((
int)thishitlist.size());
92 for (
size_t k = 0; k < thishitlist.size(); ++k) {
93 thispfp.hits.push_back(thishitlist[k]);
98 if (clustersize.size() == 3) {
99 if (!thispfp.vtx)
continue;
100 if (!thispfp.trk)
continue;
102 allpfps.push_back(thispfp);
107 std::cout <<
"pfp " << thispfp.pfp->Self() + 1 <<
" cluster sizes " << clustersize[0] <<
":" 108 << clustersize[1] <<
":" << clustersize[2] <<
" vertex " << thispfp.vtx->ID()
109 <<
" " << tick <<
":" << wire <<
" z " << thispfp.vtx->position().Z()
116 std::sort(allpfps.begin(), allpfps.end(), comparePFP);
117 std::reverse(allpfps.begin(), allpfps.end());
119 std::cout <<
"sorted pfps: ";
120 for (
size_t i = 0; i < allpfps.size(); ++i)
121 std::cout << allpfps[i].pfp->Self() + 1 <<
" ";
122 std::cout << std::endl;
124 bool showerCandidate =
false;
126 for (
size_t i = 0; i < allpfps.size(); ++i) {
132 std::vector<art::Ptr<recob::Hit>> pfphits = allpfps[i].hits;
133 std::vector<art::Ptr<recob::Cluster>> pfpcls = clspfp_fm.at(allpfps[i].pfp.key());
135 std::cout <<
"pfp " << allpfps[i].pfp->Self() + 1 <<
" hits " << pfphits.size() << std::endl;
142 if (pfptrk->
Vertex().Z() < pfptrk->
End().Z()) {
152 double pullTolerance = 0.7;
154 double minDistVert = 15;
160 if (pfphits.size() < 30)
continue;
161 if (pfphits.size() > 500)
continue;
163 if (pfphits.size() < 90) {
167 if (pfphits.size() > 400)
169 else if (pfphits.size() > 100)
173 for (
size_t ii = 0; ii < pfphits.size(); ++ii) {
178 double showerHitPull = 0;
184 std::map<geo::PlaneID, double> trk_tick1;
185 std::map<geo::PlaneID, double> trk_wire1;
188 std::map<geo::PlaneID, double> trk_tick2;
189 std::map<geo::PlaneID, double> trk_wire2;
198 for (
size_t j = 0; j < clusterlist.size(); ++j) {
199 std::vector<art::Ptr<recob::Hit>> cls_hitlist = cls_fm.at(clusterlist[j].key());
201 if (clusterlist[j]->ID() > 0 && cls_hitlist.size() > 10)
continue;
202 if (cls_hitlist.size() > 50)
continue;
204 bool isGoodCluster =
false;
208 for (
size_t k = 0; k < allpfps[i].clsIDs.size(); ++k) {
209 if (allpfps[i].clsIDs[k] == clusterlist[j]->ID()) skipit =
true;
211 if (skipit)
continue;
213 for (
size_t jj = 0; jj < cls_hitlist.size(); ++jj) {
214 int isGoodHit =
goodHit(clockData,
224 if (isGoodHit == -1) {
225 isGoodCluster =
false;
228 else if (isGoodHit == 1) {
229 isGoodCluster =
true;
236 for (
size_t jj = 0; jj < cls_hitlist.size(); ++jj) {
239 int showerHitPullAdd = 0;
250 showerHitPull += showerHitPullAdd;
259 showerHitPull /= nShowerHits;
261 std::cout <<
"shower hits " <<
showerHits.size() <<
" " << nShowerHits <<
" shower pull " 262 << showerHitPull << std::endl;
264 if (nShowerHits > tolerance &&
std::abs(showerHitPull) < pullTolerance) {
265 showerCandidate =
true;
266 std::cout <<
"SHOWER CANDIDATE" << std::endl;
268 if (nShowerHits > 400) maxDist *= 2;
269 for (
size_t k = 0; k < hitlist.size(); ++k) {
270 std::vector<art::Ptr<recob::Cluster>> hit_clslist = hitcls_fm.at(hitlist[k].key());
271 if (hit_clslist.size())
continue;
273 int isGoodHit =
goodHit(clockData,
287 for (
size_t k = 0; k < clusterlist.size(); ++k) {
288 std::vector<art::Ptr<recob::Hit>> cls_hitlist = cls_fm.at(clusterlist[k].key());
289 if (clusterlist[k]->ID() > 0 && cls_hitlist.size() > 50)
continue;
291 double thisDist = maxDist;
292 double thisMin = minDistVert;
294 if (cls_hitlist.size() < 10) {
298 else if (cls_hitlist.size() < 30)
305 for (
size_t kk = 0; kk < cls_hitlist.size(); ++kk) {
307 int isGoodHit =
goodHit(clockData,
316 if (isGoodHit == -1) {
320 else if (isGoodHit == 1) {
325 double fracGood = (double)ngoodhits / nhits;
327 bool isGoodTrack = fracGood > 0.4;
330 for (
size_t kk = 0; kk < cls_hitlist.size(); ++kk) {
339 if (showerCandidate) {
340 std::cout <<
"THIS IS THE SHOWER PFP: " << allpfps[i].pfp->Self() + 1 << std::endl;
346 if (showerCandidate)
return 1;
359 double const maxDist,
360 double const minDistVert,
361 std::map<geo::PlaneID, double>
const& trk_wire1,
362 std::map<geo::PlaneID, double>
const& trk_tick1,
363 std::map<geo::PlaneID, double>
const& trk_wire2,
364 std::map<geo::PlaneID, double>
const& trk_tick2)
const 389 double const minDistVert,
390 std::map<geo::PlaneID, double>
const& trk_wire1,
391 std::map<geo::PlaneID, double>
const& trk_tick1,
392 std::map<geo::PlaneID, double>
const& trk_wire2,
393 std::map<geo::PlaneID, double>
const& trk_tick2,
401 double UnitsPerTick = tickToDist / wirePitch;
404 double y0 = hit->
PeakTime() * UnitsPerTick;
406 double x1 = trk_wire1.at(hit->
WireID());
407 double y1 = trk_tick1.at(hit->
WireID()) * UnitsPerTick;
409 double x2 = trk_wire2.at(hit->
WireID());
410 double y2 = trk_tick2.at(hit->
WireID()) * UnitsPerTick;
412 double distToVert = std::hypot(x0 - x1, y0 - y1);
413 if (distToVert < minDistVert)
return -1;
418 else if (distToVert < 100)
423 double a = std::hypot(x2 - x1, y2 - y1);
424 double b = std::hypot(x0 - x1, y0 - y1);
425 double c = std::hypot(x0 - x2, y0 - y2);
426 double costheta = -(pow(c, 2) - pow(a, 2) - pow(b, 2)) / (2 * a * b);
427 if (costheta < 0)
return -1;
430 std::abs((y2 - y1) * x0 - (x2 - x1) * y0 + x2 * y1 - y2 * x1) / std::hypot(y2 - y1, x2 - x1);
432 if (dist < maxDist) {
433 if (((y2 - y1) * x0 - (x2 - x1) * y0 + x2 * y1 - y2 * x1) > 0)
450 for (
size_t i = 0; i < showerhits.size(); ++i) {
451 if (hit.
key() == showerhits[i].key())
return false;
details::range_type< T > Iterate() const
Initializes the specified ID with the ID of the first cryostat.
Length_t WireCoordinate(Point_t const &pos, PlaneID const &planeid) const
Returns the index of the nearest wire to the specified position.
Utilities related to art service access.
std::vector< double > totalEnergyErr
bool addShowerHit(art::Ptr< recob::Hit > hit, std::vector< art::Ptr< recob::Hit >> showerhits) const
Float_t y1[n_points_granero]
Float_t x1[n_points_granero]
std::vector< double > dEdx
The data type to uniquely identify a Plane.
double Temperature() const
In kelvin.
constexpr auto abs(T v)
Returns the absolute value of the argument.
Vector_t VertexDirection() const
Access to track direction at different points.
WireID_t Wire
Index of the wire within its plane.
geo::WireID const & WireID() const
Initial tdc tick for hit.
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)
Float_t y2[n_points_geant4]
double Efield(unsigned int planegap=0) const
kV/cm
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.
Access the description of detector geometry.
key_type key() const noexcept
Point_t const & Vertex() const
Access to track position at different points.
tick_as<> tick
Tick number, represented by std::ptrdiff_t.
double ConvertXToTicks(double X, int p, int t, int c) const
double DriftVelocity(double efield=0., double temperature=0.) const
cm/us
Declaration of cluster object.
Definition of data types for geometry description.
Detector simulation of raw signals on wires.
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
float PeakTime() const
Time of the signal peak, in tick units.
decltype(auto) get(T &&obj)
ADL-aware version of std::to_string.
Contains all timing reference information for the detector.
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
Vector_t EndDirection() const
Access to track direction at different points.
std::vector< double > totalEnergy
Point_t const & End() const
Access to track position at different points.
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
Float_t x2[n_points_geant4]
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
std::vector< art::Ptr< recob::Hit > > showerHits
std::vector< double > dEdxErr
Length_t WirePitch(PlaneID const &planeid=plane_zero) const
Returns the distance between two consecutive wires.
const Point_t & position() const
Return vertex 3D position.
art framework interface to geometry description