162 auto const det_prop =
164 auto const*
sce = lar::providerFrom<spacecharge::SpaceChargeService>();
167 std::vector<art::Ptr<recob::Track>> tracklist;
178 size_t nplanes = geom->
Nplanes();
181 std::unique_ptr<std::vector<anab::Calorimetry>> calorimetrycol(
182 new std::vector<anab::Calorimetry>);
183 std::unique_ptr<art::Assns<recob::Track, anab::Calorimetry>> assn(
187 art::FindManyP<recob::Hit, recob::TrackHitMeta> fmthm(
191 art::FindManyP<anab::T0> fmt0(trackListHandle, evt,
fT0ModuleLabel);
193 for (
size_t trkIter = 0; trkIter < tracklist.size(); ++trkIter) {
195 decltype(
auto) larEnd = tracklist[trkIter]->Trajectory().
End();
202 unsigned int cstat = 0;
203 unsigned int tpc = 0;
204 unsigned int wire = 0;
205 unsigned int plane = 0;
210 if (fmt0.isValid()) {
211 std::vector<art::Ptr<anab::T0>> allT0 = fmt0.at(trkIter);
212 if (allT0.size()) T0 = allT0[0]->Time();
216 std::vector<std::vector<unsigned int>> hits(nplanes);
219 for (
size_t ah = 0; ah < allHits.size(); ++ah) {
220 hits[allHits[ah]->WireID().Plane].push_back(ah);
223 for (
size_t ipl = 0; ipl < nplanes; ++ipl) {
240 float Trk_Length = 0.;
241 std::vector<float> vdEdx;
242 std::vector<float> vresRange;
243 std::vector<float> vdQdx;
244 std::vector<float> deadwire;
245 std::vector<TVector3> vXYZ;
248 if (hits[ipl].
size() < 2) {
249 if (hits[ipl].
size() == 1) {
251 <<
"Only one hit in plane " << ipl <<
" associated with track id " << trkIter;
267 unsigned int wire0 = 100000;
268 unsigned int wire1 = 0;
278 std::vector<double> spdelta;
280 std::vector<double> ChargeBeg;
281 std::stack<double> ChargeEnd;
284 double fTrkPitch = 0;
285 for (
size_t itp = 0; itp < tracklist[trkIter]->NumberTrajectoryPoints(); ++itp) {
287 const auto&
pos = tracklist[trkIter]->LocationAtPoint(itp);
288 const auto&
dir = tracklist[trkIter]->DirectionAtPoint(itp);
290 const double Position[3] = {
pos.X(),
pos.Y(),
pos.Z()};
300 if (
sce->EnableCalSpatialSCE() &&
fSCE)
302 if (
sce->EnableCalSpatialSCE() &&
fSCE)
303 dirOffsets =
sce->GetCalPosOffsets(
geo::Point_t{pos.X() + fTrkPitch * dir.X(),
304 pos.Y() + fTrkPitch * dir.Y(),
305 pos.Z() + fTrkPitch * dir.Z()},
307 TVector3 dir_corr = {fTrkPitch *
dir.X() - dirOffsets.X() + posOffsets.X(),
308 fTrkPitch *
dir.Y() + dirOffsets.Y() - posOffsets.Y(),
309 fTrkPitch *
dir.Z() + dirOffsets.Z() - posOffsets.Z()};
311 fTrkPitch = dir_corr.Mag();
315 <<
"caught exception " << e <<
"\n setting pitch (C) to " <<
util::kBogusD;
323 double xx = 0., yy = 0., zz = 0.;
326 std::vector<double> trkx;
327 std::vector<double> trky;
328 std::vector<double> trkz;
329 std::vector<double> trkw;
330 std::vector<double> trkx0;
331 for (
size_t i = 0; i < hits[ipl].size(); ++i) {
333 std::vector<art::Ptr<recob::SpacePoint>> sptv = fmspts.at(hits[ipl][i]);
334 for (
size_t j = 0; j < sptv.size(); ++j) {
336 double t = allHits[hits[ipl][i]]->PeakTime() -
338 double x = det_prop.ConvertTicksToX(t,
341 allHits[hits[ipl][i]]->
WireID().Cryostat);
342 double w = allHits[hits[ipl][i]]->WireID().Wire;
344 trkx.push_back(sptv[j]->XYZ()[0] -
345 det_prop.ConvertTicksToX(TickT0,
346 allHits[hits[ipl][i]]->WireID().Plane,
347 allHits[hits[ipl][i]]->WireID().TPC,
348 allHits[hits[ipl][i]]->WireID().Cryostat));
351 trkx.push_back(sptv[j]->XYZ()[0]);
353 trky.push_back(sptv[j]->XYZ()[1]);
354 trkz.push_back(sptv[j]->XYZ()[2]);
359 for (
size_t ihit = 0; ihit < hits[ipl].size();
363 plane = allHits[hits[ipl][ihit]]->WireID().Plane;
364 tpc = allHits[hits[ipl][ihit]]->WireID().TPC;
365 cstat = allHits[hits[ipl][ihit]]->WireID().Cryostat;
368 planeID.
Plane = plane;
372 wire = allHits[hits[ipl][ihit]]->WireID().Wire;
373 time = allHits[hits[ipl][ihit]]->PeakTime();
374 stime = allHits[hits[ipl][ihit]]->PeakTimeMinusRMS();
375 etime = allHits[hits[ipl][ihit]]->PeakTimePlusRMS();
376 const size_t& hitIndex = allHits[hits[ipl][ihit]].key();
378 double charge = allHits[hits[ipl][ihit]]->PeakAmplitude();
379 if (
fUseArea) charge = allHits[hits[ipl][ihit]]->Integral();
384 bool fBadhit =
false;
385 if (fmthm.isValid()) {
386 auto vhit = fmthm.at(trkIter);
387 auto vmeta = fmthm.data(trkIter);
388 for (
size_t ii = 0; ii < vhit.size(); ++ii) {
389 if (vhit[ii].
key() == allHits[hits[ipl][ihit]].key()) {
390 if (vmeta[ii]->
Index() == int_max_as_unsigned_int) {
394 if (vmeta[ii]->
Index() >= tracklist[trkIter]->NumberTrajectoryPoints()) {
396 <<
"Requested track trajectory index " << vmeta[ii]->Index()
397 <<
" exceeds the total number of trajectory points " 398 << tracklist[trkIter]->NumberTrajectoryPoints() <<
" for track index " << trkIter
399 <<
". Something is wrong with the track reconstruction. Please contact " 402 if (!tracklist[trkIter]->HasValidPoint(vmeta[ii]->
Index())) {
414 if (
sce->EnableCalSpatialSCE() &&
fSCE)
415 locOffsets =
sce->GetCalPosOffsets(loc, vhit[ii]->WireID().TPC);
416 xyz3d[0] = loc.X() - locOffsets.X();
417 xyz3d[1] = loc.Y() + locOffsets.Y();
418 xyz3d[2] = loc.Z() + locOffsets.Z();
422 vhit[ii]->
WireID().Cryostat) -
423 0.5 * ::util::pi<>();
426 std::abs(std::sin(angleToVert) * dir.Y() + std::cos(angleToVert) * dir.Z());
427 if (cosgamma) { pitch = geom->
WirePitch(vhit[ii]->
View()) / cosgamma; }
434 if (
sce->EnableCalSpatialSCE() &&
fSCE)
435 dirOffsets =
sce->GetCalPosOffsets(
geo::Point_t{loc.X() + pitch * dir.X(),
436 loc.Y() + pitch * dir.Y(),
437 loc.Z() + pitch * dir.Z()},
438 vhit[ii]->WireID().TPC);
439 const TVector3& dir_corr = {pitch * dir.X() - dirOffsets.X() + locOffsets.X(),
440 pitch * dir.Y() + dirOffsets.Y() - locOffsets.Y(),
441 pitch * dir.Z() + dirOffsets.Z() - locOffsets.Z()};
443 pitch = dir_corr.Mag();
451 allHits[hits[ipl][ihit]],
461 if (fBadhit)
continue;
463 if (pitch <= 0) pitch = fTrkPitch;
464 if (!pitch)
continue;
470 spdelta.push_back(0);
473 double dx = xyz3d[0] - xx;
474 double dy = xyz3d[1] - yy;
475 double dz = xyz3d[2] - zz;
476 spdelta.push_back(sqrt(dx * dx + dy * dy + dz * dz));
477 Trk_Length += spdelta.back();
483 ChargeBeg.push_back(charge);
484 ChargeEnd.push(charge);
486 double MIPs = charge;
487 double dQdx = MIPs / pitch;
490 dEdx =
caloAlg.
dEdx_AREA(clock_data, det_prop, *allHits[hits[ipl][ihit]], pitch, T0);
492 dEdx =
caloAlg.
dEdx_AMP(clock_data, det_prop, *allHits[hits[ipl][ihit]], pitch, T0);
494 Kin_En = Kin_En + dEdx * pitch;
496 if (allHits[hits[ipl][ihit]]->
WireID().
Wire < wire0)
497 wire0 = allHits[hits[ipl][ihit]]->WireID().Wire;
498 if (allHits[hits[ipl][ihit]]->
WireID().
Wire > wire1)
499 wire1 = allHits[hits[ipl][ihit]]->WireID().Wire;
501 fMIPs.push_back(MIPs);
502 fdEdx.push_back(dEdx);
503 fdQdx.push_back(dQdx);
504 fwire.push_back(wire);
505 ftime.push_back(time);
509 TVector3 v(xyz3d[0], xyz3d[1], xyz3d[2]);
532 for (
int isp = 0; isp <
fnsps; ++isp) {
534 USChg += ChargeBeg[isp];
537 while (!ChargeEnd.empty()) {
538 if (countsp > 3)
break;
539 DSChg += ChargeEnd.top();
545 GoingDS = (DSChg > USChg);
550 TVector3 track_start(tracklist[trkIter]->Trajectory().
Vertex().
X(),
551 tracklist[trkIter]->Trajectory().
Vertex().
Y(),
552 tracklist[trkIter]->Trajectory().
Vertex().
Z());
553 TVector3 track_end(tracklist[trkIter]->Trajectory().
End().
X(),
554 tracklist[trkIter]->Trajectory().
End().
Y(),
555 tracklist[trkIter]->Trajectory().
End().
Z());
557 if ((
fXYZ[0] - track_start).Mag() + (
fXYZ.back() - track_end).Mag() <
558 (
fXYZ[0] - track_end).Mag() + (
fXYZ.back() - track_start).Mag()) {
569 if (
fResRng.size() < 2 || spdelta.size() < 2) {
571 <<
"fResrng.size() = " <<
fResRng.size() <<
" spdelta.size() = " << spdelta.size();
574 fResRng[fnsps - 1] = spdelta[fnsps - 1] / 2;
575 for (
int isp = fnsps - 2; isp > -1; isp--) {
581 for (
int isp = 1; isp <
fnsps; isp++) {
586 MF_LOG_DEBUG(
"CaloPrtHit") <<
" pt wire time ResRng MIPs pitch dE/dx Ai X Y Z\n";
589 for (
int i = 0; i <
fnsps; ++i) {
590 vresRange.push_back(
fResRng[i]);
591 vdEdx.push_back(
fdEdx[i]);
592 vdQdx.push_back(
fdQdx[i]);
593 vXYZ.push_back(
fXYZ[i]);
594 if (i != 0 && i != fnsps - 1) {
604 << std::setiosflags(std::ios::fixed | std::ios::showpoint) <<
std::setprecision(2)
610 if (nPIDA > 0) { PIDA = PIDA / (double)nPIDA; }
615 << fTrkPitch <<
" nhits= " << fnsps <<
"\n" 616 << std::setiosflags(std::ios::fixed | std::ios::showpoint)
619 <<
" PIDA= " << PIDA <<
"\n";
622 for (
unsigned int iw = wire0; iw < wire1 + 1; ++iw) {
623 plane = allHits[hits[ipl][0]]->WireID().Plane;
624 tpc = allHits[hits[ipl][0]]->WireID().TPC;
625 cstat = allHits[hits[ipl][0]]->WireID().Cryostat;
627 if (channelStatus.
IsBad(channel)) {
628 MF_LOG_DEBUG(
"Calorimetry") <<
"Found dead wire at Plane = " << plane <<
" Wire =" << iw;
629 unsigned int closestwire = 0;
630 unsigned int endwire = 0;
631 unsigned int dwire = 100000;
632 double mindis = 100000;
633 double goodresrange = 0;
634 for (
size_t ihit = 0; ihit < hits[ipl].size(); ++ihit) {
635 channel = allHits[hits[ipl][ihit]]->Channel();
636 if (channelStatus.
IsBad(channel))
continue;
638 std::vector<art::Ptr<recob::SpacePoint>> sppv = fmspts.at(hits[ipl][ihit]);
639 if (sppv.size() < 1)
continue;
643 sppv[0]->XYZ()[0], sppv[0]->XYZ()[1], sppv[0]->XYZ()[2]};
644 double dis1 = (larEnd - xyz).Mag2();
645 if (dis1) dis1 = std::sqrt(dis1);
647 endwire = allHits[hits[ipl][ihit]]->WireID().Wire;
651 closestwire = allHits[hits[ipl][ihit]]->WireID().Wire;
658 deadwire.push_back(goodresrange + (
int(closestwire) -
int(iw)) * fTrkPitch);
661 deadwire.push_back(goodresrange + (
int(iw) -
int(closestwire)) * fTrkPitch);
virtual bool IsBad(raw::ChannelID_t channel) const =0
Returns whether the specified channel is bad in the current run.
std::vector< double > fdEdx
std::optional< double > fNotOnTrackZcut
PlaneGeo const & Plane(unsigned int const p, unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified wire.
Reconstruction base classes.
AdcChannelData::View View
The data type to uniquely identify a Plane.
bool isValid
Whether this ID points to a valid element.
CryostatID_t Cryostat
Index of cryostat.
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
std::vector< TVector3 > fXYZ
geo::Length_t WirePitch(geo::PlaneID const &planeid) const
Returns the distance between two consecutive wires.
Q_EXPORT QTSManip setprecision(int p)
std::vector< double > fResRng
geo::TPCID FindTPCAtPosition(double const worldLoc[3]) const
Returns the ID of the TPC at specified location.
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
View_t View() const
Which coordinate does this plane measure.
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
double dEdx(float dqdx, float Efield)
unsigned int Nplanes(unsigned int tpc=0, unsigned int cstat=0) const
Returns the total number of wire planes in the specified TPC.
std::string fT0ModuleLabel
double dEdx_AMP(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, recob::Hit const &hit, double pitch, double T0=0) const
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
std::vector< size_t > fHitIndex
std::string fSpacePointModuleLabel
std::vector< double > ftime
std::vector< double > fetime
ProductID put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
bool CreateAssn(PRODUCER const &prod, art::Event &evt, std::vector< T > const &a, art::Ptr< U > const &b, art::Assns< U, T > &assn, std::string a_instance, size_t indx=UINT_MAX)
Creates a single one-to-one association.
std::vector< Point_t > convertCollToPoint(std::vector< Point > const &coll)
Class providing information about the quality of channels.
The data type to uniquely identify a TPC.
PlaneID_t Plane
Index of the plane within its TPC.
std::vector< double > fMIPs
void GetPitch(detinfo::DetectorPropertiesData const &det_prop, art::Ptr< recob::Hit > const &hit, std::vector< double > const &trkx, std::vector< double > const &trky, std::vector< double > const &trkz, std::vector< double > const &trkw, std::vector< double > const &trkx0, double *xyz3d, double &pitch, double TickT0)
constexpr auto absDiff(A const &a, B const &b)
Returns the absolute value of the difference between two values.
Q_EXPORT QTSManip setw(int w)
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
raw::ChannelID_t PlaneWireToChannel(WireID const &wireid) const
Returns the ID of the TPC channel connected to the specified wire.
std::vector< double > fdQdx
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
tracking::Point_t Point_t
double TrackPitchInView(recob::Track const &track, geo::View_t view, size_t trajectory_point=0U)
Returns the projected length of track on a wire pitch step [cm].
double dEdx_AREA(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, recob::Hit const &hit, double pitch, double T0=0) const
std::string fTrackModuleLabel
boost::graph_traits< ModuleGraph >::vertex_descriptor Vertex
constexpr double kBogusD
obviously bogus double value
std::vector< double > fstime
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
TPCID_t TPC
Index of the TPC within its cryostat.
std::vector< float > fpitch
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
recob::tracking::Plane Plane
double WireAngleToVertical(geo::View_t view, geo::TPCID const &tpcid) const
Returns the angle of the wires in the specified view from vertical.
cet::coded_exception< error, detail::translate > exception