286 bool mc = !
evt.isRealData();
291 art::ServiceHandle<geo::Geometry const> geom;
295 art::PtrVector<recob::Hit> hits;
301 art::Handle<std::vector<recob::Cluster>> clusterh;
307 if (clusterh.isValid()) {
308 int nclus = clusterh->size();
310 for (
int i = 0; i < nclus; ++i) {
311 art::Ptr<recob::Cluster> pclus(clusterh, i);
312 std::vector<art::Ptr<recob::Hit>> clushits = fm.at(i);
313 int nhits = clushits.size();
314 hits.reserve(hits.size() + nhits);
316 for (
std::vector<art::Ptr<recob::Hit>>::const_iterator ihit = clushits.begin();
317 ihit != clushits.end();
319 hits.push_back(*ihit);
328 art::Handle<std::vector<recob::Hit>> hith;
330 if (hith.isValid()) {
331 int nhits = hith->size();
334 for (
int i = 0; i < nhits; ++i)
335 hits.push_back(art::Ptr<recob::Hit>(hith, i));
341 auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(
evt);
343 art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(
evt, clockData);
347 art::ServiceHandle<cheat::BackTrackerService const> bt_serv;
351 for (
auto const& hitPtr : hits) {
360 std::vector<double> hitxyz = bt_serv->HitToXYZ(clockData, hit);
364 catch (cet::exception&
x) {
369 fHDTUE->Fill(tpeak - tav);
373 fHDTVE->Fill(tpeak - tav);
377 fHDTWE->Fill(tpeak - tav);
381 throw cet::exception(
"SpacePointAna") <<
"Bad view = " << hit.
View() <<
"\n";
386 std::vector<recob::SpacePoint> spts1;
387 std::vector<recob::SpacePoint> spts2;
388 std::vector<recob::SpacePoint> spts3;
401 MF_LOG_DEBUG(
"SpacePointAna")
402 <<
"Found " << spts1.size() <<
" space points using special time cut.";
416 MF_LOG_DEBUG(
"SpacePointAna")
417 <<
"Found " << spts2.size() <<
" space points using special seperation cut.";
429 MF_LOG_DEBUG(
"SpacePointAna") <<
"Found " << spts3.size()
430 <<
" space points using default cuts.";
436 for (
auto const& spt : spts1) {
437 if (spt.XYZ()[0] <
fMinX || spt.XYZ()[0] >
fMaxX || spt.XYZ()[1] <
fMinY ||
449 unsigned int tpc1 = hit1WireID.
TPC;
450 unsigned int plane1 = hit1WireID.
Plane;
451 unsigned int wire1 = hit1WireID.
Wire;
458 unsigned int tpc2 = hit2WireID.
TPC;
459 unsigned int plane2 = hit2WireID.
Plane;
460 unsigned int wire2 = hit2WireID.
Wire;
464 if (tpc1 == tpc2 && plane1 != plane2) {
472 fHDTUVU->Fill(
double(wire1), t1 - t2);
473 fHDTUVV->Fill(
double(wire2), t1 - t2);
477 fHDTWUW->Fill(
double(wire2), t2 - t1);
478 fHDTWUU->Fill(
double(wire1), t2 - t1);
484 fHDTVWV->Fill(
double(wire1), t1 - t2);
485 fHDTVWW->Fill(
double(wire2), t1 - t2);
489 fHDTUVU->Fill(
double(wire2), t2 - t1);
490 fHDTUVV->Fill(
double(wire1), t2 - t1);
496 fHDTWUW->Fill(
double(wire1), t1 - t2);
497 fHDTWUU->Fill(
double(wire2), t1 - t2);
501 fHDTVWV->Fill(
double(wire2), t2 - t1);
502 fHDTVWW->Fill(
double(wire1), t2 - t1);
512 for (
auto const& spt : spts2) {
513 if (spt.XYZ()[0] <
fMinX || spt.XYZ()[0] >
fMaxX || spt.XYZ()[1] <
fMinY ||
530 for (
auto const& spt : spts3) {
531 if (spt.XYZ()[0] <
fMinX || spt.XYZ()[0] >
fMaxX || spt.XYZ()[1] <
fMinY ||
536 fHx->Fill(spt.XYZ()[0]);
537 fHy->Fill(spt.XYZ()[1]);
538 fHz->Fill(spt.XYZ()[2]);
542 std::vector<art::Ptr<recob::Hit>> spthits;
544 for (
auto const& ptr : av_spthits) {
545 spthits.push_back(ptr);
550 for (
auto const& hitPtr : spthits) {
574 art::ServiceHandle<cheat::BackTrackerService const> bt_serv;
577 std::vector<double> mcxyz = bt_serv->SpacePointHitsToWeightedXYZ(clockData, spthits);
578 fHMCdx->Fill(spt.XYZ()[0] - mcxyz[0]);
579 fHMCdy->Fill(spt.XYZ()[1] - mcxyz[1]);
580 fHMCdz->Fill(spt.XYZ()[2] - mcxyz[2]);
581 if (spt.ErrXYZ()[0] > 0.)
582 fHMCxpull->Fill((spt.XYZ()[0] - mcxyz[0]) / std::sqrt(spt.ErrXYZ()[0]));
583 if (spt.ErrXYZ()[2] > 0.)
584 fHMCypull->Fill((spt.XYZ()[1] - mcxyz[1]) / std::sqrt(spt.ErrXYZ()[2]));
585 if (spt.ErrXYZ()[5] > 0.)
586 fHMCzpull->Fill((spt.XYZ()[2] - mcxyz[2]) / std::sqrt(spt.ErrXYZ()[5]));
588 catch (cet::exception&) {
bool merge() const noexcept
constexpr to_element_t to_element
process_name opflash particleana ie x
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
geo::WireID WireID() const
CryostatID_t Cryostat
Index of cryostat.
Planes which measure Z direction.
float Integral() const
Integral under the calibrated signal waveform of the hit, in tick x ADC units.
geo::View_t View() const
View for the plane of the hit.
WireID_t Wire
Index of the wire within its plane.
const art::PtrVector< recob::Hit > & getAssociatedHits(const recob::SpacePoint &spt) const
const SpacePointAlg fSptalgSep
float PeakAmplitude() const
The estimated amplitude of the hit at its peak, in ADC units.
void makeMCTruthSpacePoints(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const art::PtrVector< recob::Hit > &hits, std::vector< recob::SpacePoint > &spts) const
double separation(const art::PtrVector< recob::Hit > &hits) const
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
std::string fHitModuleLabel
void makeSpacePoints(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const art::PtrVector< recob::Hit > &hits, std::vector< recob::SpacePoint > &spts) const
PlaneID_t Plane
Index of the plane within its TPC.
double correctedTime(detinfo::DetectorPropertiesData const &detProp, const recob::Hit &hit) const
float PeakTime() const
Time of the signal peak, in tick units.
const SpacePointAlg fSptalgTime
void bookHistograms(bool mc)
float SummedADC() const
The sum of calibrated ADC counts of the hit (0. by default)
const SpacePointAlg fSptalgDefault
constexpr WireID()=default
Default constructor: an invalid TPC ID.
float SigmaPeakTime() const
Uncertainty for the signal peak, in tick units.
2D representation of charge deposited in the TDC/wire plane
TPCID_t TPC
Index of the TPC within its cryostat.
std::string fClusterModuleLabel