LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
Shower2DLinearRegressionTrackHitFinder_tool.cc
Go to the documentation of this file.
1 //############################################################################
2 //### Name: Shower2DLinearRegressionTrackHitFinder ###
3 //### Author: Dominic Barker (dominic.barker@sheffield.ac.uk) ###
4 //### Date: 13.05.19 ###
5 //### Description: Tool for finding the initial shower track using a rms ###
6 //### based method to define when the shower starts to ###
7 //### shower. This methd is derived from the EMShower_module ###
8 //############################################################################
9 
10 //Framework Includes
12 
13 //LArSoft Includes
21 
22 namespace ShowerRecoTools {
23 
25  public:
27 
28  //Calculate the 2D initial track hits
29  int CalculateElement(const art::Ptr<recob::PFParticle>& pfparticle,
30  art::Event& Event,
31  reco::shower::ShowerElementHolder& ShowerEleHolder) override;
32 
33  private:
34  //Function to find the
35  std::vector<art::Ptr<recob::Hit>> FindInitialTrackHits(
36  const detinfo::DetectorPropertiesData& detProp,
38 
39  //Function to perform a weighted regression fit.
40  Int_t WeightedFit(const Int_t n,
41  const Double_t* x,
42  const Double_t* y,
43  const Double_t* w,
44  Double_t* parm);
45 
46  //fcl parameters
47  unsigned int fNfitpass; //Number of time to fit the straight
48  //line the hits.
49  std::vector<unsigned int> fNfithits; //Max number of hits to fit to.
50  std::vector<double> fToler; //Tolerance or each interaction.
51  //Defined as the perpendicualar
52  //distance from the best fit line.
53  bool fApplyChargeWeight; //Apply charge weighting to the fit.
55  int fVerbose;
61  };
62 
64  const fhicl::ParameterSet& pset)
65  : IShowerTool(pset.get<fhicl::ParameterSet>("BaseTools"))
66  , fNfitpass(pset.get<unsigned int>("Nfitpass"))
67  , fNfithits(pset.get<std::vector<unsigned int>>("Nfithits"))
68  , fToler(pset.get<std::vector<double>>("Toler"))
69  , fApplyChargeWeight(pset.get<bool>("ApplyChargeWeight"))
70  , fPFParticleLabel(pset.get<art::InputTag>("PFParticleLabel"))
71  , fVerbose(pset.get<int>("Verbose"))
72  , fHitLabel(pset.get<art::InputTag>("HitsModuleLabel"))
73  , fShowerStartPositionInputLabel(pset.get<std::string>("ShowerStartPositionInputLabel"))
74  , fShowerDirectionInputLabel(pset.get<std::string>("ShowerDirectionInputLabel"))
75  , fInitialTrackHitsOutputLabel(pset.get<std::string>("InitialTrackHitsOutputLabel"))
77  pset.get<std::string>("InitialTrackSpacePointsOutputLabel"))
78  {
79  if (fNfitpass != fNfithits.size() || fNfitpass != fToler.size()) {
81  << "Shower2DLinearRegressionTrackHitFinderEMShower: fNfithits and fToler need to have size "
82  "fNfitpass";
83  }
84  }
85 
87  const art::Ptr<recob::PFParticle>& pfparticle,
88  art::Event& Event,
89  reco::shower::ShowerElementHolder& ShowerEleHolder)
90  {
91 
92  //This is all based on the shower vertex being known. If it is not lets not do the track
93  if (!ShowerEleHolder.CheckElement(fShowerStartPositionInputLabel)) {
94  if (fVerbose)
95  mf::LogError("Shower2DLinearRegressionTrackHitFinder")
96  << "Start position not set, returning " << std::endl;
97  return 1;
98  }
99  if (!ShowerEleHolder.CheckElement(fShowerDirectionInputLabel)) {
100  if (fVerbose)
101  mf::LogError("Shower2DLinearRegressionTrackHitFinder")
102  << "Direction not set, returning " << std::endl;
103  return 1;
104  }
105 
106  geo::Point_t ShowerStartPosition{-999, -999, -999};
107  ShowerEleHolder.GetElement(fShowerStartPositionInputLabel, ShowerStartPosition);
108 
109  geo::Vector_t ShowerDirection = {-999, -999, -999};
110  ShowerEleHolder.GetElement(fShowerDirectionInputLabel, ShowerDirection);
111 
112  // Get the assocated pfParicle vertex PFParticles
113  auto const pfpHandle = Event.getValidHandle<std::vector<recob::PFParticle>>(fPFParticleLabel);
114 
115  //Get the clusters
116  auto const clusHandle = Event.getValidHandle<std::vector<recob::Cluster>>(fPFParticleLabel);
117 
118  const art::FindManyP<recob::Cluster>& fmc =
119  ShowerEleHolder.GetFindManyP<recob::Cluster>(pfpHandle, Event, fPFParticleLabel);
120  std::vector<art::Ptr<recob::Cluster>> clusters = fmc.at(pfparticle.key());
121 
122  if (clusters.size() < 2) {
123  if (fVerbose)
124  mf::LogError("Shower2DLinearRegressionTrackHitFinder")
125  << "Not enough clusters: " << clusters.size() << std::endl;
126  return 1;
127  }
128 
129  //Get the hit association
130  const art::FindManyP<recob::Hit>& fmhc =
131  ShowerEleHolder.GetFindManyP<recob::Hit>(clusHandle, Event, fPFParticleLabel);
132  std::map<geo::PlaneID, std::vector<art::Ptr<recob::Hit>>> plane_clusters;
133  //Loop over the clusters in the plane and get the hits
134  for (auto const& cluster : clusters) {
135 
136  //Get the hits
137  std::vector<art::Ptr<recob::Hit>> hits = fmhc.at(cluster.key());
138 
139  for (auto hit : hits) {
140  geo::WireID wire = hit->WireID();
141  geo::PlaneID plane = wire.asPlaneID();
142  plane_clusters[plane].push_back(hit);
143  }
144 
145  // Was having issues with clusters having hits in multiple planes breaking PMA
146  // So switched to the method above. May want to switch back when using PandoraTrack
147  //plane_clusters[plane].insert(plane_clusters[plane].end(),hits.begin(),hits.end());
148  }
149 
150  auto const clockData =
152  auto const detProp =
154 
155  std::vector<art::Ptr<recob::Hit>> InitialTrackHits;
156  //Loop over the clusters and order the hits and get the initial track hits in that plane
157  for (auto const& cluster : plane_clusters) {
158 
159  //Get the hits
160  std::vector<art::Ptr<recob::Hit>> hits = cluster.second;
161 
162  //Order the hits
164  detProp, hits, ShowerStartPosition, ShowerDirection);
165 
166  //Find the initial track hits
167  std::vector<art::Ptr<recob::Hit>> trackhits = FindInitialTrackHits(detProp, hits);
168 
169  InitialTrackHits.insert(InitialTrackHits.end(), trackhits.begin(), trackhits.end());
170  }
171 
172  //Holders for the initial track values.
173  ShowerEleHolder.SetElement(InitialTrackHits, fInitialTrackHitsOutputLabel);
174 
175  //Get the associated spacepoints
176  //Get the hits
177  auto const hitHandle = Event.getValidHandle<std::vector<recob::Hit>>(fHitLabel);
178 
179  //get the sp<->hit association
181  ShowerEleHolder.GetFindManyP<recob::SpacePoint>(hitHandle, Event, fPFParticleLabel);
182 
183  //Get the spacepoints associated to the track hit
184  std::vector<art::Ptr<recob::SpacePoint>> intitaltrack_sp;
185  for (auto const& hit : InitialTrackHits) {
186  std::vector<art::Ptr<recob::SpacePoint>> sps = fmsp.at(hit.key());
187  for (auto const sp : sps) {
188  intitaltrack_sp.push_back(sp);
189  }
190  }
191  ShowerEleHolder.SetElement(intitaltrack_sp, fInitialTrackSpacePointsOutputLabel);
192  return 0;
193  }
194 
195  //Function to calculate the what are the initial tracks hits. Adapted from EMShower FindInitialTrackHits
197  const detinfo::DetectorPropertiesData& detProp,
199  {
200 
201  std::vector<art::Ptr<recob::Hit>> trackHits;
202 
203  double parm[2];
204  int fitok = 0;
205  std::vector<double> wfit;
206  std::vector<double> tfit;
207  std::vector<double> cfit;
208 
209  for (size_t i = 0; i < fNfitpass; ++i) {
210 
211  // Fit a straight line through hits
212  unsigned int nhits = 0;
213  for (auto& hit : hits) {
214 
215  //Not sure I am a fan of doing things in wire tick space. What if id doesn't not iterate properly or the
216  //two planes in each TPC are not symmetric.
218 
219  if (i == 0 ||
220  (std::abs((coord.Y() - (parm[0] + coord.X() * parm[1])) *
221  std::cos(std::atan(parm[1]))) < fToler[i - 1]) ||
222  fitok == 1) {
223  ++nhits;
224  if (nhits == fNfithits[i] + 1) break;
225  wfit.push_back(coord.X());
226  tfit.push_back(coord.Y());
227 
228  if (fApplyChargeWeight) { cfit.push_back(hit->Integral()); }
229  else {
230  cfit.push_back(1.);
231  };
232  if (i == fNfitpass - 1) { trackHits.push_back(hit); }
233  }
234  }
235 
236  if (i < fNfitpass - 1 && wfit.size()) {
237  fitok = WeightedFit(wfit.size(), &wfit[0], &tfit[0], &cfit[0], &parm[0]);
238  }
239 
240  wfit.clear();
241  tfit.clear();
242  cfit.clear();
243  }
244  return trackHits;
245  }
246 
247  //Stolen from EMShowerAlg, a linear regression fitting function
249  const Double_t* x,
250  const Double_t* y,
251  const Double_t* w,
252  Double_t* parm)
253  {
254 
255  Double_t sumx = 0.;
256  Double_t sumx2 = 0.;
257  Double_t sumy = 0.;
258  Double_t sumxy = 0.;
259  Double_t sumw = 0.;
260  Double_t eparm[2];
261 
262  parm[0] = 0.;
263  parm[1] = 0.;
264  eparm[0] = 0.;
265  eparm[1] = 0.;
266 
267  for (Int_t i = 0; i < n; i++) {
268  sumx += x[i] * w[i];
269  sumx2 += x[i] * x[i] * w[i];
270  sumy += y[i] * w[i];
271  sumxy += x[i] * y[i] * w[i];
272  sumw += w[i];
273  }
274 
275  if (sumx2 * sumw - sumx * sumx == 0.) return 1;
276  if (sumx2 - sumx * sumx / sumw == 0.) return 1;
277 
278  parm[0] = (sumy * sumx2 - sumx * sumxy) / (sumx2 * sumw - sumx * sumx);
279  parm[1] = (sumxy - sumx * sumy / sumw) / (sumx2 - sumx * sumx / sumw);
280 
281  eparm[0] = sumx2 * (sumx2 * sumw - sumx * sumx);
282  eparm[1] = (sumx2 - sumx * sumx / sumw);
283 
284  if (eparm[0] < 0. || eparm[1] < 0.) return 1;
285 
286  eparm[0] = sqrt(eparm[0]) / (sumx2 * sumw - sumx * sumx);
287  eparm[1] = sqrt(eparm[1]) / (sumx2 - sumx * sumx / sumw);
288 
289  return 0;
290  }
291 }
292 
Float_t x
Definition: compare.C:6
#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
auto coord(Vector &v, unsigned int n) noexcept
Returns an object to manage the coordinate n of a vector.
Declaration of signal hit object.
Float_t y
Definition: compare.C:6
The data type to uniquely identify a Plane.
Definition: geo_types.h:463
void SetElement(T &dataproduct, const std::string &Name, bool checktag=false)
constexpr auto abs(T v)
Returns the absolute value of the argument.
STL namespace.
Set of hits with a 2D structure.
Definition: Cluster.h:69
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
Cluster finding and building.
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
Int_t WeightedFit(const Int_t n, const Double_t *x, const Double_t *y, const Double_t *w, Double_t *parm)
void hits()
Definition: readHits.C:15
int CalculateElement(const art::Ptr< recob::PFParticle > &pfparticle, art::Event &Event, reco::shower::ShowerElementHolder &ShowerEleHolder) override
parameter set interface
key_type key() const noexcept
Definition: Ptr.h:166
std::vector< art::Ptr< recob::Hit > > FindInitialTrackHits(const detinfo::DetectorPropertiesData &detProp, std::vector< art::Ptr< recob::Hit >> &hits)
bool CheckElement(const std::string &Name) const
constexpr PlaneID const & asPlaneID() const
Conversion to PlaneID (for convenience of notation).
Definition: geo_types.h:520
int GetElement(const std::string &Name, T &Element) const
Declaration of cluster object.
void OrderShowerHits(detinfo::DetectorPropertiesData const &detProp, std::vector< art::Ptr< recob::Hit >> &hits, geo::Point_t const &ShowerPosition, geo::Vector_t const &ShowerDirection) const
Detector simulation of raw signals on wires.
const shower::LArPandoraShowerAlg & GetLArPandoraShowerAlg() const
Definition: IShowerTool.h:82
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
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
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
decltype(auto) get(T &&obj)
ADL-aware version of std::to_string.
Definition: StdUtils.h:120
TVector2 HitCoordinates(detinfo::DetectorPropertiesData const &detProp, art::Ptr< recob::Hit > const &hit) const
Definition: MVAAlg.h:12
Char_t n[5]
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:46
Float_t w
Definition: plot.C:20
const art::FindManyP< T1 > & GetFindManyP(const art::ValidHandle< std::vector< T2 >> &handle, const art::Event &evt, const art::InputTag &moduleTag)