75 mf::LogError(
"ShowerUnidirectiondEdx") <<
"Start position not set, returning " << std::endl;
81 <<
"Initial Track Hits not set returning" << std::endl;
86 mf::LogError(
"ShowerUnidirectiondEdx") <<
"Shower Direction not set" << std::endl;
91 std::vector<art::Ptr<recob::Hit>> trackhits;
94 if (trackhits.empty()) {
96 mf::LogWarning(
"ShowerUnidirectiondEdx") <<
"Not Hits in the initial track" << std::endl;
109 std::vector<double> dEdxVec;
110 std::vector<std::vector<art::Ptr<recob::Hit>>> trackHits;
112 trackHits.resize(numPlanes);
115 for (
unsigned int hitIt = 0; hitIt < trackhits.size(); ++hitIt) {
121 if (TPC == vtxTPC) { (trackHits.at(hitWire.
Plane)).push_back(hit); }
124 int bestHitsPlane = 0;
125 int bestPlaneHits = 0;
126 int bestPlane = -999;
127 double minPitch = 999;
129 auto const clockData =
134 for (
unsigned int plane = 0; plane < numPlanes; ++plane) {
135 std::vector<art::Ptr<recob::Hit>> trackPlaneHits = trackHits.at(plane);
137 std::map<art::Ptr<recob::Hit>, std::vector<art::Ptr<recob::Hit>>> hitSnippets;
141 if (trackPlaneHits.size()) {
149 double wirepitch =
fGeom->
WirePitch(trackPlaneHits.at(0)->WireID().planeID());
152 trackPlaneHits[0]->WireID().planeID()) -
155 std::abs(std::sin(angleToVert) * showerDir.Y() + std::cos(angleToVert) * showerDir.Z());
157 pitch = wirepitch / cosgamma;
161 std::vector<float> vQ;
164 int w0 = trackPlaneHits.at(0)->WireID().Wire;
166 for (
auto const&
hit : trackPlaneHits) {
171 int w1 =
hit->WireID().Wire;
177 double q =
hit->Integral();
180 q += secondaryHit->Integral();
184 totQ +=
hit->Integral();
185 avgT +=
hit->PeakTime();
192 if (pitch < minPitch) {
198 double dQdx = TMath::Median(vQ.size(), &vQ[0]) / pitch;
200 clockData, detProp, dQdx, avgT / nhits, trackPlaneHits.at(0)->WireID().Plane);
202 if (isinf(dEdx)) { dEdx = -999; };
204 if (nhits > bestPlaneHits || ((nhits == bestPlaneHits) && (pitch < minPitch))) {
205 bestHitsPlane = plane;
206 bestPlaneHits = nhits;
209 dEdxVec.push_back(dEdx);
213 <<
"pitch is 0. I can't think how it is 0? Stopping so I can tell you" << std::endl;
217 dEdxVec.push_back(-999);
219 trackPlaneHits.clear();
223 std::vector<double> dEdxVecErr = {-999, -999, -999};
230 if (bestPlane == -999) {
231 throw cet::exception(
"ShowerUnidirectiondEdx") <<
"No best plane set";
238 std::cout <<
"Best Plane: " << bestPlane << std::endl;
239 for (
unsigned int plane = 0; plane < dEdxVec.size(); plane++) {
240 std::cout <<
"Plane: " << plane <<
" with dEdx: " << dEdxVec[plane] << std::endl;
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
The data type to uniquely identify a Plane.
void SetElement(T &dataproduct, const std::string &Name, bool checktag=false)
::geo::Point_t toPoint(Point const &p)
Convert the specified point into a geo::Point_t.
constexpr auto abs(T v)
Returns the absolute value of the argument.
geo::WireID const & WireID() const
Initial tdc tick for hit.
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
PlaneGeo const & Plane(PlaneID const &planeid) const
Returns the specified wire.
double WireAngleToVertical(View_t view, TPCID const &tpcid) const
Returns the angle of the wires in the specified view from vertical.
View_t View() const
Which coordinate does this plane measure.
TPCID FindTPCAtPosition(Point_t const &point) const
Returns the ID of the TPC at specified location.
std::map< art::Ptr< recob::Hit >, std::vector< art::Ptr< recob::Hit > > > OrganizeHits(const std::vector< art::Ptr< recob::Hit >> &hits) const
bool CheckElement(const std::string &Name) const
The data type to uniquely identify a TPC.
PlaneID_t Plane
Index of the plane within its TPC.
int GetElement(const std::string &Name, T &Element) const
float dEdx(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const TCSlice &slc, TP3D &tp3d)
Detector simulation of raw signals on wires.
constexpr TPCID const & asTPCID() const
Conversion to TPCID (for convenience of notation).
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
decltype(auto) get(T &&obj)
ADL-aware version of std::to_string.
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
double dEdx_AREA(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, recob::Hit const &hit, double pitch, double T0=0) const
unsigned int Nplanes(TPCID const &tpcid=tpc_zero) const
Returns the total number of planes in the specified TPC.
Length_t WirePitch(PlaneID const &planeid=plane_zero) const
Returns the distance between two consecutive wires.
cet::coded_exception< error, detail::translate > exception