LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
DisambigCheater_module.cc
Go to the documentation of this file.
1 // Class: DisambigCheater
3 // Module Type: producer
4 // File: DisambigCheater_module.cc
5 //
6 // tylerdalion@gmail.com
7 //
9 
16 
29 
30 namespace hit {
31 
33  public:
34  explicit DisambigCheater(fhicl::ParameterSet const& p);
35 
36  private:
37  void produce(art::Event& e) override;
38  void endJob() override;
39 
44 
45  std::string fChanHitLabel;
46  std::string fWidHitLabel;
47 
48  std::map<std::pair<double, double>, std::vector<geo::WireID>> fHitToWids;
49  void InitHitToWids(detinfo::DetectorClocksData const& clockData,
50  const std::vector<art::Ptr<recob::Hit>>& ChHits);
51 
53  geo::WireID const& wid,
54  art::Ptr<recob::Wire> const& wire,
56 
57  unsigned int fFalseChanHits =
58  0; // hits with no IDEs - could be noise or other technical problems
59  unsigned int fBadIDENearestWire = 0; // Just to count the inconsistent IDE wire returns
60 
61  std::vector<unsigned int> fMaxWireShift; // shift to account for space charge and diffusion
62  };
63 
64  //-------------------------------------------------------------------
66  {
67  fChanHitLabel = p.get<std::string>("ChanHitLabel");
68 
69  // let HitCollectionCreator declare that we are going to produce
70  // hits and associations with wires and raw digits
71  // (with no particular product label)
73 
74  // Space charge can shift true IDE postiion to far-off channels.
75  // Calculate maximum number of wires to shift hit in order to be on correct channel.
76  // Shift no more than half of the number of channels, as beyond there will
77  // have been a closer wire segment.
78  if (geom->Ncryostats() != 1 || geom->NTPC() < 1) {
79  fMaxWireShift.resize(3);
80  fMaxWireShift[0] = 1;
81  fMaxWireShift[1] = 1;
82  fMaxWireShift[2] = 1;
83  }
84  else {
85  // assume TPC 0 is typical of all in terms of number of channels
86  constexpr geo::TPCID tpcid{0, 0};
87  unsigned int np = wireReadoutGeom->Nplanes(tpcid);
88  fMaxWireShift.resize(np);
89  for (auto const& planeGeo : wireReadoutGeom->Iterate<geo::PlaneGeo>(tpcid)) {
90  unsigned int nw = planeGeo.Nwires();
91  for (unsigned int w = 0; w < nw; ++w) {
92 
93  // for vertical planes
94  if (planeGeo.View() == geo::kZ) {
95  fMaxWireShift[2] = planeGeo.Nwires();
96  break;
97  }
98 
99  auto const xyz = planeGeo.Wire(w).GetCenter();
100  auto const xyz_next = planeGeo.Wire(w + 1).GetCenter();
101 
102  if (xyz.Z() == xyz_next.Z()) {
103  fMaxWireShift[planeGeo.ID().Plane] = w;
104  break;
105  }
106  } // end wire loop
107  } // end plane loop
108 
109  for (unsigned int i = 0; i < np; i++)
110  fMaxWireShift[i] = std::floor(fMaxWireShift[i] / 2);
111  }
112  }
113 
114  //-------------------------------------------------------------------
116  {
117  // get hits on channels
119  evt.getByLabel(fChanHitLabel, ChanHits);
120  std::vector<art::Ptr<recob::Hit>> ChHits;
121  art::fill_ptr_vector(ChHits, ChanHits);
122 
123  // also get the associated wires and raw digits;
124  // we assume they have been created by the same module as the hits
125  const bool doRawDigitAssns = false;
126 
127  art::FindOneP<recob::Wire> ChannelHitWires(ChanHits, evt, fChanHitLabel);
128  const bool doWireAssns = ChannelHitWires.isValid();
129 
130  // this object contains the hit collection
131  // and its associations to wires and raw digits
132  // (if the original objects have them):
133  recob::HitCollectionCreator hits(evt, doWireAssns, doRawDigitAssns);
134 
135  // find the wireIDs each hit is on
136  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
137  InitHitToWids(clockData, ChHits);
138 
139  // make all of the hits
140  for (size_t h = 0; h < ChHits.size(); h++) {
141 
142  // get the objects associated with this hit
144  if (doWireAssns) wire = ChannelHitWires.at(h);
145 
146  // the trivial Z hits
147  if (ChHits[h]->View() == geo::kZ) {
148  MakeDisambigHit(ChHits[h], ChHits[h]->WireID(), wire, hits);
149  continue;
150  }
151 
152  // make U/V hits if any wire IDs are associated
153  // count hits without a wireID
155  std::pair<double, double> ChanTime(ChHits[h]->Channel() * 1., ChHits[h]->PeakTime() * 1.);
156  if (fHitToWids[ChanTime].size() == 1)
157  MakeDisambigHit(ChHits[h], fHitToWids[ChanTime][0], wire, hits);
158  else if (fHitToWids[ChanTime].size() == 0)
159  fFalseChanHits++;
160  else if (fHitToWids[ChanTime].size() > 1)
161  MakeDisambigHit(ChHits[h], fHitToWids[ChanTime][0], wire, hits); // same thing for now
162 
163  } // for
164 
165  // put the hit collection and associations into the event
166  hits.put_into(evt);
167 
168  fHitToWids.clear();
169  }
170 
171  //-------------------------------------------------
173  art::Ptr<recob::Hit> const& original_hit,
174  geo::WireID const& wid,
175  art::Ptr<recob::Wire> const& wire, // art::Ptr<raw::RawDigit> const& rawdigits,
177  {
178 
179  if (!wid.isValid) {
180  mf::LogWarning("InvalidWireID") << "wid is invalid, hit not being made\n";
181  return;
182  }
183 
184  // create a hit copy of the original one, but with a different wire ID
185  recob::HitCreator hit(*original_hit, wid);
186  hcol.emplace_back(hit.move(), wire);
187  }
188 
189  //-------------------------------------------------------------------
191  const std::vector<art::Ptr<recob::Hit>>& ChHits)
192  {
193  unsigned int Ucount(0), Vcount(0);
194  for (size_t h = 0; h < ChHits.size(); h++) {
195  recob::Hit const& chit = *(ChHits[h]);
196  if (ChHits[h]->View() == geo::kZ) continue;
197  if (ChHits[h]->View() == geo::kU)
198  Ucount++;
199  else if (ChHits[h]->View() == geo::kV)
200  Vcount++;
201  std::vector<geo::WireID> cwids = wireReadoutGeom->ChannelToWire(chit.Channel());
202  std::pair<double, double> ChanTime((double)chit.Channel(),
203  (double)chit.PeakTime()); // hit key value
204 
205  // get hit IDEs
206  std::vector<const sim::IDE*> ides;
207  try {
208  ides = bt_serv->HitToSimIDEs_Ps(clockData, chit);
209  }
210  catch (...) {
211  };
212 
213  // catch the hits that have no IDEs
214  bool hasIDEs = !ides.empty();
215  if (hasIDEs) {
216  try {
217  bt_serv->SimIDEsToXYZ(ides);
218  } //What is this supposed to be doing? It get a vector, but never assigns it anywhere. There has to be a better way to do this check
219  catch (...) {
220  hasIDEs = false;
221  }
222  } // if
223  if (!hasIDEs) {
224  mf::LogVerbatim("DisambigCheat")
225  << "Hit on channel " << chit.Channel() << " between " << chit.PeakTimeMinusRMS()
226  << " and " << chit.PeakTimePlusRMS() << " matches no IDEs.";
227  fHitToWids[ChanTime] = std::vector<geo::WireID>();
228  continue;
229  } // if no IDEs
230 
231  // see what wire ID(s) the hit-IDEs are on
232  // if none: hit maps to empty vector
233  // if one: we have a unique hit, vector of size 1
234  // if more than one: make a vector of all wids
235  std::vector<geo::WireID> widsWithIdes;
236  for (size_t i = 0; i < ides.size(); i++) {
237  geo::Point_t const xyzIde{ides[i]->x, ides[i]->y, ides[i]->z};
238 
239  // Occasionally, ide position is not in TPC
241  geo::TPCID tpcID = geom->FindTPCAtPosition(xyzIde);
242  if (!tpcID.isValid) {
243  mf::LogWarning("DisambigCheat")
244  << "IDE at x = " << xyzIde.X() << ", y = " << xyzIde.Y() << ", z = " << xyzIde.Z()
245  << " does not correspond to a TPC.";
246  continue;
247  }
248  unsigned int tpc = tpcID.TPC, cryo = tpcID.Cryostat;
249 
250  // NearestWire seems to be missing some correction that is applied in LArG4, and is
251  // sometimes off by one or two wires. Use it and then correct it given channel
253  geo::WireID IdeWid;
254  try {
255  IdeWid = wireReadoutGeom->Plane({tpcID, cwids[0].Plane}).NearestWireID(xyzIde);
256  }
257  catch (geo::InvalidWireError const& e) { // adopt suggestion if possible
258  if (!e.hasSuggestedWire()) throw;
259  IdeWid = e.suggestedWireID();
260  mf::LogError("DisambigCheat") << "Detected a point out of its wire plane:\n"
261  << e.what() << "\nUsing suggested wire " << IdeWid << "\n";
262  }
263  geo::WireID storethis = IdeWid; // default...
264  bool foundmatch(false);
265  for (size_t w = 0; w < cwids.size(); w++) {
266  if (cwids[w].TPC != tpc || cwids[w].Cryostat != cryo) continue;
267  if ((unsigned int)std::abs((int)(IdeWid.Wire) - (int)(cwids[w].Wire)) <=
268  fMaxWireShift[cwids[0].Plane]) {
269  storethis = cwids[w]; // ...apply correction
270  foundmatch = true;
271  break;
272  }
273  }
274  if (!foundmatch) {
275  mf::LogWarning("DisambigCheat")
276  << "IDE NearestWire return more than 1 off from channel wids: wire " << IdeWid.Wire;
278  }
279 
280  bool alreadyStored(false);
281  for (size_t wid = 0; wid < widsWithIdes.size(); wid++)
282  if (storethis == widsWithIdes[wid]) alreadyStored = true;
283  if (!alreadyStored) widsWithIdes.push_back(storethis);
284  } // end loop through ides from HitToSimIdes
285 
286  fHitToWids[ChanTime] = widsWithIdes;
287 
288  } // end U/V channel hit loop
289 
290  if (fHitToWids.size() != Ucount + Vcount) {
291  mf::LogWarning("DisambigCheat")
292  << "Nhits mismatch: " << fHitToWids.size() << " " << Ucount + Vcount;
293  }
294  }
295 
296  //-------------------------------------------------------------------
298  {
299 
300  if (fFalseChanHits > 0)
301  mf::LogWarning("DisambigCheater")
302  << fFalseChanHits << " hits had no associated IDE or WireIDs";
303  }
304 
306 
307 } // end hit namespace
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
void produce(art::Event &e) override
DisambigCheater(fhicl::ParameterSet const &p)
art::ServiceHandle< cheat::BackTrackerService const > bt_serv
WireID suggestedWireID() const
Returns a better wire ID.
Definition: Exceptions.h:85
void InitHitToWids(detinfo::DetectorClocksData const &clockData, const std::vector< art::Ptr< recob::Hit >> &ChHits)
Planes which measure V.
Definition: geo_types.h:132
unsigned int NTPC(CryostatID const &cryoid=details::cryostat_zero) const
Returns the total number of TPCs in the specified cryostat.
Definition: GeometryCore.h:416
Declaration of signal hit object.
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.cc:6
std::vector< unsigned int > fMaxWireShift
bool isValid
Whether this ID points to a valid element.
Definition: geo_types.h:194
constexpr auto abs(T v)
Returns the absolute value of the argument.
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:195
Planes which measure Z direction.
Definition: geo_types.h:134
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:430
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
unsigned int Ncryostats() const
Returns the number of cryostats in the detector.
Definition: GeometryCore.h:303
static void declare_products(art::ProducesCollector &collector, std::string instance_name="", bool doWireAssns=true, bool doRawDigitAssns=true)
Declares the hit products we are going to fill.
Definition: HitCreator.cxx:261
art::ServiceHandle< geo::Geometry const > geom
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
Class managing the creation of a new recob::Hit object.
Definition: HitCreator.h:87
Planes which measure U.
Definition: geo_types.h:131
Helper functions to create a hit.
Collection of exceptions for Geometry system.
void hits()
Definition: readHits.C:15
TPCID FindTPCAtPosition(Point_t const &point) const
Returns the ID of the TPC at specified location.
A class handling a collection of hits and its associations.
Definition: HitCreator.h:489
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
Interface for a class providing readout channel mapping to geometry.
Geometry information for a single wire plane.The plane is represented in the geometry by a solid whic...
Definition: PlaneGeo.h:67
std::map< std::pair< double, double >, std::vector< geo::WireID > > fHitToWids
void emplace_back(recob::Hit &&hit, art::Ptr< recob::Wire > const &wire=art::Ptr< recob::Wire >(), art::Ptr< raw::RawDigit > const &digits=art::Ptr< raw::RawDigit >())
Adds the specified hit to the data collection.
Definition: HitCreator.cxx:299
The data type to uniquely identify a TPC.
Definition: geo_types.h:306
Definition of data types for geometry description.
float PeakTimeMinusRMS(float sigmas=+1.) const
Returns a time sigmas RMS away from the peak time.
Definition: Hit.h:300
void put_into(art::Event &)
Moves the data into an event.
Definition: HitCreator.h:622
Detector simulation of raw signals on wires.
unsigned int Nplanes(TPCID const &tpcid=details::tpc_zero) const
Returns the total number of planes in the specified TPC.
ProducesCollector & producesCollector() noexcept
geo::WireReadoutGeom const * wireReadoutGeom
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
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:226
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
virtual std::vector< WireID > ChannelToWire(raw::ChannelID_t channel) const =0
Encapsulate the construction of a single detector plane .
Contains all timing reference information for the detector.
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
range_type< T > Iterate() const
Definition: Iterable.h:121
object containing MC truth information necessary for making RawDigits and doing back tracking ...
std::vector< double > SimIDEsToXYZ(std::vector< sim::IDE > const &ides) const
Exception thrown on invalid wire number.
Definition: Exceptions.h:37
Declaration of basic channel signal object.
std::vector< const sim::IDE * > HitToSimIDEs_Ps(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:46
PlaneGeo const & Plane(TPCID const &tpcid, View_t view) const
Returns the specified wire.
float PeakTimePlusRMS(float sigmas=+1.) const
Returns a time sigmas RMS away from the peak time.
Definition: Hit.h:295
TCEvent evt
Definition: DataStructs.cxx:8
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:306
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:315
void MakeDisambigHit(art::Ptr< recob::Hit > const &hit, geo::WireID const &wid, art::Ptr< recob::Wire > const &wire, recob::HitCollectionCreator &hcol)
Float_t e
Definition: plot.C:35
bool hasSuggestedWire() const
Returns whether we known a better wire number.
Definition: Exceptions.h:79
Float_t w
Definition: plot.C:20
recob::Hit && move()
Prepares the constructed hit to be moved away.
Definition: HitCreator.h:338
raw::ChannelID_t Channel() const
ID of the readout channel the hit was extracted from.
Definition: Hit.h:278
art framework interface to geometry description
Encapsulate the construction of a single detector plane .