137 auto showers = std::make_unique<std::vector<recob::Shower>>();
138 auto spacePoints = std::make_unique<std::vector<recob::SpacePoint>>();
139 auto clusterAssociations = std::make_unique<art::Assns<recob::Shower, recob::Cluster>>();
140 auto hitShowerAssociations = std::make_unique<art::Assns<recob::Shower, recob::Hit>>();
141 auto trackAssociations = std::make_unique<art::Assns<recob::Shower, recob::Track>>();
142 auto spShowerAssociations = std::make_unique<art::Assns<recob::Shower, recob::SpacePoint>>();
143 auto hitSpAssociations = std::make_unique<art::Assns<recob::SpacePoint, recob::Hit>>();
149 std::vector<art::Ptr<recob::Hit>>
hits;
154 std::vector<art::Ptr<recob::Track>> tracks;
159 std::vector<art::Ptr<recob::Cluster>> clusters;
165 std::vector<art::Ptr<recob::PFParticle>> pfps;
170 std::vector<art::Ptr<recob::Vertex>> vertices;
183 std::vector<std::vector<int>> newShowers;
184 std::vector<unsigned int> pfParticles;
186 std::map<int, std::vector<int>> clusterToTracks;
187 std::map<int, std::vector<int>> trackToClusters;
198 std::vector<int> clustersToIgnore;
202 if (clustersToIgnore.size() > 0) {
203 clusterToTracks.clear();
204 trackToClusters.clear();
206 clusters, fmh, fmt, clustersToIgnore, clusterToTracks, trackToClusters);
210 newShowers = initialShowers;
217 for (
size_t ipfp = 0; ipfp < pfps.size(); ++ipfp) {
224 int trkLikeIdx = hitResults->getIndex(
"track");
225 int emLikeIdx = hitResults->getIndex(
"em");
226 if ((trkLikeIdx < 0) || (emLikeIdx < 0)) {
227 throw cet::exception(
"EMShower") <<
"No em/track labeled columns in MVA data products.";
229 if (fmcp.isValid()) {
230 std::vector<art::Ptr<recob::Hit>> pfphits;
231 std::vector<art::Ptr<recob::Cluster>> clus = fmcp.at(ipfp);
232 for (
size_t iclu = 0; iclu < clus.size(); ++iclu) {
233 std::vector<art::Ptr<recob::Hit>> ClusterHits = fmh.at(clus[iclu].key());
234 pfphits.insert(pfphits.end(), ClusterHits.begin(), ClusterHits.end());
236 if (pfphits.size()) {
237 auto vout = hitResults->getOutput(pfphits);
238 double trk_like = -1, trk_or_em = vout[trkLikeIdx] + vout[emLikeIdx];
240 trk_like = vout[trkLikeIdx] / trk_or_em;
242 std::vector<int> clusters;
243 for (
size_t iclu = 0; iclu < clus.size(); ++iclu) {
244 clusters.push_back(clus[iclu].key());
246 if (clusters.size()) {
247 newShowers.push_back(clusters);
248 pfParticles.push_back(ipfp);
255 throw cet::exception(
"EMShower") <<
"Cannot get associated cluster for PFParticle " 259 else if (pfp->
PdgCode() == 11) {
260 if (fmcp.isValid()) {
261 std::vector<int> clusters;
262 std::vector<art::Ptr<recob::Cluster>> clus = fmcp.at(ipfp);
263 for (
size_t iclu = 0; iclu < clus.size(); ++iclu) {
264 clusters.push_back(clus[iclu].key());
266 if (clusters.size()) {
267 newShowers.push_back(clusters);
268 pfParticles.push_back(ipfp);
277 for (
auto const& newShower : newShowers) {
282 if (
fDebug > 0) std::cout <<
"\n\nStart shower " << showerNum <<
'\n';
290 std::vector<int> associatedTracks;
293 for (
int const showerCluster : newShower) {
300 std::vector<art::Ptr<recob::Hit>>
const& showerClusterHits = fmh.at(cluster.
key());
307 int trkLikeIdx = hitResults->getIndex(
"track");
308 int emLikeIdx = hitResults->getIndex(
"em");
309 if (trkLikeIdx < 0 || emLikeIdx < 0) {
310 throw cet::exception(
"EMShower") <<
"No em/track labeled columns in MVA data products.";
312 for (
auto const& showerHit : showerClusterHits) {
313 auto vout = hitResults->getOutput(showerHit);
314 double trk_like = -1, trk_or_em = vout[trkLikeIdx] + vout[emLikeIdx];
316 trk_like = vout[trkLikeIdx] / trk_or_em;
324 for (
auto const& showerClusterHit : showerClusterHits)
329 for (
int const clusterTrack : clusterToTracks.at(showerCluster))
330 if (not cet::search_all(associatedTracks, clusterTrack))
331 associatedTracks.push_back(clusterTrack);
337 for (
int const trackIndex : associatedTracks) {
338 showerTracks.
push_back(tracks.at(trackIndex));
343 if (fmspp.isValid()) {
344 for (
size_t ihit = 0; ihit < showerHits.
size(); ++ihit) {
345 for (
auto const& spPtr : fmspp.at(ihit))
355 std::cout <<
" ------------------ Ordering shower hits --------------------\n";
356 std::map<int, std::vector<art::Ptr<recob::Hit>>> showerHitsMap =
359 std::cout <<
" ------------------ End ordering shower hits " 360 "--------------------\n";
363 std::unique_ptr<recob::Track> initialTrack;
364 std::map<int, std::vector<art::Ptr<recob::Hit>>> initialTrackHits;
368 std::vector<std::vector<art::Ptr<recob::Hit>>> hitAssns;
369 std::vector<recob::SpacePoint> showerSpacePoints;
373 for (
auto const& trkPtr : showerTracks) {
374 for (
auto const& trackSpacePoint :
375 fmsp.at(trkPtr.key()) | ranges::views::transform(
to_element)) {
376 showerSpacePoints.push_back(trackSpacePoint);
384 size_t firstSpacePoint = spacePoints->size(), nSpacePoint = 0;
385 for (
auto const& ssp : showerSpacePoints) {
386 spacePoints->emplace_back(ssp.XYZ(), ssp.ErrXYZ(), ssp.Chisq(), spacePoints->size());
387 auto const index = spacePoints->size() - 1;
388 auto const space_point_ptr = make_space_point_ptr(index);
389 cet::for_all(hitAssns.at(nSpacePoint), AddMany{space_point_ptr, *hitSpAssociations});
391 auto const lastSpacePoint = spacePoints->size();
399 showers->push_back(shower);
402 cet::for_all(showerHits, AddMany{shower_ptr, *hitShowerAssociations});
403 cet::for_all(showerClusters, AddMany{shower_ptr, *clusterAssociations});
404 cet::for_all(showerTracks, AddMany{shower_ptr, *trackAssociations});
405 for (
size_t i = firstSpacePoint; i < lastSpacePoint; ++i) {
406 spShowerAssociations->addSingle(shower_ptr, make_space_point_ptr(i));
410 mf::LogInfo(
"EMShower") <<
"Discarding shower " << showerNum
411 <<
" due to incompleteness (SaveNonCompleteShowers == false)";
416 if (vertices.size()) {
420 for (
auto const& vtx : vertices) {
421 auto const pos = vtx->position();
422 if (pos.Z() < nuvtx.Z()) { nuvtx = pos; }
426 double mindist2 = DBL_MAX;
427 for (
auto const& sp : showerSpacePoints_p | ranges::views::transform(
to_element)) {
428 double const dist2 = cet::sum_of_squares(
429 nuvtx.X() - sp.XYZ()[0], nuvtx.Y() - sp.XYZ()[1], nuvtx.Z() - sp.XYZ()[2]);
430 if (dist2 < mindist2) {
432 shwvtx.SetXYZ(sp.XYZ()[0], sp.XYZ()[1], sp.XYZ()[2]);
438 for (
auto const& vtx : vertices) {
439 auto const pos = vtx->position();
441 cet::sum_of_squares(pos.X() - shwvtx.X(), pos.Y() - shwvtx.Y(), pos.Z() - shwvtx.Z());
442 if (dist2 < mindist2) {
452 showers->push_back(shower);
453 auto const index = showers->size() - 1;
454 showers->back().
set_id(index);
457 cet::for_all(showerHits, AddMany{shower_ptr, *hitShowerAssociations});
458 cet::for_all(showerClusters, AddMany{shower_ptr, *clusterAssociations});
459 cet::for_all(showerTracks, AddMany{shower_ptr, *trackAssociations});
460 cet::for_all(showerSpacePoints_p, AddMany{shower_ptr, *spShowerAssociations});
467 evt.
put(std::move(showers));
468 evt.
put(std::move(spacePoints));
469 evt.
put(std::move(hitShowerAssociations));
470 evt.
put(std::move(clusterAssociations));
471 evt.
put(std::move(trackAssociations));
472 evt.
put(std::move(spShowerAssociations));
473 evt.
put(std::move(hitSpAssociations));
const TVector3 & ShowerStart() const
std::vector< recob::SpacePoint > MakeSpacePoints(detinfo::DetectorPropertiesData const &detProp, std::map< int, std::vector< art::Ptr< recob::Hit >>> const &hits, std::vector< std::vector< art::Ptr< recob::Hit >>> &hitAssns) const
Makes space points from the shower hits in each plane.
bool const fUseCNNtoIDEMPFP
constexpr to_element_t to_element
art::InputTag const fVertexModuleLabel
art::InputTag const fClusterModuleLabel
recob::Shower MakeShower(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, art::PtrVector< recob::Hit > const &hits, std::unique_ptr< recob::Track > const &initialTrack, std::map< int, std::vector< art::Ptr< recob::Hit >>> const &initialTrackHits) const
Makes a recob::Shower object given the hits in the shower and the initial track-like part...
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
bool const fMakeSpacePoints
bool const fUseCNNtoIDEMHit
void AssociateClustersAndTracks(std::vector< art::Ptr< recob::Cluster >> const &clusters, art::FindManyP< recob::Hit > const &fmh, art::FindManyP< recob::Track > const &fmt, std::map< int, std::vector< int >> &clusterToTracks, std::map< int, std::vector< int >> &trackToClusters) const
Map associated tracks and clusters together given their associated hits.
int PdgCode() const
Return the type of particle as a PDG ID.
art::InputTag const fTrackModuleLabel
Cluster finding and building.
bool const fFindBadPlanes
std::vector< std::vector< int > > FindShowers(std::map< int, std::vector< int >> const &trackToClusters) const
Makes showers given a map between tracks and all clusters associated with them.
recob::tracking::Point_t Point_t
art::InputTag const fPFParticleModuleLabel
bool isValid() const noexcept
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
void set_id(const int id)
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
void push_back(Ptr< U > const &p)
void FindInitialTrack(detinfo::DetectorPropertiesData const &detProp, std::map< int, std::vector< art::Ptr< recob::Hit >>> const &hits, std::unique_ptr< recob::Track > &initialTrack, std::map< int, std::vector< art::Ptr< recob::Hit >>> &initialTrackHits) const
std::vector< int > CheckShowerPlanes(std::vector< std::vector< int >> const &initialShowers, std::vector< art::Ptr< recob::Cluster >> const &clusters, art::FindManyP< recob::Hit > const &fmh) const
Takes the initial showers found and tries to resolve issues where one bad view ruins the event...
key_type key() const noexcept
std::map< int, std::vector< art::Ptr< recob::Hit > > > OrderShowerHits(detinfo::DetectorPropertiesData const &detProp, art::PtrVector< recob::Hit > const &shower, int plane) const
Takes the hits associated with a shower and orders them so they follow the direction of the shower...
art::InputTag const fHitsModuleLabel
double const fMinTrackLikeScore
EMShowerAlg const fEMShowerAlg
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
art::InputTag const fCNNEMModuleLabel
bool const fSaveNonCompleteShowers
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< Coord_t >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space. See recob::tracking::Coord_t for more detai...
static std::unique_ptr< MVAReader > create(const art::Event &evt, const art::InputTag &tag)
cet::coded_exception< error, detail::translate > exception
1.8.11 
