LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
ShowerUnidirectiondEdx_tool.cc
Go to the documentation of this file.
1 //############################################################################
2 //### Name: ShowerUnidirectiondEdx ###
3 //### Author: Ed Tyley ###
4 //### Date: 13.05.19 ###
5 //### Description: Tool for finding the dEdx of the start track of the ###
6 //### shower using the standard calomitry module. Derived ###
7 //### from the EMShower_module.cc ###
8 //############################################################################
9 
10 //Framework Includes
12 
13 //LArSoft Includes
18 
19 namespace ShowerRecoTools {
20 
22 
23  public:
25 
26  //Generic Direction Finder
27  int CalculateElement(const art::Ptr<recob::PFParticle>& pfparticle,
28  art::Event& Event,
29  reco::shower::ShowerElementHolder& ShowerEleHolder) override;
30 
31  private:
32  //Define the services and algorithms
35 
36  //fcl parameters.
37  int fVerbose;
39  dEdxTrackLength; //Max length from a hit can be to the start point in cm.
40  bool fMaxHitPlane; //Set the best planes as the one with the most hits
41  bool fMissFirstPoint; //Do not use any hits from the first wire.
42  bool fSumHitSnippets; // Whether to treat hits individually or only one hit per snippet
48  };
49 
51  : IShowerTool(pset.get<fhicl::ParameterSet>("BaseTools"))
52  , fCalorimetryAlg(pset.get<fhicl::ParameterSet>("CalorimetryAlg"))
53  , fVerbose(pset.get<int>("Verbose"))
54  , fdEdxTrackLength(pset.get<float>("dEdxTrackLength"))
55  , fMaxHitPlane(pset.get<bool>("MaxHitPlane"))
56  , fMissFirstPoint(pset.get<bool>("MissFirstPoint"))
57  , fSumHitSnippets(pset.get<bool>("SumHitSnippets"))
58  , fShowerStartPositionInputLabel(pset.get<std::string>("ShowerStartPositionInputLabel"))
59  , fInitialTrackHitsInputLabel(pset.get<std::string>("InitialTrackHitsInputLabel"))
60  , fShowerDirectionInputLabel(pset.get<std::string>("ShowerDirectionInputLabel"))
61  , fShowerdEdxOutputLabel(pset.get<std::string>("ShowerdEdxOutputLabel"))
62  , fShowerBestPlaneOutputLabel(pset.get<std::string>("ShowerBestPlaneOutputLabel"))
63  {}
64 
66  art::Event& Event,
67  reco::shower::ShowerElementHolder& ShowerEleHolder)
68  {
69 
71 
72  // Shower dEdx calculation
73  if (!ShowerEleHolder.CheckElement(fShowerStartPositionInputLabel)) {
74  if (fVerbose)
75  mf::LogError("ShowerUnidirectiondEdx") << "Start position not set, returning " << std::endl;
76  return 1;
77  }
78  if (!ShowerEleHolder.CheckElement(fInitialTrackHitsInputLabel)) {
79  if (fVerbose)
80  mf::LogError("ShowerUnidirectiondEdx")
81  << "Initial Track Hits not set returning" << std::endl;
82  return 1;
83  }
84  if (!ShowerEleHolder.CheckElement(fShowerDirectionInputLabel)) {
85  if (fVerbose)
86  mf::LogError("ShowerUnidirectiondEdx") << "Shower Direction not set" << std::endl;
87  return 1;
88  }
89 
90  //Get the initial track hits
91  std::vector<art::Ptr<recob::Hit>> trackhits;
92  ShowerEleHolder.GetElement(fInitialTrackHitsInputLabel, trackhits);
93 
94  if (trackhits.empty()) {
95  if (fVerbose)
96  mf::LogWarning("ShowerUnidirectiondEdx") << "Not Hits in the initial track" << std::endl;
97  return 0;
98  }
99 
100  geo::Point_t ShowerStartPosition = {-999, -999, -999};
101  ShowerEleHolder.GetElement(fShowerStartPositionInputLabel, ShowerStartPosition);
102 
103  geo::Vector_t showerDir = {-999, -999, -999};
104  ShowerEleHolder.GetElement(fShowerDirectionInputLabel, showerDir);
105 
106  geo::TPCID vtxTPC = fGeom->FindTPCAtPosition(geo::vect::toPoint(ShowerStartPosition));
107 
108  // Split the track hits per plane
109  std::vector<double> dEdxVec;
110  std::vector<std::vector<art::Ptr<recob::Hit>>> trackHits;
111  unsigned int numPlanes = fGeom->Nplanes();
112  trackHits.resize(numPlanes);
113 
114  // Loop over the track hits and split into planes
115  for (unsigned int hitIt = 0; hitIt < trackhits.size(); ++hitIt) {
116  art::Ptr<recob::Hit> hit = trackhits.at(hitIt);
117  geo::PlaneID hitWire = hit->WireID();
118  geo::TPCID TPC = hitWire.asTPCID();
119 
120  //only get hits from the same TPC as the vertex
121  if (TPC == vtxTPC) { (trackHits.at(hitWire.Plane)).push_back(hit); }
122  }
123 
124  int bestHitsPlane = 0;
125  int bestPlaneHits = 0;
126  int bestPlane = -999;
127  double minPitch = 999;
128 
129  auto const clockData =
131  auto const detProp =
133 
134  for (unsigned int plane = 0; plane < numPlanes; ++plane) {
135  std::vector<art::Ptr<recob::Hit>> trackPlaneHits = trackHits.at(plane);
136 
137  std::map<art::Ptr<recob::Hit>, std::vector<art::Ptr<recob::Hit>>> hitSnippets;
138  if (fSumHitSnippets)
139  hitSnippets = IShowerTool::GetLArPandoraShowerAlg().OrganizeHits(trackPlaneHits);
140 
141  if (trackPlaneHits.size()) {
142 
143  double dEdx = -999;
144  double totQ = 0;
145  double avgT = 0;
146  double pitch = 0;
147 
148  //Calculate the pitch
149  double wirepitch = fGeom->WirePitch(trackPlaneHits.at(0)->WireID().planeID());
150  double angleToVert =
152  trackPlaneHits[0]->WireID().planeID()) -
153  0.5 * TMath::Pi();
154  double cosgamma =
155  std::abs(std::sin(angleToVert) * showerDir.Y() + std::cos(angleToVert) * showerDir.Z());
156 
157  pitch = wirepitch / cosgamma;
158 
159  if (pitch) { // Check the pitch is calculated correctly
160  int nhits = 0;
161  std::vector<float> vQ;
162 
163  //Get the first wire
164  int w0 = trackPlaneHits.at(0)->WireID().Wire;
165 
166  for (auto const& hit : trackPlaneHits) {
167 
168  if (fSumHitSnippets && !hitSnippets.count(hit)) continue;
169 
170  // Get the wire for each hit
171  int w1 = hit->WireID().Wire;
172  if (fMissFirstPoint && w0 == w1) { continue; }
173 
174  //Ignore hits that are too far away.
175  if (std::abs((w1 - w0) * pitch) < dEdxTrackLength) {
176 
177  double q = hit->Integral();
178  if (fSumHitSnippets) {
179  for (const art::Ptr<recob::Hit> secondaryHit : hitSnippets[hit])
180  q += secondaryHit->Integral();
181  }
182 
183  vQ.push_back(q);
184  totQ += hit->Integral();
185  avgT += hit->PeakTime();
186  ++nhits;
187  }
188  }
189 
190  if (totQ) {
191  // Check if the pitch is the smallest yet (best plane)
192  if (pitch < minPitch) {
193  minPitch = pitch;
194  bestPlane = plane;
195  }
196 
197  //Get the median and calculate the dEdx using the algorithm.
198  double dQdx = TMath::Median(vQ.size(), &vQ[0]) / pitch;
199  dEdx = fCalorimetryAlg.dEdx_AREA(
200  clockData, detProp, dQdx, avgT / nhits, trackPlaneHits.at(0)->WireID().Plane);
201 
202  if (isinf(dEdx)) { dEdx = -999; };
203 
204  if (nhits > bestPlaneHits || ((nhits == bestPlaneHits) && (pitch < minPitch))) {
205  bestHitsPlane = plane;
206  bestPlaneHits = nhits;
207  }
208  }
209  dEdxVec.push_back(dEdx);
210  }
211  else {
212  throw cet::exception("ShowerUnidirectiondEdx")
213  << "pitch is 0. I can't think how it is 0? Stopping so I can tell you" << std::endl;
214  }
215  }
216  else { // if not (trackPlaneHits.size())
217  dEdxVec.push_back(-999);
218  }
219  trackPlaneHits.clear();
220  } //end loop over planes
221 
222  //TODO
223  std::vector<double> dEdxVecErr = {-999, -999, -999};
224 
225  ShowerEleHolder.SetElement(dEdxVec, dEdxVecErr, fShowerdEdxOutputLabel);
226 
227  //Set The best plane
228  if (fMaxHitPlane) { bestPlane = bestHitsPlane; }
229 
230  if (bestPlane == -999) {
231  throw cet::exception("ShowerUnidirectiondEdx") << "No best plane set";
232  }
233  else {
234  ShowerEleHolder.SetElement(bestPlane, fShowerBestPlaneOutputLabel);
235  }
236 
237  if (fVerbose > 1) {
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;
241  }
242  }
243 
244  return 0;
245  }
246 }
247 
#define DEFINE_ART_CLASS_TOOL(tool)
Definition: ToolMacros.h:42
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
Definition: geo_vectors.h:160
The data type to uniquely identify a Plane.
Definition: geo_types.h:463
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.
ShowerUnidirectiondEdx(const fhicl::ParameterSet &pset)
STL namespace.
geo::WireID const & WireID() const
Initial tdc tick for hit.
Definition: Hit.h:280
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.
Definition: PlaneGeo.h:166
TPCID FindTPCAtPosition(Point_t const &point) const
Returns the ID of the TPC at specified location.
art::ServiceHandle< geo::Geometry > fGeom
parameter set interface
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.
Definition: geo_types.h:381
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:481
int GetElement(const std::string &Name, T &Element) const
int CalculateElement(const art::Ptr< recob::PFParticle > &pfparticle, art::Event &Event, reco::shower::ShowerElementHolder &ShowerEleHolder) override
float dEdx(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const TCSlice &slc, TP3D &tp3d)
Definition: PFPUtils.cxx:2675
Detector simulation of raw signals on wires.
const shower::LArPandoraShowerAlg & GetLArPandoraShowerAlg() const
Definition: IShowerTool.h:82
constexpr TPCID const & asTPCID() const
Conversion to TPCID (for convenience of notation).
Definition: geo_types.h:438
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
Definition: geo_vectors.h:180
decltype(auto) get(T &&obj)
ADL-aware version of std::to_string.
Definition: StdUtils.h:120
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.
Definition: GeometryCore.h:977
Length_t WirePitch(PlaneID const &planeid=plane_zero) const
Returns the distance between two consecutive wires.
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33