26 #include "Math/VectorUtil.h" 28 using ROOT::Math::VectorUtil::Angle;
42 void FinddEdxLength(std::vector<double>& dEdx_vec, std::vector<double>& dEdx_val);
115 <<
"Can only correct for SCE if input is already corrected" << std::endl;
130 mf::LogError(
"ShowerTrajPointdEdx") <<
"Start position not set, returning " << std::endl;
136 <<
"Initial Track Spacepoints is not set returning" << std::endl;
145 std::vector<art::Ptr<recob::SpacePoint>> tracksps;
148 if (tracksps.empty()) {
150 mf::LogWarning(
"ShowerTrajPointdEdx") <<
"no spacepointsin the initial track" << std::endl;
172 auto const pfpHandle = Event.getValidHandle<std::vector<recob::PFParticle>>(
fPFParticleLabel);
175 std::vector<art::Ptr<anab::T0>> pfpT0Vec = fmpfpt0.at(pfparticle.
key());
176 if (pfpT0Vec.size() == 1) { pfpT0Time = pfpT0Vec.front()->Time(); }
180 std::map<int, std::vector<double>> dEdx_vec;
181 std::map<int, std::vector<double>> dEdx_vecErr;
182 std::map<int, int> num_hits;
190 auto const clockData =
195 std::map<art::Ptr<recob::Hit>, std::vector<art::Ptr<recob::Hit>>> hitSnippets;
197 std::vector<art::Ptr<recob::Hit>> trackHits;
204 for (
auto const& sp : tracksps) {
207 std::vector<art::Ptr<recob::Hit>>
hits = fmsp.at(sp.key());
211 <<
"no hit for the spacepoint. This suggest the find many is wrong." << std::endl;
223 if (TPC != vtxTPC) {
continue; }
226 auto const pos = sp->position();
227 double dist_from_start = (pos - ShowerStartPosition).R();
236 unsigned int index = 999;
237 double MinDist = 999;
246 auto const dist = (pos - TrajPosition).R();
255 if (index == 999) {
continue; }
261 if ((TrajPosition - TrajPositionStart).R() == 0) {
continue; }
262 if ((TrajPosition - ShowerStartPosition).R() == 0) {
continue; }
264 if ((TrajPosition - TrajPositionStart).R() <
fMinDistCutOff * wirepitch) {
continue; }
272 geo::Vector_t const TrajDirectionYZ{0, TrajDirection.Y(), TrajDirection.Z()};
281 double velocity = detProp.DriftVelocity(detProp.Efield(), detProp.Temperature());
282 double distance_in_x = TrajDirection.X() * (wirepitch / TrajDirection.Dot(PlaneDirection));
283 double time_taken =
std::abs(distance_in_x / velocity);
291 if ((TrajPosition - TrajPositionStart).R() >
dEdxTrackLength) {
continue; }
294 ++num_hits[planeid.
Plane];
297 double trackpitch = (TrajDirection * (wirepitch / TrajDirection.Dot(PlaneDirection))).R();
301 trackpitch, pos, TrajDirection.Unit(), hit->
WireID().
TPC);
308 dQdx += secondaryHit->Integral();
313 double localEField = detProp.Efield();
318 clockData, detProp, dQdx, hit->
PeakTime(), planeid.
Plane, pfpT0Time, localEField);
321 dEdx_vec[planeid.
Plane].push_back(dEdx);
326 int best_plane = -std::numeric_limits<int>::max();
327 for (
auto const& [plane, numHits] : num_hits) {
328 if (
fVerbose > 2) std::cout <<
"Plane: " << plane <<
" with size: " << numHits << std::endl;
329 if (numHits > max_hits) {
335 if (best_plane < 0) {
337 mf::LogError(
"ShowerTrajPointdEdx") <<
"No hits in any plane, returning " << std::endl;
346 std::map<int, std::vector<double>> dEdx_vec_cut;
348 dEdx_vec_cut[plane_id.Plane] = {};
351 for (
auto& dEdx_plane : dEdx_vec) {
352 FinddEdxLength(dEdx_plane.second, dEdx_vec_cut[dEdx_plane.first]);
356 std::vector<double> dEdx_val;
357 std::vector<double> dEdx_valErr;
358 for (
auto const& dEdx_plane : dEdx_vec_cut) {
360 if ((dEdx_plane.second).empty()) {
361 dEdx_val.push_back(-999);
362 dEdx_valErr.push_back(-999);
367 dEdx_val.push_back(TMath::Median((dEdx_plane.second).size(), &(dEdx_plane.second)[0]));
371 double dEdx_mean = 0;
372 for (
auto const&
dEdx : dEdx_plane.second) {
373 if (
dEdx > 10 ||
dEdx < 0) {
continue; }
376 dEdx_val.push_back(dEdx_mean / (
float)(dEdx_plane.second).size());
381 std::cout <<
"#Best Plane: " << best_plane << std::endl;
382 for (
unsigned int plane = 0; plane < dEdx_vec.size(); plane++) {
383 std::cout <<
"#Plane: " << plane <<
" #" << std::endl;
384 std::cout <<
"#Median: " << dEdx_val[plane] <<
" #" << std::endl;
386 for (
auto const&
dEdx : dEdx_vec_cut[plane]) {
387 std::cout <<
"dEdx: " <<
dEdx << std::endl;
401 std::vector<double>& dEdx_val)
411 if (dEdx_vec.size() < 4) {
416 bool upperbound =
false;
419 int upperbound_int = 0;
420 if (dEdx_vec[0] >
fdEdxCut) { ++upperbound_int; }
421 if (dEdx_vec[1] >
fdEdxCut) { ++upperbound_int; }
422 if (dEdx_vec[2] >
fdEdxCut) { ++upperbound_int; }
423 if (upperbound_int > 1) { upperbound =
true; }
425 dEdx_val.push_back(dEdx_vec[0]);
426 dEdx_val.push_back(dEdx_vec[1]);
427 dEdx_val.push_back(dEdx_vec[2]);
429 for (
unsigned int dEdx_iter = 2; dEdx_iter < dEdx_vec.size(); ++dEdx_iter) {
435 double dEdx = dEdx_vec[dEdx_iter];
440 dEdx_val.push_back(dEdx);
441 if (
fVerbose > 1) std::cout <<
"Adding dEdx: " << dEdx << std::endl;
446 if (dEdx_iter < dEdx_vec.size() - 1) {
447 if (dEdx_vec[dEdx_iter + 1] >
fdEdxCut) {
449 std::cout <<
"Next dEdx hit is good removing hit" << dEdx << std::endl;
454 if (dEdx_iter < dEdx_vec.size() - 2) {
455 if (dEdx_vec[dEdx_iter + 2] >
fdEdxCut) {
457 std::cout <<
"Next Next dEdx hit is good removing hit" << dEdx << std::endl;
467 dEdx_val.push_back(dEdx);
468 if (
fVerbose > 1) std::cout <<
"Adding dEdx: " << dEdx << std::endl;
473 if (dEdx_iter < dEdx_vec.size() - 1) {
474 if (dEdx_vec[dEdx_iter + 1] >
fdEdxCut) {
476 std::cout <<
"Next dEdx hit is good removing hit " << dEdx << std::endl;
481 if (dEdx_iter < dEdx_vec.size() - 2) {
482 if (dEdx_vec[dEdx_iter + 2] >
fdEdxCut) {
484 std::cout <<
"Next Next dEdx hit is good removing hit " << dEdx << std::endl;
details::range_type< T > Iterate() const
Initializes the specified ID with the ID of the first cryostat.
double SCECorrectPitch(double const &pitch, geo::Point_t const &pos, geo::Vector_t const &dir, unsigned int const &TPC) const
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
static constexpr Flag_t NoPoint
The trajectory point is not defined.
Point_t const & LocationAtPoint(size_t i) const
Access to track position at different points.
The data type to uniquely identify a Plane.
void SetElement(T &dataproduct, const std::string &Name, bool checktag=false)
size_t NumberTrajectoryPoints() const
Various functions related to the presence and the number of (valid) points.
constexpr auto abs(T v)
Returns the absolute value of the argument.
float Integral() const
Integral under the calibrated signal waveform of the hit, in tick x ADC units.
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.
TPCID FindTPCAtPosition(Point_t const &point) const
Returns the ID of the TPC at specified location.
unsigned int MaxPlanes() const
Returns the largest number of planes among all TPCs in this detector.
key_type key() const noexcept
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
Provides recob::Track data product.
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.
float PeakTime() const
Time of the signal peak, in tick units.
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
decltype(auto) get(T &&obj)
ADL-aware version of std::to_string.
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
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
Vector_t const & GetIncreasingWireDirection() const
Returns the direction of increasing wires.
PointFlags_t const & FlagsAtPoint(size_t i) const
Access to i-th TrajectoryPoint or its Flags.
Vector_t DirectionAtPoint(size_t i) const
Access to track direction at different points.
double SCECorrectEField(double const &EField, geo::Point_t const &pos) const
2D representation of charge deposited in the TDC/wire plane
TPCID_t TPC
Index of the TPC within its cryostat.
Length_t WirePitch(PlaneID const &planeid=plane_zero) const
Returns the distance between two consecutive wires.
const art::FindManyP< T1 > & GetFindManyP(const art::ValidHandle< std::vector< T2 >> &handle, const art::Event &evt, const art::InputTag &moduleTag)
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a "fitted" track:
cet::coded_exception< error, detail::translate > exception