28 #include "Math/VectorUtil.h" 30 using ROOT::Math::VectorUtil::Angle;
44 void FinddEdxLength(std::vector<double>& dEdx_vec, std::vector<double>& dEdx_val);
118 <<
"Can only correct for SCE if input is already corrected" << std::endl;
133 mf::LogError(
"ShowerTrajPointdEdx") <<
"Start position not set, returning " << std::endl;
139 <<
"Initial Track Spacepoints is not set returning" << std::endl;
148 std::vector<art::Ptr<recob::SpacePoint>> tracksps;
151 if (tracksps.empty()) {
153 mf::LogWarning(
"ShowerTrajPointdEdx") <<
"no spacepointsin the initial track" << std::endl;
175 auto const pfpHandle = Event.getValidHandle<std::vector<recob::PFParticle>>(
fPFParticleLabel);
178 std::vector<art::Ptr<anab::T0>> pfpT0Vec = fmpfpt0.at(pfparticle.
key());
179 if (pfpT0Vec.size() == 1) { pfpT0Time = pfpT0Vec.front()->Time(); }
183 std::map<int, std::vector<double>> dEdx_vec;
184 std::map<int, std::vector<double>> dEdx_vecErr;
185 std::map<int, int> num_hits;
193 auto const clockData =
198 std::map<art::Ptr<recob::Hit>, std::vector<art::Ptr<recob::Hit>>> hitSnippets;
200 std::vector<art::Ptr<recob::Hit>> trackHits;
207 for (
auto const& sp : tracksps) {
210 std::vector<art::Ptr<recob::Hit>>
hits = fmsp.at(sp.key());
214 <<
"no hit for the spacepoint. This suggest the find many is wrong." << std::endl;
226 if (TPC != vtxTPC) {
continue; }
229 auto const pos = sp->position();
230 double dist_from_start = (pos - ShowerStartPosition).R();
239 unsigned int index = 999;
240 double MinDist = 999;
249 auto const dist = (pos - TrajPosition).R();
258 if (index == 999) {
continue; }
264 if ((TrajPosition - TrajPositionStart).R() == 0) {
continue; }
265 if ((TrajPosition - ShowerStartPosition).R() == 0) {
continue; }
267 if ((TrajPosition - TrajPositionStart).R() <
fMinDistCutOff * wirepitch) {
continue; }
275 geo::Vector_t const TrajDirectionYZ{0, TrajDirection.Y(), TrajDirection.Z()};
284 double velocity = detProp.DriftVelocity(detProp.Efield(), detProp.Temperature());
285 double distance_in_x = TrajDirection.X() * (wirepitch / TrajDirection.Dot(PlaneDirection));
286 double time_taken =
std::abs(distance_in_x / velocity);
294 if ((TrajPosition - TrajPositionStart).R() >
dEdxTrackLength) {
continue; }
297 ++num_hits[planeid.
Plane];
300 double trackpitch = (TrajDirection * (wirepitch / TrajDirection.Dot(PlaneDirection))).R();
304 trackpitch, pos, TrajDirection.Unit(), hit->
WireID().
TPC);
311 dQdx += secondaryHit->Integral();
316 double localEField = detProp.Efield();
321 clockData, detProp, dQdx, hit->
PeakTime(), planeid.
Plane, pfpT0Time, localEField);
324 dEdx_vec[planeid.
Plane].push_back(dEdx);
329 int best_plane = -std::numeric_limits<int>::max();
330 for (
auto const& [plane, numHits] : num_hits) {
331 if (
fVerbose > 2) std::cout <<
"Plane: " << plane <<
" with size: " << numHits << std::endl;
332 if (numHits > max_hits) {
338 if (best_plane < 0) {
340 mf::LogError(
"ShowerTrajPointdEdx") <<
"No hits in any plane, returning " << std::endl;
349 std::map<int, std::vector<double>> dEdx_vec_cut;
351 dEdx_vec_cut[plane_id.Plane] = {};
354 for (
auto& dEdx_plane : dEdx_vec) {
355 FinddEdxLength(dEdx_plane.second, dEdx_vec_cut[dEdx_plane.first]);
359 std::vector<double> dEdx_val;
360 std::vector<double> dEdx_valErr;
361 for (
auto const& dEdx_plane : dEdx_vec_cut) {
363 if ((dEdx_plane.second).empty()) {
364 dEdx_val.push_back(-999);
365 dEdx_valErr.push_back(-999);
370 dEdx_val.push_back(TMath::Median((dEdx_plane.second).size(), &(dEdx_plane.second)[0]));
374 double dEdx_mean = 0;
375 for (
auto const&
dEdx : dEdx_plane.second) {
376 if (
dEdx > 10 ||
dEdx < 0) {
continue; }
379 dEdx_val.push_back(dEdx_mean / (
float)(dEdx_plane.second).size());
384 std::cout <<
"#Best Plane: " << best_plane << std::endl;
385 for (
unsigned int plane = 0; plane < dEdx_vec.size(); plane++) {
386 std::cout <<
"#Plane: " << plane <<
" #" << std::endl;
387 std::cout <<
"#Median: " << dEdx_val[plane] <<
" #" << std::endl;
389 for (
auto const&
dEdx : dEdx_vec_cut[plane]) {
390 std::cout <<
"dEdx: " <<
dEdx << std::endl;
404 std::vector<double>& dEdx_val)
414 if (dEdx_vec.size() < 4) {
419 bool upperbound =
false;
422 int upperbound_int = 0;
423 if (dEdx_vec[0] >
fdEdxCut) { ++upperbound_int; }
424 if (dEdx_vec[1] >
fdEdxCut) { ++upperbound_int; }
425 if (dEdx_vec[2] >
fdEdxCut) { ++upperbound_int; }
426 if (upperbound_int > 1) { upperbound =
true; }
428 dEdx_val.push_back(dEdx_vec[0]);
429 dEdx_val.push_back(dEdx_vec[1]);
430 dEdx_val.push_back(dEdx_vec[2]);
432 for (
unsigned int dEdx_iter = 2; dEdx_iter < dEdx_vec.size(); ++dEdx_iter) {
438 double dEdx = dEdx_vec[dEdx_iter];
443 dEdx_val.push_back(dEdx);
444 if (
fVerbose > 1) std::cout <<
"Adding dEdx: " << dEdx << std::endl;
449 if (dEdx_iter < dEdx_vec.size() - 1) {
450 if (dEdx_vec[dEdx_iter + 1] >
fdEdxCut) {
452 std::cout <<
"Next dEdx hit is good removing hit" << dEdx << std::endl;
457 if (dEdx_iter < dEdx_vec.size() - 2) {
458 if (dEdx_vec[dEdx_iter + 2] >
fdEdxCut) {
460 std::cout <<
"Next Next dEdx hit is good removing hit" << dEdx << std::endl;
470 dEdx_val.push_back(dEdx);
471 if (
fVerbose > 1) std::cout <<
"Adding dEdx: " << dEdx << std::endl;
476 if (dEdx_iter < dEdx_vec.size() - 1) {
477 if (dEdx_vec[dEdx_iter + 1] >
fdEdxCut) {
479 std::cout <<
"Next dEdx hit is good removing hit " << dEdx << std::endl;
484 if (dEdx_iter < dEdx_vec.size() - 2) {
485 if (dEdx_vec[dEdx_iter + 2] >
fdEdxCut) {
487 std::cout <<
"Next Next dEdx hit is good removing hit " << dEdx << std::endl;
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.
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
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
unsigned int MaxPlanes() const
Returns the largest number of planes among all TPCs in this detector.
TPCID FindTPCAtPosition(Point_t const &point) const
Returns the ID of the TPC at specified location.
Interface for a class providing readout channel mapping to geometry.
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
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.
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.
double dist(const TReal *x, const TReal *y, const unsigned int dimension)
Provides recob::Track data product.
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
range_type< T > Iterate() const
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
PlaneGeo const & Plane(TPCID const &tpcid, View_t view) const
Returns the specified wire.
TPCID_t TPC
Index of the TPC within its cryostat.
const art::FindManyP< T1 > & GetFindManyP(const art::ValidHandle< std::vector< T2 >> &handle, const art::Event &evt, const art::InputTag &moduleTag)
art framework interface to geometry description
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
constexpr TPCID const & asTPCID() const
Conversion to PlaneID (for convenience of notation).