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