DisambigAlgProtoDUNESP.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 //
3 // DisambigAlgProtoDUNESP.cxx
4 //
5 // trj@fnal.gov
6 // tjyang@fnal.gov
7 //
8 // description
9 //
10 // An algorithm for disambiguation for the single-phase protoDUNE detector that is robust against missing the third wire plane.
11 // Assumes hits cannot be found outside the outer APA's. Assumes that the volumes outside of the APA's are in fact TPC's.
12 // Assumes 36-degree DUNE geometry with no more than one intersection of an induction-plane wire with a collection-plane wire.
13 // Since the U and V angles are the same, assume that the a U wire intersects a V wire in at most one place as well.
14 // Assumes there are two induction planes. kU, kV, and one collection plane kZ.
15 //
16 //
17 ////////////////////////////////////////////////////////////////////////
18 
19 
20 //Framework includes:
23 #include "art_root_io/TFileService.h"
24 #include "art_root_io/TFileDirectory.h"
26 
35 #include "DisambigAlgProtoDUNESP.h"
37 
38 #include <map>
39 #include <cmath>
40 #include <vector>
41 #include <iostream>
42 #include <fstream>
43 #include <cstdlib>
44 
45 #include "TH1D.h"
46 
47 namespace dune{
48 
50  {
51  this->reconfigure(pset);
52  }
53 
54  //----------------------------------------------------------
56  {
57 
58  fTimeCut = p.get<double>("TimeCut");
59  fDistanceCut = p.get<double>("DistanceCut");
60  fDcut2 = fDistanceCut * fDistanceCut;
61  }
62 
63 
64  //----------------------------------------------------------
65  //----------------------------------------------------------
67  const std::vector< art::Ptr<recob::Hit> > &OrigHits )
68  {
69  fDisambigHits.clear();
70 
72 
73  size_t napas = geo->NTPC()/2;
74 
75  // do only one APA at a time.
76 
77  for (size_t apa=0; apa<napas; apa++)
78  {
79 
80  int tpc = longTPC(apa); // for purposes of evaluating time offsets for time matching.
81  //tick offsets differ for the even and odd TPC's
82  int cryostat = 0; // protoDUNE-SP has only one cryostat
83 
84  // TH1D *histu = new TH1D("histu","histu",4000,0,4000);
85  // TH1D *histv = new TH1D("histv","histv",4000,0,4000);
86  // TH1D *histz = new TH1D("histz","histz",4000,0,4000);
87 
88  // copy hit ptrs to local storage
89 
90  std::vector<art::Ptr<recob::Hit> > hitsUV[2]; // index 0=U, 1=V
91  std::vector<art::Ptr<recob::Hit> > hitsZ;
92 
93  //std::cout << "Dumping hits for this APA" << std::endl;
94 
95  for (size_t i = 0; i<OrigHits.size(); ++i)
96  {
97  unsigned int hitapa(0), hitcryo(0);
98  fAPAGeo.ChannelToAPA(OrigHits[i]->Channel(), hitapa, hitcryo);
99  if (apa == hitapa)
100  {
101  //std::cout << *OrigHits[i] << std::endl;
102 
103  switch (OrigHits[i]->View())
104  {
105  case geo::kU:
106  hitsUV[0].push_back(OrigHits[i]);
107  // if (OrigHits[i]->WireID().TPC==1)
108  // histu->Fill(OrigHits[i]->PeakTime()
109  // - detProp.GetXTicksOffset(0,
110  // tpc,
111  // cryostat)
112  // ,OrigHits[i]->Charge());
113 
114  break;
115  case geo::kV:
116  hitsUV[1].push_back(OrigHits[i]);
117  // if (OrigHits[i]->WireID().TPC==1)
118  // histv->Fill(OrigHits[i]->PeakTime()
119  // - detProp.GetXTicksOffset(1,
120  // tpc,
121  // cryostat)
122  // ,OrigHits[i]->Charge());
123  break;
124  case geo::kZ:
125  hitsZ.push_back(OrigHits[i]);
126  // if (OrigHits[i]->WireID().TPC==1)
127  // histz->Fill(OrigHits[i]->PeakTime()
128  // - detProp.GetXTicksOffset(2,
129  // tpc,
130  // cryostat)
131  // ,OrigHits[i]->Charge());
132  break;
133  default:
134  throw cet::exception("DisambigAlgProtoDUNESP") <<": hit view unkonwn. \n";
135  }
136  }
137  }
138 
139  // std::cout << " DisambigAlgProtoDUNESP timing means: " <<histu->GetMean()<<" "<<histv->GetMean()<<" "<<histz->GetMean()<<std::endl;
140  // delete histu;
141  // delete histv;
142  // delete histz;
143 
144  size_t zsize = hitsZ.size();
145 
146  // loop over U and V planes to disambiguate.
147  // Identify hits matching in time in the Z and the other induction planes
148 
149  for (size_t uv=0; uv<2; uv++) // uv is the current induction plane
150  {
151  size_t other=1-uv; // the index of the other induction plane
152  size_t uvsize = hitsUV[uv].size(); // this induction plane's hits
153  size_t othersize = hitsUV[other].size(); // the other induction plane's hits
154 
155  for (size_t iuv = 0; iuv < uvsize; ++iuv)
156  {
157  //std::cout << "Attempting to Disambiguate hit: " << *hitsUV[uv][iuv] << std::endl;
158  //std::cout << "Getting time offset: " << hitsUV[uv][iuv]->WireID().Plane << " " << tpc << " " << cryostat << std::endl;
159  //std::cout << "time offset: " << detProp.GetXTicksOffset(hitsUV[uv][iuv]->WireID().Plane,tpc,cryostat) << std::endl;
160 
161  double tuv = hitsUV[uv][iuv]->PeakTime()
162  - detProp.GetXTicksOffset(hitsUV[uv][iuv]->WireID().Plane,tpc,cryostat);
163 
164  std::vector<size_t> zmatches; // list of indices of time-matched hits
165  std::vector<size_t> othermatches;
166  for (size_t z = 0; z < zsize; z++)
167  {
168  double tz = hitsZ[z]->PeakTime()
169  - detProp.GetXTicksOffset(hitsZ[z]->WireID().Plane,tpc,cryostat);
170  //std::cout << "Looking for a Z match: " << tuv << " " << tz << " Offset: " <<
171  // detProp.GetXTicksOffset(hitsZ[z]->WireID().Plane,tpc,cryostat)
172  // << " plane: " << hitsZ[z]->WireID().Plane << " TPC: " << tpc << " Cryo: " << cryostat
173  // << std::endl;
174 
175  if (std::abs(tuv-tz)<fTimeCut)
176  {
177  zmatches.push_back(z);
178  }
179  }
180  for (size_t iother = 0; iother < othersize; iother++)
181  {
182  double tother = hitsUV[other][iother]->PeakTime()
183  - detProp.GetXTicksOffset(hitsUV[other][iother]->WireID().Plane,tpc,cryostat);
184  if (std::abs(tuv-tother)<fTimeCut)
185  {
186  othermatches.push_back(iother);
187  }
188  }
189 
190  //std::cout << "Number of time-matched Z hits: " << zmatches.size() << std::endl;
191  //std::cout << "Number of time-matched Other-Ind View hits: " << othermatches.size() << std::endl;
192 
193  if (zmatches.size() == 0 && othermatches.size() == 0) continue; // hit has no matches -- possibly noise
194 
195  // find out how many of these matched-in-time hits are also matched in space.
196  // take triplets over doublets. Assume we are on the sensitive side of the APA.
197 
198  std::vector< double > zmatchz;
199  std::vector< double > zmatchy;
200  std::vector< double > othermatchz;
201  std::vector< double > othermatchy;
202 
203  std::vector<geo::WireID> wires = geo->ChannelToWire(hitsUV[uv][iuv]->Channel());
204  size_t wsize = wires.size();
205  std::vector<size_t> ndoublets(wsize,0);
206  std::vector<size_t> ntriplets(wsize,0);
207 
208  for (size_t w=0; w<wsize; w++)
209  {
210  ndoublets[w] = 0;
211  ntriplets[w] = 0;
212  geo::WireID uvwire = wires[w];
213  if ( notOuterWire(uvwire) )
214  {
215  for (size_t z=0; z<zmatches.size(); z++)
216  {
217  geo::WireID zwire = geo->ChannelToWire(hitsZ[zmatches[z]]->Channel())[0];
218  if ( notOuterWire(zwire) ) // we really shouldn't have any hits on the outer z wires
219  {
221  if (geo->WireIDsIntersect(zwire,uvwire,isect))
222  {
223  zmatchz.push_back(isect.z);
224  zmatchy.push_back(isect.y);
225  }
226  }
227  }
228  // There is at most one intersection of the u channel with a v channel. Still have to loop over possibilities though.
229  for (size_t iother=0; iother<othermatches.size(); iother++)
230  {
231  std::vector<geo::WireID> otherwires = geo->ChannelToWire(hitsUV[other][othermatches[iother]]->Channel());
232  for (size_t otherw = 0; otherw < otherwires.size(); otherw++)
233  {
234  geo::WireID otherwire = otherwires[otherw];
235  if (notOuterWire(otherwire))
236  {
238  if (geo->WireIDsIntersect(otherwire,uvwire,isect))
239  {
240  othermatchz.push_back(isect.z);
241  othermatchy.push_back(isect.y);
242  }
243  }
244  }
245  }
246 
247  // look for triplets among the matches, and then doublets
248  for (size_t izmatch=0;izmatch<zmatchz.size();izmatch++)
249  {
250  bool found_triplet = false;
251  for (size_t iothermatch=0;iothermatch<othermatchz.size();iothermatch++)
252  {
253  double disz = zmatchz[izmatch]-othermatchz[iothermatch];
254  double disy = zmatchy[izmatch]-othermatchy[iothermatch];
255  if ( disz*disz + disy*disy < fDcut2 )
256  {
257  ntriplets[w]++;
258  found_triplet = true;
259  }
260  }
261  if (!found_triplet) ndoublets[w] ++;
262  }
263  for (size_t iothermatch=0;iothermatch<othermatchz.size();iothermatch++)
264  {
265  bool found_triplet = false;
266  for (size_t izmatch=0;izmatch<zmatchz.size();izmatch++)
267  {
268  double disz = zmatchz[izmatch]-othermatchz[iothermatch];
269  double disy = zmatchy[izmatch]-othermatchy[iothermatch];
270  if ( disz*disz + disy*disy < fDcut2 )
271  {
272  found_triplet = true;
273  }
274  }
275  if (!found_triplet) ndoublets[w] ++;
276  } // end summing up doublets and triplets
277  } // end check if not an outer wire
278  } // end loop on wire candidates for U channels
279 
280  bool found_disambig = false;
281  size_t bestwire = 0;
282  for (size_t w=0; w<wsize; w++)
283  {
284  if (!found_disambig && (ntriplets[w]>0 || ndoublets[w] > 0))
285  {
286  bestwire = w;
287  found_disambig = true;
288  }
289  if (ntriplets[w] > ntriplets[bestwire])
290  {
291  bestwire = w;
292  found_disambig = true;
293  }
294  if ( (ntriplets[w] == ntriplets[bestwire]) && (ndoublets[w] > ndoublets[bestwire]))
295  {
296  bestwire = w;
297  found_disambig = true;
298  }
299  }
300  if (found_disambig)
301  {
302  fDisambigHits.push_back(std::pair<art::Ptr<recob::Hit>, geo::WireID>(hitsUV[uv][iuv],wires[bestwire]));
303  }
304  else
305  {
306  //mf::LogWarning("DisambigAlg35t")<<"Could not find disambiguated hit for "<<*hitsUV[uv][iuv]<<"\n";
307  //std::cout <<"Could not find disambiguated hit for "<<*hitsUV[uv][iuv]<<"\n";
308  }
309 
310  } // end loop over induction hits in this plane
311  } // end loop over induction planes
312  } // end loop over APA's
313  }
314 
315 // method to tell us whether a particular wire is in an outer TPC. Hardcoded outer TPC numbers.
316 // might be quicker to check the lowest two bits: 00 or 11 mean outer TPC.
317 
319  {
320  auto tpc = wireid.TPC;
321  bool result;
322  if (tpc == 0 || tpc == 3 || tpc == 4 || tpc == 7 || tpc == 8 || tpc == 11)
323  {
324  result = false;
325  }
326  else
327  {
328  result = true;
329  }
330  return(result);
331  }
332 
334  {
335  int tpc=0;
336 
337  if (apa == 0) tpc=1;
338  if (apa == 1) tpc=2;
339  if (apa == 2) tpc=5;
340  if (apa == 3) tpc=6;
341  if (apa == 4) tpc=9;
342  if (apa == 5) tpc=10;
343 
344  return(tpc);
345  }
346 }
347 
348 // todo -- look at hits nearby undisambiguated hits and assign similar disambiguation assignments
double z
z position of intersection
Definition: geo_types.h:805
static QCString result
Encapsulate the construction of a single cyostat.
AdcChannelData::View View
Planes which measure V.
Definition: geo_types.h:130
double GetXTicksOffset(int p, int t, int c) const
struct vector vector
std::vector< geo::WireID > ChannelToWire(raw::ChannelID_t const channel) const
Returns a list of wires connected to the specified TPC channel.
Planes which measure Z direction.
Definition: geo_types.h:132
void reconfigure(fhicl::ParameterSet const &p)
void RunDisambig(detinfo::DetectorPropertiesData const &detProp, const std::vector< art::Ptr< recob::Hit > > &OrigHits)
Run disambiguation as currently configured.
T abs(T value)
std::vector< std::pair< art::Ptr< recob::Hit >, geo::WireID > > fDisambigHits
The final list of hits to pass back to be made.
Planes which measure U.
Definition: geo_types.h:129
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
T get(std::string const &key) const
Definition: ParameterSet.h:271
p
Definition: test.py:223
unsigned int ChannelToAPA(uint32_t chan)
Get number of the APA containing the given channel.
Encapsulate the geometry of a wire.
unsigned int NTPC(unsigned int cstat=0) const
Returns the total number of TPCs in the specified cryostat.
bool notOuterWire(geo::WireID wireid)
bool WireIDsIntersect(WireID const &wid1, WireID const &wid2, geo::Point_t &intersection) const
Computes the intersection between two wires.
Declaration of signal hit object.
Encapsulate the construction of a single detector plane.
DisambigAlgProtoDUNESP(fhicl::ParameterSet const &pset)
double y
y position of intersection
Definition: geo_types.h:804
Declaration of basic channel signal object.
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:406
dune::apa::APAGeometryAlg fAPAGeo
recob::tracking::Plane Plane
Definition: TrackState.h:17
LArSoft geometry interface.
Definition: ChannelGeo.h:16
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
Encapsulate the construction of a single detector plane.