diff --git a/sbncode/CAFMaker/CAFMaker_module.cc b/sbncode/CAFMaker/CAFMaker_module.cc index ab33a7fa9..494de6378 100644 --- a/sbncode/CAFMaker/CAFMaker_module.cc +++ b/sbncode/CAFMaker/CAFMaker_module.cc @@ -2523,7 +2523,7 @@ void CAFMaker::produce(art::Event& evt) noexcept { assert(thisShower.size() == 1); SRShower& shw = pfp.shw; - FillShowerVars(*thisShower[0], vertex, fmShowerHit.at(iPart), wireReadout, producer, shw); + FillShowerVars(*thisShower[0], vertex, fmShowerHit.at(iPart), wireReadout, producer, shw, fDet); // We may have many residuals per shower depending on how many showers ar in the slice if (fmShowerRazzle.isValid() && fmShowerRazzle.at(iPart).size()==1) { diff --git a/sbncode/CAFMaker/FillReco.cxx b/sbncode/CAFMaker/FillReco.cxx index 49c4f7eb6..7830b48b6 100644 --- a/sbncode/CAFMaker/FillReco.cxx +++ b/sbncode/CAFMaker/FillReco.cxx @@ -391,6 +391,7 @@ namespace caf const geo::WireReadoutGeom& wireReadout, unsigned producer, caf::SRShower &srshower, + Det_t det, bool allowEmpty) { @@ -411,11 +412,56 @@ namespace caf // It's sth like this but not quite. And will need to pass a simb::MCtruth object vtx position anyway. // srshower.conversion_gap = (shower.ShowerStart() - vertex.Position()).Mag(); - if(shower.best_plane() != -999){ - srshower.bestplane = shower.best_plane(); - srshower.bestplane_dEdx = srshower.plane[shower.best_plane()].dEdx; - srshower.bestplane_energy = srshower.plane[shower.best_plane()].energy; - } + for(int p = 0; p < 3; ++p) srshower.plane[p].nHits = 0; + for (auto const& hit:hits) ++srshower.plane[hit->WireID().Plane].nHits; + + if(det == kSBND) + { + int bestplane_for_energy = -999; + int mosthits = -1; + for(int p = 0; p < 3; ++p) + { + if((int)srshower.plane[p].nHits > mosthits) + { + mosthits = srshower.plane[p].nHits; + bestplane_for_energy = p; + } + } + + if(bestplane_for_energy != -999) + { + srshower.bestplane_for_energy = bestplane_for_energy; + srshower.bestplane_energy = srshower.plane[bestplane_for_energy].energy; + } + + if(shower.best_plane() != -999 && srshower.plane[shower.best_plane()].dEdx != -999) + { + srshower.bestplane_for_dedx = shower.best_plane(); + srshower.bestplane_dEdx = srshower.plane[shower.best_plane()].dEdx; + } + else + { + for(int p = 2; p >= 0; --p) + { + if(srshower.plane[p].dEdx != -999) + { + srshower.bestplane_for_dedx = p; + srshower.bestplane_dEdx = srshower.plane[shower.best_plane()].dEdx; + break; + } + } + } + } + else + { + if(shower.best_plane() != -999) + { + srshower.bestplane_for_energy = shower.best_plane(); + srshower.bestplane_for_dedx = shower.best_plane(); + srshower.bestplane_dEdx = srshower.plane[shower.best_plane()].dEdx; + srshower.bestplane_energy = srshower.plane[shower.best_plane()].energy; + } + } if(shower.has_open_angle()) srshower.open_angle = shower.OpenAngle(); @@ -438,9 +484,6 @@ namespace caf srshower.end = shower.ShowerStart()+ (shower.Length() * shower.Direction()); } - for(int p = 0; p < 3; ++p) srshower.plane[p].nHits = 0; - for (auto const& hit:hits) ++srshower.plane[hit->WireID().Plane].nHits; - for (geo::PlaneGeo const& plane: wireReadout.Iterate()) { const double angleToVert(wireReadout.WireAngleToVertical(plane.View(), plane.ID()) - 0.5*M_PI); diff --git a/sbncode/CAFMaker/FillReco.h b/sbncode/CAFMaker/FillReco.h index 25b2bb867..b1f8d9b05 100644 --- a/sbncode/CAFMaker/FillReco.h +++ b/sbncode/CAFMaker/FillReco.h @@ -70,6 +70,7 @@ namespace caf const geo::WireReadoutGeom& wireReadout, unsigned producer, caf::SRShower& srshower, + Det_t det, bool allowEmpty = false); void FillShowerRazzle(const art::Ptr razzle,