Public Member Functions | Public Attributes | Private Member Functions | Private Attributes | List of all members
dune::DisambigAlgProtoDUNESP Class Reference

#include <DisambigAlgProtoDUNESP.h>

Public Member Functions

 DisambigAlgProtoDUNESP (fhicl::ParameterSet const &pset)
 
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. More...
 

Public Attributes

std::vector< std::pair< art::Ptr< recob::Hit >, geo::WireID > > fDisambigHits
 The final list of hits to pass back to be made. More...
 

Private Member Functions

bool notOuterWire (geo::WireID wireid)
 
int longTPC (int apa)
 

Private Attributes

dune::apa::APAGeometryAlg fAPAGeo
 
double fTimeCut
 
double fDistanceCut
 
double fDcut2
 

Detailed Description

Definition at line 44 of file DisambigAlgProtoDUNESP.h.

Constructor & Destructor Documentation

dune::DisambigAlgProtoDUNESP::DisambigAlgProtoDUNESP ( fhicl::ParameterSet const &  pset)

Definition at line 49 of file DisambigAlgProtoDUNESP.cxx.

50  {
51  this->reconfigure(pset);
52  }
dune::DisambigAlgProtoDUNESP::reconfigure
void reconfigure(fhicl::ParameterSet const &p)
Definition: DisambigAlgProtoDUNESP.cxx:55

Member Function Documentation

int dune::DisambigAlgProtoDUNESP::longTPC ( int  apa)
private

Definition at line 333 of file DisambigAlgProtoDUNESP.cxx.

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  }
bool dune::DisambigAlgProtoDUNESP::notOuterWire ( geo::WireID  wireid)
private

Definition at line 318 of file DisambigAlgProtoDUNESP.cxx.

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  }
result
static QCString result
Definition: fortranscanner.cpp:56614
geo::TPCID::TPC
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:406
void dune::DisambigAlgProtoDUNESP::reconfigure ( fhicl::ParameterSet const &  p)

Definition at line 55 of file DisambigAlgProtoDUNESP.cxx.

56  {
57 
58  fTimeCut = p.get<double>("TimeCut");
59  fDistanceCut = p.get<double>("DistanceCut");
60  fDcut2 = fDistanceCut * fDistanceCut;
61  }
test.p
p
Definition: test.py:223
void dune::DisambigAlgProtoDUNESP::RunDisambig ( detinfo::DetectorPropertiesData const &  detProp,
const std::vector< art::Ptr< recob::Hit > > &  OrigHits 
)

Run disambiguation as currently configured.

Definition at line 66 of file DisambigAlgProtoDUNESP.cxx.

68  {
69  fDisambigHits.clear();
70 
71  art::ServiceHandle<geo::Geometry> geo;
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  {
220  geo::WireIDIntersection isect;
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  {
237  geo::WireIDIntersection isect;
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  }
geo::WireIDIntersection::z
double z
z position of intersection
Definition: geo_types.h:805
View
AdcChannelData::View View
Definition: test_AdcChannelData.cxx:23
geo::kV
Planes which measure V.
Definition: geo_types.h:130
geo::GeometryCore::ChannelToWire
std::vector< geo::WireID > ChannelToWire(raw::ChannelID_t const channel) const
Returns a list of wires connected to the specified TPC channel.
Definition: GeometryCore.cxx:1107
geo::kZ
Planes which measure Z direction.
Definition: geo_types.h:132
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
dune::DisambigAlgProtoDUNESP::fDisambigHits
std::vector< std::pair< art::Ptr< recob::Hit >, geo::WireID > > fDisambigHits
The final list of hits to pass back to be made.
Definition: DisambigAlgProtoDUNESP.h:57
geo::kU
Planes which measure U.
Definition: geo_types.h:129
geo::fhicl::WireID
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
Definition: geo_types_fhicl.h:407
dune::apa::APAGeometryAlg::ChannelToAPA
unsigned int ChannelToAPA(uint32_t chan)
Get number of the APA containing the given channel.
Definition: APAGeometryAlg.cxx:141
geo::GeometryCore::NTPC
unsigned int NTPC(unsigned int cstat=0) const
Returns the total number of TPCs in the specified cryostat.
Definition: GeometryCore.h:2303
dune::DisambigAlgProtoDUNESP::notOuterWire
bool notOuterWire(geo::WireID wireid)
Definition: DisambigAlgProtoDUNESP.cxx:318
geo::GeometryCore::WireIDsIntersect
bool WireIDsIntersect(WireID const &wid1, WireID const &wid2, geo::Point_t &intersection) const
Computes the intersection between two wires.
geo::WireIDIntersection::y
double y
y position of intersection
Definition: geo_types.h:804
dune::DisambigAlgProtoDUNESP::fAPAGeo
dune::apa::APAGeometryAlg fAPAGeo
Definition: DisambigAlgProtoDUNESP.h:67
trkf::Plane
recob::tracking::Plane Plane
Definition: TrackState.h:17
geo
LArSoft geometry interface.
Definition: ChannelGeo.h:16
fhicl::exception
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33

Member Data Documentation

dune::apa::APAGeometryAlg dune::DisambigAlgProtoDUNESP::fAPAGeo
private

Definition at line 67 of file DisambigAlgProtoDUNESP.h.

double dune::DisambigAlgProtoDUNESP::fDcut2
private

Definition at line 70 of file DisambigAlgProtoDUNESP.h.

std::vector< std::pair<art::Ptr<recob::Hit>, geo::WireID> > dune::DisambigAlgProtoDUNESP::fDisambigHits

The final list of hits to pass back to be made.

Definition at line 57 of file DisambigAlgProtoDUNESP.h.

double dune::DisambigAlgProtoDUNESP::fDistanceCut
private

Definition at line 69 of file DisambigAlgProtoDUNESP.h.

double dune::DisambigAlgProtoDUNESP::fTimeCut
private

Definition at line 68 of file DisambigAlgProtoDUNESP.h.

The documentation for this class was generated from the following files:
Generated by   doxygen 1.8.11